!#################################################################
!
! tracer 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"
!
!#################################################################
!
!     Declare non-transported tracers : chem_data (type), chem_dat (var)
!     Declare     transported tracers : mass_data (type), mass_dat (var)
!     Allocate both transported and non-transported tracers : tracer_init
!     Free allocated tracers : tracer_done
!
MODULE TRACER_DATA

  use GO,            only : gol, goErr, goPr
  use dims,          only : nregions
  use tm5_distgrid,  only : dgrid, Get_DistGrid

  IMPLICIT NONE
  PRIVATE
  
  public  :: mass_data, mass_dat ! type and array for     transported tracers
  public  :: chem_data, chem_dat ! type and array for non-transported tracers

  public  :: tracer_init   ! allocate tracers arrays
  public  :: tracer_done   ! deallocate tracers arrays
  
  public  :: AdjustTracer  
  public  :: tracer_print    ! tracer and air masses in one box (location hardwire)
  public  :: par_check_mass  ! print min/max/total of Air, All-tracers, and one tracer

  public  :: init_short      ! init non-transported tracers (used when istart=4 or 5)
#ifndef __GFORTRAN__
  public  :: init_non_zero   ! init tracer to small non-zero (used when istart=2)
#endif

#ifdef with_feedback
  public  ::  Tracer_Fill_Slabs
#endif
  
  ! --- const --------------------------------------

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


  TYPE MASS_DATA    ! --- TRANSPORTED TRACERS
     !
     ! All to be allocated with halo=2
     !
     ! rm      tracer masses (kg)
     ! rxm     tracer slopes in kg/(halfgridsize)
     ! rym     tracer slopes in kg/(halfgridsize)
     ! rzm     tracer slopes in kg/(halfgridsize)
     !
     real,dimension(:,:,:,:),pointer   :: rm
#ifdef slopes
     real,dimension(:,:,:,:),pointer   :: rxm
     real,dimension(:,:,:,:),pointer   :: rym
     real,dimension(:,:,:,:),pointer   :: rzm
#ifdef secmom
     real,dimension(:,:,:,:),pointer   :: rxxm
     real,dimension(:,:,:,:),pointer   :: rxym
     real,dimension(:,:,:,:),pointer   :: rxzm
     real,dimension(:,:,:,:),pointer   :: ryym
     real,dimension(:,:,:,:),pointer   :: ryzm
     real,dimension(:,:,:,:),pointer   :: rzzm
#endif
#endif
     integer :: halo
  END TYPE MASS_DATA

  
  TYPE CHEM_DATA    ! --- NON-TRANSPORTED TRACERS
     real,dimension(:,:,:,:),pointer :: rm
     integer :: halo
  END TYPE CHEM_DATA
  
  ! --- VAR ------------------------------
  
  ! transported tracer masses
  type(mass_data),  target :: mass_dat(nregions)
  
  ! non-transported tracer masses
  ! Note expected 4th dim : (:,:,:,ntracet+1:ntracet+ntrace_chem)
  type(chem_data), target   :: chem_dat(nregions) 

CONTAINS
  !
  ! ======================================================================
  !
  SUBROUTINE TRACER_INIT

    use dims,          only : lm
    use chem_param,    only : ntracet, ntrace_chem

    integer :: i1,  i2,  j1,  j2
    integer :: n, l_halo

    l_halo = 2 ! to adapt accordingly
    
    ! allocate transported tracers
    do n=1, nregions

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

       mass_dat(n)%halo = l_halo
       
       allocate( mass_dat(n)%rm( i1-l_halo:i2+l_halo, j1-l_halo:j2+l_halo, lm(n), ntracet) )
#ifdef slopes
       allocate( mass_dat(n)%rxm( i1-l_halo:i2+l_halo, j1-l_halo:j2+l_halo, lm(n), ntracet) )
       allocate( mass_dat(n)%rym( i1-l_halo:i2+l_halo, j1-l_halo:j2+l_halo, lm(n), ntracet) )
       allocate( mass_dat(n)%rzm( i1-l_halo:i2+l_halo, j1-l_halo:j2+l_halo, lm(n), ntracet) )
#ifdef secmom
       allocate( mass_dat(n)%rxxm( i1-l_halo:i2+l_halo, j1-l_halo:j2+l_halo, lm(n), ntracet) )
       allocate( mass_dat(n)%rxym( i1-l_halo:i2+l_halo, j1-l_halo:j2+l_halo, lm(n), ntracet) )
       allocate( mass_dat(n)%rxzm( i1-l_halo:i2+l_halo, j1-l_halo:j2+l_halo, lm(n), ntracet) )
       allocate( mass_dat(n)%ryym( i1-l_halo:i2+l_halo, j1-l_halo:j2+l_halo, lm(n), ntracet) )
       allocate( mass_dat(n)%ryzm( i1-l_halo:i2+l_halo, j1-l_halo:j2+l_halo, lm(n), ntracet) )
       allocate( mass_dat(n)%rzzm( i1-l_halo:i2+l_halo, j1-l_halo:j2+l_halo, lm(n), ntracet) )
#endif
#endif
    
       mass_dat(n)%rm =0.0
#ifdef slopes
       mass_dat(n)%rxm=0.0
       mass_dat(n)%rym=0.0
       mass_dat(n)%rzm=0.0
#ifdef secmom
       mass_dat(n)%rxxm=0.0
       mass_dat(n)%rxym=0.0
       mass_dat(n)%rxzm=0.0
       mass_dat(n)%ryym=0.0
       mass_dat(n)%ryzm=0.0
       mass_dat(n)%rzzm=0.0
#endif
#endif
       ! allocate non-transported tracers, if any
       if ( ntrace_chem > 0 ) then

          chem_dat(n)%halo = 0
          allocate( chem_dat(n)%rm(i1:i2, j1:j2, lm(n), ntracet+1:ntracet+ntrace_chem) )
          chem_dat(n)%rm = 0.0

       else       
          nullify( chem_dat(n)%rm )       
       end if

    end do
       
  END SUBROUTINE TRACER_INIT
  !
  ! ======================================================================
  !
  SUBROUTINE TRACER_DONE

    USE CHEM_PARAM,    ONLY : NTRACE_CHEM

    integer :: n

    ! allocate transported tracers
    do n=1, nregions

       ! deallocate transported tracers
       if(associated( mass_dat(n)%rm)) nullify( mass_dat(n)%rm)
#ifdef slopes
       if(associated( mass_dat(n)%rxm)) nullify( mass_dat(n)%rxm)
       if(associated( mass_dat(n)%rym)) nullify( mass_dat(n)%rym)
       if(associated( mass_dat(n)%rzm)) nullify( mass_dat(n)%rzm)
#ifdef secmom
       if(associated( mass_dat(n)%rxxm)) nullify( mass_dat(n)%rxxm)
       if(associated( mass_dat(n)%rxym)) nullify( mass_dat(n)%rxym)
       if(associated( mass_dat(n)%rxzm)) nullify( mass_dat(n)%rxzm)
       if(associated( mass_dat(n)%ryym)) nullify( mass_dat(n)%ryym)
       if(associated( mass_dat(n)%ryzm)) nullify( mass_dat(n)%ryzm)
       if(associated( mass_dat(n)%rzzm)) nullify( mass_dat(n)%rzzm)
#endif
#endif

       ! deallocate non-transported tracers, if any
       if ( ntrace_chem > 0 ) then

          if(associated( chem_dat(n)%rm)) nullify( chem_dat(n)%rm)

       end if
       
    end do
    
  END SUBROUTINE TRACER_DONE  
  
  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    INIT_SHORT
  !
  ! !DESCRIPTION: Initialise short lived chemical compounds to a reasonable
  !               value e.g NO2 = NOX. It is called from initexit/start when istart=4 or 5.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE INIT_SHORT( region )
    !
    ! !USES:
    !
    use tm5_distgrid, only : dgrid, Get_DistGrid
    use chem_param,   only : ino2,  inox
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: region
    !
    ! !REVISION HISTORY: 
    !    7 May 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    integer :: i1,j1,j2,i2

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

    if ( (ino2 > 0) .and. (inox > 0)) then
       chem_dat(region)%rm(i1:i2, j1:j2, :, ino2) = mass_dat(region)%rm(i1:i2, j1:j2, :, inox)
    endif

  END SUBROUTINE INIT_SHORT
  !EOC

#ifndef __GFORTRAN__
  !------------------------------------------------------------------------------
  !                    TM5                                                      !
  !------------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    INIT_NON_ZERO
  !
  ! !DESCRIPTION: Set tracer initial values to some small value.
  !               It is called from initexit/start when ISTART = 2. 
  !\\
  !\\
  SUBROUTINE INIT_NON_ZERO
    !
    ! !USES:
    !
    use chem_param,   only : fscale, ntracet, ntrace, ntrace_chem, ino2, inox
    use meteodata,    only : m_dat
    !
    ! !REVISION HISTORY: 
    !    2 Apr 2010 - P. Le Sager - 
    !
    ! !REMARKS:
    !  - A remnant from old debugging/benchmark. Has one advantage: can be
    !     used to set to non-zero start values, which is important for chemistry.
    !  - Could be set to values based on process number for debugging. See commented code.
    !EOP
    !------------------------------------------------------------------------------
    !BOC

    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
    integer                           :: rank, i0, i1, j0, j1, n, region
    real, dimension(:,:,:), pointer   :: m

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

    do region=1,nregions

       call Get_DistGrid( dgrid(region), I_STRT=i0, I_STOP=i1, J_STRT=j0, J_STOP=j1 )

       ! pseudo-random number ( for now just retrieve rank )
       rank = ( i1 / (i1-i0) ) * ( j1 / (j1-j0) ) - 1

       m => m_dat(region)%data

       do n=1,ntracet

          mass_dat(region)%rm(:,:,:,n) = 1e-30 * m / fscale(n)

          !--for something not totally uniform and processor dependant
          !           ! fill transported tracers
          !           do n=1,ntracet
          !              mass_dat(region)%rm( i0:i1, j0:j1,:,n) = (rank+1)*1e-9*m/fscale(n)
          !           enddo

       enddo

       ! non-transported tracers
       if ( ntrace_chem > 0 ) then
          do n=ntracet+1,ntrace
             chem_dat(region)%rm(:,:,:,n) = 1e-30 * m/fscale(n)
          end do
       endif

       if ( (ino2 > 0) .and. (inox > 0)) then
          chem_dat(region)%rm(:,:,:,ino2) = mass_dat(region)%rm(:,:,:,inox)
       endif

#ifdef slopes
       mass_dat(region)%rxm = 0.0
       mass_dat(region)%rym = 0.0
       mass_dat(region)%rzm = 0.0
#ifdef secmom
       mass_dat(region)%rxxm = 0.0
       mass_dat(region)%rxym = 0.0
       mass_dat(region)%rxzm = 0.0
       mass_dat(region)%ryym = 0.0
       mass_dat(region)%ryzm = 0.0
       mass_dat(region)%rzzm = 0.0
#endif
#endif

       nullify(m)

    enddo

    !write(gol,*) ' rm initialized at mixing ratio of 1e-30' ; call goPr

  END SUBROUTINE INIT_NON_ZERO
  !EOC
#endif
  
  
  !
  ! 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       |
  !           ---------------       ---------------
  !
  !  15 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  ! 

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

    use GO      , only : gol, goErr 

    ! --- 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
    rm   => mass_dat(region)%rm
#ifdef slopes       
    rxm  => mass_dat(region)%rxm
    rym  => mass_dat(region)%rym
    rzm  => mass_dat(region)%rzm
#ifdef secmom       
    rxxm => mass_dat(region)%rxxm
    rxym => mass_dat(region)%rxym
    rxzm => mass_dat(region)%rxzm
    ryym => mass_dat(region)%ryym
    ryzm => mass_dat(region)%ryzm
    rzzm => mass_dat(region)%rzzm
#endif
#endif

    ! 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
  

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    TRACER_PRINT
  !
  ! !DESCRIPTION: Debug tool. Print value of tracer mass and air mass in one
  !               box, with message at the begining
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE TRACER_PRINT( region, message, status )
    !
    ! !USES:
    !
    use GO
    use dims,       only : im, jm, lm
    use chem_param, only : ntrace, ntracet, names
    use meteodata,  only : m_dat,  phlb_dat
    !
    ! !INPUT PARAMETERS:
    !
    character(len=*), intent(in)  :: message
    integer,          intent(in)  :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer,          intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !   15 Feb 2012 - P. Le Sager
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/tracer_print'
    integer  :: ii, jj, ll, kk, i1, j1, i2, j2
    real     :: pat, patx, paty, patz, patm, patp
    
    ! set those for the box/tracer used to debug, revert to (1,1,1,1) when you are done
    ii=1
    jj=1
    ll=1
    kk=1

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

    if((ii<i1).or.(ii>i2).or.(jj<j1).or.(jj>j2)) then 
       status=0
       return
    endif

    if (kk <= ntracet) then
       pat  = mass_dat(region)%rm (ii,jj,ll,kk)
#ifdef slopes
       patx = mass_dat(region)%rxm(ii,jj,ll,kk)
       paty = mass_dat(region)%rym(ii,jj,ll,kk)
       patz = mass_dat(region)%rzm(ii,jj,ll,kk)
#endif
    else if (kk <= ntrace) then
       pat = chem_dat(region)%rm(ii,jj,ll,kk)
#ifdef slopes
       patx = 0.                                  
       paty = 0.                                  
       patz = 0.                                  
#endif
    else
       write(gol,*) "out of range tracer index"; call goErr
       status=1
       IF_NOTOK_RETURN(status=1)
    end if
    
    patm =    m_dat(region)%data(ii,jj,ll)
    patp = phlb_dat(region)%data(ii,jj,ll)

    ! Decide what to write:
#ifdef slopes
    write (gol,*) message//" "//names(kk), pat, patx, paty, patz, patm ; call goBug
#else
    write (gol,*) message//" "//names(kk), pat, patm ; call goBug
#endif
    ! ok
    status = 0
    
  END SUBROUTINE TRACER_PRINT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PAR_CHECK_MASS
  !
  ! !DESCRIPTION: Debug tool. Print min, max, and sum of air mass and tracer
  !               mass. 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PAR_CHECK_MASS( region, text, tr_id, debug )
    !
    ! !USES:
    !
    use meteodata,    only : m_dat
    use tm5_distgrid, only : reduce
    use partools,     only : isRoot
    !
    ! !INPUT PARAMETERS:
    !
    integer,           intent(in) :: region
    character(len=*),  intent(in) :: text
    integer, optional, intent(in) :: tr_id
    logical, optional, intent(in) :: debug  ! if true, REDUCE prints min/maxloc
    !
    ! !REVISION HISTORY: 
    !   15 Feb 2012 - P. Le Sager - v0
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    real    :: min_one, max_one, tot_one
    integer :: status, trid
    logical :: dbg_

    ! check input
    trid=1
    if(present(tr_id)) trid=tr_id

    dbg_=.false.
    if(present(debug)) dbg_=debug

  
    rm => mass_dat(region)%rm
    m => m_dat(region)%data

    ! air mass
    call REDUCE( dgrid(region), m, m_dat(region)%halo, 'MAX', max_one, status, debug=dbg_)
    call REDUCE( dgrid(region), m, m_dat(region)%halo, 'MIN', min_one, status, debug=dbg_)
    call REDUCE( dgrid(region), m, m_dat(region)%halo, 'SUM', tot_one, status, debug=dbg_)

    if ( isRoot ) then
       write(gol,*) text//' [par_check_mass] AIR min, max, sum: ', min_one, max_one, tot_one ; call goPr
    endif

    ! all transported tracers
    call REDUCE( dgrid(region), rm, mass_dat(region)%halo, 'MAX', max_one, status, debug=dbg_)
    call REDUCE( dgrid(region), rm, mass_dat(region)%halo, 'MIN', min_one, status, debug=dbg_)
    call REDUCE( dgrid(region), rm, mass_dat(region)%halo, 'SUM', tot_one, status, debug=dbg_)

    if ( isRoot ) then
       write(gol,*) text//' [par_check_mass] RM min, max, sum: ', min_one, max_one, tot_one ; call goPr
    endif

    nullify(m)


    ! one selected tracer
    m => mass_dat(region)%rm(:,:,:,trid) ! no need for temp arr since slice is in last dimension
    
    call REDUCE( dgrid(region), m, mass_dat(region)%halo, 'MAX', max_one, status, debug=dbg_)
    call REDUCE( dgrid(region), m, mass_dat(region)%halo, 'MIN', min_one, status, debug=dbg_)
    call REDUCE( dgrid(region), m, mass_dat(region)%halo, 'SUM', tot_one, status, debug=dbg_)

    if ( isRoot ) then
       write(gol,*) text//' [par_check_mass] ID, min, max, sum: ', trid, min_one, max_one, tot_one ; call goPr
    endif

    nullify(m)
    
  END SUBROUTINE PAR_CHECK_MASS
  !EOC

  
  !
  ! Convert and put a slab of tracer mass (with given unit and tracer id)
  ! 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  :: i, j, l, lmr
    integer  :: i01, i02, j01, j02, lmr

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

    ! local grid size
    lmr = lm(region)

    call Get_DistGrid( dgrid(region), I_STRT=i01, I_STOP=i02, &
                                      J_STRT=j01, J_STOP=j02  )
    
    ! check input:
    if ( any(shape(rm_k)/=(/i02-i01+1,j02-j01+1,lmr/)) ) then
       write (gol,'("shape of output grid not ok:")'); call goErr
       write (gol,'("  shape(rm_k)              : ",3i4)') shape(rm_k); call goErr
       write (gol,'("  i02-i01+1,j02-j01+1,lmr  : ",3i4)') i02-i01+1,j02-j01+1,lmr; call goErr
       TRACEBACK; status=1; return
    end if

    
    ! transported or short lived tracer ?
    if ( (itr >= 1) .and. (itr <= ntracet) ) 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)
          rm_k = rm_k / plc_kg_from_tm(itr)

       case ( 'mass-mixing-ratio' )
          ! mass mixing ratio
          ! apply unit conversion from  (kg plc tracer) to (kg tm tracer)
          rm_k = rm_k / plc_kg_from_tm(itr)
          rm_k = rm_k * m_dat(region)%data(i01:i02,j01:j02,:)

       case ( 'volume-mixing-ratio' )
          ! volume mixing ratio = mass mixing ratio * fscale
          ! apply unit conversion from  (kg plc tracer) to (kg tm tracer)
          rm_k = rm_k / plc_kg_from_tm(itr)
          rm_k = rm_k * m_dat(region)%data(i01:i02,j01:j02,:) / fscale(itr)

       case default
          write (gol,'("unsupported unit for par tracer : ",a)') unit; call goErr
          TRACEBACK; status=1; return

       end select

       ! replace new mass rm_k by difference with rm :
       rm_k = rm_k - mass_dat(region)%rm_t(i01:i02,j01:j02,:,itr)   ! kg
       
       ! add change in kg to tracer arrays, adjust slopes:
       do l = 1,   lmr
       do j = j01, j02
       do i = i01, i02
                ! adjust tracer arrays for local indices:
                call AdjustTracer( rm_k(i,j,l), region, i, j, l, itr, status )
                IF_NOTOK_RETURN(status=1)
       end do
       end do
       end do


    else if ( (itr > ntracet) .and. (itr <= ntrace) ) then

       ! 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(i01:i02,j01:j02,:)

       case ( 'volume-mixing-ratio' )
          ! volume mixing ratio = mass mixing ratio * fscale
          x_k = x_k * m_dat(region)%data(i01:i02,j01:j02,:) / 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(:,:,:,itr) = rm_k  ! kg

    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

END MODULE TRACER_DATA