Thermo on fixed parallelism

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

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_parcel_merge_random.f90

@@ -90,12 +90,12 @@ program test_parcel_merge_random
                 call random_number(rn(3))
                 j = nint(n_parcels * rn(3)) + 1
                 parcels%volume(j) = 0.9d0 * vmin
-                parcels%buoyancy(j) = 1.0d0
+                parcels%theta(j) = 1.0d0
             endif
 
         enddo
 
-        n_merges = int(sum(parcels%buoyancy(1:n_parcels)))
+        n_merges = int(sum(parcels%theta(1:n_parcels)))
         call perform_integer_reduction(n_merges)
 
         if (world%rank == world%root) then
@@ -117,7 +117,7 @@ program test_parcel_merge_random
         n_parcels = n_orig
         do n = 1, n_parcels
             parcels%volume(n) = vol
-            parcels%buoyancy(n) = 0.0d0
+            parcels%theta(n) = 0.0d0
             parcels%B(:, n) = b
 
             call random_number(rn)

mpi-tests/test_parcel_split_random.f90

@@ -88,12 +88,12 @@ program test_parcel_split_random
                 call random_number(rn(3))
                 j = nint(n_parcels * rn(3)) + 1
                 parcels%B(1, j) = 5.0d0 * parcels%B(1, j)
-                parcels%buoyancy(j) = 1.0d0
+                parcels%theta(j) = 1.0d0
             endif
 
         enddo
 
-        n_splits = int(sum(parcels%buoyancy(1:n_parcels)))
+        n_splits = int(sum(parcels%theta(1:n_parcels)))
         call perform_integer_reduction(n_splits)
 
         if (world%rank == world%root) then
@@ -115,7 +115,7 @@ program test_parcel_split_random
         n_parcels = n_orig
         do n = 1, n_parcels
             parcels%volume(n) = vol
-            parcels%buoyancy(n) = 0.0d0
+            parcels%theta(n) = 0.0d0
             parcels%B(:, n) = b
 
             call random_number(rn)

python-scripts/bomex_fluxes.py

@@ -0,0 +1,152 @@
+#!/usr/bin/env python
+from tools.nc_fields import nc_fields
+import numpy as np
+import argparse
+
+try:
+    parser = argparse.ArgumentParser(
+        description="Create surface flux setup"
+    )
+
+    parser.add_argument(
+        "--nx",
+        type=int,
+        required=False,
+        default=192,
+        help="number of cells in x",
+    )
+
+    parser.add_argument(
+        "--ny",
+        type=int,
+        required=False,
+        default=192,
+        help="number of cells in y",
+    )
+
+    parser.add_argument(
+        "--nz",
+        type=int,
+        required=False,
+        default=96,
+        help="number of cells in z",
+    )
+
+    parser.add_argument(
+        "--L",
+        type=float,
+        required=False,
+        default=3200.0,
+        help="domain extent",
+    )
+
+    args = parser.parse_args()
+
+    thetafluxmag=8.0e-3
+    qvfluxmag=5.2e-5
+
+    # number of cells
+    nx = args.nx
+    ny = args.ny
+    nz = args.nz
+
+    L = args.L
+
+    ncf = nc_fields()
+
+    ncf.open('bomex_setup' + str(nx) + 'x' + str(ny) + 'x' + str(nz) + '.nc')
+
+    # domain origin
+    origin = (-1.0 * L, -1.0 * L, 0.0)
+
+    # domain extent
+    extent = (2.0 * L, 2.0 * L, L)
+
+
+    # mesh spacings
+    dx = extent[0] / nx
+    dy = extent[1] / ny
+    dz = extent[2] / nz
+
+    theta = np.zeros((nz+1, ny, nx))
+    qv = np.zeros((nz+1, ny, nx))
+    yvort = np.zeros((nz+1, ny, nx))
+    novort = np.zeros((nz+1, ny, nx))
+
+    # ranges from 0 to nz
+    for k in range(nz+1):
+        zz = k * dz
+        if(zz < 520.):
+            theta[k, :,: ] = 298.7
+        elif(zz < 1480):
+            theta[k, :, :] = 298.7 + (zz-520.)*(302.4-298.7)/(1480.-520.)
+        elif(zz < 2000):
+            theta[k, :, :] = 302.4 + (zz-1480.)*(308.2-302.4)/(2000.-1480.)
+        else:
+            theta[k, :, :] = 308.2 + (zz-2000.)*(311.85-308.2)/(3000.-2000.)
+
+        # specific humidity
+        if(zz < 520.):
+            qv[k, :, :] = 1e-3*(17.0 + zz*(16.3-17.0)/520.)
+        elif(zz < 1480):
+            qv[k, :, :] = 1.e-3*(16.3 + (zz-520.)*(10.7-16.3)/(1480.-520.))
+        elif(zz < 2000):
+            qv[k, :, :] = 1.e-3*(10.7 + (zz-1480.)*(4.2-10.7)/(2000.-1480.))
+        else:
+            qv[k, :, :] = 1.e-3*(4.2 + (zz-2000.)*(3.-4.2)/(3000.-2000.))
+
+        if(zz > 699.):
+            yvort[k, :, :] = 0.5*(-4.61+8.75)/(3000.-700.)
+        if(zz > 701.):
+            yvort[k, :, :] = (-4.61+8.75)/(3000.-700.)
+
+    # add perturbation
+    rng = np.random.default_rng(seed=42)
+    noise = rng.uniform(low=-0.05, high=0.05, size=(nz+1, ny, nx))
+    theta += noise
+
+    # write all provided fields
+    ncf.add_field('x_vorticity', novort, unit='1/s')
+    ncf.add_field('y_vorticity', yvort, unit='1/s')
+    ncf.add_field('z_vorticity', novort, unit='1/s')
+    ncf.add_field('theta', theta, unit='K', long_name='potential temperature')
+    ncf.add_field('qv', qv, unit='kg/kg', long_name='water vapour spec. hum.')
+
+    ncf.add_box(origin, extent, [nx, ny, nz])
+
+    ncf.close()
+
+    #
+    # Setup surface fluxes:
+    #
+
+    ncf_flux = nc_fields()
+
+    ncf_flux.set_dim_names(['x', 'y'])
+
+    ncf_flux.open('bomex_flux_' + str(nx) + 'x' + str(ny) + '.nc')
+
+    # domain origin
+    origin = (-1.0 * L, -1.0 * L)
+
+    # domain extent
+    extent = (2.0 * L, 2.0 * L)
+
+
+    # mesh spacings
+    dx = extent[0] / nx
+    dy = extent[1] / ny
+
+    thetaflux = np.ones((ny, nx))*thetafluxmag
+    qvflux = np.ones((ny, nx))*qvfluxmag
+
+    ncf_flux.add_field('thetaflux', thetaflux, unit='K m/s')
+    ncf_flux.add_field('qvflux', qvflux, unit='kg/kg m/s')
+
+    ncf_flux.add_box(origin, extent, [nx, ny])
+
+    ncf_flux.close()
+
+except Exception as err:
+    print(err)
+