update main with the develop

CMakeLists.txt

@@ -6,7 +6,7 @@ project(ccpp_physics
 
 #------------------------------------------------------------------------------
 set(PACKAGE "ccpp-physics")
-set(AUTHORS "Grant Firl" "Dom Heinzeller" "Man Zhang" "Mike Kavulich" "Chunxi Zhang")
+set(AUTHORS "Grant Firl" "Dustin Swales" "Man Zhang" "Mike Kavulich" )
 
 #------------------------------------------------------------------------------
 # Set OpenMP flags for C/C++/Fortran
@@ -183,7 +183,9 @@ set_target_properties(ccpp_physics PROPERTIES VERSION ${PROJECT_VERSION}
 target_include_directories(ccpp_physics PUBLIC
                            $<BUILD_INTERFACE:${CMAKE_CURRENT_BINARY_DIR}>)
 
-target_link_libraries(ccpp_physics PUBLIC w3emc::w3emc_d NetCDF::NetCDF_Fortran)
+target_link_libraries(ccpp_physics PUBLIC w3emc::w3emc_d
+                                          sp::sp_d
+                                          NetCDF::NetCDF_Fortran)
 
 # Define where to install the library
 install(TARGETS ccpp_physics

physics/GFS_rrtmg_pre.F90

@@ -20,7 +20,7 @@ module GFS_rrtmg_pre
       subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,&
         ltp, imfdeepcnv, imfdeepcnv_gf, me, ncnd, ntrac, num_p3d, npdf3d,      &
         ncnvcld3d,ntqv, ntcw,ntiw, ntlnc, ntinc, ntrnc, ntsnc, ntccn, top_at_1,&
-        ntrw, ntsw, ntgl, nthl, ntwa, ntoz,                                    &
+        ntrw, ntsw, ntgl, nthl, ntwa, ntoz, ntsmoke, ntdust, ntcoarsepm,       &
         ntclamt, nleffr, nieffr, nseffr, lndp_type, kdt,                       &
         ntdu1, ntdu2, ntdu3, ntdu4, ntdu5, ntss1, ntss2,                       &
         ntss3, ntss4, ntss5, ntsu, ntbcb, ntbcl, ntocb, ntocl, ntchm,          &
@@ -41,11 +41,10 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,&
         kd, kt, kb, mtopa, mbota, raddt, tsfg, tsfa, de_lgth, alb1d, delp, dz, & !output from here and below
         plvl, plyr, tlvl, tlyr, qlyr, olyr, gasvmr_co2, gasvmr_n2o, gasvmr_ch4,&
         gasvmr_o2, gasvmr_co, gasvmr_cfc11, gasvmr_cfc12, gasvmr_cfc22,        &
-        gasvmr_ccl4,  gasvmr_cfc113, aerodp, clouds6, clouds7, clouds8,        &
+        gasvmr_ccl4,  gasvmr_cfc113, aerodp,ext550, 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, rrfs_smoke_band, ico2, errmsg, errflg)
+        aero_dir_fdb, fdb_coef, spp_wts_rad, spp_rad, ico2, errmsg, errflg)
 
       use machine,                   only: kind_phys
 
@@ -89,6 +88,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,&
                                            ntcw, ntiw, ntlnc, ntinc,           &
                                            ntrnc, ntsnc,ntccn,                 &
                                            ntrw, ntsw, ntgl, nthl, ntwa, ntoz, &
+                                           ntsmoke, ntdust, ntcoarsepm,        &
                                            ntclamt, nleffr, nieffr, nseffr,    &
                                            lndp_type,                          &
                                            kdt, imp_physics,                   &
@@ -113,7 +113,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,&
          idcor,                            &
          idcor_hogan,                      &
          idcor_oreopoulos,                 &
-         rrfs_smoke_band,                  & ! Band number for rrfs-smoke dust and smoke
          ico2                                ! Flag for co2 source used in radiation
 
       integer, intent(in) :: ntdu1, ntdu2, ntdu3, ntdu4, ntdu5, ntss1, ntss2, ntss3,  &
@@ -126,7 +125,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,&
                                           lmfshal, lmfdeep2, pert_clds, lcrick,&
                                           lcnorm, top_at_1, lextop, mraerosol
       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
@@ -157,6 +155,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,&
                                                              clouds2, clouds3, &
                                                              clouds4, clouds5
       real(kind=kind_phys), dimension(:,:), intent(in)  :: qci_conv
+      real(kind=kind_phys), dimension(:),   intent(in)  :: fdb_coef
       real(kind=kind_phys), dimension(:),   intent(out) :: lwp_ex,iwp_ex, &
                                                            lwp_fc,iwp_fc
 
@@ -187,6 +186,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,&
                                                            gasvmr_ccl4,&
                                                            gasvmr_cfc113
       real(kind=kind_phys), dimension(:,:), intent(out) :: aerodp
+      real(kind=kind_phys), dimension(:,:), intent(out) :: ext550
       real(kind=kind_phys), dimension(:,:), intent(out) :: clouds6,   &
                                                            clouds7,   &
                                                            clouds8,   &
@@ -637,13 +637,28 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,&
           enddo
         endif
 
+!>---   add smoke and dust ---
+       if (aero_dir_fdb) then
+         do k=1,lmk
+           do i=1,im
+             aer_nm(i,k,1 )=aer_nm(i,k,1 )+ qgrs(i,k,ntdust)*fdb_coef(1)*1.e-9    ! dust bin1
+             aer_nm(i,k,2 )=aer_nm(i,k,2 )+(qgrs(i,k,ntdust)*fdb_coef(2)          &
+                           +qgrs(i,k,ntcoarsepm)*fdb_coef(3))*1.e-9               ! dust bin2
+             aer_nm(i,k,3 )=aer_nm(i,k,3 )+qgrs(i,k,ntcoarsepm)*fdb_coef(4)*1.e-9 ! dust bin3
+             aer_nm(i,k,4 )=aer_nm(i,k,4 )+qgrs(i,k,ntcoarsepm)*fdb_coef(5)*1.e-9 ! dust bin4
+             aer_nm(i,k,12)=aer_nm(i,k,12)+qgrs(i,k,ntsmoke)*fdb_coef(6)*1.e-9    ! Smoke BC
+             aer_nm(i,k,14)=aer_nm(i,k,14)+qgrs(i,k,ntsmoke)*fdb_coef(7)*1.e-9    ! Smoke OA
+            enddo
+          enddo
+       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, iaermdl, iaerflg, top_at_1, con_pi,  &
-                   con_rd, con_g, faersw, faerlw, aerodp, errflg, errmsg)         !  ---  outputs
+                   con_rd, con_g, faersw, faerlw, aerodp, ext550, errflg, errmsg)         !  ---  outputs
 
 ! CCPP
       do j = 1,NBDSW
@@ -657,16 +672,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, lextop,&
         enddo
        enddo
 
-      !> - 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
-            ! 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

physics/GFS_rrtmg_pre.meta

@@ -219,6 +219,27 @@
   dimensions = ()
   type = integer
   intent = in
+[ntsmoke]
+  standard_name = index_for_smoke_in_tracer_concentration_array
+  long_name = tracer index for smoke
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[ntdust]
+  standard_name = index_for_dust_in_tracer_concentration_array
+  long_name = tracer index for dust
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
+[ntcoarsepm]
+  standard_name = index_for_coarse_particulate_matter_in_tracer_concentration_array
+  long_name = tracer index for coarse particulate matter
+  units = index
+  dimensions = ()
+  type = integer
+  intent = in
 [iaermdl]
   standard_name = control_for_aerosol_radiation_scheme
   long_name = control of aerosol scheme in radiation
@@ -1264,6 +1285,14 @@
   type = real
   kind = kind_phys
   intent = out
+[ext550]
+  standard_name = aerosol_optical_depth_at_550nm
+  long_name = 3d optical extinction for total aerosol species
+  units = none
+  dimensions = (horizontal_loop_extent,vertical_layer_dimension)
+  type = real
+  kind = kind_phys
+  intent = out
 [clouds6]
   standard_name = cloud_rain_water_path
   long_name = cloud rain water path
@@ -1437,19 +1466,11 @@
   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)
+[fdb_coef]
+  standard_name = smoke_dust_direct_fdb_coef
+  long_name = smoke dust direct feedback coefficents
+  units = none
+  dimensions = (7)
   type = real
   kind = kind_phys
   intent = in
@@ -1468,13 +1489,6 @@
   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/cu_gf_deep.F90

@@ -25,8 +25,8 @@ module cu_gf_deep
      real(kind=kind_phys), parameter :: pgcd = 0.1
 !
 !> aerosol awareness, do not use yet!
-     integer, parameter :: autoconv=2
-     integer, parameter :: aeroevap=3
+     integer, parameter :: autoconv=1 !2
+     integer, parameter :: aeroevap=1 !3
      real(kind=kind_phys), parameter :: scav_factor = 0.5
 !> still 16 ensembles for clousres
      integer, parameter:: maxens3=16

physics/maximum_hourly_diagnostics.F90

@@ -29,14 +29,17 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics,
                                              gt0, refl_10cm, refdmax, refdmax263k, u10m, v10m,     &
                                              u10max, v10max, spd10max, pgr, t2m, q2m, t02max,      &
                                              t02min, rh02max, rh02min, dtp, rain, pratemax,        &
-                                             errmsg, errflg)
+                                             lightning_threat, ltg1_max,ltg2_max,ltg3_max,         &
+                                             wgrs, prsi, qgraupel, qsnowwat, qicewat, tgrs, con_rd,&
+                                             prsl, kdt, errmsg, errflg)
 
        ! Interface variables
-       integer, intent(in) :: im, levs
-       logical, intent(in) :: reset, lradar
+       integer, intent(in) :: im, levs, kdt
+       logical, intent(in) :: reset, lradar, lightning_threat
        integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_fer_hires, &
                               imp_physics_nssl
        real(kind_phys), intent(in   ) :: con_g
+       real(kind_phys), intent(in   ) :: con_rd
        real(kind_phys), intent(in   ) :: phil(:,:)
        real(kind_phys), intent(in   ) :: gt0(:,:)
        real(kind_phys), intent(in   ) :: refl_10cm(:,:)
@@ -55,20 +58,30 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics,
        real(kind_phys), intent(inout) :: rh02max(:)
        real(kind_phys), intent(inout) :: rh02min(:)
        real(kind_phys), intent(in   ) :: dtp
-       real(kind_phys), intent(in   ) :: rain(im)
-       real(kind_phys), intent(inout) :: pratemax(im)
+       real(kind_phys), intent(in   ) :: rain(:)
+       real(kind_phys), intent(in   ) :: tgrs(:,:)
+       real(kind_phys), intent(in   ) :: prsl(:,:)
+       real(kind_phys), intent(inout) :: pratemax(:)
+
+       real(kind_phys), intent(in), dimension(:,:) :: prsi, qgraupel, qsnowwat, qicewat, wgrs
+       real(kind_phys), intent(inout), dimension(:) :: ltg1_max, ltg2_max, ltg3_max
        character(len=*), intent(out)  :: errmsg
        integer, intent(out)           :: errflg
 
        ! Local variables
        real(kind_phys), dimension(:), allocatable :: refd, refd263k
-       real(kind_phys) :: tem, pshltr, QCQ, rh02
+       real(kind_phys) :: tem, pshltr, QCQ, rh02, dP, Q
        integer :: i
 
        ! Initialize CCPP error handling variables
        errmsg = ''
        errflg = 0
 
+!Lightning threat indices
+       if (lightning_threat) then
+         call lightning_threat_indices
+       endif
+
 !Calculate hourly max 1-km agl and -10C reflectivity
        if (lradar .and. (imp_physics == imp_physics_gfdl .or. &
           imp_physics == imp_physics_thompson .or.    &
@@ -134,6 +147,79 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics,
           pratemax(i) = max(pratemax(i),(3.6E6/dtp)*rain(i))
        enddo
 
+     contains
+
+       subroutine lightning_threat_indices
+         implicit none
+         REAL(kind_phys), PARAMETER    :: clim1=1.50
+         REAL(kind_phys), PARAMETER    :: clim2=0.40*1.22
+         REAL(kind_phys), PARAMETER    :: clim3=0.02*1.22
+         !  coef1 and coef2 are modified from the values given
+         !  in McCaul et al. 
+         !  coef1 is x 1000 x 1.22
+         !  coef2 is x 1.22
+         !  are these tuning factors, scale factors??
+         !  McCaul et al. used a 2-km WRF simulation
+         REAL(kind_phys), PARAMETER    :: coef1=0.042*1000.*1.22
+         REAL(kind_phys), PARAMETER    :: coef2=0.20*1.22
+
+         REAL(kind_phys) :: totice_colint(im), ltg1, ltg2, high_ltg1, high_wgrs, high_graupel, rho
+         LOGICAL :: ltg1_calc(im)
+         integer :: k, i, count
+
+         count = 0
+         high_ltg1 = 0
+         high_wgrs = 0
+         high_graupel = 0
+
+          totice_colint = 0
+          ltg1_calc = .false.
+          do k=1,levs-1
+             do i=1,im
+                dP = prsi(i,k) - prsi(i,k+1)
+                Q = qgraupel(i,k) + qsnowwat(i,k) + qicewat(i,k)
+                rho = prsl(i,k) / (con_rd * tgrs(i,k))
+                totice_colint(i) = totice_colint(i) + Q * rho * dP / con_g
+
+                IF ( .not.ltg1_calc(i) ) THEN
+                   IF ( 0.5*(tgrs(i,k+1) + tgrs(i,k)) < 258.15 ) THEN
+                      count = count + 1
+                      ltg1_calc(i) = .true.
+
+                      ltg1 = coef1*wgrs(i,k)* &
+                           (( qgraupel(i,k+1) + qgraupel(i,k) )*0.5 )
+                      if(ltg1 > high_ltg1) then
+                        high_ltg1 = ltg1
+                        high_graupel = qgraupel(i,k)
+                        high_wgrs = wgrs(i,k)
+                      endif
+
+                      IF ( ltg1 .LT. clim1 ) ltg1 = 0.
+
+                      IF ( ltg1 .GT. ltg1_max(i) ) THEN
+                         ltg1_max(i) = ltg1
+                      ENDIF
+                   ENDIF
+                ENDIF
+             enddo
+          enddo
+
+          do i=1,im
+             ltg2 = coef2 * totice_colint(i)
+
+             IF ( ltg2 .LT. clim2 ) ltg2 = 0.
+
+             IF ( ltg2 .GT. ltg2_max(i) ) THEN
+                ltg2_max(i) = ltg2
+             ENDIF
+
+             ltg3_max(i) = 0.95 * ltg1_max(i) + 0.05 * ltg2_max(i)
+
+             IF ( ltg3_max(i) .LT. clim3 ) ltg3_max(i) = 0.
+          enddo
+
+       end subroutine lightning_threat_indices
+
    end subroutine maximum_hourly_diagnostics_run
 
    subroutine max_fields(phil,ref3D,grav,im,levs,refd,tk,refd263k)