ecearth3/sources/tm5mp/proj/cbm4/chemistry.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:     CHEMISTRY
!
! !DESCRIPTION: perform CBM4 chemistry simulation in TM5 + M7
!\\
!\\
! !INTERFACE: 
!
MODULE CHEMISTRY
  !
  ! !USES:
  !
  use GO,              only : gol, goErr, goPr, goBug
  use tm5_distgrid,    only : dgrid, Get_DistGrid, gather
  use dims,            only : nregions
#ifdef with_budgets
  use budget_global,   only : nbudg,  nbud_vg
  use chem_param,      only : ntrace, ntracet
#endif

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public  ::  Chemistry_Init, Chemistry_Done
  public  ::  Chemie
  !
  ! !PRIVATE DATA MEMBERS:
  !
#ifdef with_budgets
  real, dimension(:,:,:), allocatable :: budchem
  real, Dimension(:,:,:), Allocatable :: budeqsam
  real :: eminox
#ifdef with_m7
  real, dimension(:,:,:), allocatable :: budm7proc
  real, dimension(:,:,:), allocatable :: budm7
#endif
#endif

  character(len=*), parameter  ::  mname = 'chemistry'
  !
  ! minimal concentration required to apply scaling for slopes
  ! (to avoid floating point overflows):
  real, parameter   ::  eps = 1.0e-20
  !
  ! !REVISION HISTORY: 
  !          2005 - Elisabetta Vignati - changed for coupling with M7
  !   22 Mar 2012 - P. Le Sager  - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CHEMISTRY_INIT
  !
  ! !DESCRIPTION: allocate and initialize budget arrays. Init M7 if needed.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CHEMISTRY_INIT( status )
    !
    ! !USES:
    !
    use Dims      , only  :  nregions, im, jm, lm
#ifdef with_budgets
    use Chem_param, only  :  ntrace, ntracet
    use ebischeme , only  :  sum_deposition, sum_wetchem, sum_chemistry
    use ebischeme , only  :  buddep_dat, budrrg, budrjg, budrwg
    !PLS    use chem_param, only  :  ch4oh, so4pg, so4pa, o3p, o3l
    use chem_param, only  :  o3p, o3l
#ifdef with_m7
    use m7_data,       only : nm7procs
#endif
#endif
#ifdef with_m7
    use mo_aero_m7,   only : m7_initialize
    use m7_data,      only : init_m7_data
#endif
    use meteoData   , only : Set
    use MeteoData   , only : phlb_dat, m_dat, temper_dat, humid_dat, gph_dat
    use MeteoData   , only : lwc_dat, iwc_dat, cc_dat
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)           ::  status
    !
    ! !REVISION HISTORY: 
    !   22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer                     ::  region, i1, i2, j1, j2
    character(len=*), parameter ::  rname = mname//'/Chemistry_Init'

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

    ! select meteo to be used:
    do region = 1, nregions
       call Set(   phlb_dat(region), status, used=.true. ) 
       call Set(      m_dat(region), status, used=.true. ) 
       call Set( temper_dat(region), status, used=.true. )
       call Set(  humid_dat(region), status, used=.true. )
       call Set(    gph_dat(region), status, used=.true. )
       call Set(    lwc_dat(region), status, used=.true. ) 
       call Set(    iwc_dat(region), status, used=.true. ) 
       call Set(     cc_dat(region), status, used=.true. ) 
    enddo

#ifdef with_budgets

    do region = 1, nregions
       sum_chemistry (region) = 0.0
       sum_deposition(region) = 0.0
       sum_wetchem   (region) = 0.0

       call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
       
       allocate(buddep_dat(region)%dry(i1:i2, j1:j2, ntrace))
       buddep_dat(region)%dry = 0.0

       allocate(o3p(region)%d3(i1:i2, j1:j2, lm(region))); o3p(region)%d3 = 0.0
       allocate(o3l(region)%d3(i1:i2, j1:j2, lm(region))); o3l(region)%d3 = 0.0
       ! >>>
       ! TvN: still in sources_sinks.F90
       !       allocate(ch4oh(region)%d3(i1:i2, j1:j2, lm(region))); ch4oh(region)%d3 = 0.0
       !       allocate(so4pg(region)%d3(i1:i2, j1:j2, lm(region))); so4pg(region)%d3 = 0.0
       !       allocate(so4pa(region)%d3(i1:i2, j1:j2, lm(region))); so4pa(region)%d3 = 0.0
       ! <<<
    enddo
    
    budrrg = 0.0
    budrjg = 0.0
    budrwg = 0.0

    Allocate(budchem(nbudg,nbud_vg,ntrace))
    budchem = 0.0

    ! Two acid-base reactions
    Allocate(budeqsam(nbudg,nbud_vg,2))
    budeqsam = 0.0
    
#ifdef with_m7
    ! M7 Bugets
    Allocate(budm7proc(nbudg,nbud_vg,nm7procs))
    Allocate(budm7(nbudg,nbud_vg,ntracet))
    budm7proc = 0.0
    budm7 = 0.0
#endif
    
    ! global lightning NOx
    eminox = 0.

#endif /* BUDGETS */

#ifdef with_m7
    CALL M7_INITIALIZE
    CALL INIT_M7_DATA( status )                   ! declare M7 output data structure
    IF_NOTOK_RETURN(status=1)
#endif

    ! ok
    status = 0

  END SUBROUTINE CHEMISTRY_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CHEMISTRY_DONE
  !
  ! !DESCRIPTION: gather and write budgets
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CHEMISTRY_DONE( status )
    !
    ! !USES:
    !
    use dims,               only : nregions, im, jm, lm
    
#ifdef with_budgets
    use dims,               only : region_name
    use file_hdf,           only : Init, Done, WriteAttribute, WriteData, SetDim, THdfFile, TSds
    use budget_global,      only : budget_file_global, budg_dat, NHAB
#ifndef without_boundary
    use boundary,           only : budstratg
#endif
#ifndef without_emission
    use emission,           only : budemig_all
#endif
    use partools,           only : isRoot, par_reduce_element, par_reduce
    use ebischeme,          only : sum_deposition, sum_wetchem, sum_chemistry
    use ebischeme,          only : buddep_dat, budrrg, budrjg, budrwg, budmarkg
    !PLS    use chem_param,         only : ch4oh, so4pg, so4pa, o3p, o3l, marknam, nmark
    use chem_param,         only : o3p, o3l, marknam, nmark
    use chem_param,         only : marknam, ntrace
    use chem_param,         only : inox
    use reaction_data,      only : nreac, ratnam, nreacw, rwnam
    use photolysis_data,    only : nj,    jnam
#ifdef with_cariolle
    use chem_param,         only :  ncar, carnam
    use chemistry_cariolle, only :  budcarg
#endif
#ifdef with_m7
    use m7_data,            only : nm7procs
#endif
#endif /* BUDGETS */

#ifdef with_m7
    use m7_data,            only : free_m7_data
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)         ::  status
    !
    ! !REVISION HISTORY: 
    !   22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter              :: rname = mname//'/Chemistry_Done'

    integer                                  ::  region
#ifdef with_budgets
    type(THdfFile)                           ::  io
    type(TSds)                               ::  sds
    real, dimension(:,:,:), allocatable      ::  collectb
    real,dimension(nbudg,nbud_vg,ntrace)     ::  buddry
    real,dimension(nbudg,nbud_vg,ntrace)     ::  budchem_result
    integer                                  ::  i, j, l, nzone, n, nsend, i1, i2, j1, j2

    ! for buggy MPI (see comment in budget_global.F90 for details)
    real, dimension(nregions) :: sum_chemistry_all, sum_wetchem_all, sum_deposition_all
    real, dimension(nbudg,nbud_vg,nreac)    :: budrrg_all
    real, dimension(nbudg,nbud_vg,nj)       :: budrjg_all
    real, dimension(nbudg,nbud_vg,nreacw)   :: budrwg_all
    real, dimension(nbudg,nbud_vg,nmark)    :: budmarkg_all
    real, dimension(nbudg,nbud_vg,ntrace)   :: buddry_all
    real, dimension(nbudg,nbud_vg,ntrace)   :: budchem_all
    real, dimension(nbudg,nbud_vg,2)        :: budeqsam_all
    real :: eminox_all
#ifndef without_boundary
    real, dimension(nbudg,nbud_vg,ntracet)  :: budstratg_all                    
#endif
#ifdef with_cariolle
    real, dimension(nbudg,nbud_vg,ncar)     :: budcarg_all
#endif
#ifdef with_m7
    real, dimension(nbudg,nbud_vg,ntracet)  :: budm7_all
    real, dimension(nbudg,nbud_vg,nm7procs) :: budm7proc_all
#endif

#ifdef with_m7
    ! Do not change this order. And if you do, look through the M7 code to pprocess (parameter process) where it is written with a fixed paramter. You can, however, change this order if it is made dynamic.
    Character(len=8), Dimension(nm7procs), Parameter :: m7procnames = (/'NucN    ','NucSU   ','Coag11N ','Coag12N ','Coag12SU',&
         'Coag13N ','Coag13SU','Coag14N ','Coag14SU','Coag15N ','Coag15SU','Coag16N ','Coag16SU','Coag17N ','Coag17SU',&
         'Coag22N ','Coag23N ','Coag23SU','Coag23BC','Coag23OC','Coag24N ','Coag24SU','Coag24BC','Coag24OC','Coag25N ',&
         'Coag25BC','Coag25OC','Coag26N ','Coag26SU','Coag26BC','Coag26OC','Coag27N ','Coag27SU','Coag27BC','Coag27OC',&
         'Coag33N ','Coag35N ','Coag35BC','Coag35OC','Coag55N ','Cond1   ','Cond2   ','Cond3   ','Cond4   ','Cond5   ',&
         'Cond6   ','Cond7   ','Age5N   ','Age5BC  ','Age5OC  ','Age6N   ','Age6DU  ','Age7N   ','Age7DU  ','Grow1N  ',&
         'Grow1SU ','Grow2N  ','Grow2SU ','Grow2BC ','Grow2OC ','Grow3N  ','Grow3SU ','Grow3BC ','Grow3OC ','Grow3SS ','Grow3DU '/)
    ! Array for left hand sides, static attribute for the M7 process budget data set.
    Character(len=8), Dimension(nm7procs), Parameter :: m7proclhs = (/'void    ','H2SO4   ','NUS_N   ','NUS_N   ','SO4NUS  ',&
         'NUS_N   ','SO4NUS  ','NUS_N   ','SO4NUS  ','NUS_N   ','SO4NUS  ','NUS_N   ','SO4NUS  ','NUS_N   ','SO4NUS  ',&
         'AIS_N   ','AIS_N   ','SO4AIS  ','BCAIS   ','POMAIS  ','AIS_N   ','SO4AIS  ','BCAIS   ','POMAIS  ','AII_N   ',&
         'BCAII   ','POMAII  ','AIS_N   ','SO4AIS  ','BCAIS   ','POMAIS  ','AIS_N   ','SO4AIS  ','BCAIS   ','POMAIS  ',&
         'ACS_N   ','AII_N   ','BCAII   ','POMAII  ','AII_N   ','H2SO4   ','H2SO4   ','H2SO4   ','H2SO4   ','H2SO4   ',&
         'H2SO4   ','H2SO4   ','AII_N   ','BCAII   ','POMAII  ','ACI_N   ','DUACI   ','COI_N   ','DUCOI   ','NUS_N   ',&
         'SO4NUS  ','AIS_N   ','SO4AIS  ','BCAIS   ','POMAIS  ','ACS_N   ','SO4ACS  ','BCACS   ','POMACS  ','SSACS   ','DUACS   '/)
    ! Array for right hand sides, static attribute for the M7 process budget data set.
    Character(len=8), Dimension(nm7procs), Parameter :: m7procrhs = (/'NUS_N   ','SO4NUS  ','void    ','void    ','SO4AIS  ',&
         'void    ','SO4ACS  ','void    ','SO4COS  ','void    ','SO4AIS  ','void    ','SO4ACS  ','void    ','SO4COS  ',&
         'void    ','void    ','SO4ACS  ','BCACS   ','POMACS  ','void    ','SO4COS  ','BCCOS   ','POMCOS  ','void    ',&
         'BCAIS   ','POMAIS  ','void    ','SO4ACS  ','BCACS   ','POMACS  ','void    ','SO4COS  ','BCCOS   ','POMCOS  ',&
         'void    ','void    ','BCACS   ','POMACS  ','void    ','SO4NUS  ','SO4AIS  ','SO4ACS  ','SO4COS  ','SO4AIS  ',&
         'SO4ACS  ','SO4COS  ','AIS_N   ','BCAIS   ','POMAIS  ','ACS_N   ','DUACS   ','COS_N   ','DUCOS   ','AIS_N   ',&
         'SO4AIS  ','ACS_N   ','SO4ACS  ','BCACS   ','POMACS  ','COS_N   ','SO4COS  ','BCCOS   ','POMCOS  ','SSCOS   ','DUCOS   '/)
#endif

#endif /* BUDGETS */

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

#ifdef with_budgets

    ! Sum up contribution from each proc into root array    
    call PAR_REDUCE_ELEMENT(budrrg, 'SUM', budrrg_all, status )

    do region = 1, nregions
    
       CALL PAR_REDUCE(sum_chemistry(region), 'SUM', sum_chemistry_all(region), status)
       IF_NOTOK_RETURN(status=1)

       CALL PAR_REDUCE(sum_deposition(region), 'SUM', sum_deposition_all(region), status)
       IF_NOTOK_RETURN(status=1)

       CALL PAR_REDUCE(sum_wetchem(region), 'SUM', sum_wetchem_all(region), status)
       IF_NOTOK_RETURN(status=1)

    end do
    
    ! Sum up contribution from each proc into root array for "total LiNOx"
    call PAR_REDUCE( eminox, 'SUM', eminox_all, status )
    eminox_all = eminox_all*1.e-9 ! kg -> Tg

    if ( isRoot ) then

       call Init(io, budget_file_global, 'write', status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(io, 'sum_chemistry' ,  sum_chemistry_all,  status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(io, 'sum_wetchem' ,    sum_wetchem_all,    status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(io, 'sum_deposition' , sum_deposition_all, status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(io, 'total_lightning_nox_emissions_Tg_N' , eminox_all, status)
       IF_NOTOK_RETURN(status=1)

       call Init(Sds, io, 'budrr', (/nbudg,nbud_vg,nreac/), 'real(8)', status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(Sds, 'ratnam', ratnam, 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, 'nreac',  'reaction number',   (/(j, j=1,nreac)/),   status)
       IF_NOTOK_RETURN(status=1)
       call WriteData(Sds, BUDRRG_all, status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)
    end if


    ! Sum up contribution from each proc into root array    
    call PAR_REDUCE_ELEMENT(budrjg, 'SUM', budrjg_all, status )    

    if ( isRoot ) then

       call Init(Sds,io, 'budrj',(/nbudg,nbud_vg,nj/), 'real(8)', status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(Sds, 'jnam', jnam, 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, 'nj','reaction number', (/(j, j=1,nj)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteData(Sds, BUDRJG_all, status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

    end if

    ! Lets us write down the WET CHEMISTRY
    
    call PAR_REDUCE_ELEMENT(budrwg, 'SUM', budrwg_all, status )

    if ( isRoot ) then

       call Init(Sds,io, 'budrw',(/nbudg,nbud_vg,nreacw/), 'real(8)', status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(Sds, 'rwnam', rwnam, 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, 'nreacw','wet reaction number', (/(j, j=1,nreacw)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteData(Sds, BUDRWG_all, status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

    end if

    
    ! Sum up contribution from each proc into root array    
    call PAR_REDUCE_ELEMENT(budmarkg, 'SUM', budmarkg_all, status )

    if ( isRoot ) then

       call Init(Sds,io, 'budmark',(/nbudg,nbud_vg,nmark/), 'real(8)', status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(Sds, 'marknam', marknam, 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, 'nmark','marked tracer number', (/(j, j=1,nmark)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteData(Sds, BUDMARKG_all, status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

    end if


    !  DRYDEP budget    
    buddry  = 0.0

    do region=1,nregions
       
       ! aggregates horizontally budget
       call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

       do n=1,ntrace
       do j=j1,j2
       do i=i1,i2
                nzone = budg_dat(region)%nzong(i,j)
                buddry(nzone,1,n) = buddry(nzone,1,n) + buddep_dat(region)%dry(i,j,n)
       end do
       end do
       end do

       call PAR_REDUCE_ELEMENT(buddry, 'SUM', buddry_all, status ) ! contribution from all procs
       
       ! gather and write Non-Horizontally-Aggregated-Budgets
       if (NHAB) then
          if ( isRoot ) then
             allocate( collectb(im(region),jm(region),ntrace ))
          else
             allocate( collectb(1,1,1))
          end if

          call GATHER( dgrid(region), buddep_dat(region)%dry, collectb, 0, status)

          if(isRoot) then
             call Init(Sds,io, 'buddep_dat_dry_'//region_name(region),(/im(region),jm(region),ntrace/), '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, 'ntrace','tracer index', (/(j, j=1,ntrace)/), status)
             IF_NOTOK_RETURN(status=1)
             call WriteData(Sds, collectb, status)
             IF_NOTOK_RETURN(status=1)
             call Done(Sds, status)
             IF_NOTOK_RETURN(status=1)
          endif
          deallocate(collectb)
          
       endif ! NHAB
       
    enddo ! regions

    !-- write more aggregated budgets
    
    call PAR_REDUCE_ELEMENT(budchem,  'SUM', budchem_all,   status )
    call PAR_REDUCE_ELEMENT(budeqsam, 'SUM', budeqsam_all,  status )
#ifndef without_boundary
    call PAR_REDUCE_ELEMENT(budstratg,'SUM', budstratg_all, status )
#endif
#ifdef with_m7
    call PAR_REDUCE_ELEMENT(budm7,    'SUM', budm7_all,     status )
    call PAR_REDUCE_ELEMENT(budm7proc,'SUM', budm7proc_all, status )
#endif
    
    
    if(isRoot) then

       call Init(Sds,io, 'buddry',(/nbudg,nbud_vg,ntrace/), '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, 'ntrace','tracer index', (/(j, j=1,ntrace)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteData(Sds,buddry_all,status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

       Call Init(Sds,io,"budchem",(/nbudg,nbud_vg,ntrace/),'real(8)',status)
       call SetDim(Sds, 0, 'nbudg','horizontal budget', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg','vertical budget', (/(j, j=1,nbud_vg)/), status)
       call SetDim(Sds, 2, 'ntrace','tracer number', (/(j, j=1,ntrace)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(Sds, 'comment', 'EBI budget stuffed together. Hopefully it does not explode', status)

       ! correct chemistry budget
       budchem_result = budchem_all+buddry_all
#ifndef without_emission
       budchem_result(:,:,inox) = budchem_result(:,:,inox) - budemig_all(:,:,inox)
#endif
       call WriteData(Sds,budchem_result,status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

       Call Init(Sds,io,"budeqsam",(/nbudg,nbud_vg,2/),'real(8)',status)
       call SetDim(Sds, 0, 'nbudg','horizontal budget', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg','vertical budget', (/(j, j=1,nbud_vg)/), status)
       call SetDim(Sds, 2, 'AB_reac','Acid-base reaction number', (/(j, j=1,2)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(Sds, 'comment', 'EQSAM evil acid base reaction. Positive forms aerosol material.', status)
       call WriteData(Sds,budeqsam_all,status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

#ifndef without_boundary
       Call Init(Sds,io,"budstrat",(/nbudg,nbud_vg,ntracet/),'real(8)',status)
       call SetDim(Sds, 0, 'nbudg','horizontal budget', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg','vertical budget', (/(j, j=1,nbud_vg)/), status)
       call SetDim(Sds, 2, 'ntracet','tracer number', (/(j, j=1,ntracet)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(Sds, 'comment', 'Stratosphere (from boundary.F90)', status)
       call WriteData(Sds,budstratg_all,status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)
#endif
    endif

    
#ifdef with_cariolle
    call PAR_REDUCE_ELEMENT(budcarg,  'SUM', budcarg_all,   status )
    
    if ( isRoot ) then

       call Init(Sds,io, 'budcar',(/nbudg,nbud_vg,ncar/), 'real(8)', status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(Sds, 'carnam', carnam, 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, 'ncar','Cariolle tracer number', (/(j, j=1,ncar)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteData(Sds,budcarg_all,status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

    end if
#endif

#ifdef with_m7
    if ( isRoot ) then
       ! M7 process budget
       Call Init(Sds,io,"budm7proc",(/nbudg,nbud_vg,nm7procs/),'real(8)',status)
       call SetDim(Sds, 0, 'nbudg','horizontal budget', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg','vertical budget', (/(j, j=1,nbud_vg)/), status)
       call SetDim(Sds, 2, 'nm7proc','M7 process number', (/(j, j=1,nm7procs)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(Sds, 'procnames',m7procnames, status)
       call WriteAttribute(Sds, 'lefthandsides',m7proclhs, status)
       call WriteAttribute(Sds, 'righthandsides',m7procrhs, status)
       call WriteData(Sds,budm7proc_all,status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)

       ! M7 budget
       Call Init(Sds,io,"budm7",(/nbudg,nbud_vg,ntracet/),'real(8)',status)
       call SetDim(Sds, 0, 'nbudg','horizontal budget', (/(j, j=1,nbudg)/), status)
       call SetDim(Sds, 1, 'nbud_vg','vertical budget', (/(j, j=1,nbud_vg)/), status)
       call SetDim(Sds, 2, 'ntracet','tracer number', (/(j, j=1,ntracet)/), status)
       IF_NOTOK_RETURN(status=1)
       call WriteAttribute(Sds, 'comment', 'M7 budget stuffed together. Hopefully the M7 tracers have now correct budgets', status)
       call WriteData(Sds,budm7_all,status)
       IF_NOTOK_RETURN(status=1)
       call Done(Sds, status)
       IF_NOTOK_RETURN(status=1)
    endif
#endif

    if ( isRoot ) then
       call Done(io, status)
       IF_NOTOK_RETURN(status=1)
    endif

    do region = 1, nregions
       deallocate(buddep_dat(region)%dry)
       deallocate(o3p(region)%d3)
       deallocate(o3l(region)%d3)
    enddo

    deallocate(budchem)
    deallocate(budeqsam)
#ifdef with_m7
    deallocate(budm7proc)
    deallocate(budm7)
#endif

#endif /* BUDGETS */

#ifdef with_m7
    call free_m7_data
#endif

    ! ok
    status = 0

  END SUBROUTINE CHEMISTRY_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CHEMIE
  !
  ! !DESCRIPTION: performs chemistry transformations of the tracers
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CHEMIE( region, tr, status )
    !
    ! !USES:
    !
    use GO,            only : TDate
    use binas,         only : Avog, pi, xmair, Rgas_air
    use dims,          only : itaur, nchem, tref, sec_month, revert, okdebug
    use dims,          only : dx, xref, dy, yref, im, jm, lm
    use dims,          only : nregions, ndyn, ndyn_max
    use dims,          only : xbeg, ybeg, adv_scheme, mlen
    use chem_param,    only : xmbc, xmpom, xmnacl, xmdust !EV
    use chem_param,    only : maxtrace, n_extra, iso4, ino3_a, inh4, ih2opart
    use chem_param,    only : irinc, ino, xmn, inox, ieno, iairm, iairn, ntrace
    use chem_param,    only : ntracet, fscale, imsa, iacid, iair, ih2o, io2, ih2on
    use chem_param,    only : iairm,i_temp, i_pres, inh3, ihno3, xmh2o,  names
#ifdef with_m7
    use chem_param,    only : iso4nus, iso4ais, iso4acs, iso4cos, ibcais, ibcacs, ibccos, ibcaii, ipomais, ipomacs
    use chem_param,    only : ipomcos, ipomaii, issacs, isscos, iduacs, iducos, iduaci,iducoi
    use chem_param,    only : inus_n, xmnumb, iais_n, iacs_n, icos_n, iaii_n, iaci_n, icoi_n
    use chem_param,    only : nh4no3_factor, nh4no3_density, msa_density
    use mo_aero_m7,    only : cmr2ram, ram2cmr
    use binas,         only : rol
#endif
    use chem_param,    only : idz
    use meteodata,     only : m_dat, gph_dat
    use reaction_data, only : nreac, aerdens, kn2o5aq
    use chem_rates,    only : calrates
    use eqsam_param,   only : eqsam_v03d
    use toolbox,       only : dumpfield
    use datetime,      only : tau2date, get_day
    use global_data,   only : region_dat, mass_dat, chem_dat
#ifndef without_photolysis
#ifdef with_optics
    use photolysis,    only : aerosol_info_m7
#else
    use photolysis,    only : aerosol_info
#endif
    use photolysis,    only : daysim, get_sza, photo_flux
    use photolysis_data
#endif
    use ebischeme,     only : ebi
#ifndef without_emission
    use dims,          only : at, bt 
    use emission_nox,  only : eminox_lightning,   nox_emis_3d
    use emission_nox,  only : nox_emis_3d_bb_app, emission_nox_bb_daily_cycle
#endif
#ifdef with_budgets
    use emission_data, only : budemi_dat, sum_emission
    use budget_global, only : budg_dat,   nzon_vg
#endif
#ifndef without_sedimentation
    use sedimentation, only : rh
#endif
#ifdef with_cariolle
    use chem_param,         only : io3
    use chemistry_cariolle, only : l_cariolle, with_trop_force
#endif

!    use tracer_data, only : tracer_print  ! for debugging

#ifdef with_m7
    use mo_aero_m7,    only : naermod, nmod, nsol
    use m7_data,       only : rw_mode, rwd_mode, dens_mode, h2o_mode
    !use emission_data, only : oc2pom !factor for conversion of OC mass to POM
#ifdef with_budgets
    Use m7_data,       only : nm7procs
#endif
    !#ifdef with_optics
    !    Use Optics,        only : nwl, aod_count, AOD, calculate_aop
    !    Use chem_param,    only : xmso4, xmno3
    !#endif
    external  ::  m7
#endif /* M7 */
    !
    ! !INPUT PARAMETERS:
    !
    integer,     intent(in)      ::  region
    type(TDate), intent(in)      ::  tr(2)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)         ::  status
    !
    ! !REVISION HISTORY: 
    !   22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/chemie'

    ! --- local ------------------------------

    real,dimension(:,:,:),pointer     :: m
    real,dimension(:,:,:,:),pointer   :: rm
#ifdef slopes
    real,dimension(:,:,:,:),pointer   :: rxm, rym, rzm
#ifdef secmom
    real,dimension(:,:,:,:),pointer   :: rxxm, rxym, rxzm, ryym, ryzm, rzzm
#endif
#endif
    real,dimension(:,:,:,:),pointer   :: rmc

    real,dimension(:,:,:),pointer     :: r_cloud

    ! store per level actinic fluxes, photolysis rates, and reaction rates
#ifndef without_photolysis
    real,dimension(:,:,:),allocatable  :: fact
#endif
    real,dimension(:,:,:),allocatable  :: rj, rr
    ! new photolysis array
    real,dimension(:,:,:,:),allocatable :: rj_new    
    ! store concentrations for chemistry ....
    real,dimension(:,:,:),allocatable  :: y, y0, yhelp
    ! store information that is transferred
    ! between chemistry modules
    real,dimension(:,:,:),allocatable  :: ye
    real,dimension(:,:),  allocatable  :: th, tha, sza, mu, reff
    real                                :: starttime,deltat
    ! eqsam:
    real,dimension(:,:), allocatable   :: yi, yo
    ! sad:
    real,dimension(:,:,:), allocatable :: sad_aer_out, sad_cld_out, sad_ice_out

#ifdef with_m7
    ! m7:
    real, dimension(:),     pointer     :: dxyp
    real, dimension(:,:),   allocatable :: presm7, rhm7, tempm7, dryairm7 ! meteo parameters
    real, dimension(:,:),   allocatable :: h2so4m7       ! h2so4 gas phase
    real, dimension(:,:,:), allocatable :: aemassm7      ! aerosol mass
    real, dimension(:,:,:), allocatable :: aenumm7       ! aerosol number

    !#ifdef with_optics
    !    Integer :: iband
    !    real,dimension(:,:,:),allocatable :: Ext   ! extinction coefficient [1/m]
    !    real,dimension(:,:,:),allocatable :: a     ! single-scattering albedo [unitless]
    !    real,dimension(:,:,:),allocatable :: g     ! asymmetry parameter [unitless]
    !#endif
    real, dimension(:,:,:), allocatable :: aedensm7     ! aerosol density
    real, dimension(:,:,:), allocatable :: aewatm7      ! aerosol water content
    real, dimension(:,:,:), allocatable :: aeradm7      ! aerosol radius in equilibrium with water
    real, dimension(:,:,:), allocatable :: aeradrm7     ! aerosol dry radius
    !
    ! To convert units from M7 to TM5 units
    real, parameter :: convnumb = 1.e3/xmnumb*Avog      ! This value is needed directly
    real, parameter :: convsu = 1.                      ! Nothing happens for sulfate.
    real, parameter :: convbc = 1e-12/xmbc*Avog         ! BC
    !>>> TvN
    ! I have the impression that OC in M7 actually refers to POM.
    ! doc seems to be the density of POM not OC, and is equal to the density of BC (dbc).
    ! Compare with GLOMAP, where the densities of POM and BC are the same (Mann et al., GMD, 2010).
    ! This would imply that also the mass of POM and not OC should be used in M7,
    ! and that "OC" in the paper by Vignati et al. (JGR, 2004) is actually referring to POM,
    ! consistent with the terminology used by Stier et al. (ACP, 2005).
    !real, parameter :: convoc = 1e-12/xmpom*Avog*oc2pom ! OC (M7) POM (TM5)
    real, parameter :: convoc = 1e-12/xmpom*Avog ! POM (in both M7 and TM5)
    !<<< TvN
    real, parameter :: convss = 1e-12/xmnacl*Avog       ! SS
    real, parameter :: convdu = 1e-12/xmdust*Avog       ! DU
#ifdef with_budgets
    real,dimension(:,:,:), allocatable   :: processm7   ! M7 processes in M7 units.
    ! Array that convertts M7 processes from M7 to TM5 units.
    real, dimension(nm7procs), parameter :: processm7totm5 = (/ &
         !  1       2         3        4        5        6       7         8       9       10 
         convnumb,convsu  ,convnumb,convnumb,convsu  ,convnumb,convsu  ,convnumb,convsu,convnumb,&  ! 0?
         convsu  ,convnumb,convsu  ,convnumb,convsu  ,convnumb,convnumb,convsu  ,convbc,convoc  ,&  ! 1?
         convnumb,convsu  ,convbc  ,convoc  ,convnumb,convbc  ,convoc  ,convnumb,convsu,convbc  ,&  ! 2? 
         convoc  ,convnumb,convsu  ,convbc  ,convoc  ,convnumb,convnumb,convbc  ,convoc,convnumb,&  ! 3?
         convsu  ,convsu  ,convsu  ,convsu  ,convsu  ,convsu  ,convsu  ,convnumb,convbc,convoc  ,&  ! 4? 
         convnumb,convdu  ,convnumb,convdu  ,convnumb,convsu  ,convnumb,convsu  ,convbc,convoc  ,&  ! 5?
         convnumb,convsu  ,convbc  ,convoc  ,convss  ,convdu  /)
#endif
    real                 :: dryvol_m7, vol_nh4no3, vol_msa, vol_h2o
    real                 :: number_conc
    integer              :: imode
#endif /* M7 */

    real                 :: dtime, gmt, dxx, dyy, lat, lon
    real                 :: px, skf, x, month2dt
    integer              :: day, level, nzone, nzone_v
    integer,dimension(6) :: idater, idate_temp
    integer              :: i, j, l, n
    integer              :: ispec
    integer              :: imax, nca, nco, iloop
    integer              :: eqsam_option
    real                 :: dry_aerosol_mass, dry_aerosol_volume
    real                 :: water_aerosol_mass, water_aerosol_volume, ccs, s
    integer              :: i1, i2, j1, j2, lmr, is3, ie3, klm

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

    if ( okdebug ) write (*,*) 'start of chemistry'

    m   => m_dat   (region)%data
    rmc => chem_dat(region)%rm
    rm  => mass_dat(region)%rm
#ifdef slopes
    rxm => mass_dat(region)%rxm
    rym => mass_dat(region)%rym
    rzm => mass_dat(region)%rzm
#ifdef secmom
    rxxm => mass_dat(region)%rxxm
    rxym => mass_dat(region)%rxym
    rxzm => mass_dat(region)%rxzm
    ryym => mass_dat(region)%ryym
    ryzm => mass_dat(region)%ryzm
    rzzm => mass_dat(region)%rzzm
#endif
#endif

#ifdef with_m7
    dxyp => region_dat(region)%dxyp
#endif

    r_cloud => phot_dat(region)%cloud_reff

    ! ------- Time stuff
    call tau2date(itaur(region), idater)   ! current time
    dtime    = nchem/(2*tref(region))      ! time step
    month2dt = dtime/sec_month    ! conversion to emission per timestep

    call tau2date(itaur(region)+revert*nint(dtime*0.5),idate_temp)     ! get date @ halfway interval.
    gmt = idate_temp(4) + idate_temp(5)/60.0 + idate_temp(6)/3600.0    ! GMT in hours
    day = get_day(idate_temp(2),idate_temp(3),mlen)

    if ( okdebug ) then
       write(gol,*)'chemistry: region ',region, ' at date: ',idater,' with time step ',dtime
       call goPr
       write(gol,*)'chemistry: day, gmt, date@halfway interval ', day, gmt, idate_temp
       call goPr
    end if

    ! ------- Resolution and bounds
    dxx = dx/xref(region)         ! gridsize at current region
    dyy = dy/yref(region)         ! gridsize at current region

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


    is3=(i1-1)*3+1   ! to cover three times oversampled
    ie3=i2*3

    !    allocate(tha   (is3:ie3,j1:j2))
    allocate(tha   (im(region)*3,j1:j2))
    allocate(sza   (i1:i2, j1:j2              ))
    allocate(mu    (i1:i2, j1:j2              ))
    allocate(reff  (i1:i2, j1:j2              ))
    allocate(fact  (i1:i2, j1:j2, nbands_trop ))
    allocate(rj    (i1:i2, j1:j2, nj          ))
    allocate(rr    (i1:i2, j1:j2, nreac       ))
    allocate(rj_new(i1:i2, j1:j2, lmr, nj     )) 
    allocate(y     (i1:i2, j1:j2, maxtrace    ))
    allocate(y0    (i1:i2, j1:j2, maxtrace    ))
    allocate(yhelp (i1:i2, j1:j2, ntrace      ))
    allocate(ye    (i1:i2, j1:j2, n_extra     ))     ! Eqsam
    
    allocate(sad_cld_out(i1:i2, j1:j2, lmr )) ; sad_cld_out = 0.
    allocate(sad_ice_out(i1:i2, j1:j2, lmr )) ; sad_ice_out = 0.
    allocate(sad_aer_out(i1:i2, j1:j2, lmr )) ; sad_aer_out = 0.

    imax = (i2-i1+1)*(j2-j1+1) !nb of grid box per layer (was im(region)*jm(region))
    nca = 11
    nco = 37 ! was 36, JadB added one for the water associated with just the nitrate.
    allocate(yi(imax,nca)) ; yi = 0.0
    allocate(yo(imax,nco)) ; yo = 0.0

#ifdef with_m7
    allocate( presm7  (imax, 1         ))
    allocate( rhm7    (imax, 1         ))
    allocate( tempm7  (imax, 1         ))
    allocate( dryairm7(imax, 1         ))
    allocate( h2so4m7 (imax, 1         ))
    allocate( aemassm7(imax, 1, naermod))
    allocate( aenumm7 (imax, 1, nmod   ))
    allocate( aedensm7(imax,1,nmod))
    allocate( aewatm7 (imax,1,nmod))
    allocate( aeradm7 (imax,1,nmod))
    allocate( aeradrm7(imax,1,nsol))

    !PLS debug
    presm7  =0.
    rhm7    =0.
    tempm7  =0.
    dryairm7=0.
    h2so4m7 =0.
    aemassm7=0.
    aenumm7 =0.
    aedensm7=0.
    aewatm7 =0.
    aeradm7 =0.
    aeradrm7=0.
    
#ifdef with_budgets
    allocate(processm7(imax,1,nm7procs))
#endif
#endif
    
    tha  = 0.0
    sza  = 0.0
    mu   = 0.0
    fact = 0.0
    rj = 0.0
    rr = 0.0
    rj_new = 0.0
    y  = 0.0
    y0 = 0.0
    ye = 0.0

    ! alternative calculation of zenith angle with higher sampling (time and space). 
    ! yields factor of two lower angles in the tropics!
    ! Difference is resolution dependent. The coarser the larger the differences
    starttime = real(idate_temp(4)) + real(idate_temp(5))/60.0 + real(idate_temp(6))/3600.
    deltat = nchem/3600.0

    call daysim(region, day, is3, i1,i2,j1, j2, tha)
    call get_sza(region, i1, j1, i2, j2, starttime, deltat, tha, sza)

    ! ascribe mu 
    do j=j1,j2
       do i=i1,i2
          mu(i,j)=max(1.e-5,COS(sza(i,j)*pi/180.))
       enddo
    enddo
     
#ifdef with_optics
    call aerosol_info_m7(region, status)
    IF_NOTOK_RETURN(status=1)
#else 
    call aerosol_info(region)
#endif

    ! calculate actinic flux outside level loop for online fluxes
    call photo_flux(region, i1, j1, sza, rj_new)

    ! scale NOx BMB if needed
#ifndef without_emission
    call emission_nox_bb_daily_cycle(status)
    IF_NOTOK_RETURN(status=1)
#endif
    
    ! compute global LiNox budget
#ifdef with_budgets
#ifndef without_emission
#ifndef without_convection
    eminox=eminox+sum(eminox_lightning(region)%d3)*month2dt  ! kg [N] month-1 * dt /(sec_month) = kg [N]
#endif
#endif
#endif

    LEVS: do level = 1, lmr

       rj(:,:,:) = rj_new (:,:,level,:)
       reff(:,:) = r_cloud(:,:,level)

#ifdef with_budgets
       nzone_v=nzon_vg(level)
#endif

       ! get/calculate  all the fields needed in chemistry
       call GET_EXTRA( region, level, ye, i1, j1)

       ! calculate the rinc: relative radius growth of aerosols due to water

       ! Eqsam wants to calculate aerosol radii
       ! irinc is used by chem_rates in calchet2
       ! calchet2 is used by no one
       ! So, this part is not so important
       ! -----------------------------------
       ! To get this code back, take a stone-age copy of chemistry.F90
       ! But, to prevent even more errors.
       ! -----------------------------------
       ye(:,:,irinc) = 1.0

       call CALRATES(region, level, i1, j1, rj, rr, reff, ye, dtime, &
                     sad_cld_out(:,:,level), sad_ice_out(:,:,level), sad_aer_out(:,:,level), status) ! also returns rh in ye
       IF_ERROR_RETURN(status=1)


       ! CALCULATE NOX EMISSIONS
       ! ---------------------------
       do j=j1,j2
       do i=i1,i2

          ! nox emis in this cell; init to zero:
          x = 0.0   ! kg(N)/month

#ifndef without_emission
             
          ! add total 3d emissions if any
          if(associated(nox_emis_3d(region)%d3)) then
             x = x + nox_emis_3d(region)%d3(i,j,level)
          end if

          if(associated(nox_emis_3d_bb_app(region)%d3)) then
             x = x + nox_emis_3d_bb_app(region)%d3(i,j,level)
          end if

#ifndef without_convection
          ! add lightning emissions:
          x = x + eminox_lightning(region)%d3(i,j,level)
#endif
#endif /* EMISS */

          ! emission budget (since nox emissions are added in chemistry)
#ifdef with_budgets
#ifndef without_emission
          budemi_dat(region)%emi(i,j,nzone_v,inox) = &
                  budemi_dat(region)%emi(i,j,nzone_v,inox) + x/xmn*month2dt*1e3   ! g
          if (inox == 1) sum_emission(region) = sum_emission(region) + x*month2dt  ! kg N
#endif
#endif
          ! put emissions in the units #/s/cm3 used in chemistry
          ye(i,j,ieno) = x/ye(i,j,iairm)*ye(i,j,iairn)* &
                  (xmair/xmn)/sec_month   !kg N/month ----> #(NO)/cm3/s
       end do
       end do


       ! INITIAL CONCENTRATIONS
       ! ---------------------------
       do n=1,ntracet
       do j=j1,j2
       do i=i1,i2
          y(i,j,n) = rm(i,j,level,n)/ye(i,j,iairm)*ye(i,j,iairn)* fscale(n)   !kg ----> #/cm3
       end do
       end do
       end do

       do n=ntracet+1,ntrace
       do j=j1,j2
       do i=i1,i2
          ! remember nontransported tracers are in rmc
          y(i,j,n) =rmc(i,j,level,n)/ye(i,j,iairm)*ye(i,j,iairn)* fscale(n)   !kg ----> #/cm3
       end do
       end do
       end do

       ! FILL EXTRA SPECIES : (from ntrace+1 to maxtrace)
       ! ---------------------------
       do j=j1,j2
       do i=i1,i2
          ! these are in #/cm3 and used in budget....should be in y....
          y(i,j,iair) = ye(i,j,iairn)
          y(i,j,ih2o) = ye(i,j,ih2on)
          
          ! introduced to calculate correct budget for JO2
          y(i,j,io2) = 1.

          ! precalculate as first order N2O5 loss.
          rr(i,j,kn2o5aq) = 0.

       end do
       end do

       ! BACKUP INIT CONCENTRATIONS
       ! ---------------------------
       y0 = y

       do klm=1,ntrace
          yhelp(:,:,klm) = y(:,:,klm)*ye(:,:,iairm)*1000./xmair/y(:,:,iair)
       end do

       ! EBI
       ! ---------------------------

       CALL EBI( region, level, i1, j1, rj, rr, y, ye, status)
       IF_ERROR_RETURN(status=1)
       
       ! EBI scheme budgets
#ifdef with_budgets
       do j=j1,j2
       do i=i1,i2
          nzone=budg_dat(region)%nzong(i,j)    ! global budget zone
          budchem(nzone,nzone_v,1:ntrace) = budchem(nzone,nzone_v,1:ntrace) + &
               y(i,j,1:ntrace)*ye(i,j,iairm)*1000./xmair/y(i,j,iair)-yhelp(i,j,:)
       end do
       end do
#endif

       !--- Begin of M7: ------------------------------------------------------
       
#ifdef with_m7
       ! M7
       iloop = 0
       do j=j1,j2
          do i=i1,i2

#ifdef with_budgets
             ! For the M7 budgets, we first subtract the old concentrations.
             nzone=budg_dat(region)%nzong(i,j)    ! global budget zone
             budm7(nzone,nzone_v,:)=budm7(nzone,nzone_v,:)-y(i,j,1:ntracet)*ye(i,j,iairm)*1000./xmair/y(i,j,iair)
#endif

             iloop = iloop + 1

             ! meteo parameters
             presm7  (iloop,1) = ye(i,j,i_pres)                             ! Pa
             rhm7    (iloop,1) = rh(region)%d3(i,j,level)           ! fraction, in cloud free part!
             tempm7  (iloop,1) = ye(i,j,i_temp)                             ! K
             dryairm7(iloop,1) = presm7(iloop,1)/(Rgas_air*tempm7(iloop,1)) ! density of dry air

             ! H2SO4 gas concentrations  molecule/cm3
             h2so4m7(iloop,1) = y(i,j,iso4)/convsu
             ! molecule/cm3 ----> molecules/cm3

             ! sulphate mass in the 4 soluble modes molecule/cm3
             aemassm7(iloop,1,1) = y(i,j,iso4nus)/convsu
             aemassm7(iloop,1,2) = y(i,j,iso4ais)/convsu
             aemassm7(iloop,1,3) = y(i,j,iso4acs)/convsu
             aemassm7(iloop,1,4) = y(i,j,iso4cos)/convsu
             ! molecule/cm3 ----> molecules/cm3

             ! BC mass in the 3 soluble modes and 1 insoluble
             aemassm7(iloop,1,5) = y(i,j,ibcais)/convbc
             aemassm7(iloop,1,6) = y(i,j,ibcacs)/convbc
             aemassm7(iloop,1,7) = y(i,j,ibccos)/convbc
             aemassm7(iloop,1,8) = y(i,j,ibcaii)/convbc
             ! molecule/cm3 -----> ug/m3

             ! was: POM (tm5) OC (m7) mass in the 3 soluble modes and 1 insoluble
             ! POM mass in the 3 soluble modes and 1 insoluble
             aemassm7(iloop,1,9)  = y(i,j,ipomais)/convoc
             aemassm7(iloop,1,10) = y(i,j,ipomacs)/convoc
             aemassm7(iloop,1,11) = y(i,j,ipomcos)/convoc
             aemassm7(iloop,1,12) = y(i,j,ipomaii)/convoc
             ! molecule/cm3 -----> ug/m3


             ! SS mass in the 2 soluble modes
             aemassm7(iloop,1,13) = y(i,j,issacs)/convss
             aemassm7(iloop,1,14) = y(i,j,isscos)/convss
             ! molecule/cm3 -----> ug/m3

             ! DUST mass in the 2 soluble modes and 2 insoluble
             aemassm7(iloop,1,15) = y(i,j,iduacs)/convdu
             aemassm7(iloop,1,16) = y(i,j,iducos)/convdu
             aemassm7(iloop,1,17) = y(i,j,iduaci)/convdu
             aemassm7(iloop,1,18) = y(i,j,iducoi)/convdu
             ! molecule/cm3 -----> ug/m3

             ! number concentrations of the 7 modes
             aenumm7(iloop,1,1) = y(i,j,inus_n)/convnumb
             aenumm7(iloop,1,2) = y(i,j,iais_n)/convnumb
             aenumm7(iloop,1,3) = y(i,j,iacs_n)/convnumb
             aenumm7(iloop,1,4) = y(i,j,icos_n)/convnumb
             aenumm7(iloop,1,5) = y(i,j,iaii_n)/convnumb
             aenumm7(iloop,1,6) = y(i,j,iaci_n)/convnumb
             aenumm7(iloop,1,7) = y(i,j,icoi_n)/convnumb

             !
             ! molecule/cm3 -----> particle/cm3:
             !                                                    1      # air
             !rm in total aerosol number (#aer) ---->     #aer * ------ * ------   is y
             !                                                   kg(air)   cm-3
             !#aer   # air                                           1
             !---- * ------- ----> 1 kg(air) = 10^3 g air = 10^3 * ----- * Avog
             !cm-3   kg(air)                                       xmair
             !

          end do
       end do

       ! The M7 processes will be calculated when budgets are used. The extra
       ! argument will appear when the compiler flag with_budgets is set.
#ifdef with_budgets
       
       call M7(iloop, imax,      1,                             &    ! TM5 indices
                      presm7,    rhm7,      tempm7,             &    ! "   thermodynamics
                      h2so4m7,   aemassm7,  aenumm7,            &    ! M7   tracers
                      aedensm7,  aewatm7,   aeradm7,  aeradrm7, dtime, &     ! M7   aerosol propertis
                      processm7  ) ! The extra argument, in which the rates of each process is stored in M7 units.

#else
       call M7(iloop, imax,      1,                             &    ! TM5 indices
                      presm7,    rhm7,      tempm7,             &    ! "   thermodynamics
                      h2so4m7,   aemassm7,  aenumm7,            &    ! M7   tracers
                      aedensm7,  aewatm7,   aeradm7,  aeradrm7, dtime)      ! M7   aerosol propertis
#endif


       iloop = 0
       do j=j1,j2
          do i=i1,i2
             iloop = iloop + 1

             ! H2SO4 gas concentrations in molecule/cm3
             y(i,j,iso4) = h2so4m7(iloop,1)*convsu
             !  molecule/cm3 ----> molecule/cm3

             ! sulphate mass in the 4 soluble modes
             y(i,j,iso4nus) = aemassm7(iloop,1,1)*convsu
             y(i,j,iso4ais) = aemassm7(iloop,1,2)*convsu
             y(i,j,iso4acs) = aemassm7(iloop,1,3)*convsu
             y(i,j,iso4cos) = aemassm7(iloop,1,4)*convsu
             !  molecule/cm3 ----> molecule/cm3

             ! BC mass in the 3 soluble modes and 1 insoluble
             y(i,j,ibcais) = aemassm7(iloop,1,5)*convbc
             y(i,j,ibcacs) = aemassm7(iloop,1,6)*convbc
             y(i,j,ibccos) = aemassm7(iloop,1,7)*convbc
             y(i,j,ibcaii) = aemassm7(iloop,1,8)*convbc
             !  ug/m3 ----> molecule/cm3


             ! was: POM (tm5) OC (m7) mass in the 3 soluble modes and 1 insoluble
             ! POM mass in the 3 soluble modes and 1 insoluble
             y(i,j,ipomais) = aemassm7(iloop,1,9)*convoc
             y(i,j,ipomacs) = aemassm7(iloop,1,10)*convoc
             y(i,j,ipomcos) = aemassm7(iloop,1,11)*convoc
             y(i,j,ipomaii) = aemassm7(iloop,1,12)*convoc
             !  ug/m3 ----> molecule/cm3

             ! SS mass in the 2 soluble modes
             y(i,j,issacs) = aemassm7(iloop,1,13)*convss
             y(i,j,isscos) = aemassm7(iloop,1,14)*convss
             !  ug/m3 ----> molecule/cm3

             ! DUST mass in the 2 soluble modes and 2 insoluble
             y(i,j,iduacs) = aemassm7(iloop,1,15)*convdu
             y(i,j,iducos) = aemassm7(iloop,1,16)*convdu
             y(i,j,iduaci) = aemassm7(iloop,1,17)*convdu
             y(i,j,iducoi) = aemassm7(iloop,1,18)*convdu
             !  ug/m3 ----> molecule/cm3

             ! number concentrations of the 7 modes
             y(i,j,inus_n) = aenumm7(iloop,1,1)*convnumb
             y(i,j,iais_n) = aenumm7(iloop,1,2)*convnumb
             y(i,j,iacs_n) = aenumm7(iloop,1,3)*convnumb
             y(i,j,icos_n) = aenumm7(iloop,1,4)*convnumb
             y(i,j,iaii_n) = aenumm7(iloop,1,5)*convnumb
             y(i,j,iaci_n) = aenumm7(iloop,1,6)*convnumb
             y(i,j,icoi_n) = aenumm7(iloop,1,7)*convnumb
             !number of particles/cm3  ---->  molecule/cm3 (see comment above)

#ifdef with_budgets
             ! For the M7 budgets, we now add the new concentrations.
             nzone=budg_dat(region)%nzong(i,j)    ! global budget zone
             budm7(nzone,nzone_v,:)=budm7(nzone,nzone_v,:)+y(i,j,1:ntracet)*ye(i,j,iairm)*1000./xmair/y(i,j,iair)
             budm7proc(nzone,nzone_v,:)=budm7proc(nzone,nzone_v,:)+processm7(iloop,1,:)*processm7totm5*ye(i,j,iairm)*&
                  1000./xmair/y(i,j,iair) ! processm7totm5 is a 66-element array with unit conversion factors.
#endif

             ! output m7
             ! kgH2O/m3  ---> kg water per gridbox
             do imode=1,nsol
                h2o_mode(region,imode)%d3(i,j,level) = aewatm7(iloop,1,imode) *dxyp(j)*ye(i,j,idz)
                rwd_mode (region,imode)%d3(i,j,level) = aeradrm7(iloop,1,imode)*1e-2   ! from cm to meter
             end Do

             do imode=1,nmod
                dens_mode(region,imode)%d3(i,j,level) = aedensm7(iloop,1,imode)*1e3   ! from g/cm^3 --> kg/m^3
                rw_mode  (region,imode)%d3(i,j,level) = aeradm7(iloop,1,imode)*1e-2   ! from cm to meter
             End Do
          end do
       end do
#endif /* M7 */

       ! ====================== EQSAM        ! Using HNO3 / NO3_A, NH4 / NH2 and SO4.
       yi = 0.0
       yo = 0.0

       iloop = 0
       do j=j1,j2
       do i=i1,i2
             iloop = iloop + 1
             yi(iloop,1)=ye(i,j,i_temp)     ! K
#ifndef without_sedimentation
#ifdef with_m7
             s = rh(region)%d3(i,j,level) ! in cloud free part!
#else
             s = 0.0
#endif
#else
             s = 0.0
#endif
             yi(iloop,2)=max(0.001,min(0.95,s))
             yi(iloop,11)=ye(i,j,i_pres)/100.                 ! Pa -> hPa
             ! yi(iloop,6)= c(nv2,ina)/Avog*1.e12
             yi(iloop,6)=0.
#ifdef with_m7
             yi(iloop,4)=(y(i,j,iso4)+y(i,j,iso4nus)+y(i,j,iso4ais)+y(i,j,iso4acs)+y(i,j,iso4cos))/Avog*1.e12      ! molec/cm3 -> umol/m3
#else
             yi(iloop,4)=(y(i,j,iso4))/Avog*1.e12  ! molec/cm3 -> umol/m3
#endif
             yi(iloop,3)=(y(i,j,inh3)+y(i,j,inh4))/Avog*1.e12
             yi(iloop,5)=(y(i,j,ihno3)+y(i,j,ino3_a))/Avog*1.e12
             !yi(iloop,7)=(c(nv2,ihcl)+c(nv2,icl))/Avog*1.e12
             yi(iloop,7)=0.
             !yi(iloop,8)=c(nv2,ik)/Avog*1.e12
             yi(iloop,8)=0.
             !yi(iloop,9)=c(nv2,ica)/Avog*1.e12
             yi(iloop,9)=0.
             !yi(iloop,10)=c(nv2,img)/Avog*1.e12
             yi(iloop,10)=0.
       end do
       end do

       eqsam_option = 1
       call EQSAM_V03D( yi, yo, nca, nco, eqsam_option, iloop, imax, 0, (/0,0,0,0,0,0/))
       iloop = 0

       do j=j1,j2
          do i=i1,i2
             iloop = iloop + 1
#ifdef with_budgets
             nzone=budg_dat(region)%nzong(i,j)    ! global budget zone
             budeqsam(nzone,nzone_v,1) = budeqsam(nzone,nzone_v,1) + &
                                         (yo(iloop,19)*Avog*1.e-12-y(i,j,inh4))  *ye(i,j,iairm)*1000./xmair/y(i,j,iair) ! Ammonium
             budeqsam(nzone,nzone_v,2) = budeqsam(nzone,nzone_v,2) + &
                                         (yo(iloop,20)*Avog*1.e-12-y(i,j,ino3_a))*ye(i,j,iairm)*1000./xmair/y(i,j,iair) ! Nitrate
#endif
             y(i,j,ihno3)  = yo(iloop, 9)*Avog*1.e-12 ! umol/m3 -> molec/cm3
             y(i,j,inh3)   = yo(iloop,10)*Avog*1.e-12 ! umol/m3 -> molec/cm3
             y(i,j,inh4)   = yo(iloop,19)*Avog*1.e-12 ! umol/m3 -> molec/cm3
             y(i,j,ino3_a) = yo(iloop,20)*Avog*1.e-12 ! umol/m3 -> molec/cm3

             ! Sulfuric acid is unchanged: 
             ! y(i,j,iso4)   = yo(iloop,21)*Avog*1.e-12

             ! The NO3_A should also get some water. M7 takes the water for the
             ! anoins: SO4-- and Cl- (Sea Salt). Cations are neglected.
             ! Just the NO3- water is used apart from the M7 water.

             ! Ignore the direct eqsam water

             ! ug/m3 -> molec/cm3
             ! y(i,j,ih2opart) = yo(iloop,12)*Avog*1.e-12/xmh2o

#ifdef with_m7
             ! The water associated with nitrate aerosol is put into 
             ! the soluble accumulation mode.
             h2o_mode(region,3)%d3(i,j,level) = h2o_mode(region,3)%d3(i,j,level) + yo(iloop,37)*dxyp(j)*ye(i,j,idz)*1.e-9 
             ! from ug/m^3 (according to eqsam), then take the M7 conversion (from kg/m^3) times 1.e-9.
             ! This EQSAM water is added to the optics water at the call of optics.
#endif
          end do
       end do

       ! END EQSAM


       ! Check for negative concentrations.
       do n=1,ntrace
       do j=j1,j2
       do i=i1,i2
                if ( y(i,j,n) < 0. ) then
                   !print *,' chemistry: negatives appeared in chemie'
                   !write(6,*)i,j,level,idater
                   !write(6,907)(y0(i,j,ispec),ispec=1,ntrace)
                   !write(6,907)(y(i,j,ispec),ispec=1,ntrace)
                   if ( n > ntracet ) then
                      do ispec=ntracet+1,ntrace
                         y(i,j,ispec)=max(0.,y(i,j,ispec))
                      end do
                   else
                      write (gol,'(a,a8,a,4i4,e12.4)') &
                           'chemistry: negatives for species:',names(n), ' at ', i,j,level,region, &
                           y(i,j,n)/y(i,j,iair); call goPr
                      y(i,j,n)=max(0.,y(i,j,n))
                   end if
                end if
       end do !i
       end do !j
       end do !n

       ! Write the tracers back to the mass_dat (rm)
       do n=1,ntrace
#ifdef with_cariolle
          if ( (n==io3) .and. ( (level.ge.l_cariolle(region)) .or. with_trop_force ) ) then
             ! O3 is changed by Cariolle chemistry
          else
#endif
             do j=j1,j2
                do i=i1,i2
                   if ( n <= ntracet ) then
                      rm(i,j,level,n)=y(i,j,n)/y(i,j,iair)*ye(i,j,iairm)/fscale(n)
                   else
                      rmc(i,j,level,n)=y(i,j,n)/y(i,j,iair)*ye(i,j,iairm)/fscale(n)
                   end if
#ifdef slopes
                   if ( n <= ntracet ) then
                      skf=0.
                      if ( y(i,j,n) > eps .and. y0(i,j,n) > eps ) skf=y(i,j,n)/y0(i,j,n)
                      rxm(i,j,level,n)=rxm(i,j,level,n)*skf
                      rym(i,j,level,n)=rym(i,j,level,n)*skf
                      rzm(i,j,level,n)=rzm(i,j,level,n)*skf
#ifdef secmom
                      rxxm(i,j,level,n)=rxxm(i,j,level,n)*skf
                      rxym(i,j,level,n)=rxym(i,j,level,n)*skf
                      rxzm(i,j,level,n)=rxzm(i,j,level,n)*skf
                      ryym(i,j,level,n)=ryym(i,j,level,n)*skf
                      ryzm(i,j,level,n)=ryzm(i,j,level,n)*skf
                      rzzm(i,j,level,n)=rzzm(i,j,level,n)*skf
#endif
                   endif
#endif
                end do
             end do   !j,i
#ifdef with_cariolle
          endif
#endif
       end do   !n

! For the optics calculations,
! it is assumed that nitrate aerosol is formed
! by condensation onto existing particles
! in the soluble accumulation mode (Adams et al., JGR, 2001).
! Therefore, the presence of nitrate does not increase
! the number of particles, nor does it have an impact
! on the processes described by M7.
! Above, the water associated with nitrate aerosol 
! is added to the total aerosol water (h2o_mode) in the accumulation mode.
! However, this only affects the calculation of the refractive index, 
! which also accounts for the nitrate mass itself.
! In order to properly account for the additional ammonium nitrate 
! and the associated water in the optics routine,
! the corresponding dry and wet particle radii have to be increased as well;
! the aerosol density (dens_mode) can be left unchanged, 
! because it is not used in the optics routine.
! We use the same refractive index for ammonium nitrate as for sulfate.
! This is a reasonable approximation, since,
! according to Tang (JGR, 1997), for sulfate and nitrate aerosols
! the effect of chemical composition on light scattering
! is outweighted by the size effect of the particles.
! Lowenthal et al. (Atmos. Environ., 2000)
! give a refractive index of 1.44 for sulfuric acid
! and 1.55 for ammonium nitrate.
! If desired, the refractive index of ammonium nitrate
! could be estimated at different wavelengths 
! by scaling the OPAC data for sulfate,
! based on these values.

! In M7 and in the optics calculations,
! sulfate aerosols are assumed to be pure sulfuric acid,
! and we do not try to account for the presence of ammonium (bi)sulfate.
! The reason is that the presence of ammonium (bi)sulfate
! has competing effects on scattering,
! which cannot be easily accounted for in a model based on M7.
! Sulfuric acid particles are extremely hygroscopic 
! and will draw significant water mass into the aerosol phase 
! at any relative humidity (RH). 
! If these particles are partially or completely neutralized 
! by drawing ammonia from the gas phase, 
! there will be an increase in particle mass due to the added ammonium
! but a decrease in particle hygroscopicity at low to moderate RH, 
! and thus a decrease in particulate water mass. 
! On the other hand, ammonium bisulfate and ammonium sulfate both 
! have a higher refractive index than sulfuric acid (see e.g. Boucher and Anderson, 1995). 
! The net result of these competing factors is that 
! one mole of sulfuric acid scatters about 25% more sunlight 
! than one mole of ammonium bisulfate at 80% relative humidity. 
! Ammonium sulfate is intermediate between the two. 
! For more information see Boucher and Anderson (1995) and Boucher et al. (JGR, 1998). 
! An excellent discussion of the optical properties of ammonium (bi)sulfate 
! versus sulfuric acid is also given by Adams et al. (JGR, 2001).
!
! In remote regions, especially outside the tropics,
! the contribution of methane sulfonate (MSA-) aerosol 
! to the total particulate sulfure burden
! cannot be neglected.
! Like nitrate, MSA- aerosols are mainly formed by
! condensation onto exisiting particles,
! in the soluble accumulation and coarse modes 
! (Saltzmann et al., JGR, 1983; Pzenny et al., J. Atmos. Chem., 1992).
! For simplicity, we assume in the model
! that MSA- is formed only in the accumulation mode.
! Moreover, it seems reasonable to use the refractive index of
! sulfuric acid (H2SO4) for MSA, as is done by Kinne et al. (2003).
! The similarity between the refractive index for MSA and H2SO4
! is also discussed by Myhre et al. (Applied Optics, 2004).
!
#ifdef with_m7
       do j=j1,j2
          do i=i1,i2
             ! prevent division for zero for extremely low aerosol concentrations
             ! calculate corresponding number concentration in #/cm3
             ! and test if this is larger than 1.e-10
             ! the same criterion is applied in m7_equil
             number_conc = rm(i,j,level,iacs_n) / ( 1.e6 * dxyp(j) * ye(i,j,idz) )
             if ( number_conc > 1.e-10 ) then
               dryvol_m7=((rwd_mode(region,3)%d3(i,j,level)*cmr2ram(3))**3.0)*pi/0.75 ! dry particle volume [m3]
                                                                                    ! in the accumulation mode from M7
               vol_nh4no3=rm(i,j,level,ino3_a)*nh4no3_factor/ &
                         (rm(i,j,level,iacs_n)*nh4no3_density)                        ! ammonium-nitrate aerosol volume per particle [m3]
               vol_msa=rm(i,j,level,imsa)/(rm(i,j,level,iacs_n)*msa_density)          ! MSA- aerosol volume per particle [m3]
               vol_h2o=h2o_mode(region,3)%d3(i,j,level)/(rm(i,j,level,iacs_n)*rol)    ! total aerosol water per particle [m3]
               rw_mode(region,3)%d3(i,j,level) = ((dryvol_m7+vol_nh4no3+vol_msa+vol_h2o)*0.75/pi)**(1./3.)*ram2cmr(3)
               rwd_mode(region,3)%d3(i,j,level) = ((dryvol_m7+vol_nh4no3+vol_msa)*0.75/pi)**(1./3.)*ram2cmr(3)
             else
               ! don't change rwd_mode and set rw_mode equal to rwd_mode
               rw_mode(region,3)%d3(i,j,level) = rwd_mode(region,3)%d3(i,j,level)
             endif
          enddo
       enddo
#endif

    end do LEVS ! loop over vertical levels...

    ! SAD
    phot_dat(region)%sad_cld=sad_cld_out
    phot_dat(region)%sad_ice=sad_ice_out
    phot_dat(region)%sad_aer=sad_aer_out
    ! perform averaging...
    phot_dat(region)%nsad_av    = phot_dat(region)%nsad_av     + float(ndyn)/float(ndyn_max)
    phot_dat(region)%sad_cld_av = phot_dat(region)%sad_cld_av  + (float(ndyn)/float(ndyn_max))*phot_dat(region)%sad_cld 
    phot_dat(region)%sad_ice_av = phot_dat(region)%sad_ice_av  + (float(ndyn)/float(ndyn_max))*phot_dat(region)%sad_ice 
    phot_dat(region)%sad_aer_av = phot_dat(region)%sad_aer_av  + (float(ndyn)/float(ndyn_max))*phot_dat(region)%sad_aer 

    ! cleanup
    deallocate(tha)
    deallocate(mu)
    deallocate(sza)
    deallocate(reff)
    deallocate(fact)
    deallocate(rj)
    deallocate(rj_new)
    deallocate(rr)
    deallocate(y)
    deallocate(y0)
    deallocate(ye, yhelp)
    ! eqsam
    deallocate(yo)
    deallocate(yi)
    ! sad
    deallocate(sad_cld_out,sad_ice_out,sad_aer_out)

#ifdef with_m7
    deallocate(presm7)
    deallocate(rhm7)
    deallocate(tempm7)
    deallocate(dryairm7)
    deallocate(h2so4m7)
    deallocate(aemassm7)
    deallocate(aenumm7)
    deallocate(aedensm7)
    deallocate(aewatm7)
    deallocate(aeradm7)
    deallocate(aeradrm7)
#ifdef with_budgets
    deallocate(processm7)
#endif
    nullify(dxyp)
#endif

    nullify(m)
    nullify(rm)
#ifdef slopes
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
#ifdef secmom
    nullify(rxxm)
    nullify(rxym)
    nullify(rxzm)
    nullify(ryym)
    nullify(ryzm)
    nullify(rzzm)
#endif
#endif
    nullify(r_cloud)


    if ( okdebug ) write (*,*) 'end of chemistry'

    ! ok
    status = 0

  END SUBROUTINE CHEMIE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    GET_EXTRA
  !
  ! !DESCRIPTION: 
  ! calculate / get 2D fields needed in the chemistry routines...
  ! some fields (like pH) are calculated later on...
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE GET_EXTRA( region, level, ye, is, js)
    !
    ! !USES:
    !
    use Binas,         only : Avog, grav
    use dims,          only : im, jm, lm
    use Dims,          only : CheckShape
    use global_data,   only : region_dat, mass_dat
    use meteodata,     only : phlb_dat, m_dat, temper_dat, humid_dat, gph_dat, cc_dat, iwc_dat, lwc_dat
    use chem_param,    only : n_extra, i_pres, i_temp, iairn, ih2on
    use chem_param,    only : xmair, xmh2o, iairm, ilwc, iiwc, icc, iph, idz,iciwc,iclwc
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)                     :: region         ! region #
    integer, intent(in)                     :: level, is, js  ! level and start indices
    !
    ! !OUTPUT PARAMETERS:
    !
    real, intent(out)                       :: ye(is:,js:,:)
    !
    ! !REVISION HISTORY: 
    !   22 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real,dimension(:,:,:),pointer          :: phlb,m,gph,t,q,lwc,iwc,cc

    integer :: i,j,l
    real    :: px,x1,x2,xice,xliq,aird

    integer ::  imr, jmr, lmr, i1, i2, j1, j2, status

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

    call Get_DistGrid( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

    ! check arguments ...
    call CheckShape( (/i2-i1+1, j2-j1+1, n_extra/), shape(ye ), status )

    lmr = lm(region)
    phlb =>   phlb_dat(region)%data
    m    =>      m_dat(region)%data
    t    => temper_dat(region)%data
    q    =>  humid_dat(region)%data
    gph  =>    gph_dat(region)%data ! (:,:,1:lm(region)+1)
    lwc  =>    lwc_dat(region)%data
    iwc  =>    iwc_dat(region)%data
    cc   =>     cc_dat(region)%data

    do j=j1,j2
    do i=i1,i2
          px = 0.5*(phlb(i,j,level)+phlb(i,j,level+1))
          ye(i,j,i_pres) = px
          ye(i,j,i_temp) = t(i,j,level)
          ye(i,j,iairn ) = 7.24e16*px/t(i,j,level)
          ye(i,j,ih2on ) = q(i,j,level)*ye(i,j,iairn)*xmair/xmh2o
          ye(i,j,iairm ) = m(i,j,level)
          ye(i,j,iclwc ) = lwc(i,j,level)
          ye(i,j,iciwc ) = iwc(i,j,level)
          ! x1: kg water to m3 water   (m3 water/ kg air)
          x1=lwc(i,j,level)*1.e-3
          aird = ye(i,j,iairn)

          x2=xmair*1e3*aird/Avog           ! kg/m3 (air)
          xliq = x1/x2                     ! dimensionless number (m^3/m^3)
          ! avoid negatives and artificial values(1e-10 is about 0.0001 g/m3)
          if ( xliq < 1e-10 ) xliq=0.
          ye(i,j,ilwc) = xliq

          x1=iwc(i,j,level)*1.e-3      ! kg ice to m3 ice
          xice=x1/x2
          ! avoid negatives and artificial values)
          if ( xice < 1e-10 ) xice=0.
          ye(i,j,iiwc) = xice
          ye(i,j,icc)  = cc(i,j,level)
          ye(i,j,iph)  = 0.0 ! only set in wetS when cloudy!!!
          ye(i,j,idz)  = gph(i,j,level+1)-gph(i,j,level)  !dz
       end do
    end do

    nullify(phlb)
    nullify(m)
    nullify(gph)
    nullify(t)
    nullify(q)
    nullify(lwc)
    nullify(iwc)
    nullify(cc)

  END SUBROUTINE GET_EXTRA
  !EOC

  
END MODULE CHEMISTRY