diff --git a/CMakeLists.txt b/CMakeLists.txt index 4d5c8eae4..e9689636d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -195,6 +195,13 @@ elseif (${CMAKE_Fortran_COMPILER_ID} STREQUAL "Intel") list(APPEND SCHEMES_SFX_OPT ${LOCAL_CURRENT_SOURCE_DIR}/physics/module_sf_mynn.F90) endif() + if (${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rrtmgp/kernels/mo_gas_optics_kernels.F90 IN_LIST SCHEMES) + # Reduce optimization for mo_gas_optics_kernels.F90 (to avoid an apparent compiler bug with Intel 19+) + SET_SOURCE_FILES_PROPERTIES(${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rrtmgp/kernels/mo_gas_optics_kernels.F90 + PROPERTIES COMPILE_FLAGS "${CMAKE_Fortran_FLAGS_OPT} -O1") + list(APPEND SCHEMES_SFX_OPT ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rrtmgp/kernels/mo_gas_optics_kernels.F90) + endif() + # Remove files with special compiler flags from list of files with standard compiler flags if (SCHEMES_SFX_OPT) list(REMOVE_ITEM SCHEMES ${SCHEMES_SFX_OPT}) diff --git a/CODEOWNERS b/CODEOWNERS index c163e7b80..5080ce409 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -3,7 +3,7 @@ # These owners will be the default owners for everything in the repo. #* @defunkt -* @climbfuji @llpcarson @grantfirl @mzhangw @panll @mkavulich +* @climbfuji @grantfirl @mzhangw @panll @mkavulich # Order is important. The last matching pattern has the most precedence. # So if a pull request only touches javascript files, only these owners diff --git a/physics/GFS_rrtmg_post.F90 b/physics/GFS_rrtmg_post.F90 index b882930bf..8584f8463 100644 --- a/physics/GFS_rrtmg_post.F90 +++ b/physics/GFS_rrtmg_post.F90 @@ -17,7 +17,7 @@ subroutine GFS_rrtmg_post_run (im, km, kmp1, lm, ltp, kt, kb, kd, nspc1, & nfxr, nday, lsswr, lslwr, lssav, fhlwr, fhswr, raddt, coszen, & coszdg, prsi, tgrs, aerodp, cldsa, mtopa, mbota, clouds1, & cldtaulw, cldtausw, sfcflw, sfcfsw, topflw, topfsw, scmpsw, & - fluxr, errmsg, errflg) + fluxr, total_albedo, errmsg, errflg) use machine, only: kind_phys use module_radsw_parameters, only: topfsw_type, sfcfsw_type, & @@ -43,6 +43,7 @@ subroutine GFS_rrtmg_post_run (im, km, kmp1, lm, ltp, kt, kb, kd, nspc1, & real(kind=kind_phys), dimension(im,lm+LTP), intent(in) :: clouds1 real(kind=kind_phys), dimension(im,lm+LTP), intent(in) :: cldtausw real(kind=kind_phys), dimension(im,lm+LTP), intent(in) :: cldtaulw + real(kind=kind_phys), dimension(im), intent(inout) :: total_albedo type(sfcflw_type), dimension(im), intent(in) :: sfcflw type(sfcfsw_type), dimension(im), intent(in) :: sfcfsw @@ -196,6 +197,12 @@ subroutine GFS_rrtmg_post_run (im, km, kmp1, lm, ltp, kt, kb, kd, nspc1, & endif endif ! end_if_lssav + +! --- The total sky (with clouds) shortwave albedo + total_albedo = 0.0 + if (lsswr) then + where(topfsw(:)%dnfxc>0) total_albedo(:) = topfsw(:)%upfxc/topfsw(:)%dnfxc + endif ! end subroutine GFS_rrtmg_post_run diff --git a/physics/GFS_rrtmg_post.meta b/physics/GFS_rrtmg_post.meta index 8677d9900..80bd5c22c 100644 --- a/physics/GFS_rrtmg_post.meta +++ b/physics/GFS_rrtmg_post.meta @@ -195,7 +195,7 @@ standard_name = total_cloud_fraction long_name = layer total cloud fraction units = frac - dimensions = (horizontal_loop_extent,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) type = real kind = kind_phys intent = in @@ -258,6 +258,14 @@ type = real kind = kind_phys intent = inout +[total_albedo] + standard_name = total_sky_albedo + long_name = total sky albedo at toa + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index dbea66985..973ac02fd 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -34,8 +34,9 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & plvl, plyr, tlvl, tlyr, qlyr, olyr, gasvmr_co2, gasvmr_n2o, gasvmr_ch4,& gasvmr_o2, gasvmr_co, gasvmr_cfc11, gasvmr_cfc12, gasvmr_cfc22, & gasvmr_ccl4, gasvmr_cfc113, aerodp, clouds6, clouds7, clouds8, & - clouds9, cldsa, cldfra, faersw1, faersw2, faersw3, faerlw1, faerlw2, & - faerlw3, alpha, errmsg, errflg) + clouds9, cldsa, cldfra, cldfra2d, lwp_ex,iwp_ex, lwp_fc,iwp_fc, & + faersw1, faersw2, faersw3, faerlw1, faerlw2, faerlw3, alpha, & + errmsg, errflg) use machine, only: kind_phys @@ -54,6 +55,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & & progcld2, & & progcld4, progcld5, & & progcld6, & + & progcld_thompson, & & progclduni, & & cal_cldfra3, & & find_cloudLayers, & @@ -125,6 +127,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & real(kind=kind_phys), dimension(:,:), intent(inout) :: clouds1, & clouds2, clouds3, & clouds4, clouds5 + real(kind=kind_phys), dimension(:), intent(out) :: lwp_ex,iwp_ex, & + & lwp_fc,iwp_fc integer, intent(out) :: kd, kt, kb @@ -158,6 +162,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & clouds8, & clouds9, & cldfra + real(kind=kind_phys), dimension(:), intent(out) :: cldfra2d real(kind=kind_phys), dimension(:,:), intent(out) :: cldsa real(kind=kind_phys), dimension(:,:,:), intent(out) :: faersw1,& @@ -191,9 +196,10 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & real(kind=kind_phys), dimension(im,lm+LTP) :: & re_cloud, re_ice, re_snow, qv_mp, qc_mp, & qi_mp, qs_mp, nc_mp, ni_mp, nwfa + real (kind=kind_phys), dimension(lm) :: cldfra1d, qv1d, & + & qc1d, qi1d, qs1d, dz1d, p1d, t1d ! for F-A MP - real(kind=kind_phys), dimension(im,lm+LTP) :: qc_save, qi_save, qs_save real(kind=kind_phys), dimension(im,lm+LTP+1) :: tem2db, hz real(kind=kind_phys), dimension(im,lm+LTP,min(4,ncnd)) :: ccnd @@ -206,6 +212,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & ! for stochastic cloud perturbations real(kind=kind_phys), dimension(im) :: cldp1d real (kind=kind_phys) :: alpha0,beta0,m,s,cldtmp,tmp_wt,cdfz + real (kind=kind_phys) :: max_relh integer :: iflag integer :: ids, ide, jds, jde, kds, kde, & @@ -228,6 +235,21 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & LP1 = LM + 1 ! num of in/out levels + gridkm = sqrt(2.0)*sqrt(dx(1)*0.001*dx(1)*0.001) + + if (imp_physics == imp_physics_thompson) then + max_relh = 1.5 + else + max_relh = 1.1 + endif + + do i = 1, IM + lwp_ex(i) = 0.0 + iwp_ex(i) = 0.0 + lwp_fc(i) = 0.0 + iwp_fc(i) = 0.0 + enddo + ! --- ... set local /level/layer indexes corresponding to in/out ! variables @@ -854,88 +876,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & enddo endif - !mz HWRF physics: icloud=3 - if(icloud == 3) then - - ! Set internal dimensions - ids = 1 - ims = 1 - its = 1 - ide = size(xlon,1) - ime = size(xlon,1) - ite = size(xlon,1) - jds = 1 - jms = 1 - jts = 1 - jde = 1 - jme = 1 - jte = 1 - kds = 1 - kms = 1 - kts = 1 - kde = lm+LTP ! should this be lmk instead of lm? no, or? - kme = lm+LTP - kte = lm+LTP - - do k = 1, LMK - do i = 1, IM - rho(i,k)=plyr(i,k)*100./(con_rd*tlyr(i,k)) - plyrpa(i,k)=plyr(i,k)*100. !hPa->Pa - end do - end do - - do i=1,im - if (slmsk(i)==1. .or. slmsk(i)==2.) then ! sea/land/ice mask (=0/1/2) in FV3 - xland(i)=1.0 ! but land/water = (1/2) in HWRF - else - xland(i)=2.0 - endif - enddo - - gridkm = sqrt(2.0)*sqrt(dx(1)*0.001*dx(1)*0.001) - - do i =1, im - do k =1, lmk - qc_save(i,k) = ccnd(i,k,1) - qi_save(i,k) = ccnd(i,k,2) - qs_save(i,k) = ccnd(i,k,4) - enddo - enddo - - - call cal_cldfra3(cldcov,qlyr,ccnd(:,:,1),ccnd(:,:,2), & - ccnd(:,:,4),plyrpa,tlyr,rho,xland,gridkm, & - ids,ide,jds,jde,kds,kde, & - ims,ime,jms,jme,kms,kme, & - its,ite,jts,jte,kts,kte) - - !mz* back to micro-only qc qi,qs - do i =1, im - do k =1, lmk - ccnd(i,k,1) = qc_save(i,k) - ccnd(i,k,2) = qi_save(i,k) - ccnd(i,k,4) = qs_save(i,k) - enddo - enddo - - endif ! icloud == 3 - - if (lextop) then - do i=1,im - cldcov(i,lyb) = cldcov(i,lya) - deltaq(i,lyb) = deltaq(i,lya) - cnvw (i,lyb) = cnvw (i,lya) - cnvc (i,lyb) = cnvc (i,lya) - enddo - if (effr_in) then - do i=1,im - effrl(i,lyb) = effrl(i,lya) - effri(i,lyb) = effri(i,lya) - effrr(i,lyb) = effrr(i,lya) - effrs(i,lyb) = effrs(i,lya) - enddo - endif - endif if (imp_physics == imp_physics_zhao_carr) then ccnd(1:IM,1:LMK,1) = ccnd(1:IM,1:LMK,1) + cnvw(1:IM,1:LMK) @@ -1012,6 +952,20 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & elseif(imp_physics == imp_physics_thompson) then ! Thompson MP if(do_mynnedmf .or. imfdeepcnv == imfdeepcnv_gf ) then ! MYNN PBL or GF conv + + if (icloud == 3) then + call progcld_thompson (plyr,plvl,tlyr,qlyr,qstl,rhly, & ! --- inputs + tracer1,xlat,xlon,slmsk,dz,delp, & + ntrac-1, ntcw-1,ntiw-1,ntrw-1, & + ntsw-1,ntgl-1, & + im, lm, lmp, uni_cld, lmfshal, lmfdeep2, & + cldcov(:,1:LM), effrl_inout, & + effri_inout, effrs_inout, & + lwp_ex, iwp_ex, lwp_fc, iwp_fc, & + dzb, xlat_d, julian, yearlen, gridkm, & + clouds, cldsa, mtopa ,mbota, de_lgth, alpha) ! --- outputs + else + !-- MYNN PBL or convective GF !-- use cloud fractions with SGS clouds do k=1,lmk @@ -1028,18 +982,35 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & effrl, effri, effrr, effrs, effr_in , & dzb, xlat_d, julian, yearlen, & clouds, cldsa, mtopa, mbota, de_lgth, alpha) ! --- outputs + endif else ! MYNN PBL or GF convective are not used - call progcld6 (plyr,plvl,tlyr,qlyr,qstl,rhly,tracer1, & ! --- inputs - xlat,xlon,slmsk,dz,delp, & + + if (icloud == 3) then + call progcld_thompson (plyr,plvl,tlyr,qlyr,qstl,rhly, & ! --- inputs + tracer1,xlat,xlon,slmsk,dz,delp, & + ntrac-1, ntcw-1,ntiw-1,ntrw-1, & + ntsw-1,ntgl-1, & + im, lm, lmp, uni_cld, lmfshal, lmfdeep2, & + cldcov(:,1:LM), effrl_inout, & + effri_inout, effrs_inout, & + lwp_ex, iwp_ex, lwp_fc, iwp_fc, & + dzb, xlat_d, julian, yearlen, gridkm, & + clouds, cldsa, mtopa ,mbota, de_lgth, alpha) ! --- outputs + + else + call progcld6 (plyr,plvl,tlyr,qlyr,qstl,rhly, & ! --- inputs + tracer1,xlat,xlon,slmsk,dz,delp, & ntrac-1, ntcw-1,ntiw-1,ntrw-1, & ntsw-1,ntgl-1, & im, lmk, lmp, uni_cld, lmfshal, lmfdeep2, & - cldcov(:,1:LMK), effrl_inout(:,:), & - effri_inout(:,:), effrs_inout(:,:), & + cldcov(:,1:LMK), effrl_inout, & + effri_inout, effrs_inout, & + lwp_ex, iwp_ex, lwp_fc, iwp_fc, & dzb, xlat_d, julian, yearlen, & clouds, cldsa, mtopa ,mbota, de_lgth, alpha) ! --- outputs + endif endif ! MYNN PBL or GF endif ! end if_imp_physics @@ -1071,7 +1042,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & enddo ! end_do_i_loop enddo ! end_do_k_loop endif - do k = 1, LMK + do k = 1, LM do i = 1, IM clouds1(i,k) = clouds(i,k,1) clouds2(i,k) = clouds(i,k,2) @@ -1085,6 +1056,12 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & cldfra(i,k) = clouds(i,k,1) enddo enddo + do i = 1, IM + cldfra2d(i) = 0.0 + do k = 1, LM-1 + cldfra2d(i) = max(cldfra2d(i), cldfra(i,k)) + enddo + enddo ! mg, sfc-perts ! --- scale random patterns for surface perturbations with diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index 09ed62f7c..85f29dc9c 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -570,7 +570,7 @@ intent = in [sppt_amp] standard_name = total_amplitude_of_sppt_perturbation - long_name = toal ampltidue of stochastic sppt perturbation + long_name = total ampltidue of stochastic sppt perturbation units = none dimensions = () type = real @@ -978,6 +978,46 @@ type = real kind = kind_phys intent = out +[cldfra2d] + standard_name = max_in_column_cloud_fraction + long_name = instantaneous 2D (max-in-column) cloud fraction + units = frac + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[lwp_ex] + standard_name = liq_water_path_from_microphysics + long_name = total liquid water path from explicit microphysics + units = kg m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[iwp_ex] + standard_name = ice_water_path_from_microphysics + long_name = total ice water path from explicit microphysics + units = kg m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[lwp_fc] + standard_name = liq_water_path_from_cloud_fraction + long_name = total liquid water path from cloud fraction scheme + units = kg m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out +[iwp_fc] + standard_name = ice_water_path_from_cloud_fraction + long_name = total ice water path from cloud fraction scheme + units = kg m-2 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = out [faersw1] standard_name = aerosol_optical_depth_for_shortwave_bands_01_16 long_name = aerosol optical depth for shortwave bands 01-16 diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index dbf554508..3183ca4bf 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -149,7 +149,7 @@ MODULE module_mp_thompson REAL, PARAMETER, PRIVATE:: fv_s = 100.0 REAL, PARAMETER, PRIVATE:: av_g = 442.0 REAL, PARAMETER, PRIVATE:: bv_g = 0.89 - REAL, PARAMETER, PRIVATE:: av_i = 1847.5 + REAL, PARAMETER, PRIVATE:: av_i = 1493.9 REAL, PARAMETER, PRIVATE:: bv_i = 1.0 REAL, PARAMETER, PRIVATE:: av_c = 0.316946E8 REAL, PARAMETER, PRIVATE:: bv_c = 2.0 @@ -214,8 +214,8 @@ MODULE module_mp_thompson REAL, PARAMETER, PRIVATE:: xm0i = 1.E-12 REAL, PARAMETER, PRIVATE:: D0c = 1.E-6 REAL, PARAMETER, PRIVATE:: D0r = 50.E-6 - REAL, PARAMETER, PRIVATE:: D0s = 200.E-6 - REAL, PARAMETER, PRIVATE:: D0g = 250.E-6 + REAL, PARAMETER, PRIVATE:: D0s = 300.E-6 + REAL, PARAMETER, PRIVATE:: D0g = 350.E-6 REAL, PRIVATE:: D0i, xm0s, xm0g !..Min and max radiative effective radius of cloud water, cloud ice, and snow; @@ -2188,7 +2188,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & ni(k) = MAX(R2, ni1d(k)*rho(k)) if (ni(k).le. R2) then lami = cie(2)/5.E-6 - ni(k) = MIN(9999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) + ni(k) = MIN(999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) endif L_qi(k) = .true. lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi @@ -2196,7 +2196,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & xDi = (bm_i + mu_i + 1.) * ilami if (xDi.lt. 5.E-6) then lami = cie(2)/5.E-6 - ni(k) = MIN(9999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) + ni(k) = MIN(999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) elseif (xDi.gt. 300.E-6) then lami = cie(2)/300.E-6 ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i @@ -2467,7 +2467,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & tau = 3.72/(rc(k)*taud) prr_wau(k) = zeta/tau prr_wau(k) = MIN(DBLE(rc(k)*odts), prr_wau(k)) - pnr_wau(k) = prr_wau(k) / (am_r*nu_c*200.*D0r*D0r*D0r) ! RAIN2M + pnr_wau(k) = prr_wau(k) / (am_r*nu_c*10.*D0r*D0r*D0r) ! RAIN2M pnc_wau(k) = MIN(DBLE(nc(k)*odts), prr_wau(k) & / (am_r*mvd_c(k)*mvd_c(k)*mvd_c(k))) ! Qc2M endif @@ -3237,7 +3237,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & xDi = (bm_i + mu_i + 1.) * ilami if (xDi.lt. 5.E-6) then lami = cie(2)/5.E-6 - xni = MIN(9999.D3, cig(1)*oig2*xri/am_i*lami**bm_i) + xni = MIN(999.D3, cig(1)*oig2*xri/am_i*lami**bm_i) niten(k) = (xni-ni1d(k)*rho(k))*odts*orho elseif (xDi.gt. 300.E-6) then lami = cie(2)/300.E-6 @@ -3248,8 +3248,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & niten(k) = -ni1d(k)*odts endif xni=MAX(0.,(ni1d(k) + niten(k)*dtsave)*rho(k)) - if (xni.gt.9999.E3) & - niten(k) = (9999.E3-ni1d(k)*rho(k))*odts*orho + if (xni.gt.999.E3) & + niten(k) = (999.E3-ni1d(k)*rho(k))*odts*orho !> - Rain tendency qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) & @@ -3835,7 +3835,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & t3_vts = Kap0*csg(1)*ils1**cse(1) t4_vts = Kap1*Mrat**mu_s*csg(7)*ils2**cse(7) vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts) - if (temp(k).gt. (T_0+0.1)) then + if (prr_sml(k) .gt. 0.0) then ! vtsk(k) = MAX(vts*vts_boost(k), & ! & vts*((vtrk(k)-vts*vts_boost(k))/(temp(k)-T_0))) SR = rs(k)/(rs(k)+rr(k)) @@ -4187,7 +4187,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & lami = cie(2)/300.E-6 endif ni1d(k) = MIN(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, & - 9999.D3/rho(k)) + 999.D3/rho(k)) endif qr1d(k) = qr1d(k) + qrten(k)*DT nr1d(k) = MAX(R2/rho(k), nr1d(k) + nrten(k)*DT) diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index dacf6e38e..f58ec8d11 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -59,7 +59,7 @@ ! ! ! 'progcld4o' --- inactive ! ! ! -! 'progcld5' --- thompson/wsm6 cloud microphysics ! +! 'progcld5' --- wsm6 cloud microphysics ! ! inputs: ! ! (plyr,plvl,tlyr,qlyr,qstl,rhly,clw, ! ! xlat,xlon,slmsk, dz, delp, ! @@ -258,8 +258,28 @@ module module_radiation_clouds integer :: llyr = 2 !< upper limit of boundary layer clouds + ! Default ice crystal sizes vs. temperature following Kristjansson and Mitchell + real (kind=kind_phys), dimension(95), parameter :: retab =(/ & + & 5.92779, 6.26422, 6.61973, 6.99539, 7.39234, & + & 7.81177, 8.25496, 8.72323, 9.21800, 9.74075, 10.2930, & + & 10.8765, 11.4929, 12.1440, 12.8317, 13.5581, 14.2319, & + & 15.0351, 15.8799, 16.7674, 17.6986, 18.6744, 19.6955, & + & 20.7623, 21.8757, 23.0364, 24.2452, 25.5034, 26.8125, & + & 27.7895, 28.6450, 29.4167, 30.1088, 30.7306, 31.2943, & + & 31.8151, 32.3077, 32.7870, 33.2657, 33.7540, 34.2601, & + & 34.7892, 35.3442, 35.9255, 36.5316, 37.1602, 37.8078, & + & 38.4720, 39.1508, 39.8442, 40.5552, 41.2912, 42.0635, & + & 42.8876, 43.7863, 44.7853, 45.9170, 47.2165, 48.7221, & + & 50.4710, 52.4980, 54.8315, 57.4898, 60.4785, 63.7898, & + & 65.5604, 71.2885, 75.4113, 79.7368, 84.2351, 88.8833, & + & 93.6658, 98.5739, 103.603, 108.752, 114.025, 119.424, & + & 124.954, 130.630, 136.457, 142.446, 148.608, 154.956, & + & 161.503, 168.262, 175.248, 182.473, 189.952, 197.699, & + & 205.728, 214.055, 222.694, 231.661, 240.971, 250.639/) + public progcld1, progcld2, progcld3, progcld4, progclduni, & - & cld_init, progcld5, progcld6, progcld4o, cal_cldfra3, & + & cld_init, progcld5, progcld4o, & + & progcld6, progcld_thompson, cal_cldfra3, & & find_cloudLayers, adjust_cloudIce, adjust_cloudH2O, & & adjust_cloudFinal, gethml @@ -934,7 +954,7 @@ subroutine progcld2 & ! ================= subprogram documentation block ================ ! ! ! ! subprogram: progcld2 computes cloud related quantities using ! -! Thompson/WSM6 cloud microphysics scheme. ! +! WSM6 cloud microphysics scheme. ! ! ! ! abstract: this program computes cloud fractions from cloud ! ! condensates, ! @@ -2644,8 +2664,8 @@ subroutine progcld5 & enddo !mz* if (uni_cld) then ! use unified sgs clouds generated outside -!mz* use unified sgs or thompson clouds generated outside - if (uni_cld .or. icloud == 3) then +!mz* use unified sgs clouds generated outside + if (uni_cld) then do k = 1, NLAY do i = 1, IX cldtot(i,k) = cldcov(i,k) @@ -2797,7 +2817,7 @@ subroutine progcld5 & clouds(i,k,3) = rew(i,k) clouds(i,k,4) = cip(i,k) clouds(i,k,5) = rei(i,k) - clouds(i,k,6) = crp(i,k) ! added for Thompson + clouds(i,k,6) = crp(i,k) clouds(i,k,7) = rer(i,k) !mz inflg .ne.5 clouds(i,k,8) = 0. @@ -2863,6 +2883,7 @@ subroutine progcld6 & & IX, NLAY, NLP1, & & uni_cld, lmfshal, lmfdeep2, cldcov, & & re_cloud,re_ice,re_snow, & + & lwp_ex, iwp_ex, lwp_fc, iwp_fc, & & dzlay, latdeg, julian, yearlen, & & clouds,clds,mtop,mbot,de_lgth,alpha & ! --- outputs: & ) @@ -2956,6 +2977,8 @@ subroutine progcld6 & real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, & & tlyr, qlyr, qstl, rhly, cldcov, delp, dz, dzlay, & & re_cloud, re_ice, re_snow + real (kind=kind_phys), dimension(:), intent(inout) :: & + & lwp_ex, iwp_ex, lwp_fc, iwp_fc real (kind=kind_phys), dimension(:,:,:), intent(in) :: clw @@ -3063,7 +3086,7 @@ subroutine progcld6 & !> - Compute cloud liquid/ice condensate path in \f$ g/m^2 \f$ . - do k = 1, NLAY + do k = 1, NLAY-1 do i = 1, IX cwp(i,k) = max(0.0, clw(i,k,ntcw) * gfac * delp(i,k)) cip(i,k) = max(0.0, clw(i,k,ntiw) * gfac * delp(i,k)) @@ -3073,6 +3096,21 @@ subroutine progcld6 & enddo enddo +!> - Sum the liquid water and ice paths that come from explicit micro + + do i = 1, IX + lwp_ex(i) = 0.0 + iwp_ex(i) = 0.0 + lwp_fc(i) = 0.0 + iwp_fc(i) = 0.0 + do k = 1, NLAY-1 + lwp_ex(i) = lwp_ex(i) + cwp(i,k) + iwp_ex(i) = iwp_ex(i) + cip(i,k) + csp(i,k) + enddo + lwp_ex(i) = lwp_ex(i)*1.E-3 + iwp_ex(i) = iwp_ex(i)*1.E-3 + enddo + if (uni_cld) then ! use unified sgs clouds generated outside do k = 1, NLAY do i = 1, IX @@ -3085,54 +3123,32 @@ subroutine progcld6 & !> - Calculate layer cloud fraction. clwmin = 0.0 - if (.not. lmfshal) then - do k = 1, NLAY - do i = 1, IX - clwt = 1.0e-6 * (plyr(i,k)*0.001) -! clwt = 2.0e-6 * (plyr(i,k)*0.001) - - if (clwf(i,k) > clwt) then + do k = 1, NLAY-1 + do i = 1, IX + clwt = 1.0e-6 * (plyr(i,k)*0.001) - onemrh= max( 1.e-10, 1.0-rhly(i,k) ) - clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) + if (clwf(i,k) > clwt) then + onemrh= max( 1.e-10, 1.0-rhly(i,k) ) + clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) + if (.not. lmfshal) then tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0) tem1 = 2000.0 / tem1 - -! tem1 = 1000.0 / tem1 - - value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) - tem2 = sqrt( sqrt(rhly(i,k)) ) - - cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 ) - endif - enddo - enddo - else - do k = 1, NLAY - do i = 1, IX - clwt = 1.0e-6 * (plyr(i,k)*0.001) -! clwt = 2.0e-6 * (plyr(i,k)*0.001) - - if (clwf(i,k) > clwt) then - onemrh= max( 1.e-10, 1.0-rhly(i,k) ) - clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) -! + else tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) !jhan if (lmfdeep2) then tem1 = xrc3 / tem1 else tem1 = 100.0 / tem1 endif -! - value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) - tem2 = sqrt( sqrt(rhly(i,k)) ) - - cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 ) endif - enddo - enddo - endif + + value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) + tem2 = sqrt( sqrt(rhly(i,k)) ) + cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 ) + endif + enddo + enddo endif ! if (uni_cld) then @@ -3148,6 +3164,18 @@ subroutine progcld6 & enddo enddo + ! 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 + lwp_fc(i) = lwp_fc(i) + cwp(i,k) + iwp_fc(i) = iwp_fc(i) + cip(i,k) + csp(i,k) + endif + enddo + lwp_fc(i) = lwp_fc(i)*1.E-3 + iwp_fc(i) = iwp_fc(i)*1.E-3 + enddo + if ( lcnorm ) then do k = 1, NLAY do i = 1, IX @@ -3222,6 +3250,374 @@ end subroutine progcld6 !mz + +! This subroutine added by G. Thompson specifically to account for +! explicit (microphysics-produced) cloud liquid water, cloud ice, and +! snow with 100% cloud fraction. Also, a parameterization for cloud +! fraction less than 1.0 but greater than 0.0 follows Mocko and Cotton +! (1996) from Sundqvist et al. (1989) with cloud fraction increasing +! as RH increases above a critical value. In locations with non-zero +! (but less than 1.0) cloud fraction, there MUST be a value assigned +! to cloud liquid water and ice or else there is zero impact in the +! RRTMG radiation scheme. + subroutine progcld_thompson & + & ( plyr,plvl,tlyr,qlyr,qstl,rhly,clw, & ! --- inputs: + & xlat,xlon,slmsk,dz,delp, & + & ntrac,ntcw,ntiw,ntrw,ntsw,ntgl, & + & IX, NLAY, NLP1, & + & uni_cld, lmfshal, lmfdeep2, cldcov, & + & re_cloud,re_ice,re_snow, & + & lwp_ex, iwp_ex, lwp_fc, iwp_fc, & + & dzlay, latdeg, julian, yearlen, gridkm, & + & clouds,clds,mtop,mbot,de_lgth,alpha & ! --- outputs: + & ) + +! ================= subprogram documentation block ================ ! +! ! +! subprogram: progcld_thompson computes cloud related quantities ! +! using Thompson cloud microphysics scheme. ! +! ! +! abstract: this program computes cloud fractions from cloud ! +! condensates, ! +! and computes the low, mid, high, total and boundary layer cloud ! +! fractions and the vertical indices of low, mid, and high cloud ! +! top and base. the three vertical cloud domains are set up in the ! +! initial subroutine "cld_init". ! +! ! +! usage: call progcld_thompson ! +! ! +! subprograms called: gethml ! +! ! +! attributes: ! +! language: fortran 90 ! +! machine: ibm-sp, sgi ! +! ! +! ! +! ==================== definition of variables ==================== ! +! ! +! ! +! input variables: ! +! plyr (IX,NLAY) : model layer mean pressure in mb (100Pa) ! +! plvl (IX,NLP1) : model level pressure in mb (100Pa) ! +! tlyr (IX,NLAY) : model layer mean temperature in k ! +! tvly (IX,NLAY) : model layer virtual temperature in k ! +! qlyr (IX,NLAY) : layer specific humidity in gm/gm ! +! qstl (IX,NLAY) : layer saturate humidity in gm/gm ! +! rhly (IX,NLAY) : layer relative humidity (=qlyr/qstl) ! +! clw (IX,NLAY,ntrac) : layer cloud condensate amount ! +! xlat (IX) : grid latitude in radians, default to pi/2 -> -pi/2! +! range, otherwise see in-line comment ! +! xlon (IX) : grid longitude in radians (not used) ! +! slmsk (IX) : sea/land mask array (sea:0,land:1,sea-ice:2) ! +! dz (ix,nlay) : layer thickness (km) ! +! delp (ix,nlay) : model layer pressure thickness in mb (100Pa) ! +! gridkm : grid length in km ! +! IX : horizontal dimention ! +! NLAY,NLP1 : vertical layer/level dimensions ! +! uni_cld : logical - true for cloud fraction from shoc ! +! lmfshal : logical - true for mass flux shallow convection ! +! lmfdeep2 : logical - true for mass flux deep convection ! +! cldcov : layer cloud fraction (used when uni_cld=.true. ! +! ! +! output variables: ! +! clouds(IX,NLAY,NF_CLDS) : cloud profiles ! +! clouds(:,:,1) - layer total cloud fraction ! +! clouds(:,:,2) - layer cloud liq water path (g/m**2) ! +! clouds(:,:,3) - mean eff radius for liq cloud (micron) ! +! clouds(:,:,4) - layer cloud ice water path (g/m**2) ! +! clouds(:,:,5) - mean eff radius for ice cloud (micron) ! +! clouds(:,:,6) - layer rain drop water path not assigned ! +! clouds(:,:,7) - mean eff radius for rain drop (micron) ! +! clouds(:,:,8) - layer snow flake water path not assigned ! +! clouds(:,:,9) - mean eff radius for snow flake (micron) ! +! *** fu's scheme need to be normalized by snow density (g/m**3/1.0e6) ! +! clds (IX,5) : fraction of clouds for low, mid, hi, tot, bl ! +! mtop (IX,3) : vertical indices for low, mid, hi cloud tops ! +! mbot (IX,3) : vertical indices for low, mid, hi cloud bases ! +! de_lgth(ix) : clouds decorrelation length (km) ! +! ! +! module variables: ! +! ivflip : control flag of vertical index direction ! +! =0: index from toa to surface ! +! =1: index from surface to toa ! +! lmfshal : mass-flux shallow conv scheme flag ! +! lmfdeep2 : scale-aware mass-flux deep conv scheme flag ! +! lcrick : control flag for eliminating CRICK ! +! =t: apply layer smoothing to eliminate CRICK ! +! =f: do not apply layer smoothing ! +! lcnorm : control flag for in-cld condensate ! +! =t: normalize cloud condensate ! +! =f: not normalize cloud condensate ! +! ! +! ==================== end of description ===================== ! +! + implicit none + +! --- inputs + integer, intent(in) :: IX, NLAY, NLP1 + integer, intent(in) :: ntrac, ntcw, ntiw, ntrw, ntsw, ntgl + + logical, intent(in) :: uni_cld, lmfshal, lmfdeep2 + + real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, & + & tlyr, qlyr, qstl, rhly, cldcov, delp, dz, dzlay, & + & re_cloud, re_ice, re_snow + real (kind=kind_phys), dimension(:), intent(inout) :: & + & lwp_ex, iwp_ex, lwp_fc, iwp_fc + + real (kind=kind_phys), dimension(:,:,:), intent(in) :: clw + + real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, & + & slmsk + + real(kind=kind_phys), dimension(:), intent(in) :: latdeg + real(kind=kind_phys), intent(in) :: julian, gridkm + integer, intent(in) :: yearlen + +! --- outputs + real (kind=kind_phys), dimension(:,:,:), intent(out) :: clouds + + real (kind=kind_phys), dimension(:,:), intent(out) :: clds + real (kind=kind_phys), dimension(:), intent(out) :: de_lgth + real (kind=kind_phys), dimension(:,:), intent(out) :: alpha + + integer, dimension(:,:), intent(out) :: mtop,mbot + +! --- local variables: + real (kind=kind_phys), dimension(IX,NLAY) :: cldtot, cldcnv, & + & cwp, cip, crp, csp, rew, rei, res, rer + + real (kind=kind_phys), dimension(NLAY) :: cldfra1d, qv1d, & + & qc1d, qi1d, qs1d, dz1d, p1d, t1d + + real (kind=kind_phys) :: ptop1(IX,NK_CLDS+1), rxlat(ix) + + real (kind=kind_phys) :: clwmin, tem1 + real (kind=kind_phys) :: corr, xland, snow_mass_factor + real (kind=kind_phys), parameter :: max_relh = 1.5 + real (kind=kind_phys), parameter :: snow_max_radius = 130.0 + + integer :: i, k, k2, id, nf, idx_rei +! +!===> ... begin here +! + + clwmin = 1.0E-9 + + do nf=1,nf_clds + do k=1,nlay + do i=1,ix + clouds(i,k,nf) = 0.0 + enddo + enddo + enddo + + do k = 1, NLAY + do i = 1, IX + cldtot(i,k) = 0.0 + cldcnv(i,k) = 0.0 + cwp (i,k) = 0.0 + cip (i,k) = 0.0 + crp (i,k) = 0.0 + csp (i,k) = 0.0 + rew (i,k) = re_cloud(i,k) + rei (i,k) = re_ice(i,k) + rer (i,k) = rrain_def ! default rain radius to 1000 micron + res (i,k) = re_snow(i,K) + enddo + enddo + +!> - Find top pressure for each cloud domain for given latitude. +!! ptopc(k,i): top presure of each cld domain (k=1-4 are sfc,L,m,h; +!! i=1,2 are low-lat (<45 degree) and pole regions) + + do i =1, IX + rxlat(i) = abs( xlat(i) / con_pi ) ! if xlat in pi/2 -> -pi/2 range +! rxlat(i) = abs(0.5 - xlat(i)/con_pi) ! if xlat in 0 -> pi range + enddo + + do id = 1, 4 + tem1 = ptopc(id,2) - ptopc(id,1) + + do i =1, IX + ptop1(i,id) = ptopc(id,1) + tem1*max( 0.0, 4.0*rxlat(i)-1.0 ) + enddo + enddo + +!> - Compute cloud liquid/ice condensate path in \f$ g/m^2 \f$ . +!> - 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.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 + snow_mass_factor = min(snow_mass_factor, & + & (snow_max_radius/re_snow(i,k)) & + & *(snow_max_radius/re_snow(i,k))) + res(i,k) = snow_max_radius + endif + csp(i,k) = max(0.,snow_mass_factor*clw(i,k,ntsw)*dz(i,k)*1.E6) + enddo + enddo + +!> - Sum the liquid water and ice paths that come from explicit micro + + do i = 1, IX + lwp_ex(i) = 0.0 + iwp_ex(i) = 0.0 + do k = 1, NLAY-1 + lwp_ex(i) = lwp_ex(i) + cwp(i,k) + iwp_ex(i) = iwp_ex(i) + cip(i,k) + csp(i,k) + enddo + enddo + +!> - Now determine the cloud fraction. Here, we will use the scheme of +!! G. Thompson that implements a variannt of Mocko and Cotton (1995) +!! based on work within HWRF and WRF. Where the bulk microphysics +!! scheme already has explicit clouds, assume cloud fraction of one, +!! but, otherwise, use a Sundqvist et al (1989) scheme and RH-critical +!! to account for sub-grid-scale clouds, include those in the water +!! and ice paths _seen_ by the radiation scheme, but do not actually +!! include these fake clouds into anything other than radiation. + + do i = 1, IX + if (slmsk(i)-0.5 .gt. 0.5 .and. slmsk(i)+0.5 .lt. 1.5) then + xland = 1.0 + else + xland = 2.0 + endif + + cldfra1d(:) = 0.0 + + if (ivflip .eq. 1) then + do k = 1, NLAY + qv1d(k) = qlyr(i,k) + qc1d(k) = max(0.0, clw(i,k,ntcw)) + qi1d(k) = max(0.0, clw(i,k,ntiw)) + qs1d(k) = max(0.0, clw(i,k,ntsw)) + dz1d(k) = dz(i,k)*1.E3 + p1d(k) = plyr(i,k)*100.0 + t1d(k) = tlyr(i,k) + enddo + else + do k = NLAY, 1, -1 + k2 = NLAY - k + 1 + qv1d(k2) = qlyr(i,k) + qc1d(k2) = max(0.0, clw(i,k,ntcw)) + qi1d(k2) = max(0.0, clw(i,k,ntiw)) + qs1d(k2) = max(0.0, clw(i,k,ntsw)) + dz1d(k2) = dz(i,k)*1.E3 + p1d(k2) = plyr(i,k)*100.0 + t1d(k2) = tlyr(i,k) + enddo + endif + + call cal_cldfra3(cldfra1d, qv1d, qc1d, qi1d, qs1d, dz1d, & + & p1d, t1d, xland, gridkm, & + & .false., max_relh, 1, nlay, .false.) + + do k = 1, NLAY + cldtot(i,k) = cldfra1d(k) + if (qc1d(k).gt.clwmin .and. cldfra1d(k).lt.ovcst) then + cwp(i,k) = qc1d(k) * dz1d(k)*1000. + if ((xland-1.5).GT.0.) then !--- Ocean + rew(i,k) = 9.5 + else !--- Land + rew(i,k) = 5.5 + endif + endif + if (qi1d(k).gt.clwmin .and. cldfra1d(k).lt.ovcst) then + cip(i,k) = qi1d(k) * dz1d(k)*1000. + idx_rei = int(t1d(k)-179.) + idx_rei = min(max(idx_rei,1),75) + corr = t1d(k) - int(t1d(k)) + rei(i,K) = max(5.0, retab(idx_rei)*(1.-corr) + & + & retab(idx_rei+1)*corr) + endif + enddo + enddo + + do k = 1, NLAY + do i = 1, IX + clouds(i,k,1) = cldtot(i,k) + clouds(i,k,2) = cwp(i,k) + clouds(i,k,3) = rew(i,k) + clouds(i,k,4) = cip(i,k) + clouds(i,k,5) = rei(i,k) + clouds(i,k,6) = crp(i,k) + clouds(i,k,7) = rer(i,k) + clouds(i,k,8) = csp(i,k) + clouds(i,k,9) = res(i,k) + enddo + enddo + +!> - Sum the liquid water and ice paths that come from fractional clouds + + do i = 1, IX + lwp_fc(i) = 0.0 + iwp_fc(i) = 0.0 + do k = 1, NLAY + lwp_fc(i) = lwp_fc(i) + cwp(i,k) + iwp_fc(i) = iwp_fc(i) + cip(i,k) + csp(i,k) + enddo + lwp_fc(i) = MAX(0.0, lwp_fc(i) - lwp_ex(i)) + iwp_fc(i) = MAX(0.0, iwp_fc(i) - iwp_ex(i)) + lwp_fc(i) = lwp_fc(i)*1.E-3 + iwp_fc(i) = iwp_fc(i)*1.E-3 + lwp_ex(i) = lwp_ex(i)*1.E-3 + iwp_ex(i) = iwp_ex(i)*1.E-3 + 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 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 and clouds top/bottom layer indices for low, mid, +!! and high clouds. +! --- compute low, mid, high, total, and boundary layer cloud fractions +! and clouds top/bottom layer indices for low, mid, and high clouds. +! The three cloud domain boundaries are defined by ptopc. The cloud +! overlapping method is defined by control flag 'iovr', which may +! be different for lw and sw radiation programs. + + call gethml & +! --- inputs: + & ( plyr, ptop1, cldtot, cldcnv, dz, de_lgth, alpha, & + & IX,NLAY, & +! --- outputs: + & clds, mtop, mbot & + & ) + +! + return + +!............................................ + end subroutine progcld_thompson +!............................................ +!mz + + !> \ingroup module_radiation_clouds !> This subroutine computes cloud related quantities using !! for unified cloud microphysics scheme. @@ -4050,6 +4446,7 @@ subroutine gethml & end subroutine gethml !----------------------------------- !! @} + !+---+-----------------------------------------------------------------+ !..Cloud fraction scheme by G. Thompson (NCAR-RAL), not intended for !.. combining with any cumulus or shallow cumulus parameterization @@ -4065,249 +4462,252 @@ end subroutine gethml ! !+---+-----------------------------------------------------------------+ - SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, & - & p,t,rho, XLAND, gridkm, & -! & rand_perturb_on, kme_stoch, rand_pert, & - & ids,ide, jds,jde, kds,kde, & - & ims,ime, jms,jme, kms,kme, & - & its,ite, jts,jte, kts,kte) + SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & + & p, t, XLAND, gridkm, & + & modify_qvapor, max_relh, & + & kts,kte, debug_flag) ! USE module_mp_thompson , ONLY : rsif, rslf IMPLICIT NONE ! - INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, & - & ims,ime, jms,jme, kms,kme, & -! & kme_stoch, & - & its,ite, jts,jte, kts,kte - -! INTEGER, INTENT(IN):: rand_perturb_on - REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN):: qv,p,t,rho - REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT):: qc,qi,qs -! REAL, DIMENSION(ims:ime,kms:kme_stoch,jms:jme), INTENT(IN):: rand_pert - REAL, DIMENSION(ims:ime,jms:jme), INTENT(IN):: XLAND - - REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT):: cldfra - REAL, INTENT(IN):: gridkm + INTEGER, INTENT(IN):: kts, kte + LOGICAL, INTENT(IN):: modify_qvapor + REAL, DIMENSION(kts:kte), INTENT(INOUT):: qv, qc, qi, cldfra + REAL, DIMENSION(kts:kte), INTENT(IN):: p, t, dz, qs + REAL, INTENT(IN):: gridkm, XLAND, max_relh + LOGICAL, INTENT(IN):: debug_flag !..Local vars. - REAL:: RH_00L, RH_00O, RH_00, RHI_max, entrmnt - REAL, DIMENSION(ims:ime,kms:kme,jms:jme):: qvsat - INTEGER:: i,j,k - REAL:: TK, TC, qvsi, qvsw, RHUM, xx, yy - REAL, DIMENSION(kts:kte):: qvs1d, cfr1d, T1d, & - & P1d, R1d, qc1d, qi1d, qs1d + REAL:: RH_00L, RH_00O, RH_00 + REAL:: entrmnt=0.5 + INTEGER:: k + REAL:: TC, qvsi, qvsw, RHUM, delz + REAL, DIMENSION(kts:kte):: qvs, rh, rhoa + integer:: ndebug = 0 character*512 dbg_msg - LOGICAL:: debug_flag !+---+ +!..Initialize cloud fraction, compute RH, and rho-air. + + DO k = kts,kte + CLDFRA(K) = 0.0 + + qvsw = rslf(P(k), t(k)) + qvsi = rsif(P(k), t(k)) + + tc = t(k) - 273.15 + if (tc .ge. -12.0) then + qvs(k) = qvsw + elseif (tc .lt. -35.0) then + qvs(k) = qvsi + else + qvs(k) = qvsw - (qvsw-qvsi)*(-12.0-tc)/(-12.0+35.) + endif + + if (modify_qvapor) then + if (qc(k).gt.1.E-8) then + qv(k) = MAX(qv(k), qvsw) + qvs(k) = qvsw + endif + if (qc(k).le.1.E-8 .and. qi(k).ge.1.E-9) then + qv(k) = MAX(qv(k), qvsi*1.005) !..To ensure a tiny bit ice supersaturated + qvs(k) = qvsi + endif + endif + + rh(k) = MAX(0.01, qv(k)/qvs(k)) + rhoa(k) = p(k)/(287.0*t(k)) + + ENDDO + + !..First cut scale-aware. Higher resolution should require closer to !.. saturated grid box for higher cloud fraction. Simple functions !.. chosen based on Mocko and Cotton (1995) starting point and desire !.. to get near 100% RH as grid spacing moves toward 1.0km, but higher !.. RH over ocean required as compared to over land. - RH_00L = 0.7 + SQRT(1./(25.0+gridkm*gridkm*gridkm)) - RH_00O = 0.81 + SQRT(1./(50.0+gridkm*gridkm*gridkm)) - - DO j = jts,jte DO k = kts,kte - DO i = its,ite - RHI_max = 0.0 - CLDFRA(I,K,J) = 0.0 - - if (qc(i,k,j).gt.1.E-6 .or. qi(i,k,j).ge.1.E-7 .or.qs(i,k,j) & - & .gt.1.E-5) then - CLDFRA(I,K,J) = 1.0 - qvsat(i,k,j) = qv(i,k,j) - else - TK = t(i,k,j) - TC = TK - 273.16 - - qvsw = rslf(P(i,k,j), TK) - qvsi = rsif(P(i,k,j), TK) - if (tc .ge. -12.0) then - qvsat(i,k,j) = qvsw - elseif (tc .lt. -20.0) then - qvsat(i,k,j) = qvsi - else - qvsat(i,k,j) = qvsw - (qvsw-qvsi)*(-12.0-tc)/(-12.0+20.) - endif - RHUM = MAX(0.01, MIN(qv(i,k,j)/qvsat(i,k,j), 0.9999)) + delz = MAX(100., dz(k)) + 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-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-5)) then + CLDFRA(K) = MIN(0.99, 0.20*(10.0 + log10(qc(k)+qi(k)))) + else - IF ((XLAND(I,J)-1.5).GT.0.) THEN !--- Ocean + IF ((XLAND-1.5).GT.0.) THEN !--- Ocean RH_00 = RH_00O - ELSE !--- Land + ELSE !--- Land RH_00 = RH_00L ENDIF - if (tc .ge. -12.0) then - RHUM = MIN(0.999, RHUM) - CLDFRA(I,K,J) = MAX(0.0, 1.0-SQRT((1.0-RHUM)/(1.-RH_00))) - elseif (tc.lt.-12..and.tc.gt.-70. .and. RHUM.gt.RH_00L) then - RHUM = MAX(0.01, MIN(qv(i,k,j)/qvsat(i,k,j), 1.0 - 1.E-6)) - CLDFRA(I,K,J) = MAX(0., 1.0-SQRT((1.0-RHUM)/(1.0-RH_00L))) + tc = t(k) - 273.15 + if (tc .lt. -12.0) RH_00 = RH_00L + + if (tc .gt. 20.0) then + CLDFRA(K) = 0.0 + elseif (tc .ge. -12.0) then + RHUM = MIN(rh(k), 1.0) + CLDFRA(K) = MAX(0., 1.0-SQRT((1.001-RHUM)/(1.001-RH_00))) + else + 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+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+112.) + CLDFRA(K) = MAX(0.,1.0-SQRT((1.06-RHUM)/(1.06-RH_00))) + endif endif - CLDFRA(I,K,J) = MIN(0.90, CLDFRA(I,K,J)) + if (CLDFRA(K).gt.0.) CLDFRA(K)=MAX(0.01,MIN(CLDFRA(K),0.99)) endif + if (cldfra(k).gt.0.0 .and. p(k).lt.7000.0) CLDFRA(K) = 0.0 ENDDO - ENDDO - ENDDO - -!..Prepare for a 1-D column to find various cloud layers. + call find_cloudLayers(qvs, cldfra, T, P, Dz, entrmnt, & + & debug_flag, qc, qi, qs, kts,kte) - DO j = jts,jte - DO i = its,ite -! if (i.gt.10.and.i.le.20 .and. j.gt.10.and.j.le.20) then -! debug_flag = .true. -! else -! debug_flag = .false. -! endif +!..Do a final total column adjustment since we may have added more than 1 mm +!.. LWP/IWP for multiple cloud decks. -! if (rand_perturb_on .eq. 1) then -! entrmnt = MAX(0.01, MIN(0.99, 0.5 + rand_pert(i,1,j)*0.5)) -! else - entrmnt = 0.5 -! endif + call adjust_cloudFinal(cldfra, qc, qi, rhoa, dz, kts,kte) - DO k = kts,kte - qvs1d(k) = qvsat(i,k,j) - cfr1d(k) = cldfra(i,k,j) - T1d(k) = t(i,k,j) - P1d(k) = p(i,k,j) - R1d(k) = rho(i,k,j) - qc1d(k) = qc(i,k,j) - qi1d(k) = qi(i,k,j) - qs1d(k) = qs(i,k,j) - ENDDO + 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 -! if (debug_flag) then -! WRITE (dbg_msg,*) 'DEBUG-GT: finding cloud layers at point (', i, ', ', j, ')' -! CALL wrf_debug (150, dbg_msg) -! endif - call find_cloudLayers(qvs1d, cfr1d, T1d, P1d, R1d, entrmnt, & - & debug_flag, qc1d, qi1d, qs1d, kts,kte) +!..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. + if (modify_qvapor) then DO k = kts,kte - cldfra(i,k,j) = cfr1d(k) - qc(i,k,j) = qc1d(k) - qi(i,k,j) = qi1d(k) + if (cldfra(k).gt.0.2 .and. cldfra(k).lt.1.0) then + qv(k) = MAX(qv(k),qvs(k)) + endif ENDDO - ENDDO - ENDDO + endif END SUBROUTINE cal_cldfra3 + !+---+-----------------------------------------------------------------+ !..From cloud fraction array, find clouds of multi-level depth and compute !.. a reasonable value of LWP or IWP that might be contained in that depth, !.. unless existing LWC/IWC is already there. - SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, R1d, entrmnt, & + SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, Dz1d, entrmnt,& & debugfl, qc1d, qi1d, qs1d, kts,kte) ! IMPLICIT NONE - +! INTEGER, INTENT(IN):: kts, kte LOGICAL, INTENT(IN):: debugfl REAL, INTENT(IN):: entrmnt - REAL, DIMENSION(kts:kte), INTENT(IN):: qvs1d,T1d,P1d,R1d - REAL, DIMENSION(kts:kte), INTENT(INOUT):: cfr1d - REAL, DIMENSION(kts:kte), INTENT(INOUT):: qc1d, qi1d, qs1d + REAL, DIMENSION(kts:kte), INTENT(IN):: qs1d,qvs1d,T1d,P1d,Dz1d + REAL, DIMENSION(kts:kte), INTENT(INOUT):: cfr1d, qc1d, qi1d !..Local vars. - REAL, DIMENSION(kts:kte):: theta, dz - REAL:: Z1, Z2, theta1, theta2, ht1, ht2 - INTEGER:: k, k2, k_tropo, k_m12C, k_m40C, k_cldb, k_cldt, kbot + REAL, DIMENSION(kts:kte):: theta + REAL:: theta1, theta2, delz + INTEGER:: k, k2, k_tropo, k_m12C, k_cldb, k_cldt, kbot, k_p200 LOGICAL:: in_cloud character*512 dbg_msg +!+---+ k_m12C = 0 - k_m40C = 0 + k_p200 = 0 DO k = kte, kts, -1 theta(k) = T1d(k)*((100000.0/P1d(k))**(287.05/1004.)) - if (T1d(k)-273.16 .gt. -40.0 .and. P1d(k).gt.7000.0) k_m40C = & - & MAX(k_m40C, k) - if (T1d(k)-273.16 .gt. -12.0 .and. P1d(k).gt.10000.0) k_m12C = & - & MAX(k_m12C, k) + if (T1d(k)-273.16 .gt. -12.0 .and. P1d(k).gt.10100.0) & + & k_m12C = MAX(k_m12C, k) + if (P1d(k).gt.19999.0 .and. k_p200.eq.0) k_p200 = k ENDDO - if (k_m40C .le. kts) k_m40C = kts if (k_m12C .le. kts) k_m12C = kts - Z2 = 44307.692 * (1.0 - (P1d(kte)/101325.)**0.190) - DO k = kte-1, kts, -1 - Z1 = 44307.692 * (1.0 - (P1d(k)/101325.)**0.190) - dz(k+1) = Z2 - Z1 - Z2 = Z1 - ENDDO - dz(kts) = dz(kts+1) - !..Find tropopause height, best surrogate, because we would not really !.. wish to put fake clouds into the stratosphere. The 10/1500 ratio !.. d(Theta)/d(Z) approximates a vertical line on typical SkewT chart !.. near typical (mid-latitude) tropopause height. Since messy data -!.. could give us a false signal of such a transition, do the check over +!.. could give us a false signal of such a transition, do the check over !.. three K-level change, not just a level-to-level check. This method !.. has potential failure in arctic-like conditions with extremely low !.. tropopause height, as would any other diagnostic, so ensure resulting -!.. k_tropo level is above 4km. +!.. k_tropo level is above 700hPa. - DO k = kte-3, kts, -1 + if ( (kte-k_p200) .lt. 3) k_p200 = kte-3 + DO k = k_p200-2, kts, -1 theta1 = theta(k) theta2 = theta(k+2) - ht1 = 44307.692 * (1.0 - (P1d(k)/101325.)**0.190) - ht2 = 44307.692 * (1.0 - (P1d(k+2)/101325.)**0.190) - if ( (((theta2-theta1)/(ht2-ht1)) .lt. 10./1500. ) .AND. & - & (ht1.lt.19000.) .and. (ht1.gt.4000.) ) then - goto 86 - endif + delz = 0.5*dz1d(k) + dz1d(k+1) + 0.5*dz1d(k+2) + if ( (((theta2-theta1)/delz).lt.10./1500.) .OR. & + & P1d(k).gt.70000.) EXIT ENDDO - 86 continue - k_tropo = MAX(kts+2, k+2) - -! if (debugfl) then -! print*, ' FOUND TROPOPAUSE ', k_tropo, ' near ', ht2, ' m' -! WRITE (dbg_msg,*) 'DEBUG-GT: FOUND TROPOPAUSE ', k_tropo, ' near ', ht2, ' m' -! CALL wrf_debug (150, dbg_msg) -! endif + k_tropo = MAX(kts+2, MIN(k+2, kte-1)) + + if (k_tropo .gt. k_p200) then + DO k = kte-3, k_p200-2, -1 + theta1 = theta(k) + theta2 = theta(k+2) + delz = 0.5*dz1d(k) + dz1d(k+1) + 0.5*dz1d(k+2) + if ( (((theta2-theta1)/delz).lt.10./1500.) .AND. & + & P1d(k).gt.9000.) EXIT + ENDDO + k_tropo = MAX(k_p200-1, MIN(k+2, kte-1)) + endif !..Eliminate possible fractional clouds above supposed tropopause. DO k = k_tropo+1, kte - if (cfr1d(k).gt.0.0 .and. cfr1d(k).lt.0.999) then + if (cfr1d(k).gt.0.0 .and. cfr1d(k).lt.1.0) then cfr1d(k) = 0. endif ENDDO -!..We would like to prevent fractional clouds below LCL in idealized -!.. situation with deep well-mixed convective PBL, that otherwise is -!.. likely to get clouds in more realistic capping inversion layer. +!..Be a bit more conservative with lower cloud fraction in scenario with +!.. well-mixed convective boundary layer below LCL. - kbot = kts+2 + kbot = kts+1 DO k = kbot, k_m12C - if ( (theta(k)-theta(k-1)) .gt. 0.05E-3*dz(k)) EXIT + if ( (theta(k)-theta(k-1)) .gt. 0.010E-3*Dz1d(k)) EXIT ENDDO kbot = MAX(kts+1, k-2) DO k = kts, kbot - if (cfr1d(k).gt.0.0 .and. cfr1d(k).lt.0.999) cfr1d(k) = 0. + if (cfr1d(k).gt.0.0 .and. cfr1d(k).lt.1.0) & + & cfr1d(k) = MAX(0.01,0.33*cfr1d(k)) + ENDDO + DO k = kts,k_tropo + if (cfr1d(k).gt.0.0) kbot = MIN(k,kbot) ENDDO - !..Starting below tropo height, if cloud fraction greater than 1 percent, -!.. compute an approximate total layer depth of cloud, determine a total -!.. liquid water/ice path (LWP/IWP), then reduce that amount with tuning +!.. compute an approximate total layer depth of cloud, determine a total +!.. liquid water/ice path (LWP/IWP), then reduce that amount with tuning !.. parameter to represent entrainment factor, then divide up LWP/IWP -!.. into delta-Z weighted amounts for individual levels per cloud layer. - +!.. into delta-Z weighted amounts for individual levels per cloud layer. k_cldb = k_tropo in_cloud = .false. k = k_tropo - DO WHILE (.not. in_cloud .AND. k.gt.k_m12C) + DO WHILE (.not. in_cloud .AND. k.gt.k_m12C+1) k_cldt = 0 if (cfr1d(k).ge.0.01) then in_cloud = .true. @@ -4324,30 +4724,20 @@ SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, R1d, entrmnt, & in_cloud = .false. endif if ((k_cldt - k_cldb + 1) .ge. 2) then -! if (debugfl) then -! print*, 'An ice cloud layer is found between ', k_cldt, -! k_cldb, P1d(k_cldt)*0.01, P1d(k_cldb)*0.01 -! WRITE (dbg_msg,*) 'DEBUG-GT: An ice cloud layer is found between -! ', k_cldt, k_cldb, P1d(k_cldt)*0.01, P1d(k_cldb)*0.01 -! CALL wrf_debug (150, dbg_msg) -! endif - call adjust_cloudIce(cfr1d, qi1d, qs1d, qvs1d, T1d,R1d,dz, & + call adjust_cloudIce(cfr1d, qi1d, qs1d, qvs1d, T1d, Dz1d, & & entrmnt, k_cldb,k_cldt,kts,kte) k = k_cldb - else - if (cfr1d(k_cldb).gt.0.and.qi1d(k_cldb).lt.1.E-6) & - & qi1d(k_cldb)=1.E-5*cfr1d(k_cldb) + elseif ((k_cldt - k_cldb + 1) .eq. 1) then + if (cfr1d(k_cldb).gt.0.and.cfr1d(k_cldb).lt.1.) & + & qi1d(k_cldb)=qi1d(k_cldb)+0.05*qvs1d(k_cldb)*cfr1d(k_cldb) + k = k_cldb endif - - k = k - 1 ENDDO - - - k_cldb = k_tropo + 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 @@ -4365,78 +4755,43 @@ SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, R1d, entrmnt, & in_cloud = .false. endif if ((k_cldt - k_cldb + 1) .ge. 2) then -! if (debugfl) then -! print*, 'A water cloud layer is found between ', k_cldt, -! k_cldb, P1d(k_cldt)*0.01, P1d(k_cldb)*0.01 -! WRITE (dbg_msg,*) 'DEBUG-GT: A water cloud layer is found -! between ', k_cldt, k_cldb, P1d(k_cldt)*0.01, P1d(k_cldb)*0.01 -! CALL wrf_debug (150, dbg_msg) -! endif - call adjust_cloudH2O(cfr1d, qc1d, qvs1d, T1d,R1d,dz, & + call adjust_cloudH2O(cfr1d, qc1d, qvs1d, T1d, Dz1d, & & entrmnt, k_cldb,k_cldt,kts,kte) k = k_cldb - else - if (cfr1d(k_cldb).gt.0.and.qc1d(k_cldb).lt.1.E-6) & - & qc1d(k_cldb)=1.E-5*cfr1d(k_cldb) + elseif ((k_cldt - k_cldb + 1) .eq. 1) then + if (cfr1d(k_cldb).gt.0.and.cfr1d(k_cldb).lt.1.) & + & qc1d(k_cldb)=qc1d(k_cldb)+0.05*qvs1d(k_cldb)*cfr1d(k_cldb) + k = k_cldb endif k = k - 1 ENDDO -!..Do a final total column adjustment since we may have added more than -!1mm -!.. LWP/IWP for multiple cloud decks. - - call adjust_cloudFinal(cfr1d, qc1d, qi1d, R1d,dz, kts,kte,k_tropo) - -! if (debugfl) then -! print*, ' Made-up fake profile of clouds' -! do k = kte, kts, -1 -! write(*,'(i3, 2x, f8.2, 2x, f9.2, 2x, f6.2, 2x, f15.7, 2x, -! f15.7)') & -! & K, T1d(k)-273.15, P1d(k)*0.01, cfr1d(k)*100., -! qc1d(k)*1000.,qi1d(k)*1000. -! enddo -! WRITE (dbg_msg,*) 'DEBUG-GT: Made-up fake profile of clouds' -! CALL wrf_debug (150, dbg_msg) -! do k = kte, kts, -1 -! write(dbg_msg,'(f8.2, 2x, f9.2, 2x, f6.2, 2x, f15.7, 2x, -! f15.7)') & -! & T1d(k)-273.15, P1d(k)*0.01, cfr1d(k)*100., -! qc1d(k)*1000.,qi1d(k)*1000. -! CALL wrf_debug (150, dbg_msg) -! enddo -! endif - END SUBROUTINE find_cloudLayers !+---+-----------------------------------------------------------------+ - SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs, T,Rho,dz, entr, k1,k2, & - & kts,kte) + SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs,T,dz,entr, k1,k2,kts,kte) ! IMPLICIT NONE ! INTEGER, INTENT(IN):: k1,k2, kts,kte REAL, INTENT(IN):: entr - REAL, DIMENSION(kts:kte), INTENT(IN):: cfr, qvs, T, Rho, dz - REAL, DIMENSION(kts:kte), INTENT(INOUT):: qi, qs - REAL:: iwc, max_iwc, tdz, this_iwc, this_dz, iwp_exists - INTEGER:: k, kmid + REAL, DIMENSION(kts:kte), INTENT(IN):: cfr, qs, qvs, T, dz + REAL, DIMENSION(kts:kte), INTENT(INOUT):: qi + REAL:: iwc, max_iwc, tdz, this_iwc, this_dz + INTEGER:: k tdz = 0. do k = k1, k2 tdz = tdz + dz(k) enddo - kmid = NINT(0.5*(k1+k2)) - max_iwc = ABS(qvs(k2-1)-qvs(k1)) -! print*, ' max_iwc = ', max_iwc, ' over DZ=',tdz + max_iwc = ABS(qvs(k2)-qvs(k1)) - iwp_exists = 0. do k = k1, k2 - iwp_exists = iwp_exists + (qi(k)+qs(k))*Rho(k)*dz(k) + max_iwc = MAX(1.E-5, max_iwc - (qi(k)+qs(k))) enddo - if (iwp_exists .gt. 1.0) RETURN + max_iwc = MIN(2.E-3, max_iwc) this_dz = 0.0 do k = k1, k2 @@ -4446,12 +4801,9 @@ SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs, T,Rho,dz, entr, k1,k2, & this_dz = this_dz + dz(k) endif this_iwc = max_iwc*this_dz/tdz - iwc = MAX(1.E-6, this_iwc*(1.-entr)) - if (cfr(k).gt.0.01.and.cfr(k).lt.0.99.and.T(k).ge.203.16) then - qi(k) = qi(k) + 0.1*cfr(k)*iwc - elseif (qi(k).lt.1.E-5.and.cfr(k).ge.0.99.and.T(k).ge.203.16) & - & then - qi(k) = qi(k) + 0.01*iwc + 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 enddo @@ -4459,30 +4811,28 @@ END SUBROUTINE adjust_cloudIce !+---+-----------------------------------------------------------------+ - SUBROUTINE adjust_cloudH2O(cfr, qc, qvs, T,Rho,dz, entr, k1,k2, & - & kts,kte) + SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte) ! IMPLICIT NONE ! INTEGER, INTENT(IN):: k1,k2, kts,kte REAL, INTENT(IN):: entr - REAL, DIMENSION(kts:kte):: cfr, qc, qvs, T, Rho, dz - REAL:: lwc, max_lwc, tdz, this_lwc, this_dz, lwp_exists - INTEGER:: k, kmid + REAL, DIMENSION(kts:kte), INTENT(IN):: cfr, qvs, T, dz + REAL, DIMENSION(kts:kte), INTENT(INOUT):: qc + REAL:: lwc, max_lwc, tdz, this_lwc, this_dz + INTEGER:: k tdz = 0. do k = k1, k2 tdz = tdz + dz(k) enddo - kmid = NINT(0.5*(k1+k2)) - max_lwc = ABS(qvs(k2-1)-qvs(k1)) + max_lwc = ABS(qvs(k2)-qvs(k1)) ! print*, ' max_lwc = ', max_lwc, ' over DZ=',tdz - lwp_exists = 0. do k = k1, k2 - lwp_exists = lwp_exists + qc(k)*Rho(k)*dz(k) + max_lwc = MAX(1.E-5, max_lwc - qc(k)) enddo - if (lwp_exists .gt. 1.0) RETURN + max_lwc = MIN(2.E-3, max_lwc) this_dz = 0.0 do k = k1, k2 @@ -4492,68 +4842,58 @@ SUBROUTINE adjust_cloudH2O(cfr, qc, qvs, T,Rho,dz, entr, k1,k2, & 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.01.and.cfr(k).lt.0.99.and.T(k).lt.298.16.and. & - & T(k).ge.253.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 - elseif (cfr(k).ge.0.99.and.qc(k).lt.1.E-5.and.T(k).lt.298.16 & - & .and.T(k).ge.253.16) then - qc(k) = qc(k) + 0.1*lwc endif enddo END SUBROUTINE adjust_cloudH2O - !+---+-----------------------------------------------------------------+ !..Do not alter any grid-explicitly resolved hydrometeors, rather only !.. the supposed amounts due to the cloud fraction scheme. - SUBROUTINE adjust_cloudFinal(cfr, qc, qi, Rho,dz, kts,kte,k_tropo) + SUBROUTINE adjust_cloudFinal(cfr, qc, qi, Rho,dz, kts,kte) ! IMPLICIT NONE ! - INTEGER, INTENT(IN):: kts,kte,k_tropo + INTEGER, INTENT(IN):: kts,kte REAL, DIMENSION(kts:kte), INTENT(IN):: cfr, Rho, dz REAL, DIMENSION(kts:kte), INTENT(INOUT):: qc, qi REAL:: lwp, iwp, xfac INTEGER:: k lwp = 0. - do k = kts, k_tropo - if (cfr(k).gt.0.0) then - lwp = lwp + qc(k)*Rho(k)*dz(k) - endif - enddo - iwp = 0. - do k = kts, k_tropo - if (cfr(k).gt.0.01 .and. cfr(k).lt.0.99) then + do k = kts, kte + if (cfr(k).gt.0.0 .and. cfr(k).lt.1.0) then + lwp = lwp + qc(k)*Rho(k)*dz(k) iwp = iwp + qi(k)*Rho(k)*dz(k) endif enddo - if (lwp .gt. 1.5) then - xfac = 1./lwp - do k = kts, k_tropo - if (cfr(k).gt.0.01 .and. cfr(k).lt.0.99) then + if (lwp .gt. 1.0) then + xfac = 1.0/lwp + do k = kts, kte + if (cfr(k).gt.0.0 .and. cfr(k).lt.1.0) then qc(k) = qc(k)*xfac endif enddo endif - if (iwp .gt. 1.5) then - xfac = 1./iwp - do k = kts, k_tropo - if (cfr(k).gt.0.01 .and. cfr(k).lt.0.99) then + if (iwp .gt. 1.0) then + xfac = 1.0/iwp + do k = kts, kte + if (cfr(k).gt.0.0 .and. cfr(k).lt.1.0) then qi(k) = qi(k)*xfac endif enddo endif END SUBROUTINE adjust_cloudFinal -! + !........................................! end module module_radiation_clouds ! !! @}