!#################################################################
!
! tracer and air mass
!
!### 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"
!
#define IF_MPI_NOTOK_RETURN(action) if (status/=MPI_SUCCESS) then; TRACEBACK; action; return; end if
!
!#################################################################

module tracer_data

  use GO, only : gol, goPr, goErr
  
  use GO          , only : gol, goErr, goPr
  use dims        , only : nregions
  use global_types, only : chem_data

#ifdef with_tendencies
  use GO          , only : TDate
  use chem_param  , only : ntrace
  use tm5_tendency, only : plc_ntr, plc_trname
  use tm5_tendency, only : plc_itr_nox
  use tm5_tendency, only : plc_npr, plc_prname
  use tm5_tendency, only : plc_ipr_pchem
  use tm5_tendency, only : plc_ipr_lchem, plc_ipr_pchem, plc_ipr_pemi , plc_ipr_conc, &
                           plc_ipr_lddep, plc_ipr_lwdep, plc_ipr_tcnvd, plc_ipr_tchem
#endif

  implicit none
  
  ! --- in/out ------------------------------------
  
  private
  
  public  ::  mass_data, mass_dat
  public  ::  assign_massdat, free_massdat
  
  public  ::  chem_data, chem_dat

  public  ::  Par_Check_Domain
  public  ::  Par_Check_Mass

  public  ::  AdjustTracer
  
  public  ::  tracer_print

#ifdef with_mpi
  public  ::  Tracer_Get_Here
#endif
#ifdef with_feedback
  public  ::  Tracer_Fill_Slabs
#endif
  
#ifdef with_tendencies
  public  ::  plc_ntr
  public  ::  plc_npr, plc_ipr_pchem, plc_ipr_lchem, plc_ipr_pemi , plc_ipr_conc, &
                       plc_ipr_lddep, plc_ipr_lwdep, plc_ipr_tcnvd, plc_ipr_tchem
  public  ::  plc_kg_from_tm
  public  ::  plc_data, plc_dat
  public  ::  PLC_Init, PLC_Done, PLC_Set, PLC_Reset
  public  ::  PLC_AddValue, PLC_AddValue_t
  public  ::  plc_reset_period
#endif


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

  character(len=*), parameter  ::  mname = 'tracer_data'


  ! --- types ------------------------------------

  type mass_data

     ! m_k     air mass (kg)   (halo = 2)
     ! p       surface pressure (Pa)   (halo = 2)
     ! phlb_k  pressure at the grid boundaries (Pa) : 
     !         l=1 means bottom lowest box (halo = 0)
     ! rm_k    tracer masses (kg)    (halo = 2)
     ! rxm_k   tracer slopes (only transported) in kg/(halfgridsize) (halo = 1)
     ! rym_k   tracer slopes (only transported) in kg/(halfgridsize) (halo = 1)
     ! rzm_k   tracer slopes (only transported) in kg/(halfgridsize) (halo = 1)
     ! m_t     air mass (kg)   (halo = 2)
     ! phlb_t  pressure at the grid boundaries (Pa) : 
     !         l=1 means bottom lowest box (halo = 0)
     ! rm_t    tracer masses (kg)    (halo = 2)
     ! rxm_t   tracer slopes (only transported) in kg/(halfgridsize) (halo = 1)
     ! rym_t   tracer slopes (only transported) in kg/(halfgridsize) (halo = 1)
     ! rzm_t   tracer slopes (only transported) in kg/(halfgridsize) (halo = 1)
     !
     real,dimension(:,:,:,:),pointer   :: rm_k
#ifdef slopes
     real,dimension(:,:,:,:),pointer   :: rxm_k
     real,dimension(:,:,:,:),pointer   :: rym_k
     real,dimension(:,:,:,:),pointer   :: rzm_k
#ifdef secmom
     real,dimension(:,:,:,:),pointer   :: rxxm_k
     real,dimension(:,:,:,:),pointer   :: rxym_k
     real,dimension(:,:,:,:),pointer   :: rxzm_k
     real,dimension(:,:,:,:),pointer   :: ryym_k
     real,dimension(:,:,:,:),pointer   :: ryzm_k
     real,dimension(:,:,:,:),pointer   :: rzzm_k
#endif
#endif
     real,dimension(:,:,:,:),pointer   :: rm_t 
#ifdef slopes
     real,dimension(:,:,:,:),pointer   :: rxm_t
     real,dimension(:,:,:,:),pointer   :: rym_t
     real,dimension(:,:,:,:),pointer   :: rzm_t
#ifdef secmom
     real,dimension(:,:,:,:),pointer   :: rxxm_t
     real,dimension(:,:,:,:),pointer   :: rxym_t
     real,dimension(:,:,:,:),pointer   :: rxzm_t
     real,dimension(:,:,:,:),pointer   :: ryym_t
     real,dimension(:,:,:,:),pointer   :: ryzm_t
     real,dimension(:,:,:,:),pointer   :: rzzm_t
#endif
#endif

  end type mass_data


#ifdef with_tendencies

  ! 3D production or loss rates, or concentrations
  type plc_data
    ! timing
    type(TDate)           ::  t
    ! slab of levels
    real, pointer         ::  rm_k(:,:,:)
    ! temporary storage for budgets if proces is parallel over tracers:
    logical               ::  par_t
    real, pointer         ::  rm_t(:,:,:)
    ! unit
    character(len=32)     ::  unit
  end type plc_data

#endif


  ! --- var ------------------------------

  !
  ! transported tracer masses and air mass:
  !
  type(mass_data), target   :: mass_dat(nregions)

  !
  ! non-transported tracer masses:
  !  o used to store emissions ?
  !  o allocated in trace0, deallocated in exitus
  !      chem_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,ntracet+1:ntracet+ntrace_chem)
  !  o read in in trace0
  !  o applied in source1
  !
  type(chem_data), target   :: chem_dat(nregions) 

#ifdef with_tendencies

  ! mapping of indices:
  integer                         ::  plc_itr_to_tm(plc_ntr)
  integer                         ::  plc_itr_from_tm(ntrace)
  
  ! mixing from tm tracer mass to tendency tracer mass:
  real                            ::  plc_kg_from_tm(ntrace)

  ! concentration, production and loss rates
  ! note: extra proces '0' is used to store concentrations at begin of a step
  type(plc_data), target, save    ::  plc_dat(nregions,plc_ntr,0:plc_npr)
  
  ! reset period in hours:
  integer                         ::  plc_reset_period  ! hours

#endif


contains


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


#ifdef with_tendencies  
  
  subroutine PLC_Init( rcfile, status )
  
    use GO          , only : AnyDate
    use GO          , only : TrcFile, Init, Done, ReadRc
    use dims        , only : nregions
    use dims        , only : im, jm, lm
    use chem_param  , only : ntracet
    use chem_param  , only : iO3, iCO, iNO2, iNOx, iCH2O, iSO2
    use chem_param  , only : xmNOx, xmNO2
    use tm5_tendency, only : plc_ipr_conc, plc_ipr_pchem, plc_ipr_lchem
    use tm5_tendency, only : plc_trname
    use partools    , only : lmloc, tracer_active
    
    ! --- in/out -------------------------------
    
    character(len=*), intent(in)   ::  rcfile
    integer, intent(out)           ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/PLC_Init'
    
    ! --- local --------------------------------
    
    type(TrcFile)   ::  rcF
    integer         ::  n
    integer         ::  imr, jmr, lmr
    integer         ::  itr, itr_tm
    integer         ::  ipr
    
    ! --- begin -------------------------------
    
    !
    ! settings
    !
    
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    call ReadRc( rcF, 'tend.reset.period', plc_reset_period, status )
    IF_ERROR_RETURN(status=1)

    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)

   
    
    !
    ! indices mapping from plc to tm and vs
    !
    
    ! initialise with dummy values:
    plc_itr_to_tm      = -1
    plc_itr_from_tm    = -1

    ! loop over plc tracers:
    do itr = 1, plc_ntr

      ! corresponding global tracer index and conversion factor:
      plc_kg_from_tm = 1.0
      select case ( plc_trname(itr) )
        case ( 'o3'   ) ; itr_tm = iO3
        case ( 'co'   ) ; itr_tm = iCO
        case ( 'no2'  ) ; itr_tm = iNO2
        case ( 'nox'  ) ; itr_tm = iNOx   ; plc_kg_from_tm(itr_tm) = xmNO2/xmNOx  ! kg N -> kg NO2
        case ( 'ch2o' ) ; itr_tm = iCH2O
        case ( 'so2'  ) ; itr_tm = iSO2
        case default
          write (gol,'("unsupported plc tracer : ",a)') trim(plc_trname(itr)); call goErr
          TRACEBACK; status=1; return
      end select
      
      ! store:
      plc_itr_to_tm  (itr   ) = itr_tm
      plc_itr_from_tm(itr_tm) = itr
      
    end do   ! plc tracers
    
    !
    ! allocate
    !
    
    ! loop over regions:
    do n = 1, nregions
    
      ! local grid sizes:
      imr = im(n)
      jmr = jm(n)
      lmr = lm(n)
      
      ! loop over tracers and proceses:
      do itr = 1, plc_ntr
        do ipr = 0, plc_npr
        
          ! no time stored yet:
          plc_dat(n,itr,ipr)%t = AnyDate()
      
          ! concentration arrays distributed over levels:
          allocate( plc_dat(n,itr,ipr)%rm_k(imr,jmr,lmloc) )
          ! init to zero:
          plc_dat(n,itr,ipr)%rm_k = 0.0
          
          ! temporary storage for some plc procs that are collected parallel over tracers ?
          ! only transported tracers could be parallel over tracers:
          if ( plc_itr_to_tm(itr) <= ntracet ) then
            ! only for some emissions and (wet) deposition:
            plc_dat(n,itr,ipr)%par_t = ( (ipr==plc_ipr_pemi) .and. (itr/=plc_itr_nox) ) .or. &
                                         (ipr==plc_ipr_lwdep) .or. &
                                         (ipr==plc_ipr_tcnvd)
          else
            plc_dat(n,itr,ipr)%par_t = .false.
          end if

          ! allocate storage ?
          if ( plc_dat(n,itr,ipr)%par_t ) then
            ! distributed to this processor ?
            if ( tracer_active(plc_itr_to_tm(itr)) ) then
               ! full 3d field:
              allocate( plc_dat(n,itr,ipr)%rm_t(imr,jmr,lmr) )
            else
              ! allocate at least one value to avoid problems in debug mode:
              allocate( plc_dat(n,itr,ipr)%rm_t(1,1,1) )
            end if
          else
            ! not used ...
            nullify( plc_dat(n,itr,ipr)%rm_t )
          end if
          
        end do
      end do
      
    end do  ! regions
    
    !
    ! done
    !
    
    ! ok
    status = 0
    
  end subroutine PLC_Init
  
  
  ! ***


  subroutine PLC_Done( status )
  
    use dims    , only : nregions
    
    ! --- in/out -------------------------------
    
    integer, intent(out)           ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/PLC_Done'
    
    ! --- local --------------------------------
    
    integer     ::  n
    integer     ::  itr, ipr
    
    ! --- begin -------------------------------
    
    ! loop over regions:
    do n = 1, nregions
    
      ! loop over tracers and proceses:
      do itr = 1, plc_ntr
        do ipr = 0, plc_npr
      
          ! clear:
          deallocate( plc_dat(n,itr,ipr)%rm_k )
          if ( associated(plc_dat(n,itr,plc_ipr_pemi)%rm_t) ) then
            deallocate( plc_dat(n,itr,plc_ipr_pemi)%rm_t )
          end if

        end do
      end do
      
    end do  ! regions
    
    ! ok
    status = 0
    
  end subroutine PLC_Done
  
  
  ! ***
  
  
  !
  ! call PLC_Set( 'action', itr, ipr, status [,ipr2=] [,fac=] )
  !
  ! action         description
  ! ----------     ---------------------------------------------------------
  !
  ! 'fill-tracer'  copy current concentrations into specified plc process
  !
  ! 'add-change'   add difference between current tracer mass and 
  !                ipr2 fields to specified process
  !
  ! 'add'          add ipr2 field to ipr
  !
  ! 'add-t'        distribute rm_t field over levels and add
  !
  !
  ! ipr
  ! -------------  ----------------------------------------------
  ! 0              temporary storage for initial concentrations
  ! plc_ipr_conc   instant concentrations
  ! plc_ipr_pchem  chemical production
  ! plc_ipr_lchem  chemical loss; now used for total tendency
  ! plc_ipr_pemi   emission production
  ! plc_ipr_ldep   deposition loss
  !
  
  
  subroutine PLC_Set( action, itr, ipr, status, ipr2, fac, t )
  
    use GO          , only : TDate, rTotal, NewDate, AnyDate, IsAnyDate, operator(-), operator(>), wrtgol
    use dims        , only : im, jm, lm
    use dims        , only : idate
    use chem_param  , only : iO3, iCO, iNO2, iNOx, iCH2O, iSO2
    use chem_param  , only : xmNOx, xmNO2
    use partools    , only : lmloc
    use partools    , only : Par_Scatter_Over_Levels, tracer_id
    use tm5_tendency, only : plc_ntr, plc_trname
    use tm5_tendency, only : plc_ipr_conc, plc_ipr_pchem, plc_ipr_lchem
    use tm5_tendency, only : plc_npr, plc_prname
  
    ! --- in/out ------------------------------------
    
    character(len=*), intent(in)    ::  action
    integer, intent(in)             ::  itr
    integer, intent(in)             ::  ipr
    integer, intent(out)            ::  status
    
    integer, intent(in), optional   ::  ipr2
    real, intent(in), optional      ::  fac
    type(TDate), optional           ::  t
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/PLC_Set'
    
    ! --- local --------------------------------
    
    integer             ::  n
    integer             ::  imr, jmr, lmr
    real, allocatable   ::  rm_k(:,:,:)
    character(len=32)   ::  rm_k_unit
    type(TDate)         ::  t_curr
    real                ::  afac
    
    ! --- begin -------------------------------
    
    ! current time:
    t_curr = NewDate(time6=idate)
    if ( present(t) )  t_curr = t

    ! loop over regions:
    do n = 1, nregions

      ! local dims:
      imr = im(n)
      jmr = jm(n)
      lmr = lm(n)

      ! what to do ?
      select case ( action )

        !
        ! concentration
        !
        case ( 'fill-tracer' )

          ! temporary storage:
          allocate( rm_k(imr,jmr,lmloc) )

          ! copy tracer fields in kg into rm_k:
          rm_k_unit = 'kg'
          call Tracer_Get_Slab( n, plc_itr_to_tm(itr), rm_k_unit, rm_k, status )
          IF_NOTOK_RETURN(status=1)

          ! apply unit conversion from (kg tm tracer) to (kg plc tracer)
          rm_k = rm_k * plc_kg_from_tm(itr)

          ! fill slab:
          plc_dat(n,itr,ipr)%rm_k = rm_k
          ! store unit:
          plc_dat(n,itr,ipr)%unit = trim(rm_k_unit)
          ! store time:
          plc_dat(n,itr,ipr)%t    = t_curr

          ! clear
          deallocate( rm_k )

        !
        ! add difference
        !
        case ( 'add-change' )

          ! requires second process ...
          if ( .not. present(ipr2) ) then
            write (gol,'("action `add-change` requires argument `ipr2` ...")'); call goErr
            TRACEBACK; status=1; return
          end if

          ! only if something has been stored ...
          if ( IsAnyDate(plc_dat(n,itr,ipr2)%t) ) then
            write (gol,'("action `add-change` requires valid time:")'); call goErr
            write (gol,'("  plc tracer name   : ",a)') trim(plc_trname(itr)); call goErr
            write (gol,'("  plc proces name   : ",a)') trim(plc_prname(ipr)); call goErr
            TRACEBACK; status=1; return
          end if

          ! add change if non-zero time step was taken since last storage:
          if ( t_curr > plc_dat(n,itr,ipr2)%t ) then

            ! temporary storage:
            allocate( rm_k(imr,jmr,lmloc) )

            ! copy tracer fields in kg into rm_k:
            rm_k_unit = 'kg'
            call Tracer_Get_Slab( n, plc_itr_to_tm(itr), rm_k_unit, rm_k, status )
            IF_NOTOK_RETURN(status=1)

            ! apply unit conversion from (kg tm tracer) to (kg plc tracer)
            rm_k = rm_k * plc_kg_from_tm(itr)

            ! add differnce between current and previous concentration:
            plc_dat(n,itr,ipr)%rm_k = plc_dat(n,itr,ipr)%rm_k + &
                                               rm_k - plc_dat(n,itr,ipr2)%rm_k

            ! store unit:
            plc_dat(n,itr,ipr)%unit = trim(rm_k_unit)

            ! clear
            deallocate( rm_k )

          else

            write (gol,'("WARNING - zero delay, thus no change in ",a," ",a," = rm - ",a)') &
                           trim(plc_trname(itr)), trim(plc_prname(ipr)), trim(plc_prname(ipr2)); call goPr
            call wrtgol( 'WARNING -   prev time : ', plc_dat(n,itr,ipr2)%t ); call goPr
            call wrtgol( 'WARNING -   curr time : ', t_curr ); call goPr

          end if


        !
        ! add two tendencies
        !
        case ( 'add' )

          ! requires second process ...
          if ( .not. present(ipr2) ) then
            write (gol,'("action `add` requires argument `ipr2` ...")'); call goErr
            TRACEBACK; status=1; return
          end if
          
          ! weight:
          afac = 1.0
          if ( present(fac) ) afac = fac

          ! add fields:
          plc_dat(n,itr,ipr)%rm_k = plc_dat(n,itr,ipr)%rm_k + &
                                             plc_dat(n,itr,ipr2)%rm_k * afac

        !
        ! add budgets collected parallel over tracers:
        !
        case ( 'add-t' )

          ! only if rm_t is setup ...
          if ( .not. plc_dat(n,itr,ipr)%par_t ) cycle

          ! temporary slab of levels:
          allocate( rm_k(imr,jmr,lmloc) )

          ! distribute over levels:
          call Par_Scatter_Over_Levels( plc_dat(n,itr,ipr)%rm_t, tracer_id(plc_itr_to_tm(itr)), rm_k, status )
          IF_NOTOK_RETURN(status=1)

          ! add to correct tendency:
          plc_dat(n,itr,ipr)%rm_k = plc_dat(n,itr,ipr)%rm_k + rm_k

          ! clear
          deallocate( rm_k )
      
        ! 
        ! unsupported ...
        !
        case default
          write (gol,'("unsupported action : ",a)') trim(action); call goErr
          TRACEBACK; status=1; return

      end select

    end do  ! regions
    
    ! ok
    status = 0
    
  end subroutine PLC_Set
  

  ! ***
  
  
  subroutine PLC_Reset( itr, ipr, status )
  
    use GO  , only : TDate, NewDate
    use dims, only : idate
  
    ! --- in/out ------------------------------------
    
    integer, intent(in)             ::  itr
    integer, intent(in)             ::  ipr
    integer, intent(out)            ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/PLC_Reset'
    
    ! --- local --------------------------------
    
    integer             ::  n
    integer             ::  itr
    type(TDate)         ::  t_curr
    
    ! --- begin -------------------------------
    
    ! current time:
    t_curr = NewDate(time6=idate)

    ! loop over regions:
    do n = 1, nregions

      ! zero concentrations:
      plc_dat(n,itr,ipr)%rm_k = 0.0
      if ( plc_dat(n,itr,ipr)%par_t ) plc_dat(n,itr,ipr)%rm_t = 0.0

      ! dummy times:
      plc_dat(n,itr,ipr)%t    = t_curr

    end do  ! regions
    
    ! ok
    status = 0
    
  end subroutine PLC_Reset
  

  ! ***
  
  
  !
  ! Add value for single cell to rm_k array
  !
  
  subroutine PLC_AddValue( region, ipr, i, j, k_loc, itrace, x, status )
  
    ! --- in/out -----------------------------------
    
    integer, intent(in)     ::  region
    integer, intent(in)     ::  ipr
    integer, intent(in)     ::  i, j, k_loc, itrace
    real, intent(in)        ::  x
    integer, intent(out)    ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/PLC_AddValue'
    
    ! --- local --------------------------------
    
    integer         ::  itr
    
    ! --- begin --------------------------------
    
    ! plc tracer index:
    itr = plc_itr_from_tm(itrace)
    
    ! add value ?
    if ( itr > 0 ) then
      plc_dat(region,itr,ipr)%rm_k(i,j,k_loc) = plc_dat(region,itr,ipr)%rm_k(i,j,k_loc) + x
    end if

    ! ok
    status = 0
    
  end subroutine PLC_AddValue
  
  
  ! ***
  
  
  !
  ! Add value for single cell to rm_k array
  !
  
  subroutine PLC_AddValue_t( region, ipr, i, j, k, itrace, x, status )
  
    ! --- in/out -----------------------------------
    
    integer, intent(in)     ::  region
    integer, intent(in)     ::  ipr
    integer, intent(in)     ::  i, j, k, itrace
    real, intent(in)        ::  x
    integer, intent(out)    ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/PLC_AddValue_t'
    
    ! --- local --------------------------------
    
    integer         ::  itr
    
    ! --- begin --------------------------------
    
    ! plc tracer index:
    itr = plc_itr_from_tm(itrace)
    
    ! not a plc tracer ?
    if ( itr <= 0 ) then
      status=0; return
    end if

    ! add value if rm_t is setup:
    if ( plc_dat(region,itr,ipr)%par_t ) then
      plc_dat(region,itr,ipr)%rm_t(i,j,k) = plc_dat(region,itr,ipr)%rm_t(i,j,k) + x
    end if

    ! ok
    status = 0
    
  end subroutine PLC_AddValue_t
  
  

#endif


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


  subroutine assign_massdat( reg, iaction )

    use dims      , only : nregions
    use dims      , only : im, jm, lm
    use chem_param, only : ntracet
    use partools  , only : ntracetloc, lmloc
    
    ! --- in/out -------------------------

    integer, intent(in)             ::  reg
    integer, intent(in)             ::  iaction
    
    ! --- local -------------------------

    integer            ::   imr, jmr, lmr, nt
    
    ! --- begin -------------------------

    select case ( iaction )
    
      !WP! swapping from tracer to k
      
      case ( 0 )
      
#ifdef MPI      

       imr = im(reg)
       jmr = jm(reg)
       lmr = lmloc
       nt  = ntracet
       allocate ( mass_dat(reg)%rm_k (-1:imr+2,-1:jmr+2,lmr,nt) )
#ifdef slopes
       allocate ( mass_dat(reg)%rxm_k(-1:imr+2,-1:jmr+2,lmr,nt) )
       allocate ( mass_dat(reg)%rym_k(-1:imr+2,-1:jmr+2,lmr,nt) )
       allocate ( mass_dat(reg)%rzm_k(-1:imr+2,-1:jmr+2,lmr,nt) )
#ifdef secmom
       allocate ( mass_dat(reg)%rxxm_k(-1:imr+2,-1:jmr+2,lmr,nt) )
       allocate ( mass_dat(reg)%rxym_k(-1:imr+2,-1:jmr+2,lmr,nt) )
       allocate ( mass_dat(reg)%rxzm_k(-1:imr+2,-1:jmr+2,lmr,nt) )
       allocate ( mass_dat(reg)%ryym_k(-1:imr+2,-1:jmr+2,lmr,nt) )
       allocate ( mass_dat(reg)%ryzm_k(-1:imr+2,-1:jmr+2,lmr,nt) )
       allocate ( mass_dat(reg)%rzzm_k(-1:imr+2,-1:jmr+2,lmr,nt) )
#endif
#endif
       mass_dat(reg)%rm_k =0.0
#ifdef slopes
       mass_dat(reg)%rxm_k=0.0
       mass_dat(reg)%rym_k=0.0
       mass_dat(reg)%rzm_k=0.0
#ifdef secmom
       mass_dat(reg)%rxxm_k=0.0
       mass_dat(reg)%rxym_k=0.0
       mass_dat(reg)%rxzm_k=0.0
       mass_dat(reg)%ryym_k=0.0
       mass_dat(reg)%ryzm_k=0.0
       mass_dat(reg)%rzzm_k=0.0
#endif
#endif

#endif

    
    ! swapping from k to tracer

    case ( 1 )
    
       imr = im(reg)
       jmr = jm(reg)
       lmr = lm(reg) 
       nt  = ntracetloc

       allocate ( mass_dat(reg)%rm_t (-1:imr+2,-1:jmr+2,lmr,nt) )
#ifdef slopes
       allocate ( mass_dat(reg)%rxm_t(-1:imr+2,-1:jmr+2,lmr,nt) )
       allocate ( mass_dat(reg)%rym_t(-1:imr+2,-1:jmr+2,lmr,nt) )
       allocate ( mass_dat(reg)%rzm_t(-1:imr+2,-1:jmr+2,lmr,nt) )
#ifdef secmom
       allocate ( mass_dat(reg)%rxxm_t(-1:imr+2,-1:jmr+2,lmr,nt) )
       allocate ( mass_dat(reg)%rxym_t(-1:imr+2,-1:jmr+2,lmr,nt) )
       allocate ( mass_dat(reg)%rxzm_t(-1:imr+2,-1:jmr+2,lmr,nt) )
       allocate ( mass_dat(reg)%ryym_t(-1:imr+2,-1:jmr+2,lmr,nt) )
       allocate ( mass_dat(reg)%ryzm_t(-1:imr+2,-1:jmr+2,lmr,nt) )
       allocate ( mass_dat(reg)%rzzm_t(-1:imr+2,-1:jmr+2,lmr,nt) )
#endif
#endif
       mass_dat(reg)%rm_t=0.0
#ifdef slopes
       mass_dat(reg)%rxm_t=0.0
       mass_dat(reg)%rym_t=0.0
       mass_dat(reg)%rzm_t=0.0
#ifdef secmom
       mass_dat(reg)%rxxm_t=0.0
       mass_dat(reg)%rxym_t=0.0
       mass_dat(reg)%rxzm_t=0.0
       mass_dat(reg)%ryym_t=0.0
       mass_dat(reg)%ryzm_t=0.0
       mass_dat(reg)%rzzm_t=0.0
#endif
#endif
  
    end select

  end subroutine assign_massdat


  ! ***
  

  subroutine free_massdat( reg, iaction )

    use dims      , only : nregions

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

    integer, intent(in)             ::  reg
    integer, intent(in)             ::  iaction
    
    ! --- local -------------------------

    integer            ::   imr, jmr, lmr, nt
    
    ! --- begin -------------------------

    select case ( iaction )
    
      ! swapping from tracer to k
      case ( 0 )
      
        deallocate ( mass_dat(reg)%rm_t)
#ifdef slopes
        deallocate ( mass_dat(reg)%rxm_t)
        deallocate ( mass_dat(reg)%rym_t)
        deallocate ( mass_dat(reg)%rzm_t)
#ifdef secmom
        deallocate ( mass_dat(reg)%rxxm_t)
        deallocate ( mass_dat(reg)%rxym_t)
        deallocate ( mass_dat(reg)%rxzm_t)
        deallocate ( mass_dat(reg)%ryym_t)
        deallocate ( mass_dat(reg)%ryzm_t)
        deallocate ( mass_dat(reg)%rzzm_t)
#endif
#endif

      !WP! swapping from k to tracer 
      case ( 1 )
    
#ifdef MPI      

        deallocate ( mass_dat(reg)%rm_k)
#ifdef slopes
        deallocate ( mass_dat(reg)%rxm_k)
        deallocate ( mass_dat(reg)%rym_k)
        deallocate ( mass_dat(reg)%rzm_k)
#ifdef secmom
        deallocate ( mass_dat(reg)%rxxm_k)
        deallocate ( mass_dat(reg)%rxym_k)
        deallocate ( mass_dat(reg)%rxzm_k)
        deallocate ( mass_dat(reg)%ryym_k)
        deallocate ( mass_dat(reg)%ryzm_k)
        deallocate ( mass_dat(reg)%rzzm_k)
#endif
#endif

#endif

    end select

  end subroutine free_massdat
  
  
  ! ===================================================================


  subroutine Par_Check_Domain( region1, others, to_check )

    use dims     , only : parent,children, okdebug
    use ParTools , only : myid, root
    use ParTools , only : previous_par
    use ParTools , only : allocate_mass
#ifdef MPI
    use mpi_comm , only : barrier
    use mpi_comm , only : escape_tm_mpi
#endif

    ! --- in/out -----------------------------------
    
    integer,intent(in)              ::  region1 
    character(len=*),intent(in)     ::  others
    character(len=*),intent(in)     ::  to_check

    ! --- local ------------------------------------
    
#ifdef MPI   
    integer :: reg,regstart,regstop,iaction
    integer :: nchild,regiop
#endif

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

#ifdef MPI   

    select case(others)
      case('p') !WP! also check parents
        regiop=parent(region1)
        if ( regiop == 0 ) then 
          regstart=region1
          regstop=region1
        else if ( regiop /= 0 ) then
          regstart=regiop
          regstop=region1
        end if
      case('c') !WP! also check child  
        nchild=children(region1,0)
        if ( nchild == 0 ) then
          regstart=region1
          regstop=region1
        else if ( nchild /= 0 ) then
          regstart=region1
          !WP! take number of last child as regstop
          regstop=children(region1,nchild)
        end if
      case('n') !WP! check no parents or children
        regstart=region1
        regstop=region1
      case default
        call escape_tm_mpi('check_domain: illegal value of *others* ')
    end select

    select case(to_check)
      case('levels')
        iaction=0
      case('tracer') 
        iaction=1
      case default
        call escape_tm_mpi('check_domain: illegal value of *to_check* ')
    end select

    do reg=regstart,regstop
      ! make sure data is in requested domain
      if ( previous_par(reg) /= to_check ) then
        if ( okdebug ) then
          write (gol,'("check_domain: data will be swapped")'); call goPr
        end if
        previous_par(reg)=to_check   ! reset previous par
        ! assign new memory
        if ( allocate_mass ) call assign_massdat(reg,iaction)
        call barrier
        ! swap data for this region
        call swap_all_mass(reg,iaction)
        call barrier
        ! free previously allocated memory
        if ( allocate_mass ) call free_massdat(reg,iaction)
      else
        if ( okdebug ) then
          write (gol,'("check_domain: data is on proper communicator, nothing to do")'); call goPr
        end if
      end if
    end do

#endif

  end subroutine Par_Check_Domain


  ! ***
  

#ifdef MPI      

  !WP! subroutine swaps the arrays m,rm,rxm,rym,and rzm 
  !WP! from levels to tracers (iaction=1)
  !WP! or from levels to tracers (iaction=0)
  !WP! for requested regions (region1=>region2)    
  !WP! after swapping, old data is set to zero to avoid mistakes 
  !WP! with 'old' data
  !WP! 23 december 2002

  subroutine swap_all_mass(region,iaction)

    use dims,       only : okdebug, im,jm ,lm
    use chem_param, only : ntracet
    use mpi_const , only : myid, root
    use mpi_const , only : ntracetloc, ntracet_ar
    use mpi_comm  , only : barrier
    use mpi_comm  , only : scatter_after_read_k_3d, gather_tracer_k_3d
    use mpi_comm  , only : tracer_to_k, k_to_tracer

    ! --- in/out ----------------------------
    
    integer,intent(in) :: iaction
    integer,intent(in) :: region
    
    ! --- local -----------------------------

    integer :: imr,jmr,lmr
    real,dimension(:,:,:,:),pointer     :: rxmsource
    real,dimension(:,:,:,:),pointer     :: rxmtarget
    integer :: nsend,xx, xx_end

    ! --- begin -----------------------------
    
    imr=im(region)
    jmr=jm(region)
    lmr=lm(region)

#ifdef slopes
#ifndef secmom
    xx_end = 4
#else
    xx_end = 10
#endif
#endif
        
    !WP! once for each x,y, and z
    do xx = 1, xx_end   
      select case(xx)
        case( 1 )
          rxmsource => mass_dat(region)%rm_t
          rxmtarget => mass_dat(region)%rm_k
#ifdef slopes
        case( 2 )
          rxmsource => mass_dat(region)%rxm_t
          rxmtarget => mass_dat(region)%rxm_k
        case( 3 )
          rxmsource => mass_dat(region)%rym_t
          rxmtarget => mass_dat(region)%rym_k
        case( 4 )
          rxmsource => mass_dat(region)%rzm_t
          rxmtarget => mass_dat(region)%rzm_k
#ifdef secmom
        case( 5 )
          rxmsource => mass_dat(region)%rxxm_t
          rxmtarget => mass_dat(region)%rxxm_k
        case( 6 )
          rxmsource => mass_dat(region)%rxym_t
          rxmtarget => mass_dat(region)%rxym_k
        case( 7 )
          rxmsource => mass_dat(region)%ryym_t
          rxmtarget => mass_dat(region)%ryym_k
        case( 8 )
          rxmsource => mass_dat(region)%rxzm_t
          rxmtarget => mass_dat(region)%rxzm_k
        case( 9 )
          rxmsource => mass_dat(region)%ryzm_t
          rxmtarget => mass_dat(region)%ryzm_k
        case( 10 )
          rxmsource => mass_dat(region)%rzzm_t
          rxmtarget => mass_dat(region)%rzzm_k
#endif
#endif
      end select
      if ( iaction == 0 ) then
        call tracer_to_k( rxmtarget, rxmsource, &
                            imr,jmr,lmr, 2,2,0, ntracet,ntracetloc, &
                            ntracet_ar )
        call barrier
      else if ( iaction == 1 ) then
        call k_to_tracer( rxmtarget, rxmsource, &
                            imr,jmr,lmr, 2,2,0, ntracet, ntracetloc, &
                            ntracet_ar)
      end if
      nullify(rxmsource)               
      nullify(rxmtarget)

      if ( okdebug ) then
        write (gol,'("swap_all_mass: Scattered tracer mass and slopes ",i4)') xx; call goPr
      end if

    end do

  end subroutine swap_all_mass

#endif


  ! ================================================================
  
  ! PLS: Brought back to life to print min, max and total sum of air and tracer
  ! masses. Used for comparison with TM6.
  
  subroutine Par_Check_Mass( region, text, tr_id )


    use dims      , only : lm, im, jm
#ifdef MPI
    use chem_param, only : ntracet
    use mpi_const , only : myid, root_t, root_k
    use mpi_const , only : ntracetloc, lmloc
    use mpi_const,   only : com_trac
    use mpi_const,   only : ierr,my_real,mpi_max,mpi_min, mpi_sum, ntracet_ar
    use mpi_comm,    only : barrier_t
    use mpi_comm  , only : dump_field3d, dump_field4d, gather_tracer_t
!    use mpi_comm  , only : summasr
#endif
    use meteodata, only : m_dat
    use partools, only : which_par, previous_par
    use partools, only : tracer_active, tracer_loc
    
    ! --- in/out -----------------------------------
    
    integer, intent(in)             ::  region
    character(len=*), intent(in)    ::  text
    integer, intent(in), optional   :: tr_id   ! default to 1

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

    real,dimension(:,:,:),pointer     :: m
    real,dimension(:,:,:,:),pointer   :: rm, rmt

    integer :: imr, jmr, lmr, trid, ijmax3(3), ijmin3(3), ijmax4(4), ijmin4(4)
    real    :: min_one, max_one, min_all, max_all, tot_one, tot_all
    ! --- begin ------------------------------------

    imr=im(region);jmr=jm(region);lmr=lm(region)

    trid=1
    if(present(tr_id)) trid=tr_id
    
#ifdef MPI
    which_par=previous_par(region)

!    if(which_par == 'tracer'.and.ntracetloc == 0) return 
    if(which_par == 'levels') then !.and.lmloc == 0) return
       print*, "check_mass not implemented for decomposition across levels"
       return
    end if
       
#endif

    rm => mass_dat(region)%rm_t
    m => m_dat(region)%data
    
    max_one=maxval(m(1:imr,1:jmr,1:lmr))
    min_one=minval(m(1:imr,1:jmr,1:lmr))
    tot_one=sum(m(1:imr,1:jmr,1:lmr))
    ijmax3=maxloc(m(1:imr,1:jmr,1:lmr) )
    ijmin3=minloc(m(1:imr,1:jmr,1:lmr) )
    
#ifdef MPI   

!    call barrier_t
!    call mpi_reduce(max_one, max_all, 1, my_real, mpi_max, root_t, com_trac, ierr)
!    call mpi_reduce(min_one, min_all, 1, my_real, mpi_min, root_t, com_trac, ierr)
!    call mpi_reduce(tot_one, tot_all, 1, my_real, mpi_max, root_t, com_trac, ierr) !same on all PEs, one is enough!
    if ( myid == root_t ) then
       print *, text//' [par_check_mass] AIR maxloc =',ijmax3
       print *, text//' [par_check_mass] AIR minloc =',ijmin3
!       print *, text//' [par_check_mass] AIR min, max, sum: ', min_all, max_all, tot_all
       print *, text//' [par_check_mass] AIR min, max, sum: ', min_one, max_one, tot_one
    end if

    allocate(rmt(-1:imr+2,-1:jmr+2,1:lmr,ntracet))
    call gather_tracer_t( rmt, im(region),jm(region),lm(region), 2,2,0, ntracet, rm, .false. ) 
    if ( myid == root_t ) then
       max_all=maxval( rmt(1:imr,1:jmr,1:lmr,:) )
       min_all=minval( rmt(1:imr,1:jmr,1:lmr,:) )
       tot_all=sum(    rmt(1:imr,1:jmr,1:lmr,:) )
       ijmax4=maxloc(  rmt(1:imr,1:jmr,1:lmr,:) )
       ijmin4=minloc(  rmt(1:imr,1:jmr,1:lmr,:) )
    end if
    deallocate(rmt)
    
    if ( myid == root_t ) then
       print *, text//' [par_check_mass] RM maxloc =',ijmax4
       print *, text//' [par_check_mass] RM minloc =',ijmin4
       print *, text//' [par_check_mass] RM min, max, sum: ', min_all, max_all, tot_all
    end if
#else
    print *,text//' [par_check_mass] AIR min, max, sum: ',min_one, max_one, tot_one
    max_one=maxval(rm(1:imr,1:jmr,1:lmr,:))
    min_one=minval(rm(1:imr,1:jmr,1:lmr,:))
    tot_one=sum(rm(1:imr,1:jmr,1:lmr,:))
    print *,text//' [par_check_mass] RM min, max, sum: ',min_one, max_one, tot_one
#endif

    ! for one tracer only
    if (tracer_active(trid)) then
       max_one=maxval( rm(1:imr,1:jmr,1:lmr,tracer_loc(trid)) )
       min_one=minval( rm(1:imr,1:jmr,1:lmr,tracer_loc(trid)) )
       tot_one=sum(    rm(1:imr,1:jmr,1:lmr,tracer_loc(trid)) )  ! THIS IS OK
       ijmax3=maxloc(rm(1:imr,1:jmr,1:lmr,tracer_loc(trid)) )
       ijmin3=minloc(rm(1:imr,1:jmr,1:lmr,tracer_loc(trid)) )
       print *, text//' [par_check_mass] ID maxloc =',ijmax3
       print *, text//' [par_check_mass] ID minloc =',ijmin3
       print *, text//' [par_check_mass] ID, min, max, sum: ', trid, min_one, max_one, tot_one
    end if
    
    ! PREVIOUSLY, but was commented anyway:
    
!#ifdef MPI
!
!    call barrier   
!    if ( previous_par(region) == 'tracer' ) then
!       if ( myid == root_t ) print*,'check_mass: checking mass in ',text
!       call barrier
!       m => mass_dat(region)%m_t
!       rm => mass_dat(region)%rm_t
!       lmr=lm(region)
!
!       call dump_field3d(region,m,lmr,(/2,2,0/),'m','mtrac.hdf')
!       call dump_field4d(region,rm,lmr,ntracetloc,(/2,2,0/),'rm','rmtrac.hdf')
!       call summasr( region, m, rm )
!    else if ( previous_par(region) == 'levels' ) then   
!       if ( myid == root_k ) print*,'check_mass: checking mass in ',text
!       call barrier
!       m => mass_dat(region)%m_k
!       rm => mass_dat(region)%rm_k
!       lmr=lmloc
!
!       call dump_field3d(region,m,lmr,(/2,2,0/),'m','mlev.hdf')
!       call dump_field4d(region,rm,lmr,ntracet,(/2,2,0/),'rm','rmlev.hdf')
!       call summasr( region, m, rm )
!    end if
!    call barrier
!    nullify(m)
!    nullify(rm)
!
!#endif
!
  end subroutine Par_Check_Mass



  ! ================================================================
  
    
  !
  ! Tracer slopes (kg/half cell) and second moments (kg/(half cell)^2)
  ! are adjusted with the same factor as applied to the total mass
  ! (background: if the concentration 'in the middle' is decreased
  ! by a factor 0.5, then the concentrations at the boundaries are
  ! also changed with a factor 0.5 :
  !
  !       ^      |       |             |       |
  ! tracer|    6 |       o             |       |
  !  mass |      |     - |             |       |
  !  (kg)      4 |   o   |     ->      |       |    
  !              | -     |           3 |     _ o
  !            2 o       |           2 | _ o   |
  !              |       |           1 o       |
  !           ---------------       ---------------
  !
  ! NOTE: not implemented for parallel yet ...
  !  ? should l be a local or a global level ?
  !  ? assume that tracer k is on this pe ?
  !

  subroutine AdjustTracer( drm, region, i, j, l, k, status )

    use GO      , only : gol, goErr 
    use ParTools, only : which_par, tracer_active, lmloc, offsetn

    ! --- in/out -----------------------------  
    
    real, intent(in)           ::  drm           ! tracer mass change (kg)
    integer, intent(in)        ::  region
    integer, intent(in)        ::  i, j           ! cell indices
    integer, intent(in)        ::  l              ! local level
    integer, intent(in)        ::  k             ! tracer index
    integer, intent(out)       ::  status
    
    ! --- const -------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/AdjustTracer'
    
    ! --- local -------------------------------
    
    real, dimension(:,:,:,:), pointer   :: rm
#ifdef slopes
    real, dimension(:,:,:,:), pointer   :: rxm, rym, rzm
#ifdef secmom
    real, dimension(:,:,:,:), pointer   :: rxxm, rxym, rxzm, ryym, ryzm, rzzm
#endif
#endif

    real   ::  rm_old
    real   ::  skf

    ! --- begin -----------------------------
    
    ! set pointers
    select case ( which_par )
      case ( 'levels' )
        ! check local level index:
        if ( l > lmloc ) then
          write (gol,'("found strange local level:")'); call goErr
          write (gol,'("  l, lmloc  :  ", 2i4)') l, lmloc; call goErr
          TRACEBACK; status=1; return
        end if
        ! set pointers:
        rm   => mass_dat(region)%rm_k
#ifdef slopes
        rxm  => mass_dat(region)%rxm_k
        rym  => mass_dat(region)%rym_k
        rzm  => mass_dat(region)%rzm_k
#ifdef secmom
        rxxm => mass_dat(region)%rxxm_k
        rxym => mass_dat(region)%rxym_k
        rxzm => mass_dat(region)%rxzm_k
        ryym => mass_dat(region)%ryym_k
        ryzm => mass_dat(region)%ryzm_k
        rzzm => mass_dat(region)%rzzm_k
#endif
#endif
      case ( 'tracer' )
        ! leave if this pe does not contain this tracer :
        if ( .not. tracer_active(k+offsetn) ) then
          status=0; return
        end if
        ! set pointers:
        rm   => mass_dat(region)%rm_t
#ifdef slopes
        rxm  => mass_dat(region)%rxm_t
        rym  => mass_dat(region)%rym_t
        rzm  => mass_dat(region)%rzm_t
#ifdef secmom
        rxxm => mass_dat(region)%rxxm_t
        rxym => mass_dat(region)%rxym_t
        rxzm => mass_dat(region)%rxzm_t
        ryym => mass_dat(region)%ryym_t
        ryzm => mass_dat(region)%ryzm_t
        rzzm => mass_dat(region)%rzzm_t
#endif
#endif
      case default
        write (gol,'("unsupported which_par : ",a)') trim(which_par); call goErr
        TRACEBACK; status=1; return
    end select

    ! store old tracer mass:
    rm_old = rm(i,j,l,k)
    
    ! store new tracer mass:
    rm(i,j,l,k) = rm_old + drm
#ifdef slopes
    ! adjust concentration gradients

    ! scale factor to be applied; trap devision by zero:
    if ( (rm(i,j,l,k) > tiny(1.0)) .and. (rm_old > tiny(0.0)) ) then
      skf = rm(i,j,l,k) / rm_old 
    else
      skf = 0.0
    end if

    ! adjust slopes:
    rxm(i,j,l,k) = rxm(i,j,l,k) * skf
    rym(i,j,l,k) = rym(i,j,l,k) * skf
    rzm(i,j,l,k) = rzm(i,j,l,k) * skf

#ifdef secmom
    ! adjust second moments:
    rxxm(i,j,l,k) = rxxm(i,j,l,k) * skf
    rxym(i,j,l,k) = rxym(i,j,l,k) * skf
    rxzm(i,j,l,k) = rxzm(i,j,l,k) * skf
    ryym(i,j,l,k) = ryym(i,j,l,k) * skf
    ryzm(i,j,l,k) = ryzm(i,j,l,k) * skf
    rzzm(i,j,l,k) = rzzm(i,j,l,k) * skf
#endif
#endif
    
    ! ok
    status = 0
    
  end subroutine AdjustTracer
  

  ! ***

  
  subroutine tracer_print( key, region, status )
  
    use GO
    use dims, only : im, jm, lm
    use chem_param, only : ntrace, ntracet, names
    !use chem_param, only : io3, ico, ino
    use partools, only : myid, root, npes
    use partools, only : lmar, lmloc, offsetl
    use partools, only : tracer_active, tracer_loc
    use partools, only : previous_par, which_par
    use partools, only : par_barrier
    use meteodata, only : m_dat, phlb_dat
  
    ! --- in/out ---------------------------------
    
    character(len=*), intent(in)  ::  key
    integer, intent(in)           ::  region
    integer, intent(out)          ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/tracer_print'
    
    ! --- local ---------------------------------
    
    integer             ::  i, j, ii, jj, ll, kk
    integer             ::  l, l_loc
    integer             ::  imr, jmr, lmr
    real, allocatable   ::  pat(:,:)
#ifdef slopes
    real, allocatable   ::  patx(:,:)
    real, allocatable   ::  paty(:,:)
    real, allocatable   ::  patz(:,:)
#endif
    real, allocatable   ::  patm(:,:)
    real, allocatable   ::  patp(:,:)
    integer             ::  itr
    logical             ::  filled

    character(len=20)   ::  pkey
    
    ! --- begin ----------------------------------

    ! set those for the box/tracer used to debug
    ii=45
    jj=41
    ll=33
    kk=18

    
    !if (region/=3) return
    
    imr = im(region)
    jmr = jm(region)
    lmr = lm(region)
    
    allocate( pat(imr,jmr) )
#ifdef slopes
    allocate( patx(imr,jmr) )
    allocate( paty(imr,jmr) )
    allocate( patz(imr,jmr) )
#endif
    allocate( patm(imr,jmr) )
    allocate( patp(imr,jmr) )

    !if (tracer_active(1)) then    
    !write (gol,'(a," region, previous_par : ",i1," ",a)') key, region, previous_par(region); call goBug
    !endif
!    call par_barrier

    ! info ...
    !write (gol,*) key, region; call goBug

    do itr = 1, ntrace

      if (itr/=kk) cycle      ! only requested tracer

      do l = 1, lmr   
      
        if (l/=ll) cycle      ! only requested level

        filled = .false.
        if ( itr <= ntracet ) then
          select case ( previous_par(region) )
            case ( 'levels' )
              l_loc = l - offsetl
              if ( (l_loc >= 1) .and. (l_loc <= lmloc) ) then
                write (pkey,'("[rm_k",2i3,"]")') l_loc, itr
                pat  = mass_dat(region)%rm_k (1:imr,1:jmr,l_loc,itr)
#ifdef slopes
                patx = mass_dat(region)%rxm_k(1:imr,1:jmr,l_loc,itr)
                paty = mass_dat(region)%rym_k(1:imr,1:jmr,l_loc,itr)
                patz = mass_dat(region)%rzm_k(1:imr,1:jmr,l_loc,itr)
#endif
                filled = .true.
              end if
            case ( 'tracer' )
              if ( tracer_active(itr) ) then
                write (pkey,'("[rm_t",2i3,"]")') l, tracer_loc(itr)
                pat  = mass_dat(region)%rm_t (1:imr,1:jmr,l,tracer_loc(itr))
#ifdef slopes
                patx = mass_dat(region)%rxm_t(1:imr,1:jmr,l,tracer_loc(itr))
                paty = mass_dat(region)%rym_t(1:imr,1:jmr,l,tracer_loc(itr))
                patz = mass_dat(region)%rzm_t(1:imr,1:jmr,l,tracer_loc(itr))
#endif
                filled = .true.
              end if
            case default
              write (gol,'("unsupported which_par : ",a)') trim(which_par); call goErr
              TRACEBACK; status=1; return
          end select
        else if ( itr <= ntrace ) then
          l_loc = l - offsetl
          if ( (l_loc >= 1) .and. (l_loc <= lmloc) ) then
            write (pkey,'("[rmck",2i3,"]")') l_loc, itr
            pat  = chem_dat(region)%rm_k(1:imr,1:jmr,l_loc,itr)    ! ntracet+1:ntrace
#ifdef slopes
            patx = 0.0
            paty = 0.0
            patz = 0.0
#endif
            filled = .true.
          end if
        end if
        
        patm =    m_dat(region)%data(1:imr,1:jmr,l)
        patp = phlb_dat(region)%data(1:imr,1:jmr,l)
        
        if ( filled ) then

          !write (gol,'(a,i3," ",a,i3,3es14.4)') key, itr, names(itr), l, maxval(pat), sum(pat)/size(pat), pat(3,3); call goBug

!#ifdef slopes
!          write (gol,*) key, region, itr, names(itr), l, &
!               sum(abs(patm)), sum(abs(patp)), &
!               sum(abs(pat)), sum(abs(patx)), sum(abs(paty)), sum(abs(patz)); call goBug
!#else
!          write (gol,*) key, region, itr, names(itr), l, &
!               sum(abs(patm)), sum(abs(patp)), &
!               sum(abs(pat)); call goBug
!#endif
!          write (gol,*) key, trim(pkey), region, itr, names(itr), l, sum(abs(pat)); call goBug
!
!           write (gol,*) key, region, itr, names(itr), l, pat(22+1,28+1), sum(abs(pat)); call goBug

#ifdef slopes           
           write (gol,*) key//" "//names(itr), pat(ii,jj), patx(ii,jj), paty(ii,jj), patz(ii,jj), patm(ii,jj) ; call goBug
#else
           write (gol,*) key//" "//names(itr), pat(ii,jj), patm(ii,jj) ; call goBug
#endif
!          do j = 1, jmr
!            do i = 1, imr
!              if ((region==1).and.((i/=28).or.(j/=26))) cycle
!              if ((region==2).and.((i/=25).or.(j/=8))) cycle
!              if ((region==3).and.((i>12).or.(j>1))) cycle
!              write (gol,*) key, region, itr, names(itr), i,j,l, &
!                patm(i,j), patp(i,j), pat(i,j), patx(i,j), paty(i,j), patz(i,j); call goBug
!            end do
!          end do

        end if

!          do j = 1, jmr
!            do i = 1, imr
!              if ((region==1).and.((i/=28).or.(j/=26))) cycle
!              if ((region==2).and.((i/=25).or.(j/=8))) cycle
!              !if ((region==3).and.((i>12).or.(j>1))) cycle
!              if ((region==3).and.((i/=7).or.(j/=1))) cycle
!              if (myid==root) then
!                write (gol,*) key, region, i,j,l, patm(i,j), patp(i,j); call goBug
!              endif
!              call par_barrier
!              if (myid==npes-1) then
!                write (gol,*) key, region, i,j,l, patm(i,j), patp(i,j); call goBug
!              endif
!              call par_barrier
!            end do
!          end do

!        call par_barrier
        
      end do  ! levels
  
!      call par_barrier
        
    end do  ! tracers

    ! info ...
!    write (gol,*) key, region; call goBug
    
    deallocate( pat )
#ifdef slopes
    deallocate( patx )
    deallocate( paty )
    deallocate( patz )
#endif
    deallocate( patm )
    deallocate( patp )

    ! ok
    status = 0
    
  end subroutine tracer_print


  ! ***
  
  
  !
  ! Return slab of tracer mass for levels on this processor.
  !  
#ifdef with_tendencies
  
  subroutine Tracer_Get_Slab( region, itr, unit, rm_k, status )
  
    use dims       , only : im, jm, lm
    use chem_param , only : ntrace, ntracet, fscale
    use partools   , only : myid, root
    use partools   , only : previous_par, tracer_loc, tracer_id, tracer_active
    use partools   , only : lmloc, lmar, offsetl
    use partools   , only : Par_Scatter_Over_Levels
    use meteodata  , only : m_dat
    
    ! --- in/out -------------------------------------
    
    integer, intent(in)             ::  region
    integer, intent(in)             ::  itr
    character(len=*), intent(in)    ::  unit
    real, intent(out)               ::  rm_k(:,:,:)
    integer, intent(out)            ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Tracer_Get_Slab'
    
    ! --- local ---------------------------------
    
    integer             ::  imr, jmr, lmr
    integer             ::  itr_loc
    
    ! --- begin ----------------------------------

    ! local grid size    
    imr = im(region)
    jmr = jm(region)
    lmr = lm(region)
    
    ! check input:
    if ( any(shape(rm_k)/=(/imr,jmr,lmloc/)) ) then
      write (gol,'("shape of output grid not ok:")'); call goErr
      write (gol,'("  shape(rm_k)    : ",3i4)') shape(rm_k); call goErr
      write (gol,'("  imr,jmr,lmloc  : ",3i4)') imr,jmr,lmloc; call goErr
      TRACEBACK; status=1; return
    end if
    
    ! transported or short lived tracer ?
    if ( (itr >= 1) .and. (itr <= ntracet) ) then
      
      ! how distributed ?
      select case ( previous_par(region) )

        case ( 'tracer' )
        
          ! local tracer index if requested tracer is on this pe,
          ! or just the first local index on other pe's to avoid array bound problems:
          if ( tracer_active(itr) ) then
            itr_loc = tracer_loc(itr)
          else
            itr_loc = 1
          end if

          ! distribute full 3D field stored in  mass_dat%rm_t  on processor  tracer_id(itr)
          ! over rm_k arrays available on each processor;
          ! distribution is over the last dimension (levels):
          call Par_Scatter_Over_Levels( &
                  mass_dat(region)%rm_t(1:imr,1:jmr,1:lmr,itr_loc), &
                  tracer_id(itr), rm_k, status )
          IF_NOTOK_RETURN(status=1)
          
          ! change unit ?
          select case ( unit )
            case ( 'kg' )
              ! rm already in kg
              ! apply unit conversion from (kg tm tracer) to (kg plc tracer)
              rm_k = rm_k * plc_kg_from_tm(itr)
            case ( 'mass-mixing-ratio' )
              ! mass mixing ratio
              ! apply unit conversion from (kg tm tracer) to (kg plc tracer)
              rm_k = rm_k * plc_kg_from_tm(itr)
              rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc)
            case ( 'volume-mixing-ratio' )
              ! volume mixing ratio = mass mixing ratio * fscale
              ! apply unit conversion from (kg tm tracer) to (kg plc tracer)
              rm_k = rm_k * plc_kg_from_tm(itr)
              rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) * fscale(itr)
            case default
              write (gol,'("unsupported unit for par tracer : ",a)') unit; call goErr
              TRACEBACK; status=1; return
          end select

        case ( 'levels' )
        
          ! array already distributed over levels; just copy:
          rm_k = mass_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,itr)

          ! change unit ?
          select case ( unit )
            case ( 'kg' )
              ! rm already in kg
              ! apply unit conversion from (kg tm tracer) to (kg plc tracer)
              rm_k = rm_k * plc_kg_from_tm(itr)
            case ( 'mass-mixing-ratio' )
              ! mass mixing ratio
              ! apply unit conversion from (kg tm tracer) to (kg plc tracer)
              rm_k = rm_k * plc_kg_from_tm(itr)
              rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc)
            case ( 'volume-mixing-ratio' )
              ! volume mixing ratio = mass mixing ratio * fscale
              ! apply unit conversion from (kg tm tracer) to (kg plc tracer)
              rm_k = rm_k * plc_kg_from_tm(itr)
              rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) * fscale(itr)
            case default
              write (gol,'("unsupported unit for par level : ",a)') unit; call goErr
              TRACEBACK; status=1; return
          end select

        case default

          write (gol,'("unsupported distribution:")'); call goPr
          write (gol,'("  region         : ",i2)') region ;call goPr
          write (gol,'("  previous_par   : ",a)') previous_par(region) ;call goPr
          TRACEBACK; status=1; return

      end select

    else if ( (itr > ntracet) .and. (itr <= ntrace) ) then
    
      ! extract short-lived tracer:
      rm_k = chem_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,itr)
      
      ! change unit ?
      select case ( unit )
        case ( 'kg' )
          ! rm already in kg
          ! apply unit conversion from (kg tm tracer) to (kg plc tracer)
          rm_k = rm_k * plc_kg_from_tm(itr)
        case ( 'mass-mixing-ratio' )
          ! mass mixing ratio
          ! apply unit conversion from (kg tm tracer) to (kg plc tracer)
          rm_k = rm_k * plc_kg_from_tm(itr)
          rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc)
        case ( 'volume-mixing-ratio' )
          ! volume mixing ratio = mass mixing ratio * fscale
          ! apply unit conversion from (kg tm tracer) to (kg plc tracer)
          rm_k = rm_k * plc_kg_from_tm(itr)
          rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) * fscale(itr)
        case default
          write (gol,'("unsupported unit for non-transported tracer slab : ",a)') unit; call goErr
          TRACEBACK; status=1; return
      end select

    else
    
      write (gol,'("unsupported tracer index : ",i4)') itr; call goErr
      write (gol,'("  ntrace, ntracet : ",2i4)') ntrace, ntracet; call goErr
      TRACEBACK; status=1; return

    end if  ! transported or chem-only

    ! ok
    status = 0
    
  end subroutine Tracer_Get_Slab

#endif

  ! ***
  
  
  !
  ! Fill slabs of tracer mass (each processor some levels) with given unit
  ! into tracer mass array (in kg); adust slopes too.
  !  
#ifdef with_feedback
  
  subroutine Tracer_Fill_Slabs( region, itr, unit, rm_k, status )
  
    use dims       , only : im, jm, lm
    use chem_param , only : ntrace, ntracet, fscale
    use partools   , only : myid, root
    use partools   , only : previous_par, tracer_loc, tracer_id, tracer_active
    use partools   , only : lmloc, lmar, offsetl
    use partools   , only : Par_Gather_From_Levels
    use meteodata  , only : m_dat
    
    ! --- in/out -------------------------------------
    
    integer, intent(in)             ::  region
    integer, intent(in)             ::  itr
    character(len=*), intent(in)    ::  unit
    real, intent(in)                ::  rm_k(:,:,:)
    integer, intent(out)            ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Tracer_Fill_Slabs'
    
    ! --- local ---------------------------------
    
    integer             ::  imr, jmr, lmr
    integer             ::  i, j, l
    integer             ::  itr_loc
    real, allocatable   ::  x_k(:,:,:)
    real, allocatable   ::  x_t(:,:,:)
    
    ! --- begin ----------------------------------

    ! local grid size    
    imr = im(region)
    jmr = jm(region)
    lmr = lm(region)
    
    ! check input:
    if ( any(shape(rm_k)/=(/imr,jmr,lmloc/)) ) then
      write (gol,'("shape of output grid not ok:")'); call goErr
      write (gol,'("  shape(rm_k)    : ",3i4)') shape(rm_k); call goErr
      write (gol,'("  imr,jmr,lmloc  : ",3i4)') imr,jmr,lmloc; call goErr
      TRACEBACK; status=1; return
    end if
    
    ! transported or short lived tracer ?
    if ( (itr >= 1) .and. (itr <= ntracet) ) then
      
      ! how distributed ?
      select case ( previous_par(region) )

        ! rm is parallel over tracers:
        case ( 'tracer' )
        
          ! local tracer index if requested tracer is on this pe,
          ! or just the first local index on other pe's to avoid array bound problems;
          ! allocate storage for 3d field on target pe:
          if ( tracer_active(itr) ) then
            itr_loc = tracer_loc(itr)
            allocate( x_t(imr,jmr,lmr) )
          else
            itr_loc = 1
            allocate( x_t(1,1,1) )
          end if
          
          ! collect rm_k arrays available on each processor
          ! into x_t on processor tracer_id(itr)
          call Par_Gather_From_Levels( rm_k, x_t, tracer_id(itr), status )
          IF_NOTOK_RETURN(status=1)
          
          ! now for processor with requested tracer only:
          if ( tracer_active(itr) ) then

            ! rm_k was in given unit; convert x_t to kg :
            select case ( unit )
              case ( 'kg' )
                ! x_t already in kg
                ! apply unit conversion from  (kg plc tracer) to (kg tm tracer)
                x_t = x_t / plc_kg_from_tm(itr)
              case ( 'mass-mixing-ratio' )
                ! mass mixing ratio
                ! apply unit conversion from  (kg plc tracer) to (kg tm tracer)
                x_t = x_t / plc_kg_from_tm(itr)
                x_t = x_t * m_dat(region)%data(1:imr,1:jmr,1:lmr)
              case ( 'volume-mixing-ratio' )
                ! volume mixing ratio = mass mixing ratio * fscale
                ! apply unit conversion from  (kg plc tracer) to (kg tm tracer)
                x_t = x_t / plc_kg_from_tm(itr)
                x_t = x_t * m_dat(region)%data(1:imr,1:jmr,1:lmr) / fscale(itr)
              case default
                write (gol,'("unsupported unit for par tracer : ",a)') unit; call goErr
                TRACEBACK; status=1; return
            end select

            ! replace new mass x_t by difference with rm :
            x_t = x_t - mass_dat(region)%rm_t(1:imr,1:jmr,1:lmr,itr_loc)   ! kg
            ! add change in kg to tracer arrays, adjust slopes:
            do l = 1, lmr
              do j = 1, jmr
                do i = 1, imr
                  ! adjust tracer arrays for local indices:
                  call AdjustTracer( x_t(i,j,l), region, i, j, l, itr_loc, status )
                  IF_NOTOK_RETURN(status=1)
                end do
              end do
            end do
            
          end if   ! tracer active

          ! clear
          deallocate( x_t )

        ! rm is parallel over levels:
        case ( 'levels' )
        
          ! storage for tracer change:
          allocate( x_k(imr,jmr,lmloc) )
          
          ! copy input:
          x_k = rm_k

          ! rm_k was in given unit; convert x_k to kg :
          select case ( unit )
            case ( 'kg' )
              ! x_k already in kg
              ! apply unit conversion from  (kg plc tracer) to (kg tm tracer)
              x_k = x_k / plc_kg_from_tm(itr)
            case ( 'mass-mixing-ratio' )
              ! from mass mixing ratio to tracer mass:
              x_k = x_k * m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc)
              ! apply unit conversion from  (kg plc tracer) to (kg tm tracer)
              x_k = x_k / plc_kg_from_tm(itr)
            case ( 'volume-mixing-ratio' )
              ! volume mixing ratio = mass mixing ratio * fscale
              x_k = x_k * m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) / fscale(itr)
              ! apply unit conversion from  (kg plc tracer) to (kg tm tracer)
              x_k = x_k / plc_kg_from_tm(itr)
            case default
              write (gol,'("unsupported unit for par level : ",a)') unit; call goErr
              TRACEBACK; status=1; return
          end select

          ! replace x_k by change between new and old tracer mass:
          x_k = x_k - mass_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,itr)
        
          ! add change in kg to tracer arrays, adjust slopes:
          do l = 1, lmloc
            do j = 1, jmr
              do i = 1, imr
                ! adjust tracer arrays for local indices:
                call AdjustTracer( x_k(i,j,l), region, i, j, l, itr, status )
                IF_NOTOK_RETURN(status=1)
              end do
            end do
          end do
            
          ! clear
          deallocate( x_k )

        ! error ...
        case default

          write (gol,'("unsupported distribution:")'); call goPr
          write (gol,'("  region         : ",i2)') region ;call goPr
          write (gol,'("  previous_par   : ",a)') previous_par(region) ;call goPr
          TRACEBACK; status=1; return

      end select

    else if ( (itr > ntracet) .and. (itr <= ntrace) ) then
    
      ! storage for tracer change:
      allocate( x_k(imr,jmr,lmloc) )
      
      ! copy input:
      x_k = rm_k

      ! rm_k was in given unit; convert x_k to kg :
      select case ( unit )
        case ( 'kg' )
          ! x_k already in kg
        case ( 'mass-mixing-ratio' )
          ! mass mixing ratio to tracer mass:
          x_k = x_k * m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc)
        case ( 'volume-mixing-ratio' )
          ! volume mixing ratio = mass mixing ratio * fscale
          x_k = x_k * m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) / fscale(itr)
        case default
          write (gol,'("unsupported unit for non-transported tracer slab : ",a)') unit; call goErr
          TRACEBACK; status=1; return
      end select

      ! replace short-lived tracer:
      chem_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,itr) = x_k  ! kg
      
      ! clear
      deallocate( x_k )

    else
    
      write (gol,'("unsupported tracer index : ",i4)') itr; call goErr
      write (gol,'("  ntrace, ntracet : ",2i4)') ntrace, ntracet; call goErr
      TRACEBACK; status=1; return

    end if  ! transported or chem-only

    ! ok
    status = 0
    
  end subroutine Tracer_Fill_Slabs

#endif


  ! ***
  
  
  !
  ! Fill array on this pe with tracer (all levels)
  !  
  
#ifdef with_mpi 
  subroutine Tracer_Get_Here( region, itr, unit, rm_t, dest_id, status )
    use mpi_const  , only : MPI_DOUBLE_PRECISION, MPI_STATUS_SIZE, MPI_SUCCESS
    use dims       , only : im, jm, lm
    use chem_param , only : ntrace, ntracet, fscale
    use partools   , only : myid, root
    use partools   , only : previous_par, tracer_loc, tracer_id, tracer_active
    use partools   , only : localComm
!    use partools   , only : lmloc, lmar, offsetl
!    use partools   , only : Par_Scatter_Over_Levels
    use meteodata  , only : m_dat
    
    ! --- in/out -------------------------------------
    
    integer, intent(in)             ::  region
    integer, intent(in)             ::  itr
    character(len=*), intent(in)    ::  unit
    real(8), intent(out)            ::  rm_t(:,:,:)
    integer, intent(in)             ::  dest_id
    integer, intent(out)            ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Tracer_Get_Here'
    
    ! --- local ---------------------------------
    
    integer             ::  imr, jmr, lmr
    integer             ::  itr_loc
    integer             ::  mpi_status(MPI_STATUS_SIZE)
    integer             ::  mpi_tag
    
    ! --- begin ----------------------------------

    ! local grid size    
    imr = im(region)
    jmr = jm(region)
    lmr = lm(region)
    
    ! check input:
    if ( any(shape(rm_t)/=(/imr,jmr,lmr/)) ) then
      write (gol,'("shape of output grid not ok:")'); call goErr
      write (gol,'("  region         : ",i4)') region; call goErr
      write (gol,'("  imr,jmr,lmr    : ",3i4)') imr,jmr,lmr; call goErr
      write (gol,'("  shape(rm_t)    : ",3i4)') shape(rm_t); call goErr
      TRACEBACK; status=1; return
    end if

    ! transported or short lived tracer ?
    if ( (itr >= 1) .and. (itr <= ntracet) ) then
      
      ! how distributed ?
      select case ( previous_par(region) )

        case ( 'tracer' )
        
          ! local tracer index:
          itr_loc = tracer_loc(itr)

          ! already on destination pe ?
          if ( tracer_id(itr) == dest_id ) then
            ! is this the destination pe ?
            if ( myid == dest_id ) then
              ! copy into destination:
              rm_t = mass_dat(region)%rm_t(1:imr,1:jmr,1:lmr,itr_loc)
            end if
          else
            ! transfer from source to destination pe ...
            ! set some tag:
            mpi_tag = 1234
            ! is this the source or destination pe ?
            if ( myid == tracer_id(itr) ) then
              ! this is the source; send to destination:
              call MPI_Send( mass_dat(region)%rm_t(1:imr,1:jmr,1:lmr,itr_loc), &
                             imr*jmr*lmr, MPI_DOUBLE_PRECISION, &
                             dest_id, mpi_tag, localComm, status )
              IF_MPI_NOTOK_RETURN(status=1)
            else if ( myid == dest_id ) then
              ! this is the destination; receive from source:
              call MPI_Recv( rm_t, imr*jmr*lmr, MPI_DOUBLE_PRECISION, &
                             tracer_id(itr), mpi_tag, localComm, mpi_status, status )
              IF_MPI_NOTOK_RETURN(status=1)
            end if
          end if

          ! rm_t now filled on destination; post process if necessary:
          if ( myid == dest_id ) then
            ! change unit ?
            select case ( unit )
              case ( 'kg' )
                ! rm already in kg
              case ( 'mass-mixing-ratio' )
                ! mass mixing ratio
                rm_t = rm_t / m_dat(region)%data(1:imr,1:jmr,1:lmr)
              case ( 'volume-mixing-ratio' )
                ! volume mixing ratio = mass mixing ratio * fscale
                rm_t = rm_t / m_dat(region)%data(1:imr,1:jmr,1:lmr) * fscale(itr)
              case default
                write (gol,'("unsupported unit for par tracer : ",a)') unit; call goErr
                TRACEBACK; status=1; return
            end select
          end if
    
!        case ( 'levels' )
!        
!          ! array already distributed over levels; just copy:
!          rm_k = mass_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,itr)
!
!          ! change unit ?
!          select case ( unit )
!            case ( 'kg' )
!              ! rm already in kg
!            case ( 'mass-mixing-ratio' )
!              ! mass mixing ratio
!              rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc)
!            case ( 'volume-mixing-ratio' )
!              ! volume mixing ratio = mass mixing ratio * fscale
!              rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) * fscale(itr)
!            case default
!              write (gol,'("unsupported unit for par level : ",a)') unit; call goErr
!              TRACEBACK; status=1; return
!          end select

        case default

          write (gol,'("unsupported distribution:")'); call goPr
          write (gol,'("  region         : ",i2)') region ;call goPr
          write (gol,'("  previous_par   : ",a)') previous_par(region) ;call goPr
          TRACEBACK; status=1; return

      end select

    else if ( (itr > ntracet) .and. (itr <= ntrace) ) then
    
      write (gol,'("routine not implemented for not-transported tracers yet ...")'); call goErr
      TRACEBACK; status=1; return
    
!      ! extract short-lived tracer:
!      rm_k = chem_dat(region)%rm_k(1:imr,1:jmr,1:lmloc,itr)
!      
!      ! change unit ?
!      select case ( unit )
!        case ( 'kg' )
!          ! rm already in kg
!        case ( 'mass-mixing-ratio' )
!          ! mass mixing ratio
!          rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc)
!        case ( 'volume-mixing-ratio' )
!          ! volume mixing ratio = mass mixing ratio * fscale
!          rm_k = rm_k / m_dat(region)%data(1:imr,1:jmr,offsetl+1:offsetl+lmloc) * fscale(itr)
!        case default
!          write (gol,'("unsupported unit for non-transported tracer slab : ",a)') unit; call goErr
!          TRACEBACK; status=1; return
!      end select

    else
    
      write (gol,'("unsupported tracer index : ",i4)') itr; call goErr
      write (gol,'("  ntrace, ntracet : ",2i4)') ntrace, ntracet; call goErr
      TRACEBACK; status=1; return

    end if  ! transported or chem-only

    ! ok
    status = 0
    
  end subroutine Tracer_Get_Here
#endif


end module tracer_data


