diff --git a/docs/Tutorials/cases.rst b/docs/Tutorials/cases.rst
index f3a3461..8657d63 100644
--- a/docs/Tutorials/cases.rst
+++ b/docs/Tutorials/cases.rst
@@ -1,8 +1,8 @@
-****************
-Ideal Test Cases
-****************
+***********
+Ideal Cases
+***********
 
-Seven test cases are described:
+Seven idealized cases are provided:
 
 * Dry neutral boundary layer
 * Dry convective boundary layer
@@ -12,7 +12,7 @@ Seven test cases are described:
 * Offshore boundary layer
 * Passive scalar transport and dispersion over an idealized hill
 
-Required tutorial resources including python utilities and Jupyter Notebooks are provided in the tutorials directory of the `FastEddy-model GitHub repository <https://github.com/NCAR/FastEddy-model>`_ with required data for the moist dynamics example available at this `Zenodo record <https://zenodo.org/records/10982246>`_. All test cases are idealized setups over flat terrain. For each case, the user will set up the input parameter file, execute FastEddy, visualize the output using a Jupyter notebook, and perform some basic analysis of the output. After examining the test cases, the user will carry out some sensitivity tests by changing various input parameters. The purpose of these tests are for the user to become more familiar with the input parameters, and how changes to those parameters affect the output. After the tutorial, the user is expected to have basic knowledge to carry out LES using FastEddy.
+Resources required for tutorial cases including python utilities and Jupyter Notebooks are provided in the tutorials directory of the `FastEddy-model GitHub repository <https://github.com/NCAR/FastEddy-model>`_ with required data for the moist dynamics example available at this `Zenodo record <https://zenodo.org/records/10982246>`_. All cases are idealized setups over flat terrain. For each case, the user will set up the input parameter file, execute FastEddy, visualize the output using a Jupyter notebook, and perform some basic analysis of the output. After examining the cases, the user will carry out some sensitivity tests by changing various input parameters. The purpose of these tests are for the user to become more familiar with the input parameters, and how changes to those parameters affect the output. After completing these tutorials, the user will have some basic knowledge and capacity to carry out LES activites using FastEddy.
 
 .. toctree::
 
@@ -23,3 +23,26 @@ Required tutorial resources including python utilities and Jupyter Notebooks are
    cases/CANOPY.rst
    cases/OFFSHORE.rst
    cases/DISPERSION.rst
+
+.. rubric:: Sensitivity Tests
+
+* Re-run the neutral case with :math:`[N_x,N_y,N_z]=[400,400,122]` and isotropic grid spacings of :math:`[dx,dy,dz]=[10,10,10]`. Adjust the model time step accordingly. Re-make all plots and discuss the differences between the control case. How much longer did it take to complete the simulation?
+* Re-run the convective case with a surface heat flux of :math:`=+0.70` Km/s. Re-make all plots and discuss the differences between the control case.
+* Re-run the neutral case with :math:`z_0=0.3` m. Re-make all plots and discuss the differences between the control case.
+* Re-run the neutral case with the first order upwind advection scheme. Re-make all plots and discuss the differences between the control case. Why is the first order scheme a bad choice?
+* Re-run the stable case with a surface cooling rate of :math:`-0.5` K/h. Re-make all plots and discuss the differences between the control case.
+* Re-run the stable case using half of the GPUs used in the control simulation. How much slower does the case run?
+
+* Re-run the BOMEX case with a higher-order advection for water vapor (moistureAdvSelectorQv = 3). What is the impact of the increased effective resolution on dynamical, thermodynamical and microphysical quantities, along with turbulence variability and fluxes? How does that change influce the comparison to the other BOMEX LES models?
+
+.. only
+
+    (Here, the user will make some modifications to the default parameters such as changing the grid spacing, stretching, model time step, advection     
+
+    scheme,    number of grid points, domain decomposition and number of GPUs, etc, etc. Here, the user will execute the sensitivity test, and visualize 
+    and analyze the     output)
+
+     x^2+y^2=z^2
+
+      frac{ sum_{t=0}^{N}f(t,k) }{N}
+
diff --git a/docs/Tutorials/cases_extensions.rst b/docs/Tutorials/cases_extensions.rst
new file mode 100644
index 0000000..97ce2d8
--- /dev/null
+++ b/docs/Tutorials/cases_extensions.rst
@@ -0,0 +1,10 @@
+*********************
+Model Extension Cases
+*********************
+
+Two model extension cases are provided excercising the Generalized Actuator Disk (GAD) parameterized model and building-resolving URBAN model extensions respectively.
+
+.. toctree::
+
+   cases_extensions/GAD.rst
+   cases_extensions/WRF_coupling_case0_FE_urban.rst
diff --git a/docs/Tutorials/cases_extensions/GAD.rst b/docs/Tutorials/cases_extensions/GAD.rst
new file mode 100644
index 0000000..1643787
--- /dev/null
+++ b/docs/Tutorials/cases_extensions/GAD.rst
@@ -0,0 +1,64 @@
+=====================================
+Turbine array flow with the GAD model
+=====================================
+
+This is an idealized scenario of wind farm (3 x 3 turbine array) flow in neutrally stratified boundary layer with an evolving wind direction spanning 90 degrees of clockwise turning over 40 minutes. This idealized scenario demonstrates the generalized actuator disk (GAD) implementation in FastEddy (*Sanchez Gomez et al., 2024*), with the inclusion of a turbine yawing capability to align with the meteorological wind direction at the turbine's nacelle. The initial and boundary conditions for this idealized case are derived from a horizontally averaged LES run of a neutral ABL with a geostrophic wind aligned in the zonal direction (:math:`[U_g,V_g]=[10.0,0.0]` m/s) with and a latitude of :math:`40.0^{\circ}` N. The required datasets to run this tutorial are provided at this Zenodo record [TO BE UPDATED].
+
+The GAD model is activated by the selector (:code:`GADSelector = 1`) in the parameters file. This case uses the parameters file **tutorials/examples/Example09_GAD.in**. The lines below correspond to additions to the FastEddy parameters file necessary for turbine-inclusive LES runs. These correspond to the turbine specification file (:code:`turbineSpecsFile`) and a parameter to choose whether or not to write GAD forces to output files (:code:`GADoutputForces`).
+
+.. code-block:: none
+
+   #--GAD
+   GADSelector = 1
+   turbineSpecsFile = ./GAD_NREL28_9WTs_tutorial.nc
+   GADoutputForces = 1
+
+A key aspect required for the GAD model to work is the specification of the aerodynamic characteristics of the simulated turbine. The turbine specifications file is a netCDF formatted file that provides the aerodynamic characteristics alongside the geometric characteristics of the wind turbine (rotor diameter: :code:`GAD_rotorD`, turbine hub height: :code:`GAD_hubHeights`, nacelle diameter: :code:`GAD_nacelleD`). The turbine specifications file also informs the location and initial orientation of the turbines (:code:`GAD_Xcoords`, :code:`GAD_Ycoords`, :code:`GAD_rotorTheta`). Aerodynamic characteristics are given as polynomial fits for lift and drag coefficient, twist, chord length, blade pitch, and rotational speed discretized over a finite number of normalized blade elements (:code:`rnorm_vect`), following the blade-element momentum theory used in the GAD formulation. For flexibility purposes, an arbitrary number of turbines (:code:`GAD_turbineType`) can be defined. This tutorial provides an example turbine specification file corresponding to the U.S. Department of Energy NREL28 turbine (**GAD_NREL28_9WTs_tutorial.nc**). All the required variables and dimensions in the turbine specifications file are listed here below.
+
+.. code-block:: none
+
+   float GAD_Xcoords(GADNumTurbines) ;
+   float GAD_Ycoords(GADNumTurbines) ;
+   float GAD_rotorTheta(GADNumTurbines) ;
+   int GAD_turbineType(GADNumTurbines) ;
+   int GADNumTurbineTypes(GADNumTurbineTypes) ;
+   int turbinePolyClCdrNormBounds(turbinePolyClCdrNormBounds) ;
+   int turbinePolyClCdrNormSegments(turbinePolyClCdrNormSegments) ;
+   int alphaBounds(alphaBounds) ;
+   int turbinePolyOrderMax(turbinePolyOrderMax) ;
+   float GAD_hubHeights(GADNumTurbineTypes) ;
+   float GAD_rotorD(GADNumTurbineTypes) ;
+   float GAD_nacelleD(GADNumTurbineTypes) ;
+   float rnorm_vect(GADNumTurbineTypes, turbinePolyClCdrNormBounds) ;
+   float alpha_minmax_vect(GADNumTurbineTypes, alphaBounds) ;
+   float turbinePolyTwist(GADNumTurbineTypes, turbinePolyOrderMax) ;
+   float turbinePolyChord(GADNumTurbineTypes, turbinePolyOrderMax) ;
+   float turbinePolyPitch(GADNumTurbineTypes, turbinePolyOrderMax) ;
+   float turbinePolyOmega(GADNumTurbineTypes, turbinePolyOrderMax) ;
+   float turbinePolyCl(GADNumTurbineTypes, turbinePolyClCdrNormSegments, turbinePolyOrderMax) ;
+   float turbinePolyCd(GADNumTurbineTypes, turbinePolyClCdrNormSegments, turbinePolyOrderMax) ;
+   int turbinePolyTwistOrder(GADNumTurbineTypes) ;
+   int turbinePolyChordOrder(GADNumTurbineTypes) ;
+   int turbinePolyPitchOrder(GADNumTurbineTypes) ;
+   int turbinePolyOmegaOrder(GADNumTurbineTypes) ;
+   int turbinePolyClOrder(GADNumTurbineTypes) ;
+   int turbinePolyCdOrder(GADNumTurbineTypes) ;
+
+The GAD turbine model capability has been implemented into FastEddy as an extension module, and is not compiled by default. The user needs to build FastEddy using the following compile flag below in order to include the GAD module:
+
+.. code-block:: none
+
+   make WITH_GAD=1
+
+The figure shows instantaneous contours of hub height (90 m) wind speed (in m/s) spatial distribution at three different times, showcasing the yawing of the turbines to align with the time-varying wind direction throughout the course of the simulation as it shifts from SW to NW. The gray areas represent the location over which GAD forces are applied (larger than the actual rotor area).
+
+.. image:: ../images/GAD_tutorial.png
+  :width: 1200
+  :alt: Alternative text
+
+.. note::
+
+   * The orientation of the turbine (:code:`GAD_rotorTheta`) is defined as the angle from the negative x axis, increasing counterclockwise. This is different from the meteorological convention for wind direction. For example, a turbine facing the west will have :code:`GAD_rotorTheta` = :math:`0.0^{\circ}`, while a turbine facing south will have :code:`GAD_rotorTheta` = :math:`90.0^{\circ}`.
+   * Application of the GAD to a real world WRF-coupled simulation does not require any additional steps besides the ones described here.
+
+Full citation references can be found in the :doc:`Publications <../../publications>` section.
diff --git a/docs/Tutorials/cases_extensions/WRF_coupling_case0_FE_urban.rst b/docs/Tutorials/cases_extensions/WRF_coupling_case0_FE_urban.rst
new file mode 100644
index 0000000..5934388
--- /dev/null
+++ b/docs/Tutorials/cases_extensions/WRF_coupling_case0_FE_urban.rst
@@ -0,0 +1,36 @@
+========================================================
+Real-world downscaled FastEddy simulation with buildings
+========================================================
+
+This tutorial involves setting up a real-world downscaled simulation that includes resolved buildings and closely follows the same procedure outlined in section 4.1. Additional steps for including buildings are described below and all required input datasets to run this tutorial are provided in this Zenodo record [TO BE UPDATED].
+
+In the *GeoSpec* preprocessing step, and additional 2d field describing building heights above ground level is required in the georeference input NetCDF file.
+
+.. code-block:: none
+
+   float BuildingHeights(y, x) ;
+
+This tutorial provides an example georeference input file with building heights for downtown Dallas, TX (:code:`Dallas_input_Oct2025_lod13.nc`). The *geospec.json* parameter file option :code:`urban_opt : 1` needs to be selected for building height information to be ingested in the reference standard-format NetCDF output file upon execution of **GeoSpec.py**.
+
+The same :code:`urban_opt : 1` option needs must be included in the subsequent *SimGrid* and *GenICBCs* stages, where building information is used for the creation of a :code:`BuildingMask` array containing gridded information of building presence, and ensuring that winds, subgrid-scale TKE and hydrometeors are set to zero within buildings for both initial and boundary conditions.
+
+After initial and boundary conditions have been created, a building-resolving FastEddy simulation can be undertaken by activating the urban model capability in the parameters file (see **tutorials/examples/Example10_REALCASE_Dallas_urban.in**).
+
+.. code-block:: none
+
+   #--URBAN
+   urbanSelector = 1 # urban selector: 0=off, 1=on
+
+The urban model capability has been implemented into FastEddy as an extension module, and is not compiled by default. To include the URBAN module in a build of FastEddy use the following compile flag:
+
+.. code-block:: none
+
+   make WITH_URBAN=1
+
+The model used to represent buildings follows the immersed body force approach described in *Muñoz-Esparza et al., 2020*, and the tutorial case corresponds to the passage of a cold front (*Muñoz-Esparza et al. (2021, 2025)*. The figure below shows instantaneous wind speed and vertical velocity fields corresponding to 30min hindcast valid at 1500 UTC on November 11th 2011 (pre-frontal conditions). These horizontal contours are from the model's third vertical level, located at approximately 23 m above ground level.
+
+.. image:: ../images/URBAN_tutorial_nz2_2panel.png
+  :width: 900
+  :alt: Alternative text
+
+Full citation references can be found in the :doc:`Publications <../../publications>` section.
diff --git a/docs/Tutorials/cases_real.rst b/docs/Tutorials/cases_real.rst
index d07c8a2..2c2499e 100644
--- a/docs/Tutorials/cases_real.rst
+++ b/docs/Tutorials/cases_real.rst
@@ -1,6 +1,6 @@
-***************
-Real Test Cases
-***************
+**********
+Real Cases
+**********
 
 Real cases are performed by dynamically downscaling from a mesoscale model such as WRF that can provide initial and boundary conditions for a one-way nested FastEddy simulation. These mesoscale-LES coupled simulations require the following preprocessing steps:
 
@@ -8,7 +8,7 @@ Real cases are performed by dynamically downscaling from a mesoscale model such
 * Step 2: **SimGrid**. Simulation grid definition step. Defines a FastEddy gridded domain at a specificed grid spacing, location and extent using the file resulting from ther previous GeoSpec step and a FastEddy input parameters file (with targeted domain configuration parameters) as inputs.
 * Step 3: **GenICBCs**. Generate initial conditions/boundary (ICBCs) conditions step. Creates ICBCs for a targeted FastEddy domain (defined in the SimGrid step) from a set of mesoscale model results. 
 
-The following tutorial provides a practical example of performing these preprocessing steps followed by a corresponding weather-driven FastEddy simulation for a real-world downscaled scenario.
+The following tutorial provides a practical example of performing these preprocessing steps followed by a corresponding weather-driven FastEddy simulation for a real-world downscaled scenario. 
 
 .. toctree::
 
diff --git a/docs/Tutorials/cases_real/WRF_coupling_case0_FE.rst b/docs/Tutorials/cases_real/WRF_coupling_case0_FE.rst
index ce0ba8b..fc7eb0f 100644
--- a/docs/Tutorials/cases_real/WRF_coupling_case0_FE.rst
+++ b/docs/Tutorials/cases_real/WRF_coupling_case0_FE.rst
@@ -1,6 +1,6 @@
-========================================================
-Running a real-world downscaled simulation with FastEddy
-========================================================
+===================================================
+Running a real-world downscaled FastEddy simulation
+===================================================
 
 After initial and boundary conditions have been created, a FastEddy simulation can undertaken with intial and boundary forcing from the mesoscale prognostic state fields. The lines below correspond to additions and modifications to the FastEddy parameters file necessary for coupled mesoscale-LES runs (corresponding to the test case from **tutorials/examples/Example08_REALCASE_FortCollins.in**)
 
diff --git a/docs/Tutorials/images/GAD_tutorial.png b/docs/Tutorials/images/GAD_tutorial.png
new file mode 100644
index 0000000..beea58b
Binary files /dev/null and b/docs/Tutorials/images/GAD_tutorial.png differ
diff --git a/docs/Tutorials/images/URBAN_tutorial_nz2_2panel.png b/docs/Tutorials/images/URBAN_tutorial_nz2_2panel.png
new file mode 100644
index 0000000..f216b3a
Binary files /dev/null and b/docs/Tutorials/images/URBAN_tutorial_nz2_2panel.png differ
diff --git a/docs/Tutorials/index.rst b/docs/Tutorials/index.rst
index fc9a5dd..ecf6283 100644
--- a/docs/Tutorials/index.rst
+++ b/docs/Tutorials/index.rst
@@ -8,7 +8,7 @@ Tutorials
       
    getting_started
    cases
-   sensitivity
    cases_real
+   cases_extensions
 
    
diff --git a/docs/Tutorials/sensitivity.rst b/docs/Tutorials/sensitivity.rst
index 74a9416..71c647a 100644
--- a/docs/Tutorials/sensitivity.rst
+++ b/docs/Tutorials/sensitivity.rst
@@ -1,3 +1,5 @@
+:orphan:
+
 *************************
 Sensitivity Tests (Ideal)
 *************************
diff --git a/docs/csv/extensions_gad.csv b/docs/csv/extensions_gad.csv
index e0d80fd..8f020fa 100644
--- a/docs/csv/extensions_gad.csv
+++ b/docs/csv/extensions_gad.csv
@@ -2,10 +2,10 @@
 **turbineSpecsFile:** |br| The path+filename to a turbine specifications file |br| *C-O if GADSelector > 0*, N/A, N/A, , C-O 
 **GADoutputForces:** |br| Flag to include GAD forces in the output: |br| 0=off; 1=on |br| *C-O if GADSelector > 0*, 0, 1, 0, C-O
 **GADofflineForces:** |br| Flag to compute GAD forces in an offline mode: |br| 0=off; 1=on |br| *C-O if GADSelector > 0*, 0, 1, 0, C-O
-**GADaxialInduction:** |br| Flag to compute axial induction factor: |br| 0=off (uses prescribed GADaxialIndVal) |br| 1=on |br| *C-O if GADSelector > 0*, 0, 1, 0, C-O
-**GADaxialIndVal:** |br| Prescribed constant axial induction factor when GADaxialInduction==1 |br| *C-O if GADSelector > 0*, 0.0, 1.0, 0.02, C-O
+**GADaxialInduction:** |br| Flag to compute axial induction factor: |br| 0=off (uses prescribed GADaxialIndVal) |br| 1=on |br| *C-O if GADSelector > 0*, 0, 1, 1, C-O
+**GADaxialIndVal:** |br| Prescribed constant axial induction factor when GADaxialInduction==0 |br| *C-O if GADSelector > 0*, 0.0, 1.0, 0.02, C-O
 **GADrefSwitch:** |br| Switch to use reference windspeed: |br| 0=off; 1=on |br| *C-O if GADSelector > 0*, 0, 1, 0, C-O
 **GADrefU:** |br| Prescribed constant reference hub-height windspeed |br| *C-M if GADSelector > 0 and* |br| *GADrefSwitch == 1*, 0.0, 50.0, 0.0, C-M
-**GADrefSampleWindow:** |br| Sample duration over which to average per-timestep values (filtering out highest frequencies) |br| *C-O if GADSelector > 0*, 1.0, 60.0, 10.0, C-O
-**GADrefSeriesLength:** |br| Number of sampling windows over which to average again for reference velocity magnitude and direction |br| *C-O if GADSelector > 0*, 1, 360, 30, C-O
+**GADrefSampleWindow:** |br| Sample duration (in seconds) over which to average per-timestep values (filtering out highest frequencies) |br| *C-O if GADSelector > 0*, 1.0, 300.0, 30.0, C-O
+**GADrefSeriesLength:** |br| Number of sampling windows over which to average again for reference velocity magnitude and direction |br| *C-O if GADSelector > 0*, 1, 360, 10, C-O
 **GADForcingSwitch:** |br| Switch to use the GADrefU-based or local windspeed in computing GAD forces: |br| 0=local |br| 1=ref |br| *C-O if GADSelector > 0*, 0, 1, 0, C-O
diff --git a/tutorials/examples/Example01_NBL.in b/tutorials/examples/Example01_NBL.in
index f940fc5..7093380 100755
--- a/tutorials/examples/Example01_NBL.in
+++ b/tutorials/examples/Example01_NBL.in
@@ -79,7 +79,7 @@ surflayer_z0 = 0.1 # roughness length (momentum) when surflayerSelector > 0
 surflayer_z0t = 0.01 # roughness length (temperature) when surflayerSelector > 0
 surflayer_wth = 0.0 # kinematic sensible heat flux at the surface when surflayerSelector = 1
 surflayer_tr = -0.0 # temperature rate at the surface when surflayerSelector = 2 (>0 for warming; <0 for cooling)
-surflayer_wq = 0.0 # 0.000052 # sensible heat flux at the surface (kg/kg m s-1) when surflayerSelector = 1
+surflayer_wq = 0.0 # latent 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
diff --git a/tutorials/examples/Example02_CBL.in b/tutorials/examples/Example02_CBL.in
index ab33a63..7600741 100755
--- a/tutorials/examples/Example02_CBL.in
+++ b/tutorials/examples/Example02_CBL.in
@@ -79,7 +79,7 @@ surflayer_z0 = 0.05 # roughness length (momentum) when surflayerSelector > 0
 surflayer_z0t = 0.05 # roughness length (temperature) when surflayerSelector > 0
 surflayer_wth = 0.35 # kinematic sensible heat flux at the surface when surflayerSelector = 1
 surflayer_tr = -0.0 # temperature rate at the surface when surflayerSelector = 2 (>0 for warming; <0 for cooling)
-surflayer_wq = 0.0 # 0.000052 # sensible heat flux at the surface (kg/kg m s-1) when surflayerSelector = 1
+surflayer_wq = 0.0 # latent 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
diff --git a/tutorials/examples/Example03_SBL.in b/tutorials/examples/Example03_SBL.in
index 92572a4..710485e 100755
--- a/tutorials/examples/Example03_SBL.in
+++ b/tutorials/examples/Example03_SBL.in
@@ -79,7 +79,7 @@ surflayer_z0 = 0.1 # roughness length (momentum) when surflayerSelector > 0
 surflayer_z0t = 0.1 # roughness length (temperature) when surflayerSelector > 0
 surflayer_wth = 0.0 # kinematic sensible heat flux at the surface when surflayerSelector = 1
 surflayer_tr = -0.25 # temperature rate at the surface when surflayerSelector = 2 (>0 for warming; <0 for cooling)
-surflayer_wq = 0.0 # 0.000052 # sensible heat flux at the surface (kg/kg m s-1) when surflayerSelector = 1
+surflayer_wq = 0.0 # latent 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
diff --git a/tutorials/examples/Example04_BOMEX.in b/tutorials/examples/Example04_BOMEX.in
index 8ecab14..3282cf8 100755
--- a/tutorials/examples/Example04_BOMEX.in
+++ b/tutorials/examples/Example04_BOMEX.in
@@ -79,7 +79,7 @@ surflayer_z0 = 0.0002 # roughness length (momentum) when surflayerSelector > 0
 surflayer_z0t = 0.0002 # roughness length (temperature) when surflayerSelector > 0
 surflayer_wth = 0.008 # kinematic sensible heat flux at the surface when surflayerSelector = 1
 surflayer_tr = -0.0 # temperature rate at the surface when surflayerSelector = 2 (>0 for warming; <0 for cooling)
-surflayer_wq = 0.052 # 0.000052 # sensible heat flux at the surface (kg/kg m s-1) when surflayerSelector = 1
+surflayer_wq = 0.052 # latent 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
diff --git a/tutorials/examples/Example05_CANOPY.in b/tutorials/examples/Example05_CANOPY.in
index 95f4226..d655d75 100755
--- a/tutorials/examples/Example05_CANOPY.in
+++ b/tutorials/examples/Example05_CANOPY.in
@@ -85,7 +85,7 @@ surflayer_z0t = 0.01 # 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.0 # kinematic sensible heat flux at the surface when surflayerSelector = 1
 surflayer_tr = -0.0 # temperature rate at the surface when surflayerSelector = 2 (>0 for warming; <0 for cooling)
-surflayer_wq = 0.0 # 0.000052 # sensible heat flux at the surface (kg/kg m s-1) when surflayerSelector = 1
+surflayer_wq = 0.0 # latent 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
diff --git a/tutorials/examples/Example06_OFFSHORE.in b/tutorials/examples/Example06_OFFSHORE.in
index d9d2f69..037395e 100755
--- a/tutorials/examples/Example06_OFFSHORE.in
+++ b/tutorials/examples/Example06_OFFSHORE.in
@@ -79,7 +79,7 @@ surflayer_z0 = 0.1 # roughness length (momentum) when surflayerSelector > 0
 surflayer_z0t = 0.01 # roughness length (temperature) when surflayerSelector > 0
 surflayer_wth = 0.0 # kinematic sensible heat flux at the surface when surflayerSelector = 1
 surflayer_tr = 0.5 # 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_wq = 0.0 # latent 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
diff --git a/tutorials/examples/Example07_DISPERSION_CBL.in b/tutorials/examples/Example07_DISPERSION_CBL.in
index 067c40a..6b6d955 100755
--- a/tutorials/examples/Example07_DISPERSION_CBL.in
+++ b/tutorials/examples/Example07_DISPERSION_CBL.in
@@ -92,7 +92,7 @@ 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_wq = 0.0 # latent 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
diff --git a/tutorials/examples/Example07_DISPERSION_SBL.in b/tutorials/examples/Example07_DISPERSION_SBL.in
index 1e977e9..99b1111 100755
--- a/tutorials/examples/Example07_DISPERSION_SBL.in
+++ b/tutorials/examples/Example07_DISPERSION_SBL.in
@@ -92,7 +92,7 @@ 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_wq = 0.0 # latent 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
diff --git a/tutorials/examples/Example08_REALCASE_FortCollins.in b/tutorials/examples/Example08_REALCASE_FortCollins.in
index 08d1af8..f37a641 100644
--- a/tutorials/examples/Example08_REALCASE_FortCollins.in
+++ b/tutorials/examples/Example08_REALCASE_FortCollins.in
@@ -97,7 +97,7 @@ surflayer_z0t = 0.1 # 0.01 # roughness length (temperature) when surflayerSelect
 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_wq = 0.0 # latent 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
diff --git a/tutorials/examples/Example09_GAD.in b/tutorials/examples/Example09_GAD.in
new file mode 100755
index 0000000..66530af
--- /dev/null
+++ b/tutorials/examples/Example09_GAD.in
@@ -0,0 +1,150 @@
+Description = This is an idealized Neutral Boundary Layer V&V scenario.  
+#--MPI_AALES
+numProcsX = 16 # 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 = ./ICBC/ # Path where initial/restart file is read in from
+inFile = FE_interp.0 # name of the input file for coordinate system and initial or restart conditions
+outPath = ./output/ # ./output_gad_10s30l/ # ./output_gad_60s5l/ # Path where output files are to be written
+outFileBase = FE_NBL # Base name of the output file series as in (outFileBase).element-in-series
+frqOutput = 6000 # 125 # 6000 # 7500 # 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 = 800 # Number of discretised domain elements in the x (zonal) direction
+Ny = 802 # Number of discretised domain elements in the y (meridional) direction
+Nz = 26 # 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 = 7.5 # Computational domain fixed resolution in the 'i' direction
+d_eta = 7.5 # Computational domain fixed resolution in the 'j' direction
+d_zeta = 15.0 # Computational domain fixed resolution in the 'k' direction
+coordHorizHalos = 1 # switch to setup coordiante halos as periodic=1 or gradient-following=0
+topoFile = # A file containing topography (surface elevation in meters ASL)
+verticalDeformSwitch = 1 # switch to use vertical coordinate deformation 0=off, 1=on
+verticalDeformFactor = 0.5 # 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 = 180000 # Number of timesteps to perform
+dt = 0.02 # timestep resolution in seconds
+NtBatch = 3000 # Number of timesteps to compute in batch launch, must have NtBatch <= Nt
+#--HYDRO_CORE
+##------------: HYDRO_CORE Submodule Selectors ---
+#----------: Boundary Conditions Set ---
+hydroBCs = 1 # Selector for hydro BC set. 2= periodicHorizVerticalAbl
+hydroBndysFileBase = ./ICBC/FE_Bndys
+hydroBndysFileStart = 0
+hydroBndysFileEnd = 60
+dtBdyPlaneBCs = 60.0
+ceilingAdvectionBC = 1 # 0 is a rigid lid
+##----------: HYDRO_IO/LOGGGING ---
+hydroForcingWrite = 0 # Switch for dumping hydroFldsFrhs for prognositic fields. 0 = off, 1=on
+hydroSubGridWrite = 1 # 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 = 40.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
+##------------: DIFFUSION ---
+diffusionSelector = 0 # diffusivity selector: 0= none, 1= const.
+nu_0 = 0.0 # constant diffusivity used when diffusionSelector = 1
+##------------: EXPLICIT FILTERS ---
+filterSelector = 1 # explicit filter selector: 0=off, 1=on
+filter_6thdiff_vert = 1 # vertical 6th-order filter on w selector: 0=off, 1=on
+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 = 100.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.01 # roughness length (temperature) when surflayerSelector > 0
+surflayer_wth = 0.0 # kinematic sensible heat flux at the surface when surflayerSelector = 1
+surflayer_tr = -0.0 # temperature rate at the surface when surflayerSelector = 2 (>0 for warming; <0 for cooling)
+surflayer_wq = 0.0 # latent 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
+#----------: CELL PERTURBATION METHOD --- 
+cellpertSelector = 1 # CP method selector: 0= off, 1= on
+cellpert_sw2b = 0 # only south and west bnds perturbed (for idealized cases): 0= off, 1= south & west, 2 = south 
+cellpert_amp = 1.0 # 0.5 # 0.1 # maximum amplitude for the potential temperature perturbations
+cellpert_nts = 750 # number of time steps after which perturbations are seeded             
+cellpert_gppc = 8 # number of grid points conforming the cell        
+cellpert_ndbc = 3 # number of cells normal to domain lateral boundaries   
+cellpert_kbottom = 1 # z-grid point where the perturbations start               
+cellpert_ktop = 30 # 20 # 30 # z-grid point where the perturbations end      
+cellpert_tvcp = 1 # time-varying CP method selector: 0= off, 1= on (hydroBCs == 1)    
+cellpert_eckert = 0.05 # 0.1 # 0.05 #0.1 # 0.2 # Eckert number for the potential temperature perturbations (hydroBCs == 1) 
+cellpert_tsfact = 0.25 # 0.5 # 1.0 # factor on the refreshing perturbation time scale (hydroBCs == 5)
+##------------ GAD  --------------------  
+GADSelector = 1 # Generalized Actuator Disk Selector: 0=off, 1=on
+turbineSpecsFile = ./GAD_NREL28_9WTs_tutorial.nc # The path+filename to a turbine specifications file
+GADoutputForces = 1 # Flag to include GAD forces in the output: 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.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 = 10.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
diff --git a/tutorials/examples/Example10_REALCASE_Dallas_urban.in b/tutorials/examples/Example10_REALCASE_Dallas_urban.in
new file mode 100644
index 0000000..b3991ce
--- /dev/null
+++ b/tutorials/examples/Example10_REALCASE_Dallas_urban.in
@@ -0,0 +1,162 @@
+Description = This is an example of multiscale WRF-FE over dowtown Dallas (TX) including resolved buildings.
+#--MPI_AALES
+numProcsX = 16 # 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 = ./ICBC/ # Path where initial/restart file is read in from
+inFile = FE_interp_143000UTC.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_Dallas # 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 = 896 # Number of discretised domain elements in the x (zonal) direction
+Ny = 898 # Number of discretised domain elements in the y (meridional) direction
+Nz = 58 # 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 = 5.0 # Computational domain fixed resolution in the 'i' direction
+d_eta = 5.0 # Computational domain fixed resolution in the 'j' direction
+d_zeta = 67.0 # Computational domain fixed resolution in the 'k' direction
+coordHorizHalos = 1 # switch to setup coordiante halos as periodic=1 or gradient-following=0
+topoFile = ./downtownDallas_urban_Topography_896x898.dat # A file containing topography (surface elevation in meters ASL)
+verticalDeformSwitch = 1 # switch to use vertical coordinate deformation 0=off, 1=on
+verticalDeformFactor = 0.14 # 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 = 180000 # 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 =  1 # Selector for hydro BC set. 2= periodicHorizVerticalAbl
+ceilingAdvectionBC = 1 #selector for normal advection at the domain ceiling (0=zero-valued w-BC, 1=non-zero w-BC)
+hydroBndysFileBase = ./ICBC/FE_Bndys # ./ICBC/FE_Bndys #
+hydroBndysFileStart = 0 #
+hydroBndysFileEnd = 16
+dtBdyPlaneBCs = 300.0
+##----------: 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 = 1 # 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 = 32.78 # 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 = 0 # 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 = 0 # Rayleigh damping layer selector: 0= off, 1= on.
+dampingLayerDepth = 600.0 # Rayleigh damping layer depth in meters
+##----------: SURFACE LAYER ---
+surflayerSelector = 3 # 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 # 0.01 # 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 # latent 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 = 1 # offshore selector: 0=off, 1=on
+surflayer_offshore_opt = 4 # offshore roughness parameterization: ==0 (Charnock), ==1 (Charnock with variable alpha), ==2 (Taylor & Yelland), ==3 (Donelan), ==4 (Drennan), ==5 (Porchetta)
+#----------: CELL PERTURBATION METHOD ---
+cellpertSelector = 1 # CP method selector: 0= off, 1= on
+cellpert_sw2b = 0 # 2 #1 # only south and west bnds perturbed (for idealized cases): 0= off, 1= south & west, 2 = south
+cellpert_amp = 2.5 # 0.1 # maximum amplitude for the potential temperature perturbations
+cellpert_nts = 250 # 100 # number of time steps after which perturbations are seeded
+cellpert_gppc = 8 # number of grid points conforming the cell
+cellpert_ndbc = 3 # number of cells normal to domain lateral boundaries
+cellpert_kbottom = 1 # z-grid point where the perturbations start
+cellpert_ktop = 40 # z-grid point where the perturbations end
+cellpert_tvcp = 1 # time-varying CP method selector: 0= off, 1= on (hydroBCs == 1)
+cellpert_eckert = 0.05 # Eckert number for the potential temperature perturbations (hydroBCs == 1)
+cellpert_tsfact = 0.333 #0.5 # 1.0 # factor on the refreshing perturbation time scale (hydroBCs == 1)
+#--URBAN
+urbanSelector = 1 # urban 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