diff --git a/physics/GFS_cloud_diagnostics.F90 b/physics/GFS_cloud_diagnostics.F90 index c62cc685d..1a7258b10 100644 --- a/physics/GFS_cloud_diagnostics.F90 +++ b/physics/GFS_cloud_diagnostics.F90 @@ -5,8 +5,9 @@ ! ######################################################################################## module GFS_cloud_diagnostics use machine, only: kind_phys - use physparam, only: iovrlw, iovrsw, ivflip, icldflg, idcor - + use physparam, only: icldflg + use module_radiation_clouds, only: gethml + ! Module parameters (imported directly from radiation_cloud.f) integer, parameter :: & NF_CLDS = 9, & ! Number of fields in cloud array @@ -23,7 +24,6 @@ module GFS_cloud_diagnostics ! Module variables integer :: & - iovr = 1, & ! Cloud overlap used for diagnostic HML cloud outputs llyr = 2 ! Upper limit of boundary layer clouds public GFS_cloud_diagnostics_run, GFS_cloud_diagnostics_init,& @@ -145,11 +145,8 @@ subroutine hml_cloud_diagnostics_initialize(imp_physics, imp_physics_fer_hires, ! Local variables integer :: iLay, kl - ! Initialize error flag - errflg = 0 - - ! Cloud overlap used for diagnostic HML cloud outputs - iovr = max(iovrsw,iovrlw) + ! Initialize error flag + errflg = 0 if (mpi_rank == 0) print *, VTAGCLD !print out version tag @@ -191,406 +188,4 @@ subroutine hml_cloud_diagnostics_initialize(imp_physics, imp_physics_fer_hires, return end subroutine hml_cloud_diagnostics_initialize - - ! ######################################################################################### - ! ######################################################################################### - subroutine gethml(plyr, ptop1, cldtot, cldcnv, dz, de_lgth, alpha, IX, NLAY, clds, mtop, mbot) - ! =================================================================== ! - ! ! - ! abstract: compute high, mid, low, total, and boundary cloud fractions ! - ! and cloud top/bottom layer indices for model diagnostic output. ! - ! the three cloud domain boundaries are defined by ptopc. the cloud ! - ! overlapping method is defined by control flag 'iovr', which is also ! - ! used by lw and sw radiation programs. ! - ! ! - ! usage: call gethml ! - ! ! - ! subprograms called: none ! - ! ! - ! attributes: ! - ! language: fortran 90 ! - ! machine: ibm-sp, sgi ! - ! ! - ! ! - ! ==================== definition of variables ==================== ! - ! ! - ! input variables: ! - ! plyr (IX,NLAY) : model layer mean pressure in mb (100Pa) ! - ! ptop1 (IX,4) : pressure limits of cloud domain interfaces ! - ! (sfc,low,mid,high) in mb (100Pa) ! - ! cldtot(IX,NLAY) : total or straiform cloud profile in fraction ! - ! cldcnv(IX,NLAY) : convective cloud (for diagnostic scheme only) ! - ! dz (ix,nlay) : layer thickness (km) ! - ! de_lgth(ix) : clouds vertical de-correlation length (km) ! - ! alpha(ix,nlay) : alpha decorrelation parameter ! - ! IX : horizontal dimention ! - ! NLAY : vertical layer dimensions ! - ! ! - ! output variables: ! - ! 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 ! - ! ! - ! external module variables: (in physparam) ! - ! ivflip : control flag of vertical index direction ! - ! =0: index from toa to surface ! - ! =1: index from surface to toa ! - ! ! - ! internal module variables: ! - ! iovr : control flag for cloud overlap ! - ! =0 random overlapping clouds ! - ! =1 max/ran overlapping clouds ! - ! =2 maximum overlapping ( for mcica only ) ! - ! =3 decorr-length ovlp ( for mcica only ) ! - ! =4 exponential cloud overlap (AER; mcica only) ! - ! =5 exponential-random overlap (AER; mcica only) ! - ! ! - ! ==================== end of description ===================== ! - ! - implicit none! - - ! --- inputs: - integer, intent(in) :: IX, NLAY - - real (kind=kind_phys), dimension(:,:), intent(in) :: plyr, ptop1, & - cldtot, cldcnv, dz - real (kind=kind_phys), dimension(:), intent(in) :: de_lgth - real (kind=kind_phys), dimension(:,:), intent(in) :: alpha - - ! --- outputs - real (kind=kind_phys), dimension(:,:), intent(out) :: clds - - integer, dimension(:,:), intent(out) :: mtop, mbot - - ! --- local variables: - real (kind=kind_phys) :: cl1(IX), cl2(IX), dz1(ix) - real (kind=kind_phys) :: pcur, pnxt, ccur, cnxt, alfa - - integer, dimension(IX):: idom, kbt1, kth1, kbt2, kth2 - integer :: i, k, id, id1, kstr, kend, kinc - - ! - !===> ... begin here - ! - clds(:,:) = 0.0 - - do i = 1, IX - cl1(i) = 1.0 - cl2(i) = 1.0 - enddo - - ! --- total and bl clouds, where cl1, cl2 are fractions of clear-sky view - ! layer processed from surface and up - - !> - Calculate total and BL cloud fractions (maximum-random cloud - !! overlapping is operational). - - if ( ivflip == 0 ) then ! input data from toa to sfc - kstr = NLAY - kend = 1 - kinc = -1 - else ! input data from sfc to toa - kstr = 1 - kend = NLAY - kinc = 1 - endif ! end_if_ivflip - - if ( iovr == 0 ) then ! random overlap - - do k = kstr, kend, kinc - do i = 1, IX - ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) )) - if (ccur >= climit) cl1(i) = cl1(i) * (1.0 - ccur) - enddo - - if (k == llyr) then - do i = 1, IX - clds(i,5) = 1.0 - cl1(i) ! save bl cloud - enddo - endif - enddo - - do i = 1, IX - clds(i,4) = 1.0 - cl1(i) ! save total cloud - enddo - - elseif ( iovr == 1 ) then ! max/ran overlap - - do k = kstr, kend, kinc - do i = 1, IX - ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) )) - if (ccur >= climit) then ! cloudy layer - cl2(i) = min( cl2(i), (1.0 - ccur) ) - else ! clear layer - cl1(i) = cl1(i) * cl2(i) - cl2(i) = 1.0 - endif - enddo - - if (k == llyr) then - do i = 1, IX - clds(i,5) = 1.0 - cl1(i) * cl2(i) ! save bl cloud - enddo - endif - enddo - - do i = 1, IX - clds(i,4) = 1.0 - cl1(i) * cl2(i) ! save total cloud - enddo - - elseif ( iovr == 2 ) then ! maximum overlap all levels - - cl1(:) = 0.0 - - do k = kstr, kend, kinc - do i = 1, IX - ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) )) - if (ccur >= climit) cl1(i) = max( cl1(i), ccur ) - enddo - - if (k == llyr) then - do i = 1, IX - clds(i,5) = cl1(i) ! save bl cloud - enddo - endif - enddo - - do i = 1, IX - clds(i,4) = cl1(i) ! save total cloud - enddo - - elseif ( iovr == 3 ) then ! random if clear-layer divided, - ! otherwise de-corrlength method - do i = 1, ix - dz1(i) = - dz(i,kstr) - enddo - - do k = kstr, kend, kinc - do i = 1, ix - ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) )) - if (ccur >= climit) then ! cloudy layer - alfa = exp( -0.5*((dz1(i)+dz(i,k)))/de_lgth(i) ) - dz1(i) = dz(i,k) - cl2(i) = alfa * min(cl2(i), (1.0 - ccur)) & ! maximum part - + (1.0 - alfa) * (cl2(i) * (1.0 - ccur)) ! random part - else ! clear layer - cl1(i) = cl1(i) * cl2(i) - cl2(i) = 1.0 - if (k /= kend) dz1(i) = -dz(i,k+kinc) - endif - enddo - - if (k == llyr) then - do i = 1, ix - clds(i,5) = 1.0 - cl1(i) * cl2(i) ! save bl cloud - enddo - endif - enddo - - do i = 1, ix - clds(i,4) = 1.0 - cl1(i) * cl2(i) ! save total cloud - enddo - - elseif ( iovr == 4 .or. iovr == 5 ) then ! exponential overlap (iovr=4), or - ! exponential-random (iovr=5); - ! distinction defined by alpha - do k = kstr, kend, kinc - do i = 1, ix - ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) )) - if (ccur >= climit) then ! cloudy layer - cl2(i) = alpha(i,k) * min(cl2(i), (1.0 - ccur)) & ! maximum part - + (1.0 - alpha(i,k)) * (cl2(i) * (1.0 - ccur)) ! random part - else ! clear layer - cl1(i) = cl1(i) * cl2(i) - cl2(i) = 1.0 - endif - enddo - if (k == llyr) then - do i = 1, ix - clds(i,5) = 1.0 - cl1(i) * cl2(i) ! save bl cloud - enddo - endif - enddo - do i = 1, ix - clds(i,4) = 1.0 - cl1(i) * cl2(i) ! save total cloud - enddo - endif ! end_if_iovr - - ! --- high, mid, low clouds, where cl1, cl2 are cloud fractions - ! layer processed from one layer below llyr and up - ! --- change! layer processed from surface to top, so low clouds will - ! contains both bl and low clouds. - - !> - Calculte high, mid, low cloud fractions and vertical indices of - !! cloud tops/bases. - if ( ivflip == 0 ) then ! input data from toa to sfc - - do i = 1, IX - cl1 (i) = 0.0 - cl2 (i) = 0.0 - kbt1(i) = NLAY - kbt2(i) = NLAY - kth1(i) = 0 - kth2(i) = 0 - idom(i) = 1 - mbot(i,1) = NLAY - mtop(i,1) = NLAY - mbot(i,2) = NLAY - 1 - mtop(i,2) = NLAY - 1 - mbot(i,3) = NLAY - 1 - mtop(i,3) = NLAY - 1 - enddo - - !org do k = llyr-1, 1, -1 - do k = NLAY, 1, -1 - do i = 1, IX - id = idom(i) - id1= id + 1 - - pcur = plyr(i,k) - ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) )) - - if (k > 1) then - pnxt = plyr(i,k-1) - cnxt = min( ovcst, max( cldtot(i,k-1), cldcnv(i,k-1) )) - else - pnxt = -1.0 - cnxt = 0.0 - endif - - if (pcur < ptop1(i,id1)) then - id = id + 1 - id1= id1 + 1 - idom(i) = id - endif - - if (ccur >= climit) then - if (kth2(i) == 0) kbt2(i) = k - kth2(i) = kth2(i) + 1 - - if ( iovr == 0 ) then - cl2(i) = cl2(i) + ccur - cl2(i)*ccur - else - cl2(i) = max( cl2(i), ccur ) - endif - - if (cnxt < climit .or. pnxt < ptop1(i,id1)) then - kbt1(i) = nint( (cl1(i)*kbt1(i) + cl2(i)*kbt2(i) ) & - / (cl1(i) + cl2(i)) ) - kth1(i) = nint( (cl1(i)*kth1(i) + cl2(i)*kth2(i) ) & - / (cl1(i) + cl2(i)) ) - cl1 (i) = cl1(i) + cl2(i) - cl1(i)*cl2(i) - - kbt2(i) = k - 1 - kth2(i) = 0 - cl2 (i) = 0.0 - endif ! end_if_cnxt_or_pnxt - endif ! end_if_ccur - - if (pnxt < ptop1(i,id1)) then - clds(i,id) = cl1(i) - mtop(i,id) = min( kbt1(i), kbt1(i)-kth1(i)+1 ) - mbot(i,id) = kbt1(i) - - cl1 (i) = 0.0 - kbt1(i) = k - 1 - kth1(i) = 0 - - if (id1 <= NK_CLDS) then - mbot(i,id1) = kbt1(i) - mtop(i,id1) = kbt1(i) - endif - endif ! end_if_pnxt - - enddo ! end_do_i_loop - enddo ! end_do_k_loop - - else ! input data from sfc to toa - - do i = 1, IX - cl1 (i) = 0.0 - cl2 (i) = 0.0 - kbt1(i) = 1 - kbt2(i) = 1 - kth1(i) = 0 - kth2(i) = 0 - idom(i) = 1 - mbot(i,1) = 1 - mtop(i,1) = 1 - mbot(i,2) = 2 - mtop(i,2) = 2 - mbot(i,3) = 2 - mtop(i,3) = 2 - enddo - - !org do k = llyr+1, NLAY - do k = 1, NLAY - do i = 1, IX - id = idom(i) - id1= id + 1 - - pcur = plyr(i,k) - ccur = min( ovcst, max( cldtot(i,k), cldcnv(i,k) )) - - if (k < NLAY) then - pnxt = plyr(i,k+1) - cnxt = min( ovcst, max( cldtot(i,k+1), cldcnv(i,k+1) )) - else - pnxt = -1.0 - cnxt = 0.0 - endif - - if (pcur < ptop1(i,id1)) then - id = id + 1 - id1= id1 + 1 - idom(i) = id - endif - - if (ccur >= climit) then - if (kth2(i) == 0) kbt2(i) = k - kth2(i) = kth2(i) + 1 - - if ( iovr == 0 ) then - cl2(i) = cl2(i) + ccur - cl2(i)*ccur - else - cl2(i) = max( cl2(i), ccur ) - endif - - if (cnxt < climit .or. pnxt < ptop1(i,id1)) then - kbt1(i) = nint( (cl1(i)*kbt1(i) + cl2(i)*kbt2(i)) & - / (cl1(i) + cl2(i)) ) - kth1(i) = nint( (cl1(i)*kth1(i) + cl2(i)*kth2(i)) & - / (cl1(i) + cl2(i)) ) - cl1 (i) = cl1(i) + cl2(i) - cl1(i)*cl2(i) - - kbt2(i) = k + 1 - kth2(i) = 0 - cl2 (i) = 0.0 - endif ! end_if_cnxt_or_pnxt - endif ! end_if_ccur - - if (pnxt < ptop1(i,id1)) then - clds(i,id) = cl1(i) - mtop(i,id) = max( kbt1(i), kbt1(i)+kth1(i)-1 ) - mbot(i,id) = kbt1(i) - - cl1 (i) = 0.0 - kbt1(i) = min(k+1, nlay) - kth1(i) = 0 - - if (id1 <= NK_CLDS) then - mbot(i,id1) = kbt1(i) - mtop(i,id1) = kbt1(i) - endif - endif ! end_if_pnxt - - enddo ! end_do_i_loop - enddo ! end_do_k_loop - - endif ! end_if_ivflip - - ! - return - !................................... - end subroutine gethml end module GFS_cloud_diagnostics diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index 2876f295d..eaa878ee7 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -3,7 +3,7 @@ type = scheme dependencies = funcphys.f90,iounitdef.f,machine.F,module_bfmicrophysics.f,module_mp_radar.F90,module_mp_thompson.F90 dependencies = module_mp_thompson_make_number_concentrations.F90,physcons.F90,physparam.f,radcons.f90,radiation_aerosols.f - dependencies = radiation_astronomy.f,radiation_clouds.f,radiation_gases.f,radlw_param.f,radsw_param.f,surface_perturbation.F90 + dependencies = radiation_astronomy.f,radiation_clouds.f,radiation_gases.f,radlw_param.f,radsw_param.f,surface_perturbation.F90,radiation_cloud_overlap.F90 ######################################################################## [ccpp-arg-table] diff --git a/physics/GFS_rrtmg_setup.F90 b/physics/GFS_rrtmg_setup.F90 index cc7e09c5e..2c18ab1e0 100644 --- a/physics/GFS_rrtmg_setup.F90 +++ b/physics/GFS_rrtmg_setup.F90 @@ -5,7 +5,7 @@ module GFS_rrtmg_setup use physparam, only : isolar , ictmflg, ico2flg, ioznflg, iaerflg,& ! & iaermdl, laswflg, lalwflg, lavoflg, icldflg, & & iaermdl, icldflg, & - & iovrsw , iovrlw , lcrick , lcnorm , lnoprec, & + & iovrRad=>iovr, lcrick , lcnorm , lnoprec, & & ialbflg, iemsflg, isubcsw, isubclw, ivflip , ipsd0, & & iswcliq, & & kind_phys @@ -45,7 +45,7 @@ module GFS_rrtmg_setup !! subroutine GFS_rrtmg_setup_init ( & si, levr, ictm, isol, ico2, iaer, ialb, iems, ntcw, num_p2d, & - num_p3d, npdf3d, ntoz, iovr_sw, iovr_lw, isubc_sw, isubc_lw, & + num_p3d, npdf3d, ntoz, iovr, isubc_sw, isubc_lw, & icliq_sw, crick_proof, ccnorm, & imp_physics, & norad_precip, idate, iflip, & @@ -131,7 +131,7 @@ subroutine GFS_rrtmg_setup_init ( & ! Stamnes(1993) \cite hu_and_stamnes_1993 method ! ! =2:cloud optical property scheme based on Hu and ! ! Stamnes(1993) -updated ! -! iovr_sw/iovr_lw : control flag for cloud overlap (sw/lw rad) ! +! iovr : control flag for cloud overlap (sw/lw rad) ! ! =0: random overlapping clouds ! ! =1: max/ran overlapping clouds ! ! =2: maximum overlap clouds (mcica only) ! @@ -178,8 +178,7 @@ subroutine GFS_rrtmg_setup_init ( & integer, intent(in) :: num_p3d integer, intent(in) :: npdf3d integer, intent(in) :: ntoz - integer, intent(in) :: iovr_sw - integer, intent(in) :: iovr_lw + integer, intent(in) :: iovr integer, intent(in) :: isubc_sw integer, intent(in) :: isubc_lw integer, intent(in) :: icliq_sw @@ -269,9 +268,10 @@ subroutine GFS_rrtmg_setup_init ( & iswcliq = icliq_sw ! optical property for liquid clouds for sw - iovrsw = iovr_sw ! cloud overlapping control flag for sw - iovrlw = iovr_lw ! cloud overlapping control flag for lw - + ! iovr comes from the model. In the RRTMG implementation this is stored in phyrparam.f, + ! it comes in from the host-model and is set here. + ! In GP, iovr is passed directly into the routines. + iovrRAD = iovr lcrick = crick_proof ! control flag for eliminating CRICK lcnorm = ccnorm ! control flag for in-cld condensate lnoprec = norad_precip ! precip effect on radiation flag (ferrier microphysics) @@ -294,8 +294,8 @@ subroutine GFS_rrtmg_setup_init ( & print *,' si =',si print *,' levr=',levr,' ictm=',ictm,' isol=',isol,' ico2=',ico2,& & ' iaer=',iaer,' ialb=',ialb,' iems=',iems,' ntcw=',ntcw - print *,' np3d=',num_p3d,' ntoz=',ntoz,' iovr_sw=',iovr_sw, & - & ' iovr_lw=',iovr_lw,' isubc_sw=',isubc_sw, & + print *,' np3d=',num_p3d,' ntoz=',ntoz, & + & ' iovr=',iovr,' isubc_sw=',isubc_sw, & & ' isubc_lw=',isubc_lw,' icliq_sw=',icliq_sw, & & ' iflip=',iflip,' me=',me print *,' crick_proof=',crick_proof, & @@ -468,8 +468,7 @@ subroutine radinit( si, NLAY, imp_physics, me ) ! =8 Thompson microphysics scheme ! ! =6 WSM6 microphysics scheme ! ! =10 MG microphysics scheme ! -! iovrsw : control flag for cloud overlap in sw radiation ! -! iovrlw : control flag for cloud overlap in lw radiation ! +! iovr : control flag for cloud overlap in radiation ! ! =0: random overlapping clouds ! ! =1: max/ran overlapping clouds ! ! isubcsw : sub-column cloud approx control flag in sw radiation ! @@ -545,10 +544,8 @@ subroutine radinit( si, NLAY, imp_physics, me ) & ' ISOLar =',isolar, ' ICO2flg=',ico2flg,' IAERflg=',iaerflg, & & ' IALBflg=',ialbflg,' IEMSflg=',iemsflg,' ICLDflg=',icldflg, & & ' IMP_PHYSICS=',imp_physics,' IOZNflg=',ioznflg - print *,' IVFLIP=',ivflip,' IOVRSW=',iovrsw,' IOVRLW=',iovrlw, & + print *,' IVFLIP=',ivflip,' IOVR=',iovrRad, & & ' ISUBCSW=',isubcsw,' ISUBCLW=',isubclw -! write(0,*)' IVFLIP=',ivflip,' IOVRSW=',iovrsw,' IOVRLW=',iovrlw,& -! & ' ISUBCSW=',isubcsw,' ISUBCLW=',isubclw print *,' LCRICK=',lcrick,' LCNORM=',lcnorm,' LNOPREC=',lnoprec print *,' LTP =',ltp,', add extra top layer =',lextop diff --git a/physics/GFS_rrtmg_setup.meta b/physics/GFS_rrtmg_setup.meta index e55147d66..b8d94db6c 100644 --- a/physics/GFS_rrtmg_setup.meta +++ b/physics/GFS_rrtmg_setup.meta @@ -112,17 +112,9 @@ type = integer intent = in optional = F -[iovr_sw] - standard_name = flag_for_cloud_overlap_method_for_shortwave_radiation - long_name = sw: max-random overlap clouds - units = flag - dimensions = () - type = integer - intent = in - optional = F -[iovr_lw] - standard_name = flag_for_cloud_overlap_method_for_longwave_radiation - long_name = lw: max-random overlap clouds +[iovr] + standard_name = flag_for_cloud_overlap_method_for_radiation + long_name = max-random overlap clouds units = flag dimensions = () type = integer diff --git a/physics/GFS_rrtmgp_gfdlmp_pre.F90 b/physics/GFS_rrtmgp_gfdlmp_pre.F90 index b67b22d41..52e1a7b74 100644 --- a/physics/GFS_rrtmgp_gfdlmp_pre.F90 +++ b/physics/GFS_rrtmgp_gfdlmp_pre.F90 @@ -4,9 +4,9 @@ ! ######################################################################################## module GFS_rrtmgp_gfdlmp_pre use machine, only: kind_phys - use physparam, only: lcnorm, lcrick, idcor, iovrlw, iovrsw use rrtmgp_aux, only: check_error_msg - use module_radiation_clouds, only: get_alpha_exp, get_alpha_dcorr + use module_radiation_cloud_overlap, only: cmp_dcorr_lgth, get_alpha_exp + ! Parameters real(kind_phys), parameter :: & reliq_def = 10.0 , & ! Default liq radius to 10 micron (used when effr_in=F) @@ -31,54 +31,64 @@ end subroutine GFS_rrtmgp_gfdlmp_pre_init !! \htmlinclude GFS_rrtmgp_gfdlmp_pre_run.html !! subroutine GFS_rrtmgp_gfdlmp_pre_run(nCol, nLev, nTracers, ncnd, i_cldliq, i_cldice, & - i_cldrain, i_cldsnow, i_cldgrpl, i_cldtot, yearlen, lsswr, lslwr, effr_in, julian,& - lat, p_lev, p_lay, tv_lay, effrin_cldliq, effrin_cldice, effrin_cldrain, & - effrin_cldsnow, tracer, con_pi, con_g, con_rd, con_epsq, & + i_cldrain, i_cldsnow, i_cldgrpl, i_cldtot, yearlen, doSWrad, doLWrad, effr_in, & + julian, lat, p_lev, p_lay, tv_lay, effrin_cldliq, effrin_cldice, effrin_cldrain, & + effrin_cldsnow, tracer, con_pi, con_g, con_rd, con_epsq, dcorr_con, idcor, iovr, & + iovr_dcorr, iovr_exprand, iovr_exp, idcor_con, idcor_hogan, idcor_oreopoulos, & cld_frac, cld_lwp, cld_reliq, cld_iwp, cld_reice, cld_swp, cld_resnow, cld_rwp, & cld_rerain, precip_frac, cloud_overlap_param, precip_overlap_param, de_lgth, & - deltaZ, errmsg, errflg) + deltaZb, errmsg, errflg) implicit none ! Inputs integer, intent(in) :: & - nCol, & ! Number of horizontal grid points - nLev, & ! Number of vertical layers - ncnd, & ! Number of cloud condensation types. - nTracers, & ! Number of tracers from model. - i_cldliq, & ! Index into tracer array for cloud liquid. - i_cldice, & ! Index into tracer array for cloud ice. - i_cldrain, & ! Index into tracer array for cloud rain. - i_cldsnow, & ! Index into tracer array for cloud snow. - i_cldgrpl, & ! Index into tracer array for cloud groupel. - i_cldtot, & ! Index into tracer array for cloud total amount. - yearlen ! Length of current year (365/366) WTF? + nCol, & ! Number of horizontal grid points + nLev, & ! Number of vertical layers + ncnd, & ! Number of cloud condensation types. + nTracers, & ! Number of tracers from model. + i_cldliq, & ! Index into tracer array for cloud liquid. + i_cldice, & ! Index into tracer array for cloud ice. + i_cldrain, & ! Index into tracer array for cloud rain. + i_cldsnow, & ! Index into tracer array for cloud snow. + i_cldgrpl, & ! Index into tracer array for cloud groupel. + i_cldtot, & ! Index into tracer array for cloud total amount. + yearlen, & ! Length of current year (365/366) WTF? + iovr, & ! Choice of cloud-overlap method + iovr_dcorr, & ! Flag for decorrelation-length cloud overlap method + iovr_exp, & ! Flag for exponential cloud overlap method + iovr_exprand, & ! Flag for exponential-random cloud overlap method + idcor, & ! Choice of method for decorrelation length computation + idcor_con, & ! Flag for decorrelation-length. Use constant value + idcor_hogan, & ! Flag for decorrelation-length. (https://rmets.onlinelibrary.wiley.com/doi/full/10.1002/qj.647) + idcor_oreopoulos ! Flag for decorrelation-length. (10.5194/acp-12-9097-2012) logical, intent(in) :: & - lsswr, & ! Call SW radiation? - lslwr, & ! Call LW radiation - effr_in ! Provide hydrometeor radii from macrophysics? + doSWrad, & ! Call SW radiation? + doLWrad, & ! Call LW radiation + effr_in ! Provide hydrometeor radii from macrophysics? real(kind_phys), intent(in) :: & - julian, & ! Julian day - con_pi, & ! Physical constant: pi - con_g, & ! Physical constant: gravitational constant - con_rd, & ! Physical constant: gas-constant for dry air - con_epsq ! Physical constant(?): Minimum value for specific humidity + julian, & ! Julian day + con_pi, & ! Physical constant: pi + con_g, & ! Physical constant: gravitational constant + con_rd, & ! Physical constant: gas-constant for dry air + con_epsq, & ! Physical constant(?): Minimum value for specific humidity + dcorr_con ! Decorrelation-length (used if idcor = idcor_con) real(kind_phys), dimension(nCol), intent(in) :: & - lat ! Latitude + lat ! Latitude real(kind_phys), dimension(nCol,nLev), intent(in) :: & - tv_lay, & ! Virtual temperature (K) - p_lay, & ! Pressure at model-layers (Pa) - effrin_cldliq, & ! Effective radius for liquid cloud-particles (microns) - effrin_cldice, & ! Effective radius for ice cloud-particles (microns) - effrin_cldrain, & ! Effective radius for rain cloud-particles (microns) - effrin_cldsnow ! Effective radius for snow cloud-particles (microns) + tv_lay, & ! Virtual temperature (K) + p_lay, & ! Pressure at model-layers (Pa) + effrin_cldliq, & ! Effective radius for liquid cloud-particles (microns) + effrin_cldice, & ! Effective radius for ice cloud-particles (microns) + effrin_cldrain, & ! Effective radius for rain cloud-particles (microns) + effrin_cldsnow ! Effective radius for snow cloud-particles (microns) real(kind_phys), dimension(nCol,nLev+1), intent(in) :: & - p_lev ! Pressure at model-level interfaces (Pa) + p_lev ! Pressure at model-level interfaces (Pa) real(kind_phys), dimension(nCol, nLev, nTracers),intent(in) :: & - tracer ! Cloud condensate amount in layer by type () + tracer ! Cloud condensate amount in layer by type () ! Outputs real(kind_phys), dimension(nCol),intent(out) :: & - de_lgth ! Decorrelation length + de_lgth ! Decorrelation length real(kind_phys), dimension(nCol,nLev),intent(out) :: & cld_frac, & ! Total cloud fraction cld_lwp, & ! Cloud liquid water path @@ -92,19 +102,22 @@ subroutine GFS_rrtmgp_gfdlmp_pre_run(nCol, nLev, nTracers, ncnd, i_cldliq, i_cld precip_frac, & ! Precipitation fraction cloud_overlap_param, & ! Cloud-overlap parameter precip_overlap_param, & ! Precipitation overlap parameter - deltaZ ! Layer thickness (km) + deltaZb ! Layer thickness (km) character(len=*), intent(out) :: & - errmsg ! Error message + errmsg ! Error message integer, intent(out) :: & - errflg ! Error flag + errflg ! Error flag ! Local variables - real(kind_phys) :: tem1 + real(kind_phys) :: tem1,pfac + real(kind_phys), dimension(nLev+1) :: hgtb + real(kind_phys), dimension(nLev) :: hgtc real(kind_phys), dimension(nCol, nLev, min(4,ncnd)) :: cld_condensate - integer :: iCol,iLay,l,ncndl,iovr - real(kind_phys), dimension(nCol,nLev) :: deltaP - - if (.not. (lsswr .or. lslwr)) return + integer :: iCol,iLay,l,ncndl,iSFC,iTOA + real(kind_phys), dimension(nCol,nLev) :: deltaP,deltaZ + logical :: top_at_1 + + if (.not. (doSWrad .or. doLWrad)) return ! Initialize CCPP error handling variables errmsg = '' @@ -117,30 +130,26 @@ subroutine GFS_rrtmgp_gfdlmp_pre_run(nCol, nLev, nTracers, ncnd, i_cldliq, i_cld call check_error_msg('GFS_rrtmgp_gfdlmp_pre_run',errmsg) return endif - ! - if (lcrick) then - errmsg = 'Namelist option lcrick is not supported.' - errflg = 1 - call check_error_msg('GFS_rrtmgp_gfdlmp_pre_run',errmsg) - return - endif - ! - if (lcnorm) then - errmsg = 'Namelist option lcnorm is not supported.' - errflg = 1 - call check_error_msg('GFS_rrtmgp_gfdlmp_pre_run',errmsg) - return + + ! What is vertical ordering? + top_at_1 = (p_lev(1,1) .lt. p_lev(1, nLev)) + if (top_at_1) then + iSFC = nLev + iTOA = 1 + else + iSFC = 1 + iTOA = nLev endif - ! Initialize outputs + ! Initialize outputs cld_lwp(:,:) = 0.0 - cld_reliq(:,:) = 0.0 + cld_reliq(:,:) = reliq_def cld_iwp(:,:) = 0.0 - cld_reice(:,:) = 0.0 + cld_reice(:,:) = reice_def cld_rwp(:,:) = 0.0 - cld_rerain(:,:) = 0.0 + cld_rerain(:,:) = rerain_def cld_swp(:,:) = 0.0 - cld_resnow(:,:) = 0.0 + cld_resnow(:,:) = resnow_def ! #################################################################################### ! Pull out cloud information for GFDL MP scheme. @@ -183,7 +192,7 @@ subroutine GFS_rrtmgp_gfdlmp_pre_run(nCol, nLev, nTracers, ncnd, i_cldliq, i_cld cld_swp(iCol,iLay) = cld_condensate(iCol,iLay,4) * tem1 endif ! Use radii provided from the macrophysics - if (effr_in) then + if (effr_in) then cld_reliq(iCol,iLay) = effrin_cldliq(iCol,iLay) cld_reice(iCol,iLay) = max(reice_min, min(reice_max,effrin_cldice(iCol,iLay))) cld_rerain(iCol,iLay) = effrin_cldrain(iCol,iLay) @@ -200,24 +209,88 @@ subroutine GFS_rrtmgp_gfdlmp_pre_run(nCol, nLev, nTracers, ncnd, i_cldliq, i_cld ! #################################################################################### ! Cloud (and precipitation) overlap ! #################################################################################### - - iovr = max(iovrsw,iovrlw) - - ! Compute layer-thickness - do iCol=1,nCol - do iLay=1,nLev - deltaZ(iCol,iLay) = ((con_rd/con_g)*0.001) * abs(log(p_lev(iCol,iLay)) - log(p_lev(iCol,iLay+1))) * tv_lay(iCol,iLay) - enddo + ! + ! Compute layer-thickness between layer boundaries (deltaZ) and layer centers (deltaZc) + ! + do iCol=1,nCol + if (top_at_1) then + ! Layer thickness (km) + do iLay=1,nLev + deltaZ(iCol,iLay) = ((con_rd/con_g)*0.001) * abs(log(p_lev(iCol,iLay+1)) - log(p_lev(iCol,iLay))) * tv_lay(iCol,iLay) + enddo + ! Height at layer boundaries + hgtb(nLev+1) = 0._kind_phys + do iLay=nLev,1,-1 + hgtb(iLay)= hgtb(iLay+1) + deltaZ(iCol,iLay) + enddo + ! Height at layer centers + do iLay = nLev, 1, -1 + pfac = abs(log(p_lev(iCol,iLay+1)) - log(p_lay(iCol,iLay))) / & + abs(log(p_lev(iCol,iLay+1)) - log(p_lev(iCol,iLay))) + hgtc(iLay) = hgtb(iLay+1) + pfac * (hgtb(iLay) - hgtb(iLay+1)) + enddo + ! Layer thickness between centers + do iLay = nLev-1, 1, -1 + deltaZb(iCol,iLay) = hgtc(iLay) - hgtc(iLay+1) + enddo + deltaZb(iCol,nLev) = hgtc(nLev) - hgtb(nLev+1) + else + do iLay=nLev,1,-1 + deltaZ(iCol,iLay) = ((con_rd/con_g)*0.001) * abs(log(p_lev(iCol,iLay)) - log(p_lev(iCol,iLay+1))) * tv_lay(iCol,iLay) + enddo + ! Height at layer boundaries + hgtb(1) = 0._kind_phys + do iLay=1,nLev + hgtb(iLay+1)= hgtb(iLay) + deltaZ(iCol,iLay) + enddo + ! Height at layer centers + do iLay = 1, nLev + pfac = abs(log(p_lev(iCol,iLay)) - log(p_lay(iCol,iLay) )) / & + abs(log(p_lev(iCol,iLay)) - log(p_lev(iCol,iLay+1))) + hgtc(iLay) = hgtb(iLay) + pfac * (hgtb(iLay+1) - hgtb(iLay)) + enddo + ! Layer thickness between centers + do iLay = 2, nLev + deltaZb(iCol,iLay) = hgtc(iLay) - hgtc(iLay-1) + enddo + deltaZb(iCol,1) = hgtc(1) - hgtb(1) + endif enddo - + + ! + ! Cloud decorrelation length + ! + if (idcor == idcor_hogan) then + call cmp_dcorr_lgth(nCol, lat, con_pi, de_lgth) + endif + if (idcor == idcor_oreopoulos) then + call cmp_dcorr_lgth(nCol, lat*(180._kind_phys/con_pi), julian, yearlen, de_lgth) + endif + if (idcor == idcor_con) then + de_lgth(:) = dcorr_con + endif + ! ! Cloud overlap parameter ! - if (iovr == 3) then - call get_alpha_dcorr(nCol, nLev, lat, con_pi, deltaZ, de_lgth, cloud_overlap_param) + if (iovr == iovr_dcorr .or. iovr == iovr_exp .or. iovr == iovr_exprand) then + call get_alpha_exp(nCol, nLev, deltaZb, de_lgth, cloud_overlap_param) + else + de_lgth(:) = 0. + cloud_overlap_param(:,:) = 0. endif - if (iovr == 4 .or. iovr == 5) then - call get_alpha_exp(nCol, nLev, deltaZ, iovr, lat, julian, yearlen, cld_frac, cloud_overlap_param) + + ! For exponential random overlap... + ! Decorrelate layers when a clear layer follows a cloudy layer to enforce + ! random correlation between non-adjacent blocks of cloudy layers + if (iovr == iovr_exprand) then + do iLay = 1, nLev + do iCol = 1, nCol + if (cld_frac(iCol,iLay) .eq. 0. .and. cld_frac(iCol,iLay-1) .gt. 0.) then + cloud_overlap_param(iCol,iLay) = 0._kind_phys + endif + enddo + enddo endif ! diff --git a/physics/GFS_rrtmgp_gfdlmp_pre.meta b/physics/GFS_rrtmgp_gfdlmp_pre.meta index 787879340..3841afc9b 100644 --- a/physics/GFS_rrtmgp_gfdlmp_pre.meta +++ b/physics/GFS_rrtmgp_gfdlmp_pre.meta @@ -39,7 +39,7 @@ type = integer intent = in optional = F -[lsswr] +[doSWrad] standard_name = flag_to_calc_sw long_name = logical flags for sw radiation calls units = flag @@ -47,7 +47,7 @@ type = logical intent = in optional = F -[lslwr] +[doLWrad] standard_name = flag_to_calc_lw long_name = logical flags for lw radiation calls units = flag @@ -245,6 +245,79 @@ kind = kind_phys intent = in optional = F +[iovr] + standard_name = flag_for_cloud_overlap_method_for_radiation + long_name = flag for cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_dcorr] + standard_name = flag_for_decorrelation_length_cloud_overlap_method + long_name = choice of decorrelation-length cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_exp] + standard_name = flag_for_exponential_cloud_overlap_method + long_name = choice of exponential cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_exprand] + standard_name = flag_for_exponential_random_cloud_overlap_method + long_name = choice of exponential-random cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F +[idcor] + standard_name = flag_for_decorrelation_length_method + long_name = flag for decorrelation length method used in cloud overlap method (iovr) + units = flag + dimensions = () + type = integer + intent = in + optional = F +[idcor_con] + standard_name = flag_for_constant_decorrelation_length_method + long_name = choice of decorrelation length computation (costant) + units = flag + dimensions = () + type = integer + intent = in + optional = F +[idcor_hogan] + standard_name = flag_for_hogan_decorrelation_length_method + long_name = choice of decorrelation length computation (hogan) + units = flag + dimensions = () + type = integer + intent = in + optional = F +[idcor_oreopoulos] + standard_name = flag_for_oreopoulos_decorrelation_length_method + long_name = choice of decorrelation length computation (oreopoulos) + units = flag + dimensions = () + type = integer + intent = in + optional = F +[dcorr_con] + standard_name = decorreltion_length_used_by_overlap_method + long_name = decorrelation length (default) used by cloud overlap method (iovr) + units = km + dimensions = () + type = real + intent = in + kind = kind_phys + optional = F [de_lgth] standard_name = cloud_decorrelation_length long_name = cloud decorrelation length @@ -362,7 +435,7 @@ kind = kind_phys intent = out optional = F -[deltaZ] +[deltaZb] standard_name = layer_thickness long_name = layer_thickness units = m diff --git a/physics/GFS_rrtmgp_setup.F90 b/physics/GFS_rrtmgp_setup.F90 index 9b503e3bc..a32f96ccf 100644 --- a/physics/GFS_rrtmgp_setup.F90 +++ b/physics/GFS_rrtmgp_setup.F90 @@ -1,634 +1,276 @@ !> \file GFS_rrtmgp_setup.f90 !! This file contains module GFS_rrtmgp_setup - - use physparam, only : & - isolar, ictmflg, ico2flg, ioznflg, iaerflg, iaermdl, icldflg, & - iovrsw, iovrlw, lcrick, lcnorm, lnoprec, ialbflg, iemsflg, & - isubcsw, isubclw, ivflip , ipsd0, iswcliq - use machine, only: & - kind_phys ! Working type - implicit none - - public GFS_rrtmgp_setup_init, GFS_rrtmgp_setup_run, GFS_rrtmgp_setup_finalize - - private - - logical :: is_initialized = .false. - - ! Version tag and last revision date - character(40), parameter :: & - VTAGRAD='NCEP-RRTMGP_driver v1.0 Sep 2019 ' - - ! Defaults - !> new data input control variables (set/reset in subroutines radinit/radupdate): - integer :: month0 = 0 - integer :: iyear0 = 0 - integer :: monthd = 0 - - !> control flag for the first time of reading climatological ozone data - !! (set/reset in subroutines radinit/radupdate, it is used only if the - !! control parameter ioznflg=0) - logical :: loz1st = .true. - - contains + use machine, only : kind_phys + use module_radiation_astronomy, only : sol_init, sol_update + use module_radiation_aerosols, only : aer_init, aer_update + use module_radiation_gases, only : gas_init, gas_update + use module_radiation_surface, only : sfc_init + use GFS_cloud_diagnostics, only : hml_cloud_diagnostics_initialize + ! *NOTE* These parameters below are required radiation_****** modules. They are not + ! directly used by the RRTMGP routines. + use physparam, only : isolar, ictmflg, ico2flg, ioznflg, iaerflg, & + iaermdl, ialbflg, iemsflg, ivflip + implicit none + + public GFS_rrtmgp_setup_init, GFS_rrtmgp_setup_run, GFS_rrtmgp_setup_finalize + + ! Version tag and last revision date + character(40), parameter :: & + VTAGRAD='NCEP-RRTMGP_driver v1.0 Sep 2019 ' + + ! Module paramaters + integer :: & + month0 = 0, & + iyear0 = 0, & + monthd = 0 + logical :: & + is_initialized = .false. + ! Control flag for the first time of reading climatological ozone data + ! (set/reset in subroutines GFS_rrtmgp_setup_init/GFS_rrtmgp_setuup_run, it is used only if + ! the control parameter ioznflg=0) + logical :: loz1st = .true. + +contains + ! ######################################################################################### + ! SUBROUTINE GFS_rrtmgp_setup_init + ! ######################################################################################### !> \defgroup GFS_rrtmgp_setup GFS RRTMGP Scheme Setup !! @{ !! \section arg_table_GFS_rrtmgp_setup_init !! \htmlinclude GFS_rrtmgp_setup_init.html !! - subroutine GFS_rrtmgp_setup_init(imp_physics, imp_physics_fer_hires, imp_physics_gfdl,& - imp_physics_thompson, imp_physics_wsm6, imp_physics_zhao_carr, & - imp_physics_zhao_carr_pdf, imp_physics_mg, si, levr, ictm, isol, ico2, iaer, & - ialb, iems, ntcw, num_p3d, ntoz, iovr_sw, iovr_lw, isubc_sw, isubc_lw, & - icliq_sw, crick_proof, ccnorm, norad_precip, idate, iflip, me, errmsg, errflg) - implicit none - - ! Inputs - integer, intent(in) :: & - imp_physics, & ! Flag for MP scheme - imp_physics_fer_hires, & ! Flag for fer-hires scheme - imp_physics_gfdl, & ! Flag for gfdl scheme - imp_physics_thompson, & ! Flag for thompsonscheme - imp_physics_wsm6, & ! Flag for wsm6 scheme - imp_physics_zhao_carr, & ! Flag for zhao-carr scheme - imp_physics_zhao_carr_pdf, & ! Flag for zhao-carr+PDF scheme - imp_physics_mg ! Flag for MG scheme - real(kind_phys), dimension(levr+1), intent(in) :: & - si - integer, intent(in) :: levr, ictm, isol, ico2, iaer, ialb, iems, & - ntcw, num_p3d, ntoz, iovr_sw, iovr_lw, isubc_sw, isubc_lw, & - icliq_sw, iflip, me - logical, intent(in) :: & - crick_proof, ccnorm, norad_precip - integer, intent(in), dimension(4) :: & - idate - ! Outputs - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg - - ! Initialize the CCPP error handling variables - errmsg = '' - errflg = 0 - if (is_initialized) return - - ! Set radiation parameters - isolar = isol ! solar constant control flag - ictmflg = ictm ! data ic time/date control flag - ico2flg = ico2 ! co2 data source control flag - ioznflg = ntoz ! ozone data source control flag - iswcliq = icliq_sw ! optical property for liquid clouds for sw - iovrsw = iovr_sw ! cloud overlapping control flag for sw - iovrlw = iovr_lw ! cloud overlapping control flag for lw - lcrick = crick_proof ! control flag for eliminating CRICK - lcnorm = ccnorm ! control flag for in-cld condensate - lnoprec = norad_precip ! precip effect on radiation flag (ferrier microphysics) - isubcsw = isubc_sw ! sub-column cloud approx flag in sw radiation - isubclw = isubc_lw ! sub-column cloud approx flag in lw radiation - ialbflg = ialb ! surface albedo control flag - iemsflg = iems ! surface emissivity control flag - ivflip = iflip ! vertical index direction control flag - - if ( ictm==0 .or. ictm==-2 ) then - iaerflg = mod(iaer, 100) ! no volcanic aerosols for clim hindcast - else - iaerflg = mod(iaer, 1000) - endif - iaermdl = iaer/1000 ! control flag for aerosol scheme selection - if ( iaermdl < 0 .or. (iaermdl>2 .and. iaermdl/=5) ) then - errmsg = trim(errmsg) // ' Error -- IAER flag is incorrect, Abort' - errflg = 1 - return - endif - - !if ( ntcw > 0 ) then - icldflg = 1 ! prognostic cloud optical prop scheme - !else - ! icldflg = 0 ! no support for diag cloud opt prop scheme - !endif - - ! Set 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) - endif - - if ( me == 0 ) then - print *,' In rad_initialize (GFS_rrtmgp_setup_init), before calling radinit' - print *,' si =',si - print *,' levr=',levr,' ictm=',ictm,' isol=',isol,' ico2=',ico2,& - ' iaer=',iaer,' ialb=',ialb,' iems=',iems,' ntcw=',ntcw - print *,' np3d=',num_p3d,' ntoz=',ntoz,' iovr_sw=',iovr_sw, & - ' iovr_lw=',iovr_lw,' isubc_sw=',isubc_sw, & - ' isubc_lw=',isubc_lw,' icliq_sw=',icliq_sw, & - ' iflip=',iflip,' me=',me - print *,' crick_proof=',crick_proof, & - ' ccnorm=',ccnorm,' norad_precip=',norad_precip - endif - - - call radinit( si, levr, imp_physics, imp_physics_fer_hires, imp_physics_gfdl, & - imp_physics_thompson, imp_physics_wsm6, imp_physics_zhao_carr, & - imp_physics_zhao_carr_pdf, imp_physics_mg, me, errflg ) - - if ( me == 0 ) then - print *,' Radiation sub-cloud initial seed =',ipsd0, & - ' IC-idate =',idate - print *,' return from rad_initialize (GFS_rrtmgp_setup_init) - after calling radinit' - endif - - is_initialized = .true. - return - end subroutine GFS_rrtmgp_setup_init - + subroutine GFS_rrtmgp_setup_init(imp_physics, imp_physics_fer_hires, imp_physics_gfdl, & + imp_physics_thompson, imp_physics_wsm6, imp_physics_zhao_carr, & + imp_physics_zhao_carr_pdf, imp_physics_mg, si, levr, ictm, isol, ico2, iaer, ialb, & + iems, ntcw, num_p3d, ntoz, iovr, isubc_sw, isubc_lw, icliq_sw, crick_proof, ccnorm, & + norad_precip, idate, iflip, me, errmsg, errflg) + + ! Inputs + integer, intent(in) :: & + imp_physics, & ! Flag for MP scheme + imp_physics_fer_hires, & ! Flag for fer-hires scheme + imp_physics_gfdl, & ! Flag for gfdl scheme + imp_physics_thompson, & ! Flag for thompsonscheme + imp_physics_wsm6, & ! Flag for wsm6 scheme + imp_physics_zhao_carr, & ! Flag for zhao-carr scheme + imp_physics_zhao_carr_pdf, & ! Flag for zhao-carr+PDF scheme + imp_physics_mg ! Flag for MG scheme + real(kind_phys), dimension(levr+1), intent(in) :: & + si + integer, intent(in) :: levr, ictm, isol, ico2, iaer, ialb, iems, & + ntcw, num_p3d, ntoz, iovr, isubc_sw, isubc_lw, & + icliq_sw, iflip, me + logical, intent(in) :: & + crick_proof, ccnorm, norad_precip + integer, intent(in), dimension(4) :: & + idate + + ! Outputs + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + + if (is_initialized) return + + ! Set radiation parameters + isolar = isol ! solar constant control flag + ictmflg = ictm ! data ic time/date control flag + ico2flg = ico2 ! co2 data source control flag + ioznflg = ntoz ! ozone data source control flag + ialbflg = ialb ! surface albedo control flag + iemsflg = iems ! surface emissivity control flag + ivflip = iflip ! vertical index direction control flag + + if ( ictm==0 .or. ictm==-2 ) then + iaerflg = mod(iaer, 100) ! no volcanic aerosols for clim hindcast + else + iaerflg = mod(iaer, 1000) + endif + iaermdl = iaer/1000 ! control flag for aerosol scheme selection + if ( iaermdl < 0 .or. (iaermdl>2 .and. iaermdl/=5) ) then + errmsg = trim(errmsg) // ' Error -- IAER flag is incorrect, Abort' + errflg = 1 + return + endif + + if ( me == 0 ) then + print *,' In rad_initialize (GFS_rrtmgp_setup_init), before calling radinit' + print *,' si = ',si + print *,' levr = ',levr, & + ' ictm = ',ictm, & + ' isol = ',isol, & + ' ico2 = ',ico2, & + ' iaer = ',iaer, & + ' ialb = ',ialb, & + ' iems = ',iems, & + ' ntcw = ',ntcw + print *,' np3d = ',num_p3d, & + ' ntoz = ',ntoz, & + ' iovr = ',iovr, & + ' isubc_sw = ',isubc_sw, & + ' isubc_lw = ',isubc_lw, & + ' icliq_sw = ',icliq_sw, & + ' iflip = ',iflip, & + ' me = ',me + endif + +#if 0 + ! GFS_radiation_driver.F90 may in the future initialize air/ground + ! temperature differently; however, this is not used at the moment + ! and as such we avoid the difficulty of dealing with exchanging + ! itsfc between GFS_rrtmgp_setup and a yet-to-be-created/-used + ! interstitial routine (or GFS_radiation_driver.F90) + itsfc = iemsflg / 10 ! sfc air/ground temp control +#endif + loz1st = (ioznflg == 0) ! first-time clim ozone data read flag + month0 = 0 + iyear0 = 0 + monthd = 0 + + ! Call initialization routines.. + call sol_init ( me ) + call aer_init ( levr, me ) + call gas_init ( me ) + call sfc_init ( me ) + call hml_cloud_diagnostics_initialize(imp_physics, imp_physics_fer_hires, & + imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, & + imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, imp_physics_mg, levr, me, si,& + errflg) + + if ( me == 0 ) then + print *,' return from rad_initialize (GFS_rrtmgp_setup_init) - after calling radinit' + endif + + is_initialized = .true. + + return + end subroutine GFS_rrtmgp_setup_init + + ! ######################################################################################### + ! SUBROUTINE GFS_rrtmgp_setup_run + ! ######################################################################################### !> \section arg_table_GFS_rrtmgp_setup_run !! \htmlinclude GFS_rrtmgp_setup_run.html !! - subroutine GFS_rrtmgp_setup_run (idate, jdate, deltsw, deltim, lsswr, me, & - slag, sdec, cdec, solcon, errmsg, errflg) - - implicit none - - ! interface variables - integer, intent(in) :: idate(:) - integer, intent(in) :: jdate(:) - real(kind=kind_phys), intent(in) :: deltsw - real(kind=kind_phys), intent(in) :: deltim - logical, intent(in) :: lsswr - integer, intent(in) :: me - real(kind=kind_phys), intent(out) :: slag - real(kind=kind_phys), intent(out) :: sdec - real(kind=kind_phys), intent(out) :: cdec - real(kind=kind_phys), intent(out) :: solcon - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg - - ! Check initialization state - if (.not.is_initialized) then - write(errmsg, fmt='((a))') 'GFS_rrtmgp_setup_run called before GFS_rrtmgp_setup_init' - errflg = 1 - return - end if - - ! Initialize the CCPP error handling variables - errmsg = '' - errflg = 0 + subroutine GFS_rrtmgp_setup_run (idate, jdate, deltsw, deltim, lsswr, me, & + slag, sdec, cdec, solcon, errmsg, errflg) - call radupdate(idate,jdate,deltsw,deltim,lsswr,me, & - slag,sdec,cdec,solcon) - - end subroutine GFS_rrtmgp_setup_run - + ! Inputs + integer, intent(in) :: idate(:) + integer, intent(in) :: jdate(:) + real(kind_phys), intent(in) :: deltsw + real(kind_phys), intent(in) :: deltim + logical, intent(in) :: lsswr + integer, intent(in) :: me + + ! Outputs + real(kind_phys), intent(out) :: slag + real(kind_phys), intent(out) :: sdec + real(kind_phys), intent(out) :: cdec + real(kind_phys), intent(out) :: solcon + character(len=*),intent(out) :: errmsg + integer, intent(out) :: errflg + + ! Locals + integer :: iyear, imon, iday, ihour + integer :: kyear, kmon, kday, khour + logical :: lmon_chg ! month change flag + logical :: lco2_chg ! cntrl flag for updating co2 data + logical :: lsol_chg ! cntrl flag for updating solar constant + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + + ! Check initialization state + if (.not.is_initialized) then + write(errmsg, fmt='((a))') 'GFS_rrtmgp_setup_run called before GFS_rrtmgp_setup_init' + errflg = 1 + return + end if + + ! Set up time stamp at fcst time and that for green house gases + iyear = jdate(1) + imon = jdate(2) + iday = jdate(3) + ihour = jdate(5) + + ! Set up time stamp used for green house gases (** currently co2 only) + ! get external data at initial condition time + if ( ictmflg==0 .or. ictmflg==-2 ) then + kyear = idate(1) + kmon = idate(2) + kday = idate(3) + khour = idate(5) + ! get external data at fcst or specified time + else + kyear = iyear + kmon = imon + kday = iday + khour = ihour + endif + + if ( month0 /= imon ) then + lmon_chg = .true. + month0 = imon + else + lmon_chg = .false. + endif + + ! Update solar forcing... + if (lsswr) then + if ( isolar == 0 .or. isolar == 10 ) then + lsol_chg = .false. + elseif ( iyear0 /= iyear ) then + lsol_chg = .true. + else + lsol_chg = ( isolar==4 .and. lmon_chg ) + endif + iyear0 = iyear + call sol_update(jdate, kyear, deltsw, deltim, lsol_chg, me, slag, sdec, cdec, solcon) + endif + + ! Update aerosols... + if ( lmon_chg ) then + call aer_update ( iyear, imon, me ) + endif + + ! Update trace gases (co2 only)... + if ( monthd /= kmon ) then + monthd = kmon + lco2_chg = .true. + else + lco2_chg = .false. + endif + call gas_update (kyear, kmon, kday, khour, loz1st, lco2_chg, me ) + + if ( loz1st ) loz1st = .false. + + return + end subroutine GFS_rrtmgp_setup_run + + ! ######################################################################################### + ! SUBROUTINE GFS_rrtmgp_setup_finalize + ! ######################################################################################### !> \section arg_table_GFS_rrtmgp_setup_finalize !! \htmlinclude GFS_rrtmgp_setup_finalize.html !! - subroutine GFS_rrtmgp_setup_finalize (errmsg, errflg) + subroutine GFS_rrtmgp_setup_finalize (errmsg, errflg) + character(len=*), intent( out) :: errmsg + integer, intent( out) :: errflg + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 - implicit none - - character(len=*), intent( out) :: errmsg - integer, intent( out) :: errflg - - ! Initialize the CCPP error handling variables - errmsg = '' - errflg = 0 - - if (.not.is_initialized) return - - ! do finalization stuff if needed - - is_initialized = .false. - - end subroutine GFS_rrtmgp_setup_finalize - - - ! Private functions - - - subroutine radinit(si, NLAY, imp_physics, imp_physics_fer_hires, imp_physics_gfdl, & - imp_physics_thompson, imp_physics_wsm6, imp_physics_zhao_carr, & - imp_physics_zhao_carr_pdf, imp_physics_mg, me, errflg ) - !................................... - -! --- inputs: -! & ( si, NLAY, imp_physics, imp_physics_fer_hires, imp_physics_gfdl, & -! & imp_physics_thompson, imp_physics_wsm6, imp_physics_zhao_carr, & -! & imp_physics_zhao_carr_pdf, imp_physics_mg, me ) -! --- outputs: -! ( errflg ) - -! ================= subprogram documentation block ================ ! -! ! -! subprogram: radinit initialization of radiation calculations ! -! ! -! usage: call radinit ! -! ! -! attributes: ! -! language: fortran 90 ! -! machine: wcoss ! -! ! -! ==================== definition of variables ==================== ! -! ! -! input parameters: ! -! si : model vertical sigma interface ! -! NLAY : number of model vertical layers ! -! imp_physics : MP identifier ! -! me : print control flag ! -! ! -! outputs: (none) ! -! ! -! external module variables: (in module physparam) ! -! isolar : solar constant cntrol flag ! -! = 0: use the old fixed solar constant in "physcon" ! -! =10: use the new fixed solar constant in "physcon" ! -! = 1: use noaa ann-mean tsi tbl abs-scale with cycle apprx! -! = 2: use noaa ann-mean tsi tbl tim-scale with cycle apprx! -! = 3: use cmip5 ann-mean tsi tbl tim-scale with cycl apprx! -! = 4: use cmip5 mon-mean tsi tbl tim-scale with cycl apprx! -! iaerflg : 3-digit aerosol flag (abc for volc, lw, sw) ! -! a:=0 use background stratospheric aerosol ! -! =1 include stratospheric vocanic aeros ! -! b:=0 no topospheric aerosol in lw radiation ! -! =1 compute tropspheric aero in 1 broad band for lw ! -! =2 compute tropspheric aero in multi bands for lw ! -! c:=0 no topospheric aerosol in sw radiation ! -! =1 include tropspheric aerosols for sw ! -! ico2flg : co2 data source control flag ! -! =0: use prescribed global mean co2 (old oper) ! -! =1: use observed co2 annual mean value only ! -! =2: use obs co2 monthly data with 2-d variation ! -! ictmflg : =yyyy#, external data ic time/date control flag ! -! = -2: same as 0, but superimpose seasonal cycle ! -! from climatology data set. ! -! = -1: use user provided external data for the ! -! forecast time, no extrapolation. ! -! = 0: use data at initial cond time, if not ! -! available, use latest, no extrapolation. ! -! = 1: use data at the forecast time, if not ! -! available, use latest and extrapolation. ! -! =yyyy0: use yyyy data for the forecast time, ! -! no further data extrapolation. ! -! =yyyy1: use yyyy data for the fcst. if needed, do ! -! extrapolation to match the fcst time. ! -! ioznflg : ozone data source control flag ! -! =0: use climatological ozone profile ! -! =1: use interactive ozone profile ! -! ialbflg : albedo scheme control flag ! -! =0: climatology, based on surface veg types ! -! =1: modis retrieval based surface albedo scheme ! -! iemsflg : emissivity scheme cntrl flag (ab 2-digit integer) ! -! a:=0 set sfc air/ground t same for lw radiation ! -! =1 set sfc air/ground t diff for lw radiation ! -! b:=0 use fixed sfc emissivity=1.0 (black-body) ! -! =1 use varying climtology sfc emiss (veg based) ! -! =2 future development (not yet) ! -! icldflg : cloud optical property scheme control flag ! -! =0: use diagnostic cloud scheme ! -! =1: use prognostic cloud scheme (default) ! -! imp_physics : cloud microphysics scheme control flag ! -! =99 zhao/carr/sundqvist microphysics scheme ! -! =98 zhao/carr/sundqvist microphysics+pdf cloud&cnvc,cnvw ! -! =11 GFDL cloud microphysics ! -! =8 Thompson microphysics scheme ! -! =6 WSM6 microphysics scheme ! -! =10 MG microphysics scheme ! -! iovrsw : control flag for cloud overlap in sw radiation ! -! iovrlw : control flag for cloud overlap in lw radiation ! -! =0: random overlapping clouds ! -! =1: max/ran overlapping clouds ! -! isubcsw : sub-column cloud approx control flag in sw radiation ! -! isubclw : sub-column cloud approx control flag in lw radiation ! -! =0: with out sub-column cloud approximation ! -! =1: mcica sub-col approx. prescribed random seed ! -! =2: mcica sub-col approx. provided random seed ! -! 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 ! -! lnoprec : precip effect in radiation flag (ferrier microphysics) ! -! =t: snow/rain has no impact on radiation ! -! =f: snow/rain has impact on radiation ! -! ivflip : vertical index direction control flag ! -! =0: index from toa to surface ! -! =1: index from surface to toa ! -! ! -! subroutines called: sol_init, aer_init, gas_init, cld_init, ! -! sfc_init, rlwinit, rswinit ! -! ! -! usage: call radinit ! -! ! -! =================================================================== ! -! - - use module_radiation_astronomy, only : sol_init - use module_radiation_aerosols, only : aer_init - use module_radiation_gases, only : gas_init - use module_radiation_surface, only : sfc_init - use GFS_cloud_diagnostics, only : hml_cloud_diagnostics_initialize - - implicit none - -! --- inputs: - integer, intent(in) :: & - imp_physics, & ! Flag for MP scheme - imp_physics_fer_hires, & ! Flag for fer-hires scheme - imp_physics_gfdl, & ! Flag for gfdl scheme - imp_physics_thompson, & ! Flag for thompsonscheme - imp_physics_wsm6, & ! Flag for wsm6 scheme - imp_physics_zhao_carr, & ! Flag for zhao-carr scheme - imp_physics_zhao_carr_pdf, & ! Flag for zhao-carr+PDF scheme - imp_physics_mg ! Flag for MG scheme - integer, intent(in) :: NLAY, me - real (kind=kind_phys), intent(in) :: si(:) - -! --- outputs: (none, to module variables) - integer, intent(out) :: & - errflg - -! --- locals: - - ! Initialize - errflg = 0 -! -!===> ... begin here -! -!> -# Set up control variables and external module variables in -!! module physparam -#if 0 - ! GFS_radiation_driver.F90 may in the future initialize air/ground - ! temperature differently; however, this is not used at the moment - ! and as such we avoid the difficulty of dealing with exchanging - ! itsfc between GFS_rrtmgp_setup and a yet-to-be-created/-used - ! interstitial routine (or GFS_radiation_driver.F90) - itsfc = iemsflg / 10 ! sfc air/ground temp control -#endif - loz1st = (ioznflg == 0) ! first-time clim ozone data read flag - month0 = 0 - iyear0 = 0 - monthd = 0 - - if (me == 0) then -! print *,' NEW RADIATION PROGRAM STRUCTURES -- SEP 01 2004' - print *,' NEW RADIATION PROGRAM STRUCTURES BECAME OPER. ', & - & ' May 01 2007' - print *, VTAGRAD !print out version tag - print *,' - Selected Control Flag settings: ICTMflg=',ictmflg, & - & ' ISOLar =',isolar, ' ICO2flg=',ico2flg,' IAERflg=',iaerflg, & - & ' IALBflg=',ialbflg,' IEMSflg=',iemsflg,' ICLDflg=',icldflg, & - & ' IMP_PHYSICS=',imp_physics,' IOZNflg=',ioznflg - print *,' IVFLIP=',ivflip,' IOVRSW=',iovrsw,' IOVRLW=',iovrlw, & - & ' ISUBCSW=',isubcsw,' ISUBCLW=',isubclw -! write(0,*)' IVFLIP=',ivflip,' IOVRSW=',iovrsw,' IOVRLW=',iovrlw,& -! & ' ISUBCSW=',isubcsw,' ISUBCLW=',isubclw - print *,' LCRICK=',lcrick,' LCNORM=',lcnorm,' LNOPREC=',lnoprec - - if ( ictmflg==0 .or. ictmflg==-2 ) then - print *,' Data usage is limited by initial condition!' - print *,' No volcanic aerosols' - endif - - if ( isubclw == 0 ) then - print *,' - ISUBCLW=',isubclw,' No McICA, use grid ', & - & 'averaged cloud in LW radiation' - elseif ( isubclw == 1 ) then - print *,' - ISUBCLW=',isubclw,' Use McICA with fixed ', & - & 'permutation seeds for LW random number generator' - elseif ( isubclw == 2 ) then - print *,' - ISUBCLW=',isubclw,' Use McICA with random ', & - & 'permutation seeds for LW random number generator' - else - print *,' - ERROR!!! ISUBCLW=',isubclw,' is not a ', & - & 'valid option ' - stop - endif - - if ( isubcsw == 0 ) then - print *,' - ISUBCSW=',isubcsw,' No McICA, use grid ', & - & 'averaged cloud in SW radiation' - elseif ( isubcsw == 1 ) then - print *,' - ISUBCSW=',isubcsw,' Use McICA with fixed ', & - & 'permutation seeds for SW random number generator' - elseif ( isubcsw == 2 ) then - print *,' - ISUBCSW=',isubcsw,' Use McICA with random ', & - & 'permutation seeds for SW random number generator' - else - print *,' - ERROR!!! ISUBCSW=',isubcsw,' is not a ', & - & 'valid option ' - stop - endif - - if ( isubcsw /= isubclw ) then - print *,' - *** Notice *** ISUBCSW /= ISUBCLW !!!', & - & isubcsw, isubclw - endif - endif - - ! Initialization - - call sol_init ( me ) ! --- ... astronomy initialization routine - call aer_init ( NLAY, me ) ! --- ... aerosols initialization routine - call gas_init ( me ) ! --- ... co2 and other gases initialization routine - call sfc_init ( me ) ! --- ... surface initialization routine - call hml_cloud_diagnostics_initialize(imp_physics, imp_physics_fer_hires, & - imp_physics_gfdl, imp_physics_thompson, imp_physics_wsm6, & - imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, imp_physics_mg, NLAY, me, si,& - errflg) - - return - !................................... - end subroutine radinit - !----------------------------------- - -!> This subroutine checks and updates time sensitive data used by -!! radiation computations. This subroutine needs to be placed inside -!! the time advancement loop but outside of the horizontal grid loop. -!! It is invoked at radiation calling frequncy but before any actual -!! radiative transfer computations. -!! \param idate NCEP absolute date and time of intial condition -!! (year,month,day,time-zone,hour,minute,second, -!! mil-second) -!! \param jdate NCEP absolute date and time at forecast time -!! (year,month,day,time-zone,hour,minute,second, -!! mil-second) -!! \param deltsw SW radiation calling time interval in seconds -!! \param deltim model advancing time-step duration in seconds -!! \param lsswr logical control flag for SW radiation calculations -!! \param me print control flag -!! \param slag equation of time in radians -!! \param sdec,cdec sine and cosine of the solar declination angle -!! \param solcon solar constant adjusted by sun-earth distance \f$(W/m^2)\f$ -!> \section gen_radupdate General Algorithm -!> @{ -!----------------------------------- - subroutine radupdate( idate,jdate,deltsw,deltim,lsswr, me, & - & slag,sdec,cdec,solcon) -!................................... - -! ================= subprogram documentation block ================ ! -! ! -! subprogram: radupdate calls many update subroutines to check and ! -! update radiation required but time varying data sets and module ! -! variables. ! -! ! -! usage: call radupdate ! -! ! -! attributes: ! -! language: fortran 90 ! -! machine: ibm sp ! -! ! -! ==================== definition of variables ==================== ! -! ! -! input parameters: ! -! idate(8) : ncep absolute date and time of initial condition ! -! (yr, mon, day, t-zone, hr, min, sec, mil-sec) ! -! jdate(8) : ncep absolute date and time at fcst time ! -! (yr, mon, day, t-zone, hr, min, sec, mil-sec) ! -! deltsw : sw radiation calling frequency in seconds ! -! deltim : model timestep in seconds ! -! lsswr : logical flags for sw radiation calculations ! -! me : print control flag ! -! ! -! outputs: ! -! slag : equation of time in radians ! -! sdec, cdec : sin and cos of the solar declination angle ! -! solcon : sun-earth distance adjusted solar constant (w/m2) ! -! ! -! external module variables: ! -! isolar : solar constant cntrl (in module physparam) ! -! = 0: use the old fixed solar constant in "physcon" ! -! =10: use the new fixed solar constant in "physcon" ! -! = 1: use noaa ann-mean tsi tbl abs-scale with cycle apprx! -! = 2: use noaa ann-mean tsi tbl tim-scale with cycle apprx! -! = 3: use cmip5 ann-mean tsi tbl tim-scale with cycl apprx! -! = 4: use cmip5 mon-mean tsi tbl tim-scale with cycl apprx! -! ictmflg : =yyyy#, external data ic time/date control flag ! -! = -2: same as 0, but superimpose seasonal cycle ! -! from climatology data set. ! -! = -1: use user provided external data for the ! -! forecast time, no extrapolation. ! -! = 0: use data at initial cond time, if not ! -! available, use latest, no extrapolation. ! -! = 1: use data at the forecast time, if not ! -! available, use latest and extrapolation. ! -! =yyyy0: use yyyy data for the forecast time, ! -! no further data extrapolation. ! -! =yyyy1: use yyyy data for the fcst. if needed, do ! -! extrapolation to match the fcst time. ! -! ! -! module variables: ! -! loz1st : first-time clim ozone data read flag ! -! ! -! subroutines called: sol_update, aer_update, gas_update ! -! ! -! =================================================================== ! -! - use module_radiation_astronomy, only : sol_update - use module_radiation_aerosols, only : aer_update - use module_radiation_gases, only : gas_update - - implicit none - -! --- inputs: - integer, intent(in) :: idate(:), jdate(:), me - logical, intent(in) :: lsswr - - real (kind=kind_phys), intent(in) :: deltsw, deltim - -! --- outputs: - real (kind=kind_phys), intent(out) :: slag, sdec, cdec, solcon - -! --- locals: - integer :: iyear, imon, iday, ihour - integer :: kyear, kmon, kday, khour - - logical :: lmon_chg ! month change flag - logical :: lco2_chg ! cntrl flag for updating co2 data - logical :: lsol_chg ! cntrl flag for updating solar constant -! -!===> ... begin here -! -!> -# Set up time stamp at fcst time and that for green house gases -!! (currently co2 only) -! --- ... time stamp at fcst time - - iyear = jdate(1) - imon = jdate(2) - iday = jdate(3) - ihour = jdate(5) - -! --- ... set up time stamp used for green house gases (** currently co2 only) - - if ( ictmflg==0 .or. ictmflg==-2 ) then ! get external data at initial condition time - kyear = idate(1) - kmon = idate(2) - kday = idate(3) - khour = idate(5) - else ! get external data at fcst or specified time - kyear = iyear - kmon = imon - kday = iday - khour = ihour - endif ! end if_ictmflg_block - - if ( month0 /= imon ) then - lmon_chg = .true. - month0 = imon - else - lmon_chg = .false. - endif - -!> -# Call module_radiation_astronomy::sol_update(), yearly update, no -!! time interpolation. - if (lsswr) then - - if ( isolar == 0 .or. isolar == 10 ) then - lsol_chg = .false. - elseif ( iyear0 /= iyear ) then - lsol_chg = .true. - else - lsol_chg = ( isolar==4 .and. lmon_chg ) - endif - iyear0 = iyear - - call sol_update & -! --- inputs: - & ( jdate,kyear,deltsw,deltim,lsol_chg, me, & -! --- outputs: - & slag,sdec,cdec,solcon & - & ) - - endif ! end_if_lsswr_block - -!> -# Call module_radiation_aerosols::aer_update(), monthly update, no -!! time interpolation - if ( lmon_chg ) then - call aer_update ( iyear, imon, me ) - endif - -!> -# Call co2 and other gases update routine: -!! module_radiation_gases::gas_update() - if ( monthd /= kmon ) then - monthd = kmon - lco2_chg = .true. - else - lco2_chg = .false. - endif - - call gas_update ( kyear,kmon,kday,khour,loz1st,lco2_chg, me ) - - if ( loz1st ) loz1st = .false. - -!> -# Call surface update routine (currently not needed) -! call sfc_update ( iyear, imon, me ) - -!> -# Call clouds update routine (currently not needed) -! call cld_update ( iyear, imon, me ) -! - return -!................................... - end subroutine radupdate -!----------------------------------- - -!! @} + if (.not.is_initialized) return + + ! do finalization stuff if needed + is_initialized = .false. + + end subroutine GFS_rrtmgp_setup_finalize end module GFS_rrtmgp_setup diff --git a/physics/GFS_rrtmgp_setup.meta b/physics/GFS_rrtmgp_setup.meta index 2e656e259..fb31f5c7a 100644 --- a/physics/GFS_rrtmgp_setup.meta +++ b/physics/GFS_rrtmgp_setup.meta @@ -160,17 +160,9 @@ type = integer intent = in optional = F -[iovr_sw] - standard_name = flag_for_cloud_overlap_method_for_shortwave_radiation - long_name = sw: max-random overlap clouds - units = flag - dimensions = () - type = integer - intent = in - optional = F -[iovr_lw] - standard_name = flag_for_cloud_overlap_method_for_longwave_radiation - long_name = lw: max-random overlap clouds +[iovr] + standard_name = flag_for_cloud_overlap_method_for_radiation + long_name = max-random overlap clouds units = flag dimensions = () type = integer diff --git a/physics/GFS_rrtmgp_sw_pre.F90 b/physics/GFS_rrtmgp_sw_pre.F90 index f6aac60b1..179c622f5 100644 --- a/physics/GFS_rrtmgp_sw_pre.F90 +++ b/physics/GFS_rrtmgp_sw_pre.F90 @@ -1,5 +1,4 @@ module GFS_rrtmgp_sw_pre - use physparam use machine, only: & kind_phys ! Working type use module_radiation_astronomy,only: & @@ -28,11 +27,11 @@ end subroutine GFS_rrtmgp_sw_pre_init !! \htmlinclude GFS_rrtmgp_sw_pre.html !! subroutine GFS_rrtmgp_sw_pre_run(me, nCol, nLev, lndp_type, n_var_lndp,lndp_var_list, & - lndp_prt_list, lsswr, solhr, & + lndp_prt_list, doSWrad, solhr, & lon, coslat, sinlat, snowd, sncovr, snoalb, zorl, tsfc, hprime, alvsf, & alnsf, alvwf, alnwf, facsf, facwf, fice, tisfc, lsmask, sfc_wts, p_lay, tv_lay, & relhum, p_lev, sw_gas_props, & - nday, idxday, alb1d, coszen, coszdg, sfc_alb_nir_dir, sfc_alb_nir_dif, & + nday, idxday, coszen, coszdg, sfc_alb_nir_dir, sfc_alb_nir_dif, & sfc_alb_uvvis_dir, sfc_alb_uvvis_dif, sfc_alb_dif, errmsg, errflg) ! Inputs @@ -47,7 +46,7 @@ subroutine GFS_rrtmgp_sw_pre_run(me, nCol, nLev, lndp_type, n_var_lndp,lndp_var_ real(kind_phys), dimension(n_var_lndp), intent(in) :: & lndp_prt_list logical,intent(in) :: & - lsswr ! Call RRTMGP SW radiation? + doSWrad ! Call RRTMGP SW radiation? real(kind_phys), intent(in) :: & solhr ! Time in hours after 00z at the current timestep real(kind_phys), dimension(nCol), intent(in) :: & @@ -86,7 +85,6 @@ subroutine GFS_rrtmgp_sw_pre_run(me, nCol, nLev, lndp_type, n_var_lndp,lndp_var_ integer, dimension(ncol), intent(out) :: & idxday ! Indices for daylit points real(kind_phys), dimension(ncol), intent(out) :: & - alb1d, & ! Surface albedo pertubation coszen, & ! Cosine of SZA coszdg, & ! Cosine of SZA, daytime sfc_alb_dif ! Mean surface diffused (nIR+uvvis) sw albedo @@ -103,65 +101,62 @@ subroutine GFS_rrtmgp_sw_pre_run(me, nCol, nLev, lndp_type, n_var_lndp,lndp_var_ ! Local variables integer :: i, j, iCol, iBand, iLay real(kind_phys), dimension(ncol, NF_ALBD) :: sfcalb + real(kind_phys), dimension(ncol) :: alb1d real(kind_phys) :: lndp_alb ! Initialize CCPP error handling variables errmsg = '' errflg = 0 - - if (.not. lsswr) return - - ! ####################################################################################### - ! Compute cosine of zenith angle (only when SW is called) - ! ####################################################################################### - call coszmn (lon, sinlat, coslat, solhr, nCol, me, coszen, coszdg) - ! ####################################################################################### - ! For SW gather daylit points - ! ####################################################################################### - nday = 0 - idxday = 0 - do i = 1, NCOL - if (coszen(i) >= 0.0001) then - nday = nday + 1 - idxday(nday) = i - endif - enddo + if (doSWrad) then - ! ####################################################################################### - ! mg, sfc-perts - ! --- scale random patterns for surface perturbations with perturbation size - ! --- turn vegetation fraction pattern into percentile pattern - ! ####################################################################################### - alb1d(:) = 0. - lndp_alb = -999. - if (lndp_type ==1) then - do k =1,n_var_lndp - if (lndp_var_list(k) == 'alb') then - do i=1,ncol - call cdfnor(sfc_wts(i,k),alb1d(i)) - lndp_alb = lndp_prt_list(k) - enddo - endif - enddo - endif - - ! ####################################################################################### - ! Call module_radiation_surface::setalb() to setup surface albedo. - ! ####################################################################################### - call setalb (lsmask, snowd, sncovr, snoalb, zorl, coszen, tsfc, tsfc, hprime, alvsf, & - alnsf, alvwf, alnwf, facsf, facwf, fice, tisfc, NCOL, alb1d, pertalb, sfcalb) + ! #################################################################################### + ! Compute cosine of zenith angle (only when SW is called) + ! #################################################################################### + call coszmn (lon, sinlat, coslat, solhr, nCol, me, coszen, coszdg) + + ! #################################################################################### + ! For SW gather daylit points + ! #################################################################################### + nday = 0 + idxday = 0 + do i = 1, NCOL + if (coszen(i) >= 0.0001) then + nday = nday + 1 + idxday(nday) = i + endif + enddo + + ! #################################################################################### + ! Call module_radiation_surface::setalb() to setup surface albedo. + ! #################################################################################### + alb1d(:) = 0. + lndp_alb = -999. + call setalb (lsmask, snowd, sncovr, snoalb, zorl, coszen, tsfc, tsfc, hprime, alvsf, & + alnsf, alvwf, alnwf, facsf, facwf, fice, tisfc, NCOL, alb1d, lndp_alb, sfcalb) - ! Approximate mean surface albedo from vis- and nir- diffuse values. - sfc_alb_dif(:) = max(0.01, 0.5 * (sfcalb(:,2) + sfcalb(:,4))) + ! Approximate mean surface albedo from vis- and nir- diffuse values. + sfc_alb_dif(:) = max(0.01, 0.5 * (sfcalb(:,2) + sfcalb(:,4))) - ! Spread across all SW bands - do iBand=1,sw_gas_props%get_nband() - sfc_alb_nir_dir(iBand,1:NCOL) = sfcalb(1:NCOL,1) - sfc_alb_nir_dif(iBand,1:NCOL) = sfcalb(1:NCOL,2) - sfc_alb_uvvis_dir(iBand,1:NCOL) = sfcalb(1:NCOL,3) - sfc_alb_uvvis_dif(iBand,1:NCOL) = sfcalb(1:NCOL,4) - enddo + ! Spread across all SW bands + do iBand=1,sw_gas_props%get_nband() + sfc_alb_nir_dir(iBand,1:NCOL) = sfcalb(1:NCOL,1) + sfc_alb_nir_dif(iBand,1:NCOL) = sfcalb(1:NCOL,2) + sfc_alb_uvvis_dir(iBand,1:NCOL) = sfcalb(1:NCOL,3) + sfc_alb_uvvis_dif(iBand,1:NCOL) = sfcalb(1:NCOL,4) + enddo + else + nday = 0 + idxday = 0 + coszen(1:nCol) = 0. + coszdg(1:nCol) = 0. + sfc_alb_nir_dir(:,1:nCol) = 0. + sfc_alb_nir_dif(:,1:nCol) = 0. + sfc_alb_uvvis_dir(:,1:nCol) = 0. + sfc_alb_uvvis_dif(:,1:nCol) = 0. + sfc_alb_dif(1:nCol) = 0. + endif + end subroutine GFS_rrtmgp_sw_pre_run diff --git a/physics/GFS_rrtmgp_sw_pre.meta b/physics/GFS_rrtmgp_sw_pre.meta index 91e875c00..b24ab5710 100644 --- a/physics/GFS_rrtmgp_sw_pre.meta +++ b/physics/GFS_rrtmgp_sw_pre.meta @@ -65,7 +65,7 @@ kind = len=3 intent = in optional = F -[lsswr] +[doSWrad] standard_name = flag_to_calc_sw long_name = logical flags for sw radiation calls units = flag @@ -297,15 +297,6 @@ type = ty_gas_optics_rrtmgp intent = in optional = F -[alb1d] - standard_name = surface_albedo_perturbation - long_name = surface albedo perturbation - units = frac - dimensions = (horizontal_loop_extent) - type = real - kind = kind_phys - intent = out - optional = F [sfc_alb_nir_dir] standard_name = surface_albedo_nearIR_direct long_name = near-IR (direct) surface albedo (sfc_alb_nir_dir) diff --git a/physics/GFS_time_vary_pre.scm.meta b/physics/GFS_time_vary_pre.scm.meta index 98d0d3b63..6241e29f1 100644 --- a/physics/GFS_time_vary_pre.scm.meta +++ b/physics/GFS_time_vary_pre.scm.meta @@ -212,7 +212,7 @@ [ipt] standard_name = index_for_diagnostic_printout long_name = horizontal index for point used for diagnostic printout - units = + units = index dimensions = () type = integer intent = out diff --git a/physics/physcons.F90 b/physics/physcons.F90 index 6a41bda44..397fee935 100644 --- a/physics/physcons.F90 +++ b/physics/physcons.F90 @@ -138,7 +138,7 @@ module physcons real(kind=kind_phys),parameter:: rhosnow = 100._kind_phys !< density of snow (kg/m^3) real(kind=kind_phys),parameter:: rhoair = 1.28_kind_phys !< density of air near surface (kg/m^3) -! Decorrelation length constant (km) for iovrlw/iovrsw = 4 or 5 and idcor = 0 +! Decorrelation length constant (km) for iovr = 4 or 5 and idcor = 0 real(kind=kind_phys),parameter:: decorr_con = 2.50_kind_phys !........................................! diff --git a/physics/physparam.f b/physics/physparam.f index c71b62e5b..5518c6163 100644 --- a/physics/physparam.f +++ b/physics/physparam.f @@ -229,25 +229,16 @@ module physparam !!\n =1:use prognostic cloud scheme for cloud cover and cloud properties integer, save :: icldflg = 1 -!> cloud overlapping control flag for SW +!> cloud overlapping control flag for Radiation !!\n =0:use random cloud overlapping method !!\n =1:use maximum-random cloud overlapping method !!\n =2:use maximum cloud overlapping method !!\n =3:use decorrelation length overlapping method !!\n =4:use exponential overlapping method !!\n =5:use exponential-random overlapping method -!!\n Opr GFS/CFS=1; see IOVR_SW in run scripts - integer, save :: iovrsw = 1 -!> cloud overlapping control flag for LW -!!\n =0:use random cloud overlapping method -!!\n =1:use maximum-random cloud overlapping method -!!\n =2:use maximum cloud overlapping method -!!\n =3:use decorrelation length overlapping method -!!\n =4:use exponential overlapping method -!!\n =5:use exponential-random overlapping method -!!\n Opr GFS/CFS=1; see IOVR_LW in run scripts - integer, save :: iovrlw = 1 -!!\n Decorrelation length type for iovrlw/iovrsw = 4 or 5 +!!\n Opr GFS/CFS=1; see IOVR in run scripts + integer, save :: iovr = 1 +!!\n Decorrelation length type for iovr = 4 or 5 !!\n =0:use constant decorrelation length defined by decorr_con (in module physcons) !!\n =1:use day-of-year and latitude-varying decorrelation length integer, save :: idcor = 1 diff --git a/physics/radiation_cloud_overlap.F90 b/physics/radiation_cloud_overlap.F90 new file mode 100644 index 000000000..a94923ba5 --- /dev/null +++ b/physics/radiation_cloud_overlap.F90 @@ -0,0 +1,116 @@ +module module_radiation_cloud_overlap + use physparam, only : kind_phys + implicit none + public :: cmp_dcorr_lgth + + interface cmp_dcorr_lgth + module procedure cmp_dcorr_lgth_hogan + module procedure cmp_dcorr_lgth_oreopoulos + end interface + +contains + ! ###################################################################################### + ! Hogan et al. (2010) + ! "Effect of improving representation of horizontal and vertical cloud structure on the + ! Earth's global radiation budget. Part I: Review and parametrization" + ! https://rmets.onlinelibrary.wiley.com/doi/full/10.1002/qj.647 + ! ###################################################################################### + subroutine cmp_dcorr_lgth_hogan(nCol, lat, con_pi, dcorr_lgth) + ! Inputs + integer, intent(in) :: & + nCol ! Number of horizontal grid-points + real(kind_phys), intent(in) :: & + con_pi ! Physical constant: Pi + real(kind_phys), dimension(nCol), intent(in) :: & + lat ! Latitude + ! Outputs + real(kind_phys), dimension(nCol),intent(out) :: & + dcorr_lgth ! Decorrelation length + + ! Local variables + integer :: iCol + + ! Parameters + real(kind_phys),parameter :: min_dcorr = 0.6 ! (see section 2.3) + + do iCol =1,nCol + dcorr_lgth(iCol) = max(min_dcorr, 2.78-4.6*abs(lat(iCol)/con_pi)) ! (eq. 13) + enddo + + end subroutine cmp_dcorr_lgth_hogan + ! ###################################################################################### + ! Oreopoulos et al. (2012) + ! "Radiative impacts of cloud heterogeneity and overlap in an + ! atmospheric General Circulation Model" + ! 10.5194/acp-12-9097-2012 + ! ###################################################################################### + subroutine cmp_dcorr_lgth_oreopoulos(nCol, lat, juldat, yearlength, dcorr_lgth) + ! Inputs + integer, intent(in) :: & + nCol, & ! Number of horizontal grid-points + yearlength ! Number of days in year + + real(kind_phys), intent(in) :: & + juldat ! Julian date + real(kind_phys), dimension(nCol), intent(in) :: & + lat ! Latitude + + ! Outputs + real(kind_phys), dimension(nCol),intent(out) :: & + dcorr_lgth ! Decorrelation length (km) + + ! Parameters for the Gaussian fits per Eqs. (10) and (11) (See Table 1) + real(kind_phys), parameter :: & ! + am1 = 1.4315_kind_phys, & ! + am2 = 2.1219_kind_phys, & ! + am4 = -25.584_kind_phys, & ! + amr = 7.0_kind_phys ! + + ! Local variables + integer :: iCol + real(kind_phys) :: am3 + + do iCol = 1, nCol + if (juldat .gt. 181._kind_phys) then + am3 = -4._kind_phys * amr * (juldat - 272._kind_phys) / yearlength ! (eq. 11a) + else + am3 = 4._kind_phys * amr * (juldat - 91._kind_phys) / yearlength ! (eq. 11b) + endif + dcorr_lgth(iCol) = am1 + am2 * exp( -(lat(iCol) - am3)**2 / am4**2) ! (eq. 10) + enddo + + end subroutine cmp_dcorr_lgth_oreopoulos + + ! ###################################################################################### + ! + ! ###################################################################################### + subroutine get_alpha_exp(nCol, nLay, dzlay, dcorr_lgth, alpha) + + ! Inputs + integer, intent(in) :: & + nCol, & ! Number of horizontal grid points + nLay ! Number of vertical grid points + real(kind_phys), dimension(nCol), intent(in) :: & + dcorr_lgth ! Decorrelation length (km) + real(kind_phys), dimension(nCol,nLay), intent(in) :: & + dzlay ! + + ! Outputs + real(kind_phys), dimension(nCol,nLay) :: & + alpha ! Cloud overlap parameter + + ! Local variables + integer :: iCol,iLay + + do iCol = 1, nCol + alpha(iCol,1) = 0.0d0 + do iLay = 2, nLay + alpha(iCol,iLay) = exp( -(dzlay(iCol,iLay)) / dcorr_lgth(iCol)) + enddo + enddo + + return + + end subroutine get_alpha_exp + +end module module_radiation_cloud_overlap diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 552037da2..056bede28 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -188,7 +188,7 @@ !!\n IMP_PHYSICS =98/99: Zhao-Carr-Sundqvist MP - Xu-Randall diagnostic cloud fraction !!\n IMP_PHYSICS =11: GFDL MP - unified diagnostic cloud fraction provided by GFDL MP !! -!! Cloud overlapping method (namelist control parameter - \b IOVR_LW, \b IOVR_SW) +!! Cloud overlapping method (namelist control parameter - \b IOVR) !!\n IOVR=0: randomly overlapping vertical cloud layers !!\n IOVR=1: maximum-random overlapping vertical cloud layers !!\n IOVR=2: maximum overlapping vertical cloud layers @@ -208,14 +208,16 @@ !> This module computes cloud related quantities for radiation computations. module module_radiation_clouds ! - use physparam, only : icldflg, iovrsw, iovrlw, & + use physparam, only : icldflg, iovr, idcor, & & lcrick, lcnorm, lnoprec, & & ivflip use physcons, only : con_fvirt, con_ttp, con_rocp, & & con_t0c, con_pi, con_g, con_rd, & - & con_thgni + & con_thgni, decorr_con use module_microphysics, only : rsipath2 use module_iounitdef, only : NICLTUN + use module_radiation_cloud_overlap, only: cmp_dcorr_lgth, & + & get_alpha_exp use machine, only : kind_phys ! implicit none @@ -255,13 +257,11 @@ module module_radiation_clouds real (kind=kind_phys), parameter :: cldasy_def = 0.84 !< default cld asymmetry factor integer :: llyr = 2 !< upper limit of boundary layer clouds - integer :: iovr = 1 !< maximum-random cloud overlapping method public progcld1, progcld2, progcld3, progcld4, progclduni, & & cld_init, progcld5, progcld6, progcld4o, cal_cldfra3, & & find_cloudLayers, adjust_cloudIce, adjust_cloudH2O, & - & adjust_cloudFinal, gethml, get_alpha_dcorr, get_alpha_exp - + & adjust_cloudFinal, gethml ! ================= contains @@ -314,7 +314,7 @@ subroutine cld_init & ! =8: Thompson microphysics ! ! =6: WSM6 microphysics ! ! =10: MG microphysics ! -! iovrsw/iovrlw : sw/lw control flag for cloud overlapping scheme ! +! iovr : control flag for cloud overlapping scheme ! ! =0: random overlapping clouds ! ! =1: max/ran overlapping clouds ! ! =2: maximum overlap clouds (mcica only) ! @@ -347,8 +347,6 @@ subroutine cld_init & ! ! --- set up module variables - iovr = max( iovrsw, iovrlw ) !cld ovlp used for diag HML cld output - if (me == 0) print *, VTAGCLD !print out version tag if ( icldflg == 0 ) then @@ -839,14 +837,23 @@ subroutine progcld1 & enddo endif -!> - Call subroutine get_alpha_exp to define alpha parameter for EXP and ER cloud overlap options - if ( iovr == 4 .or. iovr == 5 ) then - call get_alpha_exp & -! --- inputs: - & (ix, nlay, dzlay, iovr, latdeg, julian, yearlen, cldtot, & -! --- outputs: - & alpha & - & ) + ! Compute cloud decorrelation length + if (idcor == 1) then + call cmp_dcorr_lgth(ix, xlat, con_pi, de_lgth) + endif + if (idcor == 2) then + call cmp_dcorr_lgth(ix, latdeg, julian, yearlen, de_lgth) + endif + if (idcor == 0) then + de_lgth(:) = decorr_con + endif + + ! Call subroutine get_alpha_exp to define alpha parameter for exponential cloud overlap options + if (iovr == 3 .or. iovr == 4 .or. iovr == 5) then + call get_alpha_exp(ix, nLay, dzlay, de_lgth, alpha) + else + de_lgth(:) = 0. + alpha(:,:) = 0. endif !> - Call gethml() to compute low,mid,high,total, and boundary layer @@ -1226,23 +1233,23 @@ subroutine progcld2 & enddo enddo -! --- ... estimate clouds decorrelation length in km -! this is only a tentative test, need to consider change later - - if ( iovr == 3 ) then - do i = 1, ix - de_lgth(i) = max( 0.6, 2.78-4.6*rxlat(i) ) - enddo + ! Compute cloud decorrelation length + if (idcor == 1) then + call cmp_dcorr_lgth(ix, xlat, con_pi, de_lgth) + endif + if (idcor == 2) then + call cmp_dcorr_lgth(ix, latdeg, julian, yearlen, de_lgth) + endif + if (idcor == 0) then + de_lgth(:) = decorr_con endif -!> - Call subroutine get_alpha_exp to define alpha parameter for EXP and ER cloud overlap options - if ( iovr == 4 .or. iovr == 5 ) then - call get_alpha_exp & -! --- inputs: - & (ix, nlay, dzlay, iovr, latdeg, julian, yearlen, cldtot, & -! --- outputs: - & alpha & - & ) + ! Call subroutine get_alpha_exp to define alpha parameter for exponential cloud overlap options + if (iovr == 3 .or. iovr == 4 .or. iovr == 5) then + call get_alpha_exp(ix, nLay, dzlay, de_lgth, alpha) + else + de_lgth(:) = 0. + alpha(:,:) = 0. endif !> - Call gethml() to compute low,mid,high,total, and boundary layer @@ -1653,23 +1660,23 @@ subroutine progcld3 & enddo enddo -! --- ... estimate clouds decorrelation length in km -! this is only a tentative test, need to consider change later - - if ( iovr == 3 ) then - do i = 1, ix - de_lgth(i) = max( 0.6, 2.78-4.6*rxlat(i) ) - enddo + ! Compute cloud decorrelation length + if (idcor == 1) then + call cmp_dcorr_lgth(ix, xlat, con_pi, de_lgth) + endif + if (idcor == 2) then + call cmp_dcorr_lgth(ix, latdeg, julian, yearlen, de_lgth) + endif + if (idcor == 0) then + de_lgth(:) = decorr_con endif -!> - Call subroutine get_alpha_exp to define alpha parameter for EXP and ER cloud overlap options - if ( iovr == 4 .or. iovr == 5 ) then - call get_alpha_exp & -! --- inputs: - & (ix, nlay, dzlay, iovr, latdeg, julian, yearlen, cldtot, & -! --- outputs: - & alpha & - & ) + ! Call subroutine get_alpha_exp to define alpha parameter for exponential cloud overlap options + if (iovr == 3 .or. iovr == 4 .or. iovr == 5) then + call get_alpha_exp(ix, nLay, dzlay, de_lgth, alpha) + else + de_lgth(:) = 0. + alpha(:,:) = 0. endif !> -# Call gethml() to compute low,mid,high,total, and boundary layer @@ -2016,23 +2023,23 @@ subroutine progcld4 & enddo enddo -! --- ... estimate clouds decorrelation length in km -! this is only a tentative test, need to consider change later - - if ( iovr == 3 ) then - do i = 1, ix - de_lgth(i) = max( 0.6, 2.78-4.6*rxlat(i) ) - enddo + ! Compute cloud decorrelation length + if (idcor == 1) then + call cmp_dcorr_lgth(ix, xlat, con_pi, de_lgth) + endif + if (idcor == 2) then + call cmp_dcorr_lgth(ix, latdeg, julian, yearlen, de_lgth) + endif + if (idcor == 0) then + de_lgth(:) = decorr_con endif -!> - Call subroutine get_alpha_exp to define alpha parameter for EXP and ER cloud overlap options - if ( iovr == 4 .or. iovr == 5 ) then - call get_alpha_exp & -! --- inputs: - & (ix, nlay, dzlay, iovr, latdeg, julian, yearlen, cldtot, & -! --- outputs: - & alpha & - & ) + ! Call subroutine get_alpha_exp to define alpha parameter for exponential cloud overlap options + if (iovr == 3 .or. iovr == 4 .or. iovr == 5) then + call get_alpha_exp(ix, nLay, dzlay, de_lgth, alpha) + else + de_lgth(:) = 0. + alpha(:,:) = 0. endif ! --- compute low, mid, high, total, and boundary layer cloud fractions @@ -2370,23 +2377,23 @@ subroutine progcld4o & enddo enddo -! --- ... estimate clouds decorrelation length in km -! this is only a tentative test, need to consider change later - - if ( iovr == 3 ) then - do i = 1, ix - de_lgth(i) = max( 0.6, 2.78-4.6*rxlat(i) ) - enddo + ! Compute cloud decorrelation length + if (idcor == 1) then + call cmp_dcorr_lgth(ix, xlat, con_pi, de_lgth) + endif + if (idcor == 2) then + call cmp_dcorr_lgth(ix, latdeg, julian, yearlen, de_lgth) + endif + if (idcor == 0) then + de_lgth(:) = decorr_con endif -!> - Call subroutine get_alpha_exp to define alpha parameter for EXP and ER cloud overlap options - if ( iovr == 4 .or. iovr == 5 ) then - call get_alpha_exp & -! --- inputs: - & (ix, nlay, dzlay, iovr, latdeg, julian, yearlen, cldtot, & -! --- outputs: - & alpha & - & ) + ! Call subroutine get_alpha_exp to define alpha parameter for exponential cloud overlap options + if (iovr == 3 .or. iovr == 4 .or. iovr == 5) then + call get_alpha_exp(ix, nLay, dzlay, de_lgth, alpha) + else + de_lgth(:) = 0. + alpha(:,:) = 0. endif !> - Call gethml() to compute low, mid, high, total, and boundary layer cloud fractions @@ -2803,23 +2810,23 @@ subroutine progcld5 & enddo enddo -! --- ... estimate clouds decorrelation length in km -! this is only a tentative test, need to consider change later - - if ( iovr == 3 ) then - do i = 1, ix - de_lgth(i) = max( 0.6, 2.78-4.6*rxlat(i) ) - enddo + ! Compute cloud decorrelation length + if (idcor == 1) then + call cmp_dcorr_lgth(ix, xlat, con_pi, de_lgth) + endif + if (idcor == 2) then + call cmp_dcorr_lgth(ix, latdeg, julian, yearlen, de_lgth) + endif + if (idcor == 0) then + de_lgth(:) = decorr_con endif -!> - Call subroutine get_alpha_exp to define alpha parameter for EXP and ER cloud overlap options - if ( iovr == 4 .or. iovr == 5 ) then - call get_alpha_exp & -! --- inputs: - & (ix, nlay, dzlay, iovr, latdeg, julian, yearlen, cldtot, & -! --- outputs: - & alpha & - & ) + ! Call subroutine get_alpha_exp to define alpha parameter for exponential cloud overlap options + if (iovr == 3 .or. iovr == 4 .or. iovr == 5) then + call get_alpha_exp(ix, nLay, dzlay, de_lgth, alpha) + else + de_lgth(:) = 0. + alpha(:,:) = 0. endif !> - Call gethml() to compute low,mid,high,total, and boundary layer @@ -3170,23 +3177,23 @@ subroutine progcld6 & enddo enddo -! --- ... estimate clouds decorrelation length in km -! this is only a tentative test, need to consider change later - - if ( iovr == 3 ) then - do i = 1, ix - de_lgth(i) = max( 0.6, 2.78-4.6*rxlat(i) ) - enddo + ! Compute cloud decorrelation length + if (idcor == 1) then + call cmp_dcorr_lgth(ix, xlat, con_pi, de_lgth) + endif + if (idcor == 2) then + call cmp_dcorr_lgth(ix, latdeg, julian, yearlen, de_lgth) + endif + if (idcor == 0) then + de_lgth(:) = decorr_con endif -!> - Call subroutine get_alpha_exp to define alpha parameter for EXP and ER cloud overlap options - if ( iovr == 4 .or. iovr == 5 ) then - call get_alpha_exp & -! --- inputs: - & (ix, nlay, dzlay, iovr, latdeg, julian, yearlen, cldtot, & -! --- outputs: - & alpha & - & ) + ! Call subroutine get_alpha_exp to define alpha parameter for exponential cloud overlap options + if ( iovr == 4 .or. iovr == 5) then + call get_alpha_exp(ix, nLay, dzlay, de_lgth, alpha) + else + de_lgth(:) = 0. + alpha(:,:) = 0. endif !> - Call gethml() to compute low,mid,high,total, and boundary layer @@ -3567,23 +3574,23 @@ subroutine progclduni & enddo enddo -!> -# Estimate clouds decorrelation length in km -! this is only a tentative test, need to consider change later - - if ( iovr == 3 ) then - do i = 1, ix - de_lgth(i) = max( 0.6, 2.78-4.6*rxlat(i) ) - enddo + ! Compute cloud decorrelation length + if (idcor == 1) then + call cmp_dcorr_lgth(ix, xlat, con_pi, de_lgth) + endif + if (idcor == 2) then + call cmp_dcorr_lgth(ix, latdeg, julian, yearlen, de_lgth) + endif + if (idcor == 0) then + de_lgth(:) = decorr_con endif -!> - Call subroutine get_alpha_exp to define alpha parameter for EXP and ER cloud overlap options - if ( iovr == 4 .or. iovr == 5 ) then - call get_alpha_exp & -! --- inputs: - & (ix, nlay, dzlay, iovr, latdeg, julian, yearlen, cldtot, & -! --- outputs: - & alpha & - & ) + ! Call subroutine get_alpha_exp to define alpha parameter for exponential cloud overlap options + if (iovr == 3 .or. iovr == 4 .or. iovr == 5) then + call get_alpha_exp(ix, nLay, dzlay, de_lgth, alpha) + else + de_lgth(:) = 0. + alpha(:,:) = 0. endif !> - Call gethml() to compute low,mid,high,total, and boundary layer @@ -4043,219 +4050,6 @@ subroutine gethml & end subroutine gethml !----------------------------------- !! @} - ! ######################################################################################### - ! Subroutine to compute cloud-overlap parameter, alpha, for decorrelation-length cloud - ! overlap assumption. - ! ######################################################################################### - subroutine get_alpha_dcorr(nCol, nLev, lat, con_pi, deltaZ, & - & de_lgth, cloud_overlap_param) - - integer, intent(in) :: nCol, nLev - real(kind_phys), intent(in) :: con_pi - real(kind_phys), dimension(nCol), intent(in) :: lat - real(kind_phys), dimension(nCol,nLev),intent(in) :: deltaZ - real(kind_phys), dimension(nCol),intent(out) :: de_lgth - real(kind_phys), dimension(nCol,nLev),intent(out) :: & - & cloud_overlap_param - - ! Local - integer :: iCol, iLay - - do iCol =1,nCol - de_lgth(iCol) = max( 0.6, 2.78-4.6*abs(lat(iCol)/con_pi) ) - do iLay=nLev,2,-1 - if (de_lgth(iCol) .gt. 0) then - cloud_overlap_param(iCol,iLay-1) = & - & exp( -0.5 * (deltaZ(iCol,iLay)+deltaZ(iCol,iLay-1))/& - & de_lgth(iCol)) - endif - enddo - enddo - end subroutine get_alpha_dcorr - - ! ######################################################################################### -!> \ingroup module_radiation_clouds -!! This program derives the exponential transition, alpha, from maximum to -!! random overlap needed to define the fractional cloud vertical correlation -!! for the exponential (EXP, iovrlp=4) or the exponential-random (ER, iovrlp=5) -!! cloud overlap options for RRTMG/RRTMGP. For exponential, the transition from -!! maximum to random with distance through model layers occurs without regard -!! to the configuration of clear and cloudy layers. For the ER method, each -!! block of adjacent cloudy layers is treated with a separate transition from -!! maximum to random, and blocks of cloudy layers separated by one or more -!! clear layers are correlated randomly. -!> /param nlon : number of model longitude points -!> /param nlay : vertical layer dimension -!> /param dzlay(nlon,nlay) : distance between the center of model layers -!> /param iovrlp : cloud overlap method -!> : 0 = random -!> : 1 = maximum-random -!> : 2 = maximum -!> : 3 = decorrelation (NOAA/Hou) -!> : 4 = exponential (AER) -!> : 5 = exponential-random (AER) -!> /param latdeg(nlon) : latitude (in degrees 90 -> -90) -!> /param juldat : day of the year (fractional julian day) -!> /param yearlen : current length of the year (365/366 days) -!> /param cldf(nlon,nlay) : cloud fraction -!> /param idcor : decorrelation length method -!> : 0 = constant value (AER; decorr_con) -!> : 1 = latitude and day of year varying value (AER; Oreopoulos, et al., 2012) -!> /param decorr_con : decorrelation length constant -!! -!>\section detail Detailed Algorithm -!! @{ - subroutine get_alpha_exp & -! --- inputs: - & (nlon, nlay, dzlay, iovrlp, latdeg, juldat, yearlen, cldf, & -! --- outputs: - & alpha & - & ) - -! =================================================================== ! -! ! -! abstract: Derives the exponential transition, alpha, from maximum to ! -! random overlap needed to define the fractional cloud vertical ! -! correlation for the exponential (EXP, iovrlp=4) or the exponential- ! -! random (ER, iovrlp=5) cloud overlap options for RRTMG. For ! -! exponential, the transition from maximum to random with distance ! -! through model layers occurs without regard to the configuration of ! -! clear and cloudy layers. For the ER method, each block of adjacent ! -! cloudy layers is treated with a separate transition from maximum to ! -! random, and blocks of cloudy layers separated by one or more ! -! clear layers are correlated randomly. ! -! ! -! usage: call get_alpha_exp ! -! ! -! subprograms called: none ! -! ! -! attributes: ! -! language: fortran 90 ! -! machine: ibm-sp, sgi ! -! ! -! author: m.j. iacono (AER) for use with the RRTMG radiation code ! -! ! -! ==================== definition of variables ==================== ! -! ! -! Input variables: ! -! nlon : number of model longitude points ! -! nlay : vertical layer dimension ! -! dzlay(nlon,nlay) : distance between the center of model layers ! -! iovrlp : cloud overlap method ! -! : 0 = random ! -! : 1 = maximum-random ! -! : 2 = maximum ! -! : 3 = decorrelation (NOAA/Hou) ! -! : 4 = exponential (AER) ! -! : 5 = exponential-random (AER) ! -! latdeg(nlon) : latitude (in degrees 90 -> -90) ! -! juldat : day of the year (fractional julian day) ! -! yearlen : current length of the year (365/366 days) ! -! cldf(nlon,nlay) : cloud fraction ! -! ! -! output variables: ! -! alpha(nlon,nlay) : alpha exponential transition parameter for ! -! : cloud vertical correlation ! -! ! -! external module variables: (in physcons) ! -! decorr_con : decorrelation length constant (km) ! -! ! -! external module variables: (in physparam) ! -! idcor : control flag for decorrelation length method ! -! =0: constant decorrelation length (decorr_con) ! -! =1: latitude and day-of-year varying decorrelation! -! length (AER; Oreopoulos, et al., 2012) ! -! ! -! ==================== end of description ===================== ! -! - use physcons, only: decorr_con - use physparam, only: idcor - - implicit none - -! Input - integer, intent(in) :: nlon, nlay - integer, intent(in) :: iovrlp - integer, intent(in) :: yearlen - real(kind=kind_phys), dimension(:,:), intent(in) :: dzlay - real(kind=kind_phys), dimension(:,:), intent(in) :: cldf - real(kind=kind_phys), dimension(:), intent(in) :: latdeg - real(kind=kind_phys), intent(in) :: juldat - -! Output - real(kind=kind_phys), dimension(:,:), intent(out):: alpha - -! Local - integer :: i, k - real(kind=kind_phys) :: decorr_len(nlon) ! Decorrelation length (km) - -! Constants for latitude and day-of-year dependent decorrlation length (Oreopoulos et al, 2012) -! Used when idcor = 1 - real(kind=kind_phys), parameter :: am1 = 1.4315_kind_phys - real(kind=kind_phys), parameter :: am2 = 2.1219_kind_phys - real(kind=kind_phys), parameter :: am4 = -25.584_kind_phys - real(kind=kind_phys), parameter :: amr = 7.0_kind_phys - real(kind=kind_phys) :: am3 - - real(kind=kind_phys), parameter :: zero = 0.0d0 - real(kind=kind_phys), parameter :: one = 1.0d0 - -! -!===> ... begin here -! -! If exponential or exponential-random cloud overlap is used: -! derive day-of-year and latitude-varying decorrelation lendth if requested; -! otherwise use the constant decorrelation length, decorr_con, specified in physcons.F90 - do i = 1, nlon - if (iovrlp == 4 .or. iovrlp == 5) then - if (idcor .eq. 1) then - if (juldat .gt. 181._kind_phys) then - am3 = -4._kind_phys * amr * (juldat - 272._kind_phys) - & / yearlen - else - am3 = 4._kind_phys * amr * (juldat - 91._kind_phys) - & / yearlen - endif -! For latitude in degrees, decorr_len in km - decorr_len(i) = am1 + am2 * exp( -(latdeg(i) - am3)**2 - & / am4**2) - else - decorr_len(i) = decorr_con - endif - endif - enddo - -! For atmospheric data defined from surface to toa; define alpha from surface to toa -! Exponential cloud overlap - if (iovrlp == 4) then - do i = 1, nlon - alpha(i,1) = zero - do k = 2, nlay - alpha(i,k) = exp( -(dzlay(i,k)) / decorr_len(i)) - enddo - enddo - endif -! Exponential-random cloud overlap - if (iovrlp == 5) then - do i = 1, nlon - alpha(i,1) = zero - do k = 2, nlay - alpha(i,k) = exp( -(dzlay(i,k)) / decorr_len(i)) - ! Decorrelate layers when a clear layer follows a cloudy layer to enforce - ! random correlation between non-adjacent blocks of cloudy layers - if (cldf(i,k) .eq. zero .and. cldf(i,k-1) .gt. zero) then - alpha(i,k) = zero - endif - enddo - enddo - endif - - return - - end subroutine get_alpha_exp -!----------------------------------- -!! @} - !+---+-----------------------------------------------------------------+ !..Cloud fraction scheme by G. Thompson (NCAR-RAL), not intended for !.. combining with any cumulus or shallow cumulus parameterization diff --git a/physics/radlw_main.F90 b/physics/radlw_main.F90 index daa20e45d..de8d9e973 100644 --- a/physics/radlw_main.F90 +++ b/physics/radlw_main.F90 @@ -279,8 +279,7 @@ module rrtmg_lw ! use physparam, only : ilwrate, ilwrgas, ilwcliq, ilwcice, & - & isubclw, icldflg, iovrlw, ivflip, & - & kind_phys + & isubclw, icldflg, iovr, ivflip use physcons, only : con_g, con_cp, con_avgd, con_amd, & & con_amw, con_amo3 use mersenne_twister, only : random_setseed, random_number, & @@ -527,7 +526,7 @@ subroutine rrtmg_lw_run & ! =0: no sub-col cld treatment, use grid-mean cld quantities ! ! =1: mcica sub-col, prescribed seeds to get random numbers ! ! =2: mcica sub-col, providing array icseed for random numbers! -! iovrlw - cloud overlapping control flag ! +! iovr - cloud overlapping control flag ! ! =0: random overlapping clouds ! ! =1: maximum/random overlapping clouds ! ! =2: maximum overlap cloud (used for isubclw>0 only) ! @@ -849,10 +848,10 @@ subroutine rrtmg_lw_run & endif stemp = sfgtmp(iplon) ! surface ground temp - if (iovrlw == 3) delgth= de_lgth(iplon) ! clouds decorr-length + if (iovr == 3) delgth= de_lgth(iplon) ! clouds decorr-length ! mz*: HWRF - if (iovrlw == 4 ) then + if (iovr == 4 ) then !Add layer height needed for exponential (icld=4) and ! exponential-random (icld=5) overlap options @@ -876,7 +875,7 @@ subroutine rrtmg_lw_run & enddo enddo - call mcica_subcol_lw(1, iplon, nlay, iovrlw, permuteseed, & + call mcica_subcol_lw(1, iplon, nlay, iovr, permuteseed, & & irng, plyr, hgt, & & cld_cf, cld_iwp, cld_lwp,cld_swp, & & cld_ref_ice, cld_ref_liq, & @@ -910,7 +909,7 @@ subroutine rrtmg_lw_run & tavel(k)= tlyr(iplon,k1) tz(k) = tlvl(iplon,k1) dz(k) = dzlyr(iplon,k1) - if (iovrlw == 4 .or. iovrlw == 5) alph(k) = alpha(iplon,k) ! alpha decorrelation + if (iovr == 4 .or. iovr == 5) alph(k) = alpha(iplon,k) ! alpha decorrelation !> -# Set absorber amount for h2o, co2, and o3. @@ -989,7 +988,7 @@ subroutine rrtmg_lw_run & cda4(k) = cld_ref_snow(iplon,k1) enddo ! HWRF RRMTG - if (iovrlw == 4) then !mz HWRF + if (iovr == 4) then !mz HWRF do k = 1, nlay k1 = nlp1 - k do ig = 1, ngptlw @@ -1040,7 +1039,7 @@ subroutine rrtmg_lw_run & tavel(k)= tlyr(iplon,k) tz(k) = tlvl(iplon,k+1) dz(k) = dzlyr(iplon,k) - if (iovrlw == 4 .or. iovrlw == 5) alph(k) = alpha(iplon,k) ! alpha decorrelation + if (iovr == 4 .or. iovr == 5) alph(k) = alpha(iplon,k) ! alpha decorrelation ! --- ... set absorber amount !test use @@ -1112,7 +1111,7 @@ subroutine rrtmg_lw_run & cda3(k) = cld_swp(iplon,k) cda4(k) = cld_ref_snow(iplon,k) enddo - if (iovrlw == 4) then + if (iovr == 4) then !mz* Move incoming GCM cloud arrays to RRTMG cloud arrays. !For GCM input, incoming reicmcl is defined based on selected !ice parameterization (inflglw) @@ -1206,7 +1205,7 @@ subroutine rrtmg_lw_run & if ( lcf1 ) then !mz* for HWRF, save cldfmc with mcica - if (iovrlw == 4) then + if (iovr == 4) then do k = 1, nlay do ig = 1, ngptlw cldfmc_save(ig,k)=cldfmc (ig,k) @@ -1217,12 +1216,12 @@ subroutine rrtmg_lw_run & call cldprop & ! --- inputs: & ( cldfrc,clwp,relw,ciwp,reiw,cda1,cda2,cda3,cda4, & - & nlay, nlp1, ipseed(iplon), dz, delgth, iovrlw, alph, & + & nlay, nlp1, ipseed(iplon), dz, delgth, iovr, alph, & ! --- outputs: & cldfmc, taucld & & ) - if (iovrlw == 4) then + if (iovr == 4) then !mz for HWRF, still using mcica cldfmc do k = 1, nlay do ig = 1, ngptlw @@ -1251,7 +1250,7 @@ subroutine rrtmg_lw_run & endif !mz* HWRF: calculate taucmc with mcica - if (iovrlw == 4) then + if (iovr == 4) then call cldprmc(nlay, inflglw, iceflglw, liqflglw, & & cldfmc, ciwpmc, & & clwpmc, cswpmc, reicmc, relqmc, resnmc, & @@ -1344,7 +1343,7 @@ subroutine rrtmg_lw_run & if (isubclw <= 0) then - if (iovrlw <= 0) then + if (iovr <= 0) then call rtrn & ! --- inputs: @@ -1364,7 +1363,7 @@ subroutine rrtmg_lw_run & & totuflux,totdflux,htr, totuclfl,totdclfl,htrcl, htrb & & ) - endif ! end if_iovrlw_block + endif ! end if_iovr_block else @@ -1515,7 +1514,7 @@ subroutine rlwinit & ! icldflg - cloud scheme control flag ! ! =0: diagnostic scheme gives cloud tau, omiga, and g. ! ! =1: prognostic scheme gives cloud liq/ice path, etc. ! -! iovrlw - clouds vertical overlapping control flag ! +! iovr - clouds vertical overlapping control flag ! ! =0: random overlapping clouds ! ! =1: maximum/random overlapping clouds ! ! =2: maximum overlap cloud (isubcol>0 only) ! @@ -1564,19 +1563,19 @@ subroutine rlwinit & ! !===> ... begin here ! - if ( iovrlw<0 .or. iovrlw>4 ) then + if ( iovr<0 .or. iovr>4 ) then print *,' *** Error in specification of cloud overlap flag', & - & ' IOVRLW=',iovrlw,' in RLWINIT !!' + & ' IOVR=',iovr,' in RLWINIT !!' stop - elseif ( (iovrlw==2 .or. iovrlw==3) .and. isubclw==0 ) then + elseif ( (iovr==2 .or. iovr==3) .and. isubclw==0 ) then if (me == 0) then - print *,' *** IOVRLW=',iovrlw,' is not available for', & + print *,' *** IOVR=',iovr,' is not available for', & & ' ISUBCLW=0 setting!!' print *,' The program uses maximum/random overlap', & & ' instead.' endif - iovrlw = 1 + iovr = 1 endif if (me == 0) then @@ -1713,7 +1712,7 @@ end subroutine rlwinit !> @{ subroutine cldprop & & ( cfrac,cliqp,reliq,cicep,reice,cdat1,cdat2,cdat3,cdat4, & ! --- inputs - & nlay, nlp1, ipseed, dz, de_lgth, iovrlw, alpha, & + & nlay, nlp1, ipseed, dz, de_lgth, iovr, alpha, & & cldfmc, taucld & ! --- outputs & ) @@ -1814,7 +1813,7 @@ subroutine cldprop & use module_radlw_cldprlw ! --- inputs: - integer, intent(in) :: nlay, nlp1, ipseed, iovrlw + integer, intent(in) :: nlay, nlp1, ipseed, iovr real (kind=kind_phys), dimension(0:nlp1), intent(in) :: cfrac real (kind=kind_phys), dimension(nlay), intent(in) :: cliqp, & @@ -1994,7 +1993,7 @@ subroutine cldprop & ! --- ... call sub-column cloud generator !mz* - if (iovrlw .ne. 4) then + if (iovr .ne. 4) then call mcica_subcol & ! --- inputs: & ( cldf, nlay, ipseed, dz, de_lgth, alpha, & @@ -2011,7 +2010,7 @@ subroutine cldprop & endif enddo enddo - endif !iovrlw + endif !iovr endif ! end if_isubclw_block @@ -2054,7 +2053,7 @@ subroutine mcica_subcol & ! lcloudy - logical, sub-colum cloud profile flag array ngptlw*nlay! ! ! ! other control flags from module variables: ! -! iovrlw : control flag for cloud overlapping method ! +! iovr : control flag for cloud overlapping method ! ! =0:random; =1:maximum/random: =2:maximum; =3:decorr ! ! ! ! ===================== end of definitions ==================== ! @@ -2096,7 +2095,7 @@ subroutine mcica_subcol & !! - For max-random overlap, pick a random value at every level !! - For maximum overlap, pick same random numebr at every level - select case ( iovrlw ) + select case ( iovr ) case( 0 ) ! random overlap, pick a random value at every level diff --git a/physics/radsw_main.F90 b/physics/radsw_main.F90 index cf9e0e524..8ebbb3ab1 100644 --- a/physics/radsw_main.F90 +++ b/physics/radsw_main.F90 @@ -305,7 +305,7 @@ module rrtmg_sw ! use physparam, only : iswrate, iswrgas, iswcliq, iswcice, & - & isubcsw, icldflg, iovrsw, ivflip, & + & isubcsw, icldflg, iovr, ivflip, & & iswmode use physcons, only : con_g, con_cp, con_avgd, con_amd, & & con_amw, con_amo3 @@ -630,7 +630,7 @@ subroutine rrtmg_sw_run & ! =0: no sub-col cld treatment, use grid-mean cld quantities ! ! =1: mcica sub-col, prescribed seeds to get random numbers ! ! =2: mcica sub-col, providing array icseed for random numbers! -! iovrsw - cloud overlapping control flag ! +! iovr - cloud overlapping control flag ! ! =0: random overlapping clouds ! ! =1: maximum/random overlapping clouds ! ! =2: maximum overlap cloud ! @@ -972,7 +972,7 @@ subroutine rrtmg_sw_run & cosz1 = cosz(j1) sntz1 = f_one / cosz(j1) ssolar = s0fac * cosz(j1) - if (iovrsw == 3) delgth = de_lgth(j1) ! clouds decorr-length + if (iovr == 3) delgth = de_lgth(j1) ! clouds decorr-length !> -# Prepare surface albedo: bm,df - dir,dif; 1,2 - nir,uvv. albbm(1) = sfcalb_nir_dir(j1) @@ -982,7 +982,7 @@ subroutine rrtmg_sw_run & ! mz*: HWRF - if (iovrsw == 4 ) then + if (iovr == 4 ) then !Add layer height needed for exponential (icld=4) and @@ -993,7 +993,7 @@ subroutine rrtmg_sw_run & permuteseed = 1 !mz* Derive height of each layer mid-point from layer thickness. -! Needed for exponential (iovrsw=4) and exponential-random overlap +! Needed for exponential (iovr=4) and exponential-random overlap ! option (iovr=5)only. dzsum =0.0 do k = 1,nlay @@ -1012,7 +1012,7 @@ subroutine rrtmg_sw_run & enddo enddo - call mcica_subcol_sw (1, 1, nlay, iovrsw, permuteseed, & + call mcica_subcol_sw (1, 1, nlay, iovr, permuteseed, & & irng, plyr(j1:j1,:), hgt(j1:j1,:), & & cld_cf(j1:j1,:), cld_iwp(j1:j1,:), cld_lwp(j1:j1,:), & & cld_swp(j1:j1,:), cld_ref_ice(j1:j1,:), cld_ref_liq(j1:j1,:), & @@ -1040,7 +1040,7 @@ subroutine rrtmg_sw_run & tavel(k) = tlyr(j1,kk) delp (k) = delpin(j1,kk) dz (k) = dzlyr (j1,kk) - if (iovrsw == 4 .or. iovrsw == 5) alph(k) = alpha(j1,k) ! alpha decorrelation + if (iovr == 4 .or. iovr == 5) alph(k) = alpha(j1,k) ! alpha decorrelation !> -# Set absorber and gas column amount, convert from volume mixing !! ratio to molec/cm2 based on coldry (scaled to 1.0e-20) @@ -1111,7 +1111,7 @@ subroutine rrtmg_sw_run & cdat3(k) = cld_swp(j1,kk) ! cloud snow path cdat4(k) = cld_ref_snow(j1,kk) ! snow partical effctive radius enddo - if (iovrsw == 4) then !mz* HWRF + if (iovr == 4) then !mz* HWRF do k = 1, nlay kk = nlp1 - k do ig = 1, ngptsw @@ -1153,7 +1153,7 @@ subroutine rrtmg_sw_run & tavel(k) = tlyr(j1,k) delp (k) = delpin(j1,k) dz (k) = dzlyr (j1,k) - if (iovrsw == 4 .or. iovrsw == 5) alph(k) = alpha(j1,k) ! alpha decorrelation + if (iovr == 4 .or. iovr == 5) alph(k) = alpha(j1,k) ! alpha decorrelation ! --- ... set absorber amount !test use @@ -1226,7 +1226,7 @@ subroutine rrtmg_sw_run & cdat3(k) = cld_swp(j1,k) ! cloud snow path cdat4(k) = cld_ref_snow(j1,k) ! snow partical effctive radius enddo - if (iovrsw == 4) then !mz* HWRF + if (iovr == 4) then !mz* HWRF !mz* Move incoming GCM cloud arrays to RRTMG cloud arrays. !For GCM input, incoming reicmcl is defined based on selected !ice parameterization (inflglw) @@ -1269,11 +1269,11 @@ subroutine rrtmg_sw_run & zcf0 = f_one zcf1 = f_one - if (iovrsw == 0) then ! random overlapping + if (iovr == 0) then ! random overlapping do k = 1, nlay zcf0 = zcf0 * (f_one - cfrac(k)) enddo - else if (iovrsw == 1 .or. iovrsw == 4) then ! max/ran/exp overlapping + else if (iovr == 1 .or. iovr == 4) then ! max/ran/exp overlapping do k = 1, nlay if (cfrac(k) > ftiny) then ! cloudy layer zcf1 = min ( zcf1, f_one-cfrac(k) ) @@ -1283,7 +1283,7 @@ subroutine rrtmg_sw_run & endif enddo zcf0 = zcf0 * zcf1 - else if (iovrsw >= 2 .and. iovrsw /= 4) then + else if (iovr >= 2 .and. iovr /= 4) then do k = 1, nlay zcf0 = min ( zcf0, f_one-cfrac(k) ) ! used only as clear/cloudy indicator enddo @@ -1299,7 +1299,7 @@ subroutine rrtmg_sw_run & if (zcf1 > f_zero) then ! cloudy sky column !mz* for HWRF, save cldfmc with mcica - if (iovrsw == 4) then + if (iovr == 4) then do k = 1, nlay do ig = 1, ngptsw cldfmc_save(k,ig)=cldfmc (k,ig) @@ -1315,7 +1315,7 @@ subroutine rrtmg_sw_run & & taucw, ssacw, asycw, cldfrc, cldfmc & & ) - if (iovrsw == 4) then + if (iovr == 4) then !mz for HWRF, still using mcica cldfmc do k = 1, nlay do ig = 1, ngptsw @@ -1612,7 +1612,7 @@ subroutine rswinit & ! icldflg - cloud scheme control flag ! ! =0: diagnostic scheme gives cloud tau, omiga, and g. ! ! =1: prognostic scheme gives cloud liq/ice path, etc. ! -! iovrsw - clouds vertical overlapping control flag ! +! iovr - clouds vertical overlapping control flag ! ! =0: random overlapping clouds ! ! =1: maximum/random overlapping clouds ! ! =2: maximum overlap cloud ! @@ -1648,9 +1648,9 @@ subroutine rswinit & ! !===> ... begin here ! - if ( iovrsw<0 .or. iovrsw>4 ) then + if ( iovr<0 .or. iovr>4 ) then print *,' *** Error in specification of cloud overlap flag', & - & ' IOVRSW=',iovrsw,' in RSWINIT !!' + & ' IOVR=',iovr,' in RSWINIT !!' stop endif @@ -1697,15 +1697,15 @@ subroutine rswinit & stop endif - if ( isubcsw==0 .and. iovrsw>2 ) then + if ( isubcsw==0 .and. iovr>2 ) then if (me == 0) then - print *,' *** IOVRSW=',iovrsw,' is not available for', & + print *,' *** IOVR=',iovr,' is not available for', & & ' ISUBCSW=0 setting!!' print *,' The program will use maximum/random overlap', & & ' instead.' endif - iovrsw = 1 + iovr = 1 endif !> -# Setup constant factors for heating rate @@ -2116,7 +2116,7 @@ subroutine cldprop & !> -# if physparam::isubcsw > 0, call mcica_subcol() to distribute !! cloud properties to each g-point. - if ( isubcsw > 0 .and. iovrsw /= 4 ) then ! mcica sub-col clouds approx + if ( isubcsw > 0 .and. iovr /= 4 ) then ! mcica sub-col clouds approx cldf(:) = cfrac(:) where (cldf(:) < ftiny) @@ -2189,7 +2189,7 @@ subroutine mcica_subcol & ! lcloudy - logical, sub-colum cloud profile flag array nlay*ngptsw! ! ! ! other control flags from module variables: ! -! iovrsw : control flag for cloud overlapping method ! +! iovr : control flag for cloud overlapping method ! ! =0: random ! ! =1: maximum/random overlapping clouds ! ! =2: maximum overlap cloud ! @@ -2233,7 +2233,7 @@ subroutine mcica_subcol & !> -# Sub-column set up according to overlapping assumption. - select case ( iovrsw ) + select case ( iovr ) case( 0 ) ! random overlap, pick a random value at every level diff --git a/physics/rrtmg_lw_cloud_optics.F90 b/physics/rrtmg_lw_cloud_optics.F90 index ea0a703c7..02f32096a 100644 --- a/physics/rrtmg_lw_cloud_optics.F90 +++ b/physics/rrtmg_lw_cloud_optics.F90 @@ -1,6 +1,5 @@ module mo_rrtmg_lw_cloud_optics use machine, only: kind_phys - use physparam, only: ilwcliq, ilwcice, iovrlw use mersenne_twister, only: random_setseed, random_number, random_stat implicit none @@ -554,13 +553,15 @@ module mo_rrtmg_lw_cloud_optics ! subroutine rrtmg_lw_cloud_optics ! ####################################################################################### subroutine rrtmg_lw_cloud_optics(ncol, nlay, nBandsLW, cld_lwp, cld_ref_liq, cld_iwp, & - cld_ref_ice, cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow, cld_frac, tau_cld, & - tau_precip) + cld_ref_ice, cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow, cld_frac, ilwcliq, & + ilwcice, tau_cld, tau_precip) ! Inputs integer,intent(in) :: & nBandsLW, & ! Number of spectral bands ncol, & ! Number of horizontal gridpoints - nlay ! Number of vertical layers + nlay, & ! Number of vertical layers + ilwcliq, & ! + ilwcice real(kind_phys), dimension(ncol,nlay), intent(in) :: & cld_frac, & ! Cloud-fraction (1) cld_lwp, & ! Cloud liquid water path (g/m2) diff --git a/physics/rrtmg_sw_cloud_optics.F90 b/physics/rrtmg_sw_cloud_optics.F90 index 37b4e094c..01cab76e2 100644 --- a/physics/rrtmg_sw_cloud_optics.F90 +++ b/physics/rrtmg_sw_cloud_optics.F90 @@ -1,6 +1,5 @@ module mo_rrtmg_sw_cloud_optics use machine, only: kind_phys - use physparam, only: iswcliq, iswcice, iovrsw use mersenne_twister, only: random_setseed, random_number, random_stat implicit none @@ -2044,13 +2043,15 @@ module mo_rrtmg_sw_cloud_optics ! rrtmg_sw_cloud_optics ! ######################################################################################### subroutine rrtmg_sw_cloud_optics(ncol, nlay, nBandsSW, cld_lwp, cld_ref_liq, cld_iwp, & - cld_ref_ice, cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow, cld_frac, & - tau_cld, ssa_cld, asy_cld, tau_precip, ssa_precip, asy_precip) + cld_ref_ice, cld_rwp, cld_ref_rain, cld_swp, cld_ref_snow, cld_frac, iswcliq, & + iswcice, tau_cld, ssa_cld, asy_cld, tau_precip, ssa_precip, asy_precip) ! Inputs integer,intent(in) :: & nBandsSW, & ! Number of spectral bands ncol, & ! Number of horizontal gridpoints - nlay ! Number of vertical layers + nlay, & ! Number of vertical layers + iswcliq, & ! + iswcice ! real(kind_phys), dimension(ncol,nlay), intent(in) :: & cld_frac, & ! Cloud-fraction (1) cld_lwp, & ! Cloud liquid water path (g/m2) diff --git a/physics/rrtmgp_lw_cloud_optics.F90 b/physics/rrtmgp_lw_cloud_optics.F90 index 93e38994b..f45f08dd1 100644 --- a/physics/rrtmgp_lw_cloud_optics.F90 +++ b/physics/rrtmgp_lw_cloud_optics.F90 @@ -270,11 +270,11 @@ end subroutine rrtmgp_lw_cloud_optics_init !! \section arg_table_rrtmgp_lw_cloud_optics_run !! \htmlinclude rrtmgp_lw_cloud_optics.html !! - subroutine rrtmgp_lw_cloud_optics_run(doLWrad, doG_cldoptics, doGP_cldoptics_PADE, & - doGP_cldoptics_LUT, nCol, nLev, nrghice, p_lay, cld_frac, cld_lwp, cld_reliq, & - cld_iwp, cld_reice, cld_swp, cld_resnow, cld_rwp, cld_rerain, precip_frac, & - lw_cloud_props, lw_gas_props, lon, lat, cldtaulw, lw_optical_props_cloudsByBand, & - lw_optical_props_precipByBand, errmsg, errflg) + subroutine rrtmgp_lw_cloud_optics_run(doLWrad, doG_cldoptics, icliq_lw, icice_lw, & + doGP_cldoptics_PADE, doGP_cldoptics_LUT, nCol, nLev, nrghice, p_lay, cld_frac, & + cld_lwp, cld_reliq, cld_iwp, cld_reice, cld_swp, cld_resnow, cld_rwp, cld_rerain, & + precip_frac, lw_cloud_props, lw_gas_props, lon, lat, cldtaulw, & + lw_optical_props_cloudsByBand, lw_optical_props_precipByBand, errmsg, errflg) ! Inputs logical, intent(in) :: & @@ -285,7 +285,9 @@ subroutine rrtmgp_lw_cloud_optics_run(doLWrad, doG_cldoptics, doGP_cldoptics_PAD integer, intent(in) :: & nCol, & ! Number of horizontal gridpoints nLev, & ! Number of vertical levels - nrghice ! Number of ice-roughness categories + nrghice, & ! Number of ice-roughness categories + icliq_lw, & ! Choice of treatment of liquid cloud optical properties (RRTMG legacy) + icice_lw ! Choice of treatment of ice cloud optical properties (RRTMG legacy) real(kind_phys), dimension(nCol), intent(in) :: & lon, & ! Longitude lat ! Latitude @@ -378,7 +380,7 @@ subroutine rrtmgp_lw_cloud_optics_run(doLWrad, doG_cldoptics, doGP_cldoptics_PAD if (any(cld_frac .gt. 0)) then call rrtmg_lw_cloud_optics(ncol, nLev, lw_gas_props%get_nband(), cld_lwp, & cld_reliq, cld_iwp, cld_reice, cld_rwp, cld_rerain, cld_swp, cld_resnow, & - cld_frac, tau_cld, tau_precip) + cld_frac, icliq_lw, icice_lw, tau_cld, tau_precip) endif lw_optical_props_cloudsByBand%tau = tau_cld lw_optical_props_precipByBand%tau = tau_precip diff --git a/physics/rrtmgp_lw_cloud_optics.meta b/physics/rrtmgp_lw_cloud_optics.meta index d6575fa14..809e8abf0 100644 --- a/physics/rrtmgp_lw_cloud_optics.meta +++ b/physics/rrtmgp_lw_cloud_optics.meta @@ -127,6 +127,22 @@ type = logical intent = in optional = F +[icliq_lw] + standard_name = flag_for_optical_property_for_liquid_clouds_for_longwave_radiation + long_name = lw optical property for liquid clouds + units = flag + dimensions = () + type = integer + intent = in + optional = F +[icice_lw] + standard_name = flag_for_optical_property_for_ice_clouds_for_longwave_radiation + long_name = lw optical property for ice clouds + units = flag + dimensions = () + type = integer + intent = in + optional = F [doGP_cldoptics_PADE] standard_name = flag_to_calc_lw_cld_optics_using_RRTMGP_PADE long_name = logical flag to control cloud optics scheme. diff --git a/physics/rrtmgp_lw_cloud_sampling.F90 b/physics/rrtmgp_lw_cloud_sampling.F90 index 35ae3c4a8..cfb86eb3a 100644 --- a/physics/rrtmgp_lw_cloud_sampling.F90 +++ b/physics/rrtmgp_lw_cloud_sampling.F90 @@ -1,7 +1,6 @@ module rrtmgp_lw_cloud_sampling use machine, only: kind_phys use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp - use physparam, only: isubclw, iovrlw use mo_optical_props, only: ty_optical_props_1scl use rrtmgp_sampling, only: sampled_mask, draw_samples use mersenne_twister, only: random_setseed, random_number, random_stat @@ -45,7 +44,8 @@ end subroutine rrtmgp_lw_cloud_sampling_init !! \section arg_table_rrtmgp_lw_cloud_sampling_run !! \htmlinclude rrtmgp_lw_cloud_sampling_run.html !! - subroutine rrtmgp_lw_cloud_sampling_run(doLWrad, nCol, nLev, ipsdlw0, icseed_lw, & + subroutine rrtmgp_lw_cloud_sampling_run(doLWrad, nCol, nLev, ipsdlw0, icseed_lw, iovr, & + iovr_max, iovr_maxrand, iovr_rand, iovr_dcorr, iovr_exp, iovr_exprand, isubc_lw, & cld_frac, precip_frac, cloud_overlap_param, precip_overlap_param, lw_gas_props, & lw_optical_props_cloudsByBand, lw_optical_props_precipByBand, & lw_optical_props_clouds, lw_optical_props_precip, errmsg, errflg) @@ -56,12 +56,20 @@ subroutine rrtmgp_lw_cloud_sampling_run(doLWrad, nCol, nLev, ipsdlw0, icseed_lw, integer, intent(in) :: & nCol, & ! Number of horizontal gridpoints nLev, & ! Number of vertical layers - ipsdlw0 ! Initial permutation seed for McICA + iovr, & ! Choice of cloud-overlap method + iovr_max, & ! Flag for maximum cloud overlap method + iovr_maxrand, & ! Flag for maximum-random cloud overlap method + iovr_rand, & ! Flag for random cloud overlap method + iovr_dcorr, & ! Flag for decorrelation-length cloud overlap method + iovr_exp, & ! Flag for exponential cloud overlap method + iovr_exprand, & ! Flag for exponential-random cloud overlap method + ipsdlw0, & ! Initial permutation seed for McICA + isubc_lw integer,intent(in),dimension(ncol) :: & icseed_lw ! auxiliary special cloud related array when module - ! variable isubclw=2, it provides permutation seed + ! variable isubc_lw=2, it provides permutation seed ! for each column profile that are used for generating - ! random numbers. when isubclw /=2, it will not be used. + ! random numbers. when isubc_lw /=2, it will not be used. real(kind_phys), dimension(ncol,nLev),intent(in) :: & cld_frac, & ! Total cloud fraction by layer precip_frac ! Precipitation fraction by layer @@ -84,24 +92,17 @@ subroutine rrtmgp_lw_cloud_sampling_run(doLWrad, nCol, nLev, ipsdlw0, icseed_lw, lw_optical_props_precip ! RRTMGP DDT: Shortwave optical properties by spectral point (precipitation) ! Local variables - integer :: iCol + integer :: iCol, iLay integer,dimension(ncol) :: ipseed_lw type(random_stat) :: rng_stat real(kind_phys), dimension(lw_gas_props%get_ngpt(),nLev,ncol) :: rng3D,rng3D2 - real(kind_phys), dimension(lw_gas_props%get_ngpt()*nLev) :: rng1D + real(kind_phys), dimension(lw_gas_props%get_ngpt()) :: rng1D + real(kind_phys), dimension(lw_gas_props%get_ngpt()*nLev) :: rng2D logical, dimension(ncol,nLev,lw_gas_props%get_ngpt()) :: cldfracMCICA,precipfracSAMP ! Initialize CCPP error handling variables errmsg = '' errflg = 0 - - ! - if (iovrlw .ne. 1 .and. iovrlw .ne. 3 .and. iovrlw .ne. 4 .and. iovrlw .ne. 5) then - errmsg = 'Cloud overlap assumption not supported.' - errflg = 1 - call check_error_msg('rrtmgp_lw_cloud_sampling',errmsg) - return - endif if (.not. doLWrad) return @@ -113,12 +114,12 @@ subroutine rrtmgp_lw_cloud_sampling_run(doLWrad, nCol, nLev, ipsdlw0, icseed_lw, call check_error_msg('rrtmgp_lw_cloud_sampling_run',& lw_optical_props_clouds%alloc_1scl(nCol, nLev, lw_gas_props)) - ! Change random number seed value for each radiation invocation (isubclw =1 or 2). - if(isubclw == 1) then ! advance prescribed permutation seed + ! Change random number seed value for each radiation invocation (isubc_lw =1 or 2). + if(isubc_lw == 1) then ! advance prescribed permutation seed do iCol = 1, ncol ipseed_lw(iCol) = ipsdlw0 + iCol enddo - elseif (isubclw == 2) then ! use input array of permutaion seeds + elseif (isubc_lw == 2) then ! use input array of permutaion seeds do iCol = 1, ncol ipseed_lw(iCol) = icseed_lw(iCol) enddo @@ -128,34 +129,46 @@ subroutine rrtmgp_lw_cloud_sampling_run(doLWrad, nCol, nLev, ipsdlw0, icseed_lw, ! and layers. ([nGpts,nLev,nColumn]-> [nGpts*nLev]*nColumn) do iCol=1,ncol call random_setseed(ipseed_lw(icol),rng_stat) - call random_number(rng1D,rng_stat) - rng3D(:,:,iCol) = reshape(source = rng1D,shape=[lw_gas_props%get_ngpt(),nLev]) + ! Use same rng for each layer + if (iovr == iovr_max) then + call random_number(rng1D,rng_stat) + do iLay=1,nLev + rng3D(:,iLay,iCol) = rng1D + enddo + else + do iLay=1,nLev + call random_number(rng1D,rng_stat) + rng3D(:,iLay,iCol) = rng1D + enddo + endif enddo ! Cloud-overlap. - ! Maximum-random - if (iovrlw == 1) then + ! Maximum-random, random or maximum. + if (iovr == iovr_maxrand .or. iovr == iovr_rand .or. iovr == iovr_max) then call sampled_mask(rng3D, cld_frac, cldfracMCICA) endif ! Exponential decorrelation length overlap - if (iovrlw == 3) then + if (iovr == iovr_dcorr) then ! Generate second RNG do iCol=1,ncol call random_setseed(ipseed_lw(icol),rng_stat) - call random_number(rng1D,rng_stat) - rng3D2(:,:,iCol) = reshape(source = rng1D,shape=[lw_gas_props%get_ngpt(),nLev]) + call random_number(rng2D,rng_stat) + rng3D2(:,:,iCol) = reshape(source = rng2D,shape=[lw_gas_props%get_ngpt(),nLev]) enddo call sampled_mask(rng3D, cld_frac, cldfracMCICA, & overlap_param = cloud_overlap_param(:,1:nLev-1), & randoms2 = rng3D2) endif ! Exponential or Exponential-random - if (iovrlw == 4 .or. iovrlw == 5) then + if (iovr == iovr_exp .or. iovr == iovr_exprand) then call sampled_mask(rng3D, cld_frac, cldfracMCICA, & overlap_param = cloud_overlap_param(:,1:nLev-1)) endif + ! ! Sampling. Map band optical depth to each g-point using McICA + ! call check_error_msg('rrtmgp_lw_cloud_sampling_run_draw_samples',& draw_samples(cldfracMCICA, & lw_optical_props_cloudsByBand, & @@ -169,12 +182,12 @@ subroutine rrtmgp_lw_cloud_sampling_run(doLWrad, nCol, nLev, ipsdlw0, icseed_lw, call check_error_msg('rrtmgp_lw_cloud_sampling_run',& lw_optical_props_precip%alloc_1scl(nCol, nLev, lw_gas_props)) - ! Change random number seed value for each radiation invocation (isubclw =1 or 2). - if(isubclw == 1) then ! advance prescribed permutation seed + ! Change random number seed value for each radiation invocation (isubc_lw =1 or 2). + if(isubc_lw == 1) then ! advance prescribed permutation seed do iCol = 1, ncol ipseed_lw(iCol) = ipsdlw0 + iCol enddo - elseif (isubclw == 2) then ! use input array of permutaion seeds + elseif (isubc_lw == 2) then ! use input array of permutaion seeds do iCol = 1, ncol ipseed_lw(iCol) = icseed_lw(iCol) enddo @@ -190,12 +203,12 @@ subroutine rrtmgp_lw_cloud_sampling_run(doLWrad, nCol, nLev, ipsdlw0, icseed_lw, !enddo ! Precipitation overlap. - ! Maximum-random - if (iovrlw == 1) then + ! Maximum-random, random or maximum. + if (iovr == iovr_maxrand .or. iovr == iovr_rand .or. iovr == iovr_max) then call sampled_mask(rng3D, precip_frac, precipfracSAMP) endif - ! Exponential decorrelation length overlap - if (iovrlw == 3) then + ! Exponential decorrelation length overlap + if (iovr == iovr_dcorr) then ! No need to call RNG second time for now, just use the same seeds for precip as clouds. !! Generate second RNG !do iCol=1,ncol @@ -208,13 +221,14 @@ subroutine rrtmgp_lw_cloud_sampling_run(doLWrad, nCol, nLev, ipsdlw0, icseed_lw, randoms2 = rng3D2) endif ! Exponential or Exponential-random - if (iovrlw == 4 .or. iovrlw == 5) then + if (iovr == iovr_exp .or. iovr == iovr_exprand) then call sampled_mask(rng3D, precip_frac, precipfracSAMP, & overlap_param = precip_overlap_param(:,1:nLev-1)) endif - + ! ! Sampling. Map band optical depth to each g-point using McICA + ! call check_error_msg('rrtmgp_lw_precip_sampling_run_draw_samples',& draw_samples(precipfracSAMP, & lw_optical_props_precipByBand, & diff --git a/physics/rrtmgp_lw_cloud_sampling.meta b/physics/rrtmgp_lw_cloud_sampling.meta index 35699efb6..54f3c63af 100644 --- a/physics/rrtmgp_lw_cloud_sampling.meta +++ b/physics/rrtmgp_lw_cloud_sampling.meta @@ -69,6 +69,14 @@ type = integer intent = in optional = F +[isubc_lw] + standard_name = flag_for_lw_clouds_sub_grid_approximation + long_name = flag for lw clouds sub-grid approximation + units = flag + dimensions = () + type = integer + intent = in + optional = F [ipsdlw0] standard_name = initial_permutation_seed_lw long_name = initial seed for McICA LW @@ -77,6 +85,62 @@ type = integer intent = in optional = F +[iovr] + standard_name = flag_for_cloud_overlap_method_for_radiation + long_name = max-random overlap clouds + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_maxrand] + standard_name = flag_for_maximum_random_cloud_overlap_method + long_name = choice of maximum-random cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_dcorr] + standard_name = flag_for_decorrelation_length_cloud_overlap_method + long_name = choice of decorrelation-length cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_exp] + standard_name = flag_for_exponential_cloud_overlap_method + long_name = choice of exponential cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_exprand] + standard_name = flag_for_exponential_random_cloud_overlap_method + long_name = choice of exponential-random cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_rand] + standard_name = flag_for_random_cloud_overlap_method + long_name = choice of random cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_max] + standard_name = flag_for_maximum_cloud_overlap_method + long_name = choice of maximum cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F [icseed_lw] standard_name = seed_random_numbers_lw long_name = seed for random number generation for longwave radiation diff --git a/physics/rrtmgp_lw_pre.F90 b/physics/rrtmgp_lw_pre.F90 index 7ad8bd30d..caee7308e 100644 --- a/physics/rrtmgp_lw_pre.F90 +++ b/physics/rrtmgp_lw_pre.F90 @@ -1,5 +1,4 @@ module rrtmgp_lw_pre - use physparam use machine, only: & kind_phys ! Working type use module_radiation_surface, only: & diff --git a/physics/rrtmgp_sw_cloud_optics.F90 b/physics/rrtmgp_sw_cloud_optics.F90 index 7ab3c27e3..505fe7853 100644 --- a/physics/rrtmgp_sw_cloud_optics.F90 +++ b/physics/rrtmgp_sw_cloud_optics.F90 @@ -3,7 +3,6 @@ module rrtmgp_sw_cloud_optics use mo_rte_kind, only: wl use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp use mo_cloud_optics, only: ty_cloud_optics - use physparam, only: isubcsw, iovrsw use mo_optical_props, only: ty_optical_props_2str use mo_rrtmg_sw_cloud_optics, only: rrtmg_sw_cloud_optics use rrtmgp_aux, only: check_error_msg @@ -287,10 +286,10 @@ end subroutine rrtmgp_sw_cloud_optics_init !! \section arg_table_rrtmgp_sw_cloud_optics_run !! \htmlinclude rrtmgp_sw_cloud_optics.html !! - subroutine rrtmgp_sw_cloud_optics_run(doSWrad, doG_cldoptics, doGP_cldoptics_PADE, & - doGP_cldoptics_LUT, nCol, nLev, nDay, idxday, nrghice, cld_frac, cld_lwp, cld_reliq, & - cld_iwp, cld_reice, cld_swp, cld_resnow, cld_rwp, cld_rerain, precip_frac, & - sw_cloud_props, sw_gas_props, sw_optical_props_cloudsByBand, & + subroutine rrtmgp_sw_cloud_optics_run(doSWrad, doG_cldoptics, icliq_sw, icice_sw, & + doGP_cldoptics_PADE, doGP_cldoptics_LUT, nCol, nLev, nDay, idxday, nrghice, cld_frac,& + cld_lwp, cld_reliq, cld_iwp, cld_reice, cld_swp, cld_resnow, cld_rwp, cld_rerain, & + precip_frac, sw_cloud_props, sw_gas_props, sw_optical_props_cloudsByBand, & sw_optical_props_precipByBand, cldtausw, errmsg, errflg) ! Inputs @@ -303,7 +302,9 @@ subroutine rrtmgp_sw_cloud_optics_run(doSWrad, doG_cldoptics, doGP_cldoptics_PAD nCol, & ! Number of horizontal gridpoints nLev, & ! Number of vertical levels nday, & ! Number of daylit points. - nrghice ! Number of ice-roughness categories + nrghice, & ! Number of ice-roughness categories + icliq_sw, & ! Choice of treatment of liquid cloud optical properties (RRTMG legacy) + icice_sw ! Choice of treatment of ice cloud optical properties (RRTMG legacy) integer,intent(in),dimension(ncol) :: & idxday ! Indices for daylit points. real(kind_phys), dimension(ncol,nLev),intent(in) :: & @@ -417,7 +418,7 @@ subroutine rrtmgp_sw_cloud_optics_run(doSWrad, doG_cldoptics, doGP_cldoptics_PAD cld_iwp(idxday(1:nday),:), cld_reice(idxday(1:nday),:), & cld_rwp(idxday(1:nday),:), cld_rerain(idxday(1:nday),:), & cld_swp(idxday(1:nday),:), cld_resnow(idxday(1:nday),:), & - cld_frac(idxday(1:nday),:), & + cld_frac(idxday(1:nday),:), icliq_sw, icice_sw, & tau_cld, ssa_cld, asy_cld, & tau_precip, ssa_precip, asy_precip) diff --git a/physics/rrtmgp_sw_cloud_optics.meta b/physics/rrtmgp_sw_cloud_optics.meta index 0251120d3..4439a607b 100644 --- a/physics/rrtmgp_sw_cloud_optics.meta +++ b/physics/rrtmgp_sw_cloud_optics.meta @@ -143,6 +143,22 @@ type = logical intent = in optional = F +[icliq_sw] + standard_name = flag_for_optical_property_for_liquid_clouds_for_shortwave_radiation + long_name = sw optical property for liquid clouds + units = flag + dimensions = () + type = integer + intent = in + optional = F +[icice_sw] + standard_name = flag_for_optical_property_for_ice_clouds_for_shortwave_radiation + long_name = sw optical property for ice clouds + units = flag + dimensions = () + type = integer + intent = in + optional = F [doGP_cldoptics_PADE] standard_name = flag_to_calc_lw_cld_optics_using_RRTMGP_PADE long_name = logical flag to control cloud optics scheme. diff --git a/physics/rrtmgp_sw_cloud_sampling.F90 b/physics/rrtmgp_sw_cloud_sampling.F90 index 802cad840..ba4097e96 100644 --- a/physics/rrtmgp_sw_cloud_sampling.F90 +++ b/physics/rrtmgp_sw_cloud_sampling.F90 @@ -1,7 +1,6 @@ module rrtmgp_sw_cloud_sampling use machine, only: kind_phys use mo_gas_optics_rrtmgp, only: ty_gas_optics_rrtmgp - use physparam, only: isubcsw, iovrsw use mo_optical_props, only: ty_optical_props_2str use rrtmgp_sampling, only: sampled_mask, draw_samples use mersenne_twister, only: random_setseed, random_number, random_stat @@ -44,7 +43,8 @@ end subroutine rrtmgp_sw_cloud_sampling_init !! \section arg_table_rrtmgp_sw_cloud_sampling_run !! \htmlinclude rrtmgp_sw_cloud_sampling.html !! - subroutine rrtmgp_sw_cloud_sampling_run(doSWrad, nCol, nDay, nLev, ipsdsw0, idxday, & + subroutine rrtmgp_sw_cloud_sampling_run(doSWrad, nCol, nDay, nLev, ipsdsw0, idxday, iovr, & + iovr_max, iovr_maxrand, iovr_rand, iovr_dcorr, iovr_exp, iovr_exprand, isubc_sw, & icseed_sw, cld_frac, precip_frac, cloud_overlap_param, precip_overlap_param, & sw_gas_props, sw_optical_props_cloudsByBand, sw_optical_props_precipByBand, & sw_optical_props_clouds, sw_optical_props_precip, errmsg, errflg) @@ -56,14 +56,22 @@ subroutine rrtmgp_sw_cloud_sampling_run(doSWrad, nCol, nDay, nLev, ipsdsw0, idxd nCol, & ! Number of horizontal gridpoints nDay, & ! Number of daylit points. nLev, & ! Number of vertical layers - ipsdsw0 ! Initial permutation seed for McICA + ipsdsw0, & ! Initial permutation seed for McICA + iovr, & ! Choice of cloud-overlap method + iovr_max, & ! Flag for maximum cloud overlap method + iovr_maxrand, & ! Flag for maximum-random cloud overlap method + iovr_rand, & ! Flag for random cloud overlap method + iovr_dcorr, & ! Flag for decorrelation-length cloud overlap method + iovr_exp, & ! Flag for exponential cloud overlap method + iovr_exprand, & ! Flag for exponential-random cloud overlap method + isubc_sw integer,intent(in),dimension(ncol) :: & idxday ! Indices for daylit points. integer,intent(in),dimension(ncol) :: & icseed_sw ! auxiliary special cloud related array when module - ! variable isubcsw=2, it provides permutation seed + ! variable isubc_sw=2, it provides permutation seed ! for each column profile that are used for generating - ! random numbers. when isubcsw /=2, it will not be used. + ! random numbers. when isubc_sw /=2, it will not be used. real(kind_phys), dimension(ncol,nLev),intent(in) :: & cld_frac, & ! Total cloud fraction by layer precip_frac ! Precipitation fraction by layer @@ -91,21 +99,14 @@ subroutine rrtmgp_sw_cloud_sampling_run(doSWrad, nCol, nDay, nLev, ipsdsw0, idxd type(random_stat) :: rng_stat real(kind_phys) :: tauloc,asyloc,ssaloc real(kind_phys), dimension(sw_gas_props%get_ngpt(),nLev,nday) :: rng3D,rng3D2 - real(kind_phys), dimension(sw_gas_props%get_ngpt()*nLev) :: rng1D + real(kind_phys), dimension(sw_gas_props%get_ngpt()*nLev) :: rng2D + real(kind_phys), dimension(sw_gas_props%get_ngpt()) :: rng1D logical, dimension(nday,nLev,sw_gas_props%get_ngpt()) :: cldfracMCICA,precipfracSAMP ! Initialize CCPP error handling variables errmsg = '' errflg = 0 - ! Only works w/ SDFs v15p2 and v16beta - if (iovrsw .ne. 1 .and. iovrsw .ne. 3 .and. iovrsw .ne. 4 .and. iovrsw .ne. 5) then - errmsg = 'Cloud overlap assumption not supported.' - errflg = 1 - call check_error_msg('rrtmgp_sw_cloud_sampling',errmsg) - return - endif - if (.not. doSWrad) return if (nDay .gt. 0) then ! ################################################################################# @@ -115,13 +116,16 @@ subroutine rrtmgp_sw_cloud_sampling_run(doSWrad, nCol, nDay, nLev, ipsdsw0, idxd ! Allocate space RRTMGP DDTs [nday,nLev,nGpt] call check_error_msg('rrtmgp_sw_cloud_sampling_run', & sw_optical_props_clouds%alloc_2str(nday, nLev, sw_gas_props)) + sw_optical_props_clouds%tau(:,:,:) = 0._kind_phys + sw_optical_props_clouds%ssa(:,:,:) = 1._kind_phys + sw_optical_props_clouds%g(:,:,:) = 0._kind_phys - ! Change random number seed value for each radiation invocation (isubcsw =1 or 2). - if(isubcsw == 1) then ! advance prescribed permutation seed + ! Change random number seed value for each radiation invocation (isubc_sw =1 or 2). + if(isubc_sw == 1) then ! advance prescribed permutation seed do iday = 1, nday ipseed_sw(iday) = ipsdsw0 + iday enddo - elseif (isubcsw == 2) then ! use input array of permutaion seeds + elseif (isubc_sw == 2) then ! use input array of permutaion seeds do iday = 1, nday ipseed_sw(iday) = icseed_sw(iday) enddo @@ -131,33 +135,51 @@ subroutine rrtmgp_sw_cloud_sampling_run(doSWrad, nCol, nDay, nLev, ipsdsw0, idxd ! and layers. ([nGpts,nLev,nDayumn]-> [nGpts*nLev]*nDayumn) do iday=1,nday call random_setseed(ipseed_sw(iday),rng_stat) - call random_number(rng1D,rng_stat) - rng3D(:,:,iday) = reshape(source = rng1D,shape=[sw_gas_props%get_ngpt(),nLev]) + ! Use same rng for each layer + if (iovr == iovr_max) then + call random_number(rng1D,rng_stat) + do iLay=1,nLev + rng3D(:,iLay,iday) = rng1D + enddo + else + do iLay=1,nLev + call random_number(rng1D,rng_stat) + rng3D(:,iLay,iday) = rng1D + enddo + endif + enddo + + do iday=1,nday + call random_setseed(ipseed_sw(iday),rng_stat) + call random_number(rng2D,rng_stat) + rng3D(:,:,iday) = reshape(source = rng2D,shape=[sw_gas_props%get_ngpt(),nLev]) enddo ! Cloud overlap. - ! Maximum-random overlap - if (iovrsw == 1) then + ! Maximum-random, random, or maximum cloud overlap + if (iovr == iovr_maxrand .or. iovr == iovr_max .or. iovr == iovr_rand) then call sampled_mask(rng3D, cld_frac(idxday(1:nDay),:), cldfracMCICA) endif ! Decorrelation-length overlap - if (iovrsw == 3) then + if (iovr == iovr_dcorr) then do iday=1,nday call random_setseed(ipseed_sw(iday),rng_stat) - call random_number(rng1D,rng_stat) - rng3D2(:,:,iday) = reshape(source = rng1D,shape=[sw_gas_props%get_ngpt(),nLev]) + call random_number(rng2D,rng_stat) + rng3D2(:,:,iday) = reshape(source = rng2D,shape=[sw_gas_props%get_ngpt(),nLev]) enddo call sampled_mask(rng3D, cld_frac(idxday(1:nDay),:), cldfracMCICA, & overlap_param = cloud_overlap_param(idxday(1:nDay),1:nLev-1),& randoms2 = rng3D2) endif - ! Exponential overlap - if (iovrsw == 4 .or. iovrsw == 5) then + ! Exponential or exponential-random cloud overlap + if (iovr == iovr_exp .or. iovr == iovr_exprand) then call sampled_mask(rng3D, cld_frac(idxday(1:nDay),:), cldfracMCICA, & overlap_param = cloud_overlap_param(idxday(1:nDay),1:nLev-1)) endif - + + ! ! Sampling. Map band optical depth to each g-point using McICA + ! call check_error_msg('rrtmgp_sw_cloud_sampling_run_draw_samples', & draw_samples(cldfracMCICA, & sw_optical_props_cloudsByBand, & @@ -171,12 +193,12 @@ subroutine rrtmgp_sw_cloud_sampling_run(doSWrad, nCol, nDay, nLev, ipsdsw0, idxd call check_error_msg('rrtmgp_sw_cloud_sampling_run', & sw_optical_props_precip%alloc_2str( nday, nLev, sw_gas_props)) - ! Change random number seed value for each radiation invocation (isubcsw =1 or 2). - if(isubcsw == 1) then ! advance prescribed permutation seed + ! Change random number seed value for each radiation invocation (isubc_sw =1 or 2). + if(isubc_sw == 1) then ! advance prescribed permutation seed do iday = 1, nday ipseed_sw(iday) = ipsdsw0 + iday enddo - elseif (isubcsw == 2) then ! use input array of permutaion seeds + elseif (isubc_sw == 2) then ! use input array of permutaion seeds do iday = 1, nday ipseed_sw(iday) = icseed_sw(iday) enddo @@ -192,12 +214,12 @@ subroutine rrtmgp_sw_cloud_sampling_run(doSWrad, nCol, nDay, nLev, ipsdsw0, idxd !enddo ! Precipitation overlap - ! Maximum-random - if (iovrsw == 1) then + ! Maximum-random, random or maximum precipitation overlap + if (iovr == iovr_maxrand .or. iovr == iovr_max .or. iovr == iovr_rand) then call sampled_mask(rng3D, precip_frac(idxday(1:nDay),:), precipfracSAMP) endif - ! Exponential decorrelation length overlap - if (iovrsw == 3) then + ! Exponential decorrelation length overlap + if (iovr == iovr_dcorr) then !! Generate second RNG !do iday=1,nday ! call random_setseed(ipseed_sw(iday),rng_stat) @@ -208,52 +230,55 @@ subroutine rrtmgp_sw_cloud_sampling_run(doSWrad, nCol, nDay, nLev, ipsdsw0, idxd overlap_param = precip_overlap_param(idxday(1:nDay),1:nLev-1),& randoms2 = rng3D2) endif - if (iovrsw == 4 .or. iovrsw == 5) then + if (iovr == iovr_exp .or. iovr == iovr_exprand) then call sampled_mask(rng3D, precip_frac(idxday(1:nDay),:),precipfracSAMP, & overlap_param = precip_overlap_param(idxday(1:nDay),1:nLev-1)) endif - ! Map band optical depth to each g-point using McICA + ! + ! Sampling. Map band optical depth to each g-point using McICA + ! call check_error_msg('rrtmgp_sw_precip_sampling_run_draw_samples', & draw_samples(precipfracSAMP, & sw_optical_props_precipByBand, & sw_optical_props_precip)) - endif - ! #################################################################################### - ! Just add precipitation optics to cloud-optics - ! #################################################################################### - do iGpt=1,sw_gas_props%get_ngpt() - do iday=1,nDay - do iLay=1,nLev - tauloc = sw_optical_props_clouds%tau(iday,iLay,iGpt) + & - sw_optical_props_precip%tau(iday,iLay,iGpt) - if (sw_optical_props_precip%tau(iday,iLay,iGpt) > 0) then - ssaloc = (sw_optical_props_clouds%tau(iday,iLay,iGpt) * & - sw_optical_props_clouds%ssa(iday,iLay,iGpt) + & - sw_optical_props_precip%tau(iday,iLay,iGpt) * & - sw_optical_props_precip%ssa(iday,iLay,iGpt)) / & - tauloc - if (ssaloc > 0) then - asyloc = (sw_optical_props_clouds%tau(iday,iLay,iGpt) * & - sw_optical_props_clouds%ssa(iday,iLay,iGpt) * & - sw_optical_props_clouds%g(iday,iLay,iGpt) + & - sw_optical_props_precip%tau(iday,iLay,iGpt) * & - sw_optical_props_precip%ssa(iday,iLay,iGpt) * & - sw_optical_props_precip%g(iday,iLay,iGpt)) / & - (tauloc*ssaloc) - else - tauloc = sw_optical_props_clouds%tau(iday,iLay,iGpt) - ssaloc = sw_optical_props_clouds%ssa(iday,iLay,iGpt) - asyloc = sw_optical_props_clouds%g(iday,iLay,iGpt) + ! ################################################################################# + ! Just add precipitation optics to cloud-optics + ! ################################################################################# + do iGpt=1,sw_gas_props%get_ngpt() + do iday=1,nDay + do iLay=1,nLev + tauloc = sw_optical_props_clouds%tau(iday,iLay,iGpt) + & + sw_optical_props_precip%tau(iday,iLay,iGpt) + if (sw_optical_props_precip%tau(iday,iLay,iGpt) > 0) then + ssaloc = (sw_optical_props_clouds%tau(iday,iLay,iGpt) * & + sw_optical_props_clouds%ssa(iday,iLay,iGpt) + & + sw_optical_props_precip%tau(iday,iLay,iGpt) * & + sw_optical_props_precip%ssa(iday,iLay,iGpt)) / & + tauloc + if (ssaloc > 0) then + asyloc = (sw_optical_props_clouds%tau(iday,iLay,iGpt) * & + sw_optical_props_clouds%ssa(iday,iLay,iGpt) * & + sw_optical_props_clouds%g(iday,iLay,iGpt) + & + sw_optical_props_precip%tau(iday,iLay,iGpt) * & + sw_optical_props_precip%ssa(iday,iLay,iGpt) * & + sw_optical_props_precip%g(iday,iLay,iGpt)) / & + (tauloc*ssaloc) + else + tauloc = sw_optical_props_clouds%tau(iday,iLay,iGpt) + ssaloc = sw_optical_props_clouds%ssa(iday,iLay,iGpt) + asyloc = sw_optical_props_clouds%g(iday,iLay,iGpt) + endif + sw_optical_props_clouds%tau(iday,iLay,iGpt) = tauloc + sw_optical_props_clouds%ssa(iday,iLay,iGpt) = ssaloc + sw_optical_props_clouds%g(iday,iLay,iGpt) = asyloc endif - sw_optical_props_clouds%tau(iday,iLay,iGpt) = tauloc - sw_optical_props_clouds%ssa(iday,iLay,iGpt) = ssaloc - sw_optical_props_clouds%g(iday,iLay,iGpt) = asyloc - endif + enddo enddo enddo - enddo + endif + end subroutine rrtmgp_sw_cloud_sampling_run ! ######################################################################################### diff --git a/physics/rrtmgp_sw_cloud_sampling.meta b/physics/rrtmgp_sw_cloud_sampling.meta index 082704462..01a311fd4 100644 --- a/physics/rrtmgp_sw_cloud_sampling.meta +++ b/physics/rrtmgp_sw_cloud_sampling.meta @@ -77,6 +77,14 @@ type = integer intent = in optional = F +[isubc_sw] + standard_name = flag_for_sw_clouds_grid_approximation + long_name = flag for sw clouds sub-grid approximation + units = flag + dimensions = () + type = integer + intent = in + optional = F [ipsdsw0] standard_name = initial_permutation_seed_sw long_name = initial seed for McICA SW @@ -93,6 +101,62 @@ type = integer intent = in optional = F +[iovr] + standard_name = flag_for_cloud_overlap_method_for_radiation + long_name = max-random overlap clouds + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_maxrand] + standard_name = flag_for_maximum_random_cloud_overlap_method + long_name = choice of maximum-random cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_dcorr] + standard_name = flag_for_decorrelation_length_cloud_overlap_method + long_name = choice of decorrelation-length cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_exp] + standard_name = flag_for_exponential_cloud_overlap_method + long_name = choice of exponential cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_exprand] + standard_name = flag_for_exponential_random_cloud_overlap_method + long_name = choice of exponential-random cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_rand] + standard_name = flag_for_random_cloud_overlap_method + long_name = choice of random cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F +[iovr_max] + standard_name = flag_for_maximum_cloud_overlap_method + long_name = choice of maximum cloud overlap method + units = flag + dimensions = () + type = integer + intent = in + optional = F [icseed_sw] standard_name = seed_random_numbers_sw long_name = seed for random number generation for shortwave radiation diff --git a/physics/rrtmgp_sw_gas_optics.F90 b/physics/rrtmgp_sw_gas_optics.F90 index efe611e0c..ac643e71d 100644 --- a/physics/rrtmgp_sw_gas_optics.F90 +++ b/physics/rrtmgp_sw_gas_optics.F90 @@ -374,6 +374,8 @@ subroutine rrtmgp_sw_gas_optics_run(doSWrad, nCol, nLev, nday, idxday, sw_gas_pr toa_src_sw(idxday(ij),:) = toa_src_sw(idxday(ij),:)*solcon/ & sum(toa_src_sw(idxday(ij),:)) enddo + else + toa_src_sw(:,:) = 0. endif end subroutine rrtmgp_sw_gas_optics_run