Moved contents of ed_driver to the clm_fates interface

Merge branch 'rgknox-eddriver-clean'

ed_driver performed a variety of operations. Along with calling the
main ED dynamics routines, it also performed timing calculations and
the mapping and prepping of FATES/ED state variables into CLM
compliant formats for boundary conditions and IO (vi
clm_ed_inst%ed_clm_link). These commits migrate the latter functions
to live in the clm_fates interface, where CLM calls are appropriate
and the mapping functions and calls are supposed to live.

changes were also added to the udata% timing structure, which was
updated in ed_driver. Two notable changes are the migration of the
n_sub global to be included with the udata%, and removing
udata%cohort_index from regulating the cohort_fusion process. We have
found that udata%cohort_index is not threadsafe (along with the rest
of that structure), and an alternative means using
.not.associated(currentcohort, nextcohort) was implemented instead to
avoid self fusion.

Note: These changes were also built off of PR #59, if those changes
aren't up to snuff, this branch needs to update first.

Fixes: n/a

User interface changes?: No

Code review: rgknox

Test suite: ed - lawrencium-lr2 intel; yellowstone gnu, intel, pgi
Test baseline: 1aaba89
Test namelist changes: none
Test answer changes: b4b

Test summary: pass except for expected fails - #14 f19 f09 restart,
and #43 gnu f10 restart

biogeochem/EDCohortDynamicsMod.F90

@@ -68,8 +68,8 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, &
     !----------------------------------------------------------------------
 
     allocate(new_cohort)
-    udata%cohort_number = udata%cohort_number + 1 !give each cohort a unique number for checking cohort fusing routine.
-    
+    udata%cohort_number = udata%cohort_number + 1  !give each cohort a unique number for checking cohort fusing routine.
+
     call nan_cohort(new_cohort)  ! Make everything in the cohort not-a-number
     call zero_cohort(new_cohort) ! Zero things that need to be zeroed. 
 
@@ -613,7 +613,6 @@ subroutine fuse_cohorts(patchptr)
     iterate = 1
     fusion_took_place = 0   
     currentPatch => patchptr
-   ! maxcohorts =  currentPatch%NCL_p * numCohortsPerPatch
     maxcohorts = numCohortsPerPatch
 
     !---------------------------------------------------------------------!
@@ -624,8 +623,13 @@ subroutine fuse_cohorts(patchptr)
 
          currentCohort => currentPatch%tallest
 
-         !CHANGED FROM C VERSION  loop from tallest to smallest, fusing if they are similar
-         do while (currentCohort%indexnumber /= currentPatch%shortest%indexnumber)  
+         ! The following logic continues the loop while the current cohort is not the shortest cohort
+         ! if they point to the same target (ie equivalence), then the loop ends.
+         ! This loop is different than the simple "continue while associated" loop in that
+         ! it omits the last cohort (because it has already been compared by that point)
+
+         do while ( .not.associated(currentCohort,currentPatch%shortest) )
+
           nextc => currentPatch%tallest
 
           do while (associated(nextc))
@@ -636,7 +640,8 @@ subroutine fuse_cohorts(patchptr)
 
              if (diff < dynamic_fusion_tolerance) then
 
-                if (currentCohort%indexnumber /= nextc%indexnumber) then
+                ! Don't fuse a cohort with itself!
+                if (.not.associated(currentCohort,nextc) ) then
 
                    if (currentCohort%pft == nextc%pft) then              
 
@@ -982,8 +987,8 @@ subroutine copy_cohort( currentCohort,copyc )
     n => copyc
 
     udata%cohort_number = udata%cohort_number + 1
-    n%indexnumber = udata%cohort_number
-
+    n%indexnumber       = udata%cohort_number
+    
     ! VEGETATION STRUCTURE
     n%pft             = o%pft
     n%n               = o%n                         

biogeochem/EDPhysiologyMod.F90

@@ -18,7 +18,7 @@ module EDPhysiologyMod
   use EDCohortDynamicsMod , only : create_cohort, fuse_cohorts, sort_cohorts
   use EDPhenologyType     , only : ed_phenology_type
   use EDTypesMod          , only : dg_sf, dinc_ed, external_recruitment
-  use EDTypesMod          , only : ncwd, nlevcan_ed, n_sub, numpft_ed, senes
+  use EDTypesMod          , only : ncwd, nlevcan_ed, numpft_ed, senes
   use EDTypesMod          , only : ed_site_type, ed_patch_type, ed_cohort_type
 
   implicit none
@@ -250,7 +250,7 @@ subroutine phenology( currentSite, ed_phenology_inst, temperature_inst, watersta
     use EDTypesMod, only : udata
     !
     ! !ARGUMENTS:
-    type(ed_site_type)      , intent(inout), pointer:: currentSite
+    type(ed_site_type)      , intent(inout), target :: currentSite
     type(ed_phenology_type) , intent(in)            :: ed_phenology_inst
     type(temperature_type)  , intent(in)            :: temperature_inst
     type(waterstate_type)   , intent(in)            :: waterstate_inst
@@ -504,7 +504,7 @@ subroutine phenology_leafonoff(currentSite)
     ! !USES:
     !
     ! !ARGUMENTS:
-    type(ed_site_type), intent(inout), pointer:: currentSite
+    type(ed_site_type), intent(inout), target :: currentSite
     !
     ! !LOCAL VARIABLES:
     type(ed_patch_type) , pointer :: currentPatch     
@@ -771,9 +771,9 @@ subroutine Growth_Derivatives( currentCohort)
     ! NPP 
     if ( DEBUG ) write(iulog,*) 'EDphys 716 ',currentCohort%npp_acc
 
-    currentCohort%npp  = currentCohort%npp_acc  * N_SUB   !Link to CLM. convert from kgC/indiv/day into kgC/indiv/year
-    currentCohort%gpp  = currentCohort%gpp_acc  * N_SUB   !Link to CLM. convert from kgC/indiv/day into kgC/indiv/year
-    currentCohort%resp = currentCohort%resp_acc * N_SUB   !Link to CLM. convert from kgC/indiv/day into kgC/indiv/year
+    currentCohort%npp  = currentCohort%npp_acc  * udata%n_sub   !Link to CLM. convert from kgC/indiv/day into kgC/indiv/year
+    currentCohort%gpp  = currentCohort%gpp_acc  * udata%n_sub   !Link to CLM. convert from kgC/indiv/day into kgC/indiv/year
+    currentCohort%resp = currentCohort%resp_acc * udata%n_sub   !Link to CLM. convert from kgC/indiv/day into kgC/indiv/year
 
     currentSite%flux_in = currentSite%flux_in + currentCohort%npp_acc * currentCohort%n
 

main/EDMainMod.F90

@@ -5,7 +5,7 @@ module EDMainMod
   ! ============================================================================
 
   use shr_kind_mod         , only : r8 => shr_kind_r8
-  use spmdMod              , only : masterproc
+
   use decompMod            , only : bounds_type
   use clm_varctl           , only : iulog
   use atm2lndType          , only : atm2lnd_type
@@ -16,7 +16,7 @@ module EDMainMod
   use EDPatchDynamicsMod   , only : disturbance_rates, fuse_patches, spawn_patches, terminate_patches
   use EDPhysiologyMod      , only : canopy_derivs, non_canopy_derivs, phenology, recruitment, trim_canopy
   use SFMainMod            , only : fire_model
-  use EDtypesMod           , only : ncwd, n_sub, numpft_ed, udata
+  use EDtypesMod           , only : ncwd, numpft_ed, udata
   use EDtypesMod           , only : ed_site_type, ed_patch_type, ed_cohort_type
   use EDPhenologyType      , only : ed_phenology_type
   use EDCLMLinkMod         , only : ed_clm_type
@@ -26,11 +26,11 @@ module EDMainMod
 
   !
   ! !PUBLIC MEMBER FUNCTIONS:
-  public  :: ed_driver
+  public  :: ed_ecosystem_dynamics
   public  :: ed_update_site
   !
   ! !PRIVATE MEMBER FUNCTIONS:
-  private :: ed_ecosystem_dynamics
+
   private :: ed_integrate_state_variables
   private :: ed_total_balance_check
 
@@ -41,83 +41,6 @@ module EDMainMod
 
 contains
 
-  !-----------------------------------------------------------------------
-  subroutine ed_driver( bounds, ed_allsites_inst, ed_clm_inst, ed_phenology_inst, &
-       atm2lnd_inst, soilstate_inst, temperature_inst, waterstate_inst, canopystate_inst)
-    !
-    ! !DESCRIPTION:
-    ! Main ed model routine containing gridcell loop   
-    !
-    ! !USES:
-    use clm_time_manager     , only : get_days_per_year, get_curr_date
-    use clm_time_manager     , only : get_ref_date, timemgr_datediff 
-    use CanopySTateType      , only : canopystate_type
-    !
-    ! !ARGUMENTS:
-    type(bounds_type)       , intent(in)            :: bounds           
-    type(ed_site_type)      , intent(inout), target :: ed_allsites_inst( bounds%begg: )
-    type(ed_clm_type)       , intent(inout)         :: ed_clm_inst
-    type(ed_phenology_type) , intent(inout)         :: ed_phenology_inst
-    type(atm2lnd_type)      , intent(in)            :: atm2lnd_inst
-    type(soilstate_type)    , intent(in)            :: soilstate_inst
-    type(temperature_type)  , intent(in)            :: temperature_inst
-    type(waterstate_type)   , intent(inout)         :: waterstate_inst
-    type(canopystate_type)  , intent(inout)         :: canopystate_inst
-    !
-    ! !LOCAL VARIABLES:
-    type(ed_site_type), pointer :: currentSite
-    real(r8) :: dayDiff                  ! day of run
-    integer  :: dayDiffInt               ! integer of day of run
-    integer  :: g                        ! gridcell  
-    integer  :: yr                       ! year (0, ...)
-    integer  :: mon                      ! month (1, ..., 12)
-    integer  :: day                      ! day of month (1, ..., 31)
-    integer  :: sec                      ! seconds of the day
-    integer  :: ncdate                   ! current date
-    integer  :: nbdate                   ! base date (reference date)
-    !-----------------------------------------------------------------------
-
-    call ed_clm_inst%SetValues( bounds, 0._r8 )
-
-    ! timing statements. 
-    n_sub = get_days_per_year()
-    udata%deltat = 1.0_r8/n_sub !for working out age of patches in years        
-    if(udata%time_period == 0)then             
-       udata%time_period = n_sub
-    endif
-
-    call get_curr_date(yr, mon, day, sec)
-    ncdate = yr*10000 + mon*100 + day
-    call get_ref_date(yr, mon, day, sec)
-    nbdate = yr*10000 + mon*100 + day
-
-    call timemgr_datediff(nbdate, 0, ncdate, sec, dayDiff)
-
-    dayDiffInt = floor(dayDiff)
-    udata%time_period = mod( dayDiffInt , n_sub )
-
-    ! where most things happen
-    do g = bounds%begg,bounds%endg
-       if (ed_allsites_inst(g)%istheresoil) then
-          currentSite => ed_allsites_inst(g)
-          call ed_ecosystem_dynamics(currentSite, &
-               ed_clm_inst, ed_phenology_inst, atm2lnd_inst, &
-               soilstate_inst, temperature_inst, waterstate_inst)
-
-          call ed_update_site( ed_allsites_inst(g))
-       endif
-    enddo
-
-    ! link to CLM structures
-    call ed_clm_inst%ed_clm_link( bounds, ed_allsites_inst(bounds%begg:bounds%endg),  &
-         ed_phenology_inst, waterstate_inst, canopystate_inst)
-
-    if (masterproc) then
-      write(iulog, *) 'clm: leaving ED model', bounds%begg, bounds%endg, dayDiffInt
-    end if
-
-  end subroutine ed_driver
-
   !-------------------------------------------------------------------------------!
   subroutine ed_ecosystem_dynamics(currentSite, &
        ed_clm_inst, ed_phenology_inst, atm2lnd_inst, &
@@ -127,7 +50,7 @@ subroutine ed_ecosystem_dynamics(currentSite, &
     !  Core of ed model, calling all subsequent vegetation dynamics routines         
     !
     ! !ARGUMENTS:
-    type(ed_site_type)      , intent(inout), pointer :: currentSite
+    type(ed_site_type)      , intent(inout), target  :: currentSite
     type(ed_phenology_type) , intent(in)             :: ed_phenology_inst
     type(ed_clm_type)       , intent(in)             :: ed_clm_inst
     type(atm2lnd_type)      , intent(in)             :: atm2lnd_inst
@@ -221,10 +144,10 @@ subroutine ed_integrate_state_variables(currentSite, soilstate_inst, temperature
     ! !USES:
     !
     ! !ARGUMENTS:
-    type(ed_site_type)     , intent(in) :: currentSite
-    type(soilstate_type)   , intent(in) :: soilstate_inst
-    type(temperature_type) , intent(in) :: temperature_inst
-    type(waterstate_type)  , intent(in) :: waterstate_inst
+    type(ed_site_type)     , intent(in)    :: currentSite
+    type(soilstate_type)   , intent(in)    :: soilstate_inst
+    type(temperature_type) , intent(in)    :: temperature_inst
+    type(waterstate_type)  , intent(in)    :: waterstate_inst
     !
     ! !LOCAL VARIABLES:
     type(ed_patch_type)  , pointer :: currentPatch
@@ -418,7 +341,7 @@ subroutine ed_update_site( currentSite )
     enddo
 
     ! FIX(RF,032414). This needs to be monthly, not annual
-    if((udata%time_period == N_SUB-1))then 
+    if((udata%time_period == udata%n_sub-1))then 
        write(iulog,*) 'calling trim canopy' 
        call trim_canopy(currentSite)  
     endif

main/EDRestVectorMod.F90

@@ -634,7 +634,6 @@ subroutine getVectors( this, bounds, ed_allsites_inst )
     ! !USES:
     use clm_time_manager , only : get_nstep
     use EDCLMLinkMod     , only : ed_clm_type
-    use EDInitMod        , only : ed_init_sites
     use EDMainMod        , only : ed_update_site
     !
     ! !ARGUMENTS:
@@ -662,9 +661,6 @@ subroutine getVectors( this, bounds, ed_allsites_inst )
        end if
     end do
 
-!    call ed_clm_inst%ed_clm_link( bounds, ed_allsites_inst(bounds%begg:bounds%endg), &
-!         ed_phenology_inst, waterstate_inst, canopystate_inst)
-
     if (this%DEBUG) then
        call this%printIoInfoLL ( bounds, ed_allsites_inst(bounds%begg:bounds%endg) )
        call this%printDataInfoLL ( bounds, ed_allsites_inst(bounds%begg:bounds%endg) )

main/EDTypesMod.F90

@@ -14,7 +14,7 @@ module EDTypesMod
 
   ! MODEL PARAMETERS
   real(r8)            :: timestep_secs                     ! subdaily timestep in seconds (e.g. 1800 or 3600) 
-  integer             :: n_sub                             ! num of substeps in year 
+
   real(r8), parameter :: AREA                 = 10000.0_r8 ! Notional area of simulated forest m2
   integer  doy
 
@@ -438,13 +438,15 @@ module EDTypesMod
   !************************************
 
   type userdata
-     integer  ::   cohort_number            ! Counts up the number of cohorts which have been made. 
+     integer  ::   cohort_number            ! Counts up the number of cohorts which have been made.
+     integer  ::   n_sub                    ! num of substeps in year 
      real(r8) ::   deltat                   ! fraction of year used for each timestep (1/N_SUB)
      integer  ::   time_period              ! Within year timestep (1:N_SUB) day of year
      integer  ::   restart_year             ! Which year of simulation are we starting in? 
   end type userdata
 
-  type(userdata), public, target :: udata
+
+  type(userdata), public, target :: udata   ! THIS WAS NOT THREADSAFE
   !-------------------------------------------------------------------------------------!
 
   public :: ed_hist_scpfmaps