documentation and comments

Author: fjansson

README.md

@@ -18,10 +18,21 @@ and the DALES model within OMUSE. See the
 * Activate the OMUSE virtual environment
 
 ```
-python simple-sp.py # run simulations, 10 minutes on 4-core desktop
+python simple-SP.py # run simulations, 10 minutes on 4-core desktop
 python plot-lwp.py  # plot result
 ```
 
+(Instructions tested on Ubuntu 21.04, with OMUSE 2021.6.2.dev14+gf2bbc23).
+
+The script is set up to perform three simulations with a moist bubble perturbation in the initial state:
+
+* a single, wide LES
+
+* a traditional superparameterization with four LES domains
+
+* a modified superparameterization with four LES domains, where the moisture variability is coupled
+
+
 ## References
 
 DALES model description article: Formulation of the Dutch Atmospheric Large-Eddy Simulation (DALES) and overview of its applications, T. Heus et al, [Geosci. Model Dev., 3, 415-444, 2010](https://doi.org/10.5194/gmd-3-415-2010)

plot-lwp.py

@@ -11,7 +11,7 @@
 
 
 # plot results from simple-SP
-# assuming each Dales was run with one process - no merging implemented in this script yet
+# assuming each DALES was run with one process - no merging implemented in this script yet
 # Each row is one time point, columns show different DALES.
 
 mplparams = {"figure.figsize" : [3.47, 3.8],    # figure size in inches
@@ -26,15 +26,8 @@
     }
 plt.rcParams.update(mplparams)
 
-
-#   wishes
-# isolines instead of pixels
-
 C,R = 4,4 # number of columns and rows in plot
 
-field="lwp"
-#field="twp"
-
 vmin={'lwp' : 0,
       'twp' : 35}
 vmax={'lwp' : .2,
@@ -97,17 +90,6 @@ def makeplot(dirname, field, t_start, dtRow, vmin, vmax, colorbar=False):
         dirname = expname+extension
         colorbar = (extension == '-coupled-var-T')
         makeplot(dirname, field, t_start, dtRow, vmin[field], vmax[field], colorbar)
-    
-# measuring average LWP over time
-#for extension in E:
-#    dirname = expname+extension
-#    print (dirname)
-#    for ti in range(0,10):
-#        print('%2d: '%ti, end='')
-#        for i in range(0,C):
-#            print('%5.3f'%lwp_avg[(dirname, i, ti, field)], end=' ')
-#        print()
-#    print()
 
 
 # Plotting the single wide LES run
@@ -144,13 +126,11 @@ def makeplot(dirname, field, t_start, dtRow, vmin, vmax, colorbar=False):
                 lwp=dales_crossxz["qtxz"][ti]
                 im_v=axes[j,i].imshow(lwp, origin='lower', vmin=0,vmax=0.02)
 
-    #cax = fig.add_axes([0.48, 0.05, 0.03, 0.3]) #verttical,  when showing qt also
-    #cax = fig.add_axes([0.85, 0.05, 0.03, 0.3]) #vertical colorbar
     cax = fig.add_axes([0.2, 0.05, 0.3, 0.03]) #horizontal colorbar
     cbar = fig.colorbar(im, cax=cax, orientation='horizontal')
     tick_locator = ticker.MaxNLocator(nbins=4)
     cbar.locator = tick_locator
-    #cbar.ax.yaxis.set_major_locator(ticker.AutoLocator())
+
     cbar.update_ticks()
     cbar.set_label('kg/m$^2$')
     plt.subplots_adjust(left=.01, top=.99, bottom=.1, right=.99, wspace=.01, hspace=.01)

simple-SP.py

@@ -2,7 +2,7 @@
 
 # Simple SP demonstration.
 # A grid of LES models coupled to a global model consisting only of large-scale advection
-# Fredrik Jansson, CWI, 2019-2020
+# Fredrik Jansson, CWI and TU Delft, 2019-2021
 #
 
 
@@ -101,10 +101,6 @@ def evolve(grid, time):
 # upwind advection of quantity q, with velocity u
 # assuming velocity is positive and constant in x, y 
 def advect(q, u, dx, dt):
-    # boundaries are NOT periodic, and probably inconsistent. the leftmost value does not change.
-    #dq = (q[:,:-1, :] - q[:, 1:, :]) * u * dt / dx
-    #q[:, 1:, :] += dq
-
     # periodic boundaries:
     qm = numpy.roll(q.number, 1, axis=1) | q.unit # qm(j,i,k) = q(j,i-1,k), periodic
     dq = (qm - q) * u * dt / dx
@@ -139,9 +135,10 @@ def make_bubble(grid, r, center, gaussian=False):# r, center are quantities, i.e
     else:
         return numpy.where (rr < r*r, 1, 0)
 
+# Adjust the moisture varibility in an LES domain, so that the specific liquid water 
+# profile ql matches the provided ql_ref.
+# This cannot be used before the LES has been stepped - otherwise qsat and ql are not defined.
 def variability_nudge(les, ql_ref, DT, constantT=False):
-    # this cannot be used before the LES has been stepped - otherwise qsat and ql are not defined.
-
     itot, jtot = les.get_itot(), les.get_jtot()
 
     # random field to be used for additive noise when variability is very small
@@ -267,14 +264,14 @@ def get_ql_diff_additive(a):
 
     #qt_std = qt.std(axis=(0, 1))
 
-    
+
+# Performs a superparameterized simulation, where the large-scale model performs only advection using an upwind scheme.
 def run(steps=60, DT=60 | units.s, spinup = 0 | units.s, nx=4, ny=1, n=25, qt_delta=0|units.g/units.kg, name='dales',
         couple=False, bubble=False, bubbleA=1|units.g/units.kg, nudge=None, constantT=False):
 
     dx=100 | units.m     # small-scale grid size
     DX=n*dx              # large-scale grid size = horizontal size of the LES 
-
-    # note CFL on the large scale: need U * DT < DX. U ~= 10 m/s
+                         # note CFL on the large scale: need U * DT < DX. U ~= 10 m/s
 
     grid = [[initBomex(i,j,dirname=name,itot=n,jtot=n,dx=dx,nudge=nudge) for i in range(nx)] for j in range(ny)]
 
@@ -337,9 +334,6 @@ def run(steps=60, DT=60 | units.s, spinup = 0 | units.s, nx=4, ny=1, n=25, qt_de
             print ('QL target profile 2')
             print (QL[0,2,:])
 
-
-
-            
         print ('Evolving to ', DT*ti + spinup)
         evolve(grid, DT*ti + spinup)
 
@@ -354,9 +348,6 @@ def run(steps=60, DT=60 | units.s, spinup = 0 | units.s, nx=4, ny=1, n=25, qt_de
                                   ('z',   units.m)):
                     state[(ti,i,j,var)]   = getattr(grid[j][i].profiles, var).value_in(unit)
         state['time'].append((DT*ti).value_in(units.s))
-
-
-                
 
     filename='result%s.pickle'%('-coupled' if couple else '')
     with open(filename, 'wb') as f:
@@ -421,9 +412,8 @@ def run_single(steps=60, DT=60 | units.s, n=25, nx=4, ny=1, qt_delta=0|units.g/u
 
 # SP with variance nudging at constant thl
 # run       (steps=30, DT=60 | units.s, nx=4, ny=1, couple=True, name='bubble-coupled-var', bubble=True, bubbleA=A, nudge='variance')
-# with variance nudging at constant T
 
-# SP with variance nudging at constant thl
+# SP with variance nudging at constant T
 run       (steps=30, DT=60 | units.s, nx=4, ny=1, couple=True, name='bubble-coupled-var-T', bubble=True, bubbleA=A, nudge='variance', constantT=True)