Merge pull request #1321 from rgknox/twostream-restart-fixes-v2

restart fixes for two-stream

Author: rgknox

biogeochem/EDCanopyStructureMod.F90

@@ -1452,7 +1452,7 @@ subroutine canopy_summarization( nsites, sites, bc_in )
        call leaf_area_profile(sites(s))
 
        if(hlm_radiation_model.eq.twostr_solver) then
-          call FatesConstructRadElements(sites(s),bc_in(s)%fcansno_pa,bc_in(s)%coszen_pa)
+          call FatesConstructRadElements(sites(s))
        end if
 
     end do ! site loop
@@ -1858,9 +1858,10 @@ subroutine leaf_area_profile( currentSite )
        end if if_any_canopy_area
 
        cpatch => cpatch%younger
-
     enddo !patch
 
+    
+    
     return
   end subroutine leaf_area_profile
 

biogeochem/EDPatchDynamicsMod.F90

@@ -3212,6 +3212,8 @@ subroutine fuse_2_patches(csite, dp, rp)
     rp%zstar                = (dp%zstar*dp%area + rp%zstar*rp%area) * inv_sum_area
     rp%c_stomata            = (dp%c_stomata*dp%area + rp%c_stomata*rp%area) * inv_sum_area
     rp%c_lblayer            = (dp%c_lblayer*dp%area + rp%c_lblayer*rp%area) * inv_sum_area
+
+    ! Radiation
     rp%rad_error(1)         = (dp%rad_error(1)*dp%area + rp%rad_error(1)*rp%area) * inv_sum_area
     rp%rad_error(2)         = (dp%rad_error(2)*dp%area + rp%rad_error(2)*rp%area) * inv_sum_area
 

biogeochem/FatesPatchMod.F90

@@ -136,8 +136,6 @@ module FatesPatchMod
     ! RADIATION
     real(r8) :: rad_error(num_swb)                        ! radiation consv error by band [W/m2]
     real(r8) :: fcansno                                   ! fraction of canopy covered in snow [0-1]
-    logical  :: solar_zenith_flag                         ! integer flag specifying daylight (based on zenith angle)
-    real(r8) :: solar_zenith_angle                        ! solar zenith angle [radians]
     real(r8) :: gnd_alb_dif(num_swb)                      ! ground albedo for diffuse rad, both bands [0-1]
     real(r8) :: gnd_alb_dir(num_swb)                      ! ground albedo for direct rad, both bands [0-1]
 
@@ -466,8 +464,6 @@ subroutine NanValues(this)
       ! RADIATION
       this%rad_error(:)                 = nan
       this%fcansno                      = nan 
-      this%solar_zenith_flag            = .false. 
-      this%solar_zenith_angle           = nan 
       this%gnd_alb_dif(:)               = nan 
       this%gnd_alb_dir(:)               = nan
 
@@ -565,7 +561,6 @@ subroutine ZeroValues(this)
       this%c_stomata                         = 0.0_r8 
       this%c_lblayer                         = 0.0_r8
 
-      
       ! RADIATION
       this%rad_error(:)                      = 0.0_r8
       this%tr_soil_dir_dif(:)                = 0.0_r8
@@ -574,7 +569,7 @@ subroutine ZeroValues(this)
       this%fabd(:)                           = 0.0_r8
       this%sabs_dir(:)                       = 0.0_r8
       this%sabs_dif(:)                       = 0.0_r8
-
+      
       ! ROOTS
       this%btran_ft(:)                       = 0.0_r8
 
@@ -691,7 +686,7 @@ end subroutine InitLitter
     !===========================================================================
 
     subroutine Create(this, age, area, land_use_label, nocomp_pft, num_swb, num_pft,    &
-      num_levsoil, current_tod, regeneration_model) 
+      num_levsoil, current_tod, regeneration_model)
       !
       ! DESCRIPTION:
       ! create a new patch with input and default values
@@ -708,7 +703,7 @@ subroutine Create(this, age, area, land_use_label, nocomp_pft, num_swb, num_pft,
       integer,                 intent(in)    :: num_levsoil        ! number of soil layers
       integer,                 intent(in)    :: current_tod        ! time of day [seconds past 0Z]
       integer,                 intent(in)    :: regeneration_model ! regeneration model version
-    
+      
       ! initialize patch
       ! sets all values to nan, then some values to zero
       call this%Init(num_swb, num_levsoil)
@@ -1214,8 +1209,6 @@ subroutine Dump(this)
       write(fates_log(),*) 'pa%total_tree_area    = ',this%total_tree_area
       write(fates_log(),*) 'pa%total_grass_area   = ',this%total_grass_area
       write(fates_log(),*) 'pa%zstar              = ',this%zstar
-      write(fates_log(),*) 'pa%solar_zenith_flag  = ',this%solar_zenith_flag
-      write(fates_log(),*) 'pa%solar_zenith_angle = ',this%solar_zenith_angle
       write(fates_log(),*) 'pa%gnd_alb_dif        = ',this%gnd_alb_dif(:)
       write(fates_log(),*) 'pa%gnd_alb_dir        = ',this%gnd_alb_dir(:)
       write(fates_log(),*) 'pa%c_stomata          = ',this%c_stomata

biogeophys/FatesPlantRespPhotosynthMod.F90

@@ -649,7 +649,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
 
                                  else    ! Two-stream
 
-                                    if(cohort_layer_elai(iv) > nearzero .and. currentPatch%solar_zenith_flag) then
+                                    if(cohort_layer_elai(iv) > nearzero .and. sites(s)%coszen>0._r8 ) then
 
                                        call FatesGetCohortAbsRad(currentPatch, currentCohort, ipar, &
                                             cohort_vaitop(iv), cohort_vaibot(iv), cohort_elai, cohort_esai, &

main/EDInitMod.F90

@@ -282,6 +282,7 @@ subroutine zero_site( site_in )
     site_in%smp_memory(:,:)  = nan
     site_in%vegtemp_memory(:) = nan              ! record of last 10 days temperature for senescence model.
 
+    site_in%coszen = 0._r8
     site_in%phen_model_date  = fates_unset_int
 
     ! Disturbance rates tracking

main/EDTypesMod.F90

@@ -357,6 +357,9 @@ module EDTypesMod
      ! Total area of patches in each age bin [m2]
      real(r8), allocatable :: area_by_age(:)
 
+     ! Cosine of the zenith angle of the sun
+     real(r8) :: coszen
+     
 
      ! Nutrient relevant 
      real(r8), allocatable :: rec_l2fr(:,:) ! A running mean of the l2fr's for the newly

main/FatesInterfaceMod.F90

@@ -293,8 +293,6 @@ subroutine zero_bcs(fates,s)
     fates%bc_in(s)%watsat_sl(:)        = 0.0_r8
     fates%bc_in(s)%tempk_sl(:)         = 0.0_r8
     fates%bc_in(s)%h2o_liqvol_sl(:)    = 0.0_r8
-    fates%bc_in(s)%filter_vegzen_pa(:) = .false.
-    fates%bc_in(s)%coszen_pa(:)        = 0.0_r8
     fates%bc_in(s)%fcansno_pa(:)       = 0.0_r8
     fates%bc_in(s)%albgr_dir_rb(:)     = 0.0_r8
     fates%bc_in(s)%albgr_dif_rb(:)     = 0.0_r8
@@ -375,13 +373,13 @@ subroutine zero_bcs(fates,s)
     fates%bc_out(s)%rssun_pa(:)     = 0.0_r8
     fates%bc_out(s)%rssha_pa(:)     = 0.0_r8
 
-    fates%bc_out(s)%albd_parb(:,:) = 0.0_r8
-    fates%bc_out(s)%albi_parb(:,:) = 0.0_r8
-    fates%bc_out(s)%fabd_parb(:,:) = 0.0_r8
-    fates%bc_out(s)%fabi_parb(:,:) = 0.0_r8
-    fates%bc_out(s)%ftdd_parb(:,:) = 0.0_r8
-    fates%bc_out(s)%ftid_parb(:,:) = 0.0_r8
-    fates%bc_out(s)%ftii_parb(:,:) = 0.0_r8
+    fates%bc_out(s)%albd_parb(:,:) = 0.0_r8   ! zero albedo, soil absorbs all rad
+    fates%bc_out(s)%albi_parb(:,:) = 0.0_r8   ! zero albedo, soil absorbs all rad
+    fates%bc_out(s)%fabd_parb(:,:) = 0.0_r8   ! no rad absorbed by veg
+    fates%bc_out(s)%fabi_parb(:,:) = 0.0_r8   ! no rad absorbed by veg
+    fates%bc_out(s)%ftdd_parb(:,:) = 1.0_r8   ! rad flux to soil at bottom of veg is 100%
+    fates%bc_out(s)%ftid_parb(:,:) = 1.0_r8   ! rad flux to soil at bottom of veg is 100%
+    fates%bc_out(s)%ftii_parb(:,:) = 1.0_r8   ! rad flux to soil at bottom of veg is 100%
 
     fates%bc_out(s)%elai_pa(:)   = 0.0_r8
     fates%bc_out(s)%esai_pa(:)   = 0.0_r8
@@ -533,8 +531,7 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats,
       allocate(bc_in%t_soisno_sl(nlevsoil_in))
 
       ! Canopy Radiation
-      allocate(bc_in%filter_vegzen_pa(maxpatch_total))
-      allocate(bc_in%coszen_pa(maxpatch_total))
+      bc_in%coszen = nan
       allocate(bc_in%fcansno_pa(maxpatch_total))
       allocate(bc_in%albgr_dir_rb(num_swb))
       allocate(bc_in%albgr_dif_rb(num_swb))
@@ -621,13 +618,15 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in)
       allocate(bc_out%rssha_pa(maxpatch_total))
 
       ! Canopy Radiation
-      allocate(bc_out%albd_parb(fates_maxPatchesPerSite,num_swb))
-      allocate(bc_out%albi_parb(fates_maxPatchesPerSite,num_swb))
-      allocate(bc_out%fabd_parb(fates_maxPatchesPerSite,num_swb))
-      allocate(bc_out%fabi_parb(fates_maxPatchesPerSite,num_swb))
-      allocate(bc_out%ftdd_parb(fates_maxPatchesPerSite,num_swb))
-      allocate(bc_out%ftid_parb(fates_maxPatchesPerSite,num_swb))
-      allocate(bc_out%ftii_parb(fates_maxPatchesPerSite,num_swb))
+
+     
+      allocate(bc_out%albd_parb(maxpatch_total,num_swb))
+      allocate(bc_out%albi_parb(maxpatch_total,num_swb))
+      allocate(bc_out%fabd_parb(maxpatch_total,num_swb))
+      allocate(bc_out%fabi_parb(maxpatch_total,num_swb))
+      allocate(bc_out%ftdd_parb(maxpatch_total,num_swb))
+      allocate(bc_out%ftid_parb(maxpatch_total,num_swb))
+      allocate(bc_out%ftii_parb(maxpatch_total,num_swb))
 
       ! We allocate the boundary conditions to the BGC
       ! model, regardless of what scheme we use. The BGC

main/FatesInterfaceTypesMod.F90

@@ -522,17 +522,10 @@ module FatesInterfaceTypesMod
 
       ! Canopy Radiation Boundaries
       ! ---------------------------------------------------------------------------------
-      
-      ! Filter for vegetation patches with a positive zenith angle (daylight)
-      logical, allocatable :: filter_vegzen_pa(:)
-
-      ! Cosine of the zenith angle (0-1), by patch
-      ! Note RGK: It does not seem like the code would currently generate
-      !           different zenith angles for different patches (nor should it)
-      !           I am leaving it at this scale for simplicity.  Patches should
-      !           have no spacially variable information
-      real(r8), allocatable :: coszen_pa(:)
 
+      ! Cosine of the zenith angle (0-1) - site level
+      real(r8) :: coszen
+      
       ! fraction of canopy that is covered in snow
       real(r8), allocatable :: fcansno_pa(:)
 
@@ -670,6 +663,9 @@ module FatesInterfaceTypesMod
 
       ! Canopy Radiation Boundaries
       ! ---------------------------------------------------------------------------------
+
+      ! Note: We initialize and default the radiatioon balance to assume that the
+      ! canopy is invisible, and the soil absorbs all radiation.
 
       ! Surface albedo (direct) (HLMs use this for atm coupling and balance checks)
       real(r8), allocatable :: albd_parb(:,:)