Merge branch 'ufs/dev' of https://github.com/ufs-community/ccpp-physics…

… into HEAD

physics/GFS_MP_generic_post.F90

@@ -20,7 +20,7 @@ module GFS_MP_generic_post
 !> @{
       subroutine GFS_MP_generic_post_run(                                                                                 &
         im, levs, kdt, nrcm, nncl, ntcw, ntrac, imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_nssl,    &
-        imp_physics_mg, imp_physics_fer_hires, cal_pre, cplflx, cplchm, progsigma, con_g, rainmin, dtf, frain, rainc,     &
+        imp_physics_mg, imp_physics_fer_hires, cal_pre, cplflx, cplchm, cpllnd, progsigma, con_g, rainmin, dtf, frain, rainc, &
         rain1, rann, xlat, xlon, gt0, gq0, prsl, prsi, phii, tsfc, ice, snow, graupel, save_t, save_q, rain0, ice0, snow0,&
         graupel0, del, rain, domr_diag, domzr_diag, domip_diag, doms_diag, tprcp, srflag, sr, cnvprcp, totprcp, totice,   &
         totsnw, totgrp, cnvprcpb, totprcpb, toticeb, totsnwb, totgrpb, rain_cpl, rainc_cpl, snow_cpl, pwat,               &
@@ -37,7 +37,7 @@ subroutine GFS_MP_generic_post_run(
       integer, intent(in) :: im, levs, kdt, nrcm, nncl, ntcw, ntrac, num_dfi_radar, index_of_process_dfi_radar
       integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_mg, imp_physics_fer_hires
       integer, intent(in) :: imp_physics_nssl
-      logical, intent(in) :: cal_pre, lssav, ldiag3d, qdiag3d, cplflx, cplchm, progsigma
+      logical, intent(in) :: cal_pre, lssav, ldiag3d, qdiag3d, cplflx, cplchm, cpllnd, progsigma
       integer, intent(in) :: index_of_temperature,index_of_process_mp
 
       integer                                                :: dfi_radar_max_intervals
@@ -249,7 +249,7 @@ subroutine GFS_MP_generic_post_run(
       ! Conversion factor from mm per day to m per physics timestep 
       tem = dtp * con_p001 / con_day
 
-!> - For GFDL and Thompson MP scheme, determine convective snow by surface temperature;
+!> - For GFDL, Thompson and NSSL MP schemes, determine convective snow by surface temperature;
 !! and determine explicit rain/snow by snow/ice/graupel coming out directly from MP
 !! and convective rainfall from the cumulus scheme if the surface temperature is below
 !! \f$0^oC\f$.
@@ -363,7 +363,7 @@ subroutine GFS_MP_generic_post_run(
          enddo
       endif
 
-      if (cplflx .or. cplchm) then
+      if (cplflx .or. cplchm .or. cpllnd) then
         do i = 1, im
           dsnow_cpl(i)= max(zero, rain(i) * srflag(i))
           drain_cpl(i)= max(zero, rain(i) - dsnow_cpl(i))
@@ -372,7 +372,7 @@ subroutine GFS_MP_generic_post_run(
         enddo
       endif
 
-      if (cplchm) then
+      if (cplchm .or. cpllnd) then
         do i = 1, im
           rainc_cpl(i) = rainc_cpl(i) + rainc(i)
         enddo

physics/GFS_MP_generic_post.meta

@@ -128,6 +128,13 @@
   dimensions = ()
   type = logical
   intent = in
+[cpllnd]
+  standard_name = flag_for_land_coupling
+  long_name = flag controlling cpllnd collection (default off)
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
 [progsigma]
   standard_name = do_prognostic_updraft_area_fraction
   long_name = flag for prognostic sigma in cumulus scheme

physics/GFS_radiation_surface.F90

@@ -50,7 +50,7 @@ end subroutine GFS_radiation_surface_init
 !! \htmlinclude GFS_radiation_surface_run.html
 !!
       subroutine GFS_radiation_surface_run (                            &
-        im, frac_grid, lslwr, lsswr, lsm, lsm_noahmp, lsm_ruc,          &
+        im, nf_albd, frac_grid, lslwr, lsswr, lsm, lsm_noahmp, lsm_ruc, &
         xlat, xlon, slmsk, lndp_type, n_var_lndp, sfc_alb_pert,         &
         lndp_var_list, lndp_prt_list, landfrac, snodl, snodi, sncovr,   &
         sncovr_ice, fice, zorl, hprime, tsfg, tsfa, tisfc, coszen,      &
@@ -62,12 +62,12 @@ subroutine GFS_radiation_surface_run (                            &
         semisbase, semis, sfcalb, sfc_alb_dif, errmsg, errflg)
 
       use module_radiation_surface,  only: f_zero, f_one,  &
-                                           epsln, NF_ALBD, &
+                                           epsln,          &
                                            setemis, setalb
 
       implicit none
 
-      integer,               intent(in) :: im
+      integer,               intent(in) :: im, nf_albd
       logical,               intent(in) :: frac_grid, lslwr, lsswr, use_cice_alb, cplice
       integer,               intent(in) :: lsm, lsm_noahmp, lsm_ruc, lndp_type, n_var_lndp
       real(kind=kind_phys),  intent(in) :: min_seaice, min_lakeice
@@ -184,7 +184,7 @@ subroutine GFS_radiation_surface_run (                            &
                      alvsf, alnsf, alvwf, alnwf, facsf, facwf, fice, tisfc,                    &
                      albdvis_lnd, albdnir_lnd, albivis_lnd, albinir_lnd,                       &
                      albdvis_ice, albdnir_ice, albivis_ice, albinir_ice,                       &
-                     IM, sfc_alb_pert, lndp_alb, fracl, fraco, fraci, icy,                     & !  ---  inputs
+                     im, nf_albd, sfc_alb_pert, lndp_alb, fracl, fraco, fraci, icy,            & !  ---  inputs
                      sfcalb )                                                                    !  ---  outputs
 
 !> -# Approximate mean surface albedo from vis- and nir- diffuse values.

physics/GFS_radiation_surface.meta

@@ -55,6 +55,13 @@
   dimensions = ()
   type = integer
   intent = in
+[nf_albd]
+  standard_name = number_of_components_for_surface_albedo
+  long_name = number of IR/VIS/UV compinents for surface albedo
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
 [frac_grid]
   standard_name = flag_for_fractional_landmask
   long_name = flag for fractional grid

physics/GFS_rrtmg_pre.F90

@@ -73,7 +73,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp,       &
       use surface_perturbation,      only: cdfnor,ppfbet
 
       ! For Thompson MP
-      use module_mp_thompson,        only: calc_effectRad, Nt_c,     &
+      use module_mp_thompson,        only: calc_effectRad,           &
+                                           Nt_c_l, Nt_c_o,           &
                                            re_qc_min, re_qc_max,     &
                                            re_qi_min, re_qi_max,     &
                                            re_qs_min, re_qs_max
@@ -245,6 +246,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp,       &
       real (kind=kind_phys) :: alpha0,beta0,m,s,cldtmp,tmp_wt,cdfz
       real (kind=kind_phys) :: max_relh
       integer  :: iflag
+      integer  :: islmsk
 
       integer :: ids, ide, jds, jde, kds, kde, &
                  ims, ime, jms, jme, kms, kme, &
@@ -609,11 +611,10 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp,       &
 
       endif                              ! end_if_ivflip
 
-!> - Call module_radiation_aerosols::setaer(),to setup aerosols
-!! property profile for radiation.
 
 !check  print *,' in grrad : calling setaer '
 
+!> - Initialize mass mixing ratio of aerosols from NASA GOCART or NASA MERRA-2
        if (ntchm>0 .and. iaermdl==2) then
           do k=1,levs
             do i=1,im
@@ -637,6 +638,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp,       &
         endif
 
 
+!> - Call module_radiation_aerosols::setaer() to setup aerosols
+!! property profile for radiation.
       call setaer (plvl, plyr, prslk1, tvly, rhly, slmsk,    & !  ---  inputs
                    tracer1, aer_nm, xlon, xlat, IM, LMK, LMP,&
                    lsswr,lslwr,                              &
@@ -654,7 +657,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp,       &
         enddo
        enddo
 
-      !> Aerosol direct feedback effect by smoke and dust
+      !> - Add aerosol direct feedback effect by smoke and dust
       if(aero_dir_fdb) then ! add smoke/dust extinctions
         do k = 1, LMK
           do i = 1, IM
@@ -747,7 +750,11 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp,       &
                 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,k)
+                if(nint(slmsk(i)) == 1) then
+                  nc_mp (i,k) = Nt_c_l*orho(i,k)
+                else
+                  nc_mp (i,k) = Nt_c_o*orho(i,k)
+                endif
                 ni_mp (i,k) = tracer1(i,k,ntinc)/(1.-qvs)
               enddo
             enddo
@@ -876,13 +883,14 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, lextop, ltp,       &
           end do
           !> - Call Thompson's subroutine calc_effectRad() to compute effective radii
           do i=1,im
+            islmsk = nint(slmsk(i))
             ! Effective radii [m] are now intent(out), bounds applied in calc_effectRad
             !tgs: progclduni has different limits for ice radii (10.0-150.0) than
             !     calc_effectRad (4.99-125.0 for WRFv3.8.1; 2.49-125.0 for WRFv4+)
             !     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,:), &
-                                 effrl(i,:), effri(i,:), effrs(i,:), 1, lm )
+                                 effrl(i,:), effri(i,:), effrs(i,:), islmsk, 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))*1.e6

physics/GFS_rrtmg_setup.F90

@@ -325,7 +325,7 @@ end subroutine GFS_rrtmg_setup_finalize
 
 ! Private functions
 
-
+!>Initialization of radiation calculations.
    subroutine radinit( si, NLAY, imp_physics, me, ltp, lextop )
 !...................................
 

physics/GFS_rrtmgp_cloud_mp.F90

@@ -12,8 +12,8 @@ module GFS_rrtmgp_cloud_mp
   use rrtmgp_lw_cloud_optics, only: &
        radliq_lwr => radliq_lwrLW, radliq_upr => radliq_uprLW,&
        radice_lwr => radice_lwrLW, radice_upr => radice_uprLW  
-  use module_mp_thompson, only: calc_effectRad, Nt_c, re_qc_min, re_qc_max, re_qi_min, &
-       re_qi_max, re_qs_min, re_qs_max
+  use module_mp_thompson, only: calc_effectRad, Nt_c_l, Nt_c_o, re_qc_min, re_qc_max,  &
+       re_qi_min, re_qi_max, re_qs_min, re_qs_max
   use module_mp_thompson_make_number_concentrations, only: make_IceNumber,             &
        make_DropletNumber, make_RainNumber
 
@@ -254,7 +254,7 @@ subroutine GFS_rrtmgp_cloud_mp_run(nCol, nLev, nTracers, ncnd, i_cldliq, i_cldic
        ! Update particle size using modified mixing-ratios from Thompson.
        call cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice_nc,   &
             i_cldliq_nc, i_twa, q_lay, p_lay, t_lay, tracer, con_eps, con_rd, ltaerosol,&
-            mraerosol, effrin_cldliq, effrin_cldice, effrin_cldsnow)
+            mraerosol, lsmask,  effrin_cldliq, effrin_cldice, effrin_cldsnow)
        cld_reliq  = effrin_cldliq
        cld_reice  = effrin_cldice
        cld_resnow = effrin_cldsnow
@@ -820,7 +820,7 @@ function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha)
 !! \section cmp_reff_Thompson_gen General Algorithm
   subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice_nc,   &
        i_cldliq_nc, i_twa, q_lay, p_lay, t_lay, tracer, con_eps, con_rd, ltaerosol,      &
-       mraerosol, effrin_cldliq, effrin_cldice, effrin_cldsnow)
+       mraerosol, lsmask, effrin_cldliq, effrin_cldice, effrin_cldsnow)
     implicit none
 
     ! Inputs
@@ -830,6 +830,7 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice
     real(kind_phys), intent(in) :: con_eps,con_rd
     real(kind_phys), dimension(:,:),intent(in) :: q_lay, p_lay, t_lay
     real(kind_phys), dimension(:,:,:),intent(in) :: tracer
+    real(kind_phys), dimension(:), intent(in) :: lsmask
 
     ! Outputs
     real(kind_phys), dimension(:,:), intent(inout) :: effrin_cldliq, effrin_cldice,      &
@@ -840,6 +841,7 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice
     real(kind_phys) :: rho, orho
     real(kind_phys),dimension(nCol,nLev) :: qv_mp, qc_mp, qi_mp, qs_mp, ni_mp, nc_mp,    &
          nwfa, re_cloud, re_ice, re_snow
+    integer :: ilsmask 
 
     ! Prepare cloud mixing-ratios and number concentrations for calc_effectRa
     do iLay = 1, nLev
@@ -863,7 +865,11 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice
                nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho, nwfa(iCol,iLay)*rho) * orho
              endif
           else
-             nc_mp(iCol,iLay) = nt_c*orho
+             if (nint(lsmask(iCol)) == 1) then !land
+                nc_mp(iCol,iLay) = nt_c_l*orho
+             else 
+                nc_mp(iCol,iLay) = nt_c_o*orho
+             endif 
           endif
           if (qi_mp(iCol,iLay) > 1.e-12 .and. ni_mp(iCol,iLay) < 100.) then
              ni_mp(iCol,iLay) = make_IceNumber(qi_mp(iCol,iLay)*rho, t_lay(iCol,iLay)) * orho
@@ -873,9 +879,11 @@ subroutine cmp_reff_Thompson(nLev, nCol, i_cldliq, i_cldice, i_cldsnow, i_cldice
 
     ! Compute effective radii for liquid/ice/snow.
     do iCol=1,nCol
+       ilsmask = nint(lsmask(iCol))
        call calc_effectRad (t_lay(iCol,:), p_lay(iCol,:), qv_mp(iCol,:), qc_mp(iCol,:),  &
                             nc_mp(iCol,:), qi_mp(iCol,:), ni_mp(iCol,:), qs_mp(iCol,:),  &
-                            re_cloud(iCol,:), re_ice(iCol,:), re_snow(iCol,:), 1, nLev )
+                            re_cloud(iCol,:), re_ice(iCol,:), re_snow(iCol,:), ilsmask,  & 
+                            1, nLev )
        do iLay = 1, nLev
           re_cloud(iCol,iLay) = MAX(re_qc_min, MIN(re_cloud(iCol,iLay), re_qc_max))
           re_ice(iCol,iLay)   = MAX(re_qi_min, MIN(re_ice(iCol,iLay),   re_qi_max))

physics/GFS_stochastics.F90

@@ -67,7 +67,7 @@ subroutine GFS_stochastics_run (im, km, kdt, delt, do_sppt, pert_mp, use_zmtnblc
                                       gu0, gv0, gt0, gq0_wv, dtdtnp,                     &
                                       gq0_cw, gq0_rw, gq0_sw, gq0_iw, gq0_gl,            &
                                       rain, rainc, tprcp, totprcp, cnvprcp,              &
-                                      totprcpb, cnvprcpb, cplflx,                        &
+                                      totprcpb, cnvprcpb, cplflx, cpllnd,                &
                                       rain_cpl, snow_cpl, drain_cpl, dsnow_cpl,          &
                                       ntcw,ntrw,ntsw,ntiw,ntgl,                          &
                                       errmsg, errflg)
@@ -128,8 +128,10 @@ subroutine GFS_stochastics_run (im, km, kdt, delt, do_sppt, pert_mp, use_zmtnblc
          real(kind_phys), dimension(:),         intent(inout) :: totprcpb
          real(kind_phys), dimension(:),         intent(inout) :: cnvprcpb
          logical,                               intent(in)    :: cplflx
-         ! rain_cpl, snow_cpl only allocated if cplflx == .true. or cplchm == .true.
+         logical,                               intent(in)    :: cpllnd
+         ! rain_cpl only allocated if cplflx == .true. or cplchm == .true. or cpllnd == .true.
          real(kind_phys), dimension(:),         intent(inout) :: rain_cpl
+         ! snow_cpl only allocated if cplflx == .true. or cplchm == .true.
          real(kind_phys), dimension(:),         intent(inout) :: snow_cpl
          ! drain_cpl, dsnow_cpl only allocated if cplflx == .true. or cplchm == .true.
          real(kind_phys), dimension(:),         intent(in)    :: drain_cpl
@@ -240,8 +242,10 @@ subroutine GFS_stochastics_run (im, km, kdt, delt, do_sppt, pert_mp, use_zmtnblc
            totprcpb(:) = totprcpb(:) + (sppt_wts(:,15) - 1 )*rain(:)
            cnvprcpb(:) = cnvprcpb(:) + (sppt_wts(:,15) - 1 )*rainc(:)
 
-           if (cplflx) then
+           if (cplflx .or. cpllnd) then
                rain_cpl(:) = rain_cpl(:) + (sppt_wts(:,15) - 1.0)*drain_cpl(:)
+           endif
+           if (cplflx) then
                snow_cpl(:) = snow_cpl(:) + (sppt_wts(:,15) - 1.0)*dsnow_cpl(:)
            endif
            !zero out radiative heating tendency for next physics step
@@ -342,8 +346,10 @@ subroutine GFS_stochastics_run (im, km, kdt, delt, do_sppt, pert_mp, use_zmtnblc
             totprcpb(:)      = totprcpb(:)      + (ca(:,15) - 1 )*rain(:)
             cnvprcpb(:)      = cnvprcpb(:)      + (ca(:,15) - 1 )*rainc(:)
 
-            if (cplflx) then
+            if (cplflx .or. cpllnd) then
                rain_cpl(:) = rain_cpl(:) + (ca(:,15) - 1.0)*drain_cpl(:)
+            endif
+            if (cplflx) then
                snow_cpl(:) = snow_cpl(:) + (ca(:,15) - 1.0)*dsnow_cpl(:)
             endif
             !zero out radiative heating tendency for next physics step

physics/GFS_stochastics.meta

@@ -419,6 +419,13 @@
   dimensions = ()
   type = logical
   intent = in
+[cpllnd]
+  standard_name = flag_for_land_coupling
+  long_name = flag controlling cpllnd collection (default off)
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
 [rain_cpl]
   standard_name = cumulative_lwe_thickness_of_precipitation_amount_for_coupling
   long_name = total rain precipitation

physics/GFS_surface_generic_post.F90

@@ -44,7 +44,7 @@ end subroutine GFS_surface_generic_post_init
 !> \section arg_table_GFS_surface_generic_post_run Argument Table
 !! \htmlinclude GFS_surface_generic_post_run.html
 !!
-      subroutine GFS_surface_generic_post_run (im, cplflx, cplaqm, cplchm, cplwav, lssav, dry, icy, wet,                            &
+      subroutine GFS_surface_generic_post_run (im, cplflx, cplaqm, cplchm, cplwav, cpllnd, lssav, dry, icy, wet,                    &
         lsm, lsm_noahmp, dtf, ep1d, gflx, tgrs_1, qgrs_1, ugrs_1, vgrs_1,                                                           &
         adjsfcdlw, adjsfcdsw, adjnirbmd, adjnirdfd, adjvisbmd, adjvisdfd, adjsfculw, adjsfculw_wat, adjnirbmu, adjnirdfu,           &
         adjvisbmu, adjvisdfu, t2m, q2m, u10m, v10m, tsfc, tsfc_wat, pgr, xcosz, evbs, evcw, trans, sbsno, snowc, snohf, pah, pahi,  &
@@ -58,7 +58,7 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplaqm, cplchm, cplwav, lss
         implicit none
 
         integer,                                intent(in) :: im
-        logical,                                intent(in) :: cplflx, cplaqm, cplchm, cplwav, lssav
+        logical,                                intent(in) :: cplflx, cplaqm, cplchm, cplwav, cpllnd, lssav
         logical, dimension(:),                  intent(in) :: dry, icy, wet
         integer,                                intent(in) :: lsm, lsm_noahmp
         real(kind=kind_phys),                   intent(in) :: dtf
@@ -121,18 +121,24 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplaqm, cplchm, cplwav, lss
           enddo
         endif
 
-        if (cplflx .or. cplchm) then
+        if (cplflx .or. cplchm .or. cpllnd) then
           do i=1,im
             tsfci_cpl(i) = tsfc(i)
           enddo
         endif
 
+        if (cplflx .or. cpllnd) then
+          do i=1,im
+            dlwsfc_cpl  (i) = dlwsfc_cpl(i) + adjsfcdlw(i)*dtf
+            dswsfc_cpl  (i) = dswsfc_cpl(i) + adjsfcdsw(i)*dtf
+            psurfi_cpl  (i) = pgr(i)
+          enddo
+        endif
+
         if (cplflx) then
           do i=1,im
             dlwsfci_cpl (i) = adjsfcdlw(i)
             dswsfci_cpl (i) = adjsfcdsw(i)
-            dlwsfc_cpl  (i) = dlwsfc_cpl(i) + adjsfcdlw(i)*dtf
-            dswsfc_cpl  (i) = dswsfc_cpl(i) + adjsfcdsw(i)*dtf
             dnirbmi_cpl (i) = adjnirbmd(i)
             dnirdfi_cpl (i) = adjnirdfd(i)
             dvisbmi_cpl (i) = adjvisbmd(i)
@@ -148,12 +154,13 @@ subroutine GFS_surface_generic_post_run (im, cplflx, cplaqm, cplchm, cplwav, lss
             nlwsfc_cpl  (i) = nlwsfc_cpl(i) + nlwsfci_cpl(i)*dtf
             t2mi_cpl    (i) = t2m(i)
             q2mi_cpl    (i) = q2m(i)
-            psurfi_cpl  (i) = pgr(i)
           enddo
+        endif
 
 !  ---  estimate mean albedo for ocean point without ice cover and apply
 !       them to net SW heat fluxes
 
+        if (cplflx .or. cpllnd) then
           do i=1,im
 !           if (Sfcprop%landfrac(i) < one) then ! Not 100% land
             if (wet(i)) then                    ! some open water