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/EDCanopyStructureMod.F90 b/components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90 index 20ae7efe54..3b8d94e6ed 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 012926ab09..2237553ccd 100755 --- a/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 +++ b/components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90 @@ -6,14 +6,16 @@ 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 use pftconMod , only : pftcon use EDEcophysContype , only : EDecophyscon 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 ! @@ -73,7 +75,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 +83,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 +111,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 @@ -221,6 +223,7 @@ subroutine allocate_live_biomass(cc_p,mode) ! Use different proportions if the leaves are on vs off if(leaves_off_switch==0)then + new_bl = currentcohort%balive*leaf_frac new_br = pftcon%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac @@ -233,13 +236,13 @@ subroutine allocate_live_biomass(cc_p,mode) if(mode==1)then currentcohort%npp_leaf = currentcohort%npp_leaf + & - max(0.0_r8,new_bl - currentcohort%bl) / udata%deltat + max(0.0_r8,new_bl - currentcohort%bl) / freq_day currentcohort%npp_froot = currentcohort%npp_froot + & - max(0._r8,new_br - currentcohort%br) / udata%deltat - - currentcohort%npp_bsw = max(0._r8,new_bsw - currentcohort%bsw)/udata%deltat - + max(0._r8,new_br - currentcohort%br) / freq_day + + currentcohort%npp_bsw = max(0.0_r8, new_bsw - currentcohort%bsw)/freq_day + currentcohort%npp_bdead = currentCohort%dbdeaddt end if @@ -271,9 +274,9 @@ subroutine allocate_live_biomass(cc_p,mode) if(mode==1)then currentcohort%npp_froot = currentcohort%npp_froot + & - max(0.0_r8,new_br-currentcohort%br)/udata%deltat + max(0.0_r8,new_br-currentcohort%br)/freq_day - currentcohort%npp_bsw = max(0.0_r8, new_bsw-currentcohort%bsw)/udata%deltat + currentcohort%npp_bsw = max(0.0_r8, new_bsw-currentcohort%bsw)/freq_day currentcohort%npp_bdead = currentCohort%dbdeaddt @@ -1008,8 +1011,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/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 ee269e69cb..430da23a4b 100755 --- a/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 +++ b/components/clm/src/ED/biogeochem/EDPhysiologyMod.F90 @@ -6,14 +6,15 @@ module EDPhysiologyMod ! Miscellaneous physiology routines from ED. ! ============================================================================ - 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 + 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 + 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 @@ -52,7 +53,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 @@ -62,6 +63,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 @@ -72,14 +74,14 @@ 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 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 @@ -89,8 +91,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 @@ -117,11 +120,12 @@ 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) + & - (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 @@ -130,19 +134,12 @@ subroutine non_canopy_derivs( currentSite, currentPatch, temperature_inst ) 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 ! ============================================================================ @@ -183,7 +180,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 @@ -208,7 +205,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. @@ -226,7 +223,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 @@ -235,7 +232,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. @@ -247,25 +244,22 @@ subroutine trim_canopy( currentSite ) end subroutine trim_canopy ! ============================================================================ - subroutine phenology( currentSite, temperature_inst, waterstate_inst) + subroutine phenology( currentSite, bc_in ) ! ! !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 EDTypesMod, only : udata - use PatchType , only : patch + use FatesConstantsMod, only : tfrz => t_water_freeze_k_1atm + ! ! !ARGUMENTS: - type(ed_site_type) , intent(inout), target :: currentSite - type(temperature_type) , intent(in) :: temperature_inst - type(waterstate_type) , intent(in) :: waterstate_inst + type(ed_site_type), intent(inout), target :: currentSite + type(bc_in_type), intent(in) :: bc_in + ! ! !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 @@ -278,8 +272,6 @@ subroutine phenology( currentSite, temperature_inst, 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. @@ -290,27 +282,7 @@ 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 - - !------------------------------------------------------------------------ - - ! 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 - - 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 + real(r8), parameter :: mindayson = 30.0 ! Parameter of drought decid leaf loss in mm in top layer...FIX(RF,032414) ! - this is arbitrary and poorly understood. Needs work. ED_ @@ -322,15 +294,13 @@ subroutine phenology( currentSite, temperature_inst, waterstate_inst) b = 638.0_r8 c = -0.001_r8 coldday = 5.0_r8 !ed_ph_chiltemp - - mindayson = 30 !Parameters from SDGVM model of senesence ncolddayslim = 5 cold_t = 7.5_r8 ! ed_ph_coldtemp - t = udata%time_period - temp_in_C = t_veg24(patchi) - tfrz + t = day_of_year + temp_in_C = bc_in%t_veg24_si - tfrz !-----------------Cold Phenology--------------------! @@ -374,12 +344,12 @@ subroutine phenology( currentSite, 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 - timesinceleafoff = modelday - 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 @@ -388,14 +358,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' + if ( DEBUG ) write(fates_log(),*) 'leaves on' endif !ncd endif !status endif !GDD - timesinceleafon = modelday - currentSite%leafondate + timesinceleafon = model_day - currentSite%leafondate !LEAF OFF: COLD THRESHOLD @@ -409,8 +379,8 @@ 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 - if ( DEBUG ) write(iulog,*) 'leaves off' + currentSite%leafoffdate = model_day !record leaf off date + if ( DEBUG ) write(fates_log(),*) 'leaves off' endif endif endif @@ -419,8 +389,8 @@ 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 - if ( DEBUG ) write(iulog,*) 'leaves off' + currentSite%leafoffdate = model_day !record leaf off date + if ( DEBUG ) write(fates_log(),*) 'leaves off' endif endif @@ -452,7 +422,7 @@ subroutine phenology( currentSite, temperature_inst, 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 @@ -536,7 +506,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. !------------------------------------------------------------------------ @@ -565,11 +536,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 @@ -604,11 +575,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 @@ -666,7 +637,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 @@ -675,8 +647,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 @@ -732,24 +706,26 @@ 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 ! ============================================================================ - subroutine Growth_Derivatives( currentSite, currentCohort) + subroutine Growth_Derivatives( currentSite, currentCohort, bc_in) ! ! !DESCRIPTION: ! Main subroutine controlling growth and allocation derivatives ! ! !USES: use EDGrowthFunctionsMod , only : Bleaf, dDbhdBd, dhdbd, hite, mortality_rates,dDbhdBl - use EDTypesMod , only : udata + ! ! !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 @@ -794,11 +770,13 @@ 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 - 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 * 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 @@ -828,7 +806,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 @@ -840,7 +818,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 - & @@ -856,7 +834,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 @@ -867,8 +845,9 @@ subroutine Growth_Derivatives( currentSite, currentCohort) !what fraction of allocation do we divert to storage? !what is the flux into the store? currentCohort%storage_flux = currentCohort%carbon_balance * f_store + currentCohort%npp_store = 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) @@ -939,7 +918,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 @@ -953,28 +932,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, & + 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 - if ( DEBUG ) write(iulog,*) 'EDPhys dbstoredt II ',currentCohort%dbstoredt + if ( DEBUG ) write(fates_log(),*) 'EDPhys dbstoredt II ',currentCohort%dbstoredt currentCohort%dbalivedt = 0._r8 endif @@ -998,7 +977,6 @@ subroutine recruitment( t, currentSite, currentPatch ) ! ! !USES: use EDGrowthFunctionsMod, only : bdead,dbh, Bleaf - use EDTypesMod, only : udata ! ! !ARGUMENTS integer, intent(in) :: t @@ -1025,14 +1003,14 @@ 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 - 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 @@ -1051,7 +1029,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) @@ -1072,7 +1050,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 @@ -1102,7 +1080,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. @@ -1121,7 +1099,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 @@ -1132,7 +1110,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 @@ -1150,41 +1128,42 @@ 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 - ! 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 - 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 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), pointer :: t_veg24(:) + 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 )) - - t_veg24 => temperature_inst%t_veg24_patch ! Input: [real(r8) (:)] avg pft vegetation temperature for last 24 hrs - + catanf(t1) = 11.75_r8 +(29.7_r8 / pi) * atan( pi * 0.031_r8 * ( t1 - 15.4_r8 )) catanf_30 = catanf(30._r8) - p = currentPatch%clm_pno + ifp = currentPatch%patchno ! set "froz_q10" parameter froz_q10 = EDParamsShareInst%froz_q10 @@ -1193,20 +1172,22 @@ 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) - ! 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**((t_veg24(p)-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(t_veg24(p)-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 @@ -1215,7 +1196,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 @@ -1223,12 +1204,13 @@ subroutine cwd_out( currentSite, currentPatch, temperature_inst ) ! ! !USES: use SFParamsMod, only : SF_val_max_decomp - use EDTypesMod , only : udata + ! ! !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 @@ -1236,8 +1218,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. @@ -1264,19 +1246,19 @@ subroutine cwd_out( currentSite, currentPatch, temperature_inst ) 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 !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 @@ -1305,14 +1287,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 ! @@ -1371,7 +1354,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 @@ -1380,8 +1363,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 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! @@ -1505,13 +1490,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 @@ -1521,7 +1506,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 @@ -1589,12 +1574,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. @@ -1655,15 +1640,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/fire/SFMainMod.F90 b/components/clm/src/ED/fire/SFMainMod.F90 index be53100a71..b6ff07c79a 100755 --- a/components/clm/src/ED/fire/SFMainMod.F90 +++ b/components/clm/src/ED/fire/SFMainMod.F90 @@ -5,15 +5,25 @@ 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 spmdMod , only : masterproc - use clm_varctl , only : iulog - use atm2lndType , only : atm2lnd_type - use TemperatureType , only : temperature_type + use FatesConstantsMod , only : r8 => fates_r8 + +! 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 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 +52,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 +72,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 +91,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 +111,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 +184,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 ( 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(iulog,*) '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 @@ -199,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(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 ( 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 @@ -210,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(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 ( 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? @@ -227,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(iulog,*) '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) @@ -254,14 +259,14 @@ subroutine charecteristics_of_fuel ( currentSite ) if(write_SF == 1)then - if (masterproc) write(iulog,*) '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(iulog,*) '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. @@ -277,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(iulog,*) 'problem with spitfire fuel averaging' + if ( cp_masterproc == 1 ) write(fates_log(),*) 'problem with spitfire fuel averaging' endif currentPatch => currentPatch%younger @@ -288,33 +293,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 ( 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 @@ -328,7 +335,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 +347,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 +360,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 ( cp_masterproc == 1 ) write(fates_log(),*) 'grass, trees, bare',grass_fraction, tree_fraction, bare_fraction endif currentPatch=>currentSite%oldest_patch; @@ -403,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(iulog,*) 'SF - currentPatch%fuel_bulkd ',currentPatch%fuel_bulkd - if (masterproc.and.DEBUG) write(iulog,*) '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(iulog,*) 'SF - beta ',beta - if (masterproc.and.DEBUG) write(iulog,*) '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(iulog,*) '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 @@ -432,11 +439,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 ( 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 @@ -464,18 +471,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 +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(iulog,*) '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 @@ -609,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(iulog,*) '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 @@ -635,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 @@ -682,15 +686,11 @@ subroutine area_burnt ( currentSite ) ! --- calculate area burnt--- if(lb > 0.0_r8) then - 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? - 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. @@ -703,24 +703,25 @@ 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 ( cp_masterproc == 1 ) 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 ( cp_masterproc == 1 ) 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 ( 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(iulog,*) 'frac_burnt',currentPatch%frac_burnt + if ( cp_masterproc == 1 ) write(fates_log(),*) 'frac_burnt',currentPatch%frac_burnt endif endif endif! fire currentSite%frac_burnt = currentSite%frac_burnt + currentPatch%frac_burnt - currentPatch => currentPatch%younger; + currentPatch => currentPatch%younger enddo !end patch loop @@ -771,7 +772,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 ( 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 e8830e4161..76bc5ed9b2 100755 --- a/components/clm/src/ED/main/EDInitMod.F90 +++ b/components/clm/src/ED/main/EDInitMod.F90 @@ -4,23 +4,19 @@ 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 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 @@ -281,7 +277,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 9499f93d02..6882222f1e 100755 --- a/components/clm/src/ED/main/EDMainMod.F90 +++ b/components/clm/src/ED/main/EDMainMod.F90 @@ -6,17 +6,25 @@ module EDMainMod use shr_kind_mod , only : r8 => shr_kind_r8 - use clm_varctl , only : iulog + 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 - 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 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 @@ -39,24 +47,22 @@ module EDMainMod contains !-------------------------------------------------------------------------------! - subroutine ed_ecosystem_dynamics(currentSite, & - atm2lnd_inst, & - soilstate_inst, temperature_inst, waterstate_inst) + subroutine ed_ecosystem_dynamics(currentSite, bc_in) ! ! !DESCRIPTION: ! Core of ed model, calling all subsequent vegetation dynamics routines ! ! !ARGUMENTS: type(ed_site_type) , intent(inout), target :: currentSite - 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 + type(bc_in_type) , intent(in) :: bc_in ! ! !LOCAL VARIABLES: type(ed_patch_type), pointer :: currentPatch !----------------------------------------------------------------------- + if ( cp_masterproc==1 ) write(fates_log(),'(A,I4,A,I2.2,A,I2.2)') 'FATES Dynamics: ',& + current_year,'-',current_month,'-',current_day + !************************************************************************** ! Fire, growth, biogeochemistry. !************************************************************************** @@ -66,15 +72,15 @@ subroutine ed_ecosystem_dynamics(currentSite, & call ed_total_balance_check(currentSite, 0) - call phenology(currentSite, temperature_inst, waterstate_inst) + 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) ! 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 @@ -131,7 +137,7 @@ subroutine ed_ecosystem_dynamics(currentSite, & 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 @@ -139,8 +145,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 @@ -163,43 +170,43 @@ subroutine ed_integrate_state_variables(currentSite, temperature_inst ) 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(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 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(iulog,*) 'EDMainMod dbstoredt I ',currentCohort%bstore, & - currentCohort%dbstoredt,udata%deltat + write(fates_log(),*) 'EDMainMod dbstoredt I ',currentCohort%bstore, & + 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(iulog,*) 'EDMainMod dbstoredt II ',currentCohort%bstore, & - currentCohort%dbstoredt,udata%deltat + write(fates_log(),*) 'EDMainMod dbstoredt II ',currentCohort%bstore, & + currentCohort%dbstoredt,freq_day 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+ & + 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(iulog,*) 'issue with c balance in integration', abs(currentCohort%balive+currentCohort%bdead+ & - currentCohort%bstore+udata%deltat* & + write(fates_log(),*) 'issue with c balance in integration', abs(currentCohort%balive+currentCohort%bdead+ & + currentCohort%bstore+freq_day* & (currentCohort%md+currentCohort%seed_prod)-cohort_biomass_store-currentCohort%npp_acc) endif @@ -218,41 +225,43 @@ 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. ! 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 if(currentPatch%cwd_ag(c) 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 @@ -273,13 +282,13 @@ subroutine ed_integrate_state_variables(currentSite, temperature_inst ) ! 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. do ft = 1,numpft_ed if(currentSite%seed_bank(ft) currentPatch%younger @@ -344,8 +354,10 @@ 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 this is the second to last day of the year, then perform trimming + if( day_of_year == days_per_year-1) then + + write(fates_log(),*) 'calling trim canopy' call trim_canopy(currentSite) endif @@ -415,14 +427,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 07776b7b6f..1817de3e66 100755 --- a/components/clm/src/ED/main/EDTypesMod.F90 +++ b/components/clm/src/ED/main/EDTypesMod.F90 @@ -150,15 +150,17 @@ module EDTypesMod 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 + ! Module switches (this will be read in one day) ! This variable only exists now to serve as a place holder !!!!!!!!!! THIS SHOULD NOT BE SET TO TRUE !!!!!!!!!!!!!!!!! logical,parameter :: use_fates_plant_hydro = .false. - - - !************************************ !** COHORT type structure ** !************************************ @@ -555,16 +557,14 @@ module EDTypesMod !** Userdata type structure ** !************************************ - 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 +! 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 !-------------------------------------------------------------------------------------! public :: ed_hist_scpfmaps diff --git a/components/clm/src/ED/main/FatesConstantsMod.F90 b/components/clm/src/ED/main/FatesConstantsMod.F90 index 4368db1756..bf1f7e562f 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 @@ -60,7 +63,10 @@ module FatesConstantsMod real(fates_r8), parameter :: fates_tiny = tiny(g_per_kg) - + ! Geometric Constants + + ! PI + real(fates_r8), parameter :: pi_const = 3.14159265359_fates_r8 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 139dfb7fac..272bbfbc38 100644 --- a/components/clm/src/ED/main/FatesInterfaceMod.F90 +++ b/components/clm/src/ED/main/FatesInterfaceMod.F90 @@ -9,27 +9,20 @@ 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, & - maxPatchesPerCol, & - cp_nclmax, & - cp_numSWb, & - cp_numlevgrnd, & - cp_maxSWb, & - cp_numlevdecomp, & - cp_numlevdecomp_full, & - cp_hlm_name, & - cp_hio_ignore_val, & - cp_numlevsoil + use EDtypesMod , only : ed_site_type + use EDtypesMod , only : maxPatchesPerCol + 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 - use shr_kind_mod , only : r8 => shr_kind_r8 ! INTERF-TODO: REMOVE THIS - - implicit none ! ------------------------------------------------------------------------------------ @@ -41,15 +34,50 @@ 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 ! ------------------------------------------------------------------------------------ + + + type, public :: bc_in_type ! The actual number of FATES' ED patches integer :: npatches + ! Vegetation Dynamics + ! --------------------------------------------------------------------------------- + + ! 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 + + ! 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 + + ! 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 ! --------------------------------------------------------------------------------- @@ -344,6 +372,13 @@ subroutine allocate_bcin(bc_in) ! Allocate input boundaries + ! Vegetation Dynamics + allocate(bc_in%t_veg24_pa(maxPatchesPerCol)) + + 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)) allocate(bc_in%solai_parb(maxPatchesPerCol,cp_numSWb)) @@ -445,6 +480,13 @@ subroutine zero_bcs(this,s) ! Input boundaries + 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 this%bc_in(s)%smp_gl(:) = 0.0_r8 @@ -552,6 +594,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') @@ -563,6 +606,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' @@ -633,36 +684,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 0169035022..f4edad4514 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 @@ -157,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 @@ -209,6 +209,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 @@ -249,6 +250,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') @@ -449,65 +456,122 @@ 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 + 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) 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 :: jan01_curr_year + integer :: reference_date + integer :: days_per_year + real(r8) :: model_day + 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) + ! 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 + ! one day. The cost of holding site level boundary conditions is minimal + ! and it keeps all the boundaries in one location ! --------------------------------------------------------------------------------- - - ! 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 + + 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) - 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 ) + reference_date = yr*10000 + mon*100 + day + + call timemgr_datediff(reference_date, sec, current_date, current_tod, model_day) + + 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)%h2osoi_vol_si = & + waterstate_inst%h2osoi_vol_col(c,1) - ! TODO-INTEF: PROCEDURE FOR CONVERTING CLM/ALM FIELDS TO MODEL BOUNDARY - ! CONDITIONS. IE. + 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) = & + temperature_inst%t_veg24_patch(p) + + this%fates(nc)%bc_in(s)%precip24_pa(ifp) = & + atm2lnd_inst%prec24_patch(p) + + 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 + + ! --------------------------------------------------------------------------------- + ! Part II: Call the FATES model now that input boundary conditions have been + ! provided. + ! --------------------------------------------------------------------------------- - ! where most things happen do s = 1,this%fates(nc)%nsites call ed_ecosystem_dynamics(this%fates(nc)%sites(s), & - atm2lnd_inst, & - soilstate_inst, temperature_inst, waterstate_inst) + this%fates(nc)%bc_in(s)) + + call ed_update_site(this%fates(nc)%sites(s), & + this%fates(nc)%bc_in(s)) - call ed_update_site(this%fates(nc)%sites(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, & @@ -516,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, & @@ -524,7 +589,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 @@ -874,7 +939,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 ! ------------------------------------------------------------------------ @@ -936,7 +1002,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 @@ -1502,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, &