!### 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  
