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

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

  ! --- interface

  private

  public :: advectyzoom, put_yedges


contains

  subroutine advectyzoom(region)
    !
    ! set parameters for advecty 
    ! 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, xcyc
#ifdef with_budgets    
    use budget_global, only : sum_advection
#endif
    use global_data,   only : wind_dat, mass_dat
    use toolbox,       only : escape_tm
#ifdef MPI
    use mpi_const,     only : my_real,com_lev,ierr, &
                              root_t,mpi_sum, pe_first_tracer
    use mpi_comm,      only : barrier_k
#endif
    use partools     , only : myid, lmloc, ntracetloc

    implicit none

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

    ! local
    real,dimension(:,:,:),pointer               :: bm
    real,dimension(:,:,:,:),pointer             :: rm
    real,dimension(:,:),allocatable             :: bm0,bm1

    integer            :: is,ie,ls,le,n,q
    integer            :: imr,jmr,lmr,tref_,xref_,yref_,zref_
    logical            :: y_encountered
    character(len=1)   :: dir
    real               :: sum_old, sum_new,sum_old_all,sum_new_all

    ! start

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

    call put_yedges(region)

    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(bm0(0:im(region)+1,0:lmr+1))
    allocate(bm1(0:im(region)+1,0:lmr+1))

    bm => wind_dat(region)%bm_t
    rm =>  mass_dat(region)%rm_t

    sum_old = sum(rm(1:imr,1:jmr,1:lmr,1))

    ! determine the scope for advecty: 

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

    ! now track four following steps from q     !wp! changed to four steps

    do n=1,n_operators

       dir=splitorderzoom(region,q*(nsplitsteps/2)*tref_+n)

       select case(dir)
       case('x') 
          if  ((.not.y_encountered).or.(xcyc(region) == 1)) then
             is=1                             ! x-substep is before y => 
             ie=im(region)                    ! full i-scope
          else
             is=xref_+1                       ! y-substep is before x => 
             ie=im(region)-xref_              ! restricted i-scope
          end if
       case('y')
          y_encountered=.true.
       case('z')
          if (.not.y_encountered) then
             ls=1                             ! z-substep is before y => 
             le=lmr                    ! full l-scope 
          else
             ls=zref_+1                       ! y-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 *,'advectyzoom: strange value in splitorderzoom(',region,',',  &
               q*(nsplitsteps/2)*tref_+n,'): ',q
          call escape_tm('advectyzoom: Error in do_next3')
       end select

    end do

#ifdef slopes
#ifndef secmom
    if ( mod(status(region),(nsplitsteps/2)*tref_) >= n_operators ) then
       ! if (1) more than fourth substep and (2) slope scheme then
       ! at this substep no coarse fluxes will be applied 

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

       if ( okdebug ) print *,'advectyzoom: zeroing out outer ', &
            'mass fluxes (yref)',yref_
       if ( touch_sp(region) /= 1) then
          bm0 = bm(:,      yref_,:);   bm(:,      yref_,:) = zero
       else
          if ( okdebug) print *,'advectyzoom: skipping zeroing outer ', &
               'mass flux at SP' 
       end if
       if ( touch_np(region) /= 1) then
          bm1 = bm(:,jmr-yref_+2,:);   bm(:,jmr-yref_+2,:) = zero
       else
          if ( okdebug) print *,'advectyzoom: skipping zeroing outer ', &
               'mass flux at NP' 
       end if

    end if
#endif
#endif

    call advecty_work(region,is,ie,ls,le)

#ifdef slopes
#ifndef secmom
    if (mod(status(region),(nsplitsteps/2)*tref_) >= n_operators ) then

       ! for slopes only: recreate velocities via the x-edges 
       if ( touch_sp(region) /= 1 ) bm(:,      yref_,:) = bm0
       if ( touch_np(region) /= 1 ) bm(:,jmr-yref_+2,:) = bm1

    end if
#endif
#endif

#ifdef with_budgets    
    sum_new = sum(rm(1:imr,1:jmr,1:lmr,1))
#ifdef MPI
    if ( myid==pe_first_tracer ) &
         sum_advection(region) = sum_advection(region) + sum_new - sum_old
    if ( okdebug .and. myid == root_t ) &
         print *, 'advectyzoom: advect_y, 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 *, 'advectyzoom: advect_y, region, sum_advection ', &
         region,sum_advection(region),sum_new,sum_old
#endif
#endif

    deallocate(bm0)
    deallocate(bm1)
    nullify(bm)

  end subroutine advectyzoom



  subroutine advecty_work(region,is,ie,ls,le)
    !
    ! calls dynamv
    ! written by mike botchev, march-june 1999
    !
    use toolbox, only : escape_tm

    implicit none

    ! input/output
    integer,intent(in) :: region
    integer,intent(in) :: is
    integer,intent(in) :: ie
    integer,intent(in) :: ls
    integer,intent(in) :: le

    ! local

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

    call dynamv(region,is,ie,ls,le)

  end subroutine advecty_work



  subroutine dynamv(region,is,ie,ls,le)
    !-----------------------------------------------------------------------
    !
    !****   dynamv          - south-north tracer transport  v 9.1
    !
    !       programmed by           mh      mpi HH          23-feb-1994
    !
    !       purpose
    !       -------
    !       calculate amount of tracer moved in a south-north advection
    !       substep
    !
    !       interface
    !       ---------
    !       call dynamv
    !
    !       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 : xref, yref, zref, tref, im, jm, lm
    use dims,        only : okdebug, parent, nregions, xi, nxi
    use dims,        only : touch_sp, touch_np, limits, nloop_max, zero, one
    use redgridZoom, only : nred, jred, imredj, clustsize
    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, jflux1, jflux2, apply_budget_global_flux
#endif
#ifdef MPI
    use mpi_comm,    only : barrier_k
#endif
    use partools   , only : ntracetloc, myid, ntracet_ar
    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) :: is
    integer,intent(in) :: ie
    integer,intent(in) :: ls
    integer,intent(in) :: le

    ! local
    real,dimension(:,:,:,:),pointer           :: rm,rxm,rym,rzm
    real,dimension(:,:,:),  pointer           :: m,bm
    real,dimension(:,:),allocatable           :: mnew
    real,dimension(:,:),allocatable           :: f,pf,fx,fz
    integer                                   :: i,j,je,js,l,n,iee,iss
    integer                                   :: imr,imr2,jmr,lmr
    integer                                   :: tref_,xref_,yref_,zref_
    real                                      :: sfs,sfzs,sfn,sfzn
    real                                      :: beta,mxval
    integer,dimension(3)                      :: mxloc
    integer                                   :: iloop, nloop, nglob, offsetn
    integer,parameter                         :: max_nloop = 6
    real,dimension(:,:), allocatable :: mx
    logical                          :: cfl_ok
    integer        :: lrg, redfact, ixe, ixs
    real           :: summ

    ! Omp parameters
    integer        :: ns, ne, lss, lee
    real           :: l_xi, g_xi, l_nxi, g_nxi, l_nl, g_nl

    ! start

    if ( okdebug ) print *,'dynamv: region=',region,' is,ie,ls,le=',is,ie,ls,le
    if ( ( region < 0 ) .or. ( region > nregions ) )  then
       call escape_tm('dynamv: STOP, illegal number of region !!!')
    end if

    ! 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)

    call compress_yedges(region,is,ie,ls,le)

    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
    bm => wind_dat(region)%bm_t


    ! check is,ie,ls,le:
    if ( ( is /= xref_+1 ) .and. ( is /= 1 ) ) &
         call escape_tm( 'dynamv: Wrong value for IS in dynamv')
    if ( ( ie /= imr-xref_ ) .and. ( ie /= imr ) ) &
         call escape_tm( 'dynamv: Wrong value for IE in dynamv')
    if ( ( ls /= zref_+1 ) .and. ( ls /= 1 ) )  &
         call escape_tm( 'dynamv: Wrong value for LS in dynamv')
    if ( ( le /= lmr-zref_ ) .and. ( le /= lmr ) ) &
         call escape_tm( 'dynamv: Wrong value for LE in dynamv')

    ! compute js/je -- cells is:ie,js:je:ls:le will be updated
    if (region==1) then
       js = 2
       je = jmr-1
    else
       js = yref_
       je = jmr-yref_+1
       if ( touch_sp(region) == 1) js = 2    !cmk....new option.. reduced grid
       if ( touch_np(region) == 1) je = jmr-1  !exclude poles
    end if

    ! if requested limit meridional slopes such that no negative 
    if (limits) then
!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (is, ie, js, je, ls, le, ntracetloc, rm, rym) &
!$OMP  private (i, j, l, n, ns, ne)
       call my_omp_domain (1, ntracetloc, ns, ne)
       do n = ns, ne
          do l = ls, le
             do j = js - 1, je + 1
                do i = is, ie ! Otherwise negative
                   rym(i,j,l,n) = max (min (rym(i,j,l,n), rm(i,j,l,n)), &
                                  -rm(i,j,l,n))
                end do
             end do
          end do
       end do
!$OMP END PARALLEL
    end if

    ! loop over tracers and vertical layers

!$OMP PARALLEL &
!$OMP  default (none) &
#ifdef with_budgets
!$OMP  shared  (apply_budget_global_flux, budget_flux, jflux1, jflux2,ra) &
#endif
!$OMP  shared  (region, is, ie, imr, js, je, jmr, ls, le, m, bm, &
!$OMP           okdebug, ntracetloc, ntracet_ar, myid, rm, rym, rxm, &
!$OMP           rzm, limits, g_nl, g_xi, g_nxi) &
!$OMP  private (i, j, l, lss, lee, l_nl, l_xi, l_nxi, f, pf, fx, fz, &
!$OMP           mnew, nloop, beta, mx, cfl_ok, offsetn, nglob)

    l_nl  = 0.0
    l_xi  = 0.0
    l_nxi = 0.0

    allocate (f(im(region),0:jm(region)))
    allocate (pf(im(region),0:jm(region)))
    allocate (fx(im(region),0:jm(region)))
    allocate (fz(im(region),0:jm(region)))
    allocate (mnew(im(region),jm(region)))

    call my_omp_domain (ls, le, lss, lee)
    do l=lss,lee

       ! compute all inner fluxes
       ! f(*,j,*) is flux entering j-th cell from below?cmk
       
       nloop = 1   ! determine nloop per layer!
       beta = 2.0
       allocate(mx(is:ie,js-1:je+1))
       do while(abs(beta) >= one .and. nloop < max_nloop)
          mx(is:ie,js-1:je+1) = m(is:ie,js-1:je+1,l)
          
          do iloop = 1, nloop
             cfl_ok = .true.
             do j=js+1,je
                do i=is,ie
                   if (bm(i,j,l)>=zero) then
                      beta=bm(i,j,l)/mx(i,j-1)
                   else
                      beta=bm(i,j,l)/mx(i,j)
                   end if
                   if (abs(beta)>=one) then
                      cfl_ok = .false.
                   end if
                end do  !i
             end do !j
             if ( region == 1.or.touch_sp(region) == 1) then         
                do i=is,ie
                   if (bm(i,2,l)>=zero) then
                      beta=bm(i,2,l)/mx(i,1)
                   else
                      beta=bm(i,2,l)/mx(i,2)
                   end if
                   if (abs(beta)>=one) then
                      cfl_ok = .false.
                   end if
                end do
             else   !zoom region not touching south pole
                j = js
                do i=is,ie   !no reduced grid allowed
                   if (bm(i,j,l)>=zero) then
                      beta=bm(i,j,l)/mx(i,j-1)
                   else
                      beta=bm(i,j,l)/mx(i,j)
                   end if
                   if (abs(beta)>=one) then
                      cfl_ok = .false.
                   end if
                end do
             end if ! compute boundary fluxes south pole...
             if ( region == 1.or.touch_np(region) == 1) then       ! north pole  
                do i=is,ie
                   if (bm(i,jmr,l)>=zero) then
                      beta=bm(i,jmr,l)/mx(i,jmr-1)
                   else
                      beta=bm(i,jmr,l)/mx(i,jmr)
                   end if
                   if (abs(beta)>=one) then
                      cfl_ok = .false.
                   end if
                end do
             else   !zoom region not touching north pole
                j = je+1
                do i=is,ie   !no reduced grid allowed
                   if (bm(i,j,l)>=zero) then
                      beta=bm(i,j,l)/mx(i,j-1)
                   else
                      beta=bm(i,j,l)/mx(i,j)
                   end if
                   if (abs(beta)>=one) then
                      cfl_ok = .false.
                   end if
                end do
             end if ! compute boundary fluxes north pole...
             if (cfl_ok) then
                if (region==1) then
                   mx(1:imr,1:jmr) = mx(1:imr,1:jmr) + &
                        bm(1:imr,1:jmr,l) - bm(1:imr,2:jmr+1,l)
                else if ( touch_sp(region) == 1) then
                   mx(1:imr,1:je) = mx(1:imr,1:je) + &
                        bm(1:imr,1:je,l) - bm(1:imr,2:je+1,l)
                else if ( touch_np(region) == 1) then
                   mx(1:imr,js:jmr) = mx(1:imr,js:jmr) + &
                        bm(1:imr,js:jmr,l) - bm(1:imr,js+1:jmr+1,l)
                else
                   mx(is:ie,js:je) = mx(is:ie,js:je) + &
                        bm(is:ie,js  :je  ,l) - bm(is:ie,js+1:je+1,l)      
                end if
             end if
          end do
          if ( .not. cfl_ok ) then 
             beta      = 2.0
             bm(:,:,l) = bm(:,:,l)*nloop/(nloop+1) ! reduce mass flux
             if ( okdebug) print *,'dynamv: ', i,j,l, &
                  'nloop increased to', nloop + 1, beta
             nloop = nloop + 1
             if(nloop == max_nloop) call escape_tm('nloop too high in dynamv')
          end if

       end do !while 

       deallocate(mx)
       l_nl = max(l_nl, REAL(nloop))

       do iloop = 1,nloop

          ! calculate new air mass distribution

          if (region==1) then
             mnew(1:imr,1:jmr) = m(1:imr,1:jmr,l) + &
                  bm(1:imr,1:jmr,l) - bm(1:imr,2:jmr+1,l)
          else if ( touch_sp(region) == 1) then
             mnew(1:imr,1:je) = m(1:imr,1:je,l) + &
                  bm(1:imr,1:je,l) - bm(1:imr,2:je+1,l)
          else if ( touch_np(region) == 1) then
             mnew(1:imr,js:jmr) = m(1:imr,js:jmr,l) + &
                  bm(1:imr,js:jmr,l) - bm(1:imr,js+1:jmr+1,l)
          else
             mnew(is:ie,js:je) = m(is:ie,js:je,l) + & 
                  bm(is:ie,js  :je  ,l) - bm(is:ie,js+1:je+1,l)      
          end if

          do n=1,ntracetloc
             offsetn=sum(ntracet_ar(0:myid-1))
             nglob = offsetn + n
             do j=js+1,je
                do i=is,ie
                   if (bm(i,j,l)>=zero) then
                      beta=bm(i,j,l)/m(i,j-1,l)
                      f(i,j)=beta*(rm(i,j-1,l,n)+(one-beta)*rym(i,j-1,l,n))
                      pf(i,j)=bm(i,j,l)*(beta*beta*rym(i,j-1,l,n)-3.*f(i,j))
                      fx(i,j)=beta*rxm(i,j-1,l,n)
                      fz(i,j)=beta*rzm(i,j-1,l,n)
                   else
                      beta=bm(i,j,l)/m(i,j,l)
                      f(i,j)=beta*(rm(i,j,l,n)-(one+beta)*rym(i,j,l,n))
                      pf(i,j)=bm(i,j,l)*(beta*beta*rym(i,j,l,n)-3.*f(i,j))
                      fx(i,j)=beta*rxm(i,j,l,n)
                      fz(i,j)=beta*rzm(i,j,l,n)
                   end if
                   l_xi=max(l_xi,abs(beta))       
                   if (abs(beta)>=one) then
                      l_nxi=l_nxi+1
                   end if
                end do
             end do
             
             !compute boundary fluxes
             if ( region == 1.or.touch_sp(region) == 1) then         
                do i=is,ie
                   fz(i,1)=zero
                   fx(i,1)=zero
                   pf(i,1)=zero
                   f(i,1)=zero
                   if (bm(i,2,l)>=zero) then
                      beta=bm(i,2,l)/m(i,1,l)
                      f(i,2)=beta*rm(i,1,l,n)   
                      pf(i,2)=-3.*bm(i,2,l)*f(i,2)
                      fx(i,2)=zero
                      fz(i,2)=beta*rzm(i,1,l,n)  
                   else
                      beta=bm(i,2,l)/m(i,2,l)
                      f(i,2)=beta*(rm(i,2,l,n)-(one+beta)*rym(i,2,l,n))
                      pf(i,2)=bm(i,2,l)*(beta*beta*rym(i,2,l,n)-3.*f(i,2))
                      fx(i,2)=beta*rxm(i,2,l,n)
                      fz(i,2)=beta*rzm(i,2,l,n)
                   end if
                   l_xi=max(l_xi,abs(beta))       
                   if (abs(beta)>=one) then
                      l_nxi=l_nxi+1
                   end if
                end do
             else   !zoom region not touching south pole
                j = js
                do i=is,ie   !no reduced grid allowed
                   if (bm(i,j,l)>=zero) then
                      beta=bm(i,j,l)/m(i,j-1,l)
                      f(i,j)=beta*(rm(i,j-1,l,n)+(one-beta)*rym(i,j-1,l,n))
                      pf(i,j)=bm(i,j,l)*(beta*beta*rym(i,j-1,l,n)-3.*f(i,j))
                      fx(i,j)=beta*rxm(i,j-1,l,n)
                      fz(i,j)=beta*rzm(i,j-1,l,n)
                   else
                      beta=bm(i,j,l)/m(i,j,l)
                      f(i,j)=beta*(rm(i,j,l,n)-(one+beta)*rym(i,j,l,n))
                      pf(i,j)=bm(i,j,l)*(beta*beta*rym(i,j,l,n)-3.*f(i,j))
                      fx(i,j)=beta*rxm(i,j,l,n)
                      fz(i,j)=beta*rzm(i,j,l,n)
                   end if
                   l_xi=max(l_xi,abs(beta))       
                   if (abs(beta)>=one) then
                      l_nxi=l_nxi+1
                   end if
                end do
             end if ! compute boundary fluxes south pole...
             if ( region == 1.or.touch_np(region) == 1) then       ! north pole  
                do i=is,ie
                   if (bm(i,jmr,l)>=zero) then
                      beta=bm(i,jmr,l)/m(i,jmr-1,l)
                      f(i,jmr)=beta*(rm(i,jmr-1,l,n)+(one-beta)*rym(i,jmr-1,l,n))
                      pf(i,jmr)=bm(i,jmr,l)*(beta*beta*rym(i,jmr-1,l,n) - &
                           3.*f(i,jmr))
                      fx(i,jmr)=beta*rxm(i,jmr-1,l,n)
                      fz(i,jmr)=beta*rzm(i,jmr-1,l,n)
                   else
                      beta=bm(i,jmr,l)/m(i,jmr,l)
                      ! for solid-body test and polar mixing
                      f(i,jmr)=beta*rm(i,jmr,l,n)
                      pf(i,jmr)=-3.*bm(i,jmr,l)*f(i,jmr)
                      fx(i,jmr)=zero
                      ! for solid-body test and polar mixing
                      fz(i,jmr)=beta*rzm(i,jmr,l,n)
                   end if
                   l_xi=max(l_xi,abs(beta))       
                   if (abs(beta)>=one) then
                      l_nxi=l_nxi+1
                   end if
                end do
             else   !zoom region not touching north pole
                j = je+1
                do i=is,ie   !no reduced grid allowed
                   if (bm(i,j,l)>=zero) then
                      beta=bm(i,j,l)/m(i,j-1,l)
                      f(i,j)=beta*(rm(i,j-1,l,n)+(one-beta)*rym(i,j-1,l,n))
                      pf(i,j)=bm(i,j,l)*(beta*beta*rym(i,j-1,l,n)-3.*f(i,j))
                      fx(i,j)=beta*rxm(i,j-1,l,n)
                      fz(i,j)=beta*rzm(i,j-1,l,n)
                   else
                      beta=bm(i,j,l)/m(i,j,l)
                      f(i,j)=beta*(rm(i,j,l,n)-(one+beta)*rym(i,j,l,n))
                      pf(i,j)=bm(i,j,l)*(beta*beta*rym(i,j,l,n)-3.*f(i,j))
                      fx(i,j)=beta*rxm(i,j,l,n)
                      fz(i,j)=beta*rzm(i,j,l,n)
                   end if
                   l_xi=max(l_xi,abs(beta))       
                   if (abs(beta)>=one) then
                      l_nxi=l_nxi+1
                   end if
                end do
             end if ! compute boundary fluxes north pole...
             
#ifdef with_budgets    
             !cmk finished computing fluxes...now apply...
             ! first compute INCOMING fluxes from the south...
             if (apply_budget_global_flux) then
                do i=is,ie
                   budget_flux(region)%flux_y1(i,l,nglob) = budget_flux(region)%flux_y1(i,l,nglob) + &
                     f(i,jflux1(region))*1e3/ra(nglob)
                   budget_flux(region)%flux_y2(i,l,nglob) = budget_flux(region)%flux_y2(i,l,nglob) + &
                     f(i,jflux2(region))*1e3/ra(nglob)
                enddo
             endif
#endif             
             do j=js,je
                do i=is,ie
                   rm(i,j,l,n) =rm(i,j,l,n)+(f(i,j)-f(i,j+1))
                   rym(i,j,l,n)=rym(i,j,l,n)+(pf(i,j)-pf(i,j+1)            &
                        -       (bm(i,j,l)-bm(i,j+1,l))*rym(i,j,l,n)       &
                        +       3.*((bm(i,j,l)+bm(i,j+1,l))*rm(i,j,l,n)    &
                        -       (f(i,j)+f(i,j+1))*m(i,j,l)))/mnew(i,j)
                   if ( limits)  rym(i,j,l,n) = max( min(rym(i,j,l,n), &
                        rm(i,j,l,n)), -rm(i,j,l,n) )
                   rxm(i,j,l,n)=rxm(i,j,l,n)+(fx(i,j)-fx(i,j+1))
                   rzm(i,j,l,n)=rzm(i,j,l,n)+(fz(i,j)-fz(i,j+1))
                end do
             end do !forall

             ! poles were excluded.....
             if ( region==1 .or. touch_sp(region) == 1 ) then
                rm(1:imr,  1,l,n) = rm(1:imr,  1,l,n) - f(1:imr,2)
                rzm(1:imr,1  ,l,n) = rzm(1:imr,  1,l,n) - fz(1:imr,2)
             end if
             if ( region==1 .or. touch_np(region) == 1 ) then
                rm(1:imr,jmr,l,n) = rm(1:imr,jmr,l,n) + f(1:imr,jmr)
                rzm(1:imr,jmr,l,n) = rzm(1:imr,jmr,l,n) + fz(1:imr,jmr)
             end if

          end do !      & n-loop

          ! store new air mass in m array   

          if ( region == 1) then
             m(1:imr,1:jmr,l) = mnew(1:imr,1:jmr)
          else if ( touch_sp(region) == 1) then
             m(1:imr,1:je,l) = mnew(1:imr,1:je)
          else if ( touch_np(region) == 1) then
             m(1:imr,js:jmr,l) = mnew(1:imr,js:jmr)
          else
             m(is:ie,js:je,l) = mnew(is:ie,js:je) 
          end if
          !
       end do !      iloop
       bm(:,:,l) = bm(:,:,l)*nloop     ! restore old bm's
       
    end do! end of l-loop over vertical layers....

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

    deallocate (mnew)
    deallocate (f)
    deallocate (pf)
    deallocate (fx)
    deallocate (fz)

!$OMP END PARALLEL

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

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

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

    call uncompress_yedges(region,is,ie,ls,le)

    call mix_edges(region)

  end subroutine dynamv



  subroutine compress_yedges(region,is,ie,ls,le)
    !
    ! this is for slope only: condense data at the y-edges of the zoom region
    ! to allow for uniform work in dynamv
    ! written by mike botchev, march-june 1999
    ! modified by MK, dec 2002
    !
    use dims,        only : yref, parent, jm, okdebug, touch_sp, touch_np
    use tracer_data, only : mass_dat
    use meteodata  , only : m_dat
    use partools   , only : ntracetloc
    use chem_param,  only : ntracet

    ! input/output
    integer,intent(in) :: region
    integer,intent(in) :: is
    integer,intent(in) :: ie
    integer,intent(in) :: ls
    integer,intent(in) :: le

    ! local
    real,dimension(:,:,:,:),pointer           :: rm,rxm,rym,rzm
    real,dimension(:,:,:),  pointer           :: m
    integer                                   :: jmr,i,l,n,yref_

    ! start

    yref_ = yref(region)/yref(parent(region))
    jmr = jm(region)
    if ((yref_==1).or.(region==1)) then
       if ( okdebug) print *,'compress_yedges: no refinement, nothing to do'
       return
    end if
    if ( okdebug) then
       print *,'compress_yedges: region=',region,' jmr=',jmr
       print *,'compress_yedges: cells to be updated: ', &
            yref_-1,yref_,jmr-yref_+1,jmr-yref_+2
       if ( touch_sp(region) == 1) &
            print *,'compress_yedges: SP will be skipped here'
       if ( touch_np(region) == 1) &
            print *,'compress_yedges: NP will be skipped here'
    end if

    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

    do l=ls,le
       do i=is,ie

          if ( touch_sp(region) /= 1) then
             m (i,yref_,  l) = sum(m(i,1:yref_,l))
             m (i,yref_-1,l) = m (i,0,l) 
          end if
          if ( touch_np(region) /= 1) then
             m (i,jmr-yref_+1,l) = sum(m(i,jmr-yref_+1:jmr,l))
             m (i,jmr-yref_+2,l) = m (i,jmr+1,l)
          end if
          do n=1,ntracetloc
             if ( touch_sp(region) /= 1) then
                rm (i,yref_,  l,n) = sum(rm(i,1:yref_,l,n))
                rm (i,yref_-1,l,n) = rm (i,0,l,n) 
                rxm(i,yref_-1,l,n) = rxm(i,0,l,n)
                rym(i,yref_-1,l,n) = rym(i,0,l,n) 
                rzm(i,yref_-1,l,n) = rzm(i,0,l,n)
             end if
             if ( touch_np(region) /= 1) then
                rm (i,jmr-yref_+1,l,n) = sum(rm(i,jmr-yref_+1:jmr,l,n))
                rm (i,jmr-yref_+2,l,n) = rm (i,jmr+1,l,n) 
                rxm(i,jmr-yref_+2,l,n) = rxm(i,jmr+1,l,n)
                rym(i,jmr-yref_+2,l,n) = rym(i,jmr+1,l,n) 
                rzm(i,jmr-yref_+2,l,n) = rzm(i,jmr+1,l,n)
             end if
          end do !forall

       end do
    end do !forall

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

  end subroutine compress_yedges



  subroutine uncompress_yedges(region,is,ie,ls,le)
    !
    ! distribute data at the y-edges of the zoom region
    ! over the whole interface cell
    ! written by mike botchev, march-june 1999
    !
    use dims,        only : yref, parent, jm, okdebug, touch_sp, touch_np
    use tracer_data, only : mass_dat
    use meteodata  , only : m_dat
    use partools   , only : ntracetloc
    use chem_param,  only : ntracet

    ! input/output
    integer,intent(in) :: region
    integer,intent(in) :: is
    integer,intent(in) :: ie
    integer,intent(in) :: ls
    integer,intent(in) :: le

    ! local
    real,dimension(:,:,:,:),pointer           :: rm,rxm,rym,rzm
    real,dimension(:,:,:),  pointer           :: m
    integer                                   :: jmr,i,l,n,yref_
    real                                      :: m_edge,rm_edge

    ! start

    yref_ = yref(region)/yref(parent(region))
    jmr = jm(region)
    if ((yref_==1).or.(region==1)) then
       if ( okdebug) print *,'uncompress_yedges: no refinement, nothnig to do'
       return
    end if
    if ( okdebug) then
       print *,'uncompress_yedges: region=',region,' jmr=',jmr
       print *,'uncompress_yedges: cells to be updated: ',1,':', &
            yref_,' ',jmr-yref_+1,':',jmr
       if ( touch_sp(region) == 1) &
            print *,'uncompress_yedges: SP will be skipped here'
       if ( touch_np(region) == 1) &
            print *,'uncompress_yedges: NP will be skipped here'
    end if

    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

    do l=ls,le
       do i=is,ie
          if ( touch_sp(region) /= 1) then
             m_edge = m (i,yref_,l) 
             m (i,1:yref_,l) = m_edge/yref_
          end if
          if ( touch_np(region) /= 1) then
             m_edge = m (i,jmr-yref_+1,l) 
             m (i,jmr-yref_+1:jmr,l) = m_edge/yref_ 
          end if

          !forall(n=1:ntracet)
          do n=1,ntracetloc

             if ( touch_sp(region) /= 1 ) then
                rm_edge = rm (i,yref_,l,n) 
                rm (i,1:yref_,l,n) = rm_edge/yref_
                !*      rxm(i,1:yref_-1,l,n) = rxm(i,yref_,l,n) 
                !*      rym(i,1:yref_-1,l,n) = rym(i,yref_,l,n) 
                !*      rzm(i,1:yref_-1,l,n) = rzm(i,yref_,l,n) 
                rxm(i,1:yref_,l,n) = rxm(i,yref_,l,n)/yref_ 
                rym(i,1:yref_,l,n) = rym(i,yref_,l,n)/yref_ 
                rzm(i,1:yref_,l,n) = rzm(i,yref_,l,n)/yref_ 

             end if
             if ( touch_np(region) /= 1 ) then
                rm_edge = rm (i,jmr-yref_+1,l,n) 
                rm (i,jmr-yref_+1:jmr,l,n) = rm_edge/yref_ 

                !*      rxm(i,jmr-yref_+2:jmr,l,n) = rxm(i,jmr-yref_+1,l,n)
                !*      rym(i,jmr-yref_+2:jmr,l,n) = rym(i,jmr-yref_+1,l,n)
                !*      rzm(i,jmr-yref_+2:jmr,l,n) = rzm(i,jmr-yref_+1,l,n)
                rxm(i,jmr-yref_+1:jmr,l,n) = rxm(i,jmr-yref_+1,l,n)/yref_
                rym(i,jmr-yref_+1:jmr,l,n) = rym(i,jmr-yref_+1,l,n)/yref_
                rzm(i,jmr-yref_+1:jmr,l,n) = rzm(i,jmr-yref_+1,l,n)/yref_
             end if
          end do
       end do
    end do

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

  end subroutine uncompress_yedges



  subroutine put_yedges(region)
    !
    ! passes values along the y-boundaries to all the children 
    ! written by mike botchev, march-june 1999
    ! adapted for cray-run Maarten Krol, march, 2000
    ! touch_np, touch_sp implemented...MK, nov, 2000
    !

    use dims,          only : xref, yref, zref, im, jm, lm
    use dims,          only : touch_sp, touch_np, okdebug, children
    use dims,          only : ibeg, jbeg, lbeg, jend, nregions
    use tracer_data  , only : mass_dat
    use meteodata    , only : m_dat
#ifdef with_budgets    
    use budget_global, only : sum_advection
#endif
    use chem_param,    only : ra, ntracet
    use partools     , only : offsetl, lmloc
    use partools     , only : ntracetloc
    use partools     , only : which_par, previous_par
#ifdef MPI
    use mpi_const
#endif
    use tracer_data  , only : Par_Check_Domain

    implicit none

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

    ! local
    real,dimension(:,:,:,:),pointer    :: rm,rxm,rym,rzm,rmc,rxmc,rymc,rzmc
    real,dimension(:,:,:),  pointer    :: m, mc

    real,dimension(:,:),allocatable    :: toc0,toc1,tocm,tocm1

    integer :: child, ichild, j, ip, l, lp, imc, lmc, n, i
    integer :: xref_, yref_, zref_
    real    :: xzref, xyzref
    real    :: sum_old, sum_new, sum_old_all, sum_new_all
    integer :: nzne,  nzne_v, nznem
    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  !WP! 

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

    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 ( touch_sp(child) == 1) then
          if (okdebug)  & 
               print *,'put_yedges: child ',child,' sp  skipped'
       end if
       if ( touch_np(child) == 1) then
          if (okdebug)  & 
               print *,'put_yedges: child ',child,' np  skipped'
       end if

       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)

       imc = im(child)
       lmc = lm(child)
       
       allocate(toc0(imc,lmc))
       allocate(toc1(imc,lmc))
       allocate(tocm(imc,lmc))
       allocate(tocm1(imc,lmc))
       
       xzref = 1./(xref_*zref_)
       xyzref = 1./(xref_*yref_*zref_)
       
       ! loop through the cells of y-walls of the child
       do l=1,lmc
          lp = lbeg(child) + (l-1)/zref_
          do i=1,imc
             ip = mod(ibeg(child)-1 + (i-1)/xref_,im(region)) + 1
             toc0(i,l) =  m(ip,jbeg(child)-1,lp)*xzref
             toc1(i,l) =  m(ip,jbeg(child)  ,lp)*xyzref
             tocm(i,l) =  m(ip,jend(child)  ,lp)*xyzref
             tocm1(i,l) = m(ip,jend(child)+1,lp)*xzref
          end do
       end do

       if ( touch_sp(child) == 0)  mc(1:imc,0          ,1:lmc) = toc0
       if ( touch_np(child) == 0)  mc(1:imc,jm(child)+1,1:lmc) = tocm1
       do j=1,yref_
          if ( touch_sp(child) == 0)  mc( 1:imc , j,            1:lmc ) = toc1
          if ( touch_np(child) == 0)  mc( 1:imc , jm(child)+1-j,1:lmc ) = tocm
       end do

       nullify(mc)

       deallocate(toc0)
       deallocate(toc1)
       deallocate(tocm)
       deallocate(tocm1)

       ! ~~ tracers

       imc = im(child) 
       lmc = lmr
       allocate(toc0(imc,lmc))
       allocate(toc1(imc,lmc))
       allocate(tocm(imc,lmc))
       allocate(tocm1(imc,lmc))

       sum_old_all = 0.0
       sum_new = 0.0
       sum_new_all = 0.0
       sum_old = 0.0

       do n=1,nt
          do l=1,lmc
             lp = lbeg(child) + (l-1)/zref_
             do i=1,imc
                ip = mod(ibeg(child)-1 + (i-1)/xref_,im(region)) + 1
                toc0(i,l) = rm(ip,jbeg(child)-1,lp,n)*xzref
                toc1(i,l) = rm(ip,jbeg(child)  ,lp,n)*xyzref
                tocm(i,l) = rm(ip,jend(child)  ,lp,n)*xyzref
                tocm1(i,l) =rm(ip,jend(child)+1,lp,n)*xzref
             end do
          end do

          if ( touch_sp(child) == 0)  rmc(1:imc,0          ,1:lmc,n) = toc0
          if ( touch_np(child) == 0)  rmc(1:imc,jm(child)+1,1:lmc,n) = tocm1

          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 i=1,imc

                if ( touch_sp(child) == 0 .and. n == 1)  &
                     sum_new = sum_new + yref_*toc1(i,l)
                if ( touch_np(child) == 0 .and. n == 1)  &
                     sum_new = sum_new + yref_*tocm(i,l)

                do j=1,yref_
                   if ( touch_sp(child) == 0)  then
                      if ( n == 1) sum_old = sum_old + rmc( i, j, l ,1)
                      rmc(i,j,l,n) = toc1(i,l)
                   end if
                   if ( touch_np(child) == 0 )  then
                      if ( n == 1 )  sum_old = sum_old + &
                           rmc(i,jm(child)+1-j,l,1)
                      rmc(i,jm(child)+1-j,l,n) = tocm(i,l)
                   end if
                end do !j
             end do !i
          end do !l
       end do  !n

#ifdef with_budgets    

#ifdef MPI   
       !WP! only on one when in tracer  
       if ( which_par == 'tracer' .and. myid /= pe_first_tracer ) sum_new=0.0
       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 
       !WP! sum 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 MPI   
       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_yedges: region, sum_advection ', &
            child,sum_advection(child)
#else
       if ( okdebug) print *, 'put_yedges: 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 i=1,imc
                ip = mod(ibeg(child)-1 + (i-1)/xref_,im(region)) + 1
                toc0(i,l) = rxm(ip,jbeg(child)-1,lp,n)*xzref
                toc1(i,l) = rxm(ip,jbeg(child)  ,lp,n)*xzref   !note the xzref...
                tocm(i,l) = rxm(ip,jend(child)  ,lp,n)*xzref
                tocm1(i,l) =rxm(ip,jend(child)+1,lp,n)*xzref
             end do
          end do

          if ( touch_sp(child) == 0 )  rxmc(1:imc,0          ,1:lmc,n) = toc0
          if ( touch_np(child) == 0 )  rxmc(1:imc,jm(child)+1,1:lmc,n) = tocm1
          do j=1,yref_
             if ( touch_sp(child) == 0 )  &
                  rxmc( 1:imc , j            ,1:lmc ,n) = toc1
             if ( touch_np(child) == 0 )  &
                  rxmc( 1:imc , jm(child)+1-j,1:lmc ,n) = tocm
          end do

       end do  !n
       nullify(rxmc)


       do n=1,nt
          do l=1,lmc
             lp = lbeg(child) + (l-1)/zref_
             do i=1,imc
                ip = mod(ibeg(child)-1 + (i-1)/xref_,im(region)) + 1
                toc0(i,l) = rym(ip,jbeg(child)-1,lp,n)*xzref
                toc1(i,l) = rym(ip,jbeg(child)  ,lp,n)*xzref   !note the xzref...
                tocm(i,l) = rym(ip,jend(child)  ,lp,n)*xzref
                tocm1(i,l) =rym(ip,jend(child)+1,lp,n)*xzref
             end do
          end do

          if ( touch_sp(child) == 0)  rymc(1:imc,0          ,1:lmc,n) = toc0
          if ( touch_np(child) == 0)  rymc(1:imc,jm(child)+1,1:lmc,n) = tocm1
          do j=1,yref_
             if ( touch_sp(child) == 0)  &
                  rymc( 1:imc , j            ,1:lmc ,n) = toc1
             if ( touch_np(child) == 0)  &
                  rymc( 1:imc , jm(child)+1-j,1:lmc ,n) = tocm
          end do

       end do  !n
       nullify(rymc)


       do n=1,nt
          do l=1,lmc
             lp = lbeg(child) + (l-1)/zref_
             do i=1,imc
                ip = mod(ibeg(child)-1 + (i-1)/xref_,im(region)) + 1
                toc0(i,l) = rzm(ip,jbeg(child)-1,lp,n)*xzref
                toc1(i,l) = rzm(ip,jbeg(child)  ,lp,n)*xzref   !note the xzref...
                tocm(i,l) = rzm(ip,jend(child)  ,lp,n)*xzref
                tocm1(i,l) =rzm(ip,jend(child)+1,lp,n)*xzref
             end do
          end do

          if ( touch_sp(child) == 0)  rzmc(1:imc,0          ,1:lmc,n) = toc0
          if ( touch_np(child) == 0)  rzmc(1:imc,jm(child)+1,1:lmc,n) = tocm1
          do j=1,yref_
             if ( touch_sp(child) == 0)  &
                  rzmc( 1:imc , j            ,1:lmc ,n) = toc1
             if ( touch_np(child) == 0)  &
                  rzmc( 1:imc , jm(child)+1-j,1:lmc ,n) = tocm
          end do

       end do  !n

       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_yedges


end module advecty
