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/GFS_PBL_generic.F90 b/physics/GFS_PBL_generic.F90 index 4bebae589..e157013ec 100644 --- a/physics/GFS_PBL_generic.F90 +++ b/physics/GFS_PBL_generic.F90 @@ -150,11 +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,ntoz) - vdftra(i,k,9) = qgrs(i,k,ntwa) - vdftra(i,k,10) = 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 @@ -165,8 +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,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 @@ -404,11 +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,ntoz) = dvdftra(i,k,8) - dqdt(i,k,ntwa) = dvdftra(i,k,9) - dqdt(i,k,ntia) = dvdftra(i,k,10) + 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 @@ -419,8 +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,ntoz) = dvdftra(i,k,7) + 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 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) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index b179a74db..d123c9e4b 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -21,6 +21,7 @@ end subroutine GFS_rrtmg_pre_init subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input Tbd, Cldprop, Coupling, & Radtend, & ! 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 @@ -63,13 +65,20 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input & progcld1, progcld3, & & progcld2, & & progcld4, progcld5, & - & progclduni + & progcld6, progclduni use module_radsw_parameters, only: topfsw_type, sfcfsw_type, & & profsw_type, NBDSW use module_radlw_parameters, only: topflw_type, sfcflw_type, & & proflw_type, NBDLW use surface_perturbation, only: cdfnor + ! For Thompson MP + use module_mp_thompson, only: calc_effectRad, Nt_c + use module_mp_thompson_make_number_concentrations, only: & + make_IceNumber, & + make_DropletNumber, & + make_RainNumber + implicit none type(GFS_control_type), intent(in) :: Model @@ -82,6 +91,7 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input type(GFS_coupling_type), intent(in) :: Coupling 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 +133,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 +152,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, ntwa integer :: i, j, k, k1, k2, lsk, lv, n, itop, ibtc, LP1, lla, llb, lya, lyb @@ -154,7 +164,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, 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 @@ -165,6 +179,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,10 +197,13 @@ 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 ntgl = Model%ntgl + ntwa = Model%ntwa ncndl = min(Model%ncnd,4) LP1 = LM + 1 ! num of in/out levels @@ -256,6 +276,8 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input plyr(i,k1) = Statein%prsl(i,k2) * 0.01 ! pa to mb (hpa) 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 @@ -551,9 +573,36 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input ccnd(i,k,1) = tracer1(i,k,ntcw) ! liquid water 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 + ccnd(i,k,4) = tracer1(i,k,ntsw) + tracer1(i,k,ntgl) ! snow + graupel enddo enddo + ! for Thompson MP - prepare variables for calc_effr + if (Model%imp_physics == Model%imp_physics_thompson .and. Model%ltaerosol) then + do k=1,LMK + do i=1,IM + 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) + nwfa (i,k) = tracer1(i,k,ntwa) + enddo + enddo + elseif (Model%imp_physics == Model%imp_physics_thompson) then + do k=1,LMK + do i=1,IM + 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) = nt_c*orho(i,k1) + ni_mp (i,k) = tracer1(i,k,ntinc)/(1.-qvs) + enddo + enddo + endif endif do n=1,ncndl do k=1,LMK @@ -562,7 +611,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 @@ -612,7 +661,23 @@ 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 .and. Model%kdt>1) THEN + 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 + 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 +699,76 @@ 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 + ! make NC consistent with sub-grid clouds + if (Model%ltaerosol .and. 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 + ! 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 ) + 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 + + do k=1,lm + k1 = k + kd + do i=1,im + 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 +883,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 +900,126 @@ 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 progcld6 to get Xu-Randall total cloud cover (clouds(:,1:LMK,1)) at + ! --- initial time step, it takes into account subgrid PBL + ! --- clouds + call progcld6 (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 progcld6 + 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 progcld6 to get Xu-Randall total cloud cover (clouds(:,1:LMK,1)) + ! tgs: a short subroutine could be made of progcld5 to + ! compute only total cloud fraction. + call progcld6 (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 (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 progcld6 + ! 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 progcld6 + cldcov(i,k) = clouds(i,k,1) + enddo + enddo + endif + + if (.not. clduni) then + ! --- call progcld6 for interaction with the radiation with setting + ! --- uni_cld=.true. to keep precomputed cloud + ! --- fraction + call progcld6 (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..a06e718a5 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -270,6 +270,22 @@ kind = kind_phys intent = out 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 +[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 @@ -430,7 +446,7 @@ dimensions = (horizontal_dimension,vertical_dimension) type = real kind = kind_phys - intent = out + intent = inout optional = F [clouds2] standard_name = cloud_liquid_water_path @@ -439,7 +455,7 @@ dimensions = (horizontal_dimension,vertical_dimension) type = real kind = kind_phys - intent = out + intent = inout optional = F [clouds3] standard_name = mean_effective_radius_for_liquid_cloud @@ -448,7 +464,7 @@ dimensions = (horizontal_dimension,vertical_dimension) type = real kind = kind_phys - intent = out + intent = inout optional = F [clouds4] standard_name = cloud_ice_water_path @@ -457,7 +473,7 @@ dimensions = (horizontal_dimension,vertical_dimension) type = real kind = kind_phys - intent = out + intent = inout optional = F [clouds5] standard_name = mean_effective_radius_for_ice_cloud @@ -466,7 +482,7 @@ dimensions = (horizontal_dimension,vertical_dimension) type = real kind = kind_phys - intent = out + intent = inout optional = F [clouds6] standard_name = cloud_rain_water_path 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 diff --git a/physics/cu_gf_driver.F90 b/physics/cu_gf_driver.F90 index 70d1ce799..4c84621f4 100644 --- a/physics/cu_gf_driver.F90 +++ b/physics/cu_gf_driver.F90 @@ -73,7 +73,7 @@ subroutine cu_gf_driver_run(ntracer,garea,im,ix,km,dt,cactiv, & us,vs,t2di,w,qv2di_spechum,p2di,psuri, & hbot,htop,kcnv,xland,hfx2,qfx2,cliw,clcw, & pbl,ud_mf,dd_mf,dt_mf,cnvw_moist,cnvc,imfshalcnv, & - errmsg,errflg) + qci_conv,errmsg,errflg) !------------------------------------------------------------- implicit none integer, parameter :: maxiens=1 @@ -96,8 +96,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 @@ -743,6 +744,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 0733b603d..9a2e9f230 100644 --- a/physics/cu_gf_driver.meta +++ b/physics/cu_gf_driver.meta @@ -358,6 +358,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_MYNNSFC_wrapper.F90 b/physics/module_MYNNSFC_wrapper.F90 index 951d7e7c8..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 @@ -168,31 +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 - 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:" @@ -230,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) @@ -251,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, & @@ -260,23 +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 ?? - 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,*) 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_MYNNrad_post.F90 b/physics/module_MYNNrad_post.F90 deleted file mode 100644 index 1364db62e..000000000 --- a/physics/module_MYNNrad_post.F90 +++ /dev/null @@ -1,75 +0,0 @@ -!> \file module_MYNNrad_post.F90 -!! Contains the post (interstitial) work after the call to the radiation schemes: -!! 1) Restores the original qc & qi - - MODULE mynnrad_post - - contains - - subroutine mynnrad_post_init () - end subroutine mynnrad_post_init - - subroutine mynnrad_post_finalize () - end subroutine mynnrad_post_finalize - -!>\defgroup gsd_mynnrad_post GSD mynnrad_post_run Module -!>\ingroup gsd_mynn_edmf -!! This interstitial code restores the original resolved-scale clouds (qc and qi). -#if 0 -!! \section arg_table_mynnrad_post_run Argument Table -!! \htmlinclude mynnrad_post_run.html -!! -#endif -SUBROUTINE mynnrad_post_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 mynnrad_post_run - -!###================================================================= - -END MODULE mynnrad_post diff --git a/physics/module_MYNNrad_pre.F90 b/physics/module_MYNNrad_pre.F90 deleted file mode 100644 index 95dc95445..000000000 --- a/physics/module_MYNNrad_pre.F90 +++ /dev/null @@ -1,131 +0,0 @@ -!> \file module_MYNNrad_pre.F90 -!! Contains the preliminary (interstitial) work to the call to the radiation schemes: -!! 1) Backs up the original qc & qi -!! 2) Adds the subgrid clouds mixing ratio and cloud fraction to the original qc, qi and cloud fraction coming from the microphysics scheme. - - MODULE mynnrad_pre - - contains - - subroutine mynnrad_pre_init () - end subroutine mynnrad_pre_init - - subroutine mynnrad_pre_finalize () - end subroutine mynnrad_pre_finalize - -!>\defgroup gsd_mynnrad_pre GSD mynnrad_pre_run Module -!>\ingroup gsd_mynn_edmf -!! This interstitial code adds the subgrid clouds to the resolved-scale clouds if there is no resolved-scale clouds in that particular grid box. -#if 0 -!> \section arg_table_mynnrad_pre_run Argument Table -!! \htmlinclude mynnrad_pre_run.html -!! -#endif -! -! 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 ! -! -!###=================================================================== -SUBROUTINE mynnrad_pre_run( & - & ix,im,levs, & - & flag_init,flag_restart, & - & qc, qi, T3D, & - & qc_save, qi_save, & - & qc_bl,cldfra_bl, & - & delp,clouds1,clouds2,clouds3, & - & clouds4,clouds5,slmsk, & - & 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 - -!------------------------------------------------------------------- - implicit none -!------------------------------------------------------------------- - ! Interface variables - real (kind=kind_phys), parameter :: gfac=1.0e5/con_g - integer, intent(in) :: ix, im, levs - logical, intent(in) :: flag_init, flag_restart - real(kind=kind_phys), dimension(im,levs), intent(inout) :: qc, qi - 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 - character(len=*), intent(out) :: errmsg - integer, intent(out) :: errflg - ! Local variables - integer :: i, k - real :: Tc, iwc - - ! 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 - ! Add subgrid cloud information: - do k = 1, levs - do i = 1, im - - qc_save(i,k) = qc(i,k) - qi_save(i,k) = qi(i,k) - clouds1(i,k) = CLDFRA_BL(i,k) - - 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)-254.)/15.)))*CLDFRA_BL(i,k) - qi(i,k) = QC_BL(i,k)*(1. - MIN(1., MAX(0., (T3D(i,k)-254.)/15.)))*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 - - !water and ice paths - 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 - - !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 mynnrad_pre_run - -!###================================================================= - -END MODULE mynnrad_pre diff --git a/physics/module_SGSCloud_RadPost.F90 b/physics/module_SGSCloud_RadPost.F90 new file mode 100644 index 000000000..051033a26 --- /dev/null +++ b/physics/module_SGSCloud_RadPost.F90 @@ -0,0 +1,69 @@ +!> \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). +!! \section arg_table_sgscloud_radpost_run Argument Table +!! \htmlinclude sgscloud_radpost_run.html +!! + 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(inout) :: 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 sgscloud 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_MYNNrad_post.meta b/physics/module_SGSCloud_RadPost.meta similarity index 97% rename from physics/module_MYNNrad_post.meta rename to physics/module_SGSCloud_RadPost.meta index 881a19fff..f56bd30df 100644 --- a/physics/module_MYNNrad_post.meta +++ b/physics/module_SGSCloud_RadPost.meta @@ -1,5 +1,5 @@ [ccpp-arg-table] - name = mynnrad_post_run + name = sgscloud_radpost_run type = scheme [ix] standard_name = horizontal_dimension @@ -48,7 +48,7 @@ dimensions = (horizontal_dimension,vertical_dimension) type = real kind = kind_phys - intent = out + intent = inout optional = F [qi] standard_name = ice_water_mixing_ratio @@ -57,7 +57,7 @@ dimensions = (horizontal_dimension,vertical_dimension) type = real kind = kind_phys - intent = out + intent = inout optional = F [qc_save] standard_name = cloud_condensed_water_mixing_ratio_save diff --git a/physics/module_SGSCloud_RadPre.F90 b/physics/module_SGSCloud_RadPre.F90 new file mode 100644 index 000000000..544fe1004 --- /dev/null +++ b/physics/module_SGSCloud_RadPre.F90 @@ -0,0 +1,206 @@ +!>\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, imfdeepcnv_gf, & + 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 + 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, imfdeepcnv_gf, 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 + ! 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 + 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 + integer, dimension(im,3), intent(inout) :: 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) 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 == imfdeepcnv_gf) then + do k = 1, levs + do i = 1, im + if ( 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_MYNNrad_pre.meta b/physics/module_SGSCloud_RadPre.meta similarity index 57% rename from physics/module_MYNNrad_pre.meta rename to physics/module_SGSCloud_RadPre.meta index 3b5943c66..507f4ba91 100644 --- a/physics/module_MYNNrad_pre.meta +++ b/physics/module_SGSCloud_RadPre.meta @@ -1,5 +1,15 @@ [ccpp-arg-table] - name = mynnrad_pre_run + 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 @@ -68,6 +78,49 @@ 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 +[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 [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 @@ -75,7 +128,7 @@ dimensions = (horizontal_dimension,vertical_dimension) type = real kind = kind_phys - intent = out + intent = inout optional = F [qi_save] standard_name = ice_water_mixing_ratio_save @@ -84,7 +137,7 @@ dimensions = (horizontal_dimension,vertical_dimension) type = real kind = kind_phys - intent = out + intent = inout optional = F [QC_BL] standard_name = subgrid_cloud_mixing_ratio_pbl @@ -111,7 +164,7 @@ dimensions = (horizontal_dimension,adjusted_vertical_layer_dimension_for_radiation) type = real kind = kind_phys - intent = out + intent = in optional = F [clouds1] standard_name = total_cloud_fraction @@ -167,6 +220,83 @@ 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 = inout + 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 = inout + 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 = inout + 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 diff --git a/physics/module_sf_mynn.F90 b/physics/module_sf_mynn.F90 index 788ff0ace..73ef5e1fb 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)," 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)," 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)," 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) + ENDIF ENDDO ENDIF @@ -786,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. @@ -1161,8 +1192,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 +1256,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 +1323,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 +1386,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 +1452,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 +1515,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 +1660,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 +1714,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 +1752,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 +1795,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 +1930,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 +1944,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 +1958,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 +1971,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 +2686,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 diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 8c341d05b..cda06605e 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -215,294 +215,293 @@ subroutine mp_thompson_init(Data, ntqv, ntcw, ntrw, ntiw, ntsw, ntgl, & nblocks = size(Data) block_loop: do blkno=1,nblocks - ! associate_arrays: associate( & - spechum => Data(blkno)%Statein%qgrs(:,:,ntqv) !,& - qc => Data(blkno)%Statein%qgrs(:,:,ntcw) !,& - qr => Data(blkno)%Statein%qgrs(:,:,ntrw) !,& - qi => Data(blkno)%Statein%qgrs(:,:,ntiw) !,& - qs => Data(blkno)%Statein%qgrs(:,:,ntsw) !,& - qg => Data(blkno)%Statein%qgrs(:,:,ntgl) !,& - ni => Data(blkno)%Statein%qgrs(:,:,ntinc)!,& - nr => Data(blkno)%Statein%qgrs(:,:,ntrnc)!,& - nc => Data(blkno)%Statein%qgrs(:,:,ntlnc)!,& - nwfa => Data(blkno)%Statein%qgrs(:,:,ntwa) !,& - nifa => Data(blkno)%Statein%qgrs(:,:,ntia) !,& - nwfa2d => Data(blkno)%Coupling%nwfa2d !,& - nifa2d => Data(blkno)%Coupling%nifa2d !,& - tgrs => Data(blkno)%Statein%tgrs !,& - prsl => Data(blkno)%Statein%prsl !,& - phil => Data(blkno)%Statein%phil !,& - area => Data(blkno)%Grid%area !,& - re_cloud => Data(blkno)%Tbd%phy_f3d(:,:,nleffr)!,& - re_ice => Data(blkno)%Tbd%phy_f3d(:,:,nieffr)!,& - re_snow => Data(blkno)%Tbd%phy_f3d(:,:,nseffr)! ) - - ncol = size(spechum(:,1)) - nlev = size(spechum(1,:)) - allocate(qv_mp(ncol,nlev)) - allocate(qc_mp(ncol,nlev)) - allocate(qr_mp(ncol,nlev)) - allocate(qi_mp(ncol,nlev)) - allocate(qs_mp(ncol,nlev)) - allocate(qg_mp(ncol,nlev)) - allocate(ni_mp(ncol,nlev)) - allocate(nr_mp(ncol,nlev)) - allocate(nc_mp(ncol,nlev)) - allocate(hgt (ncol,nlev)) - allocate(rho (ncol,nlev)) - allocate(orho (ncol,nlev)) - - only_for_first_block: if (blkno==1) then - - ! Call Thompson init - if (is_aerosol_aware) then - call thompson_init(nwfa2d=nwfa2d, nifa2d=nifa2d, nwfa=nwfa, nifa=nifa, & - mpicomm=mpicomm, mpirank=mpirank, mpiroot=mpiroot, & - threads=threads, errmsg=errmsg, errflg=errflg) - if (errflg /= 0) return - else - call thompson_init(mpicomm=mpicomm, mpirank=mpirank, mpiroot=mpiroot, & - threads=threads, errmsg=errmsg, errflg=errflg) - if (errflg /= 0) return - end if - - ! For restart runs, the init is done here - if (restart) then - is_initialized = .true. - return - end if - - end if only_for_first_block - - ! Fix initial values of hydrometeors - where(spechum<0) spechum = 0.0 - where(qc<0) qc = 0.0 - where(qr<0) qr = 0.0 - where(qi<0) qi = 0.0 - where(qs<0) qs = 0.0 - where(qg<0) qg = 0.0 - where(ni<0) ni = 0.0 - where(nr<0) nr = 0.0 - - if (is_aerosol_aware) then - ! Fix initial values of aerosols - where(nc<0) nc = 0.0 - where(nwfa<0) nwfa = 0.0 - where(nifa<0) nifa = 0.0 - where(nwfa2d<0) nwfa2d = 0.0 - where(nifa2d<0) nifa2d = 0.0 - end if - - ! Geopotential height in m2 s-2 to height in m - hgt = phil/con_g - - ! Density of air in kg m-3 and inverse density of air - rho = prsl/(con_rd*tgrs) - orho = 1.0/rho - - ! Prior to calling the functions: make_DropletNumber, make_IceNumber, make_RainNumber, - ! the incoming mixing ratios should be converted to units of mass/num per cubic meter - ! rather than per kg of air. So, to pass back to the model state variables, - ! they also need to be switched back to mass/number per kg of air, because - ! what is returned by the functions is in units of number per cubic meter. - ! They also need to be converted to dry mixing ratios. - - !> - Convert specific humidity/moist mixing ratios to dry mixing ratios - qv_mp = spechum/(1.0_kind_phys-spechum) - qc_mp = qc/(1.0_kind_phys-spechum) - qr_mp = qr/(1.0_kind_phys-spechum) - qi_mp = qi/(1.0_kind_phys-spechum) - qs_mp = qs/(1.0_kind_phys-spechum) - qg_mp = qg/(1.0_kind_phys-spechum) - - !> - Convert number concentrations from moist to dry - ni_mp = ni/(1.0_kind_phys-spechum) - nr_mp = nr/(1.0_kind_phys-spechum) + spechum => Data(blkno)%Statein%qgrs(:,:,ntqv) + qc => Data(blkno)%Statein%qgrs(:,:,ntcw) + qr => Data(blkno)%Statein%qgrs(:,:,ntrw) + qi => Data(blkno)%Statein%qgrs(:,:,ntiw) + qs => Data(blkno)%Statein%qgrs(:,:,ntsw) + qg => Data(blkno)%Statein%qgrs(:,:,ntgl) + ni => Data(blkno)%Statein%qgrs(:,:,ntinc) + nr => Data(blkno)%Statein%qgrs(:,:,ntrnc) + if (is_aerosol_aware) then + nc => Data(blkno)%Statein%qgrs(:,:,ntlnc) + nwfa => Data(blkno)%Statein%qgrs(:,:,ntwa) + nifa => Data(blkno)%Statein%qgrs(:,:,ntia) + nwfa2d => Data(blkno)%Coupling%nwfa2d + nifa2d => Data(blkno)%Coupling%nifa2d + end if + tgrs => Data(blkno)%Statein%tgrs + prsl => Data(blkno)%Statein%prsl + phil => Data(blkno)%Statein%phil + area => Data(blkno)%Grid%area + re_cloud => Data(blkno)%Tbd%phy_f3d(:,:,nleffr) + re_ice => Data(blkno)%Tbd%phy_f3d(:,:,nieffr) + re_snow => Data(blkno)%Tbd%phy_f3d(:,:,nseffr) + + ncol = size(spechum(:,1)) + nlev = size(spechum(1,:)) + allocate(qv_mp(ncol,nlev)) + allocate(qc_mp(ncol,nlev)) + allocate(qr_mp(ncol,nlev)) + allocate(qi_mp(ncol,nlev)) + allocate(qs_mp(ncol,nlev)) + allocate(qg_mp(ncol,nlev)) + allocate(ni_mp(ncol,nlev)) + allocate(nr_mp(ncol,nlev)) + if (is_aerosol_aware) allocate(nc_mp(ncol,nlev)) + allocate(hgt (ncol,nlev)) + allocate(rho (ncol,nlev)) + allocate(orho (ncol,nlev)) + + only_for_first_block: if (blkno==1) then + + ! Call Thompson init if (is_aerosol_aware) then - nc_mp = nc/(1.0_kind_phys-spechum) + call thompson_init(nwfa2d=nwfa2d, nifa2d=nifa2d, nwfa=nwfa, nifa=nifa, & + mpicomm=mpicomm, mpirank=mpirank, mpiroot=mpiroot, & + threads=threads, errmsg=errmsg, errflg=errflg) + if (errflg /= 0) return + else + call thompson_init(mpicomm=mpicomm, mpirank=mpirank, mpiroot=mpiroot, & + threads=threads, errmsg=errmsg, errflg=errflg) + if (errflg /= 0) return end if - ! If qi is in boundary conditions but ni is not, calculate ni from qi, rho and tgrs - if (maxval(qi_mp)>0.0 .and. maxval(ni_mp)==0.0) then - ni_mp = make_IceNumber(qi_mp*rho, tgrs) * orho + ! For restart runs, the init is done here + if (restart) then + is_initialized = .true. + return end if - ! If ni is in boundary conditions but qi is not, reset ni to zero - if (maxval(ni_mp)>0.0 .and. maxval(qi_mp)==0.0) ni_mp = 0.0 - - ! If qr is in boundary conditions but nr is not, calculate nr from qr, rho and tgrs - if (maxval(qr_mp)>0.0 .and. maxval(nr_mp)==0.0) then - nr_mp = make_RainNumber(qr_mp*rho, tgrs) * orho - end if - - ! If nr is in boundary conditions but qr is not, reset nr to zero - if (maxval(nr_mp)>0.0 .and. maxval(qr_mp)==0.0) nr_mp = 0.0 - - !..Check for existing aerosol data, both CCN and IN aerosols. If missing - !.. fill in just a basic vertical profile, somewhat boundary-layer following. - if (is_aerosol_aware) then - - ! CCN - if (MAXVAL(nwfa) .lt. eps) then - if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently there are no initial CCN aerosols.' - do i = 1, ncol - if (hgt(i,1).le.1000.0) then - h_01 = 0.8 - elseif (hgt(i,1).ge.2500.0) then - h_01 = 0.01 - else - h_01 = 0.8*cos(hgt(i,1)*0.001 - 1.0) - endif - niCCN3 = -1.0*ALOG(naCCN1/naCCN0)/h_01 - nwfa(i,1) = naCCN1+naCCN0*exp(-((hgt(i,2)-hgt(i,1))/1000.)*niCCN3) - airmass = 1./orho(i,1) * (hgt(i,2)-hgt(i,1))*area(i) ! kg - nwfa2d(i) = nwfa(i,1) * 0.000196 * (airmass*2.E-10) - do k = 2, nlev - nwfa(i,k) = naCCN1+naCCN0*exp(-((hgt(i,k)-hgt(i,1))/1000.)*niCCN3) - enddo + end if only_for_first_block + + ! Fix initial values of hydrometeors + where(spechum<0) spechum = 0.0 + where(qc<0) qc = 0.0 + where(qr<0) qr = 0.0 + where(qi<0) qi = 0.0 + where(qs<0) qs = 0.0 + where(qg<0) qg = 0.0 + where(ni<0) ni = 0.0 + where(nr<0) nr = 0.0 + + if (is_aerosol_aware) then + ! Fix initial values of aerosols + where(nc<0) nc = 0.0 + where(nwfa<0) nwfa = 0.0 + where(nifa<0) nifa = 0.0 + where(nwfa2d<0) nwfa2d = 0.0 + where(nifa2d<0) nifa2d = 0.0 + end if + + ! Geopotential height in m2 s-2 to height in m + hgt = phil/con_g + + ! Density of air in kg m-3 and inverse density of air + rho = prsl/(con_rd*tgrs) + orho = 1.0/rho + + ! Prior to calling the functions: make_DropletNumber, make_IceNumber, make_RainNumber, + ! the incoming mixing ratios should be converted to units of mass/num per cubic meter + ! rather than per kg of air. So, to pass back to the model state variables, + ! they also need to be switched back to mass/number per kg of air, because + ! what is returned by the functions is in units of number per cubic meter. + ! They also need to be converted to dry mixing ratios. + + !> - Convert specific humidity/moist mixing ratios to dry mixing ratios + qv_mp = spechum/(1.0_kind_phys-spechum) + qc_mp = qc/(1.0_kind_phys-spechum) + qr_mp = qr/(1.0_kind_phys-spechum) + qi_mp = qi/(1.0_kind_phys-spechum) + qs_mp = qs/(1.0_kind_phys-spechum) + qg_mp = qg/(1.0_kind_phys-spechum) + + !> - Convert number concentrations from moist to dry + ni_mp = ni/(1.0_kind_phys-spechum) + nr_mp = nr/(1.0_kind_phys-spechum) + if (is_aerosol_aware) then + nc_mp = nc/(1.0_kind_phys-spechum) + end if + + ! If qi is in boundary conditions but ni is not, calculate ni from qi, rho and tgrs + if (maxval(qi_mp)>0.0 .and. maxval(ni_mp)==0.0) then + ni_mp = make_IceNumber(qi_mp*rho, tgrs) * orho + end if + + ! If ni is in boundary conditions but qi is not, reset ni to zero + if (maxval(ni_mp)>0.0 .and. maxval(qi_mp)==0.0) ni_mp = 0.0 + + ! If qr is in boundary conditions but nr is not, calculate nr from qr, rho and tgrs + if (maxval(qr_mp)>0.0 .and. maxval(nr_mp)==0.0) then + nr_mp = make_RainNumber(qr_mp*rho, tgrs) * orho + end if + + ! If nr is in boundary conditions but qr is not, reset nr to zero + if (maxval(nr_mp)>0.0 .and. maxval(qr_mp)==0.0) nr_mp = 0.0 + + !..Check for existing aerosol data, both CCN and IN aerosols. If missing + !.. fill in just a basic vertical profile, somewhat boundary-layer following. + if (is_aerosol_aware) then + + ! CCN + if (MAXVAL(nwfa) .lt. eps) then + if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently there are no initial CCN aerosols.' + do i = 1, ncol + if (hgt(i,1).le.1000.0) then + h_01 = 0.8 + elseif (hgt(i,1).ge.2500.0) then + h_01 = 0.01 + else + h_01 = 0.8*cos(hgt(i,1)*0.001 - 1.0) + endif + niCCN3 = -1.0*ALOG(naCCN1/naCCN0)/h_01 + nwfa(i,1) = naCCN1+naCCN0*exp(-((hgt(i,2)-hgt(i,1))/1000.)*niCCN3) + airmass = 1./orho(i,1) * (hgt(i,2)-hgt(i,1))*area(i) ! kg + nwfa2d(i) = nwfa(i,1) * 0.000196 * (airmass*2.E-10) + do k = 2, nlev + nwfa(i,k) = naCCN1+naCCN0*exp(-((hgt(i,k)-hgt(i,1))/1000.)*niCCN3) enddo - else - if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently initial CCN aerosols are present.' - if (MAXVAL(nwfa2d) .lt. eps) then + enddo + else + if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently initial CCN aerosols are present.' + if (MAXVAL(nwfa2d) .lt. eps) then ! Hard-coded switch between new (from WRFv4.0, top) and old (until WRFv3.9.1.1, bottom) surface emission rate calculations #if 0 - !+---+-----------------------------------------------------------------+ - !..Scale the lowest level aerosol data into an emissions rate. This is - !.. very far from ideal, but need higher emissions where larger amount - !.. of (climo) existing and lesser emissions where there exists fewer to - !.. begin as a first-order simplistic approach. Later, proper connection to - !.. emission inventory would be better, but, for now, scale like this: - !.. where: Nwfa=50 per cc, emit 0.875E4 aerosols per second per grid box unit - !.. that was tested as ~(20kmx20kmx50m = 2.E10 m**-3) - !+---+-----------------------------------------------------------------+ - if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently there are no initial CCN aerosol surface emission rates.' - if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Use new (WRFv4+) formula to calculate CCN surface emission rates.' - do i = 1, ncol - airmass = 1./orho(i,1) * (hgt(i,2)-hgt(i,1))*area(i) ! kg - nwfa2d(i) = nwfa(i,1) * 0.000196 * (airmass*2.E-10) - enddo + !+---+-----------------------------------------------------------------+ + !..Scale the lowest level aerosol data into an emissions rate. This is + !.. very far from ideal, but need higher emissions where larger amount + !.. of (climo) existing and lesser emissions where there exists fewer to + !.. begin as a first-order simplistic approach. Later, proper connection to + !.. emission inventory would be better, but, for now, scale like this: + !.. where: Nwfa=50 per cc, emit 0.875E4 aerosols per second per grid box unit + !.. that was tested as ~(20kmx20kmx50m = 2.E10 m**-3) + !+---+-----------------------------------------------------------------+ + if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently there are no initial CCN aerosol surface emission rates.' + if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Use new (WRFv4+) formula to calculate CCN surface emission rates.' + do i = 1, ncol + airmass = 1./orho(i,1) * (hgt(i,2)-hgt(i,1))*area(i) ! kg + nwfa2d(i) = nwfa(i,1) * 0.000196 * (airmass*2.E-10) + enddo #else - !+---+-----------------------------------------------------------------+ - !..Scale the lowest level aerosol data into an emissions rate. This is - !.. very far from ideal, but need higher emissions where larger amount - !.. of existing and lesser emissions where not already lots of aerosols - !.. for first-order simplistic approach. Later, proper connection to - !.. emission inventory would be better, but, for now, scale like this: - !.. where: Nwfa=50 per cc, emit 0.875E4 aerosols per kg per second - !.. Nwfa=500 per cc, emit 0.875E5 aerosols per kg per second - !.. Nwfa=5000 per cc, emit 0.875E6 aerosols per kg per second - !.. for a grid with 20km spacing and scale accordingly for other spacings. - !+---+-----------------------------------------------------------------+ - if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently there are no initial CCN aerosol surface emission rates.' - if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Use old (pre WRFv4) formula to calculate CCN surface emission rates.' - do i = 1, ncol - if (SQRT(area(i))/20000.0 .ge. 1.0) then - h_01 = 0.875 - else - h_01 = (0.875 + 0.125*((20000.-SQRT(area(i)))/16000.)) * SQRT(area(i))/20000. - endif - nwfa2d(i) = 10.0**(LOG10(nwfa(i,1)*1.E-6)-3.69897) - nwfa2d(i) = nwfa2d(i)*h_01 * 1.E6 - enddo -#endif - else - if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently initial CCN aerosol surface emission rates are present.' - endif - endif - - ! IN - if (MAXVAL(nifa) .lt. eps) then - if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently there are no initial IN aerosols.' + !+---+-----------------------------------------------------------------+ + !..Scale the lowest level aerosol data into an emissions rate. This is + !.. very far from ideal, but need higher emissions where larger amount + !.. of existing and lesser emissions where not already lots of aerosols + !.. for first-order simplistic approach. Later, proper connection to + !.. emission inventory would be better, but, for now, scale like this: + !.. where: Nwfa=50 per cc, emit 0.875E4 aerosols per kg per second + !.. Nwfa=500 per cc, emit 0.875E5 aerosols per kg per second + !.. Nwfa=5000 per cc, emit 0.875E6 aerosols per kg per second + !.. for a grid with 20km spacing and scale accordingly for other spacings. + !+---+-----------------------------------------------------------------+ + if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently there are no initial CCN aerosol surface emission rates.' + if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Use old (pre WRFv4) formula to calculate CCN surface emission rates.' do i = 1, ncol - if (hgt(i,1).le.1000.0) then - h_01 = 0.8 - elseif (hgt(i,1).ge.2500.0) then - h_01 = 0.01 - else - h_01 = 0.8*cos(hgt(i,1)*0.001 - 1.0) - endif - niIN3 = -1.0*ALOG(naIN1/naIN0)/h_01 - nifa(i,1) = naIN1+naIN0*exp(-((hgt(i,2)-hgt(i,1))/1000.)*niIN3) - nifa2d(i) = 0. - do k = 2, nlev - nifa(i,k) = naIN1+naIN0*exp(-((hgt(i,k)-hgt(i,1))/1000.)*niIN3) - enddo + if (SQRT(area(i))/20000.0 .ge. 1.0) then + h_01 = 0.875 + else + h_01 = (0.875 + 0.125*((20000.-SQRT(area(i)))/16000.)) * SQRT(area(i))/20000. + endif + nwfa2d(i) = 10.0**(LOG10(nwfa(i,1)*1.E-6)-3.69897) + nwfa2d(i) = nwfa2d(i)*h_01 * 1.E6 enddo +#endif else - if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently initial IN aerosols are present.' - if (MAXVAL(nifa2d) .lt. eps) then - if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently there are no initial IN aerosol surface emission rates, set to zero.' - ! calculate IN surface flux here, right now just set to zero - nifa2d = 0. + if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently initial CCN aerosol surface emission rates are present.' + endif + endif + + ! IN + if (MAXVAL(nifa) .lt. eps) then + if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently there are no initial IN aerosols.' + do i = 1, ncol + if (hgt(i,1).le.1000.0) then + h_01 = 0.8 + elseif (hgt(i,1).ge.2500.0) then + h_01 = 0.01 else - if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently initial IN aerosol surface emission rates are present.' + h_01 = 0.8*cos(hgt(i,1)*0.001 - 1.0) endif + niIN3 = -1.0*ALOG(naIN1/naIN0)/h_01 + nifa(i,1) = naIN1+naIN0*exp(-((hgt(i,2)-hgt(i,1))/1000.)*niIN3) + nifa2d(i) = 0. + do k = 2, nlev + nifa(i,k) = naIN1+naIN0*exp(-((hgt(i,k)-hgt(i,1))/1000.)*niIN3) + enddo + enddo + else + if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently initial IN aerosols are present.' + if (MAXVAL(nifa2d) .lt. eps) then + if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently there are no initial IN aerosol surface emission rates, set to zero.' + ! calculate IN surface flux here, right now just set to zero + nifa2d = 0. + else + if (mpirank==mpiroot .and. blkno==1) write(*,*) ' Apparently initial IN aerosol surface emission rates are present.' endif + endif - ! If qc is in boundary conditions but nc is not, calculate nc from qc, rho and nwfa - if (maxval(qc_mp)>0.0 .and. maxval(nc_mp)==0.0) then - nc_mp = make_DropletNumber(qc_mp*rho, nwfa) * orho - end if + ! If qc is in boundary conditions but nc is not, calculate nc from qc, rho and nwfa + if (maxval(qc_mp)>0.0 .and. maxval(nc_mp)==0.0) then + nc_mp = make_DropletNumber(qc_mp*rho, nwfa) * orho + end if - ! If nc is in boundary conditions but qc is not, reset nc to zero - if (maxval(nc_mp)>0.0 .and. maxval(qc_mp)==0.0) nc_mp = 0.0 + ! If nc is in boundary conditions but qc is not, reset nc to zero + if (maxval(nc_mp)>0.0 .and. maxval(qc_mp)==0.0) nc_mp = 0.0 - else + else - ! Constant droplet concentration for single moment cloud water as in - ! module_mp_thompson.F90, only needed for effective radii calculation - nc_mp = Nt_c/rho + ! Constant droplet concentration for single moment cloud water as in + ! module_mp_thompson.F90, only needed for effective radii calculation + nc_mp = Nt_c/rho - end if + end if - ! Calculate initial cloud effective radii if requested - do i = 1, ncol - do k = 1, nlev - re_cloud(i,k) = 2.49E-6 - re_ice(i,k) = 4.99E-6 - re_snow(i,k) = 9.99E-6 - end do - end do - do i = 1, ncol - call calc_effectRad (tgrs(i,:), prsl(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, nlev) + ! Calculate initial cloud effective radii if requested + do i = 1, ncol + do k = 1, nlev + re_cloud(i,k) = 2.49E-6 + re_ice(i,k) = 4.99E-6 + re_snow(i,k) = 9.99E-6 end do - do i = 1, ncol - do k = 1, nlev - re_cloud(i,k) = MAX(2.49E-6, MIN(re_cloud(i,k), 50.E-6)) - re_ice(i,k) = MAX(4.99E-6, MIN(re_ice(i,k), 125.E-6)) - re_snow(i,k) = MAX(9.99E-6, MIN(re_snow(i,k), 999.E-6)) - end do + end do + do i = 1, ncol + call calc_effectRad (tgrs(i,:), prsl(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, nlev) + end do + do i = 1, ncol + do k = 1, nlev + re_cloud(i,k) = MAX(2.49E-6, MIN(re_cloud(i,k), 50.E-6)) + re_ice(i,k) = MAX(4.99E-6, MIN(re_ice(i,k), 125.E-6)) + re_snow(i,k) = MAX(9.99E-6, MIN(re_snow(i,k), 999.E-6)) end do - ! Convert to micron: required for bit-for-bit identical restarts; - ! otherwise entering mp_thompson_init and converting mu to m and - ! back (without updating re_*) introduces b4b differences. - re_cloud = 1.0E6*re_cloud - re_ice = 1.0E6*re_ice - re_snow = 1.0E6*re_snow - - !> - Convert number concentrations from dry to moist - ni = ni_mp/(1.0_kind_phys+qv_mp) - nr = nr_mp/(1.0_kind_phys+qv_mp) - if (is_aerosol_aware) then - nc = nc_mp/(1.0_kind_phys+qv_mp) - end if - - deallocate(qv_mp) - deallocate(qc_mp) - deallocate(qr_mp) - deallocate(qi_mp) - deallocate(qs_mp) - deallocate(qg_mp) - deallocate(ni_mp) - deallocate(nr_mp) - deallocate(nc_mp) - deallocate(hgt ) - deallocate(rho ) - deallocate(orho ) - - !end associate associate_arrays + end do + ! Convert to micron: required for bit-for-bit identical restarts; + ! otherwise entering mp_thompson_init and converting mu to m and + ! back (without updating re_*) introduces b4b differences. + re_cloud = 1.0E6*re_cloud + re_ice = 1.0E6*re_ice + re_snow = 1.0E6*re_snow + + !> - Convert number concentrations from dry to moist + ni = ni_mp/(1.0_kind_phys+qv_mp) + nr = nr_mp/(1.0_kind_phys+qv_mp) + if (is_aerosol_aware) then + nc = nc_mp/(1.0_kind_phys+qv_mp) + end if + + deallocate(qv_mp) + deallocate(qc_mp) + deallocate(qr_mp) + deallocate(qi_mp) + deallocate(qs_mp) + deallocate(qg_mp) + deallocate(ni_mp) + deallocate(nr_mp) + if (is_aerosol_aware) deallocate(nc_mp) + deallocate(hgt ) + deallocate(rho ) + deallocate(orho ) end do block_loop @@ -548,11 +547,12 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & real(kind_phys), intent(inout) :: nr(1:ncol,1:nlev) ! Aerosols logical, intent(in) :: is_aerosol_aware - real(kind_phys), optional, intent(inout) :: nc(1:ncol,1:nlev) - real(kind_phys), optional, intent(inout) :: nwfa(1:ncol,1:nlev) - real(kind_phys), optional, intent(inout) :: nifa(1:ncol,1:nlev) - real(kind_phys), optional, intent(in ) :: nwfa2d(1:ncol) - real(kind_phys), optional, intent(in ) :: nifa2d(1:ncol) + ! The following arrays are not allocated if is_aerosol_aware is false + real(kind_phys), optional, intent(inout) :: nc(:,:) + real(kind_phys), optional, intent(inout) :: nwfa(:,:) + real(kind_phys), optional, intent(inout) :: nifa(:,:) + real(kind_phys), optional, intent(in ) :: nwfa2d(:) + real(kind_phys), optional, intent(in ) :: nifa2d(:) ! State variables and timestep information real(kind_phys), intent(inout) :: tgrs(1:ncol,1:nlev) real(kind_phys), intent(in ) :: prsl(1:ncol,1:nlev) diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 49b394fe1..176219199 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, progcld6, progcld4o, gethml ! ================= @@ -2683,6 +2683,304 @@ subroutine progcld5 & end subroutine progcld5 !................................... +!----------------------------------- +!> \ingroup module_radiation_clouds +!! This subroutine computes cloud related quantities using the Thompson +!! cloud microphysics scheme with updated microphysics-cloud-radiation +!! interaction (including subgrid clouds). Adapted from progcld5. + subroutine progcld6 & + & ( plyr,plvl,tlyr,qlyr,qstl,rhly,clw, & ! --- inputs: + & xlat,xlon,slmsk,dz,delp, & + & ntrac,ntcw,ntiw,ntrw,ntsw,ntgl, & + & IX, NLAY, NLP1, & + & uni_cld, lmfshal, lmfdeep2, cldcov, & + & re_cloud,re_ice,re_snow, & + & clouds,clds,mtop,mbot,de_lgth & ! --- outputs: + & ) + +! ================= subprogram documentation block ================ ! +! ! +! subprogram: progcld6 computes cloud related quantities using ! +! the Thompson cloud microphysics scheme with updated microphysics ! +! cloud-radiation interaction (including subgrid clouds). ! +! ! +! abstract: this program computes cloud fractions from cloud ! +! condensates, ! +! and computes the low, mid, high, total and boundary layer cloud ! +! fractions and the vertical indices of low, mid, and high cloud ! +! top and base. the three vertical cloud domains are set up in the ! +! initial subroutine "cld_init". ! +! ! +! usage: call progcld6 ! +! ! +! subprograms called: gethml ! +! ! +! attributes: ! +! language: fortran 90 ! +! machine: ibm-sp, sgi ! +! ! +! ! +! ==================== definition of variables ==================== ! +! ! +! input variables: ! +! plyr (IX,NLAY) : model layer mean pressure in mb (100Pa) ! +! plvl (IX,NLP1) : model level pressure in mb (100Pa) ! +! tlyr (IX,NLAY) : model layer mean temperature in k ! +! tvly (IX,NLAY) : model layer virtual temperature in k ! +! qlyr (IX,NLAY) : layer specific humidity in gm/gm ! +! qstl (IX,NLAY) : layer saturate humidity in gm/gm ! +! rhly (IX,NLAY) : layer relative humidity (=qlyr/qstl) ! +! clw (IX,NLAY,ntrac) : layer cloud condensate amount ! +! xlat (IX) : grid latitude in radians, default to pi/2 -> -pi/2! +! range, otherwise see in-line comment ! +! xlon (IX) : grid longitude in radians (not used) ! +! slmsk (IX) : sea/land mask array (sea:0,land:1,sea-ice:2) ! +! dz (ix,nlay) : layer thickness (km) ! +! delp (ix,nlay) : model layer pressure thickness in mb (100Pa) ! +! IX : horizontal dimention ! +! NLAY,NLP1 : vertical layer/level dimensions ! +! uni_cld : logical - true for cloud fraction from shoc ! +! lmfshal : logical - true for mass flux shallow convection ! +! lmfdeep2 : logical - true for mass flux deep convection ! +! cldcov : layer cloud fraction (used when uni_cld=.true. ! +! ! +! output variables: ! +! clouds(IX,NLAY,NF_CLDS) : cloud profiles ! +! clouds(:,:,1) - layer total cloud fraction ! +! clouds(:,:,2) - layer cloud liq water path (g/m**2) ! +! clouds(:,:,3) - mean eff radius for liq cloud (micron) ! +! clouds(:,:,4) - layer cloud ice water path (g/m**2) ! +! clouds(:,:,5) - mean eff radius for ice cloud (micron) ! +! clouds(:,:,6) - layer rain drop water path not assigned ! +! clouds(:,:,7) - mean eff radius for rain drop (micron) ! +! *** clouds(:,:,8) - layer snow flake water path not assigned ! +! clouds(:,:,9) - mean eff radius for snow flake (micron) ! +! *** fu's scheme need to be normalized by snow density (g/m**3/1.0e6) ! +! clds (IX,5) : fraction of clouds for low, mid, hi, tot, bl ! +! mtop (IX,3) : vertical indices for low, mid, hi cloud tops ! +! mbot (IX,3) : vertical indices for low, mid, hi cloud bases ! +! de_lgth(ix) : clouds decorrelation length (km) ! +! ! +! module variables: ! +! ivflip : control flag of vertical index direction ! +! =0: index from toa to surface ! +! =1: index from surface to toa ! +! lmfshal : mass-flux shallow conv scheme flag ! +! lmfdeep2 : scale-aware mass-flux deep conv scheme flag ! +! lcrick : control flag for eliminating CRICK ! +! =t: apply layer smoothing to eliminate CRICK ! +! =f: do not apply layer smoothing ! +! lcnorm : control flag for in-cld condensate ! +! =t: normalize cloud condensate ! +! =f: not normalize cloud condensate ! +! ! +! ==================== end of description ===================== ! +! + implicit none + +! --- inputs + integer, intent(in) :: IX, NLAY, NLP1 + integer, intent(in) :: ntrac, ntcw, ntiw, ntrw, ntsw, ntgl + + logical, intent(in) :: uni_cld, lmfshal, lmfdeep2 + + real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, & + & tlyr, qlyr, qstl, rhly, cldcov, delp, dz, & + & re_cloud, re_ice, re_snow + + real (kind=kind_phys), dimension(:,:,:), intent(in) :: clw + + real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, & + & slmsk + +! --- outputs + real (kind=kind_phys), dimension(:,:,:), intent(out) :: clouds + + real (kind=kind_phys), dimension(:,:), intent(out) :: clds + real (kind=kind_phys), dimension(:), intent(out) :: de_lgth + + integer, dimension(:,:), intent(out) :: mtop,mbot + +! --- local variables: + real (kind=kind_phys), dimension(IX,NLAY) :: cldtot, cldcnv, & + & cwp, cip, crp, csp, rew, rei, res, rer, tem2d, clwf + + real (kind=kind_phys) :: ptop1(IX,NK_CLDS+1), rxlat(ix) + + real (kind=kind_phys) :: clwmin, clwm, clwt, onemrh, value, & + & tem1, tem2, tem3 + + integer :: i, k, id, nf + +! --- constant values + real (kind=kind_phys), parameter :: xrc3 = 100. + +! +!===> ... begin here +! + do k = 1, NLAY + do i = 1, IX + cldtot(i,k) = 0.0 + cldcnv(i,k) = 0.0 + cwp (i,k) = 0.0 + cip (i,k) = 0.0 + crp (i,k) = 0.0 + csp (i,k) = 0.0 + rew (i,k) = re_cloud(i,k) + rei (i,k) = re_ice(i,k) + rer (i,k) = rrain_def ! default rain radius to 1000 micron + res (i,k) = re_snow(i,K) +! tem2d (i,k) = min( 1.0, max( 0.0, (con_ttp-tlyr(i,k))*0.05 ) ) + clwf(i,k) = 0.0 + enddo + enddo + + do k = 1, NLAY + do i = 1, IX + 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. +!! ptopc(k,i): top presure of each cld domain (k=1-4 are sfc,L,m,h; +!! i=1,2 are low-lat (<45 degree) and pole regions) + + do i =1, IX + rxlat(i) = abs( xlat(i) / con_pi ) ! if xlat in pi/2 -> -pi/2 range +! rxlat(i) = abs(0.5 - xlat(i)/con_pi) ! if xlat in 0 -> pi range + enddo + + do id = 1, 4 + tem1 = ptopc(id,2) - ptopc(id,1) + do i =1, IX + ptop1(i,id) = ptopc(id,1) + tem1*max( 0.0, 4.0*rxlat(i)-1.0 ) + enddo + enddo + +!> - Compute cloud liquid/ice condensate path in \f$ g/m^2 \f$ . + + do k = 1, NLAY + do i = 1, IX + cwp(i,k) = max(0.0, clw(i,k,ntcw) * gfac * delp(i,k)) + cip(i,k) = max(0.0, clw(i,k,ntiw) * gfac * delp(i,k)) + crp(i,k) = max(0.0, clw(i,k,ntrw) * gfac * delp(i,k)) + csp(i,k) = max(0.0, (clw(i,k,ntsw)+clw(i,k,ntgl)) * + & gfac * delp(i,k)) + enddo + enddo + + if (uni_cld) then ! use unified sgs clouds generated outside + do k = 1, NLAY + do i = 1, IX + cldtot(i,k) = cldcov(i,k) + enddo + enddo + + else + +!> - Calculate layer cloud fraction. + + clwmin = 0.0 + + do k = 1, NLAY + do i = 1, IX + clwt = 1.0e-6 * (plyr(i,k)*0.001) +! clwt = 2.0e-6 * (plyr(i,k)*0.001) + + if (clwf(i,k) > clwt) then + onemrh= max( 1.e-10, 1.0-rhly(i,k) ) + clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) +! + tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) !jhan + !if (lmfdeep2) then + ! tem1 = xrc3 / tem1 + !else + tem1 = 100.0 / tem1 + !endif +! + value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) + tem2 = sqrt( sqrt(rhly(i,k)) ) + + cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 ) + endif + enddo + enddo + + endif ! if (uni_cld) then + + do k = 1, NLAY + do i = 1, IX + if (cldtot(i,k) < climit) then + cwp(i,k) = 0.0 + cip(i,k) = 0.0 + crp(i,k) = 0.0 + csp(i,k) = 0.0 + endif + enddo + enddo + + if ( lcnorm ) then + do k = 1, NLAY + do i = 1, IX + if (cldtot(i,k) >= climit) then + tem1 = 1.0 / max(climit2, cldtot(i,k)) + cwp(i,k) = cwp(i,k) * tem1 + cip(i,k) = cip(i,k) * tem1 + crp(i,k) = crp(i,k) * tem1 + csp(i,k) = csp(i,k) * tem1 + endif + enddo + enddo + endif + +! + do k = 1, NLAY + do i = 1, IX + clouds(i,k,1) = cldtot(i,k) + clouds(i,k,2) = cwp(i,k) + clouds(i,k,3) = rew(i,k) + clouds(i,k,4) = cip(i,k) + clouds(i,k,5) = rei(i,k) + clouds(i,k,6) = crp(i,k) ! added for Thompson + clouds(i,k,7) = rer(i,k) + clouds(i,k,8) = csp(i,k) ! added for Thompson + clouds(i,k,9) = res(i,k) + enddo + enddo + +! --- ... estimate clouds decorrelation length in km +! this is only a tentative test, need to consider change later + + if ( iovr == 3 ) then + do i = 1, ix + de_lgth(i) = max( 0.6, 2.78-4.6*rxlat(i) ) + enddo + endif + +!> - Call gethml() to compute low,mid,high,total, and boundary layer +!! cloud fractions and clouds top/bottom layer indices for low, mid, +!! and high clouds. +! --- compute low, mid, high, total, and boundary layer cloud fractions +! and clouds top/bottom layer indices for low, mid, and high clouds. +! The three cloud domain boundaries are defined by ptopc. The cloud +! overlapping method is defined by control flag 'iovr', which may +! be different for lw and sw radiation programs. + + call gethml & +! --- inputs: + & ( plyr, ptop1, cldtot, cldcnv, dz, de_lgth, & + & IX,NLAY, & +! --- outputs: + & clds, mtop, mbot & + & ) + + +! + return +!................................... + end subroutine progcld6 +!................................... + !> \ingroup module_radiation_clouds !> This subroutine computes cloud related quantities using !! for unified cloud microphysics scheme.