Feature dispersion tutorial (#46)

* Dispersion tutorial description.

* Some edits to the description.

* Some edits to the description.

* Figs for surface meteorology

* T&D part of the tutorial

* Fix to the figure inclusion

* Small fixes to the offshore tutorial

* Notebooks and parms files for the dispersion tutorial

* Removed the need to execute FE prior to generating the terrain specification file.

* Modified DISPERSION tutorial instructions and example SBL input file.

* Couple more minor adjustments to the DISPERSION case tutroial content and SBL case parameters file.

* Fix typo in the open_dataset_line

---------

Co-authored-by: Jeremy Sauer <[email protected]>

docs/Tutorials/cases.rst

@@ -2,7 +2,7 @@
 Ideal Test Cases
 ****************
 
-Four test cases are described:
+Seven test cases are described:
 
 * Dry neutral boundary layer
 * Dry convective boundary layer
@@ -19,3 +19,4 @@ Required tutorial resources including python utilities and Jupyter Notebooks are
    cases/MBL.rst
    cases/CANOPY.rst
    cases/OFFSHORE.rst
+   cases/DISPERSION.rst

docs/Tutorials/cases/CANOPY.rst

@@ -11,7 +11,7 @@ Input parameters
 ----------------
 
 * Number of grid points: :math:`[N_x,N_y,N_z]=[252,250,90]`
-* Isotropic grid spacings in the horizontal directions: :math:`[dx,dy]=[4,4]` m, vertical grid is :math:`dz=4` m at the surface and stretched with verticalDeformFactor :math:`=0.25`
+* Isotropic grid spacings in the horizontal directions: :math:`[\Delta x,\Delta y]=[4,4]` m, vertical grid is :math:`\Delta z=4` m at the surface and stretched with verticalDeformFactor :math:`=0.25`
 * Domain size: :math:`[1.0 \times 1.0 \times 1.44]` km
 * Model time step: :math:`0.01` s
 * Advection scheme: 5th-order upwind

docs/Tutorials/cases/DISPERSION.rst

@@ -0,0 +1,100 @@
+==================================================
+Transport and dispersion over idealized topography
+==================================================
+
+Background
+----------
+
+This is an example of transport and dispersion of a scalar in the presence of idealized topography. The idealized terrain is based on the Witch of Agnesi, and a passive tracer is released on the lee side of the hill.
+
+Input parameters
+----------------
+
+* Number of grid points: :math:`[N_x,N_y,N_z]=[504,498,90]`
+* Isotropic grid spacings in the horizontal directions: :math:`[\Delta x,\Delta y]=[4,4]` m, the minimum vertical grid at the surface is :math:`\Delta z=3.6` m and stretched with verticalDeformFactor :math:`=0.26`
+* Domain size: :math:`[2.16 \times 1.99 \times 1.44]` km
+* Model time step: :math:`0.01` s
+* Advection scheme: 5th-order upwind
+* Time scheme: 3rd-order Runge Kutta
+* Geostrophic wind: :math:`[U_g,V_g]=[8,0]` m/s
+* Latitude: :math:`54.0^{\circ}` N
+* Surface potential temperature: :math:`300` K
+* Potential temperature profile:
+
+.. math::
+  \partial{\theta}/\partial z =
+    \begin{cases}
+      0 & \text{if $z$ $\le$ 500 m}\\
+      0.003 & \text{if $z$ > 500 m}
+    \end{cases} 
+
+* Rayleigh damping layer: uppermost :math:`600` m of the domain
+* Initial perturbations: :math:`\pm 0.25` K 
+* Depth of perturbations: :math:`375` m
+* Top boundary condition: free slip
+* Lateral boundary conditions: periodic
+* Time period: :math:`1` h
+
+Execute FastEddy
+----------------
+
+Note that this example requires customization of the initial condition file. Follow the sequence of steps below.
+
+#. Exectue the Jupyer notebook provided in **/tutorial/notebooks/Dispersion_Prep1.ipynb** to create the topography file *Topography_504x498.dat*. 
+#. Run FastEddy for 1 timestep using the default state of the (*Example07_DISPERSION_SBL.in*) and required binary terrain file generated in the previous step, specified as input through the FastEddy parameter file (:code:`topoFile`). This step creates an output file *FE_DISPERSION.0* that includes the topography and establishes a terrain-following vertical coordinate grid. 
+#. Exectute the Jupyer notebook provided in **/tutorial/notebooks/Dispersion_Prep2.ipynb** to modify the surface roughness distribution over the hill. 
+#. Execute the Jupyer notebook provided in **/tutorial/notebooks/Dispersion_Prep3.ipynb** to create the source specification input file. 
+#. Adjust the (*Example07_DISPERSION_SBL.in*) file to specify the targeted initial condition file (:code:`inPath`, :code:`inFile`) by removing the (:code:`#`) just to the right of the equal sign ito uncomment these parameters values. Uncomment the pre-formed passive tracer source file (:code:`srcAuxScFile`). Remove the value of :math:`0` and uncomment the value of :math:`2` for the number of source emissions (:code:`NhydroAuxScalars`).
+#. Run FastEddy for :math:`1` h of the simulation by changing :code:`frqOutput`, :code:`Nt`, and :code:`NtBatch` removing the values of :math:`1` for each and the (:code:`#`) to uncomment appropriate full-simulation parameters values. The emissions begin :math:`45` min into the simulation.  
+
+Two FastEddy simulation setups are provided for this tutorial, corresponding to weakly stable (*Example07_DISPERSION_SBL.in*) and convective conditions (*Example07_DISPERSION_CBL.in*). The initial condition and terrain preparation steps only need to be carried out once for the stable case, then can be reused directly in the convective stability case.
+
+See :ref:`run_fasteddy` for instructions on how to build and run FastEddy on NSF NCAR's High Performance Computing machines.
+
+Visualize the output
+--------------------
+
+Open the Jupyter notebook entitled *MAKE_FE_TUTORIAL_PLOTS.ipynb* and execute it using setting: :code:`case = 'dispersion'`.
+
+XY-plane views of instantaneous velocity components and potential temperature for the SBL case at :math:`t=1` h (FE_DISPERSION.360000). The contour lines in the :math:`u` panel display terrain elevation:
+
+.. image:: ../images/UVWTHETA-XY-dispersion_SBL.png
+  :width: 1200
+  :alt: Alternative text
+
+XY-plane views of instantaneous velocity components and potential temperature for the CBL case at :math:`t=1` h (FE_DISPERSION.360000). The contour lines in the :math:`u` panel display terrain elevation:
+
+.. image:: ../images/UVWTHETA-XY-dispersion_CBL.png
+  :width: 1200
+  :alt: Alternative text
+
+XY-plane views of instantaneous plume dispersion for the SBL case at :math:`z=30` m AGL and different times (:math:`t=50,55,60` min), corresponding to the windward release:
+
+.. image:: ../images/CONCENTRATION-XY-dispersion_SBL.png
+  :width: 1200
+  :alt: Alternative text
+
+XY-plane views of instantaneous plume dispersion for the CBL case at :math:`z=30` m AGL and different times (:math:`t=50,55,60` min), corresponding to the windward release:
+
+.. image:: ../images/CONCENTRATION-XY-dispersion_CBL.png
+  :width: 1200
+  :alt: Alternative text
+
+YZ-plane views of instantaneous plume dispersion for the SBL case at several downstream distances (:math:`t=1` h, FE_DISPERSION.360000), corresponding to the windward release:
+
+.. image:: ../images/CONCENTRATION-YZ-dispersion_SBL.png
+  :width: 1200
+  :alt: Alternative text
+
+YZ-plane views of instantaneous plume dispersion for the CBL case at several downstream distances (:math:`t=1` h, FE_DISPERSION.360000), corresponding to the windward release:
+
+.. image:: ../images/CONCENTRATION-YZ-dispersion_CBL.png
+  :width: 1200
+  :alt: Alternative text
+
+Analyze the output
+------------------
+
+* How does the terrain impact gets altered by the different stability conditions?
+* What are the differences in plume dispersion between stable and convective condtions?
+* How does downstream distance affect structure of the plume?

docs/Tutorials/cases/OFFSHORE.rst

@@ -11,7 +11,7 @@ Input parameters
 ----------------
 
 * Number of grid points: :math:`[N_x,N_y,N_z]=[360,362,90]`
-* Isotropic grid spacings in the horizontal directions: :math:`[dx,dy]=[15,15]` m, vertical grid is :math:`dz=15` m at the surface and stretched with verticalDeformFactor :math:`=0.50`
+* Isotropic grid spacings in the horizontal directions: :math:`[\Delta x,\Delta y]=[15,15]` m, vertical grid is :math:`\Delta z=15` m at the surface and stretched with verticalDeformFactor :math:`=0.50`
 * Domain size: :math:`[5.40 \times 5.43 \times 2.70]` km
 * Model time step: :math:`0.04` s
 * Advection scheme: 5th-order upwind
@@ -29,7 +29,7 @@ Input parameters
     \end{cases} 
 
 * Surface warming rate:  :math:`0.5` K/h
-* Surface roughness length: Donelan (1990) parameterization
+* Surface roughness length: Drennan (2003) parameterization
 * Rayleigh damping layer: uppermost :math:`400` m of the domain
 * Initial perturbations: :math:`\pm 0.50` K 
 * Depth of perturbations: :math:`825` m

tutorials/examples/Example06_OFFSHORE.in

@@ -90,7 +90,7 @@ surflayer_ideal_qte = 10.0 # end time in seconds for the idealized sinusoidal su
 surflayer_ideal_qamp = 2.0 # maximum amplitude of the idealized sinusoidal surface forcing (qv)
 surflayer_qskin_input = 1 # selector to use file input (restart) value for qskin under surflayerSelector == 2
 surflayer_offshore = 1 # offshore selector: 0=off, 1=on
-surflayer_offshore_opt = 3 # offshore roughness parameterization: ==0 (Charnock), ==1 (Charnock with variable alpha), ==2 (Taylor & Yelland), ==3 (Donelan), ==4 (Drennan), ==5 (Porchetta)
+surflayer_offshore_opt = 4 # offshore roughness parameterization: ==0 (Charnock), ==1 (Charnock with variable alpha), ==2 (Taylor & Yelland), ==3 (Donelan), ==4 (Drennan), ==5 (Porchetta)
 #------------: BASE-STATE ---
 stabilityScheme = 2 # Scheme used to set hydrostatic, stability-dependent Base-State EOS fields
 temp_grnd = 288.0 # Air Temperature (K) at the ground used to set hydrostatic Base-State EOS fields

tutorials/examples/Example07_DISPERSION_CBL.in

@@ -0,0 +1,142 @@
+Description = This is an idealized CBL over terrain with T&D of a passive scalar.
+#--MPI_AALES
+numProcsX = 8 # 4 # Number of cores to be used for horizontal domain decomposition in X
+numProcsY = 1 # Number of cores to be used for horizontal domain decomposition in Y
+#--CUDA_AALES
+tBx = 1 # Number of threads in x-dimension
+tBy = 8 # Number of threads in y-dimension
+tBz = 32 # Number of threads in z-dimension
+#--IO
+inPath = ./initial/ # Path where initial/restart file is read in from
+inFile = FE_DISPERSION.0 # name of the input file for coordinate system and initial or restart conditions
+outPath = ./output/ # ./output/ # Path where output files are to be written
+outFileBase = FE_DISPERSION # Base name of the output file series as in (outFileBase).element-in-series
+frqOutput = 30000 # frequency (in timesteps) at which to produce output
+ioOutputMode = 0 # 0: N-to-1 gather and write to a netcdf file, 1:N-to-N writes of FastEddy binary files
+#--GRID
+Nx = 504 # 500 # Number of discretised domain elements in the x (zonal) direction
+Ny = 498 # Number of discretised domain elements in the y (meridional) direction
+Nz = 90 # Number of discretised domain elements in the z (vertical) direction
+Nh = 3 # Number of halo cells to be used (dependent on largest stencil extent)
+d_xi = 4.0 # Computational domain fixed resolution in the 'i' direction
+d_eta = 4.0 # Computational domain fixed resolution in the 'j' direction
+d_zeta = 16.0 # Computational domain fixed resolution in the 'k' direction
+coordHorizHalos = 1 # switch to setup coordiante halos as periodic=1 or gradient-following=0
+topoFile = ./Topography_504x498.dat # A file containing topography (surface elevation in meters ASL)
+verticalDeformSwitch = 1 # switch to use vertical coordinate deformation 0=off, 1=on
+verticalDeformFactor = 0.26 # deformation factor (0.0=max compression,  1.0=no compression)
+verticalDeformQuadCoeff = 0.0 # deformation factor (0.0=max compression,  1.0=no compression)
+#--TIME_INTEGRATION
+timeMethod = 0 # Selector for time integration method. [0=RK3-WS2002 (default)]
+Nt = 360000 # Number of timesteps to perform
+dt = 0.01 # timestep resolution in seconds
+NtBatch = 30000 # Number of timesteps to compute in batch launch, must have NtBatch <= Nt
+#--HYDRO_CORE
+##------------: HYDRO_CORE Submodule Selectors ---
+#----------: Boundary Conditions Set ---
+hydroBCs =  2 # Selector for hydro BC set. 2= periodicHorizVerticalAbl
+##----------: HYDRO_IO/LOGGGING ---
+hydroForcingWrite = 0 # Switch for dumping hydroFldsFrhs for prognositic fields. 0 = off, 1=on
+hydroSubGridWrite = 0 # Switch for dumping Tauij fields. 0 = off, 1=on
+hydroForcingLog = 0 # switch for logging Frhs min/max, etc.
+##----------: ADVECTION ---
+advectionSelector = 3 # advection scheme selector: 0= 1st-order upwind, 1= 3rd-order QUICK, 2= hybrid 3rd-4th order, 3= hybrid 5th-6th order
+b_hyb = 0.0 # hybrid advection scheme parameter: 0.0= lower-order upwind, 1.0=higher-order cetered, 0.0 < b_hyb < 1.0 = hybrid
+##------------: MOISTURE ---
+moistureSelector = 0 # moisture selector: 0=off, 1=on
+moistureNvars = 2 # number of moisture species
+moistureAdvSelectorQv = 2 # water vapor advection scheme selector
+moistureAdvSelectorQv_b = 0.0 # hybrid advection scheme parameter for water vapor
+moistureAdvSelectorQi = 2 # moisture advection scheme selector for non-qv fields (non-oscillatory schemes)
+moistureSGSturb = 1 # selector to apply sub-grid scale diffusion to moisture fields
+moistureCond = 3 # selector to apply condensation to mositure fields
+moistureCondTscale = 1.0 # relaxation time in seconds
+moistureCondBasePres = 1 #selector to use base pressure for microphysics
+moistureMPcallTscale = 1.0 # time scale for microphysics to be called
+##------------: CORIOLIS ---
+coriolisSelector = 1  #Coriolis switch: 0 off, 1 on
+coriolisLatitude = 54.0 # Charactersitc latitude in degrees from equator of the LES domain
+## ----------: TURBULENCE ---
+turbulenceSelector = 1 # turbulence scheme selector: 0= none, 1= Lilly/Smagorinsky
+TKESelector = 1 # Prognostic TKE selector: 0= none, 1= Prognostic
+TKEAdvSelector = 3 # SGSTKE advection selector
+TKEAdvSelector_b_hyb = 0.0
+c_s = 0.18 # Smagorinsky turbulence model constant used for turbulenceSelector = 1 with TKESelector = 0
+c_k = 0.10 # Lilly turbulence model constant used for turbulenceSelector = 1 with TKESelector > 0
+##------------: CANOPY ---
+canopySelector = 0 # canopy selector: 0=off, 1=on
+canopySkinOpt = 0 # canopy selector to use additional skin friction effect on drag coefficient: 0=off, 1=on
+canopy_cd = 0.15 # non-dimensional canopy drag coefficient cd coefficient
+canopy_lf = 0.1 # representative canopy element length scale
+##------------: DIFFUSION ---
+diffusionSelector = 0 # diffusivity selector: 0= none, 1= const.
+nu_0 = 0.0 # constant diffusivity used when diffusionSelector = 1
+##----------: AUXILIARY SCALARS ---
+NhydroAuxScalars = 2 # Number of prognostic auxiliary scalar fields
+AuxScAdvSelector = 0 # advection scheme for auxiliary scalar fields
+AuxScAdvSelector_b_hyb = 0.0
+AuxScSGSturb = 1 # selector to apply sub-grid scale diffusion to auxiliary scalar fields
+## ----------: AUXILIARY SCALAR SOURCES ---
+srcAuxScFile=./Example07_sources.nc
+##------------: EXPLICIT FILTERS ---
+filterSelector = 1 # explicit filter selector: 0=off, 1=on
+filter_6thdiff_vert = 1
+filter_6thdiff_vert_coeff = 0.03 # vertical 6th-order filter factor: 0.0=off, 1.0=full
+##----------: RAYLEIGH DAMPING LAYER ---
+dampingLayerSelector = 1 # Rayleigh damping layer selector: 0= off, 1= on.
+dampingLayerDepth = 600.0 # Rayleigh damping layer depth in meters
+##----------: SURFACE LAYER ---
+surflayerSelector = 1 # surfacelayer selector: 0= off, 1= surface kinematic heat flux (surflayer_wth), 2=sking temperature rate (surflayer_tr)
+surflayer_z0 = 0.1 # roughness length (momentum) when surflayerSelector > 0
+surflayer_z0t = 0.1 # roughness length (temperature) when surflayerSelector > 0
+surflayer_z0tdyn = 1 # dynamic z0t calculation following Zilitinkevich (1995) approach: 0= off, 1= constant Zilitinkevich coeff, 2= variable Zilitinkevich coeff
+surflayer_wth = 0.25 # kinematic sensible heat flux at the surface when surflayerSelector = 1
+surflayer_tr = -0.50 # temperature rate at the surface when surflayerSelector = 2 (>0 for warming; <0 for cooling)
+surflayer_wq = 0.0 # sensible heat flux at the surface (kg/kg m s-1) when surflayerSelector = 1
+surflayer_qr = 0.0 # water vapor rate (kg/kg h-1) when surflayerSelector = 2
+surflayer_idealsine = 0 # selector for idealized sinusoidal surface heat flux or skin temperature forcing
+surflayer_ideal_ts = 0.0 # start time in seconds for the idealized sinusoidal surface forcing
+surflayer_ideal_te = 43200 # end time in seconds for the idealized sinusoidal surface forcing
+surflayer_ideal_amp = 4.0 # maximum amplitude of the idealized sinusoidal surface forcing
+surflayer_ideal_qts = 1.75 # start time in seconds for the idealized sinusoidal surface forcing (qv)
+surflayer_ideal_qte = 10.0 # end time in seconds for the idealized sinusoidal surface forcing (qv)
+surflayer_ideal_qamp = 2.0 # maximum amplitude of the idealized sinusoidal surface forcing (qv)
+surflayer_offshore = 0 # offshore selector: 0=off, 1=on
+#------------: BASE-STATE ---
+stabilityScheme = 2 # Scheme used to set hydrostatic, stability-dependent Base-State EOS fields
+temp_grnd = 300.0 # Air Temperature (K) at the ground used to set hydrostatic Base-State EOS fields
+pres_grnd = 100000.0 # Pressure (Pa) at the ground used to set hydrostatic Base-State EOS fields
+zStableBottom = 500.0 # Height (m) of the first stable upper-layer when stabilityScheme = 1 or 2
+stableGradient = 0.003 # 0.08 # Vertical gradient (K/m) of the first stable upper-layer when stabilityScheme = 1 or 2
+zStableBottom2 = 650.0 # Height (m) of the second stable upper-layer when stabilityScheme = 2
+stableGradient2 = 0.003 # Vertical gradient (K/m) of the second stable upper-layer when stabilityScheme = 2
+zStableBottom3 = 50000.0 # Height (m) of the third stable upper-layer when stabilityScheme = 2
+stableGradient3 = 0.003 # Vertical gradient (K/m) of the third stable upper-layer when stabilityScheme = 2
+thetaPerturbationSwitch = 1 # Switch to include initial theta perturbations: 0=off, 1=on
+thetaHeight = 375.0 # Height below which to include initial theta perturbations: (meters)
+thetaAmplitude = 0.25 # Maximum amplitude for theta perturbations: thetaAmplitude*[-1,+1] K
+U_g = 8.0 # Zonal (West-East) component of the geostrophic wind (m/s)
+V_g = 0.0 # Meridional (South-North) component of the geostrophic wind (m/s)
+z_Ug = 10000.0 # Height (m) above ground for linear geostrophic wind gradient (zonal component)
+z_Vg = 10000.0 # Height (m) above ground for linear geostrophic wind gradient (meridional component)
+Ug_grad = 0.0 # U_g gradient above z_Ug (ms-1/m)
+Vg_grad = 0.0 # V_g gradient above z_Vg (ms-1/m)
+# ----------: LARGE SCALE FORCINGS ---
+lsfSelector = 0 # large-scale forcings selector: 0=off, 1=on
+lsf_horMnSubTerms = 0 # large-scale subsidence terms Switch: 0= off, 1= on
+lsf_freq = 1.0 # large-scale forcing frequency (seconds)
+lsf_w_surf = 0.0 # lsf to w at the surface
+lsf_w_lev1 = -23.4 # lsf to w at the first specified level
+lsf_w_lev2 = 0.0 # lsf to w at the second specified level
+lsf_w_zlev1 = 1500.0 # lsf to w height 1
+lsf_w_zlev2 = 2100.0 # lsf to w height 2
+lsf_th_surf = -0.0833333 # lsf to theta at the surface
+lsf_th_lev1 = -0.0833333 # lsf to theta at the first specified level
+lsf_th_lev2 = 0.0 # lsf to theta at the second specified level
+lsf_th_zlev1 = 1500.0 # lsf to theta height 1
+lsf_th_zlev2 = 3000.0 # lsf to theta height 2
+lsf_qv_surf = -0.0432 # large-scale forcing to qv at the first specified level
+lsf_qv_lev1 = -0.0432 # large-scale forcing qv at height 1
+lsf_qv_lev2 = 0.0 # large-scale forcing qv at height 2
+lsf_qv_zlev1 = 300.0 # large-scale forcing qv height 1
+lsf_qv_zlev2 = 500.0 # large-scale forcing qv height 2