Merge pull request NOAA-GSL#18 from grantfirl/fix_SCM_specified_surfa…

…ce_flux_bug

Fix scm specified surface flux bug

CODEOWNERS

@@ -11,6 +11,8 @@
 #   https://docs.google.com/spreadsheets/d/14y0Th_sSpCqlssEMNfSZ_Ni9wrpPqfpPY0kRG7jCZB8/edit#gid=0
 # (Internal NOAA document.)
 
+smoke/*                                        @haiqinli                            @climbfuji @SamuelTrahanNOAA @grantfirl @mzhangw @panll @mkavulich @ChunxiZhang-NOAA
+
 physics/cs_conv_aw_adj.*                       @AnningCheng-NOAA                    @climbfuji @SamuelTrahanNOAA @grantfirl @mzhangw @panll @mkavulich @ChunxiZhang-NOAA
 physics/cs_conv.*                              @AnningCheng-NOAA                    @climbfuji @SamuelTrahanNOAA @grantfirl @mzhangw @panll @mkavulich @ChunxiZhang-NOAA
 physics/cu_gf*                                 @hannahcbarnes @haiqinli             @climbfuji @SamuelTrahanNOAA @grantfirl @mzhangw @panll @mkavulich @ChunxiZhang-NOAA

physics/GFS_debug.F90

@@ -593,6 +593,9 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling,
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Tbd%acvb'            , Tbd%acvb)
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Tbd%acvt'            , Tbd%acvt)
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Tbd%hpbl'            , Tbd%hpbl)
+                     if(Model%imfdeepcnv>0 .or. Model%imfshalcnv>0) then
+                       call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Tbd%ud_mf'         , Tbd%ud_mf)
+                     endif
                      if (Model%do_sppt) then
                        call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Tbd%dtdtnp'        , Tbd%dtdtnp)
                        call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Tbd%dtotprcp'      , Tbd%dtotprcp)
@@ -723,7 +726,7 @@ subroutine GFS_diagtoscreen_run (Model, Statein, Stateout, Sfcprop, Coupling,
                      end if
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%dkt         ',    Diag%dkt)
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%dku         ',    Diag%dku)
-                     ! CCPP/MYNNPBL only
+                     ! CCPP/MYNNEDMF only
                      if (Model%do_mynnedmf) then
                        if (Model%bl_mynn_output .ne. 0) then
                          call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Diag%edmf_a      ',  Diag%edmf_a)
@@ -1331,7 +1334,6 @@ subroutine GFS_interstitialtoscreen_run (Model, Statein, Stateout, Sfcprop, Coup
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%tsurf_ice           ', Interstitial%tsurf_ice               )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%tsurf_land          ', Interstitial%tsurf_land              )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%tsurf_water         ', Interstitial%tsurf_water             )
-                     call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%ud_mf               ', Interstitial%ud_mf                   )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%uustar_ice          ', Interstitial%uustar_ice              )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%uustar_land         ', Interstitial%uustar_land             )
                      call print_var(mpirank, omprank, blkno, Grid%xlat_d, Grid%xlon_d, 'Interstitial%uustar_water        ', Interstitial%uustar_water            )

physics/GFS_rrtmg_pre.F90

@@ -40,7 +40,8 @@ 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,           &
-        spp_wts_rad, spp_rad, errmsg, errflg)
+        aero_dir_fdb, smoke_ext, dust_ext,                                     &
+        spp_wts_rad, spp_rad, rrfs_smoke_band, errmsg, errflg)
 
       use machine,                   only: kind_phys
 
@@ -108,13 +109,16 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp,        &
          iovr_exprand,                     & ! Flag for exponential-random cloud overlap method
          idcor_con,                        &
          idcor_hogan,                      &
-         idcor_oreopoulos                            
+         idcor_oreopoulos,                 &
+         rrfs_smoke_band                     ! Band number for rrfs-smoke dust and smoke
 
       character(len=3), dimension(:), intent(in) :: lndp_var_list
 
       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
 
       logical,              intent(in) :: nssl_ccn_on, nssl_invertccn
       integer,              intent(in) :: spp_rad
@@ -616,6 +620,16 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp,        &
         enddo
        enddo
 
+      !> 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
+            ! 550nm (~18000/cm)
+            faersw1(i,k,rrfs_smoke_band) = faersw1(i,k,rrfs_smoke_band) + 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
@@ -763,21 +777,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp,        &
             enddo
           endif
         elseif (imp_physics == imp_physics_gfdl) then            ! GFDL MP
-          if ((imfdeepcnv==imfdeepcnv_gf .or. do_mynnedmf) .and. kdt>1) then
-            if (do_mynnedmf) 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,ntclamt)
-                  else
-                  ! MYNN sub-grid cloud fraction
-                    cldcov(i,k1) = clouds1(i,k1)
-                  endif
-                enddo
-              enddo
-            else ! imfdeepcnv==imfdeepcnv_gf
+          if ((imfdeepcnv==imfdeepcnv_gf) .and. kdt>1) then
               do k=1,lm
                 k1 = k + kd
                 do i=1,im
@@ -789,7 +789,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp,        &
                 endif
                 enddo
               enddo
-            endif
           else
             ! GFDL cloud fraction
             cldcov(1:IM,1+kd:LM+kd) = tracer1(1:IM,1:LM,ntclamt)

physics/GFS_rrtmg_pre.meta

@@ -1194,6 +1194,29 @@
   type = real
   kind = kind_phys
   intent = out
+[aero_dir_fdb]
+  standard_name = do_smoke_aerosol_direct_feedback
+  long_name = flag for smoke and dust radiation feedback
+  units = flag
+  dimensions = ()
+  type = logical
+  intent = in
+[smoke_ext]
+  standard_name = extinction_coefficient_in_air_due_to_smoke
+  long_name = extinction coefficient in air due to smoke
+  units = various
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = in
+[dust_ext]
+  standard_name = extinction_coefficient_in_air_due_to_dust
+  long_name = extinction coefficient in air due to dust
+  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
@@ -1209,6 +1232,13 @@
   dimensions = ()
   type = integer
   intent = in
+[rrfs_smoke_band]
+  standard_name = index_of_shortwave_band_affected_by_smoke
+  long_name = rrtmg band number that smoke and dust should affect
+  units = count
+  dimensions = ()
+  type = integer
+  intent = in
 [errmsg]
   standard_name = ccpp_error_message
   long_name = error message for error handling in CCPP

physics/GFS_surface_composites_post.F90

@@ -264,6 +264,7 @@ subroutine GFS_surface_composites_post_run (
 
         do i=1,im
           if (islmsk(i) == 1) then
+          !-- land
             zorl(i)   = zorll(i)
             cd(i)     = cd_lnd(i)
             cdq(i)    = cdq_lnd(i)
@@ -289,6 +290,7 @@ subroutine GFS_surface_composites_post_run (
             hice(i)   = zero
             cice(i)   = zero
           elseif (islmsk(i) == 0) then
+          !-- water
             zorl(i)   = zorlo(i)
             cd(i)     = cd_wat(i)
             cdq(i)    = cdq_wat(i)
@@ -315,6 +317,7 @@ subroutine GFS_surface_composites_post_run (
             hice(i)   = zero
             cice(i)   = zero
           else ! islmsk(i) == 2
+          !-- ice
             zorl(i)   = zorli(i)
             cd(i)     = cd_ice(i)
             cdq(i)    = cdq_ice(i)

physics/GFS_surface_composites_pre.F90

@@ -21,8 +21,9 @@ module GFS_surface_composites_pre
 !> \section arg_table_GFS_surface_composites_pre_run Argument Table
 !! \htmlinclude GFS_surface_composites_pre_run.html
 !!
-   subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, frac_grid,                                &
-                                 flag_cice, cplflx, cplice, cplwav2atm, landfrac, lakefrac, lakedepth, oceanfrac, frland, &
+   subroutine GFS_surface_composites_pre_run (im, lkm, frac_grid,                                                         &
+                                 flag_cice, cplflx, cplice, cplwav2atm, lsm, lsm_ruc,                                     &
+                                 landfrac, lakefrac, lakedepth, oceanfrac, frland,                                        &
                                  dry, icy, lake, use_flake, wet, hice, cice, zorlo, zorll, zorli,                         &
                                  snowd,            snowd_lnd, snowd_ice, tprcp, tprcp_wat,                                &
                                  tprcp_lnd, tprcp_ice, uustar, uustar_wat, uustar_lnd, uustar_ice,                        &
@@ -34,8 +35,8 @@ subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, fra
       implicit none
 
       ! Interface variables
-      integer,                             intent(in   ) :: im, lkm, kdt
-      logical,                             intent(in   ) :: flag_init, flag_restart, frac_grid, cplflx, cplice, cplwav2atm
+      integer,                             intent(in   ) :: im, lkm, kdt, lsm, lsm_ruc
+      logical,                             intent(in   ) :: cplflx, cplice, cplwav2atm, frac_grid
       logical, dimension(:),              intent(inout)  :: flag_cice
       logical,              dimension(:), intent(inout)  :: dry, icy, lake, use_flake, wet
       real(kind=kind_phys), dimension(:), intent(in   )  :: landfrac, lakefrac, lakedepth, oceanfrac
@@ -195,12 +196,13 @@ subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, fra
             endif
           endif
         enddo
-      endif
+      endif ! frac_grid
 
       do i=1,im
         tprcp_wat(i) = tprcp(i)
         tprcp_lnd(i) = tprcp(i)
         tprcp_ice(i) = tprcp(i)
+
         if (wet(i)) then                   ! Water
           uustar_wat(i) = uustar(i)
             tsfc_wat(i) = tsfco(i)
@@ -213,7 +215,7 @@ subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, fra
         endif
         if (dry(i)) then                   ! Land
           uustar_lnd(i) = uustar(i)
-           weasd_lnd(i) = weasd(i)
+           if(lsm /= lsm_ruc) weasd_lnd(i) = weasd(i)
            tsurf_lnd(i) = tsfcl(i)
         ! DH*
         else
@@ -224,7 +226,7 @@ subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, fra
         endif
         if (icy(i)) then                   ! Ice
           uustar_ice(i) = uustar(i)
-           weasd_ice(i) = weasd(i)
+           if(lsm /= lsm_ruc) weasd_ice(i) = weasd(i)
            tsurf_ice(i) = tisfc(i)
             ep1d_ice(i) = zero
             gflx_ice(i) = zero
@@ -272,7 +274,7 @@ subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, fra
             endif
           endif
         enddo
-      else
+      elseif(lsm /= lsm_ruc) then ! do not do snow initialization  with RUC lsm
         do i=1,im
           if (icy(i)) then
             if (kdt == 1 .or. (.not. cplflx .or. lakefrac(i) > zero)) then
@@ -290,4 +292,4 @@ subroutine GFS_surface_composites_pre_run (im, flag_init, flag_restart, lkm, fra
 
    end subroutine GFS_surface_composites_pre_run
 
-end module GFS_surface_composites_pre
+end module GFS_surface_composites_pre

physics/GFS_surface_composites_pre.meta

@@ -14,20 +14,6 @@
   dimensions = ()
   type = integer
   intent = in
-[flag_init]
-  standard_name = flag_for_first_timestep
-  long_name = flag signaling first time step for time integration loop
-  units = flag
-  dimensions = ()
-  type = logical
-  intent = in
-[flag_restart]
-  standard_name = flag_for_restart
-  long_name = flag for restart (warmstart) or coldstart
-  units = flag
-  dimensions = ()
-  type = logical
-  intent = in
 [lkm]
   standard_name = control_for_lake_surface_scheme
   long_name = flag for lake surface model
@@ -70,6 +56,20 @@
   dimensions = ()
   type = logical
   intent = in
+[lsm]
+  standard_name = control_for_land_surface_scheme
+  long_name = flag for land surface model
+  units = flag
+  dimensions = ()
+  type = integer
+  intent = in
+[lsm_ruc]
+  standard_name = identifier_for_ruc_land_surface_scheme
+  long_name = flag for RUC land surface model
+  units = flag
+  dimensions = ()
+  type = integer
+  intent = in
 [landfrac]
   standard_name = land_area_fraction
   long_name = fraction of horizontal grid area occupied by land

physics/bl_mynn_common.f90

@@ -0,0 +1,67 @@
+!>\file bl_mynn_common.f90
+!! Define Model-specific constants/parameters.
+!! This module will be used at the initialization stage
+!! where all model-specific constants are read and saved into
+!! memory. This module is then used again in the MYNN-EDMF. All
+!! MYNN-specific constants are declared globally in the main 
+!! module (module_bl_mynn) further below:
+ module bl_mynn_common
+
+!------------------------------------------
+!
+!------------------------------------------
+
+! The following 5-6 lines are the only lines in this file that are not 
+! universal for all dycores... Any ideas how to universalize it?
+! For MPAS:
+! use mpas_kind_types,only: kind_phys => RKIND
+! For CCPP:
+  use machine,  only : kind_phys
+
+ implicit none
+ save
+
+! To be specified from dycore
+ real(kind=kind_phys):: cp           != 7.*r_d/2. (J/kg/K)
+ real(kind=kind_phys):: cpv          != 4.*r_v    (J/kg/K) Spec heat H2O gas
+ real(kind=kind_phys):: cice         != 2106.     (J/kg/K) Spec heat H2O ice
+ real(kind=kind_phys):: cliq         != 4190.     (J/kg/K) Spec heat H2O liq
+ real(kind=kind_phys):: p608         != R_v/R_d-1.
+ real(kind=kind_phys):: ep_2         != R_d/R_v
+ real(kind=kind_phys):: grav         != accel due to gravity
+ real(kind=kind_phys):: karman       != von Karman constant
+ real(kind=kind_phys):: t0c          != temperature of water at freezing, 273.15 K
+ real(kind=kind_phys):: rcp          != r_d/cp
+ real(kind=kind_phys):: r_d          != 287.  (J/kg/K) gas const dry air
+ real(kind=kind_phys):: r_v          != 461.6 (J/kg/K) gas const water
+ real(kind=kind_phys):: xlf          != 0.35E6 (J/kg) fusion at 0 C
+ real(kind=kind_phys):: xlv          != 2.50E6 (J/kg) vaporization at 0 C
+ real(kind=kind_phys):: xls          != 2.85E6 (J/kg) sublimation
+ real(kind=kind_phys):: rvovrd       != r_v/r_d != 1.608
+
+! Specified locally
+ real(kind=kind_phys),parameter:: zero   = 0.0
+ real(kind=kind_phys),parameter:: half   = 0.5
+ real(kind=kind_phys),parameter:: one    = 1.0
+ real(kind=kind_phys),parameter:: two    = 2.0
+ real(kind=kind_phys),parameter:: onethird  = 1./3.
+ real(kind=kind_phys),parameter:: twothirds = 2./3.
+ real(kind=kind_phys),parameter:: tref  = 300.0   ! reference temperature (K)
+ real(kind=kind_phys),parameter:: TKmin = 253.0   ! for total water conversion, Tripoli and Cotton (1981)
+ real(kind=kind_phys),parameter:: p1000mb=100000.0
+ real(kind=kind_phys),parameter:: svp1  = 0.6112 !(kPa)
+ real(kind=kind_phys),parameter:: svp2  = 17.67  !(dimensionless)
+ real(kind=kind_phys),parameter:: svp3  = 29.65  !(K)
+ real(kind=kind_phys),parameter:: tice  = 240.0  !-33 (C), temp at saturation w.r.t. ice
+
+! To be derived in the init routine
+ real(kind=kind_phys):: ep_3         != 1.-ep_2 != 0.378
+ real(kind=kind_phys):: gtr          != grav/tref
+ real(kind=kind_phys):: rk           != cp/r_d
+ real(kind=kind_phys):: tv0          != p608*tref
+ real(kind=kind_phys):: tv1          != (1.+p608)*tref
+ real(kind=kind_phys):: xlscp        != (xlv+xlf)/cp
+ real(kind=kind_phys):: xlvcp        != xlv/cp
+ real(kind=kind_phys):: g_inv        != 1./grav
+
+ end module bl_mynn_common

physics/cu_gf_driver.F90

@@ -39,26 +39,6 @@ subroutine cu_gf_driver_init(imfshalcnv, imfshalcnv_gf, imfdeepcnv, &
          errmsg = ''
          errflg = 0
 
-         ! DH* temporary
-         ! if (mpirank==mpiroot) then
-         !    write(0,*) ' ----------------------------------------------------------'//&
-         !                '-------------------------------------------------------------------'
-         !    write(0,*) ' --- WARNING --- the CCPP Grell Freitas convection scheme is'//&
-         !                ' currently under development, use at your own risk --- WARNING ---'
-         !    write(0,*) ' --------------------------------------------------------------------'//&
-         !                '---------------------------------------------------------'
-         ! end if
-         ! *DH temporary
-
-         ! Consistency checks
-         if (.not. (imfshalcnv == imfshalcnv_gf .or.                       &
-        &        imfdeepcnv == imfdeepcnv_gf)) then
-           write(errmsg,'(*(a))') 'Logic error: namelist choice of',       &
-        &    ' convection is different from Grell-Freitas scheme'
-           errflg = 1
-           return
-         end if
-
       end subroutine cu_gf_driver_init
 
       subroutine cu_gf_driver_finalize()