From 1337aeb68955a3e1346a80f3f2cb230f7eabd89a Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 23 Nov 2021 10:02:33 -0700 Subject: [PATCH 01/18] Fix dimensions of vertical eta level variables in several metadata files --- physics/GFS_rrtmg_setup.meta | 2 +- physics/GFS_rrtmgp_setup.meta | 2 +- physics/GFS_stochastics.meta | 2 +- physics/cires_ugwp.meta | 4 ++-- physics/ugwpv1_gsldrag.meta | 4 ++-- physics/unified_ugwp.meta | 4 ++-- 6 files changed, 9 insertions(+), 9 deletions(-) diff --git a/physics/GFS_rrtmg_setup.meta b/physics/GFS_rrtmg_setup.meta index ecd849c48..d80faf8a5 100644 --- a/physics/GFS_rrtmg_setup.meta +++ b/physics/GFS_rrtmg_setup.meta @@ -12,7 +12,7 @@ standard_name = sigma_pressure_hybrid_vertical_coordinate long_name = vertical sigma coordinate for radiation initialization units = none - dimensions = (vertical_interface_dimension_for_radiation) + dimensions = (vertical_interface_dimension) type = real kind = kind_phys intent = in diff --git a/physics/GFS_rrtmgp_setup.meta b/physics/GFS_rrtmgp_setup.meta index 4043392a9..ab9b0a49c 100644 --- a/physics/GFS_rrtmgp_setup.meta +++ b/physics/GFS_rrtmgp_setup.meta @@ -76,7 +76,7 @@ standard_name = sigma_pressure_hybrid_vertical_coordinate long_name = vertical sigma coordinate for radiation initialization units = none - dimensions = (vertical_interface_dimension_for_radiation) + dimensions = (vertical_interface_dimension) type = real kind = kind_phys intent = in diff --git a/physics/GFS_stochastics.meta b/physics/GFS_stochastics.meta index 0b2c1da2f..c78dbe015 100644 --- a/physics/GFS_stochastics.meta +++ b/physics/GFS_stochastics.meta @@ -31,7 +31,7 @@ standard_name = sigma_pressure_hybrid_vertical_coordinate long_name = vertical sigma coordinate for radiation initialization units = none - dimensions = (vertical_interface_dimension_for_radiation) + dimensions = (vertical_interface_dimension) type = real kind = kind_phys intent = in diff --git a/physics/cires_ugwp.meta b/physics/cires_ugwp.meta index fe0e82390..993516e4f 100644 --- a/physics/cires_ugwp.meta +++ b/physics/cires_ugwp.meta @@ -78,7 +78,7 @@ standard_name = sigma_pressure_hybrid_coordinate_a_coefficient long_name = a parameter for sigma pressure level calculations units = Pa - dimensions = (vertical_interface_dimension_for_radiation) + dimensions = (vertical_interface_dimension) type = real kind = kind_phys intent = in @@ -86,7 +86,7 @@ standard_name = sigma_pressure_hybrid_coordinate_b_coefficient long_name = b parameter for sigma pressure level calculations units = none - dimensions = (vertical_interface_dimension_for_radiation) + dimensions = (vertical_interface_dimension) type = real kind = kind_phys intent = in diff --git a/physics/ugwpv1_gsldrag.meta b/physics/ugwpv1_gsldrag.meta index 64d6b0d64..3ffcf909c 100644 --- a/physics/ugwpv1_gsldrag.meta +++ b/physics/ugwpv1_gsldrag.meta @@ -84,7 +84,7 @@ standard_name = sigma_pressure_hybrid_coordinate_a_coefficient long_name = a parameter for sigma pressure level calculations units = Pa - dimensions = (vertical_interface_dimension_for_radiation) + dimensions = (vertical_interface_dimension) type = real kind = kind_phys intent = in @@ -92,7 +92,7 @@ standard_name = sigma_pressure_hybrid_coordinate_b_coefficient long_name = b parameter for sigma pressure level calculations units = none - dimensions = (vertical_interface_dimension_for_radiation) + dimensions = (vertical_interface_dimension) type = real kind = kind_phys intent = in diff --git a/physics/unified_ugwp.meta b/physics/unified_ugwp.meta index 547256681..246ca236e 100644 --- a/physics/unified_ugwp.meta +++ b/physics/unified_ugwp.meta @@ -86,7 +86,7 @@ standard_name = sigma_pressure_hybrid_coordinate_a_coefficient long_name = a parameter for sigma pressure level calculations units = Pa - dimensions = (vertical_interface_dimension_for_radiation) + dimensions = (vertical_interface_dimension) type = real kind = kind_phys intent = in @@ -94,7 +94,7 @@ standard_name = sigma_pressure_hybrid_coordinate_b_coefficient long_name = b parameter for sigma pressure level calculations units = none - dimensions = (vertical_interface_dimension_for_radiation) + dimensions = (vertical_interface_dimension) type = real kind = kind_phys intent = in From 98aba4a6b62e3600de76a9c74c4a09e0ce783ebe Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 23 Nov 2021 10:03:07 -0700 Subject: [PATCH 02/18] Suggested improvements for Thompson effective radii calculation and consistency checking --- physics/GFS_rrtmg_pre.F90 | 44 ++++++++++---------------- physics/module_mp_thompson.F90 | 10 ------ physics/mp_thompson.F90 | 58 ++++++++++++++++++++-------------- 3 files changed, 50 insertions(+), 62 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index ca3bf0e70..92f9d7122 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -193,9 +193,9 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & effrl, effri, effrr, effrs, rho, orho, plyrpa ! for Thompson MP - real(kind=kind_phys), dimension(im,lm+LTP) :: & - re_cloud, re_ice, re_snow, qv_mp, qc_mp, & - qi_mp, qs_mp, nc_mp, ni_mp, nwfa + real(kind=kind_phys), dimension(im,lm+LTP) :: & + qv_mp, qc_mp, qi_mp, qs_mp, & + nc_mp, ni_mp, nwfa real (kind=kind_phys), dimension(lm) :: cldfra1d, qv1d, & & qc1d, qi1d, qs1d, dz1d, p1d, t1d @@ -796,37 +796,25 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & ! it will raise the low limit from 5 to 10, but the high limit will remain 125. call calc_effectRad (tlyr(i,:), plyr(i,:)*100., 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 ) + effrl(i,:), effri(i,:), effrs(i,:), 1, lm ) do k=1,lm - re_cloud(i,k) = MAX(re_qc_min, MIN(re_cloud(i,k), re_qc_max)) - re_ice(i,k) = MAX(re_qi_min, MIN(re_ice(i,k), re_qi_max)) - re_snow(i,k) = MAX(re_qs_min, MIN(re_snow(i,k), re_qs_max)) - end do - end do - ! Scale Thompson's effective radii from meter to micron - do k=1,lm - do i=1,im - re_cloud(i,k) = re_cloud(i,k)*1.e6 - re_ice(i,k) = re_ice(i,k)*1.e6 - re_snow(i,k) = re_snow(i,k)*1.e6 + effrl(i,k) = MAX(re_qc_min, MIN(effrl(i,k), re_qc_max)) + effri(i,k) = MAX(re_qi_min, MIN(effri(i,k), re_qi_max)) + effrr(i,k) = 1000. ! rrain_def=1000. + effrs(i,k) = MAX(re_qs_min, MIN(effrs(i,k), re_qs_max)) end do + effrl(i,lmk) = re_qc_min + effri(i,lmk) = re_qi_min + effrr(i,lmk) = 1000. ! rrain_def=1000. + effrs(i,lmk) = re_qs_min end do + ! Update global arrays, scale Thompson's effective radii from meter to micron do k=1,lm k1 = k + kd do i=1,im - effrl(i,k1) = re_cloud (i,k) - effri(i,k1) = re_ice (i,k) - effrr(i,k1) = 1000. ! rrain_def=1000. - effrs(i,k1) = re_snow(i,k) - enddo - enddo - ! Update global arrays - do k=1,lm - k1 = k + kd - do i=1,im - effrl_inout(i,k) = effrl(i,k1) - effri_inout(i,k) = effri(i,k1) - effrs_inout(i,k) = effrs(i,k1) + effrl_inout(i,k) = effrl(i,k1)*1.e6 + effri_inout(i,k) = effri(i,k1)*1.e6 + effrs_inout(i,k) = effrs(i,k1)*1.e6 enddo enddo else ! all other cases diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 353f83c78..3c7e6ec2b 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1252,16 +1252,6 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & ndt = max(nint(dt_in/dt_inner),1) dt = dt_in/ndt if(dt_in .le. dt_inner) dt= dt_in - if(nsteps>1 .and. ndt>1) then - if (present(errmsg) .and. present(errflg)) then - write(errmsg, '(a)') 'Logic error in mp_gt_driver: inner loop cannot be used with subcycling' - errflg = 1 - return - else - write(*,'(a)') 'Warning: inner loop cannot be used with subcycling, resetting ndt=1' - ndt = 1 - endif - endif do it = 1, ndt diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index d5a1fcaad..d5111cf63 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -379,6 +379,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ! Reduced time step if subcycling is used real(kind_phys) :: dtstep + integer :: ndt ! Air density real(kind_phys) :: rho(1:ncol,1:nlev) !< kg m-3 ! Water vapor mixing ratio (instead of specific humidity) @@ -458,11 +459,39 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & errmsg = '' errflg = 0 - ! Check initialization state - if (.not.is_initialized) then - write(errmsg, fmt='((a))') 'mp_thompson_run called before mp_thompson_init' - errflg = 1 - return + if (first_time_step .and. istep==1 .and. blkno==1) then + ! Check initialization state + if (.not.is_initialized) then + write(errmsg, fmt='((a))') 'mp_thompson_run called before mp_thompson_init' + errflg = 1 + return + end if + ! Check forr optional arguments of aerosol-aware microphysics + if (is_aerosol_aware .and. .not. (present(nc) .and. & + present(nwfa) .and. & + present(nifa) .and. & + present(nwfa2d) .and. & + present(nifa2d) )) then + write(errmsg,fmt='(*(a))') 'Logic error in mp_thompson_run:', & + ' aerosol-aware microphysics require all of the', & + ' following optional arguments:', & + ' nc, nwfa, nifa, nwfa2d, nifa2d' + errflg = 1 + return + end if + ! Consistency cheecks - subcycling and inner loop at the same time are not supported + if (nsteps>1 .and. dt_inner < dtp) then + write(errmsg,'(*(a))') "Logic error: Subcycling and inner loop cannot be used at the same time" + errflg = 1 + return + else if (mpirank==mpiroot .and. nsteps>1) then + write(*,'(a,i0,a,a,f6.2,a)') 'Thompson MP is using ', nsteps, ' substep(s) per time step with an ', & + 'effective time step of ', dtp/real(nsteps, kind=kind_phys), ' seconds' + else if (mpirank==mpiroot .and. dt_inner < dtp) then + ndt = max(nint(dtp/dt_inner),1) + write(*,'(a,i0,a,a,f6.2,a)') 'Thompson MP is using ', ndt, ' inner loops per time step with an ', & + 'effective time step of ', dtp/real(ndt, kind=kind_phys), ' seconds' + end if end if ! Set reduced time step if subcycling is used @@ -471,25 +500,6 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & else dtstep = dtp end if - if (first_time_step .and. istep==1 .and. mpirank==mpiroot .and. blkno==1) then - write(*,'(a,i0,a,a,f6.2,a)') 'Thompson MP is using ', nsteps, ' substep(s) per time step', & - ' with an effective time step of ', dtstep, ' seconds' - end if - - if (first_time_step .and. istep==1) then - if (is_aerosol_aware .and. .not. (present(nc) .and. & - present(nwfa) .and. & - present(nifa) .and. & - present(nwfa2d) .and. & - present(nifa2d) )) then - write(errmsg,fmt='(*(a))') 'Logic error in mp_thompson_run:', & - ' aerosol-aware microphysics require all of the', & - ' following optional arguments:', & - ' nc, nwfa, nifa, nwfa2d, nifa2d' - errflg = 1 - return - end if - end if !> - Convert specific humidity to water vapor mixing ratio. !> - Also, hydrometeor variables are mass or number mixing ratio From 12ccaddd72fe83a9148eae9e4bd2f1a6282dea72 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 6 Dec 2021 08:19:30 -0700 Subject: [PATCH 03/18] Fix compile error in physics/radiation_clouds.f --- physics/radiation_clouds.f | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 3f6c54d5d..d8c6f9f59 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -3139,10 +3139,10 @@ subroutine progcld6 & clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) if (.not. lmfshal) then - tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0) - tem1 = 2000.0 / tem1 + tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0) + tem1 = 2000.0 / tem1 else - tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) !jhan + tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) !jhan if (lmfdeep2) then tem1 = xrc3 / tem1 else From 7f444eea81e5f391ab8c119b5d70f7276e1ec0c9 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Thu, 9 Dec 2021 21:12:17 -0700 Subject: [PATCH 04/18] Revert Ruiyu's changes in this branch --- physics/GFS_rrtmg_pre.F90 | 10 +------ physics/radiation_clouds.f | 53 +++++++++++++++----------------------- 2 files changed, 22 insertions(+), 41 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index f2efd9d26..b8c2edd7f 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -864,14 +864,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & enddo endif - if(imp_physics == imp_physics_thompson) then - do k=1,lm - k1 = k + kd - do i=1,im - cnvw (i,k1) = cnvw_in(i,k) - enddo - enddo - endif if (imp_physics == imp_physics_zhao_carr) then ccnd(1:IM,1:LMK,1) = ccnd(1:IM,1:LMK,1) + cnvw(1:IM,1:LMK) @@ -1001,7 +993,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & ntrac-1, ntcw-1,ntiw-1,ntrw-1, & ntsw-1,ntgl-1, & im, lmk, lmp, uni_cld, lmfshal, lmfdeep2, & - cldcov(:,1:LMK), cnvw, effrl_inout, & + cldcov(:,1:LMK), effrl_inout, & effri_inout, effrs_inout, & lwp_ex, iwp_ex, lwp_fc, iwp_fc, & dzb, xlat_d, julian, yearlen, & diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index b62f9f9e6..f58ec8d11 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -2881,7 +2881,7 @@ subroutine progcld6 & & xlat,xlon,slmsk,dz,delp, & & ntrac,ntcw,ntiw,ntrw,ntsw,ntgl, & & IX, NLAY, NLP1, & - & uni_cld, lmfshal, lmfdeep2, cldcov, cnvw, & + & uni_cld, lmfshal, lmfdeep2, cldcov, & & re_cloud,re_ice,re_snow, & & lwp_ex, iwp_ex, lwp_fc, iwp_fc, & & dzlay, latdeg, julian, yearlen, & @@ -2976,7 +2976,7 @@ subroutine progcld6 & real (kind=kind_phys), dimension(:,:), intent(in) :: plvl, plyr, & & tlyr, qlyr, qstl, rhly, cldcov, delp, dz, dzlay, & - & re_cloud, re_ice, re_snow, cnvw + & re_cloud, re_ice, re_snow real (kind=kind_phys), dimension(:), intent(inout) :: & & lwp_ex, iwp_ex, lwp_fc, iwp_fc @@ -3010,8 +3010,8 @@ subroutine progcld6 & integer :: i, k, id, nf ! --- constant values - real (kind=kind_phys), parameter :: xrc3 = 200. -! real (kind=kind_phys), parameter :: xrc3 = 100. +! real (kind=kind_phys), parameter :: xrc3 = 200. + real (kind=kind_phys), parameter :: xrc3 = 100. ! !===> ... begin here @@ -3065,7 +3065,6 @@ subroutine progcld6 & 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) + cnvw(i,k) enddo enddo !> - Find top pressure for each cloud domain for given latitude. @@ -3092,9 +3091,8 @@ subroutine progcld6 & 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)) - csp(i,k) = max(0.0, clw(i,k,ntsw) * gfac * delp(i,k)) + csp(i,k) = max(0.0, (clw(i,k,ntsw)+clw(i,k,ntgl)) * & + & gfac * delp(i,k)) enddo enddo @@ -3127,36 +3125,27 @@ subroutine progcld6 & clwmin = 0.0 do k = 1, NLAY-1 do i = 1, IX -! clwt = 1.0e-6 * (plyr(i,k)*0.001) -! clwt = 2.0e-6 * (plyr(i,k)*0.001) - clwt = 1.0e-10 * (plyr(i,k)*0.001) + clwt = 1.0e-6 * (plyr(i,k)*0.001) if (clwf(i,k) > clwt) then - if(rhly(i,k) > 1.) then - cldtot(i,k) = 1. - else - onemrh= max( 1.e-10, 1.0-rhly(i,k) ) - clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) + onemrh= max( 1.e-10, 1.0-rhly(i,k) ) + clwm = clwmin / max( 0.01, plyr(i,k)*0.001 ) - if (.not. lmfshal) then - tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0) - tem1 = 2000.0 / tem1 + if (.not. lmfshal) then + tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0) + tem1 = 2000.0 / tem1 + else + tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) !jhan + if (lmfdeep2) then + tem1 = xrc3 / tem1 else - tem1 = min(max((onemrh*qstl(i,k))**0.49,0.0001),1.0) !jhan - if (lmfdeep2) then - tem1 = xrc3 / tem1 - else - tem1 = 100.0 / tem1 - endif -! - value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) - tem2 = sqrt( sqrt(rhly(i,k)) ) - - cldtot(i,k) = max( tem2*(1.0-exp(-value)), 0.0 ) + tem1 = 100.0 / tem1 endif endif - else - cldtot(i,k) = 0.0 + + 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 From 5ebe5dc62ae85946c9ff092203a304bd2740074d Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Mon, 13 Dec 2021 14:36:41 -0700 Subject: [PATCH 05/18] following early results by Anning, make fewer clouds, especially high clouds --- physics/GFS_rrtmg_pre.F90 | 4 +++- physics/radiation_clouds.f | 12 ++++++------ 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 3a3378e15..cb20e69fb 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -216,6 +216,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & real (kind=kind_phys) :: max_relh integer :: iflag + integer :: ii_half integer :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte @@ -236,7 +237,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & LP1 = LM + 1 ! num of in/out levels - gridkm = sqrt(2.0)*sqrt(dx(1)*0.001*dx(1)*0.001) + ii_half = nint(0.5*LM) + gridkm = sqrt(dx(1)*0.001*dx(ii_half)*0.001) if (imp_physics == imp_physics_thompson) then max_relh = 1.5 diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index f58ec8d11..84dfb2667 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -3452,7 +3452,7 @@ subroutine progcld_thompson & do i = 1, IX cwp(i,k) = max(0.0, clw(i,k,ntcw) * dz(i,k)*1.E6) crp(i,k) = 0.0 - snow_mass_factor = 0.85 + snow_mass_factor = 0.90 cip(i,k) = max(0.0, (clw(i,k,ntiw) & & + (1.0-snow_mass_factor)*clw(i,k,ntsw))*dz(i,k)*1.E6) if (re_snow(i,k) .gt. snow_max_radius)then @@ -4532,8 +4532,8 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & DO k = kts,kte delz = MAX(100., dz(k)) - RH_00L = 0.74+MIN(0.25,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) - RH_00O = 0.82+MIN(0.17,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) + RH_00L = 0.77+MIN(0.22,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) + RH_00O = 0.85+MIN(0.14,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) RHUM = rh(k) if (qc(k).ge.1.E-5 .or. qi(k).ge.1.E-5 & @@ -4550,7 +4550,7 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & RH_00 = RH_00L ENDIF - tc = t(k) - 273.15 + tc = MAX(-80.0, t(k) - 273.15) if (tc .lt. -12.0) RH_00 = RH_00L if (tc .gt. 20.0) then @@ -4562,12 +4562,12 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & if (max_relh.gt.1.12 .or. (.NOT.(modify_qvapor)) ) then !..For HRRR model, the following look OK. RHUM = MIN(rh(k), 1.45) - RH_00 = RH_00 + (1.45-RH_00)*(-12.0-tc)/(-12.0+112.) + RH_00 = RH_00 + (1.45-RH_00)*(-12.0-tc)/(-12.0+85.) CLDFRA(K) = MAX(0.,1.0-SQRT((1.46-RHUM)/(1.46-RH_00))) else !..but for the GFS model, RH is way lower. RHUM = MIN(rh(k), 1.05) - RH_00 = RH_00 + (1.05-RH_00)*(-12.0-tc)/(-12.0+112.) + RH_00 = RH_00 + (1.05-RH_00)*(-12.0-tc)/(-12.0+85.) CLDFRA(K) = MAX(0.,1.0-SQRT((1.06-RHUM)/(1.06-RH_00))) endif endif From b7ddc451611c47779e5092d181d3631affbaa356 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Tue, 14 Dec 2021 16:02:27 -0700 Subject: [PATCH 06/18] add in less LWC and IWC in the partly cloudy boxes --- physics/radiation_clouds.f | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 84dfb2667..823575ddd 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -4789,9 +4789,9 @@ SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs,T,dz,entr, k1,k2,kts,kte) max_iwc = ABS(qvs(k2)-qvs(k1)) do k = k1, k2 - max_iwc = MAX(1.E-5, max_iwc - (qi(k)+qs(k))) + max_iwc = MAX(1.E-6, max_iwc - (qi(k)+qs(k))) enddo - max_iwc = MIN(2.E-3, max_iwc) + max_iwc = MIN(1.E-4, max_iwc) this_dz = 0.0 do k = k1, k2 @@ -4801,7 +4801,7 @@ SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs,T,dz,entr, k1,k2,kts,kte) this_dz = this_dz + dz(k) endif this_iwc = max_iwc*this_dz/tdz - iwc = MAX(5.E-6, this_iwc*(1.-entr)) + iwc = MAX(1.E-6, this_iwc*(1.-entr)) if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.203.16) then qi(k) = qi(k) + cfr(k)*cfr(k)*iwc endif @@ -4830,9 +4830,9 @@ SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte) ! print*, ' max_lwc = ', max_lwc, ' over DZ=',tdz do k = k1, k2 - max_lwc = MAX(1.E-5, max_lwc - qc(k)) + max_lwc = MAX(1.E-6, max_lwc - qc(k)) enddo - max_lwc = MIN(2.E-3, max_lwc) + max_lwc = MIN(1.E-4, max_lwc) this_dz = 0.0 do k = k1, k2 @@ -4842,7 +4842,7 @@ SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte) this_dz = this_dz + dz(k) endif this_lwc = max_lwc*this_dz/tdz - lwc = MAX(5.E-6, this_lwc*(1.-entr)) + lwc = MAX(1.E-6, this_lwc*(1.-entr)) if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.253.16) then qc(k) = qc(k) + cfr(k)*cfr(k)*lwc endif From a1ca10e804d68eeb8c97ff2f8c57d52ad65868d8 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Thu, 16 Dec 2021 10:18:35 -0700 Subject: [PATCH 07/18] bug fix, LM should have been IM --- physics/GFS_rrtmg_pre.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index cb20e69fb..cc68825e1 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -237,7 +237,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & LP1 = LM + 1 ! num of in/out levels - ii_half = nint(0.5*LM) + ii_half = nint(0.5*IM) gridkm = sqrt(dx(1)*0.001*dx(ii_half)*0.001) if (imp_physics == imp_physics_thompson) then From 79325a53d7ee77c9c378d59a1cfd506c5b0f78ee Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Fri, 17 Dec 2021 08:37:33 -0700 Subject: [PATCH 08/18] make gridkm even simpler --- physics/GFS_rrtmg_pre.F90 | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index cc68825e1..f08ff2752 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -216,7 +216,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & real (kind=kind_phys) :: max_relh integer :: iflag - integer :: ii_half integer :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte @@ -237,8 +236,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & LP1 = LM + 1 ! num of in/out levels - ii_half = nint(0.5*IM) - gridkm = sqrt(dx(1)*0.001*dx(ii_half)*0.001) + gridkm = dx(nint(0.5*IM))*0.001 if (imp_physics == imp_physics_thompson) then max_relh = 1.5 From 20e3b7962495858506705e77294f1c51c38dad3b Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Fri, 17 Dec 2021 09:22:19 -0700 Subject: [PATCH 09/18] make gridkm array in X-dimension --- physics/GFS_rrtmg_pre.F90 | 7 +++---- physics/radiation_clouds.f | 8 ++++---- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index f08ff2752..f6f20487e 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -183,7 +183,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & integer :: i, j, k, k1, k2, lsk, lv, n, itop, ibtc, LP1, lla, llb, lya,lyb - real(kind=kind_phys) :: es, qs, delt, tem0d, gridkm, pfac + real(kind=kind_phys) :: es, qs, delt, tem0d, pfac + real(kind=kind_phys), dimension(im) :: gridkm real(kind=kind_phys), dimension(im) :: cvt1, cvb1, tem1d, tskn, xland @@ -235,9 +236,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & LP1 = LM + 1 ! num of in/out levels - - gridkm = dx(nint(0.5*IM))*0.001 - if (imp_physics == imp_physics_thompson) then max_relh = 1.5 else @@ -245,6 +243,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & endif do i = 1, IM + gridkm(i) = dx(i)*0.001 lwp_ex(i) = 0.0 iwp_ex(i) = 0.0 lwp_fc(i) = 0.0 diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 823575ddd..87a9620b2 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -3311,7 +3311,7 @@ subroutine progcld_thompson & ! slmsk (IX) : sea/land mask array (sea:0,land:1,sea-ice:2) ! ! dz (ix,nlay) : layer thickness (km) ! ! delp (ix,nlay) : model layer pressure thickness in mb (100Pa) ! -! gridkm : grid length in km ! +! gridkm (ix) : grid length in km ! ! IX : horizontal dimention ! ! NLAY,NLP1 : vertical layer/level dimensions ! ! uni_cld : logical - true for cloud fraction from shoc ! @@ -3370,8 +3370,8 @@ subroutine progcld_thompson & real (kind=kind_phys), dimension(:), intent(in) :: xlat, xlon, & & slmsk - real(kind=kind_phys), dimension(:), intent(in) :: latdeg - real(kind=kind_phys), intent(in) :: julian, gridkm + real(kind=kind_phys), dimension(:), intent(in) :: latdeg, gridkm + real(kind=kind_phys), intent(in) :: julian integer, intent(in) :: yearlen ! --- outputs @@ -3518,7 +3518,7 @@ subroutine progcld_thompson & endif call cal_cldfra3(cldfra1d, qv1d, qc1d, qi1d, qs1d, dz1d, & - & p1d, t1d, xland, gridkm, & + & p1d, t1d, xland, gridkm(i), & & .false., max_relh, 1, nlay, .false.) do k = 1, NLAY From c6bd0ad91723771f92147c791c5859ecda5b2e4c Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 20 Dec 2021 17:08:44 -0700 Subject: [PATCH 10/18] Remove physics/rte-rrtmgp/rrtmgp/kernels-openacc/mo_gas_optics_kernels.F90 from list of RRTMGP schemes --- CMakeLists.txt | 1 - 1 file changed, 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b8cb88418..b59a8ef33 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -80,7 +80,6 @@ set(SCHEMES_OPENMP_OFF ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rrtmgp/mo_ ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rrtmgp/mo_rrtmgp_constants.F90 ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rrtmgp/mo_rrtmgp_util_reorder.F90 ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rrtmgp/mo_gas_concentrations.F90 - ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rrtmgp/kernels-openacc/mo_gas_optics_kernels.F90 ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rrtmgp/mo_rrtmgp_util_string.F90 ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rrtmgp/kernels/mo_gas_optics_kernels.F90 ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rrtmgp/kernels/mo_rrtmgp_util_reorder_kernels.F90 From 9a60c107c292101df106fcef41ad763ab1665c37 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 20 Dec 2021 17:09:31 -0700 Subject: [PATCH 11/18] Write diag messages to stdout instead of stderr, use standard _OPENMP CPP directive --- physics/GFS_debug.F90 | 20 ++++++++++---------- physics/GFS_phys_time_vary.fv3.F90 | 8 ++++---- physics/GFS_phys_time_vary.scm.F90 | 8 ++++---- 3 files changed, 18 insertions(+), 18 deletions(-) diff --git a/physics/GFS_debug.F90 b/physics/GFS_debug.F90 index 23d1be573..9d5d24aa8 100644 --- a/physics/GFS_debug.F90 +++ b/physics/GFS_debug.F90 @@ -390,7 +390,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, #ifdef MPI use mpi #endif -#ifdef OPENMP +#ifdef _OPENMP use omp_lib #endif use GFS_typedefs, only: GFS_control_type, GFS_statein_type, & @@ -437,7 +437,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, mpisize = 1 mpicomm = 0 #endif -#ifdef OPENMP +#ifdef _OPENMP omprank = OMP_GET_THREAD_NUM() ompsize = nthreads #else @@ -445,7 +445,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, ompsize = 1 #endif -#ifdef OPENMP +#ifdef _OPENMP !$OMP BARRIER #endif #ifdef MPI @@ -929,7 +929,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Grid%jindx2_tau', Grid%jindx2_tau ) endif end if -#ifdef OPENMP +#ifdef _OPENMP !$OMP BARRIER #endif end do @@ -938,7 +938,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling, #endif end do -#ifdef OPENMP +#ifdef _OPENMP !$OMP BARRIER #endif #ifdef MPI @@ -1043,7 +1043,7 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup #ifdef MPI use mpi #endif -#ifdef OPENMP +#ifdef _OPENMP use omp_lib #endif use machine, only: kind_phys @@ -1092,7 +1092,7 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup mpisize = 1 mpicomm = 0 #endif -#ifdef OPENMP +#ifdef _OPENMP omprank = OMP_GET_THREAD_NUM() ompsize = nthreads #else @@ -1100,7 +1100,7 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup ompsize = 1 #endif -#ifdef OPENMP +#ifdef _OPENMP !$OMP BARRIER #endif #ifdef MPI @@ -1451,7 +1451,7 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%precip_overlap_param', Interstitial%precip_overlap_param ) end if end if -#ifdef OPENMP +#ifdef _OPENMP !$OMP BARRIER #endif end do @@ -1460,7 +1460,7 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup #endif end do -#ifdef OPENMP +#ifdef _OPENMP !$OMP BARRIER #endif #ifdef MPI diff --git a/physics/GFS_phys_time_vary.fv3.F90 b/physics/GFS_phys_time_vary.fv3.F90 index d6155e6b1..35fe08252 100644 --- a/physics/GFS_phys_time_vary.fv3.F90 +++ b/physics/GFS_phys_time_vary.fv3.F90 @@ -7,7 +7,7 @@ !> @{ module GFS_phys_time_vary -#ifdef OPENMP +#ifdef _OPENMP use omp_lib #endif @@ -355,7 +355,7 @@ subroutine GFS_phys_time_vary_init ( !$OMP section !--- if sncovr does not exist in the restart, need to create it if (all(sncovr < zero)) then - if (me == master ) write(0,'(a)') 'GFS_phys_time_vary_init: compute sncovr from weasd and soil vegetation parameters' + if (me == master ) write(*,'(a)') 'GFS_phys_time_vary_init: compute sncovr from weasd and soil vegetation parameters' !--- compute sncovr from existing variables !--- code taken directly from read_fix.f sncovr(:) = zero @@ -376,7 +376,7 @@ subroutine GFS_phys_time_vary_init ( !--- For RUC LSM: create sncovr_ice from sncovr if (lsm == lsm_ruc) then if (all(sncovr_ice < zero)) then - if (me == master ) write(0,'(a)') 'GFS_phys_time_vary_init: fill sncovr_ice with sncovr for RUC LSM' + if (me == master ) write(*,'(a)') 'GFS_phys_time_vary_init: fill sncovr_ice with sncovr for RUC LSM' sncovr_ice(:) = sncovr(:) endif endif @@ -396,7 +396,7 @@ subroutine GFS_phys_time_vary_init ( !--- land and ice - not for restart runs lsm_init: if (.not.flag_restart) then if (lsm == lsm_noahmp .or. lsm == lsm_ruc) then - if (me == master ) write(0,'(a)') 'GFS_phys_time_vary_init: initialize albedo for land and ice' + if (me == master ) write(*,'(a)') 'GFS_phys_time_vary_init: initialize albedo for land and ice' do ix=1,im albdvis_lnd(ix) = 0.2_kind_phys albdnir_lnd(ix) = 0.2_kind_phys diff --git a/physics/GFS_phys_time_vary.scm.F90 b/physics/GFS_phys_time_vary.scm.F90 index b06e46cdc..514988a48 100644 --- a/physics/GFS_phys_time_vary.scm.F90 +++ b/physics/GFS_phys_time_vary.scm.F90 @@ -6,7 +6,7 @@ !! aerosol, IN&CCN and surface properties updates. !> @{ module GFS_phys_time_vary - + use machine, only : kind_phys use mersenne_twister, only: random_setseed, random_number @@ -313,7 +313,7 @@ subroutine GFS_phys_time_vary_init ( !--- if sncovr does not exist in the restart, need to create it if (all(sncovr < zero)) then - if (me == master ) write(0,'(a)') 'GFS_phys_time_vary_init: compute sncovr from weasd and soil vegetation parameters' + if (me == master ) write(*,'(a)') 'GFS_phys_time_vary_init: compute sncovr from weasd and soil vegetation parameters' !--- compute sncovr from existing variables !--- code taken directly from read_fix.f sncovr(:) = zero @@ -334,7 +334,7 @@ subroutine GFS_phys_time_vary_init ( !--- For RUC LSM: create sncovr_ice from sncovr if (lsm == lsm_ruc) then if (all(sncovr_ice < zero)) then - if (me == master ) write(0,'(a)') 'GFS_phys_time_vary_init: fill sncovr_ice with sncovr for RUC LSM' + if (me == master ) write(*,'(a)') 'GFS_phys_time_vary_init: fill sncovr_ice with sncovr for RUC LSM' sncovr_ice(:) = sncovr(:) endif endif @@ -351,7 +351,7 @@ subroutine GFS_phys_time_vary_init ( !--- land and ice - not for restart runs lsm_init: if (.not.flag_restart) then if (lsm == lsm_noahmp .or. lsm == lsm_ruc) then - if (me == master ) write(0,'(a)') 'GFS_phys_time_vary_init: initialize albedo for land and ice' + if (me == master ) write(*,'(a)') 'GFS_phys_time_vary_init: initialize albedo for land and ice' do ix=1,im albdvis_lnd(ix) = 0.2_kind_phys albdnir_lnd(ix) = 0.2_kind_phys From ebd9495b819fe590628d635c4cb641190c712063 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 22 Dec 2021 06:38:53 -0700 Subject: [PATCH 12/18] Remove more duplicate modules from SCHEMES_OPENMP_OFF list in CMakeLists.txt --- CMakeLists.txt | 2 -- 1 file changed, 2 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index b59a8ef33..f16014cb7 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -107,8 +107,6 @@ set(SCHEMES_OPENMP_OFF ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rrtmgp/mo_ ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rte/mo_rte_sw.F90 ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rte/mo_fluxes.F90 ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rte/mo_rte_lw.F90 - ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rte/kernels-openacc/mo_rte_solver_kernels.F90 - ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rte/kernels-openacc/mo_optical_props_kernels.F90 ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rte/mo_rte_util_array.F90 ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rte/kernels/mo_rte_solver_kernels.F90 ${LOCAL_CURRENT_SOURCE_DIR}/physics/rte-rrtmgp/rte/kernels/mo_optical_props_kernels.F90 From c73a22a507c6d480e8ed2f34ab7436eaaccdcd51 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 22 Dec 2021 13:04:22 -0700 Subject: [PATCH 13/18] Move calculation fo effrr into its own loop --- physics/GFS_rrtmg_pre.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 615a83d0f..1252418c9 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -821,14 +821,16 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & do k=1,lm effrl(i,k) = MAX(re_qc_min, MIN(effrl(i,k), re_qc_max)) effri(i,k) = MAX(re_qi_min, MIN(effri(i,k), re_qi_max)) - effrr(i,k) = 1000. ! rrain_def=1000. effrs(i,k) = MAX(re_qs_min, MIN(effrs(i,k), re_qs_max)) end do effrl(i,lmk) = re_qc_min effri(i,lmk) = re_qi_min - effrr(i,lmk) = 1000. ! rrain_def=1000. effrs(i,lmk) = re_qs_min end do + do k=1,lm + k1 = k + kd + effrr(i,k1) = 1000. ! rrain_def=1000. + end do ! Update global arrays, scale Thompson's effective radii from meter to micron do k=1,lm k1 = k + kd From b6dcb52209056db5bd0b552f2ffa81c2cc23e725 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 22 Dec 2021 13:32:37 -0700 Subject: [PATCH 14/18] Add missing i-loop around effrr in physics/GFS_rrtmg_pre.F90 --- physics/GFS_rrtmg_pre.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 1252418c9..ced2d99a4 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -829,7 +829,9 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & end do do k=1,lm k1 = k + kd - effrr(i,k1) = 1000. ! rrain_def=1000. + do i=1,im + effrr(i,k1) = 1000. ! rrain_def=1000. + end do end do ! Update global arrays, scale Thompson's effective radii from meter to micron do k=1,lm From 5a933125338806e4542f51d9912f513ee5a620c8 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Thu, 6 Jan 2022 11:09:29 -0700 Subject: [PATCH 15/18] alter aerosol surface emission based on a WRF change tested by Jimy Dudhia --- physics/mp_thompson.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index e96f0e112..d60b9f77f 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -91,7 +91,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & real(kind_phys) :: orho(1:ncol,1:nlev) ! m3 kg-1 real(kind_phys) :: nc_local(1:ncol,1:nlev) ! needed because nc is only allocated if is_aerosol_aware is true ! - real (kind=kind_phys) :: h_01, airmass, niIN3, niCCN3 + real (kind=kind_phys) :: h_01, z1, niIN3, niCCN3 integer :: i, k ! Initialize the CCPP error handling variables @@ -192,8 +192,8 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & 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*5.E-11) + z1 = hgt(i,2)-hgt(i,1) + nwfa2d(i) = nwfa(i,1) * 0.000196 * (50./z1) do k = 2, nlev nwfa(i,k) = naCCN1+naCCN0*exp(-((hgt(i,k)-hgt(i,1))/1000.)*niCCN3) enddo @@ -212,8 +212,8 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & !+---+-----------------------------------------------------------------+ if (mpirank==mpiroot) write(*,*) ' Apparently there are no initial CCN aerosol 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*5.E-11) + z1 = hgt(i,2)-hgt(i,1) + nwfa2d(i) = nwfa(i,1) * 0.000196 * (50./z1) enddo else if (mpirank==mpiroot) write(*,*) ' Apparently initial CCN aerosol surface emission rates are present.' From a45d641173850a9aa7ec398a8176254d8f587bca Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Thu, 6 Jan 2022 11:19:59 -0700 Subject: [PATCH 16/18] one more tuning for reducing cloud ice amounts in partly cloudy boxes --- physics/radiation_clouds.f | 40 ++++++++++++++++---------------------- 1 file changed, 17 insertions(+), 23 deletions(-) diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index 87a9620b2..3122a0c43 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -3164,7 +3164,8 @@ subroutine progcld6 & enddo enddo - ! What portion of water and ice contents is associated with the partly cloudy boxes + ! What portion of water and ice contents is associated with the + ! partly cloudy boxes do i = 1, IX do k = 1, NLAY-1 if (cldtot(i,k).ge.climit .and. cldtot(i,k).lt.ovcst) then @@ -3311,7 +3312,7 @@ subroutine progcld_thompson & ! slmsk (IX) : sea/land mask array (sea:0,land:1,sea-ice:2) ! ! dz (ix,nlay) : layer thickness (km) ! ! delp (ix,nlay) : model layer pressure thickness in mb (100Pa) ! -! gridkm (ix) : grid length in km ! +! gridkm (IX) : grid length in km ! ! IX : horizontal dimention ! ! NLAY,NLP1 : vertical layer/level dimensions ! ! uni_cld : logical - true for cloud fraction from shoc ! @@ -3445,14 +3446,14 @@ subroutine progcld_thompson & enddo !> - Compute cloud liquid/ice condensate path in \f$ g/m^2 \f$ . -!> - Since using Thompson MP, assume 20 percent of snow is actually in +!> - Since using Thompson MP, assume 1 percent of snow is actually in !! ice sizes. do k = 1, NLAY-1 do i = 1, IX cwp(i,k) = max(0.0, clw(i,k,ntcw) * dz(i,k)*1.E6) crp(i,k) = 0.0 - snow_mass_factor = 0.90 + snow_mass_factor = 0.99 cip(i,k) = max(0.0, (clw(i,k,ntiw) & & + (1.0-snow_mass_factor)*clw(i,k,ntsw))*dz(i,k)*1.E6) if (re_snow(i,k) .gt. snow_max_radius)then @@ -4536,12 +4537,12 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & RH_00O = 0.85+MIN(0.14,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) RHUM = rh(k) - if (qc(k).ge.1.E-5 .or. qi(k).ge.1.E-5 & - & .or. (qs(k).gt.1.E-5 .and. t(k).lt.273.)) then + if (qc(k).ge.1.E-6 .or. qi(k).ge.1.E-7 & + & .or. (qs(k).gt.1.E-6 .and. t(k).lt.273.)) then CLDFRA(K) = 1.0 elseif (((qc(k)+qi(k)).gt.1.E-10) .and. & - & ((qc(k)+qi(k)).lt.1.E-5)) then - CLDFRA(K) = MIN(0.99, 0.20*(10.0 + log10(qc(k)+qi(k)))) + & ((qc(k)+qi(k)).lt.1.E-6)) then + CLDFRA(K) = MIN(0.99, 0.1*(11.0 + log10(qc(k)+qi(k)))) else IF ((XLAND-1.5).GT.0.) THEN !--- Ocean @@ -4585,15 +4586,6 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & call adjust_cloudFinal(cldfra, qc, qi, rhoa, dz, kts,kte) - if (debug_flag .and. ndebug.lt.25) then - do k = kts,kte - write(6,'(a,i3,f9.2,f7.1,f7.2,f6.1,f6.3,f12.7,f12.7,f12.7)') & - & ' DEBUG-GT: ', k, p(k)*0.01, dz(k), t(k)-273.15, & - & rh(k)*100., cldfra(k), qc(k)*1.E3, qi(k)*1.E3, qs(k)*1.E3 - enddo - ndebug = ndebug + 1 - endif - !..Intended for cold start model runs, we use modify_qvapor to ensure that cloudy !.. areas are actually saturated such that the inserted clouds do not evaporate a !.. timestep later. @@ -4735,9 +4727,9 @@ SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, Dz1d, entrmnt,& k = k - 1 ENDDO - k_cldb = k_m12C + 5 + k_cldb = k_m12C + 3 in_cloud = .false. - k = k_m12C + 4 + k = k_m12C + 2 DO WHILE (.not. in_cloud .AND. k.gt.kbot) k_cldt = 0 if (cfr1d(k).ge.0.01) then @@ -4786,7 +4778,8 @@ SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs,T,dz,entr, k1,k2,kts,kte) do k = k1, k2 tdz = tdz + dz(k) enddo - max_iwc = ABS(qvs(k2)-qvs(k1)) +! max_iwc = ABS(qvs(k2)-qvs(k1)) + max_iwc = MAX(0.0, qvs(k1)-qvs(k2)) do k = k1, k2 max_iwc = MAX(1.E-6, max_iwc - (qi(k)+qs(k))) @@ -4826,7 +4819,8 @@ SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte) do k = k1, k2 tdz = tdz + dz(k) enddo - max_lwc = ABS(qvs(k2)-qvs(k1)) +! max_lwc = ABS(qvs(k2)-qvs(k1)) + max_lwc = MAX(0.0, qvs(k1)-qvs(k2)) ! print*, ' max_lwc = ', max_lwc, ' over DZ=',tdz do k = k1, k2 @@ -4843,7 +4837,7 @@ SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte) endif this_lwc = max_lwc*this_dz/tdz lwc = MAX(1.E-6, this_lwc*(1.-entr)) - if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.253.16) then + if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.258.16) then qc(k) = qc(k) + cfr(k)*cfr(k)*lwc endif enddo @@ -4895,6 +4889,6 @@ SUBROUTINE adjust_cloudFinal(cfr, qc, qi, Rho,dz, kts,kte) END SUBROUTINE adjust_cloudFinal !........................................! - end module module_radiation_clouds ! + end module module_radiation_clouds !! @} !========================================! From 2181d0c6e19df42ed0a50f195068519fa06b2a43 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Thu, 6 Jan 2022 11:22:31 -0700 Subject: [PATCH 17/18] reduce max ice number conc to fewer than 500 per liter of air --- physics/module_mp_thompson.F90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 3183ca4bf..c5bc99e17 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -2188,7 +2188,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & ni(k) = MAX(R2, ni1d(k)*rho(k)) if (ni(k).le. R2) then lami = cie(2)/5.E-6 - ni(k) = MIN(999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) + ni(k) = MIN(499.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) endif L_qi(k) = .true. lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi @@ -2196,7 +2196,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & xDi = (bm_i + mu_i + 1.) * ilami if (xDi.lt. 5.E-6) then lami = cie(2)/5.E-6 - ni(k) = MIN(999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) + ni(k) = MIN(499.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) elseif (xDi.gt. 300.E-6) then lami = cie(2)/300.E-6 ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i @@ -2901,7 +2901,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !> - Freezing of aqueous aerosols based on Koop et al (2001, Nature) xni = smo0(k)+ni(k) + (pni_rfz(k)+pni_wfz(k)+pni_inu(k))*dtsave - if (is_aerosol_aware .AND. homogIce .AND. (xni.le.999.E3) & + if (is_aerosol_aware .AND. homogIce .AND. (xni.le.499.E3) & & .AND.(temp(k).lt.238).AND.(ssati(k).ge.0.4) ) then xnc = iceKoop(temp(k),qv(k),qvs(k),nwfa(k), dtsave) pni_iha(k) = xnc*odts @@ -3237,7 +3237,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & xDi = (bm_i + mu_i + 1.) * ilami if (xDi.lt. 5.E-6) then lami = cie(2)/5.E-6 - xni = MIN(999.D3, cig(1)*oig2*xri/am_i*lami**bm_i) + xni = MIN(499.D3, cig(1)*oig2*xri/am_i*lami**bm_i) niten(k) = (xni-ni1d(k)*rho(k))*odts*orho elseif (xDi.gt. 300.E-6) then lami = cie(2)/300.E-6 @@ -3248,8 +3248,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & niten(k) = -ni1d(k)*odts endif xni=MAX(0.,(ni1d(k) + niten(k)*dtsave)*rho(k)) - if (xni.gt.999.E3) & - niten(k) = (999.E3-ni1d(k)*rho(k))*odts*orho + if (xni.gt.499.E3) & + niten(k) = (499.E3-ni1d(k)*rho(k))*odts*orho !> - Rain tendency qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) & @@ -4187,7 +4187,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & lami = cie(2)/300.E-6 endif ni1d(k) = MIN(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, & - 999.D3/rho(k)) + 499.D3/rho(k)) endif qr1d(k) = qr1d(k) + qrten(k)*DT nr1d(k) = MAX(R2/rho(k), nr1d(k) + nrten(k)*DT) From 1a72e9891b8580537fa13e400c7850248b1e8bf4 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Sun, 9 Jan 2022 09:08:15 -0700 Subject: [PATCH 18/18] Bugfix for cloud effective radii computation: scale local arrays from m to micron --- physics/GFS_rrtmg_pre.F90 | 37 +++++++++++++++---------------------- 1 file changed, 15 insertions(+), 22 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 80caf766e..0e398b1b9 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -818,28 +818,24 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & call calc_effectRad (tlyr(i,:), plyr(i,:)*100., qv_mp(i,:), qc_mp(i,:), & nc_mp(i,:), qi_mp(i,:), ni_mp(i,:), qs_mp(i,:), & effrl(i,:), effri(i,:), effrs(i,:), 1, lm ) + ! Scale Thompson's effective radii from meter to micron do k=1,lm - effrl(i,k) = MAX(re_qc_min, MIN(effrl(i,k), re_qc_max)) - effri(i,k) = MAX(re_qi_min, MIN(effri(i,k), re_qi_max)) - effrs(i,k) = MAX(re_qs_min, MIN(effrs(i,k), re_qs_max)) + effrl(i,k) = MAX(re_qc_min, MIN(effrl(i,k), re_qc_max))*1.e6 + effri(i,k) = MAX(re_qi_min, MIN(effri(i,k), re_qi_max))*1.e6 + effrs(i,k) = MAX(re_qs_min, MIN(effrs(i,k), re_qs_max))*1.e6 end do - effrl(i,lmk) = re_qc_min - effri(i,lmk) = re_qi_min - effrs(i,lmk) = re_qs_min + effrl(i,lmk) = re_qc_min*1.e6 + effri(i,lmk) = re_qi_min*1.e6 + effrs(i,lmk) = re_qs_min*1.e6 end do + effrr(:,:) = 1000. ! rrain_def=1000. + ! Update global arrays do k=1,lm k1 = k + kd do i=1,im - effrr(i,k1) = 1000. ! rrain_def=1000. - end do - end do - ! Update global arrays, scale Thompson's effective radii from meter to micron - do k=1,lm - k1 = k + kd - do i=1,im - effrl_inout(i,k) = effrl(i,k1)*1.e6 - effri_inout(i,k) = effri(i,k1)*1.e6 - effrs_inout(i,k) = effrs(i,k1)*1.e6 + effrl_inout(i,k) = effrl(i,k1) + effri_inout(i,k) = effri(i,k1) + effrs_inout(i,k) = effrs(i,k1) enddo enddo else ! all other cases @@ -964,8 +960,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & ntrac-1, ntcw-1,ntiw-1,ntrw-1, & ntsw-1,ntgl-1, & im, lm, lmp, uni_cld, lmfshal, lmfdeep2, & - cldcov(:,1:LM), effrl_inout, & - effri_inout, effrs_inout, & + cldcov(:,1:LM), effrl, effri, effrs, & lwp_ex, iwp_ex, lwp_fc, iwp_fc, & dzb, xlat_d, julian, yearlen, gridkm, & clouds, cldsa, mtopa ,mbota, de_lgth, alpha) ! --- outputs @@ -998,8 +993,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & ntrac-1, ntcw-1,ntiw-1,ntrw-1, & ntsw-1,ntgl-1, & im, lm, lmp, uni_cld, lmfshal, lmfdeep2, & - cldcov(:,1:LM), effrl_inout, & - effri_inout, effrs_inout, & + cldcov(:,1:LM), effrl, effri, effrs, & lwp_ex, iwp_ex, lwp_fc, iwp_fc, & dzb, xlat_d, julian, yearlen, gridkm, & clouds, cldsa, mtopa ,mbota, de_lgth, alpha) ! --- outputs @@ -1010,8 +1004,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & ntrac-1, ntcw-1,ntiw-1,ntrw-1, & ntsw-1,ntgl-1, & im, lmk, lmp, uni_cld, lmfshal, lmfdeep2, & - cldcov(:,1:LMK), cnvw, effrl_inout, & - effri_inout, effrs_inout, & + cldcov(:,1:LMK), cnvw, effrl, effri, effrs,& lwp_ex, iwp_ex, lwp_fc, iwp_fc, & dzb, xlat_d, julian, yearlen, & clouds, cldsa, mtopa ,mbota, de_lgth, alpha) ! --- outputs