Skip to content

Commit

Permalink
Merge pull request #519 from GEOS-ESM/feature/pnorris/#512-add_additi…
Browse files Browse the repository at this point in the history 
…onal_radiation_RRTMG_TAU_and_Water_Path_Diagnostics

Feature/pnorris/#512 add additional radiation rrtmg tau and water path diagnostics
  • Loading branch information
@sdrabenh
sdrabenh authored Feb 3, 2022
2 parents a70c6ab + ea5d84d commit d5ae28c
Show file tree
Hide file tree
Showing 7 changed files with 905 additions and 220 deletions.
187 changes: 148 additions & 39 deletions ...omp/GEOSphysics_GridComp/GEOSradiation_GridComp/GEOS_RadiationShared/cloud_subcol_gen.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,18 @@ module cloud_subcol_gen
! cloud output explicit by changing to a logical cldy_stoch in {true,false},
! which may permit some exterior efficiency improvements.

! PMN 2022/01: Original version only worked for layer ordering from 1=bottom
! upward. New version auto-detects vertical ordering and works for either.
! See details under generate_stochastic_clouds() and clearCounts_threeBand().
! NB: The POPULATION statistics of the subcolumn ensemble should not depend
! on the vertical ordering. But if you want an exact SAMPLE replication of
! a generation (such as the bottom-up generation used in RRTMG) then you
! should ensure the same input ordering (in the RRTMG case, by reversing the
! vertical ordering from the GEOS-5 top-down ordering BEFORE input to these
! routines). Otherwise (i.e. if POPULATION statistics equivalence is OK)
! then save time by avoiding any explicit array re-ordering and let the
! auto-detect feature do its job.

use cloud_condensate_inhomogeneity, only: &
condensate_inhomogeneous, zcw_lookup
use iso_fortran_env, only : error_unit
Expand Down Expand Up @@ -156,7 +168,7 @@ subroutine generate_stochastic_clouds( &
! of the initial seeding on each gridcolumn's PRNG stream. As such, it provides a
! way of generating a different set of random numbers for say LW and SW at the same
! timetep, when the layer pressures <play> used for seed generation are the same.
! For example, LW can be run with [1,2,3,4] and SW with [4,2,1,3].
! For example, LW can be run with [1,2,3,4] and SW with [4,3,2,1].
! NOTE: Using the same seed_order for LW and SW still does not mean using the "same
! subcolumn cloud ensemble" because the number of subcolumns (g-points) in LW and SW
! differ. Also, of course, McICA .ne. ICA, meaning that each subcolumn in McICA gets
Expand All @@ -167,13 +179,23 @@ subroutine generate_stochastic_clouds( &
! seeds for LW and SW will create an unwanted correlation of some sort, they can use
! a different seed_order for each.
!
! Vertical ordering:
! The ordering is determined by play(1,1) > play(nlay,1) for surface=1 upward, and vica
! versa. Either vertical ordering is permitted, but the generation always proceeds from
! layer 1 to layer nlay. That shouldn't matter in a statistical sense. The seed order
! does NOT need to be changed depending on the vertical ordering. The four layers nearest
! the surface are always used (in the seed_order sense). Of course, the vertical ordering
! of all profile inputs (zmid,play,cldfrac,ciwp,clwp) should be the same, and the output
! vertical ordering (of [cldy|ciwp|clwp]_stoch) is the same as the input ordering.
!---------------------------------------------------------------------------------------

integer, intent(in) :: dncol ! Dimensioned number of gridcols
integer, intent(in) :: ncol ! Actual number of gridcols
integer, intent(in) :: nsubcol ! #Subcols to generate / gridcol
integer, intent(in) :: nlay ! Number of model layers

! Note: only absolute layer separations, abs(zmid(k1,:)-zmid(k2,:)), matter
! so "zero-reference" of zmid is not important
real, intent(in) :: zmid (nlay,dncol) ! Height of midpoints [m]
real, intent(in) :: alat (dncol) ! Latitude of gridcolumn [radians]
integer, intent(in) :: doy ! Day of year
Expand Down Expand Up @@ -206,6 +228,9 @@ subroutine generate_stochastic_clouds( &

! ----- Locals -----

! vertical ordering
logical :: surface_at_one

! correlation length for cloud presence and condensate
real, dimension(dncol) :: adl, rdl

Expand Down Expand Up @@ -235,6 +260,9 @@ subroutine generate_stochastic_clouds( &
! indices
integer :: icol, isubcol, ilay, n

! detect vertical ordering
surface_at_one = (play(1,1) > play(nlay,1))

! save for speed
cond_inhomo = condensate_inhomogeneous()

Expand Down Expand Up @@ -284,12 +312,13 @@ subroutine generate_stochastic_clouds( &
do icol = 1,ncol

! exponential inter-layer correlations ...
! [alpha|rcorr](k) is correlation between layers k and k-1
do ilay = 2,nlay
alpha(ilay) = exp( -(zmid(ilay,icol)-zmid(ilay-1,icol)) / adl(icol) )
alpha(ilay) = exp( -abs(zmid(ilay,icol)-zmid(ilay-1,icol)) / adl(icol) )
enddo
if (cond_inhomo) then
do ilay = 2,nlay
rcorr(ilay) = exp( -(zmid(ilay,icol)-zmid(ilay-1,icol)) / rdl(icol) )
rcorr(ilay) = exp( -abs(zmid(ilay,icol)-zmid(ilay-1,icol)) / rdl(icol) )
enddo
endif

Expand Down Expand Up @@ -345,7 +374,12 @@ subroutine generate_stochastic_clouds( &
! pseed.

! isolate lower pressures for seeding
pseed = play(1:4,icol) * 100. ! [Pa]
if (surface_at_one) then
pseed = play(1:4,icol)
else
pseed = play(nlay:nlay-3:-1,icol)
end if
pseed = pseed * 100. ! [Pa]

! PMN 2021-11-12: Try [daPa=decaPascals] instead.
! pseed = play(1:4,icol) * 10. ! [daPa]
Expand Down Expand Up @@ -557,9 +591,24 @@ subroutine clearCounts_threeBand( &
! ----------------------------------------------------------------
! Count clear subcolumns per gridcolumn for whole column
! and for three pressure bands.
! layers [1, cloudLM] are in low pressure band
! layers [cloudLM+1, cloudMH] are in mid pressure band
! layers [cloudMH+1, nlay ] are in high pressure band
! Verical ordering is auto-detected via:
! cloudLM < cloudMH => surface at level 1
! cloudLM > cloudMH => surface at level nlay
! and concomitant layer definitions:
! layers [1, cloudLM] are in low pressure band
! layers [cloudLM+1, cloudMH] are in mid pressure band
! layers [cloudMH+1, nlay ] are in high pressure band
! for surface at level 1,
! with cloudLM = nlay - LCLDLM + 1
! and cloudMH = nlay - LCLDMH + 1,
! and
! layers [1, cloudMH-1] are in high pressure band
! layers [cloudMH, cloudLM-1] are in mid pressure band
! layers [cloudLM, nlay ] are in low pressure band
! for surface at level nlay (TOA at level 1)
! with cloudLM = LCLDLM
! and cloudMH = LCLDMH,
! all in terms of GEOS-5's LCLDLM and LCLDMH (TOA at level 1).
! ----------------------------------------------------------------

integer, intent(in) :: dncol ! Dimensioned number of gridcols
Expand Down Expand Up @@ -597,41 +646,101 @@ subroutine clearCounts_threeBand( &
if (.not. cloud_found) &
clearCnts(1,icol) = clearCnts(1,icol) + 1

! high pressure band
cloud_found = .false.
do ilay = cloudMH+1, nlay
if (cldy_stoch(ilay,isubcol,icol)) then
cloud_found = .true.
exit
endif
enddo
if (.not. cloud_found) &
clearCnts(2,icol) = clearCnts(2,icol) + 1
enddo ! subcolumns
enddo ! gridcolumns

if (cloudLM < cloudMH) then

! mid pressure band
cloud_found = .false.
do ilay = cloudLM+1, cloudMH
if (cldy_stoch(ilay,isubcol,icol)) then
cloud_found = .true.
exit
endif
enddo
if (.not. cloud_found) &
clearCnts(3,icol) = clearCnts(3,icol) + 1
! surface at level 1 ordering ...

! low pressure band
cloud_found = .false.
do ilay = 1, cloudLM
if (cldy_stoch(ilay,isubcol,icol)) then
cloud_found = .true.
exit
endif
enddo
if (.not. cloud_found) &
clearCnts(4,icol) = clearCnts(4,icol) + 1
do icol = 1,ncol
do isubcol = 1,nsubcol

enddo ! subcolumns
enddo ! gridcolumns
! low pressure band
cloud_found = .false.
do ilay = 1, cloudLM
if (cldy_stoch(ilay,isubcol,icol)) then
cloud_found = .true.
exit
endif
enddo
if (.not. cloud_found) &
clearCnts(4,icol) = clearCnts(4,icol) + 1

! mid pressure band
cloud_found = .false.
do ilay = cloudLM+1, cloudMH
if (cldy_stoch(ilay,isubcol,icol)) then
cloud_found = .true.
exit
endif
enddo
if (.not. cloud_found) &
clearCnts(3,icol) = clearCnts(3,icol) + 1

! high pressure band
cloud_found = .false.
do ilay = cloudMH+1, nlay
if (cldy_stoch(ilay,isubcol,icol)) then
cloud_found = .true.
exit
endif
enddo
if (.not. cloud_found) &
clearCnts(2,icol) = clearCnts(2,icol) + 1

enddo ! subcolumns
enddo ! gridcolumns

elseif (cloudLM > cloudMH) then

! TOA at level 1 ordering ...

do icol = 1,ncol
do isubcol = 1,nsubcol

! high pressure band
cloud_found = .false.
do ilay = 1, cloudMH-1
if (cldy_stoch(ilay,isubcol,icol)) then
cloud_found = .true.
exit
endif
enddo
if (.not. cloud_found) &
clearCnts(2,icol) = clearCnts(2,icol) + 1

! mid pressure band
cloud_found = .false.
do ilay = cloudMH, cloudLM-1
if (cldy_stoch(ilay,isubcol,icol)) then
cloud_found = .true.
exit
endif
enddo
if (.not. cloud_found) &
clearCnts(3,icol) = clearCnts(3,icol) + 1

! low pressure band
cloud_found = .false.
do ilay = cloudLM, nlay
if (cldy_stoch(ilay,isubcol,icol)) then
cloud_found = .true.
exit
endif
enddo
if (.not. cloud_found) &
clearCnts(4,icol) = clearCnts(4,icol) + 1

enddo ! subcolumns
enddo ! gridcolumns

else

write(error_unit,*) 'file:', __FILE__, ', line:', __LINE__
error stop 'invalid pressure super-layers!'

end if

end subroutine clearCounts_threeBand

Expand Down
110 changes: 65 additions & 45 deletions ...omp/GEOSphysics_GridComp/GEOSradiation_GridComp/GEOSirrad_GridComp/GEOS_IrradGridComp.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -783,33 +783,41 @@ subroutine SetServices ( GC, RC )
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, __RC__ )

call MAPL_AddExportSpec(GC, &
SHORT_NAME = 'CLDTTLW', &
LONG_NAME = 'total_cloud_area_fraction_rrtmg_lw', &
UNITS = '1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, __RC__ )

call MAPL_AddExportSpec(GC, &
SHORT_NAME = 'CLDHILW', &
LONG_NAME = 'high-level_cloud_area_fraction_rrtmg_lw', &
UNITS = '1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, __RC__ )

call MAPL_AddExportSpec(GC, &
SHORT_NAME = 'CLDMDLW', &
LONG_NAME = 'mid-level_cloud_area_fraction_rrtmg_lw', &
UNITS = '1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, __RC__ )

call MAPL_AddExportSpec(GC, &
SHORT_NAME = 'CLDLOLW', &
LONG_NAME = 'low_level_cloud_area_fraction_rrtmg_lw', &
UNITS = '1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, __RC__ )
! Note: the four CLDxxLW diagnostics below represent super-layer cloud
! fractions based on the subcolumn cloud generation called in RRTMG LW.
! They are global fields but generated only at the LW REFRESH frequency,
! NOT at the heartbeat. As such, they are useful for diagnostic comparisons
! with CLDTT above and with the full CLDxx set from the Solar GC. But they
! should NOT be used to subsample fields that are produced on the model
! heartbeat (e.g. subsampling for cloud presence).

call MAPL_AddExportSpec(GC, &
SHORT_NAME = 'CLDTTLW', &
LONG_NAME = 'total_cloud_area_fraction_rrtmg_lw_REFRESH', &
UNITS = '1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, __RC__ )

call MAPL_AddExportSpec(GC, &
SHORT_NAME = 'CLDHILW', &
LONG_NAME = 'high-level_cloud_area_fraction_rrtmg_lw_REFRESH', &
UNITS = '1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, __RC__ )

call MAPL_AddExportSpec(GC, &
SHORT_NAME = 'CLDMDLW', &
LONG_NAME = 'mid-level_cloud_area_fraction_rrtmg_lw_REFRESH', &
UNITS = '1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, __RC__ )

call MAPL_AddExportSpec(GC, &
SHORT_NAME = 'CLDLOLW', &
LONG_NAME = 'low_level_cloud_area_fraction_rrtmg_lw_REFRESH', &
UNITS = '1', &
DIMS = MAPL_DimsHorzOnly, &
VLOCATION = MAPL_VLocationNone, __RC__ )

! Irrad does not have a "real" internal state. To update the net_longwave_flux
! due to the change of surface temperature every time step, we keep
Expand Down Expand Up @@ -1701,27 +1709,39 @@ subroutine Lw_Driver(IM,JM,LM,LATS,LONS,CoresPerNode,RC)
REFF(:,:,:,KRAIN ) = RR * 1.0e6
REFF(:,:,:,KSNOW ) = RS * 1.0e6

! Determine the model level seperating high and middle clouds
!------------------------------------------------------------
! Determine the model level separating high-middle and low-middle clouds
!-----------------------------------------------------------------------

LCLDMH = 1
do K = 1, LM
if( PREF(K) >= PRS_MID_HIGH ) then
LCLDMH = K
exit
end if
end do
_ASSERT(PRS_MID_HIGH > PREF(1) , 'mid-high pressure band boundary too high!')
_ASSERT(PRS_LOW_MID > PRS_MID_HIGH, 'pressure band misordering!')
_ASSERT(PRS_LOW_MID < PREF(LM) , 'low-mid pressure band boundary too low!')

! Determine the model level seperating low and middle clouds
!-----------------------------------------------------------
! find mid-high interface level
k = 1
do while ( PREF(k) < PRS_MID_HIGH )
k=k+1
end do
LCLDMH = k
! Guaranteed that LCLDMH > 1 (by first ASSERT above)
! and that PREF(LCLDMH) >= PRS_MID_HIGH (by while loop)

LCLDLM = LM
do K = 1, LM
if( PREF(K) >= PRS_LOW_MID ) then
LCLDLM = K
exit
end if
! find low-mid interface level
do while ( PREF(k) < PRS_LOW_MID )
k=k+1
end do
LCLDLM = k
! Guaranteed that LCLDLM <= LM (by third assert above)
! and that PREF(LCLDLM) >= PRS_LOW_MID (by while loop)

! But it's still possible that LCLDLM == LCLDMH if the
! interface pressures are too close. We now ASSERT to
! prevent this.
_ASSERT(LCLDMH < LCLDLM, 'PRS_LOW_MID and PRS_MID_HIGH are too close!')

! now we have 1 < LCLDMH < LCLDLM <= LM and can use:
! layers [1, LCLDMH-1] are in high pressure band
! layers [LCLDMH, LCLDLM-1] are in mid pressure band
! layers [LCLDLM, LM ] are in low pressure band

call MAPL_GetPointer(EXPORT, CLDTTLW, 'CLDTTLW', RC=STATUS); VERIFY_(STATUS)
call MAPL_GetPointer(EXPORT, CLDHILW, 'CLDHILW', RC=STATUS); VERIFY_(STATUS)
Expand Down Expand Up @@ -3345,7 +3365,7 @@ subroutine Lw_Driver(IM,JM,LM,LATS,LONS,CoresPerNode,RC)
do I = 1,IM
IJ = IJ + 1

! convert super-band clearCounts to cloud fractions
! convert super-layer clearCounts to cloud fractions
if(associated(CLDTTLW)) then
CLDTTLW(I,J) = 1.0 - CLEARCOUNTS(IJ,1)/float(NGPTLW)
endif
Expand Down
Loading

0 comments on commit d5ae28c

Please sign in to comment.