two-stream - b4b Norman and no history changes

biogeochem/EDPatchDynamicsMod.F90

@@ -47,7 +47,6 @@ module EDPatchDynamicsMod
   use FatesConstantsMod    , only : fates_tiny
   use FatesConstantsMod    , only : nocomp_bareground
   use FatesInterfaceTypesMod    , only : hlm_use_planthydro
-  use FatesInterfaceTypesMod    , only : hlm_numSWb
   use FatesInterfaceTypesMod    , only : bc_in_type
   use FatesInterfaceTypesMod    , only : numpft
   use FatesInterfaceTypesMod    , only : hlm_stepsize
@@ -106,6 +105,7 @@ module EDPatchDynamicsMod
   use FatesRunningMeanMod,    only : ema_sdlng_mdd
   use FatesRunningMeanMod,    only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par, ema_sdlng2sap_par
   use FatesRunningMeanMod,    only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm
+  use FatesRadiationMemMod,   only : num_swb
 
   ! CIME globals
   use shr_infnan_mod       , only : nan => shr_infnan_nan, assignment(=)
@@ -603,7 +603,7 @@ subroutine spawn_patches( currentSite, bc_in)
                    allocate(newPatch)
 
                    call newPatch%Create(age, site_areadis, i_landusechange_receiverpatchlabel, i_nocomp_pft, &
-                                         hlm_numSWb, numpft, currentSite%nlevsoil, hlm_current_tod,              &
+                                         num_swb, numpft, currentSite%nlevsoil, hlm_current_tod,              &
                                          regeneration_model)
 
                    ! Initialize the litter pools to zero, these
@@ -2773,8 +2773,9 @@ 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
-    rp%radiation_error      = (dp%radiation_error*dp%area + rp%radiation_error*rp%area) * inv_sum_area
-
+    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
+
     rp%area = rp%area + dp%area !THIS MUST COME AT THE END!
 
     !insert donor cohorts into recipient patch

biogeochem/FatesAllometryMod.F90

@@ -92,12 +92,13 @@ module FatesAllometryMod
   use FatesConstantsMod, only : calloc_abs_error
   use FatesConstantsMod, only : fates_unset_r8
   use FatesConstantsMod, only : itrue
+  use FatesConstantsMod, only : nearzero
   use shr_log_mod      , only : errMsg => shr_log_errMsg
   use FatesGlobals     , only : fates_log
   use FatesGlobals     , only : endrun => fates_endrun
   use FatesGlobals     , only : FatesWarn,N2S,A2S,I2S
-  use EDParamsMod      , only : nlevleaf, dinc_vai
-  use EDParamsMod      , only : nclmax
+  use EDParamsMod      , only : nlevleaf,dinc_vai,dlower_vai
+  use EDParamsMod       , only : nclmax
   use DamageMainMod    , only : GetCrownReduction
 
   implicit none
@@ -125,7 +126,8 @@ module FatesAllometryMod
   public :: set_root_fraction  ! Generic wrapper to calculate normalized
                                ! root profiles
   public :: leafc_from_treelai ! Calculate target leaf carbon for a given treelai for SP mode
-
+  public :: VegAreaLayer
+
   logical         , parameter :: verbose_logging = .false.
   character(len=*), parameter :: sourcefile = __FILE__
 
@@ -2436,7 +2438,8 @@ real(r8) function decay_coeff_kn(pft,vcmax25top)
   end function decay_coeff_kn
 
   ! =====================================================================================
-subroutine ForceDBH( ipft, crowndamage, canopy_trim, elongf_leaf, elongf_stem, d, h, bdead, bl )
+
+  subroutine ForceDBH( ipft, crowndamage, canopy_trim, elongf_leaf, elongf_stem, d, h, bdead, bl )
 
      ! =========================================================================
      ! This subroutine estimates the diameter based on either the structural biomass
@@ -2586,6 +2589,111 @@ subroutine ForceDBH( ipft, crowndamage, canopy_trim, elongf_leaf, elongf_stem, d
      return
   end subroutine ForceDBH
 
+  ! =========================================================================
+
+  subroutine VegAreaLayer(tree_lai,tree_sai,tree_height,iv,nv,pft,snow_depth, & 
+       vai_top,vai_bot, elai_layer,esai_layer,tlai_layer,tsai_layer)
+
+    ! -----------------------------------------------------------------------------------
+    ! This routine returns the exposed leaf and stem areas (m2 of leaf and stem) per m2 of
+    ! ground inside the crown, for the leaf-layer specified.
+    ! -----------------------------------------------------------------------------------
+
+    real(r8),intent(in) :: tree_lai         ! the in-crown leaf area index for the plant
+                                            ! [m2 leaf/m2 crown footprint]
+    real(r8),intent(in) :: tree_sai         ! the in-crown stem area index for the plant
+                                            ! [m2 stem/m2 crown footprint]
+    real(r8),intent(in) :: tree_height      ! the height of the plant [m]
+    integer,intent(in)  :: iv               ! vegetation layer index
+    integer,intent(in)  :: nv               ! this plants total number of veg layers
+    integer,intent(in)  :: pft              ! plant functional type index
+    real(r8),intent(in) :: snow_depth       ! the depth of snow on the ground [m]
+    real(r8),intent(out) :: vai_top
+    real(r8),intent(out) :: vai_bot          ! the VAI of the bin top and bottom
+    real(r8),intent(out) :: elai_layer       ! exposed leaf area index of the layer
+    real(r8),intent(out) :: esai_layer       ! exposed stem area index of the layer
+    real(r8),optional,intent(out) :: tlai_layer       ! total leaf area index of the layer
+    real(r8),optional,intent(out) :: tsai_layer       ! total stem area index of the layer
+
+                                 ! [m2 of leaf in bin / m2 crown footprint]
+    real(r8) :: tree_vai         ! the in-crown veg area index for the plant
+    real(r8) :: fraction_exposed ! fraction of the veg media that is above snow
+    real(r8) :: layer_top_height ! Physical height of the layer top relative to ground [m]
+    real(r8) :: layer_bot_height ! Physical height of the layer bottom relative to ground [m]
+    real(r8) :: tlai,tsai        ! temporary total area indices [m2/m2]
+    real(r8) :: fleaf            ! fraction of biomass in layer that is leaf
+    real(r8) :: remainder        ! old-method: remainder of biomass in last bin
+    integer, parameter :: layer_height_const_depth = 1 ! constant physical depth assumption
+    integer, parameter :: layer_height_const_lad   = 2 ! constant leaf area depth assumption
+    integer, parameter :: layer_height_method = layer_height_const_depth
+
+    tree_vai = tree_lai + tree_sai
+
+    if_any_vai: if(tree_vai>0._r8)then
+
+       if(iv==0)then
+          vai_top = 0.0
+          vai_bot = tree_vai
+       else
+
+          if(iv>1)then
+             vai_top = dlower_vai(iv) - dinc_vai(iv)
+          else
+             vai_top = 0._r8
+          end if
+
+          if(iv<nv) then
+             vai_bot = dlower_vai(iv)
+          else
+             vai_bot = tree_vai
+          end if
+       end if
+
+       if(layer_height_method .eq. layer_height_const_depth)then
+          if(iv==0)then
+             layer_top_height = tree_height
+             layer_bot_height = tree_height*(1._r8 - prt_params%crown_depth_frac(pft))
+          else
+             layer_top_height = tree_height*(1._r8 - real(iv-1,r8)/real(nv,r8)*prt_params%crown_depth_frac(pft))
+             layer_bot_height = tree_height*(1._r8 - real(iv,r8)/real(nv,r8)*prt_params%crown_depth_frac(pft))
+          end if
+       else
+          layer_top_height = tree_height*(1._r8 - prt_params%crown_depth_frac(pft)*vai_top/tree_vai)
+          layer_bot_height = tree_height*(1._r8 - prt_params%crown_depth_frac(pft)*vai_bot/tree_vai)
+       end if
+
+       fraction_exposed =  1._r8 - max(0._r8,(min(1._r8, (snow_depth-layer_bot_height)/(layer_top_height-layer_bot_height))))
+
+       tlai = (vai_bot-vai_top) * tree_lai / tree_vai
+       tsai = (vai_bot-vai_top) * tree_sai / tree_vai
+
+       if(present(tlai_layer)) tlai_layer = tlai
+       if(present(tsai_layer)) tsai_layer = tsai
+
+       elai_layer = fraction_exposed * tlai
+       esai_layer = fraction_exposed * tsai
+
+       ! Update the vai at the bottom to be removed/decreased if there is no exposure
+       ! set the vai top and bottoms to the snow layer if below
+       !vai_top = min(vai_top,fraction_exposed*tree_vai)
+
+       vai_bot = vai_top + fraction_exposed*(vai_bot-vai_top)
+
+    else
+
+       if(present(tlai_layer)) tlai_layer = 0._r8
+       if(present(tsai_layer)) tsai_layer = 0._r8
+       elai_layer = 0._r8
+       esai_layer = 0._r8
+       vai_bot = 0._r8
+       vai_top = 0._r8
+
+    end if if_any_vai
+
+
+    return
+  end subroutine VegAreaLayer
+
   ! =========================================================================
 
   subroutine cspline(x1,x2,y1,y2,dydx1,dydx2,x,y,dydx)

biogeochem/FatesCohortMod.F90

@@ -200,6 +200,8 @@ module FatesCohortMod
 
     !---------------------------------------------------------------------------
 
+    integer :: twostr_col  ! The column index in the two-stream solution that this cohort is part of
+
     ! RESPIRATION COMPONENTS
     real(r8) :: rdark            ! dark respiration [kgC/indiv/s]
     real(r8) :: resp_g_tstep     ! growth respiration [kgC/indiv/timestep]

biogeochem/FatesPatchMod.F90

@@ -16,14 +16,16 @@ module FatesPatchMod
   use FatesLitterMod,      only : litter_type
   use PRTGenericMod,       only : num_elements
   use PRTGenericMod,       only : element_list
-  use EDParamsMod,         only : maxSWb, nlevleaf, nclmax, maxpft
+  use EDParamsMod,         only : nlevleaf, nclmax, maxpft
   use FatesConstantsMod,   only : n_dbh_bins, n_dist_types
-  use FatesConstantsMod,   only : n_rad_stream_types
   use FatesConstantsMod,   only : t_water_freeze_k_1atm
   use FatesRunningMeanMod, only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm
   use FatesRunningMeanMod, only : ema_sdlng_emerg_h2o, ema_sdlng_mort_par
   use FatesRunningMeanMod, only : ema_sdlng2sap_par, ema_sdlng_mdd
-
+  use TwoStreamMLPEMod,    only : twostream_type
+  use FatesRadiationMemMod,only : num_swb
+  use FatesRadiationMemMod,only : num_rad_stream_types
+  use FatesInterfaceTypesMod,only : hlm_hio_ignore_val
   use shr_infnan_mod,      only : nan => shr_infnan_nan, assignment(=)
   use shr_log_mod,         only : errMsg => shr_log_errMsg
 
@@ -111,40 +113,36 @@ module FatesPatchMod
     real(r8) :: c_stomata                                   ! mean stomatal conductance of all leaves in the patch   [umol/m2/s]
     real(r8) :: c_lblayer                                   ! mean boundary layer conductance of all leaves in the patch [umol/m2/s]
 
-    !TODO - can we delete these?
-    real(r8) :: layer_height_profile(nclmax,maxpft,nlevleaf)
     real(r8) :: psn_z(nclmax,maxpft,nlevleaf)
-    real(r8) :: nrmlzd_parprof_pft_dir_z(n_rad_stream_types,nclmax,maxpft,nlevleaf)
-    real(r8) :: nrmlzd_parprof_pft_dif_z(n_rad_stream_types,nclmax,maxpft,nlevleaf)
-    real(r8) :: nrmlzd_parprof_dir_z(n_rad_stream_types,nclmax,nlevleaf)
-    real(r8) :: nrmlzd_parprof_dif_z(n_rad_stream_types,nclmax,nlevleaf)
+    real(r8) :: nrmlzd_parprof_pft_dir_z(num_rad_stream_types,nclmax,maxpft,nlevleaf)
+    real(r8) :: nrmlzd_parprof_pft_dif_z(num_rad_stream_types,nclmax,maxpft,nlevleaf)
 
     !---------------------------------------------------------------------------
 
     ! RADIATION
-    real(r8) :: radiation_error                           ! radiation error [W/m2] 
+    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(maxSWb)                       ! ground albedo for diffuse rad, both bands [0-1]
-    real(r8) :: gnd_alb_dir(maxSWb)                       ! ground albedo for direct rad, both bands [0-1]
-
+    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]
+
+
     ! organized by canopy layer, pft, and leaf layer
     real(r8) :: fabd_sun_z(nclmax,maxpft,nlevleaf)        ! sun fraction of direct light absorbed [0-1]
     real(r8) :: fabd_sha_z(nclmax,maxpft,nlevleaf)        ! shade fraction of direct light absorbed [0-1]
     real(r8) :: fabi_sun_z(nclmax,maxpft,nlevleaf)        ! sun fraction of indirect light absorbed [0-1]
     real(r8) :: fabi_sha_z(nclmax,maxpft,nlevleaf)        ! shade fraction of indirect light absorbed [0-1]
     real(r8) :: ed_parsun_z(nclmax,maxpft,nlevleaf)       ! PAR absorbed in the sun [W/m2]   
     real(r8) :: ed_parsha_z(nclmax,maxpft,nlevleaf)       ! PAR absorbed in the shade [W/m2]
+    real(r8) :: f_sun(nclmax,maxpft,nlevleaf)             ! fraction of leaves in the sun [0-1]
     real(r8) :: ed_laisun_z(nclmax,maxpft,nlevleaf)
     real(r8) :: ed_laisha_z(nclmax,maxpft,nlevleaf)
-    real(r8) :: f_sun(nclmax,maxpft,nlevleaf)             ! fraction of leaves in the sun [0-1]
 
+
     ! radiation profiles for comparison against observations
     real(r8) :: parprof_pft_dir_z(nclmax,maxpft,nlevleaf) ! direct-beam PAR profile through canopy, by canopy, PFT, leaf level [W/m2]
     real(r8) :: parprof_pft_dif_z(nclmax,maxpft,nlevleaf) ! diffuse     PAR profile through canopy, by canopy, PFT, leaf level [W/m2]
-    real(r8) :: parprof_dir_z(nclmax,nlevleaf)            ! direct-beam PAR profile through canopy, by canopy, leaf level [W/m2]
-    real(r8) :: parprof_dif_z(nclmax,nlevleaf)            ! diffuse     PAR profile through canopy, by canopy, leaf level [W/m2]
 
     real(r8), allocatable :: tr_soil_dir(:)               ! fraction of incoming direct radiation transmitted to the soil as direct, by numSWB [0-1]
     real(r8), allocatable :: tr_soil_dif(:)               ! fraction of incoming diffuse radiation that is transmitted to the soil as diffuse [0-1]
@@ -155,6 +153,10 @@ module FatesPatchMod
     real(r8), allocatable :: sabs_dir(:)                  ! fraction of incoming direct  radiation that is absorbed by the canopy
     real(r8), allocatable :: sabs_dif(:)                  ! fraction of incoming diffuse radiation that is absorbed by the canopy
 
+    ! Twostream data structures
+    type(twostream_type) :: twostr                        ! This holds all two-stream data and procedures
+
+
     !---------------------------------------------------------------------------
 
     ! ROOTS
@@ -318,16 +320,13 @@ subroutine NanValues(this)
       this%ncan(:,:)                    = fates_unset_int
       this%c_stomata                    = nan 
       this%c_lblayer                    = nan
-      this%layer_height_profile(:,:,:)  = nan
 
       this%psn_z(:,:,:)                 = nan 
       this%nrmlzd_parprof_pft_dir_z(:,:,:,:) = nan
       this%nrmlzd_parprof_pft_dif_z(:,:,:,:) = nan
-      this%nrmlzd_parprof_dir_z(:,:,:)  = nan
-      this%nrmlzd_parprof_dir_z(:,:,:) = nan
 
       ! RADIATION
-      this%radiation_error              = nan 
+      this%rad_error(:)                 = nan
       this%fcansno                      = nan 
       this%solar_zenith_flag            = .false. 
       this%solar_zenith_angle           = nan 
@@ -336,16 +335,14 @@ subroutine NanValues(this)
       this%fabd_sun_z(:,:,:)            = nan 
       this%fabd_sha_z(:,:,:)            = nan 
       this%fabi_sun_z(:,:,:)            = nan 
-      this%fabi_sha_z(:,:,:)            = nan  
+      this%fabi_sha_z(:,:,:)            = nan
       this%ed_laisun_z(:,:,:)           = nan 
       this%ed_laisha_z(:,:,:)           = nan 
       this%ed_parsun_z(:,:,:)           = nan 
       this%ed_parsha_z(:,:,:)           = nan 
       this%f_sun(:,:,:)                 = nan
       this%parprof_pft_dir_z(:,:,:)     = nan 
       this%parprof_pft_dif_z(:,:,:)     = nan
-      this%parprof_dir_z(:,:)           = nan
-      this%parprof_dif_z(:,:)           = nan
       this%tr_soil_dir(:)               = nan    
       this%tr_soil_dif(:)               = nan    
       this%tr_soil_dir_dif(:)           = nan
@@ -418,19 +415,17 @@ subroutine ZeroValues(this)
       this%psn_z(:,:,:)                      = 0.0_r8
       this%nrmlzd_parprof_pft_dir_z(:,:,:,:) = 0.0_r8
       this%nrmlzd_parprof_pft_dif_z(:,:,:,:) = 0.0_r8
-      this%nrmlzd_parprof_dir_z(:,:,:)       = 0.0_r8
-      this%nrmlzd_parprof_dif_z(:,:,:)       = 0.0_r8
 
       ! RADIATION
-      this%radiation_error                   = 0.0_r8
+      this%rad_error(:)                      = 0.0_r8
       this%fabd_sun_z(:,:,:)                 = 0.0_r8 
       this%fabd_sha_z(:,:,:)                 = 0.0_r8 
       this%fabi_sun_z(:,:,:)                 = 0.0_r8 
       this%fabi_sha_z(:,:,:)                 = 0.0_r8  
       this%ed_parsun_z(:,:,:)                = 0.0_r8 
-      this%ed_parsha_z(:,:,:)                = 0.0_r8 
-      this%ed_laisun_z(:,:,:)                = 0.0_r8
-      this%ed_laisha_z(:,:,:)                = 0.0_r8 
+      this%ed_parsha_z(:,:,:)                = 0.0_r8
+      this%ed_laisun_z(:,:,:)           = 0._r8
+      this%ed_laisha_z(:,:,:)           = 0._r8
       this%f_sun                             = 0.0_r8
       this%tr_soil_dir_dif(:)                = 0.0_r8
       this%fab(:)                            = 0.0_r8
@@ -591,7 +586,11 @@ subroutine Create(this, age, area, label, nocomp_pft, num_swb, num_pft,    &
 
       ! initialize litter
       call this%InitLitter(num_pft, num_levsoil)
-
+
+      this%twostr%scelg => null()  ! The radiation module will check if this
+                                   ! is associated, since it is not, it will then
+                                   ! initialize and allocate
+
       ! assign known patch attributes 
       this%age       = age   
       this%age_class = 1
@@ -646,7 +645,12 @@ subroutine FreeMemory(this, regeneration_model, numpft)
         endif
         ccohort => ncohort
       end do
-
+
+      ! Deallocate Radiation scattering elements
+      if(associated(this%twostr%scelg)) then
+         call this%twostr%DeallocTwoStream()
+      end if
+
       ! deallocate all litter objects
       do el=1,num_elements
         call this%litter(el)%DeallocateLitt()