ecearth3/sources/tm5mp/base/toolbox.F90

!
#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, goBug
  use dims, only : okdebug
  
  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: zfarr                ! calculation of Arrhenius expression for rate constants
  public :: zf3bod               ! calculation of the reaction rate of 3 body reaction 
  public :: ltropo               ! find tropopause level from T gradient
  public :: ltropo_ifs           ! find tropopause level from T gradient, IFS way
  public :: lvlpress             ! return index of highest level below a pressure value
  public :: dumpfield            ! write 4D field (real) to file (debug tool)
  public :: dumpfieldi           ! write 4D field (real) to file (debug tool)
  public :: coarsen_emission     ! convert GLOBAL field to each region
  public :: escape_tm            ! error handler
  public :: distribute_emis2D    ! distribute vertically 2D field [work on MPI tiles]
  public :: distribute1x1        ! distribute vertically GLOBAL 1x1 2D field
  public :: distribute1x1b       ! alternate distribute1x1
  public :: tropospheric_columns ! integrate tropospheric ozone
  !
  ! !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; protex
  !   24 Jan 2012 - P. Le Sager - fixed coarsen_emission_2d and 3d for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !-----------------------------------------------------------------------

CONTAINS
  
  !-----------------------------------------------------------------------
  !                    TM5                                                !
  !-----------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:  ltropo
  !
  ! !DESCRIPTION: determine tropopause level from temperature gradient
  !	          reference: WMO /Bram Bregman
  !\\
  !\\
  ! !INTERFACE:
  !
  integer function ltropo(region,tx,gphx,lmx)
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in)               :: region
    integer,intent(in)               :: lmx
    real,dimension(lmx),intent(in)   :: tx
    real,dimension(lmx+1),intent(in) :: gphx 
    !
    ! !REVISION HISTORY: 
    !	programmed by		 fd 01-11-1996 
    !   changed to function      fd 12-07-2002
    !   16 Jul 2010 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
    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
  !EOC
    
  !-----------------------------------------------------------------------
  !                    TM5                                                !
  !-----------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:  ltropo_ifs
  !
  ! !DESCRIPTION: determine tropopause level from temperature gradient
  !               as it is done in IFS cmip6_strato_aero_process.F90 (EC-Earth 3.2.3)
  !\\
  !\\
  ! !INTERFACE:
  !
  integer function ltropo_ifs(region,i,j,tx,lmx)

    use binas,     only : rgas_air, grav
    use dims,      only : at, bt, lm
    use meteodata, only : phlb_dat, m_dat
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in)               :: region, i, j
    integer,intent(in)               :: lmx
    real,dimension(lmx),intent(in)   :: tx
    !
    ! !REVISION HISTORY: 
    !
    ! !REMARKS:
    !   returns the highest level *in* the troposhere, like 
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer                :: ltmin, ltmax, l, ltropo_ifs_default, jlev
    real                   :: dz, dt, lapse_rate, lapse_rate_mean
    real, dimension(lmx)   :: papf
    real, dimension(lmx+1) :: zpth
    logical :: ltest
    real, dimension(:,:,:), pointer :: phlb

    phlb => phlb_dat(region)%data
    ltest=.false.

    ! temperature at half levels
    do l=2,lmx
       zpth(l) = (tx(l) + tx(l-1))/2.
    enddo
    zpth(1)=tx(1)
    zpth(lmx+1)=tx(lmx)

    ! pressure at full-levels
    do l=1,lmx
       papf(l)=( phlb(i,j,l) + phlb(i,j,l+1) ) / 2.
    enddo

    ltropo_ifs = lmx
    ltropo_ifs_default = 1
    
    DO l=1,lmx
       
       ! register highest troposheric level
       if (papf(l) .gt. 10000.) ltropo_ifs_default = max(ltropo_ifs_default,l)

       if (papf(l)>7500. .and. papf(l)<55000.) then
          ! the lapse rate is computed for each layer (here we use the gas
          ! equation and the hydrostatic approximation):
          lapse_rate = (zpth(l)-zpth(l+1)) / (phlb(i,j,l)-phlb(i,j,l+1)) * papf(l)/tx(l)*grav/rgas_air  
          if (lapse_rate <= 0.002) then
             dz=0
             jlev=l
             do while (dz <= 2000) ! compute the delta(p) corresponding to 2km
                dz=dz+(papf(jlev)-papf(jlev+1))/phlb(i,j,jlev+1)*zpth(jlev+1)*rgas_air/grav
                jlev = jlev+1
             enddo

             lapse_rate_mean = (tx(l)-tx(jlev))/dz ! lapse rate over 2 km
             if (lapse_rate_mean <= 0.002) then
                ltropo_ifs =min(l,ltropo_ifs) 
                ltest=.true.
             endif
          endif
       endif

    enddo

    if (.not. ltest) then
       ltropo_ifs = ltropo_ifs_default
    end if

  end function ltropo_ifs
  
  !---------------------------------------------------------------------------
  !                    TM5                                                    !
  !---------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:  lvlpress
  !
  ! !DESCRIPTION: lvlpress determines the index of the level with a pressure
  !               greater than press i.e. in height below press
  !               determines level independent of vertical resolution lm 
  !               uses hybrid level coefficients to determine lvlpress
  !\\
  !\\
  ! !INTERFACE:
  !
  integer function lvlpress(region,press,press0)
    !
    ! !USES:
    !
    use dims, only: at, bt, lm
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: region
    real,    intent(in) :: press
    real,    intent(in) :: press0
    !
    ! !REVISION HISTORY: 
    !   programmed by           Peter van Velthoven, 6-november-2000
    !   16 Jul 2010 - Achim Strunk -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    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
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:  zfarr
  !
  ! !DESCRIPTION: ZFARR calculation of Arrhenius expression for rate constants
  !\\
  !\\
  ! !INTERFACE:
  !
  real function zfarr(rx1,er,ztrec) 
    !
    ! !INPUT PARAMETERS:
    !
    real,intent(in) :: rx1,er,ztrec            
    !
    ! !REVISION HISTORY: 
    !   16 Jul 2010 - Achim Strunk -
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
    zfarr=rx1*exp(er*ztrec)

  end function zfarr
  !EOC

  real function zf3bod(rx1,rx2,rx3)
  !----------------------------------------------------------------------
  !
  ! function to calculate the reaction rate of 3 body reaction
  !
  !------------------------------------------------------------------
  !
    real,intent(in) :: rx1,rx2,rx3
  !
    zf3bod=rx1/(1.+rx1/rx2)*rx3**(1./(1.+log10(rx1/rx2)**2))
  !
  end function zf3bod
    
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DumpField
  !
  ! !DESCRIPTION: Write 4D field (type real) to HDF file
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine dumpfield(flag,fname,im,jm,lm,nt,field,namfield)
    !
    ! !USES:
    !
#ifdef with_hdf4
    use io_hdf,  only : DFACC_CREATE, DFACC_WRITE, io_write
#endif
    !
    ! !INPUT PARAMETERS:
    !
    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
    !
    ! !REVISION HISTORY: 
    !   16 Jul 2010 - Achim Strunk - protex
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
#ifdef with_hdf4
    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)
#endif

  end subroutine dumpfield
  !EOC
    
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  dumpfieldi
  !
  ! !DESCRIPTION:  Write 4D field (type integer) to HDF file
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine dumpfieldi(flag,fname,im,jm,lm,nt,field,namfield)
    !
    ! !USES:
    !
#ifdef with_hdf4
    use io_hdf, only : DFACC_CREATE, DFACC_WRITE, io_write
#endif
    !
    ! !INPUT PARAMETERS:
    !
    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
    !
    ! !REVISION HISTORY: 
    !   16 Jul 2010 - Achim Strunk - protex
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
#ifdef with_hdf4
    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)
#endif

  end subroutine dumpfieldi
  !EOC
    
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    coarsen_emission_1d
  !
  ! !DESCRIPTION: Transform the 1D field available on e.g. 1 degree resolution 
  !               to the desired zoom geometry
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine coarsen_emission_1d( name_field, jm_emis, field_in, field, avg, status)
    !
    ! 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).
    !
    ! !USES:
    !
    use dims,       only: nregions, dy, yref, ybeg, jm, dxy11, idate
    use global_types, only: d2_data
    implicit none
    !
    ! !INPUT PARAMETERS:
    !
    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
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !    
    type(d2_data), dimension(nregions), target  :: field
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)          :: status
    !
    ! !REVISION HISTORY: 
    !   19 Jun 2012 - P. Le Sager - got rid of escape_tm
    !
    ! !REMARKS:
    !  (1) FIXME : should work as is for lon-lat MPI domain decomposition, but
    !       not tested.
    !  (2) This is limited to a 1 degree zonal global input [-90,+90] 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
    character(len=*), parameter ::  rname = mname//'/coarsen_emission_1d'

    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
          write (gol,'("1 degree minimum resolution in emissions")'); call goErr
          status=1
          IF_NOTOK_RETURN(status=1)          
       end if
       
       if ( jm_emis /= 180 ) then
          write (gol,'("input resolution should be 1 degree")'); call goErr
          status=1
          IF_NOTOK_RETURN(status=1)          
       end if

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

       if ( jm_start <= -1 ) then
          write (gol,'("requested emission fields outside valid region")'); call goErr
          status=1
          IF_NOTOK_RETURN(status=1)          
       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
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    COARSEN_EMISSION_2D
  !
  ! !DESCRIPTION: Transform a 2D *global* emissions, at a given resolution,  
  !                to the geometry of all model regions.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE COARSEN_EMISSION_2D(name_field, im_emis, jm_emis, field_in, field, avg, status)
    !
    ! !USES:
    !
    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 : global_lli    
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*), intent(in)  :: name_field        ! name, only for printing reasons
    integer,          intent(in)  :: im_emis, jm_emis  ! grid size of global input field
    real,             intent(in)  :: field_in(im_emis,jm_emis) ! input field
    integer,          intent(in)  :: avg  ! 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).
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !    
    type(emis_data), target       :: field(nregions)    ! per region, %surf(im,jm)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)          :: status
    !
    ! !REVISION HISTORY: 
    !   28 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer                     ::  region
    type(TllGridInfo)           ::  lli_in
    character(len=32)           ::  combkey

    ! --- 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_NOTOK_RETURN(status=1)

    ! sum or average
    status=0
    select case ( avg )
    case ( 0 ) ; combkey = 'sum'
    case ( 1 ) ; combkey = 'area-aver'
    case default
       write (gol,'("ERROR - unsupported avg = ",i6)') avg; call goErr
       status=1
    end select
    IF_NOTOK_RETURN(status=1)

    ! main loop over regions
    do region = 1, nregions

       ! convert grid:       
       call FillGrid( global_lli(region), 'n', field(region)%surf, &
                      lli_in,             'n', field_in,           trim(combkey), status)
       IF_NOTOK_RETURN(status=1)

    end do

    ! done
    call Done( lli_in, status )
    IF_NOTOK_RETURN(status=1)

  END SUBROUTINE COARSEN_EMISSION_2D
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    COARSEN_EMISSION_3D
  !
  ! !DESCRIPTION: Transform a 3D *global* emissions, at a given resolution,  
  !                to the geometry of all model regions.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE COARSEN_EMISSION_3D( name_field, im_emis, jm_emis, lm_emis, &
                                  field_in, field, avg, status )
    !
    ! !USES:
    !
    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 : global_lli    
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*), intent(in) ::  name_field
    integer,          intent(in) ::  avg
    integer,          intent(in) ::  im_emis, jm_emis, lm_emis
    real,             intent(in) ::  field_in(im_emis,jm_emis,lm_emis)
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !    
    type(d3_data), target        ::  field(nregions)  ! contains, per region, 3d field %d3(im,jm,lm)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)         :: status
    !
    ! !REVISION HISTORY: 
    !   28 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/coarsen_emission_3d'
    
    integer                     ::  region
    integer                     ::  l
    type(TllGridInfo)           ::  lli_in
    character(len=32)           ::  combkey

    ! --- 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_NOTOK_RETURN(status=1)
    
    ! sum or average
    status=0
    select case ( avg )
      case ( 0 ) ; combkey = 'sum'
      case ( 1 ) ; combkey = 'area-aver'
      case default
         write (gol,'("ERROR - unsupported avg = ",i6)') avg; call goErr
         status=1
    end select
    IF_NOTOK_RETURN(status=1)
     
    ! main loop over regions
    do region = 1, nregions

      ! convert grid:
      do l = 1, lm_emis
        call FillGrid( global_lli(region), 'n', field(region)%d3(:,:,l), &
                       lli_in,             'n', field_in(:,:,l),         trim(combkey), status)
        IF_NOTOK_RETURN(status=1)
      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_NOTOK_RETURN(status=1)

    ! ok
    status = 0

  END SUBROUTINE COARSEN_EMISSION_3D
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    COARSEN_EMISSION_D23
  !
  ! !DESCRIPTION: Transform the 2D emissions available on lat-pressure resolution 
  !               to the desired zoom geometry
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine coarsen_emission_d23( name_field, jm_emis, lm_emis,        &
                                   field_in,   field,   avg,     status )
    !
    ! 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).
    !
    ! !USES:
    !
    use dims,         only : nregions, dy, yref, ybeg, jm, lm, dxy11, idate
    use global_types, only : d23_data
    !
    ! !INPUT PARAMETERS:
    !
    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
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !    
    type(d23_data),dimension(nregions),target          :: field ! contains, per region, field %d23(jm,lm)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)         :: status
    !
    ! !REVISION HISTORY: 
    !   19 Jun 2012 - P. Le Sager - got rid of escape_tm
    !
    ! !REMARKS:
    !  (1) FIXME : should work as is for lon-lat MPI domain decomposition, but
    !       not tested.
    !  (2) First dimension of input is limited to a 1 degree global [-90,+90]
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter   :: rname = mname//'/coarsen_emission_d23'
    
    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
          write (gol,'("1 degree minimum resolution in emissions")'); call goErr
          status=1
          IF_NOTOK_RETURN(status=1)          
       end if

       if ( jm_emis /= 180 ) then
          write (gol,'("input resolution should be 1 degree")'); call goErr
          status=1
          IF_NOTOK_RETURN(status=1)          
       end if

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

       if ( jm_start <= -1 ) then
          write (gol,'("requested emission fields outside valid region")'); call goErr
          status=1
          IF_NOTOK_RETURN(status=1)          
       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
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:     ESCAPE_TM
  !
  ! !DESCRIPTION:  abnormal termination of a model run
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine escape_tm( msg )
    !
    ! !USES:
    !
    use GO,        only : GO_Print_Done
    use dims,      only : okdebug, kmain, itau
    use datetime,  only : wrtgol_tstamp
    use partools,  only : Par_StopMPI
#ifdef with_prism
#ifdef oasis3
    use mod_oasis
#endif
    use TM5_Prism, only : comp_id
#endif
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*), intent(in)  ::  msg ! character string to be printed on unit gol
    !
    ! !REVISION HISTORY: 
    !   19 Jun 2012 - P. Le Sager - use par_stopmpi instead of stop; protex
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter   ::  rname = mname//'/escape_tm'
    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 oasis_abort( comp_id, rname, 'called prism_abort' )
#endif
#endif

    call Par_StopMPI

  end subroutine escape_tm
  !EOC

  !-----------------------------------------------------------------------
  !                    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 tm5_distgrid, only : Get_DistGrid, dgrid
    use dims,         only : lm, im, jm, dy, yref, ybeg
    use meteodata,    only : gph_dat
    use global_types, only : d3_data, emis_data
    !
    ! !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 b/w HLOW and HHIGH (default=1)
    real,            intent(in), optional  :: lat_start ! lower limit (default=-90) of application domain
    real,            intent(in), optional  :: lat_end   ! upper limit (default=+90) of application domain
    !
    ! !REVISION HISTORY: 
    !   16 Jul 2010 - Achim Strunk - Added Lat_start and lat_end optional inputs.
    !   27 Mar 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !----------------------------------------------------------------------
    !BOC

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

    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
    integer                               :: i1, i2, j1, j2
    real                                  :: l_status, g_status

    ! --- begin ---------------------------------------------
    
    status = 0
    
    if (present(xfrac)) then
       fraction = xfrac
    else
       fraction = 1.0
    endif

    if (fraction < spacing(fraction)) return ! degenerated cases

    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

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

    ! get indices of application-region 
    dy1=dy/yref(region)
    j_st  = floor((lat_low  - ybeg(region))/dy1) + 1
    j_end = floor((lat_high - ybeg(region))/dy1) + 1

    ! is tile in appl. region?
    if ((j_end .lt. j1 ) .or.  (j_st  .gt. j2) ) then
       return
    end if

    ! limit range
    if (j_end .gt. j2 ) j_end = j2
    if (j_st  .lt. j1 ) j_st  = j1

    ! check inputs
    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
      TRACEBACK; status=1; return
    end if

    ! init
    nullify(gph)
    nullify(e2d)
    nullify(e3d)
    gph => gph_dat(region)%data
    e2d => emis2d%surf
    e3d => emis3d%d3

    ! get highest possible level
    ubound_e3d = ubound(e3d)
    lmmax = ubound_e3d(3)

    ! total of to-be-redistributed before redistribution
    tot2d = sum(e2d(:,j_st:j_end)*fraction)

    ! total of target array before adding redistributed data
    tot3db = sum(e3d(:,j_st:j_end,:))
    
    
    ! do work
    jlo: do j=j_st,j_end
       do i=i1, i2

          ! 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
             write (gol,'(" at i/j=",i3,1x,i3)') i,j                          ; call goErr
             do l = 1, size(height)
                write (gol,*) 'height : ', l, height(l); call goErr
             end do
             do l = 1, size(gph,3)
                write (gol,*) 'gph : ', l, gph(i,j,l); call goErr
             end do
             TRACEBACK; status=1
             exit jlo
          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
             TRACEBACK; status=1
             exit jlo
          end if
         
       enddo   !i
    enddo jlo  !j
    IF_NOTOK_RETURN(status=1)

    ! Total of target array after adding redistributed data, and check
    
    ! DO NOT FIXME ?? : the following test will fail if tot2d << tot3da (strictly: if
    ! tot2d is less than the last significant digit of tot3da). However this
    ! is good, since we actually expect that tot3da to be 0 (as used in
    ! chemistry): so we can catch some bad initialization.
    tot3da = sum(e3d(:,j_st:j_end,:))
    if (abs((tot3da-tot3db)-tot2d) > 1e-8*abs(tot2d) ) then
      write (gol,'("emissions have not been distributed mass-conserving")'); call goErr
      write(gol,*)"totals to add : ", tot2d, tot3db ; call goErr
      write(gol,*)"total after : ", tot3da ; call goErr
      write(gol,*)"with bounds : ", j1, j2, j_st, j_end; call goErr
      write(gol,*)"shapes e3d",shape(e3d) ; call goErr
      write(gol,*)"lbound e3d",lbound(e3d); call goErr
      write(gol,*)"ubound e3d",ubound(e3d); call goErr
      write(gol,*)"minmax e3d",minval(e3d), maxval(e3d); call goErr
      write(gol,*)"shapes e2d",shape(e2d) ; call goErr
      write(gol,*)"lbound e2d",lbound(e2d); call goErr
      write(gol,*)"ubound e2d",ubound(e2d); call goErr
      write(gol,*)"minmax e2d",minval(e2d), maxval(e2d); call goErr

      TRACEBACK; status=1; return
    end if

    ! done
    nullify(gph)
    nullify(e2d)
    nullify(e3d)


  END SUBROUTINE DISTRIBUTE_EMIS2D
  !EOC

  !---------------------------------------------------------------------------
  !                    TM5                                                    !
  !---------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  distribute1x1
  !
  ! !DESCRIPTION: subroutine to distribute 1x1 emissions linearly 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.
  !               
  !               To be called by one processor, with a global 2D oro input
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DISTRIBUTE1X1( emi1x1, hlow, hhigh, emis3d, oro, status, xfrac)
    !
    ! !USES:
    !
    use Binas,     only : grav
    use dims,      only : nlon360, nlat180, lm, itau, staggered, at, bt
    !
    ! !INPUT PARAMETERS:
    !
    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), intent(in)            :: oro     ! (m m/s2) 
    real,                             intent(in), optional  :: xfrac   ! fraction of emissions to put
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    real, dimension(nlon360,nlat180,lm(1)), intent(inout)   :: emis3d  ! (kg/box) distributed in height
    integer,                                intent(out)     :: status
    !
    ! !REVISION HISTORY: 
    !    8 May 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
  
    character(len=*), parameter  :: rname = mname//'/distribute1x1'
  
    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 -----------------------------------

    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(i,j)/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
                status=1
                IF_NOTOK_RETURN(status=1)
             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 ) then
                status=1
                IF_NOTOK_RETURN(status=1)
             end if
             
          endif
       enddo
    enddo
    e3da = sum(emis3d)

    if (abs(e3da-e3db-sum(fract*emi1x1)) > e3da*1e-8 ) then
       status=1
       IF_NOTOK_RETURN(status=1)
    end if

    status=0
    
  END SUBROUTINE DISTRIBUTE1X1
  !EOC
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DISTRIBUTE1X1B
  !
  ! !DESCRIPTION: same as DISTRIBUTE1X1, but with scalar for HLOW and HHIGH
  !
  !               Is it still used ???? 
  !
  !                subroutine to distribute 1x1 emissions linearly 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!
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DISTRIBUTE1X1B( emi1x1, hlow, hhigh, emis3d, oro, status, xfrac)
    !
    ! !USES:
    !
    use Binas,     only : grav
    use dims,      only : nlon360, nlat180, lm, itau, staggered, at, bt
    !
    ! !INPUT PARAMETERS:
    !
    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), intent(in)            :: oro     ! (m m/s2) 
    real,                             intent(in), optional  :: xfrac   ! fraction of emissions to put
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    real, dimension(nlon360,nlat180,lm(1)), intent(inout)   :: emis3d  ! (kg/box) distributed in height
    integer,                                intent(out)     :: status
    !
    ! !REVISION HISTORY: 
    !    8 May 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    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(i,j)/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 
                status=1
                IF_NOTOK_RETURN(status=1)
             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 ) then
                status=1
                IF_NOTOK_RETURN(status=1)
             end if
             
          endif
       enddo
    enddo
    
    e3da = sum(emis3d)
    if (abs(e3da-e3db-sum(fract*emi1x1)) > e3da*1e-8 ) then
       status=1
       IF_NOTOK_RETURN(status=1)
    end if

    status=0

  END SUBROUTINE DISTRIBUTE1X1B
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    TROPOSPHERIC_COLUMNS
  !
  ! !DESCRIPTION: 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.
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine tropospheric_columns(region, field, slope, column, thres, xmo3, status)
    !
    ! !USES:
    !
    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 PARAMETERS:
    !
    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)
    !
    ! !OUTPUT PARAMETERS:
    !
    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

    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !   19 Jun 2012 - P. Le Sager - fix calls to CheckShape
    !
    ! !REMARKS:
    !  (1) FIXME : must be adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    character(len=*), parameter  :: rname = mname//'/tropospheric_columns'
    
#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

    ! FIXME
    write (gol,'("ERROR - routine not converted to TM5-MP")') ; call goErr
    status=1
    IF_NOTOK_RETURN(status=1)

    
    ! start
    call CheckShape( (/im(region),jm(region)/),             shape(column), status )
    call CheckShape( (/im(region),jm(region), lm(region)/), shape(field),  status )
    call CheckShape( (/im(region),jm(region), lm(region)/), shape(slope),  status )

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


END MODULE TOOLBOX