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.

Sign up for GitHub

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

GFS radiation: decorrelation cloud overlap scheme #24

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

GFS radiation: decorrelation cloud overlap scheme #24

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
bc0b623
implement decorrelation length overlap
brianhenn Jun 14, 2023
cfbff72
extend decorrelation-length overlap to SW/LW computations
brianhenn Jun 15, 2023
bde3c55
fix radiation driver calls of SW/LW routines
brianhenn Jun 15, 2023
c4f3f6b
complete LW rad decorrelation length implementation
brianhenn Jun 15, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
37 changes: 32 additions & 5 deletions GFS_layer/GFS_radiation_driver.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -309,6 +309,7 @@ module module_radiation_driver !
use physcons, only: eps => con_eps, &
& epsm1 => con_epsm1, &
& fvirt => con_fvirt &
&, rog => con_rog &
&, rocp => con_rocp
use funcphys, only: fpvs

Expand Down Expand Up @@ -1212,7 +1213,7 @@ subroutine GFS_radiation_driver &
real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP) :: &
htswc, htlwc, gcice, grain, grime, htsw0, htlw0, plyr, tlyr, &
qlyr, olyr, rhly, tvly,qstl, vvel, clw, ciw, prslk1, tem2da, &
tem2db, cldcov, deltaq, cnvc, cnvw, qa, tau067, tau110
dz,delp,tem2db, cldcov, deltaq, cnvc, cnvw, qa, tau067, tau110

real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+1+LTP) :: plvl, tlvl

Expand Down Expand Up @@ -1387,7 +1388,8 @@ subroutine GFS_radiation_driver &
do i = 1, IM
tem1d (i) = QME6
tem2da(i,1) = log( plyr(i,1) )
tem2db(i,1) = 1.0
tem2db(i,1) = log( max(prsmin, plvl(i,1)) )
tem2db(i,LMP) = log( plvl(i,LMP) )
tsfa (i) = tlyr(i,LMK) ! sfc layer air temp
tlvl(i,1) = tlyr(i,1)
tlvl(i,LMP) = tskn(i)
Expand All @@ -1399,13 +1401,15 @@ subroutine GFS_radiation_driver &
qlyr(i,k1) = max( tem1d(i), Statein%qgrs(i,k,1) )
tem1d(i) = min( QME5, qlyr(i,k1) )
tvly(i,k1) = Statein%tgrs(i,k) * (1.0 + fvirt*qlyr(i,k1)) ! virtual T (K)
delp(i,lyb) = plvl(i,lla) - plvl(i,llb)
enddo
enddo

if ( lextop ) then
do i = 1, IM
qlyr(i,lyb) = qlyr(i,lya)
tvly(i,lyb) = tvly(i,lya)
delp(i,lyb) = plvl(i,lla) - plvl(i,llb)
enddo
endif

Expand All @@ -1417,12 +1421,22 @@ subroutine GFS_radiation_driver &
enddo
enddo

! --- ... level height and layer thickness (km)

tem0d = 0.001 * rog
do i = 1, IM
do k = 1, LMK
dz(i,k) = tem0d * (tem2db(i,k+1) - tem2db(i,k)) * tvly(i,k)
enddo
enddo

else ! input data from sfc to toa

do i = 1, IM
tem1d (i) = QME6
tem2da(i,1) = log( plyr(i,1) )
tem2db(i,1) = log( plvl(i,1) )
tem2db(i,LMP) = log( max(prsmin, plvl(i,LMP)) )
tsfa (i) = tlyr(i,1) ! sfc layer air temp
tlvl(i,1) = tskn(i)
tlvl(i,LMP) = tlyr(i,LMK)
Expand All @@ -1433,13 +1447,15 @@ subroutine GFS_radiation_driver &
qlyr(i,k) = max( tem1d(i), Statein%qgrs(i,k,1) )
tem1d(i) = min( QME5, qlyr(i,k) )
tvly(i,k) = Statein%tgrs(i,k) * (1.0 + fvirt*qlyr(i,k)) ! virtual T (K)
delp(i,lyb) = plvl(i,lla) - plvl(i,llb)
enddo
enddo

if ( lextop ) then
do i = 1, IM
qlyr(i,lyb) = qlyr(i,lya)
tvly(i,lyb) = tvly(i,lya)
delp(i,lyb) = plvl(i,lla) - plvl(i,llb)
enddo
endif

Expand All @@ -1451,6 +1467,15 @@ subroutine GFS_radiation_driver &
enddo
enddo

! --- ... level height and layer thickness (km)

tem0d = 0.001 * rog
do i = 1, IM
do k = LMK, 1, -1
dz(i,k) = tem0d * (tem2db(i,k) - tem2db(i,k+1)) * tvly(i,k)
enddo
enddo

endif ! end_if_ivflip

!> - Check for daytime points for SW radiation.
Expand Down Expand Up @@ -1585,7 +1610,7 @@ subroutine GFS_radiation_driver &
clw, Grid%xlat, Grid%xlon, Sfcprop%slmsk,&
tracer1(:,1:lmk,Model%ntclamt), im, lmk, &
lmp, &
clouds, cldsa, mtopa, mbota) ! --- outputs
clouds, cldsa, mtopa, mbota, de_lgth) ! --- outputs

elseif (icmphys == 5) then ! zhao/moorthi's prognostic cloud scheme + pdf cloud & cnvc and cnvw

Expand Down Expand Up @@ -1681,14 +1706,16 @@ subroutine GFS_radiation_driver &
if (Model%swhtr) then
call swrad (plyr, plvl, tlyr, tlvl, qlyr, olyr, & ! --- inputs
gasvmr, clouds, Tbd%icsdsw, faersw, &
sfcalb, Radtend%coszen, Model%solcon, &
sfcalb, dz, delp, de_lgth, &
Radtend%coszen, Model%solcon, &
nday, idxday, im, lmk, lmp, Model%lprnt,&
htswc, Diag%topfsw, Radtend%sfcfsw, & ! --- outputs
hsw0=htsw0, fdncmp=scmpsw, tau067=tau067) ! --- optional
else
call swrad (plyr, plvl, tlyr, tlvl, qlyr, olyr, & ! --- inputs
gasvmr, clouds, Tbd%icsdsw, faersw, &
sfcalb, Radtend%coszen, Model%solcon, &
sfcalb, dz, delp, de_lgth, &
Radtend%coszen, Model%solcon, &
nday, idxday, IM, LMK, LMP, Model%lprnt,&
htswc, Diag%topfsw, Radtend%sfcfsw, & ! --- outputs
FDNCMP=scmpsw, tau067=tau067) ! --- optional
Expand Down
95 changes: 83 additions & 12 deletions gsmphys/radiation_clouds.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1845,10 +1845,10 @@ subroutine progcld4 &

! --- inputs:
& ( plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw, &
& xlat,xlon,slmsk,cldtot, &
& xlat,xlon,slmsk,cldtot, dz &
& IX, NLAY, NLP1, &
! --- outputs:
& clouds,clds,mtop,mbot &
& clouds,clds,mtop,mbot,de_lgth &
& )

! ================= subprogram documentation block ================ !
Expand Down Expand Up @@ -1887,6 +1887,7 @@ subroutine progcld4 &
! range, otherwise see in-line comment !
! xlon (IX) : grid longitude in radians (not used) !
! slmsk (IX) : sea/land mask array (sea:0,land:1,sea-ice:2) !
! dz (ix,nlay) : layer thickness (km) !
! IX : horizontal dimention !
! NLAY,NLP1 : vertical layer/level dimensions !
! !
Expand All @@ -1905,6 +1906,7 @@ subroutine progcld4 &
! clds (IX,5) : fraction of clouds for low, mid, hi, tot, bl !
! mtop (IX,3) : vertical indices for low, mid, hi cloud tops !
! mbot (IX,3) : vertical indices for low, mid, hi cloud bases !
! de_lgth(ix) : clouds decorrelation length (km) !
! !
! module variables: !
! ivflip : control flag of vertical index direction !
Expand All @@ -1926,7 +1928,7 @@ subroutine progcld4 &
integer, intent(in) :: IX, NLAY, NLP1

real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, &
& tlyr, tvly, qlyr, qstl, rhly, clw, cldtot
& tlyr, tvly, qlyr, qstl, rhly, clw, cldtot, dz

real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, &
& slmsk
Expand All @@ -1935,14 +1937,14 @@ subroutine progcld4 &
real (kind=kind_phys), dimension(:,:,:), intent(out) :: clouds

real (kind=kind_phys), dimension(:,:), intent(out) :: clds

real (kind=kind_phys), dimension(:), intent(out) :: de_lgth
integer, dimension(:,:), intent(out) :: mtop,mbot

! --- local variables:
real (kind=kind_phys), dimension(IX,NLAY) :: cldcnv, &
& cwp, cip, crp, csp, rew, rei, res, rer, delp, tem2d, clwf

real (kind=kind_phys) :: ptop1(IX,NK_CLDS+1)
real (kind=kind_phys) :: ptop1(IX,NK_CLDS+1), rxlat(ix)

real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh, value, &
& tem1, tem2, tem3
Expand Down Expand Up @@ -1999,6 +2001,11 @@ subroutine progcld4 &
! ptopc(k,i): top presure of each cld domain (k=1-4 are sfc,L,m,h;
! --- i=1,2 are low-lat (<45 degree) and pole regions)

do i =1, IX
rxlat(i) = abs( xlat(i) / con_pi ) ! if xlat in pi/2 -> -pi/2 range
! rxlat(i) = abs(0.5 - xlat(i)/con_pi) ! if xlat in 0 -> pi range
enddo

do id = 1, 4
tem1 = ptopc(id,2) - ptopc(id,1)

Expand Down Expand Up @@ -2107,7 +2114,14 @@ subroutine progcld4 &
clouds(i,k,9) = rei(i,k)
enddo
enddo


! --- ... estimate clouds decorrelation length in km
! this is only a tentative test, need to consider change later
if ( iovr == 3 ) then
do i = 1, ix
de_lgth(i) = max( 0.6, 2.78-4.6*rxlat(i) )
enddo
endif

! --- compute low, mid, high, total, and boundary layer cloud fractions
! and clouds top/bottom layer indices for low, mid, and high clouds.
Expand All @@ -2117,7 +2131,7 @@ subroutine progcld4 &

call gethml &
! --- inputs:
& ( plyr, ptop1, cldtot, cldcnv, &
& ( plyr, ptop1, cldtot, cldcnv, dz, de_lgth, &
& IX,NLAY, &
! --- outputs:
& clds, mtop, mbot &
Expand Down Expand Up @@ -3611,6 +3625,8 @@ end subroutine diagcld1
!! (sfc,low,mid,high) in mb (100Pa)
!> \param cldtot (IX,NLAY), total or stratiform cloud profile in fraction
!> \param cldcnv (IX,NLAY), convective cloud (for diagnostic scheme only)
!> \param dz (IX,NLAY), layer thickness (km)
!> \param de_lgth (IX), clouds decorrelation length (km)
!> \param IX horizontal dimension
!> \param NLAY vertical layer dimensions
!> \param clds (IX,5), fraction of clouds for low, mid, hi, tot, bl
Expand All @@ -3621,7 +3637,7 @@ end subroutine diagcld1
!! @{
!----------------------------------- !
subroutine gethml &
& ( plyr, ptop1, cldtot, cldcnv, & ! --- inputs:
& ( plyr, ptop1, cldtot, cldcnv, dz, de_lgth, & ! --- inputs:
& IX, NLAY, &
& clds, mtop, mbot & ! --- outputs:
& )
Expand Down Expand Up @@ -3651,6 +3667,8 @@ subroutine gethml &
! (sfc,low,mid,high) in mb (100Pa) !
! cldtot(IX,NLAY) : total or straiform cloud profile in fraction !
! cldcnv(IX,NLAY) : convective cloud (for diagnostic scheme only) !
! dz (ix,nlay) : layer thickness (km) !
! de_lgth(ix) : clouds vertical de-correlation length (km) !
! IX : horizontal dimention !
! NLAY : vertical layer dimensions !
! !
Expand All @@ -3668,6 +3686,7 @@ subroutine gethml &
! iovr : control flag for cloud overlap !
! =0 random overlapping clouds !
! =1 max/ran overlapping clouds !
! =3 decorr-length ovlp ( for mcica only ) !
! !
! !
! ==================== end of description ===================== !
Expand All @@ -3678,16 +3697,17 @@ subroutine gethml &
integer, intent(in) :: IX, NLAY

real (kind=kind_phys), dimension(:,:), intent(in) :: plyr, ptop1, &
& cldtot, cldcnv
& cldtot, cldcnv, dz
real (kind=kind_phys), dimension(:), intent(in) :: de_lgth

! --- outputs
real (kind=kind_phys), dimension(:,:), intent(out) :: clds

integer, dimension(:,:), intent(out) :: mtop, mbot

! --- local variables:
real (kind=kind_phys) :: cl1(IX), cl2(IX)
real (kind=kind_phys) :: pcur, pnxt, ccur, cnxt
real (kind=kind_phys) :: cl1(IX), cl2(IX), dz1(ix)
real (kind=kind_phys) :: pcur, pnxt, ccur, cnxt, alfa

integer, dimension(IX):: idom, kbt1, kth1, kbt2, kth2
integer :: i, k, id, id1, kstr, kend, kinc
Expand Down Expand Up @@ -3737,7 +3757,7 @@ subroutine gethml &
clds(i,4) = 1.0 - cl1(i) ! save total cloud
enddo

else ! max/ran overlap
elseif ( iovr == 1 ) then ! max/ran overlap

do k = kstr, kend, kinc
do i = 1, IX
Expand All @@ -3761,6 +3781,57 @@ subroutine gethml &
clds(i,4) = 1.0 - cl1(i) * cl2(i) ! save total cloud
enddo

elseif ( iovr == 2 ) then ! maximum overlap all levels

cl1(:) = 0.0

do k = kstr, kend, kinc
do i = 1, IX
ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) ))
if (ccur >= climit) cl1(i) = max( cl1(i), ccur )
enddo

if (k == llyr) then
do i = 1, IX
clds(i,5) = cl1(i) ! save bl cloud
enddo
endif
enddo

do i = 1, IX
clds(i,4) = cl1(i) ! save total cloud
enddo

elseif ( iovr == 3 ) then ! random if clear-layer divided,
! otherwise de-corrlength method
do i = 1, ix
dz1(i) = - dz(i,kstr)
enddo
do k = kstr, kend, kinc
do i = 1, ix
ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) ))
if (ccur >= climit) then ! cloudy layer
alfa = exp( -0.5*((dz1(i)+dz(i,k)))/de_lgth(i) )
dz1(i) = dz(i,k)
cl2(i) = alfa * min(cl2(i), (1.0 - ccur)) & ! maximum part
& + (1.0 - alfa) * (cl2(i) * (1.0 - ccur)) ! random part
else ! clear layer
cl1(i) = cl1(i) * cl2(i)
cl2(i) = 1.0
if (k /= kend) dz1(i) = -dz(i,k+kinc)
endif
enddo
if (k == llyr) then
do i = 1, ix
clds(i,5) = 1.0 - cl1(i) * cl2(i) ! save bl cloud
enddo
endif
enddo

do i = 1, IX
clds(i,4) = 1.0 - cl1(i) * cl2(i) ! save total cloud
enddo

endif ! end_if_iovr

! --- high, mid, low clouds, where cl1, cl2 are cloud fractions
Expand Down
Loading