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

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

  ! --- interface

  private

  public :: advectzzoom



contains


  subroutine advectzzoom(region)
    !
    ! set parameters for advectz
    ! written by patrick berkvens and mike botchev, march-june 1999
    !
    use dims,        only : im, jm, lm, xref, yref, zref, tref
    use dims,        only : parent, zoom2D, touch_sp, touch_np, xcyc
    use dims,        only : splitorderzoom, n_operators, status, nsplitsteps
    use toolbox,     only : escape_tm
#ifdef MPI
    use mpi_const,   only : my_real, com_trac, ierr, mpi_sum
    use mpi_comm,    only : barrier_t
#endif
    use partools  , only : myid, root_t, ntracetloc

    implicit none

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

    ! local
    integer            :: is,ie,js,je,n,q
    integer            :: imr,jmr,lmr,tref_,xref_,yref_,zref_
    logical            :: z_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

    ! write cells along the z-walls of each child 
    ! to rm(0,...),rm(imr+1,...) and to interface cells of the child
    !
    ! this option should be there: zooming in z disabled....
    if ( .not. zoom2D ) then 
       call escape_tm('advectzzoom: ERROR zoom2D should be true:'// &
            ' zooming along z not allowed')
       !CMKCALL put_zedges(region,rm,rxm,rym,rzm,m,  &
       !CMK                        rmgl,rxmgl,rymgl,rzmgl,mgl)
    end if

    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)

    ! determine the scope for advectz: 

    if ( region == 1 ) then
       xref_ = 0; yref_ = 0  ! to have is/ie and js/je 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
    z_encountered=.false.

    do n=1,n_operators          ! now track following steps from q  

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

       select case(dir)
       case('x') 
          if ( (.not.z_encountered) .or. (xcyc(region) == 1) ) then
             is=1                          ! x-substep is before z => 
             ie=im(region)                 ! full i-scope 
          else
             is=xref_+1                    ! z-substep is before x =>
             ie=im(region)-xref_           ! restricted i-scope
          end if
       case('y') 
          if (.not.z_encountered) then
             js=1                          ! y-substep is before z => 
             je=jm(region)                 ! full j-scope 
          else
             if(touch_sp(region) == 1 ) then
                js=1                       ! sp is southern edge--->full scope.
             else
                js=yref_+1                 ! x-substep is before y =>
             end if
             if ( touch_np(region) == 1 ) then
                je=jm(region)              ! np is northern edge--->full scope
             else
                je=jm(region)-yref_        ! restricted j-scope
             end if
          end if
       case('z') 
          z_encountered=.true.
       case ('c') 
       case ('v') 
       case ('s') 
       case default
          print *,'advectzzoom: strange value in splitorderzoom(',region,',',  &
               q*(nsplitsteps/2)*tref_+n,'): ',q
          call escape_tm('advectzzoom: Error in advectz')
       end select

    end do

    call advectz_work(region,is,ie,js,je)

  end subroutine advectzzoom



  subroutine advectz_work(region,is,ie,js,je)
    !
    ! calls dynamw
    ! written by mike botchev, march-june 1999
    ! modified by MK, dec 2002
    !

    ! in/out
    integer,intent(in) :: region
    integer,intent(in) :: is
    integer,intent(in) :: ie
    integer,intent(in) :: js
    integer,intent(in) :: je

    ! start

#ifndef slopes
         stop 'advectz: Wrong advection scheme: adv_scheme/=slope'
#endif
    call dynamw(region,is,ie,js,je)
    
  end subroutine advectz_work



  subroutine dynamw(region,is,ie,js,je)
    !-----------------------------------------------------------------------
    !
    !****   dynamw          - vertical tracer transport     v 9.1
    !
    ! programmed by         mh      mpi HH          23-feb-1995
    !
    !       purpose
    !       -------
    !       calculate amount of tracer moved in a vertical advection
    !       substep
    !
    !       interface
    !       ---------
    !       call dynamw
    !
    !       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
    !       changed order of loops for increased performance
    !                               mh, 11-jun-1994
    !     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 binas,       only : xgamma=>gamma
    use dims,        only : im, jm, lm, xref, yref, zref, tref
    use dims,        only : zbeg, zend, epsz, nloop_max, zero, one
    use dims,        only : okdebug, nregions, parent, limits, xi, nxi
    use global_data, only : wind_dat, mass_dat
    use meteodata  , only : m_dat
    use chem_param,  only : ntracet, ra
    use zoom_tools,  only : mix_edges
    use toolbox,     only : escape_tm
#ifdef with_budgets    
    use budget_global, only : budget_flux, lflux1, lflux2, apply_budget_global_flux
#endif
    use redgridZoom, only : nred, clustsize, jred, imredj
#ifdef MPI
    use mpi_const,   only : my_real,com_trac,ierr, mpi_max, pe_first_tracer
    use mpi_comm,    only : barrier_t
#endif
    use partools   , only : myid, root_t
    use partools   , only : ntracetloc, ntracet_ar
    use omp_partools, only: my_omp_domain, my_omp_max

    implicit none

    ! input,output
    integer,intent(in) :: region
    integer,intent(in) :: is
    integer,intent(in) :: ie
    integer,intent(in) :: js
    integer,intent(in) :: je

    ! local
    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,imr,jmr,lmr
    integer                           :: lrg
    integer                           :: tref_,xref_,yref_,zref_
    real                              :: max_one,max_all, summ
    real                              :: mxval, gam, gamma
    integer,parameter                 :: iter_limit=15
    integer                           :: n_advz, iter, nglob, offsetn, ixs, ixe, redfact

    ! Omp parameters
    integer                           :: ns, ne, jss, jee, jje, jjs

    ! start

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

    !CMK next line: dec2003 reported by FD
    allocate(cmold(0:im(region)+1,0:jm(region)+1,0:lm(region)+1))
    allocate(mold(-1:im(region)+2,-1:jm(region)+2,lm(region)))

    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
    cm => wind_dat(region)%cm_t

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

    ! prepare z-edges to allow for uniform work:
    !CMK not any more allowed :
    !  call compress_zedges(region,is,ie,js,je,m,rm,rxm,rym,rzm)

    ! 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)
    ! check is,ie,js,je:
    if ( ( is /= xref_+1 ) .and. ( is /= 1 ) ) &
         call escape_tm('dynamw: Wrong value for IS in dynamw')
    if ( ( ie /= imr-xref_ ) .and. ( ie /= imr ) ) &
         call escape_tm('dynamw: Wrong value for IE in dynamw')
    if ( ( js /= yref_+1 ) .and. ( js /= 1 ) ) &
         call escape_tm('dynamw: Wrong value for JS in dynamw')
    if ( ( je /= jmr-yref_ ) .and. ( je /= jmr ) ) &
         call escape_tm('dynamw: Wrong value for JE in dynamw')

    ! compute ls/le -- cells is:ie,js:je:ls:le will be update
    ls = zref_; le = lmr-zref_+1

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

    !
    !     if requested limit vertical slopes such that no non-negative 
    !     tracer masses should occur
    if (limits) then
!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (is, ie, js, je, ls, le, rm, rzm) &
!$OMP  private (i, j, l, n, ns, ne)
      call my_omp_domain (1, ubound (rzm, 4), ns, ne)
      do n = ns, ne
         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
!$OMP END PARALLEL
    end if
    !    
    !     calculate new air mass distribution
    !  
    lss = ls; lee = le    
    if (abs (zbeg(region) - zbeg(1)) < epsz) lss = 2
    if (abs (zend(region) - zend(1)) < epsz) lee = lmr - 1

!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (cm, cmold, m, mold, jmr) &
!$OMP  private (jjs,jje)
     call my_omp_domain (0, jmr + 1, jjs, jje)
     cmold(:,jjs:jje,:) = cm(:,jjs:jje,:)

     call my_omp_domain (-1, jmr + 2, jjs, jje)
     mold(:,jjs:jje,:) = m(:,jjs:jje,:)
!$OMP END PARALLEL

    n_advz=1

    ! trial loop to determine # iterations needed:
    
    cfl: do

       !     calculate new air mass distribution

!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (cm, cmold, m, mold, n_advz,jmr) &
!$OMP  private (jjs,jje)
       call my_omp_domain (0, jmr + 1, jjs, jje)
       cm(:,jjs:jje,:) = cmold(:,jjs:jje,:) / float (n_advz)

       call my_omp_domain (-1, jmr + 2, jjs, jje)
       m(:,jjs:jje,:) = mold(:,jjs:jje,:)
!$OMP END PARALLEL

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

          ! reset gamma
          gam = 0.

          ! compute f, pf, fx and fy
          do l=lss-1,lee      ! 1,lmm1
             do j=1,jmr
                do i=1,imr
                   if(cm(i,j,l)>=zero) then
                      gamma=cm(i,j,l)/m(i,j,l)
                   else
                      gamma=cm(i,j,l)/m(i,j,l+1)
                   end if
                   xi(region,3)=max(xi(region,3),abs(gamma))
                   if ( abs(gamma) >= one ) then
                      nxi(region,3)=nxi(region,3)+1
                      gam=gamma
                      exit advz
                   end if
                end do
             end do
          end do

          ! no CFLs upto now.....possibly in next?

!$OMP PARALLEL &
!$OMP  default (none) &
!$OMP  shared  (region, is, ie, js, je, lmr, lss, lee, m, &
!$OMP           cm) &
!$OMP  private (jss, jee)
          call my_omp_domain (js, je, jss, jee)

          if (abs(zbeg(region)-zbeg(1))<epsz) then
             ! bottom: one-sided update of the mass
             m(is:ie,jss:jee,1)   = m(is:ie,jss:jee,1)   - cm(is:ie,jss:jee,1)
          end if
          if (abs(zend(region)-zend(1))<epsz) then
             ! top:    one-sided update of the mass
             m(is:ie,jss:jee,lmr) = m(is:ie,jss:jee,lmr) + cm(is:ie,jss:jee,lmr-1)
          end if
          m(is:ie,jss:jee,lss:lee) = m(is:ie,jss:jee,lss:lee) + &
               cm(is:ie,jss:jee,lss-1:lee-1)  - & 
               cm(is:ie,jss:jee,lss:lee)
!$OMP END PARALLEL

       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
             print *,'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: ',region,i,j,l,n
             print*,'dynamw: iterations:',n_advz
          end if
          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!______CFL

    ! n_advz is now determined: apply this number of iterations:
    ! calculate new air mass distribution

    nloop_max(region,3) = max(nloop_max(region,3), n_advz)

    ! reset m and cm to original values:

!$OMP PARALLEL &
!$OMP  default (none) &
#if defined (with_budgets)
!$OMP  shared  (apply_budget_global_flux, budget_flux, lflux1, lflux2,ra) &
#endif
!$OMP  shared  (region, limits, n_advz, is, ie, imr, js, je, &
!$OMP           jmr, le, ls, lmr, lss, lee, m, mold, cm, cmold, rm, &
!$OMP           rzm, rxm, rym, xi, ntracetloc, ntracet_ar, myid) &
!$OMP  private (i, j, l, n, iter, jjs, jje, jss, jee, gamma, &
!$OMP           mnew, f, pf, fx, fy, nglob)
    call my_omp_domain (0, jmr + 1, jjs, jje)
    cm(:,jjs:jje,:) = cmold(:,jjs:jje,:) / float (n_advz)

    call my_omp_domain (-1, jmr + 2, jjs, jje)
    m(:,jjs:jje,:) = mold(:,jjs:jje,:)

!$OMP BARRIER

    call my_omp_domain (js, je, jss, jee)

    allocate (mnew(is:ie,jss:jee,lm(region))) 
    allocate (f(is:ie,jss:jee,0:lm(region)))
    allocate (pf(is:ie,jss:jee,0:lm(region)))
    allocate (fx(is:ie,jss:jee,0:lm(region)))
    allocate (fy(is:ie,jss:jee,0:lm(region)))

    ! do the iteration
    do iter=1,n_advz

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

!$OMP BARRIER

!   0:          q1  1 3 1 3510489980.09739351
!ERROR: 0031-250  task 0: Trace/BPT trap

       !    
       !     loop over tracers
       !    
       do n=1,ntracetloc  !WP! changed scope

          ! compute f, pf, fx and fy
          do l=lss-1,lee      ! 1,lmm1
             do j=jss,jee
                !do i=1,imr
                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

!$OMP BARRIER

          ! 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
             nglob = n + sum(ntracet_ar(0:myid-1))
             do j=jss,jee
                do i=is,ie
                   budget_flux(region)%flux_z1(i,j,nglob) = budget_flux(region)%flux_z1(i,j,nglob) +   &
                       f(i,j,lflux1(region)-1)*1e3/ra(nglob)  ! moles
                   budget_flux(region)%flux_z2(i,j,nglob) = budget_flux(region)%flux_z2(i,j,nglob) +   &
                       f(i,j,lflux2(region)-1)*1e3/ra(nglob)  ! moles
                enddo
             enddo
          endif
#endif
          do l=lss,lee      !2,lmm1
             rm(is:ie,jss:jee,l,n) = rm(is:ie,jss:jee,l,n) + &
                  f(is:ie,jss:jee,l-1)-f(is:ie,jss:jee,l)
             rzm(is:ie,jss:jee,l,n) = rzm(is:ie,jss:jee,l,n) + &
                  ( pf(is:ie,jss:jee,l-1) - pf(is:ie,jss:jee,l) &
                    - (cm(is:ie,jss:jee,l-1)-cm(is:ie,jss:jee,l))* &
                            rzm(is:ie,jss:jee,l,n)               &
                    + 3.0*( (cm(is:ie,jss:jee,l-1)+cm(is:ie,jss:jee,l))* &
                                      rm(is:ie,jss:jee,l,n)            &
                            - (f(is:ie,jss:jee,l-1)+f(is:ie,jss:jee,l))* &
                                  m(is:ie,jss:jee,l) ) ) / mnew(is:ie,jss:jee,l)
             if ( limits ) then    !CMK added
                rzm(is:ie,jss:jee,l,n) = max(min(rzm(is:ie,jss:jee,l,n),  &
                     rm(is:ie,jss:jee,l,n)),  &
                     -rm(is:ie,jss:jee,l,n))
             end if
             rxm(is:ie,jss:jee,l,n)=rxm(is:ie,jss:jee,l,n)+ &
                  (fx(is:ie,jss:jee,l-1)-fx(is:ie,jss:jee,l))
             rym(is:ie,jss:jee,l,n)=rym(is:ie,jss:jee,l,n)+ &
                  (fy(is:ie,jss:jee,l-1)-fy(is:ie,jss:jee,l))
          end do

          if ( abs(zbeg(region)-zbeg(1)) < epsz ) then ! bottom
             ! update rm, rzm, rxm and rym near the ground surface
             rm(is:ie,jss:jee,1,n)=(rm(is:ie,jss:jee,1,n)-f(is:ie,jss:jee,1))
             rzm(is:ie,jss:jee,1,n)=rzm(is:ie,jss:jee,1,n)+(-pf(is:ie,jss:jee,1) & 
                  +cm(is:ie,jss:jee,1)*rzm(is:ie,jss:jee,1,n) &
                  +3.*(cm(is:ie,jss:jee,1)*rm(is:ie,jss:jee,1,n) &
                  -f(is:ie,jss:jee,1)*m(is:ie,jss:jee,1)))/mnew(is:ie,jss:jee,1)
             if ( limits ) then   !CMK added
                rzm(is:ie,jss:jee,1,n) = max(min(rzm(is:ie,jss:jee,1,n),  &
                     rm(is:ie,jss:jee,1,n)),  &
                     -rm(is:ie,jss:jee,1,n))
             end if
             rxm(is:ie,jss:jee,1,n)=rxm(is:ie,jss:jee,1,n)-fx(is:ie,jss:jee,1)
             rym(is:ie,jss:jee,1,n)=rym(is:ie,jss:jee,1,n)-fy(is:ie,jss:jee,1)
          end if

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

       end do ! loop over tracers

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

    end do !iter

!$OMP BARRIER

    call my_omp_domain (0, jmr + 1, jjs, jje)
    cm(:,jjs:jje,:) = cmold(:,jjs:jje,:)

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

!$OMP END PARALLEL

    deallocate(cmold)
    deallocate(mold)


    if ( okdebug .and. myid == root_t ) then
       max_one=0.0
       max_one=maxval(rm(1:imr,1:jmr,1:lmr,1)/m(1:imr,1:jmr,1:lmr))
       print *, 'dynamw: z Maximum value mixing ratio',max_one
    end if
    
    nullify(rm)
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
    nullify(m)
    nullify(cm)
    
    call mix_edges(region)

  end subroutine dynamw


end module advectz
