!### 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"
!
!#################################################################

module advectx

  use GO, only : gol, goPr, goErr
  
  implicit none

  ! --- interface

  private

  public :: advectxzoom, put_xedges


contains

  subroutine advectxzoom( region )

    ! set parameters for advectx 
    ! written by patrick berkvens and mike botchev, march-june 1999
    ! updated and modified by MK, dec 2002

    use dims,          only : xref, yref, zref, tref, im, jm, lm
    use dims,          only : status, nsplitsteps, n_operators, zoom2D
    use dims,          only : touch_sp, touch_np, okdebug
    use dims,          only : parent, splitorderzoom, zero
#ifdef with_budgets    
    use budget_global, only : sum_advection
#endif
    use global_data  , only : mass_dat, wind_dat
    use toolbox,       only : escape_tm
#ifdef MPI 
    use mpi_const,     only : my_real, com_lev, ierr
    use mpi_const,     only : root_t, mpi_sum, pe_first_tracer
    use mpi_comm,      only : barrier_t
#endif
    use partools     , only : myid, lmloc, ntracetloc

    implicit none

    ! input/output
    integer,intent(in)                  :: region

    ! local
    real,dimension(:,:,:),  pointer     :: am

    integer                             :: js,je,ls,le,n,q
    integer                             :: imr,jmr,lmr,tref_,xref_,yref_,zref_
    real,dimension(:,:),allocatable     :: am0,am1  ! for slopes only
    logical                             :: x_encountered
    character(len=1)                    :: dir
    real                                :: sum_old, sum_new
    real                                :: sum_old_all, sum_new_all

    ! start

    if ( ntracetloc == 0 ) return !WP! only PE's with nonzero lmloc proceed

    call put_xedges(region)   ! write BC to cildren

    tref_ = tref(region)/tref(parent(region))
    xref_ = xref(region)/xref(parent(region))
    yref_ = yref(region)/yref(parent(region))
    zref_ = zref(region)/zref(parent(region))
    imr = im(region)
    jmr = jm(region)
    lmr = lm(region)

    allocate(am0(0:jm(region)+1, 0:lmr+1))
    allocate(am1(0:jm(region)+1, 0:lmr+1))

    am => wind_dat(region)%am_t

    sum_old = sum(mass_dat(region)%rm_t(1:im(region),1:jm(region),1:lmr,1))

    ! determine the scope for advectx: 

    if ( region == 1 ) then
       yref_ = 0 
       zref_ = 0  ! to have js/je and ls/le properly computed
    end if
    q=status(region)/((nsplitsteps/2)*tref_)  
    ! find q - the place in the splitorderzoom 
    ! corresponding to the begining of the last 
    ! triple x-y-z of the parent
    x_encountered=.false.

    do n=1,n_operators    
       ! now track four following steps from q !wp! changed to n_op.. steps

       dir=splitorderzoom(region,q*(nsplitsteps/2)*tref_+n)
       select case(dir)
       case('x') 
          x_encountered=.true.
       case('y') 
          if ( .not. x_encountered) then  
             js=1                             ! y-substep is before x => 
             je=jm(region)                    ! full j-scope 
          else
             if( touch_sp(region) == 1 ) then
                js=1                             !CMK special case...touching SP
             else
                js=yref_+1                       ! x-substep is before y =>
             end if
             if( touch_np(region) == 1 ) then
                je=jm(region)                    ! CMK special case: touching NP
             else
                je=jm(region)-yref_              ! restricted j-scope
             end if
          end if
       case('z') 
          if  ( .not. x_encountered ) then
             ls=1                             ! z-substep is before x => 
             le=lmr                    ! full l-scope 
          else
             ls=zref_+1                       ! x-substep is before z =>
             le=lmr-zref_              ! restricted l-scope
          end if
          if (zoom2D) then
             ls=1; le=lmr    !WP! this is always the case
          end if
       case ('c') 
       case ('v')
       case ('s') 
       case default
          print *,'advectxzoom: strange value in splitorderzoom(',region,',',   &
               q*(nsplitsteps/2)*tref_+n,'): ',q
          call escape_tm( 'advectxzoom: error in x-advection ')
       end select

    end do  ! n=1,n_operators

#ifdef slopes
#ifndef secmom
    if ( ( mod(status(region),(nsplitsteps/2)*tref_) >= n_operators ) .and. &
         ( im(region)/xref(region) < im(1) ) ) then
       ! IF (1) more than fourth substep, (2) slope scheme, and 
       ! (3) the region is not [0;360] degrees wide THEN
       ! at this substep no coarse fluxes will be applied 

       ! for slopes only: zero out velocities via the x-edges of the region
       ! to assure that no fluxes will be applied

       if ( okdebug ) print *,'advectxzoom: zeroing out outer mass fluxes'
       am0 = am(    xref_-1,:,:)
       am(    xref_-1,:,:) = zero
       am1 = am(imr-xref_+1,:,:)
       am(imr-xref_+1,:,:) = zero

    end if
#endif
#endif

    call advectx_work(region,js,je,ls,le)

#ifdef slopes
#ifndef secmom
    if ( ( mod(status(region),(nsplitsteps/2)*tref_) >= n_operators ).and. &
         ( im(region)/xref(region) < im(1) ) ) then

       ! for slopes only: recreate velocities via the x-edges
       am(    xref_-1,:,:) = am0
       am(imr-xref_+1,:,:) = am1

    end if
#endif
#endif

    sum_new = sum(mass_dat(region)%rm_t(1:im(region),1:jm(region),1:lmr,1))

#ifdef with_budgets
#ifdef MPI
    if( myid == pe_first_tracer ) &
         sum_advection(region) = sum_advection(region) + sum_new - sum_old
    if( okdebug .and. myid == root_t ) &
         print *, 'advectxzoom: region, sum_advection ',&
         region,sum_advection(region),sum_new,sum_old
#else
    sum_advection(region) = sum_advection(region) + sum_new - sum_old
    if ( okdebug ) print *, 'advectxzoom: region, sum_advection ', &
         region,sum_advection(region),sum_new,sum_old
#endif
#endif

    deallocate(am0)
    deallocate(am1)
    nullify(am)

  end subroutine advectxzoom



  subroutine advectx_work(region,js,je,ls,le)
    !
    ! makes reduced grid pre-/postprocessing and switches between 
    ! dynamu and dynamu1
    ! written by mike botchev, march-june 1999
    !
    use dims,        only : im, jm, lm
    use redgridZoom, only : grid_reduced, nred, uni2red, uni2red_mf, red2uni
    use global_data, only : wind_dat
    use global_data, only : mass_dat
    use meteodata  , only : m_dat
    use toolbox,     only : escape_tm
    use partools   , only : lmloc
    use chem_param,  only : ntracet

    ! input
    integer,intent(in) :: region
    integer,intent(in) :: js
    integer,intent(in) :: je
    integer,intent(in) :: ls
    integer,intent(in) :: le

    ! local
    real,dimension(:,:,:),pointer     :: m
    real,dimension(:,:,:),pointer     :: am
    real,dimension(:,:,:),allocatable :: m_uni   ! for reduced grid...
    real,dimension(:,:,:),allocatable :: am_uni  ! for reduced grid...
    real                 :: mxval
    integer,dimension(3) :: mxloc
    integer              :: imr,jmr,lmr

    ! start

#ifndef slopes
    call escape_tm( 'advectx: slopes advection not defined')
#endif

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

    allocate(m_uni(-1:im(region)+2,-1:jm(region)+2, lmr))
    allocate(am_uni(0:im(region)+1, 0:jm(region)+1, 0:lmr+1))

    am => wind_dat(region)%am_t
    m  => m_dat(region)%data

    ! transform to the reduced grid:
    ! check for reduced grid in region

    if ( grid_reduced .and. (nred(region) /= 0) ) then 
       ! save non-reduced m and am in m_uni and am_uni:
       m_uni = m
       am_uni = am
       ! reduce m,rm,rxm,rym,rzm:
       call uni2red(region)
       ! reduce am:
       call uni2red_mf(region)
    end if

    call dynamu (region,js,je,ls,le)

    ! transform from the reduced grid:

    if ( grid_reduced .and. nred(region) /= 0 ) then
       ! advection on uniform grid:
       m_uni(1:imr,1:jmr,1:lmr) = m_uni(1:imr,1:jmr,1:lmr) + &
            am_uni(0:imr-1,1:jmr,1:lmr) - am_uni(1:imr,1:jmr,1:lmr)

       ! redistribute rm,rxm,rym,rzm by using m_uni and m
       call red2uni(region,m_uni)
       ! recreate rm/rxm/rym/rzm by using equal masses (not m_uni) 
       !call red2uni_em(region)
       ! recreate am:
       am = am_uni
    end if

    nullify(am)
    nullify(m)

    deallocate(am_uni)
    deallocate(m_uni)

  end subroutine advectx_work

  subroutine dynamu(region,js,je,ls,le)
    !-----------------------------------------------------------------------
    !
    !****   dynamu  - east-west tracer transport  - v 9.1
    !
    !   programmed by           mh      mpi HH          23-feb-1995
    !
    !   purpose
    !   -------
    !   calculate amount of tracer moved in an east-west advection
    !   substep
    !
    !   interface
    !   ---------
    !   call dynamu
    !
    !   method
    !   ------
    !   slopes scheme
    !
    !   externals
    !   ---------
    !   none
    !
    !   reference
    !   ---------       
    !   Russel and Lerner, 1979
    !-----------------------------------------------------------------------
    !   fixed bug in calculating maximum courant number
    !                          mh, Thu, Feb 24, 1994 12:49:27
    !     included code for limits of slopes to prevent negative tracer 
    !     masses                   mh, 20-jun-1994
    !     zoom version written by mike botchev, march-june 1999
    !-----------------------------------------------------------------------

    use dims,        only : okdebug, nregions, parent
    use dims,        only : im, jm, lm, xref, yref, zref, tref
    use dims,        only : zero, one, xi, nxi, nloop_max, limits, xcyc
    use redgridZoom, only : grid_reduced, nred, imred
    use global_data, only : wind_dat
    use tracer_data, only : mass_dat
    use meteodata  , only : m_dat
    use zoom_tools,  only : mix_edges
    use toolbox,     only : escape_tm
#ifdef with_budgets    
    use budget_global, only : budget_flux, iflux1, iflux2, apply_budget_global_flux
#endif
#ifdef MPI
    use mpi_const,   only : ntracetloc,myid,root_t,com_trac
    use mpi_const,   only : ierr,my_real,mpi_max,mpi_min, ntracet_ar
    use mpi_comm,    only : barrier_t
#endif
    use chem_param,  only : ntracet, ra
    use omp_partools,only : my_omp_domain, my_omp_max, my_omp_sum

    implicit none

    ! input/output
    integer,intent(in) :: region
    integer,intent(in) :: js
    integer,intent(in) :: je
    integer,intent(in) :: ls
    integer,intent(in) :: le

    ! local
    real,dimension(:,:,:,:),pointer    :: rm,rxm,rym,rzm
    real,dimension(:,:,:),  pointer    :: m,am
    real,dimension(:),allocatable              :: f,pf,fy,fz,mnew
    integer i,ie,is,j,l,n, nglob, offsetn
    integer imr,jmr,lmr,tref_,xref_,yref_,zref_,ic,nloop, iloop
    real :: min_one,max_one,min_all,max_all
    real :: alpha, max_alpha, x, rmold
    real,dimension(:), allocatable :: mx
    logical                        :: cfl_ok
    integer, parameter             :: max_nloop = 50

    ! Openmp parameters
    integer            :: ns, ne, jss, jee, iie, iimr
    real               :: l_xi, g_xi, l_nxi, g_nxi, l_nl, g_nl
    logical            :: special_grid


    ! start

    if ( okdebug ) print *,'dynamu: region=',region, &
         ' js,je,ls,le=',js,je,ls,le
    if ( (region < 0) .or. (region > nregions) ) &
         call escape_tm( 'dynamu: STOP, illegal number of region !!!')

    ! prepare x-edges to allow for uniform work:
    
    call compress_xedges(region,js,je,ls,le)

    ! compute refinement factors with respect to the parent
    xref_ = xref(region)/xref(parent(region))
    yref_ = yref(region)/yref(parent(region))
    zref_ = zref(region)/zref(parent(region))
    tref_ = tref(region)/tref(parent(region))
    imr=im(region);jmr=jm(region);lmr=lm(region)

    am => wind_dat(region)%am_t
    m => m_dat(region)%data
    rm => mass_dat(region)%rm_t
    rxm => mass_dat(region)%rxm_t
    rym => mass_dat(region)%rym_t
    rzm => mass_dat(region)%rzm_t

    ! check js,je,ls,le:
    if((js/=yref_+1)  .and.(js/=1))   stop 'dynamu: Wrong value for JS '
    if((je/=jmr-yref_).and.(je/=jmr)) stop 'dynamu: Wrong value for JE '
    if((ls/=zref_+1)  .and.(ls/=1))   stop 'dynamu: Wrong value for LS '
    if((le/=lmr-zref_).and.(le/=lmr)) stop 'dynamu: Wrong value for LE '

    ! compute is/ie -- cells is:ie,js:je:ls:le will be updated
    is = xref_; ie = imr-xref_+1
    
    if (im(region)/xref(region)==im(1)) then
       ! periodic boundary condition
       is = 1; ie = imr
    end if

    special_grid = grid_reduced .and. nred(region) /= 0

    ! if requested limit zonal slopes such that no negative
    ! tracer masses should occur
    if ( limits ) then
!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (special_grid, region, imred, is, ie, js, je, ls, le, &
!$OMP           rm, rxm) &
!$OMP  private (i, j, l, n, ns, ne)
       call my_omp_domain (1, ubound (rxm, 4), ns, ne)
       if (special_grid) then
          do n = ns, ne
             do l = ls, le
                do j = js, je
                   do i = is - 1, imred(j,region) + 1
                      rxm(i,j,l,n) = max (min (rxm(i,j,l,n), &
                           rm(i,j,l,n)), -rm(i,j,l,n))
                   end do
                end do
             end do
          end do
       else
          do n = ns, ne
             do l = ls, le
                do j = js, je
                   do i = is - 1, ie + 1
                      rxm(i,j,l,n) = max (min (rxm(i,j,l,n), &
                           rm(i,j,l,n)), -rm(i,j,l,n))
                   end do
                end do
             end do
          end do
       end if
!$OMP END PARALLEL
    end if

    if ( okdebug ) then

       max_one=maxval(m(1:imr,1:jmr,1:lmr))
       min_one=minval(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)
       if ( myid == root_t ) print *,'dynamu: m_min,m_max: ',min_all,max_all
       max_one=maxval(rm(1:imr,1:jmr,1:lmr,:))
       min_one=minval(rm(1:imr,1:jmr,1:lmr,:))
       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)
       if ( myid == root_t ) print *,'dynamu: rm_min,rm_max: ',min_all,max_all
#else
       print *,'dynamu: m_min,m_max: ',min_one, max_one
       max_one=maxval(rm(1:imr,1:jmr,1:lmr,:))
       min_one=minval(rm(1:imr,1:jmr,1:lmr,:))
       print *,'dynamu: rm_min,rm_max: ',min_one,max_one
#endif

    end if

    ! loop over vertical layers and latitudes

!$OMP PARALLEL &
!$OMP  default (none) &
#if defined (MPI)
!$OMP  shared  (ntracetloc, ntracet_ar, myid) &
#endif
#ifdef with_budgets    
!$OMP  shared  (apply_budget_global_flux, iflux1, iflux2, budget_flux,ra) &
#endif
!$OMP  shared  (region, is, ie, imr, js, je, ls, le, okdebug, limits, &
!$OMP           special_grid, nred, imred, m, am, rm, rxm, rym, rzm, &
!$OMP           g_nl, g_nxi, g_xi) &
!$OMP  private (i, j, l, jss, jee, iie, iimr, nloop, alpha, max_alpha, &
!$OMP           nglob, cfl_ok, mx, f, pf, fy, fz, mnew, l_nl, l_nxi, &
!$OMP           l_xi)

    iie  = ie
    iimr = imr

    l_nl  = 0.0
    l_xi  = 0.0
    l_nxi = 0.0

    allocate (f(0:im(region)))
    allocate (pf(0:im(region)))
    allocate (fy(0:im(region)))
    allocate (fz(0:im(region)))
    allocate (mnew(im(region)))

    call my_omp_domain (js, je, jss, jee)
    do j = jss, jee
       do l=ls,le
          ! reduced grid: less points will be handled
          if (special_grid) iie = imred(j,region)

          ! Determine number of loops required to avoid CFLs...   
          nloop = 1   ! default run the loop one time!
          alpha = 2.0
          allocate(mx(is-1:iie+1))
          max_alpha = 0.0

          !CMK added oct2003: nloop limit..
          do while ( abs(alpha) >= one .and. nloop < max_nloop )
             mx(is-1:iie+1) = m(is-1:iie+1,j,l)   ! copy mass to temp array
             if(xcyc(region)==1) then
                mx(0) = mx(iie)
                mx(iie+1) = mx(1)
             end if

             xloop: do iloop = 1, nloop
                cfl_ok = .true.
                do i=is-1,iie   
                   if (am(i,j,l)>=zero)  then
                      alpha=am(i,j,l)/mx(i)
                   else
                      alpha=am(i,j,l)/mx(i+1)
                   end if
                   if((abs(alpha)>=one)) then                   !PB
                      cfl_ok = .false.
!                      if ( okdebug ) &
!                           print *, 'dynamu: ',i,j,l, ' alpha > 1, since :', &
!                           am(i,j,l),mx(i)
                      exit xloop
                   end if
                   max_alpha = max(max_alpha, abs(alpha))
                end do
                ! update mass
                ! note: only is:iie updated, is-1, iie+1 are BC
                mx(is:iie) = mx(is:iie) + am(is-1:iie-1,j,l) - am(is:iie,j,l)
                if ( xcyc(region) == 1 ) then   !except for periodic BC...
                   mx(0) = mx(iie)
                   mx(iie+1) = mx(1)
                end if
             end do  xloop
             if(.not.cfl_ok) then 
                ! reduce mass flux
                am(is-1:iie,j,l) = am(is-1:iie,j,l)*nloop/(nloop+1)
                if ( okdebug ) &
                     print *, 'dynamu: ',j,l, 'nloop increased to', &
                     nloop + 1, alpha
                nloop = nloop + 1
                max_alpha = 0.0
                if (nloop == max_nloop) call escape_tm('dynuamu: nloop too high')
             end if
          end do !while alpha>1
          deallocate(mx)
          l_nl = MAX (l_nl, REAL(nloop))
          l_xi = MAX (l_xi, max_alpha)
          do iloop = 1,nloop   ! CFL loop
             ! calculate new air mass distribution
             mnew(is:iie)=m(is:iie,j,l) + am(is-1:iie-1,j,l)-am(is:iie,j,l)
             
             if (im(region)/xref(region)==im(1)) then
                m  (0    ,j,l)   = m  (iie,j,l) !periodic BCs...
                m  (iie+1,j,l)    = m  (1,  j,l)
             end if

             ! loop over tracers
#ifdef MPI   
             do n=1,ntracetloc
                nglob = n + SUM(ntracet_ar(0:myid-1))

#else
             do n=1,ntracet
                nglob = n
#endif
                if (im(region)/xref(region)==im(1)) then
                   rm (0,    j,l,n) = rm (iie,j,l,n)   !periodic BCs
                   rxm(0,    j,l,n) = rxm(iie,j,l,n)
                   rym(0,    j,l,n) = rym(iie,j,l,n)
                   rzm(0,    j,l,n) = rzm(iie,j,l,n)
                   rm (iie+1,j,l,n) = rm (1,  j,l,n)
                   rxm(iie+1,j,l,n) = rxm(1,  j,l,n)
                   rym(iie+1,j,l,n) = rym(1,  j,l,n)
                   rzm(iie+1,j,l,n) = rzm(1,  j,l,n)
                end if

                ! calculate fluxes for rm,rxm,rym,rzm
                do i=is-1,iie
                   if (am(i,j,l)>=zero)  then
                      alpha=am(i,j,l)/m(i,j,l)
                      f(i)=alpha*(rm(i,j,l,n)+(one-alpha)*rxm(i,j,l,n))
                      pf(i)=am(i,j,l)*(alpha*alpha*rxm(i,j,l,n) - 3.*f(i))
                      fy(i)=alpha*rym(i,j,l,n)
                      fz(i)=alpha*rzm(i,j,l,n)
                   else
                      alpha=am(i,j,l)/m(i+1,j,l)
                      f(i)=alpha*(rm(i+1,j,l,n)-(one+alpha)*rxm(i+1,j,l,n))
                      pf(i)=am(i,j,l)*(alpha*alpha*rxm(i+1,j,l,n) - 3.*f(i))
                      fy(i)=alpha*rym(i+1,j,l,n)
                      fz(i)=alpha*rzm(i+1,j,l,n)
                   end if
                   !xi(region,1)=max(xi(region,1),abs(alpha))
                   if(abs(alpha)>=one) then
                      if (okdebug) then
                         print *,'dynamu: CFL violation: alpha,m=',alpha, &
                              m(i,j,l),'  in region,i,j,l: ',region,i,j,l
                      end if
                      l_nxi=l_nxi+1.0
                   end if
                   
                end do
                
                if (im(region)/xref(region)==im(1)) then
                   ! redefine left fluxes using periodicity:
                   ! reduced grid
                   if (special_grid) iimr = imred(j,region)
                   f(0)  = f (iimr)
                   pf(0) = pf(iimr)
                   fy(0) = fy(iimr)
                   fz(0) = fz(iimr)
                end if

#ifdef with_budgets    
                ! 
                ! add up flux budget!
                ! note that this fails if grid is reduced (ie check)
                !
                if ( apply_budget_global_flux ) then 
                  if ( (iflux1(region)-1 >= is) .and. (iflux1(region)-1 <= iie) ) then 
                     budget_flux(region)%flux_x1(j,l,nglob) =  budget_flux(region)%flux_x1(j,l,nglob) &
                        + f(iflux1(region)-1)*1e3/ra(nglob)   ! moles
                  endif
                  if ( (iflux2(region)-1 >= is) .and. (iflux2(region)-1 <= iie) ) then 
                     budget_flux(region)%flux_x2(j,l,nglob) =  budget_flux(region)%flux_x2(j,l,nglob) &
                        + f(iflux2(region)-1)*1e3/ra(nglob)
                  endif
                end if
#endif
                !   
                ! calculate new tracer mass, and tracer mass slopes
                !
                do i=is,iie
                   rm(i,j,l,n)=rm(i,j,l,n)+(f(i-1)-f(i))
                   rxm(i,j,l,n)=rxm(i,j,l,n)+(pf(i-1)-pf(i)       &
                        -(am(i-1,j,l)-am(i,j,l))*rxm(i,j,l,n)     &
                        +3.*((am(i-1,j,l)+am(i,j,l))*rm(i,j,l,n)  &
                        -(f(i-1)+f(i))*m(i,j,l)))/mnew(i)
                   !CMK: apply limits again: might be in nloop!
                   if ( limits ) rxm(i,j,l,n) = &
                        max(min(rxm(i,j,l,n),rm(i,j,l,n)),-rm(i,j,l,n))
                   rym(i,j,l,n)=rym(i,j,l,n)+(fy(i-1)-fy(i))
                   rzm(i,j,l,n)=rzm(i,j,l,n)+(fz(i-1)-fz(i))
                end do !forall
             end do  ! end of n-loop
             ! store new air mass in m array
             m(is:iie,j,l)=mnew(is:iie)

          end do  ! iloop

          ! restore 'old' am
          am(is-1:iie,j,l) = am(is-1:iie,j,l)*nloop

       end do
    end do    ! end of l-loop/j-loop

    call my_omp_max (l_nl, g_nl)
    call my_omp_max (l_xi, g_xi)
    call my_omp_sum (l_nxi, g_nxi)

    if (im(region)/xref(region)==im(1)) then
       ! periodic boundary condition
       do j=jss,jee
          ! reduced grid
          if (special_grid) iimr = imred(j,region)
          rm (0,    j,:,:) = rm (iimr,j,:,:)
          rxm(0,    j,:,:) = rxm(iimr,j,:,:)
          rym(0,    j,:,:) = rym(iimr,j,:,:)
          rzm(0,    j,:,:) = rzm(iimr,j,:,:)
          rm (iimr+1,j,:,:) = rm (1,  j,:,:)
          rxm(iimr+1,j,:,:) = rxm(1,  j,:,:)
          rym(iimr+1,j,:,:) = rym(1,  j,:,:)
          rzm(iimr+1,j,:,:) = rzm(1,  j,:,:)
          m  (0    ,j,:)   = m  (iimr,j,:)
          m  (iimr+1,j,:)   = m  (1,  j,:)
       end do
    end if

    deallocate (f)
    deallocate (pf)
    deallocate (fy)
    deallocate (fz)
    deallocate (mnew)

!$OMP END PARALLEL

    nloop_max(region,1) = max (nloop_max(region,1), nint (g_nl))

    xi(region,1)  = max (xi(region,1), g_xi)
    nxi(region,1) = nxi(region,1) + nint (g_nxi)

    nullify(am)
    nullify(m)
    nullify(rm)
    nullify(rxm)
    nullify(rym)
    nullify(rzm)

    call uncompress_xedges(region,js,je,ls,le)
    call mix_edges(region)

  end subroutine dynamu


  subroutine compress_xedges(region,js,je,ls,le)
    !
    ! this is for slope only: condense data at the x-edges of the zoom region
    ! to allow for uniform work in dynamu/dynamu1
    ! written by mike botchev, march-june 1999
    !
    use dims,        only : xref, parent, im, okdebug
    use global_data, only : mass_dat
    use meteodata  , only : m_dat
    use partools   , only : ntracetloc
    use chem_param,  only : ntracet

    ! input
    integer,intent(in) :: region
    integer,intent(in) :: js
    integer,intent(in) :: je
    integer,intent(in) :: ls
    integer,intent(in) :: le

    ! local
    real,dimension(:,:,:,:),pointer    :: rm,rxm,rym,rzm
    real,dimension(:,:,:),  pointer    :: m
    integer :: imr,j,l,n,xref_

    ! start

    m => m_dat(region)%data
    rm => mass_dat(region)%rm_t
    rxm => mass_dat(region)%rxm_t
    rym => mass_dat(region)%rym_t
    rzm => mass_dat(region)%rzm_t
    
    xref_ = xref(region)/xref(parent(region))
    imr = im(region)
    if ( (xref_ == 1 ) .or. ( im(region)/xref(region) == im(1) ) ) then
       if ( okdebug ) print *,"compress_xedges:", &
            " no refinement or periodic bc's, nothing to do"
       return
    end if
    if ( okdebug ) then
       print *,'compress_xedges: region=',region,' imr=',imr
       print *,'compress_xedges:        cells to be updated: ', &
            xref_-1,xref_,imr-xref_+1,imr-xref_+2
    end if

    !forall(j=js:je,l=ls:le)
    do j=js,je
       do l=ls,le
          !sum_m = sum(m(1:xref_,j,l),1);  m (xref_,  j,l) = sum_m
          m (xref_,  j,l) = sum(m(1:xref_,j,l))
          m (xref_-1,j,l) = m (0,j,l) 
          m (imr-xref_+1,j,l) = sum(m(imr-xref_+1:imr,j,l))
          m (imr-xref_+2,j,l) = m (imr+1,j,l) 
          
          !forall(n=1:ntracet)
#ifdef MPI   
          do n=1,ntracetloc
#else
          do n=1,ntracet
#endif
             rm (xref_,  j,l,n) = sum(rm(1:xref_,j,l,n))
             rm (xref_-1,j,l,n) = rm (0,j,l,n) 
             rxm(xref_-1,j,l,n) = rxm(0,j,l,n)
             rym(xref_-1,j,l,n) = rym(0,j,l,n) 
             rzm(xref_-1,j,l,n) = rzm(0,j,l,n)
             
             rm (imr-xref_+1,j,l,n) = sum(rm(imr-xref_+1:imr,j,l,n))
             rm (imr-xref_+2,j,l,n) = rm (imr+1,j,l,n) 
             rxm(imr-xref_+2,j,l,n) = rxm(imr+1,j,l,n)
             rym(imr-xref_+2,j,l,n) = rym(imr+1,j,l,n) 
             rzm(imr-xref_+2,j,l,n) = rzm(imr+1,j,l,n)
          end do !forall
       end do
    end do !forall

    nullify(m)
    nullify(rm)
    nullify(rxm)
    nullify(rym)
    nullify(rzm)

  end subroutine compress_xedges



  subroutine uncompress_xedges(region,js,je,ls,le)
    !
    ! distribute data at the x-edges of the zoom region
    ! over the whole interface cell
    ! written by mike botchev, march-june 1999
    !
    use dims,        only : xref, parent, im, okdebug
    use meteodata  , only : m_dat
    use global_data, only : mass_dat
    use partools   , only : ntracetloc
    use chem_param,  only : ntracet

    implicit none

    ! input/output
    integer,intent(in) :: region
    integer,intent(in) :: js
    integer,intent(in) :: je
    integer,intent(in) :: ls
    integer,intent(in) :: le

    ! local
    real,dimension(:,:,:,:),pointer    :: rm,rxm,rym,rzm
    real,dimension(:,:,:),  pointer    :: m
    integer :: imr,j,l,n,xref_
    real    :: m_edge,rm_edge

    ! start

    m => m_dat(region)%data
    rm => mass_dat(region)%rm_t
    rxm => mass_dat(region)%rxm_t
    rym => mass_dat(region)%rym_t
    rzm => mass_dat(region)%rzm_t
    
    xref_ = xref(region)/xref(parent(region))
    imr = im(region)

    if ( (xref_==1) .or. (im(region)/xref(region)==im(1)) ) then
       if ( okdebug ) print *,"uncompress_xedges:", &
            " no refinement or periodic bc's, nothing to do"
       return
    end if
    if ( okdebug ) then
       print *,'uncompress_xedges: region=',region,' imr=',imr
       print *,'uncompress_xedges: cells to be updated: ',1,':', &
            xref_,' ',imr-xref_+1,':',imr
    end if

    ! forall(j=js:je,l=ls:le)
    do j=js,je
       do l=ls,le
          m_edge = m (xref_,j,l) 
          m (1:xref_,j,l) = m_edge/xref_
          
          m_edge = m (imr-xref_+1,j,l) 
          m (imr-xref_+1:imr,j,l) = m_edge/xref_ 
          !
          !   !forall(n=1:ntracet)
#ifdef MPI   
          do n=1,ntracetloc
#else
          do n=1,ntracet
#endif
             rm_edge = rm (xref_,j,l,n) 
             rm (1:xref_,j,l,n) = rm_edge/xref_
             !      
             !*      rxm(1:xref_-1,j,l,n) = rxm(xref_,j,l,n) 
             !*      rym(1:xref_-1,j,l,n) = rym(xref_,j,l,n) 
             !*      rzm(1:xref_-1,j,l,n) = rzm(xref_,j,l,n) 
             rxm(1:xref_,j,l,n) = rxm(xref_,j,l,n)/xref_ 
             rym(1:xref_,j,l,n) = rym(xref_,j,l,n)/xref_ 
             rzm(1:xref_,j,l,n) = rzm(xref_,j,l,n)/xref_ 
             !      
             rm_edge = rm (imr-xref_+1,j,l,n) 
             rm (imr-xref_+1:imr,j,l,n) = rm_edge/xref_ 
             !
             !*      rxm(imr-xref_+2:imr,j,l,n) = rxm(imr-xref_+1,j,l,n)
             !*      rym(imr-xref_+2:imr,j,l,n) = rym(imr-xref_+1,j,l,n)
             !*      rzm(imr-xref_+2:imr,j,l,n) = rzm(imr-xref_+1,j,l,n)
             rxm(imr-xref_+1:imr,j,l,n) = rxm(imr-xref_+1,j,l,n)/xref_
             rym(imr-xref_+1:imr,j,l,n) = rym(imr-xref_+1,j,l,n)/xref_
             rzm(imr-xref_+1:imr,j,l,n) = rzm(imr-xref_+1,j,l,n)/xref_
          end do
       end do !forall
    end do !forall

    nullify(m)
    nullify(rm)
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
    
  end subroutine uncompress_xedges



  subroutine put_xedges(region)
    !
    ! passes values along the x-boundaries to all the children 
    ! written by mike botchev, march-june 1999
    ! structure implemented by MK, dec 2002
    !
    use dims,          only : im, jm, lm, xref, yref, zref
    use dims,          only : ibeg, jbeg, lbeg, iend, nregions
    use dims,          only : children, okdebug
    use meteodata    , only : m_dat
    use global_data  , only : mass_dat
#ifdef with_budgets    
    use budget_global, only : sum_advection
#endif
    use chem_param,    only : ra, ntracet
#ifdef MPI
    use mpi_const
#endif
    use partools     , only : lmloc, offsetl
    use partools     , only : ntracetloc
    use partools     , only : which_par, previous_par
    use tracer_data  , only : Par_Check_Domain

    implicit none

    ! in/out
    integer,intent(in)   :: region

    ! local
    real,dimension(:,:,:,:),pointer     :: rm,rxm,rym,rzm,rmc,rxmc,rymc,rzmc
    real,dimension(:,:,:),  pointer     :: m,mc
    integer            :: child,ichild,j,jp,l,lp,jmc,lmc,n,i
    integer            :: xref_,yref_,zref_
    real               :: yzref,xyzref
    real               :: sum_old,sum_new,sum_old_all,sum_new_all
    integer            :: nzne, nzne_v,nznem
    real,dimension(:,:),allocatable :: toc0,toc1,tocm,tocm1
    integer :: communicator,root_id,lglob,nglob,lmr,nt
    !integer ::  mlev0

    ! start

    which_par=previous_par(region)

    !WP! make sure parents are on same domain
    call Par_Check_Domain( region, 'c', which_par )

    if ( which_par == 'tracer' .and. ntracetloc == 0 ) return 
    if ( which_par == 'levels' .and. lmloc == 0 ) return

    ichild = 0
    do while(ichild<children(region,0))
       ichild = ichild + 1
       child = children(region,ichild)
       xref_ = xref(child)/xref(region)
       yref_ = yref(child)/yref(region)
       zref_ = zref(child)/zref(region)

       if (im(child)/xref(child)==im(1)) then
          ! this child is [0,360] degrees wide - skip it since 
          ! periodic boundary conditions will be used
          if (okdebug)  & 
               print *,'put_xedges: child ',child,' skipped'
          cycle 
       endif

       m => m_dat(region)%data

       if ( which_par == 'tracer' ) then

          rm => mass_dat(region)%rm_t
          rxm => mass_dat(region)%rxm_t
          rym => mass_dat(region)%rym_t
          rzm => mass_dat(region)%rzm_t

          mc => m_dat(child)%data
          rmc => mass_dat(child)%rm_t
          rxmc => mass_dat(child)%rxm_t
          rymc => mass_dat(child)%rym_t
          rzmc => mass_dat(child)%rzm_t

          lmr = lm(region)
          nt=ntracetloc
#ifdef MPI     
          communicator=com_trac  !WP! assign com_trac as communicator
          root_id=root_t
#endif

       else if ( which_par == 'levels' ) then

          rm => mass_dat(region)%rm_k
          rxm => mass_dat(region)%rxm_k
          rym => mass_dat(region)%rym_k
          rzm => mass_dat(region)%rzm_k

          mc => m_dat(child)%data
          rmc => mass_dat(child)%rm_k
          rxmc => mass_dat(child)%rxm_k
          rymc => mass_dat(child)%rym_k
          rzmc => mass_dat(child)%rzm_k

          lmr = lmloc
          nt=ntracet
#ifdef MPI   
          communicator=com_lev  !WP! assign com_lev as communicator
          root_id=root_k
#endif
       end if
       
       ! ~~ mass (all levels on each pe)

       jmc = jm(child)
       lmc = lm(child)
       
       allocate(toc0(jmc,lmc))
       allocate(toc1(jmc,lmc))
       allocate(tocm(jmc,lmc))
       allocate(tocm1(jmc,lmc))
       
       yzref = 1./(yref_*zref_)
       xyzref = 1./(xref_*yref_*zref_)
       
       ! loop through the cells of x-walls of the child
       do l=1,lmc
          lp = lbeg(child) + (l-1)/zref_
          do j=1,jmc
             jp = jbeg(child) + (j-1)/yref_
             toc0(j,l) = m(ibeg(child)-1,jp,lp)*yzref
             toc1(j,l) = m(ibeg(child)  ,jp,lp)*xyzref
             tocm(j,l) = m(iend(child)  ,jp,lp)*xyzref
             tocm1(j,l) = m(iend(child)+1,jp,lp)*yzref
          end do
       end do
       !     copy info to the child ....
       mc(          0,1:jmc,1:lmc) = toc0
       mc(im(child)+1,1:jmc,1:lmc) = tocm1
       do i=1,xref_
          mc(i            , 1:jmc , 1:lmc ) = toc1
          mc(im(child)+1-i, 1:jmc , 1:lmc ) = tocm
       end do
       nullify(mc)
       
       deallocate(toc0)
       deallocate(toc1)
       deallocate(tocm)
       deallocate(tocm1)

       ! ~~ tracers

       jmc = jm(child) 
       lmc = lmr 

       allocate(toc0(jmc,lmc))
       allocate(toc1(jmc,lmc))
       allocate(tocm(jmc,lmc))
       allocate(tocm1(jmc,lmc))
       
       sum_old = 0.0
       sum_old_all = 0.0
       sum_new = 0.0
       sum_new_all = 0.0

       do n=1,nt
          do l=1,lmc
             lp = lbeg(child) + (l-1)/zref_
             do j=1,jmc
                jp = jbeg(child) + (j-1)/yref_
                toc0(j,l) = rm(ibeg(child)-1,jp,lp,n)*yzref
                toc1(j,l) = rm(ibeg(child)  ,jp,lp,n)*xyzref
                tocm(j,l) = rm(iend(child)  ,jp,lp,n)*xyzref
                tocm1(j,l) =rm(iend(child)+1,jp,lp,n)*yzref
             end do
          end do

          rmc(          0,1:jmc,1:lmc,n) = toc0
          rmc(im(child)+1,1:jmc,1:lmc,n) = tocm1

          ! calculate the mass that is associated with the 'implicit' 
          !advection that is done here....

          nglob=n

#ifdef MPI   
          if ( myid > 0 .and. which_par == 'tracer' ) &
               nglob=sum(ntracet_ar(0:myid-1))+n
#endif

          do l=1,lmc
             lglob=l
#ifdef MPI      
             if(myid>0.and.which_par=='levels') lglob=sum(lmar(0:myid-1))+l
#endif
             do j=1,jmc
                if (n == 1) sum_new = sum_new +  xref_*(toc1(j,l) + tocm(j,l))
                do i=1,xref_
                   if (n == 1) sum_old = sum_old +  &
                        rmc(i,j,l,1) + rmc(im(child)-i+1,j,l,1)
                   rmc(i             ,j,l,n) = toc1(j,l)
                   rmc(im(child)+1-i ,j,l,n) = tocm(j,l)
                end do
             end do
          end do

       end do

#ifdef MPI    
       !WP! only on one when in tracer  
       if ( which_par == 'tracer' .and. myid /= pe_first_tracer ) sum_new=0.0
       !WP! only on one when in tracer  
       if ( which_par == 'tracer' .and. myid /= pe_first_tracer ) sum_old=0.0
#endif

       !WP! When in levels, each sum_new and sum_old carries the full sum 
       !WP! of each level for tracer one.
       !WP! When in tracer,sum_new and sum_old is zero on all PE's 
       !WP! except PE_first_tracer

#ifdef with_budgets    

#ifdef MPI    
       !add all
       call mpi_allreduce(sum_new,sum_new_all,1,my_real,mpi_sum,communicator,ierr)
       call mpi_allreduce(sum_old,sum_old_all,1,my_real,mpi_sum,communicator,ierr)
#else
       sum_new_all = sum_new
       sum_old_all = sum_old
#endif
       sum_advection(child) = sum_advection(child) + sum_new_all - sum_old_all

#ifdef MPI    
       if ( okdebug .and. myid == root_id ) &
            print *, 'put_xedges: region, sum_advection ', &
            child,sum_advection(child)
#else
       if ( okdebug ) print *, 'put_xedges: region, sum_advection ', &
            child,sum_advection(child)
#endif

#endif

       nullify(rmc)

       do n=1,nt
          do l=1,lmc
             lp = lbeg(child) + (l-1)/zref_
             do j=1,jmc
                jp = jbeg(child) + (j-1)/yref_
                toc0(j,l) = rxm(ibeg(child)-1,jp,lp,n)*yzref
                toc1(j,l) = rxm(ibeg(child)  ,jp,lp,n)*yzref   !note the yzref...
                tocm(j,l) = rxm(iend(child)  ,jp,lp,n)*yzref
                tocm1(j,l) =rxm(iend(child)+1,jp,lp,n)*yzref
             end do
          end do

          rxmc(          0,1:jmc,1:lmc,n) = toc0
          rxmc(im(child)+1,1:jmc,1:lmc,n) = tocm1
          do i=1,xref_
             rxmc(i            , 1:jmc , 1:lmc ,n) = toc1
             rxmc(im(child)+1-i, 1:jmc , 1:lmc ,n) = tocm
          end do

       end do

       nullify(rxmc)

       do n=1,nt
          do l=1,lmc
             lp = lbeg(child) + (l-1)/zref_
             do j=1,jmc
                jp = jbeg(child) + (j-1)/yref_
                toc0(j,l) = rym(ibeg(child)-1,jp,lp,n)*yzref
                toc1(j,l) = rym(ibeg(child)  ,jp,lp,n)*yzref   !note the yzref...
                tocm(j,l) = rym(iend(child)  ,jp,lp,n)*yzref
                tocm1(j,l) =rym(iend(child)+1,jp,lp,n)*yzref
             end do
          end do
          rymc(          0,1:jmc,1:lmc,n) = toc0
          rymc(im(child)+1,1:jmc,1:lmc,n) = tocm1
          do i=1,xref_
             rymc(i            , 1:jmc , 1:lmc ,n) = toc1
             rymc(im(child)+1-i, 1:jmc , 1:lmc ,n) = tocm
          end do
       end do
       nullify(rymc)

       do n=1,nt
          do l=1,lmc
             lp = lbeg(child) + (l-1)/zref_
             do j=1,jmc
                jp = jbeg(child) + (j-1)/yref_
                toc0(j,l) = rzm(ibeg(child)-1,jp,lp,n)*yzref
                toc1(j,l) = rzm(ibeg(child)  ,jp,lp,n)*yzref   !note the yzref...
                tocm(j,l) = rzm(iend(child)  ,jp,lp,n)*yzref
                tocm1(j,l) =rzm(iend(child)+1,jp,lp,n)*yzref
             end do
          end do

          rzmc(          0,1:jmc,1:lmc,n) = toc0
          rzmc(im(child)+1,1:jmc,1:lmc,n) = tocm1
          do i=1,xref_
             rzmc(i            , 1:jmc , 1:lmc ,n) = toc1
             rzmc(im(child)+1-i, 1:jmc , 1:lmc ,n) = tocm
          end do

       end do
       nullify(rzmc)
       deallocate(toc0)
       deallocate(toc1)
       deallocate(tocm)
       deallocate(tocm1)

       nullify(rm)
       nullify(rxm)
       nullify(rym)
       nullify(rzm)
       nullify(m)
    end do


  end subroutine put_xedges



end module advectx
