From bd0421e6e7ea20105e27a6c4265a0c6690df0cce Mon Sep 17 00:00:00 2001 From: rosiealice Date: Tue, 19 Mar 2024 07:18:11 -0600 Subject: [PATCH 1/2] writing fire emission variables --- src/main/histFileMod.F90 | 8 ++++++++ src/utils/clmfates_interfaceMod.F90 | 13 +++++++++---- 2 files changed, 17 insertions(+), 4 deletions(-) diff --git a/src/main/histFileMod.F90 b/src/main/histFileMod.F90 index fb1a25db37..b67ff853ac 100644 --- a/src/main/histFileMod.F90 +++ b/src/main/histFileMod.F90 @@ -24,6 +24,7 @@ module histFileMod use PatchType , only : patch use EDParamsMod , only : nclmax use EDParamsMod , only : nlevleaf + use EDParamsMod , only : num_emission_compounds use FatesInterfaceTypesMod , only : nlevsclass, nlevage, nlevcoage use FatesInterfaceTypesMod , only : nlevheight use FatesInterfaceTypesMod , only : nlevdamage @@ -2498,6 +2499,7 @@ subroutine htape_create (t, histrest) call ncd_defdim(lnfid, 'fates_levage', nlevage, dimid) call ncd_defdim(lnfid, 'fates_levheight', nlevheight, dimid) call ncd_defdim(lnfid, 'fates_levfuel', nfsc, dimid) + call ncd_defdim(lnfid, 'fates_levemis', num_emission_compounds, dimid) call ncd_defdim(lnfid, 'fates_levcwdsc', ncwd, dimid) call ncd_defdim(lnfid, 'fates_levscpf', nlevsclass*numpft_fates, dimid) call ncd_defdim(lnfid, 'fates_levcapf', nlevcoage*numpft_fates, dimid) @@ -3045,6 +3047,7 @@ subroutine htape_timeconst(t, mode) use FatesInterfaceTypesMod, only : fates_hdim_agmap_levagepft use FatesInterfaceTypesMod, only : fates_hdim_pftmap_levagepft use FatesInterfaceTypesMod, only : fates_hdim_levfuel + use FatesInterfaceTypesMod, only : fates_hdim_levemis use FatesInterfaceTypesMod, only : fates_hdim_levdamage use FatesInterfaceTypesMod, only : fates_hdim_levcwdsc use FatesInterfaceTypesMod, only : fates_hdim_levcan @@ -3155,6 +3158,8 @@ subroutine htape_timeconst(t, mode) long_name='FATES pft number', ncid=nfid(t)) call ncd_defvar(varname='fates_levfuel',xtype=ncd_int, dim1name='fates_levfuel', & long_name='FATES fuel index', ncid=nfid(t)) + call ncd_defvar(varname='fates_levemis',xtype=ncd_int, dim1name='fates_levemis', & + long_name='FATES fire emissions index', ncid=nfid(t)) call ncd_defvar(varname='fates_levcwdsc',xtype=ncd_int, dim1name='fates_levcwdsc', & long_name='FATES cwd size class', ncid=nfid(t)) call ncd_defvar(varname='fates_levcan',xtype=ncd_int, dim1name='fates_levcan', & @@ -3227,6 +3232,7 @@ subroutine htape_timeconst(t, mode) call ncd_io(varname='fates_levheight',data=fates_hdim_levheight, ncid=nfid(t), flag='write') call ncd_io(varname='fates_levpft',data=fates_hdim_levpft, ncid=nfid(t), flag='write') call ncd_io(varname='fates_levfuel',data=fates_hdim_levfuel, ncid=nfid(t), flag='write') + call ncd_io(varname='fates_levemis',data=fates_hdim_levemis, ncid=nfid(t), flag='write') call ncd_io(varname='fates_levcdam',data=fates_hdim_levdamage, ncid=nfid(t), flag='write') call ncd_io(varname='fates_levcwdsc',data=fates_hdim_levcwdsc, ncid=nfid(t), flag='write') call ncd_io(varname='fates_levcan',data=fates_hdim_levcan, ncid=nfid(t), flag='write') @@ -5520,6 +5526,8 @@ subroutine hist_addfld2d (fname, type2d, units, avgflag, long_name, type1d_out, num2d = nlevage case ('fates_levheight') num2d = nlevheight + case ('fates_levemis') + num2d = num_emission_compounds case ('fates_levfuel') num2d = nfsc case ('fates_levcwdsc') diff --git a/src/utils/clmfates_interfaceMod.F90 b/src/utils/clmfates_interfaceMod.F90 index 7039884847..1f856859e2 100644 --- a/src/utils/clmfates_interfaceMod.F90 +++ b/src/utils/clmfates_interfaceMod.F90 @@ -1540,9 +1540,14 @@ subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump, & z0m(p) = this%fates(nc)%bc_out(s)%z0m_pa(ifp) displa(p) = this%fates(nc)%bc_out(s)%displa_pa(ifp) dleaf_patch(p) = this%fates(nc)%bc_out(s)%dleaf_pa(ifp) - end do ! veg pach - if(abs(areacheck - 1.0_r8).gt.1.e-9_r8)then + if(this%fates(nc)%bc_out(s)%fire_emissions_pa(ifp,1).gt.0_r8)then + write(*,*) 'positive fire emissions',this%fates(nc)%bc_out(s)%fire_emissions_pa(ifp,1) + endif + + end do ! veg pach + + if(abs(areacheck - 1.0_r8).gt.1.e-9_r8)then write(iulog,*) 'area wrong in interface',areacheck - 1.0_r8 call endrun(msg=errMsg(sourcefile, __LINE__)) endif @@ -3113,7 +3118,7 @@ subroutine init_history_io(this,bounds_proc) use FatesConstantsMod, only : fates_short_string_length, fates_long_string_length use FatesIOVariableKindMod, only : site_r8, site_soil_r8, site_size_pft_r8 use FatesIOVariableKindMod, only : site_size_r8, site_pft_r8, site_age_r8 - use FatesIOVariableKindMod, only : site_coage_r8, site_coage_pft_r8 + use FatesIOVariableKindMod, only : site_coage_r8, site_coage_pft_r8, site_emis_r8 use FatesIOVariableKindMod, only : site_fuel_r8, site_cwdsc_r8, site_scag_r8 use FatesIOVariableKindMod, only : site_scagpft_r8, site_agepft_r8 use FatesIOVariableKindMod, only : site_can_r8, site_cnlf_r8, site_cnlfpft_r8 @@ -3218,7 +3223,7 @@ subroutine init_history_io(this,bounds_proc) set_lake=0._r8,set_urb=0._r8) case(site_soil_r8, site_size_pft_r8, site_size_r8, site_pft_r8, & - site_age_r8, site_height_r8, site_coage_r8,site_coage_pft_r8, & + site_age_r8, site_emis_r8, site_height_r8, site_coage_r8,site_coage_pft_r8, & site_fuel_r8, site_cwdsc_r8, site_clscpf_r8, & site_can_r8,site_cnlf_r8, site_cnlfpft_r8, site_scag_r8, & site_scagpft_r8, site_agepft_r8, site_elem_r8, site_elpft_r8, & From 1470518e4aca301af76fe5a027107dc5a4b6adfb Mon Sep 17 00:00:00 2001 From: rosiealice Date: Fri, 22 Mar 2024 06:04:04 -0600 Subject: [PATCH 2/2] added unpacking of FATES fire emissions in clmfatesinterfacemod. Required also the passing of the fireemis_inst data type across theHLM-FATES interface. --- src/biogeochem/CNFireEmissionsMod.F90 | 53 +++++++++++++++++++-------- src/main/clm_driver.F90 | 2 +- src/main/clm_initializeMod.F90 | 2 +- src/main/clm_instMod.F90 | 4 +- src/utils/clmfates_interfaceMod.F90 | 39 +++++++++++++------- 5 files changed, 68 insertions(+), 32 deletions(-) diff --git a/src/biogeochem/CNFireEmissionsMod.F90 b/src/biogeochem/CNFireEmissionsMod.F90 index 5a15e138d5..791e238542 100644 --- a/src/biogeochem/CNFireEmissionsMod.F90 +++ b/src/biogeochem/CNFireEmissionsMod.F90 @@ -7,8 +7,10 @@ module CNFireEmissionsMod ! Created by F. Vitt, and revised by F. Li ! !USES: use shr_kind_mod, only : r8 => shr_kind_r8 + use clm_varctl , only : use_fates use abortutils, only : endrun use PatchType, only : patch + use EDParamsMod, only : num_emission_compounds ! for FATES emissions. use decompMod, only : bounds_type use shr_fire_emis_mod, only : shr_fire_emis_comps_n, shr_fire_emis_comp_t, shr_fire_emis_linkedlist use shr_fire_emis_mod, only : shr_fire_emis_mechcomps_n, shr_fire_emis_mechcomps @@ -27,6 +29,8 @@ module CNFireEmissionsMod ! !PUBLIC TYPES: type, public :: fireemis_type real(r8), pointer, public :: fireflx_patch(:,:) ! carbon flux from fire sources (kg/m2/sec) + real(r8), pointer, public :: fates_fire_emissions_patch(:,:) ! FATES calculated emissions from fire sources (kg/m2/sec) + real(r8), pointer, public :: fates_fire_emission_height_patch(:) ! FATES calculated emissions from fire sources (m) real(r8), pointer, public :: ztop_patch(:) ! height of the smoke plume (meters) type(emis_t), pointer, private :: comp(:) ! fire emissions component (corresponds to emis factors table input file) type(emis_t), pointer, private :: mech(:) ! cam-chem mechism species emissions @@ -58,18 +62,19 @@ subroutine Init(this, bounds) real(r8) :: molec_wght type(shr_fire_emis_comp_t), pointer :: emis_cmp - if ( shr_fire_emis_mechcomps_n < 1) return - - call fire_emis_factors_init( shr_fire_emis_factors_file ) + if ( shr_fire_emis_mechcomps_n > 0) then + + call fire_emis_factors_init( shr_fire_emis_factors_file ) - emis_cmp => shr_fire_emis_linkedlist - do while(associated(emis_cmp)) - allocate(emis_cmp%emis_factors(maxveg)) - call fire_emis_factors_get( trim(emis_cmp%name), factors, molec_wght ) - emis_cmp%emis_factors = factors*1.e-3_r8 ! convert g/kg dry fuel to kg/kg - emis_cmp%molec_weight = molec_wght - emis_cmp => emis_cmp%next_emiscomp - enddo + emis_cmp => shr_fire_emis_linkedlist + do while(associated(emis_cmp)) + allocate(emis_cmp%emis_factors(maxveg)) + call fire_emis_factors_get( trim(emis_cmp%name), factors, molec_wght ) + emis_cmp%emis_factors = factors*1.e-3_r8 ! convert g/kg dry fuel to kg/kg + emis_cmp%molec_weight = molec_wght + emis_cmp => emis_cmp%next_emiscomp + enddo + endif call this%InitAllocate(bounds) call this%InitHistory(bounds) @@ -97,14 +102,21 @@ subroutine InitAllocate(this, bounds) allocate(this%totfire%emis(beg:end)) this%totfire%emis(beg:end) = nan + if(use_fates)then + allocate(this%fates_fire_emissions_patch(beg:end,1:num_emission_compounds)) + allocate(this%fates_fire_emission_height_patch(beg:end)) + + this%fates_fire_emissions_patch(beg:end,:) = spval + this%fates_fire_emission_height_patch(beg:end) = spval + endif + if (shr_fire_emis_mechcomps_n>0) then allocate(this%ztop_patch(beg:end)) this%ztop_patch(beg:end) = spval - allocate(this%fireflx_patch(beg:end,shr_fire_emis_mechcomps_n)) this%fireflx_patch(beg:end,:) = spval + allocate(this%mech(shr_fire_emis_mechcomps_n)) - allocate(this%mech(shr_fire_emis_mechcomps_n)) do i = 1, shr_fire_emis_mechcomps_n allocate(this%mech(i)%emis(beg:end)) this%mech(i)%emis(beg:end) = nan @@ -210,6 +222,14 @@ subroutine CNFireEmisUpdate(bounds, num_bgc_vegp, filter_bgc_vegp, cnveg_cf_inst real(r8) :: epsilon ! emission factor [ug m-2 h-1] integer :: i, ii, icomp, imech, n_emis_comps, l, j + + ! initialize to zero ... + if ( use_fates) then ! always initialize the fates emission variables for now. + fireemis_inst%fates_fire_emissions_patch(bounds%begp:bounds%endp,:) = 0._r8 + fireemis_inst%fates_fire_emission_height_patch(bounds%begp:bounds%endp) = 0._r8 + endif + + if ( shr_fire_emis_mechcomps_n < 1) return associate( & @@ -221,8 +241,11 @@ subroutine CNFireEmisUpdate(bounds, num_bgc_vegp, filter_bgc_vegp, cnveg_cf_inst ) ! initialize to zero ... - fire_emis(bounds%begp:bounds%endp,:) = 0._r8 - totfire%emis(bounds%begp:bounds%endp) = 0._r8 + if ( use_fates) then + fire_emis(bounds%begp:bounds%endp,:) = 0._r8 + totfire%emis(bounds%begp:bounds%endp) = 0._r8 + endif + ztop(bounds%begp:bounds%endp) = 0._r8 do i = 1, shr_fire_emis_mechcomps_n diff --git a/src/main/clm_driver.F90 b/src/main/clm_driver.F90 index 33e9412ba9..6a2bb5ee20 100644 --- a/src/main/clm_driver.F90 +++ b/src/main/clm_driver.F90 @@ -1144,7 +1144,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro atm2lnd_inst, soilstate_inst, temperature_inst, active_layer_inst, & water_inst%waterstatebulk_inst, water_inst%waterdiagnosticbulk_inst, & water_inst%wateratm2lndbulk_inst, canopystate_inst, soilbiogeochem_carbonflux_inst, & - frictionvel_inst, soil_water_retention_curve) + frictionvel_inst, soil_water_retention_curve, fireemis_inst) ! TODO(wjs, 2016-04-01) I think this setFilters call should be replaced by a ! call to reweight_wrapup, if it's needed at all. diff --git a/src/main/clm_initializeMod.F90 b/src/main/clm_initializeMod.F90 index a53f6f2bdc..017639e76c 100644 --- a/src/main/clm_initializeMod.F90 +++ b/src/main/clm_initializeMod.F90 @@ -751,7 +751,7 @@ subroutine initialize2(ni,nj) call clm_fates%init_coldstart(water_inst%waterstatebulk_inst, & water_inst%waterdiagnosticbulk_inst, canopystate_inst, & - soilstate_inst, soilbiogeochem_carbonflux_inst) + soilstate_inst, soilbiogeochem_carbonflux_inst, fireemis_inst) end if ! topo_glc_mec was allocated in initialize1, but needed to be kept around through diff --git a/src/main/clm_instMod.F90 b/src/main/clm_instMod.F90 index 1ca450b48d..8aa4d4a59f 100644 --- a/src/main/clm_instMod.F90 +++ b/src/main/clm_instMod.F90 @@ -592,7 +592,9 @@ subroutine clm_instRest(bounds, ncid, flag, writing_finidat_interp_dest_file) soilstate_inst=soilstate_inst, & active_layer_inst=active_layer_inst, & soilbiogeochem_carbonflux_inst=soilbiogeochem_carbonflux_inst, & - soilbiogeochem_nitrogenflux_inst=soilbiogeochem_nitrogenflux_inst) + soilbiogeochem_nitrogenflux_inst=soilbiogeochem_nitrogenflux_inst,& + fireemis_inst=fireemis_inst& + ) end if diff --git a/src/utils/clmfates_interfaceMod.F90 b/src/utils/clmfates_interfaceMod.F90 index 1f856859e2..d210a1da87 100644 --- a/src/utils/clmfates_interfaceMod.F90 +++ b/src/utils/clmfates_interfaceMod.F90 @@ -43,6 +43,7 @@ module CLMFatesInterfaceMod use TemperatureType , only : temperature_type use EnergyFluxType , only : energyflux_type use SoilStateType , only : soilstate_type + use CNFireEmissionsMod, only : fireemis_type use CNProductsMod , only : cn_products_type use clm_varctl , only : iulog use clm_varctl , only : fates_parteh_mode @@ -112,7 +113,9 @@ module CLMFatesInterfaceMod use clm_varcon , only : dzsoi_decomp use FuncPedotransferMod, only: get_ipedof use CLMFatesParamInterfaceMod, only: fates_param_reader_ctsm_impl -! use SoilWaterPlantSinkMod, only : Compute_EffecRootFrac_And_VertTranSink_Default + use EDParamsMod, only : num_emission_compounds ! is this in the right place? + + ! use SoilWaterPlantSinkMod, only : Compute_EffecRootFrac_And_VertTranSink_Default ! Used FATES Modules use FatesInterfaceMod , only : fates_interface_type @@ -871,7 +874,7 @@ subroutine dynamics_driv(this, nc, bounds_clump, & atm2lnd_inst, soilstate_inst, temperature_inst, active_layer_inst, & waterstatebulk_inst, waterdiagnosticbulk_inst, wateratm2lndbulk_inst, & canopystate_inst, soilbiogeochem_carbonflux_inst, frictionvel_inst, & - soil_water_retention_curve) + soil_water_retention_curve, fireemis_inst) ! This wrapper is called daily from clm_driver ! This wrapper calls ed_driver, which is the daily dynamics component of FATES @@ -895,6 +898,7 @@ subroutine dynamics_driv(this, nc, bounds_clump, & type(waterdiagnosticbulk_type) , intent(inout) :: waterdiagnosticbulk_inst type(wateratm2lndbulk_type) , intent(inout) :: wateratm2lndbulk_inst type(canopystate_type) , intent(inout) :: canopystate_inst + type(fireemis_type), intent(inout) :: fireemis_inst type(soilbiogeochem_carbonflux_type), intent(inout) :: soilbiogeochem_carbonflux_inst type(frictionvel_type) , intent(inout) :: frictionvel_inst class(soil_water_retention_curve_type), intent(in) :: soil_water_retention_curve @@ -1138,6 +1142,7 @@ subroutine dynamics_driv(this, nc, bounds_clump, & waterdiagnosticbulk_inst, & canopystate_inst, & soilbiogeochem_carbonflux_inst, & + fireemis_inst, & .false.) ! --------------------------------------------------------------------------------- @@ -1314,7 +1319,7 @@ end subroutine UpdateCLitterFluxes subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump, & waterdiagnosticbulk_inst, canopystate_inst, & - soilbiogeochem_carbonflux_inst, is_initing_from_restart) + soilbiogeochem_carbonflux_inst, fireemis_inst, is_initing_from_restart) ! --------------------------------------------------------------------------------- ! This routine handles the updating of vegetation canopy diagnostics, (such as lai) @@ -1328,7 +1333,7 @@ subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump, & type(waterdiagnosticbulk_type) , intent(inout) :: waterdiagnosticbulk_inst type(canopystate_type) , intent(inout) :: canopystate_inst type(soilbiogeochem_carbonflux_type), intent(inout) :: soilbiogeochem_carbonflux_inst - + type(fireemis_type), intent(inout) :: fireemis_inst ! is this being called during a read from restart sequence (if so then use the restarted fates ! snow depth variable rather than the CLM variable). @@ -1340,7 +1345,8 @@ subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump, & integer :: s ! site index integer :: c ! column index integer :: g ! grid cell - + integer :: e ! emission compound + logical :: dispersal_flag ! local flag to pass to the inside of the site loop real(r8) :: areacheck call t_startf('fates_wrap_update_hlmfates_dyn') @@ -1541,10 +1547,12 @@ subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump, & displa(p) = this%fates(nc)%bc_out(s)%displa_pa(ifp) dleaf_patch(p) = this%fates(nc)%bc_out(s)%dleaf_pa(ifp) - if(this%fates(nc)%bc_out(s)%fire_emissions_pa(ifp,1).gt.0_r8)then - write(*,*) 'positive fire emissions',this%fates(nc)%bc_out(s)%fire_emissions_pa(ifp,1) - endif - + ! Unpack FATES fire emissions into CTSM side array. + fireemis_inst%fates_fire_emission_height_patch(p)=this%fates(nc)%bc_out(s)%fire_emission_height_pa(ifp) + do e = 1,num_emission_compounds + fireemis_inst%fates_fire_emissions_patch(p,e) = this%fates(nc)%bc_out(s)%fire_emissions_pa(p,e) + end do + end do ! veg pach if(abs(areacheck - 1.0_r8).gt.1.e-9_r8)then @@ -1564,7 +1572,7 @@ end subroutine wrap_update_hlmfates_dyn subroutine restart( this, bounds_proc, ncid, flag, waterdiagnosticbulk_inst, & waterstatebulk_inst, canopystate_inst, soilstate_inst, & active_layer_inst, soilbiogeochem_carbonflux_inst, & - soilbiogeochem_nitrogenflux_inst) + soilbiogeochem_nitrogenflux_inst, fireemis_inst) ! --------------------------------------------------------------------------------- ! The ability to restart the model is handled through three different types of calls @@ -1598,9 +1606,11 @@ subroutine restart( this, bounds_proc, ncid, flag, waterdiagnosticbulk_inst, & type(waterstatebulk_type) , intent(inout) :: waterstatebulk_inst type(canopystate_type) , intent(inout) :: canopystate_inst type(soilstate_type) , intent(inout) :: soilstate_inst + type(active_layer_type) , intent(in) :: active_layer_inst type(soilbiogeochem_carbonflux_type), intent(inout) :: soilbiogeochem_carbonflux_inst type(soilbiogeochem_nitrogenflux_type), intent(inout) :: soilbiogeochem_nitrogenflux_inst + type(fireemis_type), intent(inout) :: fireemis_inst ! Locals type(bounds_type) :: bounds_clump @@ -1876,7 +1886,7 @@ subroutine restart( this, bounds_proc, ncid, flag, waterdiagnosticbulk_inst, & ! ------------------------------------------------------------------------ call this%wrap_update_hlmfates_dyn(nc,bounds_clump, & waterdiagnosticbulk_inst,canopystate_inst, & - soilbiogeochem_carbonflux_inst, .true.) + soilbiogeochem_carbonflux_inst, fireemis_inst, .true.) ! ------------------------------------------------------------------------ ! Update the 3D patch level radiation absorption fractions @@ -1918,7 +1928,7 @@ end subroutine restart !===================================================================================== subroutine init_coldstart(this, waterstatebulk_inst, waterdiagnosticbulk_inst, & - canopystate_inst, soilstate_inst, soilbiogeochem_carbonflux_inst) + canopystate_inst, soilstate_inst, soilbiogeochem_carbonflux_inst, fireemis_inst) ! Arguments @@ -1928,7 +1938,8 @@ subroutine init_coldstart(this, waterstatebulk_inst, waterdiagnosticbulk_inst, & type(canopystate_type) , intent(inout) :: canopystate_inst type(soilstate_type) , intent(inout) :: soilstate_inst type(soilbiogeochem_carbonflux_type), intent(inout) :: soilbiogeochem_carbonflux_inst - + type(fireemis_type), intent(inout) :: fireemis_inst + ! locals integer :: nclumps integer :: nc @@ -2067,7 +2078,7 @@ subroutine init_coldstart(this, waterstatebulk_inst, waterdiagnosticbulk_inst, & ! ------------------------------------------------------------------------ call this%wrap_update_hlmfates_dyn(nc,bounds_clump, & waterdiagnosticbulk_inst,canopystate_inst, & - soilbiogeochem_carbonflux_inst, .false.) + soilbiogeochem_carbonflux_inst, fireemis_inst,.false.) ! ------------------------------------------------------------------------ ! Update history IO fields that depend on ecosystem dynamics