Fates two-stream restart fixes

cime_config/testdefs/ExpectedTestFails.xml

@@ -308,20 +308,6 @@
     </phase>
   </test>
 
-  <test name="ERS_D_Ld20.f45_f45_mg37.I2000Clm50FatesRs.derecho_intel.clm-FatesColdTwoStream">
-    <phase name="COMPARE_base_rest">
-      <status>FAIL</status>
-      <issue>#2325</issue>
-    </phase>
-  </test>
-
-  <test name="ERS_D_Ld20.f45_f45_mg37.I2000Clm50FatesRs.derecho_gnu.clm-FatesColdTwoStreamNoCompFixedBioGeo">
-    <phase name="COMPARE_base_rest">
-      <status>FAIL</status>
-      <issue>#2325</issue>
-    </phase>
-  </test>
-
   <test name="SMS_Ld10_D_Mmpi-serial.CLM_USRDAT.I1PtClm60Fates.derecho_intel.clm-FatesFireLightningPopDens--clm-NEON-FATES-NIWO">
     <phase name="SHAREDLIB_BUILD">
       <status>FAIL</status>

src/biogeophys/BalanceCheckMod.F90

@@ -956,6 +956,7 @@ subroutine BalanceCheck( bounds, &
 
            indexp = maxloc( abs(errsol(bounds%begp:bounds%endp)), 1 , mask = (errsol(bounds%begp:bounds%endp) /= spval) ) + bounds%begp -1
            indexg = patch%gridcell(indexp)
+           indexc = patch%column(indexp)
            write(iulog,*)'WARNING:: BalanceCheck, solar radiation balance error (W/m2)'
            write(iulog,*)'nstep         = ',nstep
            write(iulog,*)'errsol        = ',errsol(indexp)
@@ -969,7 +970,18 @@ subroutine BalanceCheck( bounds, &
                write(iulog,*)'forc_solai(1) = ',forc_solai(indexg,1)
                write(iulog,*)'forc_solai(2) = ',forc_solai(indexg,2)
                write(iulog,*)'forc_tot      = ',forc_solad(indexg,1)+forc_solad(indexg,2) &
-                  +forc_solai(indexg,1)+forc_solai(indexg,2)
+                    +forc_solai(indexg,1)+forc_solai(indexg,2)
+
+               write(iulog,*)'coszen_col:',surfalb_inst%coszen_col(indexc)
+               write(iulog,*)'fabd:',fabd(indexp,:)
+               write(iulog,*)'fabi:',fabi(indexp,:)
+               write(iulog,*)'albd:',albd(indexp,:)
+               write(iulog,*)'albi:',albi(indexp,:)
+               write(iulog,*)'ftdd:',ftdd(indexp,:)
+               write(iulog,*)'ftid:',ftid(indexp,:)
+               write(iulog,*)'ftii:',ftii(indexp,:)
+
+
                write(iulog,*)'CTSM is stopping'
                call endrun(subgrid_index=indexp, subgrid_level=subgrid_level_patch, msg=errmsg(sourcefile, __LINE__))
            end if

src/biogeophys/SurfaceAlbedoMod.F90

@@ -35,6 +35,8 @@ module SurfaceAlbedoMod
   public :: SurfaceAlbedo_readnl
   public :: SurfaceAlbedoInitTimeConst
   public :: SurfaceAlbedo  ! Surface albedo and two-stream fluxes
+  public :: UpdateZenithAngles
+
   !
   ! !PRIVATE MEMBER FUNCTIONS:
   private :: SoilAlbedo    ! Determine ground surface albedo
@@ -224,7 +226,57 @@ subroutine SurfaceAlbedoInitTimeConst(bounds)
 
   end subroutine SurfaceAlbedoInitTimeConst
 
+  ! -----------------------------------------------------------------------
+
+  subroutine UpdateZenithAngles(bounds, surfalb_inst, nextsw_cday, declinp1)
+
+    ! Incorporate surface slopes to generate column level zenith angles
+
+    use clm_varctl          , only : downscale_hillslope_meteorology
+    use shr_orb_mod
+
+    type(bounds_type)      , intent(in)            :: bounds             ! bounds
+    type(surfalb_type)     , intent(inout)         :: surfalb_inst
+    real(r8)               , intent(in)            :: nextsw_cday        ! calendar day at Greenwich (1.00, ..., days/year)
+    real(r8)               , intent(in)            :: declinp1           ! declination angle (radians) for next time step
+
+    integer :: c,g        ! indices
+
+    associate( &
+         azsun_grc     =>    surfalb_inst%azsun_grc              , & ! Output:  [real(r8) (:)   ]  cosine of solar zenith angle            
+         coszen_grc    =>    surfalb_inst%coszen_grc             , & ! Output:  [real(r8) (:)   ]  cosine of solar zenith angle            
+         coszen_col    =>    surfalb_inst%coszen_col)                ! Output:  [real(r8) (:)   ]  cosine of solar zenith angle
+
+
+      ! Cosine solar zenith angle for next time step
+      ! First calculate grid-scale zenith based on time and location, no topographic effects
+      do g = bounds%begg,bounds%endg
+         coszen_grc(g) = shr_orb_cosz (nextsw_cday, grc%lat(g), grc%lon(g), declinp1)
+         azsun_grc(g) = shr_orb_azimuth(nextsw_cday, grc%lat(g), grc%lon(g), declinp1, acos(coszen_grc(g)))
+      end do
+
+      ! calculate local incidence angle based on column slope and aspect
+      do c = bounds%begc,bounds%endc
+         g = col%gridcell(c)
+         if (col%is_hillslope_column(c) .and. downscale_hillslope_meteorology) then
+
+            ! hill_slope is [m/m], convert to radians with atan function
+            coszen_col(c) = shr_orb_cosinc(acos(coszen_grc(g)),azsun_grc(g),atan(col%hill_slope(c)),col%hill_aspect(c))
+
+            if(coszen_grc(g) > 0._r8 .and. coszen_col(c) < 0._r8) coszen_col(c) = 0._r8
+         else
+            coszen_col(c) = coszen_grc(g)
+         endif
+
+      end do
+
+    end associate
+
+    return
+  end subroutine UpdateZenithAngles
+
   !-----------------------------------------------------------------------
+
   subroutine SurfaceAlbedo(bounds,nc,  &
         num_nourbanc, filter_nourbanc, &
         num_nourbanp, filter_nourbanp, &
@@ -256,13 +308,13 @@ subroutine SurfaceAlbedo(bounds,nc,  &
     ! only computed over active points.
     !
     ! !USES:
-    use shr_orb_mod
+    !use shr_orb_mod
     use clm_time_manager   , only : get_nstep
     use abortutils         , only : endrun
     use clm_varctl         , only : use_subgrid_fluxes, use_snicar_frc, use_fates
     use CLMFatesInterfaceMod, only : hlm_fates_interface_type
     use landunit_varcon     , only : istsoil
-    use clm_varctl          , only : downscale_hillslope_meteorology
+
 
     ! !ARGUMENTS:
     type(bounds_type)      , intent(in)            :: bounds             ! bounds
@@ -427,32 +479,15 @@ subroutine SurfaceAlbedo(bounds,nc,  &
           fabi_sha_z    =>    surfalb_inst%fabi_sha_z_patch         & ! Output:  [real(r8) (:,:) ]  absorbed shaded leaf diffuse PAR (per unit lai+sai) for each canopy layer
           )
 
-    ! Cosine solar zenith angle for next time step
-
-    do g = bounds%begg,bounds%endg
-       coszen_grc(g) = shr_orb_cosz (nextsw_cday, grc%lat(g), grc%lon(g), declinp1)
-    end do
-
-    do c = bounds%begc,bounds%endc
-       g = col%gridcell(c)
-       if (col%is_hillslope_column(c) .and. downscale_hillslope_meteorology) then
-          ! calculate local incidence angle based on column slope and aspect
-          zenith_angle = acos(coszen_grc(g))
-
-          azsun_grc(g) = shr_orb_azimuth(nextsw_cday, grc%lat(g), grc%lon(g), declinp1, zenith_angle)
-          ! hill_slope is [m/m], convert to radians
-          coszen_col(c) = shr_orb_cosinc(zenith_angle,azsun_grc(g),atan(col%hill_slope(c)),col%hill_aspect(c))
-
-          if(coszen_grc(g) > 0._r8 .and. coszen_col(c) < 0._r8) coszen_col(c) = 0._r8
-
-       else
-          coszen_col(c) = coszen_grc(g)
-       endif
-    end do
+    ! Apply column level zenith angles to the patch level
     do fp = 1,num_nourbanp
        p = filter_nourbanp(fp)
        c = patch%column(p)
-       coszen_patch(p) = coszen_col(c)
+       if(patch%is_fates(p))then
+          coszen_patch(p) = spval
+       else
+          coszen_patch(p) = coszen_col(c)
+       end if
     end do
 
     ! Initialize output because solar radiation only done if coszen > 0
@@ -1045,10 +1080,7 @@ subroutine SurfaceAlbedo(bounds,nc,  &
 
     if (use_fates) then
 
-       call clm_fates%wrap_canopy_radiation(bounds, nc, &
-            num_vegsol, filter_vegsol, &
-            coszen_patch(bounds%begp:bounds%endp), &
-            fcansno(bounds%begp:bounds%endp), surfalb_inst)
+       call clm_fates%wrap_canopy_radiation(bounds, nc, fcansno(bounds%begp:bounds%endp), surfalb_inst)
 
     else
 

src/main/clm_driver.F90

@@ -52,6 +52,7 @@ module clm_driver
   use AerosolMod             , only : AerosolMasses
   use SnowSnicarMod          , only : SnowAge_grain
   use SurfaceAlbedoMod       , only : SurfaceAlbedo
+  use SurfaceAlbedoMod       , only : UpdateZenithAngles
   use UrbanAlbedoMod         , only : UrbanAlbedo
   !
   use SurfaceRadiationMod    , only : SurfaceRadiation, CanopySunShadeFracs
@@ -1224,7 +1225,15 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro
        ! Determine albedos for next time step
        ! ============================================================================
 
-       if (doalb) then
+       if (doalb .or. use_fates) then
+          ! During branch runs and non continue_run restarts, the doalb flag
+          ! does not trigger correctly for fates runs (and non-fates?), and thus
+          ! the zenith angles are not calculated and ready when radiation scattering
+          ! needs to occur.
+          call UpdateZenithAngles(bounds_clump,surfalb_inst, nextsw_cday, declinp1)
+       end if
+
+       if (doalb ) then
 
           ! Albedos for non-urban columns
           call t_startf('surfalb')

src/utils/clmfates_interfaceMod.F90

@@ -2752,28 +2752,31 @@ end subroutine wrap_WoodProducts
 
  ! ======================================================================================
 
- subroutine wrap_canopy_radiation(this, bounds_clump, nc, &
-         num_vegsol, filter_vegsol, coszen, fcansno,  surfalb_inst)
+ subroutine wrap_canopy_radiation(this, bounds_clump, nc, fcansno, surfalb_inst)
 
+
+    ! Pass boundary conditions (zenith angle, ground albedo and snow-cover)
+    ! to FATES, and have fates call normalized radiation scattering.
+    ! These methods are normalized because it is the zenith of the next
+    ! timestep, and because of atmospheric model coupling, we need to provide
+    ! an albedo. This normalized solution will be scaled by the actual
+    ! downwelling radiation during the wrap sunshade fraction call
 
+
     ! Arguments
     class(hlm_fates_interface_type), intent(inout) :: this
     type(bounds_type),  intent(in)             :: bounds_clump
-    ! filter for vegetated pfts with coszen>0
     integer            , intent(in)            :: nc ! clump index
-    integer            , intent(in)            :: num_vegsol
-    integer            , intent(in)            :: filter_vegsol(num_vegsol)
-    ! cosine solar zenith angle for next time step
-    real(r8)           , intent(in)            :: coszen( bounds_clump%begp: )
     real(r8)           , intent(in)            :: fcansno( bounds_clump%begp: )
     type(surfalb_type) , intent(inout)         :: surfalb_inst
 
     ! locals
-    integer                                    :: s,c,p,ifp,icp
+    integer                                    :: s,c,p,ifp,g
 
     call t_startf('fates_wrapcanopyradiation')
 
     associate(&
+         coszen_col   =>    surfalb_inst%coszen_col         , & !in
          albgrd_col   =>    surfalb_inst%albgrd_col         , & !in
          albgri_col   =>    surfalb_inst%albgri_col         , & !in
          albd         =>    surfalb_inst%albd_patch         , & !out
@@ -2787,56 +2790,61 @@ subroutine wrap_canopy_radiation(this, bounds_clump, nc, &
     do s = 1, this%fates(nc)%nsites
 
        c = this%f2hmap(nc)%fcolumn(s)
+       g = col%gridcell(c)
 
-         do ifp = 1,this%fates(nc)%sites(s)%youngest_patch%patchno
-           p = ifp+col%patchi(c)
+       this%fates(nc)%bc_in(s)%coszen = coszen_col(c)
 
-          if( any(filter_vegsol==p) )then
+       ! FATES only does short-wave radiation scattering when
+       ! the zenith angle is positive. In the future, FATES
+       ! may perform calculations on diffuse radiation during
+       ! dawn and dusk when the sun is below the horizon, but not yet
 
-             this%fates(nc)%bc_in(s)%filter_vegzen_pa(ifp) = .true.
-             this%fates(nc)%bc_in(s)%coszen_pa(ifp)  = coszen(p)
-             this%fates(nc)%bc_in(s)%fcansno_pa(ifp)  = fcansno(p)
-             this%fates(nc)%bc_in(s)%albgr_dir_rb(:) = albgrd_col(c,:)
-             this%fates(nc)%bc_in(s)%albgr_dif_rb(:) = albgri_col(c,:)
+       do ifp = 1,this%fates(nc)%sites(s)%youngest_patch%patchno
+          p = ifp+col%patchi(c)
+          this%fates(nc)%bc_in(s)%fcansno_pa(ifp) = fcansno(p)
+       end do
 
-          else
+
+       if(coszen_col(c) > 0._r8) then
 
-             this%fates(nc)%bc_in(s)%filter_vegzen_pa(ifp) = .false.
+          this%fates(nc)%bc_in(s)%albgr_dir_rb(:) = albgrd_col(c,:)
+          this%fates(nc)%bc_in(s)%albgr_dif_rb(:) = albgri_col(c,:)
+       else
 
-          end if
+          ! This will ensure a crash in FATES if it tries
+          this%fates(nc)%bc_in(s)%albgr_dir_rb(:) = spval
+          this%fates(nc)%bc_in(s)%albgr_dif_rb(:) = spval
+
+       end if
 
-       end do
     end do
 
-    call FatesNormalizedCanopyRadiation(this%fates(nc)%nsites,  &
-         this%fates(nc)%sites, &
-         this%fates(nc)%bc_in,  &
-         this%fates(nc)%bc_out)
+    call FatesNormalizedCanopyRadiation( &
+            this%fates(nc)%sites, &
+            this%fates(nc)%bc_in,  &
+            this%fates(nc)%bc_out)
 
     ! Pass FATES BC's back to HLM
     ! -----------------------------------------------------------------------------------
-    do icp = 1,num_vegsol
-       p = filter_vegsol(icp)
-       c = patch%column(p)
-       s = this%f2hmap(nc)%hsites(c)
-       ! do if structure here and only pass natveg columns
-       ifp = p-col%patchi(c)
+    do s = 1, this%fates(nc)%nsites
+
+       c = this%f2hmap(nc)%fcolumn(s)
+
+       do ifp = 1, this%fates(nc)%sites(s)%youngest_patch%patchno
+
+          p = ifp+col%patchi(c)
 
-       if(.not.this%fates(nc)%bc_in(s)%filter_vegzen_pa(ifp) )then
-          write(iulog,*) 's,p,ifp',s,p,ifp
-          write(iulog,*) 'Not all patches on the natveg column were passed to canrad',patch%sp_pftorder_index(p)
-!          call endrun(msg=errMsg(sourcefile, __LINE__))
-       else
           albd(p,:) = this%fates(nc)%bc_out(s)%albd_parb(ifp,:)
           albi(p,:) = this%fates(nc)%bc_out(s)%albi_parb(ifp,:)
           fabd(p,:) = this%fates(nc)%bc_out(s)%fabd_parb(ifp,:)
           fabi(p,:) = this%fates(nc)%bc_out(s)%fabi_parb(ifp,:)
           ftdd(p,:) = this%fates(nc)%bc_out(s)%ftdd_parb(ifp,:)
           ftid(p,:) = this%fates(nc)%bc_out(s)%ftid_parb(ifp,:)
           ftii(p,:) = this%fates(nc)%bc_out(s)%ftii_parb(ifp,:)
-       end if
+
+       end do
     end do
-
+          
   end associate
 
   call t_stopf('fates_wrapcanopyradiation')