From 2ec28c273bd3d02b86a83339f95c37c4c327042e Mon Sep 17 00:00:00 2001 From: brianhenn Date: Tue, 13 Jun 2023 16:27:30 +0000 Subject: [PATCH 1/2] refactor rad H&M scheme to subroutine --- FV3/gfsphysics/physics/radiation_clouds.f | 241 +++++++++------------- 1 file changed, 93 insertions(+), 148 deletions(-) diff --git a/FV3/gfsphysics/physics/radiation_clouds.f b/FV3/gfsphysics/physics/radiation_clouds.f index 99d58b677..c00628161 100644 --- a/FV3/gfsphysics/physics/radiation_clouds.f +++ b/FV3/gfsphysics/physics/radiation_clouds.f @@ -792,29 +792,11 @@ subroutine progcld1 & !! and McFarquhar (1996) \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 + call effective_ice_cloud_radius_heymsfield_mcfarquhar & + & ( tlyr, cip, plyr, delp, tvly, & ! --- inputs: + & IX, NLAY, & + & rei & ! --- outputs: + & ) endif ! @@ -1233,35 +1215,11 @@ subroutine progcld2 & !> -# Calculate effective ice cloud droplet radius. - do k = 1, NLAY - do i = 1, IX - tem1 = tlyr(i,k) - con_ttp - tem2 = cip(i,k) - - if (tem2 > 0.0) then - tem3 = tem2d(i,k) * tem2 / tvly(i,k) - - if (tem1 < -50.0) then - rei(i,k) = (1250.0/9.917) * tem3 ** 0.109 - elseif (tem1 < -40.0) then - rei(i,k) = (1250.0/9.337) * tem3 ** 0.08 - elseif (tem1 < -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 - -! if (lprnt .and. k == l) print *,' reiL=',rei(i,k),' icec=', & -! & icec,' cip=',cip(i,k),' tem=',tem,' delt=',delt - - rei(i,k) = max(10.0, min(rei(i,k), 300.0)) -! rei(i,k) = max(20.0, min(rei(i,k), 300.0)) -!!!! rei(i,k) = max(30.0, min(rei(i,k), 300.0)) -! rei(i,k) = max(50.0, min(rei(i,k), 300.0)) -! rei(i,k) = max(100.0, min(rei(i,k), 300.0)) - endif - enddo - enddo + call effective_ice_cloud_radius_heymsfield_mcfarquhar & + & ( tlyr, cip, plyr, delp, tvly, & ! --- inputs: + & IX, NLAY, & + & rei & ! --- outputs: + & ) ! do k = 1, NLAY do i = 1, IX @@ -1633,31 +1591,11 @@ subroutine progcld3 & !> -# Calculate effective ice cloud droplet radius following Heymsfield !! and McFarquhar (1996) \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) - 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 - + call effective_ice_cloud_radius_heymsfield_mcfarquhar & + & ( tlyr, cip, plyr, delp, tvly, & ! --- inputs: + & IX, NLAY, & + & rei & ! --- outputs: + & ) ! do k = 1, nlay do i = 1, ix @@ -1933,30 +1871,11 @@ subroutine progcld4 & ! --- 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 - + call effective_ice_cloud_radius_heymsfield_mcfarquhar & + & ( tlyr, cip, plyr, delp, tvly, & ! --- inputs: + & IX, NLAY, & + & rei & ! --- outputs: + & ) ! do k = 1, NLAY do i = 1, IX @@ -2217,30 +2136,11 @@ subroutine progcld4o & ! --- 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 - + call effective_ice_cloud_radius_heymsfield_mcfarquhar & + & ( tlyr, cip, plyr, delp, tvly, & ! --- inputs: + & IX, NLAY, & + & rei & ! --- outputs: + & ) ! do k = 1, NLAY do i = 1, IX @@ -2910,29 +2810,11 @@ subroutine progclduni & !> -# Compute effective ice cloud droplet radius following Heymsfield !! and McFarquhar (1996) \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)) - - 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 + call effective_ice_cloud_radius_heymsfield_mcfarquhar & + & ( tlyr, cip, plyr, delp, tvly, & ! --- inputs: + & IX, NLAY, & + & rei & ! --- outputs: + & ) endif ! do k = 1, NLAY @@ -3401,6 +3283,69 @@ end subroutine gethml !----------------------------------- !! @} + + + subroutine effective_ice_cloud_radius_heymsfield_mcfarquhar & + & ( tlyr, cip, plyr, delp, tvly, & ! --- inputs: + & IX, NLAY, & + & rei & ! --- outputs: + & ) +! =================================================================== ! +! ! +! abstract: compute cloud ice droplet effective radius, according to ! +! the Heymsfield and McFaruhar (1996) scheme +! ! +! ==================== definition of variables ==================== ! +! ! +! input variables: ! +! tlyr (IX,NLAY) : model layer mean temperature in k ! +! cip (IX,NLAY) : ice condensate path in g/m**2 ! +! plyr (IX,NLAY) : model layer mean pressure in mb (100Pa) ! +! delp (ix,nlay) : model layer pressure thickness in mb (100Pa) ! +! tvly (IX,NLAY) : model layer virtual temperature in k ! +! IX : horizontal dimention ! +! NLAY : vertical layer dimensions ! +! ! +! output variables: ! +! rei (IX,NLAY) : effective cloud ice droplet radius ! + + implicit none + + integer, intent(in) :: IX, NLAY + + real (kind=kind_phys), dimension(:,:), intent(in) :: tlyr, cip, & + & plyr, delp, tvly + + real (kind=kind_phys), dimension(:,:), intent(out) :: rei + + integer :: i, k + + real (kind=kind_phys) :: tem, tem2 + + do k = 1, NLAY + do i = 1, IX + tem = tlyr(i,k) - con_ttp + + if (cip(i,k) > 0.0) then + tem2 = gord * cip(i,k) * plyr(i,k) / (delp(i,k)*tvly(i,k)) + + if (tem < -50.0) then + rei(i,k) = (1250.0/9.917) * tem2 ** 0.109 + elseif (tem < -40.0) then + rei(i,k) = (1250.0/9.337) * tem2 ** 0.08 + elseif (tem < -30.0) then + rei(i,k) = (1250.0/9.208) * tem2 ** 0.055 + else + rei(i,k) = (1250.0/9.387) * tem2 ** 0.031 + endif + rei(i,k) = max(10.0, min(rei(i,k), 150.0)) + endif + enddo + enddo + +!................................... + end subroutine effective_ice_cloud_radius_heymsfield_mcfarquhar +!----------------------------------- ! !........................................! end module module_radiation_clouds ! From c2de2642c61b2a6a82c040347d2ad5d559f2c801 Mon Sep 17 00:00:00 2001 From: brianhenn Date: Thu, 15 Jun 2023 16:56:32 +0000 Subject: [PATCH 2/2] incomplete implementation of reiflag_rad and kristjansson scheme --- FV3/gfsphysics/GFS_layer/GFS_driver.F90 | 3 +- FV3/gfsphysics/GFS_layer/GFS_typedefs.F90 | 8 ++- FV3/gfsphysics/physics/physparam.f | 2 + FV3/gfsphysics/physics/rad_initialize.f | 10 ++- FV3/gfsphysics/physics/radiation_clouds.f | 79 +++++++++++++++++++++-- 5 files changed, 92 insertions(+), 10 deletions(-) diff --git a/FV3/gfsphysics/GFS_layer/GFS_driver.F90 b/FV3/gfsphysics/GFS_layer/GFS_driver.F90 index 71bed62b9..c7d38c861 100644 --- a/FV3/gfsphysics/GFS_layer/GFS_driver.F90 +++ b/FV3/gfsphysics/GFS_layer/GFS_driver.F90 @@ -319,7 +319,8 @@ subroutine GFS_initialize (Model, Statein, Stateout, Sfcprop, & Model%ntcw, Model%num_p2d, Model%num_p3d, Model%npdf3d, & Model%ntoz, Model%iovr_sw, Model%iovr_lw, Model%isubc_sw, & Model%isubc_lw, Model%icliq_sw, Model%crick_proof, Model%ccnorm,& - Model%imp_physics, Model%norad_precip, Model%idate, Model%iflip, Model%me) + Model%reiflag_rad, Model%imp_physics, Model%norad_precip, Model%idate, & + Model%iflip, Model%me) deallocate (si) #endif diff --git a/FV3/gfsphysics/GFS_layer/GFS_typedefs.F90 b/FV3/gfsphysics/GFS_layer/GFS_typedefs.F90 index 94892ab4d..59f83de4c 100644 --- a/FV3/gfsphysics/GFS_layer/GFS_typedefs.F90 +++ b/FV3/gfsphysics/GFS_layer/GFS_typedefs.F90 @@ -676,6 +676,7 @@ module GFS_typedefs logical :: lwhtr !< flag to output lw heating rate (Radtend%lwhc) logical :: swhtr !< flag to output sw heating rate (Radtend%swhc) logical :: do_only_clearsky_rad !< flag for whether to do only clear-sky radiation + integer :: reiflag_rad !< effective ice radius control flag !--- microphysical switch integer :: ncld !< choice of cloud scheme @@ -2855,6 +2856,10 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & logical :: lwhtr = .true. !< flag to output lw heating rate (Radtend%lwhc) logical :: swhtr = .true. !< flag to output sw heating rate (Radtend%swhc) logical :: do_only_clearsky_rad = .false. !< flag for whether to do only clear-sky radiation + integer :: reiflag_rad = 1 !< ice effective radius computation method used + !< within radiation scheme: + !< reiflag_rad=1 => Heymsfield and Mcfarquhar 1996 + !< reiflag_rad=4 => Kristjansson et al. 2000 !--- Z-C microphysical parameters integer :: ncld = 1 !< choice of cloud scheme @@ -3191,7 +3196,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & fhswr, fhlwr, levr, nfxr, aero_in, iflip, isol, ico2, ialb, & isot, iems, iaer, icliq_sw, iovr_sw, iovr_lw, ictm, isubc_sw,& isubc_lw, crick_proof, ccnorm, lwhtr, swhtr, & - do_only_clearsky_rad, & + do_only_clearsky_rad, reiflag_rad, & ! IN CCN forcing iccn, & !--- microphysical parameterizations @@ -3440,6 +3445,7 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, & Model%lwhtr = lwhtr Model%swhtr = swhtr Model%do_only_clearsky_rad = do_only_clearsky_rad + Model%reiflag_rad = reiflag_rad #ifdef CCPP ! The CCPP versions of the RRTMG lw/sw schemes are configured ! such that lw and sw heating rate are output, i.e. they rely diff --git a/FV3/gfsphysics/physics/physparam.f b/FV3/gfsphysics/physics/physparam.f index f78191278..3a6976ab9 100644 --- a/FV3/gfsphysics/physics/physparam.f +++ b/FV3/gfsphysics/physics/physparam.f @@ -270,6 +270,8 @@ module physparam ! logical, save :: lnoprec =.false. !> shallow convetion flag logical, save :: lsashal =.false. +!> ice effective radius flag + integer, save :: reiflagrad = 1 ! ............................................. ! !>\name -2.5- For module radiation_surface diff --git a/FV3/gfsphysics/physics/rad_initialize.f b/FV3/gfsphysics/physics/rad_initialize.f index 0a3d307c1..3af85fc86 100644 --- a/FV3/gfsphysics/physics/rad_initialize.f +++ b/FV3/gfsphysics/physics/rad_initialize.f @@ -4,7 +4,7 @@ subroutine rad_initialize & ! --- inputs: & ( si,levr,ictm,isol,ico2,iaer,ialb,iems,ntcw, num_p2d, & & num_p3d,npdf3d,ntoz,iovr_sw,iovr_lw,isubc_sw,isubc_lw, & - & icliq_sw,crick_proof,ccnorm, & + & icliq_sw,crick_proof,ccnorm,reiflag_rad, & & imp_physics,norad_precip,idate,iflip,me ) ! --- outputs: ( none ) @@ -99,6 +99,7 @@ subroutine rad_initialize & ! =2: mcica sub-col approx. provided random seed ! ! crick_proof : control flag for eliminating CRICK ! ! ccnorm : control flag for in-cloud condensate mixing ratio! +! reiflag_rad : control flag for ice effective radius method ! ! norad_precip : control flag for not using precip in radiation ! ! idate(4) : ncep absolute date and time of initial condition ! ! (hour, month, day, year) ! @@ -116,7 +117,7 @@ subroutine rad_initialize & & iaermdl, icldflg, & & iovrsw , iovrlw , lcrick , lcnorm , lnoprec, & & ialbflg, iemsflg, isubcsw, isubclw, ivflip , ipsd0, & - & iswcliq, & + & iswcliq, reiflagrad, & & kind_phys use module_radiation_driver, only : radinit @@ -126,7 +127,8 @@ subroutine rad_initialize & ! --- input: integer, intent(in) :: levr, ictm, isol, ico2, iaer, num_p2d, & & ntcw, ialb, iems, num_p3d, npdf3d, ntoz, iovr_sw, iovr_lw, & - & isubc_sw, isubc_lw, icliq_sw, iflip, me, idate(4) + & isubc_sw, isubc_lw, icliq_sw, iflip, me, idate(4), & + & reiflag_rad real (kind=kind_phys), intent(in) :: si(levr+1) integer, intent(in) :: imp_physics @@ -180,6 +182,8 @@ subroutine rad_initialize & ivflip = iflip ! vertical index direction control flag + reiflagrad = reiflag_rad ! ice effective radius flag + ! --- assign initial permutation seed for mcica cloud-radiation if ( isubc_sw>0 .or. isubc_lw>0 ) then ! ipsd0 = 17*idate(1)+43*idate(2)+37*idate(3)+23*idate(4) + ipsd0 diff --git a/FV3/gfsphysics/physics/radiation_clouds.f b/FV3/gfsphysics/physics/radiation_clouds.f index c00628161..327d5790e 100644 --- a/FV3/gfsphysics/physics/radiation_clouds.f +++ b/FV3/gfsphysics/physics/radiation_clouds.f @@ -246,7 +246,8 @@ module module_radiation_clouds ! ! use physparam, only : icldflg, iovrsw, iovrlw, & & lcrick, lcnorm, lnoprec, & - & ivflip, kind_phys, kind_io4 + & ivflip, kind_phys, kind_io4, & + & reiflagrad use physcons, only : con_fvirt, con_ttp, con_rocp, & & con_t0c, con_pi, con_g, con_rd, & & con_thgni @@ -293,6 +294,8 @@ module module_radiation_clouds ! real (kind=kind_phys), parameter :: cldssa_def = 0.99 !> default cld asymmetry factor real (kind=kind_phys), parameter :: cldasy_def = 0.84 +!> default min value for cloud condensates (kg/kg) + real (kind=kind_phys), parameter :: qcmin = 1.0e-15 !> upper limit of boundary layer clouds integer :: llyr = 2 @@ -792,12 +795,25 @@ subroutine progcld1 & !! and McFarquhar (1996) \cite heymsfield_and_mcfarquhar_1996. if(.not.effr_in) then - call effective_ice_cloud_radius_heymsfield_mcfarquhar & + if (reiflagrad .eq. 1) then + call effective_ice_cloud_radius_heymsfield_mcfarquhar & & ( tlyr, cip, plyr, delp, tvly, & ! --- inputs: & IX, NLAY, & & rei & ! --- outputs: & ) + elseif (reiflagrad .eq. 4) then + call effective_ice_cloud_radius_kristjansson & + & ( tlyr, cip, plyr, delp, tvly, & ! --- inputs: + & IX, NLAY, & + & rei & ! --- outputs: + & ) + else + print *,'ERROR in reiflag_rad specification, valid options ', & + & 'are 1 and 4; reiflag_rad = ',reiflagrad + stop + endif endif + ! do k = 1, NLAY @@ -3283,8 +3299,6 @@ end subroutine gethml !----------------------------------- !! @} - - subroutine effective_ice_cloud_radius_heymsfield_mcfarquhar & & ( tlyr, cip, plyr, delp, tvly, & ! --- inputs: & IX, NLAY, & @@ -3347,7 +3361,62 @@ subroutine effective_ice_cloud_radius_heymsfield_mcfarquhar & end subroutine effective_ice_cloud_radius_heymsfield_mcfarquhar !----------------------------------- ! + subroutine effective_ice_cloud_radius_kristjansson & + & ( tlyr, cip, plyr, delp, tvly, & ! --- inputs: + & IX, NLAY, & + & rei & ! --- outputs: + & ) + ! =================================================================== ! + ! ! + ! abstract: compute cloud ice droplet effective radius, according to ! + ! the Kristjansson (2000) scheme ! + ! ! + ! ==================== definition of variables ==================== ! + ! ! + ! input variables: ! + ! tlyr (IX,NLAY) : model layer mean temperature in k ! + ! cip (IX,NLAY) : ice condensate path in g/m**2 ! + ! plyr (IX,NLAY) : model layer mean pressure in mb (100Pa) ! + ! delp (ix,nlay) : model layer pressure thickness in mb (100Pa) ! + ! tvly (IX,NLAY) : model layer virtual temperature in k ! + ! IX : horizontal dimention ! + ! NLAY : vertical layer dimensions ! + ! ! + ! output variables: ! + ! rei (IX,NLAY) : effective cloud ice droplet radius ! + + implicit none + + integer, intent(in) :: IX, NLAY + + real (kind=kind_phys), dimension(:,:), intent(in) :: tlyr, cip, & + & plyr, delp, tvly + + real (kind=kind_phys), dimension(:,:), intent(out) :: rei + + integer :: i, k + + real (kind=kind_phys) :: tem, tem2 + + ! do k = 1, NLAY + ! do i = 1, IX + ! if (qmi (i, k) .gt. qcmin) then + ! qci (i, k) = dpg * qmi (i, k) * 1.0e3 + ! ind = min (max (int (tlyr (i, k) - 136.0), 44), 138 - 1) + ! cor = tlyr (i, k) - int (tlyr (i, k)) + ! rei (i, k) = retab (ind) * (1. - cor) + retab (ind + 1) * cor + ! rei (i, k) = max (reimin, min (reimax, rei (i, k))) + ! else + ! qci (i, k) = 0.0 + ! rei (i, k) = reimin + ! endif + ! enddo + ! enddo +!................................... + end subroutine effective_ice_cloud_radius_kristjansson +!----------------------------------- + !........................................! end module module_radiation_clouds ! !========================================! -!> @} +!> @} \ No newline at end of file