!### macros ##################################################### #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" !################################################################# module emission_data !----------------------------------------------------------------------- ! purpose ! ------- ! stores general data needed for emissions ! !----------------------------------------------------------------------- use Dims, only: nregions use global_types,only: emis_data #ifdef with_budgets use budget_global, only: apply_budget_global #endif implicit none private ! public variables: #ifdef with_budgets public :: budemi_dat, budemi_data #endif public :: plandr, sum_emission, emis2D public :: lm_bb, lm_bl, transcomdir, nsnap, day_of_week, cdate ! public routines: public :: msg_emis, do_add_2d, do_add_3d, flux_to_gridbox_per_month, gridbox_per_month_to_flux, normalise public :: calc_day_of_week, chardate, increase_itau #ifdef with_budgets type budemi_data real,dimension(:,:,:,:),pointer :: emi ! im,jm,nbud_vg,ntracet end type budemi_data #endif type(emis_data),dimension(nregions),target :: plandr ! field to store emissions in 2D prior to vertical redistribution in 3D type(emis_data),dimension(nregions),target :: emis2D #ifdef with_budgets ! budget of global region: type(budemi_data),dimension(nregions),target :: budemi_dat ! to collect the budget #endif real, dimension(nregions) :: sum_emission type(emis_data),dimension(nregions),target :: emis_temp ! temp array to store emission integer, parameter :: lm_bb = 15 integer, parameter :: lm_bl = 8 ! EMEP stuff integer, parameter :: nsnap = 11 integer, dimension(366) :: day_of_week character(len=130) :: transcomdir character(len=10) :: cdate ! used for character date representation contains subroutine msg_emis(year_month,msg,msg_comp,msg_mass,summed_emissions) ! ! provide output to verify the amount of emissions entering the calculations ! !************************************************ use dims, only: idate implicit none integer,intent(in) :: year_month ! 0: year,1: monthly amount character(len=*),intent(in) :: msg,msg_comp real,intent(in) :: msg_mass,summed_emissions character(len=7),dimension(2), parameter:: ym=(/'Yearly ','Monthly'/) 1111 format('Emis: ',a7,' Emission category ',a20,1x,a3,' [kg], mw:',f5.1,' month ',i2,' Sum',1x,1pe10.4) write(6,1111) ym(year_month),msg,msg_comp,msg_mass,idate(2),summed_emissions end subroutine msg_emis !-------------------------------------------------------------------------------- subroutine do_add_2d(region,i_tracer,l_level,emis_field,mol_mass,mol_mass_e, xfrac) !general routine to add a 2D emission field (given in kg/month) !of tracer i_tracer to level l_level !updtae also the budget! use toolbox, only : escape_tm, dumpfield use dims, only : CheckShape, isr, ier, jsr, jer, tref, nsrce use dims, only : sec_month, adv_scheme, im,jm, okdebug use global_data, only : mass_dat, region_dat use chem_param, only : ntracet, names #ifdef with_budgets use budget_global, only : budg_dat, nzon_vg use budget_fine, only : nzon, nzon_v, budemi #endif use ParTools, only: tracer_loc implicit none ! IO integer,intent(in) :: region integer,intent(in) :: i_tracer integer,intent(in) :: l_level real,dimension(:,:),intent(in) :: emis_field real,intent(in) :: mol_mass,mol_mass_e !first is the mass of tracer field real,intent(in), optional :: xfrac !partial addition !mol_mass_e is mass of emission field !e.g. NOx emission ar in kgN/month ! local integer :: i,j,nzone,nzone_v,itracer_loc real :: x, efrac, dtime, month2dt real,dimension(:,:,:,:),pointer :: rm,rzm integer,dimension(:,:), pointer :: zoomed real :: mbef, maft, sume zoomed => region_dat(region)%zoomed rm => mass_dat(region)%rm_t rzm => mass_dat(region)%rzm_t ! --- begin -------------------------------------- ! ! check arguments call CheckShape( (/im(region),jm(region)/), shape(emis_field) ) 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 itracer_loc = tracer_loc(i_tracer) if(itracer_loc < 0) call escape_tm('do_add_2d: itracer_loc < 0 ') if (okdebug) then print *, '--------------------------------------------------------------' print *, 'DO ADD 2D debug in region ', region print *, 'itracer ', i_tracer, names(i_tracer) print *, 'xfrac ', efrac print *, 'mol_mass ', mol_mass print *, 'mol_mass_e', mol_mass_e sume = sum(emis_field)*month2dt*(mol_mass/mol_mass_e)*efrac print *, 'sum emmis ', sume mbef = sum(rm(:,:,:,itracer_loc)) endif do j=jsr(region),jer(region) do i=isr(region),ier(region) if(zoomed(i,j)/=region) cycle x = emis_field(i,j)*month2dt*(mol_mass/mol_mass_e)*efrac rm(i,j,l_level,itracer_loc) = rm(i,j,l_level,itracer_loc) + x if (adv_scheme=='slope'.and.i_tracer.le.ntracet) & rzm(i,j,l_level,itracer_loc)=rzm(i,j,l_level,itracer_loc)-x !mk2fd: only in level1???? #ifdef with_budget ! budget___ if (apply_budget_global) then if(region==nregions) then !sub budget finest region nzone=nzon(i,j) nzone_v=nzon_v(l_level) budemi(nzone,nzone_v,i_tracer)=budemi(nzone,nzone_v,i_tracer)+x/mol_mass*1e3 !mole endif nzone=budg_dat(region)%nzong(i,j) !global budget nzone_v=nzon_vg(l_level) budemi_dat(region)%emi(i,j,nzone_v,i_tracer) = & budemi_dat(region)%emi(i,j,nzone_v,i_tracer) + x/mol_mass*1e3 endif #endif if(i_tracer ==1) sum_emission(region) = sum_emission(region) + x !in kg enddo enddo if (okdebug) then maft = sum(rm(:,:,:,itracer_loc)) print *, 'maft-mbef ', maft-mbef print *, 'END debug output' if ( sume > tiny(sume) ) then if ( (maft-mbef)/sume < 0.5) then call dumpfield(1, 'spurious.hdf', im(region), jm(region), 1, 1, emis_field, names(i_tracer)) endif endif endif nullify(zoomed) nullify(rm) nullify(rzm) end subroutine do_add_2d !--------------------------------------------------------------------------- subroutine do_add_3d(region,i_tracer,emis_field,mol_mass,mol_mass_e, xfrac) !general routine to add a 3D emission field (given in kg/month) !of tracer i_tracer use toolbox, only : escape_tm, dumpfield use dims, only : CheckShape, isr, ier, jsr, jer, tref, nsrce use dims, only : sec_month, adv_scheme, im,jm,lm, okdebug use global_data, only : mass_dat, region_dat use chem_param, only : ntracet, names #ifdef with_budget use budget_global, only : budg_dat, nzon_vg use budget_fine, only : nzon, nzon_v, budemi #endif use ParTools, only: tracer_loc implicit none ! IO integer,intent(in) :: region integer,intent(in) :: i_tracer real,dimension(:,:,:),intent(in) :: emis_field real,intent(in) :: mol_mass,mol_mass_e !first is the mass of tracer field real,intent(in), optional :: xfrac !partial addition !mol_mass_e is mass of emission field !e.g. NOx emission ar in kgN/month ! local integer :: i,j,l,nzone,nzone_v,itracer_loc integer, dimension(3) :: ubound_e real :: x, efrac, dtime, month2dt real,dimension(:,:,:,:),pointer :: rm integer,dimension(:,:), pointer :: zoomed real :: mbef, maft, sume zoomed => region_dat(region)%zoomed rm => mass_dat(region)%rm_t ! --- begin -------------------------------------- ! ! check arguments ubound_e = ubound(emis_field) if(ubound_e(1) /= im(region)) & call escape_tm('do_add_3d: im emission not consistent') if(ubound_e(2) /= jm(region)) & call escape_tm('do_add_3d: jm emission not consistent') if(ubound_e(3) > lm(region)) & call escape_tm('do_add_3d: lm emission not consistent') dtime=float(nsrce)/(2*tref(region)) !timestep emissions month2dt=dtime/sec_month ! conversion to emission per timestep itracer_loc = tracer_loc(i_tracer) if(itracer_loc < 0) call escape_tm('do_add_3d: itracer_loc < 0 ') if(present(xfrac)) then efrac = xfrac else efrac = 1.0 endif if (okdebug) then print *, '--------------------------------------------------------------' print *, 'DO ADD 3D debug in region ', region print *, 'itracer ', i_tracer, names(i_tracer) print *, 'xfrac ', efrac print *, 'mol_mass ', mol_mass print *, 'mol_mass_e', mol_mass_e sume = sum(emis_field)*month2dt*(mol_mass/mol_mass_e)*efrac print *, 'sum emmis ', sume print *, 'ubound(3) ', ubound_e(3) mbef = sum(rm(:,:,:,itracer_loc)) endif do l=1,ubound_e(3) do j=jsr(region),jer(region) do i=isr(region),ier(region) if(zoomed(i,j)/=region) cycle x = emis_field(i,j,l)*month2dt*(mol_mass/mol_mass_e)*efrac rm(i,j,l,itracer_loc) = rm(i,j,l,itracer_loc) + x #ifdef with_budget ! budget___ if(region==nregions) then !sub budget finest region nzone=nzon(i,j) nzone_v=nzon_v(l) budemi(nzone,nzone_v,i_tracer)=budemi(nzone,nzone_v,i_tracer)+x/mol_mass*1e3 !mole endif nzone=budg_dat(region)%nzong(i,j) !global budget nzone_v=nzon_vg(l) budemi_dat(region)%emi(i,j,nzone_v,i_tracer) = & budemi_dat(region)%emi(i,j,nzone_v,i_tracer) + x/mol_mass*1e3 #endif if(i_tracer ==1) sum_emission(region) = sum_emission(region) + x !in kg enddo enddo enddo !levels if (okdebug) then maft = sum(rm(:,:,:,itracer_loc)) print *, 'maft-mbef ', maft-mbef print *, 'END debug output' if ( sume > tiny(sume) ) then if ( (maft-mbef)/sume < 0.5) then call dumpfield(1, 'spurious.hdf', im(region), jm(region), ubound_e(3), 1, emis_field, names(i_tracer)) endif endif endif nullify(rm) nullify(zoomed) end subroutine do_add_3d subroutine flux_to_gridbox_per_month(field) use dims, only : nlon360, nlat180, dxy11, sec_month implicit none real, dimension(nlon360,nlat180) :: field integer :: i,j do j=1,nlat180 do i=1,nlon360 field(i,j)=field(i,j)*dxy11(j)*sec_month enddo enddo end subroutine flux_to_gridbox_per_month subroutine gridbox_per_month_to_flux(field) use dims, only : nlon360, nlat180, dxy11, sec_month implicit none real, dimension(nlon360,nlat180) :: field integer :: i,j do j=1,nlat180 do i=1,nlon360 field(i,j)=field(i,j)/dxy11(j)/sec_month enddo enddo end subroutine gridbox_per_month_to_flux subroutine normalise(tfd) use dims, only: nregions, im, tref, nsrce, dx, xref, xbeg use toolbox, only: escape_tm implicit none real, dimension(nsnap,24), intent(inout) :: tfd integer :: i, region, localt, isnap, nt real :: dtime, gmt, gmtshift, xlon, dxx real,dimension(:,:),allocatable :: fc do region = 1, nregions allocate(fc(nsnap,im(region))) fc = 0.0 dxx=dx/xref(region) !gridsize at current region dtime=float(nsrce)/(2*tref(region)) !timestep emissions gmt = dtime/2./3600.0 nt = 0 gmtloop: do nt = nt + 1 do i=1, im(region) xlon=xbeg(region)+(i-0.5)*dxx ! center of current longitude gmtshift = 24.0*xlon/360. localt = nint(gmt + gmtshift) localt = mod(24+localt,24) + 1 do isnap = 1, nsnap fc(isnap,i) = fc(isnap,i) + tfd(isnap,localt) enddo enddo gmt = gmt + dtime/3600. if (gmt > 24.0) exit gmtloop enddo gmtloop print *, 'test for the scaling factors:', region fc = abs(fc / nt - 1.0) if (maxval(fc) > 0.005) then do i=1,im(region) print '(i4, 11f10.6)' , i, fc(:,i) enddo call escape_tm(' normalise: values differ too much from 1.0' ) endif deallocate(fc) enddo end subroutine normalise subroutine calc_day_of_week(dof, year) use toolbox, only: escape_tm implicit none integer, intent(out), dimension(366) :: dof integer, intent(in) :: year integer :: i,is select case(year) case(1999) is = 4 case(2000) is = 5 case(2001) is = 0 case(2002) is = 1 case(2003) is = 2 case(2004) is = 3 case(2005) is = 5 case default call escape_tm('calc_day_of_year: invalid year') end select do i=1,366 dof(i) = is is = is + 1 if (is == 7) is = 0 enddo end subroutine calc_day_of_week subroutine chardate(chdate,idate6,cut) !----------------------------------------------------------------------- ! ! put date/time in character string ! ! on input: idate6 contains date/time ! !----------------------------------------------------------------------- implicit none integer,dimension(6),intent(in) :: idate6 ! "chardate" input date integer :: cut character(len=*) :: chdate ! character(len=2),dimension(12),parameter :: mon= & (/'01','02','03','04','05','06', & '07','08','09','10','11','12'/) ! ! put date/time in string ! if(cut.gt.8) cut=8 if(cut.eq.4) write(chdate,'(i4,a4)')idate6(1),' ' if(cut.eq.6) write(chdate,'(i4,a2,a2)')idate6(1),mon(idate6(2)),' ' if(cut.eq.8.and.idate6(3).lt.10) write(chdate,'(i4,a2,a1,i1)')idate6(1),mon(idate6(2)),'0',idate6(3) if(cut.eq.8.and.idate6(3).ge.10) write(chdate,'(i4,a2,i2)')idate6(1),mon(idate6(2)),idate6(3) ! end subroutine chardate subroutine increase_itau(itau,inc) !----------------------------------------------------------------------- implicit none integer :: itau,inc itau=itau+inc return end subroutine increase_itau end module emission_data