!### 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
#define IF_NOTOK_MPI(action) if (ierr/=MPI_SUCCESS) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      ADVECTZ
!
! !DESCRIPTION: 
!\\
!\\
! !INTERFACE: 
!
MODULE ADVECTZ
  !
  ! !USES:
  !
  USE GO,           ONLY : gol, goPr, goErr
  USE TM5_DISTGRID, ONLY : DGRID, GET_DISTGRID, UPDATE_HALO, REDUCE

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

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DYNAMW
  !
  ! !DESCRIPTION: Calculate amount of tracer moved in a vertical advection substep.
  !
  !       method    : slopes scheme
  !       reference : Russel and Lerner, 1979
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DYNAMW
    !
    ! !USES:
    !
    use dims,        only : lm
    use dims,        only : zbeg, zend, epsz, nloop_max, zero, one
    use dims,        only : okdebug, limits, xi, nxi
    use global_data, only : wind_dat, mass_dat
    use meteodata  , only : m_dat
    use chem_param,  only : ntracet, ra
    use toolbox,     only : escape_tm
#ifdef with_budgets    
    use budget_global, only : budget_flux, lflux1, lflux2, apply_budget_global_flux
#endif
    use ParTools,    only : isRoot, par_reduce
    !
    ! !REVISION HISTORY:
    ! 
    ! - programmed by mh, mpi HH, 1994-1995
    ! - zoom version written by mike botchev, march-june 1999
    ! - 3 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !    - commented method #2, which may be an option if there is CFL violation. If not, it just
    !       adds extra communication.
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real, dimension(:,:,:,:), pointer   :: rm, rxm, rym, rzm
    real, dimension(:,:,:),   pointer   :: m, cm
    real, dimension(:,:,:), allocatable :: mnew, f, pf, fx, fy, cmold, mold
    
    integer           :: i,j,l,le,lee,ls,lss,n,lmr,i1,i2,j1,j2, is, ie, js, je
    integer           :: halo
    real              :: max_one
    real              :: gam, gamma, my_gamma
    integer,parameter :: iter_limit=15
    integer           :: n_advz, iter, nglob, status
    real              :: l_xi, g_xi, l_nxi, g_nxi, g_nl
!METHOD#2 #ifdef MPI    
!METHOD#2     integer :: request, stat(MPI_STATUS_SIZE)
!METHOD#2     logical :: done
!METHOD#2 #endif
    
    !----------- START
    
    ! Indices & work arrays
    CALL GET_DISTGRID( dgrid(1), I_STRT=i1, I_STOP=i2,J_STRT=j1, J_STOP=j2)
    lmr=lm(1)

    is=i1 ; ie=i2 ; js=j1 ; je=j2
    
    halo = wind_dat(1)%halo
    allocate(cmold( i1-halo:i2+halo, j1-halo:j2+halo, 0:lmr+1) ) ! same halo as wind_dat)

    halo = m_dat(1)%halo
    allocate( mold( i1-halo:i2+halo, j1-halo:j2+halo, lmr    ) ) ! same halo as m_dat)

    rm  => mass_dat(1)%rm
    rxm => mass_dat(1)%rxm
    rym => mass_dat(1)%rym
    rzm => mass_dat(1)%rzm
    m   =>    m_dat(1)%data
    cm  => wind_dat(1)%cm

    
    ls=1
    le=lmr
    
    if(abs(zbeg(1)-zbeg(1))<epsz) ls = 1    ! bottom
    if(abs(zend(1)-zend(1))<epsz) le = lmr  ! top

    !
    !     if requested limit vertical slopes such that no non-negative 
    !     tracer masses should occur
    if (limits) then
       do n = 1, ntracet
       do l = ls, le
       do j = js, je
       do i = is, ie
          rzm(i,j,l,n) = max (min (rzm(i,j,l,n), rm(i,j,l,n)), -rm(i,j,l,n))
       end do
       end do
       end do
       end do
    end if
    !    
    !====  loop to determine number of needed iterations    
    !  
    lss = ls; lee = le    
    if (abs (zbeg(1) - zbeg(1)) < epsz) lss = 2
    if (abs (zend(1) - zend(1)) < epsz) lee = lmr - 1

    cmold  = cm
    mold   = m
    n_advz = 1
    l_xi   = 0.0
    l_nxi  = 0.0

    cfl: do

       ! calculate new air mass distribution
       cm = cmold / float (n_advz)
       m = mold

       ! do the iteration
       advz:   do iter=1,n_advz
          
          ! reset gamma
          gam = 0.

!METHOD #2 #ifdef MPI
!METHOD #2           call MPI_IRECV(gam, 1, my_real, MPI_ANY_SOURCE, iter, localComm, request, ierr)
!METHOD #2           done = .false.
!METHOD #2 #endif

          ! compute f, pf, fx and fy
          ijl: do l=lss-1,lee      ! 1,lmm1
          do j=j1,j2
          do i=i1,i2

             if(cm(i,j,l)>=zero) then
                my_gamma=cm(i,j,l)/m(i,j,l)
             else
                my_gamma=cm(i,j,l)/m(i,j,l+1)
             end if

             !aMETHOD-1
             l_xi=max(l_xi,abs(my_gamma))             
             if ( abs(my_gamma) >= one )  exit ijl
             !eMETHOD-1
             
             !aMETHOD-2
!              gam  = abs(my_gamma)
!              l_xi = max(l_xi, gam)

!              if ( gam >= one ) then
!                 do n=0,npes-1
!                    call MPI_SEND(gam, 1, my_real, n, iter, localComm, stat, ierr)
!                 end do
!              end if
!              call MPI_TEST(request, done, stat, ierr)
!              if (done) exit ijl
             !eMETHOD-2
             
          end do
          end do
          end do ijl


!aMETHOD-1
          my_gamma=abs(my_gamma)
          call Par_REDUCE( my_gamma, 'MAX', gamma, status, all=.true.)

          if (gamma>=one) then
             l_nxi=l_nxi+1
             gam=gamma             
             exit advz
          end if
!eMETHOD-1
          
          !aMETHOD #2
!           call Par_REDUCE( gam, 'MAX', gam, status, all=.true.)
!           call MPI_TEST(request, done, stat, ierr)
!           if (.not.done) then
! !            print*, "cancel ", iter, request, gam
!              call MPI_CANCEL(request, ierr)
!              call MPI_REQUEST_FREE(request, ierr)
!           end if
          
!           if (gam>=one) then
!              l_nxi=l_nxi+1
!              exit advz
!           end if
          !eMETHOD-2
          
          ! no CFLs upto now.....possibly in next?

          if (abs(zbeg(1)-zbeg(1))<epsz) then
             ! bottom: one-sided update of the mass
             m(is:ie,js:je,1)   = m(is:ie,js:je,1)   - cm(is:ie,js:je,1)
          end if
          
          if (abs(zend(1)-zend(1))<epsz) then
             ! top:    one-sided update of the mass
             m(is:ie,js:je,lmr) = m(is:ie,js:je,lmr) + cm(is:ie,js:je,lmr-1)
          end if
          
          m(is:ie,js:je,lss:lee) =  m(is:ie,js:je,lss:lee) + &
                                   cm(is:ie,js:je,lss-1:lee-1)  - & 
                                   cm(is:ie,js:je,lss:lee)

       end do advz

       ! if still CFL violation increase iteration counter and re-start cfl loop
       ! escape is called when the number of iterations becomes too large
       if(abs(gam)>=one) then
          if(okdebug)then
             if ( my_gamma >= one) then
                write(gol,*)'dynamw: CFL violation: gamma,m,cm=', gam,m(i,j,l), m(i,j,l+1),cm(i,j,l),  &
                     '    in region,i,j,l,n: ',1,i,j,l,n ; call goPr
                write(gol,*)'dynamw: iterations:',n_advz ; call goPr
             endif
          endif
          n_advz = n_advz + 1

          if ( n_advz > iter_limit ) &
               call escape_tm('dynamw: too many CFL violations!')
          cycle cfl
       else
          exit cfl
       end if

    end do cfl

    
    ! n_advz is now determined : gather and store information
    CALL PAR_REDUCE( float(n_advz), 'max', g_nl,  status, .true.)
    CALL PAR_REDUCE( l_xi,          'max', g_xi,  status, .true.)
    CALL PAR_REDUCE( l_nxi,         'sum', g_nxi, status, .true.)

    nloop_max(1,3) = max(nloop_max(1,3), nint (g_nl))
    xi(1,3)        = max (xi(1,3), g_xi)
    nxi(1,3)       = nxi(1,3) + nint (g_nxi)

    ! reset m and cm to original values:
    cm = cmold/float(n_advz)
    m = mold    

    !
    !==== Apply number of iterations to calculate new air mass distribution
    !
    allocate (mnew( is:ie, js:je, 1:lm(1) )) 
    allocate (   f( is:ie, js:je, 0:lm(1) ))
    allocate (  pf( is:ie, js:je, 0:lm(1) ))
    allocate (  fx( is:ie, js:je, 0:lm(1) ))
    allocate (  fy( is:ie, js:je, 0:lm(1) ))

    ! do the iteration
    ITERA: do iter=1,n_advz

       if (abs(zbeg(1)-zbeg(1))<epsz) then
          ! bottom: one-sided update of the mass
          mnew(is:ie,js:je,1)   = m(is:ie,js:je,1)   - cm(is:ie,js:je,1)
       end if
       
       if (abs(zend(1)-zend(1))<epsz) then
          ! top:    one-sided update of the mass
          mnew(is:ie,js:je,lmr) = m(is:ie,js:je,lmr) + cm(is:ie,js:je,lmr-1)
       end if

       mnew(is:ie,js:je,lss:lee) =  m(is:ie,js:je,lss:lee)     &
                                 + cm(is:ie,js:je,lss-1:lee-1) & 
                                 - cm(is:ie,js:je,lss:lee)

       ! loop over tracers
       TRAC : do n=1,ntracet

          ! compute f, pf, fx and fy
          do l=lss-1,lee      ! 1,lmm1
          do j=js,je
          do i=is,ie
             if(cm(i,j,l)>=zero) then
                gamma=cm(i,j,l)/m(i,j,l)
                f(i,j,l)=gamma*(rm(i,j,l,n)+(one-gamma)* rzm(i,j,l,n) )
                pf(i,j,l)=cm(i,j,l)*(gamma*gamma* rzm(i,j,l,n)-3.*f(i,j,l))
                fx(i,j,l)=gamma*rxm(i,j,l,n)
                fy(i,j,l)=gamma*rym(i,j,l,n)
             else
                gamma=cm(i,j,l)/m(i,j,l+1)
                f(i,j,l)=gamma*(rm(i,j,l+1,n)-(one+gamma)* rzm(i,j,l+1,n))
                pf(i,j,l)=cm(i,j,l)*(gamma*gamma* rzm(i,j,l+1,n)-3.*f(i,j,l))
                fx(i,j,l)=gamma*rxm(i,j,l+1,n)
                fy(i,j,l)=gamma*rym(i,j,l+1,n)
             end if
          end do
          end do
          end do
          
          ! calculate new tracer mass, and tracer mass slopes
          ! update rm, rzm, rxm and rym in interior layers of the column

#ifdef with_budgets    
          ! calculate flux flowing in from below:
          if (apply_budget_global_flux) then
             do j=js,je
                do i=is,ie
                   budget_flux(1)%flux_z1(i,j,n) = budget_flux(1)%flux_z1(i,j,n) +   &
                        f(i,j,lflux1(1)-1)*1e3/ra(n)  ! moles
                   budget_flux(1)%flux_z2(i,j,n) = budget_flux(1)%flux_z2(i,j,n) +   &
                        f(i,j,lflux2(1)-1)*1e3/ra(n)  ! moles
                enddo
             enddo
          endif
#endif
          do l=lss,lee      !2,lmm1
             rm(is:ie,js:je,l,n) = rm(is:ie,js:je,l,n) + f(is:ie,js:je,l-1)-f(is:ie,js:je,l)
             rzm(is:ie,js:je,l,n) = rzm(is:ie,js:je,l,n) + &
                  ( pf(is:ie,js:je,l-1) - pf(is:ie,js:je,l) - (cm(is:ie,js:je,l-1)-cm(is:ie,js:je,l))* rzm(is:ie,js:je,l,n) &
                  + 3.0*( (cm(is:ie,js:je,l-1)+cm(is:ie,js:je,l))* rm(is:ie,js:je,l,n) &
                  - (f(is:ie,js:je,l-1)+f(is:ie,js:je,l))* m(is:ie,js:je,l) ) ) / mnew(is:ie,js:je,l)
             if ( limits ) then    !CMK added
                rzm(is:ie,js:je,l,n) = max(min(rzm(is:ie,js:je,l,n), rm(is:ie,js:je,l,n)), -rm(is:ie,js:je,l,n))
             end if
             rxm(is:ie,js:je,l,n) = rxm(is:ie,js:je,l,n) + (fx(is:ie,js:je,l-1)-fx(is:ie,js:je,l))
             rym(is:ie,js:je,l,n) = rym(is:ie,js:je,l,n) + (fy(is:ie,js:je,l-1)-fy(is:ie,js:je,l))
          end do
          
          if ( abs(zbeg(1)-zbeg(1)) < epsz ) then ! bottom
             ! update rm, rzm, rxm and rym near the ground surface
             rm(is:ie,js:je,1,n)=(rm(is:ie,js:je,1,n)-f(is:ie,js:je,1))
             rzm(is:ie,js:je,1,n)=rzm(is:ie,js:je,1,n)+(-pf(is:ie,js:je,1) & 
                  +cm(is:ie,js:je,1)*rzm(is:ie,js:je,1,n) &
                  +3.*(cm(is:ie,js:je,1)*rm(is:ie,js:je,1,n) &
                  -f(is:ie,js:je,1)*m(is:ie,js:je,1)))/mnew(is:ie,js:je,1)
             if ( limits ) then   !CMK added
                rzm(is:ie,js:je,1,n) = max(min(rzm(is:ie,js:je,1,n),  &
                     rm(is:ie,js:je,1,n)),  &
                     -rm(is:ie,js:je,1,n))
             end if
             rxm(is:ie,js:je,1,n)=rxm(is:ie,js:je,1,n)-fx(is:ie,js:je,1)
             rym(is:ie,js:je,1,n)=rym(is:ie,js:je,1,n)-fy(is:ie,js:je,1)
          end if

          if (abs(zend(1)-zend(1))<epsz) then ! top
             ! update rm, rzm, rxm and rym at the top of the atmosphere
             rm(is:ie,js:je,lmr,n)=rm(is:ie,js:je,lmr,n)+f(is:ie,js:je,lmr-1)
             rzm(is:ie,js:je,lmr,n)=rzm(is:ie,js:je,lmr,n)+&
                  (pf(is:ie,js:je,lmr-1)    &
                  -cm(is:ie,js:je,lmr-1)*rzm(is:ie,js:je,lmr,n) &
                  +3.*(cm(is:ie,js:je,lmr-1)*rm(is:ie,js:je,lmr,n) &
                  -f(is:ie,js:je,lmr-1)*m(is:ie,js:je,lmr))) / &
                  mnew(is:ie,js:je,lmr)
             if ( limits ) then   !CMK added
                rzm(is:ie,js:je,lmr,n) = max(min(rzm(is:ie,js:je,lmr,n),  &
                     rm(is:ie,js:je,lmr,n)),  &
                     -rm(is:ie,js:je,lmr,n))
             end if
             rxm(is:ie,js:je,lmr,n) = rxm(is:ie,js:je,lmr,n) + &
                  fx(is:ie,js:je,lmr-1)
             rym(is:ie,js:je,lmr,n) = rym(is:ie,js:je,lmr,n) + &
                  fy(is:ie,js:je,lmr-1)
          end if

       end do TRAC ! loop over tracers

       ! store new air mass in m array
       m(is:ie,js:je,ls:le)=mnew(is:ie,js:je,ls:le)

    end do ITERA !iter

    ! restore cm wind
    cm = cmold

    deallocate (f)
    deallocate (pf)
    deallocate (fx)
    deallocate (fy)
    deallocate(mnew)
    deallocate(cmold)

    if ( okdebug ) then

       mold(is:ie,js:je,ls:le)=rm(is:ie,js:je,ls:le,1)/m(is:ie,js:je,ls:le)
       
       call REDUCE( dgrid(1), mold, 2, 'MAX', max_one, status)
       
       if (isRoot) then
          write(gol,*) 'dynamw: z Maximum value mixing ratio',max_one; call goPr
       end if
    end if

    deallocate(mold)
    nullify(rm)
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
    nullify(m)
    nullify(cm)

  END SUBROUTINE DYNAMW
  !EOC

END MODULE ADVECTZ