Add a scale-adaptive 3DTKE parameterization scheme as a subgrid mixing option

[zhangx-3dtke]

TYPE: new feature

KEYWORDS: 3DTKE scheme, scale-adaptive

SOURCE: Xu Zhang (Shanghai Typhoon Institute/CMA, China)

DESCRIPTION OF CHANGES:
This PR adds a new scale-adaptive 3DTKE scheme (SMS-3DTKE, Shanghai Meteorological Service)
for subgrid turbulent mixing. The new scheme is self-adaptive to the grid-size between the
large-eddy simulation (LES) and mesoscale limits. The new 3DTKE scheme is added as a new
diffusion option km_opt = 5. This option extends original LES subgrid model 1.5-order TKE closure
(km_opt=2) to the mesoscale, and can be used in LES, mesoscale, and the gray zone resolutions in
between. In the horizontal diffusion, the new scheme blends horizontal Smagorinsky first-order
closure (km_opt=4) and 1.5-order TKE-based horizontal diffusivity. In the mesoscale limit,
horizontal diffusion recovers to the 2D Smagorinsky form (km_opt=4), and the new vertical
diffusion adds a non-local term and is made implicit to allow for longer time steps and thin levels.
The new scheme can therefore replace km_opt=2 and conventional PBL schemes. In the LES limit,
km_opt=5 recovers km_opt=2. The added 3D arrays are packaged in Registry.

Several newly added subroutines:

1. SUBROUTINE nonlocal_flux: prescribes the nonlocal heat flux profile based on LES analysis
   and computes the nonlocal momentum gamma term.
2. SUBROUTINE free_atmos_length: calculates mixing length in the free atmosphere.
3. SUBROUTINE meso_length_scale: calculates mesoscale length scale.
4. SUBROUTINE vertical_diffusion_implicit: solves the vertical diffusion equations for u, v, w, th,
   tke, moist, chem, scalar, tracer using an implicit method.
5. SUBROUTINE update_tke_implicit: solves the TKE equation using implicit method. The vertical
   TKE diffusion and dissipation are implicitly treated.

From a user perspective, when the new scheme is turned on (diff_opt=2 and km_opt=5), PBL
schemes must be turned off (bl_pbl_physics = 0).

This scheme can be used with sf_sfclay_physics = 1, 5, 91.

LIST OF MODIFIED FILES:
Registry/Registry.EM_COMMON
dyn_em/module_diffusion_em.F
dyn_em/module_first_rk_step_part2.F
dyn_em/solve_em.F
run/README.namelist
test/em_real/examples.namelist
share/module_check_a_mundo.F

TESTS CONDUCTED:

1. All WTF regression tests pass.
2. Restart run test passed.
3. The real case tests give reasonable results when km_opt = 5 comparable to the conventional PBL schemes for the following:
   • Summer and winter 24-hr case with the resolution of 20km.
   • 18-hr case with resolution of 3km.
4. check_a_mundo tests correctly identified incorrect accompanying physics selections.
5. It is verified that no extra memory is being used when the scheme is not used.
6. Timing tests were conducted for the 3km resolution runs.
   • When the new option is turned on (km_opt=5), it requires 1.95s every time step.
   • When km_opt=2, it needs 1.88s every time step.
   • When turn on the YSU PBL scheme, it requires 1.42s every time step.

RELEASE NOTE:
A three dimensional (3D) scale-adaptive TKE subgrid mixing parameterization scheme (SMS-3DTKE, km_opt=5) developed by Shanghai Meteorological Service (SMS), China. (Zhang, X., Bao, J., Chen, B., and Grell, E., 2018: A Three-Dimensional Scale-Adaptive Turbulent Kinetic Energy Scheme in the WRF-ARW Model. Mon. Wea. Rev., 146(7), 2023-2045, https://doi.org/10.1175/MWR-D-17-0356.1) . The subgrid mixing parameterization extends the original 3DTKE model (km_opt=2, Deardorff 1980) that is usually used as an LES subgrid model to the mesoscale limit in the framework of the WRF-ARW. The scheme can be used in LES, mesoscale and the gray zone resolutions in between. The option must be used with diff_opt=2. The scheme can replace LES subgrid model (km_opt=2) and conventional PBL schemes. When the scheme is turned on (diff_opt=2 and km_opt=5), PBL schemes must be turned off (bl_pbl_physics = 0). This scheme can be used with sf_sfclay_physics = 1, 5, 91.

[davegill]

@smileMchen

From a user perspective, when the new scheme is turned on (diff_opt=2 and km_opt=5), PBL schemes should be turned off (bl_pbl_physics = 0).

Ming,
With these stated restrictions, we should include some IF tests in check_a_mundo

This scheme is combined with sf_sfclay_physics = 1.

Does this mean that the scheme ONLY works with sfclayrev? If so, another check_a_mundo IF test.

[davegill]

@smileMchen
Ming,
Since this is a new scheme

1. we need a release note in the PR commit message
2. add an entry in run/README.namelist for the new option
3. put use information in test/em_real/examples.namelist

Do we need to update the TechNote or the User's Guide?

[smileMchen]

@davegill
Dave,
I am sorry to answer late but this is because I was busy with Xi Zhao this afternoon. I have talked to Xu Zhang and will take care of the check_a_mundo/readme/example.namelist issue. We will also update the User's Guide.

[weiwangncar]

@zhangx-3dtke Did you test restart? If you did, and it works, please add it to 'Tests Conducted'.

[weiwangncar]

@zhangx-3dtke Do you need this variable to be in history file?

[zhangx-3dtke]

No, I output this variable for diagnosis and forgot to change back.

[zhangx-3dtke]

Restart run test passed. I added it to "Tests Conducted" @weiwangncar

[smileMchen]

I looked through this PR and I believe it is ready to merge into the development branch.

[smileMchen]

@zhangx-3dtke
Let's wait for a while for other people's opinion

[davegill]

@zhangx-3dtke
Here are a couple more tests to conduct and add to the commit message. These are related specifically to your user recommendations.

From a user perspective, when the new scheme is turned on (diff_opt=2 and km_opt=5), PBL
schemes must be turned off (bl_pbl_physics = 0).

1. Would you activate the new scheme (diff_opt=2 and km_opt=5), and INCORRECTLY leave the PBL ON. Verify that your check_a_mundo tests catch this.

This scheme can be used with sf_sfclay_physics = 1, 5, 91.

2. Try your scheme (diff_opt=2 and km_opt=5) with the WRONG surface layer option. Again, verify that check_a_mundo picks this up correctly.

3. Turning your option ON and OFF, can you give an idea of the additional amount of time per time step this option requires? Just a few time steps are required, use the first dt after the initial step (no radiation, no I/O).

4. Verify that no extra memory is being used. Please do a before (your mods) vs after (your mods). Turn off the new scheme for the after-your-mods test. The amount of memory used (as reported towards the top of the WRF model), should be quite similar.

[davegill]

@zhangx-3dtke
I assume this was closed accidentally? If so, just click the "reopen and comment" button

[zhangx-3dtke]

@davegill
1, I conducted the check_a_mundo tests, and verified it can catch the wrong setting in the namelist.
2, I did 3km resolution runs. When the new option is turned on (km_opt=5), it requires 1.95s every time step. When km_opt=2, it needs 1.88s every time step. When turn on the YSU PBL scheme, it requires 1.42s every time step.
3, I verified that no extra memory is being used. The amount of memory used are quite similar, 198126484 bytes (before mods) , 198162676 bytes (after mods).
Some test information are added to the commit message.

[davegill]

@zhangx-3dtke
Line up this new code

[davegill]

@zhangx-3dtke
NO more comments about it, but go through and fix the indenting

[weiwangncar]

@davegill @dudhia I wonder what you think of comments like '!XZ' throughout the code. I remember we talked about things like this at one point. @zhangx-3dtke Is subroutine solve_tke_implicit functionally like rk_update_scalar for tke? If it is, could a different subroutine name be used (such as rk_update_tke_implicit)? Just curious.

[zhangx-3dtke]

@weiwangncar This is a good idea! But, actually solve_tke_implicit is not using the Runge-Kutta time integration. I can change the name to update_tke_implicit.

[liujake]

This one breaks WRFPLUS compilation. Debugging …

[davegill]

@liujake
Jake,
This PR also breaks all of the km_opt==2 option idealized cases. However, we added two PRs to the top of the develop branch to provide a temporary fix. Can you get an updated develop branch and try again?

[liujake]

I cloned repo this noon, it is already the latest code.

[liujake]

@davegill Most of time, wrfplus compilation failure is caused by newly added argument in forward code, which is sometimes also called in TL and/or AD code. Usually it is easy fix by simply adding new argument in TL/AD code. But this one seems to be a bit tricky. After I added
new argument u_h_tend and v_h_tend, I got error message requesting the declaration of these new variables. Still trying. Suggestion is welcome.

s/ch/opt/netcdf/4.6.3/gnu/8.3.0/include -fdefault-real-8 module_diffusion_em_ad.f90
module_diffusion_em_ad.f90:11051:63:

CALL horizontal_diffusion_u_2(ru_tendf,config_flags,u_h_tend,defor11,defor12,div,  &
                                                           1

Error: Symbol ‘u_h_tend’ at (1) has no IMPLICIT type; did you mean ‘a_rw_tendf’?
module_diffusion_em_ad.f90:11075:63:

CALL horizontal_diffusion_v_2(rv_tendf,config_flags,v_h_tend,defor12,defor22,div,  &
                                                           1

Error: Symbol ‘v_h_tend’ at (1) has no IMPLICIT type; did you mean ‘a_rw_tendf’?
module_diffusion_em_ad.f90:7521:26:

                       cf1, cf2, cf3, msftx, msfty,             &
                      1

Error: Rank mismatch in argument ‘fnp’ at (1) (rank-1 and scalar)
module_diffusion_em_ad.f90:7528:49:

                       tke, xkhv, BN2, theta, dt,               &
                                             1

Error: Rank mismatch in argument ‘theta’ at (1) (rank-3 and scalar)
module_diffusion_em_ad.f90:7539:66:

                       its, ite, jts, jte, kts, kte           )
                                                              1

Error: Type mismatch in argument ‘hpbl’ at (1); passed INTEGER(4) to REAL(8)
module_diffusion_em_ad.f90:3803:76:

                               its, ite, jts, jte, kts, kte             )
                                                                        1

Error: Type mismatch in argument ‘hpbl’ at (1); passed INTEGER(4) to REAL(8)

[liujake]

@davegill . Resolved u_h_tend and v_h_tend lines by defining them as local variables. Move to next lines. My strategy is to simply allow code to compile, but will not try to update TL/AD code to match new modified forward code.

[weiwangncar]

@liujake u_h_tend, v_h_tend and a few others are state variables declared in Registry. Maybe they should be added to TL/AD registry?