diff --git a/components/clm/src/ED/main/EDCLMLinkMod.F90 b/components/clm/src/ED/main/EDCLMLinkMod.F90 index 009070916e..12703640da 100755 --- a/components/clm/src/ED/main/EDCLMLinkMod.F90 +++ b/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 diff --git a/components/clm/src/main/clm_driver.F90 b/components/clm/src/main/clm_driver.F90 index f968d9643e..576e39d435 100644 --- a/components/clm/src/main/clm_driver.F90 +++ b/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 ! ============================================================================ diff --git a/components/clm/src/main/clm_initializeMod.F90 b/components/clm/src/main/clm_initializeMod.F90 index 7a2274e3ab..e308d14d3f 100644 --- a/components/clm/src/main/clm_initializeMod.F90 +++ b/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. ! ------------------------------------------------------------------------ diff --git a/components/clm/src/main/clm_instMod.F90 b/components/clm/src/main/clm_instMod.F90 index ed22ed6117..7ae1d358b5 100644 --- a/components/clm/src/main/clm_instMod.F90 +++ b/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) diff --git a/components/clm/src/utils/clmfates_interfaceMod.F90 b/components/clm/src/utils/clmfates_interfaceMod.F90 index 9809eb2246..0501cd8f2b 100644 --- a/components/clm/src/utils/clmfates_interfaceMod.F90 +++ b/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,7 +243,16 @@ 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 @@ -209,6 +260,10 @@ subroutine init(this,bounds_proc, use_ed) 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)