From 630d3b1eaedf7b276cba92a4cca7b04f92da9f05 Mon Sep 17 00:00:00 2001 From: Laurie Carson Date: Thu, 7 Jun 2018 10:21:22 -0600 Subject: [PATCH 1/3] Update rrtmg radiation code from the FV3 master (6/6/2018, ~fv3gfs v1) 1. Add two new entries in the "clouds" array (diagnositc output only?) 2. surface perturbations 3. Re-order a few if-tests in radlw_main.f modified: GFS_layer/GFS_radiation_driver.F90 modified: physics/radiation_clouds.f modified: physics/radiation_surface.f modified: physics/radlw_main.f modified: physics/radsw_main.f --- GFS_layer/GFS_radiation_driver.F90 | 2 + physics/radiation_clouds.f | 1382 ++++++++++++++++++++++++---- physics/radiation_surface.f | 291 +++++- physics/radlw_main.f | 588 ++++++------ physics/radsw_main.f | 14 +- 5 files changed, 1831 insertions(+), 446 deletions(-) diff --git a/GFS_layer/GFS_radiation_driver.F90 b/GFS_layer/GFS_radiation_driver.F90 index dba3535fd..b5473b765 100644 --- a/GFS_layer/GFS_radiation_driver.F90 +++ b/GFS_layer/GFS_radiation_driver.F90 @@ -1256,6 +1256,7 @@ subroutine GFS_radiation_driver & cld_iwp=clouds(:,:,4), cld_ref_ice=clouds(:,:,5), & cld_rwp=clouds(:,:,6), cld_ref_rain=clouds(:,:,7), & cld_swp=clouds(:,:,8), cld_ref_snow=clouds(:,:,9), & + cld_od_total=clouds(:,:,10), cld_od_layer w=clouds(:,:,11), & errmsg=errmsg, errflg=errflg) ! DH* @@ -1290,6 +1291,7 @@ subroutine GFS_radiation_driver & cld_iwp=clouds(:,:,4), cld_ref_ice=clouds(:,:,5), & cld_rwp=clouds(:,:,6), cld_ref_rain=clouds(:,:,7), & cld_swp=clouds(:,:,8), cld_ref_snow=clouds(:,:,9), & + cld_od_total=clouds(:,:,10), cld_od_layer w=clouds(:,:,11), & errmsg=errmsg, errflg=errflg) !CCPP: L1718-1747 diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index f9ff28cef..109947d6e 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -13,7 +13,7 @@ ! ! ! 'cld_init' --- initialization routine ! ! inputs: ! -! ( si, NLAY, me ) ! +! ( si, NLAY, imp_physics, me ) ! ! outputs: ! ! ( none ) ! ! ! @@ -140,6 +140,7 @@ ! convective cloud cover and water for radiation ! ! ! ! jul 2014 s. moorthi - merging with gfs version ! +! feb 2017 a. cheng - add odepth output, effective radius input ! ! ! !!!!! ========================================================== !!!!! !!!!! end descriptions !!!!! @@ -231,7 +232,7 @@ module module_radiation_clouds ! !........................................! ! - use physparam, only : icldflg, icmphys, iovrsw, iovrlw, & + use physparam, only : icldflg, iovrsw, iovrlw, & & lcrick, lcnorm, lnoprec, & & ivflip, kind_phys, kind_io4 use physcons, only : con_fvirt, con_ttp, con_rocp, & @@ -252,10 +253,10 @@ module module_radiation_clouds ! ! --- set constant parameters real (kind=kind_phys), parameter :: gfac=1.0e5/con_g & &, gord=con_g/con_rd -! number of fields in cloud array - integer, parameter, public :: NF_CLDS = 9 -! number of cloud vertical domains - integer, parameter, public :: NK_CLDS = 3 +!> number of fields in cloud array + integer, parameter, public :: NF_CLDS = 11 +!> number of cloud vertical domains + integer, parameter, public :: NK_CLDS = 3 ! pressure limits of cloud domain interfaces (low,mid,high) in mb (0.1kPa) real (kind=kind_phys), save :: ptopc(NK_CLDS+1,2) @@ -319,8 +320,8 @@ module module_radiation_clouds ! ! maximum-random cloud overlapping method integer :: iovr = 1 - public progcld1, progcld2, progcld3, progclduni, diagcld1, & - & cld_init + public progcld1, progcld2, progcld3, progcld4, progclduni, & + & diagcld1, cld_init, progcld5, progcld4o ! ================= @@ -336,7 +337,7 @@ module module_radiation_clouds ! !>\section gen_cld_init cld_init General Algorithm !> @{ subroutine cld_init & - & ( si, NLAY, me ) ! --- inputs + & ( si, NLAY, imp_physics, me ) ! --- inputs ! --- outputs: ! ( none ) @@ -349,6 +350,7 @@ subroutine cld_init & ! inputs: ! ! si (L+1) : model vertical sigma layer interface ! ! NLAY : vertical layer number ! +! imp_physics : MP identifier ! ! me : print control flag ! ! ! ! outputs: (none) ! @@ -358,10 +360,13 @@ subroutine cld_init & ! icldflg : cloud optical property scheme control flag ! ! =0: model use diagnostic cloud method ! ! =1: model use prognostic cloud method ! -! icmphys : cloud microphysics scheme control flag ! -! =1: zhao/carr/sundqvist microphysics cloud ! -! =2: ferrier microphysics cloud scheme ! -! =3: zhao/carr/sundqvist microphysics cloud +pdfcld! +! imp_physics : cloud microphysics scheme control flag ! +! =99: zhao/carr/sundqvist microphysics cloud ! +! =98: zhao/carr/sundqvist microphysics cloud+pdfcld! +! =11: GFDL microphysics cloud ! +! =8: Thompson microphysics ! +! =6: WSM6 microphysics ! +! =10: MG microphysics ! ! iovrsw/iovrlw : sw/lw control flag for cloud overlapping scheme ! ! =0: random overlapping clouds ! ! =1: max/ran overlapping clouds ! @@ -377,7 +382,7 @@ subroutine cld_init & implicit none ! --- inputs: - integer, intent(in) :: NLAY, me + integer, intent(in) :: NLAY, me, imp_physics real (kind=kind_phys), intent(in) :: si(:) @@ -411,15 +416,21 @@ subroutine cld_init & else if (me == 0) then print *,' - Using Prognostic Cloud Method' - if (icmphys == 1) then + if (imp_physics == 99) then print *,' --- Zhao/Carr/Sundqvist microphysics' - elseif (icmphys == 2) then - print *,' --- Ferrier cloud microphysics' - elseif (icmphys == 3) then + elseif (imp_physics == 98) then print *,' --- zhao/carr/sundqvist + pdf cloud' + elseif (imp_physics == 11) then + print *,' --- GFDL Lin cloud microphysics' + elseif (imp_physics == 8) then + print *,' --- Thompson cloud microphysics' + elseif (imp_physics == 6) then + print *,' --- WSM6 cloud microphysics' + elseif (imp_physics == 10) then + print *,' --- MG cloud microphysics' else print *,' !!! ERROR in cloud microphysc specification!!!', & - & ' icmphys (NP3D) =',icmphys + & ' imp_physics (NP3D) =',imp_physics stop endif endif @@ -489,6 +500,7 @@ subroutine progcld1 & & ( plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw, & ! --- inputs: & xlat,xlon,slmsk, IX, NLAY, NLP1, & & uni_cld, lmfshal, lmfdeep2, cldcov, & + & effrl,effri,effrr,effrs,effr_in, & & clouds,clds,mtop,mbot & ! --- outputs: & ) @@ -571,10 +583,11 @@ subroutine progcld1 & ! --- inputs integer, intent(in) :: IX, NLAY, NLP1 - logical, intent(in) :: uni_cld, lmfshal, lmfdeep2 + logical, intent(in) :: uni_cld, lmfshal, lmfdeep2, effr_in real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, & - & tlyr, tvly, qlyr, qstl, rhly, clw, cldcov + & tlyr, tvly, qlyr, qstl, rhly, clw, cldcov, & + & effrl, effri, effrr, effrs real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, & & slmsk @@ -613,22 +626,41 @@ subroutine progcld1 & enddo ! clouds(:,:,:) = 0.0 - do k = 1, NLAY - do i = 1, IX - cldtot(i,k) = 0.0 - cldcnv(i,k) = 0.0 - cwp (i,k) = 0.0 - cip (i,k) = 0.0 - crp (i,k) = 0.0 - csp (i,k) = 0.0 - rew (i,k) = reliq_def ! default liq radius to 10 micron - rei (i,k) = reice_def ! default ice radius to 50 micron - rer (i,k) = rrain_def ! default rain radius to 1000 micron - res (i,k) = rsnow_def ! default snow radius to 250 micron - tem2d (i,k) = min( 1.0, max( 0.0, (con_ttp-tlyr(i,k))*0.05 ) ) - clwf(i,k) = 0.0 + if(effr_in) then + do k = 1, NLAY + do i = 1, IX + cldtot(i,k) = 0.0 + cldcnv(i,k) = 0.0 + cwp (i,k) = 0.0 + cip (i,k) = 0.0 + crp (i,k) = 0.0 + csp (i,k) = 0.0 + rew (i,k) = effrl (i,k) + rei (i,k) = effri (i,k) + rer (i,k) = effrr (i,k) + res (i,k) = effrs (i,k) + tem2d (i,k) = min(1.0, max(0.0,(con_ttp-tlyr(i,k))*0.05)) + clwf(i,k) = 0.0 + enddo enddo - enddo + else + do k = 1, NLAY + do i = 1, IX + cldtot(i,k) = 0.0 + cldcnv(i,k) = 0.0 + cwp (i,k) = 0.0 + cip (i,k) = 0.0 + crp (i,k) = 0.0 + csp (i,k) = 0.0 + rew (i,k) = reliq_def ! default liq radius to 10 micron + rei (i,k) = reice_def ! default ice radius to 50 micron + rer (i,k) = rrain_def ! default rain radius to 1000 micron + res (i,k) = rsnow_def ! default snow radius to 250 micron + tem2d (i,k) = min(1.0, max(0.0, (con_ttp-tlyr(i,k))*0.05)) + clwf(i,k) = 0.0 + enddo + enddo + endif ! if ( lcrick ) then do i = 1, IX @@ -687,13 +719,15 @@ subroutine progcld1 & !> -# Compute effective liquid cloud droplet radius over land. - do i = 1, IX - if (nint(slmsk(i)) == 1) then - do k = 1, NLAY - rew(i,k) = 5.0 + 5.0 * tem2d(i,k) - enddo - endif - enddo + if(.not. effr_in) then + do i = 1, IX + if (nint(slmsk(i)) == 1) then + do k = 1, NLAY + rew(i,k) = 5.0 + 5.0 * tem2d(i,k) + enddo + endif + enddo + endif if (uni_cld) then ! use unified sgs clouds generated outside do k = 1, NLAY @@ -837,32 +871,34 @@ subroutine progcld1 & enddo endif -!> -# Compute effective ice cloud droplet radius following -!! \cite heymsfield_and_mcfarquhar_1996. - - do k = 1, NLAY - do i = 1, IX - tem2 = tlyr(i,k) - con_ttp - - if (cip(i,k) > 0.0) then - tem3 = gord * cip(i,k) * plyr(i,k) / (delp(i,k)*tvly(i,k)) +!> -# Compute effective ice cloud droplet radius following Heymsfield +!! and McFarquhar (1996) \cite heymsfield_and_mcfarquhar_1996. - if (tem2 < -50.0) then - rei(i,k) = (1250.0/9.917) * tem3 ** 0.109 - elseif (tem2 < -40.0) then - rei(i,k) = (1250.0/9.337) * tem3 ** 0.08 - elseif (tem2 < -30.0) then - rei(i,k) = (1250.0/9.208) * tem3 ** 0.055 - else - rei(i,k) = (1250.0/9.387) * tem3 ** 0.031 + if(.not.effr_in) then + do k = 1, NLAY + do i = 1, IX + tem2 = tlyr(i,k) - con_ttp + + if (cip(i,k) > 0.0) then + tem3 = gord * cip(i,k) * plyr(i,k) / (delp(i,k)*tvly(i,k)) + + if (tem2 < -50.0) then + rei(i,k) = (1250.0/9.917) * tem3 ** 0.109 + elseif (tem2 < -40.0) then + rei(i,k) = (1250.0/9.337) * tem3 ** 0.08 + elseif (tem2 < -30.0) then + rei(i,k) = (1250.0/9.208) * tem3 ** 0.055 + else + rei(i,k) = (1250.0/9.387) * tem3 ** 0.031 + endif +! rei(i,k) = max(20.0, min(rei(i,k), 300.0)) +! rei(i,k) = max(10.0, min(rei(i,k), 100.0)) + rei(i,k) = max(10.0, min(rei(i,k), 150.0)) +! rei(i,k) = max(5.0, min(rei(i,k), 130.0)) endif -! rei(i,k) = max(20.0, min(rei(i,k), 300.0)) -! rei(i,k) = max(10.0, min(rei(i,k), 100.0)) - rei(i,k) = max(10.0, min(rei(i,k), 150.0)) -! rei(i,k) = max(5.0, min(rei(i,k), 130.0)) - endif + enddo enddo - enddo + endif ! do k = 1, NLAY @@ -875,7 +911,7 @@ subroutine progcld1 & ! clouds(i,k,6) = 0.0 clouds(i,k,7) = rer(i,k) ! clouds(i,k,8) = 0.0 - clouds(i,k,9) = rei(i,k) + clouds(i,k,9) = res(i,k) enddo enddo @@ -1377,7 +1413,7 @@ subroutine progcld2 & clouds(i,k,7) = rer(i,k) ! clouds(i,k,8) = csp(i,k) !ncar scheme clouds(i,k,8) = csp(i,k) * rsden(i,k) !fu's scheme - clouds(i,k,9) = rei(i,k) + clouds(i,k,9) = res(i,k) enddo enddo @@ -1769,8 +1805,8 @@ subroutine progcld3 & enddo endif -!> -# Calculate effective ice cloud droplet radius following -!! \cite heymsfield_and_mcfarquhar_1996. +!> -# Calculate effective ice cloud droplet radius following Heymsfield +!! and McFarquhar (1996) \cite heymsfield_and_mcfarquhar_1996. do k = 1, nlay do i = 1, ix @@ -1835,50 +1871,24 @@ subroutine progcld3 & !................................... end subroutine progcld3 !----------------------------------- -!> @} -!> \ingroup module_radiation_clouds -!> This subroutine computes cloud related quantities using -!! zhao/moorthi's prognostic cloud microphysics scheme. -!!\param plyr (IX,NLAY), model layer mean pressure in mb (100Pa) -!!\param plvl (IX,NLP1), model level pressure in mb (100Pa) -!!\param tlyr (IX,NLAY), model layer mean temperature in K -!!\param tvly (IX,NLAY), model layer virtual temperature in K -!!\param qlyr (IX,NLAY), layer specific humidity in gm/gm -!!\param clw (IX,NLAY), layer cloud liquid water amount -!!\param ciw (IX,NLAY), layer cloud ice water amount -!!\param xlat (IX), grid latitude in radians, default to pi/2 -> -!! -pi/2 range, otherwise see in-line comment -!!\param xlon (IX), grid longitude in radians (not used) -!!\param slmsk (IX), sea/land mask array (sea:0,land:1,sea-ice:2) -!!\param IX horizontal dimention -!!\param NLAY,NLP1 vertical layer/level dimensions -!!\param clouds (IX,NLAY,NF_CLDS), cloud profiles -!!\n (:,:,1) - layer total cloud fraction -!!\n (:,:,2) - layer cloud liq water path \f$(g/m^2)\f$ -!!\n (:,:,3) - mean eff radius for liq cloud (micron) -!!\n (:,:,4) - layer cloud ice water path \f$(g/m^2)\f$ -!!\n (:,:,5) - mean eff radius for ice cloud (micron) -!!\n (:,:,6) - layer rain drop water path (not assigned) -!!\n (:,:,7) - mean eff radius for rain drop (micron) -!!\n (:,:,8) - layer snow flake water path (not assigned) -!!\n (:,:,9) - mean eff radius for snow flake (micron) -!!\n *** fu's scheme need to be normalized by snow density \f$ (g/m^3/1.0e6)\f$ -!!\param clds (IX,5), fraction of clouds for low, mid, hi, tot, bl -!!\param mtop (IX,3), vertical indices for low, mid, hi cloud tops -!!\param mbot (IX,3), vertical indices for low, mid, hi cloud bases -!>\section gen_progclduni progclduni General Algorithm -!> @{ - subroutine progclduni & - & ( plyr,plvl,tlyr,tvly,clw,ciw, & ! --- inputs: - & xlat,xlon,slmsk, IX, NLAY, NLP1, cldcov, & - & clouds,clds,mtop,mbot & ! --- outputs: + +!----------------------------------- + subroutine progcld4 & +!................................... + +! --- inputs: + & ( plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw,cnvw,cnvc, & + & xlat,xlon,slmsk,cldtot, & + & IX, NLAY, NLP1, & +! --- outputs: + & clouds,clds,mtop,mbot & & ) ! ================= subprogram documentation block ================ ! ! ! -! subprogram: progclduni computes cloud related quantities using ! -! for unified cloud microphysics scheme. ! +! subprogram: progcld4 computes cloud related quantities using ! +! GFDL Lin MP prognostic cloud microphysics scheme. ! ! ! ! abstract: this program computes cloud fractions from cloud ! ! condensates, calculates liquid/ice cloud droplet effective radius, ! @@ -1887,7 +1897,7 @@ subroutine progclduni & ! top and base. the three vertical cloud domains are set up in the ! ! initial subroutine "cld_init". ! ! ! -! usage: call progclduni ! +! usage: call progcld4 ! ! ! ! subprograms called: gethml ! ! ! @@ -1903,8 +1913,12 @@ subroutine progclduni & ! plvl (IX,NLP1) : model level pressure in mb (100Pa) ! ! tlyr (IX,NLAY) : model layer mean temperature in k ! ! tvly (IX,NLAY) : model layer virtual temperature in k ! -! clw (IX,NLAY) : layer cloud liquid water amount ! -! ciw (IX,NLAY) : layer cloud ice water amount ! +! qlyr (IX,NLAY) : layer specific humidity in gm/gm ! +! qstl (IX,NLAY) : layer saturate humidity in gm/gm ! +! rhly (IX,NLAY) : layer relative humidity (=qlyr/qstl) ! +! clw (IX,NLAY) : layer cloud condensate amount ! +! cnvw (IX,NLAY) : layer convective cloud condensate ! +! cnvc (IX,NLAY) : layer convective cloud cover ! ! xlat (IX) : grid latitude in radians, default to pi/2 -> -pi/2! ! range, otherwise see in-line comment ! ! xlon (IX) : grid longitude in radians (not used) ! @@ -1932,8 +1946,7 @@ subroutine progclduni & ! ivflip : control flag of vertical index direction ! ! =0: index from toa to surface ! ! =1: index from surface to toa ! -! lmfshal : mass-flux shallow conv scheme flag ! -! lmfdeep2 : scale-aware mass-flux deep conv scheme flag ! +! lsashal : control flag for shallow convection ! ! lcrick : control flag for eliminating CRICK ! ! =t: apply layer smoothing to eliminate CRICK ! ! =f: do not apply layer smoothing ! @@ -1949,7 +1962,7 @@ subroutine progclduni & integer, intent(in) :: IX, NLAY, NLP1 real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, & - & tlyr, tvly, clw, ciw, cldcov + & tlyr, tvly, qlyr, qstl, rhly, clw, cldtot, cnvw, cnvc real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, & & slmsk @@ -1962,13 +1975,13 @@ subroutine progclduni & integer, dimension(:,:), intent(out) :: mtop,mbot ! --- local variables: - real (kind=kind_phys), dimension(IX,NLAY) :: cldtot, cldcnv, & - & cwp, cip, crp, csp, rew, rei, res, rer, delp, tem2d -! & cwp, cip, crp, csp, rew, rei, res, rer, delp, tem2d, clwf + 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) :: tem1, tem2, tem3 + real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh, value, & + & tem1, tem2, tem3 integer :: i, k, id, nf @@ -1986,7 +1999,6 @@ subroutine progclduni & do k = 1, NLAY do i = 1, IX - cldtot(i,k) = 0.0 cldcnv(i,k) = 0.0 cwp (i,k) = 0.0 cip (i,k) = 0.0 @@ -1997,30 +2009,30 @@ subroutine progclduni & rer (i,k) = rrain_def ! default rain radius to 1000 micron res (i,k) = rsnow_def ! default snow radius to 250 micron tem2d (i,k) = min( 1.0, max( 0.0, (con_ttp-tlyr(i,k))*0.05 ) ) -! clwf(i,k) = 0.0 + clwf(i,k) = 0.0 enddo enddo ! -! if ( lcrick ) then -! do i = 1, IX -! clwf(i,1) = 0.75*clw(i,1) + 0.25*clw(i,2) -! clwf(i,nlay) = 0.75*clw(i,nlay) + 0.25*clw(i,nlay-1) -! enddo -! do k = 2, NLAY-1 -! do i = 1, IX -! clwf(i,K) = 0.25*clw(i,k-1) + 0.5*clw(i,k) + 0.25*clw(i,k+1) -! enddo -! enddo -! else -! do k = 1, NLAY -! do i = 1, IX -! clwf(i,k) = clw(i,k) -! enddo -! enddo -! endif + if ( lcrick ) then + do i = 1, IX + clwf(i,1) = 0.75*clw(i,1) + 0.25*clw(i,2) + clwf(i,nlay) = 0.75*clw(i,nlay) + 0.25*clw(i,nlay-1) + enddo + do k = 2, NLAY-1 + do i = 1, IX + clwf(i,K) = 0.25*clw(i,k-1) + 0.5*clw(i,k) + 0.25*clw(i,k+1) + enddo + enddo + else + do k = 1, NLAY + do i = 1, IX + clwf(i,k) = clw(i,k) + enddo + enddo + endif -!> -# Find top pressure for each cloud domain for given latitude. -! ptopc(k,i): top presure of each cld domain (k=1-4 are sfc,L,m,h; +! --- find top pressure for each cloud domain for given latitude +! 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 id = 1, 4 @@ -2034,29 +2046,29 @@ subroutine progclduni & enddo enddo -!> -# Compute cloud liquid/ice condensate path in \f$ g/m^2 \f$ . +! --- compute liquid/ice condensate path in g/m**2 if ( ivflip == 0 ) then ! input data from toa to sfc do k = 1, NLAY do i = 1, IX delp(i,k) = plvl(i,k+1) - plvl(i,k) - tem1 = gfac * delp(i,k) - cip(i,k) = ciw(i,k) * tem1 - cwp(i,k) = clw(i,k) * tem1 + clwt = max(0.0,(clwf(i,k)+cnvw(i,k))) * gfac * delp(i,k) + cip(i,k) = clwt * tem2d(i,k) + cwp(i,k) = clwt - cip(i,k) enddo enddo else ! input data from sfc to toa do k = 1, NLAY do i = 1, IX delp(i,k) = plvl(i,k) - plvl(i,k+1) - tem1 = gfac * delp(i,k) - cip(i,k) = ciw(i,k) * tem1 - cwp(i,k) = clw(i,k) * tem1 + clwt = max(0.0,(clwf(i,k)+cnvw(i,k))) * gfac * delp(i,k) + cip(i,k) = clwt * tem2d(i,k) + cwp(i,k) = clwt - cip(i,k) enddo enddo endif ! end_if_ivflip -!> -# Compute effective liquid cloud droplet radius over land. +! --- effective liquid cloud droplet radius over land do i = 1, IX if (nint(slmsk(i)) == 1) then @@ -2065,16 +2077,10 @@ subroutine progclduni & enddo endif enddo - do k = 1, NLAY - do i = 1, IX - cldtot(i,k) = cldcov(i,k) - enddo - enddo do k = 1, NLAY do i = 1, IX if (cldtot(i,k) < climit) then - cldtot(i,k) = 0.0 cwp(i,k) = 0.0 cip(i,k) = 0.0 crp(i,k) = 0.0 @@ -2097,8 +2103,7 @@ subroutine progclduni & enddo endif -!> -# Compute effective ice cloud droplet radius following -!! \cite heymsfield_and_mcfarquhar_1996. +! --- effective ice cloud droplet radius do k = 1, NLAY do i = 1, IX @@ -2135,7 +2140,1070 @@ subroutine progclduni & ! clouds(i,k,6) = 0.0 clouds(i,k,7) = rer(i,k) ! clouds(i,k,8) = 0.0 - clouds(i,k,9) = rei(i,k) + clouds(i,k,9) = res(i,k) + enddo + enddo + + +! --- compute low, mid, high, total, and boundary layer cloud fractions +! and clouds top/bottom layer indices for low, mid, and high clouds. +! The three cloud domain boundaries are defined by ptopc. The cloud +! overlapping method is defined by control flag 'iovr', which may +! be different for lw and sw radiation programs. + + call gethml & +! --- inputs: + & ( plyr, ptop1, cldtot, cldcnv, & + & IX,NLAY, & +! --- outputs: + & clds, mtop, mbot & + & ) + + +! + return +!................................... + end subroutine progcld4 +!----------------------------------- + +!----------------------------------- + subroutine progcld4o & +!................................... + +! --- inputs: + & ( plyr,plvl,tlyr,tvly,qlyr,qstl,rhly,clw, & + & xlat,xlon,slmsk, & + & ntrac,ntcw,ntiw,ntrw,ntsw,ntgl,ntclamt, & + & IX, NLAY, NLP1, & +! --- outputs: + & clouds,clds,mtop,mbot & + & ) + +! ================= subprogram documentation block ================ ! +! ! +! subprogram: progcld4o computes cloud related quantities using ! +! GFDL Lin MP prognostic cloud microphysics scheme. Moist species ! +! from MP are fed into the corresponding arrays for calcuation of ! +! ! +! abstract: this program computes cloud fractions from cloud ! +! condensates, calculates liquid/ice cloud droplet effective radius, ! +! and computes the low, mid, high, total and boundary layer cloud ! +! fractions and the vertical indices of low, mid, and high cloud ! +! top and base. the three vertical cloud domains are set up in the ! +! initial subroutine "cld_init". ! +! ! +! usage: call progcld4o ! +! ! +! subprograms called: gethml ! +! ! +! attributes: ! +! language: fortran 90 ! +! machine: ibm-sp, sgi ! +! ! +! ! +! ==================== definition of variables ==================== ! +! ! +! input variables: ! +! plyr (IX,NLAY) : model layer mean pressure in mb (100Pa) ! +! plvl (IX,NLP1) : model level pressure in mb (100Pa) ! +! tlyr (IX,NLAY) : model layer mean temperature in k ! +! tvly (IX,NLAY) : model layer virtual temperature in k ! +! qlyr (IX,NLAY) : layer specific humidity in gm/gm ! +! qstl (IX,NLAY) : layer saturate humidity in gm/gm ! +! rhly (IX,NLAY) : layer relative humidity (=qlyr/qstl) ! +! clw (IX,NLAY,NTRAC) : layer cloud condensate amount ! +! xlat (IX) : grid latitude in radians, default to pi/2 -> -pi/2! +! 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) ! +! IX : horizontal dimention ! +! NLAY,NLP1 : vertical layer/level dimensions ! +! ! +! output variables: ! +! clouds(IX,NLAY,NF_CLDS) : cloud profiles ! +! clouds(:,:,1) - layer total cloud fraction ! +! clouds(:,:,2) - layer cloud liq water path (g/m**2) ! +! clouds(:,:,3) - mean eff radius for liq cloud (micron) ! +! clouds(:,:,4) - layer cloud ice water path (g/m**2) ! +! clouds(:,:,5) - mean eff radius for ice cloud (micron) ! +! clouds(:,:,6) - layer rain drop water path not assigned ! +! clouds(:,:,7) - mean eff radius for rain drop (micron) ! +! *** clouds(:,:,8) - layer snow flake water path not assigned ! +! clouds(:,:,9) - mean eff radius for snow flake (micron) ! +! *** fu's scheme need to be normalized by snow density (g/m**3/1.0e6) ! +! 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 ! +! ! +! module variables: ! +! ivflip : control flag of vertical index direction ! +! =0: index from toa to surface ! +! =1: index from surface to toa ! +! lsashal : control flag for shallow convection ! +! lcrick : control flag for eliminating CRICK ! +! =t: apply layer smoothing to eliminate CRICK ! +! =f: do not apply layer smoothing ! +! lcnorm : control flag for in-cld condensate ! +! =t: normalize cloud condensate ! +! =f: not normalize cloud condensate ! +! ! +! ==================== end of description ===================== ! +! + implicit none + +! --- inputs + integer, intent(in) :: IX, NLAY, NLP1 + integer, intent(in) :: ntrac, ntcw, ntiw, ntrw, ntsw, ntgl, & + & ntclamt + + real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, & + & tlyr, tvly, qlyr, qstl, rhly + + + real (kind=kind_phys), dimension(:,:,:), intent(in) :: clw + real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, & + & slmsk + +! --- outputs + real (kind=kind_phys), dimension(:,:,:), intent(out) :: clouds + + real (kind=kind_phys), dimension(:,:), intent(out) :: clds + + 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 + + real (kind=kind_phys) :: ptop1(IX,NK_CLDS+1) + + real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh, value, & + & tem1, tem2, tem3 + real (kind=kind_phys), dimension(IX,NLAY) :: cldtot + + integer :: i, k, id, nf + +! +!===> ... begin here +! + do nf=1,nf_clds + do k=1,nlay + do i=1,ix + clouds(i,k,nf) = 0.0 + enddo + enddo + enddo +! clouds(:,:,:) = 0.0 + + do k = 1, NLAY + do i = 1, IX + cldcnv(i,k) = 0.0 + cwp (i,k) = 0.0 + cip (i,k) = 0.0 + crp (i,k) = 0.0 + csp (i,k) = 0.0 + rew (i,k) = reliq_def ! default liq radius to 10 micron + rei (i,k) = reice_def ! default ice radius to 50 micron + rer (i,k) = rrain_def ! default rain radius to 1000 micron + res (i,k) = rsnow_def ! default snow radius to 250 micron + tem2d (i,k) = min( 1.0, max( 0.0, (con_ttp-tlyr(i,k))*0.05 ) ) + cldtot(i,k) = clw(i,k,ntclamt) + enddo + enddo + +! --- find top pressure for each cloud domain for given latitude +! 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 id = 1, 4 + tem1 = ptopc(id,2) - ptopc(id,1) + + do i =1, IX + tem2 = xlat(i) / con_pi ! if xlat in pi/2 -> -pi/2 range +! tem2 = 0.5 - xlat(i)/con_pi ! if xlat in 0 -> pi range + + ptop1(i,id) = ptopc(id,1) + tem1*max( 0.0, 4.0*abs(tem2)-1.0 ) + enddo + enddo + +! --- compute liquid/ice condensate path in g/m**2 + + if ( ivflip == 0 ) then ! input data from toa to sfc + do k = 1, NLAY + do i = 1, IX + delp(i,k) = plvl(i,k+1) - plvl(i,k) + cwp(i,k) = max(0.0, clw(i,k,ntcw) * gfac * delp(i,k)) + cip(i,k) = max(0.0, clw(i,k,ntiw) * gfac * delp(i,k)) + crp(i,k) = max(0.0, clw(i,k,ntrw) * gfac * delp(i,k)) + csp(i,k) = max(0.0, (clw(i,k,ntsw)+clw(i,k,ntgl)) * & + & gfac * delp(i,k)) + enddo + enddo + else ! input data from sfc to toa + do k = 1, NLAY + do i = 1, IX + delp(i,k) = plvl(i,k) - plvl(i,k+1) + cwp(i,k) = max(0.0, clw(i,k,ntcw) * gfac * delp(i,k)) + cip(i,k) = max(0.0, clw(i,k,ntiw) * gfac * delp(i,k)) + crp(i,k) = max(0.0, clw(i,k,ntrw) * gfac * delp(i,k)) + csp(i,k) = max(0.0, (clw(i,k,ntsw)+clw(i,k,ntgl)) * & + & gfac * delp(i,k)) + enddo + enddo + endif ! end_if_ivflip + +! --- effective liquid cloud droplet radius over land + + do i = 1, IX + if (nint(slmsk(i)) == 1) then + do k = 1, NLAY + rew(i,k) = 5.0 + 5.0 * tem2d(i,k) + enddo + endif + enddo + + do k = 1, NLAY + do i = 1, IX + if (cldtot(i,k) < climit) then + cwp(i,k) = 0.0 + cip(i,k) = 0.0 + crp(i,k) = 0.0 + csp(i,k) = 0.0 + endif + enddo + enddo + + if ( lcnorm ) then + do k = 1, NLAY + do i = 1, IX + if (cldtot(i,k) >= climit) then + tem1 = 1.0 / max(climit2, cldtot(i,k)) + cwp(i,k) = cwp(i,k) * tem1 + cip(i,k) = cip(i,k) * tem1 + crp(i,k) = crp(i,k) * tem1 + csp(i,k) = csp(i,k) * tem1 + endif + enddo + enddo + endif + +! --- effective ice cloud droplet radius + + do k = 1, NLAY + do i = 1, IX + tem2 = tlyr(i,k) - con_ttp + + if (cip(i,k) > 0.0) then + tem3 = gord * cip(i,k) * plyr(i,k) / (delp(i,k)*tvly(i,k)) + + if (tem2 < -50.0) then + rei(i,k) = (1250.0/9.917) * tem3 ** 0.109 + elseif (tem2 < -40.0) then + rei(i,k) = (1250.0/9.337) * tem3 ** 0.08 + elseif (tem2 < -30.0) then + rei(i,k) = (1250.0/9.208) * tem3 ** 0.055 + else + rei(i,k) = (1250.0/9.387) * tem3 ** 0.031 + endif +! rei(i,k) = max(20.0, min(rei(i,k), 300.0)) +! rei(i,k) = max(10.0, min(rei(i,k), 100.0)) + rei(i,k) = max(10.0, min(rei(i,k), 150.0)) +! rei(i,k) = max(5.0, min(rei(i,k), 130.0)) + endif + enddo + enddo + +! + do k = 1, NLAY + do i = 1, IX + clouds(i,k,1) = cldtot(i,k) + clouds(i,k,2) = cwp(i,k) + clouds(i,k,3) = rew(i,k) + clouds(i,k,4) = cip(i,k) + clouds(i,k,5) = rei(i,k) + clouds(i,k,6) = crp(i,k) + clouds(i,k,7) = rer(i,k) + clouds(i,k,8) = csp(i,k) + clouds(i,k,9) = rei(i,k) + enddo + enddo + + +! --- compute low, mid, high, total, and boundary layer cloud fractions +! and clouds top/bottom layer indices for low, mid, and high clouds. +! The three cloud domain boundaries are defined by ptopc. The cloud +! overlapping method is defined by control flag 'iovr', which may +! be different for lw and sw radiation programs. + + call gethml & +! --- inputs: + & ( plyr, ptop1, cldtot, cldcnv, & + & IX,NLAY, & +! --- outputs: + & clds, mtop, mbot & + & ) + + +! + return +!................................... + end subroutine progcld4o +!----------------------------------- + +!----------------------------------- + + subroutine progcld5 & + & ( plyr,plvl,tlyr,qlyr,qstl,rhly,clw, & ! --- inputs: + & xlat,xlon,slmsk, & + & ntrac,ntcw,ntiw,ntrw,ntsw,ntgl, & + & IX, NLAY, NLP1, & + & uni_cld, lmfshal, lmfdeep2, cldcov, & + & re_cloud,re_ice,re_snow, & + & clouds,clds,mtop,mbot & ! --- outputs: + & ) + +! ================= subprogram documentation block ================ ! +! ! +! subprogram: progcld5 computes cloud related quantities using ! +! Thompson/WSM6 cloud microphysics scheme. ! +! ! +! abstract: this program computes cloud fractions from cloud ! +! condensates, ! +! and computes the low, mid, high, total and boundary layer cloud ! +! fractions and the vertical indices of low, mid, and high cloud ! +! top and base. the three vertical cloud domains are set up in the ! +! initial subroutine "cld_init". ! +! ! +! usage: call progcld5 ! +! ! +! subprograms called: gethml ! +! ! +! attributes: ! +! language: fortran 90 ! +! machine: ibm-sp, sgi ! +! ! +! ! +! ==================== definition of variables ==================== ! +! ! +! input variables: ! +! plyr (IX,NLAY) : model layer mean pressure in mb (100Pa) ! +! plvl (IX,NLP1) : model level pressure in mb (100Pa) ! +! tlyr (IX,NLAY) : model layer mean temperature in k ! +! tvly (IX,NLAY) : model layer virtual temperature in k ! +! qlyr (IX,NLAY) : layer specific humidity in gm/gm ! +! qstl (IX,NLAY) : layer saturate humidity in gm/gm ! +! rhly (IX,NLAY) : layer relative humidity (=qlyr/qstl) ! +! clw (IX,NLAY,ntrac) : layer cloud condensate amount ! +! xlat (IX) : grid latitude in radians, default to pi/2 -> -pi/2! +! 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) ! +! IX : horizontal dimention ! +! NLAY,NLP1 : vertical layer/level dimensions ! +! uni_cld : logical - true for cloud fraction from shoc ! +! lmfshal : logical - true for mass flux shallow convection ! +! lmfdeep2 : logical - true for mass flux deep convection ! +! cldcov : layer cloud fraction (used when uni_cld=.true. ! +! ! +! output variables: ! +! clouds(IX,NLAY,NF_CLDS) : cloud profiles ! +! clouds(:,:,1) - layer total cloud fraction ! +! clouds(:,:,2) - layer cloud liq water path (g/m**2) ! +! clouds(:,:,3) - mean eff radius for liq cloud (micron) ! +! clouds(:,:,4) - layer cloud ice water path (g/m**2) ! +! clouds(:,:,5) - mean eff radius for ice cloud (micron) ! +! clouds(:,:,6) - layer rain drop water path not assigned ! +! clouds(:,:,7) - mean eff radius for rain drop (micron) ! +! *** clouds(:,:,8) - layer snow flake water path not assigned ! +! clouds(:,:,9) - mean eff radius for snow flake (micron) ! +! *** fu's scheme need to be normalized by snow density (g/m**3/1.0e6) ! +! 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 ! +! ! +! module variables: ! +! ivflip : control flag of vertical index direction ! +! =0: index from toa to surface ! +! =1: index from surface to toa ! +! lmfshal : mass-flux shallow conv scheme flag ! +! lmfdeep2 : scale-aware mass-flux deep conv scheme flag ! +! lcrick : control flag for eliminating CRICK ! +! =t: apply layer smoothing to eliminate CRICK ! +! =f: do not apply layer smoothing ! +! lcnorm : control flag for in-cld condensate ! +! =t: normalize cloud condensate ! +! =f: not normalize cloud condensate ! +! ! +! ==================== end of description ===================== ! +! + implicit none + +! --- inputs + integer, intent(in) :: IX, NLAY, NLP1 + integer, intent(in) :: ntrac, ntcw, ntiw, ntrw, ntsw, ntgl + + logical, intent(in) :: uni_cld, lmfshal, lmfdeep2 + + real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, & + & tlyr, qlyr, qstl, rhly, cldcov, + & re_cloud, re_ice, re_snow + + real (kind=kind_phys), dimension(:,:,:), intent(in) :: clw + + real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, & + & slmsk + +! --- outputs + real (kind=kind_phys), dimension(:,:,:), intent(out) :: clouds + + real (kind=kind_phys), dimension(:,:), intent(out) :: clds + + integer, dimension(:,:), intent(out) :: mtop,mbot + +! --- local variables: + real (kind=kind_phys), dimension(IX,NLAY) :: cldtot, 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) :: clwmin, clwm, clwt, onemrh, value, & + & tem1, tem2, tem3 + + integer :: i, k, id, nf + +! --- constant values +! real (kind=kind_phys), parameter :: xrc3 = 200. + real (kind=kind_phys), parameter :: xrc3 = 100. + +! +!===> ... begin here +! + do nf=1,nf_clds + do k=1,nlay + do i=1,ix + clouds(i,k,nf) = 0.0 + enddo + enddo + enddo +! clouds(:,:,:) = 0.0 + + do k = 1, NLAY + do i = 1, IX + cldtot(i,k) = 0.0 + cldcnv(i,k) = 0.0 + cwp (i,k) = 0.0 + cip (i,k) = 0.0 + crp (i,k) = 0.0 + csp (i,k) = 0.0 + rew (i,k) = re_cloud(i,k) + rei (i,k) = re_ice(i,k) + rer (i,k) = rrain_def ! default rain radius to 1000 micron + res (i,k) = re_snow(i,K) +! tem2d (i,k) = min( 1.0, max( 0.0, (con_ttp-tlyr(i,k))*0.05 ) ) + clwf(i,k) = 0.0 + enddo + enddo +! +! if ( lcrick ) then +! do i = 1, IX +! clwf(i,1) = 0.75*clw(i,1) + 0.25*clw(i,2) +! clwf(i,nlay) = 0.75*clw(i,nlay) + 0.25*clw(i,nlay-1) +! enddo +! do k = 2, NLAY-1 +! do i = 1, IX +! clwf(i,K) = 0.25*clw(i,k-1) + 0.5*clw(i,k) + 0.25*clw(i,k+1) +! enddo +! enddo +! else +! do k = 1, NLAY +! do i = 1, IX +! clwf(i,k) = clw(i,k) +! enddo +! enddo +! endif + + do k = 1, NLAY + do i = 1, IX + clwf(i,k) = clw(i,k,ntcw) + clw(i,k,ntiw) + clw(i,k,ntsw) + enddo + enddo +!> -# Find top pressure for each cloud domain for given latitude. +! 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 id = 1, 4 + tem1 = ptopc(id,2) - ptopc(id,1) + + do i =1, IX + tem2 = xlat(i) / con_pi ! if xlat in pi/2 -> -pi/2 range +! tem2 = 0.5 - xlat(i)/con_pi ! if xlat in 0 -> pi range + + ptop1(i,id) = ptopc(id,1) + tem1*max( 0.0, 4.0*abs(tem2)-1.0 ) + enddo + enddo + +!> -# Compute cloud liquid/ice condensate path in \f$ g/m^2 \f$ . + + if ( ivflip == 0 ) then ! input data from toa to sfc + do k = 1, NLAY + do i = 1, IX + delp(i,k) = plvl(i,k+1) - plvl(i,k) + cwp(i,k) = max(0.0, clw(i,k,ntcw) * gfac * delp(i,k)) + cip(i,k) = max(0.0, clw(i,k,ntiw) * gfac * delp(i,k)) + crp(i,k) = max(0.0, clw(i,k,ntrw) * gfac * delp(i,k)) + csp(i,k) = max(0.0, (clw(i,k,ntsw)+clw(i,k,ntgl)) * & + & gfac * delp(i,k)) + enddo + enddo + else ! input data from sfc to toa + do k = 1, NLAY + do i = 1, IX + delp(i,k) = plvl(i,k) - plvl(i,k+1) + cwp(i,k) = max(0.0, clw(i,k,ntcw) * gfac * delp(i,k)) + cip(i,k) = max(0.0, clw(i,k,ntiw) * gfac * delp(i,k)) + crp(i,k) = max(0.0, clw(i,k,ntrw) * gfac * delp(i,k)) + csp(i,k) = max(0.0, (clw(i,k,ntsw)+clw(i,k,ntgl)) * & + & gfac * delp(i,k)) + enddo + enddo + endif ! end_if_ivflip + + if (uni_cld) then ! use unified sgs clouds generated outside + do k = 1, NLAY + do i = 1, IX + cldtot(i,k) = cldcov(i,k) + enddo + enddo + + else + +!> -# Calculate layer cloud fraction. + + if ( ivflip == 0 ) then ! input data from toa to sfc + + clwmin = 0.0 + if (.not. lmfshal) then + do k = NLAY, 1, -1 + do i = 1, IX + clwt = 1.0e-6 * (plyr(i,k)*0.001) +! clwt = 2.0e-6 * (plyr(i,k)*0.001) + + if (clwf(i,k) > clwt) then + + onemrh= max( 1.e-10, 1.0-rhly(i,k) ) + clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) + + tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0) + tem1 = 2000.0 / tem1 +! tem1 = 1000.0 / tem1 + + value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) + tem2 = sqrt( sqrt(rhly(i,k)) ) + + cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 ) + endif + enddo + enddo + else + do k = NLAY, 1, -1 + do i = 1, IX + clwt = 1.0e-6 * (plyr(i,k)*0.001) +! clwt = 2.0e-6 * (plyr(i,k)*0.001) + + if (clwf(i,k) > clwt) then + onemrh= max( 1.e-10, 1.0-rhly(i,k) ) + clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) +! + tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) !jhan + if (lmfdeep2) then + tem1 = xrc3 / tem1 + else + tem1 = 100.0 / tem1 + endif +! + value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) + tem2 = sqrt( sqrt(rhly(i,k)) ) + cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 ) + endif + enddo + enddo + endif + + else ! input data from sfc to toa + + clwmin = 0.0 + if (.not. lmfshal) then + do k = 1, NLAY + do i = 1, IX + clwt = 1.0e-6 * (plyr(i,k)*0.001) +! clwt = 2.0e-6 * (plyr(i,k)*0.001) + + if (clwf(i,k) > clwt) then + + onemrh= max( 1.e-10, 1.0-rhly(i,k) ) + clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) + + tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0) + tem1 = 2000.0 / tem1 + +! tem1 = 1000.0 / tem1 + + value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) + tem2 = sqrt( sqrt(rhly(i,k)) ) + + cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 ) + endif + enddo + enddo + else + do k = 1, NLAY + do i = 1, IX + clwt = 1.0e-6 * (plyr(i,k)*0.001) +! clwt = 2.0e-6 * (plyr(i,k)*0.001) + + if (clwf(i,k) > clwt) then + onemrh= max( 1.e-10, 1.0-rhly(i,k) ) + clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) +! + tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) !jhan + if (lmfdeep2) then + tem1 = xrc3 / tem1 + else + tem1 = 100.0 / tem1 + endif +! + value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) + tem2 = sqrt( sqrt(rhly(i,k)) ) + + cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 ) + endif + enddo + enddo + endif + + endif ! end_if_flip + endif ! if (uni_cld) then + + do k = 1, NLAY + do i = 1, IX + if (cldtot(i,k) < climit) then + cldtot(i,k) = 0.0 + cwp(i,k) = 0.0 + cip(i,k) = 0.0 + crp(i,k) = 0.0 + csp(i,k) = 0.0 + endif + enddo + enddo + + if ( lcnorm ) then + do k = 1, NLAY + do i = 1, IX + if (cldtot(i,k) >= climit) then + tem1 = 1.0 / max(climit2, cldtot(i,k)) + cwp(i,k) = cwp(i,k) * tem1 + cip(i,k) = cip(i,k) * tem1 + crp(i,k) = crp(i,k) * tem1 + csp(i,k) = csp(i,k) * tem1 + endif + enddo + enddo + endif + +! + do k = 1, NLAY + do i = 1, IX + clouds(i,k,1) = cldtot(i,k) + clouds(i,k,2) = cwp(i,k) + clouds(i,k,3) = rew(i,k) + clouds(i,k,4) = cip(i,k) + clouds(i,k,5) = rei(i,k) + clouds(i,k,6) = crp(i,k) ! added for Thompson + clouds(i,k,7) = rer(i,k) + clouds(i,k,8) = csp(i,k) ! added for Thompson + clouds(i,k,9) = rei(i,k) + enddo + enddo + + +!> -# Call gethml() to compute low,mid,high,total, and boundary layer +!! cloud fractions and clouds top/bottom layer indices for low, mid, +!! and high clouds. +! --- compute low, mid, high, total, and boundary layer cloud fractions +! and clouds top/bottom layer indices for low, mid, and high clouds. +! The three cloud domain boundaries are defined by ptopc. The cloud +! overlapping method is defined by control flag 'iovr', which may +! be different for lw and sw radiation programs. + + call gethml & +! --- inputs: + & ( plyr, ptop1, cldtot, cldcnv, & + & IX,NLAY, & +! --- outputs: + & clds, mtop, mbot & + & ) + + +! + return +!................................... + end subroutine progcld5 +!................................... + +!> @} + +!> \ingroup module_radiation_clouds +!> This subroutine computes cloud related quantities using +!! zhao/moorthi's prognostic cloud microphysics scheme. +!!\param plyr (IX,NLAY), model layer mean pressure in mb (100Pa) +!!\param plvl (IX,NLP1), model level pressure in mb (100Pa) +!!\param tlyr (IX,NLAY), model layer mean temperature in K +!!\param tvly (IX,NLAY), model layer virtual temperature in K +!!\param qlyr (IX,NLAY), layer specific humidity in gm/gm +!!\param clw (IX,NLAY), layer cloud liquid water amount +!!\param ciw (IX,NLAY), layer cloud ice water amount +!!\param xlat (IX), grid latitude in radians, default to pi/2 -> +!! -pi/2 range, otherwise see in-line comment +!!\param xlon (IX), grid longitude in radians (not used) +!!\param slmsk (IX), sea/land mask array (sea:0,land:1,sea-ice:2) +!!\param IX horizontal dimention +!!\param NLAY,NLP1 vertical layer/level dimensions +!!\param clouds (IX,NLAY,NF_CLDS), cloud profiles +!!\n (:,:,1) - layer total cloud fraction +!!\n (:,:,2) - layer cloud liq water path \f$(g/m^2)\f$ +!!\n (:,:,3) - mean eff radius for liq cloud (micron) +!!\n (:,:,4) - layer cloud ice water path \f$(g/m^2)\f$ +!!\n (:,:,5) - mean eff radius for ice cloud (micron) +!!\n (:,:,6) - layer rain drop water path (not assigned) +!!\n (:,:,7) - mean eff radius for rain drop (micron) +!!\n (:,:,8) - layer snow flake water path (not assigned) +!!\n (:,:,9) - mean eff radius for snow flake (micron) +!!\n *** fu's scheme need to be normalized by snow density \f$ (g/m^3/1.0e6)\f$ +!!\param clds (IX,5), fraction of clouds for low, mid, hi, tot, bl +!!\param mtop (IX,3), vertical indices for low, mid, hi cloud tops +!!\param mbot (IX,3), vertical indices for low, mid, hi cloud bases +!>\section gen_progclduni progclduni General Algorithm +!> @{ + subroutine progclduni & + & ( plyr,plvl,tlyr,tvly,ccnd,ncnd, & ! --- inputs: + & xlat,xlon,slmsk, IX, NLAY, NLP1, cldtot, & + & effrl,effri,effrr,effrs,effr_in, & + & clouds,clds,mtop,mbot & ! --- outputs: + & ) + +! ================= subprogram documentation block ================ ! +! ! +! subprogram: progclduni computes cloud related quantities using ! +! for unified cloud microphysics scheme. ! +! ! +! abstract: this program computes cloud fractions from cloud ! +! condensates, calculates liquid/ice cloud droplet effective radius, ! +! and computes the low, mid, high, total and boundary layer cloud ! +! fractions and the vertical indices of low, mid, and high cloud ! +! top and base. the three vertical cloud domains are set up in the ! +! initial subroutine "cld_init". ! +! ! +! usage: call progclduni ! +! ! +! subprograms called: gethml ! +! ! +! attributes: ! +! language: fortran 90 ! +! machine: ibm-sp, sgi ! +! ! +! ! +! ==================== definition of variables ==================== ! +! ! +! input variables: ! +! plyr (IX,NLAY) : model layer mean pressure in mb (100Pa) ! +! plvl (IX,NLP1) : model level pressure in mb (100Pa) ! +! tlyr (IX,NLAY) : model layer mean temperature in k ! +! tvly (IX,NLAY) : model layer virtual temperature in k ! +! ccnd (IX,NLAY) : layer cloud condensate amount ! +! ncnd : number of layer cloud condensate types ! +! xlat (IX) : grid latitude in radians, default to pi/2 -> -pi/2! +! 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) ! +! IX : horizontal dimention ! +! NLAY,NLP1 : vertical layer/level dimensions ! +! cldtot : unified cloud fracrion from moist physics ! +! effrl (ix,nlay) : effective radius for liquid water ! +! effri (ix,nlay) : effective radius for ice water ! +! effrr (ix,nlay) : effective radius for rain water ! +! effrs (ix,nlay) : effective radius for snow water ! +! effr_in : logical - if .true. use input effective radii ! +! ! +! output variables: ! +! clouds(IX,NLAY,NF_CLDS) : cloud profiles ! +! clouds(:,:,1) - layer total cloud fraction ! +! clouds(:,:,2) - layer cloud liq water path (g/m**2) ! +! clouds(:,:,3) - mean eff radius for liq cloud (micron) ! +! clouds(:,:,4) - layer cloud ice water path (g/m**2) ! +! clouds(:,:,5) - mean eff radius for ice cloud (micron) ! +! clouds(:,:,6) - layer rain drop water path not assigned ! +! clouds(:,:,7) - mean eff radius for rain drop (micron) ! +! *** clouds(:,:,8) - layer snow flake water path not assigned ! +! clouds(:,:,9) - mean eff radius for snow flake (micron) ! +! *** fu's scheme need to be normalized by snow density (g/m**3/1.0e6) ! +! 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 ! +! ! +! module variables: ! +! ivflip : control flag of vertical index direction ! +! =0: index from toa to surface ! +! =1: index from surface to toa ! +! lmfshal : mass-flux shallow conv scheme flag ! +! lmfdeep2 : scale-aware mass-flux deep conv scheme flag ! +! lcrick : control flag for eliminating CRICK ! +! =t: apply layer smoothing to eliminate CRICK ! +! =f: do not apply layer smoothing ! +! lcnorm : control flag for in-cld condensate ! +! =t: normalize cloud condensate ! +! =f: not normalize cloud condensate ! +! ! +! ==================== end of description ===================== ! +! + implicit none + +! --- inputs + integer, intent(in) :: IX, NLAY, NLP1, ncnd + logical, intent(in) :: effr_in + + real (kind=kind_phys), dimension(:,:,:), intent(in) :: ccnd + real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr,& + & tlyr, tvly, cldtot, effrl, effri, effrr, effrs + + real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, & + & slmsk + +! --- outputs + real (kind=kind_phys), dimension(:,:,:), intent(out) :: clouds + + real (kind=kind_phys), dimension(:,:), intent(out) :: clds + + 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 + real (kind=kind_phys), dimension(IX,NLAY,ncnd) :: cndf + + real (kind=kind_phys) :: ptop1(IX,NK_CLDS+1) + + real (kind=kind_phys) :: tem1, tem2, tem3 + + integer :: i, k, id, nf, n + +! +!===> ... begin here +! + do nf=1,nf_clds + do k=1,nlay + do i=1,ix + clouds(i,k,nf) = 0.0 + enddo + enddo + enddo +! clouds(:,:,:) = 0.0 + + if (effr_in) then + do k = 1, NLAY + do i = 1, IX + cldcnv(i,k) = 0.0 + cwp (i,k) = 0.0 + cip (i,k) = 0.0 + crp (i,k) = 0.0 + csp (i,k) = 0.0 + rew (i,k) = effrl (i,k) + rei (i,k) = effri (i,k) + rer (i,k) = effrr (i,k) + res (i,k) = effrs (i,k) + tem2d (i,k) = min( 1.0, max( 0.0,(con_ttp-tlyr(i,k))*0.05)) + enddo + enddo + else + do k = 1, NLAY + do i = 1, IX + cldcnv(i,k) = 0.0 + cwp (i,k) = 0.0 + cip (i,k) = 0.0 + crp (i,k) = 0.0 + csp (i,k) = 0.0 + rew (i,k) = reliq_def ! default liq radius to 10 micron + rei (i,k) = reice_def ! default ice radius to 50 micron + rer (i,k) = rrain_def ! default rain radius to 1000 micron + res (i,k) = rsnow_def ! default snow radius to 250 micron + tem2d (i,k) = min(1.0, max(0.0, (con_ttp-tlyr(i,k))*0.05)) + enddo + enddo + endif +! + do n=1,ncnd + do k = 1, NLAY + do i = 1, IX + cndf(i,k,n) = ccnd(i,k,n) + enddo + enddo + enddo + if ( lcrick ) then + do n=1,ncnd + do i = 1, IX + cndf(i,1,n) = 0.75*ccnd(i,1,n) + 0.25*ccnd(i,2,n) + cndf(i,nlay,n) = 0.75*ccnd(i,nlay,n) + 0.25*ccnd(i,nlay-1,n) + enddo + do k = 2, NLAY-1 + do i = 1, IX + cndf(i,K,n) = 0.25 * (ccnd(i,k-1,n) + ccnd(i,k+1,n)) & + & + 0.5 * ccnd(i,k,n) + enddo + enddo + enddo + endif + +!> -# Find top pressure for each cloud domain for given latitude. +! 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 id = 1, 4 + tem1 = ptopc(id,2) - ptopc(id,1) + + do i =1, IX + tem2 = xlat(i) / con_pi ! if xlat in pi/2 -> -pi/2 range +! tem2 = 0.5 - xlat(i)/con_pi ! if xlat in 0 -> pi range + + ptop1(i,id) = ptopc(id,1) + tem1*max( 0.0, 4.0*abs(tem2)-1.0 ) + enddo + enddo + +!> -# Compute cloud liquid/ice condensate path in \f$ g/m^2 \f$ . + + if ( ivflip == 0 ) then ! input data from toa to sfc + if (ncnd == 2) then + do k = 1, NLAY + do i = 1, IX + delp(i,k) = plvl(i,k+1) - plvl(i,k) + tem1 = gfac * delp(i,k) + cwp(i,k) = cndf(i,k,1) * tem1 + cip(i,k) = cndf(i,k,2) * tem1 + enddo + enddo + elseif (ncnd == 4 .or. ncnd == 5) then + do k = 1, NLAY + do i = 1, IX + delp(i,k) = plvl(i,k+1) - plvl(i,k) + tem1 = gfac * delp(i,k) + cwp(i,k) = cndf(i,k,1) * tem1 + cip(i,k) = cndf(i,k,2) * tem1 + crp(i,k) = cndf(i,k,3) * tem1 + csp(i,k) = cndf(i,k,4) * tem1 + enddo + enddo + endif + else ! input data from sfc to toa + if (ncnd == 2) then + do k = 1, NLAY + do i = 1, IX + delp(i,k) = plvl(i,k) - plvl(i,k+1) + tem1 = gfac * delp(i,k) + cwp(i,k) = cndf(i,k,1) * tem1 + cip(i,k) = cndf(i,k,2) * tem1 + enddo + enddo + elseif (ncnd == 4 .or. ncnd == 5) then + do k = 1, NLAY + do i = 1, IX + delp(i,k) = plvl(i,k) - plvl(i,k+1) + tem1 = gfac * delp(i,k) + cwp(i,k) = cndf(i,k,1) * tem1 + cip(i,k) = cndf(i,k,2) * tem1 + crp(i,k) = cndf(i,k,3) * tem1 + csp(i,k) = cndf(i,k,4) * tem1 + enddo + enddo + endif + + endif ! end_if_ivflip + +!> -# Compute effective liquid cloud droplet radius over land. + + if(.not. effr_in) then + do i = 1, IX + if (nint(slmsk(i)) == 1) then + do k = 1, NLAY + rew(i,k) = 5.0 + 5.0 * tem2d(i,k) + enddo + endif + enddo + endif + + do k = 1, NLAY + do i = 1, IX + if (cldtot(i,k) < climit) then + cwp(i,k) = 0.0 + cip(i,k) = 0.0 + crp(i,k) = 0.0 + csp(i,k) = 0.0 + endif + enddo + enddo + + if ( lcnorm ) then + do k = 1, NLAY + do i = 1, IX + if (cldtot(i,k) >= climit) then + tem1 = 1.0 / max(climit2, cldtot(i,k)) + cwp(i,k) = cwp(i,k) * tem1 + cip(i,k) = cip(i,k) * tem1 + crp(i,k) = crp(i,k) * tem1 + csp(i,k) = csp(i,k) * tem1 + endif + enddo + enddo + endif + +!> -# Compute effective ice cloud droplet radius following +!! \cite heymsfield_and_mcfarquhar_1996. + + if(.not. effr_in) then + do k = 1, NLAY + do i = 1, IX + tem2 = tlyr(i,k) - con_ttp + + if (cip(i,k) > 0.0) then + tem3 = gord * cip(i,k) * plyr(i,k) / (delp(i,k)*tvly(i,k)) + + if (tem2 < -50.0) then + rei(i,k) = (1250.0/9.917) * tem3 ** 0.109 + elseif (tem2 < -40.0) then + rei(i,k) = (1250.0/9.337) * tem3 ** 0.08 + elseif (tem2 < -30.0) then + rei(i,k) = (1250.0/9.208) * tem3 ** 0.055 + else + rei(i,k) = (1250.0/9.387) * tem3 ** 0.031 + endif +! rei(i,k) = max(20.0, min(rei(i,k), 300.0)) +! rei(i,k) = max(10.0, min(rei(i,k), 100.0)) + rei(i,k) = max(10.0, min(rei(i,k), 150.0)) +! rei(i,k) = max(5.0, min(rei(i,k), 130.0)) + endif + enddo + enddo + endif + +! + do k = 1, NLAY + do i = 1, IX + clouds(i,k,1) = cldtot(i,k) + clouds(i,k,2) = cwp(i,k) + clouds(i,k,3) = rew(i,k) + clouds(i,k,4) = cip(i,k) + clouds(i,k,5) = rei(i,k) + clouds(i,k,6) = crp(i,k) + clouds(i,k,7) = rer(i,k) + clouds(i,k,8) = csp(i,k) + clouds(i,k,9) = res(i,k) enddo enddo diff --git a/physics/radiation_surface.f b/physics/radiation_surface.f index 2cd77ab9b..255c43646 100644 --- a/physics/radiation_surface.f +++ b/physics/radiation_surface.f @@ -112,7 +112,7 @@ module module_radiation_surface ! integer, allocatable :: idxems(:,:) !< global surface emissivity index array integer :: iemslw = 0 !< global surface emissivity control flag set up in 'sfc_init' ! - public sfc_init, setalb, setemis + public sfc_init, setalb, setemis, gmln, cdfbet, ppfbet ! ================= contains @@ -304,6 +304,7 @@ subroutine setalb & & ( slmsk,snowf,sncovr,snoalb,zorlf,coszf,tsknf,tairf,hprif, & ! --- inputs: & alvsf,alnsf,alvwf,alnwf,facsf,facwf,fice,tisfc, & & IMAX, & + & albPpert, pertalb, & ! sfc-perts, mgehne & sfcalb & ! --- outputs: & ) @@ -372,7 +373,8 @@ subroutine setalb & real (kind=kind_phys), dimension(:), intent(in) :: & & slmsk, snowf, zorlf, coszf, tsknf, tairf, hprif, & & alvsf, alnsf, alvwf, alnwf, facsf, facwf, fice, tisfc, & - & sncovr, snoalb + & sncovr, snoalb, albPpert ! sfc-perts, mgehne + real (kind=kind_phys), dimension(5), intent(in) :: pertalb ! sfc-perts, mgehne ! --- outputs real (kind=kind_phys), dimension(IMAX,NF_ALBD), intent(out) :: & @@ -383,11 +385,11 @@ subroutine setalb & real (kind=kind_phys) :: asnvb, asnnb, asnvd, asnnd, asevb & &, asenb, asevd, asend, fsno, fsea, rfcs, rfcw, flnd & &, asnow, argh, hrgh, fsno0, fsno1, flnd0, fsea0, csnow & - &, a1, a2, b1, b2, b3, ab1bm, ab2bm + &, a1, a2, b1, b2, b3, ab1bm, ab2bm, m, s, alpha, beta, albtmp real (kind=kind_phys) ffw, dtgd - integer :: i, k + integer :: i, k, kk, iflag ! !===> ... begin here @@ -609,6 +611,28 @@ subroutine setalb & endif ! end if_ialbflg ! + +! sfc-perts, mgehne *** +! perturb all 4 kinds of surface albedo, sfcalb(:,1:4) + if (pertalb(1)>0.0) then + do i = 1, imax + do kk=1, 4 + ! compute beta distribution parameters for all 4 albedos + m = sfcalb(i,kk) + s = pertalb(1)*m*(1.-m) + alpha = m*m*(1.-m)/(s*s)-m + beta = alpha*(1.-m)/m + ! compute beta distribution value corresponding + ! to the given percentile albPpert to use as new albedo + call ppfbet(albPpert(i),alpha,beta,iflag,albtmp) + sfcalb(i,kk) = albtmp + enddo + enddo ! end_do_i_loop + endif + +! *** sfc-perts, mgehne + + return !................................... end subroutine setalb @@ -803,6 +827,265 @@ end subroutine setemis !! @} !----------------------------------- +! mg, sfc perts **** +! --- subroutines for computing the beta distribution value that --- +! --- matches the percentile from the random pattern --- + + + subroutine ppfbet(pr,p,q,iflag,x) + use machine + implicit none + real(kind=kind_phys), intent(in) :: pr, p, q + real(kind=kind_phys), intent(out) :: x + ! local variables + integer iflag, iter, itmax + real(kind=kind_phys) tol, a, b, fa, fb, fc, cdf, tol1 + real(kind=kind_phys) c, d, e, xm, s, u, v, r, eps + data itmax, eps / 50, 1.0E-12 / + + ! Compute beta distribution value corresponding to the + ! probability and distribution parameters a,b. + ! + ! pr - a probability value in the interval [0,1] + ! p - the first parameter of the beta(p,q) distribution + ! q - the second parameter of the beta(p,q) distribution + ! iflag - erro indicator in output, 0-no errors, 1,2-error flags + ! from subroutine cdfbet, 3- pr<0 or pr>1, 4-p<=0 or + ! q<=0, 5-tol<1.E-8, 6-the cdfs at the endpoints have + ! the same sign and no value of x is defined, 7-maximum + ! iterations exceeded and current value of x returned + + tol = 1.0E-5 + + + iflag = 0 + if (pr.lt.0.0.or.pr.gt.1.) then + iflag = 3 + return + endif + if(min(p,q).le.0.) then + iflag =4 + return + endif + if (tol.lt.1.0E-8) then + iflag = 5 + return + endif + a = 0. + b = 1. + fa = -pr + fb = 1.-pr + if (fb*fa.gt.0.0) then + iflag = 6 + return + endif + + fc = fb + do iter =1,itmax + if (fb*fc.gt.0.) then + c=a + fc=fa + d = b-a + e=d + endif + if (abs(fc).lt.abs(fb)) then + a=b + b=c + c=a + fa=fb + fb=fc + fc=fa + endif + + tol1 = 2.*eps*abs(b)+0.5*tol + xm = 0.5*(c-b) + if (abs(xm).le.tol1.or.fb.eq.0.0) then + x=b + return + endif + if (abs(e).ge.tol1.and.abs(fa).gt.abs(fb)) then + s = fb/fa + if (a.eq.c) then + u = 2.0*xm*s + v = 1.0-s + else + v = fa/fc + r = fb/fc + u = s*(2.0*xm*v*(v-r)-(b-a)*(r-1.0)) + v = (v-1.0)*(r-1.0)*(s-1.0) + endif + if (u.gt.0.0) v = -v + u = abs(u) + if (2.0*u.lt.min(3.0*xm*v-ABS(tol1*v),ABS(e*v))) then + e = d + d = u/v + else + d = xm + e = d + endif + + else + + d=xm + e=d + endif + + a = b + fa = fb + if (abs(d).gt.tol1) then + b = b+d + else + b = b+sign(tol1,xm) + endif + call cdfbet(b,p,q,eps,iflag,cdf) + if (iflag.ne.0) return + fb = cdf-pr + enddo + x = b + + return + end subroutine ppfbet + + subroutine cdfbet(x,p,q,eps,iflag,cdfx) + use machine + + ! Computes the value of the cumulative beta distribution at a + ! single point x, given the distribution parameters p,q. + ! + ! x - value at which the CDF is to be computed + ! p - first parameter of the beta function + ! q - second parameter of the beta function + ! eps - desired absolute accuracy + + implicit none + real(kind=kind_phys), intent(in) :: x, p, q, eps + real(kind=kind_phys), intent(out) :: cdfx + ! local vars + integer iflag, jmax, j + logical LL + real(kind=kind_phys) dp, dq, gamln, yxeps, w, uflo + real(kind=kind_phys) xy, yx, pq, qp, pdfl, u, r, v + real(kind=kind_phys) tmp + data jmax, w, uflo / 5000, 20.0, 1.0E-30 / + + cdfx = 0.0 + + if (p.le.uflo.or.q.le.uflo.or.eps.le.uflo) then + iflag = 1 + endif + iflag = 0 + + if (x.le.0.0) return + if (x.ge.1.0) then + cdfx=1.0 + else + LL = (p+w).ge.(p+q+2.0*w)*x + if (LL) then + xy = x + yx = 1.-xy + pq = p + qp = q + else + yx = x + xy = 1.-yx + qp = p + pq = q + endif + + call gmln(pq,tmp) + dp = (pq-1.)*log(xy)-tmp + call gmln(qp,tmp) + dq = (qp-1.)*log(yx)-tmp + call gmln(pq+qp,tmp) + pdfl = tmp+dp+dq + + if (pdfl.ge.log(uflo)) then + u = exp(pdfl)*xy/pq + r = xy/yx + do while (qp.gt.1.) + if (u.le.eps*(1.-(pq+qp)*xy/(pq+1.))) then + if (.not.LL) cdfx = 1.-cdfx + return + endif + cdfx = cdfx+u + pq = pq+1. + qp = qp-1. + u = qp*r*u/pq + enddo + v = yx*u + yxeps = yx*eps + do j = 0, jmax + if (v.le.yxeps) then + if (.not.LL) cdfx = 1.-cdfx + return + endif + cdfx = cdfx + v + pq = pq+1. + v = (pq+qp-1.)*xy*v/pq + enddo + iflag = 2 + endif + if (.not.LL) cdfx = 1.-cdfx + endif + + end subroutine cdfbet + + subroutine gmln(x,y) + use machine + ! Computes the natural logarithm of the gamma distribution. Users + ! can set the absolute accuracy and corresponding xmin. + + implicit none + real(kind=kind_phys), intent(in) :: x + real(kind=kind_phys), intent(out) :: y +! local vars + integer i, n + real(kind=kind_phys) absacc, b1, b2, b3, b4, b5, b6, b7, b8 + real(kind=kind_phys) c, dx, q, r, xmin, xn +! data xmin, absacc / 6.894d0, 1.0E-15 / + data xmin, absacc / 1.357d0, 1.0E-3 / + data c / 0.918938533204672741780329736d0 / + data b1 / 0.833333333333333333333333333d-1 / + data b2 / - 0.277777777777777777777777778d-2 / + data b3 / 0.793650793650793650793650794d-3 / + data b4 / - 0.595238095238095238095238095d-3 / + data b5 / 0.841750841750841750841750842d-3 / + data b6 / - 0.191752691752691752691752692d-2 / + data b7 / 0.641025641025641025641025641d-2 / + data b8 / - 0.295506535947712418300653595d-1 / + + if (x.le.0.0) stop '*** x<=0.0 in function gamln ***' + dx = x + n = max(0,int(xmin - dx + 1.0d0) ) + xn = dx + n + r = 1.0d0/xn + q = r*r + y = r*( b1+q*( b2+q*( b3+q*( b4+q*( b5+q*( b6+q*( b7+q*b8 ) & + & )) ) ) ) ) +c + (xn-0.5d0)*log(xn)-xn + + if (n.gt.0) then + q = 1.0d0 + do i=0, n-1 + q = q*(dx+i) + enddo + y = y-log(q) + endif + + if (y + absacc.eq.y) then + print *,' ********* WARNING FROM FUNCTION GAMLN *********' + print *,' REQUIRED ABSOLUTE ACCURACY NOT ATTAINED FOR X = ',x + endif + return + end subroutine gmln + +! *** mg, sfc perts + + + + + + +!> @} ! !.........................................! end module module_radiation_surface ! diff --git a/physics/radlw_main.f b/physics/radlw_main.f index d0eeca63c..7ac2db8b4 100644 --- a/physics/radlw_main.f +++ b/physics/radlw_main.f @@ -230,7 +230,9 @@ ! apr 2012, b. ferrier and y. hou -- added conversion factor to fu's! ! cloud-snow optical property scheme. ! ! nov 2012, yu-tai hou -- modified control parameters thru ! -! module 'physparam'. ! +! module 'physparam'. ! +! FEB 2017 A.Cheng - add odpth output, effective radius input ! +! ! ! ! !!!!! ============================================================== !!!!! !!!!! end descriptions !!!!! @@ -420,6 +422,8 @@ end subroutine rrtmg_lw_init !! | cld_ref_rain | mean_effective_radius_for_rain_drop | mean effective radius for rain drop | micron | 2 | real | kind_phys | in | T | !! | cld_swp | cloud_snow_water_path | cloud snow water path | g m-2 | 2 | real | kind_phys | in | T | !! | cld_ref_snow | mean_effective_radius_for_snow_flake | mean effective radius for snow flake | micron | 2 | real | kind_phys | in | T | +!! | cld_od_weighted | cloud_optical_depth_weighted | cloud optical depth, weighted | none | 2 | real | kind_phys | in | T | +!! | cld_od_layer | cloud_optical_depth_678 | cloud optical depth, from bands 6,7,8 | none | 2 | real | kind_phys | in | T | !! | cld_od | cloud_optical_depth | cloud optical depth | none | 2 | real | kind_phys | in | T | !! | errmsg | error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | !! | errflg | error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | @@ -436,6 +440,7 @@ subroutine rrtmg_lw_run & & HLW0,HLWB,FLXPRF, & ! --- optional & cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice, & & cld_rwp,cld_ref_rain, cld_swp, cld_ref_snow, & + & cld_od_weighted, cld_od_layer, & & cld_od, errmsg, errflg & & ) @@ -630,6 +635,7 @@ subroutine rrtmg_lw_run & real (kind=kind_phys), dimension(npts,nlay),intent(in),optional:: & & cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice, & & cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow, & + & cld_od_weighted, cld_od_layer & & cld_od real (kind=kind_phys), dimension(npts), intent(in) :: sfemis, & @@ -720,12 +726,14 @@ subroutine rrtmg_lw_run & if ( .not.present(cld_lwp) .or. .not.present(cld_ref_liq) .or. & & .not.present(cld_iwp) .or. .not.present(cld_ref_ice) .or. & & .not.present(cld_rwp) .or. .not.present(cld_ref_rain) .or. & - & .not.present(cld_swp) .or. .not.present(cld_ref_snow)) then + & .not.present(cld_swp) .or. .not.present(cld_ref_snow) .or. & + & .not.present(cld_od_weighted) .or. .not.present(cld_od_layer)) then write(errmsg,'(*(a))') & & 'Logic error: ilwcliq>0 requires the following', & & ' optional arguments to be present:', & & ' cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice,', & - & ' cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow' + & ' cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow' & + & ' cld_od_weighted, cld_od_layer' errflg = 1 return end if @@ -1064,6 +1072,10 @@ subroutine rrtmg_lw_run & cldfmc = f_zero taucld = f_zero endif + do k = 1, nlay + cld_od_layer(iplon,k) = taucld(6,k) & + & + taucld(7,k) + taucld(8,k) + end do ! if (lprnt) then ! print *,' after cldprop' @@ -4720,80 +4732,43 @@ subroutine taugb05 jplp = jpl + 1 jmo3p = jmo3 + 1 - if (specparm < 0.125 .and. specparm1 < 0.125) then - p0 = fs - f_one - p40 = p0**4 + if (specparm < 0.125) then + p0 = fs - f_one + p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = fs1 - f_one - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 + 2 id210 = ind0 +11 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 + 2 - id211 = ind1 +11 - elseif (specparm > 0.875 .and. specparm1 > 0.875) then - p0 = -fs - p40 = p0**4 + elseif (specparm > 0.875) then + p0 = -fs + p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = -fs1 - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 + 1 id010 = ind0 +10 id100 = ind0 id110 = ind0 + 9 id200 = ind0 - 1 id210 = ind0 + 8 - - id001 = ind1 + 1 - id011 = ind1 +10 - id101 = ind1 - id111 = ind1 + 9 - id201 = ind1 - 1 - id211 = ind1 + 8 else fk00 = f_one - fs fk10 = fs fk20 = f_zero - fk01 = f_one - fs1 - fk11 = fs1 - fk21 = f_zero - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 id210 = ind0 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 - id211 = ind1 endif fac000 = fk00 * fac00(k) @@ -4803,6 +4778,45 @@ subroutine taugb05 fac110 = fk10 * fac10(k) fac210 = fk20 * fac10(k) + if (specparm1 < 0.125) then + p1 = fs1 - f_one + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + 2 + id211 = ind1 +11 + elseif (specparm1 > 0.875) then + p1 = -fs1 + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + 1 + id011 = ind1 +10 + id101 = ind1 + id111 = ind1 + 9 + id201 = ind1 - 1 + id211 = ind1 + 8 + else + fk01 = f_one - fs1 + fk11 = fs1 + fk21 = f_zero + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + id211 = ind1 + endif + fac001 = fk01 * fac01(k) fac101 = fk11 * fac01(k) fac201 = fk21 * fac01(k) @@ -5086,80 +5100,43 @@ subroutine taugb07 adjcolco2 = colamt(k,2) endif - if (specparm < 0.125 .and. specparm1 < 0.125) then + if (specparm < 0.125) then p0 = fs - f_one p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = fs1 - f_one - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 + 2 id210 = ind0 +11 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 + 2 - id211 = ind1 +11 - elseif (specparm > 0.875 .and. specparm1 > 0.875) then + elseif (specparm > 0.875) then p0 = -fs p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = -fs1 - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 + 1 id010 = ind0 +10 id100 = ind0 id110 = ind0 + 9 id200 = ind0 - 1 id210 = ind0 + 8 - - id001 = ind1 + 1 - id011 = ind1 +10 - id101 = ind1 - id111 = ind1 + 9 - id201 = ind1 - 1 - id211 = ind1 + 8 else fk00 = f_one - fs fk10 = fs fk20 = f_zero - fk01 = f_one - fs1 - fk11 = fs1 - fk21 = f_zero - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 id210 = ind0 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 - id211 = ind1 endif fac000 = fk00 * fac00(k) @@ -5169,6 +5146,45 @@ subroutine taugb07 fac110 = fk10 * fac10(k) fac210 = fk20 * fac10(k) + if (specparm1 < 0.125) then + p1 = fs1 - f_one + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + 2 + id211 = ind1 +11 + elseif (specparm1 > 0.875) then + p1 = -fs1 + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + 1 + id011 = ind1 +10 + id101 = ind1 + id111 = ind1 + 9 + id201 = ind1 - 1 + id211 = ind1 + 8 + else + fk01 = f_one - fs1 + fk11 = fs1 + fk21 = f_zero + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + id211 = ind1 + endif + fac001 = fk01 * fac01(k) fac101 = fk11 * fac01(k) fac201 = fk21 * fac01(k) @@ -5472,81 +5488,43 @@ subroutine taugb09 adjcoln2o = colamt(k,4) endif - if (specparm < 0.125 .and. specparm1 < 0.125) then + if (specparm < 0.125) then p0 = fs - f_one p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = fs1 - f_one - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 + 2 id210 = ind0 +11 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 + 2 - id211 = ind1 +11 - - elseif (specparm > 0.875 .and. specparm1 > 0.875) then + elseif (specparm > 0.875) then p0 = -fs p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = -fs1 - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 + 1 id010 = ind0 +10 id100 = ind0 id110 = ind0 + 9 id200 = ind0 - 1 id210 = ind0 + 8 - - id001 = ind1 + 1 - id011 = ind1 +10 - id101 = ind1 - id111 = ind1 + 9 - id201 = ind1 - 1 - id211 = ind1 + 8 else fk00 = f_one - fs fk10 = fs fk20 = f_zero - fk01 = f_one - fs1 - fk11 = fs1 - fk21 = f_zero - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 id210 = ind0 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 - id211 = ind1 endif fac000 = fk00 * fac00(k) @@ -5556,6 +5534,45 @@ subroutine taugb09 fac110 = fk10 * fac10(k) fac210 = fk20 * fac10(k) + if (specparm1 < 0.125) then + p1 = fs1 - f_one + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + 2 + id211 = ind1 +11 + elseif (specparm1 > 0.875) then + p1 = -fs1 + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + 1 + id011 = ind1 +10 + id101 = ind1 + id111 = ind1 + 9 + id201 = ind1 - 1 + id211 = ind1 + 8 + else + fk01 = f_one - fs1 + fk11 = fs1 + fk21 = f_zero + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + id211 = ind1 + endif + fac001 = fk01 * fac01(k) fac101 = fk11 * fac01(k) fac201 = fk21 * fac01(k) @@ -5862,80 +5879,43 @@ subroutine taugb12 indfp = indf + 1 jplp = jpl + 1 - if (specparm < 0.125 .and. specparm1 < 0.125) then + if (specparm < 0.125) then p0 = fs - f_one p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = fs1 - f_one - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 + 2 id210 = ind0 +11 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 + 2 - id211 = ind1 +11 - elseif (specparm > 0.875 .and. specparm1 > 0.875) then + elseif (specparm > 0.875) then p0 = -fs p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = -fs1 - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 + 1 id010 = ind0 +10 id100 = ind0 id110 = ind0 + 9 id200 = ind0 - 1 id210 = ind0 + 8 - - id001 = ind1 + 1 - id011 = ind1 +10 - id101 = ind1 - id111 = ind1 + 9 - id201 = ind1 - 1 - id211 = ind1 + 8 else fk00 = f_one - fs fk10 = fs fk20 = f_zero - fk01 = f_one - fs1 - fk11 = fs1 - fk21 = f_zero - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 id210 = ind0 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 - id211 = ind1 endif fac000 = fk00 * fac00(k) @@ -5945,6 +5925,45 @@ subroutine taugb12 fac110 = fk10 * fac10(k) fac210 = fk20 * fac10(k) + if (specparm1 < 0.125) then + p1 = fs1 - f_one + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + 2 + id211 = ind1 +11 + elseif (specparm1 > 0.875) then + p1 = -fs1 + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + 1 + id011 = ind1 +10 + id101 = ind1 + id111 = ind1 + 9 + id201 = ind1 - 1 + id211 = ind1 + 8 + else + fk01 = f_one - fs1 + fk11 = fs1 + fk21 = f_zero + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + id211 = ind1 + endif + fac001 = fk01 * fac01(k) fac101 = fk11 * fac01(k) fac201 = fk21 * fac01(k) @@ -6092,80 +6111,43 @@ subroutine taugb13 adjcolco2 = colamt(k,2) endif - if (specparm < 0.125 .and. specparm1 < 0.125) then + if (specparm < 0.125) then p0 = fs - f_one p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = fs1 - f_one - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 + 2 id210 = ind0 +11 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 + 2 - id211 = ind1 +11 - elseif (specparm > 0.875 .and. specparm1 > 0.875) then + elseif (specparm > 0.875) then p0 = -fs p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = -fs1 - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 + 1 id010 = ind0 +10 id100 = ind0 id110 = ind0 + 9 id200 = ind0 - 1 id210 = ind0 + 8 - - id001 = ind1 + 1 - id011 = ind1 +10 - id101 = ind1 - id111 = ind1 + 9 - id201 = ind1 - 1 - id211 = ind1 + 8 else fk00 = f_one - fs fk10 = fs fk20 = f_zero - fk01 = f_one - fs1 - fk11 = fs1 - fk21 = f_zero - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 id210 = ind0 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 - id211 = ind1 endif fac000 = fk00 * fac00(k) @@ -6175,6 +6157,45 @@ subroutine taugb13 fac110 = fk10 * fac10(k) fac210 = fk20 * fac10(k) + if (specparm1 < 0.125) then + p1 = fs1 - f_one + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + 2 + id211 = ind1 +11 + elseif (specparm1 > 0.875) then + p1 = -fs1 + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + 1 + id011 = ind1 +10 + id101 = ind1 + id111 = ind1 + 9 + id201 = ind1 - 1 + id211 = ind1 + 8 + else + fk01 = f_one - fs1 + fk11 = fs1 + fk21 = f_zero + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + id211 = ind1 + endif + fac001 = fk01 * fac01(k) fac101 = fk11 * fac01(k) fac201 = fk21 * fac01(k) @@ -6385,81 +6406,43 @@ subroutine taugb15 jplp = jpl + 1 jmn2p = jmn2 + 1 - - if (specparm < 0.125 .and. specparm1 < 0.125) then + if (specparm < 0.125) then p0 = fs - f_one p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = fs1 - f_one - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 + 2 id210 = ind0 +11 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 + 2 - id211 = ind1 +11 - elseif (specparm > 0.875 .and. specparm1 > 0.875) then + elseif (specparm > 0.875) then p0 = -fs p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = -fs1 - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 + 1 id010 = ind0 +10 id100 = ind0 id110 = ind0 + 9 id200 = ind0 - 1 id210 = ind0 + 8 - - id001 = ind1 + 1 - id011 = ind1 +10 - id101 = ind1 - id111 = ind1 + 9 - id201 = ind1 - 1 - id211 = ind1 + 8 else fk00 = f_one - fs fk10 = fs fk20 = f_zero - fk01 = f_one - fs1 - fk11 = fs1 - fk21 = f_zero - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 id210 = ind0 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 - id211 = ind1 endif fac000 = fk00 * fac00(k) @@ -6469,6 +6452,45 @@ subroutine taugb15 fac110 = fk10 * fac10(k) fac210 = fk20 * fac10(k) + if (specparm1 < 0.125) then + p1 = fs1 - f_one + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + 2 + id211 = ind1 +11 + elseif (specparm1 > 0.875) then + p1 = -fs1 + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + 1 + id011 = ind1 +10 + id101 = ind1 + id111 = ind1 + 9 + id201 = ind1 - 1 + id211 = ind1 + 8 + else + fk01 = f_one - fs1 + fk11 = fs1 + fk21 = f_zero + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + id211 = ind1 + endif + fac001 = fk01 * fac01(k) fac101 = fk11 * fac01(k) fac201 = fk21 * fac01(k) @@ -6577,80 +6599,43 @@ subroutine taugb16 indfp = indf + 1 jplp = jpl + 1 - if (specparm < 0.125 .and. specparm1 < 0.125) then + if (specparm < 0.125) then p0 = fs - f_one p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = fs1 - f_one - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 + 2 id210 = ind0 +11 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 + 2 - id211 = ind1 +11 - elseif (specparm > 0.875 .and. specparm1 > 0.875) then + elseif (specparm > 0.875) then p0 = -fs p40 = p0**4 fk00 = p40 fk10 = f_one - p0 - 2.0*p40 fk20 = p0 + p40 - p1 = -fs1 - p41 = p1**4 - fk01 = p41 - fk11 = f_one - p1 - 2.0*p41 - fk21 = p1 + p41 - id000 = ind0 + 1 id010 = ind0 +10 id100 = ind0 id110 = ind0 + 9 id200 = ind0 - 1 id210 = ind0 + 8 - - id001 = ind1 + 1 - id011 = ind1 +10 - id101 = ind1 - id111 = ind1 + 9 - id201 = ind1 - 1 - id211 = ind1 + 8 else fk00 = f_one - fs fk10 = fs fk20 = f_zero - fk01 = f_one - fs1 - fk11 = fs1 - fk21 = f_zero - id000 = ind0 id010 = ind0 + 9 id100 = ind0 + 1 id110 = ind0 +10 id200 = ind0 id210 = ind0 - - id001 = ind1 - id011 = ind1 + 9 - id101 = ind1 + 1 - id111 = ind1 +10 - id201 = ind1 - id211 = ind1 endif fac000 = fk00 * fac00(k) @@ -6660,6 +6645,45 @@ subroutine taugb16 fac110 = fk10 * fac10(k) fac210 = fk20 * fac10(k) + if (specparm1 < 0.125) then + p1 = fs1 - f_one + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + 2 + id211 = ind1 +11 + elseif (specparm1 > 0.875) then + p1 = -fs1 + p41 = p1**4 + fk01 = p41 + fk11 = f_one - p1 - 2.0*p41 + fk21 = p1 + p41 + + id001 = ind1 + 1 + id011 = ind1 +10 + id101 = ind1 + id111 = ind1 + 9 + id201 = ind1 - 1 + id211 = ind1 + 8 + else + fk01 = f_one - fs1 + fk11 = fs1 + fk21 = f_zero + + id001 = ind1 + id011 = ind1 + 9 + id101 = ind1 + 1 + id111 = ind1 +10 + id201 = ind1 + id211 = ind1 + endif + fac001 = fk01 * fac01(k) fac101 = fk11 * fac01(k) fac201 = fk21 * fac01(k) diff --git a/physics/radsw_main.f b/physics/radsw_main.f index 4d66815c4..91cd77fb3 100644 --- a/physics/radsw_main.f +++ b/physics/radsw_main.f @@ -554,6 +554,8 @@ end subroutine rrtmg_sw_init !! | cld_ref_rain | mean_effective_radius_for_rain_drop | mean effective radius for rain drop | micron | 2 | real | kind_phys | in | T | !! | cld_swp | cloud_snow_water_path | cloud snow water path | g m-2 | 2 | real | kind_phys | in | T | !! | cld_ref_snow | mean_effective_radius_for_snow_flake | mean effective radius for snow flake | micron | 2 | real | kind_phys | in | T | +!! | cld_od_weighted | cloud_optical_depth_weighted | cloud optical depth, weighted | none | 2 | real | kind_phys | in | T | +!! | cld_od_layer | cloud_optical_depth_678 | cloud optical depth, from bands 6,7,8 | none | 2 | real | kind_phys | in | T | !! | cld_od | cloud_optical_depth | cloud optical depth | none | 2 | real | kind_phys | in | T | !! | cld_ssa | cloud_single_scattering_albedo | cloud single scattering albedo | frac | 2 | real | kind_phys | in | T | !! | cld_asy | cloud_asymmetry_parameter | cloud asymmetry parameter | none | 2 | real | kind_phys | in | T | @@ -577,6 +579,7 @@ subroutine rrtmg_sw_run & & HSW0,HSWB,FLXPRF,FDNCMP, & ! --- optional & cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice, & & cld_rwp,cld_ref_rain, cld_swp, cld_ref_snow, & + & cld_od_weighted, cld_od_layer, & & cld_od, cld_ssa, cld_asy, errmsg, errflg & ) @@ -786,6 +789,7 @@ subroutine rrtmg_sw_run & real (kind=kind_phys), dimension(npts,nlay),intent(in),optional:: & & cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice, & & cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow, & + & cld_od_weighted, cld_od_layer, & & cld_od, cld_ssa, cld_asy real(kind=kind_phys),dimension(npts,nlay,nbdsw),intent(in)::aeraod @@ -899,12 +903,14 @@ subroutine rrtmg_sw_run & if ( .not.present(cld_lwp) .or. .not.present(cld_ref_liq) .or. & & .not.present(cld_iwp) .or. .not.present(cld_ref_ice) .or. & & .not.present(cld_rwp) .or. .not.present(cld_ref_rain) .or. & - & .not.present(cld_swp) .or. .not.present(cld_ref_snow)) then + & .not.present(cld_swp) .or. .not.present(cld_ref_snow) .or. & + & .not.present(cld_od_weighted) .or. .not.present(cld_od_layer)) then write(errmsg,'(*(a))') & & 'Logic error: ilwcliq>0 requires the following', & & ' optional arguments to be present:', & & ' cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice,', & - & ' cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow' + & ' cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow', & + & ' cld_od_weighted, cld_od_layer' errflg = 1 return end if @@ -1207,7 +1213,9 @@ subroutine rrtmg_sw_run & enddo enddo endif ! end if_zcf1_block - + do k = 1, nlay + clouds(j1,k,10) = taucw(k,10) + end do !> -# Call setcoef() to compute various coefficients needed in !! radiative transfer calculations. call setcoef & From 2d316c35f8d272a8ece514489d82524c3fa0168a Mon Sep 17 00:00:00 2001 From: Laurie Carson Date: Wed, 13 Jun 2018 15:21:29 -0600 Subject: [PATCH 2/3] update long-names for clouds(10&11) --- GFS_layer/GFS_radiation_driver.F90 | 4 ++-- physics/radlw_main.f | 12 ++++++------ physics/radsw_main.f | 12 ++++++------ 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/GFS_layer/GFS_radiation_driver.F90 b/GFS_layer/GFS_radiation_driver.F90 index b5473b765..25c7634fb 100644 --- a/GFS_layer/GFS_radiation_driver.F90 +++ b/GFS_layer/GFS_radiation_driver.F90 @@ -1256,7 +1256,7 @@ subroutine GFS_radiation_driver & cld_iwp=clouds(:,:,4), cld_ref_ice=clouds(:,:,5), & cld_rwp=clouds(:,:,6), cld_ref_rain=clouds(:,:,7), & cld_swp=clouds(:,:,8), cld_ref_snow=clouds(:,:,9), & - cld_od_total=clouds(:,:,10), cld_od_layer w=clouds(:,:,11), & + cld_od_total=clouds(:,:,10), cld_od_layer=clouds(:,:,11), & errmsg=errmsg, errflg=errflg) ! DH* @@ -1291,7 +1291,7 @@ subroutine GFS_radiation_driver & cld_iwp=clouds(:,:,4), cld_ref_ice=clouds(:,:,5), & cld_rwp=clouds(:,:,6), cld_ref_rain=clouds(:,:,7), & cld_swp=clouds(:,:,8), cld_ref_snow=clouds(:,:,9), & - cld_od_total=clouds(:,:,10), cld_od_layer w=clouds(:,:,11), & + cld_od_total=clouds(:,:,10), cld_od_layer=clouds(:,:,11), & errmsg=errmsg, errflg=errflg) !CCPP: L1718-1747 diff --git a/physics/radlw_main.f b/physics/radlw_main.f index 7ac2db8b4..61b03f40e 100644 --- a/physics/radlw_main.f +++ b/physics/radlw_main.f @@ -422,8 +422,8 @@ end subroutine rrtmg_lw_init !! | cld_ref_rain | mean_effective_radius_for_rain_drop | mean effective radius for rain drop | micron | 2 | real | kind_phys | in | T | !! | cld_swp | cloud_snow_water_path | cloud snow water path | g m-2 | 2 | real | kind_phys | in | T | !! | cld_ref_snow | mean_effective_radius_for_snow_flake | mean effective radius for snow flake | micron | 2 | real | kind_phys | in | T | -!! | cld_od_weighted | cloud_optical_depth_weighted | cloud optical depth, weighted | none | 2 | real | kind_phys | in | T | -!! | cld_od_layer | cloud_optical_depth_678 | cloud optical depth, from bands 6,7,8 | none | 2 | real | kind_phys | in | T | +!! | cld_od_total | cloud_optical_depth_total | cloud optical depth, total | none | 2 | real | kind_phys | in | T | +!! | cld_od_layer | cloud_optical_depth_layers_678 | cloud optical depth, from bands 6,7,8 | none | 2 | real | kind_phys | in | T | !! | cld_od | cloud_optical_depth | cloud optical depth | none | 2 | real | kind_phys | in | T | !! | errmsg | error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | !! | errflg | error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | @@ -440,7 +440,7 @@ subroutine rrtmg_lw_run & & HLW0,HLWB,FLXPRF, & ! --- optional & cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice, & & cld_rwp,cld_ref_rain, cld_swp, cld_ref_snow, & - & cld_od_weighted, cld_od_layer, & + & cld_od_total, cld_od_layer, & & cld_od, errmsg, errflg & & ) @@ -635,7 +635,7 @@ subroutine rrtmg_lw_run & real (kind=kind_phys), dimension(npts,nlay),intent(in),optional:: & & cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice, & & cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow, & - & cld_od_weighted, cld_od_layer & + & cld_od_total, cld_od_layer & & cld_od real (kind=kind_phys), dimension(npts), intent(in) :: sfemis, & @@ -727,13 +727,13 @@ subroutine rrtmg_lw_run & & .not.present(cld_iwp) .or. .not.present(cld_ref_ice) .or. & & .not.present(cld_rwp) .or. .not.present(cld_ref_rain) .or. & & .not.present(cld_swp) .or. .not.present(cld_ref_snow) .or. & - & .not.present(cld_od_weighted) .or. .not.present(cld_od_layer)) then + & .not.present(cld_od_total) .or. .not.present(cld_od_layer)) then write(errmsg,'(*(a))') & & 'Logic error: ilwcliq>0 requires the following', & & ' optional arguments to be present:', & & ' cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice,', & & ' cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow' & - & ' cld_od_weighted, cld_od_layer' + & ' cld_od_total, cld_od_layer' errflg = 1 return end if diff --git a/physics/radsw_main.f b/physics/radsw_main.f index 91cd77fb3..eaa199f4f 100644 --- a/physics/radsw_main.f +++ b/physics/radsw_main.f @@ -554,8 +554,8 @@ end subroutine rrtmg_sw_init !! | cld_ref_rain | mean_effective_radius_for_rain_drop | mean effective radius for rain drop | micron | 2 | real | kind_phys | in | T | !! | cld_swp | cloud_snow_water_path | cloud snow water path | g m-2 | 2 | real | kind_phys | in | T | !! | cld_ref_snow | mean_effective_radius_for_snow_flake | mean effective radius for snow flake | micron | 2 | real | kind_phys | in | T | -!! | cld_od_weighted | cloud_optical_depth_weighted | cloud optical depth, weighted | none | 2 | real | kind_phys | in | T | -!! | cld_od_layer | cloud_optical_depth_678 | cloud optical depth, from bands 6,7,8 | none | 2 | real | kind_phys | in | T | +!! | cld_od_total | cloud_optical_depth_total | cloud optical depth, weighted | none | 2 | real | kind_phys | in | T | +!! | cld_od_layer | cloud_optical_depth_layers_678 | cloud optical depth, from bands 6,7,8 | none | 2 | real | kind_phys | in | T | !! | cld_od | cloud_optical_depth | cloud optical depth | none | 2 | real | kind_phys | in | T | !! | cld_ssa | cloud_single_scattering_albedo | cloud single scattering albedo | frac | 2 | real | kind_phys | in | T | !! | cld_asy | cloud_asymmetry_parameter | cloud asymmetry parameter | none | 2 | real | kind_phys | in | T | @@ -579,7 +579,7 @@ subroutine rrtmg_sw_run & & HSW0,HSWB,FLXPRF,FDNCMP, & ! --- optional & cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice, & & cld_rwp,cld_ref_rain, cld_swp, cld_ref_snow, & - & cld_od_weighted, cld_od_layer, & + & cld_od_total, cld_od_layer, & & cld_od, cld_ssa, cld_asy, errmsg, errflg & ) @@ -789,7 +789,7 @@ subroutine rrtmg_sw_run & real (kind=kind_phys), dimension(npts,nlay),intent(in),optional:: & & cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice, & & cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow, & - & cld_od_weighted, cld_od_layer, & + & cld_od_total, cld_od_layer, & & cld_od, cld_ssa, cld_asy real(kind=kind_phys),dimension(npts,nlay,nbdsw),intent(in)::aeraod @@ -904,13 +904,13 @@ subroutine rrtmg_sw_run & & .not.present(cld_iwp) .or. .not.present(cld_ref_ice) .or. & & .not.present(cld_rwp) .or. .not.present(cld_ref_rain) .or. & & .not.present(cld_swp) .or. .not.present(cld_ref_snow) .or. & - & .not.present(cld_od_weighted) .or. .not.present(cld_od_layer)) then + & .not.present(cld_od_total) .or. .not.present(cld_od_layer)) then write(errmsg,'(*(a))') & & 'Logic error: ilwcliq>0 requires the following', & & ' optional arguments to be present:', & & ' cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice,', & & ' cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow', & - & ' cld_od_weighted, cld_od_layer' + & ' cld_od_total, cld_od_layer' errflg = 1 return end if From 7e02a215680ebda77e4e122cdf3e9bcc9051e464 Mon Sep 17 00:00:00 2001 From: Laurie Carson Date: Thu, 14 Jun 2018 11:53:45 -0600 Subject: [PATCH 3/3] syntax fixes --- physics/radlw_main.f | 11 +++++++---- physics/radsw_main.f | 11 +++++++---- 2 files changed, 14 insertions(+), 8 deletions(-) diff --git a/physics/radlw_main.f b/physics/radlw_main.f index 61b03f40e..98ea25df8 100644 --- a/physics/radlw_main.f +++ b/physics/radlw_main.f @@ -423,7 +423,7 @@ end subroutine rrtmg_lw_init !! | cld_swp | cloud_snow_water_path | cloud snow water path | g m-2 | 2 | real | kind_phys | in | T | !! | cld_ref_snow | mean_effective_radius_for_snow_flake | mean effective radius for snow flake | micron | 2 | real | kind_phys | in | T | !! | cld_od_total | cloud_optical_depth_total | cloud optical depth, total | none | 2 | real | kind_phys | in | T | -!! | cld_od_layer | cloud_optical_depth_layers_678 | cloud optical depth, from bands 6,7,8 | none | 2 | real | kind_phys | in | T | +!! | cld_od_layer | cloud_optical_depth_layers_678 | cloud optical depth, from bands 6,7,8 | none | 2 | real | kind_phys | out | T | !! | cld_od | cloud_optical_depth | cloud optical depth | none | 2 | real | kind_phys | in | T | !! | errmsg | error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | !! | errflg | error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | @@ -635,7 +635,9 @@ subroutine rrtmg_lw_run & real (kind=kind_phys), dimension(npts,nlay),intent(in),optional:: & & cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice, & & cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow, & - & cld_od_total, cld_od_layer & + & cld_od_total + real (kind=kind_phys), dimension(npts,nlay),intent(out),optional:: & + & cld_od_layer, & & cld_od real (kind=kind_phys), dimension(npts), intent(in) :: sfemis, & @@ -727,12 +729,13 @@ subroutine rrtmg_lw_run & & .not.present(cld_iwp) .or. .not.present(cld_ref_ice) .or. & & .not.present(cld_rwp) .or. .not.present(cld_ref_rain) .or. & & .not.present(cld_swp) .or. .not.present(cld_ref_snow) .or. & - & .not.present(cld_od_total) .or. .not.present(cld_od_layer)) then + & .not.present(cld_od_total) .or. & + & .not.present(cld_od_layer)) then write(errmsg,'(*(a))') & & 'Logic error: ilwcliq>0 requires the following', & & ' optional arguments to be present:', & & ' cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice,', & - & ' cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow' & + & ' cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow', & & ' cld_od_total, cld_od_layer' errflg = 1 return diff --git a/physics/radsw_main.f b/physics/radsw_main.f index eaa199f4f..6e1988a7e 100644 --- a/physics/radsw_main.f +++ b/physics/radsw_main.f @@ -554,7 +554,7 @@ end subroutine rrtmg_sw_init !! | cld_ref_rain | mean_effective_radius_for_rain_drop | mean effective radius for rain drop | micron | 2 | real | kind_phys | in | T | !! | cld_swp | cloud_snow_water_path | cloud snow water path | g m-2 | 2 | real | kind_phys | in | T | !! | cld_ref_snow | mean_effective_radius_for_snow_flake | mean effective radius for snow flake | micron | 2 | real | kind_phys | in | T | -!! | cld_od_total | cloud_optical_depth_total | cloud optical depth, weighted | none | 2 | real | kind_phys | in | T | +!! | cld_od_total | cloud_optical_depth_total | cloud optical depth, weighted | none | 2 | real | kind_phys | out | T | !! | cld_od_layer | cloud_optical_depth_layers_678 | cloud optical depth, from bands 6,7,8 | none | 2 | real | kind_phys | in | T | !! | cld_od | cloud_optical_depth | cloud optical depth | none | 2 | real | kind_phys | in | T | !! | cld_ssa | cloud_single_scattering_albedo | cloud single scattering albedo | frac | 2 | real | kind_phys | in | T | @@ -789,8 +789,10 @@ subroutine rrtmg_sw_run & real (kind=kind_phys), dimension(npts,nlay),intent(in),optional:: & & cld_lwp, cld_ref_liq, cld_iwp, cld_ref_ice, & & cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow, & - & cld_od_total, cld_od_layer, & + & cld_od_layer, & & cld_od, cld_ssa, cld_asy + real (kind=kind_phys), dimension(npts,nlay),intent(out),optional::& + & cld_od_total real(kind=kind_phys),dimension(npts,nlay,nbdsw),intent(in)::aeraod real(kind=kind_phys),dimension(npts,nlay,nbdsw),intent(in)::aerssa @@ -904,7 +906,8 @@ subroutine rrtmg_sw_run & & .not.present(cld_iwp) .or. .not.present(cld_ref_ice) .or. & & .not.present(cld_rwp) .or. .not.present(cld_ref_rain) .or. & & .not.present(cld_swp) .or. .not.present(cld_ref_snow) .or. & - & .not.present(cld_od_total) .or. .not.present(cld_od_layer)) then + & .not.present(cld_od_total) .or. & + & .not.present(cld_od_layer)) then write(errmsg,'(*(a))') & & 'Logic error: ilwcliq>0 requires the following', & & ' optional arguments to be present:', & @@ -1214,7 +1217,7 @@ subroutine rrtmg_sw_run & enddo endif ! end if_zcf1_block do k = 1, nlay - clouds(j1,k,10) = taucw(k,10) + cld_od_total(j1,k) = taucw(k,10) end do !> -# Call setcoef() to compute various coefficients needed in !! radiative transfer calculations.