!### 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:      ADVECTY
!
! !DESCRIPTION: 
!\\
!\\
! !INTERFACE: 
!
MODULE ADVECTY
  !
  ! !USES:
  !
  USE GO,           ONLY : gol, goPr, goErr
  USE TM5_DISTGRID, ONLY : DGRID, GET_DISTGRID, UPDATE_HALO, REDUCE

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  PUBLIC :: advectyzoom
  !
  ! !REVISION HISTORY: 
  !    2 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    ADVECTYZOOM
  !
  ! !DESCRIPTION: wrapper around DYNAMV
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE ADVECTYZOOM( REGION )
    !
    ! !USES:
    !
    use dims,          only : lm
    use global_data,   only : wind_dat, mass_dat
    use partools,      only : isRoot
#ifdef with_budgets    
    use budget_global, only : sum_advection
#endif
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: region
    !
    ! !REVISION HISTORY: 
    ! 
    ! written by patrick berkvens and mike botchev, march-june 1999
    ! updated and modified by MK, dec 2002
    !
    !    2 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer :: istat
    real    :: sum_old, sum_new

    ! ------ START

#ifdef with_budgets    
    ! total mass tracer #1
    call REDUCE( dgrid(region), mass_dat(region)%rm(:,:,:,1), mass_dat(region)%halo, 'SUM', sum_old, istat)
#endif

    CALL DYNAMV(region)

#ifdef with_budgets    
    call REDUCE( dgrid(region), mass_dat(region)%rm(:,:,:,1), mass_dat(region)%halo, 'SUM', sum_new, istat)

    if ( isRoot ) sum_advection(region) = sum_advection(region) + sum_new - sum_old
#endif

  END SUBROUTINE ADVECTYZOOM
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DYNAMV
  !
  ! !DESCRIPTION: south-north tracer transport: calculate amount of tracer moved in a south-north
  !               advection substep.
  !
  !       method    : slopes scheme
  !       reference : Russel and Lerner, 1979
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DYNAMV( region )!, is, ie, ls, le)
    !
    ! !USES:
    !
    use dims,        only : im, jm, lm
    use dims,        only : okdebug, xi, nxi
    use dims,        only : limits, nloop_max, zero, one
    use global_data, only : wind_dat, mass_dat
    use meteodata  , only : m_dat
    use toolbox,     only : escape_tm
#ifdef with_budgets    
    use budget_global, only : budget_flux, jflux1, jflux2, apply_budget_global_flux
#endif
    use partools   , only : par_reduce
    use chem_param,  only : ntracet, ra
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in) :: region
    !
    ! !REVISION HISTORY: 
    !    programmed by : mh, mpi HH, feb-jun 1994
    !    zoom version written by mike botchev, march-june 1999
    !    2 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real, dimension(:,:,:,:), pointer :: rm, rxm, rym, rzm
    real, dimension(:,:,:),   pointer :: m, bm, mx
    real, dimension(:,:), allocatable :: mnew
    real, dimension(:,:), allocatable :: f,pf,fx,fz

    integer               :: i,j,l,n, is,ie, js,je, ls,le, iee,iss
    integer               :: i1, i2, j1, j2, jss, jee, status
    integer               :: imr,imr2,jmr,lmr
    real                  :: sfs,sfzs,sfn,sfzn
    real                  :: my_beta, beta,mxval
    integer,dimension(3)  :: mxloc
    integer               :: iloop, nloop(lm(region)), nglob, offsetn
    logical        :: cfl_ok, SouthPole, NorthPole, SouthBorder, NorthBorder
    integer        :: lrg, redfact, ixe, ixs
    real           :: summ
    integer        :: ns, ne, lss, lee
    real           :: l_xi, g_xi, l_nxi, g_nxi, g_nl

    integer,parameter     :: max_nloop = 6

    !------ Start

    ! Indices
    CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2,&
                       hasSouthPole=SouthPole,     hasNorthPole=NorthPole,       & 
                       hasSouthBorder=SouthBorder, hasNorthBorder=NorthBorder    )
    lmr = lm(region)
    jmr = jm(region)
    is=i1 ; ie=i2 ; js=j1 ; je=j2 ; ls=1 ; le=lmr

    ! Work arrays
    m   => m_dat(region)%data     ! halo = 2
    rm  => mass_dat(region)%rm    ! halo = 2
    rxm => mass_dat(region)%rxm   ! halo = 2
    rym => mass_dat(region)%rym   ! halo = 2
    rzm => mass_dat(region)%rzm   ! halo = 2
    bm  => wind_dat(region)%bm    ! halo = 1

    CALL UPDATE_HALO( dgrid(region),  bm, wind_dat(region)%halo, status)
    CALL UPDATE_HALO( dgrid(region),   m,    m_dat(region)%halo, status)
    CALL UPDATE_HALO( dgrid(region),  rm, mass_dat(region)%halo, status)
    CALL UPDATE_HALO( dgrid(region), rxm, mass_dat(region)%halo, status)
    CALL UPDATE_HALO( dgrid(region), rym, mass_dat(region)%halo, status)
    CALL UPDATE_HALO( dgrid(region), rzm, mass_dat(region)%halo, status)

    ! Reduce work domain if border of zoom region, and take care of poles    
    if (NorthBorder) then ! j2=jmr by design
       if (NorthPole) je=j2-1
    endif

    if (SouthBorder) then ! j1=1 by design
       if (SouthPole) js=j1+1
    endif

    ! indices for mass update
    jss=js
    jee=je
    if (SouthPole) jss=js-1  ! =j1=1
    if (NorthPole) jee=je+1  ! =j2=jmr

    
    ! if requested limit meridional slopes such that no negative 
    if (limits) then
       do n = 1,  ntracet
       do l = ls, le
       do j = js - 1, je + 1
       do i = is, ie 
          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
    end if

    g_nl  = 0.0
    l_xi  = 0.0
    l_nxi = 0.0

    allocate (   f( i1-1:i2+1, j1-1:j2+1 ))   ! halo=1
    allocate (  pf( i1-1:i2+1, j1-1:j2+1 ))   ! halo=1
    allocate (  fx( i1-1:i2+1, j1-1:j2+1 ))   ! halo=1
    allocate (  fz( i1-1:i2+1, j1-1:j2+1 ))   ! halo=1

    allocate ( mnew( i1:i2, j1:j2 )) ! halo=0
    
    allocate ( mx(i1-1:i2+1, j1-1:j2+1, ls:le)) ! halo=1

    ! ================= Find number of iterations needed (per layer)
    LEVELS: do l=ls,le

       nloop(l) = 1
       beta     = 2.0

       do while(abs(beta) >= one .and. nloop(l) < max_nloop)

          mx(is:ie,js-1:je+1,1) = m(is:ie,js-1:je+1,l) ! use 1st layer of max as a 2D work array

          xloop : do iloop = 1, nloop(l)
             cfl_ok = .true.
             
             IJ: do j=js,je+1
             do i=is,ie
                   if (bm(i,j,l)>=zero) then
                      my_beta=bm(i,j,l)/mx(i,j-1,1)
                   else
                      my_beta=bm(i,j,l)/mx(i,j,1)
                   end if
                   if (abs(my_beta)>=one) exit IJ
             end do
             end do IJ

             my_beta = abs(my_beta)
             call Par_Reduce(my_beta, 'MAX', beta, status, all=.true.)
             if(status==1) then
                write(gol,*) "error par_reduce";call goPr
             end if

             if (abs(beta)>=one) then
                cfl_ok = .false.
                exit xloop
             end if
             ! else...
             if (cfl_ok.and.(iloop/=nloop(l))) then
                mx(is:ie,jss:jee,1) = mx(is:ie,jss:jee,1) + bm(is:ie,jss:jee,l) - bm(is:ie,jss+1:jee+1,l)
                CALL UPDATE_HALO( dgrid(region), mx(:,:,1), 1, status)
             end if
             
          end do xloop
          
          if ( .not. cfl_ok ) then 
             beta      = 2.0
             bm(:,:,l) = bm(:,:,l)*nloop(l)/(nloop(l)+1) ! reduce mass flux
             if ( okdebug) print *,'dynamv: ', i,j,l, &
                  'nloop increased to', nloop(l) + 1, beta
             nloop(l) = nloop(l) + 1
             if(nloop(l) == max_nloop) call escape_tm('nloop too high in dynamv')
          end if

       end do !while

    end do LEVELS

    
    ! advect info
    g_nl = max(g_nl, REAL(maxval(nloop)))

    
    ! ================= Loop over tracers and vertical layers to update masses
    
    ! Store init M, so we can get the loop over tracer outside the CFL.
    ! Thus we can limit the mpi comm, and switch Trac and levels loops => faster
    mx(i1-1:i2+1, j1-1:j2+1, ls:le) = m(i1-1:i2+1, j1-1:j2+1, ls:le)

    TRAC: do n=1,ntracet

       m(i1-1:i2+1, j1-1:j2+1, ls:le) = mx(i1-1:i2+1, j1-1:j2+1, ls:le) !reset

       LEV2: do l=ls,le

          CFL: DO iloop = 1,nloop(l)

          !-- calculate new air mass distribution
          mnew(is:ie,jss:jee) = m(is:ie,jss:jee,l) + &
                               bm(is:ie,jss:jee,l) - bm(is:ie,jss+1:jee+1,l)

          !-- compute all inner fluxes
          !   f(*,j,*) is flux entering j-th cell from below?cmk    
          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

          ! Case of JS
          if ( SouthPole ) then ! js=2         
             do i=is,ie

                ! this is the specific for js-1:
                fz(i,1)=zero
                fx(i,1)=zero
                pf(i,1)=zero
                f(i,1)=zero

                if (bm(i,2,l)>=zero) then ! this is the specific case
                   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 ! like above
                   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   ! do like above
             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...

          ! Case of JE+1
          if ( NorthPole ) then ! je+1=jmr
             do i=is,ie
                if (bm(i,jmr,l)>=zero) then !like above

                   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 ! this is the specific case

                   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
          !
          ! Finished computing fluxes. Now apply
          ! 
          !-- first compute INCOMING fluxes from the south...
          if (apply_budget_global_flux) then

             do i=is,ie
                if ((jflux1(region) <= j2 ) .and. (jflux1(region) >= j1 )) then
                   budget_flux(region)%flux_y1(i,l,n) = budget_flux(region)%flux_y1(i,l,n) + &
                        f(i,jflux1(region))*1e3/ra(n)
                end if

                if ((jflux2(region) <= j2 ) .and. (jflux2(region) >= j1 )) then

                   budget_flux(region)%flux_y2(i,l,n) = budget_flux(region)%flux_y2(i,l,n) + &
                        f(i,jflux2(region))*1e3/ra(n)
                end if
             enddo
          endif
#endif             
          !   
          !-- Calculate new tracer mass, and tracer mass slopes
          !
          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

           ! case of JS-1 if south pole (js=2)
           if ( SouthPole ) then
              rm (is:ie, 1,l,n) =  rm(is:ie, 1,l,n) -  f(is:ie,2)
              rzm(is:ie, 1,l,n) = rzm(is:ie, 1,l,n) - fz(is:ie,2)
           end if

           ! case of JE+1, only if north pole (je+1=jmr)
           if ( NorthPole ) then
              rm (is:ie,jmr,l,n) =  rm(is:ie,jmr,l,n) +  f(is:ie,jmr)
              rzm(is:ie,jmr,l,n) = rzm(is:ie,jmr,l,n) + fz(is:ie,jmr)
           end if

           ! store new air mass in m array   
           m(is:ie,jss:jee,l) = mnew(is:ie,jss:jee)

           ! update anything that changes and may be used with +/- index shift,
           ! but only if not the last step
           if (iloop/=nloop(l)) then
              !print*, "PLS-test-halo" (FIXME: j_only is not yet implemented but should reduce comm by a factor of 2 here)
              CALL UPDATE_HALO( dgrid(region),   m(:,:,l),      m_dat(region)%halo, status, j_only=.true.)
              CALL UPDATE_HALO( dgrid(region),  rm(:,:,l,n), mass_dat(region)%halo, status, j_only=.true.)
              CALL UPDATE_HALO( dgrid(region), rxm(:,:,l,n), mass_dat(region)%halo, status, j_only=.true.)
              CALL UPDATE_HALO( dgrid(region), rym(:,:,l,n), mass_dat(region)%halo, status, j_only=.true.)
              CALL UPDATE_HALO( dgrid(region), rzm(:,:,l,n), mass_dat(region)%halo, status, j_only=.true.)
           end if

        end do CFL ! iloop

    end do LEV2 ! end of l-loop over vertical layers....
    end do TRAC ! n-loop

    ! restore old bm's
    do l=ls,le
       bm(:,:,l) = bm(:,:,l)*nloop(l)
    end do
    
    ! Gather more alorithm information
    CALL PAR_REDUCE( l_xi, 'max', g_xi, status, .true.)
    CALL PAR_REDUCE( l_nxi, 'sum', g_nxi, status, .true.)

    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)

    ! Done
    deallocate (mx, mnew)
    deallocate (f)
    deallocate (pf)
    deallocate (fx)
    deallocate (fz)

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

  END SUBROUTINE DYNAMV
  !EOC

END MODULE ADVECTY