rectifying how to allocate ed-sites with the potential of dynamic col…

…umn status. 1x1br cold-starts running with what appears to be non-nonsense results.

components/clm/src/ED/main/EDCLMLinkMod.F90

@@ -1409,7 +1409,6 @@ subroutine ed_update_history_variables( this, bounds, sites, nsites, fcolumn, ca
 
             ! INTERF-TODO: THIS LOGIC SHOULDN'T BE NECESSARY, SHOULD BE CHECKED AT THE BEGINNING
             ! OF LINKING, ONCE
-
             ! %patchno is the local index of the ED/FATES patches, starting at 1
             if(currentPatch%patchno  <= numpft - numcft)then !don't expand into crop patches.   
 
@@ -1421,6 +1420,7 @@ subroutine ed_update_history_variables( this, bounds, sites, nsites, fcolumn, ca
                currentCohort => currentPatch%shortest
                do while(associated(currentCohort))
                   !accumulate into history variables. 
+
                   ft = currentCohort%pft
 
                   ed_ncohorts(c) = ed_ncohorts(c) + 1._r8
@@ -1517,9 +1517,9 @@ subroutine ed_update_history_variables( this, bounds, sites, nsites, fcolumn, ca
                !Patch specific variables that are already calculated
 
                !These things are all duplicated. Should they all be converted to LL or array structures RF? 
-
                ! define scalar to counteract the patch albedo scaling logic for conserved quantities
-               if (currentPatch%area .gt. 0._r8) then
+
+               if (currentPatch%area .gt. 0._r8 .and. currentPatch%total_canopy_area .gt.0 ) then
                   patch_scaling_scalar  = min(1._r8, currentPatch%area / currentPatch%total_canopy_area)
                else
                   patch_scaling_scalar = 0._r8
@@ -1545,7 +1545,11 @@ subroutine ed_update_history_variables( this, bounds, sites, nsites, fcolumn, ca
                seed_germination(p)     = sum(currentPatch%seed_germination) * 1.e3_r8 * 365.0_r8 * SHR_CONST_CDAY * patch_scaling_scalar
                canopy_spread(p)        = currentPatch%spread(1) 
                area_plant(p)           = 1._r8
-               area_trees(p)           = currentPatch%total_tree_area   / min(currentPatch%total_canopy_area,currentPatch%area)
+               if (min(currentPatch%total_canopy_area,currentPatch%area)>0.0_r8) then
+                  area_trees(p)           = currentPatch%total_tree_area   / min(currentPatch%total_canopy_area,currentPatch%area)
+               else
+                  area_trees(p)           = 0.0_r8
+               end if
                if(associated(currentPatch%tallest))then
                   trimming(p)          = currentPatch%tallest%canopy_trim                
                else

components/clm/src/main/clm_driver.F90

@@ -839,9 +839,9 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate)
              write(iulog,*)  'clm: calling ED model ', get_nstep()
           end if
 
-          ! INTERF-TODO: FATES(NC) SHOULD ONLY BE VISIBLE TO THE INTERFACE
-          ! AND ONLY FATES API DEFINED TYPES SHOULD BE PASSED TO IT
-          ! NEEDS A WRAPPER
+
+          call clm_fates%check_hlm_active(nc, bounds_clump)
+
           call clm_fates%dynamics_driv( nc, bounds_clump,                        &
                atm2lnd_inst, soilstate_inst, temperature_inst,                   &
                waterstate_inst, canopystate_inst)
@@ -944,7 +944,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate)
                lakestate_inst, temperature_inst, surfalb_inst)
 
 
-          ! INTERF-TOD: THIS ACTUALLY WON'T BE TO BAD TO PULL OUT
+          ! INTERF-TOD: THIS ACTUALLY WON'T BE TO HARD TO PULL OUT
           ! ED_Norman_Radiation() is the last thing called
           ! in SurfaceAlbedo, we can simply remove it
           ! The clm_fates interfac called below will split
@@ -1098,10 +1098,6 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate)
        call t_stopf('d2dgvm')
     end if
 
-    ! ============================================================================
-    ! Call ED model on daily timestep
-    ! ============================================================================
-
     ! ============================================================================
     ! History/Restart output
     ! ============================================================================

components/clm/src/main/clm_initializeMod.F90

@@ -481,6 +481,9 @@ subroutine initialize2( )
        call SatellitePhenologyInit(bounds_proc)
     end if
 
+
+
+
     ! ------------------------------------------------------------------------
     ! On restart only - process the history namelist. 
     ! ------------------------------------------------------------------------

components/clm/src/main/clm_instMod.F90

@@ -428,7 +428,7 @@ subroutine clm_instInit(bounds)
     !                          Incrementally changing to ED names to FATES
 
     call clm_fates%Init(bounds,use_ed)
-
+    call clm_fates%init_allocate()
 
     deallocate (h2osno_col)
     deallocate (snow_depth_col)

components/clm/src/utils/clmfates_interfaceMod.F90

@@ -53,6 +53,8 @@ module CLMFatesInterfaceMod
    use ColumnType        , only : col
    use LandunitType      , only : lun
    use landunit_varcon   , only : istsoil
+   use abortutils      , only : endrun
+   use shr_log_mod     , only : errMsg => shr_log_errMsg    
 
    ! Used FATES Modules
    use FatesInterfaceMod     , only : fates_interface_type
@@ -118,14 +120,23 @@ module CLMFatesInterfaceMod
       !      type(fates_bounds_type) :: bound_fate
 
 
+
+
    contains
 
       procedure, public :: init
+      procedure, public :: init_allocate
+      procedure, public :: check_hlm_active
       procedure, public :: init_restart
       procedure, public :: dynamics_driv
 
+
+
    end type hlm_fates_interface_type
 
+
+   logical :: DEBUG  = .true.
+
 contains
 
    subroutine init(this,bounds_proc, use_ed)
@@ -158,13 +169,6 @@ subroutine init(this,bounds_proc, use_ed)
                                                                  ! ONLY PART OF THIS MAY BE OPERATIVE
       ! local variables
       integer                                        :: nclumps   ! Number of threads
-      integer                                        :: nc        ! thread index
-      integer                                        :: s         ! FATES site index
-      integer                                        :: c         ! HLM column index
-      integer                                        :: l         ! HLM LU index
-      integer, allocatable                           :: collist (:)
-      type(bounds_type)                              :: bounds_clump
-      integer                                        :: nmaxcol
 
       if (use_ed) then
 
@@ -178,15 +182,53 @@ subroutine init(this,bounds_proc, use_ed)
 
       end if
 
+      if(DEBUG)then
+         write(iulog,*) 'Entering clm_fates%init'
+      end if
+
 
       nclumps = get_proc_clumps()
       allocate(this%fates(nclumps))
       allocate(this%f2hmap(nclumps))
 
+      if(DEBUG)then
+         write(iulog,*) 'clm_fates%init():  allocating for ',nclumps,' threads'
+      end if
+
+
+   end subroutine init
+
+
+   subroutine init_allocate(this)
+
+      implicit none
+
+      ! Input Arguments
+      class(hlm_fates_interface_type), intent(inout) :: this
+      ! local variables
+      integer                                        :: nclumps   ! Number of threads
+      integer                                        :: nc        ! thread index
+      integer                                        :: s         ! FATES site index
+      integer                                        :: c         ! HLM column index
+      integer                                        :: l         ! HLM LU index
+      integer, allocatable                           :: collist (:)
+      type(bounds_type)                              :: bounds_clump
+      integer                                        :: nmaxcol
+
+      if(DEBUG)then
+         write(iulog,*) 'Entering clm_fates%init_allocate'
+      end if
+
+      nclumps = get_proc_clumps()
       do nc = 1,nclumps
 
          call get_clump_bounds(nc, bounds_clump)
          nmaxcol = bounds_clump%endc - bounds_clump%begc + 1
+
+         if(DEBUG)then
+            write(iulog,*) 'clm_fates%init(): thread',nc,': allocating ',nmaxcol,' column space'
+         end if
+
          allocate(collist(1:nmaxcol))
 
          ! Allocate the mapping that points columns to FATES sites, 0 is NA
@@ -201,14 +243,27 @@ subroutine init(this,bounds_proc, use_ed)
 
             ! These are the key constraints that determine if this column
             ! will have a FATES site associated with it
-            if (col%active(c) .and. lun%itype(l) == istsoil ) then 
+
+            if(DEBUG)then
+               write(iulog,*) 'clm_fates%init(): thread',nc,': found column',c,'with lu',l
+               write(iulog,*) '  LU type:', lun%itype(l)
+            end if
+
+            ! INTERF-TODO: WE HAVE NOT FILTERED OUT FATES SITES ON INACTIVE COLUMNS.. YET
+            ! NEED A RUN-TIME ROUTINE THAT CLEARS AND REWRITES THE SITE LIST
+
+            if (lun%itype(l) == istsoil ) then
                s = s + 1
                collist(s) = c
                this%f2hmap(nc)%hsites(c) = s
             endif
 
          enddo
 
+         if(DEBUG)then
+            write(iulog,*) 'clm_fates%init(): thread',nc,': allocated ',s,' sites'
+         end if
+
          ! Allocate vectors that match FATES sites with HLM columns
          allocate(this%f2hmap(nc)%fcolumn(s))
 
@@ -228,11 +283,50 @@ subroutine init(this,bounds_proc, use_ed)
       end do
 
 
-   end subroutine init
+   end subroutine init_allocate
 
 
    ! ------------------------------------------------------------------------------------
 
+
+   subroutine check_hlm_active(this, nc, bounds_clump)
+
+
+      implicit none
+      class(hlm_fates_interface_type), intent(inout) :: this
+      integer                                        :: nc
+      type(bounds_type),intent(in)                   :: bounds_clump
+
+      ! local variables
+      integer :: c
+
+      ! FATES-TODO: THIS SHOULD BE CHANGED TO DO RE-ALLOCATION
+      ! INSTEAD OF FAILURE
+
+      do c = bounds_clump%begc,bounds_clump%endc
+
+         ! FATES ACTIVE BUT HLM IS NOT
+         if(this%f2hmap(nc)%hsites(c)>0 .and. .not.col%active(c)) then
+
+
+            write(iulog,*) 'INACTIVE COLUMN WITH ACTIVE FATES SITE'
+            write(iulog,*) 'c = ',c
+            call endrun(msg=errMsg(__FILE__, __LINE__))
+
+         elseif (this%f2hmap(nc)%hsites(c)==0 .and. col%active(c)) then
+
+            write(iulog,*) 'ACTIVE COLUMN WITH INACTIVE FATES SITE'
+            write(iulog,*) 'c = ',c
+            call endrun(msg=errMsg(__FILE__, __LINE__))
+         end if
+      end do
+
+
+
+   end subroutine check_hlm_active
+
+   ! ------------------------------------------------------------------------------------
+
    subroutine dynamics_driv(this, nc, bounds_clump,      &
          atm2lnd_inst, soilstate_inst, temperature_inst, &
          waterstate_inst, canopystate_inst)