Changes for BOMEX..PR in progress.

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

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
+/

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

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)
 

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)

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)