cleaning up nocomp patch indexing

biogeochem/EDCanopyStructureMod.F90

@@ -49,7 +49,7 @@ module EDCanopyStructureMod
   use PRTGenericMod,          only : carbon12_element
   use FatesTwoStreamUtilsMod, only : FatesConstructRadElements
   use FatesRadiationMemMod  , only : twostr_solver
-  use FatesRadiationMemMod  , only : num_rad_stream_types
+
 
   ! CIME Globals
   use shr_log_mod           , only : errMsg => shr_log_errMsg
@@ -1909,127 +1909,120 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
        c = fcolumn(s)
        do while(associated(currentPatch))
 
-          if(currentPatch%nocomp_pft_label.ne.nocomp_bareground)then  ! ignore the bare-ground-PFT patch entirely for these BC outs
+          ifp = ifp+1
 
-             ifp = ifp+1
+          if ( currentPatch%total_canopy_area-currentPatch%area > 0.000001_r8 ) then
+             if(debug)then
+                write(fates_log(),*) 'ED: canopy area bigger than area', &
+                     currentPatch%total_canopy_area ,currentPatch%area
+             end if
+             currentPatch%total_canopy_area = currentPatch%area
+          endif
 
-             if ( currentPatch%total_canopy_area-currentPatch%area > 0.000001_r8 ) then
-                if(debug)then
-                   write(fates_log(),*) 'ED: canopy area bigger than area', &
-                        currentPatch%total_canopy_area ,currentPatch%area
-                end if
-                currentPatch%total_canopy_area = currentPatch%area
-             endif
+          if (associated(currentPatch%tallest)) then
+             bc_out(s)%htop_pa(ifp) = currentPatch%tallest%height
+          else
+             ! FIX(RF,040113) - should this be a parameter for the minimum possible vegetation height?
+             bc_out(s)%htop_pa(ifp) = 0.1_r8
+          endif
 
-             if (associated(currentPatch%tallest)) then
-                bc_out(s)%htop_pa(ifp) = currentPatch%tallest%height
-             else
-                ! FIX(RF,040113) - should this be a parameter for the minimum possible vegetation height?
-                bc_out(s)%htop_pa(ifp) = 0.1_r8
-             endif
+          bc_out(s)%hbot_pa(ifp) = max(0._r8, min(0.2_r8, bc_out(s)%htop_pa(ifp)- 1.0_r8))
 
-             bc_out(s)%hbot_pa(ifp) = max(0._r8, min(0.2_r8, bc_out(s)%htop_pa(ifp)- 1.0_r8))
+          ! Use canopy-only crown area weighting for all cohorts in the patch to define the characteristic
+          ! Roughness length and displacement height used by the HLM
+          ! use total LAI + SAI to weight the leaft characteristic dimension
+          ! Avoid this if running in satellite phenology mode
+          ! ----------------------------------------------------------------------------
 
-             ! Use canopy-only crown area weighting for all cohorts in the patch to define the characteristic
-             ! Roughness length and displacement height used by the HLM
-             ! use total LAI + SAI to weight the leaft characteristic dimension
-             ! Avoid this if running in satellite phenology mode
-             ! ----------------------------------------------------------------------------
+          if (currentPatch%total_canopy_area > nearzero) then
+             currentCohort => currentPatch%shortest
+             do while(associated(currentCohort))
+                if (currentCohort%canopy_layer .eq. 1) then
+                   weight = min(1.0_r8,currentCohort%c_area/currentPatch%total_canopy_area)
+                   bc_out(s)%z0m_pa(ifp) = bc_out(s)%z0m_pa(ifp) + &
+                        EDPftvarcon_inst%z0mr(currentCohort%pft) * currentCohort%height * weight
+                   bc_out(s)%displa_pa(ifp) = bc_out(s)%displa_pa(ifp) + &
+                        EDPftvarcon_inst%displar(currentCohort%pft) * currentCohort%height * weight
+                endif
+                currentCohort => currentCohort%taller
+             end do
 
-             if (currentPatch%total_canopy_area > nearzero) then
-                currentCohort => currentPatch%shortest
-                do while(associated(currentCohort))
-                   if (currentCohort%canopy_layer .eq. 1) then
-                      weight = min(1.0_r8,currentCohort%c_area/currentPatch%total_canopy_area)
-                      bc_out(s)%z0m_pa(ifp) = bc_out(s)%z0m_pa(ifp) + &
-                           EDPftvarcon_inst%z0mr(currentCohort%pft) * currentCohort%height * weight
-                      bc_out(s)%displa_pa(ifp) = bc_out(s)%displa_pa(ifp) + &
-                           EDPftvarcon_inst%displar(currentCohort%pft) * currentCohort%height * weight
-                   endif
-                   currentCohort => currentCohort%taller
-                end do
+             ! for lai, scale to total LAI + SAI in patch.  first add up all the LAI and SAI in the patch
+             total_patch_leaf_stem_area = 0._r8
+             currentCohort => currentPatch%shortest
+             do while(associated(currentCohort))
+                total_patch_leaf_stem_area = total_patch_leaf_stem_area + &
+                     (currentCohort%treelai + currentCohort%treesai) * currentCohort%c_area
+                currentCohort => currentCohort%taller
+             end do
 
-                ! for lai, scale to total LAI + SAI in patch.  first add up all the LAI and SAI in the patch
-                total_patch_leaf_stem_area = 0._r8
+             ! make sure there is some leaf and stem area
+             if (total_patch_leaf_stem_area > nearzero) then
                 currentCohort => currentPatch%shortest
                 do while(associated(currentCohort))
-                   total_patch_leaf_stem_area = total_patch_leaf_stem_area + &
-                        (currentCohort%treelai + currentCohort%treesai) * currentCohort%c_area
+                   ! weight dleaf by the relative totals of leaf and stem area
+                   weight = (currentCohort%treelai + currentCohort%treesai) * currentCohort%c_area / total_patch_leaf_stem_area
+                   bc_out(s)%dleaf_pa(ifp) = bc_out(s)%dleaf_pa(ifp) + &
+                        EDPftvarcon_inst%dleaf(currentCohort%pft) * weight
                    currentCohort => currentCohort%taller
                 end do
-
-                ! make sure there is some leaf and stem area
-                if (total_patch_leaf_stem_area > nearzero) then
-                   currentCohort => currentPatch%shortest
-                   do while(associated(currentCohort))
-                      ! weight dleaf by the relative totals of leaf and stem area
-                      weight = (currentCohort%treelai + currentCohort%treesai) * currentCohort%c_area / total_patch_leaf_stem_area
-                      bc_out(s)%dleaf_pa(ifp) = bc_out(s)%dleaf_pa(ifp) + &
-                           EDPftvarcon_inst%dleaf(currentCohort%pft) * weight
-                      currentCohort => currentCohort%taller
-                   end do
-                else
-                   ! dummy case
-                   bc_out(s)%dleaf_pa(ifp)  = EDPftvarcon_inst%dleaf(1)
-                endif
              else
-                ! if no canopy, then use dummy values (first PFT) of aerodynamic properties
-                bc_out(s)%z0m_pa(ifp)    = EDPftvarcon_inst%z0mr(1) * bc_out(s)%htop_pa(ifp)
-                bc_out(s)%displa_pa(ifp) = EDPftvarcon_inst%displar(1) * bc_out(s)%htop_pa(ifp)
+                ! dummy case
                 bc_out(s)%dleaf_pa(ifp)  = EDPftvarcon_inst%dleaf(1)
              endif
-             ! -----------------------------------------------------------------------------
-
-             ! We are assuming here that grass is all located underneath tree canopies.
-             ! The alternative is to assume it is all spatial distinct from tree canopies.
-             ! In which case, the bare area would have to be reduced by the grass area...
-             ! currentPatch%total_canopy_area/currentPatch%area is fraction of this patch cover by plants
-             ! currentPatch%area/AREA is the fraction of the soil covered by this patch.
-
-             if(currentPatch%area.gt.0.0_r8)then
-                bc_out(s)%canopy_fraction_pa(ifp) = &
-                     min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)*(currentPatch%area/AREA)
-             else
-                bc_out(s)%canopy_fraction_pa(ifp) = 0.0_r8
-             endif
-
-             bare_frac_area = (1.0_r8 - min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)) * &
-                  (currentPatch%area/AREA)
+          else
+             ! if no canopy, then use dummy values (first PFT) of aerodynamic properties
+             bc_out(s)%z0m_pa(ifp)    = EDPftvarcon_inst%z0mr(1) * bc_out(s)%htop_pa(ifp)
+             bc_out(s)%displa_pa(ifp) = EDPftvarcon_inst%displar(1) * bc_out(s)%htop_pa(ifp)
+             bc_out(s)%dleaf_pa(ifp)  = EDPftvarcon_inst%dleaf(1)
+          endif
+          ! -----------------------------------------------------------------------------
 
-             total_patch_area = total_patch_area + bc_out(s)%canopy_fraction_pa(ifp) + bare_frac_area
+          ! We are assuming here that grass is all located underneath tree canopies.
+          ! The alternative is to assume it is all spatial distinct from tree canopies.
+          ! In which case, the bare area would have to be reduced by the grass area...
+          ! currentPatch%total_canopy_area/currentPatch%area is fraction of this patch cover by plants
+          ! currentPatch%area/AREA is the fraction of the soil covered by this patch.
 
-             total_canopy_area = total_canopy_area + bc_out(s)%canopy_fraction_pa(ifp)
+          if(currentPatch%area.gt.0.0_r8)then
+             bc_out(s)%canopy_fraction_pa(ifp) = &
+                  min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)*(currentPatch%area/AREA)
+          else
+             bc_out(s)%canopy_fraction_pa(ifp) = 0.0_r8
+          endif
 
-             bc_out(s)%nocomp_pft_label_pa(ifp) = currentPatch%nocomp_pft_label
+          bare_frac_area = (1.0_r8 - min(1.0_r8,currentPatch%total_canopy_area/currentPatch%area)) * &
+               (currentPatch%area/AREA)
 
-             ! Calculate area indices for output boundary to HLM
-             ! It is assumed that cpatch%canopy_area_profile and cpat%xai_profiles
-             ! have been updated (ie ed_leaf_area_profile has been called since dynamics has been called)
+          total_patch_area = total_patch_area + bc_out(s)%canopy_fraction_pa(ifp) + bare_frac_area
 
-             bc_out(s)%elai_pa(ifp) = calc_areaindex(currentPatch,'elai')
-             bc_out(s)%tlai_pa(ifp) = calc_areaindex(currentPatch,'tlai')
-             bc_out(s)%esai_pa(ifp) = calc_areaindex(currentPatch,'esai')
-             bc_out(s)%tsai_pa(ifp) = calc_areaindex(currentPatch,'tsai')
+          total_canopy_area = total_canopy_area + bc_out(s)%canopy_fraction_pa(ifp)
 
-             ! Fraction of vegetation free of snow. This is used to flag those
-             ! patches which shall under-go photosynthesis
-             ! INTERF-TODO: we may want to stop using frac_veg_nosno_alb and let
-             ! FATES internal variables decide if photosynthesis is possible
-             ! we are essentially calculating it inside FATES to tell the
-             ! host to tell itself when to do things (circuitous). Just have
-             ! to determine where else it is used
+          bc_out(s)%nocomp_pft_label_pa(ifp) = currentPatch%nocomp_pft_label
 
-             if ((bc_out(s)%elai_pa(ifp) + bc_out(s)%esai_pa(ifp)) > 0._r8) then
-                bc_out(s)%frac_veg_nosno_alb_pa(ifp) = 1.0_r8
-             else
-                bc_out(s)%frac_veg_nosno_alb_pa(ifp) = 0.0_r8
-             end if
+          ! Calculate area indices for output boundary to HLM
+          ! It is assumed that cpatch%canopy_area_profile and cpat%xai_profiles
+          ! have been updated (ie ed_leaf_area_profile has been called since dynamics has been called)
 
-          else  ! nocomp or SP, and currentPatch%nocomp_pft_label .eq. 0
+          bc_out(s)%elai_pa(ifp) = calc_areaindex(currentPatch,'elai')
+          bc_out(s)%tlai_pa(ifp) = calc_areaindex(currentPatch,'tlai')
+          bc_out(s)%esai_pa(ifp) = calc_areaindex(currentPatch,'esai')
+          bc_out(s)%tsai_pa(ifp) = calc_areaindex(currentPatch,'tsai')
 
-             total_patch_area = total_patch_area + currentPatch%area/AREA
+          ! Fraction of vegetation free of snow. This is used to flag those
+          ! patches which shall under-go photosynthesis
+          ! INTERF-TODO: we may want to stop using frac_veg_nosno_alb and let
+          ! FATES internal variables decide if photosynthesis is possible
+          ! we are essentially calculating it inside FATES to tell the
+          ! host to tell itself when to do things (circuitous). Just have
+          ! to determine where else it is used
 
+          if ((bc_out(s)%elai_pa(ifp) + bc_out(s)%esai_pa(ifp)) > 0._r8) then
+             bc_out(s)%frac_veg_nosno_alb_pa(ifp) = 1.0_r8
+          else
+             bc_out(s)%frac_veg_nosno_alb_pa(ifp) = 0.0_r8
           end if
+
           currentPatch => currentPatch%younger
        end do
 
@@ -2051,11 +2044,8 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
           currentPatch => sites(s)%oldest_patch
           ifp = 0
           do while(associated(currentPatch))
-             if(currentPatch%nocomp_pft_label.ne.nocomp_bareground)then ! for vegetated patches only
-                ifp = ifp+1
-                bc_out(s)%canopy_fraction_pa(ifp) = bc_out(s)%canopy_fraction_pa(ifp)/total_patch_area
-             endif ! veg patch
-
+             ifp = ifp+1
+             bc_out(s)%canopy_fraction_pa(ifp) = bc_out(s)%canopy_fraction_pa(ifp)/total_patch_area
              currentPatch => currentPatch%younger
           end do
 

biogeochem/EDPatchDynamicsMod.F90

@@ -1899,22 +1899,6 @@ subroutine set_patchno( currentSite )
        currentPatch => currentPatch%younger
     enddo
 
-    if(hlm_use_fixed_biogeog.eq.itrue .and. hlm_use_nocomp.eq.itrue)then
-      patchno = 1
-      currentPatch => currentSite%oldest_patch
-      do while(associated(currentPatch))
-        if(currentPatch%nocomp_pft_label.eq.nocomp_bareground)then
-         ! for bareground patch, we make the patch number 0
-         ! we also do not count this in the veg. patch numbering scheme.
-          currentPatch%patchno = 0
-        else
-         currentPatch%patchno = patchno
-         patchno = patchno + 1
-        endif
-        currentPatch => currentPatch%younger
-       enddo
-    endif
-
   end subroutine set_patchno
 
   ! ============================================================================

biogeochem/EDPhysiologyMod.F90

@@ -34,7 +34,6 @@ module EDPhysiologyMod
   use FatesConstantsMod, only    : mpa_per_mm_suction
   use FatesConstantsMod, only    : g_per_kg
   use FatesConstantsMod, only    : ndays_per_year
-  use FatesConstantsMod, only    : nocomp_bareground
   use FatesConstantsMod, only    : nocomp_bareground_land
   use FatesConstantsMod, only    : is_crop
   use FatesConstantsMod, only    : area_error_2
@@ -3175,44 +3174,40 @@ subroutine fragmentation_scaler( currentPatch, bc_in)
     catanf(t1) = 11.75_r8 +(29.7_r8 / pi) * atan( pi * 0.031_r8  * ( t1 - 15.4_r8 ))
     catanf_30 = catanf(30._r8)
 
-    if(currentPatch%nocomp_pft_label.ne.nocomp_bareground)then
+    ! Use the hlm temp and moisture decomp fractions by default
+    if ( use_hlm_soil_scalar ) then
 
-       ! Use the hlm temp and moisture decomp fractions by default
-       if ( use_hlm_soil_scalar ) then
+       ! Calculate the fragmentation_scaler
+       currentPatch%fragmentation_scaler =  min(1.0_r8,max(0.0_r8,bc_in%t_scalar_sisl * bc_in%w_scalar_sisl))
 
-         ! Calculate the fragmentation_scaler
-         currentPatch%fragmentation_scaler =  min(1.0_r8,max(0.0_r8,bc_in%t_scalar_sisl * bc_in%w_scalar_sisl))
+    else
 
+       if ( .not. use_century_tfunc ) then
+          !calculate rate constant scalar for soil temperature,assuming that the base rate constants
+          !are assigned for non-moisture limiting conditions at 25C.
+          if (currentPatch%tveg24%GetMean()  >=  tfrz) then
+             t_scalar = q10_mr**((currentPatch%tveg24%GetMean()-(tfrz+25._r8))/10._r8)
+             !  Q10**((t_soisno(c,j)-(tfrz+25._r8))/10._r8)
+          else
+             t_scalar = (q10_mr**(-25._r8/10._r8))*(q10_froz**((currentPatch%tveg24%GetMean()-tfrz)/10._r8))
+             !  Q10**(-25._r8/10._r8))*(froz_q10**((t_soisno(c,j)-tfrz)/10._r8)
+          endif
        else
+          ! original century uses an arctangent function to calculate the
+          ! temperature dependence of decomposition
+          t_scalar = max(catanf(currentPatch%tveg24%GetMean()-tfrz)/catanf_30,0.01_r8)
+       endif
 
-         if ( .not. use_century_tfunc ) then
-            !calculate rate constant scalar for soil temperature,assuming that the base rate constants
-            !are assigned for non-moisture limiting conditions at 25C.
-            if (currentPatch%tveg24%GetMean()  >=  tfrz) then
-               t_scalar = q10_mr**((currentPatch%tveg24%GetMean()-(tfrz+25._r8))/10._r8)
-               !  Q10**((t_soisno(c,j)-(tfrz+25._r8))/10._r8)
-            else
-               t_scalar = (q10_mr**(-25._r8/10._r8))*(q10_froz**((currentPatch%tveg24%GetMean()-tfrz)/10._r8))
-               !  Q10**(-25._r8/10._r8))*(froz_q10**((t_soisno(c,j)-tfrz)/10._r8)
-            endif
-         else
-            ! original century uses an arctangent function to calculate the
-            ! temperature dependence of decomposition
-            t_scalar = max(catanf(currentPatch%tveg24%GetMean()-tfrz)/catanf_30,0.01_r8)
-         endif
-
-         !Moisture Limitations
-         !BTRAN APPROACH - is quite simple, but max's out decomp at all unstressed
-         !soil moisture values, which is not realistic.
-         !litter decomp is proportional to water limitation on average...
-         w_scalar = sum(currentPatch%btran_ft(1:numpft))/real(numpft,r8)
+       !Moisture Limitations
+       !BTRAN APPROACH - is quite simple, but max's out decomp at all unstressed
+       !soil moisture values, which is not realistic.
+       !litter decomp is proportional to water limitation on average...
+       w_scalar = sum(currentPatch%btran_ft(1:numpft))/real(numpft,r8)
 
        ! Calculate the fragmentation_scaler
-         currentPatch%fragmentation_scaler(:) =  min(1.0_r8,max(0.0_r8,t_scalar * w_scalar))
-
-      endif ! scalar
+       currentPatch%fragmentation_scaler(:) =  min(1.0_r8,max(0.0_r8,t_scalar * w_scalar))
 
-    endif ! not bare ground
+    endif ! scalar
 
   end subroutine fragmentation_scaler