! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" ! !----------------------------------------------------------------------------- ! TM5 ! !----------------------------------------------------------------------------- !BOP ! ! !MODULE: EMISSION_DATA ! ! !DESCRIPTION: Provides general methods and fields needed for the emission ! routines. ! ! emission_vdist_by_sector : distribute vertically emissions according to source sector ! msg_emis : print formated emission total ! do_add_2d : add 2D emission to tracer mass array ! do_add_3d : add 3D emission to tracer mass array !\\ !\\ ! !INTERFACE: ! MODULE EMISSION_DATA ! ! !USES: ! use GO, only : gol, goPr, goErr use Dims, only : nregions, lm use global_types, only : emis_data, d2_data, d23_data, d3_data #ifdef with_budgets use budget_global,only : budg_dat, nzon_vg #endif IMPLICIT NONE ! ! !PUBLIC DATA MEMBERS & TYPES: ! public :: emis_data, d2_data, d23_data ! type(emis_data), target :: plandr(nregions) ! Land fraction type(emis_data), target :: emis2D(nregions) ! 2D emiss prior vertical redistribution type(emis_data), target :: emis_temp(nregions) ! temp array to store emission ! CMIP6 - AR5 logical :: LCMIP6 logical :: LAR5 character(len=256) :: emis_input_dir character(len=256) :: emis_input_dir_cmip6 character(len=256) :: emis_input_dir_ar5 integer :: emis_input_year_co2 ! biomass burning levels: integer, parameter :: bb_nlev = 5 integer, parameter :: bb_lm = lm(1) ! to be reduced to strip stratosphere! integer, parameter :: bl_lm = lm(1) ! to be reduced to the BL (around 8) real, parameter :: bb_hlow (0:bb_nlev) = (/ 0., 100., 500., 1000., 2000., 3000./) real, parameter :: bb_hhigh(0:bb_nlev) = (/ 100., 500., 1000., 2000., 3000., 6000./) ! biomass burning levels in tropics (-20 - 20):: integer, parameter :: bb_nlev_trop = 3 real, parameter :: bb_hlow_trop (0:bb_nlev_trop) = (/ 0., 100., 500., 1000./) real, parameter :: bb_hhigh_trop(0:bb_nlev_trop) = (/ 100., 500., 1000., 2000./) ! ----------------------------------------------------------------------------------- ! Budgets #ifdef with_budgets real, dimension(nregions) :: sum_emission type budemi_data real,dimension(:,:,:,:),pointer :: emi ! [i, j, nbud_vg, ntracet] end type budemi_data type(budemi_data),dimension(nregions),target :: budemi_dat #endif ! ---------------------------------------------------- ! Vertical splitting : used now for all sectors ! (tables in Dentener, 2006, ACP) ! ! ATTENTION: changes here have impact on the routine emission_vdist_by_sector (contained) ! Note: vdist_* variable could be hold private ! integer, parameter :: vd_class_name_len = 12 ! biomassburning integer, parameter :: vdist_bb_nlev = 6 real,dimension(vdist_bb_nlev),parameter :: vdist_bb_hlow = (/ 0., 100., 500., 1000., 2000., 3000./) real,dimension(vdist_bb_nlev),parameter :: vdist_bb_hhigh = (/ 100., 500., 1000., 2000., 3000., 6000./) ! other sources (related to surface) integer, parameter :: vdist_nn_nlev = 6 real,dimension(vdist_nn_nlev),parameter :: vdist_nn_hlow = (/ 0., 30. , 100., 300., 600., 1000. /) real,dimension(vdist_nn_nlev),parameter :: vdist_nn_hhigh = (/ 30., 100., 300., 600., 1000., 2000. /) ! ! !PRIVATE TYPES: ! character(len=*), parameter, private :: mname = 'emission_data' ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - adapted for AR5 ! 28 Jun 2012 - Ph. Le Sager - made everything public by default ! - adapted for lon-lat MPI domain decomposition ! !REMARKS: ! !EOP !------------------------------------------------------------------------ CONTAINS !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: EMISSION_VDIST_BY_SECTOR ! ! !DESCRIPTION: Return vertically distributed emissions depending on given ! sector and constituent. ! Distribution is done as 'compromise' between ! 1) Dentener et al., 2006, ACP, slightly modified for ! emissions supposed to enter the surface layer. ! 2) De Meij et al., ACP, 2006 ! 3) EMEP model (docu, model code, publications) ! 4) Bieser et al., GMDD, 2010 ! 5) Small updates to biomassburning following ! Huijnen et al., GMDD, 2010 !\\ !\\ ! !INTERFACE: ! SUBROUTINE EMISSION_VDIST_BY_SECTOR( splitname, constituent, region, emis2D, emis3D, status ) ! ! !USES: ! use toolbox, only : distribute_emis2D ! ! !INPUT PARAMETERS: ! character(len=vd_class_name_len), intent(in) :: splitname character(len=*), intent(in) :: constituent integer, intent(in) :: region type(emis_data), intent(in) :: emis2D ! ! !INPUT/OUTPUT PARAMETERS: ! type(d3_data), intent(inout) :: emis3D ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - v0 ! ! !REMARKS: Splitting depending on constituent is still to code. ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/emission_vdist_by_sector' integer :: ilev character(len=2) :: splittype integer, parameter :: maxvdist_nlev = max(vdist_bb_nlev,vdist_nn_nlev) real, dimension(maxvdist_nlev) :: sfglobal, sfintrop, sfintemp ! --- begin -------------------------------------- ! control implicit settings here in this routine if( vdist_bb_nlev /= 6 .or. vdist_nn_nlev /= 6 ) then write(gol, '("ERROR: wrong number of layers: bb_nlev /= 6 or nn_lev /= 6")'); call goErr status=1; return end if ! initialise sfglobal = 0.0 sfintrop = 0.0 sfintemp = 0.0 ! ------------------------------------------------- ! This here is selection by emission source sectors ! ------------------------------------------------- select case( trim(splitname) ) case( 'combenergy' ) ! Combustion from energy production and transformation, power-plants ! no zonal differences, no dependency on species ! -- ! assumed to be related to EMEP SNAP category [1] ! AEROCOM equivalent: power-plants ! -- splittype = 'nn' sfglobal = (/ 0.0, 0.1, 0.7, 0.2, 0.0, 0.0 /) case( 'combrescom' ) ! Residential and commercial combustion ! no zonal differences, no dependency on species ! -- ! assumed to be related to EMEP SNAP category [2] ! AEROCOM equivalent: domestic/industry ! -- splittype = 'nn' sfglobal = (/ 0.4, 0.4, 0.2, 0.0, 0.0, 0.0 /) case( 'industry' ) ! Industrial processes and combustion ! no zonal differences, no dependency on species ! -- ! assumed to be an aggregate of EMEP SNAP categories [3,4,5] ! AEROCOM equivalent: industry ! -- splittype = 'nn' sfglobal = (/ 0.1, 0.2, 0.6, 0.1, 0.0, 0.0 /) case( 'waste' ) ! Waste treatment and disposal ! no zonal differences, no dependency on species ! -- ! assumed to be related to EMEP SNAP category [9] ! AEROCOM equivalent: - ! -- splittype = 'nn' sfglobal = (/ 0.1, 0.2, 0.4, 0.3, 0.0, 0.0 /) case( 'nearsurface' ) ! Near surface emissions ! no zonal differences, no dependency on species ! -- ! AR5: solvent production and use, agricultural waste burning, & ! ships, grassland fire, mineral dust (AEROCOM or Tegen-Vignati) ! EMEP equivalents: SNAP categories [6,8] ! AEROCOM equivalents: ships, agricultural waste ! -- splittype = 'nn' sfglobal = (/ 0.8, 0.2, 0.0, 0.0, 0.0, 0.0 /) case( 'surface' ) ! Surface emissions ! no zonal differences, no dependency on species ! -- ! AR5: agriculture, transport, soil, oceanic, biogenic ! EMEP equivalents: SNAP categories [7,8,10] ! AEROCOM equivalents: surface emissions from road/off-road, ships, agricultural waste ! -- splittype = 'nn' sfglobal = (/ 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 /) case( 'forestfire' ) ! forest fires: 6 layer model, different splitting for different regions ! dont distinguish between boreal Europe and boreal Canada, use Europe globally ! Dentener et al., modified in tropics (up to 2 km) with respect to Huijnen et al., GMD 2010, \ ! who are citing Labonne et al., 2007 ofr new information based on satellite data. splittype = 'bb' sfglobal = (/ 0.1, 0.1, 0.2, 0.2, 0.4, 0.0 /) sfintemp = (/ 0.2, 0.2, 0.2, 0.4, 0.0, 0.0 /) sfintrop = (/ 0.2, 0.2, 0.2, 0.4, 0.0, 0.0 /) case default write(gol, '("ERROR: wrong specification of emission sector ",a,"in vdist_by_sector called for constituent ",a)') & trim(splitname), trim(constituent); call goErr status=1; return end select ! -------------------- ! do the splitting ! -------------------- select case( splittype ) case( 'nn' ) do ilev = 1, vdist_nn_nlev call distribute_emis2D( region, emis2D, emis3D, vdist_nn_hlow(ilev), vdist_nn_hhigh(ilev), status, sfglobal(ilev) ) IF_NOTOK_RETURN(status=1) end do case( 'bb' ) do ilev = 1, vdist_bb_nlev ! boreal call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, sfglobal(ilev), -90.,-61.) IF_NOTOK_RETURN(status=1) call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, sfglobal(ilev), 60., 90.) IF_NOTOK_RETURN(status=1) ! temperate call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, sfintemp(ilev), -60.,-31.) IF_NOTOK_RETURN(status=1) call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, sfintemp(ilev), 30., 59.) IF_NOTOK_RETURN(status=1) ! tropics call distribute_emis2D( region, emis2D, emis3D, vdist_bb_hlow(ilev), vdist_bb_hhigh(ilev), status, sfintrop(ilev), -30., 29.) IF_NOTOK_RETURN(status=1) end do case default write(gol, '("ERROR: wrong specification of splitting type in vdist_by_sector.")'); call goErr status=1; return end select ! OK status=0 END SUBROUTINE EMISSION_VDIST_BY_SECTOR !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: MSG_EMIS ! ! !DESCRIPTION: Provide output to verify the amount of emissions ! entering the calculations !\\ !\\ ! !INTERFACE: ! SUBROUTINE MSG_EMIS(year_month, provider, sector, msg_comp, msg_mass, summed_emissions) ! ! !USES: ! use dims, only: idate ! ! !INPUT PARAMETERS: ! integer, intent(in) :: year_month ! 1: year, 2: monthly amount character(len=*), intent(in) :: provider, sector, msg_comp real, intent(in) :: msg_mass, summed_emissions ! ! !REVISION HISTORY: ! 1 Oct 2010 - Achim Strunk - protex; increase space for output ! !EOP !------------------------------------------------------------------------ !BOC character(len=7),dimension(2), parameter:: ym=(/'Yearly ','Monthly'/) 1111 format(' EMISS-INFO - ',a7,' Emission from ',a8,' -> ',a18,1x,a6,' : mw=',f5.1,', month=',i2.2,', Sum [kg]=',1x,1pe11.4) 1112 format(' EMISS-INFO - ',a7,' Emission from ',a8,' -> ',a18,1x,a6,' : mw=',f5.1,', --------, Sum [kg]=',1x,1pe11.4) if (year_month==1) then write (gol,1112) ym(year_month), provider, sector, msg_comp, msg_mass, summed_emissions; call goPr else write (gol,1111) ym(year_month), provider, sector, msg_comp, msg_mass, idate(2), summed_emissions; call goPr endif END SUBROUTINE MSG_EMIS !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: DO_ADD_2D ! ! !DESCRIPTION: General routine to add a 2D emission field (given in kg ! /month) into tracer mass array. The mass increase is done ! at level L_LEVEL, and for tracer I_TRACER. ! ! Update accordingly the budget. !\\ !\\ ! !INTERFACE: ! SUBROUTINE DO_ADD_2D( region, i_tracer, l_level, is, js, emis_field, & mol_mass, mol_mass_e, status, xfrac ) ! ! !USES: ! use dims, only : CheckShape, tref, nsrce use dims, only : sec_month, im, jm, okdebug use global_data, only : mass_dat, region_dat use chem_param, only : ntracet, names use tm5_distgrid, only : get_distgrid, dgrid #ifdef with_budgets use budget_global, only : budg_dat, nzon_vg #endif ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region integer, intent(in) :: i_tracer integer, intent(in) :: l_level integer, intent(in) :: is, js real, intent(in) :: emis_field(is:,js:) real, intent(in) :: mol_mass ! mol mass of tracer field real, intent(in) :: mol_mass_e ! mol mass of emission field (e.g. NOx emission ar in kgN/month) real, intent(in), optional :: xfrac ! partial addition ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/do_add_2d' integer :: i, j, nzone, nzone_v, i1, i2, j1, j2 real :: x, efrac, dtime, month2dt real :: mbef, maft, sume ! mass BEFore, AFTer, SUM Emiss (debug) real, dimension(:,:,:,:), pointer :: rm, rzm ! --- begin -------------------------------------- call GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) call CheckShape( (/i2-i1+1,j2-j1+1/), shape(emis_field), status ) IF_NOTOK_RETURN(status=1) dtime=float(nsrce)/(2*tref(region)) ! timestep emissions month2dt=dtime/sec_month ! conversion to emission per timestep rm => mass_dat(region)%rm rzm => mass_dat(region)%rzm if(present(xfrac)) then efrac = xfrac else efrac = 1.0 endif if (okdebug) then write (gol, '("--------------------------------------------------------------")') ; call goPr write (gol, '(" DO ADD 2D debug in region:",i3) ' ) region ; call goPr write (gol, '(" itracer :",i3,a8)' ) i_tracer, ' :'//names(i_tracer) ; call goPr write (gol, '(" xfrac :",f6.1)' ) efrac ; call goPr write (gol, '(" mol_mass :",f6.1 ) ' ) mol_mass ; call goPr write (gol, '(" mol_mass_e :",f6.1 )' ) mol_mass_e ; call goPr write (gol, '(" emisfield :",2e12.4 )') minval(emis_field),maxval(emis_field) ; call goPr sume = sum(emis_field)*month2dt*(mol_mass/mol_mass_e)*efrac write (gol, '(" sum emmis :",e12.4 )' ) sume ; call goPr mbef = sum(rm(i1:i2,j1:j2,:,i_tracer)) endif do j=j1,j2 do i=i1,i2 x = emis_field(i,j)*month2dt*(mol_mass/mol_mass_e)*efrac rm(i,j,l_level,i_tracer) = rm(i,j,l_level,i_tracer) + x #ifdef slopes ! injection of emissions at surface, thus vertical slope more negative ! (keep concentration at top of layer the same as before emission) if ( i_tracer <= ntracet ) then rzm(i,j,l_level,i_tracer) = rzm(i,j,l_level,i_tracer) - x end if #endif #ifdef with_budgets nzone = budg_dat(region)%nzong(i,j) !global budget nzone_v = nzon_vg(l_level) #ifndef without_emission budemi_dat(region)%emi(i,j,nzone_v,i_tracer) = & budemi_dat(region)%emi(i,j,nzone_v,i_tracer) + x/mol_mass*1e3 if(i_tracer == 1) sum_emission(region) = sum_emission(region) + x !in kg #endif #endif enddo enddo if (okdebug) then maft = sum(rm(i1:i2,j1:j2,:,i_tracer)) write (gol, '(" Added amount :",e12.4 )' ) maft-mbef ; call goPr write (gol, '(" Added amount can be different when inner zoom is present!" )' ) ; call goPr if (maft-mbef-sume > 1e-8*sume) then write (gol, '(" Warning: error in emission ")' ) ; call goErr endif endif nullify(rm) nullify(rzm) status=0 END SUBROUTINE DO_ADD_2D !EOC !-------------------------------------------------------------------------- ! TM5 ! !-------------------------------------------------------------------------- !BOP ! ! !IROUTINE: DO_ADD_3D ! ! !DESCRIPTION: General routine to add a 3D emission field (given in kg ! /month) into tracer mass array. The mass increase is done ! for tracer I_TRACER. ! ! Update accordingly the budget. !\\ !\\ ! !INTERFACE: ! SUBROUTINE DO_ADD_3D( region, i_tracer, is, js, emis_field, mol_mass, mol_mass_e, status, xfrac) ! ! !USES: ! use dims, only : CheckShape, nsrce, tref use dims, only : sec_month, im, jm, lm, okdebug use global_data, only : mass_dat, region_dat use chem_param, only : ntracet, names use tm5_distgrid, only : get_distgrid, dgrid #ifdef with_budgets use budget_global, only : budg_dat, nzon_vg #endif ! ! !INPUT PARAMETERS: ! integer, intent(in) :: region ! region # integer, intent(in) :: i_tracer ! id of tracer to increase integer, intent(in) :: is, js real, intent(in) :: emis_field(is:,js:,:) real, intent(in) :: mol_mass, mol_mass_e real, intent(in), optional :: xfrac !partial addition ! ! !OUTPUT PARAMETERS: ! integer, intent(out) :: status ! ! !REVISION HISTORY: ! 27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition ! !EOP !------------------------------------------------------------------------ !BOC character(len=*), parameter :: rname = mname//'/do_add_3d' integer :: i, j, l, nzone, nzone_v, i1,j1,i2,j2 integer, dimension(3) :: ubound_e real :: x, efrac, dtime, month2dt real, dimension(:,:,:,:), pointer :: rm real :: mbef, maft, sume ! debug ! --- begin -------------------------------------- call GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 ) status = 0 ! check arguments ubound_e = ubound(emis_field) if(ubound_e(1) /= i2) then write(gol,'(A)') 'do_add_3d: im emission not consistent' ; call goPr status = 1 endif if(ubound_e(2) /= j2) then write(gol,'(A)') 'do_add_3d: jm emission not consistent' ; call goPr status = 1 endif if(ubound_e(3) > lm(region)) then write(gol,'(A)') 'do_add_3d: lm emission not consistent' ; call goPr status = 1 endif IF_NOTOK_RETURN(status=1) rm => mass_dat(region)%rm dtime=float(nsrce)/(2*tref(region)) ! timestep emissions month2dt=dtime/sec_month ! conversion to emission per timestep if(present(xfrac)) then efrac = xfrac else efrac = 1.0 endif ! if (okdebug) then ! write (gol, '("--------------------------------------------------------------")') ; call goPr ! write (gol, '(" DO ADD 3D debug in region:",i3) ' ) region ; call goPr ! write (gol, '(" itracer :",i3,a8)' ) i_tracer, ' :'//names(i_tracer) ; call goPr ! write (gol, '(" xfrac :",f6.1)' ) efrac ; call goPr ! write (gol, '(" mol_mass :",f6.1 ) ' ) mol_mass ; call goPr ! write (gol, '(" mol_mass_e :",f6.1 )' ) mol_mass_e ; call goPr ! write (gol, '(" emisfield :",2e12.4 )') minval(emis_field), maxval(emis_field) ; call goPr ! write (gol, '(" ubound(3) :",i3 )' ) ubound_e(3) ; call goPr ! sume = sum(emis_field)*month2dt*(mol_mass/mol_mass_e)*efrac ! write (gol, '(" sume :", e12.4 )') sume ; call goPr ! mbef = sum(rm(i1:i2,j1:j2,:,i_tracer)) ! endif do l=1,ubound_e(3) do j=j1,j2 do i=i1,i2 x = emis_field(i,j,l)*month2dt*(mol_mass/mol_mass_e)*efrac rm(i,j,l,i_tracer) = rm(i,j,l,i_tracer) + x #ifdef with_budgets ! budget___ nzone=budg_dat(region)%nzong(i,j) !global budget nzone_v=nzon_vg(l) #ifndef without_emission budemi_dat(region)%emi(i,j,nzone_v,i_tracer) = & budemi_dat(region)%emi(i,j,nzone_v,i_tracer) + x/mol_mass*1e3 if(i_tracer ==1) sum_emission(region) = sum_emission(region) + x !in kg #endif #endif enddo enddo enddo !levels ! if (okdebug) then ! maft = sum(rm(i1:i2,j1:j2,:,i_tracer)) ! write (gol, '(" after-before :", e12.4 )' ) maft-mbef ; call goPr ! write (gol, '(" END debug output ")' ) ; call goPr ! endif nullify(rm) END SUBROUTINE DO_ADD_3D !EOC END MODULE EMISSION_DATA