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

  use GO, only : gol, goPr, goErr
  
use go, only : gol, gobug
use partools, only : tracer_active
use tracer_data, only : tracer_print
use meteodata, only : m_dat

  implicit none

  ! --- interface ---

  private

  public :: mix_edges, coarsen_region, update_parent


contains



  subroutine mix_edges(region)
    !
    ! distribute data at the x-edges of the zoom region
    ! over the whole interface cell
    ! written by mike botchev, march-june 1999
    ! adapted for the CRAY by Maarten Krol, 3-2000
    !
    use dims
    use global_data, only : mass_dat, region_dat
    use chem_param,  only : ntracet
    use toolbox,     only : escape_tm
    use partools   , only : myid
    use partools   , only : ntracetloc
    use partools   , only : lmloc, offsetl
    use partools   , only : which_par, previous_par
    use meteodata  , only : m_dat

    ! input
    integer,intent(in)            ::  region

    ! local
    real,dimension(:,:,:,:),pointer    :: rm
#ifdef slopes
    real,dimension(:,:,:,:),pointer    :: rxm,rym,rzm
#endif
    real,dimension(:,:,:),  pointer    :: m
    integer,dimension(:,:), pointer    :: edge

    integer    :: imr, i, j, l, n, jmr, lmr, ip, jp, imp, jmp, nt
    integer    :: xref_, yref_, zref_, xyzref_, xyref_
    real       :: m_edge, rm_edge, rxm_edge, rym_edge, rzm_edge
    logical    :: xedges , ynp,ysp

    ! relaxation factor for slopes; relax=0.0 means upwinding
    real,parameter:: relax=1.0
    real,dimension(:,:),allocatable   :: to_parent
#ifdef slopes
    real,dimension(:,:),allocatable   :: to_parentx
    real,dimension(:,:),allocatable   :: to_parenty
    real,dimension(:,:),allocatable   :: to_parentz
#endif

    ! start

    xedges = .true.
    ysp = .true.  !mix also NP/SP?
    ynp = .true.

    which_par=previous_par(region)
    if ( which_par == 'tracer' .and. ntracetloc == 0 ) return 
    if ( which_par == 'levels' .and. lmloc == 0 ) return  !WP! 

    imr = im(region)
    jmr = jm(region)

    m => m_dat(region)%data
    
    if ( which_par == 'tracer') then

       rm => mass_dat(region)%rm_t
#ifdef slopes
       rxm => mass_dat(region)%rxm_t
       rym => mass_dat(region)%rym_t
       rzm => mass_dat(region)%rzm_t
#endif

       lmr = lm(region)
       nt=ntracetloc

    else if ( which_par == 'levels' ) then

       rm => mass_dat(region)%rm_k
#ifdef slopes
       rxm => mass_dat(region)%rxm_k
       rym => mass_dat(region)%rym_k
       rzm => mass_dat(region)%rzm_k
#endif

       lmr = lmloc
       nt=ntracet
    end if

    edge => region_dat(region)%edge

    xref_ = xref(region)/xref(parent(region))
    yref_ = yref(region)/yref(parent(region))
    zref_ = zref(region)/zref(parent(region))
    if ( zref_ /= 1 ) call escape_tm( &
         'mix_edges: zooming in z direction not supported!')

    xyref_ = xref_*yref_

    if ( okdebug ) then
       write (gol,*) 'mix_edges: region=',region; call goPr
    end if
    if ( (xref_ == 1) .or. (im(region)/xref(region) == im(1)) ) then
       if ( okdebug ) then
         write (gol,*) 'mix_edges: no refinement or periodic bc, skipping x-walls'; call goPr
       end if
       xedges = .false. 
       if ( region == 1 ) then
          if ( okdebug ) then
            write (gol,*) 'mix_edges: region = 1, returning'; call goPr
          end if
          return
       end if
    end if

    if ( touch_sp(region) == 1) then
       ysp = .false.
       if ( okdebug ) then 
         write (gol,*) 'mix_edges: SP touching...skip SP y-walls'; call goPr
       end if
    end if
    if(touch_np(region) == 1) then
       ynp = .false.
       if ( okdebug ) then
         write (gol,*) 'mix_edges: NP touching...skip NP y-walls'; call goPr
       end if
    end if
    imp = imr/xref_   !in 'coarse' resolution
    jmp = jmr/yref_

    ! ~~ mass (all levels!)

    allocate(to_parent(imp,jmp))

    !do l=1,lmr !WP! depends on tracers or levels
    do l=1,lm(region)
       to_parent = 0.0
       do j=1,jmr
          jp = 1 + (j-1)/yref_
          do i=1,imr
             if ( edge(i,j) /= -1 ) cycle   !do only for edge of child...
             ip = 1 + (i-1)/xref_
             to_parent(ip,jp) = to_parent(ip,jp) + m(i,j,l)/xyref_
          end do
       end do
       if ( xedges ) then     !mix the x-wall of the zoom region.
          do j = 1,jmr
             jp = 1 + (j-1)/yref_
             do i = 1, xref_
                m(i,j,l) = to_parent(1,jp)
             end do
             do i = imr-xref_+1,imr
                m(i,j,l) = to_parent(imp,jp)
             end do
          end do
       end if    !xedges
       if ( ysp ) then
          do j = 1,yref_
             do i = 1, imr
                ip = 1 + (i-1)/xref_
                m(i,j,l) = to_parent(ip,1)
             end do
          end do
       end if
       if ( ynp ) then
          do j = jmr-yref_+1,jmr
             do i = 1, imr
                ip = 1 + (i-1)/xref_
                m(i,j,l) = to_parent(ip,jmp)
             end do
          end do
       end if
    end do  !l

    deallocate(to_parent)
    
    ! ~~ tracers

    allocate(to_parent(imp,jmp))
#ifdef slopes
    allocate(to_parentx(imp,jmp))
    allocate(to_parenty(imp,jmp))
    allocate(to_parentz(imp,jmp))
#endif

    do n=1,nt  !WP! depends on tracers or levels
       do l=1,lmr !WP! depends on tracers or levels
          to_parent = 0.0
#ifdef slopes
          to_parentx = 0.0
          to_parenty = 0.0
          to_parentz = 0.0
#endif
          do j=1,jmr
             jp = 1 + (j-1)/yref_
             do i=1,imr
                if ( edge(i,j) /= -1 ) cycle   !do only for edge of child...
                ip = 1 + (i-1)/xref_
                to_parent(ip,jp) = to_parent(ip,jp) + rm(i,j,l,n)/xyref_
#ifdef slopes
                to_parentx(ip,jp) = to_parentx(ip,jp) + rxm(i,j,l,n)/xyref_
                to_parenty(ip,jp) = to_parenty(ip,jp) + rym(i,j,l,n)/xyref_
                to_parentz(ip,jp) = to_parentz(ip,jp) + rzm(i,j,l,n)/xyref_
#endif
             end do
          end do
          if ( xedges ) then     !mix the x-wall of the zoom region.
             do j = 1,jmr
                jp = 1 + (j-1)/yref_
                do i = 1, xref_
                   rm(i,j,l,n) = to_parent(1,jp)
#ifdef slopes
                   rxm(i,j,l,n) = to_parentx(1,jp)
                   rym(i,j,l,n) = to_parenty(1,jp)
                   rzm(i,j,l,n) = to_parentz(1,jp)
#endif
                end do
                do i = imr-xref_+1,imr
                   rm(i,j,l,n) = to_parent(imp,jp)
#ifdef slopes
                   rxm(i,j,l,n) = to_parentx(imp,jp)
                   rym(i,j,l,n) = to_parenty(imp,jp)
                   rzm(i,j,l,n) = to_parentz(imp,jp)
#endif
                end do
             end do
          end if !xedges
          if ( ysp ) then
             do j = 1,yref_
                do i = 1, imr
                   ip = 1 + (i-1)/xref_
                   rm(i,j,l,n) = to_parent(ip,1)
#ifdef slopes
                   rxm(i,j,l,n) = to_parentx(ip,1)
                   rym(i,j,l,n) = to_parenty(ip,1)
                   rzm(i,j,l,n) = to_parentz(ip,1)
#endif
                end do
             end do
          end if
          if ( ynp ) then
             do j = jmr-yref_+1,jmr
                do i = 1, imr
                   ip = 1 + (i-1)/xref_
                   rm(i,j,l,n) = to_parent(ip,jmp)
#ifdef slopes
                   rxm(i,j,l,n) = to_parentx(ip,jmp)
                   rym(i,j,l,n) = to_parenty(ip,jmp)
                   rzm(i,j,l,n) = to_parentz(ip,jmp)
#endif
                end do
             end do
          end if
       end do   !l
    end do   !nt

    deallocate(to_parent)
#ifdef slopes
    deallocate(to_parentx)
    deallocate(to_parenty)
    deallocate(to_parentz)
#endif

    nullify(m)
    nullify(rm)
#ifdef slopes
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
#endif
    nullify(edge)

  end subroutine mix_edges



  recursive subroutine coarsen_region(region)
    !
    ! calculates values in the parent identical to the child. 
    ! Avoids small numerical differences that are caused by the 
    ! limited accuracy of the stored pressure and wind-data.
    ! written by Maarten Krol, december 1999
    !
    use dims
    use global_data,   only : wind_dat
    use meteodata    , only : m_dat
    use advect_tools,  only : m2phlb, m2phlb1
    use toolbox,       only : escape_tm
#ifdef MPI
    use mpi_comm, only : barrier,barrier_t
#endif
    use partools   , only : myid, root_t, root_k
    use partools   , only : ntracetloc, lmar, lmloc, ntracet_ar
    use tracer_data, only : Par_Check_Domain

    implicit none

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

    ! local:
    real,dimension(:,:,:),pointer     :: am,bm,cm,m,amp,bmp,cmp,mp
    real,dimension(:,:,:),allocatable :: top
    integer i,j,l,ip,jp,lp,regiop,imp,jmp,lmp,imr,jmr,lmr
    integer xref_,yref_,zref_,tref_
    real,dimension(0:im(region)+1,0:jm(region)+1,0:lm(region) ) :: fieldglob

    ! start

    regiop = parent(region) 
#ifdef MPI
    !WP! check whether region and
    call Par_Check_Domain( region, 'p', 'tracer' ) 
#endif
    !parents are on tracer domain
    imr = im(region)
    jmr = jm(region)
    lmr = lm(region)
    if ( regiop > 0 ) then 
#ifdef MPI 
       if ( ntracetloc /= 0 ) then !only PE's with nonzero ntracetloc continue
#endif
          if ( okdebug .and. myid == root_t ) then
            write (gol,*) 'coarsen_region: region:',region; call goPr
          end if
          xref_ = xref(region)/xref(regiop)
          yref_ = yref(region)/yref(regiop)
          zref_ = zref(region)/zref(regiop)
          tref_ = tref(region)/tref(regiop)
          imp = imr/xref_
          jmp = jmr/yref_
          lmp = lmr/zref_

          if ( ibeg(region) < iend(region) .and. &
               imp /= iend(region)-ibeg(region)+1 ) then
            write (gol,'("wrong region i dimensions:")'); call goErr
            write (gol,'("  ibeg, iend, imp : ",3i4)') ibeg(region), iend(region), imp; call goErr
            call escape_tm('coarsen_region: stop')
          end if
          if ( jmp /= jend(region)-jbeg(region)+1 ) then
            write (gol,'("wrong region j dimensions:")'); call goErr
            write (gol,'("  jbeg, jend, jmp : ",3i4)') jbeg(region), jend(region), jmp; call goErr
            call escape_tm('coarsen_region: stop')
          end if
          if ( lmp /= lend(region)-lbeg(region)+1 ) then
            write (gol,'("wrong region l dimensions:")'); call goErr
            write (gol,'("  lbeg, lend, lmp : ",3i4)') lbeg(region), lend(region), lmp; call goErr
            call escape_tm('coarsen_region: stop')
          end if
          write (gol,*) 'coarsen_region: coarsening parent',regiop,'with child',region; call goPr

          allocate(top(0:imp,jmp,lmp))
          amp => wind_dat(regiop)%am_t
          am => wind_dat(region)%am_t
          top = 0.0
          do l=1,lmr
             lp = 1 + (l-1)/zref_
             do j=1,jmr
                jp = 1 + (j-1)/yref_
                do ip=0,imp
                   top(ip,jp,lp) = top(ip,jp,lp) + am(ip*xref_,j,l)*tref_
                end do
             end do
          end do
          !   copy top to parent
          do l=1,lmp
             lp = lbeg(region)-1+l
             do j=1,jmp
                jp = jbeg(region)-1+j
                do i=0,imp
                   ip = mod(ibeg(region)-2+i,im(regiop))+1
                   amp(ip,jp,lp) = top(i,j,l)
                end do
             end do
          end do

          deallocate(top)
          nullify(am)
          nullify(amp)

          bm => wind_dat(region)%bm_t
          bmp => wind_dat(regiop)%bm_t
          allocate(top(imp,jmp+1,lmp))
          top = 0.0
          do l=1,lmr
             lp = 1 + (l-1)/zref_
             do i=1,imr
                ip = 1 + (i-1)/xref_
                do jp = 1,jmp+1
                   top(ip,jp,lp) = top(ip,jp,lp) + bm(i,1 + &
                        (jp-1)*yref_,l)*tref_
                end do
             end do
          end do
          do l=1,lmp
             lp = lbeg(region)-1+l
             do j=1,jmp+1
                jp = jbeg(region)-1+j
                do i=1,imp
                   ip = mod(ibeg(region)-2+i,im(regiop))+1
                   bmp(ip,jp,lp) = top(i,j,l)
                end do
             end do
          end do

          deallocate(top)
          nullify(bm)
          nullify(bmp)

          allocate(top(imp,jmp,lmp))
          cm => wind_dat(region)%cm_t
          cmp => wind_dat(regiop)%cm_t
          top = 0.0
          do i=1,imr
             ip = 1 + (i-1)/xref_
             do j = 1,jmr
                jp = 1 + (j-1)/yref_  
                do lp=1,lmp
                   top(ip,jp,lp) = top(ip,jp,lp) + cm(i,j,lp*zref_)*tref_
                end do
             end do
          end do
          do l=1,lmp
             lp = lbeg(region)-1+l
             do j=1,jmp
                jp = jbeg(region)-1+j
                do i=1,imp
                   ip = mod(ibeg(region)-2+i,im(regiop))+1
                   cmp(ip,jp,lp) = top(i,j,l)
                end do
             end do
          end do

          nullify(cm)
          nullify(cmp)

          top = 0.0

          m  => m_dat(region)%data
          mp => m_dat(regiop)%data

          do l=1,lmr
             lp = 1 + (l-1)/zref_
             do j = 1,jmr
                jp = 1 + (j-1)/yref_  
                do i=1,imr
                   ip = 1 + (i-1)/xref_
                   top(ip,jp,lp) = top(ip,jp,lp) + m(i,j,l)
                end do
             end do
          end do

          do l=1,lmp
             lp = lbeg(region)-1+l
             do j=1,jmp
                jp = jbeg(region)-1+j
                do i=1,imp
                   ip = mod(ibeg(region)-2+i,im(regiop))+1
                   mp(ip,jp,lp) = top(i,j,l)
                end do
             end do
          end do

          deallocate(top)
          nullify(m)
          nullify(mp)

#ifdef MPI 
       end if ! all PE's from here
#endif

       if ( regiop /= 1 ) then
          call m2phlb(regiop)
       else
          call m2phlb1(regiop)
       end if
       ! call grand-parent recursively...
       call coarsen_region(regiop)
#ifdef MPI
       call barrier   
#endif

    end if !regiop>0


  end subroutine coarsen_region



  subroutine update_parent(region)
    !
    ! update rm of its parent in correspondence with itself
    ! written by mike botchev, march-june 1999
    !
    use dims
    use global_data  ,only : mass_dat
    use meteodata   , only : m_dat
#ifdef with_budgets
    use budget_global,only : sum_update
#endif
    use chem_param,   only : ntracet
    use toolbox,      only : escape_tm
#ifdef MPI
    use mpi_const
#endif
    use partools   , only : which_par, previous_par
    use partools   , only : lmloc, offsetl
    use partools   , only : ntracetloc
    use tracer_data, only : Par_Check_Domain

    implicit none

    ! input
    integer,intent(in)   :: region

    ! local
    real,dimension(:,:,:,:),pointer          :: rm, rmp
#ifdef slopes
    real,dimension(:,:,:,:),pointer          :: rxm,rym,rzm, rxmp,rymp,rzmp
#endif
    real,dimension(:,:,:),pointer            :: m,mp
    real,dimension(:,:,:),allocatable        :: to_parent

    real,dimension(xref(region)/xref(parent(region)), &
         yref(region)/yref(parent(region)), &
         zref(region)/zref(parent(region)))   :: slope
    real    :: dummy
    real    :: sum_parent, sum_to_parent, sum_parent_all, sum_to_parent_all
    integer :: i, iend_loop, iip, ip, j, jp, l, lp, n
    integer :: lb, le, ib, ie, jb, je
    integer :: my_parent, xref_, yref_, zref_
    integer :: imp, jmp, lmp, imr, jmr, lmr, nt
    real    :: m_value, mxval, sum_mass
    integer,dimension(3) :: mxloc
    integer :: communicator, root_id
    integer :: lmpa

    ! start

    if (region == 1) return

    which_par=previous_par(region)

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

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

    imr = im(region)
    jmr = jm(region)
    my_parent = parent(region)
    xref_ = xref(region)/xref(my_parent)
    yref_ = yref(region)/yref(my_parent)
    zref_ = zref(region)/zref(my_parent)

    m  => m_dat(region   )%data
    mp => m_dat(my_parent)%data
       
    if ( which_par == 'tracer' ) then

       rm => mass_dat(region)%rm_t
#ifdef slopes
       rxm => mass_dat(region)%rxm_t
       rym => mass_dat(region)%rym_t
       rzm => mass_dat(region)%rzm_t
#endif

       rmp => mass_dat(my_parent)%rm_t
#ifdef slopes
       rxmp => mass_dat(my_parent)%rxm_t
       rymp => mass_dat(my_parent)%rym_t
       rzmp => mass_dat(my_parent)%rzm_t
#endif

       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
#ifdef slopes
       rxm => mass_dat(region)%rxm_k
       rym => mass_dat(region)%rym_k
       rzm => mass_dat(region)%rzm_k
#endif

       rmp => mass_dat(my_parent)%rm_k
#ifdef slopes
       rxmp => mass_dat(my_parent)%rxm_k
       rymp => mass_dat(my_parent)%rym_k
       rzmp => mass_dat(my_parent)%rzm_k
#endif

       lmr = lmloc
       nt=ntracet
#ifdef MPI  
       communicator=com_lev  !WP! assign com_lev as communicator
       root_id=root_k
#endif
    end if

    if ( okdebug ) then
       write (gol,*) 'update_parent: my_parent=',my_parent, &
            ' x-,y-,zref_: ',xref_,yref_,zref_ ; call goPr
    end if

    imp = im(region)/xref_
    jmp = jm(region)/yref_
    lmp = lmr

    if ( ibeg(region) < iend(region) .and. imp /= iend(region)-ibeg(region)+1 )  &
         call escape_tm('update_parent: stop')
    if ( jmp /= jend(region)-jbeg(region)+1 ) &
         call escape_tm('update_parent: stop')

    ! ~~ mass (all levels!)

    lmpa = lm(region)
    allocate(to_parent(imp,jmp,lmpa))

    to_parent = 0.0
    do i=1,imr
       ip = 1 + (i-1)/xref_
       do j=1,jmr
          jp = 1 + (j-1)/yref_
          do l=1,lmpa
             to_parent(ip,jp,l) = to_parent(ip,jp,l) + m(i,j,l)
          end do
       end do
    end do
    if ( ibeg(region) < iend(region) ) then   !no dateline crossing!
       mp(ibeg(region):iend(region),jbeg(region):jend(region),1:lmpa) = &
            to_parent 
    else
       mp(ibeg(region):im(region),jbeg(region):jend(region),1:lmpa)  = &
            to_parent(1:im(region)-ibeg(region)+1,:,:)
       mp(1:iend(region),jbeg(region):jend(region),1:lmpa)  = &
            to_parent(im(region)-ibeg(region)+2:,:,:)
    end if

    deallocate(to_parent)
    
    ! ~~ tracers

    allocate(to_parent(imp,jmp,lmp))

    sum_parent = 0.0
    sum_to_parent = 0.0
    sum_parent_all = 0.0
    sum_to_parent_all = 0.0
    do n=1,nt
       to_parent = 0.0
       do i=1,imr
          ip = 1 + (i-1)/xref_
          do j=1,jmr
             jp = 1 + (j-1)/yref_
             do l=1,lmr
                lp = 1 + (l-1)/zref_
                to_parent(ip,jp,lp) = to_parent(ip,jp,lp) + rm(i,j,l,n)
                if ( n == 1 ) sum_to_parent = sum_to_parent + rm(i,j,l,1)
#ifdef MPI         
                if ( which_par == 'tracer' .and. myid/=pe_first_tracer) &
                     sum_to_parent=0.0 !WP! only on one when in tracer
#endif
             end do
          end do
       end do
       if ( ibeg(region) < iend(region) ) then   !no dateline crossing!
          if ( n == 1 ) then
             sum_parent = sum(rmp(ibeg(region):iend(region), &
                  jbeg(region):jend(region), 1:lmr,1) )
#ifdef MPI     
             if ( which_par == 'tracer' .and. myid /= pe_first_tracer ) &
                  sum_parent=0.0 !WP! only on one when in tracer
#endif
          end if
          ! CMK this line gave a problem 1:lmr /= to_parent dimension!
          rmp(ibeg(region):iend(region),jbeg(region):jend(region),1:lmr,n) = &
               to_parent 
       else
          if (n == 1) then
             sum_parent = sum(rmp(ibeg(region):im(region), &
                  jbeg(region):jend(region),   &
                  1:lmr,1) ) + &
                  sum(rmp(1:iend(region), &
                  jbeg(region):jend(region),   &
                  1:lmr,1) ) 
#ifdef MPI     
             if ( which_par == 'tracer' .and. myid /= pe_first_tracer ) &
                  sum_parent=0.0 !WP! only on one when in tracer
#endif
          end if
          rmp(ibeg(region):im(region),jbeg(region):jend(region),1:lmr,n)  = &
               to_parent(1:im(region)-ibeg(region)+1,:,:)
          rmp(1:iend(region),jbeg(region):jend(region),1:lmr,n)  = &
               to_parent(im(region)-ibeg(region)+2:,:,:)
       end if
    end do   ! ntracet

#ifdef with_budgets

#ifdef MPI
    call mpi_allreduce(sum_parent,sum_parent_all,1,my_real, &
         mpi_sum,communicator,ierr)
    call mpi_allreduce(sum_to_parent,sum_to_parent_all,1,my_real, &
         mpi_sum,communicator,ierr)
#else
    sum_parent_all = sum_parent
    sum_to_parent_all = sum_to_parent
#endif

    sum_update(my_parent) = sum_update(my_parent) &
         + sum_to_parent_all - sum_parent_all

#endif

    deallocate(to_parent)

    nullify(m)
    nullify(rm)
#ifdef slopes
    nullify(rxm)
    nullify(rym)
    nullify(rzm)
#endif

    nullify(mp)
    nullify(rmp)
#ifdef slopes
    nullify(rxmp)
    nullify(rymp)
    nullify(rzmp)
#endif

    ! mk still to do : limit the slopes after rm update!

  end subroutine update_parent



end module zoom_tools