Merge pull request #147 from haiqinli/gsl/develop-smoke

gsl/develop RRFS-Smoke

Author: SamuelTrahanNOAA

physics/GFS_rrtmg_pre.F90

@@ -36,6 +36,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp,        &
         gasvmr_ccl4,  gasvmr_cfc113, aerodp, clouds6, clouds7, clouds8,        &
         clouds9, cldsa, cldfra, cldfra2d, lwp_ex,iwp_ex, lwp_fc,iwp_fc,        &
         faersw1, faersw2, faersw3, faerlw1, faerlw2, faerlw3, alpha,           &
+        aero_dir_fdb, smoke_ext, dust_ext,                                     &
         spp_wts_rad, spp_rad, errmsg, errflg)
 
       use machine,                   only: kind_phys
@@ -103,6 +104,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp,        &
       logical,              intent(in) :: lsswr, lslwr, ltaerosol, lgfdlmprad, &
                                           uni_cld, effr_in, do_mynnedmf,       &
                                           lmfshal, lmfdeep2, pert_clds
+      logical,              intent(in) :: aero_dir_fdb
+      real(kind=kind_phys), dimension(:,:), intent(in) :: smoke_ext, dust_ext
 
       integer,    intent(in) :: spp_rad
       real(kind_phys),              intent(in) :: spp_wts_rad(:,:)
@@ -601,6 +604,15 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp,        &
         enddo
        enddo
 
+      if(aero_dir_fdb) then ! add smoke/dust extinctions
+        do k = 1, LMK
+          do i = 1, IM
+            ! 550nm (~18000/cm)
+            faersw1(i,k,10) = faersw1(i,k,10) + MIN(4.,smoke_ext(i,k) + dust_ext(i,k))
+          enddo
+        enddo
+      endif
+
       do j = 1,NBDLW
         do k = 1, LMK
           do i = 1, IM

physics/GFS_rrtmg_pre.meta

@@ -1082,6 +1082,29 @@
   type = real
   kind = kind_phys
   intent = out
+[aero_dir_fdb]
+  standard_name = rrfs_smoke_dust_rad_fdb_opt
+  long_name = flag for rrfs smoke dust rad feedback
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+[smoke_ext]
+  standard_name = smoke_ext
+  long_name = smoke optical extinction
+  units = various
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[dust_ext]
+  standard_name = dust_ext
+  long_name = dust optical extinction
+  units = various
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
 [spp_wts_rad]
   standard_name = spp_weights_for_radiation_scheme
   long_name = spp weights for radiation scheme

physics/module_MYNNPBL_wrapper.F90

@@ -159,6 +159,7 @@ SUBROUTINE mynnedmf_wrapper_run(        &
      &  icloud_bl, do_mynnsfclay,                          &
      &  imp_physics, imp_physics_gfdl,                     &
      &  imp_physics_thompson, imp_physics_wsm6,            &
+     &  chem3d, frp, mix_chem, fire_turb, nchem, ndvel,    &
      &  ltaerosol, spp_wts_pbl, spp_pbl, lprnt, huge, errmsg, errflg  )
 
 ! should be moved to inside the mynn:
@@ -179,11 +180,8 @@ SUBROUTINE mynnedmf_wrapper_run(        &
      logical, intent(in) :: cplflx
 
      !smoke/chem
-     !logical, intent(in) :: mix_chem, fire_turb
-     !integer, intent(in) :: nchem, ndvel, kdvel
-     !for testing only:
-     logical, parameter :: mix_chem=.false., fire_turb=.false.
-     integer, parameter :: nchem=2, ndvel=2, kdvel=1
+     integer, intent(in) :: nchem, ndvel
+     integer, parameter  :: kdvel=1
 
 ! NAMELIST OPTIONS (INPUT):
      logical, intent(in) ::                                 &
@@ -285,17 +283,11 @@ SUBROUTINE mynnedmf_wrapper_run(        &
       real(kind=kind_phys), allocatable :: old_ozone(:,:)
 
 !smoke/chem arrays
-   !   real(kind=kind_phys), dimension(:,:), intent(inout), optional ::   &                                     
-   !  &         qgrs_smoke_conc, qgrs_dust_conc
-   !   real(kind=kind_phys), allocatable, dimension(:,:,:) :: chem3d
-   !   real(kind=kind_phys), dimension(:,:), intent(in), optional :: vdep
-   !   real(kind=kind_phys), dimension(:),   intent(in), optional :: frp, emis_ant_no
-!for testing only
-      real(kind=kind_phys), dimension(im,levs)                      ::   &
-     &         qgrs_smoke_conc, qgrs_dust_conc
-      real(kind=kind_phys), allocatable, dimension(:,:,:) :: chem3d
-      real(kind=kind_phys), dimension(im,ndvel)      :: vdep    !not passed in yet???
-      real(kind=kind_phys), dimension(im)            :: frp, emis_ant_no
+      real(kind_phys), dimension(:), intent(inout) :: frp
+      logical, intent(in) :: mix_chem, fire_turb
+      real(kind=kind_phys), dimension(:,:,:), intent(inout) :: chem3d
+      real(kind=kind_phys), dimension(im)   :: emis_ant_no
+      real(kind=kind_phys), dimension(im,ndvel) :: vdep
 
 !MYNN-2D
       real(kind=kind_phys), dimension(:), intent(in) ::                  &
@@ -361,20 +353,8 @@ SUBROUTINE mynnedmf_wrapper_run(        &
       endif
 
       !initialize arrays for test
-      qgrs_smoke_conc = 1.0
-      qgrs_dust_conc  = 1.0
-      FRP         = 0.
       EMIS_ANT_NO = 0.
       vdep = 0. ! hli for chem dry deposition, 0 temporarily
-      if (mix_chem) then
-         allocate ( chem3d(im,levs,nchem) )
-         do k=1,levs
-         do i=1,im
-            chem3d(i,k,1)=qgrs_smoke_conc(i,k)
-            chem3d(i,k,2)=qgrs_dust_conc (i,k)
-         enddo
-         enddo
-      endif
 
   ! Check incoming moist species to ensure non-negative values
   ! First, create height (dz) and pressure differences (delp) 
@@ -966,10 +946,6 @@ SUBROUTINE mynnedmf_wrapper_run(        &
          deallocate(save_qke_adv)
        endif
 
-       if(allocated(chem3d))then
-         deallocate(chem3d)
-       endif
-
   CONTAINS
 
     SUBROUTINE dtend_helper(itracer,field,mult)

physics/module_MYNNPBL_wrapper.meta

@@ -1313,6 +1313,50 @@
   dimensions = ()
   type = integer
   intent = in
+[chem3d]
+  standard_name = chem3d_mynn_pbl_transport
+  long_name = mynn pbl transport of smoke and dust
+  units = various
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension,2)
+  type = real
+  kind = kind_phys
+  intent = inout
+[frp]
+  standard_name = frp_hourly
+  long_name = hourly frp
+  units = mw
+  dimensions = (horizontal_loop_extent)
+  type = real
+  kind = kind_phys
+  intent = inout
+[mix_chem]
+  standard_name = rrfs_smoke_mynn_tracer_mixing
+  long_name = flag for rrfs smoke mynn tracer mixing
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+[fire_turb]
+  standard_name = rrfs_smoke_mynn_enh_vermix
+  long_name = flag for rrfs smoke mynn enh vermix
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+[nchem]
+  standard_name = number_of_chemical_species_vertically_mixed
+  long_name = number of chemical vertically mixed
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
+[ndvel]
+  standard_name = number_of_chemical_species_deposited
+  long_name = number of chemical pbl deposited
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
 [ltaerosol]
   standard_name = flag_for_aerosol_physics
   long_name = flag for aerosol physics

physics/module_bl_mynn.F90

@@ -5348,7 +5348,7 @@ SUBROUTINE mynn_mix_chem(kts,kte,i,     &
     REAL, DIMENSION(kts:kte) :: rhoinv
     REAL, DIMENSION(kts:kte+1) :: rhoz,khdz
     REAL, PARAMETER :: no_threshold    = 0.1
-    REAL, PARAMETER :: frp_threshold   = 1.0     ! JLS 12/21/21
+    REAL, PARAMETER :: frp_threshold   = 10.0     ! RAR 02/11/22: I increased the frp threshold to enhance mixing
     REAL, PARAMETER :: pblh_threshold  = 250.0
 
     dztop=.5*(dz(kte)+dz(kte-1))
@@ -5362,6 +5362,7 @@ SUBROUTINE mynn_mix_chem(kts,kte,i,     &
     rhoz(kts)  =rho(kts)
     rhoinv(kts)=1./rho(kts)
     khdz(kts)  =rhoz(kts)*dfh(kts)
+
     DO k=kts+1,kte
        rhoz(k)  =(rho(k)*dz(k-1) + rho(k-1)*dz(k))/(dz(k-1)+dz(k))
        rhoz(k)  =  MAX(rhoz(k),1E-4)

physics/module_mp_thompson.F90

@@ -968,7 +968,7 @@ END SUBROUTINE thompson_init
 !> @{
       SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc,     &
                               nwfa, nifa, nwfa2d, nifa2d,             &
-                              tt, th, pii,                            &
+                              aero_ind_fdb, tt, th, pii,              &
                               p, w, dz, dt_in, dt_inner,              &
                               sedi_semi, decfl,                       &
                               RAINNC, RAINNCV,                        &
@@ -1023,6 +1023,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc,     &
       REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: &
                           nc, nwfa, nifa
       REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d, nifa2d
+      LOGICAL, OPTIONAL, INTENT(IN):: aero_ind_fdb
       REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: &
                           re_cloud, re_ice, re_snow
       INTEGER, INTENT(IN) :: rand_perturb_on, kme_stoch
@@ -1446,8 +1447,10 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc,     &
 !.. Changed 13 May 2013 to fake emissions in which nwfa2d is aerosol
 !.. number tendency (number per kg per second).
          if (is_aerosol_aware) then
+            if ( .not. aero_ind_fdb) then
             nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt
             nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt
+            endif
 
             do k = kts, kte
                nc(i,k,j) = nc1d(k)

physics/mp_thompson.F90

@@ -300,7 +300,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd,        &
                               con_eps, convert_dry_rho,            &
                               spechum, qc, qr, qi, qs, qg, ni, nr, &
                               is_aerosol_aware, nc, nwfa, nifa,    &
-                              nwfa2d, nifa2d,                      &
+                              nwfa2d, nifa2d, aero_ind_fdb,        &
                               tgrs, prsl, phii, omega,             &
                               sedi_semi, decfl, dtp, dt_inner,     & 
                               first_time_step, istep, nsteps,      &
@@ -339,6 +339,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd,        &
          real(kind_phys), optional, intent(inout) :: nifa(:,:)
          real(kind_phys), optional, intent(in   ) :: nwfa2d(:)
          real(kind_phys), optional, intent(in   ) :: nifa2d(:)
+         logical,         optional, intent(in   ) :: aero_ind_fdb
          ! State variables and timestep information
          real(kind_phys),           intent(inout) :: tgrs(:,:)
          real(kind_phys),           intent(in   ) :: prsl(:,:)
@@ -634,6 +635,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd,        &
          if (is_aerosol_aware) then
             call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr,        &
                               nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d,     &
+                              aero_ind_fdb=aero_ind_fdb,                                     &
                               tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner,  &
                               sedi_semi=sedi_semi, decfl=decfl,                              &
                               rainnc=rain_mp, rainncv=delta_rain_mp,                         &

physics/mp_thompson.meta

@@ -435,6 +435,13 @@
   type = real
   kind = kind_phys
   intent = in
+[aero_ind_fdb]
+  standard_name = rrfs_smoke_aero_ind_fdb_opt
+  long_name = flag for rrfs wfa ifa emission
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
 [tgrs]
   standard_name = air_temperature_of_new_state
   long_name = model layer mean temperature