diff --git a/convention.txt b/convention.txt index 1ea29ad9..d79d67cd 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/examples/fomex.config b/examples/fomex.config new file mode 100755 index 00000000..48a7e374 --- /dev/null +++ b/examples/fomex.config @@ -0,0 +1,40 @@ +``` +&EPIC + + field_file = 'bomex_setup192x192x96.nc' ! input field file + flux_file = 'bomex_flux_192x192.nc' ! input field file + + ! + ! output info + ! + output%field_freq = 300 ![s] write after these many seconds to the field NetCDF file + output%parcel_freq = 3600 ![s] write after these many seconds to the parcel NetCDF file + output%parcel_stats_freq = 100 ![s] write after these many seconds to parcel stats NetCDF file + output%field_stats_freq = 100 ![s] write after these many seconds to the field stats NetCDF file + output%write_fields = .true. ! enable / disable field dump + output%write_parcels = .true. ! enable / disable parcel dump + output%write_parcel_stats = .true. ! enable / disable parcel statistics + output%write_field_stats = .true. ! enable / disable field statistics + output%overwrite = .true. ! replace existing NetCDF files + output%basename = 'fomex' ! NetCDF output base name + + ! + ! parcel info + ! + parcel%n_per_cell = 8 ! initial number of parcels per cell + parcel%lambda_max = 4.0 ! maximum parcel aspect ratio + parcel%min_vratio = 20.0 ! minimum ratio of grid cell volume / parcel volume + parcel%correction_iters = 2 ! how many parcel correction iterations + parcel%gradient_pref = 1.8 ! gradient correction prefactor + parcel%max_compression = 0.5 ! gradient correction maximum compression + + ! BOMEX ls forcings + forcings%l_ls_forcings = .true. ! use BOMEX large-scale forcings + + ! + ! stepper info + ! + time%limit = 21600.0 ! time limit (s) + time%alpha = 0.2 ! scaling factor for the strain and buoyancy gradient time step + time%precise_stop = .false. ! time limit exact +/ diff --git a/models/3d/moist_3d.f90 b/models/3d/moist_3d.f90 index e7946b1d..dad16e51 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/mpi-tests/parcel_merge_serial.f90 b/mpi-tests/parcel_merge_serial.f90 index f3f45e78..79b62884 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/mpi-tests/test_utils.f90 b/mpi-tests/test_utils.f90 index 5248150b..066b39b6 100644 --- a/mpi-tests/test_utils.f90 +++ b/mpi-tests/test_utils.f90 @@ -170,7 +170,7 @@ subroutine setup_parcels(xlen, ylen, zlen, l_shuffle, l_variable_nppc) parcels%vorticity(3, l) = 20.0d0 * rn(6) - 10.d0 ! buoyancy between -1 and 1: y = 2 * x - 1 - parcels%buoyancy(l) = 2.0d0 * rn(7) - 1.d0 + parcels%theta(l) = 2.0d0 * rn(7) - 1.d0 if ((x >= xlo) .and. (x <= xhi) .and. & (y >= ylo) .and. (y <= yhi) .and. & @@ -242,9 +242,9 @@ subroutine shuffleall parcels%vorticity(:, rand_target) = parcels%vorticity(:, shuffle_index) parcels%vorticity(:, shuffle_index) = tmp_vec - tmp_var = parcels%buoyancy(rand_target) - parcels%buoyancy(rand_target) = parcels%buoyancy(shuffle_index) - parcels%buoyancy(shuffle_index) = tmp_var + tmp_var = parcels%theta(rand_target) + parcels%theta(rand_target) = parcels%theta(shuffle_index) + parcels%theta(shuffle_index) = tmp_var tmp_var = parcels%volume(rand_target) parcels%volume(rand_target) = parcels%volume(shuffle_index) diff --git a/python-scripts/bomex_with_fluxes.py b/python-scripts/bomex_with_fluxes.py new file mode 100644 index 00000000..289c00a4 --- /dev/null +++ b/python-scripts/bomex_with_fluxes.py @@ -0,0 +1,164 @@ +#!/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)) + xvort = 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.0: + theta[k, :, :] = 298.7 + elif zz < 1480: + theta[k, :, :] = 298.7 + (zz - 520.0) * (302.4 - 298.7) / (1480.0 - 520.0) + elif zz < 2000: + theta[k, :, :] = 302.4 + (zz - 1480.0) * (308.2 - 302.4) / (2000.0 - 1480.0) + else: + theta[k, :, :] = 308.2 + (zz - 2000.0) * (311.85 - 308.2) / ( + 3000.0 - 2000.0 + ) + + # specific humidity + if zz < 520.0: + qv[k, :, :] = 1e-3 * (17.0 + zz * (16.3 - 17.0) / 520.0) + elif zz < 1480: + qv[k, :, :] = 1.0e-3 * ( + 16.3 + (zz - 520.0) * (10.7 - 16.3) / (1480.0 - 520.0) + ) + elif zz < 2000: + qv[k, :, :] = 1.0e-3 * ( + 10.7 + (zz - 1480.0) * (4.2 - 10.7) / (2000.0 - 1480.0) + ) + else: + qv[k, :, :] = 1.0e-3 * ( + 4.2 + (zz - 2000.0) * (3.0 - 4.2) / (3000.0 - 2000.0) + ) + + if zz < 300.0: + xvort[k, :, :] = 0.0043 * np.cos(np.pi * (zz + 600) / 1200.0) + 0.0023 + yvort[k, :, :] = 0.0062 * np.cos((np.pi * (zz - 300) / 600.0) ** 2) - 0.0067 + elif zz < 600.0: + xvort[k, :, :] = 0.0043 * np.cos(np.pi * (zz + 600) / 1200.0) + 0.0023 + yvort[k, :, :] = 0.0027 * np.cos(np.pi * (zz - 300) / 600.0) - 0.0032 + elif zz < 1400.0: + xvort[k, :, :] = -0.001 - 0.001 * np.cos(np.pi * (zz - 600.0) / 800.0) + yvort[k, :, :] = 0.005 * np.cos(np.pi * (zz - 1400.0) / 1600.0) - 0.0032 + else: + xvort[k, :, :] = 0.0 + yvort[k, :, :] = 0.0018 + + # 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", xvort, 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/python-scripts/surface_fluxes.py b/python-scripts/surface_fluxes.py new file mode 100644 index 00000000..e4b786f5 --- /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) + diff --git a/src/3d/Makefile.am b/src/3d/Makefile.am index 62b18565..f9ca509e 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_bc.f90 \ @@ -26,9 +23,13 @@ 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 \ + parcels/parcel_ls_forcings.f90 \ parcels/parcel_netcdf.f90 \ parcels/parcel_diagnostics_netcdf.f90 \ boundary_layer/bndry_fluxes.f90 \ diff --git a/src/3d/boundary_layer/bndry_fluxes.f90 b/src/3d/boundary_layer/bndry_fluxes.f90 index 5cf8c2f0..cffff959 100644 --- a/src/3d/boundary_layer/bndry_fluxes.f90 +++ b/src/3d/boundary_layer/bndry_fluxes.f90 @@ -9,6 +9,7 @@ module bndry_fluxes use parcel_container, only : n_parcels, parcels use netcdf_reader use omp_lib + use physics, only : gravity, theta_0 use options, only : time use mpi_utils, only : mpi_stop use field_mpi, only : field_mpi_alloc & @@ -27,9 +28,9 @@ module bndry_fluxes integer :: bndry_flux_timer ! Spatial form of the buoyancy and humidity fluxes through lower surface: - double precision, dimension(:, :), allocatable :: binc + double precision, dimension(:, :), allocatable :: thetaflux #ifndef ENABLE_DRY_MODE - double precision, dimension(:, :), allocatable :: hinc + double precision, dimension(:, :), allocatable :: qvflux #endif ! Denotes height below which surface fluxes are applied: @@ -52,9 +53,9 @@ subroutine bndry_fluxes_allocate l_bndry_flux_allocated = .true. - allocate(binc(box%hlo(2):box%hhi(2), box%hlo(1):box%hhi(1))) + allocate(thetaflux(box%hlo(2):box%hhi(2), box%hlo(1):box%hhi(1))) #ifndef ENABLE_DRY_MODE - allocate(hinc(box%hlo(2):box%hhi(2), box%hlo(1):box%hhi(1))) + allocate(qvflux(box%hlo(2):box%hhi(2), box%hlo(1):box%hhi(1))) #endif end subroutine bndry_fluxes_allocate @@ -69,9 +70,9 @@ subroutine bndry_fluxes_deallocate l_bndry_flux_allocated = .false. - deallocate(binc) + deallocate(thetaflux) #ifndef ENABLE_DRY_MODE - deallocate(hinc) + deallocate(qvflux) #endif call stop_timer(bndry_flux_timer) @@ -85,9 +86,9 @@ subroutine bndry_fluxes_default call bndry_fluxes_allocate - binc = zero + thetaflux = zero #ifndef ENABLE_DRY_MODE - hinc = zero + qvflux = zero #endif call stop_timer(bndry_flux_timer) @@ -123,27 +124,27 @@ subroutine read_bndry_fluxes(fname) call bndry_fluxes_default - if (has_dataset(ncid, 'bflux')) then + if (has_dataset(ncid, 'thetaflux')) then call read_netcdf_dataset(ncid, & - 'bflux', & - binc(lo(2):hi(2), & - lo(1):hi(1)), & + 'thetaflux', & + thetaflux(lo(2):hi(2), & + lo(1):hi(1)), & start, & cnt) else - call mpi_stop("No buoyancy flux field 'bflux' found in file.") + call mpi_stop("No theta flux field 'thetaflux' found in file.") endif #ifndef ENABLE_DRY_MODE - if (has_dataset(ncid, 'hflux')) then + if (has_dataset(ncid, 'qvflux')) then call read_netcdf_dataset(ncid, & - 'hflux', & - hinc(lo(2):hi(2), & - lo(1):hi(1)), & + 'qvflux', & + qvflux(lo(2):hi(2), & + lo(1):hi(1)), & start, & cnt) else - call mpi_stop("No humidity flux field 'hflux' found in file.") + call mpi_stop("No qv flux field 'qvflux' found in file.") endif #endif @@ -164,9 +165,9 @@ subroutine read_bndry_fluxes(fname) !$omp parallel !$omp workshare - binc = two * dxi(3) * binc + thetaflux = two * dxi(3) * thetaflux #ifndef ENABLE_DRY_MODE - hinc = two * dxi(3) * hinc + qvflux = two * dxi(3) * qvflux #endif !$omp end workshare !$omp end parallel @@ -188,7 +189,7 @@ subroutine bndry_fluxes_time_step(dt) endif ! local maximum of absolute value (units: m/s**3) - abs_max = maxval(abs(binc(box%lo(2):box%hi(2), box%lo(1):box%hi(1)))) + abs_max = (gravity/theta_0)*maxval(abs(thetaflux(box%lo(2):box%hi(2), box%lo(1):box%hi(1)))) ! get global abs_max call MPI_Allreduce(MPI_IN_PLACE, & @@ -236,11 +237,11 @@ subroutine apply_bndry_fluxes(dt) ! The multiplication by dt is necessary to provide the amount of b or h ! entering through the bottom surface over a time interval of dt. - parcels%buoyancy(n) = parcels%buoyancy(n) & - + fac * sum(weights * binc(js:js+1, is:is+1)) + parcels%theta(n) = parcels%theta(n) & + + fac * sum(weights * thetaflux(js:js+1, is:is+1)) #ifndef ENABLE_DRY_MODE - parcels%humidity(n) = parcels%humidity(n) & - + fac * sum(weights * hinc(js:js+1, is:is+1)) + parcels%qv(n) = parcels%qv(n) & + + fac * sum(weights * qvflux(js:js+1, is:is+1)) #endif endif enddo @@ -261,16 +262,16 @@ subroutine bndry_fluxes_fill_halo #endif call field_mpi_alloc(nfluxes, ndim=2) - call field_interior_to_buffer(binc, 1) + call field_interior_to_buffer(thetaflux, 1) #ifndef ENABLE_DRY_MODE - call field_interior_to_buffer(hinc, 2) + call field_interior_to_buffer(qvflux, 2) #endif call interior_to_halo_communication - call field_buffer_to_halo(binc, 1, .false.) + call field_buffer_to_halo(thetaflux, 1, .false.) #ifndef ENABLE_DRY_MODE - call field_buffer_to_halo(hinc, 2, .false.) + call field_buffer_to_halo(qvflux, 2, .false.) #endif call field_mpi_dealloc diff --git a/src/3d/epic3d.f90 b/src/3d/epic3d.f90 index fbac6931..9b8930d9 100644 --- a/src/3d/epic3d.f90 +++ b/src/3d/epic3d.f90 @@ -37,6 +37,7 @@ program epic3d use inversion_utils, only : init_inversion, finalise_inversion use parcel_interpl, only : grid2par_timer, & par2grid_timer, & + saturation_adjustment, & halo_swap_timer use parcel_init, only : init_timer use ls_rk, only : ls_rk_step, rk_timer, ls_rk_setup @@ -133,7 +134,8 @@ subroutine run do while (t < time%limit) - call apply_vortcor + ! no longer apply vorticity correction for real cases + !call apply_vortcor call ls_rk_step(t) @@ -149,11 +151,14 @@ subroutine run call apply_gradient(parcel%gradient_pref, parcel%max_compression, .true.) enddo + ! call saturation adjustment again after all corrections + call saturation_adjustment + enddo - ! write final step (we only write if we really advanced in time) + ! write final step (we only write if we really advanced in time) if (t > time%initial) then - call apply_vortcor + !call apply_vortcor call write_last_step(t) endif diff --git a/src/3d/fields/field_diagnostics.f90 b/src/3d/fields/field_diagnostics.f90 index b4425736..8e5747eb 100644 --- a/src/3d/fields/field_diagnostics.f90 +++ b/src/3d/fields/field_diagnostics.f90 @@ -11,9 +11,7 @@ module field_diagnostics use mpi_collectives, only : mpi_blocking_reduce use physics, only : ape_calculation use ape_density, only : ape_den -#ifdef ENABLE_BUOYANCY_PERTURBATION_MODE - use physics, only : bfsq -#endif + use parcel_interpl, only: par2grid_diag implicit none integer :: field_stats_timer @@ -50,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 5147a6fa..40b8888a 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 495de9e4..7fb8dd69 100644 --- a/src/3d/fields/field_netcdf.f90 +++ b/src/3d/fields/field_netcdf.f90 @@ -38,17 +38,17 @@ module field_netcdf , NC_Z_VTEND = 9 & , NC_NPARG = 10 & , NC_NSPARG = 11 & - , NC_TBUOY = 12 & - , NC_VOL = 13 + , NC_THETA = 12 & + , NC_TBUOY = 13 & + , NC_VOL = 14 #ifndef ENABLE_DRY_MODE - integer, parameter :: NC_DBUOY = 14 & - , NC_HUM = 15 & - , NC_LBUOY = 16 + integer, parameter :: NC_QV = 15 & + , NC_QL = 16 type(netcdf_info) :: nc_dset(16) #else - type(netcdf_info) :: nc_dset(13) + type(netcdf_info) :: nc_dset(14) #endif integer :: n_writes @@ -237,13 +237,12 @@ subroutine write_netcdf_fields(t) #endif call write_field_double(NC_TBUOY, tbuoyg, start, cnt) + call write_field_double(NC_THETA, thetag, start, cnt) #ifndef ENABLE_DRY_MODE - call write_field_double(NC_DBUOY, dbuoyg, start, cnt) + call write_field_double(NC_QV, qvg, start, cnt) - call write_field_double(NC_LBUOY, glati * (tbuoyg - dbuoyg), start, cnt) - - call write_field_double(NC_HUM, humg, start, cnt) + call write_field_double(NC_QL, qlg, start, cnt) #endif call write_field_double(NC_VOL, volg, start, cnt) @@ -370,11 +369,31 @@ subroutine read_netcdf_fields(fname, step) cnt) endif + if (has_dataset(lid, nc_dset(NC_THETA)%name)) then + call read_netcdf_dataset(lid, & + nc_dset(NC_THETA)%name, & + thetag(box%lo(3):box%hi(3), & + box%lo(2):box%hi(2), & + box%lo(1):box%hi(1)), & + start, & + cnt) + endif + #ifndef ENABLE_DRY_MODE - if (has_dataset(lid, nc_dset(NC_HUM)%name)) then + if (has_dataset(lid, nc_dset(NC_QV)%name)) then + call read_netcdf_dataset(lid, & + nc_dset(NC_QV)%name, & + 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, nc_dset(NC_QL)%name)) then call read_netcdf_dataset(lid, & - nc_dset(NC_HUM)%name, & - humg(box%lo(3):box%hi(3), & + nc_dset(NC_QL)%name, & + qlg(box%lo(3):box%hi(3), & box%lo(2):box%hi(2), & box%lo(1):box%hi(1)), & start, & @@ -402,11 +421,11 @@ subroutine set_netcdf_field_output nc_dset(NC_X_VEL)%l_enabled = .true. nc_dset(NC_Y_VEL)%l_enabled = .true. nc_dset(NC_Z_VEL)%l_enabled = .true. + nc_dset(NC_THETA)%l_enabled = .true. nc_dset(NC_TBUOY)%l_enabled = .true. #ifndef ENABLE_DRY_MODE - nc_dset(NC_DBUOY)%l_enabled = .true. - nc_dset(NC_HUM)%l_enabled = .true. - nc_dset(NC_LBUOY)%l_enabled = .true. + nc_dset(NC_QV)%l_enabled = .true. + nc_dset(NC_QL)%l_enabled = .true. #endif nc_dset(NC_VOL)%l_enabled = .true. else @@ -436,11 +455,11 @@ subroutine set_netcdf_field_output nc_dset(NC_X_VEL)%l_enabled = .true. nc_dset(NC_Y_VEL)%l_enabled = .true. nc_dset(NC_Z_VEL)%l_enabled = .true. + nc_dset(NC_THETA)%l_enabled = .true. nc_dset(NC_TBUOY)%l_enabled = .true. #ifndef ENABLE_DRY_MODE - nc_dset(NC_DBUOY)%l_enabled = .true. - nc_dset(NC_HUM)%l_enabled = .true. - nc_dset(NC_LBUOY)%l_enabled = .true. + nc_dset(NC_QV)%l_enabled = .true. + nc_dset(NC_QL)%l_enabled = .true. #endif nc_dset(NC_VOL)%l_enabled = .true. endif @@ -528,6 +547,12 @@ subroutine set_netcdf_field_info unit='1/s', & dtype=NF90_DOUBLE) + call nc_dset(NC_THETA)%set_info(name='theta', & + long_name='potential temperature', & + std_name='', & + unit='K', & + dtype=NF90_DOUBLE) + call nc_dset(NC_TBUOY)%set_info(name='buoyancy', & long_name='total buoyancy', & std_name='', & @@ -535,23 +560,17 @@ subroutine set_netcdf_field_info dtype=NF90_DOUBLE) #ifndef ENABLE_DRY_MODE - call nc_dset(NC_DBUOY)%set_info(name='dry_buoyancy', & - long_name='dry buoyancy', & + call nc_dset(NC_QV)%set_info(name='qv', & + long_name='water vapor mixing ratio', & std_name='', & - unit='m/s^2', & + unit='kg/kg', & dtype=NF90_DOUBLE) - call nc_dset(NC_HUM)%set_info(name='humidity', & - long_name='specific humidity', & + call nc_dset(NC_QL)%set_info(name='ql', & + long_name='liquid water mixing ratio', & std_name='', & unit='kg/kg', & dtype=NF90_DOUBLE) - - call nc_dset(NC_LBUOY)%set_info(name='liquid_water_content', & - long_name='liquid-water content', & - std_name='', & - unit='1', & - dtype=NF90_DOUBLE) #endif call nc_dset(NC_VOL)%set_info(name='volume', & diff --git a/src/3d/fields/fields.f90 b/src/3d/fields/fields.f90 index 8a82c624..2e9a2e67 100644 --- a/src/3d/fields/fields.f90 +++ b/src/3d/fields/fields.f90 @@ -36,9 +36,10 @@ module fields double precision, allocatable, dimension(:, :, :) :: & #ifndef ENABLE_DRY_MODE - dbuoyg, & ! dry buoyancy (or liquid-water buoyancy) - humg, & ! humidity + qvg, & ! humidity + qlg, & ! liquid water #endif + thetag, & ! dry buoyancy (or liquid-water buoyancy) tbuoyg, & ! buoyancy #ifndef NDEBUG sym_volg, & ! symmetry volume (debug mode only) @@ -80,8 +81,8 @@ subroutine field_alloc allocate(velog(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1), n_dim)) allocate(velgradg(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1), 8)) - allocate(volg(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1))) allocate(strain_mag(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1))) + allocate(volg(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1))) #ifndef NDEBUG allocate(sym_volg(hlo(3):hhi(3), hlo(2):hhi(2), hlo(1):hhi(1))) @@ -92,10 +93,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))) @@ -116,9 +117,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 @@ -135,14 +137,15 @@ subroutine field_dealloc if (allocated(velog)) then deallocate(velog) deallocate(velgradg) - deallocate(volg) deallocate(strain_mag) + deallocate(volg) 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 bf273fd3..0718e6e5 100644 --- a/src/3d/parcels/parcel_container.f90 +++ b/src/3d/parcels/parcel_container.f90 @@ -38,9 +38,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 #ifdef ENABLE_LABELS IDX_LABEL, & ! label @@ -81,12 +82,13 @@ module parcel_container double precision, allocatable, dimension(:) :: & volume, & #ifndef ENABLE_DRY_MODE - humidity, & + qv, & + ql, & #endif #ifdef ENABLE_LABELS dilution, & #endif - buoyancy + theta ! LS-RK4 variables double precision, allocatable, dimension(:, :) :: & @@ -125,12 +127,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 #ifdef ENABLE_LABELS IDX_LABEL = i @@ -265,9 +268,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 #ifdef ENABLE_LABELS parcels%label(n) = parcels%label(m) @@ -306,9 +310,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 #ifdef ENABLE_LABELS call resize_array(parcels%label, new_size, n_parcels) @@ -339,9 +344,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 #ifdef ENABLE_LABELS allocate(parcels%label(num)) @@ -373,9 +379,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 #ifdef ENABLE_LABELS deallocate(parcels%label) @@ -402,9 +409,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 #ifdef ENABLE_LABELS buffer(IDX_LABEL) = parcels%label(n) @@ -430,9 +438,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 #ifdef ENABLE_LABELS parcels%label(n) = nint(buffer(IDX_LABEL)) diff --git a/src/3d/parcels/parcel_correction.f90 b/src/3d/parcels/parcel_correction.f90 index 8ab7afe4..b35541e4 100644 --- a/src/3d/parcels/parcel_correction.f90 +++ b/src/3d/parcels/parcel_correction.f90 @@ -128,7 +128,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 diff --git a/src/3d/parcels/parcel_damping.f90 b/src/3d/parcels/parcel_damping.f90 index 1e0b9677..97606aa9 100644 --- a/src/3d/parcels/parcel_damping.f90 +++ b/src/3d/parcels/parcel_damping.f90 @@ -60,14 +60,13 @@ subroutine parcel_damp(dt) vortg(nz+1, :, :, :) = vortg(nz-1, :, :, :) #ifndef ENABLE_DRY_MODE - humg(-1, :, :) = humg(1, :, :) - humg(nz+1, :, :) = humg(nz-1, :, :) - dbuoyg(-1, :, :) = dbuoyg(1, :, :) - dbuoyg(nz+1, :, :) = dbuoyg(nz-1, :, :) -#else - tbuoyg(-1, :, :) = tbuoyg(1, :, :) - tbuoyg(nz+1, :, :) = tbuoyg(nz-1, :, :) + qvg(-1, :, :) = qvg(1, :, :) + qvg(nz+1, :, :) = qvg(nz-1, :, :) + qlg(-1, :, :) = qlg(1, :, :) + qlg(nz+1, :, :) = qlg(nz-1, :, :) #endif + thetag(-1, :, :) = thetag(1, :, :) + thetag(nz+1, :, :) = thetag(nz-1, :, :) !$omp end parallel workshare call get_strain_magnitude_field @@ -91,11 +90,11 @@ subroutine perturbation_damping(dt, l_reuse) ! before modifying the parcel attribute double precision :: vortend(3) #ifndef ENABLE_DRY_MODE - double precision :: dbuoytend - double precision :: humtend -#else - double precision :: tbuoytend + double precision :: qvtend + double precision :: qltend #endif + double precision :: thetatend + call start_timer(damping_timer) ! This is only here to allow debug compilation @@ -107,9 +106,9 @@ subroutine perturbation_damping(dt, l_reuse) !$omp parallel default(shared) !$omp do private(n, p, l, points, pvol, weight, surface_index) & #ifndef ENABLE_DRY_MODE - !$omp& private(is, js, ks, weights, vortend, humtend, dbuoytend, time_fact) + !$omp& private(is, js, ks, weights, vortend, qvtend, qltend, thetatend, time_fact) #else - !$omp& private(is, js, ks, weights, vortend, tbuoytend, time_fact) + !$omp& private(is, js, ks, weights, vortend, thetatend, time_fact) #endif do n = 1, n_parcels ! check if only surface damping applies and we are far from surfaces @@ -131,11 +130,10 @@ subroutine perturbation_damping(dt, l_reuse) #endif vortend = zero #ifndef ENABLE_DRY_MODE - humtend = zero - dbuoytend = zero -#else - tbuoytend = zero + qvtend = zero + qltend = zero #endif + thetatend = zero ! we have 4 points per ellipsoid do p = 1, n_points_p2g @@ -154,11 +152,10 @@ subroutine perturbation_damping(dt, l_reuse) if (damping%l_scalars) then time_fact = one - exp(-damping%scalars_prefactor * strain_mag(ks:ks+1, js:js+1, is:is+1) * dt) #ifndef ENABLE_DRY_MODE - humtend = humtend + sum(weight * time_fact * (humg(ks:ks+1, js:js+1, is:is+1) - parcels%humidity(n))) - dbuoytend = dbuoytend + sum(weight * time_fact * (dbuoyg(ks:ks+1, js:js+1, is:is+1) - parcels%buoyancy(n))) -#else - tbuoytend = tbuoytend + sum(weight * time_fact * (tbuoyg(ks:ks+1, js:js+1, is:is+1) - parcels%buoyancy(n))) + qvtend = qvtend + sum(weight * time_fact * (qvg(ks:ks+1, js:js+1, is:is+1) - parcels%qv(n))) + qltend = qltend + sum(weight * time_fact * (qlg(ks:ks+1, js:js+1, is:is+1) - parcels%ql(n))) #endif + thetatend = thetatend + sum(weight * time_fact * (thetag(ks:ks+1, js:js+1, is:is+1) - parcels%theta(n))) endif if (damping%l_surface_vorticity .or. damping%l_surface_scalars) then @@ -189,14 +186,13 @@ subroutine perturbation_damping(dt, l_reuse) ! Note this exponential factor can be different for vorticity/scalars time_fact = one - exp(-damping%scalars_prefactor * strain_mag(ks:ks+1, js:js+1, is:is+1) * dt) #ifndef ENABLE_DRY_MODE - humtend = humtend + sum(weight(surface_index, :, :) * time_fact(surface_index, :, :) * & - (humg(ks+surface_index, js:js+1, is:is+1) - parcels%humidity(n))) - dbuoytend = dbuoytend + sum(weight(surface_index, :, :) * time_fact(surface_index, :, :) * & - (dbuoyg(ks+surface_index, js:js+1, is:is+1) - parcels%buoyancy(n))) -#else - tbuoytend = tbuoytend + sum(weight(surface_index, :, :) * time_fact(surface_index, :, :) * & - (tbuoyg(ks+surface_index, js:js+1, is:is+1) - parcels%buoyancy(n))) + qvtend = qvtend + sum(weight(surface_index, :, :) * time_fact(surface_index, :, :) * & + (qvg(ks+surface_index, js:js+1, is:is+1) - parcels%qv(n))) + qltend = qltend + sum(weight(surface_index, :, :) * time_fact(surface_index, :, :) * & + (qlg(ks+surface_index, js:js+1, is:is+1) - parcels%ql(n))) #endif + thetatend = thetatend + sum(weight(surface_index, :, :) * time_fact(surface_index, :, :) * & + (thetag(ks+surface_index, js:js+1, is:is+1) - parcels%theta(n))) endif enddo ! Add all the tendencies only at the end @@ -207,11 +203,10 @@ subroutine perturbation_damping(dt, l_reuse) endif if (damping%l_scalars .or. damping%l_surface_scalars) then #ifndef ENABLE_DRY_MODE - parcels%humidity(n) = parcels%humidity(n) + humtend - parcels%buoyancy(n) = parcels%buoyancy(n) + dbuoytend -#else - parcels%buoyancy(n) = parcels%buoyancy(n) + tbuoytend + parcels%qv(n) = parcels%qv(n) + qvtend + parcels%ql(n) = parcels%ql(n) + qltend #endif + parcels%theta(n) = parcels%theta(n) + thetatend endif enddo !$omp end do diff --git a/src/3d/parcels/parcel_diagnostics.f90 b/src/3d/parcels/parcel_diagnostics.f90 index 6334a07c..c9811b36 100644 --- a/src/3d/parcels/parcel_diagnostics.f90 +++ b/src/3d/parcels/parcel_diagnostics.f90 @@ -10,7 +10,7 @@ module parcel_diagnostics use parcel_split_mod, only : n_parcel_splits use parcel_merging, only : n_parcel_merges, n_big_close, n_way_parcel_mergers 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_environment @@ -73,7 +73,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_init.f90 b/src/3d/parcels/parcel_init.f90 index b1faf75d..94bef68a 100644 --- a/src/3d/parcels/parcel_init.f90 +++ b/src/3d/parcels/parcel_init.f90 @@ -98,9 +98,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 @@ -195,11 +196,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 @@ -213,7 +216,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 @@ -222,9 +225,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 @@ -232,9 +236,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_interpl.f90 b/src/3d/parcels/parcel_interpl.f90 index a1508d87..9b6dc3c7 100644 --- a/src/3d/parcels/parcel_interpl.f90 +++ b/src/3d/parcels/parcel_interpl.f90 @@ -24,7 +24,7 @@ module parcel_interpl , field_buffer_to_interior_integer & , field_interior_to_buffer_integer & , field_buffer_to_halo_integer - use physics, only : glat, lambda_c, q_0 + use physics, only : gravity, theta_0, qv_dens_coeff, r_d, c_p, L_v #ifdef ENABLE_BUOYANCY_PERTURBATION_MODE use physics, only : bfsq #endif @@ -52,13 +52,17 @@ module parcel_interpl , IDX_TBUOY_SWAP = 5 & , IDX_NPARG_SWAP = 6 & , IDX_NSPARG_SWAP = 7 -#ifndef ENABLE_DRY_MODE - integer, parameter :: IDX_DBUOY_SWAP = 8 & - , IDX_HUM_SWAP = 9 - integer, parameter :: n_field_swap = 9 -#else integer, parameter :: n_field_swap = 7 + + ! restart indices for par2grid_diag + integer, parameter :: IDX_THETA_SWAP = 2 +#ifndef ENABLE_DRY_MODE + integer, parameter :: IDX_QV_SWAP = 3 & + , IDX_QL_SWAP = 4 + integer, parameter :: n_field_swap_diag = 4 +#else + integer, parameter :: n_field_swap_diag = 2 #endif #ifndef ENABLE_G2P_1POINT @@ -78,12 +82,14 @@ module parcel_interpl #endif public :: par2grid & + , par2grid_diag & , vol2grid & , grid2par & , par2grid_timer & , grid2par_timer & , halo_swap_timer & , trilinear & + , saturation_adjustment & , bilinear & , n_points_p2g & , point_weight_p2g @@ -153,44 +159,30 @@ subroutine par2grid(l_reuse) double precision :: points(3, n_points_p2g) 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 * exp(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_dens_coeff*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 #ifdef ENABLE_BUOYANCY_PERTURBATION_MODE @@ -223,12 +215,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) & @@ -280,19 +266,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, :, :) @@ -312,10 +285,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 @@ -328,6 +297,129 @@ 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(l_reuse) + logical, optional :: l_reuse + double precision :: points(3, 4) + integer :: n, p, i, j, k + double precision :: pvol, weight(0:1,0:1,0:1) + + thetag = zero + volg = zero +#ifndef ENABLE_DRY_MODE + qvg = zero + qlg = zero +#endif + !$omp parallel default(shared) +#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, 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, volg) +#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) + 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) + 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 + + !$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, :, :) + 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, :, :) + 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, :, :) +#endif + !$omp end parallel workshare + + + 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 @@ -347,10 +439,6 @@ subroutine par2grid_halo_swap call field_halo_to_buffer(tbuoyg, IDX_TBUOY_SWAP) call field_halo_to_buffer_integer(nparg, IDX_NPARG_SWAP) call field_halo_to_buffer_integer(nsparg, IDX_NSPARG_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 @@ -364,10 +452,6 @@ subroutine par2grid_halo_swap call field_buffer_to_interior(tbuoyg, IDX_TBUOY_SWAP, .true.) call field_buffer_to_interior_integer(nparg, IDX_NPARG_SWAP, .true.) call field_buffer_to_interior_integer(nsparg, IDX_NSPARG_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: @@ -379,10 +463,6 @@ subroutine par2grid_halo_swap call field_interior_to_buffer(tbuoyg, IDX_TBUOY_SWAP) call field_interior_to_buffer_integer(nparg, IDX_NPARG_SWAP) call field_interior_to_buffer_integer(nsparg, IDX_NSPARG_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 @@ -393,15 +473,60 @@ subroutine par2grid_halo_swap call field_buffer_to_halo(tbuoyg, IDX_TBUOY_SWAP, .false.) call field_buffer_to_halo_integer(nparg, IDX_NPARG_SWAP, .false.) call field_buffer_to_halo_integer(nsparg, IDX_NSPARG_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, ndim=3) + + !------------------------------------------------------------------ + ! 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) + 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(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.) + call field_buffer_to_interior(qlg, IDX_QL_SWAP, .true.) +#endif + + !------------------------------------------------------------------ + ! 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) + call field_interior_to_buffer(qlg, IDX_QL_SWAP) +#endif + + 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.) + 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. @@ -524,7 +649,65 @@ pure subroutine trilinear(pos, ii, jj, kk, ww) end subroutine trilinear - !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: + + 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 = 8619.0 ! Scale height for bomex + double precision, parameter :: surf_press = 100000.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 + double precision :: err_at_temp, err_at_temp_inv_deriv,efact,divfact + integer :: n, iter + + ! TO DO: ADD OPENMP + do n = 1, n_parcels + 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 + 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/(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 + 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 + ! Use Newton-Raphson to converge + else + do iter=1,3 + efact=0.01*press*exp(qsa2*(temp - tk0c)/(temp - qsa3)) + qsat=qsa1/(efact - qsa4) + ql_iter=max(qt_start-qsat,0.0) + err_at_temp=temp-(temp_start-(L_v/c_p)*(ql_start-ql_iter)) + if(ql_iter>0.0) then + !calculate 1/(d err/ dt) to save a division latet on + divfact=((efact - qsa4)*(efact - qsa4)*(temp - qsa3)*(temp - qsa3)) + err_at_temp_inv_deriv=divfact/(divfact+(L_v/c_p)*(qsa1*qsa2*efact*(qsa3-tk0c))) + else + err_at_temp_inv_deriv=1.0 + endif + temp=temp-err_at_temp*err_at_temp_inv_deriv + enddo + 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 + parcels%ql(n)=ql_iter + end if + end do + + end subroutine saturation_adjustment + + !:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ! Bi-linear interpolation ! @param[in] pos position of the parcel diff --git a/src/3d/parcels/parcel_ls_forcings.f90 b/src/3d/parcels/parcel_ls_forcings.f90 new file mode 100644 index 00000000..8fdcf5eb --- /dev/null +++ b/src/3d/parcels/parcel_ls_forcings.f90 @@ -0,0 +1,248 @@ +! ============================================================================= +! This module contains the subroutines to do parcel-to-grid and grid-to-parcel +! interpolation. +! ============================================================================= +module parcel_ls_forcings + use constants, only : zero, one, pi + use mpi_timer, only : start_timer, stop_timer + use options, only : parcel, forcings + use parameters, only : nx,ny,nz, dxi, dx + use parcel_container, only : parcels, n_parcels + use parcel_ellipsoid + use fields + use omp_lib + use inversion_utils, only : field_decompose_physical + use mpi_utils, only : mpi_exit_on_error + use mpi_layout + use mpi_environment + use mpi_collectives, only : mpi_blocking_reduce + use parcel_interpl, only : saturation_adjustment, par2grid, par2grid_diag, trilinear + + implicit none + + private + + public :: apply_ls_forcings + + contains + + subroutine apply_ls_forcings(dt) + double precision, intent(in) :: dt + if(forcings%l_ls_forcings) then + call apply_subsidence_and_vorticity_adjustment(dt) + call apply_ls_tendencies(dt) + call saturation_adjustment + end if + + end subroutine apply_ls_forcings + ! + ! @pre + subroutine apply_subsidence_and_vorticity_adjustment(dt) + double precision, intent(in) :: dt + double precision :: weights(0:1,0:1,0:1) + double precision :: xf, yf, zf, xfc, yfc, zfc + double precision :: thetagradz, qvgradz, qlgradz, ww + integer :: n, kk, is, js, ks + double precision :: xvortbar(0:nz), yvortbar(0:nz), nudgefac + + ! use par2grid_diag theta + ! use par2grid_diag ql+qv + ! to be replaced by special par2grid? + call par2grid(.false.) + call par2grid_diag(.false.) + + ! code to apply subsidence, use local gradients + + !$omp parallel default(shared) + !$omp do private(n, is, js, ks, weights, xf, yf, zf, xfc, yfc, zfc, thetagradz, qvgradz, qlgradz, ww) + do n = 1, n_parcels + + call trilinear(parcels%position(:, n), is, js, ks, weights) + + xf = weights(0,0,1) + weights(0,1,1) + weights(1,0,1) + weights(1,1,1) ! fractional position along x + yf = weights(0,1,0) + weights(0,1,1) + weights(1,1,0) + weights(1,1,1) ! fractional position along y + zf = weights(1,0,0) + weights(1,0,1) + weights(1,1,0) + weights(1,1,1) ! fractional position along z + xfc=one-xf + yfc=one-yf + zfc=one-zf + + thetagradz = dxi(3)*(& + xfc * (& + yfc * (thetag(ks+1, js, is ) - thetag(ks, js, is)) & + + yf * (thetag(ks+1, js+1, is) - thetag(ks, js+1, is))) & + + xf * (& + yfc * (thetag(ks+1, js, is+1) - thetag(ks, js, is+1)) & + + yf * (thetag(ks+1, js+1, is+1) - thetag(ks, js+1, is+1)))) + + qvgradz = dxi(3)*(& + xfc * (& + yfc * (qvg(ks+1, js, is ) - qvg(ks, js, is)) & + + yf * (qvg(ks+1, js+1, is) - qvg(ks, js+1, is))) & + + xf * (& + yfc * (qvg(ks+1, js, is+1) - qvg(ks, js, is+1)) & + + yf * (qvg(ks+1, js+1, is+1) - qvg(ks, js+1, is+1)))) + + qlgradz = dxi(3)*(& + xfc * (& + yfc * (qlg(ks+1, js, is ) - qlg(ks, js, is)) & + + yf * (qlg(ks+1, js+1, is) - qlg(ks, js+1, is))) & + + xf * (& + yfc * (qlg(ks+1, js, is+1) - qlg(ks, js, is+1)) & + + yf * (qlg(ks+1, js+1, is+1) - qlg(ks, js+1, is+1)))) + + ww=get_bomex_subsidence_velocity(parcels%position(3, n)) + + parcels%theta(n) = parcels%theta(n)-dt*ww*thetagradz + parcels%qv(n) = parcels%qv(n)-dt*ww*(qvgradz+qlgradz) + enddo + !$omp end do + !$omp end parallel + + do kk=0,nz + xvortbar(kk)= sum(vortg(kk,box%lo(2):box%hi(2),box%lo(1):box%hi(1),1)) + yvortbar(kk)= sum(vortg(kk,box%lo(2):box%hi(2),box%lo(1):box%hi(1),2)) + end do + + call MPI_Allreduce(MPI_IN_PLACE, & + xvortbar, & + nz, & + MPI_DOUBLE_PRECISION, & + MPI_SUM, & + world%comm, & + world%err) + call MPI_Allreduce(MPI_IN_PLACE, & + yvortbar, & + nz, & + MPI_DOUBLE_PRECISION, & + MPI_SUM, & + world%comm, & + world%err) + + do kk=0,nz + xvortbar(kk)= xvortbar(kk)/(nx*ny)-get_bomex_xvort(lower(3)+kk*dx(3)) + yvortbar(kk)= yvortbar(kk)/(nx*ny)-get_bomex_yvort(lower(3)+kk*dx(3)) + enddo + + !if(cart%rank == cart%root) then + ! write(*,*) 'xvortbar' + ! write(*,*) xvortbar + ! write(*,*) 'yvortbar' + ! write(*,*) yvortbar + !end if + + ! nudge velocity profiles with 120 seconds time scale + nudgefac=(1.0-exp(-dt/120.0)) + + !$omp parallel default(shared) + !$omp do private(n, is, js, ks, weights, zf, zfc) + do n = 1, n_parcels + + call trilinear(parcels%position(:, n), is, js, ks, weights) + + zf = weights(1,0,0) + weights(1,0,1) + weights(1,1,0) + weights(1,1,1) ! fractional position along z + zfc= one-zf + + parcels%vorticity(1, n) = parcels%vorticity(1, n) - zfc*xvortbar(ks)*nudgefac & + - zf*xvortbar(ks+1)*nudgefac + parcels%vorticity(2, n) = parcels%vorticity(2, n) - zfc*yvortbar(ks)*nudgefac & + - zf*yvortbar(ks+1)*nudgefac + enddo + !$omp end do + !$omp end parallel + + end subroutine apply_subsidence_and_vorticity_adjustment + + ! + ! @pre + subroutine apply_ls_tendencies(dt) + double precision, intent(in) :: dt + integer :: n + + !$omp parallel default(shared) + !$omp do private(n) + do n = 1, n_parcels + parcels%theta(n) = parcels%theta(n) + get_bomex_theta_ls(parcels%position(3, n))*dt + parcels%qv(n) = parcels%qv(n) + get_bomex_qt_ls(parcels%position(3, n))*dt + enddo + !$omp end do + !$omp end parallel + + end subroutine apply_ls_tendencies + ! + ! @pre + elemental function get_bomex_subsidence_velocity(zz) result (ww) + double precision, intent(in) :: zz + double precision :: ww + + if (zz < 1500.0) then + ww = zz*(-0.0065)/1500.0 + else if (zz < 2100.0) then + ww = -0.0065 + (zz-1500.0)*(0.0065)/(2100.0-1500.0) + else + ww = 0.0 + endif + end function get_bomex_subsidence_velocity + + ! + ! @pre + elemental function get_bomex_theta_ls(zz) result (theta_ls) + double precision, intent(in) :: zz + double precision :: theta_ls + double precision, parameter :: iday=1./86400. + + if (zz < 1500.0) then + theta_ls = -2.0*iday + else + theta_ls = (-2.0 + (zz-1500.0)*(2.0/(3000.0-1500.0)))*iday + endif + end function get_bomex_theta_ls + + ! + ! @pre + elemental function get_bomex_qt_ls(zz) result (qt_ls) + double precision, intent(in) :: zz + double precision :: qt_ls + + if(zz <= 300.0) then + qt_ls = -1.2e-8 + else if(zz <= 500.0) then + qt_ls = -1.2e-8 + (zz-300.)*(1.2e-8)/(500.0-300.0) + else + qt_ls = 0.0 + endif + end function get_bomex_qt_ls + + ! + ! @pre + elemental function get_bomex_xvort(zz) result (xvort) + double precision, intent(in) :: zz + double precision :: xvort + + if (zz<600.0) then + xvort=0.0043*cos(pi*(zz+600)/1200.0)+0.0023 + else if (zz<1400.0) then + xvort=-0.001-0.001*cos(pi*(zz-600)/800.0) + else + xvort=0.0 + end if + + end function get_bomex_xvort + + ! @pre + elemental function get_bomex_yvort(zz) result (yvort) + double precision, intent(in) :: zz + double precision :: yvort + + if(zz<300.0) then + yvort=0.0062*cos((pi*(zz-300.0)/600.0)**2)-0.0067 + else if (zz<600.0) then + yvort=0.0027*cos(pi*(zz-300.0)/600.0)-0.0032 + else if (zz<1400.0) then + yvort=0.005*cos(pi*(zz-1400.0)/1600.0)-0.0032 + else + yvort=0.0018 + end if + + end function get_bomex_yvort + +end module parcel_ls_forcings diff --git a/src/3d/parcels/parcel_merge.f90 b/src/3d/parcels/parcel_merge.f90 index 9ec55394..6ae9a9ab 100644 --- a/src/3d/parcels/parcel_merge.f90 +++ b/src/3d/parcels/parcel_merge.f90 @@ -128,9 +128,10 @@ subroutine do_group_merge(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 #ifdef ENABLE_LABELS double precision :: labelm(n_merge) @@ -175,9 +176,10 @@ subroutine do_group_merge(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) @@ -220,9 +222,10 @@ subroutine do_group_merge(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 @@ -246,9 +249,10 @@ subroutine do_group_merge(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 @@ -286,9 +290,10 @@ subroutine do_group_merge(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 #ifdef ENABLE_LABELS parcels%label(ic) = labelm(l) diff --git a/src/3d/parcels/parcel_netcdf.f90 b/src/3d/parcels/parcel_netcdf.f90 index a4609fbe..477771ad 100644 --- a/src/3d/parcels/parcel_netcdf.f90 +++ b/src/3d/parcels/parcel_netcdf.f90 @@ -47,7 +47,7 @@ module parcel_netcdf , NC_X_POS = 7 & , NC_Y_POS = 8 & , NC_Z_POS = 9 & - , NC_BUOY = 10 & + , NC_THETA = 10 & , NC_X_VOR = 11 & , NC_Y_VOR = 12 & , NC_Z_VOR = 13 & @@ -61,15 +61,16 @@ module parcel_netcdf logical :: l_unable = .false. #ifndef ENABLE_DRY_MODE - integer, parameter :: NC_HUM = 19 + integer, parameter :: NC_QV = 19 & + , NC_QL = 20 #ifdef ENABLE_LABELS - integer, parameter :: NC_LABEL = 20 & - , NC_DILUTION = 21 + integer, parameter :: NC_LABEL = 21 & + , NC_DILUTION = 22 type(netcdf_info) :: nc_dset(NC_DILUTION) #else - type(netcdf_info) :: nc_dset(NC_HUM) + type(netcdf_info) :: nc_dset(NC_QL) #endif #else @@ -302,10 +303,11 @@ subroutine write_netcdf_parcels(t) call write_parcel_attribute(NC_Y_VOR, parcels%vorticity(2, :), start, cnt) call write_parcel_attribute(NC_Z_VOR, parcels%vorticity(3, :), start, cnt) - call write_parcel_attribute(NC_BUOY, parcels%buoyancy, start, cnt) + call write_parcel_attribute(NC_THETA, parcels%theta, start, cnt) #ifndef ENABLE_DRY_MODE - call write_parcel_attribute(NC_HUM, parcels%humidity, start, cnt) + call write_parcel_attribute(NC_QV, parcels%qv, start, cnt) + call write_parcel_attribute(NC_QL, parcels%ql, start, cnt) #endif #ifdef ENABLE_LABELS call write_parcel_attribute_int(NC_LABEL, parcels%label, start, cnt) @@ -656,17 +658,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, 'humidity', & - parcels%humidity(pfirst:plast), start, cnt) + 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, 'ql', & + parcels%ql(pfirst:plast), start, cnt) endif #endif #ifdef ENABLE_LABELS @@ -681,7 +688,7 @@ subroutine read_chunk(first, last, pfirst) if (.not. l_valid) then call mpi_exit_on_error(& - "Either the parcel buoyancy or vorticity must be present! Exiting.") + "Either the parcel theta or vorticity must be present! Exiting.") endif end subroutine read_chunk @@ -703,9 +710,10 @@ subroutine set_netcdf_parcel_output nc_dset(NC_X_VOR)%l_enabled = .true. nc_dset(NC_Y_VOR)%l_enabled = .true. nc_dset(NC_Z_VOR)%l_enabled = .true. - nc_dset(NC_BUOY)%l_enabled = .true. + nc_dset(NC_THETA)%l_enabled = .true. #ifndef ENABLE_DRY_MODE - nc_dset(NC_HUM)%l_enabled = .true. + nc_dset(NC_QV)%l_enabled = .true. + nc_dset(NC_QL)%l_enabled = .true. #endif #ifdef ENABLE_LABELS nc_dset(NC_LABEL)%l_enabled = .true. @@ -754,9 +762,10 @@ subroutine set_netcdf_parcel_output nc_dset(NC_X_VOR)%l_enabled = .true. nc_dset(NC_Y_VOR)%l_enabled = .true. nc_dset(NC_Z_VOR)%l_enabled = .true. - nc_dset(NC_BUOY)%l_enabled = .true. + nc_dset(NC_THETA)%l_enabled = .true. #ifndef ENABLE_DRY_MODE - nc_dset(NC_HUM)%l_enabled = .true. + nc_dset(NC_QV)%l_enabled = .true. + nc_dset(NC_QL)%l_enabled = .true. #endif #ifdef ENABLE_LABELS nc_dset(NC_LABEL)%l_enabled = .true. @@ -799,13 +808,13 @@ subroutine set_netcdf_parcel_output ((nc_dset(NC_X_VOR)%l_enabled .and. & nc_dset(NC_Y_VOR)%l_enabled .and. & nc_dset(NC_Z_VOR)%l_enabled) .or. & - nc_dset(NC_BUOY)%l_enabled)) + nc_dset(NC_THETA)%l_enabled)) if ((.not. l_enabled_restart) .and. (world%rank == world%root)) then print *, "WARNING: EPIC will not be able to restart from a parcel file." print *, " You must at least write the B-shape matrix, parcel position" - print *, " parcel volume and parcel vorticity or buoyancy to enable a" + print *, " parcel volume and parcel vorticity or theta to enable a" print *, " restart. If you intend to restart from a parcel file later," print *, " you must stop the simulation immediately. Furthermore, you can" print *, " write the MPI 'start_index' to speed up the restart." @@ -920,17 +929,22 @@ subroutine set_netcdf_parcel_info unit='1/s', & dtype=NF90_DOUBLE) - call nc_dset(NC_BUOY)%set_info(name='buoyancy', & - long_name='parcel buoyancy', & - std_name='', & - unit='m/s^2', & + call nc_dset(NC_THETA)%set_info(name='theta' , & + long_name='parcel potential temperature', & + std_name='', & + unit='K', & dtype=NF90_DOUBLE) #ifndef ENABLE_DRY_MODE - call nc_dset(NC_HUM)%set_info(name='humidity', & - long_name='parcel humidity', & - std_name='', & - unit='1', & + call nc_dset(NC_QV)%set_info(name='qv', & + long_name='parcel water vapour mixing ratio', & + std_name='', & + unit='1', & + dtype=NF90_DOUBLE) + call nc_dset(NC_QL)%set_info(name='ql', & + long_name='parcel liquid water mixing ratio', & + std_name='', & + unit='1', & dtype=NF90_DOUBLE) #endif #ifdef ENABLE_LABELS diff --git a/src/3d/parcels/parcel_split.f90 b/src/3d/parcels/parcel_split.f90 index 22060f9c..e70625f0 100644 --- a/src/3d/parcels/parcel_split.f90 +++ b/src/3d/parcels/parcel_split.f90 @@ -144,9 +144,10 @@ subroutine parcel_split 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 #ifdef ENABLE_LABELS parcels%label(n_thread_loc) = parcels%label(n) diff --git a/src/3d/stepper/ls_rk.f90 b/src/3d/stepper/ls_rk.f90 index f0767eb9..cb266923 100644 --- a/src/3d/stepper/ls_rk.f90 +++ b/src/3d/stepper/ls_rk.f90 @@ -8,8 +8,11 @@ module ls_rk use parcel_container use parcel_bc use parcel_mpi, only : parcel_communicate - use rk_utils, only: get_dBdt, get_time_step + use rk_utils, only: get_dBdt, get_time_step, get_strain_magnitude_field use utils, only : write_step + use parcel_interpl, only : par2grid, grid2par, saturation_adjustment + use parcel_damping, only : parcel_damp + use parcel_ls_forcings, only : apply_ls_forcings use parcel_interpl, only : par2grid, grid2par use parcel_damping, only : parcel_damp use fields, only : velgradg, velog, vortg, vtend, tbuoyg @@ -127,6 +130,7 @@ subroutine ls_rk_step(t) call stop_timer(rk_timer) call parcel_damp(dt) + call apply_ls_forcings(dt) ! we need to subtract 14 calls since we start and stop ! the timer multiple times which increments n_calls @@ -195,8 +199,10 @@ subroutine ls_rk_substep(dt, step) enddo !$omp end parallel do - call stop_timer(rk_timer) + ! call saturation adjustment after each integration + call saturation_adjustment + call stop_timer(rk_timer) if (step == n_stages) then call parcel_communicate return diff --git a/src/3d/utils/options.f90 b/src/3d/utils/options.f90 index af37a4de..61825861 100644 --- a/src/3d/utils/options.f90 +++ b/src/3d/utils/options.f90 @@ -77,6 +77,17 @@ module options ! logical :: allow_larger_anisotropy = .false. + + ! + ! + ! + type forcings_type + ! use large-scale forcing + logical :: l_ls_forcings = .false. + end type forcings_type + + type(forcings_type) :: forcings + ! ! parcel options ! @@ -131,7 +142,7 @@ subroutine read_config_file logical :: exists = .false. ! namelist definitions - namelist /EPIC/ field_file, flux_file, rk_order, boundary, output, parcel, time, damping + namelist /EPIC/ field_file, flux_file, rk_order, boundary, output, parcel, forcings, time, damping ! check whether file exists inquire(file=filename, exist=exists) @@ -188,6 +199,8 @@ subroutine write_netcdf_options(ncid) call write_netcdf_attribute(ncid, "correction_iters", parcel%correction_iters) call write_netcdf_attribute(ncid, "gradient_pref", parcel%gradient_pref) call write_netcdf_attribute(ncid, "max_compression", parcel%max_compression) + + call write_netcdf_attribute(ncid, "l_ls_forcings", forcings%l_ls_forcings) call write_netcdf_attribute(ncid, "parcel_freq", output%parcel_freq) call write_netcdf_attribute(ncid, "field_freq", output%field_freq) diff --git a/src/utils/physics.f90 b/src/utils/physics.f90 index 643091bb..9cff4338 100644 --- a/src/utils/physics.f90 +++ b/src/utils/physics.f90 @@ -70,7 +70,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: ! @@ -171,7 +176,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 diff --git a/unit-tests/3d/test_ellipsoid_split.f90 b/unit-tests/3d/test_ellipsoid_split.f90 index 4cda52b3..e6978b7d 100644 --- a/unit-tests/3d/test_ellipsoid_split.f90 +++ b/unit-tests/3d/test_ellipsoid_split.f90 @@ -76,9 +76,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 @@ -133,9 +134,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 @@ -148,9 +150,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 9827bf8b..8d4a5cff 100644 --- a/unit-tests/3d/test_mpi_nearest_1.f90 +++ b/unit-tests/3d/test_mpi_nearest_1.f90 @@ -110,8 +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 + world%rank * 100 - + parcels%theta(l) = l + world%rank * 100 l = l + 1 do m = -1, 1, 2 @@ -119,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 enddo @@ -128,7 +127,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 + world%rank * 100 + parcels%theta(l) = l + world%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 12310b5a..398b1736 100644 --- a/unit-tests/3d/test_mpi_nearest_10.f90 +++ b/unit-tests/3d/test_mpi_nearest_10.f90 @@ -111,7 +111,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel c @@ -119,7 +119,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -127,7 +127,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 + world%rank * 100 + parcels%theta(l) = l + world%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 94b6aadd..6e16074c 100644 --- a/unit-tests/3d/test_mpi_nearest_11.f90 +++ b/unit-tests/3d/test_mpi_nearest_11.f90 @@ -111,7 +111,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel c @@ -119,7 +119,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -127,7 +127,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 + world%rank * 100 + parcels%theta(l) = l + world%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 373a2f67..2b961cb9 100644 --- a/unit-tests/3d/test_mpi_nearest_12.f90 +++ b/unit-tests/3d/test_mpi_nearest_12.f90 @@ -154,7 +154,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -162,7 +162,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -170,7 +170,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_1 @@ -195,7 +195,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -203,7 +203,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -211,7 +211,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 + world%rank * 100 + parcels%theta(l) = l + world%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.44d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -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) = 1.1d0 * vmin - parcels%buoyancy(l) = l + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -252,7 +252,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_3 @@ -277,7 +277,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -285,7 +285,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -293,7 +293,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_4 @@ -322,7 +322,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -330,7 +330,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -338,7 +338,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! @@ -350,7 +350,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -358,7 +358,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -366,7 +366,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_5 @@ -395,7 +395,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -403,7 +403,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -411,7 +411,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! @@ -423,7 +423,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -431,7 +431,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -439,7 +439,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 + world%rank * 100 + parcels%theta(l) = l + world%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 12597a01..29e4307b 100644 --- a/unit-tests/3d/test_mpi_nearest_13.f90 +++ b/unit-tests/3d/test_mpi_nearest_13.f90 @@ -180,7 +180,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -188,7 +188,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -196,7 +196,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_1 @@ -221,7 +221,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -229,7 +229,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -237,7 +237,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_2 @@ -262,7 +262,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -270,7 +270,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -278,7 +278,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_3 @@ -303,7 +303,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -311,7 +311,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -319,7 +319,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_4 @@ -348,7 +348,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 + world%rank * 100 + parcels%theta(l) = l + world%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 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + world%rank * 100 + parcels%theta(l) = l + world%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 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! @@ -376,7 +376,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -384,7 +384,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -392,7 +392,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_5 @@ -421,7 +421,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 + world%rank * 100 + parcels%theta(l) = l + world%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 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + world%rank * 100 + parcels%theta(l) = l + world%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 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! @@ -449,7 +449,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -457,7 +457,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -465,7 +465,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_6 @@ -494,7 +494,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -502,7 +502,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -510,7 +510,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! @@ -522,7 +522,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -530,7 +530,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -538,7 +538,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_7 @@ -567,7 +567,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -575,7 +575,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -583,7 +583,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! @@ -595,7 +595,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel C @@ -603,7 +603,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel a @@ -611,7 +611,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 + world%rank * 100 + parcels%theta(l) = l + world%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 b757e207..c0bf6b4e 100644 --- a/unit-tests/3d/test_mpi_nearest_14.f90 +++ b/unit-tests/3d/test_mpi_nearest_14.f90 @@ -154,7 +154,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel B @@ -162,7 +162,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel c @@ -170,7 +170,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_1 @@ -195,7 +195,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel B @@ -203,7 +203,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel c @@ -211,7 +211,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 + world%rank * 100 + parcels%theta(l) = l + world%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 + world%rank * 100 + parcels%theta(l) = l + world%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.44d0 parcels%position(3, l) = z parcels%volume(l) = 1.1d0 * vmin - parcels%buoyancy(l) = l + world%rank * 100 + parcels%theta(l) = l + world%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.34d0 parcels%position(3, l) = z parcels%volume(l) = 0.9d0 * vmin - parcels%buoyancy(l) = l + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_3 @@ -277,7 +277,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel B @@ -285,7 +285,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel c @@ -293,7 +293,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_4 @@ -322,7 +322,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel B @@ -330,7 +330,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel c @@ -338,7 +338,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! @@ -350,7 +350,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel B @@ -358,7 +358,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel c @@ -366,7 +366,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 end subroutine cell_placement_5 @@ -395,7 +395,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel B @@ -403,7 +403,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel c @@ -411,7 +411,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! @@ -423,7 +423,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! big parcel B @@ -431,7 +431,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 ! small parcel c @@ -439,7 +439,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 + world%rank * 100 + parcels%theta(l) = l + world%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 df68326a..51e773b0 100644 --- a/unit-tests/3d/test_mpi_nearest_15.f90 +++ b/unit-tests/3d/test_mpi_nearest_15.f90 @@ -141,7 +141,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -149,7 +150,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -157,7 +159,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -165,7 +168,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_1 @@ -190,7 +194,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -198,7 +203,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -206,7 +212,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -214,7 +221,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_2 @@ -239,7 +247,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -247,7 +256,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -255,7 +265,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -263,7 +274,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_3 @@ -292,7 +304,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -300,7 +313,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -308,7 +322,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -316,7 +331,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 @@ -329,7 +345,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -337,7 +354,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -345,7 +363,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -353,7 +372,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_4 @@ -378,7 +398,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -386,7 +407,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -394,7 +416,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -402,7 +425,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 @@ -415,7 +439,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -423,7 +448,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -431,7 +457,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -439,7 +466,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%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 bd776b67..311efe1f 100644 --- a/unit-tests/3d/test_mpi_nearest_16.f90 +++ b/unit-tests/3d/test_mpi_nearest_16.f90 @@ -141,7 +141,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -149,7 +150,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -157,7 +159,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -165,7 +168,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_1 @@ -190,7 +194,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -198,7 +203,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -206,7 +212,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -214,7 +221,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_2 @@ -238,7 +246,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -246,7 +255,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -254,7 +264,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -262,7 +273,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_3 @@ -291,7 +303,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -299,7 +312,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -307,7 +321,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -315,7 +330,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 @@ -328,7 +344,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -336,7 +353,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -344,7 +362,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -352,7 +371,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_4 @@ -381,7 +401,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -389,7 +410,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -397,7 +419,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -405,7 +428,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 @@ -418,7 +442,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -426,7 +451,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -434,7 +460,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -442,7 +469,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%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 1994f053..f2451178 100644 --- a/unit-tests/3d/test_mpi_nearest_17.f90 +++ b/unit-tests/3d/test_mpi_nearest_17.f90 @@ -141,7 +141,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel B @@ -149,7 +150,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -157,7 +159,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -165,7 +168,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_1 @@ -190,7 +194,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel B @@ -198,7 +203,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -206,7 +212,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -214,7 +221,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_2 @@ -238,7 +246,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel B @@ -246,7 +255,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -254,7 +264,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -262,7 +273,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_3 @@ -291,7 +303,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel B @@ -299,7 +312,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -307,7 +321,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -315,7 +330,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 @@ -328,7 +344,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel B @@ -336,7 +353,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -344,7 +362,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -352,7 +371,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 @@ -382,7 +402,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel B @@ -390,7 +411,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -398,7 +420,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -406,7 +429,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 @@ -419,7 +443,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel B @@ -427,7 +452,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -435,7 +461,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -443,7 +470,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%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 2628f82b..8cde4332 100644 --- a/unit-tests/3d/test_mpi_nearest_18.f90 +++ b/unit-tests/3d/test_mpi_nearest_18.f90 @@ -128,7 +128,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -136,7 +137,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -144,7 +146,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -152,7 +155,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_1 @@ -177,7 +181,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -185,7 +190,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -193,7 +199,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -201,7 +208,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_2 @@ -230,7 +238,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -238,7 +247,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -246,7 +256,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -254,7 +265,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 @@ -267,7 +279,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -275,7 +288,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -283,7 +297,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -291,7 +306,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_3 @@ -320,7 +336,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -328,7 +345,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -336,7 +354,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -344,7 +363,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 @@ -357,7 +377,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -365,7 +386,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -373,7 +395,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -381,7 +404,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%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 60adf5ef..6177ddf5 100644 --- a/unit-tests/3d/test_mpi_nearest_19.f90 +++ b/unit-tests/3d/test_mpi_nearest_19.f90 @@ -141,7 +141,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -149,7 +150,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -157,7 +159,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel D @@ -165,7 +168,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_1 @@ -190,7 +194,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -198,7 +203,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -206,7 +212,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel D @@ -214,7 +221,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_2 @@ -239,7 +247,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -247,7 +256,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -255,7 +265,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel D @@ -263,7 +274,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_3 @@ -292,7 +304,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -300,7 +313,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -308,7 +322,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel D @@ -316,7 +331,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 @@ -329,7 +345,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -337,7 +354,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -345,7 +363,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel D @@ -353,7 +372,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_4 @@ -382,7 +402,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -390,7 +411,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -398,7 +420,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel D @@ -406,7 +429,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 @@ -419,7 +443,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -427,7 +452,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -435,7 +461,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel D @@ -443,7 +470,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%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 a73f7c52..3a7e5749 100644 --- a/unit-tests/3d/test_mpi_nearest_2.f90 +++ b/unit-tests/3d/test_mpi_nearest_2.f90 @@ -110,7 +110,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 @@ -119,7 +120,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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 l = l + 1 enddo @@ -128,7 +129,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 + world%rank * 100 + parcels%theta(l) = l + world%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 d69ce794..1c8b143a 100644 --- a/unit-tests/3d/test_mpi_nearest_20.f90 +++ b/unit-tests/3d/test_mpi_nearest_20.f90 @@ -110,8 +110,8 @@ program test_mpi_nearest_20 ! do n = 1, n_merge ! is = isma(n) ! ic = iclo(n) -! print *, world%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 *, world%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(world%comm, world%err) @@ -163,7 +163,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -171,7 +172,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel c @@ -179,7 +181,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -187,7 +190,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel e @@ -195,7 +199,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 end subroutine cell_placement_1 @@ -220,7 +225,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -228,7 +234,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel C @@ -236,7 +243,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel d @@ -244,7 +252,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel e @@ -252,7 +261,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%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 f69a688c..71ade79c 100644 --- a/unit-tests/3d/test_mpi_nearest_3.f90 +++ b/unit-tests/3d/test_mpi_nearest_3.f90 @@ -112,7 +112,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel @@ -120,7 +121,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel @@ -128,7 +130,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%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 d601e58c..7cb1ba71 100644 --- a/unit-tests/3d/test_mpi_nearest_4.f90 +++ b/unit-tests/3d/test_mpi_nearest_4.f90 @@ -112,7 +112,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel @@ -120,7 +121,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! big parcel @@ -128,7 +130,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%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 736e0506..63c54858 100644 --- a/unit-tests/3d/test_mpi_nearest_5.f90 +++ b/unit-tests/3d/test_mpi_nearest_5.f90 @@ -111,7 +111,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -119,7 +120,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel a @@ -127,7 +129,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%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 eed620bf..6bd2cc1b 100644 --- a/unit-tests/3d/test_mpi_nearest_6.f90 +++ b/unit-tests/3d/test_mpi_nearest_6.f90 @@ -111,7 +111,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -119,7 +120,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel a @@ -127,7 +129,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%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 e2f6c15f..9977d283 100644 --- a/unit-tests/3d/test_mpi_nearest_7.f90 +++ b/unit-tests/3d/test_mpi_nearest_7.f90 @@ -111,7 +111,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -119,7 +120,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel a @@ -127,7 +129,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%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 4a23bfdd..677f4b83 100644 --- a/unit-tests/3d/test_mpi_nearest_8.f90 +++ b/unit-tests/3d/test_mpi_nearest_8.f90 @@ -111,7 +111,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -119,7 +120,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel a @@ -127,7 +129,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%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 a33e0e7d..fb2fe7af 100644 --- a/unit-tests/3d/test_mpi_nearest_9.f90 +++ b/unit-tests/3d/test_mpi_nearest_9.f90 @@ -111,7 +111,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel b @@ -119,7 +120,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%rank * 100 + l = l + 1 ! small parcel a @@ -127,7 +129,8 @@ 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 + world%rank * 100 + parcels%theta(l) = l + world%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 35993b40..3663bbff 100644 --- a/unit-tests/3d/test_mpi_parcel_communicate.f90 +++ b/unit-tests/3d/test_mpi_parcel_communicate.f90 @@ -96,8 +96,7 @@ program test_mpi_parcel_communicate parcels%volume(1:n_parcels) = cart%rank + 1 parcels%B(:, 1:n_parcels) = cart%rank + 1 parcels%vorticity(:, 1:n_parcels) = cart%rank + 1 - parcels%buoyancy(1:n_parcels) = cart%rank + 1 - + parcels%theta(1:n_parcels) = cart%rank + 1 call parcel_communicate n_total_verify = n_parcels diff --git a/unit-tests/3d/test_mpi_parcel_delete.f90 b/unit-tests/3d/test_mpi_parcel_delete.f90 index 84d5ca66..40190803 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 7814a585..e9b064fc 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 57dbb52b..c69a3524 100644 --- a/unit-tests/3d/test_mpi_parcel_pack.f90 +++ b/unit-tests/3d/test_mpi_parcel_pack.f90 @@ -48,9 +48,10 @@ 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 + parcels%qv(n) = 14.0d0 + a + parcels%ql(n) = 15.0d0 + a #endif enddo @@ -93,9 +94,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 ce9e26f0..f8d4f00e 100644 --- a/unit-tests/3d/test_mpi_parcel_read.f90 +++ b/unit-tests/3d/test_mpi_parcel_read.f90 @@ -47,10 +47,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 @@ -70,9 +71,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') @@ -94,9 +96,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 d1b5da8f..c530e8b4 100644 --- a/unit-tests/3d/test_mpi_parcel_read_rejection.f90 +++ b/unit-tests/3d/test_mpi_parcel_read_rejection.f90 @@ -30,14 +30,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 @@ -102,10 +102,11 @@ program test_mpi_parcel_read_rejection parcels%B(:, 1:n_parcels) = world%rank + 1 parcels%volume(1:n_parcels) = world%rank + 1 parcels%vorticity(:, 1:n_parcels) = world%rank + 1 - parcels%buoyancy(1:n_parcels) = world%rank + 1 + parcels%theta(1:n_parcels) = world%rank + 1 #ifndef ENABLE_DRY_MODE - parcels%humidity(1:n_parcels) = world%rank + 1 + parcels%qv(1:n_parcels) = world%rank + 1 + parcels%ql(1:n_parcels) = world%rank + 1 #endif call create_file('nctest') @@ -124,9 +125,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') @@ -148,9 +150,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 @@ -320,23 +323,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) @@ -389,10 +401,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 86086d2b..a86217a3 100644 --- a/unit-tests/3d/test_mpi_parcel_split.f90 +++ b/unit-tests/3d/test_mpi_parcel_split.f90 @@ -112,7 +112,7 @@ program test_mpi_parcel_split parcels%volume(1:n_parcels) = f12 * vcell parcels%vorticity(:, 1:n_parcels) = world%rank + 1 - parcels%buoyancy(1:n_parcels) = world%rank + 1 + parcels%theta(1:n_parcels) = world%rank + 1 call parcel_split diff --git a/unit-tests/3d/test_mpi_parcel_unpack.f90 b/unit-tests/3d/test_mpi_parcel_unpack.f90 index f45dcd65..878bfe5f 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 a654c941..5b134473 100644 --- a/unit-tests/3d/test_mpi_parcel_write.f90 +++ b/unit-tests/3d/test_mpi_parcel_write.f90 @@ -30,10 +30,11 @@ program test_mpi_parcel_write parcels%B(:, 1:n_parcels) = world%rank parcels%volume(1:n_parcels) = world%rank parcels%vorticity(:, 1:n_parcels) = world%rank - parcels%buoyancy(1:n_parcels) = world%rank + parcels%theta(1:n_parcels) = world%rank #ifndef ENABLE_DRY_MODE - parcels%humidity(1:n_parcels) = world%rank + parcels%qv(1:n_parcels) = world%rank + parcels%ql(1:n_parcels) = world%rank #endif call create_netcdf_parcel_file('nctest', .true., .false.)