module FatesInventoryInitMod !----------------------------------------------------------------------------------------------- ! This module contains the majority of the code used to read forest inventory data and ! initialize a simulation from that data. Some important points: ! - This procedure is called through the host model's "cold-start" initialization and not a ! restart type of simulation. ! - This procedure, if called from a massive grid, is both probably inappropriate, and probably ! time consuming. ! - This procedure is assumed to be called over a small subset of sites, for instance a single ! site, or a small collection of sparse/irregularly spaced group of sites ! ! Created: Ryan Knox June 2017 ! This code borrows heavily in concept from what is done in ED2. We will also do our best to ! maintain compatibility with the PSS/CSS file formats that were used in ED2. ! See: https://github.com/EDmodel/ED2/blob/master/ED/src/io/ed_read_ed10_20_history.f90 ! At the time of writing this ED2 is unlicensed, and only concepts were borrowed with no direct ! code copied. !----------------------------------------------------------------------------------------------- ! CIME GLOBALS use shr_log_mod , only : errMsg => shr_log_errMsg ! FATES GLOBALS use FatesConstantsMod, only : r8 => fates_r8 use FatesConstantsMod, only : pi_const use FatesConstantsMod, only : itrue use FatesGlobals , only : endrun => fates_endrun use FatesGlobals , only : fates_log use FatesInterfaceMod, only : bc_in_type use FatesInterfaceMod, only : hlm_inventory_ctrl_file use FatesInterfaceMod, only : nleafage use FatesLitterMod , only : litter_type 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 : leaves_on use EDTypesMod , only : leaves_off use EDTypesMod , only : num_elements use EDTypesMod , only : element_list use EDTypesMod , only : phen_cstat_nevercold use EDTypesMod , only : phen_cstat_iscold use EDTypesMod , only : phen_dstat_timeoff use EDTypesMod , only : phen_dstat_moistoff use EDPftvarcon , only : EDPftvarcon_inst use FatesInterfaceMod, only : hlm_parteh_mode use EDCohortDynamicsMod, only : InitPRTObject use PRTGenericMod, only : prt_carbon_allom_hyp use PRTGenericMod, only : prt_cnp_flex_allom_hyp use PRTGenericMod, only : prt_vartypes use PRTGenericMod, only : leaf_organ use PRTGenericMod, only : fnrt_organ use PRTGenericMod, only : sapw_organ use PRTGenericMod, only : store_organ use PRTGenericMod, only : struct_organ use PRTGenericMod, only : repro_organ use PRTGenericMod, only : carbon12_element use PRTGenericMod, only : nitrogen_element use PRTGenericMod, only : phosphorus_element use PRTGenericMod, only : SetState use FatesConstantsMod, only : primaryforest implicit none private ! This derived type is to allow an array of pointers to the LL patch structure ! This is different than allocating a vector of patches. This is needed for ! quickly matching cohort string identifiers, the indices that match thos identifiers ! with a patch. BY having a vector of patch pointers that lines up with the string ! identifier array, this can be done quickly. type pp_array type(ed_patch_type), pointer :: cpatch end type pp_array character(len=*), parameter, private :: sourcefile = __FILE__ logical, parameter :: debug_inv = .false. ! Debug flag for devs ! String length specifiers integer, parameter :: patchname_strlen = 64 integer, parameter :: cohortname_strlen = 64 integer, parameter :: line_strlen = 512 integer, parameter :: path_strlen = 256 real(r8), parameter :: max_site_adjacency_deg = 0.05_r8 ! The maximum distance in degrees ! allowed between a site's coordinate ! defined in model memory and a physical ! site listed in the file logical, parameter :: do_inventory_out = .false. public :: initialize_sites_by_inventory contains ! ============================================================================================== subroutine initialize_sites_by_inventory(nsites,sites,bc_in) ! !USES: use shr_file_mod, only : shr_file_getUnit use shr_file_mod, only : shr_file_freeUnit use EDPatchDynamicsMod, only : create_patch use EDPatchDynamicsMod, only : fuse_patches use EDCohortDynamicsMod, only : fuse_cohorts use EDCohortDynamicsMod, only : sort_cohorts use EDcohortDynamicsMod, only : count_cohorts ! Arguments integer, intent(in) :: nsites type(ed_site_type), intent(inout), target :: sites(nsites) type(bc_in_type), intent(in) :: bc_in(nsites) ! Locals type(ed_site_type), pointer :: currentSite type(ed_patch_type), pointer :: currentpatch type(ed_cohort_type), pointer :: currentcohort type(ed_patch_type), pointer :: newpatch type(ed_patch_type), pointer :: olderpatch integer :: sitelist_file_unit ! fortran file unit for site list integer :: pss_file_unit ! fortran file unit for the pss file integer :: css_file_unit ! fortran file unit for the css file integer :: nfilesites ! number of sites in file list logical :: lopen ! logical, file is open logical :: lexist ! logical, file exists integer :: ios ! integer, "IO" status character(len=line_strlen) :: header_str ! large string for whole lines real(r8) :: age_init ! dummy value for creating a patch real(r8) :: area_init ! dummy value for creating a patch integer :: s ! site index integer :: ipa ! patch index integer :: total_cohorts ! cohort counter for error checking integer, allocatable :: inv_format_list(:) ! list of format specs character(len=path_strlen), allocatable :: inv_css_list(:) ! list of css file names character(len=path_strlen), allocatable :: inv_pss_list(:) ! list of pss file names real(r8), allocatable :: inv_lat_list(:) ! list of lat coords real(r8), allocatable :: inv_lon_list(:) ! list of lon coords integer :: invsite ! index of inventory site ! closest to actual site integer :: el ! loop counter for number of elements character(len=patchname_strlen) :: patch_name ! patch ID string in the file integer :: npatches ! number of patches found in PSS type(pp_array), allocatable :: patch_pointer_vec(:) ! vector of pointers to patch LL character(len=patchname_strlen), allocatable :: patch_name_vec(:) ! vector of patch ID strings real(r8) :: basal_area_postf ! basal area before fusion (m2/ha) real(r8) :: basal_area_pref ! basal area after fusion (m2/ha) ! I. Load the inventory list file, do some file handle checks ! ------------------------------------------------------------------------------------------ sitelist_file_unit = shr_file_getUnit() inquire(file=trim(hlm_inventory_ctrl_file),exist=lexist,opened=lopen) if( .not.lexist ) then ! The inventory file list DNE write(fates_log(), *) 'An inventory Initialization was requested.' write(fates_log(), *) 'However the inventory file: ',trim(hlm_inventory_ctrl_file),' DNE' write(fates_log(), *) 'Aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end if if( lopen ) then ! The inventory file should not be open write(fates_log(), *) 'The inventory list file is open but should not be.' write(fates_log(), *) 'Aborting.' call endrun(msg=errMsg(sourcefile, __LINE__)) end if open(unit=sitelist_file_unit,file=trim(hlm_inventory_ctrl_file),status='OLD',action='READ',form='FORMATTED') rewind(sitelist_file_unit) ! There should be at least 1 line read(sitelist_file_unit,fmt='(A)',iostat=ios) header_str read(sitelist_file_unit,fmt='(A)',iostat=ios) header_str if( ios /= 0 ) then write(fates_log(), *) 'The inventory file does not contain at least two lines' write(fates_log(), *) 'of data, ie a header and 1 site. Aborting.' call endrun(msg=errMsg(sourcefile, __LINE__)) end if rewind(unit=sitelist_file_unit) ! Count the number of sites that are listed in this file, and allocate storage arrays ! ------------------------------------------------------------------------------------------ nfilesites = count_inventory_sites(sitelist_file_unit) allocate(inv_format_list(nfilesites)) allocate(inv_pss_list(nfilesites)) allocate(inv_css_list(nfilesites)) allocate(inv_lat_list(nfilesites)) allocate(inv_lon_list(nfilesites)) ! Check through the sites that are listed and do some sanity checks ! ------------------------------------------------------------------------------------------ call assess_inventory_sites(sitelist_file_unit, nfilesites, inv_format_list, & inv_pss_list, inv_css_list, & inv_lat_list, inv_lon_list) ! We can close the list file now close(sitelist_file_unit, iostat = ios) if( ios /= 0 ) then write(fates_log(), *) 'The inventory file needed to be closed, but was still open' write(fates_log(), *) 'aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end if call shr_file_freeUnit(sitelist_file_unit) ! For each site, identify the most proximal PSS/CSS couplet, read-in the data ! allocate linked lists and assign to memory do s = 1, nsites invsite = & minloc( (sites(s)%lat-inv_lat_list(:))**2.0_r8 + & (sites(s)%lon-inv_lon_list(:))**2.0_r8 , dim=1) ! Do a sanity check on the distance separation between physical site and model site if ( sqrt( (sites(s)%lat-inv_lat_list(invsite))**2.0_r8 + & (sites(s)%lon-inv_lon_list(invsite))**2.0_r8 ) > max_site_adjacency_deg ) then write(fates_log(), *) 'Model site at lat:',sites(s)%lat,' lon:',sites(s)%lon write(fates_log(), *) 'has no reasonably proximal site in the inventory site list.' write(fates_log(), *) 'Closest is at lat:',inv_lat_list(invsite),' lon:',inv_lon_list(invsite) write(fates_log(), *) 'Separation must be less than ',max_site_adjacency_deg,' degrees' write(fates_log(), *) 'Exiting' call endrun(msg=errMsg(sourcefile, __LINE__)) end if ! Open the PSS/CSS couplet and initialize the ED data structures. ! Lets start withe the PSS ! --------------------------------------------------------------------------------------- pss_file_unit = shr_file_getUnit() open(unit=pss_file_unit,file=trim(inv_pss_list(invsite)), & status='OLD',action='READ',form='FORMATTED') rewind(pss_file_unit) read(pss_file_unit,fmt=*) header_str ! Do one quick pass through just to count lines ipa = 0 countpatchloop: do read(pss_file_unit,fmt=*,iostat=ios) header_str if(ios/=0) exit ipa = ipa + 1 end do countpatchloop rewind(pss_file_unit) read(pss_file_unit,fmt=*) header_str npatches = ipa allocate(patch_name_vec(npatches)) allocate(patch_pointer_vec(npatches)) do ipa=1,npatches allocate(newpatch) newpatch%patchno = ipa newpatch%younger => null() newpatch%older => null() ! This call doesn't do much asside from initializing the patch with ! nominal values, NaNs, zero's and allocating some vectors. We should ! be able to get the following values from the patch files. But on ! the patch creation step, we don't have that information. age_init = 0.0_r8 area_init = 0.0_r8 call create_patch(sites(s), newpatch, age_init, area_init, primaryforest ) if( inv_format_list(invsite) == 1 ) then call set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name) end if ! Add it to the site's patch list ! ------------------------------------------------------------------------------------ patch_name_vec(ipa) = trim(patch_name) patch_pointer_vec(ipa)%cpatch => newpatch if(ipa == 1) then ! This is the first patch to be added ! It starts off as the oldest and youngest patch in the list sites(s)%youngest_patch => newpatch sites(s)%oldest_patch => newpatch else ! Insert this patch into the patch LL ! First check the two end-cases ! Youngest Patch if(newpatch%age<=sites(s)%youngest_patch%age)then newpatch%older => sites(s)%youngest_patch newpatch%younger => null() sites(s)%youngest_patch%younger => newpatch sites(s)%youngest_patch => newpatch ! Oldest Patch else if(newpatch%age>sites(s)%oldest_patch%age)then newpatch%older => null() newpatch%younger => sites(s)%oldest_patch sites(s)%oldest_patch%older => newpatch sites(s)%oldest_patch => newpatch ! Somewhere in the middle else currentpatch => sites(s)%youngest_patch do while(associated(currentpatch)) olderpatch => currentpatch%older if(associated(currentpatch%older)) then if(newpatch%age >= currentpatch%age .and. & newpatch%age < olderpatch%age) then ! Set the new patches pointers newpatch%older => currentpatch%older newpatch%younger => currentpatch ! Fix the patch's older pointer currentpatch%older => newpatch ! Fix the older patch's younger pointer olderpatch%younger => newpatch end if end if currentPatch => olderpatch enddo end if end if end do close(pss_file_unit,iostat=ios) if( ios /= 0 ) then write(fates_log(), *) 'The pss file: ',inv_pss_list(invsite),' could not be closed' write(fates_log(), *) 'aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end if call shr_file_freeUnit(pss_file_unit) if(debug_inv) then write(fates_log(),*) 'Raw List of Inventory Patches, Age Sorted:' currentpatch => sites(s)%youngest_patch do while(associated(currentpatch)) write(fates_log(),*) ' AGE: ',currentpatch%age,' AREA: ',currentpatch%area currentPatch => currentpatch%older enddo end if ! OPEN THE CSS FILE ! --------------------------------------------------------------------------------------- css_file_unit = shr_file_getUnit() open(unit=css_file_unit,file=trim(inv_css_list(invsite)), & status='OLD',action='READ',form='FORMATTED') rewind(css_file_unit) read(css_file_unit,fmt=*) header_str ! Read in each cohort line. Each line is associated with a patch from the PSS ! file via a patch name identification string. We pass the whole site pointer ! to this routine, because inside the routine we identify the patch by making ! comparisons with patch_name_vec and identifying the patch pointer ! from patch_pointer_vec invcohortloop: do if ( inv_format_list(invsite) == 1 ) then call set_inventory_edcohort_type1(sites(s),bc_in(s),css_file_unit, & npatches, patch_pointer_vec,patch_name_vec, ios) end if if ( ios/=0 ) exit end do invcohortloop close(css_file_unit,iostat=ios) if( ios/=0 ) then write(fates_log(), *) 'The css file: ',inv_css_list(invsite),' could not be closed' write(fates_log(), *) 'aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end if call shr_file_freeUnit(css_file_unit) deallocate(patch_pointer_vec,patch_name_vec) ! Report Basal Area (as a check on if things were read in) ! ------------------------------------------------------------------------------ basal_area_pref = 0.0_r8 currentpatch => sites(s)%youngest_patch do while(associated(currentpatch)) currentcohort => currentpatch%tallest do while(associated(currentcohort)) basal_area_pref = basal_area_pref + & currentcohort%n*0.25*((currentcohort%dbh/100.0_r8)**2.0_r8)*pi_const currentcohort => currentcohort%shorter end do currentPatch => currentpatch%older enddo write(fates_log(),*) '-------------------------------------------------------' write(fates_log(),*) 'Basal Area from inventory, BEFORE fusion' write(fates_log(),*) 'Lat: ',sites(s)%lat,' Lon: ',sites(s)%lon write(fates_log(),*) basal_area_pref,' [m2/ha]' write(fates_log(),*) '-------------------------------------------------------' ! Update the patch index numbers and fuse the cohorts in the patches ! ---------------------------------------------------------------------------------------- ipa=1 total_cohorts = 0 currentpatch => sites(s)%youngest_patch do while(associated(currentpatch)) currentpatch%patchno = ipa ipa=ipa+1 ! Perform Cohort Fusion call fuse_cohorts(sites(s), currentpatch,bc_in(s)) call sort_cohorts(currentpatch) ! This calculates %countcohorts call count_cohorts(currentpatch) total_cohorts = total_cohorts + currentPatch%countcohorts currentPatch => currentpatch%older enddo if(total_cohorts .eq. 0)then write(fates_log(), *) 'The inventory initialization produced no cohorts.' write(fates_log(), *) 'aborting' call endrun(msg=errMsg(sourcefile, __LINE__)) end if ! Fuse patches ! ---------------------------------------------------------------------------------------- call fuse_patches(sites(s), bc_in(s) ) ! Report Basal Area (as a check on if things were read in) ! ---------------------------------------------------------------------------------------- basal_area_postf = 0.0_r8 currentpatch => sites(s)%youngest_patch do while(associated(currentpatch)) currentcohort => currentpatch%tallest do while(associated(currentcohort)) basal_area_postf = basal_area_postf + & currentcohort%n*0.25*((currentcohort%dbh/100.0_r8)**2.0_r8)*pi_const currentcohort => currentcohort%shorter end do currentPatch => currentpatch%older enddo write(fates_log(),*) '-------------------------------------------------------' write(fates_log(),*) 'Basal Area from inventory, AFTER fusion' write(fates_log(),*) 'Lat: ',sites(s)%lat,' Lon: ',sites(s)%lon write(fates_log(),*) basal_area_postf,' [m2/ha]' write(fates_log(),*) '-------------------------------------------------------' ! If this is flagged as true, the post-fusion inventory will be written to file ! in the run directory. if(do_inventory_out)then call write_inventory_type1(sites(s)) end if end do deallocate(inv_format_list, inv_pss_list, inv_css_list, inv_lat_list, inv_lon_list) return end subroutine initialize_sites_by_inventory ! ============================================================================================== function count_inventory_sites(sitelist_file_unit) result(nsites) ! Simple function that counts the number of lines in the inventory descriptor file ! Arguments integer, intent(in) :: sitelist_file_unit ! Locals character(len=line_strlen) :: header_str character(len=line_strlen) :: site_str integer :: ios integer :: nsites ! Set the file position to the top of the file ! Read in the header line ! Read through sites and check coordinates and file existence rewind(unit=sitelist_file_unit) read(sitelist_file_unit,fmt='(A)') header_str nsites = 0 do read(sitelist_file_unit,fmt='(A)',iostat=ios) site_str if (ios/=0) exit nsites = nsites + 1 end do return end function count_inventory_sites ! ============================================================================================== subroutine assess_inventory_sites(sitelist_file_unit,nsites, inv_format_list, & inv_pss_list,inv_css_list, & inv_lat_list,inv_lon_list) ! ------------------------------------------------------------------------------------------- ! This subroutine looks through the inventory descriptor file ! and line by line reads information about the available inventory ! sites, and saves their information (such as location and file path) ! to arrays. This routine also does some simple checks to make ! sure it is not reading nonsense ! ! File Format for the inventory site file: ! 1 line header ! 1 line listing each available inventory site with the following fields: ! type latitude longitude pss-name css-name ! ! The fields for each site are described as follows: ! ! <short-name> <value-kind> <description> ! ! type integer We will accomodate different file format with different ! field values as the need arises. format 1 will read in ! datasets via "set_inventory_edpatch_type1()", ! "set_inventory_edcohort_type1()" ! ! latitude float The geographic latitude coordinate of the site ! longitude float The geogarphic longitude coordinate of the site ! pss-name string The full path to the patch descriptor file (PSS) ! css-name string The full path to the cohort descriptor file (CSS) ! ------------------------------------------------------------------------------------------- ! Arguments integer, intent(in) :: sitelist_file_unit ! file unit for sitelist integer, intent(in) :: nsites ! number of inventory sites integer, intent(inout) :: inv_format_list(nsites) ! array of formats for each inventory site character(len=path_strlen),intent(inout) :: inv_pss_list(nsites) ! array of pss file paths for each site character(len=path_strlen),intent(inout) :: inv_css_list(nsites) ! array of css file paths for each site real(r8),intent(inout) :: inv_lat_list(nsites) ! array of latitudes for each site real(r8),intent(inout) :: inv_lon_list(nsites) ! array of longitudes for each site ! Locals character(len=line_strlen) :: header_str ! a string to hold the header information character(len=line_strlen) :: site_str ! a string to hold each site-line in the file integer :: isite ! site index integer :: ios ! fortran read status flag character(len=path_strlen) :: pss_file character(len=path_strlen) :: css_file real(r8) :: site_lat ! inventory site latitude real(r8) :: site_lon ! site longitude integer :: iblnk ! Index used for string parsing integer :: file_format ! format type (1=legacy ED pss/css) logical :: lexist ! file existence flag rewind(unit=sitelist_file_unit) read(sitelist_file_unit,fmt='(4A)') header_str do isite=1,nsites ! Read in the whole line read(sitelist_file_unit,fmt='(a)',iostat=ios) site_str ! Parse the format identifier read(site_str,*) file_format ! Parse the latitude site_str = adjustl(site_str) iblnk = index(site_str,' ') site_str = adjustl(site_str(iblnk:)) read(site_str,*) site_lat ! Parse the longitude site_str = adjustl(site_str) iblnk = index(site_str,' ') site_str = adjustl(site_str(iblnk:)) read(site_str,*) site_lon ! Parse the pss file name site_str = adjustl(site_str) iblnk = index(site_str,' ') site_str = adjustl(site_str(iblnk:)) iblnk = index(site_str,' ') read(site_str(:iblnk),fmt='(1A)') pss_file ! Parse the css file name site_str = adjustl(site_str) iblnk = index(site_str,' ') site_str = adjustl(site_str(iblnk:)) read(site_str,fmt='(1A)') css_file if ( site_lat < -90.0_r8 .or. site_lat > 90.0_r8 ) then write(fates_log(), *) 'read invalid latitude: ',site_lat,' from inventory site list' call endrun(msg=errMsg(sourcefile, __LINE__)) end if ! Longitude should be converted to positive coordinate if( site_lon<0.0_r8 ) site_lon = 360.0_r8 + site_lon if ( site_lon < 0.0_r8 .or. site_lon > 360.0_r8 ) then write(fates_log(), *) 'read invalid longitude: ',site_lon,' from inventory site list' call endrun(msg=errMsg(sourcefile, __LINE__)) end if inquire(file=trim(pss_file),exist=lexist) if( .not.lexist ) then write(fates_log(), *) 'the following pss file could not be found:' write(fates_log(), *) trim(pss_file) call endrun(msg=errMsg(sourcefile, __LINE__)) end if inquire(file=trim(css_file),exist=lexist) if( .not.lexist ) then write(fates_log(), *) 'the following css file could not be found:' write(fates_log(), *) trim(css_file) call endrun(msg=errMsg(sourcefile, __LINE__)) end if ! If we have made it to this point, then in all likelihood, the PSS/CSS ! File has probably been correctly identified inv_format_list(isite) = file_format inv_pss_list(isite) = pss_file inv_css_list(isite) = css_file inv_lat_list(isite) = site_lat inv_lon_list(isite) = site_lon end do return end subroutine assess_inventory_sites ! ============================================================================================== subroutine set_inventory_edpatch_type1(newpatch,pss_file_unit,ipa,ios,patch_name) ! -------------------------------------------------------------------------------------------- ! This subroutine reads in a line of an inventory patch file (pss) ! And populates a new patch with that information. ! This routine specifically reads PSS files that are "Type 1" formatted ! ! The file is formatted text, which contains 1 header line to label columns ! and then 1 line for each patch containing the following fields: ! ! time (year) year of measurement ! patch (string) patch id string ! trk (integer) LU type index (0 non-forest, 1 secondary, 2 primary ! age (years) Time since this patch was disturbed (created) ! area (fraction) Fraction of the site occupied by this patch ! water (NA) Water content of soil (NOT USED) ! fsc (kg/m2) Fast Soil Carbon ! stsc (kg/m2) Structural Soil Carbon ! stsl (kg/m2) Structural Soil Lignan ! ssc (kg/m2) Slow Soil Carbon ! psc (NA) Passive Soil Carbon (NOT USED) ! msn (kg/m2) Mineralized Soil Nitrogen ! fsn (kg/m2) Fast Soil Nitrogen ! -------------------------------------------------------------------------------------------- use FatesSizeAgeTypeIndicesMod, only: get_age_class_index use EDtypesMod, only: AREA use SFParamsMod , only : SF_val_CWD_frac ! Arguments type(ed_patch_type),intent(inout), target :: newpatch ! Patch structure integer,intent(in) :: pss_file_unit ! Self explanatory integer,intent(in) :: ipa ! Patch index (line number) integer,intent(out) :: ios ! Return flag character(len=patchname_strlen),intent(out) :: patch_name ! unique string identifier of patch ! Locals type(litter_type),pointer :: litt integer :: el ! index for elements real(r8) :: p_time ! Time patch was recorded real(r8) :: p_trk ! Land Use index (see above descriptions) character(len=patchname_strlen) :: p_name ! unique string identifier of patch real(r8) :: p_age ! Patch age [years] real(r8) :: p_area ! Patch area [fraction] real(r8) :: p_water ! Patch water (unused) real(r8) :: p_fsc ! Patch fast soil carbon real(r8) :: p_stsc ! Patch structural soil carbon real(r8) :: p_stsl ! Patch structural soil lignans real(r8) :: p_ssc ! Patch slow soil carbon real(r8) :: p_psc ! Patch P soil carbon real(r8) :: p_msn ! Patch mean soil nitrogen real(r8) :: p_fsn ! Patch fast soil nitrogen integer :: icwd ! index for counting CWD pools integer :: ipft ! index for counting PFTs real(r8) :: pftfrac ! the inverse of the total number of PFTs character(len=128),parameter :: wr_fmt = & '(F5.2,2X,A4,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2)' read(pss_file_unit,fmt=*,iostat=ios) p_time, p_name, p_trk, p_age, p_area, & p_water,p_fsc, p_stsc, p_stsl, p_ssc, & p_psc, p_msn, p_fsn if (ios/=0) return patch_name = trim(p_name) if( debug_inv) then write(*,fmt=wr_fmt) & p_time, p_name, p_trk, p_age, p_area, & p_water,p_fsc, p_stsc, p_stsl, p_ssc, & p_psc, p_msn, p_fsn end if ! Fill in the patch's memory structures newpatch%age = p_age newpatch%age_class = get_age_class_index(newpatch%age) newpatch%area = p_area*AREA ! The non-litter patch soil variables need to be sent to the HLM ! This is not being sent as of this message (RGK 06-2017) ! --------------------------------------------------------------------- ! The litter and CWD could be estimated from at least two methods ! 1) after reading in the cohort data, assuming a SS turnover rate ! 2) again assuming SS, back out the CWD and Litter flux rates into ! the SSC, STSC and FSC pools that balance with their decomp rates ! and then use those flux rates, to calculate the CWD and litter ! pool sizes based on another SS model where flux out balances with ! mortality and litter fluxes into the non-decomposing pool ! ! This is significant science modeling and does not have a simple ! first hack solution. (RGK 06-2017) ! ---------------------------------------------------------------------- do el=1,num_elements litt => newpatch%litter(el) call litt%InitConditions(init_leaf_fines=0._r8, & init_root_fines=0._r8, & init_ag_cwd=0._r8, & init_bg_cwd=0._r8, & init_seed=0._r8, & init_seed_germ=0._r8) end do return end subroutine set_inventory_edpatch_type1 ! ============================================================================================== subroutine set_inventory_edcohort_type1(csite,bc_in,css_file_unit,npatches, & patch_pointer_vec,patch_name_vec,ios) ! -------------------------------------------------------------------------------------------- ! This subroutine reads in a line of an inventory cohort file (css) ! And populates a new cohort with that information. ! This routine specifically reads CSS files that are "Type 1" formatted ! ! The file formatted text, which contains 1 header line to label columns ! and then 1 line for each cohort containing the following fields: ! ! FILE FORMAT: ! time (year) year of measurement ! patch (string) patch id string associated with this cohort ! index (integer) cohort index ! dbh (cm) diameter at breast height ! height (m) height of the tree ! pft (integer) the plant functional type index (must be consistent with param file) ! n (/m2) The plant number density ! bdead (kgC/plant)The dead biomass per indiv of this cohort (NOT USED) ! balive (kgC/plant)The live biomass per indiv of this cohort (NOT USED) ! avgRG (cm/yr?) Average Radial Growth (NOT USED) ! -------------------------------------------------------------------------------------------- use FatesAllometryMod , only : h_allom use FatesAllometryMod , only : h2d_allom use FatesAllometryMod , only : bagw_allom use FatesAllometryMod , only : bbgw_allom use FatesAllometryMod , only : bleaf use FatesAllometryMod , only : bfineroot use FatesAllometryMod , only : bsap_allom use FatesAllometryMod , only : bdead_allom use FatesAllometryMod , only : bstore_allom use EDCohortDynamicsMod , only : create_cohort use FatesInterfaceMod , only : numpft ! Arguments type(ed_site_type),intent(inout), target :: csite ! current site type(bc_in_type),intent(in) :: bc_in ! boundary conditions integer, intent(in) :: css_file_unit ! Self explanatory integer, intent(in) :: npatches ! number of patches type(pp_array), intent(in) :: patch_pointer_vec(npatches) character(len=patchname_strlen), intent(in) :: patch_name_vec(npatches) integer,intent(out) :: ios ! Return flag ! Locals class(prt_vartypes), pointer :: prt_obj real(r8) :: c_time ! Time patch was recorded character(len=patchname_strlen) :: p_name ! The patch associated with this cohort character(len=cohortname_strlen) :: c_name ! cohort index real(r8) :: c_dbh ! diameter at breast height (cm) real(r8) :: c_height ! tree height (m) integer :: c_pft ! plant functional type index real(r8) :: c_nplant ! plant density (/m2) real(r8) :: c_bdead ! dead biomass (kg) real(r8) :: c_balive ! live biomass (kg) real(r8) :: c_avgRG ! avg radial growth (NOT USED) real(r8) :: site_spread ! initial guess of site spread ! should be quickly re-calculated integer :: cstatus ! cohort status integer,parameter :: rstatus = 0 ! recruit status type(ed_patch_type), pointer :: cpatch ! current patch pointer type(ed_cohort_type), pointer :: temp_cohort ! temporary patch (needed for allom funcs) integer :: ipa ! patch idex integer :: iage integer :: el integer :: element_id logical :: matched_patch ! check if cohort was matched w/ patch real(r8) :: b_agw ! biomass above ground non-leaf [kgC] real(r8) :: b_bgw ! biomass below ground non-leaf [kgC] real(r8) :: b_leaf ! biomass in leaves [kgC] real(r8) :: b_fnrt ! biomass in fine roots [kgC] real(r8) :: b_sapw ! biomass in sapwood [kgC] real(r8) :: b_struct real(r8) :: b_store real(r8) :: a_sapw ! area of sapwood at reference height [m2] real(r8) :: m_struct ! Generic (any element) mass for structure [kg] real(r8) :: m_leaf ! Generic mass for leaf [kg] real(r8) :: m_fnrt ! Generic mass for fine-root [kg] real(r8) :: m_sapw ! Generic mass for sapwood [kg] real(r8) :: m_store ! Generic mass for storage [kg] real(r8) :: m_repro ! Generic mass for reproductive tissues [kg] integer :: i_pft, ncohorts_to_create character(len=128),parameter :: wr_fmt = & '(F7.1,2X,A20,2X,A20,2X,F5.2,2X,F5.2,2X,I4,2X,F5.2,2X,F5.2,2X,F5.2,2X,F5.2)' real(r8), parameter :: abnormal_large_nplant = 1000.0_r8 ! Used to catch bad values real(r8), parameter :: abnormal_large_dbh = 500.0_r8 ! I've never heard of a tree > 3m integer, parameter :: recruitstatus = 0 read(css_file_unit,fmt=*,iostat=ios) c_time, p_name, c_name, c_dbh, c_height, & c_pft, c_nplant, c_bdead, c_balive, c_avgRG if( debug_inv) then write(*,fmt=wr_fmt) & c_time, p_name, c_name, c_dbh, c_height, & c_pft, c_nplant, c_bdead, c_balive, c_avgRG end if if (ios/=0) return ! Identify the patch based on the patch_name matched_patch = .false. do ipa=1,npatches if( trim(p_name) == trim(patch_name_vec(ipa))) then cpatch => patch_pointer_vec(ipa)%cpatch matched_patch = .true. end if end do if(.not.matched_patch)then write(fates_log(), *) 'could not match a cohort with a patch' write(fates_log(),fmt=wr_fmt) & c_time, p_name, c_name, c_dbh, c_height, & c_pft, c_nplant, c_bdead, c_balive, c_avgRG call endrun(msg=errMsg(sourcefile, __LINE__)) end if ! Run some sanity checks on the input data ! pft, nplant and dbh are the critical ones in this format specification ! ------------------------------------------------------------------------------------------- if (c_pft > numpft ) then write(fates_log(), *) 'inventory pft: ',c_pft write(fates_log(), *) 'An inventory cohort file specified a pft index' write(fates_log(), *) 'greater than the maximum specified pfts numpft' call endrun(msg=errMsg(sourcefile, __LINE__)) end if if (c_pft < 0 ) then write(fates_log(), *) 'inventory pft: ',c_pft write(fates_log(), *) 'The inventory produced a cohort with <0 pft index' call endrun(msg=errMsg(sourcefile, __LINE__)) end if if (c_dbh <=0 ) then write(fates_log(), *) 'inventory dbh: ', c_dbh write(fates_log(), *) 'The inventory produced a cohort with <= 0 dbh' call endrun(msg=errMsg(sourcefile, __LINE__)) end if if (c_dbh > abnormal_large_dbh ) then write(fates_log(), *) 'inventory dbh: ', c_nplant write(fates_log(), *) 'The inventory produced a cohort with very large diameter [cm]' call endrun(msg=errMsg(sourcefile, __LINE__)) end if if (c_nplant <=0 ) then write(fates_log(), *) 'inventory nplant: ', c_nplant write(fates_log(), *) 'The inventory produced a cohort with <= 0 density /m2' call endrun(msg=errMsg(sourcefile, __LINE__)) end if if (c_nplant > abnormal_large_nplant ) then write(fates_log(), *) 'inventory nplant: ', c_nplant write(fates_log(), *) 'The inventory produced a cohort with very large density /m2' call endrun(msg=errMsg(sourcefile, __LINE__)) end if if (c_pft .eq. 0 ) then write(fates_log(), *) 'inventory pft: ',c_pft write(fates_log(), *) 'SPECIAL CASE TRIGGERED: PFT == 0 and therefore this subroutine' write(fates_log(), *) 'will assign a cohort with n = n_orig/numpft to every cohort in range 1 to numpft' ncohorts_to_create = numpft else ncohorts_to_create = 1 end if do i_pft = 1,ncohorts_to_create allocate(temp_cohort) ! A temporary cohort is needed because we want to make if (c_pft .ne. 0 ) then ! normal case: assign each cohort to its specified PFT temp_cohort%pft = c_pft else ! special case, make an identical cohort for each PFT temp_cohort%pft = i_pft endif temp_cohort%n = c_nplant * cpatch%area / real(ncohorts_to_create,r8) temp_cohort%dbh = c_dbh call h_allom(c_dbh,temp_cohort%pft,temp_cohort%hite) temp_cohort%canopy_trim = 1.0_r8 call bagw_allom(temp_cohort%dbh,temp_cohort%pft,b_agw) ! Calculate coarse root biomass from allometry call bbgw_allom(temp_cohort%dbh,temp_cohort%pft,b_bgw) ! Calculate the leaf biomass (calculates a maximum first, then applies canopy trim ! and sla scaling factors) call bleaf(temp_cohort%dbh,temp_cohort%pft,temp_cohort%canopy_trim,b_leaf) ! Calculate fine root biomass call bfineroot(temp_cohort%dbh,temp_cohort%pft,temp_cohort%canopy_trim,b_fnrt) ! Calculate sapwood biomass call bsap_allom(temp_cohort%dbh,temp_cohort%pft,temp_cohort%canopy_trim, a_sapw, b_sapw) call bdead_allom( b_agw, b_bgw, b_sapw, temp_cohort%pft, b_struct ) call bstore_allom(temp_cohort%dbh, temp_cohort%pft, temp_cohort%canopy_trim, b_store) temp_cohort%laimemory = 0._r8 cstatus = leaves_on if( EDPftvarcon_inst%season_decid(temp_cohort%pft) == itrue .and. & any(csite%cstatus == [phen_cstat_nevercold,phen_cstat_iscold])) then temp_cohort%laimemory = b_leaf b_leaf = 0._r8 cstatus = leaves_off endif if ( EDPftvarcon_inst%stress_decid(temp_cohort%pft) == itrue .and. & any(csite%dstatus == [phen_dstat_timeoff,phen_dstat_moistoff])) then temp_cohort%laimemory = b_leaf b_leaf = 0._r8 cstatus = leaves_off endif prt_obj => null() call InitPRTObject(prt_obj) do el = 1,num_elements element_id = element_list(el) ! If this is carbon12, then the initialization is straight forward ! otherwise, we use stoichiometric ratios select case(element_id) case(carbon12_element) m_struct = b_struct m_leaf = b_leaf m_fnrt = b_fnrt m_sapw = b_sapw m_store = b_store m_repro = 0._r8 case(nitrogen_element) m_struct = b_struct*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,struct_organ) m_leaf = b_leaf*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,leaf_organ) m_fnrt = b_fnrt*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,fnrt_organ) m_sapw = b_sapw*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,sapw_organ) m_store = b_store*EDPftvarcon_inst%prt_nitr_stoich_p1(temp_cohort%pft,store_organ) m_repro = 0._r8 case(phosphorus_element) m_struct = b_struct*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,struct_organ) m_leaf = b_leaf*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,leaf_organ) m_fnrt = b_fnrt*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,fnrt_organ) m_sapw = b_sapw*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,sapw_organ) m_store = b_store*EDPftvarcon_inst%prt_phos_stoich_p1(temp_cohort%pft,store_organ) m_repro = 0._r8 end select select case(hlm_parteh_mode) case (prt_carbon_allom_hyp,prt_cnp_flex_allom_hyp ) ! Equally distribute leaf mass into available age-bins do iage = 1,nleafage call SetState(prt_obj,leaf_organ, element_id,m_leaf/real(nleafage,r8),iage) end do call SetState(prt_obj,fnrt_organ, element_id, m_fnrt) call SetState(prt_obj,sapw_organ, element_id, m_sapw) call SetState(prt_obj,store_organ, element_id, m_store) call SetState(prt_obj,struct_organ, element_id, m_struct) call SetState(prt_obj,repro_organ, element_id, m_repro) case default write(fates_log(),*) 'Unspecified PARTEH module during inventory intitialization' call endrun(msg=errMsg(sourcefile, __LINE__)) end select end do call prt_obj%CheckInitialConditions() ! Since spread is a canopy level calculation, we need to provide an initial guess here. call create_cohort(csite, cpatch, temp_cohort%pft, temp_cohort%n, temp_cohort%hite, temp_cohort%dbh, & prt_obj, temp_cohort%laimemory, cstatus, rstatus, temp_cohort%canopy_trim, & 1, csite%spread, bc_in) deallocate(temp_cohort) ! get rid of temporary cohort end do return end subroutine set_inventory_edcohort_type1 ! ==================================================================================== subroutine write_inventory_type1(currentSite) ! -------------------------------------------------------------------------------- ! This subroutine writes the cohort/patch inventory type files in the "type 1" ! format. Note that for compatibility with ED2, we chose an old type that has ! both extra unused fields and is missing fields from FATES. THis is not ! a recommended file type for restarting a run. ! The files will have a lat/long tag added to their name, and will be ! generated in the run folder. ! -------------------------------------------------------------------------------- use shr_file_mod, only : shr_file_getUnit use shr_file_mod, only : shr_file_freeUnit ! Arguments type(ed_site_type), target :: currentSite ! Locals type(ed_patch_type), pointer :: currentpatch type(ed_cohort_type), pointer :: currentcohort character(len=128) :: pss_name_out ! output file string character(len=128) :: css_name_out ! output file string integer :: pss_file_out integer :: css_file_out integer :: ilat_int,ilat_dec ! for output string parsing integer :: ilon_int,ilon_dec ! for output string parsing character(len=32) :: patch_str character(len=32) :: cohort_str integer :: ipatch integer :: icohort character(len=1) :: ilat_sign,ilon_sign ! Generate pss/css file name based on the location of the site ilat_int = abs(int(currentSite%lat)) ilat_dec = int(100000*(abs(currentSite%lat) - real(ilat_int,r8))) ilon_int = abs(int(currentSite%lon)) ilon_dec = int(100000*(abs(currentSite%lon) - real(ilon_int,r8))) if(currentSite%lat>=0._r8)then ilat_sign = 'N' else ilat_sign = 'S' end if if(currentSite%lon>=0._r8)then ilon_sign = 'E' else ilon_sign = 'W' end if write(pss_name_out,'(A8,I2.2,A1,I5.5,A1,A1,I3.3,A1,I5.5,A1,A4)') & 'pss_out_',ilat_int,'.',ilat_dec,ilat_sign,'_',ilon_int,'.',ilon_dec,ilon_sign,'.txt' write(css_name_out,'(A8,I2.2,A1,I5.5,A1,A1,I3.3,A1,I5.5,A1,A4)') & 'css_out_',ilat_int,'.',ilat_dec,ilat_sign,'_',ilon_int,'.',ilon_dec,ilon_sign,'.txt' pss_file_out = shr_file_getUnit() css_file_out = shr_file_getUnit() open(unit=pss_file_out,file=trim(pss_name_out), status='UNKNOWN',action='WRITE',form='FORMATTED') open(unit=css_file_out,file=trim(css_name_out), status='UNKNOWN',action='WRITE',form='FORMATTED') write(pss_file_out,*) 'time patch trk age area water fsc stsc stsl ssc psc msn fsn' write(css_file_out,*) 'time patch cohort dbh hite pft nplant bdead alive Avgrg' ipatch=0 currentpatch => currentSite%youngest_patch do while(associated(currentpatch)) ipatch=ipatch+1 write(patch_str,'(A7,i4.4,A)') '<patch_',ipatch,'>' write(pss_file_out,*) '0000 ',trim(patch_str),' 2 ',currentPatch%age,currentPatch%area/AREA, & '0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000' icohort=0 currentcohort => currentpatch%tallest do while(associated(currentcohort)) icohort=icohort+1 write(cohort_str,'(A7,i4.4,A)') '<coh_',icohort,'>' write(css_file_out,*) '0000 ',trim(patch_str),' ',trim(cohort_str), & currentCohort%dbh,0.0,currentCohort%pft,currentCohort%n/currentPatch%area,0.0,0.0,0.0 currentcohort => currentcohort%shorter end do currentPatch => currentpatch%older enddo close(css_file_out) close(pss_file_out) call shr_file_freeUnit(css_file_out) call shr_file_freeUnit(pss_file_out) end subroutine write_inventory_type1 end module FatesInventoryInitMod