From 29db559173428a05c6509a7d2f886bf521dfc8ef Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Tue, 16 Nov 2021 13:44:09 -0700 Subject: [PATCH 01/18] added a new cloud fraction scheme (actually upgraded cal_cldfra3) designed with Thompson-MP for GFSv17-prototype8 --- physics/GFS_rrtmg_post.F90 | 9 +- physics/GFS_rrtmg_post.meta | 9 + physics/GFS_rrtmg_pre.F90 | 66 ++- physics/GFS_rrtmg_pre.meta | 67 ++- physics/module_mp_thompson.F90 | 4 +- physics/radiation_clouds.f | 919 ++++++++++++++++++++++----------- 6 files changed, 743 insertions(+), 331 deletions(-) diff --git a/physics/GFS_rrtmg_post.F90 b/physics/GFS_rrtmg_post.F90 index b882930bf..e0278c45e 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(out) :: 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 + + do i=1,im + total_albedo(i) = topfsw(i)%upfxc/topfsw(i)%dnfxc + enddo ! end subroutine GFS_rrtmg_post_run diff --git a/physics/GFS_rrtmg_post.meta b/physics/GFS_rrtmg_post.meta index 6564f5025..7a5144739 100644 --- a/physics/GFS_rrtmg_post.meta +++ b/physics/GFS_rrtmg_post.meta @@ -275,6 +275,15 @@ type = topfsw_type intent = in optional = F +[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 = out + optional = F [scmpsw] standard_name = components_of_surface_downward_shortwave_fluxes long_name = derived type for special components of surface downward shortwave fluxes diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index dbea66985..17c8fa2e7 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,& @@ -228,6 +233,15 @@ 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) + + 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 @@ -637,7 +651,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & enddo ! for Thompson MP - prepare variables for calc_effr if_thompson: if (imp_physics == imp_physics_thompson .and. ltaerosol) then - do k=1,LMK + do k=1,LM do i=1,IM qvs = qlyr(i,k) qv_mp (i,k) = qvs/(1.-qvs) @@ -652,7 +666,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & enddo enddo elseif (imp_physics == imp_physics_thompson) then - do k=1,LMK + do k=1,LM do i=1,IM qvs = qlyr(i,k) qv_mp (i,k) = qvs/(1.-qvs) @@ -892,8 +906,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & 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) @@ -903,11 +915,11 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & 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) +! 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 @@ -1031,14 +1043,26 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & else ! MYNN PBL or GF convective are not used - call progcld6 (plyr,plvl,tlyr,qlyr,qstl,rhly,tracer1, & ! --- inputs - xlat,xlon,slmsk,dz,delp, & +! call progcld6 (plyr,plvl,tlyr,qlyr,qstl,rhly,tracer1, & ! --- inputs +! 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, & +! lwp_ex, iwp_ex, lwp_fc, iwp_fc, & +! dzb, xlat_d, julian, yearlen, & +! clouds, cldsa, mtopa ,mbota, de_lgth, alpha) ! --- outputs + + 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, lmk, lmp, uni_cld, lmfshal, lmfdeep2, & - cldcov(:,1:LMK), effrl_inout(:,:), & - effri_inout(:,:), effrs_inout(:,:), & - dzb, xlat_d, julian, yearlen, & + 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 endif ! MYNN PBL or GF @@ -1071,7 +1095,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 +1109,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 48ddc586d..6fac1ef2d 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -646,7 +646,7 @@ optional = F [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 @@ -755,7 +755,7 @@ standard_name = total_cloud_fraction long_name = layer total cloud fraction units = frac - dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = inout @@ -764,7 +764,7 @@ standard_name = cloud_liquid_water_path long_name = layer cloud liquid water path units = g m-2 - dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = inout @@ -773,7 +773,7 @@ standard_name = mean_effective_radius_for_liquid_cloud long_name = mean effective radius for liquid cloud units = um - dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = inout @@ -782,7 +782,7 @@ standard_name = cloud_ice_water_path long_name = layer cloud ice water path units = g m-2 - dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = inout @@ -791,7 +791,7 @@ standard_name = mean_effective_radius_for_ice_cloud long_name = mean effective radius for ice cloud units = um - dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = inout @@ -1056,7 +1056,7 @@ standard_name = cloud_rain_water_path long_name = cloud rain water path units = g m-2 - dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = out @@ -1065,7 +1065,7 @@ standard_name = mean_effective_radius_for_rain_drop long_name = mean effective radius for rain drop units = um - dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = out @@ -1074,7 +1074,7 @@ standard_name = cloud_snow_water_path long_name = cloud snow water path units = g m-2 - dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = out @@ -1083,7 +1083,7 @@ standard_name = mean_effective_radius_for_snow_flake long_name = mean effective radius for snow flake units = um - dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys intent = out @@ -1101,7 +1101,52 @@ standard_name = instantaneous_3d_cloud_fraction long_name = instantaneous 3D cloud fraction for all MPs units = frac - dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[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 + optional = F +[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 + optional = F +[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 + optional = F +[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 + optional = F +[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 diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index f05aa8ba2..028395285 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -2458,7 +2458,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 @@ -3826,7 +3826,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)) diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index dacf6e38e..39a40ed67 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,7 +2664,7 @@ subroutine progcld5 & enddo !mz* if (uni_cld) then ! use unified sgs clouds generated outside -!mz* use unified sgs or thompson clouds generated outside +!mz* use unified sgs clouds generated outside if (uni_cld .or. icloud == 3) then do k = 1, NLAY do i = 1, IX @@ -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,19 @@ 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 + enddo + if (uni_cld) then ! use unified sgs clouds generated outside do k = 1, NLAY do i = 1, IX @@ -3085,54 +3121,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 @@ -3222,6 +3236,364 @@ 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: progcld6 computes cloud related quantities using ! +! Thompson/WSM6 cloud microphysics scheme. ! +! ! +! abstract: this program computes cloud fractions from cloud ! +! condensates, ! +! and computes the low, mid, high, total and boundary layer cloud ! +! fractions and the vertical indices of low, mid, and high cloud ! +! top and base. the three vertical cloud domains are set up in the ! +! initial subroutine "cld_init". ! +! ! +! usage: call progcld6 ! +! ! +! 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, id, nf, idx_rei +! +!===> ... begin here +! + + if (ivflip .ne. 1) then + STOP ' K must be bottom to top oriented by this point.' + endif + + 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 25 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) = max(0.0, clw(i,k,ntrw) * dz(i,k)*1.E6) + snow_mass_factor = 0.75 + 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 + do k = 1, NLAY-1 + 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 + + call cal_cldfra3(cldfra1d, qv1d, qc1d, qi1d, qs1d, dz1d, & + & p1d, t1d, xland, gridkm, & + & .false., max_relh, 1, nlay-1, .false.) + + do k = 1, NLAY-1 + 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-1 + 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 +4422,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 +4438,257 @@ 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) + delz = MAX(100., dz(k)) + RH_00L = 0.63+MIN(0.36,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) + RH_00O = 0.79+MIN(0.20,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) + RHUM = rh(k) - 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)) + if (qc(k).ge.1.E-8 .or. qi(k).ge.1.E-9 & + & .or. (qs(k).gt.1.E-6 .and. t(k).lt.273.)) then + CLDFRA(K) = 1.0 + 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 ENDDO - ENDDO - ENDDO + call find_cloudLayers(qvs, cldfra, T, P, Dz, entrmnt, & + & debug_flag, qc, qi, qs, kts,kte) -!..Prepare for a 1-D column to find various cloud layers. +!..Do a final total column adjustment since we may have added more than 1 mm +!.. LWP/IWP for multiple cloud decks. - 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 + call adjust_cloudFinal(cldfra, qc, qi, rhoa, dz, kts,kte) -! 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 +!..Last adjustment to cloud fraction already set to 1.0 when the explicit +!.. clouds are present but extremely low mixing ratios. Also, no way in this +!.. world should we permit clouds above the 70 hPa level. - 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 + DO k = kts,kte + if (cldfra(k).eq.1.0 .and. ((qc(k)+qi(k)).gt.1.E-10) .and. & + & ((qc(k)+qi(k)).lt.1.E-6)) then + CLDFRA(K) = MIN(0.99, 0.25*(10.0 + log10(qc(k)+qi(k)))) + endif + if (cldfra(k).gt.0.0 .and. p(k).gt.7000.0) CLDFRA(K) = 0.0 + if (debug_flag .and. ndebug.lt.25) then + write(6,'(a,x,i3,x,f8.2,f7.1,f7.2,f6.1,x,f5.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 + if (k.eq.kte) ndebug = ndebug + 1 + endif + ENDDO -! 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 +4705,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 +4736,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 +4782,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 +4792,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 +4823,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 SUBROUTINE adjust_cloudFinal! + !........................................! end module module_radiation_clouds ! !! @} From 052145c70dd0d95fa93d349be29c224a57b264d3 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Tue, 16 Nov 2021 14:16:38 -0700 Subject: [PATCH 02/18] correct a dumb mistake --- physics/radiation_clouds.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 39a40ed67..0c31623fd 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -4566,7 +4566,7 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & & ((qc(k)+qi(k)).lt.1.E-6)) then CLDFRA(K) = MIN(0.99, 0.25*(10.0 + log10(qc(k)+qi(k)))) endif - if (cldfra(k).gt.0.0 .and. p(k).gt.7000.0) CLDFRA(K) = 0.0 + if (cldfra(k).gt.0.0 .and. p(k).lt.7000.0) CLDFRA(K) = 0.0 if (debug_flag .and. ndebug.lt.25) then write(6,'(a,x,i3,x,f8.2,f7.1,f7.2,f6.1,x,f5.3,f12.7,f12.7, & f12.7)') ' DEBUG-GT: ', k, p(k)*0.01, dz(k), t(k)-273.15, & From 8da0705fd02df720f80cb4a1ee2330582ab3789c Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Tue, 16 Nov 2021 15:23:28 -0700 Subject: [PATCH 03/18] make icloud=3 call cldfra3 compatible with updated subroutine --- physics/GFS_rrtmg_pre.F90 | 42 +++++++++++++++++++++------------------ 1 file changed, 23 insertions(+), 19 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 17c8fa2e7..ca3bf0e70 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -196,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 @@ -211,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, & @@ -906,27 +908,29 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & endif enddo + if (imp_physics == imp_physics_thompson) then + max_relh = 1.5 + else + max_relh = 1.1 + endif + 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) + cldfra1d(:) = 0.0 + do k = 1, lm-1 + qv1d(k) = qlyr(i,k) + qc1d(k) = max(0.0, tracer1(i,k,ntcw)) + qi1d(k) = max(0.0, tracer1(i,k,ntiw)) + qs1d(k) = max(0.0, tracer1(i,k,ntsw)) + dz1d(k) = dz(i,k)*1.E3 + p1d(k) = plyr(i,k)*100.0 + t1d(k) = tlyr(i,k) 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) + call cal_cldfra3(cldfra1d, qv1d, qc1d, qi1d, qs1d, dz1d, & + & p1d, t1d, xland(i), gridkm, & + & .false., max_relh, 1, lm-1, .false.) + do k = 1, lm-1 + cldcov(i,k) = cldfra1d(k) enddo enddo From 3787899c6e84edc1e739a91fc8a1bd773ed6f49f Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Wed, 17 Nov 2021 10:55:39 -0700 Subject: [PATCH 04/18] bug fix lwp_ex and iwp_ex inside progcld6; factor of 1000 --- physics/radiation_clouds.f | 2 ++ 1 file changed, 2 insertions(+) diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 0c31623fd..13f0ae76b 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -3107,6 +3107,8 @@ subroutine progcld6 & 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 From 4742485d4e279c7d97f00f557ec0ba03552e5cd0 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 19 Nov 2021 15:04:38 -0700 Subject: [PATCH 05/18] Fix compilation error with gfortran in physics/radiation_clouds.f --- physics/radiation_clouds.f | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 13f0ae76b..39296aec0 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -4570,8 +4570,8 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & endif if (cldfra(k).gt.0.0 .and. p(k).lt.7000.0) CLDFRA(K) = 0.0 if (debug_flag .and. ndebug.lt.25) then - write(6,'(a,x,i3,x,f8.2,f7.1,f7.2,f6.1,x,f5.3,f12.7,f12.7, - & f12.7)') ' DEBUG-GT: ', k, p(k)*0.01, dz(k), t(k)-273.15, & + write(6,'(a,x,i3,x,f8.2,f7.1,f7.2,f6.1,x,f5.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 if (k.eq.kte) ndebug = ndebug + 1 endif From 30ae919d7dda5d816c2a6acf05c7254bd371c982 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 22 Nov 2021 06:36:21 -0700 Subject: [PATCH 06/18] Avoid dividing by zero in physics/GFS_rrtmg_post.F90 --- physics/GFS_rrtmg_post.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/physics/GFS_rrtmg_post.F90 b/physics/GFS_rrtmg_post.F90 index e0278c45e..f8a63789a 100644 --- a/physics/GFS_rrtmg_post.F90 +++ b/physics/GFS_rrtmg_post.F90 @@ -199,10 +199,10 @@ subroutine GFS_rrtmg_post_run (im, km, kmp1, lm, ltp, kt, kb, kd, nspc1, & endif ! end_if_lssav ! --- The total sky (with clouds) shortwave albedo - - do i=1,im - total_albedo(i) = topfsw(i)%upfxc/topfsw(i)%dnfxc - enddo + total_albedo = 0.0 + if (lsswr) then + where(topfsw(:)%dnfxc>0) total_albedo(:) = topfsw(:)%upfxc/topfsw(:)%dnfxc + endif ! end subroutine GFS_rrtmg_post_run From 74cae0c8fd5ec3907b45ff5b8edf2bf274ac86f2 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 22 Nov 2021 06:44:22 -0700 Subject: [PATCH 07/18] Move total albedo to right place in physics/GFS_rrtmg_post.meta, make intent(inout) --- physics/GFS_rrtmg_post.F90 | 2 +- physics/GFS_rrtmg_post.meta | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/physics/GFS_rrtmg_post.F90 b/physics/GFS_rrtmg_post.F90 index f8a63789a..8584f8463 100644 --- a/physics/GFS_rrtmg_post.F90 +++ b/physics/GFS_rrtmg_post.F90 @@ -43,7 +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(out) :: total_albedo + 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 diff --git a/physics/GFS_rrtmg_post.meta b/physics/GFS_rrtmg_post.meta index 3332589de..761affd53 100644 --- a/physics/GFS_rrtmg_post.meta +++ b/physics/GFS_rrtmg_post.meta @@ -243,14 +243,6 @@ dimensions = (horizontal_loop_extent) type = topfsw_type intent = in -[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 = out [scmpsw] standard_name = components_of_surface_downward_shortwave_fluxes long_name = derived type for special components of surface downward shortwave fluxes @@ -266,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 From fe6677fe1f9a1ac115d838c7880a9eedc92d5303 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Mon, 22 Nov 2021 09:03:53 -0700 Subject: [PATCH 08/18] minor adjustments to reduce coverage of partly cloudy conditions --- physics/radiation_clouds.f | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 13f0ae76b..3384b1f72 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -3437,14 +3437,14 @@ subroutine progcld_thompson & enddo !> - Compute cloud liquid/ice condensate path in \f$ g/m^2 \f$ . -!> - Since using Thompson MP, assume 25 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) = max(0.0, clw(i,k,ntrw) * dz(i,k)*1.E6) - snow_mass_factor = 0.75 + snow_mass_factor = 0.80 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 @@ -4510,11 +4510,11 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & DO k = kts,kte delz = MAX(100., dz(k)) - RH_00L = 0.63+MIN(0.36,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) + RH_00L = 0.68+MIN(0.31,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) RH_00O = 0.79+MIN(0.20,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) RHUM = rh(k) - if (qc(k).ge.1.E-8 .or. qi(k).ge.1.E-9 & + if (qc(k).ge.1.E-7 .or. qi(k).ge.1.E-7 & & .or. (qs(k).gt.1.E-6 .and. t(k).lt.273.)) then CLDFRA(K) = 1.0 else @@ -4569,9 +4569,10 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & CLDFRA(K) = MIN(0.99, 0.25*(10.0 + log10(qc(k)+qi(k)))) endif if (cldfra(k).gt.0.0 .and. p(k).lt.7000.0) CLDFRA(K) = 0.0 + if (debug_flag .and. ndebug.lt.25) then - write(6,'(a,x,i3,x,f8.2,f7.1,f7.2,f6.1,x,f5.3,f12.7,f12.7, - & f12.7)') ' DEBUG-GT: ', k, p(k)*0.01, dz(k), t(k)-273.15, & + 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 if (k.eq.kte) ndebug = ndebug + 1 endif From b4d7ab0694eba3a49bad228cc5807bf2ac0722c2 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 22 Nov 2021 13:49:29 -0700 Subject: [PATCH 09/18] Use correct vertical dimension for several cloud arrays in physics/GFS_rrtmg_post.meta and physics/GFS_rrtmg_pre.meta --- physics/GFS_rrtmg_post.meta | 2 +- physics/GFS_rrtmg_pre.meta | 20 ++++++++++---------- 2 files changed, 11 insertions(+), 11 deletions(-) diff --git a/physics/GFS_rrtmg_post.meta b/physics/GFS_rrtmg_post.meta index 761affd53..6cff0cad4 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 diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index 3debc33e4..85f29dc9c 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -667,7 +667,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 = inout @@ -675,7 +675,7 @@ standard_name = cloud_liquid_water_path long_name = layer cloud liquid water path units = g m-2 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) type = real kind = kind_phys intent = inout @@ -683,7 +683,7 @@ standard_name = mean_effective_radius_for_liquid_cloud long_name = mean effective radius for liquid cloud units = um - dimensions = (horizontal_loop_extent,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) type = real kind = kind_phys intent = inout @@ -691,7 +691,7 @@ standard_name = cloud_ice_water_path long_name = layer cloud ice water path units = g m-2 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) type = real kind = kind_phys intent = inout @@ -699,7 +699,7 @@ standard_name = mean_effective_radius_for_ice_cloud long_name = mean effective radius for ice cloud units = um - dimensions = (horizontal_loop_extent,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) type = real kind = kind_phys intent = inout @@ -934,7 +934,7 @@ standard_name = cloud_rain_water_path long_name = cloud rain water path units = g m-2 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) type = real kind = kind_phys intent = out @@ -942,7 +942,7 @@ standard_name = mean_effective_radius_for_rain_drop long_name = mean effective radius for rain drop units = um - dimensions = (horizontal_loop_extent,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) type = real kind = kind_phys intent = out @@ -950,7 +950,7 @@ standard_name = cloud_snow_water_path long_name = cloud snow water path units = g m-2 - dimensions = (horizontal_loop_extent,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) type = real kind = kind_phys intent = out @@ -958,7 +958,7 @@ standard_name = mean_effective_radius_for_snow_flake long_name = mean effective radius for snow flake units = um - dimensions = (horizontal_loop_extent,vertical_layer_dimension) + dimensions = (horizontal_loop_extent,adjusted_vertical_layer_dimension_for_radiation) type = real kind = kind_phys intent = out @@ -974,7 +974,7 @@ standard_name = instantaneous_3d_cloud_fraction long_name = instantaneous 3D cloud fraction for all MPs 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 = out From 6f91aea381633b45dea5d172354713a2079f1fc0 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Tue, 23 Nov 2021 09:26:49 -0700 Subject: [PATCH 10/18] silly stray comment/exclamation point --- physics/radiation_clouds.f | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 3384b1f72..114bd3108 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -4876,7 +4876,7 @@ SUBROUTINE adjust_cloudFinal(cfr, qc, qi, Rho,dz, kts,kte) enddo endif - END SUBROUTINE adjust_cloudFinal! + END SUBROUTINE adjust_cloudFinal !........................................! end module module_radiation_clouds ! From 211a41311a9ef07e86e84d157e9cb569d729f4b0 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Tue, 23 Nov 2021 09:29:46 -0700 Subject: [PATCH 11/18] fix a few comment lines per review of pull request 781 --- physics/radiation_clouds.f | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 114bd3108..be45b0b50 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -3262,8 +3262,8 @@ subroutine progcld_thompson & ! ================= subprogram documentation block ================ ! ! ! -! subprogram: progcld6 computes cloud related quantities using ! -! Thompson/WSM6 cloud microphysics scheme. ! +! subprogram: progcld_thompson computes cloud related quantities ! +! using Thompson cloud microphysics scheme. ! ! ! ! abstract: this program computes cloud fractions from cloud ! ! condensates, ! @@ -3272,7 +3272,7 @@ subroutine progcld_thompson & ! top and base. the three vertical cloud domains are set up in the ! ! initial subroutine "cld_init". ! ! ! -! usage: call progcld6 ! +! usage: call progcld_thompson ! ! ! ! subprograms called: gethml ! ! ! From 0a54eecbebb14140685c8c9e0203cf30ca2b0a7e Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Tue, 23 Nov 2021 09:39:33 -0700 Subject: [PATCH 12/18] per discussions with Ruiyu, increasing min snow and graupel size seems advantageous --- physics/module_mp_thompson.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 271535581..353f83c78 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; From 112077914442e3555effbba0a4a231483206ce8a Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Tue, 23 Nov 2021 10:03:08 -0700 Subject: [PATCH 13/18] tiny adjustment of indent/format consistency --- physics/GFS_rrtmg_pre.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index ca3bf0e70..3ffc0d1be 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -35,8 +35,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & gasvmr_o2, gasvmr_co, gasvmr_cfc11, gasvmr_cfc12, gasvmr_cfc22, & gasvmr_ccl4, gasvmr_cfc113, aerodp, clouds6, clouds7, clouds8, & clouds9, cldsa, cldfra, cldfra2d, lwp_ex,iwp_ex, lwp_fc,iwp_fc, & - faersw1, faersw2, faersw3, & - faerlw1, faerlw2, faerlw3, alpha, errmsg, errflg) + faersw1, faersw2, faersw3, faerlw1, faerlw2, faerlw3, alpha, & + errmsg, errflg) use machine, only: kind_phys @@ -653,7 +653,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & enddo ! for Thompson MP - prepare variables for calc_effr if_thompson: if (imp_physics == imp_physics_thompson .and. ltaerosol) then - do k=1,LM + do k=1,LMK do i=1,IM qvs = qlyr(i,k) qv_mp (i,k) = qvs/(1.-qvs) @@ -668,7 +668,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & enddo enddo elseif (imp_physics == imp_physics_thompson) then - do k=1,LM + do k=1,LMK do i=1,IM qvs = qlyr(i,k) qv_mp (i,k) = qvs/(1.-qvs) From 7e119fc859eb7c3fe705538e2da701b5a9dad2f5 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Tue, 23 Nov 2021 11:00:23 -0700 Subject: [PATCH 14/18] re-use icloud=3 option for progcld_thompson --- physics/GFS_rrtmg_pre.F90 | 131 +++++++++++--------------------------- 1 file changed, 37 insertions(+), 94 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 3ffc0d1be..722ada210 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -237,6 +237,12 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & 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 @@ -870,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 - - if (imp_physics == imp_physics_thompson) then - max_relh = 1.5 - else - max_relh = 1.1 - endif - - do i =1, im - cldfra1d(:) = 0.0 - do k = 1, lm-1 - qv1d(k) = qlyr(i,k) - qc1d(k) = max(0.0, tracer1(i,k,ntcw)) - qi1d(k) = max(0.0, tracer1(i,k,ntiw)) - qs1d(k) = max(0.0, tracer1(i,k,ntsw)) - dz1d(k) = dz(i,k)*1.E3 - p1d(k) = plyr(i,k)*100.0 - t1d(k) = tlyr(i,k) - enddo - - call cal_cldfra3(cldfra1d, qv1d, qc1d, qi1d, qs1d, dz1d, & - & p1d, t1d, xland(i), gridkm, & - & .false., max_relh, 1, lm-1, .false.) - do k = 1, lm-1 - cldcov(i,k) = cldfra1d(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) @@ -1028,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 .eq. 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 @@ -1044,21 +982,13 @@ 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, & -! 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, & -! lwp_ex, iwp_ex, lwp_fc, iwp_fc, & -! dzb, xlat_d, julian, yearlen, & -! clouds, cldsa, mtopa ,mbota, de_lgth, alpha) ! --- outputs - - call progcld_thompson (plyr,plvl,tlyr,qlyr,qstl,rhly, & ! --- inputs + + if (icloud .eq. 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, & @@ -1068,6 +998,19 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & 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, & + 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 From bc8c26d8796767d0b2b2c648a8fa31200e3d0e81 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Tue, 23 Nov 2021 16:50:44 -0700 Subject: [PATCH 15/18] per review by climbfuji, make suggested changes --- physics/GFS_rrtmg_pre.F90 | 4 ++-- physics/radiation_clouds.f | 46 +++++++++++++++++++++++--------------- 2 files changed, 30 insertions(+), 20 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 722ada210..973ac02fd 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -953,7 +953,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & if(do_mynnedmf .or. imfdeepcnv == imfdeepcnv_gf ) then ! MYNN PBL or GF conv - if (icloud .eq. 3) then + 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, & @@ -987,7 +987,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & else ! MYNN PBL or GF convective are not used - if (icloud .eq. 3) then + 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, & diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index be45b0b50..e6e11131a 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -2665,7 +2665,7 @@ subroutine progcld5 & !mz* if (uni_cld) then ! use unified sgs clouds generated outside !mz* use unified sgs clouds generated outside - if (uni_cld .or. icloud == 3) then + if (uni_cld) then do k = 1, NLAY do i = 1, IX cldtot(i,k) = cldcov(i,k) @@ -3385,15 +3385,11 @@ subroutine progcld_thompson & real (kind=kind_phys), parameter :: max_relh = 1.5 real (kind=kind_phys), parameter :: snow_max_radius = 130.0 - integer :: i, k, id, nf, idx_rei + integer :: i, k, k2, id, nf, idx_rei ! !===> ... begin here ! - if (ivflip .ne. 1) then - STOP ' K must be bottom to top oriented by this point.' - endif - clwmin = 1.0E-9 do nf=1,nf_clds @@ -3485,21 +3481,35 @@ subroutine progcld_thompson & endif cldfra1d(:) = 0.0 - do k = 1, NLAY-1 - 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 + + 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-1, .false.) + & .false., max_relh, 1, nlay, .false.) - do k = 1, NLAY-1 + 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. @@ -3539,7 +3549,7 @@ subroutine progcld_thompson & do i = 1, IX lwp_fc(i) = 0.0 iwp_fc(i) = 0.0 - do k = 1, NLAY-1 + 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 From 2b2e3b195d47ed2bdea7ffc99e5183d467ad2b62 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Wed, 24 Nov 2021 12:04:23 -0700 Subject: [PATCH 16/18] disallow rain in Thompson to be used in radiation scheme; also set fractional cloudy LWP and IWP inside progcld6 --- physics/radiation_clouds.f | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index e6e11131a..c5b94053c 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -3164,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 @@ -3439,7 +3451,7 @@ subroutine progcld_thompson & 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) = max(0.0, clw(i,k,ntrw) * dz(i,k)*1.E6) + crp(i,k) = 0.0 snow_mass_factor = 0.80 cip(i,k) = max(0.0, (clw(i,k,ntiw) & & + (1.0-snow_mass_factor)*clw(i,k,ntsw))*dz(i,k)*1.E6) From 7d0d7364007310a28afc3ae6af97df73a73c0eb3 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Fri, 26 Nov 2021 09:45:52 -0700 Subject: [PATCH 17/18] minor rearrangement of cloud fraction with very low mixing ratios --- physics/radiation_clouds.f | 26 ++++++++++---------------- 1 file changed, 10 insertions(+), 16 deletions(-) diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index c5b94053c..e73f91ea2 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -4536,9 +4536,12 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & RH_00O = 0.79+MIN(0.20,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) RHUM = rh(k) - if (qc(k).ge.1.E-7 .or. qi(k).ge.1.E-7 & + if (qc(k).ge.1.E-6 .or. qi(k).ge.1.E-6 & & .or. (qs(k).gt.1.E-6 .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.25*(10.0 + log10(qc(k)+qi(k)))) else IF ((XLAND-1.5).GT.0.) THEN !--- Ocean @@ -4571,6 +4574,7 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & 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 call find_cloudLayers(qvs, cldfra, T, P, Dz, entrmnt, & @@ -4581,24 +4585,14 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & call adjust_cloudFinal(cldfra, qc, qi, rhoa, dz, kts,kte) -!..Last adjustment to cloud fraction already set to 1.0 when the explicit -!.. clouds are present but extremely low mixing ratios. Also, no way in this -!.. world should we permit clouds above the 70 hPa level. - - DO k = kts,kte - if (cldfra(k).eq.1.0 .and. ((qc(k)+qi(k)).gt.1.E-10) .and. & - & ((qc(k)+qi(k)).lt.1.E-6)) then - CLDFRA(K) = MIN(0.99, 0.25*(10.0 + log10(qc(k)+qi(k)))) - endif - if (cldfra(k).gt.0.0 .and. p(k).lt.7000.0) CLDFRA(K) = 0.0 - - if (debug_flag .and. ndebug.lt.25) then + 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 - if (k.eq.kte) ndebug = ndebug + 1 - endif - ENDDO + 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 From 472987f8df3ce53850801d7531923306fc7b6849 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Tue, 7 Dec 2021 15:06:49 -0700 Subject: [PATCH 18/18] permit few max ice number conc; scale back cloud fraction a little; reduce assumed snow as ice in IWP for radiation --- physics/module_mp_thompson.F90 | 12 ++++++------ physics/radiation_clouds.f | 14 +++++++------- 2 files changed, 13 insertions(+), 13 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index fb1a4a5f2..3183ca4bf 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -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 @@ -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) & @@ -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 e73f91ea2..f58ec8d11 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -3452,7 +3452,7 @@ subroutine progcld_thompson & 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.80 + 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 @@ -4532,16 +4532,16 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & DO k = kts,kte delz = MAX(100., dz(k)) - RH_00L = 0.68+MIN(0.31,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) - RH_00O = 0.79+MIN(0.20,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-6 & - & .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.25*(10.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