From 1a943a035bd881938bdd1b7c447a2a11dadb7dd6 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 6 Jan 2017 16:05:00 -0800 Subject: [PATCH 01/11] Initial changes to timing boundary conditions, mostly directed towards the calculation of model-day for phenology. --- .../clm/src/ED/biogeochem/EDPhysiologyMod.F90 | 39 ++++++++---------- components/clm/src/ED/main/EDMainMod.F90 | 9 ++++- components/clm/src/ED/main/EDTypesMod.F90 | 1 + .../clm/src/ED/main/FatesInterfaceMod.F90 | 19 +++++++++ .../clm/src/utils/clmfates_interfaceMod.F90 | 40 +++++++++++++++++-- 5 files changed, 81 insertions(+), 27 deletions(-) diff --git a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 index fccd8c0843..7ff0274efd 100755 --- a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 +++ b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 @@ -8,7 +8,7 @@ module EDPhysiologyMod use shr_kind_mod , only : r8 => shr_kind_r8 use clm_varctl , only : iulog - use spmdMod , only : masterproc + use TemperatureType , only : temperature_type use SoilStateType , only : soilstate_type use WaterstateType , only : waterstate_type @@ -240,20 +240,21 @@ subroutine trim_canopy( currentSite ) end subroutine trim_canopy ! ============================================================================ - subroutine phenology( currentSite, temperature_inst, waterstate_inst) + subroutine phenology( currentSite, bc_in, temperature_inst, waterstate_inst) ! ! !DESCRIPTION: ! Phenology. ! ! !USES: use clm_varcon, only : tfrz - use clm_time_manager, only : get_curr_date - use clm_time_manager, only : get_ref_date, timemgr_datediff + use FatesInterfaceMod, only : bc_in_type use EDTypesMod, only : udata - use PatchType , only : patch + use PatchType , only : patch ! ! !ARGUMENTS: - type(ed_site_type) , intent(inout), target :: currentSite + type(ed_site_type), intent(inout), target :: currentSite + type(bc_in_type), intent(in) :: bc_in + type(temperature_type) , intent(in) :: temperature_inst type(waterstate_type) , intent(in) :: waterstate_inst ! @@ -283,8 +284,9 @@ subroutine phenology( currentSite, temperature_inst, waterstate_inst) real(r8) :: drought_threshold real(r8) :: off_time ! minimum number of days between leaf off and leaf on for drought phenology real(r8) :: temp_in_C ! daily averaged temperature in celcius - real(r8) :: mindayson - real(r8) :: modelday + + real(r8), parameter :: mindayson = 30.0 + !------------------------------------------------------------------------ @@ -294,16 +296,9 @@ subroutine phenology( currentSite, temperature_inst, waterstate_inst) patchi = currentSite%oldest_patch%clm_pno-1 coli = patch%column(patchi) - t_veg24 => temperature_inst%t_veg24_patch ! Input: [real(r8) (:)] avg pft vegetation temperature for last 24 hrs + t_veg24 => temperature_inst%t_veg24_patch ! Input: [real(r8) (:)] avg pft vegetation temperature for last 24 hrs - call get_curr_date(yr, mon, day, sec) - curdate = yr*10000 + mon*100 + day - call get_ref_date(yr, mon, day, sec) - refdate = yr*10000 + mon*100 + day - - call timemgr_datediff(refdate, 0, curdate, sec, modelday) - if ( masterproc ) write(iulog,*) 'modelday',modelday ! Parameter of drought decid leaf loss in mm in top layer...FIX(RF,032414) ! - this is arbitrary and poorly understood. Needs work. ED_ @@ -316,7 +311,7 @@ subroutine phenology( currentSite, temperature_inst, waterstate_inst) c = -0.001_r8 coldday = 5.0_r8 !ed_ph_chiltemp - mindayson = 30 + !Parameters from SDGVM model of senesence ncolddayslim = 5 @@ -372,7 +367,7 @@ subroutine phenology( currentSite, temperature_inst, waterstate_inst) endif - timesinceleafoff = modelday - currentSite%leafoffdate + timesinceleafoff = bc_in%model_day - currentSite%leafoffdate !LEAF ON: COLD DECIDUOUS. Needs to !1) have exceeded the growing degree day threshold !2) The leaves should not be on already @@ -381,14 +376,14 @@ subroutine phenology( currentSite, temperature_inst, waterstate_inst) if (currentSite%status == 1) then if (currentSite%ncd >= 1) then currentSite%status = 2 !alter status of site to 'leaves on' - ! NOTE(bja, 2015-01) should leafondate = modelday to be consistent with leaf off? + ! NOTE(bja, 2015-01) should leafondate = model_day to be consistent with leaf off? currentSite%leafondate = t !record leaf on date if ( DEBUG ) write(iulog,*) 'leaves on' endif !ncd endif !status endif !GDD - timesinceleafon = modelday - currentSite%leafondate + timesinceleafon = bc_in%model_day - currentSite%leafondate !LEAF OFF: COLD THRESHOLD @@ -402,7 +397,7 @@ subroutine phenology( currentSite, temperature_inst, waterstate_inst) if (timesinceleafon > mindayson)then if (currentSite%status == 2)then currentSite%status = 1 !alter status of site to 'leaves on' - currentSite%leafoffdate = modelday !record leaf off date + currentSite%leafoffdate = bc_in%model_day !record leaf off date if ( DEBUG ) write(iulog,*) 'leaves off' endif endif @@ -412,7 +407,7 @@ subroutine phenology( currentSite, temperature_inst, waterstate_inst) if(timesinceleafoff > 400)then !remove leaves after a whole year when there is no 'off' period. if(currentSite%status == 2)then currentSite%status = 1 !alter status of site to 'leaves on' - currentSite%leafoffdate = modelday !record leaf off date + currentSite%leafoffdate = bc_in%model_day !record leaf off date if ( DEBUG ) write(iulog,*) 'leaves off' endif endif diff --git a/components/clm/src/ED/main/EDMainMod.F90 b/components/clm/src/ED/main/EDMainMod.F90 index 9499f93d02..f496fb4044 100755 --- a/components/clm/src/ED/main/EDMainMod.F90 +++ b/components/clm/src/ED/main/EDMainMod.F90 @@ -17,6 +17,8 @@ module EDMainMod use SFMainMod , only : fire_model use EDtypesMod , only : ncwd, numpft_ed, udata use EDtypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type + use FatesInterfaceMod , only : bc_in_type + use spmdMod , only : masterproc implicit none private @@ -39,7 +41,7 @@ module EDMainMod contains !-------------------------------------------------------------------------------! - subroutine ed_ecosystem_dynamics(currentSite, & + subroutine ed_ecosystem_dynamics(currentSite, bc_in, & atm2lnd_inst, & soilstate_inst, temperature_inst, waterstate_inst) ! @@ -48,6 +50,7 @@ subroutine ed_ecosystem_dynamics(currentSite, & ! ! !ARGUMENTS: type(ed_site_type) , intent(inout), target :: currentSite + type(bc_in_type) , intent(in) :: bc_in type(atm2lnd_type) , intent(in) :: atm2lnd_inst type(soilstate_type) , intent(in) :: soilstate_inst type(temperature_type) , intent(in) :: temperature_inst @@ -57,6 +60,8 @@ subroutine ed_ecosystem_dynamics(currentSite, & type(ed_patch_type), pointer :: currentPatch !----------------------------------------------------------------------- + if ( masterproc ) write(iulog,*) 'modelday',bc_in%model_day + !************************************************************************** ! Fire, growth, biogeochemistry. !************************************************************************** @@ -66,7 +71,7 @@ subroutine ed_ecosystem_dynamics(currentSite, & call ed_total_balance_check(currentSite, 0) - call phenology(currentSite, temperature_inst, waterstate_inst) + call phenology(currentSite, bc_in, temperature_inst, waterstate_inst) call fire_model(currentSite, atm2lnd_inst, temperature_inst) diff --git a/components/clm/src/ED/main/EDTypesMod.F90 b/components/clm/src/ED/main/EDTypesMod.F90 index 6de2f1ea2c..4d80cf314c 100755 --- a/components/clm/src/ED/main/EDTypesMod.F90 +++ b/components/clm/src/ED/main/EDTypesMod.F90 @@ -551,6 +551,7 @@ module EDTypesMod 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? + integer :: modelday ! Number of days since reference end type userdata diff --git a/components/clm/src/ED/main/FatesInterfaceMod.F90 b/components/clm/src/ED/main/FatesInterfaceMod.F90 index 356951bcf0..46f1268aa2 100644 --- a/components/clm/src/ED/main/FatesInterfaceMod.F90 +++ b/components/clm/src/ED/main/FatesInterfaceMod.F90 @@ -50,6 +50,18 @@ module FatesInterfaceMod ! The actual number of FATES' ED patches integer :: npatches + ! Timing Variables + integer :: current_year ! Current year + integer :: current_month ! month of year + integer :: current_day ! day of month + integer :: current_tod ! time of day (seconds past 0Z) + integer :: current_date ! time of day (seconds past 0Z) + integer :: reference_date ! YYYYMMDD + real(r8) :: model_day ! elapsed days between current date and reference + ! uses ESMF functions, so prefered to pass it in as + ! argument rather than calculate directly + + ! Radiation variables for calculating sun/shade fractions ! --------------------------------------------------------------------------------- @@ -447,6 +459,13 @@ subroutine zero_bcs(this,s) integer, intent(in) :: s ! Input boundaries + this%bc_in(s)%current_year = 0 + this%bc_in(s)%current_month = 0 + this%bc_in(s)%current_day = 0 + this%bc_in(s)%current_tod = 0 + this%bc_in(s)%current_date = 0 + this%bc_in(s)%reference_date = 0 + this%bc_in(s)%model_day = 0.0_r8 this%bc_in(s)%solad_parb(:,:) = 0.0_r8 this%bc_in(s)%solai_parb(:,:) = 0.0_r8 diff --git a/components/clm/src/utils/clmfates_interfaceMod.F90 b/components/clm/src/utils/clmfates_interfaceMod.F90 index 963d013e37..1432277d0b 100644 --- a/components/clm/src/utils/clmfates_interfaceMod.F90 +++ b/components/clm/src/utils/clmfates_interfaceMod.F90 @@ -447,6 +447,13 @@ subroutine dynamics_driv(this, nc, bounds_clump, & integer :: sec ! seconds of the day integer :: ncdate ! current date integer :: nbdate ! base date (reference date) + integer :: current_year + integer :: current_month + integer :: current_day + integer :: current_tod + integer :: current_date + integer :: reference_date + real(r8) :: model_day !----------------------------------------------------------------------- @@ -463,7 +470,7 @@ subroutine dynamics_driv(this, nc, bounds_clump, & ! timing statements. udata%n_sub = get_days_per_year() - udata%deltat = 1.0_r8/dble(udata%n_sub) !for working out age of patches in years + udata%deltat = 1.0_r8/dble(udata%n_sub) !for working out age of patches in years if(udata%time_period == 0)then udata%time_period = udata%n_sub endif @@ -474,12 +481,38 @@ subroutine dynamics_driv(this, nc, bounds_clump, & nbdate = yr*10000 + mon*100 + day call timemgr_datediff(nbdate, 0, ncdate, sec, dayDiff) - + udata%modelday = dayDiff dayDiffInt = floor(dayDiff) udata%time_period = mod( dayDiffInt , udata%n_sub ) + ! --------------------------------------------------------------------------------- + ! Prepare input boundary conditions for FATES dynamics + ! Note that timing information is the same across all sites, this may + ! seem redundant, but it is possible that we may have asynchronous site simulations + ! one day. The cost of holding site level boundary conditions is minimal + ! and it keeps all the boundaries in one location + ! --------------------------------------------------------------------------------- + call get_curr_date(current_year,current_month,current_day,current_tod) + current_date = current_year*10000 + current_month*100 + current_day + + call get_ref_date(yr, mon, day, sec) + reference_date = yr*10000 + mon*100 + day + + call timemgr_datediff(reference_date, sec, current_date, current_tod, model_day) + if ( masterproc ) write(iulog,*) 'modelday',model_day + + do s=1,this%fates(nc)%nsites + this%fates(nc)%bc_in(s)%current_year = current_year + this%fates(nc)%bc_in(s)%current_month = current_month + this%fates(nc)%bc_in(s)%current_day = current_day + this%fates(nc)%bc_in(s)%current_tod = current_tod + this%fates(nc)%bc_in(s)%current_tod = current_date + this%fates(nc)%bc_in(s)%reference_date = reference_date + this%fates(nc)%bc_in(s)%model_day = model_day + end do + - ! TODO-INTEF: PROCEDURE FOR CONVERTING CLM/ALM FIELDS TO MODEL BOUNDARY + ! TODO-INTEF: PROCEDURE FOR CONVERTING CLM/ALM FIELDS TO MODEL BOUNDARY ! CONDITIONS. IE. @@ -487,6 +520,7 @@ subroutine dynamics_driv(this, nc, bounds_clump, & do s = 1,this%fates(nc)%nsites call ed_ecosystem_dynamics(this%fates(nc)%sites(s), & + this%fates(nc)%bc_in(s), & atm2lnd_inst, & soilstate_inst, temperature_inst, waterstate_inst) From f5fb7126bcfb39ea72f2e28cca37f7a78489a4ba Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Sun, 8 Jan 2017 21:39:55 -0800 Subject: [PATCH 02/11] Converted 24 vegetation temperatures used in phenology to bc_in --- .../clm/src/ED/biogeochem/EDPhysiologyMod.F90 | 28 +++++-------------- components/clm/src/ED/main/EDMainMod.F90 | 2 +- .../clm/src/ED/main/FatesInterfaceMod.F90 | 21 ++++++++++++++ .../clm/src/utils/clmfates_interfaceMod.F90 | 15 ++++++++-- 4 files changed, 42 insertions(+), 24 deletions(-) diff --git a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 index 7ff0274efd..02eedb26e1 100755 --- a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 +++ b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 @@ -240,26 +240,25 @@ subroutine trim_canopy( currentSite ) end subroutine trim_canopy ! ============================================================================ - subroutine phenology( currentSite, bc_in, temperature_inst, waterstate_inst) + subroutine phenology( currentSite, bc_in, waterstate_inst) ! ! !DESCRIPTION: ! Phenology. ! ! !USES: - use clm_varcon, only : tfrz + use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm use FatesInterfaceMod, only : bc_in_type use EDTypesMod, only : udata - use PatchType , only : patch + ! ! !ARGUMENTS: type(ed_site_type), intent(inout), target :: currentSite type(bc_in_type), intent(in) :: bc_in - type(temperature_type) , intent(in) :: temperature_inst type(waterstate_type) , intent(in) :: waterstate_inst ! ! !LOCAL VARIABLES: - real(r8), pointer :: t_veg24(:) + integer :: t ! day of year integer :: ncolddays ! no days underneath the threshold for leaf drop integer :: ncolddayslim ! critical no days underneath the threshold for leaf drop @@ -287,19 +286,6 @@ subroutine phenology( currentSite, bc_in, temperature_inst, waterstate_inst) real(r8), parameter :: mindayson = 30.0 - - !------------------------------------------------------------------------ - - ! INTERF-TODO: THIS IS A BAND-AID, AS I WAS HOPING TO REMOVE CLM_PNO - ! ALREADY REMOVED currentSite%clmcolumn, hence the need for these - - patchi = currentSite%oldest_patch%clm_pno-1 - coli = patch%column(patchi) - - t_veg24 => temperature_inst%t_veg24_patch ! Input: [real(r8) (:)] avg pft vegetation temperature for last 24 hrs - - - ! Parameter of drought decid leaf loss in mm in top layer...FIX(RF,032414) ! - this is arbitrary and poorly understood. Needs work. ED_ drought_threshold = 0.15 @@ -318,7 +304,7 @@ subroutine phenology( currentSite, bc_in, temperature_inst, waterstate_inst) cold_t = 7.5_r8 ! ed_ph_coldtemp t = udata%time_period - temp_in_C = t_veg24(patchi) - tfrz + temp_in_C = bc_in%t_veg24_si - tfrz !-----------------Cold Phenology--------------------! @@ -362,8 +348,8 @@ subroutine phenology( currentSite, bc_in, temperature_inst, waterstate_inst) endif ! ! accumulate the GDD using daily mean temperatures - if (t_veg24(patchi) .gt. tfrz) then - currentSite%ED_GDD_site = currentSite%ED_GDD_site + t_veg24(currentSite%oldest_patch%clm_pno-1) - tfrz + if (bc_in%t_veg24_si .gt. tfrz) then + currentSite%ED_GDD_site = currentSite%ED_GDD_site + bc_in%t_veg24_si - tfrz endif diff --git a/components/clm/src/ED/main/EDMainMod.F90 b/components/clm/src/ED/main/EDMainMod.F90 index f496fb4044..ffd7948d31 100755 --- a/components/clm/src/ED/main/EDMainMod.F90 +++ b/components/clm/src/ED/main/EDMainMod.F90 @@ -71,7 +71,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, & call ed_total_balance_check(currentSite, 0) - call phenology(currentSite, bc_in, temperature_inst, waterstate_inst) + call phenology(currentSite, bc_in, waterstate_inst ) call fire_model(currentSite, atm2lnd_inst, temperature_inst) diff --git a/components/clm/src/ED/main/FatesInterfaceMod.F90 b/components/clm/src/ED/main/FatesInterfaceMod.F90 index 46f1268aa2..e61a6dda5d 100644 --- a/components/clm/src/ED/main/FatesInterfaceMod.F90 +++ b/components/clm/src/ED/main/FatesInterfaceMod.F90 @@ -61,6 +61,20 @@ module FatesInterfaceMod ! uses ESMF functions, so prefered to pass it in as ! argument rather than calculate directly + ! Vegetation Dynamics + ! --------------------------------------------------------------------------------- + + ! This 24 hour vegetation temperature is used for various purposes during vegetation + ! dynamics. However, we are currently using the bare-ground patch's value [K] + ! TO-DO: Get some consensus on the correct vegetation temperature used for phenology. + ! It is possible that the bare-ground value is where the average is being stored. + ! (RGK-01-2017) + + real(r8) :: t_veg24_si ! See above [K] + + ! Patch 24 hour vegetation temperature [K] + real(r8),allocatable :: t_veg24_pa(:) + ! Radiation variables for calculating sun/shade fractions ! --------------------------------------------------------------------------------- @@ -357,6 +371,10 @@ subroutine allocate_bcin(bc_in) ! Allocate input boundaries + ! Vegetation Dynamics + allocate(bc_in%t_veg24_pa(numPatchesPerCol)) + + ! Radiation allocate(bc_in%solad_parb(numPatchesPerCol,cp_numSWb)) allocate(bc_in%solai_parb(numPatchesPerCol,cp_numSWb)) @@ -467,6 +485,9 @@ subroutine zero_bcs(this,s) this%bc_in(s)%reference_date = 0 this%bc_in(s)%model_day = 0.0_r8 + this%bc_in(s)%t_veg24_pa(:) = 0.0_r8 + this%bc_in(s)%t_veg24_si = 0.0_r8 + this%bc_in(s)%solad_parb(:,:) = 0.0_r8 this%bc_in(s)%solai_parb(:,:) = 0.0_r8 this%bc_in(s)%smp_gl(:) = 0.0_r8 diff --git a/components/clm/src/utils/clmfates_interfaceMod.F90 b/components/clm/src/utils/clmfates_interfaceMod.F90 index 1432277d0b..747427312b 100644 --- a/components/clm/src/utils/clmfates_interfaceMod.F90 +++ b/components/clm/src/utils/clmfates_interfaceMod.F90 @@ -440,7 +440,10 @@ subroutine dynamics_driv(this, nc, bounds_clump, & ! !LOCAL VARIABLES: real(r8) :: dayDiff ! day of run integer :: dayDiffInt ! integer of day of run - integer :: s ! site + integer :: s ! site index + integer :: c ! column index (HLM) + integer :: ifp ! patch index + integer :: p ! HLM patch index integer :: yr ! year (0, ...) integer :: mon ! month (1, ..., 12) integer :: day ! day of month (1, ..., 31) @@ -502,6 +505,7 @@ subroutine dynamics_driv(this, nc, bounds_clump, & if ( masterproc ) write(iulog,*) 'modelday',model_day do s=1,this%fates(nc)%nsites + c = this%f2hmap(nc)%fcolumn(s) this%fates(nc)%bc_in(s)%current_year = current_year this%fates(nc)%bc_in(s)%current_month = current_month this%fates(nc)%bc_in(s)%current_day = current_day @@ -509,10 +513,17 @@ subroutine dynamics_driv(this, nc, bounds_clump, & this%fates(nc)%bc_in(s)%current_tod = current_date this%fates(nc)%bc_in(s)%reference_date = reference_date this%fates(nc)%bc_in(s)%model_day = model_day + this%fates(nc)%bc_in(s)%t_veg24_si = & + temperature_inst%t_veg24_patch(col%patchi(c)-1) + do ifp = 1, this%fates(nc)%sites(s)%youngest_patch%patchno + p = ifp+col%patchi(c) + this%fates(nc)%bc_in(s)%t_veg24_pa(ifp) = & + temperature_inst%t_veg24_patch(p) + end do end do - ! TODO-INTEF: PROCEDURE FOR CONVERTING CLM/ALM FIELDS TO MODEL BOUNDARY + ! TODO-INTER: PROCEDURE FOR CONVERTING CLM/ALM FIELDS TO MODEL BOUNDARY ! CONDITIONS. IE. From 96e680e920d2a19c01eb95796f4217d199902488 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 10 Jan 2017 17:23:35 -0800 Subject: [PATCH 03/11] Purging waterstate_inst, soilstate_inst from calls to dynamics and replacing it with bc_in. Temperature_inst for all calls but fire are also purged. Compiles and runs, not regression tested. --- .../cesm/machines/config_compilers.xml | 4 +- .../clm/src/ED/biogeochem/EDPhysiologyMod.F90 | 58 ++++++++----------- components/clm/src/ED/main/EDMainMod.F90 | 18 +++--- .../clm/src/ED/main/FatesInterfaceMod.F90 | 18 +++--- .../clm/src/utils/clmfates_interfaceMod.F90 | 9 ++- 5 files changed, 52 insertions(+), 55 deletions(-) diff --git a/cime/cime_config/cesm/machines/config_compilers.xml b/cime/cime_config/cesm/machines/config_compilers.xml index 318ac8ed9a..ee13c73857 100644 --- a/cime/cime_config/cesm/machines/config_compilers.xml +++ b/cime/cime_config/cesm/machines/config_compilers.xml @@ -556,8 +556,8 @@ for mct, etc. $(NETCDF_HOME) -DLinux -DCPRGNU - -ffree-line-length-132 - -ffpe-trap=invalid,zero,overflow + -ffree-line-length-132 + -g -fbacktrace -fbounds-check -ffpe-trap=invalid,zero,overflow -Wline-truncation -L$(NETCDF_HOME)/lib/ -lnetcdff -lnetcdf -lcurl -llapack -lblas -DHAVE_VPRINTF -DHAVE_GETTIMEOFDAY diff --git a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 index 02eedb26e1..9d84507865 100755 --- a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 +++ b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 @@ -10,10 +10,9 @@ module EDPhysiologyMod use clm_varctl , only : iulog use TemperatureType , only : temperature_type - use SoilStateType , only : soilstate_type - use WaterstateType , only : waterstate_type use pftconMod , only : pftcon use EDEcophysContype , only : EDecophyscon + use FatesInterfaceMod, only : bc_in_type use EDCohortDynamicsMod , only : allocate_live_biomass, zero_cohort use EDCohortDynamicsMod , only : create_cohort, fuse_cohorts, sort_cohorts use EDTypesMod , only : dg_sf, dinc_ed, external_recruitment @@ -72,7 +71,7 @@ subroutine canopy_derivs( currentSite, currentPatch ) end subroutine canopy_derivs ! ============================================================================ - subroutine non_canopy_derivs( currentSite, currentPatch, temperature_inst ) + subroutine non_canopy_derivs( currentSite, currentPatch, bc_in ) ! ! !DESCRIPTION: ! Returns time differentials of the state vector @@ -82,8 +81,9 @@ subroutine non_canopy_derivs( currentSite, currentPatch, temperature_inst ) ! ! !ARGUMENTS type(ed_site_type), intent(inout), target :: currentSite - type(ed_patch_type) , intent(inout) :: currentPatch - type(temperature_type) , intent(in) :: temperature_inst + type(ed_patch_type), intent(inout) :: currentPatch + type(bc_in_type), intent(in) :: bc_in + ! ! !LOCAL VARIABLES: integer c,p @@ -110,7 +110,7 @@ subroutine non_canopy_derivs( currentSite, currentPatch, temperature_inst ) ! update fragmenting pool fluxes call cwd_input(currentPatch) - call cwd_out( currentSite, currentPatch, temperature_inst) + call cwd_out( currentSite, currentPatch, bc_in) do p = 1,numpft_ed currentSite%dseed_dt(p) = currentSite%dseed_dt(p) + & @@ -240,14 +240,13 @@ subroutine trim_canopy( currentSite ) end subroutine trim_canopy ! ============================================================================ - subroutine phenology( currentSite, bc_in, waterstate_inst) + subroutine phenology( currentSite, bc_in ) ! ! !DESCRIPTION: ! Phenology. ! ! !USES: use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm - use FatesInterfaceMod, only : bc_in_type use EDTypesMod, only : udata ! @@ -255,7 +254,6 @@ subroutine phenology( currentSite, bc_in, waterstate_inst) type(ed_site_type), intent(inout), target :: currentSite type(bc_in_type), intent(in) :: bc_in - type(waterstate_type) , intent(in) :: waterstate_inst ! ! !LOCAL VARIABLES: @@ -271,8 +269,6 @@ subroutine phenology( currentSite, bc_in, waterstate_inst) integer :: mon ! month (1, ..., 12) integer :: day ! day of month (1, ..., 31) integer :: sec ! seconds of the day - integer :: patchi ! the first CLM/ALM patch index of the associated column - integer :: coli ! the CLM/ALM column index of the associated site real(r8) :: gdd_threshold real(r8) :: a,b,c ! params of leaf-pn model from botta et al. 2000. @@ -283,7 +279,6 @@ subroutine phenology( currentSite, bc_in, waterstate_inst) real(r8) :: drought_threshold real(r8) :: off_time ! minimum number of days between leaf off and leaf on for drought phenology real(r8) :: temp_in_C ! daily averaged temperature in celcius - real(r8), parameter :: mindayson = 30.0 ! Parameter of drought decid leaf loss in mm in top layer...FIX(RF,032414) @@ -296,8 +291,6 @@ subroutine phenology( currentSite, bc_in, waterstate_inst) b = 638.0_r8 c = -0.001_r8 coldday = 5.0_r8 !ed_ph_chiltemp - - !Parameters from SDGVM model of senesence ncolddayslim = 5 @@ -426,7 +419,7 @@ subroutine phenology( currentSite, bc_in, waterstate_inst) ! distinction actually matter??).... !Accumulate surface water memory of last 10 days. - currentSite%water_memory(1) = waterstate_inst%h2osoi_vol_col(coli,1) + currentSite%water_memory(1) = bc_in%h2osoi_vol_si !waterstate_inst%h2osoi_vol_col(coli,1) do i = 1,9 !shift memory along one currentSite%water_memory(11-i) = currentSite%water_memory(10-i) enddo @@ -1121,7 +1114,7 @@ subroutine CWD_Input( currentPatch) end subroutine CWD_Input ! ============================================================================ - subroutine fragmentation_scaler( currentPatch, temperature_inst ) + subroutine fragmentation_scaler( currentPatch, bc_in) ! ! !DESCRIPTION: ! Simple CWD fragmentation Model @@ -1133,12 +1126,14 @@ subroutine fragmentation_scaler( currentPatch, temperature_inst ) ! ! !ARGUMENTS - type(ed_patch_type) , intent(inout) :: currentPatch - type(temperature_type) , intent(in) :: temperature_inst + type(ed_patch_type), intent(inout) :: currentPatch + type(bc_in_type), intent(in) :: bc_in + ! ! !LOCAL VARIABLES: logical :: use_century_tfunc = .false. - integer :: p,j + integer :: j + integer :: ifp ! Index of a FATES Patch "ifp" real(r8) :: t_scalar real(r8) :: w_scalar real(r8) :: catanf ! hyperbolic temperature function from CENTURY @@ -1146,16 +1141,12 @@ subroutine fragmentation_scaler( currentPatch, temperature_inst ) real(r8) :: t1 ! temperature argument real(r8) :: Q10 ! temperature dependence real(r8) :: froz_q10 ! separate q10 for frozen soil respiration rates. default to same as above zero rates - real(r8), pointer :: t_veg24(:) !---------------------------------------------------------------------- catanf(t1) = 11.75_r8 +(29.7_r8 / SHR_CONST_PI) * atan( SHR_CONST_PI * 0.031_r8 * ( t1 - 15.4_r8 )) - - t_veg24 => temperature_inst%t_veg24_patch ! Input: [real(r8) (:)] avg pft vegetation temperature for last 24 hrs - catanf_30 = catanf(30._r8) - p = currentPatch%clm_pno + ifp = currentPatch%patchno ! set "froz_q10" parameter froz_q10 = EDParamsShareInst%froz_q10 @@ -1164,16 +1155,16 @@ subroutine fragmentation_scaler( currentPatch, temperature_inst ) 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 (t_veg24(p) >= SHR_CONST_TKFRZ) then - t_scalar = Q10**((t_veg24(p)-(SHR_CONST_TKFRZ+25._r8))/10._r8) + if (bc_in%t_veg24_pa(ifp) >= SHR_CONST_TKFRZ) then + t_scalar = Q10**((bc_in%t_veg24_pa(ifp)-(SHR_CONST_TKFRZ+25._r8))/10._r8) ! Q10**((t_soisno(c,j)-(SHR_CONST_TKFRZ+25._r8))/10._r8) else - t_scalar = (Q10**(-25._r8/10._r8))*(froz_q10**((t_veg24(p)-SHR_CONST_TKFRZ)/10._r8)) + t_scalar = (Q10**(-25._r8/10._r8))*(froz_q10**((bc_in%t_veg24_pa(ifp)-SHR_CONST_TKFRZ)/10._r8)) !Q10**(-25._r8/10._r8))*(froz_q10**((t_soisno(c,j)-SHR_CONST_TKFRZ)/10._r8) endif else ! original century uses an arctangent function to calculate the temperature dependence of decomposition - t_scalar = max(catanf(t_veg24(p)-SHR_CONST_TKFRZ)/catanf_30,0.01_r8) + t_scalar = max(catanf(bc_in%t_veg24_pa(ifp)-SHR_CONST_TKFRZ)/catanf_30,0.01_r8) endif !Moisture Limitations @@ -1186,7 +1177,7 @@ subroutine fragmentation_scaler( currentPatch, temperature_inst ) end subroutine fragmentation_scaler ! ============================================================================ - subroutine cwd_out( currentSite, currentPatch, temperature_inst ) + subroutine cwd_out( currentSite, currentPatch, bc_in ) ! ! !DESCRIPTION: ! Simple CWD fragmentation Model @@ -1198,8 +1189,9 @@ subroutine cwd_out( currentSite, currentPatch, temperature_inst ) ! ! !ARGUMENTS type(ed_site_type), intent(inout), target :: currentSite - type(ed_patch_type) , intent(inout), target :: currentPatch - type(temperature_type) , intent(in) :: temperature_inst + type(ed_patch_type), intent(inout), target :: currentPatch + type(bc_in_type), intent(in) :: bc_in + ! ! !LOCAL VARIABLES: integer :: c,ft @@ -1207,8 +1199,8 @@ subroutine cwd_out( currentSite, currentPatch, temperature_inst ) currentPatch%root_litter_out(:) = 0.0_r8 currentPatch%leaf_litter_out(:) = 0.0_r8 - - call fragmentation_scaler(currentPatch, temperature_inst) + + call fragmentation_scaler(currentPatch, bc_in) !Flux of coarse woody debris into decomposing litter pool. diff --git a/components/clm/src/ED/main/EDMainMod.F90 b/components/clm/src/ED/main/EDMainMod.F90 index ffd7948d31..d43c2be6e4 100755 --- a/components/clm/src/ED/main/EDMainMod.F90 +++ b/components/clm/src/ED/main/EDMainMod.F90 @@ -10,7 +10,6 @@ module EDMainMod use atm2lndType , only : atm2lnd_type use SoilStateType , only : soilstate_type use TemperatureType , only : temperature_type - use WaterStateType , only : waterstate_type use EDCohortDynamicsMod , only : allocate_live_biomass, terminate_cohorts, fuse_cohorts, sort_cohorts, count_cohorts use EDPatchDynamicsMod , only : disturbance_rates, fuse_patches, spawn_patches, terminate_patches use EDPhysiologyMod , only : canopy_derivs, non_canopy_derivs, phenology, recruitment, trim_canopy @@ -43,7 +42,7 @@ module EDMainMod !-------------------------------------------------------------------------------! subroutine ed_ecosystem_dynamics(currentSite, bc_in, & atm2lnd_inst, & - soilstate_inst, temperature_inst, waterstate_inst) + temperature_inst) ! ! !DESCRIPTION: ! Core of ed model, calling all subsequent vegetation dynamics routines @@ -52,9 +51,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, & type(ed_site_type) , intent(inout), target :: currentSite type(bc_in_type) , intent(in) :: bc_in 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(in) :: waterstate_inst ! ! !LOCAL VARIABLES: type(ed_patch_type), pointer :: currentPatch @@ -71,7 +68,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, & call ed_total_balance_check(currentSite, 0) - call phenology(currentSite, bc_in, waterstate_inst ) + call phenology(currentSite, bc_in ) call fire_model(currentSite, atm2lnd_inst, temperature_inst) @@ -79,7 +76,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, & call disturbance_rates(currentSite) ! Integrate state variables from annual rates to daily timestep - call ed_integrate_state_variables(currentSite, temperature_inst ) + call ed_integrate_state_variables(currentSite, bc_in ) !****************************************************************************** ! Reproduction, Recruitment and Cohort Dynamics : controls cohort organisation @@ -136,7 +133,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, & end subroutine ed_ecosystem_dynamics !-------------------------------------------------------------------------------! - subroutine ed_integrate_state_variables(currentSite, temperature_inst ) + subroutine ed_integrate_state_variables(currentSite, bc_in ) ! ! !DESCRIPTION: ! FIX(SPM,032414) refactor so everything goes through interface @@ -144,8 +141,9 @@ subroutine ed_integrate_state_variables(currentSite, temperature_inst ) ! !USES: ! ! !ARGUMENTS: - type(ed_site_type) , intent(inout) :: currentSite - type(temperature_type) , intent(in) :: temperature_inst + type(ed_site_type) , intent(inout) :: currentSite + type(bc_in_type) , intent(in) :: bc_in + ! ! !LOCAL VARIABLES: type(ed_patch_type) , pointer :: currentPatch @@ -223,7 +221,7 @@ subroutine ed_integrate_state_variables(currentSite, temperature_inst ) write(6,*)'DEBUG18: calling non_canopy_derivs with pno= ',currentPatch%clm_pno endif - call non_canopy_derivs( currentSite, currentPatch, temperature_inst ) + call non_canopy_derivs( currentSite, currentPatch, bc_in) !update state variables simultaneously according to derivatives for this time period. diff --git a/components/clm/src/ED/main/FatesInterfaceMod.F90 b/components/clm/src/ED/main/FatesInterfaceMod.F90 index e61a6dda5d..afc39fcdd6 100644 --- a/components/clm/src/ED/main/FatesInterfaceMod.F90 +++ b/components/clm/src/ED/main/FatesInterfaceMod.F90 @@ -64,17 +64,21 @@ module FatesInterfaceMod ! Vegetation Dynamics ! --------------------------------------------------------------------------------- - ! This 24 hour vegetation temperature is used for various purposes during vegetation - ! dynamics. However, we are currently using the bare-ground patch's value [K] + ! The site level 24 hour vegetation temperature is used for various purposes during vegetation + ! dynamics. However, we are currently using the bare ground patch's value [K] ! TO-DO: Get some consensus on the correct vegetation temperature used for phenology. ! It is possible that the bare-ground value is where the average is being stored. ! (RGK-01-2017) - - real(r8) :: t_veg24_si ! See above [K] + real(r8) :: t_veg24_si ! Patch 24 hour vegetation temperature [K] real(r8),allocatable :: t_veg24_pa(:) - + + ! NOTE: h2osoi_vol_si is used to update surface water memory + ! CLM/ALM may be using "waterstate%h2osoi_vol_col" on the first index (coli,1) + ! to inform this. I think this should be re-evaluated (RGK 01/2017) + ! Site volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3] + real(r8) :: h2osoi_vol_si ! Radiation variables for calculating sun/shade fractions ! --------------------------------------------------------------------------------- @@ -484,9 +488,9 @@ subroutine zero_bcs(this,s) this%bc_in(s)%current_date = 0 this%bc_in(s)%reference_date = 0 this%bc_in(s)%model_day = 0.0_r8 - - this%bc_in(s)%t_veg24_pa(:) = 0.0_r8 this%bc_in(s)%t_veg24_si = 0.0_r8 + this%bc_in(s)%t_veg24_pa(:) = 0.0_r8 + this%bc_in(s)%h2osoi_vol_si = 0.0_r8 this%bc_in(s)%solad_parb(:,:) = 0.0_r8 this%bc_in(s)%solai_parb(:,:) = 0.0_r8 diff --git a/components/clm/src/utils/clmfates_interfaceMod.F90 b/components/clm/src/utils/clmfates_interfaceMod.F90 index 747427312b..7400a69d99 100644 --- a/components/clm/src/utils/clmfates_interfaceMod.F90 +++ b/components/clm/src/utils/clmfates_interfaceMod.F90 @@ -513,8 +513,11 @@ subroutine dynamics_driv(this, nc, bounds_clump, & this%fates(nc)%bc_in(s)%current_tod = current_date this%fates(nc)%bc_in(s)%reference_date = reference_date this%fates(nc)%bc_in(s)%model_day = model_day - this%fates(nc)%bc_in(s)%t_veg24_si = & - temperature_inst%t_veg24_patch(col%patchi(c)-1) + this%fates(nc)%bc_in(s)%h2osoi_vol_si = & + waterstate_inst%h2osoi_vol_col(c,1) + + this%fates(nc)%bc_in(s)%t_veg24_si = & + temperature_inst%t_veg24_patch(col%patchi(c)) do ifp = 1, this%fates(nc)%sites(s)%youngest_patch%patchno p = ifp+col%patchi(c) this%fates(nc)%bc_in(s)%t_veg24_pa(ifp) = & @@ -533,7 +536,7 @@ subroutine dynamics_driv(this, nc, bounds_clump, & call ed_ecosystem_dynamics(this%fates(nc)%sites(s), & this%fates(nc)%bc_in(s), & atm2lnd_inst, & - soilstate_inst, temperature_inst, waterstate_inst) + temperature_inst ) call ed_update_site(this%fates(nc)%sites(s)) From 76662a1456365f85353e7d750c686801dde31ad4 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 11 Jan 2017 14:26:36 -0800 Subject: [PATCH 04/11] Clean up of globals in EDPhysiology --- .../clm/src/ED/biogeochem/EDPhysiologyMod.F90 | 183 +++++++++--------- .../clm/src/ED/main/FatesConstantsMod.F90 | 11 +- 2 files changed, 104 insertions(+), 90 deletions(-) diff --git a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 index 9d84507865..918065809b 100755 --- a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 +++ b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 @@ -6,10 +6,8 @@ module EDPhysiologyMod ! Miscellaneous physiology routines from ED. ! ============================================================================ - use shr_kind_mod , only : r8 => shr_kind_r8 - use clm_varctl , only : iulog - - use TemperatureType , only : temperature_type + use FatesGlobals, only : fates_log + use FatesConstantsMod, only : r8 => fates_r8 use pftconMod , only : pftcon use EDEcophysContype , only : EDecophyscon use FatesInterfaceMod, only : bc_in_type @@ -114,7 +112,8 @@ subroutine non_canopy_derivs( currentSite, currentPatch, bc_in ) do p = 1,numpft_ed currentSite%dseed_dt(p) = currentSite%dseed_dt(p) + & - (currentPatch%seeds_in(p) - currentPatch%seed_decay(p) - currentPatch%seed_germination(p)) * currentPatch%area/AREA + (currentPatch%seeds_in(p) - currentPatch%seed_decay(p) - & + currentPatch%seed_germination(p)) * currentPatch%area/AREA enddo do c = 1,ncwd @@ -123,19 +122,12 @@ subroutine non_canopy_derivs( currentSite, currentPatch, bc_in ) enddo do p = 1,numpft_ed - currentPatch%dleaf_litter_dt(p) = currentPatch%leaf_litter_in(p) - currentPatch%leaf_litter_out(p) - currentPatch%droot_litter_dt(p) = currentPatch%root_litter_in(p) - currentPatch%root_litter_out(p) + currentPatch%dleaf_litter_dt(p) = currentPatch%leaf_litter_in(p) - & + currentPatch%leaf_litter_out(p) + currentPatch%droot_litter_dt(p) = currentPatch%root_litter_in(p) - & + currentPatch%root_litter_out(p) enddo - ! currentPatch%leaf_litter_in(:) = 0.0_r8 - ! currentPatch%root_litter_in(:) = 0.0_r8 - ! currentPatch%leaf_litter_out(:) = 0.0_r8 - ! currentPatch%root_litter_out(:) = 0.0_r8 - ! currentPatch%CWD_AG_in(:) = 0.0_r8 - ! currentPatch%cwd_bg_in(:) = 0.0_r8 - ! currentPatch%CWD_AG_out(:) = 0.0_r8 - ! currentPatch%cwd_bg_out(:) = 0.0_r8 - end subroutine non_canopy_derivs ! ============================================================================ @@ -176,7 +168,7 @@ subroutine trim_canopy( currentSite ) currentCohort%treelai = tree_lai(currentCohort) currentCohort%nv = ceiling((currentCohort%treelai+currentCohort%treesai)/dinc_ed) if (currentCohort%nv > cp_nlevcan)then - write(iulog,*) 'nv > cp_nlevcan',currentCohort%nv,currentCohort%treelai,currentCohort%treesai, & + write(fates_log(),*) 'nv > cp_nlevcan',currentCohort%nv,currentCohort%treelai,currentCohort%treesai, & currentCohort%c_area,currentCohort%n,currentCohort%bl endif @@ -201,7 +193,7 @@ subroutine trim_canopy( currentSite ) if (currentCohort%canopy_trim > trim_limit)then if ( DEBUG ) then - write(iulog,*) 'trimming leaves',currentCohort%canopy_trim,currentCohort%leaf_cost + write(fates_log(),*) 'trimming leaves',currentCohort%canopy_trim,currentCohort%leaf_cost endif ! keep trimming until none of the canopy is in negative carbon balance. @@ -219,7 +211,7 @@ subroutine trim_canopy( currentSite ) if (currentCohort%NV.gt.2)then ! leaf_cost may be uninitialized, removing its diagnostic from the log ! to allow checking with fpe_traps (RGK) - write(iulog,*) 'nv>4',currentCohort%year_net_uptake(1:6),currentCohort%canopy_trim + write(fates_log(),*) 'nv>4',currentCohort%year_net_uptake(1:6),currentCohort%canopy_trim endif currentCohort%year_net_uptake(:) = 999.0_r8 @@ -228,7 +220,7 @@ subroutine trim_canopy( currentSite ) endif if ( DEBUG ) then - write(iulog,*) 'trimming',currentCohort%canopy_trim + write(fates_log(),*) 'trimming',currentCohort%canopy_trim endif ! currentCohort%canopy_trim = 1.0_r8 !FIX(RF,032414) this turns off ctrim for now. @@ -357,7 +349,7 @@ subroutine phenology( currentSite, bc_in ) currentSite%status = 2 !alter status of site to 'leaves on' ! NOTE(bja, 2015-01) should leafondate = model_day to be consistent with leaf off? currentSite%leafondate = t !record leaf on date - if ( DEBUG ) write(iulog,*) 'leaves on' + if ( DEBUG ) write(fates_log(),*) 'leaves on' endif !ncd endif !status endif !GDD @@ -377,7 +369,7 @@ subroutine phenology( currentSite, bc_in ) if (currentSite%status == 2)then currentSite%status = 1 !alter status of site to 'leaves on' currentSite%leafoffdate = bc_in%model_day !record leaf off date - if ( DEBUG ) write(iulog,*) 'leaves off' + if ( DEBUG ) write(fates_log(),*) 'leaves off' endif endif endif @@ -387,7 +379,7 @@ subroutine phenology( currentSite, bc_in ) if(currentSite%status == 2)then currentSite%status = 1 !alter status of site to 'leaves on' currentSite%leafoffdate = bc_in%model_day !record leaf off date - if ( DEBUG ) write(iulog,*) 'leaves off' + if ( DEBUG ) write(fates_log(),*) 'leaves off' endif endif @@ -503,7 +495,8 @@ subroutine phenology_leafonoff(currentSite) type(ed_patch_type) , pointer :: currentPatch type(ed_cohort_type), pointer :: currentCohort - real(r8) :: store_output ! the amount of the store to put into leaves - is a barrier against negative storage and C starvation. + real(r8) :: store_output ! the amount of the store to put into leaves - + ! is a barrier against negative storage and C starvation. !------------------------------------------------------------------------ @@ -532,11 +525,11 @@ subroutine phenology_leafonoff(currentSite) ! Add deployed carbon to alive biomass pool currentCohort%balive = currentCohort%balive + currentCohort%bl - if ( DEBUG ) write(iulog,*) 'EDPhysMod 1 ',currentCohort%bstore + if ( DEBUG ) write(fates_log(),*) 'EDPhysMod 1 ',currentCohort%bstore currentCohort%bstore = currentCohort%bstore - currentCohort%bl ! Drain store - if ( DEBUG ) write(iulog,*) 'EDPhysMod 2 ',currentCohort%bstore + if ( DEBUG ) write(fates_log(),*) 'EDPhysMod 2 ',currentCohort%bstore currentCohort%laimemory = 0.0_r8 @@ -571,11 +564,11 @@ subroutine phenology_leafonoff(currentSite) endif currentCohort%balive = currentCohort%balive + currentCohort%bl - if ( DEBUG ) write(iulog,*) 'EDPhysMod 3 ',currentCohort%bstore + if ( DEBUG ) write(fates_log(),*) 'EDPhysMod 3 ',currentCohort%bstore currentCohort%bstore = currentCohort%bstore - currentCohort%bl ! empty store - if ( DEBUG ) write(iulog,*) 'EDPhysMod 4 ',currentCohort%bstore + if ( DEBUG ) write(fates_log(),*) 'EDPhysMod 4 ',currentCohort%bstore currentCohort%laimemory = 0.0_r8 @@ -633,7 +626,8 @@ subroutine seeds_in( currentSite, cp_pnt ) currentCohort => currentPatch%tallest do while (associated(currentCohort)) p = currentCohort%pft - currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + currentCohort%seed_prod * currentCohort%n/currentPatch%area + currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + & + currentCohort%seed_prod * currentCohort%n/currentPatch%area currentCohort => currentCohort%shorter enddo !cohort loop @@ -642,8 +636,10 @@ subroutine seeds_in( currentSite, cp_pnt ) do while(associated(currentPatch)) if (EXTERNAL_RECRUITMENT == 1) then !external seed rain - needed to prevent extinction do p = 1,numpft_ed - currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + EDecophyscon%seed_rain(p) !KgC/m2/year - currentSite%seed_rain_flux(p) = currentSite%seed_rain_flux(p) + EDecophyscon%seed_rain(p) * currentPatch%area/AREA !KgC/m2/year + currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + & + EDecophyscon%seed_rain(p) !KgC/m2/year + currentSite%seed_rain_flux(p) = currentSite%seed_rain_flux(p) + & + EDecophyscon%seed_rain(p) * currentPatch%area/AREA !KgC/m2/year enddo endif currentPatch => currentPatch%younger @@ -699,7 +695,8 @@ subroutine seed_germination( currentSite, currentPatch ) max_germination = 1.0_r8 !this is arbitrary do p = 1,numpft_ed - currentPatch%seed_germination(p) = min(currentSite%seed_bank(p) * germination_timescale,max_germination) + currentPatch%seed_germination(p) = min(currentSite%seed_bank(p) * & + germination_timescale,max_germination) enddo end subroutine seed_germination @@ -761,7 +758,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort) endif ! NPP - if ( DEBUG ) write(iulog,*) 'EDphys 716 ',currentCohort%npp_acc + if ( DEBUG ) write(fates_log(),*) 'EDphys 716 ',currentCohort%npp_acc currentCohort%npp_acc_hold = currentCohort%npp_acc * udata%n_sub ! convert from kgC/indiv/day into kgC/indiv/year currentCohort%gpp_acc_hold = currentCohort%gpp_acc * udata%n_sub ! convert from kgC/indiv/day into kgC/indiv/year @@ -795,7 +792,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort) if (pftcon%stress_decid(currentCohort%pft) /= 1.and.pftcon%season_decid(currentCohort%pft) /= 1.and. & pftcon%evergreen(currentCohort%pft) /= 1)then - write(iulog,*) 'problem with phenology definitions',currentCohort%pft,pftcon%stress_decid(currentCohort%pft), & + write(fates_log(),*) 'problem with phenology definitions',currentCohort%pft,pftcon%stress_decid(currentCohort%pft), & pftcon%season_decid(currentCohort%pft),pftcon%evergreen(currentCohort%pft) endif @@ -807,7 +804,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort) ! Calculate carbon balance ! this is the fraction of maintenance demand we -have- to do... - if ( DEBUG ) write(iulog,*) 'EDphys 760 ',currentCohort%npp_acc_hold, currentCohort%md, & + if ( DEBUG ) write(fates_log(),*) 'EDphys 760 ',currentCohort%npp_acc_hold, currentCohort%md, & EDecophyscon%leaf_stor_priority(currentCohort%pft) currentCohort%carbon_balance = currentCohort%npp_acc_hold - & @@ -823,7 +820,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort) if (Bleaf(currentCohort) > 0._r8)then - if ( DEBUG ) write(iulog,*) 'EDphys A ',currentCohort%carbon_balance + if ( DEBUG ) write(fates_log(),*) 'EDphys A ',currentCohort%carbon_balance if (currentCohort%carbon_balance > 0._r8)then !spend C on growing and storing @@ -835,7 +832,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort) !what is the flux into the store? currentCohort%storage_flux = currentCohort%carbon_balance * f_store - if ( DEBUG ) write(iulog,*) 'EDphys B ',f_store + if ( DEBUG ) write(fates_log(),*) 'EDphys B ',f_store !what is the tax on the carbon available for growth? currentCohort%carbon_balance = currentCohort%carbon_balance * (1.0_r8 - f_store) @@ -848,7 +845,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort) currentCohort%storage_flux = 0._r8 currentCohort%carbon_balance = 0._r8 - write(iulog,*) 'ED: no leaf area in gd', currentCohort%indexnumber,currentCohort%n,currentCohort%bdead, & + write(fates_log(),*) 'ED: no leaf area in gd', currentCohort%indexnumber,currentCohort%n,currentCohort%bdead, & currentCohort%dbh,currentCohort%balive endif @@ -902,7 +899,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort) !FIX(RF,032414) - to fix high bl's. needed to prevent numerical errors without the ODEINT. if (currentCohort%balive > target_balive*1.1_r8)then va = 0.0_r8; vs = 1._r8 - write(iulog,*) 'using high bl cap',target_balive,currentCohort%balive + write(fates_log(),*) 'using high bl cap',target_balive,currentCohort%balive endif else @@ -916,28 +913,28 @@ subroutine Growth_Derivatives( currentSite, currentCohort) currentCohort%dbdeaddt = gr_fract * vs * currentCohort%carbon_balance currentCohort%dbstoredt = currentCohort%storage_flux - if ( DEBUG ) write(iulog,*) 'EDPhys dbstoredt I ',currentCohort%dbstoredt + if ( DEBUG ) write(fates_log(),*) 'EDPhys dbstoredt I ',currentCohort%dbstoredt currentCohort%seed_prod = (1.0_r8 - gr_fract) * currentCohort%carbon_balance if (abs(currentCohort%npp_acc_hold-(currentCohort%dbalivedt+currentCohort%dbdeaddt+currentCohort%dbstoredt+ & currentCohort%seed_prod+currentCohort%md)) > 0.0000000001_r8)then - write(iulog,*) 'error in carbon check growth derivs',currentCohort%npp_acc_hold- & + write(fates_log(),*) 'error in carbon check growth derivs',currentCohort%npp_acc_hold- & (currentCohort%dbalivedt+currentCohort%dbdeaddt+currentCohort%dbstoredt+currentCohort%seed_prod+currentCohort%md) - write(iulog,*) 'cohort fluxes',currentCohort%pft,currentCohort%canopy_layer,currentCohort%n, & + write(fates_log(),*) 'cohort fluxes',currentCohort%pft,currentCohort%canopy_layer,currentCohort%n, & currentCohort%npp_acc_hold,currentCohort%dbalivedt,balive_loss, & currentCohort%dbdeaddt,currentCohort%dbstoredt,currentCohort%seed_prod,currentCohort%md * & EDecophyscon%leaf_stor_priority(currentCohort%pft) - write(iulog,*) 'proxies' ,target_balive,currentCohort%balive,currentCohort%dbh,va,vs,gr_fract + write(fates_log(),*) 'proxies' ,target_balive,currentCohort%balive,currentCohort%dbh,va,vs,gr_fract endif ! prevent negative leaf pool (but not negative store pool). This is also a numerical error prevention, ! but it shouldn't happen actually... if (-1.0_r8*currentCohort%dbalivedt * udata%deltat > currentCohort%balive*0.99)then - write(iulog,*) 'using non-neg leaf mass cap',currentCohort%balive , currentCohort%dbalivedt,currentCohort%dbstoredt, & + write(fates_log(),*) 'using non-neg leaf mass cap',currentCohort%balive , currentCohort%dbalivedt,currentCohort%dbstoredt, & currentCohort%carbon_balance currentCohort%dbstoredt = currentCohort%dbstoredt + currentCohort%dbalivedt - if ( DEBUG ) write(iulog,*) 'EDPhys dbstoredt II ',currentCohort%dbstoredt + if ( DEBUG ) write(fates_log(),*) 'EDPhys dbstoredt II ',currentCohort%dbstoredt currentCohort%dbalivedt = 0._r8 endif @@ -993,10 +990,10 @@ subroutine recruitment( t, currentSite, currentPatch ) / (temp_cohort%bdead+temp_cohort%balive+temp_cohort%bstore) if (t == 1)then - write(iulog,*) 'filling in cohorts where there are none left; this will break carbon balance', & + write(fates_log(),*) 'filling in cohorts where there are none left; this will break carbon balance', & currentPatch%patchno,currentPatch%area temp_cohort%n = 0.1_r8*currentPatch%area - write(iulog,*) 'cohort n',ft,temp_cohort%n + write(fates_log(),*) 'cohort n',ft,temp_cohort%n endif temp_cohort%laimemory = 0.0_r8 @@ -1015,7 +1012,7 @@ subroutine recruitment( t, currentSite, currentPatch ) endif if (temp_cohort%n > 0.0_r8 )then - if ( DEBUG ) write(iulog,*) 'EDPhysiologyMod.F90 call create_cohort ' + if ( DEBUG ) write(fates_log(),*) 'EDPhysiologyMod.F90 call create_cohort ' call create_cohort(currentPatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, & temp_cohort%balive, temp_cohort%bdead, temp_cohort%bstore, & temp_cohort%laimemory, cohortstatus, temp_cohort%canopy_trim, currentPatch%NCL_p) @@ -1096,7 +1093,7 @@ subroutine CWD_Input( currentPatch) SF_val_CWD_frac(c) * dead_n * (1.0_r8-ED_val_ag_biomass) if (currentPatch%cwd_AG_in(c) < 0.0_r8)then - write(iulog,*) 'negative CWD in flux',currentPatch%cwd_AG_in(c), & + write(fates_log(),*) 'negative CWD in flux',currentPatch%cwd_AG_in(c), & (currentCohort%bdead+currentCohort%bsw), dead_n endif enddo @@ -1118,12 +1115,14 @@ subroutine fragmentation_scaler( currentPatch, bc_in) ! ! !DESCRIPTION: ! Simple CWD fragmentation Model - ! FIX(SPM, 091914) this should be a function as it returns a value in currentPatch%fragmentation_scaler + ! FIX(SPM, 091914) this should be a function as it returns a value in + ! currentPatch%fragmentation_scaler ! ! !USES: - use shr_const_mod , only : SHR_CONST_PI, SHR_CONST_TKFRZ - use EDSharedParamsMod , only : EDParamsShareInst + use EDSharedParamsMod , only : EDParamsShareInst + use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm + use FatesConstantsMod, only : pi => pi_const ! ! !ARGUMENTS type(ed_patch_type), intent(inout) :: currentPatch @@ -1140,10 +1139,11 @@ subroutine fragmentation_scaler( currentPatch, bc_in) real(r8) :: catanf_30 ! hyperbolic temperature function from CENTURY real(r8) :: t1 ! temperature argument real(r8) :: Q10 ! temperature dependence - real(r8) :: froz_q10 ! separate q10 for frozen soil respiration rates. default to same as above zero rates + real(r8) :: froz_q10 ! separate q10 for frozen soil respiration rates. + ! default to same as above zero rates !---------------------------------------------------------------------- - catanf(t1) = 11.75_r8 +(29.7_r8 / SHR_CONST_PI) * atan( SHR_CONST_PI * 0.031_r8 * ( t1 - 15.4_r8 )) + catanf(t1) = 11.75_r8 +(29.7_r8 / pi) * atan( pi * 0.031_r8 * ( t1 - 15.4_r8 )) catanf_30 = catanf(30._r8) ifp = currentPatch%patchno @@ -1155,20 +1155,22 @@ subroutine fragmentation_scaler( currentPatch, bc_in) 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 (bc_in%t_veg24_pa(ifp) >= SHR_CONST_TKFRZ) then - t_scalar = Q10**((bc_in%t_veg24_pa(ifp)-(SHR_CONST_TKFRZ+25._r8))/10._r8) - ! Q10**((t_soisno(c,j)-(SHR_CONST_TKFRZ+25._r8))/10._r8) + if (bc_in%t_veg24_pa(ifp) >= tfrz) then + t_scalar = Q10**((bc_in%t_veg24_pa(ifp)-(tfrz+25._r8))/10._r8) + ! Q10**((t_soisno(c,j)-(tfrz+25._r8))/10._r8) else - t_scalar = (Q10**(-25._r8/10._r8))*(froz_q10**((bc_in%t_veg24_pa(ifp)-SHR_CONST_TKFRZ)/10._r8)) - !Q10**(-25._r8/10._r8))*(froz_q10**((t_soisno(c,j)-SHR_CONST_TKFRZ)/10._r8) + t_scalar = (Q10**(-25._r8/10._r8))*(froz_q10**((bc_in%t_veg24_pa(ifp)-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(bc_in%t_veg24_pa(ifp)-SHR_CONST_TKFRZ)/catanf_30,0.01_r8) + ! original century uses an arctangent function to calculate the + ! temperature dependence of decomposition + t_scalar = max(catanf(bc_in%t_veg24_pa(ifp)-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. + !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_ed))/numpft_ed @@ -1227,7 +1229,7 @@ subroutine cwd_out( currentSite, currentPatch, bc_in ) currentPatch%root_litter_out(ft) = max(0.0_r8,currentPatch%root_litter(ft)* SF_val_max_decomp(dg_sf) * & currentPatch%fragmentation_scaler ) if ( currentPatch%leaf_litter_out(ft)<0.0_r8.or.currentPatch%root_litter_out(ft)<0.0_r8)then - write(iulog,*) 'root or leaf out is negative?',SF_val_max_decomp(dg_sf),currentPatch%fragmentation_scaler + write(fates_log(),*) 'root or leaf out is negative?',SF_val_max_decomp(dg_sf),currentPatch%fragmentation_scaler endif enddo @@ -1268,14 +1270,15 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) !use EDCLMLinkMod, only: cwd_fcel_ed, cwd_flig use pftconMod, only : pftcon - use shr_const_mod, only: SHR_CONST_CDAY + use FatesConstantsMod, only : sec_per_day use clm_varcon, only : zisoi, dzsoi_decomp, zsoi use EDParamsMod, only : ED_val_ag_biomass use FatesInterfaceMod, only : bc_in_type, bc_out_type use clm_varctl, only : use_vertsoilc use abortutils , only : endrun - ! INTERF-TODO: remove the control parameters: exponential_rooting_profile, pftspecific_rootingprofile, rootprof_exp, surfprof_exp, zisoi, dzsoi_decomp, zsoi + ! INTERF-TODO: remove the control parameters: exponential_rooting_profile, + ! pftspecific_rootingprofile, rootprof_exp, surfprof_exp, zisoi, dzsoi_decomp, zsoi ! implicit none ! @@ -1334,7 +1337,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) delta = 0.001_r8 !no of seconds in a year. - time_convert = 365.0_r8*SHR_CONST_CDAY + time_convert = 365.0_r8*sec_per_day ! number of grams in a kilogram mass_convert = 1000._r8 @@ -1343,8 +1346,10 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! first calculate vertical profiles ! define two types of profiles: - ! (1) a surface profile, for leaves and stem inputs, which is the same for each pft but differs from one site to the next to avoid inputting any C into permafrost or bedrock - ! (2) a fine root profile, which is indexed by both site and pft, differs for each pft and also from one site to the next to avoid inputting any C into permafrost or bedrock + ! (1) a surface profile, for leaves and stem inputs, which is the same for each + ! pft but differs from one site to the next to avoid inputting any C into permafrost or bedrock + ! (2) a fine root profile, which is indexed by both site and pft, differs for + ! each pft and also from one site to the next to avoid inputting any C into permafrost or bedrock ! (3) a coarse root profile, which is the root-biomass=weighted average of the fine root profiles !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1468,13 +1473,13 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) stem_prof_sum = stem_prof_sum + stem_prof(s,j) * dzsoi_decomp(j) end do if ( ( abs(stem_prof_sum - 1._r8) > delta ) .or. ( abs(leaf_prof_sum - 1._r8) > delta ) ) then - write(iulog, *) 'profile sums: ', leaf_prof_sum, stem_prof_sum - write(iulog, *) 'surface_prof: ', surface_prof - write(iulog, *) 'surface_prof_tot: ', surface_prof_tot - write(iulog, *) 'leaf_prof: ', leaf_prof(s,:) - write(iulog, *) 'stem_prof: ', stem_prof(s,:) - write(iulog, *) 'max_rooting_depth_index_col: ', bc_in(s)%max_rooting_depth_index_col - write(iulog, *) 'dzsoi_decomp: ', dzsoi_decomp + write(fates_log(), *) 'profile sums: ', leaf_prof_sum, stem_prof_sum + write(fates_log(), *) 'surface_prof: ', surface_prof + write(fates_log(), *) 'surface_prof_tot: ', surface_prof_tot + write(fates_log(), *) 'leaf_prof: ', leaf_prof(s,:) + write(fates_log(), *) 'stem_prof: ', stem_prof(s,:) + write(fates_log(), *) 'max_rooting_depth_index_col: ', bc_in(s)%max_rooting_depth_index_col + write(fates_log(), *) 'dzsoi_decomp: ', dzsoi_decomp call endrun() endif ! now check each fine root profile @@ -1484,7 +1489,7 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) froot_prof_sum = froot_prof_sum + froot_prof(s,ft,j) * dzsoi_decomp(j) end do if ( ( abs(froot_prof_sum - 1._r8) > delta ) ) then - write(iulog, *) 'profile sums: ', froot_prof_sum + write(fates_log(), *) 'profile sums: ', froot_prof_sum call endrun() endif end do @@ -1552,12 +1557,12 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) ! now disaggregate, vertically and by decomposition substrate type, the actual fluxes from CWD and litter pools ! ! do c = 1, ncwd - ! write(iulog,*)'cdk CWD_AG_out', c, currentpatch%CWD_AG_out(c), cwd_fcel, currentpatch%area/AREA - ! write(iulog,*)'cdk CWD_BG_out', c, currentpatch%CWD_BG_out(c), cwd_fcel, currentpatch%area/AREA + ! write(fates_log(),*)'cdk CWD_AG_out', c, currentpatch%CWD_AG_out(c), cwd_fcel, currentpatch%area/AREA + ! write(fates_log(),*)'cdk CWD_BG_out', c, currentpatch%CWD_BG_out(c), cwd_fcel, currentpatch%area/AREA ! end do ! do ft = 1,numpft_ed - ! write(iulog,*)'cdk leaf_litter_out', ft, currentpatch%leaf_litter_out(ft), cwd_fcel, currentpatch%area/AREA - ! write(iulog,*)'cdk root_litter_out', ft, currentpatch%root_litter_out(ft), cwd_fcel, currentpatch%area/AREA + ! write(fates_log(),*)'cdk leaf_litter_out', ft, currentpatch%leaf_litter_out(ft), cwd_fcel, currentpatch%area/AREA + ! write(fates_log(),*)'cdk root_litter_out', ft, currentpatch%root_litter_out(ft), cwd_fcel, currentpatch%area/AREA ! end do ! ! ! CWD pools fragmenting into decomposing litter pools. @@ -1618,15 +1623,15 @@ subroutine flux_into_litter_pools(nsites, sites, bc_in, bc_out) end do end do - ! write(iulog,*)'cdk FATES_c_to_litr_lab_c: ', FATES_c_to_litr_lab_c - ! write_col(iulog,*)'cdk FATES_c_to_litr_cel_c: ', FATES_c_to_litr_cel_c - ! write_col(iulog,*)'cdk FATES_c_to_litr_lig_c: ', FATES_c_to_litr_lig_c - ! write_col(iulog,*)'cdk cp_numlevdecomp_full, bounds%begc, bounds%endc: ', cp_numlevdecomp_full, bounds%begc, bounds%endc - ! write(iulog,*)'cdk leaf_prof: ', leaf_prof - ! write(iulog,*)'cdk stem_prof: ', stem_prof - ! write(iulog,*)'cdk froot_prof: ', froot_prof - ! write(iulog,*)'cdk croot_prof_perpatch: ', croot_prof_perpatch - ! write(iulog,*)'cdk croot_prof: ', croot_prof + ! write(fates_log(),*)'cdk FATES_c_to_litr_lab_c: ', FATES_c_to_litr_lab_c + ! write_col(fates_log(),*)'cdk FATES_c_to_litr_cel_c: ', FATES_c_to_litr_cel_c + ! write_col(fates_log(),*)'cdk FATES_c_to_litr_lig_c: ', FATES_c_to_litr_lig_c + ! write_col(fates_log(),*)'cdk cp_numlevdecomp_full, bounds%begc, bounds%endc: ', cp_numlevdecomp_full, bounds%begc, bounds%endc + ! write(fates_log(),*)'cdk leaf_prof: ', leaf_prof + ! write(fates_log(),*)'cdk stem_prof: ', stem_prof + ! write(fates_log(),*)'cdk froot_prof: ', froot_prof + ! write(fates_log(),*)'cdk croot_prof_perpatch: ', croot_prof_perpatch + ! write(fates_log(),*)'cdk croot_prof: ', croot_prof end subroutine flux_into_litter_pools diff --git a/components/clm/src/ED/main/FatesConstantsMod.F90 b/components/clm/src/ED/main/FatesConstantsMod.F90 index 3df36d6b56..b7bf5edb9b 100644 --- a/components/clm/src/ED/main/FatesConstantsMod.F90 +++ b/components/clm/src/ED/main/FatesConstantsMod.F90 @@ -39,8 +39,11 @@ module FatesConstantsMod real(fates_r8), parameter :: umol_per_mol = 1.0E6_fates_r8 - ! Conversion: secons per minute + ! Conversion: seconds per minute real(fates_r8), parameter :: sec_per_min = 60.0_fates_r8 + + ! Conversion: seconds per day + real(fates_r8), parameter :: sec_per_day = 86400.0_fates_r8 ! Physical constants @@ -55,4 +58,10 @@ module FatesConstantsMod real(fates_r8), parameter :: t_water_freeze_k_triple = 273.16_fates_r8 + ! Geometric Constants + + ! PI + real(fates_r8), parameter :: pi_const = 3.14159265359_fates_r8 + + end module FatesConstantsMod From 78ae15a9eacceb738a47c7fc41e9f792942b12da Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 11 Jan 2017 15:13:44 -0800 Subject: [PATCH 05/11] Removal of CLM boundary conditions to fire model and some CLM globals. ldomain%area(g) is still called in fire model, as well as masterproc and use_spitfire. --- components/clm/src/ED/fire/SFMainMod.F90 | 210 +++++++++--------- components/clm/src/ED/main/EDMainMod.F90 | 8 +- .../clm/src/ED/main/FatesInterfaceMod.F90 | 18 ++ .../clm/src/utils/clmfates_interfaceMod.F90 | 14 +- 4 files changed, 138 insertions(+), 112 deletions(-) diff --git a/components/clm/src/ED/fire/SFMainMod.F90 b/components/clm/src/ED/fire/SFMainMod.F90 index be53100a71..1f5f14f235 100755 --- a/components/clm/src/ED/fire/SFMainMod.F90 +++ b/components/clm/src/ED/fire/SFMainMod.F90 @@ -5,15 +5,24 @@ module SFMainMod ! Code originally developed by Allan Spessa & Rosie Fisher as part of the NERC-QUEST project. ! ============================================================================ - use shr_kind_mod , only : r8 => shr_kind_r8; + use FatesConstantsMod , only : r8 => fates_r8 + use spmdMod , only : masterproc - use clm_varctl , only : iulog - use atm2lndType , only : atm2lnd_type - use TemperatureType , only : temperature_type + use FatesGlobals , only : fates_log + + use FatesInterfaceMod , only : bc_in_type use pftconMod , only : pftcon use EDEcophysconType , only : EDecophyscon - use EDtypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type, AREA, DG_SF, FIRE_THRESHOLD - use EDtypesMod , only : LB_SF, LG_SF, NCWD, TR_SF + use EDtypesMod , only : ed_site_type + use EDtypesMod , only : ed_patch_type + use EDtypesMod , only : ed_cohort_type + use EDtypesMod , only : AREA + use EDtypesMod , only : DG_SF + use EDtypesMod , only : FIRE_THRESHOLD + use EDtypesMod , only : LB_SF + use EDtypesMod , only : LG_SF + use EDtypesMod , only : NCWD + use EDtypesMod , only : TR_SF implicit none private @@ -42,13 +51,13 @@ module SFMainMod ! ============================================================================ ! Area of site burned by fire ! ============================================================================ - subroutine fire_model( currentSite, atm2lnd_inst, temperature_inst) + subroutine fire_model( currentSite, bc_in) use clm_varctl, only : use_ed_spit_fire type(ed_site_type) , intent(inout), target :: currentSite - type(atm2lnd_type) , intent(in) :: atm2lnd_inst - type(temperature_type) , intent(in) :: temperature_inst + type(bc_in_type) , intent(in) :: bc_in + type (ed_patch_type), pointer :: currentPatch @@ -62,12 +71,12 @@ subroutine fire_model( currentSite, atm2lnd_inst, temperature_inst) enddo if(write_SF==1)then - write(iulog,*) 'use_ed_spit_fire',use_ed_spit_fire + write(fates_log(),*) 'use_ed_spit_fire',use_ed_spit_fire endif if(use_ed_spit_fire)then - call fire_danger_index(currentSite, temperature_inst, atm2lnd_inst) - call wind_effect(currentSite, atm2lnd_inst) + call fire_danger_index(currentSite, bc_in) + call wind_effect(currentSite, bc_in) call charecteristics_of_fuel(currentSite) call rate_of_spread(currentSite) call ground_fuel_consumption(currentSite) @@ -81,20 +90,19 @@ subroutine fire_model( currentSite, atm2lnd_inst, temperature_inst) end subroutine fire_model - !***************************************************************** - subroutine fire_danger_index ( currentSite, temperature_inst, atm2lnd_inst) + !***************************************************************** + subroutine fire_danger_index ( currentSite, bc_in) - !***************************************************************** + !***************************************************************** ! currentSite%acc_NI is the accumulated Nesterov fire danger index - use clm_varcon , only : tfrz - use SFParamsMod, only : SF_val_fdi_a, SF_val_fdi_b - + use FatesConstantsMod , only : tfrz => t_water_freeze_k_1atm + use FatesConstantsMod , only : sec_per_day + type(ed_site_type) , intent(inout), target :: currentSite - type(temperature_type) , intent(in) :: temperature_inst - type(atm2lnd_type) , intent(in) :: atm2lnd_inst - + type(bc_in_type) , intent(in) :: bc_in + real(r8) :: temp_in_C ! daily averaged temperature in celcius real(r8) :: rainfall ! daily precip real(r8) :: rh ! daily rh @@ -102,35 +110,31 @@ subroutine fire_danger_index ( currentSite, temperature_inst, atm2lnd_inst) real yipsolon; !intermediate varable for dewpoint calculation real dewpoint; !dewpoint in K real d_NI; !daily change in Nesterov Index. C^2 + integer :: iofp ! index of oldest the fates patch - associate( & - t_veg24 => temperature_inst%t_veg24_patch , & ! Input: [real(r8) (:)] avg pft vegetation temperature for last 24 hrs - - prec24 => atm2lnd_inst%prec24_patch , & ! Input: [real(r8) (:)] avg pft rainfall for last 24 hrs - rh24 => atm2lnd_inst%rh24_patch & ! Input: [real(r8) (:)] avg pft relative humidity for last 24 hrs - ) - - ! NOTE: t_veg24(:), prec24(:) and rh24(:) are p level temperatures, precipitation and RH, - ! which probably won't have much inpact, unless we decide to ever calculated the NI for each patch. - - temp_in_C = t_veg24(currentSite%oldest_patch%clm_pno) - tfrz - rainfall = prec24(currentSite%oldest_patch%clm_pno) *24.0_r8*3600._r8 - rh = rh24(currentSite%oldest_patch%clm_pno) - - if (rainfall > 3.0_r8) then !rezero NI if it rains... - d_NI = 0.0_r8 - currentSite%acc_NI = 0.0_r8 - else - yipsolon = (SF_val_fdi_a* temp_in_C)/(SF_val_fdi_b+ temp_in_C)+log(rh/100.0_r8) - dewpoint = (SF_val_fdi_b*yipsolon)/(SF_val_fdi_a-yipsolon) !Standard met. formula - d_NI = ( temp_in_C-dewpoint)* temp_in_C !follows Nesterov 1968. Equation 5. Thonicke et al. 2010. - if (d_NI < 0.0_r8) then !Change in NI cannot be negative. - d_NI = 0.0_r8 !check - endif - endif - currentSite%acc_NI = currentSite%acc_NI + d_NI !Accumulate Nesterov index over the fire season. - - end associate + ! NOTE that the boundary conditions of temperature, precipitation and relative humidity + ! are available at the patch level. We are currently using a simplification where the whole site + ! is simply using the values associated with the first patch. + ! which probably won't have much inpact, unless we decide to ever calculated the NI for each patch. + + iofp = currentSite%oldest_patch%patchno + + temp_in_C = bc_in%t_veg24_pa(iofp) - tfrz + rainfall = bc_in%precip24_pa(iofp)*sec_per_day + rh = bc_in%relhumid24_pa(iofp) + + if (rainfall > 3.0_r8) then !rezero NI if it rains... + d_NI = 0.0_r8 + currentSite%acc_NI = 0.0_r8 + else + yipsolon = (SF_val_fdi_a* temp_in_C)/(SF_val_fdi_b+ temp_in_C)+log(rh/100.0_r8) + dewpoint = (SF_val_fdi_b*yipsolon)/(SF_val_fdi_a-yipsolon) !Standard met. formula + d_NI = ( temp_in_C-dewpoint)* temp_in_C !follows Nesterov 1968. Equation 5. Thonicke et al. 2010. + if (d_NI < 0.0_r8) then !Change in NI cannot be negative. + d_NI = 0.0_r8 !check + endif + endif + currentSite%acc_NI = currentSite%acc_NI + d_NI !Accumulate Nesterov index over the fire season. end subroutine fire_danger_index @@ -179,15 +183,15 @@ subroutine charecteristics_of_fuel ( currentSite ) currentPatch%fuel_frac = 0.0_r8 if(write_sf == 1)then - if (masterproc) write(iulog,*) ' leaf_litter1 ',currentPatch%leaf_litter - if (masterproc) write(iulog,*) ' leaf_litter2 ',sum(currentPatch%CWD_AG) - if (masterproc) write(iulog,*) ' leaf_litter3 ',currentPatch%livegrass - if (masterproc) write(iulog,*) ' sum fuel', currentPatch%sum_fuel + if (masterproc) write(fates_log(),*) ' leaf_litter1 ',currentPatch%leaf_litter + if (masterproc) write(fates_log(),*) ' leaf_litter2 ',sum(currentPatch%CWD_AG) + if (masterproc) write(fates_log(),*) ' leaf_litter3 ',currentPatch%livegrass + if (masterproc) write(fates_log(),*) ' sum fuel', currentPatch%sum_fuel endif currentPatch%sum_fuel = sum(currentPatch%leaf_litter) + sum(currentPatch%CWD_AG) + currentPatch%livegrass if(write_SF == 1)then - if (masterproc) write(iulog,*) 'sum fuel', currentPatch%sum_fuel,currentPatch%area + if (masterproc) write(fates_log(),*) 'sum fuel', currentPatch%sum_fuel,currentPatch%area endif ! =============================================== ! Average moisture, bulk density, surface area-volume and moisture extinction of fuel @@ -199,9 +203,9 @@ subroutine charecteristics_of_fuel ( currentSite ) currentPatch%fuel_frac(dg_sf+1:tr_sf) = currentPatch%CWD_AG / currentPatch%sum_fuel if(write_sf == 1)then - if (masterproc) write(iulog,*) 'ff1 ',currentPatch%fuel_frac - if (masterproc) write(iulog,*) 'ff2 ',currentPatch%fuel_frac - if (masterproc) write(iulog,*) 'ff2a ',lg_sf,currentPatch%livegrass,currentPatch%sum_fuel + if (masterproc) write(fates_log(),*) 'ff1 ',currentPatch%fuel_frac + if (masterproc) write(fates_log(),*) 'ff2 ',currentPatch%fuel_frac + if (masterproc) write(fates_log(),*) 'ff2a ',lg_sf,currentPatch%livegrass,currentPatch%sum_fuel endif currentPatch%fuel_frac(lg_sf) = currentPatch%livegrass / currentPatch%sum_fuel @@ -210,10 +214,10 @@ subroutine charecteristics_of_fuel ( currentSite ) !Equation 6 in Thonicke et al. 2010. fuel_moisture(dg_sf+1:tr_sf) = exp(-1.0_r8 * SF_val_alpha_FMC(dg_sf+1:tr_sf) * currentSite%acc_NI) if(write_SF == 1)then - if (masterproc) write(iulog,*) 'ff3 ',currentPatch%fuel_frac - if (masterproc) write(iulog,*) 'fm ',fuel_moisture - if (masterproc) write(iulog,*) 'csa ',currentSite%acc_NI - if (masterproc) write(iulog,*) 'sfv ',SF_val_alpha_FMC + if (masterproc) write(fates_log(),*) 'ff3 ',currentPatch%fuel_frac + if (masterproc) write(fates_log(),*) 'fm ',fuel_moisture + if (masterproc) write(fates_log(),*) 'csa ',currentSite%acc_NI + if (masterproc) write(fates_log(),*) 'sfv ',SF_val_alpha_FMC endif ! FIX(RF,032414): needs refactoring. ! average water content !is this the correct metric? @@ -227,7 +231,7 @@ subroutine charecteristics_of_fuel ( currentSite ) currentPatch%fuel_mef = sum(currentPatch%fuel_frac(dg_sf:lb_sf) * MEF(dg_sf:lb_sf)) currentPatch%fuel_eff_moist = sum(currentPatch%fuel_frac(dg_sf:lb_sf) * fuel_moisture(dg_sf:lb_sf)) if(write_sf == 1)then - if (masterproc) write(iulog,*) 'ff4 ',currentPatch%fuel_eff_moist + if (masterproc) write(fates_log(),*) 'ff4 ',currentPatch%fuel_eff_moist endif ! Add on properties of live grass multiplied by grass fraction. (6) currentPatch%fuel_bulkd = currentPatch%fuel_bulkd + currentPatch%fuel_frac(lg_sf) * SF_val_FBD(lg_sf) @@ -254,14 +258,14 @@ subroutine charecteristics_of_fuel ( currentSite ) if(write_SF == 1)then - if (masterproc) write(iulog,*) 'no litter fuel at all',currentPatch%patchno, & + if (masterproc) write(fates_log(),*) 'no litter fuel at all',currentPatch%patchno, & currentPatch%sum_fuel,sum(currentPatch%cwd_ag), & sum(currentPatch%cwd_bg),sum(currentPatch%leaf_litter) endif currentPatch%fuel_sav = sum(SF_val_SAV(1:ncwd+2))/(ncwd+2) ! make average sav to avoid crashing code. - if (masterproc) write(iulog,*) 'problem with spitfire fuel averaging' + if (masterproc) write(fates_log(),*) 'problem with spitfire fuel averaging' ! FIX(SPM,032414) refactor...should not have 0 fuel unless everything is burnt ! off. @@ -277,7 +281,7 @@ subroutine charecteristics_of_fuel ( currentSite ) ! FIX(SPM,032414) refactor... if(write_SF == 1.and.currentPatch%fuel_sav <= 0.0_r8.or.currentPatch%fuel_bulkd <= & 0.0_r8.or.currentPatch%fuel_mef <= 0.0_r8.or.currentPatch%fuel_eff_moist <= 0.0_r8)then - if (masterproc) write(iulog,*) 'problem with spitfire fuel averaging' + if (masterproc) write(fates_log(),*) 'problem with spitfire fuel averaging' endif currentPatch => currentPatch%younger @@ -288,33 +292,35 @@ end subroutine charecteristics_of_fuel !***************************************************************** - subroutine wind_effect ( currentSite, atm2lnd_inst) + subroutine wind_effect ( currentSite, bc_in) !*****************************************************************. ! Routine called daily from within ED within a site loop. ! Calculates the effective windspeed based on vegetation charecteristics. + use FatesConstantsMod, only : sec_per_min + type(ed_site_type) , intent(inout), target :: currentSite - type(atm2lnd_type) , intent(in) :: atm2lnd_inst + type(bc_in_type) , intent(in) :: bc_in type(ed_patch_type) , pointer :: currentPatch type(ed_cohort_type), pointer :: currentCohort - ! note - this is a p level temperature, which probably won't have much inpact, - ! unless we decide to ever calculated the NI for each patch. - real(r8), pointer :: wind24(:) - real(r8) :: wind ! daily wind real(r8) :: total_grass_area ! per patch,in m2 real(r8) :: tree_fraction ! site level. no units real(r8) :: grass_fraction ! site level. no units real(r8) :: bare_fraction ! site level. no units + integer :: iofp ! index of oldest fates patch - wind24 => atm2lnd_inst%wind24_patch ! Input: [real(r8) (:)] avg pft windspeed (m/s) + ! note - this is a patch level temperature, which probably won't have much inpact, + ! unless we decide to ever calculated the NI for each patch. + + iofp = currentSite%oldest_patch%patchno + wind = bc_in%wind24_pa(iofp) * sec_per_min ! Convert to m/min for SPITFIRE units. - wind = wind24(currentSite%oldest_patch%clm_pno) * 60._r8 ! Convert to m/min for SPITFIRE units. if(write_SF == 1)then - if (masterproc) write(iulog,*) 'wind24', wind24(currentSite%oldest_patch%clm_pno) + if (masterproc) write(fates_log(),*) 'wind24', wind endif ! --- influence of wind speed, corrected for surface roughness---- ! --- averaged over the whole grid cell to prevent extreme divergence @@ -328,7 +334,7 @@ subroutine wind_effect ( currentSite, atm2lnd_inst) currentCohort => currentPatch%tallest do while(associated(currentCohort)) - write(iulog,*) 'SF currentCohort%c_area ',currentCohort%c_area + write(fates_log(),*) 'SF currentCohort%c_area ',currentCohort%c_area if(pftcon%woody(currentCohort%pft) == 1)then currentPatch%total_tree_area = currentPatch%total_tree_area + currentCohort%c_area else @@ -340,10 +346,10 @@ subroutine wind_effect ( currentSite, atm2lnd_inst) grass_fraction = grass_fraction + min(currentPatch%area,total_grass_area)/AREA if(DEBUG)then - !write(iulog,*) 'SF currentPatch%area ',currentPatch%area - !write(iulog,*) 'SF currentPatch%total_area ',currentPatch%total_tree_area - !write(iulog,*) 'SF total_grass_area ',tree_fraction,grass_fraction - !write(iulog,*) 'SF AREA ',AREA + !write(fates_log(),*) 'SF currentPatch%area ',currentPatch%area + !write(fates_log(),*) 'SF currentPatch%total_area ',currentPatch%total_tree_area + !write(fates_log(),*) 'SF total_grass_area ',tree_fraction,grass_fraction + !write(fates_log(),*) 'SF AREA ',AREA endif currentPatch => currentPatch%younger @@ -353,7 +359,7 @@ subroutine wind_effect ( currentSite, atm2lnd_inst) grass_fraction = min(grass_fraction,1.0_r8-tree_fraction) bare_fraction = 1.0 - tree_fraction - grass_fraction if(write_sf == 1)then - if (masterproc) write(iulog,*) 'grass, trees, bare',grass_fraction, tree_fraction, bare_fraction + if (masterproc) write(fates_log(),*) 'grass, trees, bare',grass_fraction, tree_fraction, bare_fraction endif currentPatch=>currentSite%oldest_patch; @@ -403,18 +409,18 @@ subroutine rate_of_spread ( currentSite ) currentPatch%sum_fuel = currentPatch%sum_fuel * (1.0_r8 - SF_val_miner_total) !net of minerals ! ----start spreading--- - if (masterproc.and.DEBUG) write(iulog,*) 'SF - currentPatch%fuel_bulkd ',currentPatch%fuel_bulkd - if (masterproc.and.DEBUG) write(iulog,*) 'SF - SF_val_part_dens ',SF_val_part_dens + if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - currentPatch%fuel_bulkd ',currentPatch%fuel_bulkd + if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - SF_val_part_dens ',SF_val_part_dens beta = (currentPatch%fuel_bulkd / 0.45_r8) / SF_val_part_dens ! Equation A6 in Thonicke et al. 2010 beta_op = 0.200395_r8 *(currentPatch%fuel_sav**(-0.8189_r8)) - if (masterproc.and.DEBUG) write(iulog,*) 'SF - beta ',beta - if (masterproc.and.DEBUG) write(iulog,*) 'SF - beta_op ',beta_op + if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - beta ',beta + if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - beta_op ',beta_op bet = beta/beta_op if(write_sf == 1)then - if (masterproc) write(iulog,*) 'esf ',currentPatch%fuel_eff_moist + if (masterproc) write(fates_log(),*) 'esf ',currentPatch%fuel_eff_moist endif ! ---heat of pre-ignition--- ! Equation A4 in Thonicke et al. 2010 @@ -432,11 +438,11 @@ subroutine rate_of_spread ( currentSite ) ! Equation A5 in Thonicke et al. 2010 if (DEBUG) then - if (masterproc.and.DEBUG) write(iulog,*) 'SF - c ',c - if (masterproc.and.DEBUG) write(iulog,*) 'SF - currentPatch%effect_wspeed ',currentPatch%effect_wspeed - if (masterproc.and.DEBUG) write(iulog,*) 'SF - b ',b - if (masterproc.and.DEBUG) write(iulog,*) 'SF - bet ',bet - if (masterproc.and.DEBUG) write(iulog,*) 'SF - e ',e + if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - c ',c + if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - currentPatch%effect_wspeed ',currentPatch%effect_wspeed + if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - b ',b + if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - bet ',bet + if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - e ',e endif ! convert from m/min to ft/min for Rothermel ROS eqn @@ -464,18 +470,18 @@ subroutine rate_of_spread ( currentSite ) ! FIX(SPM, 040114) ask RF if this should be an endrun ! if(write_SF == 1)then - ! write(iulog,*) 'moist_damp' ,moist_damp,mw_weight,currentPatch%fuel_eff_moist,currentPatch%fuel_mef + ! write(fates_log(),*) 'moist_damp' ,moist_damp,mw_weight,currentPatch%fuel_eff_moist,currentPatch%fuel_mef ! endif ir = gamma_aptr*(currentPatch%sum_fuel/0.45_r8)*SF_val_fuel_energy*moist_damp*SF_val_miner_damp ! currentPatch%sum_fuel needs to be converted from kgC/m2 to kgBiomass/m2 - ! write(iulog,*) 'ir',gamma_aptr,moist_damp,SF_val_fuel_energy,SF_val_miner_damp + ! write(fates_log(),*) 'ir',gamma_aptr,moist_damp,SF_val_fuel_energy,SF_val_miner_damp if (((currentPatch%fuel_bulkd/0.45_r8) <= 0.0_r8).or.(eps <= 0.0_r8).or.(q_ig <= 0.0_r8)) then currentPatch%ROS_front = 0.0_r8 else ! Equation 9. Thonicke et al. 2010. currentPatch%ROS_front = (ir*xi*(1.0_r8+phi_wind)) / (currentPatch%fuel_bulkd/0.45_r8*eps*q_ig) - ! write(iulog,*) 'ROS',currentPatch%ROS_front,phi_wind,currentPatch%effect_wspeed - ! write(iulog,*) 'ros calcs',currentPatch%fuel_bulkd,ir,xi,eps,q_ig + ! write(fates_log(),*) 'ROS',currentPatch%ROS_front,phi_wind,currentPatch%effect_wspeed + ! write(fates_log(),*) 'ros calcs',currentPatch%fuel_bulkd,ir,xi,eps,q_ig endif ! Equation 10 in Thonicke et al. 2010 ! Can FBP System in m/min @@ -598,7 +604,7 @@ subroutine fire_intensity ( currentSite ) W = currentPatch%TFC_ROS / 0.45_r8 !kgC/m2 to kgbiomass/m2 currentPatch%FI = SF_val_fuel_energy * W * ROS !kj/m/s, or kW/m if(write_sf == 1)then - if(masterproc) write(iulog,*) 'fire_intensity',currentPatch%fi,W,currentPatch%ROS_front + if(masterproc) write(fates_log(),*) 'fire_intensity',currentPatch%fi,W,currentPatch%ROS_front endif !'decide_fire' subroutine shortened and put in here... if (currentPatch%FI >= fire_threshold) then ! 50kW/m is the threshold for a self-sustaining fire @@ -609,7 +615,7 @@ subroutine fire_intensity ( currentSite ) ! Equation 14 in Thonicke et al. 2010 currentPatch%FD = SF_val_max_durat / (1.0_r8 + SF_val_max_durat * exp(SF_val_durat_slope*d_FDI)) if(write_SF == 1)then - if (masterproc) write(iulog,*) 'fire duration minutes',currentPatch%fd + if (masterproc) write(fates_log(),*) 'fire duration minutes',currentPatch%fd endif !equation 15 in Arora and Boer CTEM model.Average fire is 1 day long. !currentPatch%FD = 60.0_r8 * 24.0_r8 !no minutes in a day @@ -684,7 +690,7 @@ subroutine area_burnt ( currentSite ) p = currentPatch%clm_pno g = patch%gridcell(p) - ! g = currentSite%clmgcell (DEPRECATED VARIABLE) + ! INTERF-TODO: ! THIS SHOULD HAVE THE COLUMN AND LU AREA WEIGHT ALSO, NO? @@ -703,18 +709,18 @@ subroutine area_burnt ( currentSite ) currentPatch%AB = size_of_fire * currentPatch%nf if (currentPatch%AB > gridarea*currentPatch%area/area) then !all of patch burnt. - if (masterproc) write(iulog,*) 'burnt all of patch',currentPatch%patchno, & + if (masterproc) write(fates_log(),*) 'burnt all of patch',currentPatch%patchno, & currentPatch%area/area,currentPatch%ab,currentPatch%area/area*gridarea - if (masterproc) write(iulog,*) 'ros',currentPatch%ROS_front,currentPatch%FD, & + if (masterproc) write(fates_log(),*) 'ros',currentPatch%ROS_front,currentPatch%FD, & currentPatch%NF,currentPatch%FI,size_of_fire - if (masterproc) write(iulog,*) 'litter',currentPatch%sum_fuel,currentPatch%CWD_AG,currentPatch%leaf_litter + if (masterproc) write(fates_log(),*) 'litter',currentPatch%sum_fuel,currentPatch%CWD_AG,currentPatch%leaf_litter ! turn km2 into m2. work out total area burnt. currentPatch%AB = currentPatch%area * gridarea/AREA endif currentPatch%frac_burnt = currentPatch%AB / (gridarea*currentPatch%area/area) if(write_SF == 1)then - if (masterproc) write(iulog,*) 'frac_burnt',currentPatch%frac_burnt + if (masterproc) write(fates_log(),*) 'frac_burnt',currentPatch%frac_burnt endif endif endif! fire @@ -771,7 +777,7 @@ subroutine crown_scorching ( currentSite ) currentCohort%bdead))*currentCohort%n)/tree_ag_biomass !equation 16 in Thonicke et al. 2010 if(write_SF == 1)then - if (masterproc) write(iulog,*) 'currentPatch%SH',currentPatch%SH,f_ag_bmass + if (masterproc) write(fates_log(),*) 'currentPatch%SH',currentPatch%SH,f_ag_bmass endif !2/3 Byram (1959) currentPatch%SH = currentPatch%SH + f_ag_bmass * SF_val_alpha_SH * (currentPatch%FI**0.667_r8) diff --git a/components/clm/src/ED/main/EDMainMod.F90 b/components/clm/src/ED/main/EDMainMod.F90 index d43c2be6e4..81c7fac44e 100755 --- a/components/clm/src/ED/main/EDMainMod.F90 +++ b/components/clm/src/ED/main/EDMainMod.F90 @@ -40,9 +40,7 @@ module EDMainMod contains !-------------------------------------------------------------------------------! - subroutine ed_ecosystem_dynamics(currentSite, bc_in, & - atm2lnd_inst, & - temperature_inst) + subroutine ed_ecosystem_dynamics(currentSite, bc_in) ! ! !DESCRIPTION: ! Core of ed model, calling all subsequent vegetation dynamics routines @@ -50,8 +48,6 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, & ! !ARGUMENTS: type(ed_site_type) , intent(inout), target :: currentSite type(bc_in_type) , intent(in) :: bc_in - type(atm2lnd_type) , intent(in) :: atm2lnd_inst - type(temperature_type) , intent(in) :: temperature_inst ! ! !LOCAL VARIABLES: type(ed_patch_type), pointer :: currentPatch @@ -70,7 +66,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in, & call phenology(currentSite, bc_in ) - call fire_model(currentSite, atm2lnd_inst, temperature_inst) + call fire_model(currentSite, bc_in) ! Calculate disturbance and mortality based on previous timestep vegetation. call disturbance_rates(currentSite) diff --git a/components/clm/src/ED/main/FatesInterfaceMod.F90 b/components/clm/src/ED/main/FatesInterfaceMod.F90 index afc39fcdd6..210aece16d 100644 --- a/components/clm/src/ED/main/FatesInterfaceMod.F90 +++ b/components/clm/src/ED/main/FatesInterfaceMod.F90 @@ -80,6 +80,18 @@ module FatesInterfaceMod ! Site volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3] real(r8) :: h2osoi_vol_si + ! Fire Model + + ! Average precipitation over the last 24 hours [mm/s] + real(r8), allocatable :: precip24_pa(:) + + ! Average relative humidity over past 24 hours [-] + real(r8), allocatable :: relhumid24_pa(:) + + ! Patch 24-hour running mean of wind (m/s ?) + real(r8), allocatable :: wind24_pa(:) + + ! Radiation variables for calculating sun/shade fractions ! --------------------------------------------------------------------------------- @@ -378,6 +390,9 @@ subroutine allocate_bcin(bc_in) ! Vegetation Dynamics allocate(bc_in%t_veg24_pa(numPatchesPerCol)) + allocate(bc_in%wind24_pa(numPatchesPerCol)) + allocate(bc_in%relhumid24_pa(numPatchesPerCol)) + allocate(bc_in%precip24_pa(numPatchesPerCol)) ! Radiation allocate(bc_in%solad_parb(numPatchesPerCol,cp_numSWb)) @@ -491,6 +506,9 @@ subroutine zero_bcs(this,s) this%bc_in(s)%t_veg24_si = 0.0_r8 this%bc_in(s)%t_veg24_pa(:) = 0.0_r8 this%bc_in(s)%h2osoi_vol_si = 0.0_r8 + this%bc_in(s)%precip24_pa(:) = 0.0_r8 + this%bc_in(s)%relhumid24_pa(:) = 0.0_r8 + this%bc_in(s)%wind24_pa(:) = 0.0_r8 this%bc_in(s)%solad_parb(:,:) = 0.0_r8 this%bc_in(s)%solai_parb(:,:) = 0.0_r8 diff --git a/components/clm/src/utils/clmfates_interfaceMod.F90 b/components/clm/src/utils/clmfates_interfaceMod.F90 index 7400a69d99..97f32087be 100644 --- a/components/clm/src/utils/clmfates_interfaceMod.F90 +++ b/components/clm/src/utils/clmfates_interfaceMod.F90 @@ -518,17 +518,23 @@ subroutine dynamics_driv(this, nc, bounds_clump, & this%fates(nc)%bc_in(s)%t_veg24_si = & temperature_inst%t_veg24_patch(col%patchi(c)) + do ifp = 1, this%fates(nc)%sites(s)%youngest_patch%patchno p = ifp+col%patchi(c) this%fates(nc)%bc_in(s)%t_veg24_pa(ifp) = & temperature_inst%t_veg24_patch(p) - end do - end do + this%fates(nc)%bc_in(s)%precip24_pa(ifp) = & + atm2lnd_inst%prec24_patch(p) - ! TODO-INTER: PROCEDURE FOR CONVERTING CLM/ALM FIELDS TO MODEL BOUNDARY - ! CONDITIONS. IE. + this%fates(nc)%bc_in(s)%relhumid24_pa(ifp) = & + atm2lnd_inst%rh24_patch(p) + this%fates(nc)%bc_in(s)%wind24_pa(ifp) = & + atm2lnd_inst%wind24_patch(p) + + end do + end do ! where most things happen do s = 1,this%fates(nc)%nsites From ff8cf21ae4e7cc56feae80b3fd499e7cd27b9f4d Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 11 Jan 2017 17:15:18 -0800 Subject: [PATCH 06/11] Added masterproc (used in firemod for log messaging) to the control parameters passed during initialization. Tracked down other instances of use, and also fixed calls to clm version of r8 or iulog when found. --- components/clm/src/ED/fire/SFMainMod.F90 | 74 ++++++++++--------- components/clm/src/ED/main/EDInitMod.F90 | 12 +-- components/clm/src/ED/main/EDMainMod.F90 | 4 +- components/clm/src/ED/main/EDTypesMod.F90 | 6 ++ .../clm/src/ED/main/FatesInterfaceMod.F90 | 51 ++++++++----- .../clm/src/utils/clmfates_interfaceMod.F90 | 11 ++- 6 files changed, 89 insertions(+), 69 deletions(-) diff --git a/components/clm/src/ED/fire/SFMainMod.F90 b/components/clm/src/ED/fire/SFMainMod.F90 index 1f5f14f235..e5557d1ff2 100755 --- a/components/clm/src/ED/fire/SFMainMod.F90 +++ b/components/clm/src/ED/fire/SFMainMod.F90 @@ -7,7 +7,8 @@ module SFMainMod use FatesConstantsMod , only : r8 => fates_r8 - use spmdMod , only : masterproc +! use spmdMod , only : masterproc + use EDTypesMod , only : cp_masterproc ! 1= master process, 0=not master process use FatesGlobals , only : fates_log use FatesInterfaceMod , only : bc_in_type @@ -183,15 +184,15 @@ subroutine charecteristics_of_fuel ( currentSite ) currentPatch%fuel_frac = 0.0_r8 if(write_sf == 1)then - if (masterproc) write(fates_log(),*) ' leaf_litter1 ',currentPatch%leaf_litter - if (masterproc) write(fates_log(),*) ' leaf_litter2 ',sum(currentPatch%CWD_AG) - if (masterproc) write(fates_log(),*) ' leaf_litter3 ',currentPatch%livegrass - if (masterproc) write(fates_log(),*) ' sum fuel', currentPatch%sum_fuel + if ( cp_masterproc == 1 ) write(fates_log(),*) ' leaf_litter1 ',currentPatch%leaf_litter + if ( cp_masterproc == 1 ) write(fates_log(),*) ' leaf_litter2 ',sum(currentPatch%CWD_AG) + if ( cp_masterproc == 1 ) write(fates_log(),*) ' leaf_litter3 ',currentPatch%livegrass + if ( cp_masterproc == 1 ) write(fates_log(),*) ' sum fuel', currentPatch%sum_fuel endif currentPatch%sum_fuel = sum(currentPatch%leaf_litter) + sum(currentPatch%CWD_AG) + currentPatch%livegrass if(write_SF == 1)then - if (masterproc) write(fates_log(),*) 'sum fuel', currentPatch%sum_fuel,currentPatch%area + if ( cp_masterproc == 1 ) write(fates_log(),*) 'sum fuel', currentPatch%sum_fuel,currentPatch%area endif ! =============================================== ! Average moisture, bulk density, surface area-volume and moisture extinction of fuel @@ -203,9 +204,9 @@ subroutine charecteristics_of_fuel ( currentSite ) currentPatch%fuel_frac(dg_sf+1:tr_sf) = currentPatch%CWD_AG / currentPatch%sum_fuel if(write_sf == 1)then - if (masterproc) write(fates_log(),*) 'ff1 ',currentPatch%fuel_frac - if (masterproc) write(fates_log(),*) 'ff2 ',currentPatch%fuel_frac - if (masterproc) write(fates_log(),*) 'ff2a ',lg_sf,currentPatch%livegrass,currentPatch%sum_fuel + if ( cp_masterproc == 1 ) write(fates_log(),*) 'ff1 ',currentPatch%fuel_frac + if ( cp_masterproc == 1 ) write(fates_log(),*) 'ff2 ',currentPatch%fuel_frac + if ( cp_masterproc == 1 ) write(fates_log(),*) 'ff2a ',lg_sf,currentPatch%livegrass,currentPatch%sum_fuel endif currentPatch%fuel_frac(lg_sf) = currentPatch%livegrass / currentPatch%sum_fuel @@ -214,10 +215,10 @@ subroutine charecteristics_of_fuel ( currentSite ) !Equation 6 in Thonicke et al. 2010. fuel_moisture(dg_sf+1:tr_sf) = exp(-1.0_r8 * SF_val_alpha_FMC(dg_sf+1:tr_sf) * currentSite%acc_NI) if(write_SF == 1)then - if (masterproc) write(fates_log(),*) 'ff3 ',currentPatch%fuel_frac - if (masterproc) write(fates_log(),*) 'fm ',fuel_moisture - if (masterproc) write(fates_log(),*) 'csa ',currentSite%acc_NI - if (masterproc) write(fates_log(),*) 'sfv ',SF_val_alpha_FMC + if ( cp_masterproc == 1 ) write(fates_log(),*) 'ff3 ',currentPatch%fuel_frac + if ( cp_masterproc == 1 ) write(fates_log(),*) 'fm ',fuel_moisture + if ( cp_masterproc == 1 ) write(fates_log(),*) 'csa ',currentSite%acc_NI + if ( cp_masterproc == 1 ) write(fates_log(),*) 'sfv ',SF_val_alpha_FMC endif ! FIX(RF,032414): needs refactoring. ! average water content !is this the correct metric? @@ -231,7 +232,7 @@ subroutine charecteristics_of_fuel ( currentSite ) currentPatch%fuel_mef = sum(currentPatch%fuel_frac(dg_sf:lb_sf) * MEF(dg_sf:lb_sf)) currentPatch%fuel_eff_moist = sum(currentPatch%fuel_frac(dg_sf:lb_sf) * fuel_moisture(dg_sf:lb_sf)) if(write_sf == 1)then - if (masterproc) write(fates_log(),*) 'ff4 ',currentPatch%fuel_eff_moist + if ( cp_masterproc == 1 ) write(fates_log(),*) 'ff4 ',currentPatch%fuel_eff_moist endif ! Add on properties of live grass multiplied by grass fraction. (6) currentPatch%fuel_bulkd = currentPatch%fuel_bulkd + currentPatch%fuel_frac(lg_sf) * SF_val_FBD(lg_sf) @@ -258,14 +259,14 @@ subroutine charecteristics_of_fuel ( currentSite ) if(write_SF == 1)then - if (masterproc) write(fates_log(),*) 'no litter fuel at all',currentPatch%patchno, & + if ( cp_masterproc == 1 ) write(fates_log(),*) 'no litter fuel at all',currentPatch%patchno, & currentPatch%sum_fuel,sum(currentPatch%cwd_ag), & sum(currentPatch%cwd_bg),sum(currentPatch%leaf_litter) endif currentPatch%fuel_sav = sum(SF_val_SAV(1:ncwd+2))/(ncwd+2) ! make average sav to avoid crashing code. - if (masterproc) write(fates_log(),*) 'problem with spitfire fuel averaging' + if ( cp_masterproc == 1 ) write(fates_log(),*) 'problem with spitfire fuel averaging' ! FIX(SPM,032414) refactor...should not have 0 fuel unless everything is burnt ! off. @@ -281,7 +282,7 @@ subroutine charecteristics_of_fuel ( currentSite ) ! FIX(SPM,032414) refactor... if(write_SF == 1.and.currentPatch%fuel_sav <= 0.0_r8.or.currentPatch%fuel_bulkd <= & 0.0_r8.or.currentPatch%fuel_mef <= 0.0_r8.or.currentPatch%fuel_eff_moist <= 0.0_r8)then - if (masterproc) write(fates_log(),*) 'problem with spitfire fuel averaging' + if ( cp_masterproc == 1 ) write(fates_log(),*) 'problem with spitfire fuel averaging' endif currentPatch => currentPatch%younger @@ -320,7 +321,7 @@ subroutine wind_effect ( currentSite, bc_in) wind = bc_in%wind24_pa(iofp) * sec_per_min ! Convert to m/min for SPITFIRE units. if(write_SF == 1)then - if (masterproc) write(fates_log(),*) 'wind24', wind + if ( cp_masterproc == 1 ) write(fates_log(),*) 'wind24', wind endif ! --- influence of wind speed, corrected for surface roughness---- ! --- averaged over the whole grid cell to prevent extreme divergence @@ -359,7 +360,7 @@ subroutine wind_effect ( currentSite, bc_in) grass_fraction = min(grass_fraction,1.0_r8-tree_fraction) bare_fraction = 1.0 - tree_fraction - grass_fraction if(write_sf == 1)then - if (masterproc) write(fates_log(),*) 'grass, trees, bare',grass_fraction, tree_fraction, bare_fraction + if ( cp_masterproc == 1 ) write(fates_log(),*) 'grass, trees, bare',grass_fraction, tree_fraction, bare_fraction endif currentPatch=>currentSite%oldest_patch; @@ -409,18 +410,18 @@ subroutine rate_of_spread ( currentSite ) currentPatch%sum_fuel = currentPatch%sum_fuel * (1.0_r8 - SF_val_miner_total) !net of minerals ! ----start spreading--- - if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - currentPatch%fuel_bulkd ',currentPatch%fuel_bulkd - if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - SF_val_part_dens ',SF_val_part_dens + if ( cp_masterproc == 1 .and.DEBUG) write(fates_log(),*) 'SF - currentPatch%fuel_bulkd ',currentPatch%fuel_bulkd + if ( cp_masterproc == 1 .and.DEBUG) write(fates_log(),*) 'SF - SF_val_part_dens ',SF_val_part_dens beta = (currentPatch%fuel_bulkd / 0.45_r8) / SF_val_part_dens ! Equation A6 in Thonicke et al. 2010 beta_op = 0.200395_r8 *(currentPatch%fuel_sav**(-0.8189_r8)) - if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - beta ',beta - if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - beta_op ',beta_op + if ( cp_masterproc == 1 .and.DEBUG) write(fates_log(),*) 'SF - beta ',beta + if ( cp_masterproc == 1 .and.DEBUG) write(fates_log(),*) 'SF - beta_op ',beta_op bet = beta/beta_op if(write_sf == 1)then - if (masterproc) write(fates_log(),*) 'esf ',currentPatch%fuel_eff_moist + if ( cp_masterproc == 1 ) write(fates_log(),*) 'esf ',currentPatch%fuel_eff_moist endif ! ---heat of pre-ignition--- ! Equation A4 in Thonicke et al. 2010 @@ -438,11 +439,11 @@ subroutine rate_of_spread ( currentSite ) ! Equation A5 in Thonicke et al. 2010 if (DEBUG) then - if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - c ',c - if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - currentPatch%effect_wspeed ',currentPatch%effect_wspeed - if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - b ',b - if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - bet ',bet - if (masterproc.and.DEBUG) write(fates_log(),*) 'SF - e ',e + if ( cp_masterproc == 1 .and.DEBUG) write(fates_log(),*) 'SF - c ',c + if ( cp_masterproc == 1 .and.DEBUG) write(fates_log(),*) 'SF - currentPatch%effect_wspeed ',currentPatch%effect_wspeed + if ( cp_masterproc == 1 .and.DEBUG) write(fates_log(),*) 'SF - b ',b + if ( cp_masterproc == 1 .and.DEBUG) write(fates_log(),*) 'SF - bet ',bet + if ( cp_masterproc == 1 .and.DEBUG) write(fates_log(),*) 'SF - e ',e endif ! convert from m/min to ft/min for Rothermel ROS eqn @@ -604,7 +605,7 @@ subroutine fire_intensity ( currentSite ) W = currentPatch%TFC_ROS / 0.45_r8 !kgC/m2 to kgbiomass/m2 currentPatch%FI = SF_val_fuel_energy * W * ROS !kj/m/s, or kW/m if(write_sf == 1)then - if(masterproc) write(fates_log(),*) 'fire_intensity',currentPatch%fi,W,currentPatch%ROS_front + if( cp_masterproc == 1 ) write(fates_log(),*) 'fire_intensity',currentPatch%fi,W,currentPatch%ROS_front endif !'decide_fire' subroutine shortened and put in here... if (currentPatch%FI >= fire_threshold) then ! 50kW/m is the threshold for a self-sustaining fire @@ -615,7 +616,7 @@ subroutine fire_intensity ( currentSite ) ! Equation 14 in Thonicke et al. 2010 currentPatch%FD = SF_val_max_durat / (1.0_r8 + SF_val_max_durat * exp(SF_val_durat_slope*d_FDI)) if(write_SF == 1)then - if (masterproc) write(fates_log(),*) 'fire duration minutes',currentPatch%fd + if ( cp_masterproc == 1 ) write(fates_log(),*) 'fire duration minutes',currentPatch%fd endif !equation 15 in Arora and Boer CTEM model.Average fire is 1 day long. !currentPatch%FD = 60.0_r8 * 24.0_r8 !no minutes in a day @@ -709,18 +710,19 @@ subroutine area_burnt ( currentSite ) currentPatch%AB = size_of_fire * currentPatch%nf if (currentPatch%AB > gridarea*currentPatch%area/area) then !all of patch burnt. - if (masterproc) write(fates_log(),*) 'burnt all of patch',currentPatch%patchno, & + if ( cp_masterproc == 1 ) write(fates_log(),*) 'burnt all of patch',currentPatch%patchno, & currentPatch%area/area,currentPatch%ab,currentPatch%area/area*gridarea - if (masterproc) write(fates_log(),*) 'ros',currentPatch%ROS_front,currentPatch%FD, & + if ( cp_masterproc == 1 ) write(fates_log(),*) 'ros',currentPatch%ROS_front,currentPatch%FD, & currentPatch%NF,currentPatch%FI,size_of_fire - if (masterproc) write(fates_log(),*) 'litter',currentPatch%sum_fuel,currentPatch%CWD_AG,currentPatch%leaf_litter + if ( cp_masterproc == 1 ) write(fates_log(),*) 'litter', & + currentPatch%sum_fuel,currentPatch%CWD_AG,currentPatch%leaf_litter ! turn km2 into m2. work out total area burnt. currentPatch%AB = currentPatch%area * gridarea/AREA endif currentPatch%frac_burnt = currentPatch%AB / (gridarea*currentPatch%area/area) if(write_SF == 1)then - if (masterproc) write(fates_log(),*) 'frac_burnt',currentPatch%frac_burnt + if ( cp_masterproc == 1 ) write(fates_log(),*) 'frac_burnt',currentPatch%frac_burnt endif endif endif! fire @@ -777,7 +779,7 @@ subroutine crown_scorching ( currentSite ) currentCohort%bdead))*currentCohort%n)/tree_ag_biomass !equation 16 in Thonicke et al. 2010 if(write_SF == 1)then - if (masterproc) write(fates_log(),*) 'currentPatch%SH',currentPatch%SH,f_ag_bmass + if ( cp_masterproc == 1 ) write(fates_log(),*) 'currentPatch%SH',currentPatch%SH,f_ag_bmass endif !2/3 Byram (1959) currentPatch%SH = currentPatch%SH + f_ag_bmass * SF_val_alpha_SH * (currentPatch%FI**0.667_r8) diff --git a/components/clm/src/ED/main/EDInitMod.F90 b/components/clm/src/ED/main/EDInitMod.F90 index 18a52c16d6..0b8f0ef693 100755 --- a/components/clm/src/ED/main/EDInitMod.F90 +++ b/components/clm/src/ED/main/EDInitMod.F90 @@ -4,16 +4,12 @@ module EDInitMod ! Contains all modules to set up the ED structure. ! ============================================================================ - use shr_kind_mod , only : r8 => shr_kind_r8; - use spmdMod , only : masterproc - use decompMod , only : bounds_type + use FatesConstantsMod , only : r8 => fates_r8 use abortutils , only : endrun use EDTypesMod , only : cp_nclmax - use clm_varctl , only : iulog, use_ed_spit_fire + use FatesGlobals , only : fates_log + use clm_varctl , only : use_ed_spit_fire use clm_time_manager , only : is_restart - use CanopyStateType , only : canopystate_type - use WaterStateType , only : waterstate_type - use GridcellType , only : grc use pftconMod , only : pftcon use EDEcophysConType , only : EDecophyscon use EDGrowthFunctionsMod , only : bdead, bleaf, dbh @@ -285,7 +281,7 @@ subroutine init_cohorts( patch_in ) cstatus = patch_in%siteptr%dstatus endif - if ( DEBUG ) write(iulog,*) 'EDInitMod.F90 call create_cohort ' + if ( DEBUG ) write(fates_log(),*) 'EDInitMod.F90 call create_cohort ' call create_cohort(patch_in, pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, & temp_cohort%balive, temp_cohort%bdead, temp_cohort%bstore, & diff --git a/components/clm/src/ED/main/EDMainMod.F90 b/components/clm/src/ED/main/EDMainMod.F90 index 81c7fac44e..e82dc7b512 100755 --- a/components/clm/src/ED/main/EDMainMod.F90 +++ b/components/clm/src/ED/main/EDMainMod.F90 @@ -17,7 +17,7 @@ module EDMainMod use EDtypesMod , only : ncwd, numpft_ed, udata use EDtypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use FatesInterfaceMod , only : bc_in_type - use spmdMod , only : masterproc + use EDTypesMod , only : cp_masterproc implicit none private @@ -53,7 +53,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in) type(ed_patch_type), pointer :: currentPatch !----------------------------------------------------------------------- - if ( masterproc ) write(iulog,*) 'modelday',bc_in%model_day + if ( cp_masterproc==1 ) write(iulog,*) 'modelday',bc_in%model_day !************************************************************************** ! Fire, growth, biogeochemistry. diff --git a/components/clm/src/ED/main/EDTypesMod.F90 b/components/clm/src/ED/main/EDTypesMod.F90 index 4d80cf314c..842e5a8a63 100755 --- a/components/clm/src/ED/main/EDTypesMod.F90 +++ b/components/clm/src/ED/main/EDTypesMod.F90 @@ -149,6 +149,12 @@ module EDTypesMod ! HLM will interpret that the value should not be included in the average. real(r8) :: cp_hio_ignore_val + + ! Is this the master processor, typically useful for knowing if + ! the current machine should be printing out messages to the logs or terminals + ! 1 = TRUE (is master) 0 = FALSE (is not master) + integer :: cp_masterproc + !************************************ !** COHORT type structure ** !************************************ diff --git a/components/clm/src/ED/main/FatesInterfaceMod.F90 b/components/clm/src/ED/main/FatesInterfaceMod.F90 index 210aece16d..ed66e4f9d9 100644 --- a/components/clm/src/ED/main/FatesInterfaceMod.F90 +++ b/components/clm/src/ED/main/FatesInterfaceMod.F90 @@ -15,20 +15,19 @@ module FatesInterfaceMod ! PUBLIC API!!!! ! ------------------------------------------------------------------------------------ - use EDtypesMod , only : ed_site_type, & - numPatchesPerCol, & - cp_nclmax, & - cp_numSWb, & - cp_numlevgrnd, & - cp_maxSWb, & - cp_numlevdecomp, & - cp_numlevdecomp_full, & - cp_hlm_name, & - cp_hio_ignore_val, & - cp_numlevsoil - - use shr_kind_mod , only : r8 => shr_kind_r8 ! INTERF-TODO: REMOVE THIS - + use EDtypesMod , only : ed_site_type + use EDtypesMod , only : numPatchesPerCol + use EDtypesMod , only : cp_nclmax + use EDtypesMod , only : cp_numSWb + use EDtypesMod , only : cp_numlevgrnd + use EDtypesMod , only : cp_maxSWb + use EDtypesMod , only : cp_numlevdecomp + use EDtypesMod , only : cp_numlevdecomp_full + use EDtypesMod , only : cp_hlm_name + use EDtypesMod , only : cp_hio_ignore_val + use EDtypesMod , only : cp_numlevsoil + use EDtypesMod , only : cp_masterproc + use FatesConstantsMod , only : r8 => fates_r8 implicit none @@ -41,6 +40,7 @@ module FatesInterfaceMod ! (Intel-Forum Post), ALLOCATABLES are better perfomance wise as long as they point ! to contiguous memory spaces and do not alias other variables, the case here. ! Naming conventions: _gl means ground layer dimensions + ! _si means site dimensions (scalar in that case) ! _pa means patch dimensions ! _rb means radiation band ! ------------------------------------------------------------------------------------ @@ -620,6 +620,7 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) cp_numlevdecomp = unset_int cp_hlm_name = 'unset' cp_hio_ignore_val = unset_double + cp_masterproc = unset_int case('check_allset') @@ -631,6 +632,14 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) ! end_run('MESSAGE') end if + if(cp_masterproc .eq. unset_int) then + if (fates_global_verbose()) then + write(fates_log(), *) 'FATES parameter unset: cp_masterproc' + end if + ! INTERF-TODO: FATES NEEDS INTERNAL end_run + ! end_run('MESSAGE') + end if + if(cp_numSWb > cp_maxSWb) then if (fates_global_verbose()) then write(fates_log(), *) 'FATES sets a maximum number of shortwave bands' @@ -701,36 +710,38 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval) if(present(ival))then select case (trim(tag)) - - case('num_sw_bbands') + case('masterproc') + cp_masterproc = ival + if (fates_global_verbose()) then + write(fates_log(),*) 'Transfering masterproc = ',ival,' to FATES' + end if + + case('num_sw_bbands') cp_numSwb = ival if (fates_global_verbose()) then write(fates_log(),*) 'Transfering num_sw_bbands = ',ival,' to FATES' end if case('num_lev_ground') - cp_numlevgrnd = ival if (fates_global_verbose()) then write(fates_log(),*) 'Transfering num_lev_ground = ',ival,' to FATES' end if + case('num_lev_soil') - cp_numlevsoil = ival if (fates_global_verbose()) then write(fates_log(),*) 'Transfering num_lev_ground = ',ival,' to FATES' end if case('num_levdecomp_full') - cp_numlevdecomp_full = ival if (fates_global_verbose()) then write(fates_log(),*) 'Transfering num_levdecomp_full = ',ival,' to FATES' end if case('num_levdecomp') - cp_numlevdecomp = ival if (fates_global_verbose()) then write(fates_log(),*) 'Transfering num_levdecomp = ',ival,' to FATES' diff --git a/components/clm/src/utils/clmfates_interfaceMod.F90 b/components/clm/src/utils/clmfates_interfaceMod.F90 index 97f32087be..95e4999fbc 100644 --- a/components/clm/src/utils/clmfates_interfaceMod.F90 +++ b/components/clm/src/utils/clmfates_interfaceMod.F90 @@ -205,6 +205,7 @@ subroutine init(this, bounds_proc, use_ed) ! local variables integer :: nclumps ! Number of threads logical :: verbose_output + integer :: pass_masterproc if (use_ed) then @@ -246,6 +247,12 @@ subroutine init(this, bounds_proc, use_ed) call set_fates_ctrlparms('num_levdecomp_full',ival=nlevdecomp_full) call set_fates_ctrlparms('hlm_name',cval='CLM') call set_fates_ctrlparms('hio_ignore_val',rval=spval) + if(masterproc)then + pass_masterproc = 1 + else + pass_masterproc = 0 + end if + call set_fates_ctrlparms('masterproc',ival=pass_masterproc) ! Check through FATES parameters to see if all have been set call set_fates_ctrlparms('check_allset') @@ -540,9 +547,7 @@ subroutine dynamics_driv(this, nc, bounds_clump, & do s = 1,this%fates(nc)%nsites call ed_ecosystem_dynamics(this%fates(nc)%sites(s), & - this%fates(nc)%bc_in(s), & - atm2lnd_inst, & - temperature_inst ) + this%fates(nc)%bc_in(s)) call ed_update_site(this%fates(nc)%sites(s)) From 390aeb646c93a6c01200115d96b25e5008f8da8f Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 11 Jan 2017 17:25:50 -0800 Subject: [PATCH 07/11] Missed a fix on the merge with master: in FatesInterfaceMod, allocation of boundary conditions was using numPatchesPerCohort while master changed that to maxPatchesPerCohort. --- components/clm/src/ED/main/FatesInterfaceMod.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/components/clm/src/ED/main/FatesInterfaceMod.F90 b/components/clm/src/ED/main/FatesInterfaceMod.F90 index ec9ae4ac71..a8bf0cc6a2 100644 --- a/components/clm/src/ED/main/FatesInterfaceMod.F90 +++ b/components/clm/src/ED/main/FatesInterfaceMod.F90 @@ -388,11 +388,11 @@ subroutine allocate_bcin(bc_in) ! Allocate input boundaries ! Vegetation Dynamics - allocate(bc_in%t_veg24_pa(numPatchesPerCol)) + allocate(bc_in%t_veg24_pa(maxPatchesPerCol)) - allocate(bc_in%wind24_pa(numPatchesPerCol)) - allocate(bc_in%relhumid24_pa(numPatchesPerCol)) - allocate(bc_in%precip24_pa(numPatchesPerCol)) + allocate(bc_in%wind24_pa(maxPatchesPerCol)) + allocate(bc_in%relhumid24_pa(maxPatchesPerCol)) + allocate(bc_in%precip24_pa(maxPatchesPerCol)) ! Radiation allocate(bc_in%solad_parb(maxPatchesPerCol,cp_numSWb)) From d04ba37484c9e524f2751af7ec6eef9058c40fdf Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 11 Jan 2017 19:04:35 -0800 Subject: [PATCH 08/11] Partial re-working of time-control from host to FATES. --- .../ED/biogeochem/EDCanopyStructureMod.F90 | 6 +- .../src/ED/biogeochem/EDCohortDynamicsMod.F90 | 10 ++-- .../ED/biogeochem/EDGrowthFunctionsMod.F90 | 34 ++++++------ .../clm/src/ED/biogeochem/EDPhysiologyMod.F90 | 18 +++--- components/clm/src/ED/main/EDMainMod.F90 | 55 +++++++++++-------- components/clm/src/ED/main/EDTypesMod.F90 | 3 +- .../clm/src/ED/main/FatesInterfaceMod.F90 | 10 +++- .../clm/src/utils/clmfates_interfaceMod.F90 | 26 +++++++-- 8 files changed, 98 insertions(+), 64 deletions(-) diff --git a/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 b/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 index f5419cedd7..00a969a78c 100755 --- a/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 +++ b/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 @@ -465,12 +465,12 @@ subroutine canopy_structure( currentSite ) endif !call terminate_cohorts(currentPatch) if(promswitch == 1)then - ! write(fates_log(),*) 'cohort loop',currentCohort%pft,currentCohort%indexnumber,currentPatch%patchno + ! write(fates_log(),*) 'cohort loop',currentCohort%pft,currentPatch%patchno endif !----------- End of cohort splitting ------------------------------! else if(promswitch == 1)then - ! write(fates_log(),*) 'cohort list',currentCohort%pft,currentCohort%indexnumber, & + ! write(fates_log(),*) 'cohort list',currentCohort%pft, & ! currentCohort%canopy_layer,currentCohort%c_area endif endif @@ -485,7 +485,7 @@ subroutine canopy_structure( currentSite ) !currentPatch%patchno,z,i,lower_cohort_switch endif if(promswitch == 1.and.associated(currentPatch%tallest))then - ! write(fates_log(),*) 'cohorts',currentCohort%pft,currentCohort%indexnumber,currentPatch%patchno, & + ! write(fates_log(),*) 'cohorts',currentCohort%pft,currentPatch%patchno, & !currentCohort%c_area endif enddo !arealayer loop diff --git a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 index 2ccea2cad1..70b00c3307 100755 --- a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 +++ b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 @@ -7,6 +7,7 @@ module EDCohortDynamicsMod use abortutils , only : endrun use FatesGlobals , only : fates_log use FatesConstantsMod , only : r8 => fates_r8 + use FatesConstantsMod , only : fates_unset_int use shr_log_mod , only : errMsg => shr_log_errMsg use pftconMod , only : pftcon use EDEcophysContype , only : EDecophyscon @@ -73,7 +74,6 @@ 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. 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. @@ -82,7 +82,8 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, & ! Define cohort state variable !**********************/ - new_cohort%indexnumber = udata%cohort_number + new_cohort%indexnumber = fates_unset_int ! Cohort indexing was not thread-safe, setting + ! bogus value for the time being (RGK-012017) new_cohort%siteptr => patchptr%siteptr new_cohort%patchptr => patchptr new_cohort%pft = pft @@ -109,7 +110,7 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, & if (new_cohort%dbh <= 0.0_r8 .or. new_cohort%n == 0._r8 .or. new_cohort%pft == 0 ) then write(fates_log(),*) 'ED: something is zero in create_cohort', & - new_cohort%indexnumber,new_cohort%dbh,new_cohort%n, & + new_cohort%dbh,new_cohort%n, & new_cohort%pft call endrun(msg=errMsg(sourcefile, __LINE__)) endif @@ -1006,8 +1007,7 @@ subroutine copy_cohort( currentCohort,copyc ) o => currentCohort n => copyc - udata%cohort_number = udata%cohort_number + 1 - n%indexnumber = udata%cohort_number + n%indexnumber = fates_unset_int ! VEGETATION STRUCTURE n%pft = o%pft diff --git a/components/clm/src/ED/biogeochem/EDGrowthFunctionsMod.F90 b/components/clm/src/ED/biogeochem/EDGrowthFunctionsMod.F90 index a400f46ab9..12a46c79b1 100755 --- a/components/clm/src/ED/biogeochem/EDGrowthFunctionsMod.F90 +++ b/components/clm/src/ED/biogeochem/EDGrowthFunctionsMod.F90 @@ -7,7 +7,7 @@ module EDGrowthFunctionsMod ! ============================================================================ use shr_kind_mod , only : r8 => shr_kind_r8 - use clm_varctl , only : iulog + use FatesGlobals , only : fates_log use pftconMod , only : pftcon use EDEcophysContype , only : EDecophyscon use EDTypesMod , only : ed_cohort_type, cp_nlevcan, dinc_ed @@ -76,7 +76,7 @@ real(r8) function Hite( cohort_in ) c = 0.37_r8 if(cohort_in%dbh <= 0._r8)then - write(iulog,*) 'ED: dbh less than zero problem!',cohort_in%indexnumber + write(fates_log(),*) 'ED: dbh less than zero problem!' cohort_in%dbh = 0.1_r8 endif @@ -106,7 +106,7 @@ real(r8) function Bleaf( cohort_in ) real(r8) :: slascaler ! changes the target biomass according to the SLA if(cohort_in%dbh < 0._r8.or.cohort_in%pft == 0.or.cohort_in%dbh > 1000.0_r8)then - write(iulog,*) 'problems in bleaf',cohort_in%dbh,cohort_in%pft + write(fates_log(),*) 'problems in bleaf',cohort_in%dbh,cohort_in%pft endif if(cohort_in%dbh <= EDecophyscon%max_dbh(cohort_in%pft))then @@ -117,7 +117,7 @@ real(r8) function Bleaf( cohort_in ) slascaler = 0.03_r8/pftcon%slatop(cohort_in%pft) bleaf = bleaf * slascaler - !write(iulog,*) 'bleaf',bleaf, slascaler,cohort_in%pft + !write(fates_log(),*) 'bleaf',bleaf, slascaler,cohort_in%pft !Adjust for canopies that have become so deep that their bottom layer is not producing any carbon... !nb this will change the allometry and the effects of this remain untested. RF. April 2014 @@ -141,7 +141,7 @@ real(r8) function tree_lai( cohort_in ) real(r8) :: slat ! the sla of the top leaf layer. m2/kgC if( cohort_in%bl < 0._r8 .or. cohort_in%pft == 0 ) then - write(iulog,*) 'problem in treelai',cohort_in%bl,cohort_in%pft + write(fates_log(),*) 'problem in treelai',cohort_in%bl,cohort_in%pft endif if( cohort_in%status_coh == 2 ) then ! are the leaves on? @@ -162,7 +162,7 @@ real(r8) function tree_lai( cohort_in ) ! at the moments cp_nlevcan default is 40, which is very large, so exceeding this would clearly illustrate a ! huge error if(cohort_in%treelai > cp_nlevcan*dinc_ed)then - write(iulog,*) 'too much lai' , cohort_in%treelai , cohort_in%pft , cp_nlevcan * dinc_ed + write(fates_log(),*) 'too much lai' , cohort_in%treelai , cohort_in%pft , cp_nlevcan * dinc_ed endif return @@ -186,7 +186,7 @@ real(r8) function tree_sai( cohort_in ) sai_scaler = 0.05_r8 ! here, a high biomass of 20KgC per m2 gives us a high SAI of 1.0. if( cohort_in%bdead < 0._r8 .or. cohort_in%pft == 0 ) then - write(iulog,*) 'problem in treesai',cohort_in%bdead,cohort_in%pft + write(fates_log(),*) 'problem in treesai',cohort_in%bdead,cohort_in%pft endif cohort_in%c_area = c_area(cohort_in) ! call the tree area @@ -199,7 +199,7 @@ real(r8) function tree_sai( cohort_in ) ! at the moments cp_nlevcan default is 40, which is very large, so exceeding this would clearly illustrate a ! huge error if(cohort_in%treesai > cp_nlevcan*dinc_ed)then - write(iulog,*) 'too much sai' , cohort_in%treesai , cohort_in%pft , cp_nlevcan * dinc_ed + write(fates_log(),*) 'too much sai' , cohort_in%treesai , cohort_in%pft , cp_nlevcan * dinc_ed endif return @@ -223,13 +223,13 @@ real(r8) function c_area( cohort_in ) real(r8) :: dbh ! Tree diameter at breat height. cm. if (DEBUG_growth) then - write(iulog,*) 'z_area 1',cohort_in%dbh,cohort_in%pft - write(iulog,*) 'z_area 2',EDecophyscon%max_dbh - write(iulog,*) 'z_area 3',pftcon%woody - write(iulog,*) 'z_area 4',cohort_in%n - write(iulog,*) 'z_area 5',cohort_in%patchptr%spread - write(iulog,*) 'z_area 6',cohort_in%canopy_layer - write(iulog,*) 'z_area 7',ED_val_grass_spread + write(fates_log(),*) 'z_area 1',cohort_in%dbh,cohort_in%pft + write(fates_log(),*) 'z_area 2',EDecophyscon%max_dbh + write(fates_log(),*) 'z_area 3',pftcon%woody + write(fates_log(),*) 'z_area 4',cohort_in%n + write(fates_log(),*) 'z_area 5',cohort_in%patchptr%spread + write(fates_log(),*) 'z_area 6',cohort_in%canopy_layer + write(fates_log(),*) 'z_area 7',ED_val_grass_spread end if dbh = min(cohort_in%dbh,EDecophyscon%max_dbh(cohort_in%pft)) @@ -371,8 +371,8 @@ subroutine mortality_rates( cohort_in,cmort,hmort,bmort ) endif else - write(iulog,*) 'dbh problem in mortality_rates', & - cohort_in%dbh,cohort_in%pft,cohort_in%n,cohort_in%canopy_layer,cohort_in%indexnumber + write(fates_log(),*) 'dbh problem in mortality_rates', & + cohort_in%dbh,cohort_in%pft,cohort_in%n,cohort_in%canopy_layer endif !mortality_rates = bmort + hmort + cmort diff --git a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 index 918065809b..b771c890da 100755 --- a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 +++ b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 @@ -42,7 +42,7 @@ module EDPhysiologyMod contains ! ============================================================================ - subroutine canopy_derivs( currentSite, currentPatch ) + subroutine canopy_derivs( currentSite, currentPatch, bc_in ) ! ! !DESCRIPTION: ! spawn new cohorts of juveniles of each PFT @@ -52,6 +52,7 @@ subroutine canopy_derivs( currentSite, currentPatch ) ! !ARGUMENTS type(ed_site_type), intent(inout), target :: currentSite type(ed_patch_type) , intent(inout), target :: currentPatch + type(bc_in_type), intent(in) :: bc_in ! ! !LOCAL VARIABLES: type(ed_cohort_type), pointer ::currentCohort @@ -62,7 +63,7 @@ subroutine canopy_derivs( currentSite, currentPatch ) currentCohort => currentPatch%shortest do while(associated(currentCohort)) - call Growth_Derivatives(currentSite, currentCohort) + call Growth_Derivatives(currentSite, currentCohort, bc_in ) currentCohort => currentCohort%taller enddo @@ -702,7 +703,7 @@ subroutine seed_germination( currentSite, currentPatch ) end subroutine seed_germination ! ============================================================================ - subroutine Growth_Derivatives( currentSite, currentCohort) + subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) ! ! !DESCRIPTION: ! Main subroutine controlling growth and allocation derivatives @@ -714,6 +715,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort) ! !ARGUMENTS type(ed_site_type), intent(inout), target :: currentSite type(ed_cohort_type),intent(inout), target :: currentCohort + type(bc_in_type), intent(in) :: bc_in ! ! !LOCAL VARIABLES: real(r8) :: dbldbd !rate of change of dead biomass per unit dbh @@ -760,9 +762,11 @@ subroutine Growth_Derivatives( currentSite, currentCohort) ! NPP if ( DEBUG ) write(fates_log(),*) 'EDphys 716 ',currentCohort%npp_acc - currentCohort%npp_acc_hold = currentCohort%npp_acc * udata%n_sub ! convert from kgC/indiv/day into kgC/indiv/year - currentCohort%gpp_acc_hold = currentCohort%gpp_acc * udata%n_sub ! convert from kgC/indiv/day into kgC/indiv/year - currentCohort%resp_acc_hold = currentCohort%resp_acc * udata%n_sub ! convert from kgC/indiv/day into kgC/indiv/year + ! convert from kgC/indiv/day into kgC/indiv/year + ! TODO: CONVERT DAYS_PER_YEAR TO DBLE (HOLDING FOR B4B COMPARISONS, RGK-01-2017) + currentCohort%npp_acc_hold = currentCohort%npp_acc * bc_in%days_per_year + currentCohort%gpp_acc_hold = currentCohort%gpp_acc * bc_in%days_per_year + currentCohort%resp_acc_hold = currentCohort%resp_acc * bc_in%days_per_year currentSite%flux_in = currentSite%flux_in + currentCohort%npp_acc * currentCohort%n @@ -845,7 +849,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort) currentCohort%storage_flux = 0._r8 currentCohort%carbon_balance = 0._r8 - write(fates_log(),*) 'ED: no leaf area in gd', currentCohort%indexnumber,currentCohort%n,currentCohort%bdead, & + write(fates_log(),*) 'ED: no leaf area in gd',currentCohort%n,currentCohort%bdead, & currentCohort%dbh,currentCohort%balive endif diff --git a/components/clm/src/ED/main/EDMainMod.F90 b/components/clm/src/ED/main/EDMainMod.F90 index e82dc7b512..e18676e261 100755 --- a/components/clm/src/ED/main/EDMainMod.F90 +++ b/components/clm/src/ED/main/EDMainMod.F90 @@ -6,7 +6,7 @@ module EDMainMod use shr_kind_mod , only : r8 => shr_kind_r8 - use clm_varctl , only : iulog + use FatesGlobals , only : fates_log use atm2lndType , only : atm2lnd_type use SoilStateType , only : soilstate_type use TemperatureType , only : temperature_type @@ -53,7 +53,8 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in) type(ed_patch_type), pointer :: currentPatch !----------------------------------------------------------------------- - if ( cp_masterproc==1 ) write(iulog,*) 'modelday',bc_in%model_day + if ( cp_masterproc==1 ) write(fates_log(),'(A,I4,A,I2.2,A,I2.2)') 'FATES Dynamics: ',& + bc_in%current_year,'-',bc_in%current_month,'-',bc_in%current_day !************************************************************************** ! Fire, growth, biogeochemistry. @@ -165,12 +166,12 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) currentPatch%age = currentPatch%age + udata%deltat ! FIX(SPM,032414) valgrind 'Conditional jump or move depends on uninitialised value' if( currentPatch%age < 0._r8 )then - write(iulog,*) 'negative patch age?',currentPatch%age, & + write(fates_log(),*) 'negative patch age?',currentPatch%age, & currentPatch%patchno,currentPatch%area endif ! Find the derivatives of the growth and litter processes. - call canopy_derivs(currentSite, currentPatch) + call canopy_derivs(currentSite, currentPatch, bc_in) ! Update Canopy Biomass Pools currentCohort => currentPatch%shortest @@ -181,23 +182,23 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) currentCohort%balive = currentCohort%balive + currentCohort%dbalivedt * udata%deltat currentCohort%bdead = max(small_no,currentCohort%bdead + currentCohort%dbdeaddt * udata%deltat ) if ( DEBUG ) then - write(iulog,*) 'EDMainMod dbstoredt I ',currentCohort%bstore, & + write(fates_log(),*) 'EDMainMod dbstoredt I ',currentCohort%bstore, & currentCohort%dbstoredt,udata%deltat end if currentCohort%bstore = currentCohort%bstore + currentCohort%dbstoredt * udata%deltat if ( DEBUG ) then - write(iulog,*) 'EDMainMod dbstoredt II ',currentCohort%bstore, & + write(fates_log(),*) 'EDMainMod dbstoredt II ',currentCohort%bstore, & currentCohort%dbstoredt,udata%deltat end if if( (currentCohort%balive+currentCohort%bdead+currentCohort%bstore)*currentCohort%n<0._r8)then - write(iulog,*) 'biomass is negative', currentCohort%n,currentCohort%balive, & + write(fates_log(),*) 'biomass is negative', currentCohort%n,currentCohort%balive, & currentCohort%bdead,currentCohort%bstore endif if(abs((currentCohort%balive+currentCohort%bdead+currentCohort%bstore+udata%deltat*(currentCohort%md+ & currentCohort%seed_prod)-cohort_biomass_store)-currentCohort%npp_acc) > 1e-8_r8)then - write(iulog,*) 'issue with c balance in integration', abs(currentCohort%balive+currentCohort%bdead+ & + write(fates_log(),*) 'issue with c balance in integration', abs(currentCohort%balive+currentCohort%bdead+ & currentCohort%bstore+udata%deltat* & (currentCohort%md+currentCohort%seed_prod)-cohort_biomass_store-currentCohort%npp_acc) endif @@ -234,23 +235,25 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) do c = 1,ncwd if(currentPatch%cwd_ag(c) currentPatch%younger @@ -343,8 +347,13 @@ subroutine ed_update_site( currentSite ) enddo ! FIX(RF,032414). This needs to be monthly, not annual - if((udata%time_period == udata%n_sub-1))then - write(iulog,*) 'calling trim canopy' +! if((udata%time_period == udata%n_sub-1))then + + ! If this is the second to last day of the year, then perform trimming + + if( bc_in%day_of_year == bc_in%days_per_year-1) then + + write(fates_log(),*) 'calling trim canopy' call trim_canopy(currentSite) endif @@ -414,14 +423,14 @@ subroutine ed_total_balance_check (currentSite, call_index ) error = abs(net_flux - change_in_stock) if ( abs(error) > 10e-6 ) then - write(iulog,*) 'total error: call index: ',call_index, & + write(fates_log(),*) 'total error: call index: ',call_index, & 'in: ',currentSite%flux_in, & 'out: ',currentSite%flux_out, & 'net: ',net_flux, & 'dstock: ',change_in_stock, & 'error=net_flux-dstock:', error - write(iulog,*) 'biomass,litter,seeds', biomass_stock,litter_stock,seed_stock - write(iulog,*) 'lat lon',currentSite%lat,currentSite%lon + write(fates_log(),*) 'biomass,litter,seeds', biomass_stock,litter_stock,seed_stock + write(fates_log(),*) 'lat lon',currentSite%lat,currentSite%lon endif currentSite%flux_in = 0.0_r8 diff --git a/components/clm/src/ED/main/EDTypesMod.F90 b/components/clm/src/ED/main/EDTypesMod.F90 index 821ed08b52..ec47eba03b 100755 --- a/components/clm/src/ED/main/EDTypesMod.F90 +++ b/components/clm/src/ED/main/EDTypesMod.F90 @@ -552,12 +552,11 @@ 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? - integer :: modelday ! Number of days since reference end type userdata diff --git a/components/clm/src/ED/main/FatesInterfaceMod.F90 b/components/clm/src/ED/main/FatesInterfaceMod.F90 index a8bf0cc6a2..f24a415ea3 100644 --- a/components/clm/src/ED/main/FatesInterfaceMod.F90 +++ b/components/clm/src/ED/main/FatesInterfaceMod.F90 @@ -60,7 +60,15 @@ module FatesInterfaceMod real(r8) :: model_day ! elapsed days between current date and reference ! uses ESMF functions, so prefered to pass it in as ! argument rather than calculate directly - + integer :: day_of_year ! The integer day of the year + integer :: days_per_year ! The HLM controls time, some HLMs may include a leap + ! day, some actually don't. This is the number of + ! days in the current year + real(r8) :: deltat_day ! fraction of year for each time-step (1/days_per_year) + + + + ! Vegetation Dynamics ! --------------------------------------------------------------------------------- diff --git a/components/clm/src/utils/clmfates_interfaceMod.F90 b/components/clm/src/utils/clmfates_interfaceMod.F90 index 493105fbad..f05ed166d3 100644 --- a/components/clm/src/utils/clmfates_interfaceMod.F90 +++ b/components/clm/src/utils/clmfates_interfaceMod.F90 @@ -473,8 +473,11 @@ subroutine dynamics_driv(this, nc, bounds_clump, & integer :: current_day integer :: current_tod integer :: current_date + integer :: jan01_curr_year integer :: reference_date + integer :: days_per_year real(r8) :: model_day + real(r8) :: day_of_year !----------------------------------------------------------------------- @@ -498,14 +501,14 @@ subroutine dynamics_driv(this, nc, bounds_clump, & 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) - udata%modelday = dayDiff dayDiffInt = floor(dayDiff) udata%time_period = mod( dayDiffInt , udata%n_sub ) - + ! --------------------------------------------------------------------------------- ! Prepare input boundary conditions for FATES dynamics ! Note that timing information is the same across all sites, this may @@ -513,14 +516,19 @@ subroutine dynamics_driv(this, nc, bounds_clump, & ! one day. The cost of holding site level boundary conditions is minimal ! and it keeps all the boundaries in one location ! --------------------------------------------------------------------------------- + + days_per_year = get_days_per_year() call get_curr_date(current_year,current_month,current_day,current_tod) current_date = current_year*10000 + current_month*100 + current_day + jan01_curr_year = current_year*10000 + 100 + 1 call get_ref_date(yr, mon, day, sec) reference_date = yr*10000 + mon*100 + day call timemgr_datediff(reference_date, sec, current_date, current_tod, model_day) - if ( masterproc ) write(iulog,*) 'modelday',model_day + + ! Calculate the day of the year + call timemgr_datediff(jan01_curr_year,0,ncdate,sec,day_of_year) do s=1,this%fates(nc)%nsites c = this%f2hmap(nc)%fcolumn(s) @@ -531,6 +539,9 @@ subroutine dynamics_driv(this, nc, bounds_clump, & this%fates(nc)%bc_in(s)%current_tod = current_date this%fates(nc)%bc_in(s)%reference_date = reference_date this%fates(nc)%bc_in(s)%model_day = model_day + this%fates(nc)%bc_in(s)%days_per_year = days_per_year + this%fates(nc)%bc_in(s)%day_of_year = floor(day_of_year) + this%fates(nc)%bc_in(s)%deltat_day = 1.0_r8/dble(day_of_year) this%fates(nc)%bc_in(s)%h2osoi_vol_si = & waterstate_inst%h2osoi_vol_col(c,1) @@ -560,7 +571,8 @@ subroutine dynamics_driv(this, nc, bounds_clump, & call ed_ecosystem_dynamics(this%fates(nc)%sites(s), & this%fates(nc)%bc_in(s)) - call ed_update_site(this%fates(nc)%sites(s)) + call ed_update_site(this%fates(nc)%sites(s), & + this%fates(nc)%bc_in(s)) enddo @@ -933,7 +945,8 @@ subroutine restart( this, bounds_proc, ncid, flag, waterstate_inst, canopystate_ ! I think ed_update_site and update_hlmfates_dyn are doing some similar ! update type stuff, should consolidate (rgk 11-2016) do s = 1,this%fates(nc)%nsites - call ed_update_site( this%fates(nc)%sites(s) ) + call ed_update_site( this%fates(nc)%sites(s), & + this%fates(nc)%bc_in(s) ) end do ! ------------------------------------------------------------------------ @@ -995,7 +1008,8 @@ subroutine init_coldstart(this, waterstate_inst, canopystate_inst) call init_patches(this%fates(nc)%nsites, this%fates(nc)%sites) do s = 1,this%fates(nc)%nsites - call ed_update_site(this%fates(nc)%sites(s)) + call ed_update_site(this%fates(nc)%sites(s), & + this%fates(nc)%bc_in(s)) end do From 69b4c24ebdb3a1b5e9b30f4aebf5e481a94d6957 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 12 Jan 2017 16:35:04 -0800 Subject: [PATCH 09/11] Finished first pass of time-control refactor. --- .../src/ED/biogeochem/EDCohortDynamicsMod.F90 | 13 ++-- .../src/ED/biogeochem/EDPatchDynamicsMod.F90 | 14 ++--- .../clm/src/ED/biogeochem/EDPhysiologyMod.F90 | 44 ++++++------- components/clm/src/ED/main/EDInitMod.F90 | 2 +- components/clm/src/ED/main/EDMainMod.F90 | 48 ++++++++------- components/clm/src/ED/main/EDTypesMod.F90 | 16 +++-- components/clm/src/ED/main/FatesGlobals.F90 | 61 ++++++++++++++++++- .../clm/src/ED/main/FatesInterfaceMod.F90 | 36 ++--------- .../clm/src/utils/clmfates_interfaceMod.F90 | 58 ++++-------------- 9 files changed, 145 insertions(+), 147 deletions(-) diff --git a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 index 70b00c3307..7b40aa4203 100755 --- a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 +++ b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 @@ -6,6 +6,7 @@ module EDCohortDynamicsMod ! !USES: use abortutils , only : endrun use FatesGlobals , only : fates_log + use FatesGlobals , only : freq_day use FatesConstantsMod , only : r8 => fates_r8 use FatesConstantsMod , only : fates_unset_int use shr_log_mod , only : errMsg => shr_log_errMsg @@ -14,7 +15,7 @@ module EDCohortDynamicsMod use EDGrowthFunctionsMod , only : c_area, tree_lai use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : fusetol, cp_nclmax - use EDtypesMod , only : ncwd, maxcohortsperpatch, udata + use EDtypesMod , only : ncwd, maxcohortsperpatch use EDtypesMod , only : sclass_ed,nlevsclass_ed,AREA use EDtypesMod , only : min_npm2, min_nppatch, min_n_safemath ! @@ -223,7 +224,7 @@ subroutine allocate_live_biomass(cc_p,mode) if(mode==1)then ! it will not be able to put out as many leaves as it had previous timestep currentcohort%npp_leaf = currentcohort%npp_leaf + & - max(0.0_r8,currentcohort%balive*leaf_frac - currentcohort%bl)/udata%deltat + max(0.0_r8,currentcohort%balive*leaf_frac - currentcohort%bl)/freq_day end if currentcohort%bl = currentcohort%balive*leaf_frac @@ -234,10 +235,10 @@ subroutine allocate_live_biomass(cc_p,mode) currentcohort%npp_froot = currentcohort%npp_froot + & max(0._r8,pftcon%froot_leaf(ft)*(currentcohort%balive+currentcohort%laimemory)*leaf_frac - currentcohort%br) / & - udata%deltat + freq_day currentcohort%npp_bsw = max(0._r8,EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + & - currentcohort%laimemory)*leaf_frac - currentcohort%bsw)/udata%deltat + currentcohort%laimemory)*leaf_frac - currentcohort%bsw)/freq_day currentcohort%npp_bdead = currentCohort%dbdeaddt @@ -273,10 +274,10 @@ subroutine allocate_live_biomass(cc_p,mode) currentcohort%npp_froot = currentcohort%npp_froot + & max(0.0_r8,pftcon%froot_leaf(ft)*(ideal_balive + & - currentcohort%laimemory)*leaf_frac*ratio_balive-currentcohort%br)/udata%deltat + currentcohort%laimemory)*leaf_frac*ratio_balive-currentcohort%br)/freq_day currentcohort%npp_bsw = max(0.0_r8,EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(ideal_balive + & - currentcohort%laimemory)*leaf_frac*ratio_balive - currentcohort%bsw)/udata%deltat + currentcohort%laimemory)*leaf_frac*ratio_balive - currentcohort%bsw)/freq_day currentcohort%npp_bdead = currentCohort%dbdeaddt diff --git a/components/clm/src/ED/biogeochem/EDPatchDynamicsMod.F90 b/components/clm/src/ED/biogeochem/EDPatchDynamicsMod.F90 index 5fae1a783f..6d7c84ebb1 100755 --- a/components/clm/src/ED/biogeochem/EDPatchDynamicsMod.F90 +++ b/components/clm/src/ED/biogeochem/EDPatchDynamicsMod.F90 @@ -7,10 +7,11 @@ module EDPatchDynamicsMod use shr_kind_mod , only : r8 => shr_kind_r8; use shr_infnan_mod , only : nan => shr_infnan_nan, assignment(=) use clm_varctl , only : iulog + use FatesGlobals , only : freq_day use pftconMod , only : pftcon use EDCohortDynamicsMod , only : fuse_cohorts, sort_cohorts, insert_cohort use EDtypesMod , only : ncwd, n_dbh_bins, ntol, numpft_ed, area, dbhmax, maxPatchesPerCol - use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type, udata + use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use EDTypesMod , only : min_patch_area, cp_numlevgrnd, cp_numSWb ! implicit none @@ -45,7 +46,6 @@ subroutine disturbance_rates( site_in) ! ! !USES: use EDGrowthFunctionsMod , only : c_area, mortality_rates - use EDTypesMod , only : udata ! ! !ARGUMENTS: type(ed_site_type) , intent(inout), target :: site_in @@ -85,7 +85,7 @@ subroutine disturbance_rates( site_in) if(currentCohort%canopy_layer == 1)then currentPatch%disturbance_rates(1) = currentPatch%disturbance_rates(1) + & - min(1.0_r8,currentCohort%dmort)*udata%deltat*currentCohort%c_area/currentPatch%area + min(1.0_r8,currentCohort%dmort)*freq_day*currentCohort%c_area/currentPatch%area endif @@ -271,7 +271,7 @@ subroutine spawn_patches( currentSite ) ! because this is the part of the original patch where no trees have actually fallen ! The diagnostic cmort,bmort and hmort rates have already been saved - currentCohort%n = currentCohort%n * (1.0_r8 - min(1.0_r8,currentCohort%dmort * udata%deltat)) + currentCohort%n = currentCohort%n * (1.0_r8 - min(1.0_r8,currentCohort%dmort * freq_day)) nc%n = 0.0_r8 ! kill all of the trees who caused the disturbance. nc%cmort = nan ! The mortality diagnostics are set to nan because the cohort should dissappear @@ -298,7 +298,7 @@ subroutine spawn_patches( currentSite ) ! so with the number density must come the effective mortality rates. nc%fmort = 0.0_r8 ! Should had also been zero in the donor - nc%imort = ED_val_understorey_death/udata%deltat ! This was zero in the donor + nc%imort = ED_val_understorey_death/freq_day ! This was zero in the donor nc%cmort = currentCohort%cmort nc%hmort = currentCohort%hmort nc%bmort = currentCohort%bmort @@ -336,7 +336,7 @@ subroutine spawn_patches( currentSite ) ! loss of individual from fire in new patch. nc%n = nc%n * (1.0_r8 - currentCohort%fire_mort) - nc%fmort = currentCohort%fire_mort/udata%deltat + nc%fmort = currentCohort%fire_mort/freq_day nc%imort = 0.0_r8 nc%cmort = currentCohort%cmort nc%hmort = currentCohort%hmort @@ -716,7 +716,7 @@ subroutine mortality_litter_fluxes(cp_target, new_patch_target, patch_site_aread !currentCohort%dmort = mortality_rates(currentCohort) !the disturbance calculations are done with the previous n, c_area and d_mort. So it's probably & !not right to recalcualte dmort here. - canopy_dead = currentCohort%n * min(1.0_r8,currentCohort%dmort * udata%deltat) + canopy_dead = currentCohort%n * min(1.0_r8,currentCohort%dmort * freq_day) currentPatch%canopy_mortality_woody_litter = currentPatch%canopy_mortality_woody_litter + & canopy_dead*(currentCohort%bdead+currentCohort%bsw) diff --git a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 index b771c890da..a496767038 100755 --- a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 +++ b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 @@ -7,6 +7,10 @@ module EDPhysiologyMod ! ============================================================================ use FatesGlobals, only : fates_log + use FatesGlobals, only : days_per_year + use FatesGlobals, only : model_day + use FatesGlobals, only : freq_day + use FatesGlobals, only : day_of_year use FatesConstantsMod, only : r8 => fates_r8 use pftconMod , only : pftcon use EDEcophysContype , only : EDecophyscon @@ -240,7 +244,6 @@ subroutine phenology( currentSite, bc_in ) ! ! !USES: use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm - use EDTypesMod, only : udata ! ! !ARGUMENTS: @@ -289,7 +292,7 @@ subroutine phenology( currentSite, bc_in ) ncolddayslim = 5 cold_t = 7.5_r8 ! ed_ph_coldtemp - t = udata%time_period + t = day_of_year temp_in_C = bc_in%t_veg24_si - tfrz !-----------------Cold Phenology--------------------! @@ -339,7 +342,7 @@ subroutine phenology( currentSite, bc_in ) endif - timesinceleafoff = bc_in%model_day - currentSite%leafoffdate + timesinceleafoff = model_day - currentSite%leafoffdate !LEAF ON: COLD DECIDUOUS. Needs to !1) have exceeded the growing degree day threshold !2) The leaves should not be on already @@ -355,7 +358,7 @@ subroutine phenology( currentSite, bc_in ) endif !status endif !GDD - timesinceleafon = bc_in%model_day - currentSite%leafondate + timesinceleafon = model_day - currentSite%leafondate !LEAF OFF: COLD THRESHOLD @@ -369,7 +372,7 @@ subroutine phenology( currentSite, bc_in ) if (timesinceleafon > mindayson)then if (currentSite%status == 2)then currentSite%status = 1 !alter status of site to 'leaves on' - currentSite%leafoffdate = bc_in%model_day !record leaf off date + currentSite%leafoffdate = model_day !record leaf off date if ( DEBUG ) write(fates_log(),*) 'leaves off' endif endif @@ -379,7 +382,7 @@ subroutine phenology( currentSite, bc_in ) if(timesinceleafoff > 400)then !remove leaves after a whole year when there is no 'off' period. if(currentSite%status == 2)then currentSite%status = 1 !alter status of site to 'leaves on' - currentSite%leafoffdate = bc_in%model_day !record leaf off date + currentSite%leafoffdate = model_day !record leaf off date if ( DEBUG ) write(fates_log(),*) 'leaves off' endif endif @@ -710,7 +713,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) ! ! !USES: use EDGrowthFunctionsMod , only : Bleaf, dDbhdBd, dhdbd, hite, mortality_rates,dDbhdBl - use EDTypesMod , only : udata + ! ! !ARGUMENTS type(ed_site_type), intent(inout), target :: currentSite @@ -764,9 +767,9 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) ! convert from kgC/indiv/day into kgC/indiv/year ! TODO: CONVERT DAYS_PER_YEAR TO DBLE (HOLDING FOR B4B COMPARISONS, RGK-01-2017) - currentCohort%npp_acc_hold = currentCohort%npp_acc * bc_in%days_per_year - currentCohort%gpp_acc_hold = currentCohort%gpp_acc * bc_in%days_per_year - currentCohort%resp_acc_hold = currentCohort%resp_acc * bc_in%days_per_year + currentCohort%npp_acc_hold = currentCohort%npp_acc * days_per_year + currentCohort%gpp_acc_hold = currentCohort%gpp_acc * days_per_year + currentCohort%resp_acc_hold = currentCohort%resp_acc * days_per_year currentSite%flux_in = currentSite%flux_in + currentCohort%npp_acc * currentCohort%n @@ -933,7 +936,7 @@ subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) ! prevent negative leaf pool (but not negative store pool). This is also a numerical error prevention, ! but it shouldn't happen actually... - if (-1.0_r8*currentCohort%dbalivedt * udata%deltat > currentCohort%balive*0.99)then + if (-1.0_r8*currentCohort%dbalivedt * freq_day > currentCohort%balive*0.99)then write(fates_log(),*) 'using non-neg leaf mass cap',currentCohort%balive , currentCohort%dbalivedt,currentCohort%dbstoredt, & currentCohort%carbon_balance currentCohort%dbstoredt = currentCohort%dbstoredt + currentCohort%dbalivedt @@ -963,7 +966,6 @@ subroutine recruitment( t, currentSite, currentPatch ) ! ! !USES: use EDGrowthFunctionsMod, only : bdead,dbh, Bleaf - use EDTypesMod, only : udata ! ! !ARGUMENTS integer, intent(in) :: t @@ -990,7 +992,7 @@ subroutine recruitment( t, currentSite, currentPatch ) + EDecophyscon%sapwood_ratio(ft)*temp_cohort%hite) temp_cohort%bstore = EDecophyscon%cushion(ft)*(temp_cohort%balive/ (1.0_r8 + pftcon%froot_leaf(ft) & + EDecophyscon%sapwood_ratio(ft)*temp_cohort%hite)) - temp_cohort%n = currentPatch%area * currentPatch%seed_germination(ft)*udata%deltat & + temp_cohort%n = currentPatch%area * currentPatch%seed_germination(ft)*freq_day & / (temp_cohort%bdead+temp_cohort%balive+temp_cohort%bstore) if (t == 1)then @@ -1037,7 +1039,7 @@ subroutine CWD_Input( currentPatch) ! !USES: use SFParamsMod , only : SF_val_CWD_frac use EDParamsMod , only : ED_val_ag_biomass - use EDTypesMod , only : udata + ! ! !ARGUMENTS type(ed_patch_type),intent(inout), target :: currentPatch @@ -1067,7 +1069,7 @@ subroutine CWD_Input( currentPatch) currentPatch%root_litter_in(pft) = currentPatch%root_litter_in(pft) + & currentCohort%root_md * currentCohort%n/currentPatch%area !turnover currentPatch%leaf_litter_in(pft) = currentPatch%leaf_litter_in(pft) + & - currentCohort%leaf_litter * currentCohort%n/currentPatch%area/udata%deltat + currentCohort%leaf_litter * currentCohort%n/currentPatch%area/freq_day !daily leaf loss needs to be scaled up to the annual scale here. @@ -1086,7 +1088,7 @@ subroutine CWD_Input( currentPatch) dead_n = -1.0_r8 * currentCohort%dndt / currentPatch%area currentPatch%leaf_litter_in(pft) = currentPatch%leaf_litter_in(pft) + & - (currentCohort%bl+currentCohort%leaf_litter/udata%deltat)* dead_n + (currentCohort%bl+currentCohort%leaf_litter/freq_day)* dead_n currentPatch%root_litter_in(pft) = currentPatch%root_litter_in(pft) + & (currentCohort%br+currentCohort%bstore) * dead_n @@ -1191,7 +1193,7 @@ subroutine cwd_out( currentSite, currentPatch, bc_in ) ! ! !USES: use SFParamsMod, only : SF_val_max_decomp - use EDTypesMod , only : udata + ! ! !ARGUMENTS type(ed_site_type), intent(inout), target :: currentSite @@ -1239,13 +1241,13 @@ subroutine cwd_out( currentSite, currentPatch, bc_in ) !add up carbon going into fragmenting pools currentSite%flux_out = currentSite%flux_out + sum(currentPatch%leaf_litter_out) * & - currentPatch%area *udata%deltat!kgC/site/day + currentPatch%area *freq_day!kgC/site/day currentSite%flux_out = currentSite%flux_out + sum(currentPatch%root_litter_out) * & - currentPatch%area *udata%deltat!kgC/site/day + currentPatch%area *freq_day!kgC/site/day currentSite%flux_out = currentSite%flux_out + sum(currentPatch%cwd_ag_out) * & - currentPatch%area *udata%deltat!kgC/site/day + currentPatch%area *freq_day!kgC/site/day currentSite%flux_out = currentSite%flux_out + sum(currentPatch%cwd_bg_out) * & - currentPatch%area *udata%deltat!kgC/site/day + currentPatch%area *freq_day!kgC/site/day end subroutine cwd_out diff --git a/components/clm/src/ED/main/EDInitMod.F90 b/components/clm/src/ED/main/EDInitMod.F90 index 55a0fd8a48..76bc5ed9b2 100755 --- a/components/clm/src/ED/main/EDInitMod.F90 +++ b/components/clm/src/ED/main/EDInitMod.F90 @@ -16,7 +16,7 @@ module EDInitMod use EDCohortDynamicsMod , only : create_cohort, fuse_cohorts, sort_cohorts use EDPatchDynamicsMod , only : create_patch use EDTypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type, area - use EDTypesMod , only : cohorts_per_col, ncwd, numpft_ed, udata + use EDTypesMod , only : cohorts_per_col, ncwd, numpft_ed implicit none private diff --git a/components/clm/src/ED/main/EDMainMod.F90 b/components/clm/src/ED/main/EDMainMod.F90 index e18676e261..6882222f1e 100755 --- a/components/clm/src/ED/main/EDMainMod.F90 +++ b/components/clm/src/ED/main/EDMainMod.F90 @@ -7,6 +7,12 @@ module EDMainMod use shr_kind_mod , only : r8 => shr_kind_r8 use FatesGlobals , only : fates_log + use FatesGlobals , only : freq_day + use FatesGlobals , only : day_of_year + use FatesGlobals , only : days_per_year + use FatesGlobals , only : current_year + use FatesGlobals , only : current_month + use FatesGlobals , only : current_day use atm2lndType , only : atm2lnd_type use SoilStateType , only : soilstate_type use TemperatureType , only : temperature_type @@ -14,10 +20,11 @@ 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, numpft_ed, udata + use EDtypesMod , only : ncwd, numpft_ed use EDtypesMod , only : ed_site_type, ed_patch_type, ed_cohort_type use FatesInterfaceMod , only : bc_in_type use EDTypesMod , only : cp_masterproc + implicit none private @@ -54,7 +61,7 @@ subroutine ed_ecosystem_dynamics(currentSite, bc_in) !----------------------------------------------------------------------- if ( cp_masterproc==1 ) write(fates_log(),'(A,I4,A,I2.2,A,I2.2)') 'FATES Dynamics: ',& - bc_in%current_year,'-',bc_in%current_month,'-',bc_in%current_day + current_year,'-',current_month,'-',current_day !************************************************************************** ! Fire, growth, biogeochemistry. @@ -163,7 +170,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) do while(associated(currentPatch)) - currentPatch%age = currentPatch%age + udata%deltat + currentPatch%age = currentPatch%age + freq_day ! FIX(SPM,032414) valgrind 'Conditional jump or move depends on uninitialised value' if( currentPatch%age < 0._r8 )then write(fates_log(),*) 'negative patch age?',currentPatch%age, & @@ -178,17 +185,17 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) do while(associated(currentCohort)) cohort_biomass_store = (currentCohort%balive+currentCohort%bdead+currentCohort%bstore) - currentCohort%dbh = max(small_no,currentCohort%dbh + currentCohort%ddbhdt * udata%deltat ) - currentCohort%balive = currentCohort%balive + currentCohort%dbalivedt * udata%deltat - currentCohort%bdead = max(small_no,currentCohort%bdead + currentCohort%dbdeaddt * udata%deltat ) + currentCohort%dbh = max(small_no,currentCohort%dbh + currentCohort%ddbhdt * freq_day ) + currentCohort%balive = currentCohort%balive + currentCohort%dbalivedt * freq_day + currentCohort%bdead = max(small_no,currentCohort%bdead + currentCohort%dbdeaddt * freq_day ) if ( DEBUG ) then write(fates_log(),*) 'EDMainMod dbstoredt I ',currentCohort%bstore, & - currentCohort%dbstoredt,udata%deltat + currentCohort%dbstoredt,freq_day end if - currentCohort%bstore = currentCohort%bstore + currentCohort%dbstoredt * udata%deltat + currentCohort%bstore = currentCohort%bstore + currentCohort%dbstoredt * freq_day if ( DEBUG ) then write(fates_log(),*) 'EDMainMod dbstoredt II ',currentCohort%bstore, & - currentCohort%dbstoredt,udata%deltat + currentCohort%dbstoredt,freq_day end if if( (currentCohort%balive+currentCohort%bdead+currentCohort%bstore)*currentCohort%n<0._r8)then @@ -196,10 +203,10 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) currentCohort%bdead,currentCohort%bstore endif - if(abs((currentCohort%balive+currentCohort%bdead+currentCohort%bstore+udata%deltat*(currentCohort%md+ & + if(abs((currentCohort%balive+currentCohort%bdead+currentCohort%bstore+freq_day*(currentCohort%md+ & currentCohort%seed_prod)-cohort_biomass_store)-currentCohort%npp_acc) > 1e-8_r8)then write(fates_log(),*) 'issue with c balance in integration', abs(currentCohort%balive+currentCohort%bdead+ & - currentCohort%bstore+udata%deltat* & + currentCohort%bstore+freq_day* & (currentCohort%md+currentCohort%seed_prod)-cohort_biomass_store-currentCohort%npp_acc) endif @@ -224,13 +231,13 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) ! first update the litter variables that are tracked at the patch level do c = 1,ncwd - currentPatch%cwd_ag(c) = currentPatch%cwd_ag(c) + currentPatch%dcwd_ag_dt(c)* udata%deltat - currentPatch%cwd_bg(c) = currentPatch%cwd_bg(c) + currentPatch%dcwd_bg_dt(c)* udata%deltat + currentPatch%cwd_ag(c) = currentPatch%cwd_ag(c) + currentPatch%dcwd_ag_dt(c)* freq_day + currentPatch%cwd_bg(c) = currentPatch%cwd_bg(c) + currentPatch%dcwd_bg_dt(c)* freq_day enddo do ft = 1,numpft_ed - currentPatch%leaf_litter(ft) = currentPatch%leaf_litter(ft) + currentPatch%dleaf_litter_dt(ft)* udata%deltat - currentPatch%root_litter(ft) = currentPatch%root_litter(ft) + currentPatch%droot_litter_dt(ft)* udata%deltat + currentPatch%leaf_litter(ft) = currentPatch%leaf_litter(ft) + currentPatch%dleaf_litter_dt(ft)* freq_day + currentPatch%root_litter(ft) = currentPatch%root_litter(ft) + currentPatch%droot_litter_dt(ft)* freq_day enddo do c = 1,ncwd @@ -254,7 +261,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) endif if(currentPatch%root_litter(ft) currentPatch%shortest do while(associated(currentCohort)) - currentCohort%n = max(small_no,currentCohort%n + currentCohort%dndt * udata%deltat ) + currentCohort%n = max(small_no,currentCohort%n + currentCohort%dndt * freq_day ) currentCohort => currentCohort%taller enddo @@ -275,7 +282,7 @@ subroutine ed_integrate_state_variables(currentSite, bc_in ) ! at the site level, update the seed bank mass do ft = 1,numpft_ed - currentSite%seed_bank(ft) = currentSite%seed_bank(ft) + currentSite%dseed_dt(ft)*udata%deltat + currentSite%seed_bank(ft) = currentSite%seed_bank(ft) + currentSite%dseed_dt(ft)*freq_day enddo ! Check for negative values. Write out warning to show carbon balance. @@ -347,11 +354,8 @@ subroutine ed_update_site( currentSite, bc_in ) enddo ! FIX(RF,032414). This needs to be monthly, not annual -! if((udata%time_period == udata%n_sub-1))then - ! If this is the second to last day of the year, then perform trimming - - if( bc_in%day_of_year == bc_in%days_per_year-1) then + if( day_of_year == days_per_year-1) then write(fates_log(),*) 'calling trim canopy' call trim_canopy(currentSite) diff --git a/components/clm/src/ED/main/EDTypesMod.F90 b/components/clm/src/ED/main/EDTypesMod.F90 index ec47eba03b..a6a34ba71f 100755 --- a/components/clm/src/ED/main/EDTypesMod.F90 +++ b/components/clm/src/ED/main/EDTypesMod.F90 @@ -551,16 +551,14 @@ module EDTypesMod !** Userdata type structure ** !************************************ - type userdata +! type userdata ! 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 ! THIS WAS NOT THREADSAFE +! 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 ! THIS WAS NOT THREADSAFE !-------------------------------------------------------------------------------------! public :: ed_hist_scpfmaps diff --git a/components/clm/src/ED/main/FatesGlobals.F90 b/components/clm/src/ED/main/FatesGlobals.F90 index 9ae06e207c..0b4e11e7f3 100644 --- a/components/clm/src/ED/main/FatesGlobals.F90 +++ b/components/clm/src/ED/main/FatesGlobals.F90 @@ -4,14 +4,36 @@ module FatesGlobals ! global data that needs to be dealt with, but doesn't have an ! immediately obvious home. + use FatesConstantsMod , only : r8 => fates_r8 + implicit none - integer, private :: fates_log_ - logical, private :: fates_global_verbose_ + public :: FatesGlobalsInit public :: fates_log public :: fates_global_verbose + public :: SetFatesTime + + ! ------------------------------------------------------------------------------------- + ! Timing Variables + ! It is assumed that all of the sites on a given machine will be synchronous. + ! It is also assumed that the HLM will control time. + ! ------------------------------------------------------------------------------------- + integer, protected :: current_year ! Current year + integer, protected :: current_month ! month of year + integer, protected :: current_day ! day of month + integer, protected :: current_tod ! time of day (seconds past 0Z) + integer, protected :: current_date ! time of day (seconds past 0Z) + integer, protected :: reference_date ! YYYYMMDD + real(r8), protected :: model_day ! elapsed days between current date and reference + integer, protected :: day_of_year ! The integer day of the year + integer, protected :: days_per_year ! The HLM controls time, some HLMs may include a leap + real(r8), protected :: freq_day ! fraction of year for daily time-step (1/days_per_year) + ! this is a frequency + + integer, private :: fates_log_ + logical, private :: fates_global_verbose_ contains @@ -35,4 +57,39 @@ logical function fates_global_verbose() fates_global_verbose = fates_global_verbose_ end function fates_global_verbose + ! ===================================================================================== + + subroutine SetFatesTime(current_year_in, current_month_in, & + current_day_in, current_tod_in, & + current_date_in, reference_date_in, & + model_day_in, day_of_year_in, & + days_per_year_in, freq_day_in) + + ! This subroutine should be called directly from the HLM + + integer, intent(in) :: current_year_in + integer, intent(in) :: current_month_in + integer, intent(in) :: current_day_in + integer, intent(in) :: current_tod_in + integer, intent(in) :: current_date_in + integer, intent(in) :: reference_date_in + real(r8), intent(in) :: model_day_in + integer, intent(in) :: day_of_year_in + integer, intent(in) :: days_per_year_in + real(r8), intent(in) :: freq_day_in + + current_year = current_year_in + current_month = current_month_in + current_day = current_day_in + current_tod = current_tod_in + current_date = current_date_in + reference_date = reference_date_in + model_day = model_day_in + day_of_year = day_of_year_in + days_per_year = days_per_year_in + freq_day = freq_day_in + + end subroutine SetFatesTime + + end module FatesGlobals diff --git a/components/clm/src/ED/main/FatesInterfaceMod.F90 b/components/clm/src/ED/main/FatesInterfaceMod.F90 index f24a415ea3..0d44975be5 100644 --- a/components/clm/src/ED/main/FatesInterfaceMod.F90 +++ b/components/clm/src/ED/main/FatesInterfaceMod.F90 @@ -9,12 +9,6 @@ module FatesInterfaceMod ! which is allocated by thread ! ------------------------------------------------------------------------------------ - ! ------------------------------------------------------------------------------------ - ! Used CLM Modules - ! INTERF-TODO: NO CLM MODULES SHOULD BE ACCESSIBLE BY THE FATES - ! PUBLIC API!!!! - ! ------------------------------------------------------------------------------------ - use EDtypesMod , only : ed_site_type use EDtypesMod , only : maxPatchesPerCol use EDtypesMod , only : cp_nclmax @@ -45,30 +39,14 @@ module FatesInterfaceMod ! _rb means radiation band ! ------------------------------------------------------------------------------------ + + + type, public :: bc_in_type ! The actual number of FATES' ED patches integer :: npatches - ! Timing Variables - integer :: current_year ! Current year - integer :: current_month ! month of year - integer :: current_day ! day of month - integer :: current_tod ! time of day (seconds past 0Z) - integer :: current_date ! time of day (seconds past 0Z) - integer :: reference_date ! YYYYMMDD - real(r8) :: model_day ! elapsed days between current date and reference - ! uses ESMF functions, so prefered to pass it in as - ! argument rather than calculate directly - integer :: day_of_year ! The integer day of the year - integer :: days_per_year ! The HLM controls time, some HLMs may include a leap - ! day, some actually don't. This is the number of - ! days in the current year - real(r8) :: deltat_day ! fraction of year for each time-step (1/days_per_year) - - - - ! Vegetation Dynamics ! --------------------------------------------------------------------------------- @@ -504,13 +482,7 @@ subroutine zero_bcs(this,s) integer, intent(in) :: s ! Input boundaries - this%bc_in(s)%current_year = 0 - this%bc_in(s)%current_month = 0 - this%bc_in(s)%current_day = 0 - this%bc_in(s)%current_tod = 0 - this%bc_in(s)%current_date = 0 - this%bc_in(s)%reference_date = 0 - this%bc_in(s)%model_day = 0.0_r8 + this%bc_in(s)%t_veg24_si = 0.0_r8 this%bc_in(s)%t_veg24_pa(:) = 0.0_r8 this%bc_in(s)%h2osoi_vol_si = 0.0_r8 diff --git a/components/clm/src/utils/clmfates_interfaceMod.F90 b/components/clm/src/utils/clmfates_interfaceMod.F90 index f05ed166d3..150b85a159 100644 --- a/components/clm/src/utils/clmfates_interfaceMod.F90 +++ b/components/clm/src/utils/clmfates_interfaceMod.F90 @@ -82,11 +82,12 @@ module CLMFatesInterfaceMod allocate_bcin, & allocate_bcout + use FatesGlobals , only : SetFatesTime + use FatesHistoryInterfaceMod, only : fates_history_interface_type use FatesRestartInterfaceMod, only : fates_restart_interface_type use ChecksBalancesMod , only : SummarizeNetFluxes, FATES_BGC_Carbon_BalanceCheck - use EDTypesMod , only : udata use EDTypesMod , only : ed_patch_type use EDtypesMod , only : cp_numlevgrnd use EDMainMod , only : ed_ecosystem_dynamics @@ -456,8 +457,6 @@ subroutine dynamics_driv(this, nc, bounds_clump, & type(soilbiogeochem_carbonflux_type), intent(inout) :: soilbiogeochem_carbonflux_inst ! !LOCAL VARIABLES: - real(r8) :: dayDiff ! day of run - integer :: dayDiffInt ! integer of day of run integer :: s ! site index integer :: c ! column index (HLM) integer :: ifp ! patch index @@ -466,8 +465,6 @@ subroutine dynamics_driv(this, nc, bounds_clump, & 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) integer :: current_year integer :: current_month integer :: current_day @@ -480,35 +477,6 @@ subroutine dynamics_driv(this, nc, bounds_clump, & real(r8) :: day_of_year !----------------------------------------------------------------------- - - ! --------------------------------------------------------------------------------- - ! INTERF-TODO: REMOVE ED_DRIVER ARGUMENTS OF CLM STUCTURED TYPES AND - ! REPLACE THEM WITH FATES_BC TYPES WITH ITS OWN MAPPING SCHEME - ! ALSO, NOTE THAT THE ED_DYNAMICS IS A MODULE OF FATES NOW - ! ie: - ! fates(nc)%fatesbc%leaf_temp <=> canopystate_inst% - ! - ! call this%fates(nc)%ed_driver(this%fates(nc)%site, & - ! this%fates(nc)%fatesbc) - ! --------------------------------------------------------------------------------- - - ! timing statements. - udata%n_sub = get_days_per_year() - udata%deltat = 1.0_r8/dble(udata%n_sub) !for working out age of patches in years - if(udata%time_period == 0)then - udata%time_period = udata%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 , udata%n_sub ) - ! --------------------------------------------------------------------------------- ! Prepare input boundary conditions for FATES dynamics ! Note that timing information is the same across all sites, this may @@ -527,21 +495,17 @@ subroutine dynamics_driv(this, nc, bounds_clump, & call timemgr_datediff(reference_date, sec, current_date, current_tod, model_day) - ! Calculate the day of the year - call timemgr_datediff(jan01_curr_year,0,ncdate,sec,day_of_year) + call timemgr_datediff(jan01_curr_year,0,current_date,sec,day_of_year) + + call SetFatesTime(current_year, current_month, & + current_day, current_tod, & + current_date, reference_date, & + model_day, floor(day_of_year), & + days_per_year, 1.0_r8/dble(days_per_year)) + do s=1,this%fates(nc)%nsites c = this%f2hmap(nc)%fcolumn(s) - this%fates(nc)%bc_in(s)%current_year = current_year - this%fates(nc)%bc_in(s)%current_month = current_month - this%fates(nc)%bc_in(s)%current_day = current_day - this%fates(nc)%bc_in(s)%current_tod = current_tod - this%fates(nc)%bc_in(s)%current_tod = current_date - this%fates(nc)%bc_in(s)%reference_date = reference_date - this%fates(nc)%bc_in(s)%model_day = model_day - this%fates(nc)%bc_in(s)%days_per_year = days_per_year - this%fates(nc)%bc_in(s)%day_of_year = floor(day_of_year) - this%fates(nc)%bc_in(s)%deltat_day = 1.0_r8/dble(day_of_year) this%fates(nc)%bc_in(s)%h2osoi_vol_si = & waterstate_inst%h2osoi_vol_col(c,1) @@ -595,7 +559,7 @@ subroutine dynamics_driv(this, nc, bounds_clump, & if (masterproc) then write(iulog, *) 'clm: leaving ED model', bounds_clump%begg, & - bounds_clump%endg, dayDiffInt + bounds_clump%endg end if From db109adf28c773f0c4416c5185ab563849ce15cb Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 12 Jan 2017 17:30:34 -0800 Subject: [PATCH 10/11] Removed the use of the HLMs grid area global ldomain(g)%area from SF. --- components/clm/src/ED/fire/SFMainMod.F90 | 13 +++---------- 1 file changed, 3 insertions(+), 10 deletions(-) diff --git a/components/clm/src/ED/fire/SFMainMod.F90 b/components/clm/src/ED/fire/SFMainMod.F90 index e5557d1ff2..b6ff07c79a 100755 --- a/components/clm/src/ED/fire/SFMainMod.F90 +++ b/components/clm/src/ED/fire/SFMainMod.F90 @@ -642,12 +642,9 @@ subroutine area_burnt ( currentSite ) !currentPatch%AB daily area burnt (m2) !currentPatch%NF !Daily number of ignitions (lightning and human-caused), adjusted for size of patch. - use domainMod, only : ldomain use EDParamsMod, only : ED_val_nfires - use PatchType, only : patch type(ed_site_type), intent(inout), target :: currentSite - type(ed_patch_type), pointer :: currentPatch real lb !length to breadth ratio of fire ellipse @@ -689,15 +686,11 @@ subroutine area_burnt ( currentSite ) ! --- calculate area burnt--- if(lb > 0.0_r8) then - p = currentPatch%clm_pno - g = patch%gridcell(p) - - ! INTERF-TODO: ! THIS SHOULD HAVE THE COLUMN AND LU AREA WEIGHT ALSO, NO? - gridarea = ldomain%area(g) *1000000.0_r8 !convert from km2 into m2 - currentPatch%NF = ldomain%area(g) * ED_val_nfires * currentPatch%area/area /365 + gridarea = 1000000.0_r8 ! 1M m2 in a km2 + currentPatch%NF = ED_val_nfires * currentPatch%area/area /365 ! If there are 15 lightening strickes per year, per km2. (approx from NASA product) ! then there are 15/365 s/km2 each day. @@ -728,7 +721,7 @@ subroutine area_burnt ( currentSite ) endif! fire currentSite%frac_burnt = currentSite%frac_burnt + currentPatch%frac_burnt - currentPatch => currentPatch%younger; + currentPatch => currentPatch%younger enddo !end patch loop From e2262c482def34bdd6735fc20cdfc5abc1beea6e Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 12 Jan 2017 18:33:58 -0800 Subject: [PATCH 11/11] Simplified the dynamics call sequence by pulling out the litter flux wrapper. Essentially it was just a wrapper that appeared only once within another wrapper. --- .../clm/src/utils/clmfates_interfaceMod.F90 | 78 ++++++++----------- 1 file changed, 34 insertions(+), 44 deletions(-) diff --git a/components/clm/src/utils/clmfates_interfaceMod.F90 b/components/clm/src/utils/clmfates_interfaceMod.F90 index 150b85a159..fd92a13dfe 100644 --- a/components/clm/src/utils/clmfates_interfaceMod.F90 +++ b/components/clm/src/utils/clmfates_interfaceMod.F90 @@ -158,7 +158,6 @@ module CLMFatesInterfaceMod procedure, public :: wrap_accumulatefluxes procedure, public :: prep_canopyfluxes procedure, public :: wrap_canopy_radiation - procedure, private :: wrap_litter_fluxout procedure, public :: wrap_bgc_summary procedure, private :: init_history_io procedure, private :: wrap_update_hlmfates_dyn @@ -478,6 +477,7 @@ subroutine dynamics_driv(this, nc, bounds_clump, & !----------------------------------------------------------------------- ! --------------------------------------------------------------------------------- + ! Part I. ! Prepare input boundary conditions for FATES dynamics ! Note that timing information is the same across all sites, this may ! seem redundant, but it is possible that we may have asynchronous site simulations @@ -512,6 +512,8 @@ subroutine dynamics_driv(this, nc, bounds_clump, & this%fates(nc)%bc_in(s)%t_veg24_si = & temperature_inst%t_veg24_patch(col%patchi(c)) + this%fates(nc)%bc_in(s)%max_rooting_depth_index_col = canopystate_inst%altmax_lastyear_indx_col(c) + do ifp = 1, this%fates(nc)%sites(s)%youngest_patch%patchno p = ifp+col%patchi(c) this%fates(nc)%bc_in(s)%t_veg24_pa(ifp) = & @@ -529,7 +531,11 @@ subroutine dynamics_driv(this, nc, bounds_clump, & end do end do - ! where most things happen + ! --------------------------------------------------------------------------------- + ! Part II: Call the FATES model now that input boundary conditions have been + ! provided. + ! --------------------------------------------------------------------------------- + do s = 1,this%fates(nc)%nsites call ed_ecosystem_dynamics(this%fates(nc)%sites(s), & @@ -537,12 +543,35 @@ subroutine dynamics_driv(this, nc, bounds_clump, & call ed_update_site(this%fates(nc)%sites(s), & this%fates(nc)%bc_in(s)) - + enddo + + ! call subroutine to aggregate ED litter output fluxes and + ! package them for handing across interface + call flux_into_litter_pools(this%fates(nc)%nsites, & + this%fates(nc)%sites, & + this%fates(nc)%bc_in, & + this%fates(nc)%bc_out) - call this%wrap_litter_fluxout(nc, bounds_clump, canopystate_inst, soilbiogeochem_carbonflux_inst) ! --------------------------------------------------------------------------------- + ! Part III: Process FATES output into the dimensions and structures that are part + ! of the HLMs API. (column, depth, and litter fractions) + ! --------------------------------------------------------------------------------- + + do s = 1, this%fates(nc)%nsites + c = this%f2hmap(nc)%fcolumn(s) + soilbiogeochem_carbonflux_inst%FATES_c_to_litr_lab_c_col(c,:) = & + this%fates(nc)%bc_out(s)%FATES_c_to_litr_lab_c_col(:) + soilbiogeochem_carbonflux_inst%FATES_c_to_litr_cel_c_col(c,:) = & + this%fates(nc)%bc_out(s)%FATES_c_to_litr_cel_c_col(:) + soilbiogeochem_carbonflux_inst%FATES_c_to_litr_lig_c_col(c,:) = & + this%fates(nc)%bc_out(s)%FATES_c_to_litr_lig_c_col(:) + end do + + + ! --------------------------------------------------------------------------------- + ! Part III.2 (continued). ! Update diagnostics of the FATES ecosystem structure that are used in the HLM. ! --------------------------------------------------------------------------------- call this%wrap_update_hlmfates_dyn(nc, & @@ -551,6 +580,7 @@ subroutine dynamics_driv(this, nc, bounds_clump, & canopystate_inst) ! --------------------------------------------------------------------------------- + ! Part IV: ! Update history IO fields that depend on ecosystem dynamics ! --------------------------------------------------------------------------------- call this%fates_hist%update_history_dyn( nc, & @@ -1539,46 +1569,6 @@ subroutine wrap_canopy_radiation(this, bounds_clump, nc, & end subroutine wrap_canopy_radiation - ! ====================================================================================== - - subroutine wrap_litter_fluxout(this, nc, bounds_clump, canopystate_inst, soilbiogeochem_carbonflux_inst) - - implicit none - - ! Arguments - class(hlm_fates_interface_type), intent(inout) :: this - integer , intent(in) :: nc - type(bounds_type),intent(in) :: bounds_clump - type(canopystate_type) , intent(inout) :: canopystate_inst - type(soilbiogeochem_carbonflux_type), intent(inout) :: soilbiogeochem_carbonflux_inst - - ! local variables - integer :: s, c - - - ! process needed input boundary conditions to define rooting profiles - ! call subroutine to aggregate ED litter output fluxes and package them for handing across interface - ! process output into the dimensions that the BGC model wants (column, depth, and litter fractions) - - do s = 1, this%fates(nc)%nsites - c = this%f2hmap(nc)%fcolumn(s) - this%fates(nc)%bc_in(s)%max_rooting_depth_index_col = canopystate_inst%altmax_lastyear_indx_col(c) - end do - - call flux_into_litter_pools(this%fates(nc)%nsites, & - this%fates(nc)%sites, & - this%fates(nc)%bc_in, & - this%fates(nc)%bc_out) - - do s = 1, this%fates(nc)%nsites - c = this%f2hmap(nc)%fcolumn(s) - soilbiogeochem_carbonflux_inst%FATES_c_to_litr_lab_c_col(c,:) = this%fates(nc)%bc_out(s)%FATES_c_to_litr_lab_c_col(:) - soilbiogeochem_carbonflux_inst%FATES_c_to_litr_cel_c_col(c,:) = this%fates(nc)%bc_out(s)%FATES_c_to_litr_cel_c_col(:) - soilbiogeochem_carbonflux_inst%FATES_c_to_litr_lig_c_col(c,:) = this%fates(nc)%bc_out(s)%FATES_c_to_litr_lig_c_col(:) - end do - - end subroutine wrap_litter_fluxout - ! ====================================================================================== subroutine wrap_bgc_summary(this, nc, soilbiogeochem_carbonflux_inst, &