!### macro's #####################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!-------------------------------------------------------------------------
!                    TM5                                                 !
!-------------------------------------------------------------------------
!BOP
!
! !MODULE:  toolbox
!
! !DESCRIPTION: General tools for chemistry/emission treatment. 
!\\
!\\
! !INTERFACE: 
!
module toolbox
  !
  ! !USES:
  !
  use GO,   only : gol, goErr, goPr, goLabel
  use dims, only : okdebug
  
  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: zfarr
  public :: ltropo, lvlpress, print_totalmass
  public :: dumpfield, dumpfieldi
  public :: coarsen_emission, coarsen_emission_2d
  public :: escape_tm
  public :: distribute_emis2D
  public :: distribute1x1
  public :: distribute1x1b
  public :: tropospheric_columns
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter     ::  mname = 'toolbox'
  !
  ! !INTERFACE:
  !
  interface coarsen_emission
     module procedure coarsen_emission_1d
     module procedure coarsen_emission_2d
     module procedure coarsen_emission_3d
     module procedure coarsen_emission_d23
  end interface
  !
  ! !REVISION HISTORY: 
  !   16 Jul 2010 - Achim Strunk - optional input to distribute_emis2D
  !
  ! !REMARKS:
  !
  !EOP
  !-----------------------------------------------------------------------

contains


  ! ============================================================================


  !-----------------------------------------------------------------------
  !
  !****	tropo determine tropopause level
  !
  !	programmed by		 fd 01-11-1996 
  !       changed to function      fd 12-07-2002
  !
  !	purpose
  !	-------
  !	
  !	interface
  !	---------
  !	i=ltropo(region,tx,gphx,lmx)
  !
  !	method
  !	------
  !	determine tropopause level
  !       from temperature gradient
  !
  !	externals
  !	---------
  !       function lvlpress
  !
  !	reference
  !	---------	
  !	WMO /Bram Bregman
  !----------------------------------------------------------------------

  integer function ltropo(region,tx,gphx,lmx)

    ! in/out
    integer,intent(in)               :: region
    integer,intent(in)               :: lmx
    real,dimension(lmx),intent(in)   :: tx
    real,dimension(lmx+1),intent(in) :: gphx 

    ! local
    integer :: ltmin,ltmax,l
    real    :: dz,dt

    !   ltropo is highest tropospheric level
    !   above is defined as stratosphere

    ! min tropopause level 450 hPa (ca. 8 km)
    ltmin=lvlpress(region,45000.,98400.)
    ! max tropopause level 70 hPa (ca. 20 km)
    ltmax=lvlpress(region,7000.,98400.)
    ltropo=ltmin
    do l=ltmin,ltmax
       dz=(gphx(l)-gphx(l-2))/2. 
       dt=(tx(l)-tx(l-1))/dz
       if ( dt < -0.002) ltropo=l !wmo 85 criterium for stratosphere
       ! increase upper tropospheric level
    end do !l

  end function ltropo


  !-----------------------------------------------------------------------
  !
  !****   lvlpress determines the index of the level with a pressure
  !       greater than press i.e. in height below press
  !
  !       programmed by           Peter van Velthoven, 6-november-2000
  !
  !       purpose
  !       -------
  !       determines level independent of vertical resolution lm 
  !
  !       interface
  !       ---------
  !       i = lvlpress(press,press0)
  !
  !       method
  !       ------
  !       use hybrid level coefficients to determine lvlpress
  !
  !       input
  !       ---------
  !       press  : pressure in Pa
  !       press0 : surface pressure in Pa
  !
  !----------------------------------------------------------------------

  integer function lvlpress(region,press,press0)

    use dims, only: at, bt, lm

    ! in/out
    integer, intent(in) :: region
    real, intent(in)    :: press
    real, intent(in)    :: press0

    ! local
    real             :: zpress0, zpress
    integer          :: l,lvl

    if ( press0 == 0.0 ) then
       zpress0 = 98400.
    else
       zpress0 = press0
    endif
    lvl = 1
    ! l increases from bottom (l=1) to top (l=lm)
    do l = 1, lm(region)
       zpress =  (at(l)+at(l+1) + (bt(l)+bt(l+1))*zpress0)*0.5
       if ( press < zpress ) then
          lvl = l
       endif
    end do
    lvlpress = lvl

  end function lvlpress



!  subroutine calc_region_coordinates(region,xlat,xlon)
!    !
!    !
!    !
!    use dims, only: im, jm, dx, dy, xref, yref, xbeg, ybeg
!    implicit none
!
!    ! in/out
!    integer, intent(in)                       :: region
!    real, dimension(im(region)+1),intent(out) :: xlon
!    real, dimension(jm(region)+1),intent(out) :: xlat
!
!    ! local
!    integer :: i,j,l,n
!    real    :: dx1,dy1   ! gridsize in degrees longitude and latitude
!
!    dx1=dx/xref(region) !longitude
!    dy1=dy/yref(region) !latitude
!
!
!    do j=1,jm(region)+1
!       ! southern border,last value is North Border
!       xlat(j)=ybeg(region)+dy1*(j-1)
!    enddo
!
!    do i=1,im(region)+1
!       ! western border, last value is East Border
!       xlon(i)=xbeg(region)+dx1*(i-1)
!    enddo
!
!  end subroutine calc_region_coordinates



  real function zfarr(rx1,er,ztrec) 
    !------------------------------------------------------------------
    !     
    !****  ZFARR calculation of Arrhenius expression for rate constants
    !
    !------------------------------------------------------------------
    !
    implicit none
    real,intent(in) :: rx1,er,ztrec            
    !
    zfarr=rx1*exp(er*ztrec)
    !
  end function zfarr



  subroutine print_totalmass(region)

    use dims,        only : im,jm,lm
    use global_data, only : mass_dat

    implicit none

    ! in/out
    integer,intent(in)            :: region

    ! local
    integer       :: i
    real          :: value
    logical,save  :: first=.true.  

    value = &
         sum(mass_dat(region)%rm_t(1:im(region),1:jm(region),1:lm(region),1))

    if ( first ) then
       first = .false.
       write (gol,*) 'print_totalmass: Initial Total mass tracer 1 ',value; call goPr 
    else
       write (gol,*) 'print_totalmass: Total mass tracer 1 ',value; call goPr 
    end if

  end subroutine print_totalmass



  subroutine dumpfield(flag,fname,im,jm,lm,nt,field,namfield)
    !
    ! Write 4D field (type real) to HDF file
    !
    use io_hdf,  only : DFACC_CREATE, DFACC_WRITE, io_write

    implicit none

    ! in/out
    integer,intent(in)                     :: im, jm, lm, nt
    integer,intent(in)                     :: flag
    real,dimension(im,jm,lm,nt),intent(in) :: field     ! 4D field
    character(len=*),intent(in)            :: fname     ! file name
    character(len=*),intent(in)            :: namfield  ! field name

    ! local
    integer              :: sfstart, io, sfend

    if ( flag==0 ) then
       io = sfstart(fname,DFACC_CREATE)
    else
       io = sfstart(fname,DFACC_WRITE)
       if ( io == -1 ) io = sfstart(fname,DFACC_CREATE)
    end if
    call io_write(io,im,'X',jm,'Y',lm,'Z',nt,'N',field,namfield)
    io = sfend(io)

  end subroutine dumpfield



  subroutine dumpfieldi(flag,fname,im,jm,lm,nt,field,namfield)
    !
    ! Write 4D field (type integer) to HDF file
    !
    use io_hdf, only : DFACC_CREATE, DFACC_WRITE, io_write

    implicit none

    ! in/out
    integer,intent(in)                        :: im, jm, lm, nt
    integer,intent(in)                        :: flag
    integer,dimension(im,jm,lm,nt),intent(in) :: field     ! 4D field
    character(len=*),intent(in)               :: fname     ! file name
    character(len=*),intent(in)               :: namfield  ! field name

    ! local
    integer              :: sfstart,io,sfend

    if ( flag == 0 ) then
       io = sfstart(fname,DFACC_CREATE)
    else
       io = sfstart(fname,DFACC_WRITE)
       if (io==-1) io = sfstart(fname,DFACC_CREATE)
    end if
    call io_write(io,im,'X',jm,'Y',lm,'Z',nt,'N',field,namfield)
    io = sfend(io)

  end subroutine dumpfieldi



  subroutine coarsen_emission_1d(name_field,jm_emis,field_in,field,avg)
    !
    ! Purpose: Transform the 1D field available on e.g. 1 degree resolution 
    !          to the desired zoom geometry
    ! name_field : name, only for printing reasons
    ! jm_emis    : the dimension of the input field 
    ! field_in   : the input field
    ! field      : type d2_data  (defined in chem_param)
    ! avg        : flags the mode: avg = 1 means that a surface area 
    !              weighted area is calulated (e.g. land fraction)
    !              avg = 0 means that the sum of the underlying field 
    !              is taken (e.g. emissions in kg/month).
    !
    use dims,       only: nregions, dy, yref, ybeg, jm, dxy11, idate
    use global_types, only: d2_data
    implicit none

    ! in/out
    character(len=*),intent(in)               :: name_field
    integer,intent(in)                        :: avg  ! 0=no 1=yes
    integer,intent(in)                        :: jm_emis
    real,dimension(jm_emis),intent(in)        :: field_in
    type(d2_data),dimension(nregions),target  :: field
    ! contains per region field d2(jm)

    ! local
    real,dimension(:),pointer    :: coarse_field
    integer :: region
    real    :: scale,w,wtot
    integer :: ystart
    integer :: yres
    integer :: j,j_index,jj
    integer :: jm_start

    ! start

    do region=1,nregions      ! main loop over regions

       yres=dy/yref(region)

       if ( yres < 1 ) then
          call escape_tm('coarsen_emission_1d:'//&
               ' 1 degree minimum resolution in emissions')
       end if
       if ( jm_emis /= 180 ) then
          call escape_tm( &
               'coarsen_emission_1d: input resolution should be 1 degree')
       end if

       !WP! find index of southmost gridpoint based on 1x1 degree
       jm_start=ybeg(region)+91

       if ( jm_start <= -1 ) then
          call escape_tm('coarsen_emission_1d:'//&
               ' requested emission fields outside valid region')
       end if

       coarse_field => field(region)%d2
       do j=1,jm(region)
          !cycle through 1x1 grid with steps for this region
          j_index=jm_start+(j-1)*yres
          coarse_field(j) = 0.0
          wtot = 0.0
          do jj=0,yres-1
             if(avg==1) then
                w = dxy11(j_index+jj)
                wtot = wtot + w
                coarse_field(j)=coarse_field(j) + field_in(j_index+jj)*w
             else
                coarse_field(j)=coarse_field(j) + field_in(j_index+jj)
             end if
          end do
          if ( avg == 1 )  coarse_field(j) = coarse_field(j)/wtot
       end do

       !if ( avg == 0 ) then
       !   write(gol,'(a,i1,x,a12,a,i2,a,1pe14.3)') &
       !        ' coarsen_emission_1d: region ',region,name_field, &
       !        ' Field coarsened month :',idate(2),'total:',&
       !        sum(coarse_field) ; call goPr
       !else
       !   write(gol,'(a,i1,x,a12,a,i2,a,1pe14.3)') &
       !        ' coarsen_emission_1d: region:',region,name_field, &
       !        ' Field averaged month :',idate(2) ; call goPr
       !end if

       nullify(coarse_field)

    end do ! loop over regions....

  end subroutine coarsen_emission_1d


  !
  ! Purpose: Transform the 2D emissions available on 
  !          e.g. 1x1 degree resolution 
  !          to the desired zoom geometry
  ! name_field : name, only for printing reasons
  ! im_emis    : the x-dimension of the input field 
  ! jm_emis    : the y-dimension of the input field 
  ! field_in   : the input field
  ! field      : type emis_data  (defined in chem_param)
  ! avg        : flags the mode: avg = 1 means that a surface area 
  !              weighted area is calulated (e.g. land fraction)
  !              avg = 0 means that the sum of the underlying field 
  !              is taken (e.g. emissions in kg/month).
  !

  subroutine coarsen_emission_2d(name_field,im_emis,jm_emis,field_in,field,avg)

    use GO          , only : gol, goPr, goErr
    use Grid        , only : TllGridInfo, Init, Done, FillGrid
    use dims        , only : nregions
    use global_types, only : emis_data
    use meteodata   , only : lli    

    ! --- in/out ------------------------------------
    
    character(len=*),intent(in)   :: name_field
    integer,intent(in)            :: avg  ! 0=no 1=yes
    integer,intent(in)            :: im_emis,jm_emis
    real, intent(in)              :: field_in(im_emis,jm_emis)
    type(emis_data), target       :: field(nregions)
    ! contains per region field surf(im,jm)

    ! --- local -------------------------------------
    
    integer                     ::  region
    type(TllGridInfo)           ::  lli_in
    character(len=32)           ::  combkey
    integer                     ::  status

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

    ! define input grid
    call Init( lli_in, -180.0+0.5*360.0/im_emis, 360.0/im_emis, im_emis, &
                        -90.0+0.5*180.0/jm_emis, 180.0/jm_emis, jm_emis, status )
    if (status/=0) call escape_tm('toolbox/coarsen_emission_2d: error from init lli')
    
    ! sum or average ?
    select case ( avg )
      case ( 0 ) ; combkey = 'sum'
      case ( 1 ) ; combkey = 'area-aver'
      case default
        write (gol,'("ERROR - unsupported avg = ",i6)') avg; call goErr
        call escape_tm('ERROR - in toolbox/coarsen_emission_2d')
    end select

    ! main loop over regions
    do region = 1, nregions
    
      ! convert grid:
      call FillGrid( lli(region), 'n', field(region)%surf, &
                     lli_in     , 'n', field_in          , trim(combkey), status )
      if (status/=0) call escape_tm('toolbox/coarsen_emission_2d: error from fillgrid')
      
      !! info ...
      !select case ( combkey )
      !  case ( 'aver', 'area-aver' )
      !    write (gol,'("coarsen_emission : region ",i1," ",a," field averaged")') &
      !             region, trim(name_field); call goPr
      !  case ( 'sum' )
      !    write (gol,'("coarsen_emission : region ",i1," ",a," field coarsened/distr.; total ",es12.3)') &
      !             region, trim(name_field), sum(field(region)%surf); call goPr
      !  case default
      !    write (gol,'("coarsen_emission : region ",i1," ",a," field coarsened/distributed")') &
      !             region, trim(name_field); call goPr
      !end select

    end do ! regions
    
    ! done
    call Done( lli_in, status )
    if (status/=0) call escape_tm('toolbox/coarsen_emission_2d: error from done lli')

    ! ok
    status = 0
    
  end subroutine coarsen_emission_2d



   !
   ! Purpose: Transform the 3D emissions available on 
   !          e.g. 1x1 degree resolution 
   !          to the desired zoom geometry
   ! name_field : name, only for printing reasons
   ! im_emis    : the x-dimension of the input field 
   ! jm_emis    : the y-dimension of the input field 
   ! lm_emis    : the z-dimension of the input field 
   ! field_in   : the input field
   ! field      : type d3_data  (defined in chem_param)
   ! avg        : flags the mode: avg = 1 means that a surface area 
   !              weighted area is calulated (e.g. land fraction)
   !              avg = 0 means that the sum of the underlying field 
   !              is taken (e.g. emissions in kg/month).
   !

  subroutine coarsen_emission_3d( name_field, im_emis, jm_emis, lm_emis, &
                                  field_in, field, avg )

    use GO          , only : gol, goPr, goErr
    use Grid        , only : TllGridInfo, Init, Done, FillGrid
    use dims        , only : nregions
    use global_types, only : d3_data
    use meteodata   , only : lli    

    ! --- in/out ------------------------------------
    
    character(len=*), intent(in)    ::  name_field
    integer, intent(in)             ::  avg  !0=no 1=yes
    integer, intent(in)             ::  im_emis, jm_emis, lm_emis
    real, intent(in)                ::  field_in(im_emis,jm_emis,lm_emis)
    type(d3_data), target           ::  field(nregions)
    ! contains per region 3d field d3(im,jm,lm)

    ! --- local -----------------------------
    
    integer                     ::  region
    integer                     ::  l
    type(TllGridInfo)           ::  lli_in
    character(len=32)           ::  combkey
    integer                     ::  status

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

    ! define input grid
    call Init( lli_in, -180.0+0.5*360.0/im_emis, 360.0/im_emis, im_emis, &
                        -90.0+0.5*180.0/jm_emis, 180.0/jm_emis, jm_emis, status )
    if (status/=0) call escape_tm('toolbox/coarsen_emission_3d: error from init lli')
    
    ! sum or average ?
    select case ( avg )
      case ( 0 ) ; combkey = 'sum'
      case ( 1 ) ; combkey = 'area-aver'
      case default
        write (gol,'("ERROR - unsupported avg = ",i6)') avg; call goErr
        call escape_tm('ERROR - in toolbox/coarsen_emission_2d')
    end select

    ! main loop over regions
    do region = 1, nregions

      ! convert grid:
      do l = 1, lm_emis
        call FillGrid( lli(region), 'n', field(region)%d3(:,:,l), &
                       lli_in     , 'n', field_in(:,:,l)        , trim(combkey), status )
        if (status/=0) call escape_tm('toolbox/coarsen_emission_3d: error from fillgrid')
      end do
      
      !! info ...
      !select case ( combkey )
      !  case ( 'aver', 'area-aver' )
      !    write (gol,'("coarsen_emission : region ",i1," ",a," field averaged")') &
      !             region, trim(name_field); call goPr
      !  case ( 'sum' )
      !    write (gol,'("coarsen_emission : region ",i1," ",a," field coarsened/distr.; total ",es12.3)') &
      !             region, trim(name_field), sum(field(region)%d3); call goPr
      !  case default
      !    write (gol,'("coarsen_emission : region ",i1," ",a," field coarsened/distributed")') &
      !             region, trim(name_field); call goPr
      !end select

    end do ! regions
    
    ! done
    call Done( lli_in, status )
    if (status/=0) call escape_tm('toolbox/coarsen_emission_3d: error from done lli')

    ! ok
    status = 0
    
  end subroutine coarsen_emission_3d



  subroutine coarsen_emission_d23( name_field, jm_emis, lm_emis, &
                                   field_in, field, avg )
    !
    ! Purpose: Transform the 2D emissions available on lat-pressure resolution 
    !          to the desired zoom geometry
    ! name_field : name, only for printing reasons
    ! jm_emis    : the y-dimension of the input field 
    ! lm_emis    : the z-dimension of the input field 
    ! field_in   : the input field
    ! field      : type d23_data (defined in chem_param)
    ! avg        : flags the mode: avg = 1 means that a surface area 
    !              weighted area is calulated (e.g. land fraction)
    !              avg = 0 means that the sum of the underlying field 
    !              is taken (e.g. emissions in kg/month).
    !
    use dims,       only: nregions, dy, yref, ybeg, jm, lm, dxy11, idate
    use global_types, only: d23_data
    implicit none

    ! in/out
    character(len=*),intent(in)                        :: name_field
    integer,intent(in)                                 :: avg  !0=no 1=yes
    integer,intent(in)                                 :: jm_emis,lm_emis
    real,dimension(jm_emis,lm_emis),intent(in)         :: field_in
    type(d23_data),dimension(nregions),target          :: field
    !                contains per region field d23(jm,lm)

    ! local
    real,dimension(:,:),pointer    :: coarse_field
    integer :: region
    real    :: scale,w,wtot
    integer :: ystart
    integer :: yres
    integer :: j,j_index,jj,l
    integer :: jm_start

    ! start

    do region=1,nregions      ! main loop over regions

       yres=dy/yref(region)

       if ( yres < 1 ) then
          call escape_tm('coarsen_emission_d23:'//&
               ' 1 degree minimum resolution in emissions')
       end if
       if ( jm_emis /= 180 ) then
          call escape_tm( 'coarsen_emission_d23:'//&
               'input resolution should be 1 degree')
       end if

       !WP! find index of southmost gridpoint based on 1x1 degree
       jm_start=ybeg(region)+91

       if(jm_start<=(-1)) then
          call escape_tm( 'coarsen_emission_d23:'//&
               'requested emission fields outside valid region')
       end if

       coarse_field => field(region)%d23
       do l=1,lm(region)
          do j=1,jm(region)
             !cycle through 1x1 grid with steps for this region
             j_index=jm_start+(j-1)*yres
             coarse_field(j,l) = 0.0
             wtot = 0.0
             do jj=0,yres-1
                if(avg==1) then
                   w = dxy11(j_index+jj)
                   wtot = wtot + w
                   coarse_field(j,l)=coarse_field(j,l) + &
                        field_in(j_index+jj,l)*w
                else
                   coarse_field(j,l)=coarse_field(j,l) + &
                        field_in(j_index+jj,l)
                end if
             end do
             if ( avg == 1 )  coarse_field(j,l) = coarse_field(j,l)/wtot
          end do  !j
       end do  !l

       !if ( avg == 0 ) then
       !
       !   write(gol,'(a,i1,x,a12,a,i2,a,1pe14.3)') &
       !        ' coarsen_emission_d23: region ',region,name_field// &
       !        ' Field coarsened month ',idate(2),'  total ',&
       !        sum(coarse_field); call goPr
       !else
       !   write(gol,'(a,i1,x,a12,a,i2,a,1pe14.3)') &
       !        ' coarsen_emission_d23: region ',region,name_field// &
       !        ' Field averaged month ',idate(2); call goPr
       !end if

       nullify(coarse_field)

    end do ! loop over regions....

  end subroutine coarsen_emission_d23



  !--------------------------------------------------------------
  !
  ! abnormal termination of a model run
  ! 
  ! msg       character string to be printed on unit kmain
  !
  !-------------------------------------------------------------

  subroutine escape_tm( msg )

    use GO       , only : GO_Print_Done
    use dims     , only : okdebug, kmain, itau
    use datetime , only : wrtgol_tstamp
#ifdef MPI      
    use mpi_const, only : myid, root
    use mpi_comm , only : abortmpi
#endif
#ifdef with_prism
#ifdef oasis3
    use mod_prism_proto
#endif
#ifdef oasis4
    use PRISM    , only : prism_abort
#endif
    use TM5_Prism, only : comp_id
#endif

    ! --- in/out -----------------------------
    
    character(len=*), intent(in)  ::  msg
    
    ! --- const ------------------------------

    character(len=*), parameter   ::  rname = mname//'/escape_tm'
    
    ! --- local ------------------------------
    
    integer         ::  status

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

    write (gol,'(1x,39("--"))'); call goErr
    write (gol,'(1x,39("--"))'); call goErr
    call wrtgol_tstamp(itau,'ERROR - ESCAPE_TM'); call goErr
    write (gol,'(1x,a)') trim(msg); call goErr
    write (gol,'(1x,39("--"))'); call goErr
    write (gol,'(1x,39("--"))'); call goErr

    ! try to display at least the messages in a proper way ...
    call GO_Print_Done( status )
    if (status/=0) write (*,'("ERROR status from GO_Print_Done : ",i6)') status

#ifdef with_prism
#ifdef oasis3
    call PRISM_Abort_Proto( comp_id, rname, 'called prism_abort' )
#endif
#ifdef oasis4
    call prism_abort( comp_id, rname, 'called prism_abort' )
#endif
#endif

#ifdef MPI
    call abortmpi
#endif

    stop 'Fortran stop in escape_tm'

  end subroutine escape_tm



  ! 
  !-----------------------------------------------------------------------
  !                    TM5                                               !
  !-----------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  distribute_emis2D
  !
  ! !DESCRIPTION: subroutine to distribute the emissions given as a TM5 2D field
  !               into a TM5 3D emission field. Hlow and Hhigh are the bounds of
  !               the 2D emission fields. They give the height RELATIVE to oro
  !               From that the distribution is calculated
  !               employing the geopotential height (relative to oro)
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine distribute_emis2D( region, emis2D, emis3D, hlow, hhigh, status, xfrac,lat_start,lat_end)
    !
    ! !USES:
    !
    use dims,         only : lm, im, jm, dy, yref, ybeg
    use meteodata,    only : gph_dat
    use global_types, only : d3_data, emis_data
    use OMP_partools, only : my_omp_domain, my_omp_max
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)                   :: region
    type(emis_data), intent(in), target   :: emis2D    ! 2D emission field (kg/gridbox/month)
    type(d3_data), intent(inout), target  :: emis3D    ! 3D emission field (kg/gridbox/month)
    real,intent(in)                       :: hlow      ! lowest level (m)
    real,intent(in)                       :: hhigh     ! highest level (m)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)                  :: status
    real,intent(in), optional             :: xfrac     ! fraction of emissions to put
    real,intent(in), optional             :: lat_start
    real,intent(in), optional	          :: lat_end     ! latitude range for which the evaluation should be done.
    !
    ! !REVISION HISTORY: 
    !   16 Jul 2010 - Achim Strunk - Added Lat_start and lat_end optional inputs.
    !
    ! !REMARKS:
    !
    !EOP
    !----------------------------------------------------------------------
    !BOC


    ! --- const ------------------------------

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

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

    real, dimension(:,:,:),pointer        :: gph       ! geopotential height (m)
    real, dimension(:,:,:),pointer        :: e3d       ! 3D emissions
    real, dimension(:,:),pointer          :: e2d       ! 2D emissions
    integer                               :: i,j,l,j_st,j_end
    real, dimension(lm(1)+1)              :: height   
    real, dimension(lm(1))                :: dz
    real                                  :: dy1,dze
    real                                  :: totw, f, tot2d, tot3db, tot3da, fraction
    integer, dimension(3)                 :: ubound_e3d
    integer                               :: lmmax
    real                                  :: hhighb
    real                                  :: lat_low,lat_high
    ! Openmp paramters
    integer                               :: js, je
    real                                  :: l_status, g_status

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

    if (present(xfrac)) then
       fraction = xfrac
    else
       fraction = 1.0
    endif

    if (present(lat_start)) then
       lat_low = lat_start
    else
       lat_low = -90.
    endif

    if (present(lat_end)) then
       lat_high = lat_end
    else
       lat_high = 90.
    endif

    dy1=dy/yref(region)
    j_st  = floor ((lat_low  - ybeg(region))/dy1) + 1
    j_end = floor((lat_high - ybeg(region))/dy1) + 1
    if ((j_end .lt. 1 )         .or.  (j_st  .gt. jm(region)) )return ! outside requested range
    if ((j_end .gt. jm(region)) .and. (j_st  .lt. jm(region)))  j_end = jm(region) ! limit region
    if ((j_st  .lt. 1 )         .and. (j_end .gt. 1         ))  j_st = 1           ! limit region

    dze = hhigh - hlow  
    if ( (hhigh <= 0.0) .or. (hlow < 0.0) .or. (dze < 0.0) ) then
      write (gol,'("found strange emission heights:")'); call goErr
      write (gol,'("  hhigh      (>  0.0 ?) : ",es12.4)') hhigh; call goErr
      write (gol,'("  hlow       (>= 0.0 ?) : ",es12.4)') hlow; call goErr
      write (gol,'("  hhigh-hlow (>  0.0 ?) : ",es12.4)') dze; call goErr
      write (gol,'("in ",a)') rname; call goErr; status=1; return
    end if

    nullify(gph)
    nullify(e2d)
    nullify(e3d)
    gph => gph_dat(region)%data
    e2d => emis2d%surf
    e3d => emis3d%d3
    ubound_e3d = ubound(e3d)
    lmmax = ubound_e3d(3)
    tot2d = sum(e2d(:,j_st:j_end)*fraction)
    tot3db = sum(e3d(:,j_st:j_end,:))

    !$OMP PARALLEL &
    !$OMP  default (none) &
    !$OMP  shared  (region, g_status, gph, hlow, hhigh, lmmax, &
    !$OMP           dze, e2d, e3d, fraction, j_st, j_end) &
    !$OMP  private (i, j, js, je, l_status, height, dz, totw, gol, f)
    call my_omp_domain (j_st,j_end,js,je)

    l_status = 0.0
    do j=js,je
       do i=1,im(region)

          ! zero based height:
          height(1) = 0.0 
          do l=1, lm(region)
             dz(l) = gph(i,j,l+1) - gph(i,j,l)
             height(l+1) = height(l) + dz(l)
          enddo
          totw = 0.0
          if( hhigh > height(lmmax+1) ) then
            write (gol,'("try to put emissions higher than array allows:")'); call goErr
            write (gol,'("  hhigh            : ",es12.4)') hhigh; call goErr
            write (gol,'("  height(lmmax+1)  : ",es12.4)') height(lmmax+1); call goErr
            do l = 1, size(height)
              write (gol,*) 'height : ', l, height(l); call goErr
            end do
            write (gol,'("in ",a)') rname; call goErr; l_status=1.0
          endif   
          zz: do l=1, lmmax
              if (hhigh > height(l+1)) then 
                 if ( hlow < height(l) ) then
                    f =  dz(l)/dze
                    totw = totw + f
                    e3d(i,j,l) = e3d(i,j,l) + f*fraction*e2d(i,j)
                 else if( hlow < height(l+1)) then
                    f = (height(l+1)-hlow)/dze
                    totw = totw + f
                    e3d(i,j,l) = e3d(i,j,l) + f*fraction*e2d(i,j)
                 endif
              else
                 if ( hlow < height(l)) then
                    f = (hhigh - height(l))/dze
                    totw = totw + f
                    e3d(i,j,l) = e3d(i,j,l) + f*fraction*e2d(i,j)
                 else
                    totw = totw + 1.0
                    e3d(i,j,l) = e3d(i,j,l) + fraction*e2d(i,j)
                 endif
                 exit zz
              endif
          enddo zz
          if ( abs(totw-1.0) > 1e-14 ) then
            write (gol,'("sum weight /= 1 : ",es12.4)') totw-1.0; call goErr
            write (gol,'("in ",a)') rname; call goErr; l_status=1.0
          end if
       enddo   !i
    enddo   !j
    call my_omp_max(l_status,g_status)
    !$OMP END PARALLEL

    if (g_status > 0.5) then
       status = 1
       return
    end if

    tot3da = sum(e3d(:,j_st:j_end,:))
    if (abs((tot3da-tot3db)-tot2d) > 1e-8*tot2d ) then
      write (gol,'("emissions have not been distributed mass-conserving")'); call goErr
      write (gol,'("in ",a)') rname; call goErr; status=1; return
    end if

    nullify(gph)
    nullify(e2d)
    nullify(e3d)

    ! ok
    status = 0

  end subroutine distribute_emis2D
  !EOC



  !
  !

  !---------------------------------------------------------------------------
  !                    TM5                                                    !
  !---------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  distribute1x1
  !
  ! !DESCRIPTION: subroutine to distribute 1x1 emissions liniairly between
  !               hlow and hhigh. The vertical level is determined by
  !               the orography which is read from the surface file...
  !               A simple scale height vertical structure is assumed.
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine distribute1x1(emi1x1, hlow, hhigh, emis3d, xfrac)

  use Binas    , only : grav
  use dims     , only : nlon360, nlat180, lm, itau, staggered, at, bt
  use dims     , only : iglbsfc
  use MeteoData, only : oro_dat

  ! --- in/out ----------------------------------
  
  real, dimension(nlon360,nlat180), intent(in)            :: emi1x1  ! (kg/1x1 gridbox) 2D field of emissions
  real, dimension(nlon360,nlat180), intent(in)            :: hlow    ! (m) lower bound of emission
  real, dimension(nlon360,nlat180), intent(in)            :: hhigh   ! (m) higher bound of emission
  real, dimension(nlon360,nlat180,lm(1)), intent(inout)   :: emis3d  ! (kg/box) distributed in height
  real,intent(in), optional                               :: xfrac   ! fraction of emissions to put

  ! --- const ----------------------------------
  
  character(len=*), parameter  ::  rname = mname//'/distribute1x1'
  
  !--- local ------------------------------------
  
  real, parameter                                :: scalh = 8000.0
  real                                           :: p0, pt
  real,dimension(lm(1))                          :: height, dz
  integer                                        :: i,j,l
  real                                           :: hh,hl,dze,totw,f,e3da,e3db, fract
  
  ! --- begin -----------------------------------

  ! check ...
  if ( .not. oro_dat(iglbsfc)%used ) then
    write (gol,'("oro_dat(iglbsfc) not in use ...")'); call goErr
    TRACEBACK; call escape_tm('escape')
  end if
  
  if (present(xfrac)) then
     fract = xfrac
  else
     fract = 1.0
  endif
  e3db = sum(emis3d)
  do j=1,nlat180
    do i=1,nlon360
       if (fract*emi1x1(i,j) > 1e-14) then 
          height(1) = oro_dat(iglbsfc)%data(i,j,1)/grav
          p0 = 1e5*exp(-height(1)/scalh)
          do l=1,lm(1)-1   ! bug reported by FD: alog(0) crashes!
             pt = at(l+1) + bt(l+1)*p0
             height(l+1) = height(1)-scalh*alog(pt/p0)
             dz(l) = height(l+1)-height(l)
          enddo
          hl = max(height(1),hlow(i,j))
          hh = max(hhigh(i,j),height(1))
          dze = hh-hl
          if(dze < 0.0) then 
             call escape_tm('In distribute1x1: dze <= 0')
          else if ( dze == 0.0) then   ! this somehow happens!
             hh = height(1)+1.0
             hl = height(1)
          endif   
          totw = 0.0
          zz: do l=1, lm(1)
              if (hh > height(l+1)) then 
                 if ( hl < height(l) ) then
                    f =  dz(l)/dze
                    totw = totw + f
                    emis3d(i,j,l) = emis3d(i,j,l) + f*fract*emi1x1(i,j)
                 else if( hl < height(l+1)) then
                    f = (height(l+1)-hl)/dze
                    totw = totw + f
                    emis3d(i,j,l) = emis3d(i,j,l) + f*fract*emi1x1(i,j)
                 endif
              else
                 if ( hl < height(l)) then
                    f = (hh - height(l))/dze
                    totw = totw + f
                    emis3d(i,j,l) = emis3d(i,j,l) + f*fract*emi1x1(i,j)
                 else
                    totw = totw + 1.0
                    emis3d(i,j,l) = emis3d(i,j,l) + fract*emi1x1(i,j)
                 endif
                 exit zz
              endif
          enddo zz
          if ( abs(totw-1.0) > 1e-14 ) & 
               call escape_tm(' Routine distribute1x1 : sum weight /= 1')
       endif 
    enddo
  enddo
  e3da = sum(emis3d)
  if (abs(e3da-e3db-sum(fract*emi1x1)) > e3da*1e-8 ) &
                call escape_tm(' Routine distribute1x1 : distributed amount differs!')
  end subroutine distribute1x1


  !
  ! subroutine to distribute 1x1 emissions liniairly between
  ! hlow and hhigh. The vertical level is determined by
  ! the orography which is read from the surface file...
  ! A simple scale height vertical structure is assumed.
  ! same as distribute1x1 but hlow, hhigh now scalar
  ! ALSO: the height is now defined relative to the orography!!!
  !

  subroutine distribute1x1b(emi1x1, hlow, hhigh, emis3d, xfrac)

  use Binas    , only : grav
  use dims     , only : nlon360, nlat180, lm, itau, staggered, at, bt
  use dims     , only : iglbsfc
  use MeteoData, only : oro_dat

  ! --- in/out ----------------------------------

  real, dimension(nlon360,nlat180), intent(in)            :: emi1x1  ! (kg/1x1 gridbox) 2D field of emissions
  real,  intent(in)                                       :: hlow    ! (m) lower bound of emission
  real,  intent(in)                                       :: hhigh   ! (m) higher bound of emission
  real, dimension(nlon360,nlat180,lm(1)), intent(inout)   :: emis3d  ! (kg/box) distributed in height
  real,intent(in), optional                               :: xfrac   ! fraction of emissions to put

  ! --- const ----------------------------------
  
  character(len=*), parameter  ::  rname = mname//'/distribute1x1b'
  
  !--- local ------------------------------------
  
  real, parameter                                :: scalh = 8000.0
  real                                           :: p0, pt
  real,dimension(lm(1))                          :: height, dz
  integer                                        :: i,j,l
  real                                           :: hh,hl,dze,totw,f,e3da,e3db, fract, hlow_oro, hhigh_oro

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

  ! check ...
  if ( .not. oro_dat(iglbsfc)%used ) then
    write (gol,'("oro_dat(iglbsfc) not in use ...")'); call goErr
    TRACEBACK; call escape_tm('escape')
  end if
  
  if (present(xfrac)) then
     fract = xfrac
  else
     fract = 1.0
  endif
  e3db = sum(emis3d)
  do j=1,nlat180
    do i=1,nlon360
       if (fract*emi1x1(i,j) > 1e-14) then 
          height(1) = oro_dat(iglbsfc)%data(i,j,1)/grav
          hlow_oro = hlow + height(1)
          hhigh_oro = hhigh + height(1)
          p0 = 1e5*exp(-height(1)/scalh)
          do l=1,lm(1)-1   ! bug reported by FD: alog(0) crashes!
             pt = at(l+1) + bt(l+1)*p0
             height(l+1) = height(1)-scalh*alog(pt/p0)
             dz(l) = height(l+1)-height(l)
          enddo
          hl = max(height(1),hlow_oro)
          hh = max(hhigh_oro,height(1))
          dze = hh-hl
          if(dze < 0.0) then 
             call escape_tm('In distribute1x1b: dze <= 0')
          else if ( dze == 0.0) then   ! this somehow happens!
             hh = height(1)+1.0
             hl = height(1)
          endif   
          totw = 0.0
          zz: do l=1, lm(1)
              if (hh > height(l+1)) then 
                 if ( hl < height(l) ) then
                    f =  dz(l)/dze
                    totw = totw + f
                    emis3d(i,j,l) = emis3d(i,j,l) + f*fract*emi1x1(i,j)
                 else if( hl < height(l+1)) then
                    f = (height(l+1)-hl)/dze
                    totw = totw + f
                    emis3d(i,j,l) = emis3d(i,j,l) + f*fract*emi1x1(i,j)
                 endif
              else
                 if ( hl < height(l)) then
                    f = (hh - height(l))/dze
                    totw = totw + f
                    emis3d(i,j,l) = emis3d(i,j,l) + f*fract*emi1x1(i,j)
                 else
                    totw = totw + 1.0
                    emis3d(i,j,l) = emis3d(i,j,l) + fract*emi1x1(i,j)
                 endif
                 exit zz
              endif
          enddo zz
          if ( abs(totw-1.0) > 1e-14 ) & 
               call escape_tm(' Routine distribute1x1 : sum weight /= 1')
       endif 
    enddo
  enddo
  e3da = sum(emis3d)
  if (abs(e3da-e3db-sum(fract*emi1x1)) > e3da*1e-8 ) &
                call escape_tm(' Routine distribute1x1 : distributed amount differs!')
  end subroutine distribute1x1b


  ! routine to integrate tropospheric (ozone) 
  ! note: routine is now written for Ozone, but may be changed
  ! to be more general. The definition of tropopause is critical for ozone.
  ! here, the slope is used in the interpolation.
  ! multiple tropopause values are ignored, but are known to 
  ! occur. The lowest crossing of thres is used.
  ! In the lower atmosphere > 600 hPa, values > thres are ignored.

  subroutine tropospheric_columns(region,field,slope, column,thres, xmo3)

    use binas,       only : Dobs, xmair, Avog 
    use global_data, only : region_dat
    use meteodata  , only : phlb_dat, m_dat
    use Dims,        only : lm, isr, ier, jsr, jer, CheckShape, im, jm

    !input 
    integer, intent(in)               :: region  ! region in the TM5 zoom definition
    real, dimension(:,:,:),intent(in) :: field   ! 3D field of tracer (O3)
    real, dimension(:,:,:),intent(in) :: slope   ! 3D field of tracer z -slope (O3)
    real, dimension(:,:),intent(out)  :: column  ! output: tropospheric column in DU
    real, intent(in)                  :: thres   ! ppb threshold for tropospheric column
    real, intent(in)                  :: xmo3    ! mol mass tracer
    ! locals and pointers
#ifdef with_zoom
    integer,dimension(:,:),pointer               :: zoomed
#endif
    real,dimension(:,:,:),pointer                :: m
    real,dimension(:,:,:),pointer                :: phlb
    real,dimension(:),pointer                    :: dxyp

    real, dimension(2*lm(1)+1)                   :: o3a
    real     :: o3mm, o3mmu, o3mml, o3trop, o3ml, o3mu, frac
    integer  :: i,j,l,ip
    ! start
      call CheckShape( (/im(region),jm(region)/), shape(column) )
      call CheckShape( (/im(region),jm(region), lm(region)/), shape(field) )
      call CheckShape( (/im(region),jm(region), lm(region)/), shape(slope) )

      dxyp=> region_dat(region)%dxyp
#ifdef with_zoom
      zoomed => region_dat(region)%zoomed
#endif
      phlb => phlb_dat(region)%data
      m => m_dat(region)%data

      ! collect ozone on mid levels and at boundaries
      ! average the estimates from upper/lower gridboxes
      ! except for bottom and top.
      ! 
      column(:,:) = 0.0
      do j=jsr(region), jer(region)
         do i=isr(region), ier(region)
#ifdef with_zoom
           if(zoomed(i,j)/=region) cycle
#endif
           do l=1,lm (region)
              ip = 2*l-1
              o3mm = xmair/xmo3*1e9*field(i,j,l)/m(i,j,l)    ! level
              o3mmu = xmair/xmo3*1e9*(field(i,j,l)+slope(i,j,l))/m(i,j,l) !  top
              o3mml = xmair/xmo3*1e9*(field(i,j,l)-slope(i,j,l))/m(i,j,l) !  bottom
              if(l == 1) then 
                 o3a(ip) = o3mml
              else
                 o3a(ip) = o3a(ip) + 0.5*o3mml
              endif
              o3a(ip+1) = o3mm
              if(l == lm(region)) then 
                 o3a(ip+2) = o3mmu
              else
                 o3a(ip+2) = 0.5*o3mmu
              endif
           enddo

           o3trop = 0.0
           height: do l=1,lm (region)
              ip = 2*l-1
              ! split gridbox in upper and lower part
              ! for more accurate interpolation
              o3ml = 0.5*field(i,j,l) - 0.25*slope(i,j,l)
              o3mu = 0.5*field(i,j,l) + 0.25*slope(i,j,l)
              if (phlb(i,j,l) > 60000.0) then  ! avoid surface o3>150ppb
                 o3trop = o3trop + field(i,j,l)
                 cycle height
              endif
              if (phlb(i,j,l+1) < 7000.0) exit height  ! now about time 70 hPa at top
              if(o3a(ip)  < thres) then   ! bottom value less than thres
                 if(o3a(ip+1) >= thres) then  ! but central value is not
                    frac = (thres - o3a(ip))/(o3a(ip+1)-o3a(ip))
                    o3trop = o3trop + frac*o3ml
                    exit height
                 else if (o3a(ip+2) >= thres) then  ! but upper is not
                    frac = (thres - o3a(ip+1))/(o3a(ip+2)-o3a(ip+1))
                    o3trop = o3trop + o3ml + frac*o3mu
                    exit height
                 else  ! entire layer is not
                    o3trop = o3trop + field(i,j,l)
                 endif  
              else
                 exit height
              endif
           enddo height
           column(i,j) = o3trop*1e3/xmo3*Avog*1e-4/dxyp(j)/Dobs  ! kg ->dobs
        enddo ! i
      enddo !j
      call dumpfield(0,'dump.hdf', im(region), jm(region), lm(region), 1, field, 'o3')
      call dumpfield(1,'dump.hdf', im(region), jm(region), lm(region), 1, slope, 'o3slope')
      call dumpfield(1,'dump.hdf', im(region), jm(region), 1, 1, column, 'o3column')
      call dumpfield(1,'dump.hdf', im(region), jm(region), lm(region)+1, 1, phlb, 'phlb')
      call dumpfield(1,'dump.hdf', im(region), jm(region), lm(region), 1, m, 'm')
      call dumpfield(1,'dump.hdf', jm(region),1, 1,  1, dxyp, 'dxyp')
      nullify(dxyp)
#ifdef with_zoom
      nullify(zoomed)
#endif
      nullify(phlb)
      nullify(m)
  end subroutine tropospheric_columns

end module toolbox
