From 8da1fd1aff10fa4deb424b4393c7d5b55a5f3b23 Mon Sep 17 00:00:00 2001 From: Steven Boeing Date: Wed, 10 May 2023 08:22:51 +0100 Subject: [PATCH 01/17] Initial work on new thermodynamics --- convention.txt | 8 +-- models/3d/moist_3d.f90 | 2 +- src/3d/fields/field_netcdf.f90 | 76 ++++++++++++++++------------- src/3d/fields/fields.f90 | 21 ++++---- src/3d/parcels/parcel_container.f90 | 44 ++++++++++------- src/3d/parcels/parcel_init.f90 | 27 +++++----- src/3d/parcels/parcel_merge.f90 | 25 ++++++---- src/3d/parcels/parcel_netcdf.f90 | 55 +++++++++++++-------- 8 files changed, 152 insertions(+), 106 deletions(-) diff --git a/convention.txt b/convention.txt index 1ea29ad9f..d79d67cd2 100644 --- a/convention.txt +++ b/convention.txt @@ -4,8 +4,8 @@ volg -- volume velog -- velocity velgradg -- velocity gradient tensor vortg -- vorticity -tbuoyg -- total buoyancy -dbuoyg -- dry buoyancy -humg -- specific humidity -humlig -- condensed humidity +tbuoyg -- total buoyancy (in m/s^2) +thetag -- potential temperature +qvg -- water vapour specific humidity +qlg -- liquid water specific humidity nparg -- number of parcels per grid box diff --git a/models/3d/moist_3d.f90 b/models/3d/moist_3d.f90 index 4e2c18a27..fceeb88b5 100644 --- a/models/3d/moist_3d.f90 +++ b/models/3d/moist_3d.f90 @@ -30,7 +30,7 @@ module moist_3d double precision :: e_values(3) ! To create asymmetry, we vary the buoyancy in the plume ! according to b = b_pl*[1 + (e1*x*y+e2*x*z+e3*yz)/R^2]. double precision :: r_smooth_frac ! Fraction of radius where smooth transition starts - + end type plume_type type(plume_type) :: moist diff --git a/src/3d/fields/field_netcdf.f90 b/src/3d/fields/field_netcdf.f90 index dbd763bb3..290c09610 100644 --- a/src/3d/fields/field_netcdf.f90 +++ b/src/3d/fields/field_netcdf.f90 @@ -23,7 +23,7 @@ module field_netcdf integer :: x_vel_id, y_vel_id, z_vel_id, & x_vor_id, y_vor_id, z_vor_id, & - tbuoy_id, vol_id, n_writes + tbuoy_id, theta_id, vol_id, n_writes #ifdef ENABLE_DIAGNOSE integer :: x_vtend_id, y_vtend_id, z_vtend_id, & @@ -31,7 +31,7 @@ module field_netcdf #endif #ifndef ENABLE_DRY_MODE - integer :: dbuoy_id, hum_id, lbuoy_id + integer :: qv_id, ql_id #endif double precision :: restart_time @@ -41,7 +41,7 @@ module field_netcdf coord_ids, t_axis_id, & x_vel_id, y_vel_id, z_vel_id, & x_vor_id, y_vor_id, z_vor_id, & - tbuoy_id, vol_id, & + tbuoy_id, theta_id, vol_id, & n_writes, restart_time #ifdef ENABLE_DIAGNOSE @@ -50,7 +50,7 @@ module field_netcdf #endif #ifndef ENABLE_DRY_MODE - private :: dbuoy_id, hum_id, lbuoy_id + private :: qv_id, ql_id #endif contains @@ -222,33 +222,33 @@ subroutine create_netcdf_field_file(basename, overwrite, l_restart) dimids=dimids, & varid=tbuoy_id) -#ifndef ENABLE_DRY_MODE call define_netcdf_dataset(ncid=ncid, & - name='dry_buoyancy', & - long_name='dry buoyancy', & + name='theta', & + long_name='potential temperature', & std_name='', & - unit='m/s^2', & + unit='K', & dtype=NF90_DOUBLE, & dimids=dimids, & - varid=dbuoy_id) + varid=theta_id) +#ifndef ENABLE_DRY_MODE call define_netcdf_dataset(ncid=ncid, & - name='humidity', & - long_name='specific humidity', & + name='qv', & + long_name='water vapour spec. hum.', & std_name='', & unit='kg/kg', & dtype=NF90_DOUBLE, & dimids=dimids, & - varid=hum_id) + varid=qv_id) call define_netcdf_dataset(ncid=ncid, & - name='liquid_water_content', & - long_name='liquid-water content', & + name='ql', & + long_name='liquid water spec. hum.', & std_name='', & - unit='1', & + unit='kg/kg', & dtype=NF90_DOUBLE, & dimids=dimids, & - varid=lbuoy_id) + varid=ql_id) #endif call define_netcdf_dataset(ncid=ncid, & @@ -312,13 +312,14 @@ subroutine read_netcdf_field_content call get_var_id(ncid, 'buoyancy', tbuoy_id) -#ifndef ENABLE_DRY_MODE - call get_var_id(ncid, 'dry_buoyancy', dbuoy_id) + call get_var_id(ncid, 'theta', theta_id) - call get_var_id(ncid, 'humidity', hum_id) +#ifndef ENABLE_DRY_MODE + call get_var_id(ncid, 'qv', qv_id) - call get_var_id(ncid, 'liquid_water_content', lbuoy_id) + call get_var_id(ncid, 'ql', ql_id) #endif + call get_var_id(ncid, 'volume', vol_id) end subroutine read_netcdf_field_content @@ -392,15 +393,14 @@ subroutine write_netcdf_fields(t) call write_netcdf_dataset(ncid, tbuoy_id, tbuoyg(lo(3):hi(3), lo(2):hi(2), lo(1):hi(1)), & start, cnt) -#ifndef ENABLE_DRY_MODE - call write_netcdf_dataset(ncid, dbuoy_id, dbuoyg(lo(3):hi(3), lo(2):hi(2), lo(1):hi(1)), & + call write_netcdf_dataset(ncid, theta_id, theta(lo(3):hi(3), lo(2):hi(2), lo(1):hi(1)), & start, cnt) - call write_netcdf_dataset(ncid, lbuoy_id, glati * (tbuoyg(lo(3):hi(3), lo(2):hi(2), lo(1):hi(1)) & - - dbuoyg(lo(3):hi(3), lo(2):hi(2), lo(1):hi(1))), & +#ifndef ENABLE_DRY_MODE + call write_netcdf_dataset(ncid, qv_id, qvg(lo(3):hi(3), lo(2):hi(2), lo(1):hi(1)), & start, cnt) - call write_netcdf_dataset(ncid, hum_id, humg(lo(3):hi(3), lo(2):hi(2), lo(1):hi(1)), & + call write_netcdf_dataset(ncid, ql_id, qlg(lo(3):hi(3), lo(2):hi(2), lo(1):hi(1)), & start, cnt) #endif @@ -472,10 +472,10 @@ subroutine read_netcdf_fields(fname, step) cnt) endif - if (has_dataset(lid, 'buoyancy')) then + if (has_dataset(lid, 'theta')) then call read_netcdf_dataset(lid, & - 'buoyancy', & - tbuoyg(box%lo(3):box%hi(3), & + 'theta', & + thetag(box%lo(3):box%hi(3), & box%lo(2):box%hi(2), & box%lo(1):box%hi(1)), & start, & @@ -483,12 +483,22 @@ subroutine read_netcdf_fields(fname, step) endif #ifndef ENABLE_DRY_MODE - if (has_dataset(lid, 'humidity')) then + if (has_dataset(lid, 'qv')) then + call read_netcdf_dataset(lid, & + 'qv', & + qvg(box%lo(3):box%hi(3), & + box%lo(2):box%hi(2), & + box%lo(1):box%hi(1)), & + start, & + cnt) + endif + + if (has_dataset(lid, 'ql')) then call read_netcdf_dataset(lid, & - 'humidity', & - humg(box%lo(3):box%hi(3), & - box%lo(2):box%hi(2), & - box%lo(1):box%hi(1)), & + 'ql', & + qlg(box%lo(3):box%hi(3), & + box%lo(2):box%hi(2), & + box%lo(1):box%hi(1)), & start, & cnt) endif diff --git a/src/3d/fields/fields.f90 b/src/3d/fields/fields.f90 index f307d90d5..f54ba2f3a 100644 --- a/src/3d/fields/fields.f90 +++ b/src/3d/fields/fields.f90 @@ -36,8 +36,9 @@ module fields double precision, allocatable, dimension(:, :, :) :: & #ifndef ENABLE_DRY_MODE - dbuoyg, & ! dry buoyancy (or liquid-water buoyancy) - humg, & ! humidity + thetag, & ! dry buoyancy (or liquid-water buoyancy) + qvg, & ! humidity + qlg, & ! liquid water #endif tbuoyg, & ! buoyancy #ifndef NDEBUG @@ -87,10 +88,10 @@ subroutine field_alloc allocate(vtend(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1), n_dim)) allocate(tbuoyg(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1))) - + allocate(thetag(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1))) #ifndef ENABLE_DRY_MODE - allocate(dbuoyg(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1))) - allocate(humg(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1))) + allocate(qvg(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1))) + allocate(qlg(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1))) #endif allocate(nparg(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1))) @@ -110,9 +111,10 @@ subroutine field_default vortg = zero vtend = zero tbuoyg = zero + thetag = zero #ifndef ENABLE_DRY_MODE - dbuoyg = zero - humg = zero + qvg = zero + qlg = zero #endif nparg = zero nsparg = zero @@ -133,9 +135,10 @@ subroutine field_dealloc deallocate(vortg) deallocate(vtend) deallocate(tbuoyg) + deallocate(thetag) #ifndef ENABLE_DRY_MODE - deallocate(dbuoyg) - deallocate(humg ) + deallocate(qvg) + deallocate(qlg) #endif deallocate(nparg) deallocate(nsparg) diff --git a/src/3d/parcels/parcel_container.f90 b/src/3d/parcels/parcel_container.f90 index 21440b400..3ed3e2894 100644 --- a/src/3d/parcels/parcel_container.f90 +++ b/src/3d/parcels/parcel_container.f90 @@ -27,9 +27,10 @@ module parcel_container IDX_B22, & ! B22 shape matrix element IDX_B23, & ! B23 shape matrix element IDX_VOL, & ! volume - IDX_BUO, & ! buoyancy + IDX_THETA, & ! buoyancy #ifndef ENABLE_DRY_MODE - IDX_HUM, & ! humidity + IDX_QV, & ! water vapour specific humidity + IDX_QL, & ! liquid vapour specific humidity #endif IDX_RK4_X_DPOS, & ! RK4 variable delta x-position IDX_RK4_Y_DPOS, & ! RK4 variable delta y-position @@ -63,9 +64,10 @@ module parcel_container double precision, allocatable, dimension(:) :: & volume, & #ifndef ENABLE_DRY_MODE - humidity, & + qv, & + ql, & #endif - buoyancy + theta ! LS-RK4 variables @@ -100,12 +102,13 @@ subroutine set_buffer_indices IDX_B22 = 10 ! B22 shape matrix element IDX_B23 = 11 ! B23 shape matrix element IDX_VOL = 12 ! volume - IDX_BUO = 13 ! buoyancy + IDX_THETA = 13 ! potential temperature - i = IDX_BUO + 1 + i = IDX_THETA + 1 #ifndef ENABLE_DRY_MODE - IDX_HUM = i - i = i + 1 + IDX_QV = i + IDX_QL = i + 1 + i = i + 2 #endif ! LS-RK4 variables @@ -230,9 +233,10 @@ subroutine parcel_replace(n, m) parcels%position(:, n) = parcels%position(:, m) parcels%vorticity(:, n) = parcels%vorticity(:, m) parcels%volume(n) = parcels%volume(m) - parcels%buoyancy(n) = parcels%buoyancy(m) + parcels%theta(n) = parcels%theta(m) #ifndef ENABLE_DRY_MODE - parcels%humidity(n) = parcels%humidity(m) + parcels%qv(n) = parcels%qv(m) + parcels%ql(n) = parcels%ql(m) #endif parcels%B(:, n) = parcels%B(:, m) @@ -257,9 +261,10 @@ subroutine parcel_alloc(num) allocate(parcels%vorticity(3, num)) allocate(parcels%B(5, num)) allocate(parcels%volume(num)) - allocate(parcels%buoyancy(num)) + allocate(parcels%theta(num)) #ifndef ENABLE_DRY_MODE - allocate(parcels%humidity(num)) + allocate(parcels%qv(num)) + allocate(parcels%ql(num)) #endif call parcel_ellipsoid_allocate(num) @@ -286,9 +291,10 @@ subroutine parcel_dealloc deallocate(parcels%vorticity) deallocate(parcels%B) deallocate(parcels%volume) - deallocate(parcels%buoyancy) + deallocate(parcels%theta) #ifndef ENABLE_DRY_MODE - deallocate(parcels%humidity) + deallocate(parcels%qv) + deallocate(parcels%ql) #endif call parcel_ellipsoid_deallocate @@ -311,9 +317,10 @@ subroutine parcel_serialize(n, buffer) buffer(IDX_X_VOR:IDX_Z_VOR) = parcels%vorticity(:, n) buffer(IDX_B11:IDX_B23) = parcels%B(:, n) buffer(IDX_VOL) = parcels%volume(n) - buffer(IDX_BUO) = parcels%buoyancy(n) + buffer(IDX_THETA) = parcels%theta(n) #ifndef ENABLE_DRY_MODE - buffer(IDX_HUM) = parcels%humidity(n) + buffer(IDX_QV) = parcels%qv(n) + buffer(IDX_QL) = parcels%ql(n) #endif ! LS-RK4 variables: buffer(IDX_RK4_X_DPOS:IDX_RK4_Z_DPOS) = parcels%delta_pos(:, n) @@ -333,9 +340,10 @@ subroutine parcel_deserialize(n, buffer) parcels%vorticity(:, n) = buffer(IDX_X_VOR:IDX_Z_VOR) parcels%B(:, n) = buffer(IDX_B11:IDX_B23) parcels%volume(n) = buffer(IDX_VOL) - parcels%buoyancy(n) = buffer(IDX_BUO) + parcels%theta(n) = buffer(IDX_THETA) #ifndef ENABLE_DRY_MODE - parcels%humidity(n) = buffer(IDX_HUM) + parcels%qv(n) = buffer(IDX_QV) + parcels%ql(n) = buffer(IDX_QL) #endif ! LS-RK4 variables: parcels%delta_pos(:, n) = buffer(IDX_RK4_X_DPOS:IDX_RK4_Z_DPOS) diff --git a/src/3d/parcels/parcel_init.f90 b/src/3d/parcels/parcel_init.f90 index ef30266a4..abe667750 100644 --- a/src/3d/parcels/parcel_init.f90 +++ b/src/3d/parcels/parcel_init.f90 @@ -99,9 +99,10 @@ subroutine parcel_default !$omp do private(n) do n = 1, n_parcels parcels%vorticity(:, n) = zero - parcels%buoyancy(n) = zero + parcels%theta(n) = zero #ifndef ENABLE_DRY_MODE - parcels%humidity(n) = zero + parcels%qv(n) = zero + parcels%ql(n) = zero #endif enddo !$omp end do @@ -190,11 +191,13 @@ subroutine init_parcels_from_grids parcels%vorticity(l, n) = parcels%vorticity(l, n) & + sum(weights * vortg(ks:ks+1, js:js+1, is:is+1, l)) enddo - parcels%buoyancy(n) = parcels%buoyancy(n) & - + sum(weights * tbuoyg(ks:ks+1, js:js+1, is:is+1)) + parcels%theta(n) = parcels%theta(n) & + + sum(weights * thetag(ks:ks+1, js:js+1, is:is+1)) #ifndef ENABLE_DRY_MODE - parcels%humidity(n) = parcels%humidity(n) & - + sum(weights * humg(ks:ks+1, js:js+1, is:is+1)) + parcels%qv(n) = parcels%qv(n) & + + sum(weights * qvg(ks:ks+1, js:js+1, is:is+1)) + parcels%ql(n) = parcels%ql(n) & + + sum(weights * qlg(ks:ks+1, js:js+1, is:is+1)) #endif enddo !$omp end do @@ -208,7 +211,7 @@ end subroutine init_parcels_from_grids subroutine init_fill_halo #ifndef ENABLE_DRY_MODE - integer, parameter :: n_fields = 5 + integer, parameter :: n_fields = 6 #else integer, parameter :: n_fields = 4 #endif @@ -217,9 +220,10 @@ subroutine init_fill_halo call field_interior_to_buffer(vortg(:, :, :, I_X), 1) call field_interior_to_buffer(vortg(:, :, :, I_Y), 2) call field_interior_to_buffer(vortg(:, :, :, I_Z), 3) - call field_interior_to_buffer(tbuoyg, 4) + call field_interior_to_buffer(thetag, 4) #ifndef ENABLE_DRY_MODE - call field_interior_to_buffer(humg, 5) + call field_interior_to_buffer(qvg, 5) + call field_interior_to_buffer(qlg, 6) #endif call interior_to_halo_communication @@ -227,9 +231,10 @@ subroutine init_fill_halo call field_buffer_to_halo(vortg(:, :, :, I_X), 1, .false.) call field_buffer_to_halo(vortg(:, :, :, I_Y), 2, .false.) call field_buffer_to_halo(vortg(:, :, :, I_Z), 3, .false.) - call field_buffer_to_halo(tbuoyg, 4, .false.) + call field_buffer_to_halo(thetag, 4, .false.) #ifndef ENABLE_DRY_MODE - call field_buffer_to_halo(humg, 5, .false.) + call field_buffer_to_halo(qvg, 5, .false.) + call field_buffer_to_halo(qlg, 6, .false.) #endif call field_mpi_dealloc diff --git a/src/3d/parcels/parcel_merge.f90 b/src/3d/parcels/parcel_merge.f90 index 7fc2f7bae..af6584a4d 100644 --- a/src/3d/parcels/parcel_merge.f90 +++ b/src/3d/parcels/parcel_merge.f90 @@ -106,9 +106,10 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, Bm, vm) double precision :: x0(n_merge), y0(n_merge) double precision :: posm(3, n_merge) double precision :: delx, vmerge, dely, delz, B33, mu - double precision :: buoym(n_merge), vortm(3, n_merge) + double precision :: thetam(n_merge), vortm(3, n_merge) #ifndef ENABLE_DRY_MODE - double precision :: hum(n_merge) + double precision :: qvm(n_merge) + double precision :: qlm(n_merge) #endif double precision, intent(out) :: Bm(6, n_merge) ! B11, B12, B13, B22, B23, B33 double precision, intent(out) :: vm(n_merge) @@ -148,9 +149,10 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, Bm, vm) posm(3, l) = parcels%volume(ic) * parcels%position(3, ic) ! buoyancy and humidity - buoym(l) = parcels%volume(ic) * parcels%buoyancy(ic) + thetam(l) = parcels%volume(ic) * parcels%theta(ic) #ifndef ENABLE_DRY_MODE - hum(l) = parcels%volume(ic) * parcels%humidity(ic) + qvm(l) = parcels%volume(ic) * parcels%qv(ic) + qlm(l) = parcels%volume(ic) * parcels%ql(ic) #endif vortm(:, l) = parcels%volume(ic) * parcels%vorticity(:, ic) @@ -175,9 +177,10 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, Bm, vm) posm(3, n) = posm(3, n) + parcels%volume(is) * parcels%position(3, is) ! Accumulate buoyancy and humidity - buoym(n) = buoym(n) + parcels%volume(is) * parcels%buoyancy(is) + thetam(n) = thetam(n) + parcels%volume(is) * parcels%theta(is) #ifndef ENABLE_DRY_MODE - hum(n) = hum(n) + parcels%volume(is) * parcels%humidity(is) + qvm(n) = qvm(n) + parcels%volume(is) * parcels%qv(is) + qlm(n) = qlm(n) + parcels%volume(is) * parcels%ql(is) #endif vortm(:, n) = vortm(:, n) + parcels%volume(is) * parcels%vorticity(:, is) enddo @@ -204,9 +207,10 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, Bm, vm) call apply_periodic_bc(posm(:, m)) ! buoyancy and humidity - buoym(m) = vmerge * buoym(m) + thetam(m) = vmerge * thetam(m) #ifndef ENABLE_DRY_MODE - hum(m) = vmerge * hum(m) + qvm(m) = vmerge * qvm(m) + qlm(m) = vmerge * qlm(m) #endif vortm(:, m) = vmerge * vortm(:, m) enddo @@ -244,9 +248,10 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, Bm, vm) parcels%position(2, ic) = posm(2, l) parcels%position(3, ic) = posm(3, l) - parcels%buoyancy(ic) = buoym(l) + parcels%theta(ic) = thetam(l) #ifndef ENABLE_DRY_MODE - parcels%humidity(ic) = hum(l) + parcels%qv(ic) = qvm(l) + parcels%ql(ic) = qlm(l) #endif parcels%vorticity(:, ic) = vortm(:, l) diff --git a/src/3d/parcels/parcel_netcdf.f90 b/src/3d/parcels/parcel_netcdf.f90 index cae365b41..e7e830766 100644 --- a/src/3d/parcels/parcel_netcdf.f90 +++ b/src/3d/parcels/parcel_netcdf.f90 @@ -25,7 +25,7 @@ module parcel_netcdf character(len=512) :: ncfname integer :: ncid - integer :: npar_dim_id, vol_id, buo_id, & + integer :: npar_dim_id, vol_id, theta_id, & x_pos_id, y_pos_id, z_pos_id, & x_vor_id, y_vor_id, z_vor_id, & b11_id, b12_id, b13_id, & @@ -34,19 +34,19 @@ module parcel_netcdf double precision :: restart_time #ifndef ENABLE_DRY_MODE - integer :: hum_id + integer :: qv_id, ql_id #endif private :: ncid, ncfname, n_writes, npar_dim_id, & x_pos_id, y_pos_id, z_pos_id, start_id, & x_vor_id, y_vor_id, z_vor_id, & b11_id, b12_id, b13_id, b22_id, b23_id, & - vol_id, buo_id, t_dim_id, t_axis_id, & + vol_id, theta_id, t_dim_id, t_axis_id, & restart_time, mpi_dim_id, & read_chunk #ifndef ENABLE_DRY_MODE - private :: hum_id + private :: qv_id, ql_id #endif private :: ncbasename @@ -238,23 +238,32 @@ subroutine create_netcdf_parcel_file(basename, overwrite, l_restart) varid=z_vor_id) call define_netcdf_dataset(ncid=ncid, & - name='buoyancy', & - long_name='parcel buoyancy', & + name='theta', & + long_name='parcel potential temperature',& std_name='', & - unit='m/s^2', & + unit='K', & dtype=NF90_DOUBLE, & dimids=dimids, & - varid=buo_id) + varid=theta_id) #ifndef ENABLE_DRY_MODE call define_netcdf_dataset(ncid=ncid, & - name='humidity', & - long_name='parcel humidity', & + name='qv', & + long_name='water vapour spec. hum.', & std_name='', & - unit='1', & + unit='kg/kg', & dtype=NF90_DOUBLE, & dimids=dimids, & - varid=hum_id) + varid=qv_id) + + call define_netcdf_dataset(ncid=ncid, & + name='ql', & + long_name='liquid water spec. hum.', & + std_name='', & + unit='kg/kg', & + dtype=NF90_DOUBLE, & + dimids=dimids, & + varid=ql_id) #endif call close_definition(ncid) @@ -330,10 +339,11 @@ subroutine write_netcdf_parcels(t) call write_netcdf_dataset(ncid, y_vor_id, parcels%vorticity(2, 1:n_parcels), start, cnt) call write_netcdf_dataset(ncid, z_vor_id, parcels%vorticity(3, 1:n_parcels), start, cnt) - call write_netcdf_dataset(ncid, buo_id, parcels%buoyancy(1:n_parcels), start, cnt) + call write_netcdf_dataset(ncid, theta_id, parcels%theta(1:n_parcels), start, cnt) #ifndef ENABLE_DRY_MODE - call write_netcdf_dataset(ncid, hum_id, parcels%humidity(1:n_parcels), start, cnt) + call write_netcdf_dataset(ncid, qv_id, parcels%qv(1:n_parcels), start, cnt) + call write_netcdf_dataset(ncid, ql_id, parcels%ql(1:n_parcels), start, cnt) #endif ! increment counter n_writes = n_writes + 1 @@ -566,17 +576,22 @@ subroutine read_chunk(first, last, pfirst) parcels%vorticity(3, pfirst:plast), start, cnt) endif - if (has_dataset(ncid, 'buoyancy')) then + if (has_dataset(ncid, 'theta')) then l_valid = .true. - call read_netcdf_dataset(ncid, 'buoyancy', & - parcels%buoyancy(pfirst:plast), start, cnt) + call read_netcdf_dataset(ncid, 'theta', & + parcels%theta(pfirst:plast), start, cnt) endif #ifndef ENABLE_DRY_MODE - if (has_dataset(ncid, 'humidity')) then + if (has_dataset(ncid, 'qv')) then + l_valid = .true. + call read_netcdf_dataset(ncid, 'qv', & + parcels%qv(pfirst:plast), start, cnt) + endif + if (has_dataset(ncid, 'ql')) then l_valid = .true. - call read_netcdf_dataset(ncid, 'humidity', & - parcels%humidity(pfirst:plast), start, cnt) + call read_netcdf_dataset(ncid, 'ql', & + parcels%ql(pfirst:plast), start, cnt) endif #endif From 1d788d115bc50e7e22aa21c958fa20e75a32d0e3 Mon Sep 17 00:00:00 2001 From: Steven Boeing Date: Thu, 11 May 2023 07:45:18 +0100 Subject: [PATCH 02/17] Ongoing work on thermodynamics --- mpi-tests/parcel_merge_serial.f90 | 25 ++++--- src/3d/fields/field_diagnostics.f90 | 2 + src/3d/fields/field_diagnostics_netcdf.f90 | 1 + src/3d/fields/field_netcdf.f90 | 4 ++ src/3d/parcels/parcel_container.f90 | 2 +- src/3d/parcels/parcel_netcdf.f90 | 4 +- src/3d/parcels/parcel_split.f90 | 9 +-- unit-tests/3d/test_ellipsoid_split.f90 | 15 ++-- unit-tests/3d/test_mpi_nearest_1.f90 | 6 +- unit-tests/3d/test_mpi_nearest_10.f90 | 6 +- unit-tests/3d/test_mpi_nearest_11.f90 | 6 +- unit-tests/3d/test_mpi_nearest_12.f90 | 48 ++++++------- unit-tests/3d/test_mpi_nearest_13.f90 | 72 +++++++++---------- unit-tests/3d/test_mpi_nearest_14.f90 | 48 ++++++------- unit-tests/3d/test_mpi_nearest_15.f90 | 56 +++++++-------- unit-tests/3d/test_mpi_nearest_16.f90 | 56 +++++++-------- unit-tests/3d/test_mpi_nearest_17.f90 | 56 +++++++-------- unit-tests/3d/test_mpi_nearest_18.f90 | 48 ++++++------- unit-tests/3d/test_mpi_nearest_19.f90 | 56 +++++++-------- unit-tests/3d/test_mpi_nearest_2.f90 | 6 +- unit-tests/3d/test_mpi_nearest_20.f90 | 24 +++---- unit-tests/3d/test_mpi_nearest_3.f90 | 6 +- unit-tests/3d/test_mpi_nearest_4.f90 | 6 +- unit-tests/3d/test_mpi_nearest_5.f90 | 6 +- unit-tests/3d/test_mpi_nearest_6.f90 | 6 +- unit-tests/3d/test_mpi_nearest_7.f90 | 6 +- unit-tests/3d/test_mpi_nearest_8.f90 | 6 +- unit-tests/3d/test_mpi_nearest_9.f90 | 6 +- unit-tests/3d/test_mpi_parcel_communicate.f90 | 2 +- unit-tests/3d/test_mpi_parcel_delete.f90 | 10 +-- unit-tests/3d/test_mpi_parcel_diagnostics.f90 | 2 +- unit-tests/3d/test_mpi_parcel_pack.f90 | 7 +- unit-tests/3d/test_mpi_parcel_read.f90 | 15 ++-- .../3d/test_mpi_parcel_read_rejection.f90 | 49 ++++++++----- unit-tests/3d/test_mpi_parcel_split.f90 | 2 +- unit-tests/3d/test_mpi_parcel_unpack.f90 | 15 ++-- unit-tests/3d/test_mpi_parcel_write.f90 | 5 +- 37 files changed, 369 insertions(+), 330 deletions(-) diff --git a/mpi-tests/parcel_merge_serial.f90 b/mpi-tests/parcel_merge_serial.f90 index 8f500b9e1..ad1302cea 100644 --- a/mpi-tests/parcel_merge_serial.f90 +++ b/mpi-tests/parcel_merge_serial.f90 @@ -87,9 +87,10 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, Bm, vm) double precision :: x0(n_merge), y0(n_merge) double precision :: posm(3, n_merge) double precision :: delx, vmerge, dely, delz, B33, mu - double precision :: buoym(n_merge), vortm(3, n_merge) + double precision :: thetam(n_merge), vortm(3, n_merge) #ifndef ENABLE_DRY_MODE - double precision :: hum(n_merge) + double precision :: qvm(n_merge) + double precision :: qlm(n_merge) #endif double precision, intent(out) :: Bm(6, n_merge) ! B11, B12, B13, B22, B23, B33 double precision, intent(out) :: vm(n_merge) @@ -124,9 +125,10 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, Bm, vm) posm(3, l) = parcels%volume(ic) * parcels%position(3, ic) ! buoyancy and humidity - buoym(l) = parcels%volume(ic) * parcels%buoyancy(ic) + thetam(l) = parcels%volume(ic) * parcels%theta(ic) #ifndef ENABLE_DRY_MODE - hum(l) = parcels%volume(ic) * parcels%humidity(ic) + qvm(l) = parcels%volume(ic) * parcels%qv(ic) + qlm(l) = parcels%volume(ic) * parcels%ql(ic) #endif vortm(:, l) = parcels%volume(ic) * parcels%vorticity(:, ic) @@ -151,9 +153,10 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, Bm, vm) posm(3, n) = posm(3, n) + parcels%volume(is) * parcels%position(3, is) ! Accumulate buoyancy and humidity - buoym(n) = buoym(n) + parcels%volume(is) * parcels%buoyancy(is) + thetam(n) = thetam(n) + parcels%volume(is) * parcels%theta(is) #ifndef ENABLE_DRY_MODE - hum(n) = hum(n) + parcels%volume(is) * parcels%humidity(is) + qvm(n) = qvm(n) + parcels%volume(is) * parcels%qv(is) + qlm(n) = qlm(n) + parcels%volume(is) * parcels%ql(is) #endif vortm(:, n) = vortm(:, n) + parcels%volume(is) * parcels%vorticity(:, is) enddo @@ -180,9 +183,10 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, Bm, vm) call apply_periodic_bc(posm(:, m)) ! buoyancy and humidity - buoym(m) = vmerge * buoym(m) + thetam(m) = vmerge * thetam(m) #ifndef ENABLE_DRY_MODE - hum(m) = vmerge * hum(m) + qvm(m) = vmerge * qvm(m) + qlm(m) = vmerge * qlm(m) #endif vortm(:, m) = vmerge * vortm(:, m) enddo @@ -220,9 +224,10 @@ subroutine do_group_merge(parcels, isma, iclo, n_merge, Bm, vm) parcels%position(2, ic) = posm(2, l) parcels%position(3, ic) = posm(3, l) - parcels%buoyancy(ic) = buoym(l) + parcels%theta(ic) = thetam(l) #ifndef ENABLE_DRY_MODE - parcels%humidity(ic) = hum(l) + parcels%qv(ic) = qvm(l) + parcels%ql(ic) = qlm(l) #endif parcels%vorticity(:, ic) = vortm(:, l) diff --git a/src/3d/fields/field_diagnostics.f90 b/src/3d/fields/field_diagnostics.f90 index b932063dc..352a9efc1 100644 --- a/src/3d/fields/field_diagnostics.f90 +++ b/src/3d/fields/field_diagnostics.f90 @@ -11,6 +11,7 @@ module field_diagnostics use mpi_collectives, only : mpi_blocking_reduce use physics, only : ape_calculation use ape_density, only : ape_den + use parcel_interpl, only: par2grid_diag implicit none integer :: field_stats_timer @@ -47,6 +48,7 @@ subroutine calculate_field_diagnostics ! ! calculate locally ! + call par2grid_diag ! do not take halo cells into account field_stats(IDX_RMS_V) = sum((volg(lo(3):hi(3), lo(2):hi(2), lo(1):hi(1)) - vcell) ** 2) diff --git a/src/3d/fields/field_diagnostics_netcdf.f90 b/src/3d/fields/field_diagnostics_netcdf.f90 index 6ef467526..8318f1767 100644 --- a/src/3d/fields/field_diagnostics_netcdf.f90 +++ b/src/3d/fields/field_diagnostics_netcdf.f90 @@ -14,6 +14,7 @@ module field_diagnostics_netcdf use mpi_timer, only : start_timer, stop_timer use options, only : write_netcdf_options use physics, only : write_physical_quantities, ape_calculation + implicit none private diff --git a/src/3d/fields/field_netcdf.f90 b/src/3d/fields/field_netcdf.f90 index 290c09610..97e195de2 100644 --- a/src/3d/fields/field_netcdf.f90 +++ b/src/3d/fields/field_netcdf.f90 @@ -11,6 +11,8 @@ module field_netcdf use mpi_layout, only : box use parameters, only : write_zeta_boundary_flag use mpi_utils, only : mpi_stop + use parcel_interpl, only: par2grid_diag + implicit none integer :: field_io_timer @@ -338,6 +340,8 @@ subroutine write_netcdf_fields(t) return endif + call par2grid_diag + call open_netcdf_file(ncfname, NF90_WRITE, ncid) lo = box%lo diff --git a/src/3d/parcels/parcel_container.f90 b/src/3d/parcels/parcel_container.f90 index 3ed3e2894..5b30732f4 100644 --- a/src/3d/parcels/parcel_container.f90 +++ b/src/3d/parcels/parcel_container.f90 @@ -317,7 +317,7 @@ subroutine parcel_serialize(n, buffer) buffer(IDX_X_VOR:IDX_Z_VOR) = parcels%vorticity(:, n) buffer(IDX_B11:IDX_B23) = parcels%B(:, n) buffer(IDX_VOL) = parcels%volume(n) - buffer(IDX_THETA) = parcels%theta(n) + buffer(IDX_THETA) = parcels%theta(n) #ifndef ENABLE_DRY_MODE buffer(IDX_QV) = parcels%qv(n) buffer(IDX_QL) = parcels%ql(n) diff --git a/src/3d/parcels/parcel_netcdf.f90 b/src/3d/parcels/parcel_netcdf.f90 index e7e830766..9efa41417 100644 --- a/src/3d/parcels/parcel_netcdf.f90 +++ b/src/3d/parcels/parcel_netcdf.f90 @@ -249,7 +249,7 @@ subroutine create_netcdf_parcel_file(basename, overwrite, l_restart) #ifndef ENABLE_DRY_MODE call define_netcdf_dataset(ncid=ncid, & name='qv', & - long_name='water vapour spec. hum.', & + long_name='parcel water vapour spec. hum.',& std_name='', & unit='kg/kg', & dtype=NF90_DOUBLE, & @@ -258,7 +258,7 @@ subroutine create_netcdf_parcel_file(basename, overwrite, l_restart) call define_netcdf_dataset(ncid=ncid, & name='ql', & - long_name='liquid water spec. hum.', & + long_name='parcel liquid water spec. hum.',& std_name='', & unit='kg/kg', & dtype=NF90_DOUBLE, & diff --git a/src/3d/parcels/parcel_split.f90 b/src/3d/parcels/parcel_split.f90 index 365a2d334..ee727230b 100644 --- a/src/3d/parcels/parcel_split.f90 +++ b/src/3d/parcels/parcel_split.f90 @@ -88,9 +88,10 @@ subroutine parcel_split(parcels, threshold) parcels%vorticity(:, n_thread_loc) = parcels%vorticity(:, n) parcels%volume(n_thread_loc) = parcels%volume(n) - parcels%buoyancy(n_thread_loc) = parcels%buoyancy(n) + parcels%theta(n_thread_loc) = parcels%theta(n) #ifndef ENABLE_DRY_MODE - parcels%humidity(n_thread_loc) = parcels%humidity(n) + parcels%qv(n_thread_loc) = parcels%qv(n) + parcels%ql(n_thread_loc) = parcels%ql(n) #endif V(:, 1) = V(:, 1) * dh * dsqrt(D(1)) parcels%position(:, n_thread_loc) = parcels%position(:, n) - V(:, 1) @@ -104,11 +105,11 @@ subroutine parcel_split(parcels, threshold) call apply_reflective_bc(parcels%position(:, n), parcels%B(:, n)) ! save parcel indices of child parcels for the - ! halo swap routine +!~ ! halo swap routine pid(n) = n pid(n_thread_loc) = n_thread_loc enddo - !$omp end do + !$omp end dosrc/2d/parcels/parcel_diagnostics.f90: !$omp end parallel n_parcel_splits = n_parcel_splits + n_parcels - last_index diff --git a/unit-tests/3d/test_ellipsoid_split.f90 b/unit-tests/3d/test_ellipsoid_split.f90 index 4aaeedfdb..401015241 100644 --- a/unit-tests/3d/test_ellipsoid_split.f90 +++ b/unit-tests/3d/test_ellipsoid_split.f90 @@ -72,9 +72,10 @@ subroutine setup_parcels parcels%position(:, 1) = zero parcels%volume(1) = four / three * abc * pi - parcels%buoyancy(1) = one + parcels%theta(1) = one #ifndef ENABLE_DRY_MODE - parcels%humidity(1) = one + parcels%qv(1) = one + parcels%ql(1) = one #endif ! 7 Nov 2021 ! https://mathworld.wolfram.com/SphericalCoordinates.html @@ -129,9 +130,10 @@ subroutine check_result error = max(error, abs(get_B33(parcels%B(:, 1), parcels%volume(1)) - B33)) error = max(error, sum(abs(pos(:, 1) - parcels%position(:, 1)))) error = max(error, abs(f12 * four / three * abc * pi - parcels%volume(1))) - error = max(error, abs(parcels%buoyancy(1) - one)) + error = max(error, abs(parcels%theta(1) - one)) #ifndef ENABLE_DRY_MODE - error = max(error, abs(parcels%humidity(1) - one)) + error = max(error, abs(parcels%qv(1) - one)) + error = max(error, abs(parcels%ql(1) - one)) #endif ! second parcel @@ -144,9 +146,10 @@ subroutine check_result error = max(error, sum(abs(pos(:, 2) - parcels%position(:, 2)))) error = max(error, abs(f12 * four / three * abc * pi - parcels%volume(2))) error = max(error, dble(abs(n_parcels - 2))) - error = max(error, abs(parcels%buoyancy(2) - one)) + error = max(error, abs(parcels%theta(2) - one)) #ifndef ENABLE_DRY_MODE - error = max(error, abs(parcels%humidity(2) - one)) + error = max(error, abs(parcels%qv(2) - one)) + error = max(error, abs(parcels%ql(2) - one)) #endif end subroutine check_result diff --git a/unit-tests/3d/test_mpi_nearest_1.f90 b/unit-tests/3d/test_mpi_nearest_1.f90 index c48e0c927..9ea00b38a 100644 --- a/unit-tests/3d/test_mpi_nearest_1.f90 +++ b/unit-tests/3d/test_mpi_nearest_1.f90 @@ -108,7 +108,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 @@ -117,7 +117,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 enddo @@ -126,7 +126,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y + dble(m) * dx(2) * 0.45 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 enddo diff --git a/unit-tests/3d/test_mpi_nearest_10.f90 b/unit-tests/3d/test_mpi_nearest_10.f90 index 22e33b8ef..35e086f6a 100644 --- a/unit-tests/3d/test_mpi_nearest_10.f90 +++ b/unit-tests/3d/test_mpi_nearest_10.f90 @@ -109,7 +109,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -117,7 +117,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -125,7 +125,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.33d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement diff --git a/unit-tests/3d/test_mpi_nearest_11.f90 b/unit-tests/3d/test_mpi_nearest_11.f90 index 94850e359..f2501fe8f 100644 --- a/unit-tests/3d/test_mpi_nearest_11.f90 +++ b/unit-tests/3d/test_mpi_nearest_11.f90 @@ -109,7 +109,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -117,7 +117,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -125,7 +125,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.33d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement diff --git a/unit-tests/3d/test_mpi_nearest_12.f90 b/unit-tests/3d/test_mpi_nearest_12.f90 index a31ef9b36..7e034fd70 100644 --- a/unit-tests/3d/test_mpi_nearest_12.f90 +++ b/unit-tests/3d/test_mpi_nearest_12.f90 @@ -152,7 +152,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -160,7 +160,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -168,7 +168,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.28d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_1 @@ -193,7 +193,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -201,7 +201,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -209,7 +209,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.34d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_2 @@ -234,7 +234,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -242,7 +242,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -250,7 +250,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.28d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_3 @@ -275,7 +275,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -283,7 +283,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -291,7 +291,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.34d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_4 @@ -320,7 +320,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -328,7 +328,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -336,7 +336,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! @@ -348,7 +348,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -356,7 +356,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -364,7 +364,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.30d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_5 @@ -393,7 +393,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -401,7 +401,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -409,7 +409,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! @@ -421,7 +421,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -429,7 +429,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.32d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -437,7 +437,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.33d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_6 diff --git a/unit-tests/3d/test_mpi_nearest_13.f90 b/unit-tests/3d/test_mpi_nearest_13.f90 index 7e2332358..8015d225c 100644 --- a/unit-tests/3d/test_mpi_nearest_13.f90 +++ b/unit-tests/3d/test_mpi_nearest_13.f90 @@ -178,7 +178,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -186,7 +186,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -194,7 +194,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.34d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_1 @@ -219,7 +219,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -227,7 +227,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -235,7 +235,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.26d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_2 @@ -260,7 +260,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -268,7 +268,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -276,7 +276,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.34d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_3 @@ -301,7 +301,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -309,7 +309,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.42d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -317,7 +317,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.28d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_4 @@ -346,7 +346,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -354,7 +354,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -362,7 +362,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! @@ -374,7 +374,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -382,7 +382,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -390,7 +390,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.34d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_5 @@ -419,7 +419,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -427,7 +427,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -435,7 +435,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! @@ -447,7 +447,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -455,7 +455,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -463,7 +463,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.34d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_6 @@ -492,7 +492,7 @@ subroutine cell_placement_7(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -500,7 +500,7 @@ subroutine cell_placement_7(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -508,7 +508,7 @@ subroutine cell_placement_7(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! @@ -520,7 +520,7 @@ subroutine cell_placement_7(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -528,7 +528,7 @@ subroutine cell_placement_7(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.28d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -536,7 +536,7 @@ subroutine cell_placement_7(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_7 @@ -565,7 +565,7 @@ subroutine cell_placement_8(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -573,7 +573,7 @@ subroutine cell_placement_8(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -581,7 +581,7 @@ subroutine cell_placement_8(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! @@ -593,7 +593,7 @@ subroutine cell_placement_8(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -601,7 +601,7 @@ subroutine cell_placement_8(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.28d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -609,7 +609,7 @@ subroutine cell_placement_8(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_8 diff --git a/unit-tests/3d/test_mpi_nearest_14.f90 b/unit-tests/3d/test_mpi_nearest_14.f90 index 6bcf177f7..1752f11dd 100644 --- a/unit-tests/3d/test_mpi_nearest_14.f90 +++ b/unit-tests/3d/test_mpi_nearest_14.f90 @@ -152,7 +152,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -160,7 +160,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -168,7 +168,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.34d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_1 @@ -193,7 +193,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -201,7 +201,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -209,7 +209,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.34d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_2 @@ -234,7 +234,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -242,7 +242,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -250,7 +250,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.34d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_3 @@ -275,7 +275,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -283,7 +283,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -291,7 +291,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.34d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_4 @@ -320,7 +320,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -328,7 +328,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -336,7 +336,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! @@ -348,7 +348,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.28d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -356,7 +356,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -364,7 +364,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_5 @@ -393,7 +393,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -401,7 +401,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -409,7 +409,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! @@ -421,7 +421,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.28d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -429,7 +429,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -437,7 +437,7 @@ subroutine cell_placement_6(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_6 diff --git a/unit-tests/3d/test_mpi_nearest_15.f90 b/unit-tests/3d/test_mpi_nearest_15.f90 index d0add1759..2afec529c 100644 --- a/unit-tests/3d/test_mpi_nearest_15.f90 +++ b/unit-tests/3d/test_mpi_nearest_15.f90 @@ -139,7 +139,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -147,7 +147,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -155,7 +155,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -163,7 +163,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_1 @@ -188,7 +188,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -196,7 +196,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -204,7 +204,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -212,7 +212,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_2 @@ -237,7 +237,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -245,7 +245,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -253,7 +253,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -261,7 +261,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_3 @@ -290,7 +290,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -298,7 +298,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -306,7 +306,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -314,7 +314,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 @@ -327,7 +327,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -335,7 +335,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.45d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -343,7 +343,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.35d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -351,7 +351,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.2d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_4 @@ -376,7 +376,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -384,7 +384,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -392,7 +392,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -400,7 +400,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 @@ -413,7 +413,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.2d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -421,7 +421,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -429,7 +429,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.45d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -437,7 +437,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.2d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_5 diff --git a/unit-tests/3d/test_mpi_nearest_16.f90 b/unit-tests/3d/test_mpi_nearest_16.f90 index 7b4d5b94a..a2d803e15 100644 --- a/unit-tests/3d/test_mpi_nearest_16.f90 +++ b/unit-tests/3d/test_mpi_nearest_16.f90 @@ -139,7 +139,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -147,7 +147,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -155,7 +155,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -163,7 +163,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.43d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_1 @@ -188,7 +188,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -196,7 +196,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -204,7 +204,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -212,7 +212,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.43d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_2 @@ -236,7 +236,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -244,7 +244,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -252,7 +252,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -260,7 +260,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.43d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_3 @@ -289,7 +289,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -297,7 +297,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -305,7 +305,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -313,7 +313,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 @@ -326,7 +326,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.3d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -334,7 +334,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.45d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -342,7 +342,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.3d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -350,7 +350,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.2d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_4 @@ -379,7 +379,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -387,7 +387,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -395,7 +395,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -403,7 +403,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 @@ -416,7 +416,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.2d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -424,7 +424,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -432,7 +432,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.45d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -440,7 +440,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.35d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_5 diff --git a/unit-tests/3d/test_mpi_nearest_17.f90 b/unit-tests/3d/test_mpi_nearest_17.f90 index dd6ccb501..f9e3c8df5 100644 --- a/unit-tests/3d/test_mpi_nearest_17.f90 +++ b/unit-tests/3d/test_mpi_nearest_17.f90 @@ -139,7 +139,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -147,7 +147,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -155,7 +155,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -163,7 +163,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.43d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_1 @@ -188,7 +188,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -196,7 +196,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -204,7 +204,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -212,7 +212,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.43d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_2 @@ -236,7 +236,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.35d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -244,7 +244,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -252,7 +252,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -260,7 +260,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.35d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_3 @@ -289,7 +289,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -297,7 +297,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -305,7 +305,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -313,7 +313,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 @@ -326,7 +326,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -334,7 +334,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -342,7 +342,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.3d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -350,7 +350,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.15d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 @@ -380,7 +380,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -388,7 +388,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -396,7 +396,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -404,7 +404,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 @@ -417,7 +417,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.25d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel B @@ -425,7 +425,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -433,7 +433,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.45d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -441,7 +441,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.25d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_5 diff --git a/unit-tests/3d/test_mpi_nearest_18.f90 b/unit-tests/3d/test_mpi_nearest_18.f90 index 678468c3f..749be5eff 100644 --- a/unit-tests/3d/test_mpi_nearest_18.f90 +++ b/unit-tests/3d/test_mpi_nearest_18.f90 @@ -126,7 +126,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.38d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -134,7 +134,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.38d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -142,7 +142,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.38d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -150,7 +150,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.38d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_1 @@ -175,7 +175,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.3d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -183,7 +183,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.38d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -191,7 +191,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.3d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -199,7 +199,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.38d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_2 @@ -228,7 +228,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -236,7 +236,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -244,7 +244,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -252,7 +252,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 @@ -265,7 +265,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.45d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -273,7 +273,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.45d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -281,7 +281,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.3d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -289,7 +289,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.2d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_3 @@ -318,7 +318,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -326,7 +326,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -334,7 +334,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -342,7 +342,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 @@ -355,7 +355,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.45d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -363,7 +363,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.45d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -371,7 +371,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.3d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -379,7 +379,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.2d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_4 diff --git a/unit-tests/3d/test_mpi_nearest_19.f90 b/unit-tests/3d/test_mpi_nearest_19.f90 index 60c74f733..6b1bc009e 100644 --- a/unit-tests/3d/test_mpi_nearest_19.f90 +++ b/unit-tests/3d/test_mpi_nearest_19.f90 @@ -139,7 +139,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.23d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -147,7 +147,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -155,7 +155,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.43d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel D @@ -163,7 +163,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.43d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_1 @@ -188,7 +188,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.23d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -196,7 +196,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -204,7 +204,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.43d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel D @@ -212,7 +212,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.33d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_2 @@ -237,7 +237,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.43d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -245,7 +245,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -253,7 +253,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.43d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel D @@ -261,7 +261,7 @@ subroutine cell_placement_3(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_3 @@ -290,7 +290,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -298,7 +298,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -306,7 +306,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel D @@ -314,7 +314,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 @@ -327,7 +327,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -335,7 +335,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -343,7 +343,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.25d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel D @@ -351,7 +351,7 @@ subroutine cell_placement_4(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.13d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_4 @@ -380,7 +380,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -388,7 +388,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -396,7 +396,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel D @@ -404,7 +404,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 @@ -417,7 +417,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.2d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -425,7 +425,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.35d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -433,7 +433,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.46d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel D @@ -441,7 +441,7 @@ subroutine cell_placement_5(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.46d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_5 diff --git a/unit-tests/3d/test_mpi_nearest_2.f90 b/unit-tests/3d/test_mpi_nearest_2.f90 index 1f33b85b9..72815e153 100644 --- a/unit-tests/3d/test_mpi_nearest_2.f90 +++ b/unit-tests/3d/test_mpi_nearest_2.f90 @@ -108,7 +108,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 @@ -117,7 +117,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 enddo @@ -126,7 +126,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 enddo diff --git a/unit-tests/3d/test_mpi_nearest_20.f90 b/unit-tests/3d/test_mpi_nearest_20.f90 index eae7d50d7..91e8c10d2 100644 --- a/unit-tests/3d/test_mpi_nearest_20.f90 +++ b/unit-tests/3d/test_mpi_nearest_20.f90 @@ -106,8 +106,8 @@ program test_mpi_nearest_20 ! do n = 1, n_merge ! is = isma(n) ! ic = iclo(n) -! print *, comm%rank, parcels%position(1, is), parcels%position(2, is), int(parcels%buoyancy(is)), & -! parcels%position(1, ic), parcels%position(2, ic), int(parcels%buoyancy(ic)) +! print *, comm%rank, parcels%position(1, is), parcels%position(2, is), int(parcels%theta(is)), & +! parcels%position(1, ic), parcels%position(2, ic), int(parcels%theta(ic)) ! enddo ! endif ! call MPI_Barrier(comm%world, comm%err) @@ -159,7 +159,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.35d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -167,7 +167,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel c @@ -175,7 +175,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -183,7 +183,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel e @@ -191,7 +191,7 @@ subroutine cell_placement_1(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.38d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_1 @@ -216,7 +216,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.35d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -224,7 +224,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel C @@ -232,7 +232,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel d @@ -240,7 +240,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel e @@ -248,7 +248,7 @@ subroutine cell_placement_2(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.38d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement_2 diff --git a/unit-tests/3d/test_mpi_nearest_3.f90 b/unit-tests/3d/test_mpi_nearest_3.f90 index 558372fa5..d653b8d8c 100644 --- a/unit-tests/3d/test_mpi_nearest_3.f90 +++ b/unit-tests/3d/test_mpi_nearest_3.f90 @@ -110,7 +110,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel @@ -118,7 +118,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel @@ -126,7 +126,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement diff --git a/unit-tests/3d/test_mpi_nearest_4.f90 b/unit-tests/3d/test_mpi_nearest_4.f90 index 24e0a762e..12411a694 100644 --- a/unit-tests/3d/test_mpi_nearest_4.f90 +++ b/unit-tests/3d/test_mpi_nearest_4.f90 @@ -110,7 +110,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y + 0.1d0 * dx(2) parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel @@ -118,7 +118,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.35d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! big parcel @@ -126,7 +126,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y - 0.42d0 * dx(2) parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement diff --git a/unit-tests/3d/test_mpi_nearest_5.f90 b/unit-tests/3d/test_mpi_nearest_5.f90 index 19756d931..91c7147d6 100644 --- a/unit-tests/3d/test_mpi_nearest_5.f90 +++ b/unit-tests/3d/test_mpi_nearest_5.f90 @@ -109,7 +109,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -117,7 +117,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -125,7 +125,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement diff --git a/unit-tests/3d/test_mpi_nearest_6.f90 b/unit-tests/3d/test_mpi_nearest_6.f90 index 16412b8cb..137374665 100644 --- a/unit-tests/3d/test_mpi_nearest_6.f90 +++ b/unit-tests/3d/test_mpi_nearest_6.f90 @@ -109,7 +109,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -117,7 +117,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -125,7 +125,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement diff --git a/unit-tests/3d/test_mpi_nearest_7.f90 b/unit-tests/3d/test_mpi_nearest_7.f90 index a0cfc0ff3..8107e1241 100644 --- a/unit-tests/3d/test_mpi_nearest_7.f90 +++ b/unit-tests/3d/test_mpi_nearest_7.f90 @@ -109,7 +109,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -117,7 +117,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y - dx(1) * 0.4d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -125,7 +125,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y - dx(1) * 0.3d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement diff --git a/unit-tests/3d/test_mpi_nearest_8.f90 b/unit-tests/3d/test_mpi_nearest_8.f90 index d679b8ead..f6c4dc358 100644 --- a/unit-tests/3d/test_mpi_nearest_8.f90 +++ b/unit-tests/3d/test_mpi_nearest_8.f90 @@ -109,7 +109,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.24d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -117,7 +117,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.42d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -125,7 +125,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.46d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement diff --git a/unit-tests/3d/test_mpi_nearest_9.f90 b/unit-tests/3d/test_mpi_nearest_9.f90 index 7e35b7719..d97c36335 100644 --- a/unit-tests/3d/test_mpi_nearest_9.f90 +++ b/unit-tests/3d/test_mpi_nearest_9.f90 @@ -109,7 +109,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel b @@ -117,7 +117,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y - dx(2) * 0.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 ! small parcel a @@ -125,7 +125,7 @@ subroutine cell_placement(l, i, j, k) parcels%position(2, l) = y + dx(2) * 0.28d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + comm%rank * 100 + parcels%theta(l) = l + comm%rank * 100 l = l + 1 end subroutine cell_placement diff --git a/unit-tests/3d/test_mpi_parcel_communicate.f90 b/unit-tests/3d/test_mpi_parcel_communicate.f90 index abceefa22..c724110a1 100644 --- a/unit-tests/3d/test_mpi_parcel_communicate.f90 +++ b/unit-tests/3d/test_mpi_parcel_communicate.f90 @@ -96,7 +96,7 @@ program test_mpi_parcel_communicate parcels%volume(1:n_parcels) = comm%rank + 1 parcels%B(:, 1:n_parcels) = comm%rank + 1 parcels%vorticity(:, 1:n_parcels) = comm%rank + 1 - parcels%buoyancy(1:n_parcels) = comm%rank + 1 + parcels%theta(1:n_parcels) = comm%rank + 1 call parcel_communicate diff --git a/unit-tests/3d/test_mpi_parcel_delete.f90 b/unit-tests/3d/test_mpi_parcel_delete.f90 index 82dfb3c5f..a58d9f3e7 100644 --- a/unit-tests/3d/test_mpi_parcel_delete.f90 +++ b/unit-tests/3d/test_mpi_parcel_delete.f90 @@ -50,9 +50,10 @@ program test_mpi_parcel_delete parcels%B(4, n) = 10.0d0 + a parcels%B(5, n) = 11.0d0 + a parcels%volume(n) = 12.0d0 + a - parcels%buoyancy(n) = 13.0d0 + a + parcels%theta(n) = 13.0d0 + a #ifndef ENABLE_DRY_MODE - parcels%humidity(n) = 14.0d0 + a + parcels%qv(n) = 14.0d0 + a + parcels%ql(n) = 15.0d0 + a #endif enddo @@ -107,9 +108,10 @@ program test_mpi_parcel_delete passed = (passed .and. (parcels%B(4, ii(n)) - (10.0d0 + a) == zero)) passed = (passed .and. (parcels%B(5, ii(n)) - (11.0d0 + a) == zero)) passed = (passed .and. (parcels%volume(ii(n)) - (12.0d0 + a) == zero)) - passed = (passed .and. (parcels%buoyancy(ii(n)) - (13.0d0 + a) == zero)) + passed = (passed .and. (parcels%theta(ii(n)) - (13.0d0 + a) == zero)) #ifndef ENABLE_DRY_MODE - passed = (passed .and. (parcels%humidity(ii(n)) - (14.0d0 + a) == zero)) + passed = (passed .and. (parcels%qv(ii(n)) - (14.0d0 + a) == zero)) + passed = (passed .and. (parcels%ql(ii(n)) - (15.0d0 + a) == zero)) #endif i = i + 1 endif diff --git a/unit-tests/3d/test_mpi_parcel_diagnostics.f90 b/unit-tests/3d/test_mpi_parcel_diagnostics.f90 index b4c5525d9..adf4e241b 100644 --- a/unit-tests/3d/test_mpi_parcel_diagnostics.f90 +++ b/unit-tests/3d/test_mpi_parcel_diagnostics.f90 @@ -53,7 +53,7 @@ program test_mpi_parcel_diagnostics parcels%position(1, l) = corner(1) + dx(1) * (dble(i) - f12) * im parcels%position(2, l) = corner(2) + dx(2) * (dble(j) - f12) * im parcels%position(3, l) = corner(3) + dx(3) * (dble(k) - f12) * im - parcels%buoyancy(l) = parcels%position(3, l) + parcels%theta(l) = parcels%position(3, l) l = l + 1 enddo enddo diff --git a/unit-tests/3d/test_mpi_parcel_pack.f90 b/unit-tests/3d/test_mpi_parcel_pack.f90 index ad8d3c3b9..58c94adc8 100644 --- a/unit-tests/3d/test_mpi_parcel_pack.f90 +++ b/unit-tests/3d/test_mpi_parcel_pack.f90 @@ -48,7 +48,7 @@ program test_mpi_parcel_pack parcels%B(4, n) = 10.0d0 + a parcels%B(5, n) = 11.0d0 + a parcels%volume(n) = 12.0d0 + a - parcels%buoyancy(n) = 13.0d0 + a + parcels%theta(n) = 13.0d0 + a #ifndef ENABLE_DRY_MODE parcels%humidity(n) = 14.0d0 + a #endif @@ -93,9 +93,10 @@ program test_mpi_parcel_pack passed = (passed .and. (parcels%B(4, l) - buffer(i + IDX_B22) == zero)) passed = (passed .and. (parcels%B(5, l) - buffer(i + IDX_B23) == zero)) passed = (passed .and. (parcels%volume(l) - buffer(i + IDX_VOL) == zero)) - passed = (passed .and. (parcels%buoyancy(l) - buffer(i + IDX_BUO) == zero)) + passed = (passed .and. (parcels%theta(l) - buffer(i + IDX_THETA) == zero)) #ifndef ENABLE_DRY_MODE - passed = (passed .and. (parcels%humidity(l) - buffer(i + IDX_HUM) == zero)) + passed = (passed .and. (parcels%qv(l) - buffer(i + IDX_QV) == zero)) + passed = (passed .and. (parcels%ql(l) - buffer(i + IDX_QL) == zero)) #endif enddo diff --git a/unit-tests/3d/test_mpi_parcel_read.f90 b/unit-tests/3d/test_mpi_parcel_read.f90 index 904e2314b..8f8c93dee 100644 --- a/unit-tests/3d/test_mpi_parcel_read.f90 +++ b/unit-tests/3d/test_mpi_parcel_read.f90 @@ -46,10 +46,11 @@ program test_mpi_parcel_read parcels%B(:, n) = start_index + n parcels%volume(n) = start_index + n parcels%vorticity(:, n) = start_index + n - parcels%buoyancy(n) = start_index + n + parcels%theta(n) = start_index + n #ifndef ENABLE_DRY_MODE - parcels%humidity(n) = start_index + n + parcels%qv(n) = start_index + n + parcels%ql(n) = start_index + n #endif enddo @@ -69,9 +70,10 @@ program test_mpi_parcel_read parcels%B = 0 parcels%volume = 0 parcels%vorticity = 0 - parcels%buoyancy = 0 + parcels%theta = 0 #ifndef ENABLE_DRY_MODE - parcels%humidity = 0 + parcels%qv = 0 + parcels%ql = 0 #endif call read_netcdf_parcels('nctest_0000000001_parcels.nc') @@ -93,9 +95,10 @@ program test_mpi_parcel_read passed = (passed .and. (abs(parcels%vorticity(1, n) - res) == zero)) passed = (passed .and. (abs(parcels%vorticity(2, n) - res) == zero)) passed = (passed .and. (abs(parcels%vorticity(3, n) - res) == zero)) - passed = (passed .and. (abs(parcels%buoyancy(n) - res) == zero)) + passed = (passed .and. (abs(parcels%theta(n) - res) == zero)) #ifndef ENABLE_DRY_MODE - passed = (passed .and. (abs(parcels%humidity(n) - res) == zero)) + passed = (passed .and. (abs(parcels%qv(n) - res) == zero)) + passed = (passed .and. (abs(parcels%ql(n) - res) == zero)) #endif enddo endif diff --git a/unit-tests/3d/test_mpi_parcel_read_rejection.f90 b/unit-tests/3d/test_mpi_parcel_read_rejection.f90 index cc5449302..74965a626 100644 --- a/unit-tests/3d/test_mpi_parcel_read_rejection.f90 +++ b/unit-tests/3d/test_mpi_parcel_read_rejection.f90 @@ -25,14 +25,14 @@ program test_mpi_parcel_read_rejection double precision :: x_sum, y_sum, z_sum integer :: ncid - integer :: npar_dim_id, vol_id, buo_id, & + integer :: npar_dim_id, vol_id, theta_id, & x_pos_id, y_pos_id, z_pos_id, & x_vor_id, y_vor_id, z_vor_id, & b11_id, b12_id, b13_id, & b22_id, b23_id, & t_axis_id, t_dim_id, mpi_dim_id #ifndef ENABLE_DRY_MODE - integer :: hum_id + integer :: qv_id, ql_id #endif character(len=512) :: ncbasename @@ -97,10 +97,11 @@ program test_mpi_parcel_read_rejection parcels%B(:, 1:n_parcels) = comm%rank + 1 parcels%volume(1:n_parcels) = comm%rank + 1 parcels%vorticity(:, 1:n_parcels) = comm%rank + 1 - parcels%buoyancy(1:n_parcels) = comm%rank + 1 + parcels%theta(1:n_parcels) = comm%rank + 1 #ifndef ENABLE_DRY_MODE - parcels%humidity(1:n_parcels) = comm%rank + 1 + parcels%qv(1:n_parcels) = comm%rank + 1 + parcels%ql(1:n_parcels) = comm%rank + 1 #endif call create_file('nctest') @@ -119,9 +120,10 @@ program test_mpi_parcel_read_rejection parcels%B = 0 parcels%volume = 0 parcels%vorticity = 0 - parcels%buoyancy = 0 + parcels%theta = 0 #ifndef ENABLE_DRY_MODE - parcels%humidity = 0 + parcels%qv = 0 + parcels%ql = 0 #endif call read_netcdf_parcels('nctest_0000000001_parcels.nc') @@ -143,9 +145,10 @@ program test_mpi_parcel_read_rejection passed = (passed .and. (maxval(abs(parcels%vorticity(1, 1:n_parcels) - res)) == zero)) passed = (passed .and. (maxval(abs(parcels%vorticity(2, 1:n_parcels) - res)) == zero)) passed = (passed .and. (maxval(abs(parcels%vorticity(3, 1:n_parcels) - res)) == zero)) - passed = (passed .and. (maxval(abs(parcels%buoyancy(1:n_parcels) - res)) == zero)) + passed = (passed .and. (maxval(abs(parcels%theta(1:n_parcels) - res)) == zero)) #ifndef ENABLE_DRY_MODE - passed = (passed .and. (maxval(abs(parcels%humidity(1:n_parcels) - res)) == zero)) + passed = (passed .and. (maxval(abs(parcels%qv(1:n_parcels) - res)) == zero)) + passed = (passed .and. (maxval(abs(parcels%ql(1:n_parcels) - res)) == zero)) #endif endif @@ -315,23 +318,32 @@ subroutine create_file(basename) varid=z_vor_id) call define_netcdf_dataset(ncid=ncid, & - name='buoyancy', & - long_name='parcel buoyancy', & + name='theta', & + long_name='parcel potential temperature',& std_name='', & - unit='m/s^2', & + unit='K', & dtype=NF90_DOUBLE, & dimids=dimids, & - varid=buo_id) + varid=theta_id) #ifndef ENABLE_DRY_MODE call define_netcdf_dataset(ncid=ncid, & - name='humidity', & - long_name='parcel humidity', & + name='qv', & + long_name='parcel water vapour spec. hum.',& std_name='', & - unit='1', & + unit='kg/kg', & dtype=NF90_DOUBLE, & dimids=dimids, & - varid=hum_id) + varid=qv_id) + + call define_netcdf_dataset(ncid=ncid, & + name='ql', & + long_name='parcel liquid water spec. hum.',& + std_name='', & + unit='kg/kg', & + dtype=NF90_DOUBLE, & + dimids=dimids, & + varid=ql_id) #endif call close_definition(ncid) @@ -384,10 +396,11 @@ subroutine write_parcels(t) call write_netcdf_dataset(ncid, y_vor_id, parcels%vorticity(2, 1:n_parcels), start, cnt) call write_netcdf_dataset(ncid, z_vor_id, parcels%vorticity(3, 1:n_parcels), start, cnt) - call write_netcdf_dataset(ncid, buo_id, parcels%buoyancy(1:n_parcels), start, cnt) + call write_netcdf_dataset(ncid, theta_id, parcels%theta(1:n_parcels), start, cnt) #ifndef ENABLE_DRY_MODE - call write_netcdf_dataset(ncid, hum_id, parcels%humidity(1:n_parcels), start, cnt) + call write_netcdf_dataset(ncid, qv_id, parcels%qv(1:n_parcels), start, cnt) + call write_netcdf_dataset(ncid, ql_id, parcels%ql(1:n_parcels), start, cnt) #endif call close_netcdf_file(ncid) diff --git a/unit-tests/3d/test_mpi_parcel_split.f90 b/unit-tests/3d/test_mpi_parcel_split.f90 index b92bc525c..70c4ea0d7 100644 --- a/unit-tests/3d/test_mpi_parcel_split.f90 +++ b/unit-tests/3d/test_mpi_parcel_split.f90 @@ -106,7 +106,7 @@ program test_mpi_parcel_split parcels%volume(1:n_parcels) = f12 * vcell parcels%vorticity(:, 1:n_parcels) = comm%rank + 1 - parcels%buoyancy(1:n_parcels) = comm%rank + 1 + parcels%theta(1:n_parcels) = comm%rank + 1 call parcel_split(parcels, threshold=four) diff --git a/unit-tests/3d/test_mpi_parcel_unpack.f90 b/unit-tests/3d/test_mpi_parcel_unpack.f90 index 4bee8835e..95f8a1bd6 100644 --- a/unit-tests/3d/test_mpi_parcel_unpack.f90 +++ b/unit-tests/3d/test_mpi_parcel_unpack.f90 @@ -40,9 +40,10 @@ program test_mpi_parcel_unpack parcels%B(4, n) = 10.0d0 + a parcels%B(5, n) = 11.0d0 + a parcels%volume(n) = 12.0d0 + a - parcels%buoyancy(n) = 13.0d0 + a + parcels%theta(n) = 13.0d0 + a #ifndef ENABLE_DRY_MODE - parcels%humidity(n) = 14.0d0 + a + parcels%qv(n) = 14.0d0 + a + parcels%ql(n) = 15.0d0 + a #endif enddo @@ -65,9 +66,10 @@ program test_mpi_parcel_unpack buffer(i + IDX_B22) = 10.0d0 + a buffer(i + IDX_B23) = 11.0d0 + a buffer(i + IDX_VOL) = 12.0d0 + a - buffer(i + IDX_BUO) = 13.0d0 + a + buffer(i + IDX_THETA) = 13.0d0 + a #ifndef ENABLE_DRY_MODE - buffer(i + IDX_HUM) = 14.0d0 + a + buffer(i + IDX_QV) = 14.0d0 + a + buffer(i + IDX_QL) = 15.0d0 + a #endif enddo @@ -92,9 +94,10 @@ program test_mpi_parcel_unpack passed = (passed .and. (parcels%B(4, n) == 10.0d0 + a)) passed = (passed .and. (parcels%B(5, n) == 11.0d0 + a)) passed = (passed .and. (parcels%volume(n) == 12.0d0 + a)) - passed = (passed .and. (parcels%buoyancy(n) == 13.0d0 + a)) + passed = (passed .and. (parcels%theta(n) == 13.0d0 + a)) #ifndef ENABLE_DRY_MODE - passed = (passed .and. (parcels%humidity(n) == 14.0d0 + a)) + passed = (passed .and. (parcels%qv(n) == 14.0d0 + a)) + passed = (passed .and. (parcels%ql(n) == 15.0d0 + a)) #endif enddo diff --git a/unit-tests/3d/test_mpi_parcel_write.f90 b/unit-tests/3d/test_mpi_parcel_write.f90 index 02329683b..27236cd87 100644 --- a/unit-tests/3d/test_mpi_parcel_write.f90 +++ b/unit-tests/3d/test_mpi_parcel_write.f90 @@ -29,10 +29,11 @@ program test_mpi_parcel_write parcels%B(:, 1:n_parcels) = comm%rank parcels%volume(1:n_parcels) = comm%rank parcels%vorticity(:, 1:n_parcels) = comm%rank - parcels%buoyancy(1:n_parcels) = comm%rank + parcels%theta(1:n_parcels) = comm%rank #ifndef ENABLE_DRY_MODE - parcels%humidity(1:n_parcels) = comm%rank + parcels%qv(1:n_parcels) = comm%rank + parcels%ql(1:n_parcels) = comm%rank #endif call create_netcdf_parcel_file('nctest', .true., .false.) From 061400fe7cc847548270300404d9eaa0e108f2b3 Mon Sep 17 00:00:00 2001 From: Steven Boeing Date: Thu, 11 May 2023 20:49:57 +0100 Subject: [PATCH 03/17] More ongoing work on thermodynamics --- src/3d/parcels/parcel_interpl.f90 | 271 ++++++++++++++++++++++-------- src/utils/physics.f90 | 8 +- 2 files changed, 210 insertions(+), 69 deletions(-) diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index 58bb2ca71..0272e584d 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -20,7 +20,7 @@ module parcel_interpl , interior_to_halo_communication & , halo_to_interior_communication & , field_halo_swap - use physics, only : glat, lambda_c, q_0 + use physics, only : gravity, theta_0, qv_dens_coeff, r_d, c_p, L_v use omp_lib use mpi_utils, only : mpi_exit_on_error implicit none @@ -43,22 +43,27 @@ module parcel_interpl , IDX_VOR_Y_SWAP = 3 & , IDX_VOR_Z_SWAP = 4 & , IDX_TBUOY_SWAP = 5 -#ifndef ENABLE_DRY_MODE - integer, parameter :: IDX_DBUOY_SWAP = 6 & - , IDX_HUM_SWAP = 7 - integer, parameter :: n_field_swap = 7 -#else integer, parameter :: n_field_swap = 5 + + ! restart indices for par2grid_diag + integer, parameter :: IDX_THETA_SWAP = 1 +#ifndef ENABLE_DRY_MODE + integer, parameter :: IDX_QV_SWAP = 2 & + , IDX_QL_SWAP = 3 + integer, parameter :: n_field_swap_diag = 3 +#else + integer, parameter :: n_field_swap_diag = 1 #endif public :: par2grid & + , par2grid_diag & , vol2grid & , grid2par & , par2grid_timer & , grid2par_timer & , halo_swap_timer & - , trilinear + , trilinear contains @@ -93,7 +98,7 @@ subroutine vol2grid(l_reuse) enddo !$omp end do !$omp end parallel - + call start_timer(halo_swap_timer) call field_halo_swap(volg) call stop_timer(halo_swap_timer) @@ -125,44 +130,30 @@ subroutine par2grid(l_reuse) double precision :: points(3, 4) integer :: n, p, l, i, j, k double precision :: pvol, weight(0:1,0:1,0:1), btot -#ifndef ENABLE_DRY_MODE - double precision :: q_c -#endif call start_timer(par2grid_timer) +#ifndef ENABLE_DRY_MODE + call saturation_adjustment +#endif + vortg = zero volg = zero nparg = zero nsparg = zero -#ifndef ENABLE_DRY_MODE - dbuoyg = zero - humg = zero -#endif tbuoyg = zero !$omp parallel default(shared) -#ifndef ENABLE_DRY_MODE - !$omp do private(n, p, l, i, j, k, points, pvol, weight, btot, q_c) & - !$omp& private( is, js, ks, weights) & - !$omp& reduction(+:nparg, nsparg, vortg, dbuoyg, humg, tbuoyg, volg) -#else !$omp do private(n, p, l, i, j, k, points, pvol, weight, btot) & !$omp& private( is, js, ks, weights) & !$omp& reduction(+:nparg, nsparg, vortg, tbuoyg, volg) -#endif do n = 1, n_parcels pvol = parcels%volume(n) #ifndef ENABLE_DRY_MODE - ! liquid water content - q_c = parcels%humidity(n) & - - q_0 * dexp(lambda_c * (lower(3) - parcels%position(3, n))) - q_c = max(zero, q_c) - ! total buoyancy (including effects of latent heating) - btot = parcels%buoyancy(n) + glat * q_c + btot = gravity*(parcels%theta(n)*(1.0+qv_theta_coefficient*parcels%qv(n)-parcels%ql(n))/theta_0-1.0) #else - btot = parcels%buoyancy(n) + btot = gravity*(parcels%theta(n)/theta_0-1.0) #endif points = get_ellipsoid_points(parcels%position(:, n), & pvol, parcels%B(:, n), n, l_reuse) @@ -186,12 +177,6 @@ subroutine par2grid(l_reuse) vortg(ks:ks+1, js:js+1, is:is+1, l) = vortg(ks:ks+1, js:js+1, is:is+1, l) & + weight * parcels%vorticity(l, n) enddo -#ifndef ENABLE_DRY_MODE - dbuoyg(ks:ks+1, js:js+1, is:is+1) = dbuoyg(ks:ks+1, js:js+1, is:is+1) & - + weight * parcels%buoyancy(n) - humg(ks:ks+1, js:js+1, is:is+1) = humg(ks:ks+1, js:js+1, is:is+1) & - + weight * parcels%humidity(n) -#endif tbuoyg(ks:ks+1, js:js+1, is:is+1) = tbuoyg(ks:ks+1, js:js+1, is:is+1) & + weight * btot volg(ks:ks+1, js:js+1, is:is+1) = volg(ks:ks+1, js:js+1, is:is+1) & @@ -226,19 +211,6 @@ subroutine par2grid(l_reuse) tbuoyg(nz-1, :, :) = tbuoyg(nz-1, :, :) + tbuoyg(nz+1, :, :) !$omp end parallel workshare -#ifndef ENABLE_DRY_MODE - !$omp parallel workshare - dbuoyg(0, :, :) = two * dbuoyg(0, :, :) - dbuoyg(nz, :, :) = two * dbuoyg(nz, :, :) - dbuoyg(1, :, :) = dbuoyg(1, :, :) + dbuoyg(-1, :, :) - dbuoyg(nz-1, :, :) = dbuoyg(nz-1, :, :) + dbuoyg(nz+1, :, :) - humg(0, :, :) = two * humg(0, :, :) - humg(nz, :, :) = two * humg(nz, :, :) - humg(1, :, :) = humg(1, :, :) + humg(-1, :, :) - humg(nz-1, :, :) = humg(nz-1, :, :) + humg(nz+1, :, :) - !$omp end parallel workshare -#endif - ! exclude halo cells to avoid division by zero do p = 1, 3 vortg(0:nz, :, :, p) = vortg(0:nz, :, :, p) / volg(0:nz, :, :) @@ -258,10 +230,6 @@ subroutine par2grid(l_reuse) vortg(-1, :, :, :) = two * vortg(0, :, :, :) - vortg(1, :, :, :) vortg(nz+1, :, :, :) = two * vortg(nz, :, :, :) - vortg(nz-1, :, :, :) -#ifndef ENABLE_DRY_MODE - dbuoyg(0:nz, :, :) = dbuoyg(0:nz, :, :) / volg(0:nz, :, :) - humg(0:nz, :, :) = humg(0:nz, :, :) / volg(0:nz, :, :) -#endif tbuoyg(0:nz, :, :) = tbuoyg(0:nz, :, :) / volg(0:nz, :, :) ! extrapolate to halo grid points (needed to compute @@ -288,6 +256,99 @@ subroutine par2grid(l_reuse) end subroutine par2grid + ! Interpolate parcel quantities to the grid, these consist of the parcel + ! - vorticity + ! - buoyancy + ! - volume + ! It also updates the scalar fields: + ! - nparg, that is the number of parcels per grid cell + ! - nsparg, that is the number of small parcels per grid cell + ! @pre The parcel must be assigned to the correct MPI process. + subroutine par2grid_diag + double precision :: points(3, 4) + integer :: n, p, l, i, j, k + double precision :: pvol, weight(0:1,0:1,0:1) + + thetag = zero +#ifndef ENABLE_DRY_MODE + qcg = zero + qvg = zero +#endif + !$omp parallel default(shared) +#ifndef ENABLE_DRY_MODE + !$omp do private(n, p, l, i, j, k, points, pvol, weight) & + !$omp& private( is, js, ks, weights) & + !$omp& reduction(+:nparg, nsparg, vortg, thetag, qvg, qlg) +#else + !$omp do private(n, p, l, i, j, k, points, pvol, weight) & + !$omp& private( is, js, ks, weights) & + !$omp& reduction(+:nparg, nsparg, vortg, thetag) +#endif + + do n = 1, n_parcels + pvol = parcels%volume(n) + + points = get_ellipsoid_points(parcels%position(:, n), & + pvol, parcels%B(:, n), n, l_reuse) + + call get_index(parcels%position(:, n), i, j, k) + nparg(k, j, i) = nparg(k, j, i) + 1 + if (parcels%volume(n) <= vmin) then + nsparg(k, j, i) = nsparg(k, j, i) + 1 + endif + + ! we have 4 points per ellipsoid + do p = 1, 4 + + call trilinear(points(:, p), is, js, ks, weights) + + ! loop over grid points which are part of the interpolation + ! the weight is a quarter due to 4 points per ellipsoid + weight = f14 * pvol* weights + + thetag(ks:ks+1, js:js+1, is:is+1) = thetag(ks:ks+1, js:js+1, is:is+1) & + + weight * parcels%theta(n) +#ifndef ENABLE_DRY_MODE + qvg(ks:ks+1, js:js+1, is:is+1) = qvg(ks:ks+1, js:js+1, is:is+1) & + + weight * parcels%qv(n) + qlg(ks:ks+1, js:js+1, is:is+1) = qlg(ks:ks+1, js:js+1, is:is+1) & + + weight * parcels%ql(n) +#endif + enddo + enddo + !$omp end do + !$omp end parallel + + call start_timer(halo_swap_timer) + call par2grid_halo_swap_diag + call stop_timer(halo_swap_timer) + + !$omp parallel workshare + thetag(0:nz, :, :) = thetag(0:nz, :, :) / volg(0:nz, :, :) + + ! extrapolate to halo grid points (needed to compute + ! z derivative used for the time step) + thetag(-1, :, :) = two * thetag(0, :, :) - thetag(1, :, :) + thetag(nz+1, :, :) = two * thetag(nz, :, :) - thetag(nz-1, :, :) + +#ifndef ENABLE_DRY_MODE + qvg(0:nz, :, :) = qvg(0:nz, :, :) / volg(0:nz, :, :) + + ! extrapolate to halo grid points (needed to compute + ! z derivative used for the time step) + qvg(-1, :, :) = two * qvg(0, :, :) - qvg(1, :, :) + qvg(nz+1, :, :) = two * qvg(nz, :, :) - qvg(nz-1, :, :) + + qlg(0:nz, :, :) = qlg(0:nz, :, :) / volg(0:nz, :, :) + + ! extrapolate to halo grid points (needed to compute + ! z derivative used for the time step) + qlg(-1, :, :) = two * qlg(0, :, :) - qlg(1, :, :) + qlg(nz+1, :, :) = two * qlg(nz, :, :) - qlg(nz-1, :, :) +#endif + !$omp end parallel workshare + + end subroutine par2grid_diag subroutine par2grid_halo_swap ! we must first fill the interior grid points @@ -305,10 +366,6 @@ subroutine par2grid_halo_swap call field_halo_to_buffer(vortg(:, :, :, I_Y), IDX_VOR_Y_SWAP) call field_halo_to_buffer(vortg(:, :, :, I_Z), IDX_VOR_Z_SWAP) call field_halo_to_buffer(tbuoyg, IDX_TBUOY_SWAP) -#ifndef ENABLE_DRY_MODE - call field_halo_to_buffer(dbuoyg, IDX_DBUOY_SWAP) - call field_halo_to_buffer(humg, IDX_HUM_SWAP) -#endif ! send halo data to valid regions of other processes call halo_to_interior_communication @@ -320,10 +377,7 @@ subroutine par2grid_halo_swap call field_buffer_to_interior(vortg(:, :, :, I_Y), IDX_VOR_Y_SWAP, .true.) call field_buffer_to_interior(vortg(:, :, :, I_Z), IDX_VOR_Z_SWAP, .true.) call field_buffer_to_interior(tbuoyg, IDX_TBUOY_SWAP, .true.) -#ifndef ENABLE_DRY_MODE - call field_buffer_to_interior(dbuoyg, IDX_DBUOY_SWAP, .true.) - call field_buffer_to_interior(humg, IDX_HUM_SWAP, .true.) -#endif + !------------------------------------------------------------------ ! Fill halo: @@ -333,10 +387,6 @@ subroutine par2grid_halo_swap call field_interior_to_buffer(vortg(:, :, :, I_Y), IDX_VOR_Y_SWAP) call field_interior_to_buffer(vortg(:, :, :, I_Z), IDX_VOR_Z_SWAP) call field_interior_to_buffer(tbuoyg, IDX_TBUOY_SWAP) -#ifndef ENABLE_DRY_MODE - call field_interior_to_buffer(dbuoyg, IDX_DBUOY_SWAP) - call field_interior_to_buffer(humg, IDX_HUM_SWAP) -#endif call interior_to_halo_communication @@ -345,15 +395,55 @@ subroutine par2grid_halo_swap call field_buffer_to_halo(vortg(:, :, :, I_Y), IDX_VOR_Y_SWAP, .false.) call field_buffer_to_halo(vortg(:, :, :, I_Z), IDX_VOR_Z_SWAP, .false.) call field_buffer_to_halo(tbuoyg, IDX_TBUOY_SWAP, .false.) -#ifndef ENABLE_DRY_MODE - call field_buffer_to_halo(dbuoyg, IDX_DBUOY_SWAP, .false.) - call field_buffer_to_halo(humg, IDX_HUM_SWAP, .false.) -#endif call field_mpi_dealloc end subroutine par2grid_halo_swap + subroutine par2grid_halo_swap_diag + ! we must first fill the interior grid points + ! correctly, and then the halo; otherwise + ! halo grid points do not have correct values at + ! corners where multiple processes share grid points. + + call field_mpi_alloc(n_field_swap_diag) + + !------------------------------------------------------------------ + ! Accumulate interior: + call field_halo_to_buffer(thetag, IDX_THETA_SWAP) +#ifndef ENABLE_DRY_MODE + call field_halo_to_buffer(qvg, IDX_QV_SWAP) + call field_halo_to_buffer(qlg, IDX_QL_SWAP) +#endif + ! send halo data to valid regions of other processes + call halo_to_interior_communication + + ! accumulate interior; after this operation + ! all interior grid points have the correct value + call field_buffer_to_interior(thetag, IDX_THETA_SWAP, .true.) + call field_buffer_to_interior(qvg, IDX_QV_SWAP, .true.) + call field_buffer_to_interior(qlg, IDX_QL_SWAP, .true.) +#endif + + !------------------------------------------------------------------ + ! Fill halo: + + call field_interior_to_buffer(thetag, IDX_THETA_SWAP) +#ifndef ENABLE_DRY_MODE + call field_interior_to_buffer(qvg, IDX_QV_SWAP) + call field_interior_to_buffer(qlg, IDX_QL_SWAP) +#endif + + call interior_to_halo_communication + + call field_buffer_to_halo(thetag, IDX_THETA_SWAP, .false.) +#ifndef ENABLE_DRY_MODE + call field_buffer_to_halo(qvg, IDX_QV_SWAP, .false.) + call field_buffer_to_halo(qlg, IDX_QL_SWAP, .false.) +#endif + call field_mpi_dealloc + + end subroutine par2grid_halo_swap_diag ! Interpolate the gridded quantities to the parcels ! @param[in] add contributions, i.e. do not reset parcel quantities to zero before doing grid2par. @@ -484,5 +574,52 @@ pure subroutine get_weights(xyz, i, j, k, ww) end subroutine get_weights + subroutine saturation_adjustment + double precision, parameter :: tk0c = 273.15 &! Temperature of freezing in Kelvin + double precision, parameter :: qsa1 = 3.8, &! Top in equation to calculate qsat + double precision, parameter :: qsa2 = -17.2693882, &! Constant in qsat equation + double precision, parameter :: qsa3 = 35.86, &! Constant in qsat equation + double precision, parameter :: qsa4 = 6.109, &! Constant in qsat equation + double precision, parameter :: pressure_scale_height = 8000. ! Constant in qsat equation + double precision, parameter :: ref_press = 100000 + double precision :: press, exn, temp, temp_low, qsat_low, qt_start, ql_start, ql_iter, temp_start, qsat + integer :: n, iter + + !$omp parallel default(shared) + !$omp do private(n, iter, press, exn, temp, temp_low, qsat_low, qt_start, ql_start, ql_iter, temp_start, qsat) + do n = 1, n_parcels + press=ref_press*exp(-parcels%position(3, n)/pressure_scale_height) + exn=(press/ref_press)**(r_d/c_p) + temp=parcels%theta(n)*exn + temp_start=temp + ql_start=parcels%ql(n) + qt_start=ql_start+parcels%qv(n) + ! Test unsaturated case first + temp_low=temp-(L_v/c_p)*ql_start + qsat_low = qsa1/(press*exp(qsa2*(temp_low - tk0c)/(temp_low - qsa3)) - qsa4) + if(qt_start < qsat_low) then ! Evaporate everything, if needed at all + if(ql_start>0.) then + parcels%theta(n)=parcels%theta(n)-(L_v/(c_p*exn))*ql_start + parcels%qv(n)=parcels%qv(n)+ql_start + parcels%ql(n)=0. + end if + ! Moist case: iterate a few times, start from temp instead of temp_low + else + do iter=1,4 + qsat=qsa1/(press*exp(qsa2*(temp - tk0c)/(temp - qsa3)) - qsa4) + ql_iter=max(qt_start-qsat,0.0) + temp=temp_start-(L_v/c_p)*(ql_start-ql_iter) + enddo + qsat=qsa1/(press*exp(qsa2*(temp - tk0c)/(temp - qsa3)) - qsa4) + ql_iter=max(qt_start-qsat,0.0) + parcels%theta(n)=parcels%theta(n)-(L_v/(c_p*exn))*(ql_start-ql_iter) + parcels%qv(n)=qt_start-ql_iter + parcels%ql(n)=ql_iter + end if + end do + !$omp end do + !$omp end parallel + + end subroutine saturation_adjustment end module parcel_interpl diff --git a/src/utils/physics.f90 b/src/utils/physics.f90 index 35166167e..7726d0e6d 100644 --- a/src/utils/physics.f90 +++ b/src/utils/physics.f90 @@ -69,7 +69,12 @@ module physics ![m] MPIC specific, scale-height, H double precision, protected :: height_c = 1000.0d0 - ! + ![-] Molecular weight of dry air/ molecular weight of water - 1 + double precision, protected :: qv_dens_coeff= 0.608 + + ! gas constant of gry air + double precision, protected :: r_d=287.05 + ! The following quantities are calculated: ! @@ -159,7 +164,6 @@ subroutine update_physical_quantities glat = gravity * L_v / (c_p * theta_0) glati = one / glat - lambda_c = one / height_c if (l_planet_vorticity) then From 16da9d5b915330624987b5efd75f57be4e623c52 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Thu, 11 May 2023 21:30:02 +0100 Subject: [PATCH 04/17] Compiling code --- mpi-tests/test_parcel_merge_random.f90 | 6 +++--- mpi-tests/test_parcel_split_random.f90 | 6 +++--- src/3d/Makefile.am | 6 +++--- src/3d/fields/field_netcdf.f90 | 2 +- src/3d/parcels/parcel_diagnostics.f90 | 9 +++++++-- src/3d/parcels/parcel_interpl.f90 | 24 +++++++++++++----------- unit-tests/3d/test_mpi_parcel_pack.f90 | 3 ++- 7 files changed, 32 insertions(+), 24 deletions(-) diff --git a/mpi-tests/test_parcel_merge_random.f90 b/mpi-tests/test_parcel_merge_random.f90 index 3d5e1488f..d324a42e5 100644 --- a/mpi-tests/test_parcel_merge_random.f90 +++ b/mpi-tests/test_parcel_merge_random.f90 @@ -90,12 +90,12 @@ program test_parcel_merge_random call random_number(rn(3)) j = nint(n_parcels * rn(3)) + 1 parcels%volume(j) = 0.9d0 * vmin - parcels%buoyancy(j) = 1.0d0 + parcels%theta(j) = 1.0d0 endif enddo - n_merges = int(sum(parcels%buoyancy(1:n_parcels))) + n_merges = int(sum(parcels%theta(1:n_parcels))) call perform_integer_reduction(n_merges) if (comm%rank == comm%master) then @@ -117,7 +117,7 @@ program test_parcel_merge_random n_parcels = n_orig do n = 1, n_parcels parcels%volume(n) = vol - parcels%buoyancy(n) = 0.0d0 + parcels%theta(n) = 0.0d0 parcels%B(:, n) = b call random_number(rn) diff --git a/mpi-tests/test_parcel_split_random.f90 b/mpi-tests/test_parcel_split_random.f90 index b807e291f..f8feea551 100644 --- a/mpi-tests/test_parcel_split_random.f90 +++ b/mpi-tests/test_parcel_split_random.f90 @@ -87,12 +87,12 @@ program test_parcel_split_random call random_number(rn(3)) j = nint(n_parcels * rn(3)) + 1 parcels%volume(j) = 1.1d0 * vmax - parcels%buoyancy(j) = 1.0d0 + parcels%theta(j) = 1.0d0 endif enddo - n_splits = int(sum(parcels%buoyancy(1:n_parcels))) + n_splits = int(sum(parcels%theta(1:n_parcels))) call perform_integer_reduction(n_splits) if (comm%rank == comm%master) then @@ -114,7 +114,7 @@ program test_parcel_split_random n_parcels = n_orig do n = 1, n_parcels parcels%volume(n) = vol - parcels%buoyancy(n) = 0.0d0 + parcels%theta(n) = 0.0d0 parcels%B(:, n) = b call random_number(rn) diff --git a/src/3d/Makefile.am b/src/3d/Makefile.am index 3f2ca37de..aaf6e6242 100644 --- a/src/3d/Makefile.am +++ b/src/3d/Makefile.am @@ -13,9 +13,6 @@ epic3d_SOURCES = \ fields/field_mpi.f90 \ fields/fields.f90 \ fields/field_ops.f90 \ - fields/field_diagnostics.f90 \ - fields/field_netcdf.f90 \ - fields/field_diagnostics_netcdf.f90 \ parcels/parcel_ellipsoid.f90 \ parcels/parcel_container.f90 \ parcels/parcel_mpi.f90 \ @@ -26,6 +23,9 @@ epic3d_SOURCES = \ parcels/parcel_diagnostics.f90 \ parcels/parcel_interpl.f90 \ parcels/parcel_init.f90 \ + fields/field_diagnostics.f90 \ + fields/field_netcdf.f90 \ + fields/field_diagnostics_netcdf.f90 \ inversion/inversion_utils.f90 \ inversion/inversion.f90 \ parcels/parcel_correction.f90 \ diff --git a/src/3d/fields/field_netcdf.f90 b/src/3d/fields/field_netcdf.f90 index 97e195de2..f7b9a067f 100644 --- a/src/3d/fields/field_netcdf.f90 +++ b/src/3d/fields/field_netcdf.f90 @@ -397,7 +397,7 @@ subroutine write_netcdf_fields(t) call write_netcdf_dataset(ncid, tbuoy_id, tbuoyg(lo(3):hi(3), lo(2):hi(2), lo(1):hi(1)), & start, cnt) - call write_netcdf_dataset(ncid, theta_id, theta(lo(3):hi(3), lo(2):hi(2), lo(1):hi(1)), & + call write_netcdf_dataset(ncid, theta_id, thetag(lo(3):hi(3), lo(2):hi(2), lo(1):hi(1)), & start, cnt) #ifndef ENABLE_DRY_MODE diff --git a/src/3d/parcels/parcel_diagnostics.f90 b/src/3d/parcels/parcel_diagnostics.f90 index 8b312bff1..555c057ec 100644 --- a/src/3d/parcels/parcel_diagnostics.f90 +++ b/src/3d/parcels/parcel_diagnostics.f90 @@ -9,7 +9,7 @@ module parcel_diagnostics use parcel_split_mod, only : n_parcel_splits use parcel_merge, only : n_parcel_merges use omp_lib - use physics, only : ape_calculation + use physics, only : ape_calculation, gravity, theta_0, qv_dens_coeff use ape_density, only : ape_den use mpi_timer, only : start_timer, stop_timer use mpi_communicator @@ -70,7 +70,12 @@ subroutine calculate_parcel_diagnostics vel = parcels%delta_pos(:, n) vor = parcels%vorticity(:, n) vol = parcels%volume(n) - b = parcels%buoyancy(n) +#ifndef ENABLE_DRY_MODE + ! total buoyancy (including effects of latent heating) + b = gravity*(parcels%theta(n)*(1.0+qv_dens_coeff*parcels%qv(n)-parcels%ql(n))/theta_0-1.0) +#else + b = gravity*(parcels%theta(n)/theta_0-1.0) +#endif z = parcels%position(3, n) ! kinetic energy diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index 0272e584d..64fc3cc43 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -151,7 +151,7 @@ subroutine par2grid(l_reuse) #ifndef ENABLE_DRY_MODE ! total buoyancy (including effects of latent heating) - btot = gravity*(parcels%theta(n)*(1.0+qv_theta_coefficient*parcels%qv(n)-parcels%ql(n))/theta_0-1.0) + btot = gravity*(parcels%theta(n)*(1.0+qv_dens_coeff*parcels%qv(n)-parcels%ql(n))/theta_0-1.0) #else btot = gravity*(parcels%theta(n)/theta_0-1.0) #endif @@ -264,23 +264,24 @@ end subroutine par2grid ! - nparg, that is the number of parcels per grid cell ! - nsparg, that is the number of small parcels per grid cell ! @pre The parcel must be assigned to the correct MPI process. - subroutine par2grid_diag + subroutine par2grid_diag(l_reuse) + logical, optional :: l_reuse double precision :: points(3, 4) - integer :: n, p, l, i, j, k + integer :: n, p, i, j, k double precision :: pvol, weight(0:1,0:1,0:1) thetag = zero #ifndef ENABLE_DRY_MODE - qcg = zero qvg = zero + qlg = zero #endif !$omp parallel default(shared) #ifndef ENABLE_DRY_MODE - !$omp do private(n, p, l, i, j, k, points, pvol, weight) & + !$omp do private(n, p, i, j, k, points, pvol, weight) & !$omp& private( is, js, ks, weights) & !$omp& reduction(+:nparg, nsparg, vortg, thetag, qvg, qlg) #else - !$omp do private(n, p, l, i, j, k, points, pvol, weight) & + !$omp do private(n, p, i, j, k, points, pvol, weight) & !$omp& private( is, js, ks, weights) & !$omp& reduction(+:nparg, nsparg, vortg, thetag) #endif @@ -421,6 +422,7 @@ subroutine par2grid_halo_swap_diag ! accumulate interior; after this operation ! all interior grid points have the correct value call field_buffer_to_interior(thetag, IDX_THETA_SWAP, .true.) +#ifndef ENABLE_DRY_MODE call field_buffer_to_interior(qvg, IDX_QV_SWAP, .true.) call field_buffer_to_interior(qlg, IDX_QL_SWAP, .true.) #endif @@ -575,11 +577,11 @@ pure subroutine get_weights(xyz, i, j, k, ww) end subroutine get_weights subroutine saturation_adjustment - double precision, parameter :: tk0c = 273.15 &! Temperature of freezing in Kelvin - double precision, parameter :: qsa1 = 3.8, &! Top in equation to calculate qsat - double precision, parameter :: qsa2 = -17.2693882, &! Constant in qsat equation - double precision, parameter :: qsa3 = 35.86, &! Constant in qsat equation - double precision, parameter :: qsa4 = 6.109, &! Constant in qsat equation + double precision, parameter :: tk0c = 273.15 ! Temperature of freezing in Kelvin + double precision, parameter :: qsa1 = 3.8 ! Top in equation to calculate qsat + double precision, parameter :: qsa2 = -17.2693882 ! Constant in qsat equation + double precision, parameter :: qsa3 = 35.86 ! Constant in qsat equation + double precision, parameter :: qsa4 = 6.109 ! Constant in qsat equation double precision, parameter :: pressure_scale_height = 8000. ! Constant in qsat equation double precision, parameter :: ref_press = 100000 double precision :: press, exn, temp, temp_low, qsat_low, qt_start, ql_start, ql_iter, temp_start, qsat diff --git a/unit-tests/3d/test_mpi_parcel_pack.f90 b/unit-tests/3d/test_mpi_parcel_pack.f90 index 58c94adc8..b96216b20 100644 --- a/unit-tests/3d/test_mpi_parcel_pack.f90 +++ b/unit-tests/3d/test_mpi_parcel_pack.f90 @@ -50,7 +50,8 @@ program test_mpi_parcel_pack parcels%volume(n) = 12.0d0 + a parcels%theta(n) = 13.0d0 + a #ifndef ENABLE_DRY_MODE - parcels%humidity(n) = 14.0d0 + a + parcels%qv(n) = 14.0d0 + a + parcels%ql(n) = 15.0d0 + a #endif enddo From 54189c1482f66d215648b46ad151d8cbac3bc921 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Thu, 11 May 2023 21:53:43 +0100 Subject: [PATCH 05/17] Small fixes in openmp --- src/3d/parcels/parcel_interpl.f90 | 8 ++++---- src/3d/parcels/parcel_split.f90 | 4 ++-- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index 64fc3cc43..1d2fc83ca 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -279,11 +279,11 @@ subroutine par2grid_diag(l_reuse) #ifndef ENABLE_DRY_MODE !$omp do private(n, p, i, j, k, points, pvol, weight) & !$omp& private( is, js, ks, weights) & - !$omp& reduction(+:nparg, nsparg, vortg, thetag, qvg, qlg) + !$omp& reduction(+:nparg, nsparg, thetag, qvg, qlg) #else !$omp do private(n, p, i, j, k, points, pvol, weight) & !$omp& private( is, js, ks, weights) & - !$omp& reduction(+:nparg, nsparg, vortg, thetag) + !$omp& reduction(+:nparg, nsparg, thetag) #endif do n = 1, n_parcels @@ -582,8 +582,8 @@ subroutine saturation_adjustment double precision, parameter :: qsa2 = -17.2693882 ! Constant in qsat equation double precision, parameter :: qsa3 = 35.86 ! Constant in qsat equation double precision, parameter :: qsa4 = 6.109 ! Constant in qsat equation - double precision, parameter :: pressure_scale_height = 8000. ! Constant in qsat equation - double precision, parameter :: ref_press = 100000 + double precision, parameter :: pressure_scale_height = 8000.0 ! Constant in qsat equation + double precision, parameter :: ref_press = 100000.0 double precision :: press, exn, temp, temp_low, qsat_low, qt_start, ql_start, ql_iter, temp_start, qsat integer :: n, iter diff --git a/src/3d/parcels/parcel_split.f90 b/src/3d/parcels/parcel_split.f90 index ee727230b..b7d7699fe 100644 --- a/src/3d/parcels/parcel_split.f90 +++ b/src/3d/parcels/parcel_split.f90 @@ -105,11 +105,11 @@ subroutine parcel_split(parcels, threshold) call apply_reflective_bc(parcels%position(:, n), parcels%B(:, n)) ! save parcel indices of child parcels for the -!~ ! halo swap routine + ! halo swap routine pid(n) = n pid(n_thread_loc) = n_thread_loc enddo - !$omp end dosrc/2d/parcels/parcel_diagnostics.f90: + !$omp end do !$omp end parallel n_parcel_splits = n_parcel_splits + n_parcels - last_index From 7977b2a23ea7f1e4e5148adb2adb34f7ad36d88e Mon Sep 17 00:00:00 2001 From: sjboeing Date: Thu, 11 May 2023 22:41:29 +0100 Subject: [PATCH 06/17] Restoring double correction factor, some name fixes --- src/3d/boundary_layer/bndry_fluxes.f90 | 10 +++++----- src/3d/parcels/parcel_interpl.f90 | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/3d/boundary_layer/bndry_fluxes.f90 b/src/3d/boundary_layer/bndry_fluxes.f90 index c296ccb4b..0225fc6a6 100644 --- a/src/3d/boundary_layer/bndry_fluxes.f90 +++ b/src/3d/boundary_layer/bndry_fluxes.f90 @@ -25,9 +25,9 @@ module bndry_fluxes private ! Spatial form of the buoyancy and humidity fluxes through lower surface: - double precision, dimension(:, :), allocatable :: bflux + double precision, dimension(:, :), allocatable :: thetaflux #ifndef ENABLE_DRY_MODE - double precision, dimension(:, :), allocatable :: hflux + double precision, dimension(:, :), allocatable :: qvflux #endif double precision :: fluxfac @@ -86,9 +86,9 @@ subroutine bndry_fluxes_default call bndry_fluxes_allocate - bflux = zero + thetaflux = zero #ifndef ENABLE_DRY_MODE - hflux = zero + qvflux = zero #endif end subroutine bndry_fluxes_default @@ -183,7 +183,7 @@ subroutine apply_bndry_fluxes(dt) call bilinear(xy, is, js, weights) - fac = 2.0*(zdepth - z)/(dx(3)*dx(3)) * dt + fac = (zdepth - z) * dt do l = 1, ngp ! The multiplication by dt is necessary to provide the amount of b or h diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index c4029f63c..f6eda039d 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -20,7 +20,7 @@ module parcel_interpl , interior_to_halo_communication & , halo_to_interior_communication & , field_halo_swap_scalar - , + use physics, only : gravity, theta_0, qv_dens_coeff, r_d, c_p, L_v use omp_lib use mpi_utils, only : mpi_exit_on_error @@ -409,7 +409,7 @@ subroutine par2grid_halo_swap_diag ! halo grid points do not have correct values at ! corners where multiple processes share grid points. - call field_mpi_alloc(n_field_swap_diag) + call field_mpi_alloc(n_field_swap_diag, ndim=3) !------------------------------------------------------------------ ! Accumulate interior: From 667d52e7bea2f6a4fc7df410610abdddf8b07c45 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Fri, 12 May 2023 14:50:01 +0100 Subject: [PATCH 07/17] Adding back multiplication at boundaries in par2grid_diag (needed) --- src/3d/parcels/parcel_interpl.f90 | 71 ++++++++++++++++++++++++++++--- 1 file changed, 65 insertions(+), 6 deletions(-) diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index 1d2fc83ca..1f6d84f4b 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -19,7 +19,8 @@ module parcel_interpl , field_interior_to_buffer & , interior_to_halo_communication & , halo_to_interior_communication & - , field_halo_swap + , field_halo_swap_scalar + use physics, only : gravity, theta_0, qv_dens_coeff, r_d, c_p, L_v use omp_lib use mpi_utils, only : mpi_exit_on_error @@ -62,8 +63,9 @@ module parcel_interpl , grid2par & , par2grid_timer & , grid2par_timer & - , halo_swap_timer & - , trilinear + , halo_swap_timer & + , trilinear & + , bilinear contains @@ -100,7 +102,7 @@ subroutine vol2grid(l_reuse) !$omp end parallel call start_timer(halo_swap_timer) - call field_halo_swap(volg) + call field_halo_swap_scalar(volg) call stop_timer(halo_swap_timer) ! apply free slip boundary condition @@ -320,6 +322,25 @@ subroutine par2grid_diag(l_reuse) !$omp end do !$omp end parallel + !$omp parallel workshare + ! apply free slip boundary condition + thetag(0, :, :) = two * thetag(0, :, :) + thetag(nz, :, :) = two * thetag(nz, :, :) + thetag(1, :, :) = thetag(1, :, :) + thetag(-1, :, :) + thetag(nz-1, :, :) = thetag(nz-1, :, :) + thetag(nz+1, :, :) + + qvg(0, :, :) = two * qvg(0, :, :) + qvg(nz, :, :) = two * qvg(nz, :, :) + qvg(1, :, :) = qvg(1, :, :) + qvg(-1, :, :) + qvg(nz-1, :, :) = qvg(nz-1, :, :) + qvg(nz+1, :, :) + + qlg(0, :, :) = two * qlg(0, :, :) + qlg(nz, :, :) = two * qlg(nz, :, :) + qlg(1, :, :) = qlg(1, :, :) + qlg(-1, :, :) + qlg(nz-1, :, :) = qlg(nz-1, :, :) + qlg(nz+1, :, :) + !$omp end parallel workshare + + call start_timer(halo_swap_timer) call par2grid_halo_swap_diag call stop_timer(halo_swap_timer) @@ -357,7 +378,7 @@ subroutine par2grid_halo_swap ! halo grid points do not have correct values at ! corners where multiple processes share grid points. - call field_mpi_alloc(n_field_swap) + call field_mpi_alloc(n_field_swap, ndim=3) !------------------------------------------------------------------ ! Accumulate interior: @@ -407,7 +428,7 @@ subroutine par2grid_halo_swap_diag ! halo grid points do not have correct values at ! corners where multiple processes share grid points. - call field_mpi_alloc(n_field_swap_diag) + call field_mpi_alloc(n_field_swap_diag, ndim=3) !------------------------------------------------------------------ ! Accumulate interior: @@ -624,4 +645,42 @@ subroutine saturation_adjustment end subroutine saturation_adjustment + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + ! Bi-linear interpolation + ! @param[in] pos position of the parcel + ! @param[out] ii horizontal grid points for interoplation + ! @param[out] jj meridional grid points for interpolation + ! @param[out] ww interpolation weights + subroutine bilinear(pos, ii, jj, ww) + double precision, intent(in) :: pos(2) + integer, intent(out) :: ii(4), jj(4) + double precision, intent(out) :: ww(4) + double precision :: xy(2) + + ! (i, j) + call get_horizontal_index(pos, ii(1), jj(1)) + call get_horizontal_position(ii(1), jj(1), xy) + ww(1) = product(one - abs(pos - xy) * dxi(1:2)) + + ! (i+1, j) + ii(2) = ii(1) + 1 + jj(2) = jj(1) + call get_horizontal_position(ii(2), jj(2), xy) + ww(2) = product(one - abs(pos - xy) * dxi(1:2)) + + ! (i, j+1) + ii(3) = ii(1) + jj(3) = jj(1) + 1 + call get_horizontal_position(ii(3), jj(3), xy) + ww(3) = product(one - abs(pos - xy) * dxi(1:2)) + + ! (i+1, j+1) + ii(4) = ii(2) + jj(4) = jj(3) + call get_horizontal_position(ii(4), jj(4), xy) + ww(4) = product(one - abs(pos - xy) * dxi(1:2)) + + end subroutine bilinear + end module parcel_interpl From 5b2d499b50669d4a72647e22f98629fd8ecd4c3d Mon Sep 17 00:00:00 2001 From: sjboeing Date: Fri, 12 May 2023 15:00:05 +0100 Subject: [PATCH 08/17] Adding in pre-compiler directives for dry case --- src/3d/parcels/parcel_interpl.f90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index 1f6d84f4b..58f0b6724 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -329,6 +329,7 @@ subroutine par2grid_diag(l_reuse) thetag(1, :, :) = thetag(1, :, :) + thetag(-1, :, :) thetag(nz-1, :, :) = thetag(nz-1, :, :) + thetag(nz+1, :, :) +#ifndef ENABLE_DRY_MODE qvg(0, :, :) = two * qvg(0, :, :) qvg(nz, :, :) = two * qvg(nz, :, :) qvg(1, :, :) = qvg(1, :, :) + qvg(-1, :, :) @@ -338,6 +339,7 @@ subroutine par2grid_diag(l_reuse) qlg(nz, :, :) = two * qlg(nz, :, :) qlg(1, :, :) = qlg(1, :, :) + qlg(-1, :, :) qlg(nz-1, :, :) = qlg(nz-1, :, :) + qlg(nz+1, :, :) +#endif !$omp end parallel workshare From 14b927a6854a737e055658b3913bd1d20c5408e4 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Fri, 12 May 2023 22:54:21 +0100 Subject: [PATCH 09/17] Bugfixes in 1) Dry mode and 2) Diagnostic par2grid somehow needs volg recalculated --- src/3d/fields/fields.f90 | 2 +- src/3d/parcels/parcel_interpl.f90 | 28 +++++++++++++++++++++------- 2 files changed, 22 insertions(+), 8 deletions(-) diff --git a/src/3d/fields/fields.f90 b/src/3d/fields/fields.f90 index 6b13d8654..030f97c1d 100644 --- a/src/3d/fields/fields.f90 +++ b/src/3d/fields/fields.f90 @@ -36,10 +36,10 @@ module fields double precision, allocatable, dimension(:, :, :) :: & #ifndef ENABLE_DRY_MODE - thetag, & ! dry buoyancy (or liquid-water buoyancy) qvg, & ! humidity qlg, & ! liquid water #endif + thetag, & ! dry buoyancy (or liquid-water buoyancy) tbuoyg, & ! buoyancy #ifndef NDEBUG sym_volg, & ! symmetry volume (debug mode only) diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index 58f0b6724..01fe1566a 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -48,13 +48,13 @@ module parcel_interpl integer, parameter :: n_field_swap = 5 ! restart indices for par2grid_diag - integer, parameter :: IDX_THETA_SWAP = 1 + integer, parameter :: IDX_THETA_SWAP = 2 #ifndef ENABLE_DRY_MODE - integer, parameter :: IDX_QV_SWAP = 2 & - , IDX_QL_SWAP = 3 - integer, parameter :: n_field_swap_diag = 3 + integer, parameter :: IDX_QV_SWAP = 3 & + , IDX_QL_SWAP = 4 + integer, parameter :: n_field_swap_diag = 4 #else - integer, parameter :: n_field_swap_diag = 1 + integer, parameter :: n_field_swap_diag = 2 #endif public :: par2grid & @@ -273,6 +273,7 @@ subroutine par2grid_diag(l_reuse) double precision :: pvol, weight(0:1,0:1,0:1) thetag = zero + volg = zero #ifndef ENABLE_DRY_MODE qvg = zero qlg = zero @@ -281,11 +282,11 @@ subroutine par2grid_diag(l_reuse) #ifndef ENABLE_DRY_MODE !$omp do private(n, p, i, j, k, points, pvol, weight) & !$omp& private( is, js, ks, weights) & - !$omp& reduction(+:nparg, nsparg, thetag, qvg, qlg) + !$omp& reduction(+:nparg, nsparg, thetag, volg, qvg, qlg) #else !$omp do private(n, p, i, j, k, points, pvol, weight) & !$omp& private( is, js, ks, weights) & - !$omp& reduction(+:nparg, nsparg, thetag) + !$omp& reduction(+:nparg, nsparg, thetag, volg) #endif do n = 1, n_parcels @@ -311,6 +312,8 @@ subroutine par2grid_diag(l_reuse) thetag(ks:ks+1, js:js+1, is:is+1) = thetag(ks:ks+1, js:js+1, is:is+1) & + weight * parcels%theta(n) + volg(ks:ks+1, js:js+1, is:is+1) = volg(ks:ks+1, js:js+1, is:is+1) & + + weight #ifndef ENABLE_DRY_MODE qvg(ks:ks+1, js:js+1, is:is+1) = qvg(ks:ks+1, js:js+1, is:is+1) & + weight * parcels%qv(n) @@ -324,6 +327,11 @@ subroutine par2grid_diag(l_reuse) !$omp parallel workshare ! apply free slip boundary condition + volg(0, :, :) = two * volg(0, :, :) + volg(nz, :, :) = two * volg(nz, :, :) + volg(1, :, :) = volg(1, :, :) + volg(-1, :, :) + volg(nz-1, :, :) = volg(nz-1, :, :) + volg(nz+1, :, :) + thetag(0, :, :) = two * thetag(0, :, :) thetag(nz, :, :) = two * thetag(nz, :, :) thetag(1, :, :) = thetag(1, :, :) + thetag(-1, :, :) @@ -434,6 +442,7 @@ subroutine par2grid_halo_swap_diag !------------------------------------------------------------------ ! Accumulate interior: + call field_halo_to_buffer(volg, IDX_VOL_SWAP) call field_halo_to_buffer(thetag, IDX_THETA_SWAP) #ifndef ENABLE_DRY_MODE call field_halo_to_buffer(qvg, IDX_QV_SWAP) @@ -444,6 +453,7 @@ subroutine par2grid_halo_swap_diag ! accumulate interior; after this operation ! all interior grid points have the correct value + call field_buffer_to_interior(volg, IDX_VOL_SWAP, .true.) call field_buffer_to_interior(thetag, IDX_THETA_SWAP, .true.) #ifndef ENABLE_DRY_MODE call field_buffer_to_interior(qvg, IDX_QV_SWAP, .true.) @@ -453,6 +463,7 @@ subroutine par2grid_halo_swap_diag !------------------------------------------------------------------ ! Fill halo: + call field_interior_to_buffer(volg, IDX_VOL_SWAP) call field_interior_to_buffer(thetag, IDX_THETA_SWAP) #ifndef ENABLE_DRY_MODE call field_interior_to_buffer(qvg, IDX_QV_SWAP) @@ -461,6 +472,7 @@ subroutine par2grid_halo_swap_diag call interior_to_halo_communication + call field_buffer_to_halo(volg, IDX_VOL_SWAP, .false.) call field_buffer_to_halo(thetag, IDX_THETA_SWAP, .false.) #ifndef ENABLE_DRY_MODE call field_buffer_to_halo(qvg, IDX_QV_SWAP, .false.) @@ -610,6 +622,7 @@ subroutine saturation_adjustment double precision :: press, exn, temp, temp_low, qsat_low, qt_start, ql_start, ql_iter, temp_start, qsat integer :: n, iter +#ifndef ENABLE_DRY_MODE !$omp parallel default(shared) !$omp do private(n, iter, press, exn, temp, temp_low, qsat_low, qt_start, ql_start, ql_iter, temp_start, qsat) do n = 1, n_parcels @@ -644,6 +657,7 @@ subroutine saturation_adjustment end do !$omp end do !$omp end parallel +#endif end subroutine saturation_adjustment From e1091a74719a6881de592132a2e4d9fef3a63730 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Fri, 12 May 2023 22:54:21 +0100 Subject: [PATCH 10/17] Bugfixes in 1) Dry mode and 2) Diagnostic par2grid somehow needs volg recalculated --- src/3d/fields/fields.f90 | 2 +- src/3d/parcels/parcel_interpl.f90 | 28 +++++++++++++++++++++------- 2 files changed, 22 insertions(+), 8 deletions(-) diff --git a/src/3d/fields/fields.f90 b/src/3d/fields/fields.f90 index f54ba2f3a..2dc61d804 100644 --- a/src/3d/fields/fields.f90 +++ b/src/3d/fields/fields.f90 @@ -36,10 +36,10 @@ module fields double precision, allocatable, dimension(:, :, :) :: & #ifndef ENABLE_DRY_MODE - thetag, & ! dry buoyancy (or liquid-water buoyancy) qvg, & ! humidity qlg, & ! liquid water #endif + thetag, & ! dry buoyancy (or liquid-water buoyancy) tbuoyg, & ! buoyancy #ifndef NDEBUG sym_volg, & ! symmetry volume (debug mode only) diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index 58f0b6724..01fe1566a 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -48,13 +48,13 @@ module parcel_interpl integer, parameter :: n_field_swap = 5 ! restart indices for par2grid_diag - integer, parameter :: IDX_THETA_SWAP = 1 + integer, parameter :: IDX_THETA_SWAP = 2 #ifndef ENABLE_DRY_MODE - integer, parameter :: IDX_QV_SWAP = 2 & - , IDX_QL_SWAP = 3 - integer, parameter :: n_field_swap_diag = 3 + integer, parameter :: IDX_QV_SWAP = 3 & + , IDX_QL_SWAP = 4 + integer, parameter :: n_field_swap_diag = 4 #else - integer, parameter :: n_field_swap_diag = 1 + integer, parameter :: n_field_swap_diag = 2 #endif public :: par2grid & @@ -273,6 +273,7 @@ subroutine par2grid_diag(l_reuse) double precision :: pvol, weight(0:1,0:1,0:1) thetag = zero + volg = zero #ifndef ENABLE_DRY_MODE qvg = zero qlg = zero @@ -281,11 +282,11 @@ subroutine par2grid_diag(l_reuse) #ifndef ENABLE_DRY_MODE !$omp do private(n, p, i, j, k, points, pvol, weight) & !$omp& private( is, js, ks, weights) & - !$omp& reduction(+:nparg, nsparg, thetag, qvg, qlg) + !$omp& reduction(+:nparg, nsparg, thetag, volg, qvg, qlg) #else !$omp do private(n, p, i, j, k, points, pvol, weight) & !$omp& private( is, js, ks, weights) & - !$omp& reduction(+:nparg, nsparg, thetag) + !$omp& reduction(+:nparg, nsparg, thetag, volg) #endif do n = 1, n_parcels @@ -311,6 +312,8 @@ subroutine par2grid_diag(l_reuse) thetag(ks:ks+1, js:js+1, is:is+1) = thetag(ks:ks+1, js:js+1, is:is+1) & + weight * parcels%theta(n) + volg(ks:ks+1, js:js+1, is:is+1) = volg(ks:ks+1, js:js+1, is:is+1) & + + weight #ifndef ENABLE_DRY_MODE qvg(ks:ks+1, js:js+1, is:is+1) = qvg(ks:ks+1, js:js+1, is:is+1) & + weight * parcels%qv(n) @@ -324,6 +327,11 @@ subroutine par2grid_diag(l_reuse) !$omp parallel workshare ! apply free slip boundary condition + volg(0, :, :) = two * volg(0, :, :) + volg(nz, :, :) = two * volg(nz, :, :) + volg(1, :, :) = volg(1, :, :) + volg(-1, :, :) + volg(nz-1, :, :) = volg(nz-1, :, :) + volg(nz+1, :, :) + thetag(0, :, :) = two * thetag(0, :, :) thetag(nz, :, :) = two * thetag(nz, :, :) thetag(1, :, :) = thetag(1, :, :) + thetag(-1, :, :) @@ -434,6 +442,7 @@ subroutine par2grid_halo_swap_diag !------------------------------------------------------------------ ! Accumulate interior: + call field_halo_to_buffer(volg, IDX_VOL_SWAP) call field_halo_to_buffer(thetag, IDX_THETA_SWAP) #ifndef ENABLE_DRY_MODE call field_halo_to_buffer(qvg, IDX_QV_SWAP) @@ -444,6 +453,7 @@ subroutine par2grid_halo_swap_diag ! accumulate interior; after this operation ! all interior grid points have the correct value + call field_buffer_to_interior(volg, IDX_VOL_SWAP, .true.) call field_buffer_to_interior(thetag, IDX_THETA_SWAP, .true.) #ifndef ENABLE_DRY_MODE call field_buffer_to_interior(qvg, IDX_QV_SWAP, .true.) @@ -453,6 +463,7 @@ subroutine par2grid_halo_swap_diag !------------------------------------------------------------------ ! Fill halo: + call field_interior_to_buffer(volg, IDX_VOL_SWAP) call field_interior_to_buffer(thetag, IDX_THETA_SWAP) #ifndef ENABLE_DRY_MODE call field_interior_to_buffer(qvg, IDX_QV_SWAP) @@ -461,6 +472,7 @@ subroutine par2grid_halo_swap_diag call interior_to_halo_communication + call field_buffer_to_halo(volg, IDX_VOL_SWAP, .false.) call field_buffer_to_halo(thetag, IDX_THETA_SWAP, .false.) #ifndef ENABLE_DRY_MODE call field_buffer_to_halo(qvg, IDX_QV_SWAP, .false.) @@ -610,6 +622,7 @@ subroutine saturation_adjustment double precision :: press, exn, temp, temp_low, qsat_low, qt_start, ql_start, ql_iter, temp_start, qsat integer :: n, iter +#ifndef ENABLE_DRY_MODE !$omp parallel default(shared) !$omp do private(n, iter, press, exn, temp, temp_low, qsat_low, qt_start, ql_start, ql_iter, temp_start, qsat) do n = 1, n_parcels @@ -644,6 +657,7 @@ subroutine saturation_adjustment end do !$omp end do !$omp end parallel +#endif end subroutine saturation_adjustment From c6b5dff8b6001d344a0199c361ea7aefe9e54d5b Mon Sep 17 00:00:00 2001 From: sjboeing Date: Sat, 13 May 2023 08:09:37 +0100 Subject: [PATCH 11/17] Surface fluxes setup added --- python-scripts/surface_fluxes.py | 125 +++++++++++++++++++++++++++++++ 1 file changed, 125 insertions(+) create mode 100644 python-scripts/surface_fluxes.py diff --git a/python-scripts/surface_fluxes.py b/python-scripts/surface_fluxes.py new file mode 100644 index 000000000..e4b786f5f --- /dev/null +++ b/python-scripts/surface_fluxes.py @@ -0,0 +1,125 @@ +#!/usr/bin/env python +from tools.nc_fields import nc_fields +import numpy as np +import argparse + +#The surface flux problem has few parameters, and it is not worth trying to make an exact comparison with Sam's results because he made some odd choices for some of the parameters. The cleanest thing you can do is start with b = z and omega = 0, with f = 0 (rotation would be interesting, later). Then apply a surface flux of 1- units do not matter and db/dz and the flux set the time and length scales to be O(1). So, make the domain big enough, maybe 4L x 4L x L (in x, y and z) with L > 1, maybe 5. Use a uniform flux but slightly disturb the initial parcel positions - this is enough. I'll ask Sam to send his thesis to us. + +try: + parser = argparse.ArgumentParser( + description="Create surface flux setup" + ) + + parser.add_argument( + "--nx", + type=int, + required=False, + default=256, + help="number of cells in x", + ) + + parser.add_argument( + "--ny", + type=int, + required=False, + default=256, + help="number of cells in y", + ) + + parser.add_argument( + "--nz", + type=int, + required=False, + default=64, + help="number of cells in z", + ) + + parser.add_argument( + "--L", + type=float, + required=False, + default=3200.0, + help="domain extent", + ) + + args = parser.parse_args() + + lapse=0.003 + thetafluxmag=0.1 + + # number of cells + nx = args.nx + ny = args.ny + nz = args.nz + + L = args.L + + ncf = nc_fields() + + ncf.open('boundary_layer_setup' + str(nx) + 'x' + str(ny) + 'x' + str(nz) + '.nc') + + # domain origin + origin = (-2.0 * L, -2.0 * L, 0.0) + + # domain extent + extent = (4.0 * L, 4.0 * L, L) + + + # mesh spacings + dx = extent[0] / nx + dy = extent[1] / ny + dz = extent[2] / nz + + buoy = np.zeros((nz+1, ny, nx)) + + # ranges from 0 to nz + for k in range(nz+1): + zrel = k * dz + buoy[k, :, :] = 300.+lapse*zrel + + # add perturbation + rng = np.random.default_rng(seed=42) + noise = rng.uniform(low=0.0, high=0.01*(buoy.max()-300.), size=(nz+1, ny, nx)) + buoy += noise + + # write all provided fields + ncf.add_field('theta', buoy, unit='K', long_name='potential temperature') + + ncf.add_box(origin, extent, [nx, ny, nz]) + + ncf.close() + + # + # Setup surface fluxes: + # + + ncf_flux = nc_fields() + + ncf_flux.set_dim_names(['x', 'y']) + + ncf_flux.open('surface_flux_' + str(nx) + 'x' + str(ny) + '.nc') + + # domain origin + origin = (-2.0 * L, -2.0 * L) + + # domain extent + extent = (4.0 * L, 4.0 * L) + + + # mesh spacings + dx = extent[0] / nx + dy = extent[1] / ny + + thetaflux = np.ones((ny, nx))*thetafluxmag + qvflux = np.zeros((ny, nx)) + + ncf_flux.add_field('thetaflux', thetaflux, unit='K m/s') + ncf_flux.add_field('qvflux', qvflux, unit='kg/kg m/s') + + ncf_flux.add_box(origin, extent, [nx, ny]) + + ncf_flux.close() + +except Exception as err: + print(err) + From 02d27734e7ae4df4d12d221896458e99ea6356f1 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Sat, 13 May 2023 16:11:39 +0100 Subject: [PATCH 12/17] Fixes to run BOMEX, temporarily disable vorticity correction --- python-scripts/bomex_fluxes.py | 152 ++++++++++++++++++++++++++++++ src/3d/epic3d.f90 | 2 +- src/3d/parcels/parcel_interpl.f90 | 13 +-- 3 files changed, 160 insertions(+), 7 deletions(-) create mode 100644 python-scripts/bomex_fluxes.py diff --git a/python-scripts/bomex_fluxes.py b/python-scripts/bomex_fluxes.py new file mode 100644 index 000000000..e633fe1e2 --- /dev/null +++ b/python-scripts/bomex_fluxes.py @@ -0,0 +1,152 @@ +#!/usr/bin/env python +from tools.nc_fields import nc_fields +import numpy as np +import argparse + +try: + parser = argparse.ArgumentParser( + description="Create surface flux setup" + ) + + parser.add_argument( + "--nx", + type=int, + required=False, + default=192, + help="number of cells in x", + ) + + parser.add_argument( + "--ny", + type=int, + required=False, + default=192, + help="number of cells in y", + ) + + parser.add_argument( + "--nz", + type=int, + required=False, + default=96, + help="number of cells in z", + ) + + parser.add_argument( + "--L", + type=float, + required=False, + default=3200.0, + help="domain extent", + ) + + args = parser.parse_args() + + thetafluxmag=8.0e-3 + qvfluxmag=5.2e-5 + + # number of cells + nx = args.nx + ny = args.ny + nz = args.nz + + L = args.L + + ncf = nc_fields() + + ncf.open('bomex_setup' + str(nx) + 'x' + str(ny) + 'x' + str(nz) + '.nc') + + # domain origin + origin = (-1.0 * L, -1.0 * L, 0.0) + + # domain extent + extent = (2.0 * L, 2.0 * L, L) + + + # mesh spacings + dx = extent[0] / nx + dy = extent[1] / ny + dz = extent[2] / nz + + theta = np.zeros((nz+1, ny, nx)) + qv = np.zeros((nz+1, ny, nx)) + yvort = np.zeros((nz+1, ny, nx)) + novort = np.zeros((nz+1, ny, nx)) + + # ranges from 0 to nz + for k in range(nz+1): + zz = k * dz + if(zz < 520.): + theta[k, :,: ] = 298.7 + elif(zz < 1480): + theta[k, :, :] = 298.7 + (zz-520.)*(302.4-298.7)/(1480.-520.) + elif(zz < 2000): + theta[k, :, :] = 302.4 + (zz-1480.)*(308.2-302.4)/(2000.-1480.) + else: + theta[k, :, :] = 308.2 + (zz-2000.)*(311.85-308.2)/(3000.-2000.) + + # specific humidity + if(zz < 520.): + qv[k, :, :] = 1e-3*(17.0 + zz*(16.3-17.0)/520.) + elif(zz < 1480): + qv[k, :, :] = 1.e-3*(16.3 + (zz-520.)*(10.7-16.3)/(1480.-520.)) + elif(zz < 2000): + qv[k, :, :] = 1.e-3*(10.7 + (zz-1480.)*(4.2-10.7)/(2000.-1480.)) + else: + qv[k, :, :] = 1.e-3*(4.2 + (zz-2000.)*(3.-4.2)/(3000.-2000.)) + + if(zz > 699.): + yvort[k, :, :] = 0.5*(-4.61+8.75)/(3000.-700.) + if(zz > 701.): + yvort[k, :, :] = (-4.61+8.75)/(3000.-700.) + + # add perturbation + rng = np.random.default_rng(seed=42) + noise = rng.uniform(low=-0.05, high=0.05, size=(nz+1, ny, nx)) + theta += noise + + # write all provided fields + ncf.add_field('x_vorticity', novort, unit='1/s') + ncf.add_field('y_vorticity', yvort, unit='1/s') + ncf.add_field('z_vorticity', novort, unit='1/s') + ncf.add_field('theta', theta, unit='K', long_name='potential temperature') + ncf.add_field('qv', qv, unit='kg/kg', long_name='water vapour spec. hum.') + + ncf.add_box(origin, extent, [nx, ny, nz]) + + ncf.close() + + # + # Setup surface fluxes: + # + + ncf_flux = nc_fields() + + ncf_flux.set_dim_names(['x', 'y']) + + ncf_flux.open('bomex_flux_' + str(nx) + 'x' + str(ny) + '.nc') + + # domain origin + origin = (-1.0 * L, -1.0 * L) + + # domain extent + extent = (2.0 * L, 2.0 * L) + + + # mesh spacings + dx = extent[0] / nx + dy = extent[1] / ny + + thetaflux = np.ones((ny, nx))*thetafluxmag + qvflux = np.ones((ny, nx))*qvfluxmag + + ncf_flux.add_field('thetaflux', thetaflux, unit='K m/s') + ncf_flux.add_field('qvflux', qvflux, unit='kg/kg m/s') + + ncf_flux.add_box(origin, extent, [nx, ny]) + + ncf_flux.close() + +except Exception as err: + print(err) + diff --git a/src/3d/epic3d.f90 b/src/3d/epic3d.f90 index 714c73bca..1b3d672aa 100644 --- a/src/3d/epic3d.f90 +++ b/src/3d/epic3d.f90 @@ -122,7 +122,7 @@ subroutine run print "(a15, f0.4)", "time: ", t endif #endif - call apply_vortcor + !call apply_vortcor call ls_rk4_step(t) diff --git a/src/3d/parcels/parcel_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index 01fe1566a..a65145e9a 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -617,8 +617,9 @@ subroutine saturation_adjustment double precision, parameter :: qsa2 = -17.2693882 ! Constant in qsat equation double precision, parameter :: qsa3 = 35.86 ! Constant in qsat equation double precision, parameter :: qsa4 = 6.109 ! Constant in qsat equation - double precision, parameter :: pressure_scale_height = 8000.0 ! Constant in qsat equation - double precision, parameter :: ref_press = 100000.0 + double precision, parameter :: pressure_scale_height = 8500.0 + double precision, parameter :: surf_press = 101500.0 + double precision, parameter :: ref_press = 100000.0 double precision :: press, exn, temp, temp_low, qsat_low, qt_start, ql_start, ql_iter, temp_start, qsat integer :: n, iter @@ -626,7 +627,7 @@ subroutine saturation_adjustment !$omp parallel default(shared) !$omp do private(n, iter, press, exn, temp, temp_low, qsat_low, qt_start, ql_start, ql_iter, temp_start, qsat) do n = 1, n_parcels - press=ref_press*exp(-parcels%position(3, n)/pressure_scale_height) + press=surf_press*exp(-parcels%position(3, n)/pressure_scale_height) exn=(press/ref_press)**(r_d/c_p) temp=parcels%theta(n)*exn temp_start=temp @@ -634,7 +635,7 @@ subroutine saturation_adjustment qt_start=ql_start+parcels%qv(n) ! Test unsaturated case first temp_low=temp-(L_v/c_p)*ql_start - qsat_low = qsa1/(press*exp(qsa2*(temp_low - tk0c)/(temp_low - qsa3)) - qsa4) + qsat_low = qsa1/(0.01*press*exp(qsa2*(temp_low - tk0c)/(temp_low - qsa3)) - qsa4) if(qt_start < qsat_low) then ! Evaporate everything, if needed at all if(ql_start>0.) then parcels%theta(n)=parcels%theta(n)-(L_v/(c_p*exn))*ql_start @@ -644,11 +645,11 @@ subroutine saturation_adjustment ! Moist case: iterate a few times, start from temp instead of temp_low else do iter=1,4 - qsat=qsa1/(press*exp(qsa2*(temp - tk0c)/(temp - qsa3)) - qsa4) + qsat=qsa1/(0.01*press*exp(qsa2*(temp - tk0c)/(temp - qsa3)) - qsa4) ql_iter=max(qt_start-qsat,0.0) temp=temp_start-(L_v/c_p)*(ql_start-ql_iter) enddo - qsat=qsa1/(press*exp(qsa2*(temp - tk0c)/(temp - qsa3)) - qsa4) + qsat=qsa1/(0.01*press*exp(qsa2*(temp - tk0c)/(temp - qsa3)) - qsa4) ql_iter=max(qt_start-qsat,0.0) parcels%theta(n)=parcels%theta(n)-(L_v/(c_p*exn))*(ql_start-ql_iter) parcels%qv(n)=qt_start-ql_iter From 41c80a94742e508d0f0cd8fd9ac78be50d8e4fb5 Mon Sep 17 00:00:00 2001 From: sjboeing Date: Sat, 13 May 2023 19:53:59 +0100 Subject: [PATCH 13/17] Fix to parcel correction --- src/3d/epic3d.f90 | 2 +- src/3d/parcels/parcel_correction.f90 | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/3d/epic3d.f90 b/src/3d/epic3d.f90 index 1b3d672aa..714c73bca 100644 --- a/src/3d/epic3d.f90 +++ b/src/3d/epic3d.f90 @@ -122,7 +122,7 @@ subroutine run print "(a15, f0.4)", "time: ", t endif #endif - !call apply_vortcor + call apply_vortcor call ls_rk4_step(t) diff --git a/src/3d/parcels/parcel_correction.f90 b/src/3d/parcels/parcel_correction.f90 index a2cac6676..5f3cebf4f 100644 --- a/src/3d/parcels/parcel_correction.f90 +++ b/src/3d/parcels/parcel_correction.f90 @@ -91,7 +91,7 @@ subroutine apply_vortcor call start_timer(vort_corr_timer) vsum = zero - dvor = - vor_bar + dvor = zero !$omp parallel default(shared) !$omp do private(n) reduction(+: dvor, vsum) @@ -121,7 +121,7 @@ subroutine apply_vortcor !$omp parallel default(shared) !$omp do private(n) do n = 1, n_parcels - parcels%vorticity(:, n) = parcels%vorticity(:, n) - dvor + parcels%vorticity(:, n) = parcels%vorticity(:, n) - (dvor - vor_bar) enddo !$omp end do !$omp end parallel From f8c16b94b98ecdde31ad3a91c591f47c792346d0 Mon Sep 17 00:00:00 2001 From: Steven Boeing Date: Sun, 30 Jul 2023 16:59:54 +0100 Subject: [PATCH 14/17] Cleanup --- src/3d/parcels/parcel_container.f90 | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/3d/parcels/parcel_container.f90 b/src/3d/parcels/parcel_container.f90 index cd51025b4..4e4ae38b7 100644 --- a/src/3d/parcels/parcel_container.f90 +++ b/src/3d/parcels/parcel_container.f90 @@ -281,9 +281,10 @@ subroutine parcel_resize(new_size) call resize_array(parcels%vorticity, new_size, n_parcels) call resize_array(parcels%B, new_size, n_parcels) call resize_array(parcels%volume, new_size, n_parcels) - call resize_array(parcels%buoyancy, new_size, n_parcels) + call resize_array(parcels%theta, new_size, n_parcels) #ifndef ENABLE_DRY_MODE - call resize_array(parcels%humidity, new_size, n_parcels) + call resize_array(parcels%qv, new_size, n_parcels) + call resize_array(parcels%ql, new_size, n_parcels) #endif call parcel_ellipsoid_resize(new_size, n_parcels) From 3314ea951cbc49cab1168f606841e65bcaa11775 Mon Sep 17 00:00:00 2001 From: Steven Boeing Date: Sun, 30 Jul 2023 17:03:26 +0100 Subject: [PATCH 15/17] Fix mistake in random merge test --- mpi-tests/test_parcel_split_random.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mpi-tests/test_parcel_split_random.f90 b/mpi-tests/test_parcel_split_random.f90 index 5d5f7e36f..ddf8552a5 100644 --- a/mpi-tests/test_parcel_split_random.f90 +++ b/mpi-tests/test_parcel_split_random.f90 @@ -87,7 +87,7 @@ program test_parcel_split_random if (rn(3) > 0.5d0) then call random_number(rn(3)) j = nint(n_parcels * rn(3)) + 1 - parcels%volume(j) = 1.1d0 * vmax + parcels%B(1, j) = 5.0d0 * parcels%B(1, j) parcels%theta(j) = 1.0d0 endif From d6e9e693a0dea3c20384840d16a2559ca6385c00 Mon Sep 17 00:00:00 2001 From: Matthias Frey Date: Mon, 31 Jul 2023 15:25:37 +0100 Subject: [PATCH 16/17] adjust amax --- src/3d/utils/parameters.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/3d/utils/parameters.f90 b/src/3d/utils/parameters.f90 index da55eec0c..f1c0a1004 100644 --- a/src/3d/utils/parameters.f90 +++ b/src/3d/utils/parameters.f90 @@ -145,7 +145,7 @@ subroutine update_parameters vmin = vcell / parcel%min_vratio - amax = minval(dx) + amax = (f34 * fpi) ** f13 * minval(dx) max_num_parcels = int(box%halo_ncell * parcel%min_vratio * parcel%size_factor) From 9c0c0080a00be1862b25df9b60c3f057db11674b Mon Sep 17 00:00:00 2001 From: sjboeing Date: Tue, 1 Aug 2023 08:27:04 +0100 Subject: [PATCH 17/17] Replace remaining binc --- src/3d/boundary_layer/bndry_fluxes.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/3d/boundary_layer/bndry_fluxes.f90 b/src/3d/boundary_layer/bndry_fluxes.f90 index 6a441e8d4..ae32c92b4 100644 --- a/src/3d/boundary_layer/bndry_fluxes.f90 +++ b/src/3d/boundary_layer/bndry_fluxes.f90 @@ -187,8 +187,8 @@ subroutine bndry_fluxes_time_step(dt) return endif - ! local maximum of absolute value (units: m/s**3) - abs_max = maxval(dabs(binc(box%lo(2):box%hi(2), box%lo(1):box%hi(1)))) + ! local maximum of absolute value + abs_max = maxval(dabs(thetaflux(box%lo(2):box%hi(2), box%lo(1):box%hi(1)))) ! get global abs_max call MPI_Allreduce(MPI_IN_PLACE, &