addition of HWRF Noah LSM and GFDL surface layer

[grantfirl]

This PR supercedes #431 and #441 (and combines changes for HWRF Noah and GFDL surface layer).

HWRF Noah changes:

• added "working" LSM routines module_sf_noahlsm.F90 and module_sf_noahlsm_glacial_only.F90
  • modified working routines to use CCPP error messaging
  • these files originated with the HWRF v4.0a codebase, obtained via https://dtcenter.org/community-code/hurricane-wrf-hwrf/hwrf-system-components-v4-0a
• added CCPP-compliant interfaces to the working subroutines in sfc_noah_wrfv4.F90/meta
  • note that this scheme is "pure" in that it only calls the working subroutines for each point in the horizontal "chunk" of grid points (it should be extremely portable)
• added CCPP-compliant interfaces for pre/post interstitial schemes that handle all data preparation and diagnostics as part of a GFS-based suite (sfc_noah_wrfv4_interstitial.F90/meta)
• remove inoperable set_source_file_properties statement in CMakeLists.txt
• edit GFS_surface_generic.F90 to handle more ivegsrc values
• edit sfc_drv.f, sfc_drv_ruc.F90, sfc_noahmp_drv.F to check for valid values of ivegsrc and isot
• note that this scheme does not utilize a set_soilveg routine or a namelist_soilveg module (it uses module variables in module_sf_noahlsm though)

for GFDL surface layer

• add module_sf_exchcoef.f90 for internal subroutines
• add gfdl_sfc_layer.F90/meta as CCPP-compliant interface to GFDL surface layer scheme
  • edited the working routine (MFLUX2) to operate on 1-element arrays and use CCPP error messaging
  • added checks for proper values of flags in initialization routine
  • scheme now looks for previously defined values of variables drysmc and maxsmc instead of hard-coding them (see TODO note below)
  • define local flags that turn on/off features that may not be used in suites (calculating quantities that are done in different schemes within a GFS-based suite)
  • added logic for working with fractional land masks (using dry, icy, wet variables and potentially calling the working routine up to 3 times for each grid point)
  • added checks on roughness lengths from GFS surface layer and example scheme from Chunxi Zhang
  • added bounds on output terms (like surface exchange coefficients, following other surface layer schemes)

TODO:

• The use of namelist_soilveg modules and set_soilveg subroutines is now a problem because more than one scheme needs access to that data. The module-passing of variables was originally allowed in the CCPP via this mechanism because it was deemed internal to a scheme. This is no longer the case. I think to fix this requires re-working this mechanism or at least exposing variables that are in namelist_soilveg to the CCPP and passing in data through argument lists. This was deemed beyond the scope of this work, so a workaround was written in to the GFDL surface layer scheme that grabs the drysmc and maxsmc variables depending on the value of the 'lsm' variable. This is clunky and not portable if more LSMs join the CCPP, for example.
• Many "features" in both of these schemes were not ported in this work (urban physics, FASDAS, and even restarts for the HWRF LSM; using the HWRF preprocessor flag, coupling to waves, using ocean currents, pertubations, or running with fhcyc /= 0 in the GFDL surface layer). If any of these features are needed, they can be addressed in followup PRs.

Regression testing and simple reality-checking plots are included in the ufs-weather-model PR.

[grantfirl]

Associated PRs:
#452
NCAR/fv3atm#49
NCAR/ufs-weather-model#47

[grantfirl]

The following corresponds to commits b850fe, 6d3ce4, b8629e, and 0298ce:

1. for the GFDL surface layer scheme to use precalculated values of the lowest model layer wind speed that includes a convective or gustiness component (like the GFS surface layer, instead of recalculating it without this component)
2. tighter bounds on surface exchange coefficients following Chunxi's advice and GFS-based schemes
3. maintaining the same ratio relationship between M-O stability functions for enthalpy at the lowest model layer and 2m (for diagnostic 2m purposes)
4. introducing a time-averaged component to the calculation of surface skin temperature and soil temperatures in the HWRF Noah LSM (following the GFS version).

The first change saves a bit of computation and allows for better consistency amongst physics schemes, and primarily affects low wind speed areas over land (despite the change technically being "active" over all surface types). This might affect tropical cyclogenesis or at least African Easterly waves over land? This can likely be reverted, but surface layer schemes can produce funny results with low wind speeds (even with an enforced minimum of 1 m s-1) and adding the gustiness component is a somewhat artificial, yet physically-based way to boost surface wind speeds in quiescent but convectively active areas over land that has been widely used in many models/schemes.

The second change is primarily for computational stability. Low values of exchange coefficients lead to spurious values of several surface-related quantities since they are found in the denominator of their formulas. Although these bounds were taken from either Chunxi's suggestions of GFS-based schemes, I'm not totally confident that they are "correct". I'm sure this would benefit from sensitivity studies or input from surface layer experts.

The third change was necessary because bounding only the surface exchange coefficients but not their 2m counterparts can lead to unphysical values of 2m diagnostics. Maintaining the calculated ratio between, e.g. c_h and c_h_2m, from the output of the "working" part of the GFDL surface layer scheme when the surface exchange coefficient is bounded by change 2) above fixes this. I'm not 100% confident that this is the "correct" thing to do either, but it at least seems like a reasonable thing to do.

Even with changes 1-3, the HWRF Noah and GFDL surface layer schemes could still produce spuriously high values of T2m and q2m. I noticed in the GFS version of the Noah LSM that they introduced time-averaging of the calculation of surface skin temperature and soil temperatures (averaging old, input values and new, output values inside the LSM prior to returning them). Surface skin temperature is used to diagnose T2m and in the calculation of q2m when certain conditions are met. Adding this "feature" to the HWRF Noah LSM seemed to do the trick of eliminating spuriously large values of T2m and q2m.

Plots comparing the control GFS v15p2 suite with the commits listed above at the 12h mark (where the largest spurious values were) are attached.

HWRF_sfc_comp_plots.pptx

[climbfuji]

Thanks for adding these two schemes. We are ready to start the merge process. Please merge ccpp-physics. Then, in your local working directory, checkout the head of dtc/hwrf-develop in ccpp-physics, step up to fv3atm, revert the change to .gitmodules and update the submodule pointer for ccpp-physics.

[grantfirl]

Thanks for adding these two schemes. We are ready to start the merge process. Please merge ccpp-physics. Then, in your local working directory, checkout the head of dtc/hwrf-develop in ccpp-physics, step up to fv3atm, revert the change to .gitmodules and update the submodule pointer for ccpp-physics.

Finished.

[ChunxiZhang-NOAA]

Hi Grant,

The GFDL scheme layer scheme looks good. In terms of your recent changes:
 1) It is good to include the convective component w* in the total wind speed calculation.
 2) It may help to set bounds on surface exchange coefficients.
 
 My experience with my surface layer scheme is that some very large q2m grid points occur over large roughness length (z0). The z0s are usually larger than 2 meters, which are too high.  As a result, the integrated stability functions are very small and the exchange coefficients are very large.  You can do some tests with the upper limit of z0 to 0.5 m. 

Chunxi

[grantfirl]

Thanks @ChunxiZhang-NOAA for taking the time to look at the code (and to offer helpful suggestions in the past!).

I also noticed the large q2m values over large roughness lengths. I had done some tests with limiting z0 to 0.5m in the past, but that alone wasn't sufficient to prevent spurious values with this code.

I think that the plan is to go ahead with the code as-is for hurricane testing this summer, since our charge (as DTC) was to simply get the HWRF suite working in the CCPP. IMO, it could definitely benefit from some fine-tuning in all aspects should the CCPP-based HWRF suite to go into operations in the future. Although the values that "stood out" were 2m T and q, the GFS surface diagnostics scheme bases their calculation on their respective skin values, which can also be large, and those come from the LSM, which of course use surface layer exchange coefficients as input. Regardless, I think that if tuning is performed with the suite, it would be better performed by those more familiar with surface layer, LSM, ocean, and ice schemes, rather than DTC, who are positioning themselves as physics "generalists" and/or CCPP framework experts.

[ChunxiZhang-NOAA]

I agree with you. DTC's task is to make the HWRF suite work in the CCPP. The further tuning and improving is not an easy task.