ecearth3/sources/tm5mp/proj/cb05/chemistry_cariolle.F90

#define TRACEBACK write (gol,'("in ",a," (",a,i6,")")') rname, __FILE__, __LINE__ ; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!#################################################################
!
! Parameterized ozone chemistry:
!  o Cariolle, Cold Tracer
!  o relaxation to climatology
!
!  Relaxation time scale is altitude depended:
! 
!      |                           ------ 0
!      |  ------                    ^
!      |   ^                        |
!      |   | ColdTracer             |
!      |   |                        |Cariolle
!      |   v                        | 
!      |  ------                    |
!      |                            v
!      |                      *    ------ cariolle_plev (230 hPa)
!      |              o  
!      |       *                   ------ tropo_plev (500 hPa)
!      |       o  
!      |       o
!    ----------+--------------+------> tau
!        tropo_tau_lower tropo_tau_upper
!           (weeks)        (months)
!  
!
!   22 Jun 2011 - P. Le Sager - Bug fix in O3 overhead. Plus added two
!                               checks on negative mixing ratio from Vincent.
!
!### macro's #####################################################
!
  
module Chemistry_Cariolle

  use GO          , only : gol, goPr, goErr, goLabel
  use dims        , only : nregions
  use global_types, only : chem_data
  use Climat      , only : TClimat
#ifdef with_budgets
  use budget_global,only : nbudg, nbud_vg
  use chem_param,   only : ncar
#endif
  
  implicit none
  
  
  ! --- in/out ----------------------------

  private

  public  ::  Cariolle_Init, Cariolle_Done, Cariolle_Chemie
  
  public  ::  o3clim_vmr, l_cariolle, with_trop_force
  
  
  ! --- const -----------------------------
  
  character(len=*), parameter  ::  mname = 'cariolle'
  
  !
  ! overhead fields
  !
  
  ! 3d fields with total column ABOVE each level !
  !
  !  rm_ovh(:,:,lm  ,iovh(k)) = 0.0
  !  rm_ovh(:,:,lm-1,iovh(k)) =      rm(:,:,lm,k)
  !             :
  !  rm_ovh(:,:,l   ,iovh(k)) = sum( rm(:,:,l+1:lm,k) , 3 )
  !
  ! Total column is thus :  
  !
  !  rm(:,:,1,k) + rm_ovh(:,:,1,iovh(k))
  !

  ! number of overhead fields:
  integer, parameter         ::  novh = 1
  ! indices:
  integer, parameter         ::  iovh_o3  = 1
  
  ! --- var ------------------------------
  
  ! overhead fields (4d arrays, parallel over layers)
  type(chem_data)              ::  ovh_dat(nregions)
  
  ! ozone climatology
  character(len=256)           ::  climat_fname
  type(TClimat),save           ::  o3clim_vmr(nregions)

  ! cariolle or linoz  coefficients
  character(len=256)           ::  cariolle_fname
  type(TClimat),save           ::  o3coef(nregions)
  type(TClimat),save           ::  o3coefp(nregions)
  integer, parameter           ::  nc = 7
  
  ! with cold tracer
  logical                      ::  with_cold_tracer

  ! with forcing troposphere to climatology
  logical                      ::  with_trop_force

  ! compute temper climat from model temperatures ?
  logical                      ::  apply_taver
  
  ! minimum time scale for slower chemistry:
  real                         ::  cariolle_tau_min

  ! pressure range over which Cariolle is applied:
  real                         ::  cariolle_plev 

  ! lowest model level for Cariolle chemistry:
  integer                      ::  l_cariolle(nregions)
  
  ! Pressure range over which cold tracer is active:
  !  o production of cold tracer
  !  o hetrochem 
  ! Loss of cold tracer is applied everywhere ...
  real                         ::  coldtracer_plevs(2) !  35.0e2,228.0e2  Pa

  ! tropospheric forcing to climatology:
  real                         ::  tropo_plev   ! Pa
  real                         ::  tropo_tau    ! sec

#ifdef with_budgets
  real,dimension(nbudg, nbud_vg, ncar ),public    :: budcarg
#endif


contains


  subroutine Cariolle_Init( status )
  
    use GO           , only : TrcFile, Init, Done, ReadRc
    use dims         , only : im, jm, lm
    use global_data  , only : rcfile
    use Climat       , only : Init
    use file_cariolle, only : Cariolle_Info
    use toolbox      , only : lvlpress
  
    ! --- in/out --------------------------------
    
    integer, intent(out)           ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Cariolle_Init'
    
    ! --- local ------------------------------
    
    type(TrcFile)   ::  rcF
    integer         ::  region
    integer         ::  imr, jmr, lmr
    real            ::  tau_day, plev_hPa
    integer         ::  nlevs

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

    call ReadRc( rcF, 'chem.o3tracer.climat', climat_fname, status )
    IF_NOTOK_RETURN(status=1)

    call ReadRc( rcF, 'chem.o3tracer.cariolle.coeff', cariolle_fname, status )
    IF_NOTOK_RETURN(status=1)

    call ReadRc( rcF, 'chem.o3tracer.cariolle.taver', apply_taver, status )
    IF_NOTOK_RETURN(status=1)

    ! pressure level below which Cariolle is applied:
    call ReadRc( rcF, 'chem.o3tracer.cariolle.strato.plev', plev_hPa, status )
    IF_NOTOK_RETURN(status=1)
    cariolle_plev    = plev_hPa * 1e2            ! Pa

    call ReadRc( rcF, 'chem.o3tracer.cariolle.strato.tau_min', tau_day, status )
    IF_NOTOK_RETURN(status=1)
    cariolle_tau_min = tau_day * 24.0 * 3600.0   ! sec


    ! with cold tracer?
    call ReadRc( rcF, 'chem.o3tracer.coldtracer', with_cold_tracer, status )
    IF_NOTOK_RETURN(status=1)

    ! pressure levels between which cold tracer scheme is applied:
    call ReadRc( rcF, 'chem.o3tracer.coldtracer.plevtop', coldtracer_plevs(1), status ) ! hPa
    IF_NOTOK_RETURN(status=1)
    call ReadRc( rcF, 'chem.o3tracer.coldtracer.plevbottom', coldtracer_plevs(2), status ) ! hPa
    IF_NOTOK_RETURN(status=1)
    coldtracer_plevs = coldtracer_plevs * 1.0e2   ! Pa


    call ReadRc( rcF, 'chem.o3tracer.tropforce', with_trop_force, status )
    IF_NOTOK_RETURN(status=1)

    ! tropospheric forcing to climatology:
    !   on and below plev0 (hPa), use relax timescale tau0 (days)
    !   between plev0 and plev1 (hPa), linear interpol between tau0 and tau1
    !
    call ReadRc( rcF, 'chem.o3tracer.cariolle.tropo.plev', plev_hPa, status )
    IF_NOTOK_RETURN(status=1)
    tropo_plev = plev_hPa * 1e2            ! Pa
    call ReadRc( rcF, 'chem.o3tracer.cariolle.tropo.tau', tau_day, status )
    IF_NOTOK_RETURN(status=1)
    tropo_tau  = tau_day  * 24.0 * 3600.0  ! sec

    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)
    
    !
    ! other
    !
    
    ! loop over regions:
    do region = 1, nregions
    
      ! regional grid sizes:
      imr = im(region)
      jmr = jm(region)
      lmr = lm(region)

      ! setup climatology; linear time interpolation
      call Init( o3clim_vmr(region), 'O3', 'ppv', 'linear', 1, jmr, lmr, status )
      IF_NOTOK_RETURN(status=1)

      ! setup cariole coeff and pressure; constant within month
      ! grid sizes:
      call Cariolle_Info( status, nlevs=nlevs )
      IF_NOTOK_RETURN(status=1)
      call Init( o3coef(region), 'coeff', 'unit', 'constant', nc, jmr, nlevs, status )
      IF_NOTOK_RETURN(status=1)
      call Init( o3coefp(region), 'pressure', 'Pa', 'constant', 1, jmr, nlevs, status )
      IF_NOTOK_RETURN(status=1)

      ! overhead fields (concentration arrays, thus 2 halo cells):
      allocate( ovh_dat(region)%rm_k(-1:imr+2,-1:jmr+2,lmr,novh) )

      ! lowest model level for Cariolle chemistry
      ! the use of 1013.25 hPa allows reading full level pressures from ECMWF tables
      l_cariolle=lvlpress(region,cariolle_plev, 101325.)
     
    end do   ! regions

    ! ok
    status = 0
    
  end subroutine Cariolle_Init
  
  
  ! ***
  
  
  subroutine Cariolle_Done( status )
  
    use Climat, only : Done

    ! --- in/out --------------------------------
    
    integer, intent(out)           ::  status
    
    ! --- const ------------------------------

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

    ! --- var ----------------------------------
    
    integer    ::  region

    ! --- begin --------------------------------
   
    ! loop over regions:
    do region = 1, nregions

      ! clear climatology:
      call Done( o3clim_vmr(region), status )
      IF_NOTOK_RETURN(status=1)

      ! clear cariolle coeff:
      call Done( o3coef(region), status )
      IF_NOTOK_RETURN(status=1)
      call Done( o3coefp(region), status )
      IF_NOTOK_RETURN(status=1)
     
      ! overhead fields:
      deallocate( ovh_dat(region)%rm_k )
     
    end do   ! regions

    ! ok
    status = 0
    
  end subroutine Cariolle_Done

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

  subroutine Cariolle_Chemie( region, tr, status )

    use binas       , only : Avog, xm_o3,xmair
    use GO          , only : TDate, Get, DayNumber, rTotal, wrtgol
    use GO          , only : operator(+), operator(-), operator(/)
    use toolbox     , only : dumpfield,escape_tm
    use Num         , only : Interp_Lin
    use Grid        , only : AreaOper
    use Phys        , only : cos_sza
    use dims        , only : im, jm, lm
    use dims        , only : isr, ier, jsr, jer
    use chem_param  , only : fscale, io3, ipsc,xmo3
    use global_data , only : mass_dat, region_dat
    use tracer_data , only : AdjustTracer
    use meteodata   , only : phlb_dat, m_dat
    use meteo       , only : lli, levi
    use meteo       , only : temper_dat
    use chemistry_0d, only : Cariolle_0d, ColdTracer_0d
    use ParTools    , only : tracer_active, lmloc, offsetl
    use ParTools    , only : which_par, previous_par, myid
    use photolysis_data, only : phot_dat
    use mpi_comm    , only : dump_field3d
    use boundary    , only : use_o3du
#ifdef with_budgets
    use budget_global,only : budg_dat,nzon_vg
    use chem_param  , only : ra
#endif


    implicit none

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

    integer, intent(in)          ::  region
    type(TDate), intent(in)      ::  tr(2)
    integer, intent(out)         ::  status

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

    character(len=*), parameter ::  rname = mname//'/Chemie'
!    real, parameter  :: conv = 0.1*Avog/xm_o3   !  from kg/m2 --> #/cm2
    real, parameter  :: conv = 0.1*Avog/xmo3   !  from kg/m2 --> #/cm2
    real, parameter  :: todu = 3.767e-17        !  from #/cm2 --> du

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

    type(TDate)         ::  tmid

    integer             ::  imr, jmr, lmr
    integer             ::  i, j, l, ll, lll
    integer             ::  level, lglob
#ifdef with_budgets   
    integer             ::  nzone,nzone_v
#endif

    real, pointer       ::  rm(:,:,:,:)
    real, pointer       ::  vo3(:,:,:)
    real, pointer       ::  ozone_klimtop(:)
    real, pointer       ::  phlb(:,:,:)
    real, pointer       ::  m(:,:,:)
    real, pointer       ::  temper(:,:,:)
#ifdef with_zoom
    integer, pointer    ::  zoomed(:,:)
#endif
    real, pointer       ::  area(:)        ! cell area (m2) for each lat band

    real, allocatable   ::  Tclim(:,:)
!    real,allocatable    ::  o3_overhead(:,:,:) ! in mlc/cm2
    real,allocatable    ::  o3_overhead(:,:) ! in mlc/cm2

    real                ::  lon, lat
    integer             ::  daynr, hour, minu, sec
    real                ::  csza

    real                ::  pf    ! full level pressure (Pa)

    real                ::  dt_sec
    real                ::  o3_vmr    ! ozone mixing ratio (ppv)
    real                ::  o3_ohc    ! overhead column (mlc/cm2)
    real                ::  o3_ohc_fac
    real                ::  o3c_vmr   ! ozone mixing ratio (ppv)
    real                ::  o3c_ohc   ! overhead column (mlc/cm2)
    real                ::  kg_TO_mlc_m2
    real                ::  Xpsc, dXpsc, kXpsc
    real                ::  do3_vmr, do3

    real                ::  cc(nc)
    integer             ::  ilast
    real                ::  k_h

    real                ::  rm_new,A,B,sgn,timesc

!type(TTimeSeriesHDF)   ::  F

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

    if ( lmloc == 0 ) then
      status=0; return
    end if

#ifdef MPI
    which_par = previous_par(region)
    if (which_par == 'tracers') &
       call escape_tm('Model should be parallel over levels in Cariolle chemistry')
#endif

    call goLabel( rname )

    if (myid == 0) then
    
      call wrtgol( 'ozone parchem : ', tr(1), ' - ', tr(2) ); call goPr

      if (      with_cold_tracer ) then
        write (gol,*) 'ozone parchem with cold tracer' ; call goPr
      else
        write (gol,*) 'ozone parchem without cold tracer' ; call goPr
      endif


      if (      with_trop_force ) then
        write (gol,*) 'ozone parchem with nudged tropospheric climatology' ; call goPr
      else
        write (gol,*) 'ozone parchem without nudged tropospheric climatology' ; call goPr
      endif 

    end if

    ! regional grid sizes:
    imr = im(region)
    jmr = jm(region)
    lmr = lm(region)

      ! set pointers
#ifdef MPI
    rm     => mass_dat(region)%rm_k       ! kg tracer, parallel over levels
#else
    rm     => mass_dat(region)%rm_t       ! kg tracer, parallel over tracers
#endif
    m      => m_dat(region)%data          ! kg air
    phlb   => phlb_dat(region)%data       ! air pressure at boundaries (1:lm+1)
    vo3    => phot_dat(region)%vo3        ! overhead O3, #/cm2, parallel over levels
    ozone_klimtop => phot_dat(region)%o3klim_top ! top layers overhead O3, #/cm2
    temper => temper_dat(region)%data     ! K
    area   => region_dat(region)%dxyp     ! m2
#ifdef with_zoom
    zoomed => region_dat(region)%zoomed
#endif

    allocate(o3_overhead(jmr,lmr)) ; o3_overhead = 0.0

    ! mid of interval:
    tmid = tr(1) + (tr(2)-tr(1))/2

    !
    ! ** setup climatology
    !

    ! interpolate to this time, eventually read new fields:
    call Setup_O3Clim( lli(region), levi, O3clim_vmr(region), tmid, status )
    IF_NOTOK_RETURN(status=1)

    if (use_o3du) then  
    ! overhead ozone column for model top is calculated (in photolysis routine)
    ! based on Fortuin-Kelder climatology scaled by observations
    ! same should be done here for climatological overhead column 
      do j=1,jmr 
        o3_overhead(j,lmr) = ozone_klimtop(j)
      enddo
    else
       ! convert climatology from ppv to mlc/cm2 using standard airmass/m2 (p=1e5 Pa)

!PRIOR
!       o3_overhead(:,:,lmr) = 0.5 * O3clim_vmr(region)%data(:,:,lmr)*levi%m0(lmr) * conv / fscale(io3)
!NOW
       o3_overhead(:,lmr) = 0.5 * O3clim_vmr(region)%data(1,:,lmr)*levi%m0(lmr) * conv / fscale(io3)
       
    endif
   
    do l = lmr-1,1,-1
       ! convert climatology from ppv to mlc/cm2 using standard airmass/m2 (p=1e5 Pa)
       
!PRIOR
!        o3_overhead(:,:,l) = o3_overhead(:,:,l+1) + &
!                             0.5 * ( O3clim_vmr(region)%data(:,:,l)*levi%m0(l) + &
!                                     O3clim_vmr(region)%data(:,:,l+1)*levi%m0(l+1) ) * conv / fscale(io3)
!NOW        
       o3_overhead(:,l) = o3_overhead(:,l+1) + &
                            0.5 * ( O3clim_vmr(region)%data(1,:,l)*levi%m0(l) + &
                                    O3clim_vmr(region)%data(1,:,l+1)*levi%m0(l+1) ) * conv / fscale(io3)
    end do

    !
    ! ** setup coefficients
    !

    ! tables:
    !  bsf15k:/usr/people/eskes/fdscia/input/cariolle/OZONE_64x26.ASCII
    !                                        linoz/linoz_coeff.dat
    ! input routines:
    !   /usr/people/segers/work/tm/TM5/proj/chem/ozone.bak/src/chem_ozonepar.F90

    ! check if current coeff are valid for this month, 
    ! eventually (re)load:
    call Setup_Cariolle( lli(region), o3coef(region), o3coefp(region), tmid, status )
    IF_NOTOK_RETURN(status=1)

    !
    ! zonal average temperature
    !

    ! compute smoothed temperature climatology ?
    if ( apply_taver ) then

      ! setup zonal field:
      allocate( Tclim(jmr,lmr) )

      ! fill smoothed temperature:
      call Setup_Tclim( lli(region), temper, Tclim, status )
      IF_NOTOK_RETURN(status=1)

    end if

    !
    ! apply chemistry
    !

    ! time step
    dt_sec = rTotal( tr(2) - tr(1), 'sec' )

    ! loop over all grid cells:
    do j = jsr(region), jer(region)
      do i = isr(region), ier(region)

        ! skip if this cell is calculated in an overlying zoom region:
#ifdef with_zoom
        if ( region_dat(region)%zoomed(i,j) /= region ) cycle
#endif

        ! grid cell location:
        lon = lli(region)%lon_deg(i)   ! deg
        lat = lli(region)%lat_deg(j)   ! deg

        ! compute cos(solar_zenith_angle)
        daynr = DayNumber( tmid )
        call Get( tmid, hour=hour, min=minu, sec=sec )
        csza = cos_sza( daynr, hour, minu, sec, lon, lat )

        ! conversion from (kg o3) to (mlc o3/m2) :
        !    ( kg O3 / m2    ) / (kg o3 / mol) * (mlc/mol) * (1e-4 m2/cm2) = mlc/cm2
        kg_TO_mlc_m2 =   1.0  / area(j) /     xm_o3     *    Avog   *     1.0e-4

        ! loop over levels
        ilast = 1
        do level = 1,lmloc

          lglob = offsetl + level

          ! Cariolle chemistry not applied below l_cariolle
          ! except when tropospheric nudging is active
          if ( (.not.with_trop_force) .and. (lglob .lt. l_cariolle(region)) ) exit

          ! full level pressure:
          pf = ( phlb(i,j,lglob) + phlb(i,j,lglob+1) )/2.0    ! Pa

          ! ozone volume mixing ratios:

          o3_vmr  = fscale(io3) * rm(i,j,level,io3) / m(i,j,lglob)    ! ppv
	  o3_vmr  = max(o3_vmr,0.0) ! can become negative...
          o3c_vmr = o3clim_vmr(region)%data(1,j,lglob)

          ! overhead column (including contribution of upper-half of current layer):
          o3_ohc  = vo3(i,j,level)  ! mlc/cm2

          o3c_ohc = o3_overhead(j,lglob) ! mlc/cm2

          o3_ohc_fac = 0.5*kg_TO_mlc_m2 / fscale(io3) * m(i,j,lglob)
          o3_ohc_fac = max(o3_ohc_fac,0.0)

!if (i.eq.10.and.j.eq.10)then
! write(*,'a15,3i3,5e10.2')'DEBUG CARIOLLE',i,j,lglob,o3_ohc*todu,o3c_ohc*todu,o3_vmr*o3_ohc_fac*todu,o3_vmr,o3_ohc_fac
!endif
          ! ** cold tracer scheme
          ! standard put values to zero, overwritten if with_cold_tracer = T

          Xpsc = 0. 

          ! convert to mass again:
          dXpsc = 0. 

          ! hetero chem (might be zero):
          k_h = 0. !kXpsc

          if ( with_cold_tracer ) then 

            ! cold tracer (marked air mass)
            Xpsc = rm(i,j,level,ipsc) / m(i,j,lglob)      ! ~ 0-1
            Xpsc = min( max( 0.0, Xpsc ), 1.0 )   ! 0-1

            ! hetrogeneous chemistry (cold tracer);
            ! fills kXpsc to be used as hetrochem reaction rate in Cariolle scheme:

            call ColdTracer_0d( Xpsc, pf, coldtracer_plevs, &
                                 temper(i,j,lglob), csza, lat, dt_sec, &
                                 dXpsc, kXpsc, status )
            IF_NOTOK_RETURN(status=1)

            ! convert to mass again:
            dXpsc = dXpsc * m(i,j,lglob)   ! kg

            ! add tracer mass changes:
            if (dXpsc .gt. 0.) then
              call AdjustTracer( dXpsc, region, i, j, level, ipsc, status )
              IF_NOTOK_RETURN(status=1)
            end if

#ifdef with_budgets
            nzone=budg_dat(region)%nzong(i,j)    ! global budget
            nzone_v=nzon_vg(lglob)
            budcarg(nzone,nzone_v,2)=budcarg(nzone,nzone_v,2)+ dXpsc*1000./ra(ipsc) ! mol psc per cell
#endif

           ! hetero chem (might be zero):
            k_h = kXpsc

          end if

          ! ** end applying cold tracer

          ! ** Cariolle scheme

          ! interpolate to pressure level:
          call Interp_Lin( o3coefp(region)%data(1,j,:), &
                             o3coef(region)%data(:,j,:), pf, cc, ilast, status )
          IF_NOTOK_RETURN(status=1)

          ! replace some coeffs by climatological data
          ! (have no impact if corresponding reaction rates are zero)
          cc(3) = o3c_vmr     ! zonal aver ozone mixing ratio
          cc(7) = o3c_ohc     ! zonal aver ozone overhead column

          ! use zonal aver temperature ?
          ! (has no impact if corresponding reaction rate is zero)
          if ( apply_taver ) cc(5) = Tclim(j,lglob)

          ! apply parameterized ozone chemistry for current cell:
          call Cariolle_0d( o3_vmr, o3_ohc, o3_ohc_fac, temper(i,j,lglob), &
                              cc, k_h, cariolle_tau_min, dt_sec, do3_vmr )
          IF_NOTOK_RETURN(status=1)

!          if (i.eq.5.and. (j.eq.1 .or. j.eq.11 .or. j.eq. 80 .or. j.eq.90))then
! write(*,'a18,i3,i3,8e11.3')'DEBUG CARIOLLE cc',lglob,j,cc,dt_sec
! write(*,'a18,i3,i3,7e11.3')'DEBUG CARIOLLE o3',lglob,j,o3_vmr,o3_ohc,o3_ohc_fac,temper(i,j,lglob),pf,k_h,do3_vmr
! 
! write(*,'a18,i3,i3,7e11.3')'DEBUG CAR timesc', lglob, j, 1./abs(cc(2) + cc(6)*o3_ohc_fac - k_h)
!       endif
          
          ! ozone mass change:
          do3 = do3_vmr * m(i,j,lglob) / fscale(io3)   ! kg

          ! add tracer mass changes:

          call AdjustTracer(   do3, region, i, j, level,  io3, status )
          IF_NOTOK_RETURN(status=1)

#ifdef with_budgets
          nzone=budg_dat(region)%nzong(i,j)    ! global budget
          nzone_v=nzon_vg(lglob)
          budcarg(nzone,nzone_v,1)=budcarg(nzone,nzone_v,1)+do3*1000./ra(io3) ! mol o3 per cell
#endif

        end do  ! levels

      end do  ! i
    end do  ! j

    ! 
    ! done
    !

    nullify( rm     )
    nullify( m      )
    nullify( phlb   )
    nullify( vo3    )
    nullify( ozone_klimtop )
    nullify( temper )
#ifdef with_zoom
    nullify( zoomed )
#endif
    nullify( area   )

#ifdef MPI
    rm     => mass_dat(region)%rm_k       ! kg tracer, parallel over levels
#else
    rm     => mass_dat(region)%rm_t       ! kg tracer, parallel over tracers
#endif
    m      => m_dat(region)%data          ! kg air

    nullify( rm     )
    nullify( m      )

    if ( apply_taver ) deallocate( Tclim )
    deallocate (o3_overhead)

    ! ok
    call goLabel(); status=0

  end subroutine Cariolle_Chemie


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


  subroutine Setup_Cariolle( lli, carcoef, carcoefp, t, status )
  
    use Binas        , only : p0    ! standard pressure (1e5 Pa)
    use GO           , only : TDate, Get, NewDate
    use Num          , only : Interp_Lin
    use Grid         , only : TllGridInfo
    use Climat       , only : TClimat, Set, Setup
    use file_cariolle, only : Cariolle_Info, CARIOLLE_READ
  
    ! --- in/out --------------------------------
    
    type(TllGridInfo), intent(in)   ::  lli
    type(TClimat), intent(inout)    ::  carcoef    ! 1-7
    type(TClimat), intent(inout)    ::  carcoefp   ! Pa
    type(TDate), intent(in)         ::  t
    integer, intent(out)            ::  status
    
    ! --- const ------------------------------

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

    ! --- local ------------------------------
    
    integer             ::  year, month
    integer             ::  nlats, nlevs
    real, allocatable   ::  lats(:)
    real, allocatable   ::  pp(:,:)
    real, allocatable   ::  ppX(:,:,:)
    real, allocatable   ::  cc(:,:,:)
    real, allocatable   ::  ccX(:,:,:)
    type(TDate)         ::  tr_in(2)
    integer             ::  j, l, ic
    
    ! --- begin ------------------------------
    
    ! try to setup for requested time;
    ! return status=-1 if not filled yet or filled for wrong time:
    call Setup( carcoef, t, status )
    if (status==0) return
    
    ! data needs to be (re)setup ...
    
    ! grid sizes:
    call Cariolle_Info( status, nlats=nlats, nlevs=nlevs )
    IF_NOTOK_RETURN(status=1)
    
    ! setup storage:
    allocate( lats(nlats) )
    allocate( pp(nlats,nlevs) )
    allocate( ppX(1,lli%nlat,nlevs) )
    allocate( cc(nlats,nlevs,nc) )
    allocate( ccX(nc,lli%nlat,nlevs) )
    
    ! extract month:
    call Get( t, year=year, month=month )

    ! load coeff for this month (ppmv)
    call Cariolle_Read( trim(cariolle_fname), month, lats, pp, cc, status )
    IF_NOTOK_RETURN(status=1)
    
    ! interpolate pressure field to latitudes
    do l = 1, nlevs
      call Interp_Lin( lats, pp(:,l), lli%lat_deg, ppX(1,:,l), status  )
      IF_NOTOK_RETURN(status=1)
    end do

    ! spatial interpolation;
    ! loop over coeffs:
    do ic = 1, nc

      ! interpolate to latitudes
      do l = 1, nlevs
        call Interp_Lin( lats, cc(:,l,ic), lli%lat_deg, ccX(ic,:,l), status )
        IF_NOTOK_RETURN(status=1)
      end do

      !! interpolate to standard pressures
      !do j = 1, lli%nlat
      !  ! Cariolle coefs from top->surf, thus no negation required:
      !  call Interp_Lin( ppX(j,:), ccX(j,:), levi%fp0, cc_in(ic,j,:) )
      !end do
      
    end do   ! coeff loop
    
    ! set time range for which this data is valid:
    !  o begin of current month:
    tr_in(1) = NewDate( year, month, 01, 00, 00 )
    !  o begin of next month:
    month = month + 1
    if ( month > 12 ) then
      year = year + 1
      month = 1
    end if
    tr_in(2) = NewDate( year, month, 01, 00, 00 )
    
    ! store:
    call Set( carcoef, status, data=ccX, tr=tr_in )
    IF_NOTOK_RETURN(status=1)
    call Set( carcoefp, status, data=ppX, tr=tr_in )
    IF_NOTOK_RETURN(status=1)
      
    ! set coeff such that in troposphere only relaxation to climatology remains:
    
    if ( with_trop_force ) then    

      call Setup_TropForce( carcoef, carcoefp, status )
      IF_NOTOK_RETURN(status=1)
    
    end if 

    ! clear storage:
    deallocate( lats )
    deallocate( pp   )
    deallocate( ppX  )
    deallocate( cc   )
    deallocate( ccX  )
    
    ! ok
    status = 0

  end subroutine Setup_Cariolle


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

  subroutine Setup_TropForce( carcoef, carcoefp, status )
  
    ! --- in/out --------------------------------
    
    type(TClimat), intent(inout)    ::  carcoef     ! c1-c7
    type(TClimat), intent(inout)    ::  carcoefp    ! Pa
    integer, intent(out)            ::  status
    
    ! --- const ------------------------------

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

    ! --- local ------------------------------
    
    integer             ::  nlat, nlev
    integer             ::  j, l
    integer             ::  ilev0
    real                ::  plev0, tau0
    real                ::  plev, tau
    real                ::  cc(7)
    real                ::  alfa
    
    ! --- begin ------------------------------
    
    ! grid size
    nlat = size(carcoefp%data,2)
    nlev = size(carcoefp%data,3)
    
    ! loop over latitudes
    do j = 1, nlat
    
      ! first level above troposphere:
      if ( carcoefp%data(1,j,1) > carcoefp%data(1,j,nlev) ) then
        ! upwards
        ilev0 = nlev   ! model top
        do l = 1, nlev
          if ( carcoefp%data(1,j,l) <= cariolle_plev ) then
            ilev0 = l
            exit
          end if
        end do
      else
        ! downwards
        ilev0 = 1    ! model top
        do l = nlev, 1, -1
          if ( carcoefp%data(1,j,l) <= cariolle_plev ) then
            ilev0 = l
            exit
          end if
        end do
      end if
      ! pressure in lowest cariolle level:
      plev0 = carcoefp%data(1,j,ilev0)   ! Pa
      tau0 = - 1.0/carcoef%data(2,j,ilev0)     ! sec
    
      ! loop over levels
      do l = 1, nlev
      
        ! current pressure and coeffs:
        cc   = carcoef%data(:,j,l)
        plev = carcoefp%data(1,j,l)
        
        ! different pressure ranges:
        if ( plev >= tropo_plev ) then

          ! ~~ relaxation to climatology, slow

          ! relaxation towards climatology via Cariolle coeff:
          cc    = 0.0
          cc(2) = -1.0/tropo_tau

        else if ( plev >= cariolle_plev ) then

          ! ~~ relaxation to climatology, slower downwards

          ! interpolate tau between lower constant value and cariolle value:
          alfa = ( plev - plev0 ) / ( tropo_plev - plev0 )
          tau = tau0 * (1.0-alfa ) + tropo_tau * alfa

          ! relaxation towards climatology via Cariolle coeff:
          cc    = 0.0
          cc(2) = -1.0/tau

        end if

        ! restore:
        carcoef%data(:,j,l) = cc

      end do   ! levs
      
    end do   ! lats  
    
    ! ok
    status = 0

  end subroutine Setup_TropForce


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

  subroutine Setup_O3clim( lli, levi, o3c, t, status )
  
    use GO           , only : TDate, Get, NewDate
    use Num          , only : Interp_Lin
    use Grid         , only : TllGridInfo, TLevelInfo
    use Climat       , only : TClimat, Set, Setup
    use file_fortkeld, only : FortuinKelder_Info, FortuinKelder_Read
  
    ! --- in/out --------------------------------
    
    type(TllGridInfo), intent(in)   ::  lli
    type(TLevelInfo), intent(in)    ::  levi
    type(TClimat), intent(inout)    ::  o3c
    type(TDate), intent(in)         ::  t
    integer, intent(out)            ::  status
    
    ! --- const ------------------------------

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

    ! --- local ------------------------------
    
    integer             ::  year, month, day
    type(TDate)         ::  t_in(2)
    integer             ::  it
    
    integer             ::  nlats, nlevs
    real, allocatable   ::  lats_deg(:)
    real, allocatable   ::  pres_Pa(:)
    real, allocatable   ::  o3_ppv(:,:)
    real, allocatable   ::  o3X(:,:)
    real, allocatable   ::  sp(:)
    real, allocatable   ::  o3_in(:,:,:)
    integer             ::  j, l
    
    ! --- begin ------------------------------
    
    ! try to interpolate in time to t;
    ! return status=-1 if not filled yet or filled for wrong time:
    call Setup( o3c, t, status )
    if (status==0) return
    
    ! data needs to be (re)setup ...
    
    ! extract grid sizes:
    call FortuinKelder_Info( status, nlats=nlats, nlevs=nlevs )
    IF_NOTOK_RETURN(status=1)
    
    ! setup arrays:
    allocate( lats_deg(nlats)     )
    allocate( pres_Pa(nlevs)      )
    allocate( o3_ppv(nlats,nlevs) )
    allocate( o3X(lli%nlat,nlevs) )
    allocate(  sp(lli%nlat)       )
    allocate( o3_in(1,lli%nlat,levi%nlev) )
    
    ! loop over interpolation times:
    do it = 1, 2
    
      ! set times to 16th 00:00
      call Get( t, year=year, month=month, day=day )
      if ( it == 1 ) then
        ! previous mid of month
        if ( day <= 15 ) then
          ! mid of previous month
          month = month - 1
          if ( month < 1 ) then
            year = year - 1
            month = 12
          end if
        end if
      else
        ! next mid of month
        if ( day >= 16 ) then
          ! mid of next month
          month = month + 1
          if ( month > 12 ) then
            year = year + 1
            month = 1
          end if
        end if
      end if
      
      ! load climatology for this month (ppmv)
      call FortuinKelder_Read( trim(climat_fname), month, &
                                 lats_deg, pres_Pa, o3_ppv, status )
      IF_NOTOK_RETURN(status=1)
      
      ! convert from ppmv to ppv:
      o3_ppv = o3_ppv * 1.0e-6

      ! interpolate to latitudes
      do l = 1, nlevs
        call Interp_Lin( lats_deg, o3_ppv(:,l), lli%lat_deg, o3X(:,l), status )
	IF_NOTOK_RETURN(status=1)
      end do

      ! interpolate to standard pressures
      do j = 1, lli%nlat
        ! negate pressure ax; should be increasing ...
        call Interp_Lin( -pres_Pa, o3X(j,:), -levi%fp0, o3_in(1,j,:), status )
	IF_NOTOK_RETURN(status=1)
      end do
      
      ! time is mid of month:
      t_in(it) = NewDate( year, month, 16, 00, 00 )
      
      ! store:
      call Set( o3c, status, t_in=t_in(it), data_in=o3_in, it_in=it )
      IF_NOTOK_RETURN(status=1)
    
    end do  ! times
    
    ! clear:
    deallocate( lats_deg )
    deallocate( pres_Pa  )
    deallocate( o3_ppv   )
    deallocate( o3X      )
    deallocate( sp       )
    deallocate( o3_in    )
    
    ! set time range for which data could be interpolated:
    call Set( o3c, status, tr=t_in )
    IF_NOTOK_RETURN(status=1)
    
    ! apply interpolation: should succeed now:
    call Setup( o3c, t, status )
    if ( status < 0 ) then
      write (gol,'("time interpolation failed, while data has just been read.")'); call goErr
      write (gol,'("in ",a)') rname; call goErr; status=1; return
    else if ( status > 0 ) then
      write (gol,'("in ",a)') rname; call goErr; status=1; return
    end if

    ! ok
    status = 0

  end subroutine Setup_O3clim

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

  subroutine Setup_Tclim( lli, T, Taver, status )
  
    use Grid, only : TllGridInfo
  
    ! --- in/out --------------------------------
    
    type(TllGridInfo), intent(in)   ::  lli
    real, intent(in)                ::  T(:,:,:)
    real, intent(out)               ::  Taver(:,:)
    integer, intent(out)            ::  status
    
    ! --- const ------------------------------

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

    ! latitude range for model temperature average
    real, parameter  ::  latrange = 6.2

    ! --- local ------------------------------
    
    integer                   ::  im, jm, lm
    integer                   ::  i, j, l, k
    integer                   ::  kmin, kmax
    real                      ::  wk
    real                      ::  taccu, waccu
    
    ! --- begin ------------------------------
    
    im = size(T,1)
    jm = size(T,2)
    lm = size(T,3)
    
    do l = 1, lm
      do j = 1, jm

        taccu = 0.0
        waccu = 0.0

        ! cells within lat+[-latrange,+latrange]
        kmin = max(j-int(latrange/lli%dlat_deg),1)
        kmax = min(j+int(latrange/lli%dlat_deg),jm)

        ! loop over smoothing cells:
        do k = kmin, kmax
           ! weight decreasing with distance:
           wk = 1.0 - abs(real(k-j))*lli%dlat_deg/latrange
           if ( wk <= 0.0 ) then
             write (gol,'("internal error : wk = ",es12.4)') wk
             write (gol,'("in ",a)') rname; call goErr; status=1; return
           end if
           ! accumulate
           do i = 1, im
             taccu = taccu + wk * T(i,k,l)
             waccu = waccu + wk
           end do   ! lon
        end do   ! smoothing lats

        ! average:
        Taver(j,l) = taccu / waccu
        
      end do   ! lats
    end do   ! levs

    ! ok
    status = 0

  end subroutine Setup_Tclim


end module Chemistry_Cariolle