Radiation Improvements for prototype 8 #871
Conversation
…agnostics.F90 and related meta files
…mpson MP plus radiatively active convective cloud.
…cpp-physics into enhanced_GP2cld_coupling
…nd pbl (MYNN) clouds to RRTMGP.
dustinswales
requested review from
mjiacono,
Qingfu-Liu,
climbfuji,
SamuelTrahanNOAA,
grantfirl,
mzhangw,
panll and
mkavulich
as code owners
|
@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 | |||
There was a problem hiding this comment.
Strongly recommend to change 'saturation vapor pressure' to 'model layer mean saturation vapor pressure'. Same for q_lay and relhum.
There was a problem hiding this comment.
@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.
There was a problem hiding this comment.
@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.
There was a problem hiding this comment.
@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.
There was a problem hiding this comment.
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'.
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
@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.
There was a problem hiding this comment.
@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?
There was a problem hiding this comment.
@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) | |||
There was a problem hiding this comment.
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
There was a problem hiding this comment.
@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.
There was a problem hiding this comment.
@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.
There was a problem hiding this comment.
@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.
There was a problem hiding this comment.
@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] |
There was a problem hiding this comment.
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 () |
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
@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) |
There was a problem hiding this comment.
Please double check if p_lev is necessary. It would be great if it can use 'prsi' directly.
| call check_error_msg('rrtmgp_lw_rte_run',lw_optical_props_cnvcloudsByBand%increment(lw_optical_props_clrsky)) | ||
| endif | ||
|
|
||
| ! Include MYNN-EDMF PBL clouds? |
There was a problem hiding this comment.
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) |
There was a problem hiding this comment.
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. |
| ! **Additionally, Conditioned on relative-humidity** | ||
| ! | ||
| ! ###################################################################################### | ||
| subroutine cloud_mp_thompson(nCol, nLev, nTracers, ncnd, i_cldliq, i_cldice, i_cldrain,& |
There was a problem hiding this comment.
I recommend that all no CCPP entry subroutines be documented in Doxygen format.
…nced_GP2cld_coupling_tight
This reverts commit 2617af6.
physics/GFS_rrtmgp_pre.F90
Outdated
| @@ -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 | |||
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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.
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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:
There was a problem hiding this comment.
If a compiler can optimize it, it is ok.
There was a problem hiding this comment.
@dustinswales Thanks for sharing the statistics.
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:
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