Drydeposition

.gitmodules

@@ -27,8 +27,8 @@
 #
 [submodule "fates"]
 path = src/fates
-url = https://github.com/NGEET/fates
-fxtag = sci.1.78.2_api.36.0.0
+url = https://github.com/rosiealice/fates
+fxtag = drydep_main
 fxrequired = AlwaysRequired
 # Standard Fork to compare to with "git fleximod test" to ensure personal forks aren't committed
 fxDONOTUSEurl = https://github.com/NCAR/fates-release

src/biogeochem/DryDepVelocity.F90

@@ -208,6 +208,7 @@ subroutine depvel_compute( bounds, &
     use pftconMod      , only : nc3_nonarctic_grass, nc4_grass, nc3crop
     use pftconMod      , only : nc3irrig, npcropmin, npcropmax
     use clm_varcon     , only : spval
+    use clm_varctl     , only : use_fates
 
     !
     ! !ARGUMENTS:
@@ -358,7 +359,20 @@ subroutine depvel_compute( bounds, &
                     msg='ERROR: Not able to determine Wesley vegetation type'//&
                     errMsg(sourcefile, __LINE__))
             end if
-
+
+             if(use_fates)then
+                if(patch%is_fates(pi))then
+                   wesveg = canopystate_inst%wesley_veg_index_patch(pi)
+                else
+                   wesveg = 8 !make bare ground for non-fates patches. Some of these are overwritten below.
+                endif
+             endif   
+
+             if(wesveg<0 .or. wesveg>1 )then
+                call endrun(subgrid_index=pi, subgrid_level=subgrid_level_patch, &
+                    msg='ERROR: No sensible Wesley vegetation type'//&
+                    errMsg(sourcefile, __LINE__))
+             endif
             ! create seasonality index used to index wesely data tables from LAI,  Bascially
             !if elai is between max lai from input data and half that max the index_season=1
 
@@ -376,12 +390,13 @@ subroutine depvel_compute( bounds, &
             ! 4 - Winter, snow on ground and subfreezing
             ! 5 - Transitional spring with partially green short annuals
 
-
-            !mlaidiff=jan-feb
-            minlai=minval(annlai(:,pi))
-            maxlai=maxval(annlai(:,pi))
-
-            index_season = -1
+            if(.not. use_fates)then !for non-FATES runs we use satellite phenology to choose the season for the ressistance parameters. 
+               !mlaidiff=jan-feb
+               minlai=minval(annlai(:,pi))
+               maxlai=maxval(annlai(:,pi))
+            endif
+
+               index_season = -1
 
             if ( lun%itype(l) /= istsoil )then
                if ( lun%itype(l) == istice ) then
@@ -399,27 +414,38 @@ subroutine depvel_compute( bounds, &
                end if
             else if ( snow_depth(c) > 0 ) then
                index_season = 4
-            else if(elai(pi) > 0.5_r8*maxlai) then
+            else if(.not.use_fates .and. elai(pi) > 0.5_r8*maxlai) then
                index_season = 1
             endif
-
-            if (index_season<0) then
-               if (elai(pi) < (minlai+0.05_r8*(maxlai-minlai))) then
-                  index_season = 3
+
+            if(use_fates.and.index_season<1)then 
+               if(patch%is_fates(pi))then
+                  index_season = canopystate_inst%wesley_season_index_patch(pi)
+               else
+                  index_season = 2 !set intermediate spring seson for bare ground. (as for urban)
                endif
-            endif
+
+            else ! not fates
 
-            if (index_season<0) then
-               if (mlaidiff(pi) > 0.0_r8) then
-                  index_season = 2
-               elseif (mlaidiff(pi) < 0.0_r8) then
-                  index_season = 5
-               elseif (mlaidiff(pi).eq.0.0_r8) then
-                  index_season = 3
+               if (index_season<0) then
+                  if (elai(pi) < (minlai+0.05_r8*(maxlai-minlai))) then
+                     index_season = 3
+                  endif
                endif
-            endif
 
-            if (index_season<0) then
+               if (index_season<0) then
+                  if (mlaidiff(pi) > 0.0_r8) then
+                     index_season = 2
+                  elseif (mlaidiff(pi) < 0.0_r8) then
+                     index_season = 5
+                  elseif (mlaidiff(pi).eq.0.0_r8) then
+                     index_season = 3
+                  endif
+               endif
+
+            endif ! use_fates
+
+            if (index_season<0.or.index_season>5) then
                call endrun('ERROR: not able to determine season'//errmsg(sourcefile, __LINE__))
             endif
 
@@ -472,14 +498,14 @@ subroutine depvel_compute( bounds, &
 
                ! correction for frost
                cts = 1000._r8*exp( -tc - 4._r8 )
-
                !ground resistance
                rgsx(ispec) = 1._r8/((heff(ispec)/(1.e5_r8*(rgss(index_season,wesveg)+cts))) + &
                     (foxd(ispec)/(rgso(index_season,wesveg)+cts)))
 
                !-------------------------------------------------------------------------------------
                ! special case for H2 and CO;; CH4 is set ot a fraction of dv(H2)
                !-------------------------------------------------------------------------------------
+
                if( ispec == index_h2 .or. ispec == index_co .or. ispec == index_ch4 ) then
 
                   if( ispec == index_co ) then

src/biogeophys/CanopyStateType.F90

@@ -52,6 +52,9 @@ module CanopyStateType
 
      real(r8) , pointer :: dleaf_patch              (:)   ! patch characteristic leaf width (diameter) [m]
                                                           ! for non-ED/FATES this is the same as pftcon%dleaf()
+     integer ,  pointer :: wesley_veg_index_patch   (:)   ! Wesley PFT index for FATES dry deposition calculations
+     integer ,  pointer :: wesley_season_index_patch (:)   ! Season for dry deposition calculations
+
      real(r8) , pointer :: rscanopy_patch           (:)   ! patch canopy stomatal resistance (s/m) (ED specific)
 
      real(r8) , pointer :: vegwp_patch              (:,:) ! patch vegetation water matric potential (mm)
@@ -144,6 +147,8 @@ subroutine InitAllocate(this, bounds)
     allocate(this%fsun240_patch            (begp:endp))           ; this%fsun240_patch            (:)   = nan
 
     allocate(this%dleaf_patch              (begp:endp))           ; this%dleaf_patch              (:)   = nan
+    allocate(this%wesley_veg_index_patch   (begp:endp))           ; ! is an integer so can't be nan
+    allocate(this%wesley_season_index_patch   (begp:endp))           ; ! is an integer so can't be nan    
     allocate(this%rscanopy_patch           (begp:endp))           ; this%rscanopy_patch           (:)   = nan
 !    allocate(this%gccanopy_patch           (begp:endp))           ; this%gccanopy_patch           (:)   = 0.0_r8
     allocate(this%vegwp_patch              (begp:endp,1:nvegwcs)) ; this%vegwp_patch              (:,:) = nan

src/main/clm_initializeMod.F90

@@ -481,7 +481,7 @@ subroutine initialize2(ni,nj)
        ! Also for FATES with Dry-Deposition on as well (see above)
        ! For now don't allow for dry-deposition with full fates 
        ! because of issues in #1044 EBK Jun/17/2022
-       if( use_fates_nocomp .or. (.not. use_fates )) then
+       if( use_fates_sp .or. (.not. use_fates )) then
           if (masterproc) then
              write(iulog,'(a)')'Initializing Satellite Phenology'
           end if
@@ -698,12 +698,13 @@ subroutine initialize2(ni,nj)
     ! Read monthly vegetation
     ! Even if CN or FATES is on, and dry-deposition is active, read CLMSP annual vegetation
     ! to get estimates of monthly LAI
-    if ( n_drydep > 0 ) then
+    if ( n_drydep > 0 .and. .not. use_fates ) then
        call readAnnualVegetation(bounds_proc, canopystate_inst)
        ! Call interpMonthlyVeg for dry-deposition so that mlaidiff will be calculated
        ! This needs to be done even if FATES, CN or CNDV is on!
        call interpMonthlyVeg(bounds_proc, canopystate_inst)
-    elseif ( use_fates_sp ) then
+    endif 
+    if ( use_fates_sp ) then
        ! If fates has satellite phenology enabled, get the monthly veg values
        ! prior to the first call to SatellitePhenology()
        call interpMonthlyVeg(bounds_proc, canopystate_inst)

src/utils/clmfates_interfaceMod.F90

@@ -1485,6 +1485,8 @@ subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump,      &
          z0m  => canopystate_inst%z0m_patch  , & ! Output: [real(r8) (:)   ] momentum roughness length (m)
          displa => canopystate_inst%displa_patch, &
          dleaf_patch => canopystate_inst%dleaf_patch, &
+         wesley_veg_index_patch => canopystate_inst%wesley_veg_index_patch, &
+         wesley_season_index_patch => canopystate_inst%wesley_season_index_patch, &
          snow_depth => waterdiagnosticbulk_inst%snow_depth_col, &
          frac_sno_eff => waterdiagnosticbulk_inst%frac_sno_eff_col, &
          frac_veg_nosno_alb => canopystate_inst%frac_veg_nosno_alb_patch)
@@ -1608,7 +1610,8 @@ subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump,      &
           z0m(col%patchi(c)+1:col%patchf(c)) = 0.0_r8
           displa(col%patchi(c)+1:col%patchf(c)) = 0.0_r8
           dleaf_patch(col%patchi(c)+1:col%patchf(c)) = 0.0_r8
-
+          wesley_veg_index_patch(col%patchi(c)+1:col%patchf(c)) = 0
+          wesley_season_index_patch(col%patchi(c)+1:col%patchf(c)) = 0
           frac_veg_nosno_alb(col%patchi(c):col%patchf(c)) = 0.0_r8
 
           ! Set the bareground patch indicator
@@ -1670,6 +1673,8 @@ subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump,      &
              z0m(p)    = this%fates(nc)%bc_out(s)%z0m_pa(ifp)
              displa(p) = this%fates(nc)%bc_out(s)%displa_pa(ifp)
              dleaf_patch(p) = this%fates(nc)%bc_out(s)%dleaf_pa(ifp)
+             wesley_veg_index_patch(p) = this%fates(nc)%bc_out(s)%wesley_pft_label_pa(ifp)
+             wesley_season_index_patch(p) = this%fates(nc)%bc_out(s)%drydep_season_pa(ifp)
           end do ! veg pach
 
           if(abs(areacheck - 1.0_r8).gt.1.e-9_r8)then