Implement NOAA/ARL FENGSHA dust scheme

CHANGELOG.md

@@ -23,6 +23,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ### Added
 
 - Compute and export PM2.5 and PM10 diagnostic tracers in UFS interface
+- Add NOAA/ARL FENGSHA dust scheme
 
 ## [1.0.1] - 2021-03-22
 

ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_GridCompMod.F90

@@ -38,15 +38,27 @@ module DU2G_GridCompMod
 
    integer, parameter         :: NHRES = 6
 
+!  !Supported dust schemes
+   enum, bind(C)
+     enumerator :: DUST_SCHEME_DATA    = 0
+     enumerator :: DUST_SCHEME_GOCART
+     enumerator :: DUST_SCHEME_FENGSHA
+   end enum
+
 !  !Dust state
    type, extends(GA_GridComp) :: DU2G_GridComp
-       real, allocatable      :: rlow(:)        ! particle effective radius lower bound [um]
-       real, allocatable      :: rup(:)         ! particle effective radius upper bound [um]
+       real, allocatable      :: rlow(:)        ! particle radius lower bound [um]
+       real, allocatable      :: rup(:)         ! particle radius upper bound [um]
        real, allocatable      :: sfrac(:)       ! fraction of total source
+       real, allocatable      :: sdist(:)       ! aerosol fractional size distribution [1]
+       real                   :: alpha          ! FENGSHA scaling factor
+       real                   :: gamma          ! FENGSHA tuning exponent
+       real                   :: kvhmax         ! FENGSHA max. vertical/horizontal mass flux ratio [1]
        real                   :: Ch_DU_res(NHRES) ! resolutions used for Ch_DU
        real                   :: Ch_DU          ! dust emission tuning coefficient [kg s2 m-5].
        logical                :: maringFlag=.false.  ! maring settling velocity correction
        integer                :: day_save = -1      
+       character(len=:), allocatable :: emission_scheme     ! emission scheme selector
 !      !Workspae for point emissions
        logical                :: doing_point_emissions = .false.
        character(len=255)     :: point_emissions_srcfilen   ! filename for pointwise emissions
@@ -98,6 +110,7 @@ subroutine SetServices (GC, RC)
     type (DU2G_GridComp), pointer      :: self
 
     character (len=ESMF_MAXSTR)        :: field_name
+    character (len=ESMF_MAXSTR)        :: emission_scheme
     integer                            :: i
     real                               :: DEFVAL
     logical                            :: data_driven = .true.
@@ -137,6 +150,8 @@ subroutine SetServices (GC, RC)
     call ESMF_ConfigGetAttribute (cfg, self%Ch_DU_res,  label='Ch_DU:', __RC__)
     call ESMF_ConfigGetAttribute (cfg, self%rlow,       label='radius_lower:', __RC__)
     call ESMF_ConfigGetAttribute (cfg, self%rup,        label='radius_upper:', __RC__)
+    call ESMF_ConfigGetAttribute (cfg, emission_scheme, label='emission_scheme:', default='ginoux', __RC__)
+    self%emission_scheme = ESMF_UtilStringLowerCase(trim(emission_scheme), __RC__)
     call ESMF_ConfigGetAttribute (cfg, self%point_emissions_srcfilen, &
                                   label='point_emissions_srcfilen:', default='/dev/null', __RC__)
     if ( (index(self%point_emissions_srcfilen,'/dev/null')>0) ) then
@@ -145,6 +160,14 @@ subroutine SetServices (GC, RC)
        self%doing_point_emissions = .true.  ! we are good to go
     end if
 
+!   read FENGSHA-specific parameters
+!   --------------------------------
+    if (self%emission_scheme == 'fengsha') then
+      call ESMF_ConfigGetAttribute (cfg, self%alpha,  label='alpha:', __RC__)
+      call ESMF_ConfigGetAttribute (cfg, self%gamma,  label='gamma:', __RC__)
+      call ESMF_ConfigGetAttribute (cfg, self%kvhmax, label='vertical_to_horizontal_flux_ratio_limit:', __RC__)
+    end if
+
 !   Is DU data driven?
 !   ------------------
     call determine_data_driven (COMP_NAME, data_driven, __RC__)
@@ -247,8 +270,10 @@ subroutine SetServices (GC, RC)
 !   ----------------------
     if (.not. data_driven) then
 #include "DU2G_Export___.h"
-#include "DU2G_Import___.h"
 #include "DU2G_Internal___.h"
+      associate (scheme => self%emission_scheme)
+#include "DU2G_Import___.h"
+      end associate
     end if
 
 !   This state holds fields needed by radiation
@@ -371,6 +396,13 @@ subroutine Initialize (GC, import, export, clock, RC)
     self%Ch_DU = Chem_UtilResVal(dims(1), dims(2), self%Ch_DU_res(:), __RC__)
     self%Ch_DU = self%Ch_DU * 1.00E-09
 
+!   Dust emission size distribution for FENGSHA
+!   ---------------------------------------------------------------
+    if (self%emission_scheme == 'fengsha') then
+      allocate(self%sdist(self%nbins), __STAT__)
+      call DustAerosolDistributionKok(self%radius, self%rup, self%rlow, self%sdist)
+    end if
+
 !   Get dimensions
 !   ---------------
     km = dims(3)
@@ -676,16 +708,6 @@ subroutine Run1 (GC, import, export, clock, RC)
 !   -----------------------------------
     call MAPL_Get (mapl, INTERNAL_ESMF_STATE=internal, __RC__)
 
-#include "DU2G_GetPointer___.h"
-
-do n=1,5
-   if(mapl_am_i_root()) print*,'n = ', n,' : Run1 B DU2G sum(du00n) = ',sum(DU(:,:,:,n))
-end do
-
-!   Set du_src to 0 where undefined
-!   --------------------------------
-    where (1.01*du_src > MAPL_UNDEF) du_src = 0.
-
 !   Get my private internal state
 !   ------------------------------
     call ESMF_UserCompGetInternalState(GC, 'DU2G_GridComp', wrap, STATUS)
@@ -699,6 +721,14 @@ subroutine Run1 (GC, import, export, clock, RC)
     call MAPL_PackTime (nymd, iyr, imm , idd)
     call MAPL_PackTime (nhms, ihr, imn, isc)
 
+    associate (scheme => self%emission_scheme)
+#include "DU2G_GetPointer___.h"
+    end associate
+
+do n=1,5
+   if(mapl_am_i_root()) print*,'n = ', n,' : Run1 B DU2G sum(du00n) = ',sum(DU(:,:,:,n))
+end do
+
 !   Get dimensions
 !   ---------------
     import_shape = shape(wet1)
@@ -714,9 +744,23 @@ subroutine Run1 (GC, import, export, clock, RC)
 
 !   Get surface gridded emissions
 !   -----------------------------
-    call DustEmissionGOCART2G(self%radius*1.e-6, frlake, wet1, lwi, u10m, v10m, &
-                              self%Ch_DU, du_src, MAPL_GRAV, &
-                              emissions_surface, __RC__)
+    select case (self%emission_scheme)
+      case ('fengsha')
+        call DustEmissionFENGSHA (frlake, frsnow, lwi, slc, du_clay, du_sand, du_silt,       &
+                                  du_ssm, du_rdrag, airdens(:,:,self%km), ustar, du_uthres,  &
+                                  self%alpha, self%gamma, self%kvhmax, MAPL_GRAV,   &
+                                  self%rhop, self%sdist, emissions_surface, __RC__)
+      case ('ginoux')
+!       Set du_src to 0 where undefined
+!       --------------------------------
+        where (1.01*du_src > MAPL_UNDEF) du_src = 0.
+
+        call DustEmissionGOCART2G(self%radius*1.e-6, frlake, wet1, lwi, u10m, v10m, &
+                                  self%Ch_DU, du_src, MAPL_GRAV, &
+                                  emissions_surface, __RC__)
+      case default
+        _ASSERT_RC(.false.,'missing dust emission scheme',ESMF_RC_NOT_IMPL)
+    end select
 
 !   Read point emissions file once per day
 !   --------------------------------------
@@ -828,18 +872,20 @@ subroutine Run2 (GC, import, export, clock, RC)
 !   -----------------------------------
     call MAPL_Get (MAPL, INTERNAL_ESMF_STATE=internal, __RC__)
 
-#include "DU2G_GetPointer___.h"
-
-do n=1,5
-   if(mapl_am_i_root()) print*,'n = ', n,' : Run2 B DU2G sum(du00n) = ',sum(DU(:,:,:,n))
-end do
-
 !   Get my private internal state
 !   ------------------------------
     call ESMF_UserCompGetInternalState(GC, 'DU2G_GridComp', wrap, STATUS)
     VERIFY_(STATUS)
     self => wrap%ptr
 
+    associate (scheme => self%emission_scheme)
+#include "DU2G_GetPointer___.h"
+    end associate
+
+do n=1,5
+   if(mapl_am_i_root()) print*,'n = ', n,' : Run2 B DU2G sum(du00n) = ',sum(DU(:,:,:,n))
+end do
+
     allocate(dqa, mold=wet1, __STAT__)
     allocate(drydepositionfrequency, mold=wet1, __STAT__)
 

ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_GridComp_DU.rc

@@ -37,5 +37,13 @@ nbins: 5
 
 pressure_lid_in_hPa: 0.01
 
+# Emissions methods
+emission_scheme: Ginoux       # available schemes: Ginoux, FENGSHA
+
+# FENGSHA settings
+alpha: 0.3
+gamma: 1.3
+vertical_to_horizontal_flux_ratio_limit: 2.e-04
+
 #point_emissions_srcfilen: /gpfsm/dnb32/esherman/gocartRefactor/RC/CA2G_point_src_test.rc
 

ESMF/GOCART2G_GridComp/DU2G_GridComp/DU2G_StateSpecs.rc

@@ -7,10 +7,18 @@ category: IMPORT
 #----------------------------------------------------------------------------------------
      NAME   | UNITS  | DIMS | VLOC| COND | LONG NAME
 #----------------------------------------------------------------------------------------
- DU_SRC     | 1      | xy   | N   |      | erod - dust emissions
- FRLAKE     | 1      | xy   | N   |      | fraction_of_lake
+ DU_SRC     | 1      | xy   | N   | scheme == 'ginoux'  | erod - dust emissions
+ DU_CLAY    | 1      | xy   | N   | scheme == 'fengsha' | volume_fraction_of_clay_in_soil
+ DU_SAND    | 1      | xy   | N   | scheme == 'fengsha' | volume_fraction_of_sand_in_soil
+ DU_SILT    | 1      | xy   | N   | scheme == 'fengsha' | volume_fraction_of_silt_in_soil
+ DU_RDRAG   | m-1    | xy   | N   | scheme == 'fengsha' | drag_partition
+ DU_SSM     | 1      | xy   | N   | scheme == 'fengsha' | sediment_supply_map
+ DU_UTHRES  | m s-1  | xy   | N   | scheme == 'fengsha' | surface_dry_threshold_velocity
+ FRSNOW     | 1      | xy   | N   | scheme == 'fengsha' | surface_snow_area_fraction
+ SLC        | 1      | xy   | N   | scheme == 'fengsha' | liquid_water_content_of_soil_layer
  WET1       | 1      | xy   | N   |      | surface_soil_wetness
  LWI        | 1      | xy   | N   |      | land-ocean-ice_mask
+ FRLAKE     | 1      | xy   | N   |      | fraction_of_lake
  U10M       | m s-1  | xy   | N   |      | 10-meter_eastward_wind
  V10M       | m s-1  | xy   | N   |      | 10-meter_northward_wind
  AREA       | m^2    | xy   | N   |      | agrid_cell_area

ESMF/UFS/Aerosol_Cap.F90

@@ -24,9 +24,8 @@ module Aerosol_Cap
   implicit none
 
   ! -- import fields
-  integer, parameter :: importFieldCount = 26
-  character(len=*), dimension(importFieldCount), parameter :: &
-    importFieldNames = (/ &
+  character(len=*), dimension(*), parameter :: &
+    importFieldNames = [ &
       "inst_pres_interface                  ", &
       "inst_pres_levels                     ", &
       "inst_geop_interface                  ", &
@@ -50,18 +49,19 @@ module Aerosol_Cap
       "inst_sensi_heat_flx                  ", &
       "inst_surface_roughness               ", &
       "inst_surface_soil_wetness            ", &
+      "inst_soil_moisture_content           ", &
       "ice_fraction                         ", &
       "lake_fraction                        ", &
-      "ocean_fraction                       "  &
-    /)
+      "ocean_fraction                       ", &
+      "surface_snow_area_fraction           "  &
+    ]
   ! -- export fields
-  integer, parameter :: exportFieldCount = 3
-  character(len=*), dimension(exportFieldCount), parameter :: &
-    exportFieldNames = (/ &
+  character(len=*), dimension(*), parameter :: &
+    exportFieldNames = [ &
       "inst_tracer_mass_frac                ", &
       "inst_tracer_up_surface_flx           ", &
       "inst_tracer_down_surface_flx         "  &
-    /)
+    ]
 
   type Aerosol_InternalState_T
     type(MAPL_Cap), pointer :: maplCap
@@ -187,28 +187,24 @@ subroutine ModelInitializeP1(model, importState, exportState, clock, rc)
     rc = ESMF_SUCCESS
 
     ! -- advertise imported fields
-    if (importFieldCount > 0) then
-      call NUOPC_Advertise(importState, importFieldNames, &
-        TransferOfferGeomObject="cannot provide", &
-        SharePolicyField="share", &
-        rc=rc)
-      if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
-        line=__LINE__, &
-        file=__FILE__)) &
-        return  ! bail out
-    end if
+    call NUOPC_Advertise(importState, importFieldNames, &
+      TransferOfferGeomObject="cannot provide", &
+      SharePolicyField="share", &
+      rc=rc)
+    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+      line=__LINE__, &
+      file=__FILE__)) &
+      return  ! bail out
 
     ! -- advertise exported fields
-    if (exportFieldCount > 0) then
-      call NUOPC_Advertise(exportState, exportFieldNames, &
-        TransferOfferGeomObject="cannot provide", &
-        SharePolicyField="share", &
-        rc=rc)
-      if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
-        line=__LINE__, &
-        file=__FILE__)) &
-        return  ! bail out
-    end if
+    call NUOPC_Advertise(exportState, exportFieldNames, &
+      TransferOfferGeomObject="cannot provide", &
+      SharePolicyField="share", &
+      rc=rc)
+    if (ESMF_LogFoundError(rcToCheck=rc, msg=ESMF_LOGERR_PASSTHRU, &
+      line=__LINE__, &
+      file=__FILE__)) &
+      return  ! bail out
 
   end subroutine ModelInitializeP1
 

ESMF/UFS/Aerosol_Comp_Mod.F90

@@ -6,12 +6,12 @@ module Aerosol_Comp_Mod
 
   implicit none
 
-  integer, parameter :: fieldMapSize = 19
-  character(len=*), dimension(fieldMapSize, 2), parameter :: &
-    fieldMap = reshape((/ &
+  character(len=*), dimension(*), parameter :: &
+    fieldPairList = [ &
       "FROCEAN                         ", "ocean_fraction                  ", &
       "FRACI                           ", "ice_fraction                    ", &
       "FRLAKE                          ", "lake_fraction                   ", &
+      "FRSNOW                          ", "surface_snow_area_fraction      ", &
       "LWI                             ", "inst_land_sea_mask              ", &
       "WET1                            ", "inst_surface_soil_wetness       ", &
       "U10M                            ", "inst_zonal_wind_height10m       ", &
@@ -28,7 +28,10 @@ module Aerosol_Comp_Mod
       "PFI_LSAN                        ", "inst_ice_nonconv_tendency_levels", &
       "PFL_LSAN                        ", "inst_liq_nonconv_tendency_levels", &
       "FCLD                            ", "inst_cloud_frac_levels          "  &
-      /), (/fieldMapSize, 2/), order=(/2,1/))
+    ]
+
+  character(len=*), dimension(*,2), parameter :: &
+    fieldMap = reshape(fieldPairList, [size(fieldPairList)/2,2], order=[2,1])
 
   ! gravity (m/s2)
   real(ESMF_KIND_R8), parameter :: con_g = 9.80665e+0_ESMF_KIND_R8
@@ -98,6 +101,7 @@ subroutine AerosolStateUpdate(model, cap, action, rc)
     integer :: localrc, stat
     integer :: imap, item, itemCount
     integer :: rank
+    integer :: fieldMapSize
     integer :: i,j, k, kk, n, ni, nj, nk, nlev, offset, v
     integer, dimension(1) :: plb, pub, rlb, rub
     real(ESMF_KIND_R4) :: blkevap, blkesat
@@ -107,7 +111,7 @@ subroutine AerosolStateUpdate(model, cap, action, rc)
     real(ESMF_KIND_R4), dimension(:,:,:,:), pointer :: fp4dr
     real(ESMF_KIND_R8), dimension(:,:),     pointer :: fp2dp
     real(ESMF_KIND_R8), dimension(:,:),     pointer :: rain, rainc, zorl
-    real(ESMF_KIND_R8), dimension(:,:,:),   pointer :: fp3dp
+    real(ESMF_KIND_R8), dimension(:,:,:),   pointer :: fp3dp, slc
     real(ESMF_KIND_R8), dimension(:,:,:),   pointer :: phii, phil, prsi, prsl, temp
     real(ESMF_KIND_R8), dimension(:,:,:,:), pointer :: q
     type(ESMF_Clock) :: clock
@@ -226,6 +230,12 @@ subroutine AerosolStateUpdate(model, cap, action, rc)
             file=__FILE__,  &
             rcToReturn=rc)) return  ! bail out
 
+          call AerosolGetPtr(state, "inst_soil_moisture_content", slc, rc=localrc)
+          if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
+            line=__LINE__,  &
+            file=__FILE__,  &
+            rcToReturn=rc)) return  ! bail out
+
           call AerosolGetPtr(state, "inst_tracer_mass_frac",  q, rc=localrc)
           if (ESMF_LogFoundError(rcToCheck=localrc, msg=ESMF_LOGERR_PASSTHRU, &
             line=__LINE__,  &
@@ -235,6 +245,7 @@ subroutine AerosolStateUpdate(model, cap, action, rc)
           ni = size(q, 1)
           nj = size(q, 2)
           nk = size(q, 3)
+          fieldMapSize = size(fieldMap, 1)
 
           do item = 1, itemCount
             if (itemTypeList(item) == ESMF_STATEITEM_FIELD) then
@@ -371,6 +382,8 @@ subroutine AerosolStateUpdate(model, cap, action, rc)
                       end do
                     end do
                   end do
+                case ("SLC")
+                  fp2dr = slc(:,:,1)
                 case ("Z0H")
                   fp2dr = 0.01_ESMF_KIND_R4 * zorl
                 case ("ZLE")