From 2ee8e48f45cc10c05a78ffabaae0bfdf034cc515 Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Fri, 31 Jan 2020 20:35:33 +0000 Subject: [PATCH 01/12] Changes in cloud/radiation interaction in GSD physics suite that uses Thompson MP, MYNN pbl and GF convection: 1. Switch the order of calls, first MYNNrad_pre (or SGSCloud_RadPre), then rrtmg_pre. This will add sub-grid clouds from MYNN PBL (or MYNN PBL and GF) to QC and QI, and these updated hydrometeors will be used to compute cloud paths and effective radii. 2. In rrtmg_pre with the use of THompson MP: - use Thompson's subroutines make_IceNumber and make_DropletNumber to compute number concentrations for subgrid clouds. - use calc_effectRad to compute effective radii for QC and QI with sub-grid clouds. - added option (clduni) to use the same subroutine to compute water paths as with the GFDL MP. For this input.nl should set effr_in=.true. - the progcld5 is used mostly to compute Xu-Randall cloud fraction. 3. Added *SGSCloud_* modules to replace *MYNNrad* to add all subgrid clouds to QC and QI (from MYNN PBL and GF conv). 4. Added convective clouds qci_conv to GF scheme and SGSCloud_RadPre. 5. Computation of total cloud fraction in progcld5 is change not to depend on shallow/deep convection. Not needed in the current version of GSD suite. --- physics/GFS_rrtmg_pre.F90 | 306 ++++++++++++++++++++++++-- physics/GFS_rrtmg_pre.meta | 61 ++++++ physics/cu_gf_driver.F90 | 8 +- physics/cu_gf_driver.meta | 9 + physics/module_MYNNrad_pre.F90 | 7 + physics/module_SGSCloud_RadPost.F90 | 75 +++++++ physics/module_SGSCloud_RadPost.meta | 96 +++++++++ physics/module_SGSCloud_RadPre.F90 | 211 ++++++++++++++++++ physics/module_SGSCloud_RadPre.meta | 308 +++++++++++++++++++++++++++ physics/radiation_clouds.f | 63 +++--- 10 files changed, 1095 insertions(+), 49 deletions(-) create mode 100644 physics/module_SGSCloud_RadPost.F90 create mode 100644 physics/module_SGSCloud_RadPost.meta create mode 100644 physics/module_SGSCloud_RadPre.F90 create mode 100644 physics/module_SGSCloud_RadPre.meta diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index b179a74db..6b5382e65 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -20,7 +20,8 @@ end subroutine GFS_rrtmg_pre_init ! in the CCPP version - they are defined in the interstitial_create routine subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input Tbd, Cldprop, Coupling, & - Radtend, & ! input/output + Radtend, qc, qi, nc, ni, nwfa, & ! input/output + imfdeepcnv, imfdeepcnv_gf, & f_ice, f_rain, f_rimef, flgmin, cwm, & ! F-A mp scheme only lm, im, lmk, lmp, & ! input kd, kt, kb, raddt, delp, dz, plvl, plyr, & ! output @@ -50,7 +51,8 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input & epsm1 => con_epsm1, & & fvirt => con_fvirt & &, rog => con_rog & - &, rocp => con_rocp + &, rocp => con_rocp & + &, con_rd use radcons, only: itsfc,ltp, lextop, qmin, & qme5, qme6, epsq, prsmin use funcphys, only: fpvs @@ -70,6 +72,10 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input & proflw_type, NBDLW use surface_perturbation, only: cdfnor + !tgs for Thompson MP + use module_mp_thompson, only : calc_effectRad + use module_mp_thompson_make_number_concentrations, only: make_IceNumber, make_DropletNumber, make_RainNumber + implicit none type(GFS_control_type), intent(in) :: Model @@ -81,7 +87,15 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input type(GFS_cldprop_type), intent(in) :: Cldprop type(GFS_coupling_type), intent(in) :: Coupling + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: qc + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: qi + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: nc + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: ni + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: nwfa + + integer, intent(in) :: im, lm, lmk, lmp + integer, intent(in) :: imfdeepcnv, imfdeepcnv_gf integer, intent(out) :: kd, kt, kb ! F-A mp scheme only @@ -123,11 +137,11 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP,NBDLW), intent(out) :: faerlw3 real(kind=kind_phys), dimension(size(Grid%xlon,1),NSPC1), intent(out) :: aerodp - real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(out) :: clouds1 - real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(out) :: clouds2 - real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(out) :: clouds3 - real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(out) :: clouds4 - real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(out) :: clouds5 + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(inout) :: clouds1 + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(inout) :: clouds2 + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(inout) :: clouds3 + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(inout) :: clouds4 + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(inout) :: clouds5 real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(out) :: clouds6 real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(out) :: clouds7 real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(out) :: clouds8 @@ -142,7 +156,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input integer, intent(out) :: errflg ! Local variables - integer :: me, nfxr, ntrac, ntcw, ntiw, ncld, ntrw, ntsw, ntgl, ncndl + integer :: me, nfxr, ntrac, ntcw, ntiw, ncld, ntrw, ntsw, ntgl, ncndl, ntlnc, ntinc integer :: i, j, k, k1, k2, lsk, lv, n, itop, ibtc, LP1, lla, llb, lya, lyb @@ -154,7 +168,11 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input htswc, htlwc, gcice, grain, grime, htsw0, htlw0, & rhly, tvly,qstl, vvel, clw, ciw, prslk1, tem2da, & cldcov, deltaq, cnvc, cnvw, & - effrl, effri, effrr, effrs + effrl, effri, effrr, effrs, rho, orho + ! for Thompson MP + real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP) :: & + re_cloud, re_ice, re_snow, qv_mp, qc_mp, & + qi_mp, qs_mp, nc_mp, ni_mp real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP+1) :: tem2db ! real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP+1) :: hz @@ -165,6 +183,9 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP,NF_VGAS) :: gasvmr real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP,NBDSW,NF_AESW)::faersw real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP,NBDLW,NF_AELW)::faerlw + + logical :: clduni + real(kind=kind_phys) :: qvs ! !===> ... begin here ! @@ -180,6 +201,8 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input NTRAC = Model%ntrac ! tracers in grrad strip off sphum - start tracer1(2:NTRAC) ntcw = Model%ntcw ntiw = Model%ntiw + ntlnc = Model%ntlnc + ntinc = Model%ntinc ncld = Model%ncld ntrw = Model%ntrw ntsw = Model%ntsw @@ -257,6 +280,9 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input tlyr(i,k1) = Statein%tgrs(i,k2) prslk1(i,k1) = Statein%prslk(i,k2) + rho(i,k1) = plyr(i,k1)/(con_rd*tlyr(i,k1)) + orho(i,k1) = 1.0/rho(i,k1) + !> - Compute relative humidity. es = min( Statein%prsl(i,k2), fpvs( Statein%tgrs(i,k2) ) ) ! fpvs and prsl in pa qs = max( QMIN, eps * es / (Statein%prsl(i,k2) + epsm1*es) ) @@ -273,6 +299,15 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input tracer1(:,k1,j) = max(0.0, Statein%qgrs(:,k2,j)) enddo enddo + if ((Model%do_mynnedmf.or. (imfdeepcnv == imfdeepcnv_gf)) .and. Model%kdt > 1) then + ! for MYNN PBL and GF convective include subgrid clouds into tracer1 + do k = 1, LM + k1 = k + kd + k2 = k + lsk + tracer1(:,k1,ntcw) = max(0.0, qc(:,k2)) + tracer1(:,k1,ntiw) = max(0.0, qi(:,k2)) + enddo + endif ! if (ivflip == 0) then ! input data from toa to sfc do i = 1, IM @@ -552,6 +587,17 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input ccnd(i,k,2) = tracer1(i,k,ntiw) ! ice water ccnd(i,k,3) = tracer1(i,k,ntrw) ! rain water ccnd(i,k,4) = tracer1(i,k,ntsw) + tracer1(i,k,ntgl) ! snow + grapuel + + ! for Thompson MP - prepare variables for calc_effr + if (Model%imp_physics == Model%imp_physics_thompson) then + qvs = Statein%qgrs(i,k2,1) + qv_mp (i,k) = qvs/(1.-qvs) + qc_mp (i,k) = tracer1(i,k,ntcw)/(1.-qvs) + qi_mp (i,k) = tracer1(i,k,ntiw)/(1.-qvs) + qs_mp (i,k) = tracer1(i,k,ntsw)/(1.-qvs) + nc_mp (i,k) = tracer1(i,k,ntlnc)/(1.-qvs) + ni_mp (i,k) = tracer1(i,k,ntinc)/(1.-qvs) + endif enddo enddo endif @@ -562,7 +608,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input enddo enddo enddo - if (Model%imp_physics == 11 ) then + if (Model%imp_physics == Model%imp_physics_gfdl ) then if (.not. Model%lgfdlmprad) then @@ -583,7 +629,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input ! do j=1,Model%ncld ! ccnd(:,:,1) = ccnd(:,:,1) + tracer1(:,1:LMK,ntcw+j-1) ! cloud condensate amount ! enddo - endif + endif ! imp_physics == 11 do k=1,LMK do i=1,IM if (ccnd(i,k,1) < EPSQ ) ccnd(i,k,1) = 0.0 @@ -612,7 +658,29 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input enddo endif elseif (Model%imp_physics == Model%imp_physics_gfdl) then ! GFDL MP - cldcov(1:IM,1+kd:LM+kd) = tracer1(1:IM,1:LM,Model%ntclamt) + IF (Model%do_mynnedmf) THEN + if(Model%kdt == 1) then + ! GFDL cloud fraction + cldcov(1:IM,1+kd:LM+kd) = tracer1(1:IM,1:LM,Model%ntclamt) + else ! kdt > 1 + do k=1,lm + k1 = k + kd + do i=1,im + IF (tracer1(i,k1,ntrw)>1.0e-7 .OR. tracer1(i,k1,ntsw)>1.0e-7) then + ! GFDL cloud fraction + cldcov(i,k1) = tracer1(I,k1,Model%ntclamt) + ELSE + ! MYNN sub-grid cloud fraction + cldcov(i,k1) = clouds1(i,k1) + ENDIF + enddo + enddo + endif + ELSE + ! GFDL cloud fraction + cldcov(1:IM,1+kd:LM+kd) = tracer1(1:IM,1:LM,Model%ntclamt) + ENDIF + if(Model%effr_in) then do k=1,lm k1 = k + kd @@ -634,6 +702,103 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input enddo enddo endif + elseif (Model%imp_physics == Model%imp_physics_thompson) then ! Thompson MP + if(Model%kdt == 1 ) then + do k=1,lm + k1 = k + kd + do i=1,im + effrl(i,k1) = Tbd%phy_f3d(i,k,Model%nleffr) + effri(i,k1) = Tbd%phy_f3d(i,k,Model%nieffr) + effrr(i,k1) = 1000. ! rrain_def=1000. + effrs(i,k1) = Tbd%phy_f3d(i,k,Model%nseffr) + enddo + enddo + else ! kdt>1 + if(Model%do_mynnedmf .or. & + Model%imfdeepcnv == Model%imfdeepcnv_gf ) then + !tgs - take into account sub-grid clouds from GF or MYNN PBL + + ! Compute effective radii for QC and QI with sub-grid clouds + do k=1,lm + do i=1,im + re_cloud(i,k) = 2.49E-6 + re_ice(i,k) = 4.99E-6 + re_snow(i,k) = 999.E-6 + ! make NC consistent with sub-grid clouds + if (qc_mp(i,k)>1.e-12 .and. nc_mp(i,k)<100.) then + nc_mp(i,k) = make_DropletNumber(qc_mp(i,k)*rho(i,k), nwfa(i,k)) * orho(i,k) + endif + if (qi_mp(i,k)>1.e-12 .and. ni_mp(i,k)<100.) then + ni_mp(i,k) = make_IceNumber(qi_mp(i,k)*rho(i,k), tlyr(i,k)) * orho(i,k) + endif + end do + end do + do i = 1, im + call calc_effectRad (tlyr(i,:), plyr(i,:), qv_mp(i,:), qc_mp(i,:), & + nc_mp(i,:), qi_mp(i,:), ni_mp(i,:), qs_mp(i,:), & + re_cloud(i,:), re_ice(i,:), re_snow(i,:), 1, lm ) + end do + do k=1,lm + do i=1,im + re_cloud(i,k) = MAX(2.49, MIN(re_cloud(i,k)*1.e6, 50.)) + re_ice(i,k) = MAX(4.99, MIN(re_ice(i,k)*1.e6, 125.)) + !tgs: clduni has different limits for ice radii: 10.0-150.0 + ! it will raise the low limit from 5 to 10, but the + ! high limit will remain 125. + re_snow(i,k) = MAX(9.99, MIN(re_snow(i,k)*1.e6, 999.)) + end do + end do + if(1==2) then + write(0,'(a,3e16.7)') " before progclduni: re_cloud min/mean/max =", & + minval(re_cloud), & + sum(re_cloud)/real(size(re_cloud)), & + maxval(re_cloud) + write(0,'(a,3e16.7)') " before progclduni: re_ice min/mean/max =", & + minval(re_ice), & + sum(re_ice)/real(size(re_ice)), & + maxval(re_ice) + write(0,'(a,3e16.7)') " before progclduni: clouds3 min/mean/max =", & + minval(clouds3), & + sum(clouds3)/real(size(clouds3)), & + maxval(clouds3) + write(0,'(a,3e16.7)') " before progclduni: clouds5 min/mean/max =", & + minval(clouds5), & + sum(clouds5)/real(size(clouds5)), & + maxval(clouds5) + write(0,'(a,3e16.7)') " before progcld5: phy_f3d cl min/mean/max =", & + minval(Tbd%phy_f3d(:,:,Model%nleffr)), & + sum(Tbd%phy_f3d(:,:,Model%nleffr))/real(size(Tbd%phy_f3d(:,:,Model%nleffr))), & + maxval(Tbd%phy_f3d(:,:,Model%nleffr)) + write(0,'(a,3e16.7)')" before progcld5: phy_f3d ice min/mean/max =", & + minval(Tbd%phy_f3d(:,:,Model%nieffr)), & + sum(Tbd%phy_f3d(:,:,Model%nieffr))/real(size(Tbd%phy_f3d(:,:,Model%nieffr))), & + maxval(Tbd%phy_f3d(:,:,Model%nieffr)) + endif + + do k=1,lm + k1 = k + kd + do i=1,im + !effrl(i,k1) = clouds3 (i,k) ! Tbd%phy_f3d(i,k,Model%nleffr) + !effri(i,k1) = clouds5 (i,k) ! Tbd%phy_f3d(i,k,Model%nieffr) + effrl(i,k1) = re_cloud (i,k) ! Tbd%phy_f3d(i,k,Model%nleffr) + effri(i,k1) = re_ice (i,k) ! Tbd%phy_f3d(i,k,Model%nieffr) + effrr(i,k1) = 1000. ! rrain_def=1000. + effrs(i,k1) = Tbd%phy_f3d(i,k,Model%nseffr) + enddo + enddo + else ! not MYNN or not GF + do k=1,lm + k1 = k + kd + do i=1,im + effrl(i,k1) = Tbd%phy_f3d(i,k,Model%nleffr) + effri(i,k1) = Tbd%phy_f3d(i,k,Model%nieffr) + effrr(i,k1) = 1000. ! rrain_def=1000. + effrs(i,k1) = Tbd%phy_f3d(i,k,Model%nseffr) + enddo + enddo + endif ! MYNN PBL or GF conv + endif ! kdt + else ! neither of the other two cases cldcov = 0.0 endif @@ -748,9 +913,8 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input ! clouds, cldsa, mtopa, mbota, de_lgth) ! --- outputs endif - elseif(Model%imp_physics == 8 .or. Model%imp_physics == 6 .or. & - Model%imp_physics == 15) then - if (Model%kdt == 1 .and. .not.Model%imp_physics == 8) then + elseif(Model%imp_physics == 6 .or. Model%imp_physics == 15) then + if (Model%kdt == 1 ) then Tbd%phy_f3d(:,:,Model%nleffr) = 10. Tbd%phy_f3d(:,:,Model%nieffr) = 50. Tbd%phy_f3d(:,:,Model%nseffr) = 250. @@ -766,6 +930,118 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input Tbd%phy_f3d(:,:,2), Tbd%phy_f3d(:,:,3), & clouds,cldsa,mtopa,mbota, de_lgth) ! --- outputs + + elseif(Model%imp_physics == Model%imp_physics_thompson) then ! Thompson MP + + clduni = .true. + + if(Model%do_mynnedmf .or. & + Model%imfdeepcnv == Model%imfdeepcnv_gf ) then ! MYNN PBL or GF conv + ! MYNN PBL or convective GF + + if (Model%kdt == 1 ) then + ! --- call progcld5 to get Xu-Randall total cloud cover (clouds(:,1:LMK,1)) at + ! --- initial time step, it takes into account subgrid PBL + ! --- clouds + call progcld5 (plyr,plvl,tlyr,qlyr,qstl,rhly,tracer1, & ! --- inputs + Grid%xlat,Grid%xlon,Sfcprop%slmsk,dz,delp, & + ntrac-1, ntcw-1,ntiw-1,ntrw-1, & + ntsw-1,ntgl-1, & + im, lmk, lmp, Model%uni_cld, & + Model%lmfshal,Model%lmfdeep2, & + cldcov(:,1:LMK),Tbd%phy_f3d(:,:,Model%nleffr), & + Tbd%phy_f3d(:,:,Model%nieffr), & + Tbd%phy_f3d(:,:,Model%nseffr), & + clouds,cldsa,mtopa,mbota, de_lgth) ! --- outputs + if ( clduni) then + ! use progclduni for interaction with radiation, + ! overwrites 'clouds' from progcld5 + call progclduni (plyr, plvl, tlyr, tvly, ccnd, ncndl, & ! --- inputs + Grid%xlat, Grid%xlon, Sfcprop%slmsk, dz,delp, & + IM, LMK, LMP, clouds(:,1:LMK,1), & + effrl, effri, effrr, effrs, Model%effr_in , & + clouds, cldsa, mtopa, mbota, de_lgth) ! --- outputs + endif + + else ! kdt > 1 + + do k=1,lm + k1 = k + kd + do i=1,im + Tbd%phy_f3d(i,k,Model%nleffr) = effrl(i,k1) + Tbd%phy_f3d(i,k,Model%nieffr) = effri(i,k1) + Tbd%phy_f3d(i,k,Model%nseffr) = effrs(i,k1) + enddo + enddo + + ! --- call progcld5 to get Xu-Randall total cloud cover (clouds(:,1:LMK,1)) + ! tgs: a short subroutine could be made of progcld5 only to + ! compute total cloud fraction. + call progcld5 (plyr,plvl,tlyr,qlyr,qstl,rhly,tracer1, & ! --- inputs + Grid%xlat,Grid%xlon,Sfcprop%slmsk,dz,delp, & + ntrac-1, ntcw-1,ntiw-1,ntrw-1, & + ntsw-1,ntgl-1, & + im, lmk, lmp, Model%uni_cld, & + Model%lmfshal,Model%lmfdeep2, & + cldcov(:,1:LMK),Tbd%phy_f3d(:,:,Model%nleffr), & + Tbd%phy_f3d(:,:,Model%nieffr), & + Tbd%phy_f3d(:,:,Model%nseffr), & + clouds,cldsa,mtopa,mbota, de_lgth) ! --- outputs + + do k=1,lmk + do i=1,im + !IF (tracer1(i,k,ntrw) > 1.0e-7 .OR. tracer1(i,k,ntsw) > 1.0e-7) then + ! ! Xu-Randall cloud fraction computed in progcld5 + ! cldcov(i,k) = clouds(i,k,1) + ! clouds(i,k,1) = clouds(i,k,1) + !ELSE + ! MYNN sub-grid cloud fraction + !tgs - let's use only PBL cloud fraction + cldcov(i,k) = clouds1(i,k) + clouds(i,k,1) = clouds1(i,k) + !ENDIF + enddo + enddo + if( .not. clduni) then + ! --- call progcld5 for interaction with the radiation with setting + ! --- uni_cld=.true. to keep precomputed cloud + ! --- fraction + call progcld5 (plyr,plvl,tlyr,qlyr,qstl,rhly,tracer1, & ! --- inputs + Grid%xlat,Grid%xlon,Sfcprop%slmsk,dz,delp, & + ntrac-1, ntcw-1,ntiw-1,ntrw-1, & + ntsw-1,ntgl-1, & + im, lmk, lmp, .true., & ! Model%uni_cld, + Model%lmfshal,Model%lmfdeep2, & + cldcov(:,1:LMK),Tbd%phy_f3d(:,:,Model%nleffr), & + Tbd%phy_f3d(:,:,Model%nieffr), & + Tbd%phy_f3d(:,:,Model%nseffr), & + clouds,cldsa,mtopa,mbota, de_lgth) ! --- outputs + + else ! clduni + ! --- use clduni as with the GFDL microphysics. + ! --- make sure that effr_in=.true. in the input.nml! + call progclduni (plyr, plvl, tlyr, tvly, ccnd, ncndl, & ! --- inputs + Grid%xlat, Grid%xlon, Sfcprop%slmsk, dz,delp, & + IM, LMK, LMP, clouds(:,1:LMK,1), & + effrl, effri, effrr, effrs, Model%effr_in , & + clouds, cldsa, mtopa, mbota, de_lgth) ! --- outputs + endif ! clduni + + endif ! kdt + + else + ! MYNN PBL or GF convective are not used + call progcld5 (plyr,plvl,tlyr,qlyr,qstl,rhly,tracer1, & ! --- inputs + Grid%xlat,Grid%xlon,Sfcprop%slmsk,dz,delp, & + ntrac-1, ntcw-1,ntiw-1,ntrw-1, & + ntsw-1,ntgl-1, & + im, lmk, lmp, Model%uni_cld, & + Model%lmfshal,Model%lmfdeep2, & + cldcov(:,1:LMK),Tbd%phy_f3d(:,:,1), & + Tbd%phy_f3d(:,:,2), Tbd%phy_f3d(:,:,3), & + clouds,cldsa,mtopa,mbota, de_lgth) ! --- outputs + endif ! MYNN PBL or GF + endif ! end if_imp_physics ! endif ! end_if_ntcw diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index 7b40e2c1d..423a50ff0 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -270,6 +270,67 @@ kind = kind_phys intent = out optional = F +[qc] + standard_name = cloud_condensed_water_mixing_ratio + long_name = moist (dry+vapor, no condensates) mixing ratio of cloud water (condensate) + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qi] + standard_name = ice_water_mixing_ratio + long_name = moist (dry+vapor, no condensates) mixing ratio of ice water + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[nc] + standard_name = cloud_droplet_number_concentration + long_name = cloud droplet number concentration + units = kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = T +[ni] + standard_name = ice_number_concentration + long_name = ice number concentration + units = kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[nwfa] + standard_name = water_friendly_aerosol_number_concentration + long_name = number concentration of water-friendly aerosols + units = kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = T +[imfdeepcnv] + standard_name = flag_for_mass_flux_deep_convection_scheme + long_name = flag for mass-flux deep convection scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F +[imfdeepcnv_gf] + standard_name = flag_for_gf_deep_convection_scheme + long_name = flag for Grell-Freitas deep convection scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F [gasvmr_co2] standard_name = volume_mixing_ratio_co2 long_name = CO2 volume mixing ratio diff --git a/physics/cu_gf_driver.F90 b/physics/cu_gf_driver.F90 index 53e26fb46..1d21a7f4e 100644 --- a/physics/cu_gf_driver.F90 +++ b/physics/cu_gf_driver.F90 @@ -75,7 +75,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,ix,km,dt,cactiv, & hbot,htop,kcnv,xland,hfx2,qfx2,cliw,clcw, & pbl,ud_mf,dd_mf,dt_mf,cnvw_moist,cnvc,imfshalcnv, & nwfa,con_rd,gq0,ntinc,ntlnc,imp_physics,imp_physics_thompson, & - errmsg,errflg) + qci_conv,errmsg,errflg) !------------------------------------------------------------- implicit none integer, parameter :: maxiens=1 @@ -98,8 +98,9 @@ subroutine cu_gf_driver_run(ntracer,garea,im,ix,km,dt,cactiv, & integer :: its,ite, jts,jte, kts,kte integer, intent(in ) :: im,ix,km,ntracer - real(kind=kind_phys), dimension( ix , km ), intent(in ) :: forcet,forceqv_spechum,w,phil - real(kind=kind_phys), dimension( ix , km ), intent(inout ) :: t,us,vs + real(kind=kind_phys), dimension( ix , km ), intent(in ) :: forcet,forceqv_spechum,w,phil + real(kind=kind_phys), dimension( ix , km ), intent(inout ) :: t,us,vs + real(kind=kind_phys), dimension( ix , km ), intent(inout ) :: qci_conv real(kind=kind_phys), dimension( ix ) :: rand_mom,rand_vmas real(kind=kind_phys), dimension( ix,4 ) :: rand_clos real(kind=kind_phys), dimension( ix , km, 11 ) :: gdc,gdc2 @@ -751,6 +752,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,ix,km,dt,cactiv, & gdc(i,k,1)= max(0.,tun_rad_shall(i)*cupclws(i,k)*cutens(i)) ! my mod gdc2(i,k,1)=max(0.,tun_rad_deep(i)*(cupclwm(i,k)*cutenm(i)+cupclw(i,k)*cuten(i))) + qci_conv(i,k)=gdc2(i,k,1) gdc(i,k,2)=(outt(i,k))*86400. gdc(i,k,3)=(outtm(i,k))*86400. gdc(i,k,4)=(outts(i,k))*86400. diff --git a/physics/cu_gf_driver.meta b/physics/cu_gf_driver.meta index d3687a352..3966c1eba 100644 --- a/physics/cu_gf_driver.meta +++ b/physics/cu_gf_driver.meta @@ -417,6 +417,15 @@ type = integer intent = in optional = F +[qci_conv] + standard_name = convective_cloud_condesate_after_rainout + long_name = convective cloud condesate after rainout + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/module_MYNNrad_pre.F90 b/physics/module_MYNNrad_pre.F90 index 95dc95445..54c47f681 100644 --- a/physics/module_MYNNrad_pre.F90 +++ b/physics/module_MYNNrad_pre.F90 @@ -85,6 +85,13 @@ SUBROUTINE mynnrad_pre_run( & qi_save(i,k) = qi(i,k) clouds1(i,k) = CLDFRA_BL(i,k) + !IF (qr(i,k) > 1.0e-7 .OR. qs(i,k) > 1.0e-7) then ! .OR. & + !(Model%imfdeepcnv == Model%imfdeepcnv_gf .AND. qci_conv(i,k)>1.0e-7)) THEN + !Keep Xu-RandalL clouds fraction + !ELSE + ! clouds1(i,k) = CLDFRA_BL(i,k) + !ENDIF + IF (qc(i,k) < 1.E-6 .AND. qi(i,k) < 1.E-8 .AND. CLDFRA_BL(i,k)>0.001) THEN !Partition the BL clouds into water & ice according to a linear !approximation of Hobbs et al. (1974). This allows us to only use diff --git a/physics/module_SGSCloud_RadPost.F90 b/physics/module_SGSCloud_RadPost.F90 new file mode 100644 index 000000000..810c3bcd3 --- /dev/null +++ b/physics/module_SGSCloud_RadPost.F90 @@ -0,0 +1,75 @@ +!> \file module_SGSCloud_RadPost.F90 +!! Contains the post (interstitial) work after the call to the radiation schemes: +!! 1) Restores the original qc & qi + + MODULE sgscloud_radpost + + contains + + subroutine sgscloud_radpost_init () + end subroutine sgscloud_radpost_init + + subroutine sgscloud_radpost_finalize () + end subroutine sgscloud_radpost_finalize + +!>\defgroup sgscloud_radpost GSD sgscloud_radpost_run Module +!>\ingroup gsd_mynn_edmf +!! This interstitial code restores the original resolved-scale clouds (qc and qi). +#if 0 +!! \section arg_table_sgscloud_radpost_run Argument Table +!! \htmlinclude sgscloud_radpost_run.html +!! +#endif +SUBROUTINE sgscloud_radpost_run( & + & ix,im,levs, & + & flag_init,flag_restart, & + & qc,qi, & + & qc_save, qi_save, & + & errmsg, errflg ) + +! should be moved to inside the mynn: + use machine , only : kind_phys + +!------------------------------------------------------------------- + implicit none +!------------------------------------------------------------------- + + integer, intent(in) :: ix, im, levs + logical, intent(in) :: flag_init, flag_restart + real(kind=kind_phys), dimension(im,levs), intent(out) :: qc, qi + real(kind=kind_phys), dimension(im,levs), intent(in) :: qc_save, qi_save + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + ! Local variable + integer :: i, k + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + !write(0,*)"==============================================" + !write(0,*)"in mynn rad post" + + if (flag_init .and. (.not. flag_restart)) then + !write (0,*) 'Skip MYNNrad_post flag_init = ', flag_init + return + endif + + ! Add subgrid cloud information: + do k = 1, levs + do i = 1, im + + qc(i,k) = qc_save(i,k) + qi(i,k) = qi_save(i,k) + + enddo + enddo + + ! print*,"===Finished restoring the resolved-scale clouds" + ! print*,"qc_save:",qc_save(1,1)," qc:",qc(1,1) + + END SUBROUTINE sgscloud_radpost_run + +!###================================================================= + +END MODULE sgscloud_radpost diff --git a/physics/module_SGSCloud_RadPost.meta b/physics/module_SGSCloud_RadPost.meta new file mode 100644 index 000000000..0318aa231 --- /dev/null +++ b/physics/module_SGSCloud_RadPost.meta @@ -0,0 +1,96 @@ +[ccpp-arg-table] + name = sgscloud_radpost_run + type = scheme +[ix] + standard_name = horizontal_dimension + long_name = horizontal dimension + units = count + dimensions = () + type = integer + intent = in + optional = F +[im] + standard_name = horizontal_loop_extent + long_name = horizontal loop extent + units = count + dimensions = () + type = integer + intent = in + optional = F +[levs] + standard_name = vertical_dimension + long_name = vertical layer dimension + units = count + dimensions = () + type = integer + intent = in + optional = F +[flag_init] + standard_name = flag_for_first_time_step + long_name = flag signaling first time step for time integration loop + units = flag + dimensions = () + type = logical + intent = in + optional = F +[flag_restart] + standard_name = flag_for_restart + long_name = flag for restart (warmstart) or coldstart + units = flag + dimensions = () + type = logical + intent = in + optional = F +[qc] + standard_name = cloud_condensed_water_mixing_ratio + long_name = moist (dry+vapor, no condensates) mixing ratio of cloud water (condensate) + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[qi] + standard_name = ice_water_mixing_ratio + long_name = moist (dry+vapor, no condensates) mixing ratio of ice water + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[qc_save] + standard_name = cloud_condensed_water_mixing_ratio_save + long_name = moist (dry+vapor, no condensates) mixing ratio of cloud water (condensate) before entering a physics scheme + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[qi_save] + standard_name = ice_water_mixing_ratio_save + long_name = moist (dry+vapor, no condensates) mixing ratio of ice water before entering a physics scheme + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out + optional = F +[errflg] + standard_name = ccpp_error_flag + long_name = error flag for error handling in CCPP + units = flag + dimensions = () + type = integer + intent = out + optional = F diff --git a/physics/module_SGSCloud_RadPre.F90 b/physics/module_SGSCloud_RadPre.F90 new file mode 100644 index 000000000..91617b06b --- /dev/null +++ b/physics/module_SGSCloud_RadPre.F90 @@ -0,0 +1,211 @@ +!>\file module_SGSCloud_RadPre.F90 +!! Contains the preliminary (interstitial) work to the call to the radiation schemes: +!! 1) Backs up the original qc & qi +!! 2) Adds the partioning of convective condensate into liqice/ice for effective radii +!! 3) Adds the subgrid clouds mixing ratio and cloud fraction to the original qc, qi and cloud fraction coming from the microphysics scheme. +!! 4) Recompute the diagnostic high, mid, low, total and bl clouds to be consistent with radiation + + MODULE sgscloud_radpre + + contains + + subroutine sgscloud_radpre_init () + end subroutine sgscloud_radpre_init + + subroutine sgscloud_radpre_finalize () + end subroutine sgscloud_radpre_finalize + +!> \defgroup sgsrad_group GSD sgscloud_radpre_run Module +!> \ingroup sgscloud_radpre +!! This interstitial code adds the subgrid clouds to the resolved-scale clouds if there is no resolved-scale clouds in that particular grid box. +!> \section arg_table_sgscloud_radpre_run Argument Table +!! \htmlinclude sgscloud_radpre_run.html +!! +!! +!! cloud array description: ! +!! clouds(:,:,1) - layer total cloud fraction ! +!! clouds(:,:,2) - layer cloud liq water path ! +!! clouds(:,:,3) - mean effective radius for liquid cloud ! +!! clouds(:,:,4) - layer cloud ice water path ! +!! clouds(:,:,5) - mean effective radius for ice cloud ! +!! +!>\section sgscloud_radpre GSD SGS Scheme General Algorithm +!> @{ +SUBROUTINE sgscloud_radpre_run( & + & ix,im,levs, & + & flag_init,flag_restart, & + & do_mynnedmf, & + & qc, qi, T3D, & + & qr, qs, & + & qci_conv, & + & imfdeepcnv, & + & qc_save, qi_save, & + & qc_bl,cldfra_bl, & + & delp,clouds1,clouds2,clouds3, & + & clouds4,clouds5,slmsk, & + & nlay, plyr, xlat, dz,de_lgth, & + & cldsa,mtopa,mbota, & + & errmsg, errflg ) + +! should be moved to inside the mynn: + use machine , only : kind_phys + ! DH* TODO - input argument, not constant + use physcons, only : con_g, con_pi + use module_radiation_clouds, only : gethml + +!------------------------------------------------------------------- + implicit none +!------------------------------------------------------------------- + ! Interface variables + real (kind=kind_phys), parameter :: gfac=1.0e5/con_g + integer, intent(in) :: ix, im, levs, imfdeepcnv, nlay + logical, intent(in) :: flag_init, flag_restart, do_mynnedmf + real(kind=kind_phys), dimension(im,levs), intent(inout) :: qc, qi + real(kind=kind_phys), dimension(im,levs), intent(inout) :: qr, qs + real(kind=kind_phys), dimension(im,levs), intent(inout) :: qci_conv + real(kind=kind_phys), dimension(im,levs), intent(in) :: T3D,delp + real(kind=kind_phys), dimension(im,levs), intent(inout) :: & + & clouds1,clouds2,clouds3,clouds4,clouds5 + real(kind=kind_phys), dimension(im,levs), intent(out) :: qc_save, qi_save + real(kind=kind_phys), dimension(im,levs), intent(in) :: qc_bl, cldfra_bl + ! DH* TODO add intent() information for delp,clouds1,clouds2,clouds3,clouds4,clouds5 + real(kind=kind_phys), dimension(im), intent(in) :: slmsk, xlat, de_lgth + real(kind=kind_phys), dimension(im,nlay), intent(in) :: plyr, dz + real(kind=kind_phys), dimension(im,5), intent(out) :: cldsa + integer, dimension(im,3), intent(out) :: mbota, mtopa + character(len=*), intent(out) :: errmsg + integer, intent(out) :: errflg + ! Local variables + ! pressure limits of cloud domain interfaces (low,mid,high) in mb (0.1kPa) + real (kind=kind_phys) :: ptop1(im,3+1) !< pressure limits of cloud domain interfaces + real (kind=kind_phys) :: ptopc(3+1,2 ) !< pressure limits of cloud domain interfaces + !! (low, mid, high) in mb (0.1kPa) + data ptopc / 1050., 650., 400., 0.0, 1050., 750., 500., 0.0 / + real(kind=kind_phys), dimension(im,nlay) :: cldcnv + real(kind=kind_phys), dimension(im) :: rxlat + real (kind=kind_phys):: Tc, iwc, tem1 + integer :: i, k, id + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 + + !write(0,*)"==============================================" + !write(0,*)"in mynn rad pre" + + if (flag_init .and. (.not. flag_restart)) then + !write (0,*) 'Skip MYNNrad_pre flag_init = ', flag_init + return + endif + ! Back-up microphysics cloud information: + do k = 1, levs + do i = 1, im + qc_save(i,k) = qc(i,k) + qi_save(i,k) = qi(i,k) + end do + end do + + ! add boundary layer clouds + IF (do_mynnedmf == .true.) THEN + do k = 1, levs + do i = 1, im + + clouds1(i,k) = CLDFRA_BL(i,k) + + !IF( qr(i,k) > 1.0e-7 .OR. qs(i,k) > 1.0e-7.or.qci_conv(i,k)>1.0e-7)THEN + !Keep Xu-RandalL clouds fraction - do not overwrite + !ELSE + ! clouds1(i,k) = CLDFRA_BL(i,k) + !ENDIF + + IF (qc(i,k) < 1.E-6 .AND. qi(i,k) < 1.E-8 .AND. CLDFRA_BL(i,k)>0.001) THEN + !Partition the BL clouds into water & ice according to a linear + !approximation of Hobbs et al. (1974). This allows us to only use + !one 3D array for both cloud water & ice. +! Wice = 1. - MIN(1., MAX(0., (t(i,k)-254.)/15.)) +! Wh2o = 1. - Wice + !clouds1(i,k)=MAX(clouds1(i,k),CLDFRA_BL(i,k)) + !clouds1(i,k)=MAX(0.0,MIN(1.0,clouds1(i,k))) + qc(i,k) = QC_BL(i,k)*(MIN(1., MAX(0., (T3D(i,k)-244.)/25.)))*CLDFRA_BL(i,k) + qi(i,k) = QC_BL(i,k)*(1. - MIN(1., MAX(0., (T3D(i,k)-244.)/25.)))*CLDFRA_BL(i,k) + + Tc = T3D(i,k) - 273.15 + !iwc = qi(i,k)*1.0e6*rho(i,k) + + IF (nint(slmsk(i)) == 1) then !land + IF(qc(i,k)>1.E-8)clouds3(i,k)=5.4 !eff radius cloud water (microns) + !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos) + IF(qi(i,k)>1.E-8)clouds5(i,k)=MAX(173.45 + 2.14*Tc, 20.) + ELSE + !eff radius cloud water (microns), from Miles et al. + IF(qc(i,k)>1.E-8)clouds3(i,k)=9.6 + !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 6b) + IF(qi(i,k)>1.E-8)clouds5(i,k)=MAX(173.45 + 2.14*Tc, 20.) + !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 8b) + !IF(qi(i,k)>1.E-8)clouds5(i,k)=MAX(139.7 + 1.76*Tc + 13.49*LOG(iwc), 20.) + ENDIF + !calculate water and ice paths for additional BL clouds + clouds2(i,k) = max(0.0, qc(i,k) * gfac * delp(i,k)) + clouds4(i,k) = max(0.0, qi(i,k) * gfac * delp(i,k)) + ENDIF + enddo + enddo + ENDIF ! do_mynnedmf + + ! add convective clouds + IF (imfdeepcnv == 3) THEN + do k = 1, levs + do i = 1, im + IF ( qci_conv(i,k) > 0.) THEN + !IF (qc(i,k) < 1.E-6 .AND. qi(i,k) < 1.E-8 .AND. qci_conv(i,k) > 0.) THEN + !Partition the convective clouds into water & ice according to a linear + qc(i,k) = qc(i,k)+qci_conv(i,k)*(MIN(1., MAX(0., (T3D(i,k)-244.)/25.))) + qi(i,k) = qi(i,k)+qci_conv(i,k)*(1. - MIN(1., MAX(0., (T3D(i,k)-244.)/25.))) + + Tc = T3D(i,k) - 273.15 + + IF (nint(slmsk(i)) == 1) then !land + IF(qc(i,k)>1.E-8)clouds3(i,k)=5.4 !eff radius cloud water (microns) + !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos) + IF(qi(i,k)>1.E-8)clouds5(i,k)=MAX(173.45 + 2.14*Tc, 20.) + ELSE + !eff radius cloud water (microns), from Miles et al. + IF(qc(i,k)>1.E-8)clouds3(i,k)=9.6 + !eff radius cloud ice (microns), from Mishra et al. (2014, JGR Atmos, fig 6b) + IF(qi(i,k)>1.E-8)clouds5(i,k)=MAX(173.45 + 2.14*Tc, 20.) + ENDIF + ENDIF + enddo + enddo + ENDIF +!> - Compute SFC/low/middle/high cloud top pressure for each cloud domain for given latitude. + + do i =1, im + 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, im + ptop1(i,id) = ptopc(id,1) + tem1*max( 0.0, 4.0*rxlat(i)-1.0 ) + enddo + enddo + + cldcnv = 0. + +!> - Recompute the diagnostic high, mid, low, total and bl cloud fraction + call gethml & +! --- inputs: + & ( plyr, ptop1, clouds1, cldcnv, dz, de_lgth, im, nlay, & +! --- outputs: + & cldsa, mtopa, mbota) + + !print*,"===Finished adding subgrid clouds to the resolved-scale clouds" + !print*,"qc_save:",qc_save(1,1)," qi_save:",qi_save(1,1) + + END SUBROUTINE sgscloud_radpre_run + +!###================================================================= + +END MODULE sgscloud_radpre diff --git a/physics/module_SGSCloud_RadPre.meta b/physics/module_SGSCloud_RadPre.meta new file mode 100644 index 000000000..349d37885 --- /dev/null +++ b/physics/module_SGSCloud_RadPre.meta @@ -0,0 +1,308 @@ +[ccpp-arg-table] + name = sgscloud_radpre_init + type = scheme + +######################################################################## +[ccpp-arg-table] + name = sgscloud_radpre_finalize + type = scheme + +######################################################################## +[ccpp-arg-table] + name = sgscloud_radpre_run + type = scheme +[ix] + standard_name = horizontal_dimension + long_name = horizontal dimension + units = count + dimensions = () + type = integer + intent = in + optional = F +[im] + standard_name = horizontal_loop_extent + long_name = horizontal loop extent + units = count + dimensions = () + type = integer + intent = in + optional = F +[levs] + standard_name = vertical_dimension + long_name = vertical layer dimension + units = count + dimensions = () + type = integer + intent = in + optional = F +[flag_init] + standard_name = flag_for_first_time_step + long_name = flag signaling first time step for time integration loop + units = flag + dimensions = () + type = logical + intent = in + optional = F +[flag_restart] + standard_name = flag_for_restart + long_name = flag for restart (warmstart) or coldstart + units = flag + dimensions = () + type = logical + intent = in + optional = F +[qc] + standard_name = cloud_condensed_water_mixing_ratio + long_name = moist (dry+vapor, no condensates) mixing ratio of cloud water (condensate) + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qi] + standard_name = ice_water_mixing_ratio + long_name = moist (dry+vapor, no condensates) mixing ratio of ice water + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[T3D] + standard_name = air_temperature + long_name = layer mean air temperature + units = K + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[qr] + standard_name = rain_water_mixing_ratio + long_name = moist (dry+vapor, no condensates) mixing ratio of rain water + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qs] + standard_name = snow_water_mixing_ratio + long_name = moist (dry+vapor, no condensates) mixing ratio of snow water + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qci_conv] + standard_name = convective_cloud_condesate_after_rainout + long_name = convective cloud condesate after rainout + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[imfdeepcnv] + standard_name = flag_for_mass_flux_deep_convection_scheme + long_name = flag for mass-flux deep convection scheme + units = flag + dimensions = () + type = integer + intent = in + optional = F +[qc_save] + standard_name = cloud_condensed_water_mixing_ratio_save + long_name = moist (dry+vapor, no condensates) mixing ratio of cloud water (condensate) before entering a physics scheme + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[qi_save] + standard_name = ice_water_mixing_ratio_save + long_name = moist (dry+vapor, no condensates) mixing ratio of ice water before entering a physics scheme + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = out + optional = F +[QC_BL] + standard_name = subgrid_cloud_mixing_ratio_pbl + long_name = subgrid cloud cloud mixing ratio from PBL scheme + units = kg kg-1 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[CLDFRA_BL] + standard_name = subgrid_cloud_fraction_pbl + long_name = subgrid cloud fraction from PBL scheme + units = frac + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[delp] + standard_name = layer_pressure_thickness_for_radiation + long_name = layer pressure thickness on radiation levels + units = hPa + dimensions = (horizontal_dimension,adjusted_vertical_layer_dimension_for_radiation) + type = real + kind = kind_phys + intent = out + optional = F +[clouds1] + standard_name = total_cloud_fraction + long_name = layer total cloud fraction + units = frac + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[clouds2] + standard_name = cloud_liquid_water_path + long_name = layer cloud liquid water path + units = g m-2 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[clouds3] + standard_name = mean_effective_radius_for_liquid_cloud + long_name = mean effective radius for liquid cloud + units = micron + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[clouds4] + standard_name = cloud_ice_water_path + long_name = layer cloud ice water path + units = g m-2 + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[clouds5] + standard_name = mean_effective_radius_for_ice_cloud + long_name = mean effective radius for ice cloud + units = micron + dimensions = (horizontal_dimension,vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[slmsk] + standard_name = sea_land_ice_mask_real + long_name = landmask: sea/land/ice=0/1/2 + units = flag + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[nlay] + standard_name = adjusted_vertical_layer_dimension_for_radiation + long_name = number of vertical layers for radiation + units = count + dimensions = () + type = integer + intent = in + optional = F +[plyr] + standard_name = air_pressure_at_layer_for_radiation_in_hPa + long_name = air pressure at vertical layer for radiation calculation + units = hPa + dimensions = (horizontal_dimension,adjusted_vertical_layer_dimension_for_radiation) + type = real + kind = kind_phys + intent = in + optional = F +[xlat] + standard_name = latitude + long_name = grid latitude in radians + units = radians + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[dz] + standard_name = layer_thickness_for_radiation + long_name = layer thickness on radiation levels + units = km + dimensions = (horizontal_dimension,adjusted_vertical_layer_dimension_for_radiation) + type = real + kind = kind_phys + intent = in + optional = F +[de_lgth] + standard_name = cloud_decorrelation_length + long_name = cloud decorrelation length + units = km + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[cldsa] + standard_name = cloud_area_fraction_for_radiation + long_name = fraction of clouds for low, middle,high, total and BL + units = frac + dimensions = (horizontal_dimension,5) + type = real + kind = kind_phys + intent = out + optional = F +[mtopa] + standard_name = model_layer_number_at_cloud_top + long_name = vertical indices for low, middle and high cloud tops + units = index + dimensions = (horizontal_dimension,3) + type = integer + intent = out + optional = F +[mbota] + standard_name = model_layer_number_at_cloud_base + long_name = vertical indices for low, middle and high cloud bases + units = index + dimensions = (horizontal_dimension,3) + type = integer + intent = out + optional = F +[do_mynnedmf] + standard_name = do_mynnedmf + long_name = flag to activate MYNN-EDMF + units = flag + dimensions = () + type = logical + intent = in + optional = F +[errmsg] + standard_name = ccpp_error_message + long_name = error message for error handling in CCPP + units = none + dimensions = () + type = character + kind = len=* + intent = out + optional = F +[errflg] + standard_name = ccpp_error_flag + long_name = error flag for error handling in CCPP + units = flag + dimensions = () + type = integer + intent = out + optional = F diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 49b394fe1..8e5c099aa 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -243,7 +243,7 @@ module module_radiation_clouds integer :: iovr = 1 !< maximum-random cloud overlapping method public progcld1, progcld2, progcld3, progcld4, progclduni, & - & cld_init, progcld5, progcld4o + & cld_init, progcld5, progcld4o, gethml ! ================= @@ -2468,13 +2468,13 @@ subroutine progcld5 & ! !===> ... begin here ! - do nf=1,nf_clds - do k=1,nlay - do i=1,ix - clouds(i,k,nf) = 0.0 - enddo - enddo - enddo + !do nf=1,nf_clds + ! do k=1,nlay + ! do i=1,ix + ! clouds(i,k,nf) = 0.0 + ! enddo + ! enddo + !enddo ! clouds(:,:,:) = 0.0 do k = 1, NLAY @@ -2514,7 +2514,8 @@ subroutine progcld5 & do k = 1, NLAY do i = 1, IX - clwf(i,k) = clw(i,k,ntcw) + clw(i,k,ntiw) + clw(i,k,ntsw) + clwf(i,k) = clw(i,k,ntcw) + clw(i,k,ntiw) + clw(i,k,ntsw) & + & + clw(i,k,ntrw) + clw(i,k,ntgl) enddo enddo !> - Find top pressure for each cloud domain for given latitude. @@ -2558,30 +2559,30 @@ subroutine progcld5 & !> - 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) + !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 + !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 ) + ! onemrh= max( 1.e-10, 1.0-rhly(i,k) ) + ! clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) - tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0) - tem1 = 2000.0 / tem1 + ! 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)) ) + ! 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 + ! 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) @@ -2592,11 +2593,11 @@ subroutine progcld5 & clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) ! tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) !jhan - if (lmfdeep2) then - tem1 = xrc3 / tem1 - else + !if (lmfdeep2) then + ! tem1 = xrc3 / tem1 + !else tem1 = 100.0 / tem1 - endif + !endif ! value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) tem2 = sqrt( sqrt(rhly(i,k)) ) @@ -2605,14 +2606,14 @@ subroutine progcld5 & endif enddo enddo - endif + !endif endif ! if (uni_cld) then do k = 1, NLAY do i = 1, IX if (cldtot(i,k) < climit) then - cldtot(i,k) = 0.0 + !cldtot(i,k) = 0.0 cwp(i,k) = 0.0 cip(i,k) = 0.0 crp(i,k) = 0.0 From 63303d37527213f1a786ec1832494074c2f74468 Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Tue, 4 Feb 2020 18:33:43 +0000 Subject: [PATCH 02/12] Several changes in the comments. --- physics/GFS_rrtmg_pre.F90 | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 6b5382e65..950ea3d5d 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -733,6 +733,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input endif end do end do + ! Call Thompson's subroutine to compoute effective radii do i = 1, im call calc_effectRad (tlyr(i,:), plyr(i,:), qv_mp(i,:), qc_mp(i,:), & nc_mp(i,:), qi_mp(i,:), ni_mp(i,:), qs_mp(i,:), & @@ -975,8 +976,8 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input enddo ! --- call progcld5 to get Xu-Randall total cloud cover (clouds(:,1:LMK,1)) - ! tgs: a short subroutine could be made of progcld5 only to - ! compute total cloud fraction. + ! tgs: a short subroutine could be made of progcld5 to + ! compute only total cloud fraction. call progcld5 (plyr,plvl,tlyr,qlyr,qstl,rhly,tracer1, & ! --- inputs Grid%xlat,Grid%xlon,Sfcprop%slmsk,dz,delp, & ntrac-1, ntcw-1,ntiw-1,ntrw-1, & @@ -988,6 +989,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input Tbd%phy_f3d(:,:,Model%nseffr), & clouds,cldsa,mtopa,mbota, de_lgth) ! --- outputs + !tgs - let's use the PBL cloud fraction do k=1,lmk do i=1,im !IF (tracer1(i,k,ntrw) > 1.0e-7 .OR. tracer1(i,k,ntsw) > 1.0e-7) then @@ -996,7 +998,6 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input ! clouds(i,k,1) = clouds(i,k,1) !ELSE ! MYNN sub-grid cloud fraction - !tgs - let's use only PBL cloud fraction cldcov(i,k) = clouds1(i,k) clouds(i,k,1) = clouds1(i,k) !ENDIF From 4ae3591c2ab5c3f8ad912028e0385ffda0655433 Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Wed, 5 Feb 2020 22:06:40 +0000 Subject: [PATCH 03/12] 1.The unnecessary arays NI and NC are removed. 2. Bug fix for the case when GF scheme is used without MYNN. In this case always use Xu-Randall cloud fraction. --- physics/GFS_rrtmg_pre.F90 | 41 ++++++++++++++++++++++---------------- physics/GFS_rrtmg_pre.meta | 18 ----------------- 2 files changed, 24 insertions(+), 35 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 950ea3d5d..351862cf5 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -20,7 +20,7 @@ end subroutine GFS_rrtmg_pre_init ! in the CCPP version - they are defined in the interstitial_create routine subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input Tbd, Cldprop, Coupling, & - Radtend, qc, qi, nc, ni, nwfa, & ! input/output + Radtend, qc, qi, nwfa, & ! input/output imfdeepcnv, imfdeepcnv_gf, & f_ice, f_rain, f_rimef, flgmin, cwm, & ! F-A mp scheme only lm, im, lmk, lmp, & ! input @@ -89,8 +89,6 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: qc real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: qi - real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: nc - real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: ni real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: nwfa @@ -989,20 +987,29 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input Tbd%phy_f3d(:,:,Model%nseffr), & clouds,cldsa,mtopa,mbota, de_lgth) ! --- outputs - !tgs - let's use the PBL cloud fraction - do k=1,lmk - do i=1,im - !IF (tracer1(i,k,ntrw) > 1.0e-7 .OR. tracer1(i,k,ntsw) > 1.0e-7) then - ! ! Xu-Randall cloud fraction computed in progcld5 - ! cldcov(i,k) = clouds(i,k,1) - ! clouds(i,k,1) = clouds(i,k,1) - !ELSE - ! MYNN sub-grid cloud fraction - cldcov(i,k) = clouds1(i,k) - clouds(i,k,1) = clouds1(i,k) - !ENDIF - enddo - enddo + if(Model%do_mynnedmf) then + !tgs - let's use the PBL cloud fraction for now + do k=1,lmk + do i=1,im + !IF (tracer1(i,k,ntrw) > 1.0e-7 .OR. tracer1(i,k,ntsw) > 1.0e-7) then + ! ! Xu-Randall cloud fraction computed in progcld5 + ! cldcov(i,k) = clouds(i,k,1) + !ELSE + ! MYNN sub-grid cloud fraction + cldcov(i,k) = clouds1(i,k) + clouds(i,k,1) = clouds1(i,k) + !ENDIF + enddo + enddo + elseif (Model%imfdeepcnv == Model%imfdeepcnv_gf ) then ! GF conv + do k=1,lmk + do i=1,im + ! Xu-Randall cloud fraction computed in progcld5 + cldcov(i,k) = clouds(i,k,1) + enddo + enddo + endif + if( .not. clduni) then ! --- call progcld5 for interaction with the radiation with setting ! --- uni_cld=.true. to keep precomputed cloud diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index 423a50ff0..9a46ae3d9 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -288,24 +288,6 @@ kind = kind_phys intent = inout optional = F -[nc] - standard_name = cloud_droplet_number_concentration - long_name = cloud droplet number concentration - units = kg-1 - dimensions = (horizontal_dimension,vertical_dimension) - type = real - kind = kind_phys - intent = inout - optional = T -[ni] - standard_name = ice_number_concentration - long_name = ice number concentration - units = kg-1 - dimensions = (horizontal_dimension,vertical_dimension) - type = real - kind = kind_phys - intent = inout - optional = F [nwfa] standard_name = water_friendly_aerosol_number_concentration long_name = number concentration of water-friendly aerosols From 4261b1554689bd5faad1370ef3f2ebf670dfb916 Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Fri, 14 Feb 2020 19:31:17 +0000 Subject: [PATCH 04/12] QC, Qi and NWFA are not needed in the parameters list as they come into this subroutine as the qgrs entries. Results before/after this change are identical. --- physics/GFS_rrtmg_pre.F90 | 75 ++++++++++++++++---------------------- physics/GFS_rrtmg_pre.meta | 27 -------------- 2 files changed, 31 insertions(+), 71 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 351862cf5..b5055757c 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -20,7 +20,7 @@ end subroutine GFS_rrtmg_pre_init ! in the CCPP version - they are defined in the interstitial_create routine subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input Tbd, Cldprop, Coupling, & - Radtend, qc, qi, nwfa, & ! input/output + Radtend, & ! input/output imfdeepcnv, imfdeepcnv_gf, & f_ice, f_rain, f_rimef, flgmin, cwm, & ! F-A mp scheme only lm, im, lmk, lmp, & ! input @@ -87,11 +87,6 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input type(GFS_cldprop_type), intent(in) :: Cldprop type(GFS_coupling_type), intent(in) :: Coupling - real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: qc - real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: qi - real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP), intent(in) :: nwfa - - integer, intent(in) :: im, lm, lmk, lmp integer, intent(in) :: imfdeepcnv, imfdeepcnv_gf integer, intent(out) :: kd, kt, kb @@ -154,7 +149,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input integer, intent(out) :: errflg ! Local variables - integer :: me, nfxr, ntrac, ntcw, ntiw, ncld, ntrw, ntsw, ntgl, ncndl, ntlnc, ntinc + integer :: me, nfxr, ntrac, ntcw, ntiw, ncld, ntrw, ntsw, ntgl, ncndl, ntlnc, ntinc, ntwa integer :: i, j, k, k1, k2, lsk, lv, n, itop, ibtc, LP1, lla, llb, lya, lyb @@ -170,7 +165,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input ! for Thompson MP real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP) :: & re_cloud, re_ice, re_snow, qv_mp, qc_mp, & - qi_mp, qs_mp, nc_mp, ni_mp + qi_mp, qs_mp, nc_mp, ni_mp, nwfa real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP+1) :: tem2db ! real(kind=kind_phys), dimension(size(Grid%xlon,1),Model%levr+LTP+1) :: hz @@ -205,6 +200,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input ntrw = Model%ntrw ntsw = Model%ntsw ntgl = Model%ntgl + ntwa = Model%ntwa ncndl = min(Model%ncnd,4) LP1 = LM + 1 ! num of in/out levels @@ -297,15 +293,6 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input tracer1(:,k1,j) = max(0.0, Statein%qgrs(:,k2,j)) enddo enddo - if ((Model%do_mynnedmf.or. (imfdeepcnv == imfdeepcnv_gf)) .and. Model%kdt > 1) then - ! for MYNN PBL and GF convective include subgrid clouds into tracer1 - do k = 1, LM - k1 = k + kd - k2 = k + lsk - tracer1(:,k1,ntcw) = max(0.0, qc(:,k2)) - tracer1(:,k1,ntiw) = max(0.0, qi(:,k2)) - enddo - endif ! if (ivflip == 0) then ! input data from toa to sfc do i = 1, IM @@ -595,6 +582,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input qs_mp (i,k) = tracer1(i,k,ntsw)/(1.-qvs) nc_mp (i,k) = tracer1(i,k,ntlnc)/(1.-qvs) ni_mp (i,k) = tracer1(i,k,ntinc)/(1.-qvs) + nwfa (i,k) = tracer1(i,k,ntwa) endif enddo enddo @@ -731,7 +719,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input endif end do end do - ! Call Thompson's subroutine to compoute effective radii + ! Call Thompson's subroutine to compute effective radii do i = 1, im call calc_effectRad (tlyr(i,:), plyr(i,:), qv_mp(i,:), qc_mp(i,:), & nc_mp(i,:), qi_mp(i,:), ni_mp(i,:), qs_mp(i,:), & @@ -747,32 +735,31 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input re_snow(i,k) = MAX(9.99, MIN(re_snow(i,k)*1.e6, 999.)) end do end do - if(1==2) then - write(0,'(a,3e16.7)') " before progclduni: re_cloud min/mean/max =", & - minval(re_cloud), & - sum(re_cloud)/real(size(re_cloud)), & - maxval(re_cloud) - write(0,'(a,3e16.7)') " before progclduni: re_ice min/mean/max =", & - minval(re_ice), & - sum(re_ice)/real(size(re_ice)), & - maxval(re_ice) - write(0,'(a,3e16.7)') " before progclduni: clouds3 min/mean/max =", & - minval(clouds3), & - sum(clouds3)/real(size(clouds3)), & - maxval(clouds3) - write(0,'(a,3e16.7)') " before progclduni: clouds5 min/mean/max =", & - minval(clouds5), & - sum(clouds5)/real(size(clouds5)), & - maxval(clouds5) - write(0,'(a,3e16.7)') " before progcld5: phy_f3d cl min/mean/max =", & - minval(Tbd%phy_f3d(:,:,Model%nleffr)), & - sum(Tbd%phy_f3d(:,:,Model%nleffr))/real(size(Tbd%phy_f3d(:,:,Model%nleffr))), & - maxval(Tbd%phy_f3d(:,:,Model%nleffr)) - write(0,'(a,3e16.7)')" before progcld5: phy_f3d ice min/mean/max =", & - minval(Tbd%phy_f3d(:,:,Model%nieffr)), & - sum(Tbd%phy_f3d(:,:,Model%nieffr))/real(size(Tbd%phy_f3d(:,:,Model%nieffr))), & - maxval(Tbd%phy_f3d(:,:,Model%nieffr)) - endif + + !write(0,'(a,3e16.7)') " before progclduni: re_cloud min/mean/max =", & + ! minval(re_cloud), & + ! sum(re_cloud)/real(size(re_cloud)), & + ! maxval(re_cloud) + !write(0,'(a,3e16.7)') " before progclduni: re_ice min/mean/max =", & + ! minval(re_ice), & + ! sum(re_ice)/real(size(re_ice)), & + ! maxval(re_ice) + !write(0,'(a,3e16.7)') " before progclduni: clouds3 min/mean/max =", & + ! minval(clouds3), & + ! sum(clouds3)/real(size(clouds3)), & + ! maxval(clouds3) + !write(0,'(a,3e16.7)') " before progclduni: clouds5 min/mean/max =", & + ! minval(clouds5), & + ! sum(clouds5)/real(size(clouds5)), & + ! maxval(clouds5) + !write(0,'(a,3e16.7)') " before progcld5: phy_f3d cl min/mean/max =", & + ! minval(Tbd%phy_f3d(:,:,Model%nleffr)), & + ! sum(Tbd%phy_f3d(:,:,Model%nleffr))/real(size(Tbd%phy_f3d(:,:,Model%nleffr))), & + ! maxval(Tbd%phy_f3d(:,:,Model%nleffr)) + !write(0,'(a,3e16.7)')" before progcld5: phy_f3d ice min/mean/max =", & + ! minval(Tbd%phy_f3d(:,:,Model%nieffr)), & + ! sum(Tbd%phy_f3d(:,:,Model%nieffr))/real(size(Tbd%phy_f3d(:,:,Model%nieffr))), & + ! maxval(Tbd%phy_f3d(:,:,Model%nieffr)) do k=1,lm k1 = k + kd diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index 9a46ae3d9..901015f04 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -270,33 +270,6 @@ kind = kind_phys intent = out optional = F -[qc] - standard_name = cloud_condensed_water_mixing_ratio - long_name = moist (dry+vapor, no condensates) mixing ratio of cloud water (condensate) - units = kg kg-1 - dimensions = (horizontal_dimension,vertical_dimension) - type = real - kind = kind_phys - intent = inout - optional = F -[qi] - standard_name = ice_water_mixing_ratio - long_name = moist (dry+vapor, no condensates) mixing ratio of ice water - units = kg kg-1 - dimensions = (horizontal_dimension,vertical_dimension) - type = real - kind = kind_phys - intent = inout - optional = F -[nwfa] - standard_name = water_friendly_aerosol_number_concentration - long_name = number concentration of water-friendly aerosols - units = kg-1 - dimensions = (horizontal_dimension,vertical_dimension) - type = real - kind = kind_phys - intent = inout - optional = T [imfdeepcnv] standard_name = flag_for_mass_flux_deep_convection_scheme long_name = flag for mass-flux deep convection scheme From 3a852e8a3cd016571a5b08ddffda28585b2347f9 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 21 Feb 2020 08:12:13 -0700 Subject: [PATCH 05/12] physics/GFS_PBL_generic.F90: add missing tracers to vertical diffusion array for Thompson MP --- physics/GFS_PBL_generic.F90 | 40 ++++++++++++++++++++----------------- 1 file changed, 22 insertions(+), 18 deletions(-) diff --git a/physics/GFS_PBL_generic.F90 b/physics/GFS_PBL_generic.F90 index 7e28d2cec..e157013ec 100644 --- a/physics/GFS_PBL_generic.F90 +++ b/physics/GFS_PBL_generic.F90 @@ -150,12 +150,13 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, vdftra(i,k,3) = qgrs(i,k,ntiw) vdftra(i,k,4) = qgrs(i,k,ntrw) vdftra(i,k,5) = qgrs(i,k,ntsw) - vdftra(i,k,6) = qgrs(i,k,ntlnc) - vdftra(i,k,7) = qgrs(i,k,ntinc) - vdftra(i,k,8) = qgrs(i,k,ntrnc) - vdftra(i,k,9) = qgrs(i,k,ntoz) - vdftra(i,k,10) = qgrs(i,k,ntwa) - vdftra(i,k,11) = qgrs(i,k,ntia) + vdftra(i,k,6) = qgrs(i,k,ntgl) + vdftra(i,k,7) = qgrs(i,k,ntlnc) + vdftra(i,k,8) = qgrs(i,k,ntinc) + vdftra(i,k,9) = qgrs(i,k,ntrnc) + vdftra(i,k,10) = qgrs(i,k,ntoz) + vdftra(i,k,11) = qgrs(i,k,ntwa) + vdftra(i,k,12) = qgrs(i,k,ntia) enddo enddo else @@ -166,9 +167,10 @@ subroutine GFS_PBL_generic_pre_run (im, levs, nvdiff, ntrac, vdftra(i,k,3) = qgrs(i,k,ntiw) vdftra(i,k,4) = qgrs(i,k,ntrw) vdftra(i,k,5) = qgrs(i,k,ntsw) - vdftra(i,k,6) = qgrs(i,k,ntinc) - vdftra(i,k,7) = qgrs(i,k,ntrnc) - vdftra(i,k,8) = qgrs(i,k,ntoz) + vdftra(i,k,6) = qgrs(i,k,ntgl) + vdftra(i,k,7) = qgrs(i,k,ntinc) + vdftra(i,k,8) = qgrs(i,k,ntrnc) + vdftra(i,k,9) = qgrs(i,k,ntoz) enddo enddo endif @@ -406,12 +408,13 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, dqdt(i,k,ntiw) = dvdftra(i,k,3) dqdt(i,k,ntrw) = dvdftra(i,k,4) dqdt(i,k,ntsw) = dvdftra(i,k,5) - dqdt(i,k,ntlnc) = dvdftra(i,k,6) - dqdt(i,k,ntinc) = dvdftra(i,k,7) - dqdt(i,k,ntrnc) = dvdftra(i,k,8) - dqdt(i,k,ntoz) = dvdftra(i,k,9) - dqdt(i,k,ntwa) = dvdftra(i,k,10) - dqdt(i,k,ntia) = dvdftra(i,k,11) + dqdt(i,k,ntgl) = dvdftra(i,k,6) + dqdt(i,k,ntlnc) = dvdftra(i,k,7) + dqdt(i,k,ntinc) = dvdftra(i,k,8) + dqdt(i,k,ntrnc) = dvdftra(i,k,9) + dqdt(i,k,ntoz) = dvdftra(i,k,10) + dqdt(i,k,ntwa) = dvdftra(i,k,11) + dqdt(i,k,ntia) = dvdftra(i,k,12) enddo enddo else @@ -422,9 +425,10 @@ subroutine GFS_PBL_generic_post_run (im, levs, nvdiff, ntrac, dqdt(i,k,ntiw) = dvdftra(i,k,3) dqdt(i,k,ntrw) = dvdftra(i,k,4) dqdt(i,k,ntsw) = dvdftra(i,k,5) - dqdt(i,k,ntinc) = dvdftra(i,k,6) - dqdt(i,k,ntrnc) = dvdftra(i,k,7) - dqdt(i,k,ntoz) = dvdftra(i,k,8) + dqdt(i,k,ntgl) = dvdftra(i,k,6) + dqdt(i,k,ntinc) = dvdftra(i,k,7) + dqdt(i,k,ntrnc) = dvdftra(i,k,8) + dqdt(i,k,ntoz) = dvdftra(i,k,9) enddo enddo endif From 2edeeadb4416937dfb18eb32f3bf327449d47de2 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Thu, 12 Mar 2020 15:14:23 -0600 Subject: [PATCH 06/12] Bugfixes: uninitialized data before entering effective radii calculation; array qci_conv may not be allocated, thus use assumed-size declaration --- physics/GFS_rrtmg_pre.F90 | 4 ++++ physics/module_SGSCloud_RadPre.F90 | 4 ++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 170cb707a..d123c9e4b 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -729,6 +729,10 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input end do ! Call Thompson's subroutine to compute effective radii do i=1,im + ! Initialize to default in units m as in module_mp_thompson.F90 + re_cloud(i,:) = 2.49E-6 + re_ice(i,:) = 4.99E-6 + re_snow(i,:) = 9.99E-6 call calc_effectRad (tlyr(i,:), plyr(i,:), qv_mp(i,:), qc_mp(i,:), & nc_mp(i,:), qi_mp(i,:), ni_mp(i,:), qs_mp(i,:), & re_cloud(i,:), re_ice(i,:), re_snow(i,:), 1, lm ) diff --git a/physics/module_SGSCloud_RadPre.F90 b/physics/module_SGSCloud_RadPre.F90 index 15ac383f5..544fe1004 100644 --- a/physics/module_SGSCloud_RadPre.F90 +++ b/physics/module_SGSCloud_RadPre.F90 @@ -61,13 +61,13 @@ subroutine sgscloud_radpre_run( & logical, intent(in) :: flag_init, flag_restart, do_mynnedmf real(kind=kind_phys), dimension(im,levs), intent(inout) :: qc, qi real(kind=kind_phys), dimension(im,levs), intent(inout) :: qr, qs - real(kind=kind_phys), dimension(im,levs), intent(inout) :: qci_conv + ! qci_conv only allocated if GF is used + real(kind=kind_phys), dimension(:,:), intent(inout) :: qci_conv real(kind=kind_phys), dimension(im,levs), intent(in) :: T3D,delp real(kind=kind_phys), dimension(im,levs), intent(inout) :: & & clouds1,clouds2,clouds3,clouds4,clouds5 real(kind=kind_phys), dimension(im,levs), intent(inout) :: qc_save, qi_save real(kind=kind_phys), dimension(im,levs), intent(in) :: qc_bl, cldfra_bl - ! DH* TODO add intent() information for delp,clouds1,clouds2,clouds3,clouds4,clouds5 real(kind=kind_phys), dimension(im), intent(in) :: slmsk, xlat, de_lgth real(kind=kind_phys), dimension(im,nlay), intent(in) :: plyr, dz real(kind=kind_phys), dimension(im,5), intent(inout) :: cldsa From 5960e5e8a1e15ebb7040656ebe1ea2c62253ae87 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 17 Mar 2020 14:30:47 -0600 Subject: [PATCH 07/12] Cosmetic changes to physics/GFS_debug.F90 --- physics/GFS_debug.F90 | 2 ++ 1 file changed, 2 insertions(+) diff --git a/physics/GFS_debug.F90 b/physics/GFS_debug.F90 index 486ee604e..3bb50d9ef 100644 --- a/physics/GFS_debug.F90 +++ b/physics/GFS_debug.F90 @@ -310,6 +310,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank,omprank, blkno, 'Diag%tdomzr ', Diag%tdomzr) call print_var(mpirank,omprank, blkno, 'Diag%tdomip ', Diag%tdomip) call print_var(mpirank,omprank, blkno, 'Diag%tdoms ', Diag%tdoms) + ! CCPP/RUC only if (Model%lsm == Model%lsm_ruc) then call print_var(mpirank,omprank, blkno, 'Diag%wet1 ', Sfcprop%wetness) else @@ -345,6 +346,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, if(Model%lradar) then call print_var(mpirank,omprank, blkno, 'Diag%refl_10cm ', Diag%refl_10cm) end if + ! CCPP/MYNNPBL only if (Model%do_mynnedmf) then call print_var(mpirank,omprank, blkno, 'Diag%edmf_a ', Diag%edmf_a) call print_var(mpirank,omprank, blkno, 'Diag%edmf_w ', Diag%edmf_w) From 10e357f8f03d03ec2a704529685e2b047cec9af0 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 20 Mar 2020 09:34:37 -0600 Subject: [PATCH 08/12] physics/module_MYNNSFC_wrapper.F90: add comment about CCPP being able to do automatic unit conversions --- physics/module_MYNNSFC_wrapper.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/physics/module_MYNNSFC_wrapper.F90 b/physics/module_MYNNSFC_wrapper.F90 index 951d7e7c8..9fd71c37d 100644 --- a/physics/module_MYNNSFC_wrapper.F90 +++ b/physics/module_MYNNSFC_wrapper.F90 @@ -186,9 +186,10 @@ SUBROUTINE mynnsfc_wrapper_run( & endif qgh(i)=0.0 !snowh(i)=snowd(i)*800. !mm -> m + ! DH* note - this could be automated (CCPP knows how to convert cm to m) znt_lnd(i)=znt_lnd(i)*0.01 !cm -> m znt_ocn(i)=znt_ocn(i)*0.01 !cm -> m - znt_ice(i)=znt_ice(i)*0.01 !cm -> m + znt_ice(i)=znt_ice(i)*0.01 !cm -> m ts(i)=tskin_ocn(i)/exner(i,1) !theta mavail(i)=1.0 !???? cpm(i)=cp @@ -272,6 +273,7 @@ SUBROUTINE mynnsfc_wrapper_run( & !NOTE: evap & qflx will be solved for later !qflx(i)=QFX(i)/ !evap(i)=QFX(i) !or /rho ?? + ! DH* note - this could be automated (CCPP knows how to convert m to cm) znt_lnd(i)=znt_lnd(i)*100. !m -> cm znt_ocn(i)=znt_ocn(i)*100. znt_ice(i)=znt_ice(i)*100. From 92b6ee80c7ee93b91b8072c1f6b4b7a619d4af44 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 23 Mar 2020 10:32:06 -0600 Subject: [PATCH 09/12] physics/module_MYNNSFC_wrapper.F90: perform unit conversion m <-> cm only for valid data --- physics/module_MYNNSFC_wrapper.F90 | 98 ++++++++++++++++-------------- 1 file changed, 54 insertions(+), 44 deletions(-) diff --git a/physics/module_MYNNSFC_wrapper.F90 b/physics/module_MYNNSFC_wrapper.F90 index 9fd71c37d..42d0108a1 100644 --- a/physics/module_MYNNSFC_wrapper.F90 +++ b/physics/module_MYNNSFC_wrapper.F90 @@ -168,32 +168,38 @@ SUBROUTINE mynnsfc_wrapper_run( & ! write(0,*)"iter=",iter ! endif - !prep MYNN-only variables - do k=1,2 !levs - do i=1,im - dz(i,k)=(phii(i,k+1) - phii(i,k))*g_inv - th(i,k)=t3d(i,k)/exner(i,k) - !qc(i,k)=MAX(qgrs(i,k,ntcw),0.0) - qv(i,k)=qvsh(i,k)/(1.0 - qvsh(i,k)) - pattern_spp_pbl(i,k)=0.0 - enddo - enddo - 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 SFCLAY_mynn - else - xland(i)=2.0 - endif - qgh(i)=0.0 - !snowh(i)=snowd(i)*800. !mm -> m - ! DH* note - this could be automated (CCPP knows how to convert cm to m) - znt_lnd(i)=znt_lnd(i)*0.01 !cm -> m - znt_ocn(i)=znt_ocn(i)*0.01 !cm -> m - znt_ice(i)=znt_ice(i)*0.01 !cm -> m - ts(i)=tskin_ocn(i)/exner(i,1) !theta - mavail(i)=1.0 !???? - cpm(i)=cp - enddo + ! prep MYNN-only variables + do k=1,2 !levs + do i=1,im + dz(i,k)=(phii(i,k+1) - phii(i,k))*g_inv + th(i,k)=t3d(i,k)/exner(i,k) + !qc(i,k)=MAX(qgrs(i,k,ntcw),0.0) + qv(i,k)=qvsh(i,k)/(1.0 - qvsh(i,k)) + pattern_spp_pbl(i,k)=0.0 + enddo + enddo + 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 SFCLAY_mynn + else + xland(i)=2.0 + endif + qgh(i)=0.0 + !snowh(i)=snowd(i)*800. !mm -> m + !znt_lnd(i)=znt_lnd(i)*0.01 !cm -> m + !znt_ocn(i)=znt_ocn(i)*0.01 !cm -> m + !znt_ice(i)=znt_ice(i)*0.01 !cm -> m + ! DH* do the following line only if wet(i)? + ts(i)=tskin_ocn(i)/exner(i,1) !theta + ! *DH + mavail(i)=1.0 !???? + cpm(i)=cp + enddo + + ! cm -> m + where (dry) znt_lnd=znt_lnd*0.01 + where (wet) znt_ocn=znt_ocn*0.01 + where (icy) znt_ice=znt_ice*0.01 ! if (lprnt) then ! write(0,*)"CALLING SFCLAY_mynn; input:" @@ -261,24 +267,28 @@ SUBROUTINE mynnsfc_wrapper_run( & its=1,ite=im, jts=1,jte=1, kts=1,kte=levs ) - ! POST MYNN SURFACE LAYER (INTERSTITIAL) WORK: - do i = 1, im - !* Taken from sfc_nst.f - !* ch = surface exchange coeff heat & moisture(m/s) im - !* rch(i) = rho_a(i) * cp * ch(i) * wind(i) - !* hflx(i) = rch(i) * (tsurf(i) - theta1(i)) !K m s-1 - !* hflx(i)=hfx(i)/(rho(i,1)*cp) - now calculated inside module_sf_mynn.F90 - !* Taken from sfc_nst.f - !* evap(i) = elocp * rch(i) * (qss(i) - q0(i)) !kg kg-1 m s-1 - !NOTE: evap & qflx will be solved for later - !qflx(i)=QFX(i)/ - !evap(i)=QFX(i) !or /rho ?? - ! DH* note - this could be automated (CCPP knows how to convert m to cm) - znt_lnd(i)=znt_lnd(i)*100. !m -> cm - znt_ocn(i)=znt_ocn(i)*100. - znt_ice(i)=znt_ice(i)*100. - enddo - + !! POST MYNN SURFACE LAYER (INTERSTITIAL) WORK: + !do i = 1, im + ! !* Taken from sfc_nst.f + ! !* ch = surface exchange coeff heat & moisture(m/s) im + ! !* rch(i) = rho_a(i) * cp * ch(i) * wind(i) + ! !* hflx(i) = rch(i) * (tsurf(i) - theta1(i)) !K m s-1 + ! !* hflx(i)=hfx(i)/(rho(i,1)*cp) - now calculated inside module_sf_mynn.F90 + ! !* Taken from sfc_nst.f + ! !* evap(i) = elocp * rch(i) * (qss(i) - q0(i)) !kg kg-1 m s-1 + ! !NOTE: evap & qflx will be solved for later + ! !qflx(i)=QFX(i)/ + ! !evap(i)=QFX(i) !or /rho ?? + ! ! DH* note - this could be automated (CCPP knows how to convert m to cm) + ! znt_lnd(i)=znt_lnd(i)*100. !m -> cm + ! znt_ocn(i)=znt_ocn(i)*100. + ! znt_ice(i)=znt_ice(i)*100. + !enddo + + ! m -> cm + where (dry) znt_lnd=znt_lnd*100. + where (wet) znt_ocn=znt_ocn*100. + where (icy) znt_ice=znt_ice*100. ! if (lprnt) then ! write(0,*) From 65da24ee5def2e86a19614bf2536af76b9074d99 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 23 Mar 2020 16:52:00 -0600 Subject: [PATCH 10/12] physics/GFS_surface_composites.*: initialize composites uustar_*, qss_*, hflx_* --- physics/GFS_surface_composites.F90 | 21 +++++--- physics/GFS_surface_composites.meta | 81 +++++++++++++++++++++++++++++ 2 files changed, 96 insertions(+), 6 deletions(-) diff --git a/physics/GFS_surface_composites.F90 b/physics/GFS_surface_composites.F90 index 9636eb384..0060e1a7b 100644 --- a/physics/GFS_surface_composites.F90 +++ b/physics/GFS_surface_composites.F90 @@ -28,10 +28,11 @@ end subroutine GFS_surface_composites_pre_finalize subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, landfrac, lakefrac, oceanfrac, & frland, dry, icy, lake, ocean, wet, cice, cimin, zorl, zorlo, zorll, zorl_ocn, & zorl_lnd, zorl_ice, snowd, snowd_ocn, snowd_lnd, snowd_ice, tprcp, tprcp_ocn, & - tprcp_lnd, tprcp_ice, uustar, uustar_lnd, uustar_ice, weasd, weasd_ocn, & - weasd_lnd, weasd_ice, ep1d_ice, tsfc, tsfco, tsfcl, tsfc_ocn, tsfc_lnd, & - tsfc_ice, tisfc, tice, tsurf, tsurf_ocn, tsurf_lnd, tsurf_ice, gflx_ice, & - tgice, islmsk, semis_rad, semis_ocn, semis_lnd, semis_ice, & + tprcp_lnd, tprcp_ice, uustar, uustar_ocn, uustar_lnd, uustar_ice, & + weasd, weasd_ocn, weasd_lnd, weasd_ice, ep1d_ice, tsfc, tsfco, tsfcl, tsfc_ocn,& + tsfc_lnd, tsfc_ice, tisfc, tice, tsurf, tsurf_ocn, tsurf_lnd, tsurf_ice, & + gflx_ice, tgice, islmsk, semis_rad, semis_ocn, semis_lnd, semis_ice, & + qss, qss_ocn, qss_lnd, qss_ice, hflx, hflx_ocn, hflx_lnd, hflx_ice, & min_lakeice, min_seaice, errmsg, errflg) implicit none @@ -45,12 +46,13 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, lan real(kind=kind_phys), dimension(im), intent(in ) :: landfrac, lakefrac, oceanfrac real(kind=kind_phys), dimension(im), intent(inout) :: cice real(kind=kind_phys), dimension(im), intent( out) :: frland - real(kind=kind_phys), dimension(im), intent(in ) :: zorl, snowd, tprcp, uustar, weasd + real(kind=kind_phys), dimension(im), intent(in ) :: zorl, snowd, tprcp, uustar, weasd, qss, hflx real(kind=kind_phys), dimension(im), intent(inout) :: zorlo, zorll, tsfc, tsfco, tsfcl, tisfc, tsurf real(kind=kind_phys), dimension(im), intent(inout) :: snowd_ocn, snowd_lnd, snowd_ice, tprcp_ocn, & tprcp_lnd, tprcp_ice, zorl_ocn, zorl_lnd, zorl_ice, tsfc_ocn, tsfc_lnd, tsfc_ice, tsurf_ocn, & - tsurf_lnd, tsurf_ice, uustar_lnd, uustar_ice, weasd_ocn, weasd_lnd, weasd_ice, ep1d_ice, gflx_ice + tsurf_lnd, tsurf_ice, uustar_ocn, uustar_lnd, uustar_ice, weasd_ocn, weasd_lnd, weasd_ice, & + qss_ocn, qss_lnd, qss_ice, hflx_ocn, hflx_lnd, hflx_ice, ep1d_ice, gflx_ice real(kind=kind_phys), dimension(im), intent( out) :: tice real(kind=kind_phys), intent(in ) :: tgice integer, dimension(im), intent(in ) :: islmsk @@ -145,6 +147,7 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, lan tprcp_lnd(i) = tprcp(i) tprcp_ice(i) = tprcp(i) if (wet(i)) then ! Water + uustar_ocn(i) = uustar(i) zorl_ocn(i) = zorlo(i) tsfc_ocn(i) = tsfco(i) tsurf_ocn(i) = tsfco(i) @@ -153,6 +156,8 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, lan weasd_ocn(i) = zero snowd_ocn(i) = zero semis_ocn(i) = 0.984d0 + qss_ocn(i) = qss(i) + hflx_ocn(i) = hflx(i) endif if (dry(i)) then ! Land uustar_lnd(i) = uustar(i) @@ -162,6 +167,8 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, lan tsurf_lnd(i) = tsfcl(i) snowd_lnd(i) = snowd(i) semis_lnd(i) = semis_rad(i) + qss_lnd(i) = qss(i) + hflx_lnd(i) = hflx(i) end if if (icy(i)) then ! Ice uustar_ice(i) = uustar(i) @@ -173,6 +180,8 @@ subroutine GFS_surface_composites_pre_run (im, frac_grid, flag_cice, cplflx, lan ep1d_ice(i) = zero gflx_ice(i) = zero semis_ice(i) = 0.95d0 + qss_ice(i) = qss(i) + hflx_ice(i) = hflx(i) end if enddo diff --git a/physics/GFS_surface_composites.meta b/physics/GFS_surface_composites.meta index 74c6b9575..bf613e160 100644 --- a/physics/GFS_surface_composites.meta +++ b/physics/GFS_surface_composites.meta @@ -262,6 +262,15 @@ kind = kind_phys intent = in optional = F +[uustar_ocn] + standard_name = surface_friction_velocity_over_ocean + long_name = surface friction velocity over ocean + units = m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F [uustar_lnd] standard_name = surface_friction_velocity_over_land long_name = surface friction velocity over land @@ -495,6 +504,78 @@ kind = kind_phys intent = inout optional = F +[qss] + standard_name = surface_specific_humidity + long_name = surface air saturation specific humidity + units = kg kg-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[qss_ocn] + standard_name = surface_specific_humidity_over_ocean + long_name = surface air saturation specific humidity over ocean + units = kg kg-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qss_lnd] + standard_name = surface_specific_humidity_over_land + long_name = surface air saturation specific humidity over land + units = kg kg-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qss_ice] + standard_name = surface_specific_humidity_over_ice + long_name = surface air saturation specific humidity over ice + units = kg kg-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[hflx] + standard_name = kinematic_surface_upward_sensible_heat_flux + long_name = kinematic surface upward sensible heat flux + units = K m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[hflx_ocn] + standard_name = kinematic_surface_upward_sensible_heat_flux_over_ocean + long_name = kinematic surface upward sensible heat flux over ocean + units = K m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[hflx_lnd] + standard_name = kinematic_surface_upward_sensible_heat_flux_over_land + long_name = kinematic surface upward sensible heat flux over land + units = K m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[hflx_ice] + standard_name = kinematic_surface_upward_sensible_heat_flux_over_ice + long_name = kinematic surface upward sensible heat flux over ice + units = K m s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F [min_lakeice] standard_name = lake_ice_minimum long_name = minimum lake ice value From a6f3dedf3311863e4a3fd66086d7c53472b11fdc Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Thu, 26 Mar 2020 13:06:25 -0600 Subject: [PATCH 11/12] Updates from @joeolson42 for physics/module_MYNNSFC_wrapper.F90, physics/module_MYNNSFC_wrapper.meta, physics/module_sf_mynn.F90 --- physics/module_MYNNSFC_wrapper.F90 | 14 +- physics/module_MYNNSFC_wrapper.meta | 17 ++ physics/module_sf_mynn.F90 | 256 ++++++++++++++++++---------- 3 files changed, 189 insertions(+), 98 deletions(-) diff --git a/physics/module_MYNNSFC_wrapper.F90 b/physics/module_MYNNSFC_wrapper.F90 index 42d0108a1..82cdbca76 100644 --- a/physics/module_MYNNSFC_wrapper.F90 +++ b/physics/module_MYNNSFC_wrapper.F90 @@ -28,7 +28,7 @@ end subroutine mynnsfc_wrapper_finalize SUBROUTINE mynnsfc_wrapper_run( & & ix,im,levs, & & itimestep,iter, & - & flag_init,flag_restart, & + & flag_init,flag_restart,lsm, & & delt,dx, & & u, v, t3d, qvsh, qc, prsl, phii, & & exner, ps, PBLH, slmsk, & @@ -47,8 +47,8 @@ SUBROUTINE mynnsfc_wrapper_run( & & fh_ocn, fh_lnd, fh_ice, & !intent(inout) & fm10_ocn, fm10_lnd, fm10_ice, & !intent(inout) & fh2_ocn, fh2_lnd, fh2_ice, & !intent(inout) - & QSFC, USTM, ZOL, MOL, RMOL, & - & WSPD, ch, HFLX, QFLX, LH, & + & QSFC, qsfc_ruc, USTM, ZOL, MOL, & + & RMOL, WSPD, ch, HFLX, QFLX, LH, & & FLHC, FLQC, & & U10, V10, TH2, T2, Q2, & & wstar, CHS2, CQS2, & @@ -106,7 +106,7 @@ SUBROUTINE mynnsfc_wrapper_run( & !MYNN-1D REAL :: delt INTEGER :: im, ix, levs - INTEGER :: iter, k, i, itimestep + INTEGER :: iter, k, i, itimestep, lsm LOGICAL :: flag_init,flag_restart,lprnt INTEGER :: IDS,IDE,JDS,JDE,KDS,KDE, & & IMS,IME,JMS,JME,KMS,KME, & @@ -146,7 +146,7 @@ SUBROUTINE mynnsfc_wrapper_run( & & dx, pblh, slmsk, ps real(kind=kind_phys), dimension(im), intent(inout) :: & - & ustm, hflx, qflx, wspd, qsfc, & + & ustm, hflx, qflx, wspd, qsfc, qsfc_ruc, & & FLHC, FLQC, U10, V10, TH2, T2, Q2, & & CHS2, CQS2, rmol, zol, mol, ch, & & lh, wstar @@ -237,7 +237,7 @@ SUBROUTINE mynnsfc_wrapper_run( & CP=cp,G=g,ROVCP=rcp,R=r_d,XLV=xlv, & SVP1=svp1,SVP2=svp2,SVP3=svp3,SVPT0=svpt0, & EP1=ep_1,EP2=ep_2,KARMAN=karman, & - ISFFLX=isfflx,isftcflx=isftcflx, & + ISFFLX=isfflx,isftcflx=isftcflx,LSM=lsm, & iz0tlnd=iz0tlnd,itimestep=itimestep,iter=iter, & wet=wet, dry=dry, icy=icy, & !intent(in) tskin_ocn=tskin_ocn, tskin_lnd=tskin_lnd, tskin_ice=tskin_ice, & !intent(in) @@ -258,7 +258,7 @@ SUBROUTINE mynnsfc_wrapper_run( & ZNT=znt,USTM=ustm,ZOL=zol,MOL=mol,RMOL=rmol, & psim=psim,psih=psih, & HFLX=hflx,HFX=hfx,QFLX=qflx,QFX=qfx,LH=lh,FLHC=flhc,FLQC=flqc, & - QGH=qgh,QSFC=qsfc, & + QGH=qgh,QSFC=qsfc,QSFC_RUC=qsfc_ruc, & U10=u10,V10=v10,TH2=th2,T2=t2,Q2=q2, & GZ1OZ0=GZ1OZ0,WSPD=wspd,wstar=wstar, & spp_pbl=spp_pbl,pattern_spp_pbl=pattern_spp_pbl, & diff --git a/physics/module_MYNNSFC_wrapper.meta b/physics/module_MYNNSFC_wrapper.meta index 0a988f575..b12837233 100644 --- a/physics/module_MYNNSFC_wrapper.meta +++ b/physics/module_MYNNSFC_wrapper.meta @@ -57,6 +57,14 @@ type = logical intent = in optional = F +[lsm] + standard_name = flag_for_land_surface_scheme + long_name = flag for land surface model + units = flag + dimensions = () + type = integer + intent = in + optional = F [delt] standard_name = time_step_for_physics long_name = time step for physics @@ -585,6 +593,15 @@ kind = kind_phys intent = inout optional = F +[qsfc_ruc] + standard_name = water_vapor_mixing_ratio_at_surface + long_name = water vapor mixing ratio at surface + units = kg kg-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F [ustm] standard_name = surface_friction_velocity_drag long_name = friction velocity isolated for momentum only diff --git a/physics/module_sf_mynn.F90 b/physics/module_sf_mynn.F90 index 788ff0ace..2ac9d832c 100644 --- a/physics/module_sf_mynn.F90 +++ b/physics/module_sf_mynn.F90 @@ -106,6 +106,9 @@ MODULE module_sf_mynn !1: some step-by-step output !2: everything - heavy I/O LOGICAL, PARAMETER :: compute_diag = .false. + LOGICAL, PARAMETER :: compute_flux = .false. !shouldn't need compute + ! these in FV3. They will be written over anyway. + ! Computing the fluxes here is leftover from the WRF world. REAL, DIMENSION(0:1000 ),SAVE :: psim_stab,psim_unstab, & psih_stab,psih_unstab @@ -137,7 +140,8 @@ SUBROUTINE SFCLAY_mynn( & PSFCPA,PBLH,MAVAIL,XLAND,DX, & !in CP,G,ROVCP,R,XLV, & !in SVP1,SVP2,SVP3,SVPT0,EP1,EP2,KARMAN, & !in - ISFFLX,isftcflx,iz0tlnd,itimestep,iter,& !in + ISFFLX,isftcflx,lsm,iz0tlnd, & !in + itimestep,iter, & !in wet, dry, icy, & !intent(in) tskin_ocn, tskin_lnd, tskin_ice, & !intent(in) tsurf_ocn, tsurf_lnd, tsurf_ice, & !intent(in) @@ -157,7 +161,7 @@ SUBROUTINE SFCLAY_mynn( & ZNT,USTM,ZOL,MOL,RMOL, & PSIM,PSIH, & HFLX,HFX,QFLX,QFX,LH,FLHC,FLQC, & - QGH,QSFC, & + QGH,QSFC,QSFC_RUC, & U10,V10,TH2,T2,Q2, & GZ1OZ0,WSPD,WSTAR, & spp_pbl,pattern_spp_pbl, & @@ -268,8 +272,8 @@ SUBROUTINE SFCLAY_mynn( & REAL, INTENT(IN) :: EP1,EP2,KARMAN REAL, INTENT(IN) :: CP,G,ROVCP,R,XLV !,DX !NAMELIST OPTIONS: - INTEGER, INTENT(IN) :: ISFFLX - INTEGER, OPTIONAL, INTENT(IN) :: ISFTCFLX, IZ0TLND + INTEGER, INTENT(IN) :: ISFFLX, LSM + INTEGER, OPTIONAL, INTENT(IN) :: ISFTCFLX, IZ0TLND INTEGER, OPTIONAL, INTENT(IN) :: spp_pbl !=================================== @@ -306,7 +310,8 @@ SUBROUTINE SFCLAY_mynn( & QFLX,QFX, & LH, & MOL,RMOL, & - QSFC, QGH, & + QSFC, & + QGH, & ZNT, & ZOL, & USTM, & @@ -339,7 +344,8 @@ SUBROUTINE SFCLAY_mynn( & & fh_ocn, fh_lnd, fh_ice, & & fm10_ocn, fm10_lnd, fm10_ice, & & fh2_ocn, fh2_lnd, fh2_ice, & - & qsfc_ocn, qsfc_lnd, qsfc_ice + & qsfc_ocn, qsfc_lnd, qsfc_ice, & + & qsfc_ruc !ADDITIONAL OUTPUT !JOE-begin @@ -402,10 +408,21 @@ SUBROUTINE SFCLAY_mynn( & UST_ICE(i)=MAX(0.04*SQRT(U1D(i)*U1D(i) + V1D(i)*V1D(i)),0.001) MOL(i,j)=0. ! Tstar QSFC(i,j)=QV3D(i,kts,j)/(1.+QV3D(i,kts,j)) + QSFC_OCN(i)=QSFC(i,j) + QSFC_LND(i)=QSFC(i,j) + QSFC_ICE(i)=QSFC(i,j) qstar(i,j)=0.0 QFX(i,j)=0. HFX(i,j)=0. + QFLX(i,j)=0. + HFLX(i,j)=0. ENDDO + ELSE + IF (LSM == 3) THEN + DO i=its,ite + QSFC_LND(i)=QSFC_RUC(i) + ENDDO + ENDIF ENDIF CALL SFCLAY1D_mynn( & @@ -453,7 +470,10 @@ END SUBROUTINE SFCLAY_MYNN !------------------------------------------------------------------- !>\ingroup module_sf_mynn_mod -!! This subroutine calculates +!! This subroutine calculates u*, z/L, and the exchange coefficients +!! which are passed to subsequent scheme to calculate the fluxes. +!! This scheme has options to calculate the fluxes and near-surface +!! diagnostics, as was needed in WRF, but these are skipped for FV3. SUBROUTINE SFCLAY1D_mynn( & J,U1D,V1D,T1D,QV1D,P1D,dz8w1d,U1D2,V1D2,dz2w1d, & PSFCPA,PBLH,MAVAIL,XLAND,DX, & @@ -621,20 +641,27 @@ SUBROUTINE SFCLAY1D_mynn( & !------------------------------------------------------------------- IF (debug_code >= 1) THEN - write(*,*)"ITIMESTEP=",ITIMESTEP," iter=",iter + write(0,*)"ITIMESTEP=",ITIMESTEP," iter=",iter DO I=its,ite - write(*,*)"=== input to mynnsfclayer, i:", i - !write(*,*)" land, ice, water" - write(*,*)"dry=",dry(i)," icy=",icy(i)," wet=",wet(i) - write(*,*)"tsk=", tskin_lnd(i),tskin_ice(i),tskin_ocn(i) - write(*,*)"tsurf=", tsurf_lnd(i),tsurf_ice(i),tsurf_ocn(i) - write(*,*)"qsfc=", qsfc_lnd(i),qsfc_ice(i),qsfc_ocn(i) - write(*,*)"znt=", znt_lnd(i),znt_ice(i),znt_ocn(i) - write(*,*)"ust=", ust_lnd(i),ust_ice(i),ust_ocn(i) - write(*,*)"snowh=", snowh_lnd(i),snowh_ice(i),snowh_ocn(i) - write(*,*)"psfcpa=",PSFCPA(i)," dz=",dz8w1d(i) - write(*,'(A5,F0.8,A6,F0.6,A6,F5.0)') & - "qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + write(0,*)"=== imortant input to mynnsfclayer, i:", i + IF (dry(i)) THEN + write(0,*)"dry=",dry(i)," tsk=", tskin_lnd(i),& + " tsurf=", tsurf_lnd(i)," qsfc=", qsfc_lnd(i)," znt=", znt_lnd(i),& + " ust=", ust_lnd(i)," snowh=", snowh_lnd(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF + IF (icy(i)) THEN + write(0,*)"icy=",icy(i)," tsk=", tskin_ice(i),& + " tsurf=", tsurf_ice(i)," qsfc=", qsfc_ice(i)," znt=", znt_ice(i),& + " ust=", ust_ice(i)," snowh=", snowh_ice(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF + IF (wet(i)) THEN + write(0,*)"wet=",wet(i)," tsk=", tskin_ocn(i),& + " tsurf=", tsurf_ocn(i)," qsfc=", qsfc_ocn(i)," znt=", znt_ocn(i),& + " ust=", ust_ocn(i)," snowh=", snowh_ocn(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF ENDDO ENDIF @@ -1161,8 +1188,14 @@ SUBROUTINE SFCLAY1D_mynn( & ENDIF IF (debug_code >= 1) THEN - write(0,*)"===(wet) capture bad input in mynn sfc layer, i=:",i - write(0,*)"rb=", rb_ocn(I)," ZNT=", ZNTstoch_ocn(i)," ZT=",Zt_ocn(i) + IF (ZNTstoch_ocn(i) < 1E-8 .OR. Zt_ocn(i) < 1E-10) THEN + write(0,*)"===(wet) capture bad input in mynn sfc layer, i=:",i + write(0,*)"rb=", rb_ocn(I)," ZNT=", ZNTstoch_ocn(i)," ZT=",Zt_ocn(i) + write(0,*)" tsk=", tskin_ocn(i)," prev z/L=",ZOL(I),& + " tsurf=", tsurf_ocn(i)," qsfc=", qsfc_ocn(i)," znt=", znt_ocn(i),& + " ust=", ust_ocn(i)," snowh=", snowh_ocn(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF ENDIF !Use Pedros iterative function to find z/L zol(I)=zolri(rb_ocn(I),ZA(I),ZNTstoch_ocn(I),ZT_ocn(I),ZOL(I)) @@ -1219,8 +1252,14 @@ SUBROUTINE SFCLAY1D_mynn( & ENDIF IF (debug_code >= 1) THEN - write(0,*)"===(wet) capture bad input in mynn sfc layer, i=:",i - write(0,*)"rb=", rb_ocn(I)," ZNT=", ZNTstoch_ocn(i)," ZT=",Zt_ocn(i) + IF (ZNTstoch_ocn(i) < 1E-8 .OR. Zt_ocn(i) < 1E-10) THEN + write(0,*)"===(wet) capture bad input in mynn sfc layer, i=:",i + write(0,*)"rb=", rb_ocn(I)," ZNT=", ZNTstoch_ocn(i)," ZT=",Zt_ocn(i) + write(0,*)" tsk=", tskin_ocn(i)," wstar=",wstar(i)," prev z/L=",ZOL(I),& + " tsurf=", tsurf_ocn(i)," qsfc=", qsfc_ocn(i)," znt=", znt_ocn(i),& + " ust=", ust_ocn(i)," snowh=", snowh_ocn(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF ENDIF !Use Pedros iterative function to find z/L zol(I)=zolri(rb_ocn(I),ZA(I),ZNTstoch_ocn(I),ZT_ocn(I),ZOL(I)) @@ -1280,8 +1319,14 @@ SUBROUTINE SFCLAY1D_mynn( & ENDIF IF (debug_code >= 1) THEN - write(0,*)"===(dry) capture bad input in mynn sfc layer, i=:",i - write(0,*)"rb=", rb_lnd(I)," ZNT=", ZNTstoch_lnd(i)," ZT=",Zt_lnd(i) + IF (ZNTstoch_lnd(i) < 1E-8 .OR. Zt_lnd(i) < 1E-10) THEN + write(0,*)"===(land) capture bad input in mynn sfc layer, i=:",i + write(0,*)"rb=", rb_lnd(I)," ZNT=", ZNTstoch_lnd(i)," ZT=",Zt_lnd(i) + write(0,*)" tsk=", tskin_lnd(i)," prev z/L=",ZOL(I),& + " tsurf=", tsurf_lnd(i)," qsfc=", qsfc_lnd(i)," znt=", znt_lnd(i),& + " ust=", ust_lnd(i)," snowh=", snowh_lnd(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF ENDIF !Use Pedros iterative function to find z/L zol(I)=zolri(rb_lnd(I),ZA(I),ZNTstoch_lnd(I),ZT_lnd(I),ZOL(I)) @@ -1337,8 +1382,14 @@ SUBROUTINE SFCLAY1D_mynn( & ENDIF IF (debug_code >= 1) THEN - write(0,*)"===(dry) capture bad input in mynn sfc layer, i=:",i - write(0,*)"rb=", rb_lnd(I)," ZNT=", ZNTstoch_lnd(i)," ZT=",Zt_lnd(i) + IF (ZNTstoch_lnd(i) < 1E-8 .OR. Zt_lnd(i) < 1E-10) THEN + write(0,*)"===(land) capture bad input in mynn sfc layer, i=:",i + write(0,*)"rb=", rb_lnd(I)," ZNT=", ZNTstoch_lnd(i)," ZT=",Zt_lnd(i) + write(0,*)" tsk=", tskin_lnd(i)," wstar=",wstar(i)," prev z/L=",ZOL(I),& + " tsurf=", tsurf_lnd(i)," qsfc=", qsfc_lnd(i)," znt=", znt_lnd(i),& + " ust=", ust_lnd(i)," snowh=", snowh_lnd(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF ENDIF !Use Pedros iterative function to find z/L zol(I)=zolri(rb_lnd(I),ZA(I),ZNTstoch_lnd(I),ZT_lnd(I),ZOL(I)) @@ -1397,8 +1448,14 @@ SUBROUTINE SFCLAY1D_mynn( & ENDIF IF (debug_code >= 1) THEN - write(0,*)"===(ice) capture bad input in mynn sfc layer, i=:",i - write(0,*)"rb=", rb_ice(I)," ZNT=", ZNTstoch_ice(i)," ZT=",Zt_ice(i) + IF (ZNTstoch_ice(i) < 1E-8 .OR. Zt_ice(i) < 1E-10) THEN + write(0,*)"===(ice) capture bad input in mynn sfc layer, i=:",i + write(0,*)"rb=", rb_ice(I)," ZNT=", ZNTstoch_ice(i)," ZT=",Zt_ice(i) + write(0,*)" tsk=", tskin_ice(i)," prev z/L=",ZOL(I),& + " tsurf=", tsurf_ice(i)," qsfc=", qsfc_ice(i)," znt=", znt_ice(i),& + " ust=", ust_ice(i)," snowh=", snowh_ice(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF ENDIF !Use Pedros iterative function to find z/L zol(I)=zolri(rb_ice(I),ZA(I),ZNTstoch_ice(I),ZT_ice(I),ZOL(I)) @@ -1454,8 +1511,14 @@ SUBROUTINE SFCLAY1D_mynn( & ENDIF IF (debug_code >= 1) THEN - write(0,*)"===(ice) capture bad input in mynn sfc layer, i=:",i - write(0,*)"rb=", rb_ice(I)," ZNT=", ZNTstoch_ice(i)," ZT=",Zt_ice(i) + IF (ZNTstoch_ice(i) < 1E-8 .OR. Zt_ice(i) < 1E-10) THEN + write(0,*)"===(ice) capture bad input in mynn sfc layer, i=:",i + write(0,*)"rb=", rb_ice(I)," ZNT=", ZNTstoch_ice(i)," ZT=",Zt_ice(i) + write(0,*)" tsk=", tskin_ice(i)," wstar=",wstar(i)," prev z/L=",ZOL(I),& + " tsurf=", tsurf_ice(i)," qsfc=", qsfc_ice(i)," znt=", znt_ice(i),& + " ust=", ust_ice(i)," snowh=", snowh_ice(i),"psfcpa=",PSFCPA(i), & + " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) + ENDIF ENDIF !Use Pedros iterative function to find z/L zol(I)=zolri(rb_ice(I),ZA(I),ZNTstoch_ice(I),ZT_ice(I),ZOL(I)) @@ -1593,9 +1656,9 @@ SUBROUTINE SFCLAY1D_mynn( & IF (debug_code == 2) THEN DO I=its,ite - IF(wet(i))write(*,*)"==== AT END OF ITER LOOP, i=",i, "(wet)" - IF(dry(i))write(*,*)"==== AT END OF ITER LOOP, i=",i, "(land)" - IF(icy(i))write(*,*)"==== AT END OF ITER LOOP, i=",i, "(ice)" + IF(wet(i))write(*,*)"==== AT END OF MAIN LOOP, i=",i, "(wet)" + IF(dry(i))write(*,*)"==== AT END OF MAIN LOOP, i=",i, "(land)" + IF(icy(i))write(*,*)"==== AT END OF MAIN LOOP, i=",i, "(ice)" write(*,*)"z/L:",ZOL(I)," wspd:",wspd(I)," Tstar:",MOL(I) IF(wet(i))write(*,*)"PSIM:",PSIM(I)," PSIH:",PSIH(I)," W*:",WSTAR(I),& " DTHV:",THV1D(I)-THVSK_ocn(I) @@ -1647,20 +1710,23 @@ SUBROUTINE SFCLAY1D_mynn( & FLQC(I)=RHO1D(I)*MAVAIL(I)*UST_lnd(I)*KARMAN/PSIQ_lnd(i) FLHC(I)=RHO1D(I)*CPM(I)*UST_lnd(I)*KARMAN/PSIT_lnd(I) - !---------------------------------- - ! COMPUTE SURFACE MOISTURE FLUX: - !---------------------------------- - QFX(I)=FLQC(I)*(QSFCMR_lnd(I)-QV1D(I)) - QFX(I)=MAX(QFX(I),-0.02) !allows small neg QFX - LH(i)=XLV*QFX(i) - QFLX(i)=QFX(i)/RHO1D(i) - - !---------------------------------- - ! COMPUTE SURFACE HEAT FLUX: - !---------------------------------- - HFX(I)=FLHC(I)*(THSK_lnd(I)-TH1D(I)) - HFX(I)=MAX(HFX(I),-250.) - HFLX(I)=HFX(I)/(RHO1D(I)*cpm(I)) + IF (compute_flux) THEN + !---------------------------------- + ! COMPUTE SURFACE MOISTURE FLUX: + !---------------------------------- + !QFX(I)=FLQC(I)*(QSFCMR_lnd(I)-QV1D(I)) + QFX(I)=FLQC(I)*(QSFC_lnd(I)-QV1D(I)) + QFX(I)=MAX(QFX(I),-0.02) !allows small neg QFX + LH(i)=XLV*QFX(i) + QFLX(i)=QFX(i)/RHO1D(i) + + !---------------------------------- + ! COMPUTE SURFACE HEAT FLUX: + !---------------------------------- + HFX(I)=FLHC(I)*(THSK_lnd(I)-TH1D(I)) + HFX(I)=MAX(HFX(I),-250.) + HFLX(I)=HFX(I)/(RHO1D(I)*cpm(I)) + ENDIF !TRANSFER COEFF FOR SOME LSMs: !CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) & @@ -1682,25 +1748,28 @@ SUBROUTINE SFCLAY1D_mynn( & FLQC(I)=RHO1D(I)*MAVAIL(I)*UST_ocn(I)*KARMAN/PSIQ_ocn(i) FLHC(I)=RHO1D(I)*CPM(I)*UST_ocn(I)*KARMAN/PSIT_ocn(I) - !---------------------------------- - ! COMPUTE SURFACE MOISTURE FLUX: - !---------------------------------- - QFX(I)=FLQC(I)*(QSFCMR_ocn(I)-QV1D(I)) - QFX(I)=MAX(QFX(I),-0.02) !allows small neg QFX - LH(I)=XLV*QFX(I) - QFLX(i)=QFX(i)/RHO1D(i) - - !---------------------------------- - ! COMPUTE SURFACE HEAT FLUX: - !---------------------------------- - HFX(I)=FLHC(I)*(THSK_ocn(I)-TH1D(I)) - IF ( PRESENT(ISFTCFLX) ) THEN - IF ( ISFTCFLX.NE.0 ) THEN - ! AHW: add dissipative heating term - HFX(I)=HFX(I)+RHO1D(I)*USTM(I)*USTM(I)*WSPDI(I) + IF (compute_flux) THEN + !---------------------------------- + ! COMPUTE SURFACE MOISTURE FLUX: + !---------------------------------- + !QFX(I)=FLQC(I)*(QSFCMR_ocn(I)-QV1D(I)) + QFX(I)=FLQC(I)*(QSFC_ocn(I)-QV1D(I)) + QFX(I)=MAX(QFX(I),-0.02) !allows small neg QFX + LH(I)=XLV*QFX(I) + QFLX(i)=QFX(i)/RHO1D(i) + + !---------------------------------- + ! COMPUTE SURFACE HEAT FLUX: + !---------------------------------- + HFX(I)=FLHC(I)*(THSK_ocn(I)-TH1D(I)) + IF ( PRESENT(ISFTCFLX) ) THEN + IF ( ISFTCFLX.NE.0 ) THEN + ! AHW: add dissipative heating term + HFX(I)=HFX(I)+RHO1D(I)*USTM(I)*USTM(I)*WSPDI(I) + ENDIF ENDIF + HFLX(I)=HFX(I)/(RHO1D(I)*cpm(I)) ENDIF - HFLX(I)=HFX(I)/(RHO1D(I)*cpm(I)) !TRANSFER COEFF FOR SOME LSMs: !CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) & @@ -1722,20 +1791,23 @@ SUBROUTINE SFCLAY1D_mynn( & FLQC(I)=RHO1D(I)*MAVAIL(I)*UST_ice(I)*KARMAN/PSIQ_ice(i) FLHC(I)=RHO1D(I)*CPM(I)*UST_ice(I)*KARMAN/PSIT_ice(I) - !---------------------------------- - ! COMPUTE SURFACE MOISTURE FLUX: - !---------------------------------- - QFX(I)=FLQC(I)*(QSFCMR_ice(I)-QV1D(I)) - QFX(I)=MAX(QFX(I),-0.02) !allows small neg QFX - LH(I)=XLF*QFX(I) - QFLX(i)=QFX(i)/RHO1D(i) - - !---------------------------------- - ! COMPUTE SURFACE HEAT FLUX: - !---------------------------------- - HFX(I)=FLHC(I)*(THSK_ice(I)-TH1D(I)) - HFX(I)=MAX(HFX(I),-250.) - HFLX(I)=HFX(I)/(RHO1D(I)*cpm(I)) + IF (compute_flux) THEN + !---------------------------------- + ! COMPUTE SURFACE MOISTURE FLUX: + !---------------------------------- + !QFX(I)=FLQC(I)*(QSFCMR_ice(I)-QV1D(I)) + QFX(I)=FLQC(I)*(QSFC_ice(I)-QV1D(I)) + QFX(I)=MAX(QFX(I),-0.02) !allows small neg QFX + LH(I)=XLF*QFX(I) + QFLX(i)=QFX(i)/RHO1D(i) + + !---------------------------------- + ! COMPUTE SURFACE HEAT FLUX: + !---------------------------------- + HFX(I)=FLHC(I)*(THSK_ice(I)-TH1D(I)) + HFX(I)=MAX(HFX(I),-250.) + HFLX(I)=HFX(I)/(RHO1D(I)*cpm(I)) + ENDIF !TRANSFER COEFF FOR SOME LSMs: !CHS(I)=UST(I)*KARMAN/(ALOG(KARMAN*UST(I)*ZA(I) & @@ -1854,8 +1926,8 @@ SUBROUTINE SFCLAY1D_mynn( & ENDIF T2(I)=TH2(I)*(PSFCPA(I)/100000.)**ROVCP - Q2(I)=QSFCMR_lnd(I)+(QV1D(I)-QSFCMR_lnd(I))*PSIQ2_lnd(i)/PSIQ_lnd(i) - Q2(I)= MAX(Q2(I), MIN(QSFCMR_lnd(I), QV1D(I))) + Q2(I)=QSFC_lnd(I)+(QV1D(I)-QSFC_lnd(I))*PSIQ2_lnd(i)/PSIQ_lnd(i) + Q2(I)= MAX(Q2(I), MIN(QSFC_lnd(I), QV1D(I))) Q2(I)= MIN(Q2(I), 1.05*QV1D(I)) ELSEIF (wet(i)) THEN DTG=TH1D(I)-THSK_ocn(I) @@ -1868,8 +1940,8 @@ SUBROUTINE SFCLAY1D_mynn( & ENDIF T2(I)=TH2(I)*(PSFCPA(I)/100000.)**ROVCP - Q2(I)=QSFCMR_ocn(I)+(QV1D(I)-QSFCMR_ocn(I))*PSIQ2_ocn(i)/PSIQ_ocn(i) - Q2(I)= MAX(Q2(I), MIN(QSFCMR_ocn(I), QV1D(I))) + Q2(I)=QSFC_ocn(I)+(QV1D(I)-QSFC_ocn(I))*PSIQ2_ocn(i)/PSIQ_ocn(i) + Q2(I)= MAX(Q2(I), MIN(QSFC_ocn(I), QV1D(I))) Q2(I)= MIN(Q2(I), 1.05*QV1D(I)) ELSEIF (icy(i)) THEN DTG=TH1D(I)-THSK_ice(I) @@ -1882,8 +1954,8 @@ SUBROUTINE SFCLAY1D_mynn( & ENDIF T2(I)=TH2(I)*(PSFCPA(I)/100000.)**ROVCP - Q2(I)=QSFCMR_ice(I)+(QV1D(I)-QSFCMR_ice(I))*PSIQ2_ice(i)/PSIQ_ice(i) - Q2(I)= MAX(Q2(I), MIN(QSFCMR_ice(I), QV1D(I))) + Q2(I)=QSFC_ice(I)+(QV1D(I)-QSFC_ice(I))*PSIQ2_ice(i)/PSIQ_ice(i) + Q2(I)= MAX(Q2(I), MIN(QSFC_ice(I), QV1D(I))) Q2(I)= MIN(Q2(I), 1.05*QV1D(I)) ENDIF ENDDO @@ -1895,15 +1967,17 @@ SUBROUTINE SFCLAY1D_mynn( & IF ( debug_code == 2) THEN DO I=its,ite yesno = 0 - IF (HFX(I) > 1200. .OR. HFX(I) < -700.)THEN + IF (compute_flux) THEN + IF (HFX(I) > 1200. .OR. HFX(I) < -700.)THEN print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",& I,J, "HFX: ",HFX(I) yesno = 1 - ENDIF - IF (LH(I) > 1200. .OR. LH(I) < -700.)THEN + ENDIF + IF (LH(I) > 1200. .OR. LH(I) < -700.)THEN print*,"SUSPICIOUS VALUES IN MYNN SFCLAYER",& I,J, "LH: ",LH(I) yesno = 1 + ENDIF ENDIF IF (wet(i)) THEN IF (UST_ocn(I) < 0.0 .OR. UST_ocn(I) > 4.0 )THEN @@ -2608,9 +2682,9 @@ SUBROUTINE PSI_CB2005(psim1,psih1,zL,z0L) REAL, INTENT(IN) :: zL,z0L REAL, INTENT(OUT) :: psim1,psih1 - psim1 = -6.1*LOG(zL + (1.+ zL**2.5)**0.4) - & + psim1 = -6.1*LOG(zL + (1.+ zL**2.5)**0.4) & -6.1*LOG(z0L + (1.+ z0L**2.5)**0.4) - psih1 = -5.5*log(zL + (1.+ zL**1.1)**0.90909090909) - & + psih1 = -5.5*log(zL + (1.+ zL**1.1)**0.90909090909) & -5.5*log(z0L + (1.+ z0L**1.1)**0.90909090909) return From afd6481120054a5531bdaa3a8ccbf24c884b7f88 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 1 Apr 2020 15:25:28 -0600 Subject: [PATCH 12/12] Compile physics/module_sf_mynn.F90 with -O1 instead of -O2 to avoid a bug with Intel 18 on hera; add a corresponding note in physics/module_sf_mynn.F90 --- CMakeLists.txt | 8 +++++--- physics/module_sf_mynn.F90 | 12 ++++++++---- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b8d3c3e18..8e6785c71 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -191,15 +191,17 @@ elseif (${CMAKE_Fortran_COMPILER_ID} STREQUAL "Intel") ${CMAKE_CURRENT_SOURCE_DIR}/physics/cu_gf_sh.F90 ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_bl_mynn.F90 ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_MYNNPBL_wrapper.F90 - ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_sf_mynn.F90 ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_MYNNSFC_wrapper.F90 - ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_MYNNrad_pre.F90 - ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_MYNNrad_post.F90 ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_mp_thompson_make_number_concentrations.F90 ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_SF_JSFC.F90 ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_BL_MYJPBL.F90 PROPERTIES COMPILE_FLAGS "-r8 -ftz") + # Reduce optimization for module_sf_mynn.F90 (to avoid an apparent compiler bug with Intel 18 on Hera) + SET_SOURCE_FILES_PROPERTIES(${CMAKE_CURRENT_SOURCE_DIR}/physics/module_sf_mynn.F90 + PROPERTIES COMPILE_FLAGS "${CMAKE_Fortran_FLAGS_OPT} -O1") + list(APPEND SCHEMES_SFX_OPT ${CMAKE_CURRENT_SOURCE_DIR}/physics/module_sf_mynn.F90) + # Replace -xHost or -xCORE-AVX2 with -xCORE-AVX-I for certain files set(CMAKE_Fortran_FLAGS_LOPT1 ${CMAKE_Fortran_FLAGS_OPT}) string(REPLACE "-xHOST" "-xCORE-AVX-I" diff --git a/physics/module_sf_mynn.F90 b/physics/module_sf_mynn.F90 index 2ac9d832c..73ef5e1fb 100644 --- a/physics/module_sf_mynn.F90 +++ b/physics/module_sf_mynn.F90 @@ -645,19 +645,19 @@ SUBROUTINE SFCLAY1D_mynn( & DO I=its,ite write(0,*)"=== imortant input to mynnsfclayer, i:", i IF (dry(i)) THEN - write(0,*)"dry=",dry(i)," tsk=", tskin_lnd(i),& + write(0,*)"dry=",dry(i)," pblh=",pblh(i)," tsk=", tskin_lnd(i),& " tsurf=", tsurf_lnd(i)," qsfc=", qsfc_lnd(i)," znt=", znt_lnd(i),& " ust=", ust_lnd(i)," snowh=", snowh_lnd(i),"psfcpa=",PSFCPA(i), & " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) ENDIF IF (icy(i)) THEN - write(0,*)"icy=",icy(i)," tsk=", tskin_ice(i),& + write(0,*)"icy=",icy(i)," pblh=",pblh(i)," tsk=", tskin_ice(i),& " tsurf=", tsurf_ice(i)," qsfc=", qsfc_ice(i)," znt=", znt_ice(i),& " ust=", ust_ice(i)," snowh=", snowh_ice(i),"psfcpa=",PSFCPA(i), & " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) ENDIF IF (wet(i)) THEN - write(0,*)"wet=",wet(i)," tsk=", tskin_ocn(i),& + write(0,*)"wet=",wet(i)," pblh=",pblh(i)," tsk=", tskin_ocn(i),& " tsurf=", tsurf_ocn(i)," qsfc=", qsfc_ocn(i)," znt=", znt_ocn(i),& " ust=", ust_ocn(i)," snowh=", snowh_ocn(i),"psfcpa=",PSFCPA(i), & " dz=",dz8w1d(i)," qflx=",qflx(i)," hflx=",hflx(i)," hpbl=",pblh(i) @@ -813,7 +813,11 @@ SUBROUTINE SFCLAY1D_mynn( & ENDIF DO I=its,ite - WSPD(I)=SQRT(U1D(I)*U1D(I)+V1D(I)*V1D(I)) + ! DH* 20200401 - note. A weird bug in Intel 18 on hera prevents using the + ! normal -O2 optimization in REPRO and PROD mode for this file. Not reproducible + ! by every user, the bug manifests itself in the resulting wind speed WSPD(I) + ! being -99.0 despite the assignments in lines 932 and 933. *DH + WSPD(I)=SQRT(U1D(I)*U1D(I)+V1D(I)*V1D(I)) WSPD_ocn = -99. WSPD_ice = -99. WSPD_lnd = -99.