ecearth3/sources/tm5mp/base/redgridZoom.F90

!### macro's #####################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:     REDGRIDZOOM
!
! !DESCRIPTION: 
!\\
!\\
! !INTERFACE: 
!
MODULE REDGRIDZOOM
  !
  ! !USES:
  !
  use GO   , only : gol, goErr, goPr
  use binas, only : pi
  use dims,  only : nregions, im, jm, okdebug

  implicit none
  private
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: uni2red_mf, redgrid_init, redgrid_done, calc_pdiff
  public :: uni2red, red2uni, red2uni_em
#ifdef secmom
  public :: uni2red_2nd, red2uni_em_2nd
#endif
  !
  ! !PUBLIC DATA MEMBERS:
  !
  public :: nred, nredmax, clustsize, grid_reduced, idle_xadv
  public :: jred, imredj
  public :: mfl_alt1, m_alt1, rm_alt1, rxm_alt1, rym_alt1, rzm_alt1

  integer, parameter                   :: nredmax=60
  logical                              :: grid_reduced=.true.
  logical                              :: idle_xadv=.false. ! proc is idle during x-advection (possible if reduced grid with npe_lon/=1)
  ! the number of reduced latitude circles per region on the current proc
  integer, dimension(nregions)         :: nred=0
  ! number of joined cells per latitude circle and region
  integer, allocatable                 :: clustsize(:,:)  ! j1:j2,nregions
  ! reduced dimension......
  integer, allocatable                 :: imredj(:,:)  ! j1:j2,nregions
  ! latitaude numbers where the reduction applies
  integer, dimension(nredmax,nregions) :: jred

  ! alternate mass array (used on row_roots, to deal with reduced grid). With and w/o halo.
  real, dimension(:,:,:),   target, allocatable :: mfl_alt1, m_alt1
  real, dimension(:,:,:,:), target, allocatable :: rm_alt1, rxm_alt1, rym_alt1, rzm_alt1
  real, dimension(:,:,:),   allocatable :: mfl_alt2, m_alt2
  real, dimension(:,:,:,:), allocatable :: rm_alt2

  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter :: mname = 'redgridZoom'
  real,             parameter :: dl = 2.0*pi/im(1)
  real,             parameter :: dp = pi/jm(1)
  real,             parameter :: epsx = epsilon(0.0)
  !
  ! !REVISION HISTORY: 
  !   20 Dec 2012 - Ph Le Sager - TM5-MP version.
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !  Fortran code for Split-rg (reduced grid) scheme for advection on the 
  !  globe.
  !
  !  VERSION:        1.1
  !  DATE:           November 12, 1998
  !
  !  Version 1.0     Dated October 19, 1998.
  !  Version 1.1     Corrected red2uni routine and added new uni2red_m 
  !                  routine.
  !                  (Orography effects were erroneously not taken into 
  !                  account in the reduced to uniform grid conversion)
  !
  !  This code solves the discrete advection equation for the tracer
  !  mass, with the wind field specified in air mass fluxes:
  !
  !  dtracm/dt = (
  !   mu chi (i-half) - mu chi (i+half)
  !   mv chi (j-half) - mv chi (j+half)
  !   mw chi (k-half) - mw chi (k+half) )
  !
  !  where
  !
  !  tracm = tracer mass
  !  mu, mv, mw = air mass fluxes
  !  chi = tracm / mass = tracer mixing ratio
  !  mass = air mass
  !
  !  The integer values of indices i, j, and k correspond to cell centers.
  !
  !  See:
  !
  !  Petersen, A.C., E.J. Spee, H. van Dop, and W. Hundsdorfer,
  !  An evaluation and intercomparison of four new advection schemes
  !  for use in global chemistry models,
  !  Journal of Geophysical Research, 103, 19253-19269, 1998.
  !
  !  Written by Edwin Spee, CWI, Amsterdam, The Netherlands
  !
  !  Parts of the code are written by
  !  Joke Blom and Willem Hundsdorfer (CWI) and
  !
  !  The advection routine is called in the following way:
  !
  !     call advect_split(tracm, mu, mv, mw, mass, dt)
  !
  !  tracm = tracer mass            real(im,jm,lm,ntracet) [kg]
  !
  !  mu = longitudinal mass flux    real(im,jm,lm) [kg/s]
  !                                 (defined on eastward side of grid
  !                                 cells; for all grid cells)
  !
  !  mv = latitudinal mass flux     real(im,jm,lm) [kg/s]
  !                                 (defined on southward side of grid
  !                                 cells; only j=2,..,jm are used)
  !
  !  mw = vertical massflux         real(im,jm,lm) [kg/s]
  !                                 (defined on upper side of grid cell;
  !                                 only l=1,..,lm-1 are used since it is 
  !                                 assumed that there is no flux across 
  !                                 the top of the model)
  !
  !  The number of latitudinal elements actually used for mu, mv, and mw
  !  depends on the advection grid; for mv this number is the determined
  !  by the latitude row with the largest number of cells bordering a
  !  certain cell-bounding latitude circle. 
  !
  !  mass = air mass                real(im,jm,lm) [kg]
  !
  !  dt   = timestep                real [s]
  !
  !  The physical units were chosen for compatibility with the TM3 model
  !  but can be changed for other global tracer models, provided that
  !  the calculation of the fluxes is done with an equivalent of tracer
  !  mixing ratios and not of tracer concentrations (make for other tracer
  !  models, if necessary, the appropriate changes in the routines 
  !  split_lab, split_phi, and split_eta, where now tracer mixing ratios 
  !  are calculated as tracm / mass; also the grid transformations as done
  !  in advect_split are different for tracer mass and concentration).
  !
  !  The longitude is counted from west to east, the latitude from south 
  !  to north, and the height from surface to top.
  !
  !  The variables tracm, mu, mv, mw and mass should be declared in
  !  the calling program under any name; dt is a  parameter. The 
  !  parameters im, jm, and lm denote the number of grid cells in 
  !  longitudinal, latitudinal, and vertical direction, respectively. 
  !  ntracet is the number of transported tracers. The letters 'u', 'v', 
  !  and 'w' for the three dimensions are in the code often replaced by 
  !  'l', 'p', and 'e', denoting the coordinates lambda, phi, and eta, 
  !  respectively (see Petersen et al. 1998). The advection scheme works 
  !  in principle with any grid and coordinate system of the main model. 
  !  However, in the current implementation (made for the TM3 model) the 
  !  variables tracm and mass are assumed to be defined in the main 
  !  program on a uniform grid with polar cap, and the variables mu, mv, 
  !  and mw on the advection grid (i.e., reduced grid with polar cap). 
  !  Locally in the advection routines a transformation of tracm and mass,
  !  if necessary, is made to the advection grid. Please contact us for 
  !  information on the changes that have to be made to work with other 
  !  main model grids. Depending on the interests of users, we intend to 
  !  generalize the code with respect to the main model grid in the next 
  !  version.
  !
  !  The following parameters and variables, related to the specification
  !  of the advection grid are used by the routines of the advection 
  !  scheme:
  !
  !     integer, parameter :: nredmax = 10
  !     real, parameter    :: dl = 2.0 * pi / im
  !     real, parameter    :: dp = pi / jm
  !     integer               nred
  !     integer               clustsize(0:nredmax),
  !                           jred(-nredmax:nredmax+1),imredj(jm)
  !     real                  maxtau
  !     integer               advsub
  !
  !   nredmax   = maximum number of grid reductions per hemisphere
  !               (in this implementation the minimum number of grid
  !               reductions per hemisphere equals one, since use is
  !               made of polar caps)
  !   dl        = longitudinal uniform grid cell size in lambda 
  !               coordinates
  !   dp        = latitudinal grid cell size in phi coordinates
  !   nred      = number of grid reductions
  !   clustsize = array with the ratio im / imred, as a function of the 
  !               absolute value of the latitude zone index, where imred 
  !               is the actually used number of cells in the longitudinal
  !               direction for the given latitude zone
  !   jred      = latitude number where grid reduction starts (seen from 
  !               SP), as a function of the reduction zone index
  !   imredj(j) = number of cells for latitude #j
  !   maxtau    = time step such that CFL <= 1 in all directions (given
  !               the directional splitting, see routine advect_split)
  !   advsub    = largest integer such that dtadv / advsub <= maxtau,
  !               with dtadv the advection time step in operator splitting
  !
  !  The advection grid variables for the reduced grid are automatically 
  !  set by the following line:
  !
  !     call redgrid_init(clustsize, jred, imredj, nred, nredmax, im, jm,
  !   .                  pi, .false.)
  !
  !  The variables maxtau and advsub should be set each time new mass 
  !  fluxes (wind fields) are available in the model through
  !
  !     call split_maxtau(mu, mv, mw, mass, dtadv)
  !  
  !  where dtadv is the advection time step in the main model.
  !
  !  Furthermore, some routines in the advection scheme need information
  !  on the machine precision. The following variable is used for this
  !  purpose:
  !
  !     real, parameter    :: epsx = epsilon(0.0)
  !
  !  The value of epsx is dependent on the computing platform that is used.
  !  epsilon is a Fortran 90 intrinsic function. For example, for our Cray
  !  C90 eps is about 1e-14.
  !
  !  The above-mentioned parameters and variables are assumed to be 
  !  declared in the module global_data (Fortran 90) or in some common 
  !  block (Fortran 77). The current implementation (Fortran 90) uses the 
  !  statements
  !
  !     use global_data
  !     implicit none
  !
  !  If Fortran 77 is preferred one should define a common block or
  !  use an existing common block to include the above-given parameters
  !  and variables, e.g.
  !
  !     implicit none
  !     (declarations as given in the above)
  !     common/advgrid/nredmax,dl,dp,nred,clustsize,jred,maxtau,advsub
  !
  !  Also the parameters im, jm, lm, ntracet, pi, twopi (= 2 * pi), and 
  !  eps are assumed to be declared in the module global_data.
  !
  !  --------------------------
  !  Example of a reduced grid
  !  --------------------------
  ! 
  !  NP=North Pole
  !  SP=South Pole
  !  EQ=Equator
  !  | grid boundary in longitudal direction
  !  o cell center
  !                                 j   imredj   clustsize   latitude zone
  !   NP|           o           |  13     1        12             3
  !     |   o   |   o   |   o   |  12     3         4             2
  !     | o | o | o | o | o | o |  11     6         2             1
  !     | o | o | o | o | o | o |  10     6         2             1
  !     |o|o|o|o|o|o|o|o|o|o|o|o|   9    12         1             0
  !     |o|o|o|o|o|o|o|o|o|o|o|o|   8    12         1             0
  !   EQ|o|o|o|o|o|o|o|o|o|o|o|o|   7    12         1             0
  !     |o|o|o|o|o|o|o|o|o|o|o|o|   6    12         1             0
  !     |o|o|o|o|o|o|o|o|o|o|o|o|   5    12         1             0
  !     | o | o | o | o | o | o |   4     6         2            -1
  !     | o | o | o | o | o | o |   3     6         2            -1
  !     |   o   |   o   |   o   |   2     3         4            -2
  !   SP|           o           |   1     1        12            -3
  !
  !  This grid would have
  !  im = 12, jm = 13, nred = 3
  !
  !  jred(-3) = 1 by definition
  !  jred(-2) = 2
  !  jred(-1) = 3
  !  jred( 0) = 5
  !  jred( 1) =10
  !  jred( 2) =12
  !  jred( 3) =13
  !  jred( 4) =14=jm+1 by definition
  !
  !  clustsize(3) = 12
  !  clustsize(2) =  4
  !  clustsize(1) =  2
  !  clustsize(0) =  1

  !CMK total new definition to allow zoom regios..
  !  nred(nregions)    :: number of reduced latitudes circles   < jm(region)
  !  nredmax           :: maximum number of reduced latitude circles...
  !  clustsize(nredmax.nregions)  :: number of joined cells for n'th circle
  !  jred(nredmax.nregions)       :: corresponding j value in 1...jm(region) array
  !  imredj(nredmax.nregions)     :: 'new' im value for the reduced grid...
  !CMK

  !PLS new definitions for TM5-MP
  !  nred(nregions)    :: number of reduced latitudes circles on current tile (ie proc) and region
  !  nredmax           :: maximum number of reduced latitude circles
  !  clustsize(j1:j2,   nregions) :: number of joined cells for J'th circle (default=1=no_cluster=no_reduced_grid)
  !  jred     (nredmax, nregions) :: j value for 1...nred(region) in 1..jmr
  !  imredj   (j1:j2,   nregions) :: 'new' im value for the reduced grid, default to im(region)
  !PLS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  uni2red_mf
  !
  ! !DESCRIPTION:  Converts mass fluxes mfl from uniform grid to reduced grid
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine uni2red_mf(region,i1,i2,j1,j2, status)
    !
    ! !USES:
    !
    use dims,         only : lm, im
    use global_data,  only : wind_dat
    use partools,     only : isRowRoot, npe_lon
    use tm5_distgrid, only : dgrid, gather_j_band
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in) :: region,i1,i2,j1,j2
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)  :: status
    !
    ! !REVISION HISTORY: 
    !   Written by Edwin Spee, CWI, Amsterdam, The Netherlands
    !   cmk changed.....
    !   20 Dec 2012 - Ph Le Sager - TM5-MP version
    !
    ! !REMARKS:
    !  - no need to check on tile bounds, zero-trip loop for jreg does the job
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/uni2red_mf'
    
    real,dimension(:,:,:),pointer :: mfl, mfl_alt
    integer :: clust_size, i, j, l, jreg
    integer :: lmr, imr, halo

    ! -- start
    
    mfl => wind_dat(region)%am
    halo=wind_dat(region)%halo

    imr=im(region)
    lmr=lm(region)

    ! Two scenarios
    if (npe_lon /= 1) then
       !
       ! decomposition along longitudes => work on row_root
       !
       if (nred(region)/=0) then  ! some reduce grid on this tile

          ! array to gather (remove halo)
          allocate(mfl_alt( i1:i2, j1:j2, 0:lmr+1))
          mfl_alt = mfl(i1:i2,j1:j2,:)

          CALL GATHER_J_BAND(dgrid(region), mfl_alt, mfl_alt2, status, jref=j1, rowcom=.true.)
          IF_NOTOK_RETURN(status=1)
          
          if (isRowRoot) then

             mfl_alt1(1:imr,:,:) = mfl_alt2
             mfl_alt1(    0,:,:) = mfl_alt1(imr,:,:)
             mfl_alt1(   -1,:,:) = mfl_alt1(imr-1,:,:)
             mfl_alt1(imr+1,:,:) = mfl_alt1(1,:,:)
             !PLS mfl_alt1(imr+2,:,:) = mfl_alt1(2,:,:)  ! full on

             do l=1,lmr
                do jreg = 1,nred(region) ! loop over the reduced latitudes

                   j = jred(jreg,region)            ! latitude
                   clust_size = clustsize(j,region) ! clustersize

                   do i=1,imredj(j,region)
                      mfl_alt1(i,j,l) = mfl_alt1(i*clust_size,j,l) ! compress mfl
                   end do
                   do i=imredj(j,region)+1,im(region)
                      mfl_alt1(i,j,l) = -1.0
                   end do

                   ! Update boundary condition. Only meaningful for x-cyclic regions (which is the
                   !  only possibility in TM5-MP).
                   mfl_alt1( 0,j,l)=mfl_alt1(imredj(j,region),j,l)
                   mfl_alt1(-1,j,l)=mfl_alt1(imredj(j,region)-1,j,l) ! needed to limit mpi comm in advectx__slopes.F90
                   mfl_alt1(imredj(j,region)+1,j,l)=mfl_alt1(1,j,l)  ! needed to limit mpi comm in advectx__slopes.F90
                end do
             end do             
             
          endif
          deallocate( mfl_alt )
       endif

    else
       !
       ! Reduced grid without longitudinal decomposition
       !
       do l=1,lmr
          do jreg = 1,nred(region) ! loop over the reduced latitudes
             
             j = jred(jreg,region)            ! latitude
             clust_size = clustsize(j,region) ! clustersize

             do i=1,imredj(j,region)
                mfl(i,j,l) = mfl(i*clust_size,j,l) ! compress mfl
             end do
             do i=imredj(j,region)+1,im(region)
                mfl(i,j,l) = -1.0
             end do

             ! Update boundary condition. Only meaningful for x-cyclic regions (which is the
             !  only possibility in TM5-MP).
             mfl( 0,j,l)=mfl(imredj(j,region),j,l)
             mfl(-1,j,l)=mfl(imredj(j,region)-1,j,l) ! needed to limit mpi comm in advectx__slopes.F90
             mfl(imredj(j,region)+1,j,l)=mfl(1,j,l)  ! needed to limit mpi comm in advectx__slopes.F90
          end do
       end do

    endif

    nullify(mfl)
    status=0
    
  end subroutine uni2red_mf
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    REDGRID_INIT
  !
  ! !DESCRIPTION: allocate work arrays, read settings from rc file and fill
  !               reduced grid parameters (clustsize, jred, and imredj) for
  !               current region and tile.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE REDGRID_INIT( region, status )
    !
    ! !USES:
    !
    use GO,           only : TrcFile, Init, Done, ReadRc
    use dims,         only : lm, im
    use global_data,  only : wind_dat
    use tracer_data,  only : mass_dat
    use meteodata  ,  only : m_dat
    use chem_param,   only : ntracet
    use tm5_distgrid, only : dgrid, Get_DistGrid
    use partools,     only : isRowRoot, npe_lon
    use global_data,  only : rcfile
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !   Written by Edwin Spee, CWI, Amsterdam, The Netherlands, based on work of Joke Blom.
    !   20 Dec 2012 - Ph Le Sager - TM5-MP
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer             :: j, i1, i2, j1, j2, imr, lmr, halo
    character(len=256)  :: fname, line
    type(TrcFile)       :: rcF

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

    ! Open rcfile
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    ! grid reduction ?
    call ReadRc( rcF, 'proces.advection.reduced', grid_reduced, status )
    IF_NOTOK_RETURN(status=1)
    
    if (grid_reduced) then 
       
       ! Arrays
       CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

       allocate( clustsize(j1:j2, nregions))
       allocate(    imredj(j1:j2, nregions))
    
       call Read_RedGrid_rc( region, status )
       IF_NOTOK_RETURN(status=1)
       
       ! Extra arrays to work on row_root (only if really needed)
       if (nred(region)/=0 .and. (npe_lon/=1)) then
          if (isRowRoot) then
             imr = im(region)
             lmr = lm(region)
       
             halo= wind_dat(1)%halo
             allocate(mfl_alt1( 1-halo:imr+halo, j1:j2, 0:lmr+1))
             allocate(mfl_alt2(             imr, j1:j2, 0:lmr+1))

             halo= m_dat(1)%halo
             allocate(  m_alt1( 1-halo:imr+halo, j1:j2, lmr))
             allocate(  m_alt2(             imr, j1:j2, lmr))

             halo= mass_dat(1)%halo
             allocate(rm_alt1( 1-halo:imr+halo, j1:j2, lmr, ntracet))
             allocate(rm_alt2(             imr, j1:j2, lmr, ntracet))
#ifdef slopes
             allocate(rxm_alt1( 1-halo:imr+halo, j1:j2, lmr, ntracet))
             allocate(rym_alt1( 1-halo:imr+halo, j1:j2, lmr, ntracet))
             allocate(rzm_alt1( 1-halo:imr+halo, j1:j2, lmr, ntracet))
#endif
          else
             idle_xadv=.true.
             allocate(mfl_alt2(1,1,1)  )
             allocate(  m_alt2(1,1,1)  )
             allocate( rm_alt2(1,1,1,1))
          endif
          
       endif
       
    endif
    
    ! close
    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)

    status = 0

  END SUBROUTINE REDGRID_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    REDGRID_DONE
  !
  ! !DESCRIPTION: deallocate work arrays
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE REDGRID_DONE( status )
    !
    ! !USES:
    !
    use partools,     only : isRowRoot, npe_lon
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: status
    !
    ! !REMARKS: used only for region=1, in a hardcoded way 
    ! 
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

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

    if (grid_reduced) then 

       deallocate( clustsize)
       deallocate(    imredj)

       if (nred(1)/=0 .and. (npe_lon/=1)) then          
          if (isRowRoot) then
             deallocate(mfl_alt1)
             deallocate(mfl_alt2)
             deallocate(  m_alt1)
             deallocate(  m_alt2)
             deallocate( rm_alt1)
             deallocate( rm_alt2)
#ifdef slopes
             deallocate(rxm_alt1)
             deallocate(rym_alt1)
             deallocate(rzm_alt1)
#endif
          else
             deallocate(mfl_alt2)
             deallocate(  m_alt2)
             deallocate( rm_alt2)
          endif
       endif
    endif

    status = 0

  END SUBROUTINE REDGRID_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    READ_REDGRID_RC
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Read_RedGrid_rc( region, status )
    !
    ! !USES:
    !
    use GO,           only : TrcFile, Init, Done, ReadRc
    use global_data,  only : rcfile
    use dims,         only : region_name, ybeg, yend, dy, yref
    use tm5_distgrid, only : dgrid, Get_DistGrid
    use partools,     only : isRoot
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)    ::  region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)   ::  status
    !
    ! !REVISION HISTORY: 
    !   20 Dec 2012 - Ph Le Sager - TM5-MP version
    !    6 Jun 2013 - P. Le Sager - check for consecutive 1-cell rings
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    type(TrcFile)         ::  rcF
    integer               ::  imr, jmr, xam, nim
    integer               ::  i,j, i1, i2, j1, j2
    integer               ::  jband
    integer               ::  nred_nh,      nred_sh
    integer, allocatable  ::  ncombs_nh(:), ncombs_sh(:)
    real                  ::  y0, y1

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

    ! number of lats
    imr = im(region)
    jmr = jm(region)

    ! open settings:
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    ! number of reduced lat bands:
    call ReadRc( rcF, 'region.'//trim(region_name(region))//'.redgrid.nh.n', nred_nh, status )
    IF_NOTOK_RETURN(status=1)

    ! number of reduced lat bands:
    call ReadRc( rcF, 'region.'//trim(region_name(region))//'.redgrid.sh.n', nred_sh, status )
    IF_NOTOK_RETURN(status=1)

    ! total number of reduced bands
    nred(region) = nred_nh + nred_sh

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

    
    ! reduced grid ?
    if ( nred(region) > 0 ) then

       if ( nred(region) > nredmax ) then
          write (gol,'("problem with reduced grid. nred > nredmax")'); call goErr
          TRACEBACK; status=1; return
       end if

       ! test done, reset
       nred(region)=0

       ! to read number of combined cells per lat band
       allocate( ncombs_sh(jmr) )
       ncombs_sh=1
       allocate( ncombs_nh(jmr) )

       ! southern hemisphere
       if ( nred_sh > 0 ) then
          
          call ReadRc( rcF, 'region.'//trim(region_name(region))//'.redgrid.sh.comb', ncombs_sh(1:nred_sh), status )
          IF_NOTOK_RETURN(status=1)

       end if

       ! northern hemisphere
       if ( nred_nh > 0 ) then

          call ReadRc( rcF, 'region.'//trim(region_name(region))//'.redgrid.nh.comb', ncombs_nh(1:nred_nh), status )
          IF_NOTOK_RETURN(status=1)

          do jband= 1, nred_nh
             ncombs_sh(jmr-jband+1) = ncombs_nh(jband)
          end do

       end if

       ! fill work arrays
       clustsize(:,region) = ncombs_sh(j1:j2)
       nred(region)        = count(clustsize(:,region)/=1)
       jred(1:nred(region), region) = pack((/(j,j=j1,j2)/), clustsize(:,region)/=1)

       imredj(:,region) = imr ! default

       do i=1,nred(region)
          j = jred(i, region)
          imredj(j,region) = imr/clustsize(j,region)
       end do
       

       ! Testing (global)
       do j = 1, jmr
          ! check if number of combined cells matches with grid:
          if ( modulo(imr,ncombs_sh(j)) /= 0 ) then
             write (gol,'("number of combined cells not ok:")')                               ; call goErr
             write (gol,'("  region         : ",i2," ",a)') region, trim(region_name(region)) ; call goErr
             write (gol,'("  lat band       : ",i4)') j                                       ; call goErr
             write (gol,'("  combined cells : ",i4)') ncombs_sh(j)                            ; call goErr
             write (gol,'("  im             : ",i4)') imr                                     ; call goErr
             TRACEBACK; status=1; return
          end if
          ! check with previous ...
          if ( j > 1 ) then
             xam=max(ncombs_sh(j-1),ncombs_sh(j))
             nim=min(ncombs_sh(j-1),ncombs_sh(j))
             if ( modulo(xam,nim) /= 0 ) then
                write (gol,'("number of combined cells does match with previous:")')             ; call goErr
                write (gol,'("  region         : ",i2," ",a)') region, trim(region_name(region)) ; call goErr
                write (gol,'("  lat band       : ",i4)') j                                       ; call goErr
                write (gol,'("  combined cells : ",i4)') ncombs_sh(j)                            ; call goErr
                write (gol,'("  previous       : ",i4)') ncombs_sh(j-1)                          ; call goErr
                TRACEBACK; status=1; return
             end if
          end if
          ! check for consecutive rings with a single reduced grid cell
          if ( j > 1 ) then
             if ( ncombs_sh(j-1)==imr .and. ncombs_sh(j)==imr ) then
                write (gol,'("Consecutive rings with a single reduced grid cell detected:")')    ; call goErr
                write (gol,'("  region         : ",i2," ",a)') region, trim(region_name(region)) ; call goErr
                write (gol,'("  lat band       : ",i4)') j                                       ; call goErr
                write (gol,'("  combined cells : ",i4)') ncombs_sh(j)                            ; call goErr
                write (gol,'("  previous       : ",i4)') ncombs_sh(j-1)                          ; call goErr
                TRACEBACK; status=1; return
             end if
          end if
          
       end do

       ! Verbose
       if (isRoot) then
          write (gol,'("   [Reduced grid] Region """,a,""":")') trim(region_name(region)); call goPr
          write (gol,'("   [Reduced grid]    - Uniform grid: no. of cells in each latitude band: ",i3,".")') imr; call goPr
          do j=1,jmr
             if ( ncombs_sh(j) /= 1 ) then
                y0=ybeg(region)+(dy/yref(region))*(j-1)
                y1=ybeg(region)+(dy/yref(region))*j
                write (gol,'("   [Reduced grid]    - from ",f6.2," to ",f6.2," degrees latitude: ",i3," cells.")') &
                     y0,y1,ncombs_sh(j); call goPr
             end if
          end do
       end if

       deallocate( ncombs_sh )
       deallocate( ncombs_nh )

    else
       write (gol,'("No reduced grid for region :",a)') trim(region_name(region)); call goPr
       nred(region) = 0
    end if
    
    
    ! Done
    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)

    status = 0

  END SUBROUTINE READ_REDGRID_RC
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CALC_PDIFF
  !
  ! !DESCRIPTION: calculate max pressure difference, accounting for reduce
  !               grid if used.
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CALC_PDIFF( region, pdiffmax )
    !
    ! !USES:
    !
    use meteodata,    only : sp_dat, sp1_dat
    use tm5_distgrid, only : dgrid, Get_DistGrid
    use partools,     only : isRowRoot, npe_lon
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    real,    intent(out) :: pdiffmax   ! max pressure difference 
    !
    ! !REVISION HISTORY: 
    !   20 Dec 2012 - Ph Le Sager - TM5-MP version
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer :: lrg, i, j, ratio, idx, iu, i1, i2, j1, j2
    real    :: work, work1
    
    real, pointer :: p(:,:,:), pold(:,:,:)

    !---- start
    p    =>  sp_dat(region)%data
    pold => sp1_dat(region)%data
    
    pdiffmax = 0.0

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

    ! if no reduced grid
    if ( .not. grid_reduced ) then

       pdiffmax = maxval(abs(p(i1:i2,j1:j2,1)-pold(i1:i2,j1:j2,1)))

    else
       if (npe_lon==1) then  
          
          ! first maximum over the reduced grid
          do lrg = 1,nred(region)
             j = jred(lrg,region)
             ratio = clustsize(j,region)
             do i =  1,imredj(j,region)
                work = 0.0
                work1 = 0.0
                do iu = 1,ratio
                   idx  = (i-1)*ratio+iu
                   work = work + p(idx,j,1)
                   work1 = work1 + pold(idx,j,1)
                end do
                pdiffmax=max(pdiffmax,abs(work-work1)/ratio)
             end do
          end do

          ! now the pressure difference in the non-reduced part
          j2loop: do j=j1,j2
             if(imredj(j,region).ne.im(region)) cycle j2loop   !if reduced...skip
             i2loop: do i=i1,i2  ! should be same as 1,im(region) by design
                pdiffmax = max(pdiffmax,abs(p(i,j,1)-pold(i,j,1)))
             end do i2loop
          end do j2loop
          
       else
          if (nred(region)/=0) then
             
             jloop:do j=j1,j2
                if(imredj(j,region).ne.im(region)) then
                   ! this requires gathering of sp_dat and sp1_dat across mpi tasks - just skip it (FIXME?)
                   cycle jloop
                else
                   pdiffmax = max( pdiffmax, maxval(abs(p(i1:i2,j,1)-pold(i1:i2,j,1))) )
                endif
             enddo jloop
          else
             pdiffmax = maxval(abs(p(i1:i2,j1:j2,1)-pold(i1:i2,j1:j2,1)))
          endif
       endif
    end if
    
  end subroutine calc_pdiff
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    UNI2RED
  !
  ! !DESCRIPTION: transforms data (air mass & tracers) from uniform grid to reduced grid
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine uni2red( region, i1, i2, j1, j2, status)
    !
    ! !USES:
    ! 
    use dims,         only : im,jm,lm
    use meteodata  ,  only : m_dat
    use tracer_data,  only : mass_dat
    use chem_param,   only : ntracet
    use partools,     only : isRowRoot, npe_lon
    use tm5_distgrid, only : dgrid, gather_j_band
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region,i1,i2,j1,j2
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)  :: status
    !
    ! !REVISION HISTORY: 
    !   written by mike botchev, march-june 1999
    !   modified by Maarten Krol, dec 2002
    !   20 Dec 2012 - Ph Le Sager - TM5-MP
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC
    
    character(len=*), parameter  ::  rname = mname//'_MOD_uni2red'

    real,dimension(:,:,:,:),pointer         :: rm, rm_alt
#ifdef slopes
    real,dimension(:,:,:,:),pointer         :: rxm,rym,rzm
#endif    
    real,dimension(:,:,:),  pointer         :: m, m_alt
    
    integer :: i,ie,is,j,l,lrg,redfact,n, imr, lmr, halo
    real    :: summ,sumrm
    
    ! -------------- start

    m   => m_dat(region)%data
    rm  => mass_dat(region)%rm
#ifdef slopes
    rxm => mass_dat(region)%rxm
    rym => mass_dat(region)%rym
    rzm => mass_dat(region)%rzm
#endif
    
    imr=im(region)
    lmr=lm(region)

    ! Two scenarios
    if (npe_lon /= 1) then
       !
       ! decomposition along longitudes => work on row_root
       !
       if (nred(region)/=0) then  ! some reduce grid on this tile

          halo=m_dat(region)%halo
          
          ! local array to gather (remove i-halo)
          allocate( m_alt( i1:i2, j1:j2, lmr))
          allocate(rm_alt( i1:i2, j1:j2, lmr, ntracet))

          m_alt = m(i1:i2,j1:j2,:)
          CALL GATHER_J_BAND(dgrid(region), m_alt, m_alt2, status, jref=j1, rowcom=.true.)
          IF_NOTOK_RETURN(status=1)
          if (isRowRoot) then
             m_alt1( 1:imr,:,:) = m_alt2
             m_alt1(     0,:,:) = m_alt1(  imr,:,:)
             m_alt1(    -1,:,:) = m_alt1(imr-1,:,:)
             !PLS m_alt1( imr+1,:,:) = m_alt1(1,:,:)
             !PLS m_alt1( imr+2,:,:) = m_alt1(2,:,:)
          endif
          
          rm_alt = rm(i1:i2,j1:j2,:,:)
          CALL GATHER_J_BAND(dgrid(region), rm_alt, rm_alt2, status, jref=j1, rowcom=.true.)
          IF_NOTOK_RETURN(status=1)
          if (isRowRoot) then
             rm_alt1( 1:imr,:,:,:) = rm_alt2
             rm_alt1(     0,:,:,:) = rm_alt1(  imr,:,:,:)
             rm_alt1(    -1,:,:,:) = rm_alt1(imr-1,:,:,:)
             !PLS rm_alt1( imr+1,:,:,:) = rm_alt1(    1,:,:,:)
             !PLS rm_alt1( imr+2,:,:,:) = rm_alt1(    2,:,:,:)
          endif

#ifdef slopes
          rm_alt = rxm(i1:i2,j1:j2,:,:)
          CALL GATHER_J_BAND(dgrid(region), rm_alt, rm_alt2, status, jref=j1, rowcom=.true.)
          IF_NOTOK_RETURN(status=1)
          if (isRowRoot) then
             rxm_alt1( 1:imr,:,:,:) = rm_alt2
             rxm_alt1(     0,:,:,:) = rm_alt1(  imr,:,:,:)
             rxm_alt1(    -1,:,:,:) = rm_alt1(imr-1,:,:,:)
             !PLS rxm_alt1( imr+1,:,:,:) = rm_alt1(    1,:,:,:)
             !PLS rxm_alt1( imr+2,:,:,:) = rm_alt1(    2,:,:,:)
          endif
          rm_alt = rym(i1:i2,j1:j2,:,:)
          CALL GATHER_J_BAND(dgrid(region), rm_alt, rm_alt2, status, jref=j1, rowcom=.true.)
          IF_NOTOK_RETURN(status=1)
          if (isRowRoot) then
             rym_alt1( 1:imr,:,:,:) = rm_alt2
             rym_alt1(     0,:,:,:) = rm_alt1(  imr,:,:,:)
             rym_alt1(    -1,:,:,:) = rm_alt1(imr-1,:,:,:)
             !PLS rym_alt1( imr+1,:,:,:) = rm_alt1(    1,:,:,:)
             !PLS rym_alt1( imr+2,:,:,:) = rm_alt1(    2,:,:,:)
          endif
          rm_alt = rzm(i1:i2,j1:j2,:,:)
          CALL GATHER_J_BAND(dgrid(region), rm_alt, rm_alt2, status, jref=j1, rowcom=.true.)
          IF_NOTOK_RETURN(status=1)
          if (isRowRoot) then
             rzm_alt1( 1:imr,:,:,:) = rm_alt2
             rzm_alt1(     0,:,:,:) = rm_alt1(  imr,:,:,:)
             rzm_alt1(    -1,:,:,:) = rm_alt1(imr-1,:,:,:)
             !PLS rzm_alt1( imr+1,:,:,:) = rm_alt1(    1,:,:,:)
             !PLS rzm_alt1( imr+2,:,:,:) = rm_alt1(    2,:,:,:)
          endif
#endif
          
          if (isRowRoot) then

             do lrg=1,nred(region)

                j = jred(lrg,region)
                redfact=clustsize(j,region)

                do l=1,lmr
                   ! air mass
                   do i = 1,imredj(j,region)
                      ! the is:ie  array section will be reduced to i
                      is = (i-1)*redfact + 1  
                      ie = i*redfact
                      
                      summ = sum(m_alt1(is:ie,j,l))
                      m_alt1(i,j,l) = summ
                   end do
                   ! tracers
                   do n=1, ntracet
                      do i =  1,imredj(j,region)
                         is = (i-1)*redfact + 1  
                         ie = i*redfact
                         sumrm = sum(rm_alt1(is:ie,j,l,n))
                         rm_alt1(i,j,l,n)  = sumrm
#ifdef slopes
                         sumrm = sum(rxm_alt1(is:ie,j,l,n))
                         rxm_alt1(i,j,l,n) = sumrm
                         sumrm = sum(rym_alt1(is:ie,j,l,n))
                         rym_alt1(i,j,l,n) = sumrm
                         sumrm = sum(rzm_alt1(is:ie,j,l,n))
                         rzm_alt1(i,j,l,n) = sumrm
#endif
                      end do
                      ! JFM: set remaining masses to zero
                      !      for consistency with adjoint
                      do i = imredj(j,region)+1, im(region)
                         rm_alt1(i,j,l,n) = 0.
#ifdef slopes
                         rxm_alt1(i,j,l,n) = 0.
                         rym_alt1(i,j,l,n) = 0.
                         rzm_alt1(i,j,l,n) = 0.
#endif
                      end do
                   end do

                end do  !l
             enddo
          endif
          deallocate(m_alt, rm_alt)
       endif

    else
       !
       ! Reduced grid without longitudinal decomposition
       !
       do lrg=1,nred(region)

          j = jred(lrg,region)
          redfact=clustsize(j,region)

          do l=1,lmr
             ! air mass
             do i =  1,imredj(j,region)
                ! the is:ie  array section will be reduced to i
                is = (i-1)*redfact + 1  
                ie = i*redfact
                
                summ = sum(m(is:ie,j,l))
                m(i,j,l) = summ
             end do
             !cmkm_uni(is:ie,j,l) = m_uni(is:ie,j,l)/summ   
             !   use as distribution function
             ! when transferring back from reduced--->uniform grid...
             ! these summations mean that mixing ratio and the 
             ! the slopes are averaged out within the is:ie section
             ! with m(is:ie,...) taken as the weights

             ! tracers
             do n=1, ntracet
                do i =  1,imredj(j,region)
                   is = (i-1)*redfact + 1  
                   ie = i*redfact
                   sumrm = sum(rm(is:ie,j,l,n))
                   rm(i,j,l,n)  = sumrm
#ifdef slopes
                   sumrm = sum(rxm(is:ie,j,l,n))
                   rxm(i,j,l,n) = sumrm
                   sumrm = sum(rym(is:ie,j,l,n))
                   rym(i,j,l,n) = sumrm
                   sumrm = sum(rzm(is:ie,j,l,n))
                   rzm(i,j,l,n) = sumrm
#endif
                end do
                ! JFM: set remaining masses to zero
                !      for consistency with adjoint
                do i = imredj(j,region)+1, im(region)
                   rm(i,j,l,n) = 0.
#ifdef slopes
                   rxm(i,j,l,n) = 0.
                   rym(i,j,l,n) = 0.
                   rzm(i,j,l,n) = 0.
#endif
                end do
             end do

             !put periodic boundary...
          end do  !l

       end do   !redgrid...

    endif
 
    nullify(m)
    nullify(rm)
#ifdef slopes
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
#endif

    status=0
    
  end subroutine uni2red
  !EOC


#ifdef secmom
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  uni2red_2nd
  !
  ! !DESCRIPTION: Transforms data from uniform grid to reduced grid. @nd
  !               moments version.
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine uni2red_2nd(region)
    !
    ! !USES:
    !
    use dims
    use meteodata,   only : m_dat
    use global_data, only : mass_dat
    use chem_param,  only : ntracet
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in)  :: region
    !
    ! !REVISION HISTORY: 
    !   written by mike botchev, march-june 1999
    !   modified by Maarten Krol, dec 2002
    !   20 Dec 2012 - Ph Le Sager -
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real, dimension(:,:,:,:), pointer  :: rm, rxm,rym,rzm, rxxm, rxym, rxzm, ryym, ryzm, rzzm
    real, dimension(:,:,:),   pointer  :: m

    integer i,ie,is,j,l,lrg,redfact,n,lmr
    real summ,sumrm

    ! start
    m    => m_dat(region)%data
    rm   => mass_dat(region)%rm
    rxm  => mass_dat(region)%rxm
    rym  => mass_dat(region)%rym
    rzm  => mass_dat(region)%rzm
    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

    lmr=lm(region)

    do lrg=1,nred(region)
       j = jred(lrg,region)
       redfact=clustsize(j,region)
       do l=1,lmr 
          do i =  1,imredj(j,region)
             ! the is:ie  array section will be reduced to i
             is = (i-1)*redfact + 1  
             ie = i*redfact
             summ = sum(m(is:ie,j,l))
             m(i,j,l) = summ
             !cmkm_uni(is:ie,j,l) = m_uni(is:ie,j,l)/summ   
             !   use as distribution function
             ! when transferring back from reduced--->uniform grid...
             ! these summations mean that mixing ratio and the 
             ! the slopes are averaged out within the is:ie section
             ! with m(is:ie,...) taken as the weights
             do n=1,ntracet
                !#endif
                sumrm = sum(rm(is:ie,j,l,n))
                rm(i,j,l,n)  = sumrm
                sumrm = sum(rxm(is:ie,j,l,n))
                rxm(i,j,l,n) = sumrm
                sumrm = sum(rym(is:ie,j,l,n))
                rym(i,j,l,n) = sumrm
                sumrm = sum(rzm(is:ie,j,l,n))
                rzm(i,j,l,n) = sumrm
                sumrm = sum(rxxm(is:ie,j,l,n))
                rxxm(i,j,l,n) = sumrm
                sumrm = sum(rxym(is:ie,j,l,n))
                rxym(i,j,l,n) = sumrm
                sumrm = sum(rxzm(is:ie,j,l,n))
                rxzm(i,j,l,n) = sumrm
                sumrm = sum(ryym(is:ie,j,l,n))
                ryym(i,j,l,n) = sumrm
                sumrm = sum(ryzm(is:ie,j,l,n))
                ryzm(i,j,l,n) = sumrm
                sumrm = sum(rzzm(is:ie,j,l,n))
                rzzm(i,j,l,n) = sumrm
             end do   !n
          end do  !i
          !put periodic boundary...
       end do  !l

    end do   !redgrid...

    nullify(m)
    nullify(rm)
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
    nullify(rxxm)
    nullify(rxym)
    nullify(rxzm)
    nullify(ryym)
    nullify(ryzm)
    nullify(rzzm)

  end subroutine uni2red_2nd
  !EOC

#endif


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RED2UNI
  !
  ! !DESCRIPTION: transforms data from reduced grid back to uniform grid
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine red2uni(region, i1, i2, j1, j2, halo, m_uni, status)
    !
    ! !USES:
    !
    use dims,         only : im, lm
    use tracer_data,  only : mass_dat
    use meteodata,    only : m_dat
    use chem_param,   only : ntracet
    use partools,     only : isRowRoot, npe_lon
    use tm5_distgrid, only : dgrid, scatter_j_band, gather_j_band, update_halo_jband
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: region, i1, i2, j1, j2, halo
    real,    intent(in) :: m_uni(i1-halo:,j1-halo:,1:)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !   written by mike botchev, march-june 1999
    !   20 Dec 2012 - Ph Le Sager - TM5-MP version
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    ! local
    real,dimension(:,:,:,:),pointer    :: rm, rm_
#ifdef slopes
    real,dimension(:,:,:,:),pointer    :: rxm,rym,rzm
#endif
    real,dimension(:,:,:)  ,pointer    :: m, m_, m_alt
    integer :: i, ie, is, j, l, lrg, n, redfact, imr, lmr, nt, halo_r
    real    :: hi, mass, mass_coord, rmm, slope, m_old
    character(len=5) :: distr_mode
        
    ! start
    nt=ntracet
    imr=im(region)
    lmr=lm(region)

    m   => m_dat(region)%data
    rm  => mass_dat(region)%rm
#ifdef slopes
    rxm => mass_dat(region)%rxm
    rym => mass_dat(region)%rym
    rzm => mass_dat(region)%rzm
#endif

    ! Two scenarios
    if (npe_lon /= 1) then
       !
       ! decomposition along longitudes => work on row_root
       !
       if (nred(region)/=0) then
          
          ! gather m_uni on row_root, where we redistribute the masses
          allocate(m_(i1:i2,j1:j2,1:lmr))
          m_ = m_uni(i1:i2,j1:j2,1:lmr)
                    
          CALL GATHER_J_BAND(dgrid(region), m_, m_alt2, status, jref=j1, rowcom=.true.)
          IF_NOTOK_RETURN(status=1)
                    
          if (isRowRoot) then
             
             do lrg=1,nred(region)
                j = jred(lrg,region)
                redfact=clustsize(j,region)
                !m(1:imr,j,:)= m_alt2(:,j,:)

                do l=1,lmr
                   do i=imredj(j,region),1,-1

                      is = (i-1)*redfact + 1  
                      ie = i*redfact
                                            
                      mass=m_alt1(i,j,l)
                      
                      do n=1,nt
                         rmm = rm_alt1(i,j,l,n) 
                         rm_alt1(is:ie,j,l,n)= m_alt2(is:ie,j,l)/mass* rmm
#ifdef slopes
                         rmm = rxm_alt1(i,j,l,n) 
                         rxm_alt1(is:ie,j,l,n)= m_alt2(is:ie,j,l)/mass * rmm
                         ! rym and rzm are always distributed uniformly:
                         rmm = rym_alt1(i,j,l,n) 
                         rym_alt1(is:ie,j,l,n)= m_alt2(is:ie,j,l)/mass * rmm
                         rmm = rzm_alt1(i,j,l,n) 
                         rzm_alt1(is:ie,j,l,n)= m_alt2(is:ie,j,l)/mass * rmm
#endif
                      end do
                      
                   end do
                enddo
             enddo
          endif
                   
          ! scatter results

          if (nred(region)< (j2-j1+1)) then
             
             if (isRowRoot) m_alt2 = m_alt1(1:imr,:,:)

             CALL SCATTER_J_BAND(dgrid(region), m_, m_alt2, status, jref=j1, rowcom=.true.)
             IF_NOTOK_RETURN(status=1)
             
             m(i1:i2,j1:j2,:) = m_
          endif
          
          allocate(rm_( i1:i2, j1:j2, lmr, ntracet))
          
          if (isRowRoot) rm_alt2 = rm_alt1(1:imr,:,:,:)

          CALL SCATTER_J_BAND(dgrid(region), rm_, rm_alt2, status, jref=j1, rowcom=.true.)
          IF_NOTOK_RETURN(status=1)

          rm(i1:i2,j1:j2,:,:) = rm_
          
#ifdef slopes
          if (isRowRoot) rm_alt2 = rxm_alt1(1:imr,:,:,:)
          CALL SCATTER_J_BAND(dgrid(region), rm_, rm_alt2, status, jref=j1, rowcom=.true.)
          IF_NOTOK_RETURN(status=1)
          rxm(i1:i2,j1:j2,:,:) = rm_

          if (isRowRoot) rm_alt2 = rym_alt1(1:imr,:,:,:)
          CALL SCATTER_J_BAND(dgrid(region), rm_, rm_alt2, status, jref=j1, rowcom=.true.)
          IF_NOTOK_RETURN(status=1)
          rym(i1:i2,j1:j2,:,:) = rm_

          if (isRowRoot) rm_alt2 = rzm_alt1(1:imr,:,:,:)
          CALL SCATTER_J_BAND(dgrid(region), rm_, rm_alt2, status, jref=j1, rowcom=.true.)
          IF_NOTOK_RETURN(status=1)
          rzm(i1:i2,j1:j2,:,:) = rm_
#endif

          ! (2) bands with reduced grid, just use m_uni
          do lrg=1,nred(region)
             j = jred(lrg,region)
             m(:,j,:)= m_uni(:,j,:)
          end do
          
          ! update halo
          halo_r = mass_dat(region)%halo
          call UPDATE_HALO_JBAND(dgrid(region), m, m_dat(region)%halo, status )
          call UPDATE_HALO_JBAND(dgrid(region), rm(:,:,:,1:nt), halo_r, status )
#ifdef slopes          
          call UPDATE_HALO_JBAND(dgrid(region), rxm, halo_r, status )
          call UPDATE_HALO_JBAND(dgrid(region), rym, halo_r, status )
          call UPDATE_HALO_JBAND(dgrid(region), rzm, halo_r, status )
#endif
          deallocate(m_,rm_)

       endif
       
    else
       !
       ! Reduced grid without longitudinal decomposition
       !
       do lrg=1,nred(region)
          j = jred(lrg,region)
          redfact=clustsize(j,region)
          
          do l=1,lmr 
             do i =  imredj(j,region),1,-1
                ! the i cell will be distributed within the is:ie array section
                is = (i-1)*redfact + 1  
                ie = i*redfact

                !m_uni is the mass-distribution in the non-reduced grid/divided by
                !the reduced_grid mass. This is used as distribution function!....
                mass=m(i,j,l)
                m(is:ie,j,l)= m_uni(is:ie,j,l)
                
                do n=1,nt
                   rmm = rm(i,j,l,n) 
                   rm(is:ie,j,l,n)= m_uni(is:ie,j,l)/mass* rmm
#ifdef slopes
                   rmm = rxm(i,j,l,n) 
                   rxm(is:ie,j,l,n)= m_uni(is:ie,j,l)/mass * rmm
                   ! rym and rzm are always distributed uniformly:
                   rmm = rym(i,j,l,n) 
                   rym(is:ie,j,l,n)= m_uni(is:ie,j,l)/mass * rmm
                   rmm = rzm(i,j,l,n) 
                   rzm(is:ie,j,l,n)= m_uni(is:ie,j,l)/mass * rmm
#endif
                end do
             end do

             ! update cell(0,...) according to the periodic bc's:
             rm(0,j,l,1:nt) = rm(im(region),j,l,1:nt)
#ifdef slopes
             rxm(0,j,l,:) = rxm(im(region),j,l,:)
             rym(0,j,l,:) = rym(im(region),j,l,:)
             rzm(0,j,l,:) = rzm(im(region),j,l,:)
#endif
             m(0,j,l) = m(im(region),j,l)

          end do
       end do

    endif

    nullify(m)
    nullify(rm)
#ifdef slopes
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
#endif

  end subroutine red2uni
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RED2UNI_EM
  !
  ! !DESCRIPTION: transforms data from reduced grid back to uniform grid,
  !               using EQUAL MASSES instead of reduced mass array (m_uni).
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine red2uni_em(region)
    ! 
    !USES:
    !
    use dims
    use meteodata  , only : m_dat
    use tracer_data, only : mass_dat
    use chem_param,  only : ntracet
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in) :: region
    !
    ! !REVISION HISTORY: 
    !  written by mike botchev, march-june 1999
    !  20 Dec 2012 - Ph Le Sager - TM5-MP version
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real,dimension(:,:,:,:),pointer    :: rm
#ifdef slopes
    real,dimension(:,:,:,:),pointer    :: rxm,rym,rzm
#endif
    real,dimension(:,:,:)  ,pointer    :: m
    integer i,ie,ii,is,j,l,lrg,n,redfact,lmr
    real hi,mass,mass_coord,rmm,slope,m_old
    character(len=5) distr_mode

    !------------ start

    lmr=lm(region)

    m   => m_dat(region)%data
    rm  => mass_dat(region)%rm
#ifdef slopes
    rxm => mass_dat(region)%rxm
    rym => mass_dat(region)%rym
    rzm => mass_dat(region)%rzm
#endif

    distr_mode = 'unfrm' ! 'slope' or 'unfrm' 

    do lrg=1,nred(region)
       j = jred(lrg,region)
       redfact=clustsize(j,region)
       do l=1,lmr 
          do i =  imredj(j,region),1,-1
             ! the i cell will be distributed within the is:ie array section
             is = (i-1)*redfact + 1  
             ie = i*redfact

             !m_uni is the mass-distribution in the non-reduced grid/divided by
             !the reduced_grid mass. This is used as distribution function!....
             mass=m(i,j,l); m(is:ie,j,l)= mass/(ie-is+1)

             if (distr_mode=='unfrm') then
                ! mixing ratio and x-slope will be UNiFoRMly distributed 
                ! within is:ie
                do n=1,ntracet
                   rmm = rm(i,j,l,n) 
                   rm(is:ie,j,l,n)= m(is:ie,j,l)/mass* rmm
#ifdef slopes
                   rmm = rxm(i,j,l,n) 
                   rxm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm
                   ! rym and rzm are always distributed uniformly:
                   rmm = rym(i,j,l,n) 
                   rym(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm
                   rmm = rzm(i,j,l,n) 
                   rzm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm
#endif
                end do
             end if
             !cmkelseif(distr_mode=='slope') then
          end do
          ! update cell(0,...) according to the periodic bc's:
          rm(0,j,l,1:ntracet) = rm(im(region),j,l,1:ntracet)
#ifdef slopes
          rxm(0,j,l,:) = rxm(im(region),j,l,:)
          rym(0,j,l,:) = rym(im(region),j,l,:)
          rzm(0,j,l,:) = rzm(im(region),j,l,:)
#endif
          m(0,j,l) = m(im(region),j,l)
       end do
    end do

    nullify(m)
    nullify(rm)
#ifdef slopes
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
#endif

  end subroutine red2uni_em
  !EOC

#ifdef secmom
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  red2uni_em_2nd
  !
  ! !DESCRIPTION: transforms data from reduced grid back to uniform grid
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine red2uni_em_2nd(region)
    !
    ! !USES:
    !
    use dims
    use meteodata,   only : m_dat
    use global_data, only : mass_dat
    use chem_param,  only : ntracet
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in) :: region
    !
    ! !REVISION HISTORY: 
    !   written by mike botchev, march-june 1999
    !   20 Dec 2012 - Ph Le Sager - TM5-MP version
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real,dimension(:,:,:,:),pointer    :: rm,rxm,rym,rzm, rxxm, rxym, rxzm, ryym, ryzm, rzzm
    real,dimension(:,:,:)  ,pointer    :: m
    integer i,ie,ii,is,j,l,lrg,n,redfact,lmr
    real hi,mass,mass_coord,rmm,slope,m_old
    character(len=5) distr_mode

    ! ----- start

    lmr=lm(region)

    m    => m_dat(region)%data
    rm   => mass_dat(region)%rm
    rxm  => mass_dat(region)%rxm
    rym  => mass_dat(region)%rym
    rzm  => mass_dat(region)%rzm
    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

    distr_mode = 'unfrm' ! 'slope' or 'unfrm' 

    do lrg=1,nred(region)
       j = jred(lrg,region)
       redfact=clustsize(j,region)
       do l=1,lmr 
          do i =  imredj(j,region),1,-1
             ! the i cell will be distributed within the is:ie array section
             is = (i-1)*redfact + 1  
             ie = i*redfact

             !m_uni is the mass-distribution in the non-reduced grid/divided by
             !the reduced_grid mass. This is used as distribution function!....
             mass=m(i,j,l); m(is:ie,j,l)= mass/(ie-is+1)

             if (distr_mode=='unfrm') then
                ! mixing ratio and x-slope will be UNiFoRMly distributed 
                ! within is:ie
                !#ifdef MPI                  
                !                do n=1,ntracetloc
                !#else
                do n=1,ntracet
                   !#endif
                   rmm = rm(i,j,l,n) 
                   rm(is:ie,j,l,n)= m(is:ie,j,l)/mass* rmm
                   rmm = rxm(i,j,l,n) 
                   rxm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm
                   ! rym and rzm are always distributed uniformly:
                   rmm = rym(i,j,l,n) 
                   rym(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm
                   rmm = rzm(i,j,l,n) 
                   rzm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm
                   rmm = rxxm(i,j,l,n)
                   rxxm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm
                   rmm = rxym(i,j,l,n)
                   rxym(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm
                   rmm = rxzm(i,j,l,n)
                   rxzm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm
                   rmm = ryym(i,j,l,n)
                   ryym(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm
                   rmm = ryzm(i,j,l,n)
                   ryzm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm
                   rmm = rzzm(i,j,l,n)
                   rzzm(is:ie,j,l,n)= m(is:ie,j,l)/mass * rmm
                end do
             end if
             !cmkelseif(distr_mode=='slope') then
          end do
          ! update cell(0,...) according to the periodic bc's:
          rm(0,j,l,1:ntracet) = rm(im(region),j,l,1:ntracet)
          rxm(0,j,l,:) = rxm(im(region),j,l,:)
          rym(0,j,l,:) = rym(im(region),j,l,:)
          rzm(0,j,l,:) = rzm(im(region),j,l,:)
          m(0,j,l) = m(im(region),j,l)

       end do
    end do

    nullify(m)
    nullify(rm)
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
    nullify(rxxm)
    nullify(rxym)
    nullify(rxzm)
    nullify(ryym)
    nullify(ryzm)
    nullify(rzzm)

  end subroutine red2uni_em_2nd
  !EOC

#endif

end module redgridZoom