Update the MPAS-A Kain-Fritsch CPS to WRF 4.4.1

[jherwehe]

This EPA_KF PR updates the MPAS-A release version of the KF convection
parameterization scheme to a kfeta module based on the latest version
included with WRF 4.4.1, thereby adding a convective trigger option,
subgrid-scale cloud feedbacks (utilizing cloud fraction, and water and
ice condensates) to RRTMG SW and LW radiation (Alapaty et al., 2012;
Herwehe et al., 2014), and a dynamic convective time scale (Bullock et
al., 2015) to MPAS-A. Other KF enhancements include extending the KF
look-up table (to prevent "OUT OF BOUNDS" errors being written to a
fort.98 file) and improving the DQEVAP and QICE calculations. This
version of module_cu_kfeta.F has also been generalized to work with
either MPAS or WRF.

New namelist.atmosphere options (with default values) added under &physics:
   config_kfeta_trigger = 1
   config_cu_rad_feedback = false

New variables available for output:
   cldfrac_dp ; sugbgrid-scale cloud fraction from KF deep convection
   cldfrac_sh ; sugbgrid-scale cloud fraction from KF shallow convection
   cldfracwcu ; horizontal cloud fraction, resolved + KF Cu
   cldfracwcut ; total cloud fraction, "vert-integrated" resolved + KF Cu

NOTE: These EPA_KF code changes to MPAS-A are based on the 22 August 2022
"develop" branch of MPAS v7.3.

References:
Alapaty, K., J. A. Herwehe, T. L. Otte, C. G. Nolte, O. R. Bullock, M.
   S. Mallard, J. S. Kain, and J. Dudhia, 2012: Introducing subgrid-
   scale cloud feedbacks to radiation for regional meteorological and
   climate modeling. Geophysical Research Letters, 39, L24808.
   https://doi.org/10.1029/2012GL054031
Bullock, O. R., K. Alapaty, J. A. Herwehe, and J. S. Kain, 2015: A
   dynamically computed convective time scale for the Kain–Fritsch
   convective parameterization scheme. Monthly Weather Review, 143,
   2105-2120. https://doi.org/10.1175/MWR-D-14-00251.1
Herwehe, J. A., K. Alapaty, T. L. Spero, and C. G. Nolte, 2014:
   Increasing the credibility of regional climate simulations by
   introducing subgrid-scale cloud-radiation interactions. Journal of
   Geophysical Research: Atmospheres, 119, 5317-5330.
   https://doi.org/10.1002/2014JD021504

[ldfowler58]

Here a few comments that I have regarding the PR:

• In Registry.xml:

1. what about adding snowncv and graupelncv in the output stream since you added the two arrays in the
   restart streams (I am actually adding them as well in some of my updates)?
2. rename config_cu_rad_feedback to config_cu_radt_feedback.
3. cldfracwcu: description - change "including Cu" with "including convective cloud fraction".
4. cldfracwcut: description - change "including Cu" with "including convective cloud fraction".
5. cldfract: description - add using "random overlap assumption".

• In mpas_atmphys_driver_radiation_lw.F and mpas_atmphys_driver_radiation_sw.F:

1. in subroutine radiation_lw_from_MPAS, change cu_ra_feedback to cu_radt_feedback.
2. change call mpas_pool_get_array(diag_physics,'cldfrac' ,cldfrac ) with
   call mpas_pool_get_array(diag_physics,'cldfrac',cldfrac) in two places.

• In mpas_atmphys_driver_cloudiness.F:
  use cu_radt_feedback instead of cu_rad_feedback.
  in subroutine cloudiness_to_MPAS, rephrase "!compute cloud diagnosis (random overlapping) by ZCX (from WRF)".
  We do not know who ZCX is anyway. No need to reference WRF.
  removed definition 351-356 (already defined twice).

• In mpas_atmphys_vars.F:
  add unit to cldfracrad_p.

• In module_cu_kfeta.F:

1. have you actually tested the option trigger = 2 and = 3? I do not think that it will work? I think that we should abort the run when those trigger options are used? What do you think?

[jherwehe]

Reply for PR #994 comments from ldfowler58 on 11 April 2024:

• In Registry.xml:

1. Probably do not need to add snowncv and graupelncv to the output stream since snownc and graupelnc are not currently in the output stream; snowncv and graupelncv were added to the restart stream for completeness (and as a precaution) since snownc and graupelnc were already there.
2. If "radt" is being used as the abbreviation for "radiation", then changing the option to config_cu_radt_feedback is okay; we had retained the "cu_rad_feedback" naming for consistency with the same option in WRF.
3. Yes; "including convective cloud fraction" is a better description.
4. Yes; "including convective cloud fraction" is a better description.
5. Yes; adding "using random overlap assumption" is a better description.

• In mpas_atmphys_driver_radiation_lw.F and mpas_atmphys_driver_radiation_sw.F:

1. Yes, changing to "cu_radt_feedback" everywhere is necessary, as well as changing to "config_cu_radt_feedback", in mpas_atmphys_driver_radiation_lw.F and mpas_atmphys_driver_radiation_sw.F (see 2. above).
2. Yes, if only extra spaces are being omitted.

• In mpas_atmphys_driver_cloudiness.F:

1. Yes, please use "cu_radt_feedback" and "config_cu_radt_feedback" for consistency with the above points.
2. Agreed; the ZCX and WRF comment can be removed.
3. Not sure what repeated definition at 351-356 was removed; Can you elaborate by showing the code or indicating where I can find the change in GitHub?

• In mpas_atmphys_vars.F:

1. Yes, adding a unit of [-] to cldfracrad_p is good.

• In module_cu_kfeta.F:

1. Yes; though we usually use config_kfeta_trigger=1, we have successfully tested config_kfeta_trigger=2 and config_kfeta_trigger=3 using our EPA-enhanced version of the MPAS v7.3 develop branch code; the resulting accumulated precipitation fields for a July 2016 simulation on the 92-25km mesh were compared against 4-km PRISM data (CONUS) and 0.1-degree IMERG data (global) revealing that trigger3 is similar to trigger1, and trigger2 is too dry (global and CONUS precipitation plots are available and can be provided upon request).
   Triggers 2 and 3 do seem to work, but code could be added to abort the MPAS run if those triggers are selected; have those additional triggers been useful in WRF? If so, maybe their functionality should be retained.

[ldfowler58]

Can you please resolve the conflicts in mpas_atmphys_driver.F, mpas_atmphys_driver_radiation_lw.F, and mpas_atmphys_driver_radiation_sw.F. I think that you need to resolve those conflicts first so that I can see the updated files with a new hash. Thanks.

[jherwehe]

The code conflicts in these modules have been resolved as of 25 July 2024.

[ldfowler58]

Thanks for merging the sourcecode with v8.2.0.

Makefile: Michael will need to review changes made to Makefile (see lines 363-367 for intel-mpi and line 782). you need to replace line 979 with that in develop.

./src/core_atmosphere/Registry.xml:

1. please move 2175-2185 below line 2228 for consistency with the organization of the config options. In the develop branch, the config options between lines 2161 and 2204 list the names of the different parameterizations whereas options for the individual parameterizations start at line 2206. To preserve the original organization, I suggest to move lines 2175-2185 below the possible values for config_gfconv_closure_shallow. Thanks.

2. As suggested earlier:
   -> rename config_cu_rad_feedback to config_cu_radt_feedback.
   -> cldfracwcu: description - change "including Cu" with "including convective cloud fraction".
   -> cldfracwcut: description - change "including Cu" with "including convective cloud fraction".
   -> cldfract: description - add using "randowm overlap assumption".

[ldfowler58]

Jerry: I did not see that you modified the sourcecode in v8.2.0 as I suggested earlier. Hopefully, I am not overlooking something.

• mpas_atmphys_driver.F:
  -> please remove lines 102-103 since the actual sourcecode in subroutine physics_driver is not modified relative to the develop branch.

• mpas_atmphys_vars.F:
  line 764: add unit to cldfracrad_p.

• mpas_atmphys_driver_radiation_lw.F and mpas_atmphys_driver_radiation_sw.F:

1. in subroutine radiation_lw_from_MPAS (radiation_sw_from_MPAS), change cu_ra_feedback to cu_radt_feedback.
2. change call mpas_pool_get_array(diag_physics,'cldfrac' ,cldfrac ) with
   call mpas_pool_get_array(diag_physics,'cldfrac',cldfrac) in two places.

• mpas_atmphys_driver_cloudiness.F:
  use cu_radt_feedback instead of cu_rad_feedback.
  in subroutine cloudiness_to_MPAS, rephrase "!compute cloud diagnosis (random overlapping) by ZCX (from WRF)".
  We do not know who ZCX is anyway. No need to reference WRF.
  removed definition 351-356 (already defined twice).

[jherwehe]

Laura, I may have misunderstood my responsibilities in the process of preparing our PRs for merging. I am not very experienced with Git or GitHub, so I may not be working this process correctly according to NCAR's practices. Please note the following points:

• The changes for yesterday's code conflict resolution commit (3d9de31) went to into my original branch for this PR (Update the MPAS-A Kain-Fritsch CPS to WRF 4.4.1 #994) called atmosphere/develop-EPA_KF, which is based on MPAS v7.3 develop branch code.
• My master branch at https://github.com/jherwehe/MPAS-Model/tree/master is just a pure clone or fork of NCAR's release code for MPAS v8.2.0, unmodified. All of our EPA-modified codes start from the release code clone/fork and are developed, tested, and reside on EPA's HPC system outside of GitHub. The only time I have pushed our EPA-modified codes to GitHub has been for the separate branches for each of our PRs (Update the MPAS-A Kain-Fritsch CPS to WRF 4.4.1 #994, Add 3-D grid analysis nudging FDDA to MPAS-A #995, Add the ACM2 PBL scheme to MPAS-A #1002, Add the Pleim surface layer scheme to MPAS-A #1003, Add options to core_init_atmosphere for P-X LSM #1010, Add the Pleim-Xiu (P-X) LSM scheme to MPAS-A #1025, and Add the Pleim-Xiu (P-X) LSM and FDDA to MPAS-A #1027), and all of these modifications were originally based on the MPAS v7.3 develop branch.
• In my response to your original April 11 suggestions for this PR, I thought that these were changes that you were planning to make to the code for this PR in whatever development branch you were working in and were seeking my confirmation or feedback before proceeding. So on second thought, was I supposed to make those suggested changes to the code in my atmosphere/develop-EPA_KF branch? If so, sorry; that was not done. Guess I did not understand what my responsibilities were on this end. Please clarify.
• I have all of our EPA physics modifications working in MPAS v8.1.0, but have had problems getting our EPA-modified version of MPAS v8.2.0 to run more than a few hours of simulation. MPAS v8.2.0 seems to be set up to run in single precision (the new default), but for our air quality applications, we prefer to keep MPAS running in double precision (though it is about 25% slower; and I had to make a modification to the tnccn_act declaration in module_mp_thompson.F to even get the release code to compile in double precision). [On a side note: There is a problem (and this is something I should probably submit to the Forum if it isn't already noted there) introduced in MPAS v8.1.0 that does not allow the model to run with config_frac_seaice turned off; the fix is to include additional array allocations in mpas_atmphys_driver_seaice.F.] I am still working on stabilizing our EPA modifications to MPAS v8.2.0.

Please tell me how to proceed to meet NCAR's requirements and facilitate your work with these EPA PRs.

• Can you work with my MPAS v7.3-based branches (i.e., atmosphere/develop-EPA_*) for our PRs that I currently have in my GitHub if I resolve the current code conflicts? And would implementing your suggested code changes into the existing MPAS v7.3-based branches (and resolving any resulting code conflicts) be sufficient for you to then implement our codes into whatever branch of MPAS v8.2.0 that you are currently working with?
• Or do I need to implement your suggested code changes into MPAS v8.2.0 code, then retest, redo, and resubmit all 7 of our PRs as new branches based on MPAS v8.2.0 code now?

So to conclude, all of our EPA development work is conducted on our EPA HPC system, and code is only pushed back to my GitHub fork/clone as needed as new branches for PRs to the MPAS repo.

I'll wait to hear back from you before taking any more steps. Thank you for your patience with my GitHub inexperience.

[ldfowler58]

1. I rebased your original sourcecode (that using MPAS version 7.3) using the latest develop branch so that the modifications that you made now show at the top of the logfile, or:

commit a618d57770614909a12b6e423342982bdf15c78f
Author: Jerold Herwehe herwehe.jerry@epa.gov
Date: Thu Sep 29 15:09:37 2022 -0400

Update the MPAS-A Kain-Fritsch CPS to WRF 4.4.1

This EPA_KF commit will update the MPAS-A release version of the KF

.
Alapaty, K., J. A. Herwehe, T. L. Otte, C. G. Nolte, O. R. Bullock, M.
S. Mallard, J. S. Kain, and J. Dudhia, 2012: Introducing subgrid-
scale cloud feedbacks to radiation for regional meteorological and
climate modeling. Geophysical Research Letters, 39, L24808.
https://doi.org/10.1029/2012GL054031

commit 66d1ab8
Merge: 02e2dea b566fc8
Author: Michael Duda duda@ucar.edu
Date: Thu Aug 8 13:55:30 2024 -0600

Merge branch 'master' into develop

This merge introduces fixes from MPAS v8.2.1, and it connects the v8.2.1 tag to
future commits on 'develop'.

2. Using the intel compiler after adding the option "-fp-model precise", I compared differences in the restart file between your sourcecode and the develop branch at the end of a 3-hour forecast when running the convection-permitting suite. Differences result because of changes that were made in module_ra_rrtmg_sw.F, or:

< !Balwinder.Singh@pnnl.gov: Code added to avoid 'divide by zero' error in zwo computation
< denom = max((1._rb - (1._rb - zw) * (zg / (1._rb - zg))**2),1.0E-30_rb)

        denom = (1._rb - (1._rb - zw) * (zg / (1._rb - zg))**2)
        ! Code added to avoid 'divide by zero' error in zwo computation
        if (abs(denom).eq.0.0_rb) then
           denom = 1.0E-30_rb
        end if

Do you actually need those changes in the sourcecode? It seems to me that the formulation of Balwinder.Singh does take care of dividing by zero.

I would really welcome the input from Michael, of course, but I would prefer to keep module_ra_rrtmg_sw.F so that changes only relate to the changes made to the KF scheme and calculation of the cloud fraction when KF is used.

Thanks.

[jherwehe]

Regarding item 1 in the previous post: Thank you. Hope this helps the forward movement on this PR.

Regarding item 2: Yes, please go ahead and retain the Balwinder.Singh coding, which accomplishes the same as the changed code. That code change that you encountered in module_ra_rrtmg_sw.F was carried over from a supposed fix to an issue reported on the WRF Model Users website and introduced in the same module back in WRF-4.1.3. Later, apparently a newer coding to avoid divide-by-0 was introduced in module_ra_rrtmg_sw.F in WRF-4.2, but those changes did not get transferred to module_ra_rrtmg_sw.F in MPAS.

Not sure why you would see different results or restart files due to this source code change to module_ra_rrtmg_sw.F.

Aside: At the time of compiling the MPASv7.3-based code for this PR, a slightly modified version of the Makefile's intel-mpi section (which included "-fp-model precise" among other additional compile options) was used with the Intel 2021.4 compilers.