ecearth3/sources/tm5mp/proj/cbm4/emission.F90

#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') 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
!
! !DESCRIPTION: wrappers around various emissions (init/declare/apply/done)
!                routines, needed for TM5 CBM4 version.
!               Also hold emissions budget variables.
! 
!\\
!\\
! !INTERFACE: 
!
MODULE EMISSION
  !
  ! !USES:
  !
  USE GO,             ONLY : gol, goErr, goPr
  use GO,             ONLY : GO_Timer_Def, GO_Timer_End, GO_Timer_Start
  USE TM5_DISTGRID,   ONLY : dgrid, Get_DistGrid, scatter, gather
  USE dims,           ONLY : nregions, okdebug
  
  USE emission_data,  ONLY : plandr, emis2D ! , bmbcycle, bb_lm

#ifdef with_m7
  USE emission_data,  ONLY : emis_number, emis_mass, emis_temp
#endif

#ifdef with_budgets
  USE emission_data,  ONLY : budemi_dat, budemi_data, sum_emission
  USE budget_global,  ONLY : nbud_vg,nbudg
  USE chem_param,     ONLY : ntracet
#endif

  use emission_nox  , only : Emission_NOx_Init  , Emission_NOx_Done   , emission_nox_declare
  use emission_co   , only : Emission_CO_Init   , Emission_CO_Done    , emission_co_declare   , emission_co_apply
  use emission_nmvoc, only : Emission_NMVOC_Init, Emission_NMVOC_Done , emission_nmvoc_declare, emission_nmvoc_apply
  use emission_ch4 ,  only : Emission_CH4_Init  , Emission_CH4_Done   , emission_ch4_declare  , emission_ch4_apply
  use emission_nh3 ,  only : Emission_nh3_Init  , Emission_nh3_Done   , emission_nh3_declare  , emission_nh3_apply
  use emission_sox ,  only : Emission_SOx_Init  , Emission_SOx_Done   , emission_sox_declare  , emission_sox_apply
  use emission_dms ,  only : Emission_DMS_Init  , emission_dms_done   , emission_dms_declare  , emission_dms_apply
  use emission_rn222, only : Emission_rn222_Init, emission_rn222_done , emission_rn222_declare, emission_rn222_apply

#ifdef with_online_bvoc
  USE emission_bvoc,  ONLY : declare_emission_bvoc,  free_emission_bvoc,  emission_apply_bvoc
#else
  use emission_isop , only:  Emission_isop_Init , Emission_isop_Done , emission_isop_declare , emission_isop_apply
#endif
#ifdef with_m7 
  use emission_dust,  only :                    emission_dust_done , emission_dust_declare
  use emission_pom ,  only : Emission_POM_Init, Emission_POM_Done  , emission_pom_declare
  use emission_bc  ,  only : Emission_BC_Init , Emission_BC_Done   , emission_bc_declare
!!$  use emission_ss,    only: declare_emission_ss  , free_emission_ss   
#endif

  IMPLICIT NONE  
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  PUBLIC  :: Emission_Init     ! allocate/init budget var; call other emiss-related init
  PUBLIC  :: Emission_Done     ! gather/write final budget
  PUBLIC  :: Declare_Emission  ! allocate emiss var (new run), read emiss data (new month)
  PUBLIC  :: Emission_Apply    ! 
  !
  ! !PRIVATE DATA MEMBERS:
  !
  ! budemig_all (used to sum budemig from all processors) is used in chemistry for its NOx values.
#ifdef with_budgets
  REAL, DIMENSION(nbudg, nbud_vg, ntracet)         :: budemig
  REAL, DIMENSION(nbudg, nbud_vg, ntracet), PUBLIC :: budemig_all  ! for buggy MPI (see budget_global.F90)
#endif

  integer ::  itim_appl, itim_co, itim_voc, itim_nh3, itim_sox, itim_dms, itim_ch4, itim_isop, itim_rn222
  
  CHARACTER(len=*), PARAMETER  ::  mname = 'emission'  
  !
  ! !REVISION HISTORY: 
  !   16 Jul 2010 - P. Le Sager - fix for m7 with GFED_8days 
  !   20 Aug 2010 - A. Strunk   - Adapted to AR5 emissions + various other changes. 
  !   27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  !EOP
  !----------------------------------------------------------------------

CONTAINS

  !----------------------------------------------------------------------
  !                    TM5                                              !
  !----------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_INIT
  !
  ! !DESCRIPTION: Initialise emission fields and parameters by reading 
  !               the rc-file. Allocate and initialize budget variables.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_INIT( status )
    !
    ! !USES:
    !
    use GO,              only : TrcFile, Init, Done, ReadRc
    use meteodata,       only : Set,     oro_dat
    use dims,            only : iglbsfc, nregions, lm
    use dims,            only : idate,   ndyn_max, tref
    use global_data,     only : inputdir, rcfile
    use emission_data,   only : emis_input_dir,  use_tiedkte
    use emission_data,   only : emis_input_dir_gfed
    use emission_data,   only : emis_input_dir_retro
    use emission_data,   only : emis_input_year
    use emission_data,   only : LAR5, LAR5BMB, LEDGAR4,  LRETROF, LGFED3
    use emission_data,   only : LLPJ, LHYMN,   LMACCITY, LMEGAN,  LMACC
#ifdef with_ch4_emis
    use emission_data,   only : emis_input_dir_natch4
#endif
#ifdef with_m7
    use chem_param,      only : mode_nm
    use emission_data,   only : emis_input_dir_aerocom
    use emission_data,   only : emis_input_dir_dust, emis_input_dust
    use mo_aero_m7,      only : nmod
#endif
    use emission_data,   only : emis_input_dir_ar5
    use emission_data,   only : emis_input_dir_mac
    use emission_data,   only : emis_input_dir_ed4
    use emission_data,   only : emis_input_dir_dms
    use emission_data,   only : emis_input_dir_rn222
    use emission_data,   only : emis_input_dir_megan
    use emission_data,   only : emis_ch4_single, emis_ch4_fix3d
    use emission_data,   only : emis_ch4_fixed_ppb, emis_zch4_fname
    use emission_data,   only : emis_bb_trop_cycle, bb_cycle, scale_cycle
    use emission_read,   only : emission_read_init
#ifdef with_online_nox
    use online_nox_data, only : input_nox_dir
#endif
    !
    ! !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//'/Emission_Init'
    INTEGER                        ::  region, i1, i2, j1, j2, ntim, lmr, imode
    TYPE(TrcFile)                  ::  rcF
    REAL                           ::  dtime

    ! -----------------------------------
    ! read settings from rcfile
    ! -----------------------------------
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    if (okdebug) then
       write(gol,*) "EMISS-INFO - running year : ", idate(1) ; call goPr
    end if

    ! emission base year (assumption: no run overlapping more than one year)
    call ReadRc( rcF, 'input.emis.year', emis_input_year, status, default=idate(1) )
    IF_ERROR_RETURN(status=1)

    write(gol,*) 'EMISS-INFO - Emissions base year : ', emis_input_year; call goPr

    ! default directory for emissions data is "standard input files" dir
    emis_input_dir=trim(inputdir)

    ! directory of each data provider

    call ReadRc( rcF, 'input.emis.dir.AR5',   emis_input_dir_ar5,   status, default=emis_input_dir )
    IF_ERROR_RETURN(status=1)

    call ReadRc( rcF, 'input.emis.dir.MACC',  emis_input_dir_mac,   status, default=emis_input_dir )
    IF_ERROR_RETURN(status=1)

    call ReadRc( rcF, 'input.emis.dir.ED41',  emis_input_dir_ed4,   status, default=emis_input_dir )
    IF_ERROR_RETURN(status=1)

    call ReadRc( rcF, 'input.emis.dir.gfed',  emis_input_dir_gfed,  status, default=emis_input_dir )
    IF_ERROR_RETURN(status=1)

    call ReadRc( rcF, 'input.emis.dir.retro', emis_input_dir_retro, status, default=emis_input_dir )
    IF_ERROR_RETURN(status=1)

    call ReadRc( rcF, 'input.emis.dir.MEGAN', emis_input_dir_megan, status, default=emis_input_dir )
    IF_ERROR_RETURN(status=1)

#ifdef with_ch4_emis
    ! for both HYMN and LPJ datasets
    call ReadRc( rcF, 'input.emis.dir.natch4',  emis_input_dir_natch4, status, default=emis_input_dir )
    IF_ERROR_RETURN(status=1)
#endif

    ! Flags
    call ReadRc( rcF, 'use_ar5',         LAR5,    status, default=.false. )
    IF_ERROR_RETURN(status=1) 

    call ReadRc( rcF, 'use_ar5_fires',   LAR5BMB, status, default=.false. )
    IF_ERROR_RETURN(status=1) 

    call ReadRc( rcF, 'use_edgar4',      LEDGAR4, status, default=.false. )
    IF_ERROR_RETURN(status=1) 

    call ReadRc( rcF, 'use_retro_fires', LRETROF, status, default=.false. )
    IF_ERROR_RETURN(status=1) 

    call ReadRc( rcF, 'use_gfed3',       LGFED3,  status, default=.false. )
    IF_ERROR_RETURN(status=1) 

    call ReadRc( rcF, 'use_macc',        LMACC,   status, default=.false. )
    IF_ERROR_RETURN(status=1) 

    call ReadRc( rcF, 'use_lpj',         LLPJ,    status, default=.false. )
    IF_ERROR_RETURN(status=1) 

    call ReadRc( rcF, 'use_hymn',        LHYMN,   status, default=.false. )
    IF_ERROR_RETURN(status=1) 

    call ReadRc( rcF, 'use_maccity',     LMACCITY,status, default=.false. )
    IF_ERROR_RETURN(status=1)

    IF(LMACCITY) LMACC=.true.   ! ensure that MACC data are read. LMACCITY just add anthro sector to MACC's sector list.
    
    call ReadRc( rcF, 'use_megan',       LMEGAN,  status, default=.false. )
    IF_ERROR_RETURN(status=1) 

    ! very basic checks
    if (count((/ LAR5, LEDGAR4, LMACCITY /)) > 1) then
       write(gol,*) 'You use more than one ANTHROPOGENIC inventory'; call goErr
       status=1; TRACEBACK; return
    end if
    
    if (count((/ LAR5BMB, LRETROF, LGFED3 /)) > 1) then
       write(gol,*) 'You use more than one BIOMASS BURNING inventory'; call goErr
       status=1; TRACEBACK; return
    end if

    
    ! init providers info 
    call emission_read_init( rcF, status )

#ifdef with_online_nox
    call ReadRc( rcF, 'input.onlinenox.dir', input_nox_dir, status )
    IF_NOTOK_RETURN(status=1)
#endif

    ! are convection fluxes computed (Tiedkte) or read?
    call ReadRc( rcF, 'tiedtke', use_tiedkte, status )
    IF_NOTOK_RETURN(status=1)

    
#ifdef with_m7
    call ReadRc( rcF, 'input.emis.dir.aerocom', emis_input_dir_aerocom, status, default=emis_input_dir )
    IF_NOTOK_RETURN(status=1)

    call ReadRc( rcF, 'input.emis.dust', emis_input_dust, status, default="AEROCOM" )
    IF_NOTOK_RETURN(status=1) 

    call ReadRc( rcF, 'input.emis.dir.dust', emis_input_dir_dust, status, default=emis_input_dir )
    IF_NOTOK_RETURN(status=1) 
#endif
    call ReadRc( rcF, 'input.emis.dir.dms', emis_input_dir_dms, status, default=emis_input_dir )
    IF_NOTOK_RETURN(status=1)

    call ReadRc( rcF, 'input.emis.dir.rn222', emis_input_dir_rn222, status, default=emis_input_dir )
    IF_ERROR_RETURN(status=1)

    ! Get biomassburning time splitting factors (same for all constituents)
    ! -----------------------------------------------------------------------------------
    call ReadRc( rcF, 'input.emis.bb.dailycycle', emis_bb_trop_cycle, status, default=.false. )
    IF_ERROR_RETURN(status=1) 

    if (emis_bb_trop_cycle) then 
       do region = 1, nregions
          dtime = float(ndyn_max)/(2*tref(region))    ! timestep emissions  (CMK changed 5/2006)
          ntim  = 86400/nint(dtime)                   ! number of timesteps in 24 hours for this region
          allocate(bb_cycle(region)%scalef(ntim))
          call scale_cycle(ntim, bb_cycle(region)%scalef)
       end do
    end if
    
    ! CH4
    call ReadRc( rcF, 'input.emis.ch4.single', emis_ch4_single, status )
    IF_NOTOK_RETURN(status=1)

    if( emis_ch4_single ) then
       call ReadRc( rcF, 'input.emis.ch4.fixed_ppb', emis_ch4_fixed_ppb, status )
       IF_NOTOK_RETURN(status=1)
       call ReadRc( rcF, 'input.emis.ch4.fix3d', emis_ch4_fix3d, status )
       IF_NOTOK_RETURN(status=1)
       if( emis_ch4_fix3d ) then
          write (gol,*) 'EMISS-INFO - 3-D CH4 field will be fixed to ', emis_ch4_fixed_ppb, ' ppb'; call goPr
       else
          write (gol,*) 'EMISS-INFO - surface CH4 will be fixed to ', emis_ch4_fixed_ppb, ' ppb'; call goPr
       endif
    else
       ! root name
       call ReadRc( rcF, 'input.emis.ch4.surf', emis_zch4_fname, status )
       IF_NOTOK_RETURN(status=1)
       write (gol,*) 'EMISS-INFO - surface CH4 will be fixed by a zonal mean surface field'; call goPr
    endif ! ch4_single

    ! used by vertical distribution:
    CALL Set( oro_dat(iglbsfc), status, used=.TRUE. )


    ! Allocate data
    ! -------------
    DO region=1,nregions

       CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )          
       lmr = lm(region)

       ALLOCATE( plandr(region)%surf(i1:i2, j1:j2))
       ALLOCATE( emis2D(region)%surf(i1:i2, j1:j2))    

#ifdef with_m7
       ! aerosols (up to lmr instead of bb_lm)
       DO imode=1,nmod
          ALLOCATE(emis_number(region,imode)%d4( i1:i2, j1:j2, lmr, mode_nm(imode)))
          ALLOCATE(emis_mass  (region,imode)%d4( i1:i2, j1:j2, lmr, mode_nm(imode)))
       ENDDO
       ALLOCATE(emis_temp(region)%surf(i1:i2, j1:j2))
#endif

#ifdef with_budgets
       ALLOCATE( budemi_dat(region)%emi(i1:i2, j1:j2, nbud_vg, ntracet) )
       budemi_dat(region)%emi = 0.0
       sum_emission(region)   = 0.0
#endif       
    ENDDO
        
    ! Done
    ! -------------
    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)


    ! define timers:
    call GO_Timer_Def( itim_appl, 'emission appl', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_co, 'emission co', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_voc, 'emission voc', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_dms, 'emission dms', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_nh3, 'emission nh3', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_sox, 'emission sox', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_ch4, 'emission ch4 ', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_isop, 'emission isop', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_rn222, 'emission rn222', status )
    IF_NOTOK_RETURN(status=1)

    
    status = 0

  END SUBROUTINE EMISSION_INIT
  !EOC
  
  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_DONE
  !
  ! !DESCRIPTION: calculate and write final budgets
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_DONE( status )
    !
    ! !USES:
    !
    USE dims,          ONLY : nregions, im, jm

#ifdef with_budgets
    USE chem_param,    ONLY : ntracet, names
    USE budget_global, ONLY : budget_file_global, nbud_vg, budg_dat, nbudg, NHAB
    USE file_hdf,      ONLY : THdfFile, TSds
    USE file_hdf,      ONLY : Init, Done, WriteAttribute, WriteData, SetDim
    USE Dims,          ONLY : region_name
    USE partools,      ONLY : isRoot, par_reduce, par_reduce_element
#endif
    use emission_data, only : bb_cycle
    use emission_data, only : emis_bb_trop_cycle
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)           ::  status
    !
    ! !REVISION HISTORY: 
    !   16 Jul 2010 - A. Strunk - 
    !   27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER :: rname = mname//'/Emission_Done'
    INTEGER                     :: region, i1, i2, j1, j2

#ifdef with_budgets
    TYPE(THdfFile)                        :: io
    TYPE(TSds)                            :: sds
    REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: collect_emissions
    INTEGER                               :: nsend,j,i,n,nzone,nzone_v
    real, dimension(nregions) :: sum_emission_all
#endif

    ! --- begin ---------------------------------

#ifdef with_budgets

    ! add up contribution from all proc
    DO region = 1, nregions
       
       CALL PAR_REDUCE(sum_emission(region), 'SUM', sum_emission_all(region), status)
       IF_NOTOK_RETURN(status=1)
       
    END DO
    
    ! Write global budget of tracer #1 
    IF ( isRoot ) THEN 

       write (gol,'("EMISS-INFO - ----------------------------------------------")'); call goPr
       write (gol,'("EMISS-INFO -  Budget of tracer ",a," (kg) ")') trim(names(1))  ; call goPr
       write (gol,'("EMISS-INFO - ----------------------------------------------")'); call goPr

       do region = 1, nregions
          write (gol,'(A,E13.6)') 'EMISS-INFO - mass emitted   :            ',sum_emission_all(region); call goPr
       enddo

       CALL Init(io, budget_file_global, 'write', status)
       IF_NOTOK_RETURN(status=1)
       
       CALL WriteAttribute(io, 'sum_emission', sum_emission_all, status)
       IF_NOTOK_RETURN(status=1)

       budemig = 0.0
    END IF

    ! Gather budgets
    REG: DO region = 1, nregions

       CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

       if (isRoot) then 
          ALLOCATE(collect_emissions(im(region), jm(region), nbud_vg, ntracet))
       else
          ALLOCATE(collect_emissions(1,1,1,1) )
       end if

       CALL GATHER( dgrid(region), budemi_dat(region)%emi, collect_emissions, 0, status)
       IF_NOTOK_RETURN(status=1)

       ! write Not-Horizontally-Aggregated-Budgets
       IF (isRoot.and.NHAB) THEN
          CALL Init(Sds,io, 'budemi_dat_'//region_name(region),(/im(region),jm(region),nbud_vg,ntracet/), 'real(4)', status)
          CALL SetDim(Sds, 0, 'im_'//region_name(region),'longitude', (/(j, j=1,im(region))/), status)
          CALL SetDim(Sds, 1, 'jm_'//region_name(region),'latitude', (/(j, j=1,jm(region))/), status)
          CALL SetDim(Sds, 2, 'nbud_vg','vertical layer', (/(j, j=1,nbud_vg)/), status)
          CALL SetDim(Sds, 3, 'ntracet','tracer number', (/(j, j=1,ntracet)/), status)
          IF_NOTOK_RETURN(status=1)
          CALL WriteData(Sds,collect_emissions,status)
          IF_NOTOK_RETURN(status=1)
          CALL Done(Sds, status)
          IF_NOTOK_RETURN(status=1)
       ENDIF

       ! horizontally aggregates budgets
       DO n=1,ntracet
       DO nzone_v=1,nbud_vg
       DO j=j1,j2
       DO i=i1,i2
          nzone = budg_dat(region)%nzong(i,j)
          budemig(nzone,nzone_v,n) = budemig(nzone,nzone_v,n) + budemi_dat(region)%emi(i,j,nzone_v,n)
       END DO
       END DO !j
       END DO !nzone_v
       END DO !nt
       
       DEALLOCATE( collect_emissions      )
       DEALLOCATE( budemi_dat(region)%emi )

    ENDDO REG


    CALL PAR_REDUCE_ELEMENT( budemig, 'SUM', budemig_all, status)
    IF_NOTOK_RETURN(status=1)
    
    ! Write horizontally aggregated budget
    IF ( isRoot ) THEN 
       CALL Init(Sds,io, 'budemi',(/nbudg,nbud_vg,ntracet/), 'real(8)', status)
       IF_NOTOK_RETURN(status=1)
       CALL SetDim(Sds, 0, 'nbudg','horizontal region', (/(j, j=1,nbudg)/), status)
       CALL SetDim(Sds, 1, 'nbud_vg','vertical layer', (/(j, j=1,nbud_vg)/), status)
       CALL SetDim(Sds, 2, 'ntracet','tracer number', (/(j, j=1,ntracet)/), status)
       IF_NOTOK_RETURN(status=1)
       CALL WriteData(Sds,budemig_all,status)
       IF_NOTOK_RETURN(status=1)
       CALL Done(Sds, status)
       IF_NOTOK_RETURN(status=1)
       CALL Done(io, status)
       IF_NOTOK_RETURN(status=1)
    ENDIF
    
#endif /* BUDGETS */

    ! call other emission_*_done routines
    CALL FREE_EMISSION(status)
    IF_NOTOK_RETURN(status=1)
    
    ! -----------------------------------------------------------------------------------
    ! Free biomassburning time splitting factors (now globally, instead of by constituent)
    if( emis_bb_trop_cycle ) then 
       do region = 1, nregions
          deallocate(bb_cycle(region)%scalef)
       end do
    end if
    ! -----------------------------------------------------------------------------------

    ! ok
    status = 0

  END SUBROUTINE EMISSION_DONE
  !EOC

  !---------------------------------------------------------------------------
  !                    TM5                                                   !
  !---------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DECLARE_EMISSION
  !
  ! !DESCRIPTION: Called at run start (init/allocate emiss data) and then at
  !                beginning of every month to just read in data.
  !               Called from SS_MONTHLY_UPDATE.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DECLARE_EMISSION( status )
    !
    ! !USES:
    !
    USE Grid,            ONLY : FillGrid
    USE MDF,             ONLY : MDF_Open, MDF_HDF4, MDF_READ, MDF_Inq_VarID, MDF_Get_Var, MDF_Close
    USE dims,            ONLY : im, jm, lm, newsrun
    USE dims,            ONLY : nregions, iglbsfc, nlat180, nlon360
    USE chem_param        
#ifdef with_m7            
    USE mo_aero_m7,      ONLY : nmod
#endif                    
    USE partools,        ONLY : isRoot
    USE global_data,     ONLY : emis_data
    USE meteodata,       ONLY : global_lli
#ifdef with_online_nox
    use Online_NOx,      only : Online_NOx_Init
!    use online_nox_data, only : mlai2d_onlinenox
#endif
#ifdef with_online_bvoc
    USE Emission_BVOC,   ONLY : Online_BVOC_Init
#endif
    ! AR5/EDGAR4
    use emission_data,   only : emis_input_dir
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)      ::  status
    !
    ! !REVISION HISTORY: 
    !   16 Jul 2010 - A. Strunk - Adapted to revised emission_*.F90 routines
    !   27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !   - anything that is done only if it's a newrun, and that does not require meteo data, should go in INIT
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER          :: rname = mname//'/declare_emission'
    INTEGER                              :: region, imode, hid
    REAL, DIMENSION(:,:), ALLOCATABLE    :: pland
    type(emis_data), dimension(nregions) :: wrk
    
    ! -----------------
    ! Reset M7 emission
    ! -----------------
#ifdef with_m7
    DO region=1,nregions
       DO imode=1,nmod
          emis_number(region,imode)%d4 = 0.0
          emis_mass  (region,imode)%d4 = 0.0
       END DO
    END DO
#endif
    
    ! ---------------------------------------------------------------
    ! ** land fraction
    ! ---------------------------------------------------------------
    IF( newsrun ) THEN

       if(isRoot)then
          ALLOCATE( pland(nlon360,nlat180) )
          DO region=1,nregions
             allocate( wrk(region)%surf(im(region),jm(region)) )
          end DO
       else
          ALLOCATE( pland(1,1) )
          DO region=1,nregions
             allocate( wrk(region)%surf(1,1))
          end DO
       end if

       if (isRoot)then
          CALL MDF_Open( TRIM(emis_input_dir)//'/land/landfraction.hdf', MDF_HDF4, MDF_READ, hid, status )
          IF_NOTOK_RETURN(status=1)
          !       CALL MDF_Inq_VarID( hid, TRIM(filemon), varid, status )  ! more than one 'landfraction' sds
          !       IF_NOTOK_RETURN(status=1)
          CALL MDF_Get_Var( hid, 1, pland, status )
          IF_NOTOK_RETURN(status=1)
          CALL MDF_Close( hid, status )
          IF_NOTOK_RETURN(status=1)

          ! coarsen or distribute to zoom regions:
          DO region = 1, nregions
             ! convert grid:
             CALL FillGrid( global_lli(region),  'n', wrk(region)%surf,                     &
                            global_lli(iglbsfc), 'n', pland,            'area-aver', status )
             IF_NOTOK_RETURN(status=1)
          END DO
       end if

       DO region = 1, nregions
          call scatter( dgrid(region), plandr(region)%surf, wrk(region)%surf, 0, status)
          IF_NOTOK_RETURN(status=1)
          DEALLOCATE( wrk(region)%surf )
       END DO
       
       DEALLOCATE( pland )       

    ENDIF
    

    ! ---------------------------------------------------------------
    ! ** init each constituent
    ! ---------------------------------------------------------------

    ! ** 1st time, initialise emissions

    IF ( newsrun ) THEN
       CALL Emission_NOx_Init( status )
       IF_NOTOK_RETURN(status=1)
#ifdef with_online_nox
       ! init online nox module:
       CALL Online_NOx_Init( status )
       IF_NOTOK_RETURN(status=1)
#endif
       CALL Emission_NMVOC_Init( status )
       IF_NOTOK_RETURN(status=1)

       CALL Emission_NH3_Init( status )
       IF_NOTOK_RETURN(status=1)

       CALL Emission_CO_Init( status )
       IF_NOTOK_RETURN(status=1)

       CALL Emission_CH4_Init( status )
       IF_NOTOK_RETURN(status=1)

       CALL Emission_SOx_Init( status )
       IF_NOTOK_RETURN(status=1)

       CALL Emission_DMS_Init( status )
       IF_NOTOK_RETURN(status=1)

       CALL Emission_rn222_Init( status )
       IF_NOTOK_RETURN(status=1)

#ifdef with_m7
       call emission_bc_init( status )
       IF_NOTOK_RETURN(status=1)

       call emission_pom_init( status )
       IF_NOTOK_RETURN(status=1)
#endif
       
    END IF
    

    ! ** every month, read and re-grid
    
    CALL emission_nox_declare( status )
    IF_NOTOK_RETURN(status=1)
    
    CALL emission_nmvoc_declare( status )
    IF_NOTOK_RETURN(status=1)

    CALL emission_nh3_declare( status )
    IF_NOTOK_RETURN(status=1)

    CALL emission_co_declare( status )
    IF_NOTOK_RETURN(status=1)

    CALL emission_ch4_declare( status )
    IF_NOTOK_RETURN(status=1)

    CALL emission_sox_declare( status )
    IF_NOTOK_RETURN(status=1)

    CALL emission_dms_declare( status )
    IF_NOTOK_RETURN(status=1)

    
    ! ** special case of bio voc/isoprene
#ifdef with_online_bvoc
    CALL Online_BVOC_Init( status )
    IF_NOTOK_RETURN(status=1)

    CALL declare_emission_bvoc( status )
#else
    IF ( newsrun ) THEN
       CALL Emission_isop_Init( status )
       IF_NOTOK_RETURN(status=1)
    END IF

    CALL emission_isop_declare( status )
    IF_NOTOK_RETURN(status=1)
#endif

    CALL emission_rn222_declare( status )
    IF_NOTOK_RETURN(status=1)

#ifdef with_m7
    CALL emission_bc_declare( status)
    IF_NOTOK_RETURN(status=1)

    CALL emission_pom_declare( status )
    IF_NOTOK_RETURN(status=1)
  
    CALL emission_dust_declare(status)
    IF_NOTOK_RETURN(status=1)
    
!!$    CALL declare_emission_ss(status)
!!$    IF_NOTOK_RETURN(status=1)
#endif

    ! ok
    status = 0

  END SUBROUTINE DECLARE_EMISSION
  !EOC


  !-------------------------------------------------------------------
  !                    TM5                                           !
  !-------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    EMISSION_APPLY
  !
  ! !DESCRIPTION: Call emission_apply methods of constituent modules. 
  !               --> add current emissions to tracers array. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE EMISSION_APPLY( region, status )
    !
    ! !USES:
    !
    USE chem_param
    !
    ! !INPUT PARAMETERS:
    !
    INTEGER, INTENT(in)        ::  region
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)       ::  status
    !
    ! !REVISION HISTORY: 
    !   16 Jul 2010 - A. Strunk   - Adapted to revised emission_*.F90 routines
    !   27 Mar 2012 - P. Le Sager - cleanup for lat-lon mpi decomposition
    !
    !EOP
    !-----------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER ::  rname = mname//'/emission_apply'

    ! --- begin -------------------------------------- 

    ! start timing:
    call GO_Timer_Start( itim_appl, status )
    IF_NOTOK_RETURN(status=1)

    IF (okdebug) then
       WRITE(gol,*) 'start of emission_apply for region:',region ; call goPr
    END IF

    ! CO emissions 
    call GO_Timer_Start( itim_co, status )
    IF_NOTOK_RETURN(status=1)
    CALL emission_co_apply( region, status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_End( itim_co, status )
    IF_NOTOK_RETURN(status=1)

    ! CH4 emissions 
    call GO_Timer_Start( itim_ch4, status )
    IF_NOTOK_RETURN(status=1)
    CALL emission_ch4_apply(region, status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_End( itim_ch4, status )
    IF_NOTOK_RETURN(status=1)

    !  biogenic NMHC emissions (isoprene)
    call GO_Timer_Start( itim_isop, status )
    IF_NOTOK_RETURN(status=1)
#ifdef with_online_bvoc
    CALL emission_apply_bvoc( region, status )
#else
    CALL emission_isop_apply( region, status )
    IF_NOTOK_RETURN(status=1)
#endif
    call GO_Timer_End( itim_isop, status )
    IF_NOTOK_RETURN(status=1)

    ! add di-methyl sulfide emissions:
    call GO_Timer_Start( itim_dms, status )
    IF_NOTOK_RETURN(status=1)
    CALL emission_dms_apply(region, status)
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_End( itim_dms, status )
    IF_NOTOK_RETURN(status=1)

    ! add SOx emissions:
    call GO_Timer_Start( itim_sox, status )
    IF_NOTOK_RETURN(status=1)
    CALL emission_sox_apply( region, status ) 
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_End( itim_sox, status )
    IF_NOTOK_RETURN(status=1)

    ! add NH3 emissions:                  
    call GO_Timer_Start( itim_nh3, status )
    IF_NOTOK_RETURN(status=1)
    CALL emission_nh3_apply(region, status)
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_End( itim_nh3, status )
    IF_NOTOK_RETURN(status=1)

    ! add Rn222 emissions:
    call GO_Timer_Start( itim_rn222, status )
    IF_NOTOK_RETURN(status=1)
    CALL emission_rn222_apply(region, status)
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_End( itim_rn222, status )
    IF_NOTOK_RETURN(status=1)

    ! black carbon and particulate organic matter emissions are added in the sedimentation routine...
    
    ! seasalt and dust (which also sediment) are added in tracer_after_read...

    ! add non-methane voc emissions:
    call GO_Timer_Start( itim_voc, status )
    IF_NOTOK_RETURN(status=1)
    CALL emission_nmvoc_apply( region, status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_End( itim_voc, status )
    IF_NOTOK_RETURN(status=1)

    IF(okdebug) then
       WRITE(gol,*) 'End of adding emission '; call goPr
    END IF

    ! end timing:
    call GO_Timer_End( itim_appl, status )
    IF_NOTOK_RETURN(status=1)
    ! ok
    status = 0

  END SUBROUTINE EMISSION_APPLY
  !EOC

  !----------------------------------------------------------------------------
  !                    TM5                                                    !
  !----------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    FREE_EMISSION
  !
  ! !DESCRIPTION: Deallocate space needed to handle the emissions  by calling 
  !               *done methods of constituents' modules.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE FREE_EMISSION( status )
    !
    ! !USES:
    !
    USE dims,          ONLY : nregions
#ifdef with_m7
    USE mo_aero_m7,    ONLY : nmod
#endif
#ifdef with_online_nox
    USE Online_NOx,    ONLY : Online_NOx_Done
#endif
#ifdef with_online_bvoc
    USE Emission_BVOC, ONLY : Online_BVOC_Done
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    INTEGER, INTENT(out)        ::  status
    !
    ! !REVISION HISTORY: 
    !   16 Jul 2010 - A. Strunk   - Adapted to revised emission_*.F90 routines
    !   27 Mar 2012 - P. Le Sager - Adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !---------------------------------------------------------------------------
    !BOC

    CHARACTER(len=*), PARAMETER ::  rname = mname//'/free_emission'

    INTEGER                     ::  region, imode

    ! --- begin -----------------------------------

    DO region = 1, nregions

       DEALLOCATE(plandr(region)%surf)
       DEALLOCATE(emis2D(region)%surf)
#ifdef with_m7
       ! aerosols:
       DO imode = 1, nmod
          DEALLOCATE(emis_number(region,imode)%d4)
          DEALLOCATE(emis_mass(region,imode)%d4)
       ENDDO
       DEALLOCATE(emis_temp(region)%surf)
#endif
    ENDDO

    CALL Emission_NOx_Done( status )
    IF_NOTOK_RETURN(status=1)

#ifdef with_online_nox
    CALL Online_NOx_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif

    CALL emission_nh3_done( status )
    IF_NOTOK_RETURN(status=1)

    CALL Emission_CO_Done( status )
    IF_NOTOK_RETURN(status=1)

    CALL Emission_CH4_Done( status )
    IF_NOTOK_RETURN(status=1)

    CALL Emission_SOx_Done( status )
    IF_NOTOK_RETURN(status=1)

    CALL emission_dms_done( status )
    IF_NOTOK_RETURN(status=1)

#ifdef with_online_bvoc
    CALL free_emission_bvoc ( status )
    ! done with online bvoc module:
    CALL Online_BVOC_Done( status )
    IF_NOTOK_RETURN(status=1)
#else
    CALL Emission_isop_Done ( status )
    IF_NOTOK_RETURN(status=1)
#endif

    CALL emission_rn222_done( status )
    IF_NOTOK_RETURN(status=1)

#ifdef with_m7
    CALL emission_bc_done( status )
    IF_NOTOK_RETURN(status=1)

    CALL emission_pom_done( status )
    IF_NOTOK_RETURN(status=1)

    CALL emission_dust_done
!!$    CALL free_emission_ss
#endif

    CALL Emission_NMVOC_Done( status )
    IF_NOTOK_RETURN(status=1)

    ! done
    status = 0

  END SUBROUTINE FREE_EMISSION
  !EOC


END MODULE EMISSION