diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 90b7460a0..888087e56 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -876,6 +876,19 @@ subroutine progcld1 & alpha(:,:) = 0. endif + ! Revise alpha for exponential-random cloud overlap + ! Decorrelate layers when a clear layer follows a cloudy layer to enforce + ! random correlation between non-adjacent blocks of cloudy layers + if (iovr == 5) then + do k = 2, nLay + do i = 1, ix + if (clouds(i,k,1) == 0.0 .and. clouds(i,k-1,1) > 0.0) then + alpha(i,k) = 0.0 + endif + enddo + enddo + endif + !> - Call gethml() to compute low,mid,high,total, and boundary layer !! cloud fractions and clouds top/bottom layer indices for low, mid, !! and high clouds. The three cloud domain boundaries are defined by @@ -1272,6 +1285,19 @@ subroutine progcld2 & alpha(:,:) = 0. endif + ! Revise alpha for exponential-random cloud overlap + ! Decorrelate layers when a clear layer follows a cloudy layer to enforce + ! random correlation between non-adjacent blocks of cloudy layers + if (iovr == 5) then + do k = 2, nLay + do i = 1, ix + if (clouds(i,k,1) == 0.0 .and. clouds(i,k-1,1) > 0.0) then + alpha(i,k) = 0.0 + endif + enddo + enddo + endif + !> - Call gethml() to compute low,mid,high,total, and boundary layer !! cloud fractions and clouds top/bottom layer indices for low, mid, !! and high clouds. @@ -1699,6 +1725,19 @@ subroutine progcld3 & alpha(:,:) = 0. endif + ! Revise alpha for exponential-random cloud overlap + ! Decorrelate layers when a clear layer follows a cloudy layer to enforce + ! random correlation between non-adjacent blocks of cloudy layers + if (iovr == 5) then + do k = 2, nLay + do i = 1, ix + if (clouds(i,k,1) == 0.0 .and. clouds(i,k-1,1) > 0.0) then + alpha(i,k) = 0.0 + endif + enddo + enddo + endif + !> -# Call gethml() to compute low,mid,high,total, and boundary layer !! cloud fractions and clouds top/bottom layer indices for low, mid, !! and high clouds. @@ -2062,6 +2101,19 @@ subroutine progcld4 & alpha(:,:) = 0. endif + ! Revise alpha for exponential-random cloud overlap + ! Decorrelate layers when a clear layer follows a cloudy layer to enforce + ! random correlation between non-adjacent blocks of cloudy layers + if (iovr == 5) then + do k = 2, nLay + do i = 1, ix + if (clouds(i,k,1) == 0.0 .and. clouds(i,k-1,1) > 0.0) then + alpha(i,k) = 0.0 + endif + enddo + enddo + endif + ! --- compute low, mid, high, total, and boundary layer cloud fractions ! and clouds top/bottom layer indices for low, mid, and high clouds. ! The three cloud domain boundaries are defined by ptopc. The cloud @@ -2416,6 +2468,19 @@ subroutine progcld4o & alpha(:,:) = 0. endif + ! Revise alpha for exponential-random cloud overlap + ! Decorrelate layers when a clear layer follows a cloudy layer to enforce + ! random correlation between non-adjacent blocks of cloudy layers + if (iovr == 5) then + do k = 2, nLay + do i = 1, ix + if (clouds(i,k,1) == 0.0 .and. clouds(i,k-1,1) > 0.0) then + alpha(i,k) = 0.0 + endif + enddo + enddo + endif + !> - Call gethml() to compute low, mid, high, total, and boundary layer cloud fractions !! and clouds top/bottom layer indices for low, mid, and high clouds. !! The three cloud domain boundaries are defined by ptopc. The cloud @@ -2792,6 +2857,19 @@ subroutine progcld5 & alpha(:,:) = 0. endif + ! Revise alpha for exponential-random cloud overlap + ! Decorrelate layers when a clear layer follows a cloudy layer to enforce + ! random correlation between non-adjacent blocks of cloudy layers + if (iovr == 5) then + do k = 2, nLay + do i = 1, ix + if (clouds(i,k,1) == 0.0 .and. clouds(i,k-1,1) > 0.0) then + alpha(i,k) = 0.0 + endif + enddo + enddo + endif + !> - Call gethml() to compute low,mid,high,total, and boundary layer !! cloud fractions and clouds top/bottom layer indices for low, mid, !! and high clouds. @@ -3128,8 +3206,7 @@ subroutine progcld6 & enddo enddo - ! What portion of water and ice contents is associated with the - ! partly cloudy boxes + ! What portion of water and ice contents is associated with the partly cloudy boxes do i = 1, IX do k = 1, NLAY-1 if (cldtot(i,k).ge.climit .and. cldtot(i,k).lt.ovcst) then @@ -3188,6 +3265,19 @@ subroutine progcld6 & alpha(:,:) = 0. endif + ! Revise alpha for exponential-random cloud overlap + ! Decorrelate layers when a clear layer follows a cloudy layer to enforce + ! random correlation between non-adjacent blocks of cloudy layers + if (iovr == 5) then + do k = 2, nLay + do i = 1, ix + if (clouds(i,k,1) == 0.0 .and. clouds(i,k-1,1) > 0.0) then + alpha(i,k) = 0.0 + endif + enddo + enddo + endif + !> - Call gethml() to compute low,mid,high,total, and boundary layer !! cloud fractions and clouds top/bottom layer indices for low, mid, !! and high clouds. @@ -3274,7 +3364,7 @@ subroutine progcld_thompson & ! slmsk (IX) : sea/land mask array (sea:0,land:1,sea-ice:2) ! ! dz (ix,nlay) : layer thickness (km) ! ! delp (ix,nlay) : model layer pressure thickness in mb (100Pa) ! -! gridkm (IX) : grid length in km ! +! gridkm : grid length in km ! ! IX : horizontal dimention ! ! NLAY,NLP1 : vertical layer/level dimensions ! ! uni_cld : logical - true for cloud fraction from shoc ! @@ -3333,8 +3423,8 @@ subroutine progcld_thompson & real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, & & slmsk - real(kind=kind_phys), dimension(:), intent(in) :: latdeg, gridkm - real(kind=kind_phys), intent(in) :: julian + real(kind=kind_phys), dimension(:), intent(in) :: latdeg + real(kind=kind_phys), intent(in) :: julian, gridkm integer, intent(in) :: yearlen ! --- outputs @@ -3408,14 +3498,14 @@ subroutine progcld_thompson & enddo !> - Compute cloud liquid/ice condensate path in \f$ g/m^2 \f$ . -!> - Since using Thompson MP, assume 1 percent of snow is actually in +!> - Since using Thompson MP, assume 20 percent of snow is actually in !! ice sizes. do k = 1, NLAY-1 do i = 1, IX cwp(i,k) = max(0.0, clw(i,k,ntcw) * dz(i,k)*1.E6) crp(i,k) = 0.0 - snow_mass_factor = 0.99 + snow_mass_factor = 0.85 cip(i,k) = max(0.0, (clw(i,k,ntiw) & & + (1.0-snow_mass_factor)*clw(i,k,ntsw))*dz(i,k)*1.E6) if (re_snow(i,k) .gt. snow_max_radius)then @@ -3481,7 +3571,7 @@ subroutine progcld_thompson & endif call cal_cldfra3(cldfra1d, qv1d, qc1d, qi1d, qs1d, dz1d, & - & p1d, t1d, xland, gridkm(i), & + & p1d, t1d, xland, gridkm, & & .false., max_relh, 1, nlay, .false.) do k = 1, NLAY @@ -3555,6 +3645,19 @@ subroutine progcld_thompson & alpha(:,:) = 0. endif + ! Revise alpha for exponential-random cloud overlap + ! Decorrelate layers when a clear layer follows a cloudy layer to enforce + ! random correlation between non-adjacent blocks of cloudy layers + if (iovr == 5) then + do k = 2, nLay + do i = 1, ix + if (clouds(i,k,1) == 0.0 .and. clouds(i,k-1,1) > 0.0) then + alpha(i,k) = 0.0 + endif + enddo + enddo + endif + !> - Call gethml() to compute low,mid,high,total, and boundary layer !! cloud fractions and clouds top/bottom layer indices for low, mid, !! and high clouds. @@ -3952,6 +4055,19 @@ subroutine progclduni & alpha(:,:) = 0. endif + ! Revise alpha for exponential-random cloud overlap + ! Decorrelate layers when a clear layer follows a cloudy layer to enforce + ! random correlation between non-adjacent blocks of cloudy layers + if (iovr == 5) then + do k = 2, nLay + do i = 1, ix + if (clouds(i,k,1) == 0.0 .and. clouds(i,k-1,1) > 0.0) then + alpha(i,k) = 0.0 + endif + enddo + enddo + endif + !> - Call gethml() to compute low,mid,high,total, and boundary layer !! cloud fractions and clouds top/bottom layer indices for low, mid, !! and high clouds. @@ -4495,16 +4611,16 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & DO k = kts,kte delz = MAX(100., dz(k)) - RH_00L = 0.77+MIN(0.22,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) - RH_00O = 0.85+MIN(0.14,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) + RH_00L = 0.74+MIN(0.25,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) + RH_00O = 0.82+MIN(0.17,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) RHUM = rh(k) - if (qc(k).ge.1.E-6 .or. qi(k).ge.1.E-7 & - & .or. (qs(k).gt.1.E-6 .and. t(k).lt.273.)) then + if (qc(k).ge.1.E-5 .or. qi(k).ge.1.E-5 & + & .or. (qs(k).gt.1.E-5 .and. t(k).lt.273.)) then CLDFRA(K) = 1.0 elseif (((qc(k)+qi(k)).gt.1.E-10) .and. & - & ((qc(k)+qi(k)).lt.1.E-6)) then - CLDFRA(K) = MIN(0.99, 0.1*(11.0 + log10(qc(k)+qi(k)))) + & ((qc(k)+qi(k)).lt.1.E-5)) then + CLDFRA(K) = MIN(0.99, 0.20*(10.0 + log10(qc(k)+qi(k)))) else IF ((XLAND-1.5).GT.0.) THEN !--- Ocean @@ -4513,7 +4629,7 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & RH_00 = RH_00L ENDIF - tc = MAX(-80.0, t(k) - 273.15) + tc = t(k) - 273.15 if (tc .lt. -12.0) RH_00 = RH_00L if (tc .gt. 20.0) then @@ -4525,12 +4641,12 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & if (max_relh.gt.1.12 .or. (.NOT.(modify_qvapor)) ) then !..For HRRR model, the following look OK. RHUM = MIN(rh(k), 1.45) - RH_00 = RH_00 + (1.45-RH_00)*(-12.0-tc)/(-12.0+85.) + RH_00 = RH_00 + (1.45-RH_00)*(-12.0-tc)/(-12.0+112.) CLDFRA(K) = MAX(0.,1.0-SQRT((1.46-RHUM)/(1.46-RH_00))) else !..but for the GFS model, RH is way lower. RHUM = MIN(rh(k), 1.05) - RH_00 = RH_00 + (1.05-RH_00)*(-12.0-tc)/(-12.0+85.) + RH_00 = RH_00 + (1.05-RH_00)*(-12.0-tc)/(-12.0+112.) CLDFRA(K) = MAX(0.,1.0-SQRT((1.06-RHUM)/(1.06-RH_00))) endif endif @@ -4548,6 +4664,15 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & call adjust_cloudFinal(cldfra, qc, qi, rhoa, dz, kts,kte) + if (debug_flag .and. ndebug.lt.25) then + do k = kts,kte + write(6,'(a,i3,f9.2,f7.1,f7.2,f6.1,f6.3,f12.7,f12.7,f12.7)') & + & ' DEBUG-GT: ', k, p(k)*0.01, dz(k), t(k)-273.15, & + & rh(k)*100., cldfra(k), qc(k)*1.E3, qi(k)*1.E3, qs(k)*1.E3 + enddo + ndebug = ndebug + 1 + endif + !..Intended for cold start model runs, we use modify_qvapor to ensure that cloudy !.. areas are actually saturated such that the inserted clouds do not evaporate a !.. timestep later. @@ -4689,9 +4814,9 @@ SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, Dz1d, entrmnt,& k = k - 1 ENDDO - k_cldb = k_m12C + 3 + k_cldb = k_m12C + 5 in_cloud = .false. - k = k_m12C + 2 + k = k_m12C + 4 DO WHILE (.not. in_cloud .AND. k.gt.kbot) k_cldt = 0 if (cfr1d(k).ge.0.01) then @@ -4740,13 +4865,12 @@ SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs,T,dz,entr, k1,k2,kts,kte) do k = k1, k2 tdz = tdz + dz(k) enddo -! max_iwc = ABS(qvs(k2)-qvs(k1)) - max_iwc = MAX(0.0, qvs(k1)-qvs(k2)) + max_iwc = ABS(qvs(k2)-qvs(k1)) do k = k1, k2 - max_iwc = MAX(1.E-6, max_iwc - (qi(k)+qs(k))) + max_iwc = MAX(1.E-5, max_iwc - (qi(k)+qs(k))) enddo - max_iwc = MIN(1.E-4, max_iwc) + max_iwc = MIN(2.E-3, max_iwc) this_dz = 0.0 do k = k1, k2 @@ -4756,7 +4880,7 @@ SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs,T,dz,entr, k1,k2,kts,kte) this_dz = this_dz + dz(k) endif this_iwc = max_iwc*this_dz/tdz - iwc = MAX(1.E-6, this_iwc*(1.-entr)) + iwc = MAX(5.E-6, this_iwc*(1.-entr)) if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.203.16) then qi(k) = qi(k) + cfr(k)*cfr(k)*iwc endif @@ -4781,14 +4905,13 @@ SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte) do k = k1, k2 tdz = tdz + dz(k) enddo -! max_lwc = ABS(qvs(k2)-qvs(k1)) - max_lwc = MAX(0.0, qvs(k1)-qvs(k2)) + max_lwc = ABS(qvs(k2)-qvs(k1)) ! print*, ' max_lwc = ', max_lwc, ' over DZ=',tdz do k = k1, k2 - max_lwc = MAX(1.E-6, max_lwc - qc(k)) + max_lwc = MAX(1.E-5, max_lwc - qc(k)) enddo - max_lwc = MIN(1.E-4, max_lwc) + max_lwc = MIN(2.E-3, max_lwc) this_dz = 0.0 do k = k1, k2 @@ -4798,8 +4921,8 @@ SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte) this_dz = this_dz + dz(k) endif this_lwc = max_lwc*this_dz/tdz - lwc = MAX(1.E-6, this_lwc*(1.-entr)) - if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.258.16) then + lwc = MAX(5.E-6, this_lwc*(1.-entr)) + if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.253.16) then qc(k) = qc(k) + cfr(k)*cfr(k)*lwc endif enddo @@ -4851,6 +4974,6 @@ SUBROUTINE adjust_cloudFinal(cfr, qc, qi, Rho,dz, kts,kte) END SUBROUTINE adjust_cloudFinal !........................................! - end module module_radiation_clouds + end module module_radiation_clouds ! !! @} !========================================! diff --git a/physics/radlw_main.F90 b/physics/radlw_main.F90 index 14a28cf6b..95bc0b059 100644 --- a/physics/radlw_main.F90 +++ b/physics/radlw_main.F90 @@ -1363,7 +1363,8 @@ subroutine rlwinit & ! =1: maximum/random overlapping clouds ! ! =2: maximum overlap cloud (isubcol>0 only) ! ! =3: decorrelation-length overlap (for isubclw>0 only) ! -! =4: exponential overlap cloud +! =4: exponential cloud overlap (AER) ! +! =5: exponential-random cloud overlap (AER) ! ! ! ! ******************************************************************* ! ! original code description ! @@ -1407,7 +1408,7 @@ subroutine rlwinit & ! !===> ... begin here ! - if ( iovr<0 .or. iovr>4 ) then + if ( iovr<0 .or. iovr>5 ) then print *,' *** Error in specification of cloud overlap flag', & & ' IOVR=',iovr,' in RLWINIT !!' stop @@ -1896,6 +1897,7 @@ subroutine mcica_subcol & ! other control flags from module variables: ! ! iovr : control flag for cloud overlapping method ! ! =0:random; =1:maximum/random: =2:maximum; =3:decorr ! +! =4:exponential; =5:exponential-random ! ! ! ! ===================== end of definitions ==================== ! @@ -2084,39 +2086,39 @@ subroutine mcica_subcol & ! --- setup 2 sets of random numbers -! call random_number ( rand2d, stat ) + call random_number ( rand2d, stat ) -! k1 = 0 -! do k = 1, nlay -! do n = 1, ngptlw -! k1 = k1 + 1 -! cdfunc(n,k) = rand2d(k1) -! enddo -! enddo + k1 = 0 + do k = 1, nlay + do n = 1, ngptlw + k1 = k1 + 1 + cdfunc(n,k) = rand2d(k1) + enddo + enddo -! call random_number ( rand2d, stat ) + call random_number ( rand2d, stat ) -! k1 = 0 -! do k = 1, nlay -! do n = 1, ngptlw -! k1 = k1 + 1 -! cdfun2(n,k) = rand2d(k1) -! enddo -! enddo + k1 = 0 + do k = 1, nlay + do n = 1, ngptlw + k1 = k1 + 1 + cdfun2(n,k) = rand2d(k1) + enddo + enddo ! --- then working upward from the surface: ! if a random number (from an independent set: cdfun2) is smaller than ! alpha, then use the previous layer's number, otherwise use a new random ! number (keep the originally assigned one in cdfunc for that layer). -! do k = 2, nlay -! k1 = k - 1 -! do n = 1, ngptlw -! if ( cdfun2(n,k) < alpha(k) ) then -! cdfunc(n,k) = cdfunc(n,k1) -! endif -! enddo -! enddo + do k = 2, nlay + k1 = k - 1 + do n = 1, ngptlw + if ( cdfun2(n,k) < alpha(k) ) then + cdfunc(n,k) = cdfunc(n,k1) + endif + enddo + enddo end select diff --git a/physics/radsw_main.F90 b/physics/radsw_main.F90 index 44de9848c..d09f586a3 100644 --- a/physics/radsw_main.F90 +++ b/physics/radsw_main.F90 @@ -1120,7 +1120,7 @@ subroutine rrtmg_sw_run & endif enddo zcf0 = zcf0 * zcf1 - else if (iovr >= 2 .and. iovr /= 4) then + else if (iovr >= 2) then do k = 1, nlay zcf0 = min ( zcf0, f_one-cfrac(k) ) ! used only as clear/cloudy indicator enddo @@ -1436,6 +1436,8 @@ subroutine rswinit & ! =1: maximum/random overlapping clouds ! ! =2: maximum overlap cloud ! ! =3: decorrelation-length overlap clouds ! +! =4: exponential cloud overlap (AER) ! +! =5: exponential-random cloud overlap (AER) ! ! iswmode - control flag for 2-stream transfer scheme ! ! =1; delta-eddington (joseph et al., 1976) ! ! =2: pifm (zdunkowski et al., 1980) ! @@ -1467,7 +1469,7 @@ subroutine rswinit & ! !===> ... begin here ! - if ( iovr<0 .or. iovr>4 ) then + if ( iovr<0 .or. iovr>5 ) then print *,' *** Error in specification of cloud overlap flag', & & ' IOVR=',iovr,' in RSWINIT !!' stop @@ -1935,7 +1937,7 @@ subroutine cldprop & !> -# if physparam::isubcsw > 0, call mcica_subcol() to distribute !! cloud properties to each g-point. - if ( isubcsw > 0 .and. iovr /= 4 ) then ! mcica sub-col clouds approx + if ( isubcsw > 0 ) then ! mcica sub-col clouds approx cldf(:) = cfrac(:) where (cldf(:) < ftiny)