Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Radiation Improvements for prototype 8 #871

Merged
merged 41 commits into from
Apr 8, 2022

Conversation

dustinswales
Copy link
Collaborator

@dustinswales dustinswales commented Mar 8, 2022

This PR contains code changes to both the RRTMG and RRTMGP radiation schemes.

RRTMGP: Many of the changes are organizational, or added infrastructure for new features. This includes:

  • New single module for all of the RRTMGP/cloudMP coupling.
  • Updates to RRTMGP Thompson MP.
  • Option (default=F) for explicit sgs clouds coupling in GP (GF/SAMF convective clouds: MYNN-EDMF PBL cloud).

RRTMG: A request has been made by many physics developers and users to rewrite
the cloud computation routines (progcldXXX) in the program
radiation_clouds.f. Each microphysics scheme has one subroutine to
calculate the cloud radiation properties, and those similar subroutines
have many lines of common code. In the new code, we wrapped the calculation
of cloud radiation properties in one single module “radiation_clouds.f” and
moved all the common code from subroutines progcldXXX to the new subroutine
“radiation_clouds_prop”. This subroutine can be called by RRTMG and RRTMGP.
A single call to the subroutine “radiation_clouds_prop” can connect to the
calculations of the cloud radiation properties for all the microphysics
schemes. There are also many other changes in the following listed program
based on the reviewer’s suggestions/comments.
See https://docs.google.com/document/d/1y-0K3GS6ZwwKipwNNfBWE-k-TQaN3NwrgIbxRkdInYY/edit

physics/GFS_rrtmg_pre.F90 Outdated Show resolved Hide resolved
@climbfuji
Copy link
Collaborator

@mzhangw @ChunxiZhang-NOAA You had previously commented on this PR, which is now approved by two of the CODEOWNERS. Can you please check if your comments have been addressed and approve as well? Thanks.

@@ -151,6 +281,30 @@
type = real
kind = kind_phys
intent = in
[qs_lay]
standard_name = saturation_vapor_pressure
long_name = saturation vapor pressure
Copy link
Collaborator

@ChunxiZhang-NOAA ChunxiZhang-NOAA Mar 25, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Strongly recommend to change 'saturation vapor pressure' to 'model layer mean saturation vapor pressure'. Same for q_lay and relhum.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

done.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChunxiZhang-NOAA @dustinswales I'm sorry that I didn't catch this until now, but why do we need to add "model layer mean" to the standard names? According to the naming rules (https://github.com/ESCOMP/CCPPStandardNames/blob/main/StandardNamesRules.rst, see rules 3 and 4) all variables are located at model layer means by default and therefore it should not need to be specified.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChunxiZhang-NOAA @dustinswales I'm sorry that I didn't catch this until now, but why do we need to add "model layer mean" to the standard names? According to the naming rules (https://github.com/ESCOMP/CCPPStandardNames/blob/main/StandardNamesRules.rst, see rules 3 and 4) all variables are located at model layer means by default and therefore it should not need to be specified.

@grantfirl @dustinswales I am sorry I didn't follow the Standard Names Rules. We should follow the rules then.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChunxiZhang-NOAA @dustinswales I'm sorry that I didn't catch this until now, but why do we need to add "model layer mean" to the standard names? According to the naming rules (https://github.com/ESCOMP/CCPPStandardNames/blob/main/StandardNamesRules.rst, see rules 3 and 4) all variables are located at model layer means by default and therefore it should not need to be specified.

@grantfirl @dustinswales I am sorry I didn't follow the Standard Names Rules. We should follow the rules then.

@grantfirl @dustinswales I mean add change the long name to "model layer mean ...". The rule is only for standard name. So the long name can have a better description.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Strongly recommend to change 'saturation vapor pressure' to 'model layer mean saturation vapor pressure'. Same for q_lay and relhum.

@dustinswales @grantfirl At the beginning here, I meant to change the long name. That's why I wrote 'saturation vapor pressure' instead of 'saturation_vapor_pressure'.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Strongly recommend to change 'saturation vapor pressure' to 'model layer mean saturation vapor pressure'. Same for q_lay and relhum.

@dustinswales @grantfirl At the beginning here, I meant to change the long name. That's why I wrote 'saturation vapor pressure' instead of 'saturation_vapor_pressure'.

@dustinswales If possible, please change it ASAP. Your PR is scheduled today.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChunxiZhang-NOAA
This will have to wait. I only have access to one machine, which is down today, so I can't test any of the changes.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChunxiZhang-NOAA This will have to wait. I only have access to one machine, which is down today, so I can't test any of the changes.

@dustinswales Just want to confirm: you only have access to Hera?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChunxiZhang-NOAA
Yes, Hera only for me :(
well for now...

@@ -256,6 +261,57 @@ subroutine GFS_rrtmgp_pre_run(me, nCol, nLev, nTracers, i_o3, lsswr, lslwr, fhsw
enddo
enddo

!
! Compute layer-thickness between layer boundaries (deltaZ) and layer centers (deltaZc)
Copy link
Collaborator

@ChunxiZhang-NOAA ChunxiZhang-NOAA Mar 25, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why don't use phii and phil to calculate thicknesses? Also notice that layer thickness is calculated in subroutine setaer, which is probably not necessary.
!> -# Compute level height and layer thickness.

  if ( laswflg .or. lalwflg ) then

    lab_do_IMAX : do i = 1, IMAX

      lab_if_flip : if (ivflip == 1) then       ! input from sfc to toa

        do k = 1, NLAY
          prsln(k) = log(prsi(i,k))
        enddo
        prsln(NLP1)= log(prsl(i,NLAY))

        do k = NLAY, 1, -1
          dz(i,k) = rovg * (prsln(k) - prsln(k+1)) * tvly(i,k)
        enddo
        dz(i,NLAY)  = 2.0 * dz(i,NLAY)

        hz(i,1) = f_zero
        do k = 1, NLAY
          hz(i,k+1) = hz(i,k) + dz(i,k)
        enddo

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChunxiZhang-NOAA
This is how the layer thickness is computed for both RRTMG and RRTMGP. In RRTMGP, there are many interstitial variables that were internal variables in RRTMG.

GP is organized differently than G, so some fields are needed by different GP scheme components.
Layer-thickness is one field that is interstitial in GP, but internal in RRTMG (GFS_rrtmg_pre.F90)
"tracer" is another example.
There are many many more Interstitial flat-fields and DDTs in GP, most of which are analogous to some field in the RRTMG implementation.

GP isn't allocating any more fields than G, just more of them are of the Interstitial type.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChunxiZhang-NOAA This is how the layer thickness is computed for both RRTMG and RRTMGP. In RRTMGP, there are many interstitial variables that were internal variables in RRTMG.

GP is organized differently than G, so some fields are needed by different GP scheme components. Layer-thickness is one field that is interstitial in GP, but internal in RRTMG (GFS_rrtmg_pre.F90) "tracer" is another example. There are many many more Interstitial flat-fields and DDTs in GP, most of which are analogous to some field in the RRTMG implementation.

GP isn't allocating any more fields than G, just more of them are of the Interstitial type.
@dustinswales My suggestion is to use the heights for both mid-layer and interface layer directly from the Statein DDT. By this way, the code to calculate the heights by using hydrostatic relation can be deleted. And virtual temperature (tvly) is probably no longer needed. You do need the layer thickness as an interstitial variable. Same for the 'tracer' variable, if change a few lines of the code, this variable is not needed. Due to the limited time for P8 implementation, those can be modified later I think.

Copy link
Collaborator Author

@dustinswales dustinswales Apr 5, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChunxiZhang-NOAA
RRTMGP is much more picky about the inputs than RRTMG.

For example, both radiation schemes use LUTs to calculate the clear-sky gaseous optics. These LUTs have temperature and pressure ranges for which they are valid.

In RRTMG, if a pressure or temperature was out-of-range of the LUT, it would just use the last/first index in the LUTs. "Essentially bounding the problem internally to the limits of the LUTs"

In RRTMGP, if you provide an out-of-range value, it reports an error. So to avoid this we need to "bound the inputs to RRTMGP."

This is why we can't use the statein DDT fields directly.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChunxiZhang-NOAA RRTMGP is much more picky about the inputs than RRTMG.

For example, both radiation schemes use LUTs to calculate the clear-sky gaseous optics. These LUTs have temperature and pressure ranges for which they are valid.

In RRTMG, if a pressure or temperature was out-of-range of the LUT, it would just use the last/first index in the LUTs. "Essentially bounding the problem internally to the limits of the LUTs"

In RRTMGP, if you provide an out-of-range value, it reports an error. So to avoid this we need to "bound the inputs to RRTMGP."

This is why we can't use the statein DDT files directly.

@dustinswales Yes, you told me last time. I agreed to keep some variables like air temperature that need to set boundaries or some variables used across RRTMGP subroutines. But some variables can be avoided by just using the Statein DDT and modifying your interface a little bit. We don't have time to change these this time. Keep it as a future work.

units = flag
dimensions = ()
type = logical
type = integer
intent = in
[i_cldliq]
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Recommend to use the same local names for tracer indexes, e.g., ntcw, ntiw, etc..

real(kind_phys), dimension(:,:), intent(in) :: &
p_lev ! Pressure at model-level interfaces (Pa)
real(kind_phys), dimension(:,:,:),intent(in) :: &
tracer ! Cloud condensate amount in layer by type ()
Copy link
Collaborator

@ChunxiZhang-NOAA ChunxiZhang-NOAA Mar 25, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please double check if you really need this huge array 'tracer'. I checked the code and found that 'tracer' was used by the following code:
GFS_rrtmgp_pre.meta: standard_name = chemical_tracers
GFS_rrtmgp_gfdlmp_pre.meta: standard_name = chemical_tracers
GFS_rrtmgp_thompsonmp_pre.meta: standard_name = chemical_tracers
GFS_rrtmgp_zhaocarr_pre.meta: standard_name = chemical_tracers
rrtmgp_lw_aerosol_optics.meta: standard_name = chemical_tracers
rrtmgp_sw_aerosol_optics.meta: standard_name = chemical_tracers

Let's take a look one by one:
in GFS_rrtmgp_pre.F90, you defined 'tracer' by assigned them to >=0. values.
in GFS_rrtmgp_MP_pre.F90, you assigned a few tracers to some variables:
cld_condensate(1:nCol,1:nLev,1) = tracer(1:nCol,1:nLev,i_cldliq) ! -liquid water
cld_condensate(1:nCol,1:nLev,2) = tracer(1:nCol,1:nLev,i_cldice) ! -ice water
cld_condensate(1:nCol,1:nLev,3) = tracer(1:nCol,1:nLev,i_cldrain) ! -rain water
cld_condensate(1:nCol,1:nLev,4) = tracer(1:nCol,1:nLev,i_cldsnow) + &! -snow + grapuel
tracer(1:nCol,1:nLev,i_cldgrpl)
qc_mp(iCol,iLay) = tracer(iCol,iLay,i_cldliq) / (1.-q_lay(iCol,iLay))
qi_mp(iCol,iLay) = tracer(iCol,iLay,i_cldice) / (1.-q_lay(iCol,iLay))
qs_mp(iCol,iLay) = tracer(iCol,iLay,i_cldsnow) / (1.-q_lay(iCol,iLay))
ni_mp(iCol,iLay) = tracer(iCol,iLay,i_cldice_nc) / (1.-q_lay(iCol,iLay))
if (ltaerosol) then
nc_mp(iCol,iLay) = tracer(iCol,iLay,i_cldliq_nc) / (1.-q_lay(iCol,iLay))
nwfa(iCol,iLay) = tracer(iCol,iLay,i_twa)
in rrtmgp_lw_aerosol_optics.F90 and rrtmgp_sw_aerosol_optics.F90, it calls subroutine setaer, and setaer calls aer_property and aer_perperty_gocart, but none of them actually used 'tracer'.

So I think you can use the original 'qgrs' directly. Maybe I missed something. Please double check.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@ChunxiZhang-NOAA
See comment above.
I'm willing to revisit this at a later time, these changes you suggest are a bit out of scope for this PR and aren't changes I'm comfortable pushing in with so little time. The deadline for p8c/ccpp6 is approaching and this is in the queue to get merged tomorrow.

real(kind_phys), dimension(:,:), intent(in) :: &
effrin_cldrain ! Effective radius for stratiform rain cloud-particles (microns)
real(kind_phys), dimension(:,:), intent(in) :: &
p_lev ! Pressure at model-level interfaces (Pa)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please double check if p_lev is necessary. It would be great if it can use 'prsi' directly.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ditto

call check_error_msg('rrtmgp_lw_rte_run',lw_optical_props_cnvcloudsByBand%increment(lw_optical_props_clrsky))
endif

! Include MYNN-EDMF PBL clouds?
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Similar to convective cloud, I believe the flag here could be generic. It could refer to PBL cloud from any PBL scheme. So it would be great if doGP_sgs_mynn can be changed to doGP_sgs_pbl

cldtausw ! Approx 10.mu band layer cloud optical depth
sw_optical_props_cloudsByBand, & ! RRTMGP DDT: Shortwave optical properties in each band (clouds)
sw_optical_props_cnvcloudsByBand, & ! RRTMGP DDT: Shortwave optical properties in each band (convective cloud)
sw_optical_props_MYNNcloudsByBand,& ! RRTMGP DDT: Shortwave optical properties in each band (MYNN PBL cloud)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Again, it would be great if MYNN can be changed to pbl.

do iDay=1,nDay
! Near IR
scmpsw(idxday(iDay))%nirbm = sum(flux_allsky%bnd_flux_dn_dir(iDay,iSFC,1:ibd-1)) + &
flux_allsky%bnd_flux_dn_dir(iDay,iSFC,ibd)/2.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

/2 -> *0.5

Copy link
Collaborator

@ChunxiZhang-NOAA ChunxiZhang-NOAA left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Some interstitial or global variables need to be reconsidered.

! **Additionally, Conditioned on relative-humidity**
!
! ######################################################################################
subroutine cloud_mp_thompson(nCol, nLev, nTracers, ncnd, i_cldliq, i_cldice, i_cldrain,&
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I recommend that all no CCPP entry subroutines be documented in Doxygen format.

@@ -244,6 +246,12 @@ subroutine GFS_rrtmgp_pre_run(me, nCol, nLev, nTracers, i_o3, lsswr, lslwr, fhsw

! Temperature at layer-interfaces
call cmp_tlev(nCol,nLev,minGPpres,p_lay,t_lay,p_lev,tsfc,t_lev)
do iCol=1,nCol
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

FORTRAN arrays are stored in column-major order. So it should be:
do iLev=1,nlev+1
do iCol=1,nCol
....
end do
end do

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I just found out that you need to change many loop orders in the code. I believe the wrong orders can increase the computational costs.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As Dom mentioned, the optimization is probably catching these.
Additionally, the profiling of the code does not point to these _pre routines, they take a trivial amount of time.
The major contributors to GP's timing are from rrtmgp_lw(sw)_gas_optics and rrtmgp_lw(sw)_rte, ~60% of GP's total runtime

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't recall the configuration, it's been a long time, but here is where the time is being spent in the GP suite:
Screen Shot 2022-04-06 at 2 27 06 PM

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If a compiler can optimize it, it is ok.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@dustinswales Thanks for sharing the statistics.

@grantfirl grantfirl merged commit 0ccc8ac into NCAR:main Apr 8, 2022
@dustinswales dustinswales deleted the enhanced_GP2cld_coupling_tight branch April 15, 2022 15:20
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
CCPP v6 Needed for CCPP v6 public release
Projects
None yet
Development

Successfully merging this pull request may close these issues.

8 participants