Skip to content
Merged
4 changes: 3 additions & 1 deletion physics/GFS_rrtmg_pre.F90
Original file line number Diff line number Diff line change
Expand Up @@ -216,6 +216,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, &
real (kind=kind_phys) :: max_relh
integer :: iflag

integer :: ii_half
integer :: ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte
Expand All @@ -236,7 +237,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, &
LP1 = LM + 1 ! num of in/out levels


gridkm = sqrt(2.0)*sqrt(dx(1)*0.001*dx(1)*0.001)
ii_half = nint(0.5*IM)
gridkm = sqrt(dx(1)*0.001*dx(ii_half)*0.001)

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

@gthompsnWRF I'm confused by this change. The dx variable is an array over the horizontal index but you're referencing with ii_half which is a function of the vertical domain size. Also, the original formulation had sqrt(2) presumably to calculate the diagonal length of the grid, which is now missing. Would you mind explaining this change? I realize that gridkm is only used in the progcld_thompson, so any effects are confined there, but in the spirit of good code reviewing, I'm trying to understand this.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Doggone! That should have been variable IM not LM. Thanks for catching it! Yes, I dropped out the diagonal (sqrt(2)). I really do not know the DX spacing of a cube-sphere grid, but does DX vary a whole lot across 1..IM? Perhaps I don't even need to bother with the half-point and first-point average, I just thought it might better represent avg. grid spacing.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Thanks @gthompsnWRF. Regarding your question, I'm definitely not the authority to answer it, but my guess would be that it depends quite a bit on how big the "chunks" (IM) are and whether the columns are near a cubed sphere edge or not. Also, I don't think that we can assume that there is a specific horizontal relationship amongst the columns. I.e., using the midpoint index of the array (IM/2) doesn't necessarily mean that it will be in the "center" of the current chunk of columns in any physically-consistent sense. My guess would be that if you want one "average" gridkm value, you would actually need to calculate an average, which may be overkill for your use, I don't know.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

I am curious, will this impact the decomposition reproducibility? At the single point, it seems the gridkm will change if we have different decompositions (different mpi tasks) or different chunksize (IM).

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Yes, it likely will. UGWP had exactly the same problem. I think the best solution would be to make gridkm a gridkm(i) with i=1...IM. That will also help with non-uniform grids such as stretched nests and regional grids that are not using the ESG (Extended Schmidt Gnomonic) grid. It's a small overhead but necessary.

Copy link
Copy Markdown

Choose a reason for hiding this comment

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

As others have already mentioned, it is best to use the exact grid spacing - especially as CCPP has application beyond the UFS/GFS with the FV3 dycore.

As I haven't delved very far into the inner workings of CCPP metadata and variable definitions, is it always the case that dx will be provided/stored in meters?

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

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

Not necessarily, it could be in km as well. Meters of course is the standard unit and makes a lot of sense, and it is also the unit in which the host model (fv3atm) stores the variable. If host model and scheme use different units, then CCPP can implement automatic unit conversions, as long as these unit conversions are implemented (for my example above, CCPP knows how to convert m to km and back). These conversions obviously introduce an overhead in the auto-generated caps, therefore it is best to use the same unit as the host model. And, if all host models used the same standard units, then there would be no unit conversions at all.


if (imp_physics == imp_physics_thompson) then
max_relh = 1.5
Expand Down
24 changes: 12 additions & 12 deletions physics/radiation_clouds.f
Original file line number Diff line number Diff line change
Expand Up @@ -3452,7 +3452,7 @@ subroutine progcld_thompson &
do i = 1, IX
cwp(i,k) = max(0.0, clw(i,k,ntcw) * dz(i,k)*1.E6)
crp(i,k) = 0.0
snow_mass_factor = 0.85
snow_mass_factor = 0.90
cip(i,k) = max(0.0, (clw(i,k,ntiw) &
& + (1.0-snow_mass_factor)*clw(i,k,ntsw))*dz(i,k)*1.E6)
if (re_snow(i,k) .gt. snow_max_radius)then
Expand Down Expand Up @@ -4532,8 +4532,8 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, &
DO k = kts,kte

delz = MAX(100., dz(k))
RH_00L = 0.74+MIN(0.25,SQRT(1./(50.0+gridkm*gridkm*delz*0.01)))
RH_00O = 0.82+MIN(0.17,SQRT(1./(50.0+gridkm*gridkm*delz*0.01)))
RH_00L = 0.77+MIN(0.22,SQRT(1./(50.0+gridkm*gridkm*delz*0.01)))
RH_00O = 0.85+MIN(0.14,SQRT(1./(50.0+gridkm*gridkm*delz*0.01)))
RHUM = rh(k)

if (qc(k).ge.1.E-5 .or. qi(k).ge.1.E-5 &
Expand All @@ -4550,7 +4550,7 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, &
RH_00 = RH_00L
ENDIF

tc = t(k) - 273.15
tc = MAX(-80.0, t(k) - 273.15)
if (tc .lt. -12.0) RH_00 = RH_00L

if (tc .gt. 20.0) then
Expand All @@ -4562,12 +4562,12 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, &
if (max_relh.gt.1.12 .or. (.NOT.(modify_qvapor)) ) then
!..For HRRR model, the following look OK.
RHUM = MIN(rh(k), 1.45)
RH_00 = RH_00 + (1.45-RH_00)*(-12.0-tc)/(-12.0+112.)
RH_00 = RH_00 + (1.45-RH_00)*(-12.0-tc)/(-12.0+85.)
CLDFRA(K) = MAX(0.,1.0-SQRT((1.46-RHUM)/(1.46-RH_00)))
else
!..but for the GFS model, RH is way lower.
RHUM = MIN(rh(k), 1.05)
RH_00 = RH_00 + (1.05-RH_00)*(-12.0-tc)/(-12.0+112.)
RH_00 = RH_00 + (1.05-RH_00)*(-12.0-tc)/(-12.0+85.)
CLDFRA(K) = MAX(0.,1.0-SQRT((1.06-RHUM)/(1.06-RH_00)))
endif
endif
Expand Down Expand Up @@ -4789,9 +4789,9 @@ SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs,T,dz,entr, k1,k2,kts,kte)
max_iwc = ABS(qvs(k2)-qvs(k1))

do k = k1, k2
max_iwc = MAX(1.E-5, max_iwc - (qi(k)+qs(k)))
max_iwc = MAX(1.E-6, max_iwc - (qi(k)+qs(k)))
enddo
max_iwc = MIN(2.E-3, max_iwc)
max_iwc = MIN(1.E-4, max_iwc)

this_dz = 0.0
do k = k1, k2
Expand All @@ -4801,7 +4801,7 @@ SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs,T,dz,entr, k1,k2,kts,kte)
this_dz = this_dz + dz(k)
endif
this_iwc = max_iwc*this_dz/tdz
iwc = MAX(5.E-6, this_iwc*(1.-entr))
iwc = MAX(1.E-6, this_iwc*(1.-entr))
if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.203.16) then
qi(k) = qi(k) + cfr(k)*cfr(k)*iwc
endif
Expand Down Expand Up @@ -4830,9 +4830,9 @@ SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte)
! print*, ' max_lwc = ', max_lwc, ' over DZ=',tdz

do k = k1, k2
max_lwc = MAX(1.E-5, max_lwc - qc(k))
max_lwc = MAX(1.E-6, max_lwc - qc(k))
enddo
max_lwc = MIN(2.E-3, max_lwc)
max_lwc = MIN(1.E-4, max_lwc)

this_dz = 0.0
do k = k1, k2
Expand All @@ -4842,7 +4842,7 @@ SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte)
this_dz = this_dz + dz(k)
endif
this_lwc = max_lwc*this_dz/tdz
lwc = MAX(5.E-6, this_lwc*(1.-entr))
lwc = MAX(1.E-6, this_lwc*(1.-entr))
if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.253.16) then
qc(k) = qc(k) + cfr(k)*cfr(k)*lwc
endif
Expand Down