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

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

  ! --- in/out ------------------------------

  private

  public :: advect_m, m2phlb1, m2phlb, coarsen_ambm


contains


  ! ********************************************************************


  subroutine m2phlb1(region)
    !-----------------------------------------
    ! routine to convert m(im,jm,lm) into pressure layers.
    ! First, calculate the surface pressure in p(im,jm), THEN
    ! using the hybrid coordinate system, calculate phlb(im,jm,lmp1)
    ! special routine for region 1
    !   mk december 1998
    !-----------------------------------------
    use binas,       only : grav
    use dims
    use global_data, only : region_dat
    use meteo      , only : sp_dat
    use meteo      , only : phlb_dat, m_dat
#ifdef MPI 
    use mpi_const
#endif

    implicit none

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

    ! local
    real,dimension(:,:,:), pointer :: m, phlb
    real,dimension(:,:,:), pointer :: p
    real,dimension(:),     pointer :: dxyp
    integer :: i,j,l,imax,xref_,yref_,imr,jmr
    real    :: maxdiff,m_edge,maxval_one,maxdiff_all,maxval_all
    real    :: mnew
    real,dimension(-1:im(region)+2,-1:jm(region)+2) :: pold,phelp
    integer                                         :: communicator,ll,root_id,lmr

    ! start

    dxyp => region_dat(region)%dxyp
    p    =>   sp_dat(region)%data
    m    =>    m_dat(region)%data
    phlb => phlb_dat(region)%data
    lmr = lm(region)
#ifdef MPI   
       communicator=com_trac  !WP! assign com_trac as communicator
       root_id=root_t
#endif
    
    imr=im(region)
    jmr=jm(region)

    pold = p(:,:,1)
    phelp = zero
    p = zero

    do l=1,lmr
       do j=1,jmr
          do i=1,imr
             phelp(i,j) = phelp(i,j) + m(i,j,l)*grav/dxyp(j)
          end do
       end do
    end do

    ! surface pressure
    p(:,:,1) = phelp

    ! pressure at half level boundaries:
    do  l = 1, lmr+1
      phlb(:,:,l) = at(l) + bt(l)*p(:,:,1)
    end do

    maxdiff = 0.0

    do l = 1,lmr
       do j=1,jmr
          do  i=1,imr
             mnew=(phlb(i,j,l)-phlb(i,j,l+1))*dxyp(j)/grav
             maxdiff = max(maxdiff,abs(mnew-m(i,j,l))/mnew)
          end do
       end do
    end do

    maxval_one=maxval(abs(p(1:im(region),1:jm(region),1)-  &
         pold(1:im(region),1:jm(region)))/ &
         pold(1:im(region),1:jm(region)) )
#ifdef MPI 
    call mpi_allreduce(maxdiff,maxdiff_all,1,my_real,mpi_max,communicator,ierr)
    call mpi_allreduce(maxval_one,maxval_all,1,my_real,mpi_max,communicator,ierr)
    if ( okdebug .and. myid == root_id ) &
         print*,'m2phlb1: maxdiff (mnew-m)/mnew on all pes:',maxdiff_all
    if ( okdebug .and. myid == root_id ) &
         print*,'m2phlb1: max pressure difference (%) on all pes :',maxval_all
#else
    if ( okdebug ) print*,'m2phlb1: maxdiff (mnew-m)/mnew :',maxdiff
    if ( okdebug ) print*,'m2phlb1: max pressure difference (%) :',maxval_one
#endif

    nullify(p)
    nullify(m)
    nullify(phlb)
    nullify(dxyp)

  end subroutine m2phlb1



  subroutine m2phlb(region)
    !-----------------------------------------
    !routine to convert m(im,jm,lm) into pressure layers.
    !First, calculate the surface pressure in p(im,jm), then
    !using the hybrid coordinate system, calculate phlb(im,jm,lmp1)
    ! mk december 1998
    ! special routine for region > 1
    !-----------------------------------------
    use binas,       only : grav
    use dims
    use io_hdf
    use global_data, only : region_dat
    use meteodata  , only : sp_dat
    use meteodata  , only : phlb_dat, m_dat
    use toolbox,     only : escape_tm
#ifdef MPI
    use mpi_const
#endif

    implicit none

    ! input
    integer,intent(in)   ::region 

    ! local
    real,dimension(:,:,:), pointer    :: m, phlb, mpar
    real,dimension(:,:,:), pointer    :: p
    real,dimension(:),     pointer    :: dxyp, dxy
    real,dimension(-1:im(region)+2,-1:jm(region)+2)  :: pold,phelp
    integer                           :: i, j, l, imax, xref_, yref_
    integer                           :: imr, jmr, ipar, jpar, lmr
    real                              :: maxdiff, m_edge, maxdiff_all
    real                              :: mnew
    real,dimension(:,:,:),allocatable :: mn, mo   
    integer                           :: io, sfstart, sfend
    integer                           :: communicator, ll, root_id

    ! start

    p    =>   sp_dat(region)%data
    m    =>    m_dat(region)%data
    phlb => phlb_dat(region)%data
    dxy  => region_dat(parent(region))%dxyp
    
    mpar =>    m_dat(parent(region))%data
    dxyp => region_dat(region)%dxyp

    lmr = lm(region)

    pold = p(:,:,1)

    yref_ = yref(region)/yref(parent(region))
    xref_ = xref(region)/xref(parent(region))
    imr = im(region)
    jmr = jm(region)
    ! check mass in the interface region....
    if ( okdebug ) then
       maxdiff = zero
       do l=1,lmr
          do j=1,jmr-yref_+1,yref_
             do i=1,imr-xref_+1,xref_
                if ((i > xref_ .and. i < imr-xref_+1) .and.  &
                     (j > yref_ .and. j < jmr-yref_+1) ) cycle  ! not on edge...
                if ((xcyc(region) == 1) .and.  &
                     (j > yref_ .and. j < jmr-yref_+1) ) cycle  ! periodic boundary
                if ((touch_np(region) == 1).and.(j > jmr-yref_)) cycle ! np_touching
                if ((touch_sp(region) == 1).and.(j < yref_+1)  ) cycle ! sp_touching
                ipar = ibeg(region) + (i-1)/xref_
                jpar = jbeg(region) + (j-1)/yref_   !to get mass in the parent cell
                m_edge = sum(m(i:i+xref_-1,j:j+yref_-1,l)) !get mass in zoom region
                if ( abs((m_edge-mpar(ipar,jpar,l))/m_edge) > 1e-8) then
                   print*,'m2phlb: Encountered mass-difference at',i,j,l
                   print*,'m2phlb: Parent mass comes from        ',ipar,jpar,l
                   print*,'m2phlb: Parent mass is                ',mpar(ipar,jpar,l)
                   print*,'m2phlb: Child  mass is                ',m_edge
                   print*,'m2phlb: Made out of                   ', &
                        m(i:i+xref_-1,j:j+yref_-1,l)
                end if
                maxdIFf = max(maxdIFf,100*abs((m_edge-mpar(ipar,jpar,l))/m_edge))
             end do
          end do
       end do
       print*,'m2phlb: maxdiff mass  on edges (%)',maxdIFf, 'region:', region
    end if

    phelp = zero
    p = zero
    do l=1,lmr
       do j=1,jmr
          do i=1,imr
             phelp(i,j) = phelp(i,j) + m(i,j,l)*grav/dxyp(j)
          end do
       end do
    end do

    ! surface pressure
    p(:,:,1) = phelp

    ! correct pressure edges for coarse grid area...
    ! MKADJ: one exception: in adjoint pressure is restored
    ! and is already corrected for this if newsrun=.TRUE.
    if ( (revert /= -1) .or. (.not. newsrun) ) then
      do j = 1, jmr
         do i = 1, imr
            if ((i > xref_ .and. i < imr-xref_+1) .and.  &
                 (j > yref_ .and. j < jmr-yref_+1) ) cycle    !not on edge...
            if ((xcyc(region) == 1) .and.  &
                 (j > yref_ .and. j < jmr-yref_+1) ) cycle    !periodic boundary
            if ((touch_np(region) == 1).and.(j > jmr-yref_)) cycle    !np_touching
            if ((touch_sp(region) == 1).and.(j < yref_+1)  ) cycle    !sp_touching
            jpar = jbeg(region) + (j-1)/yref_   !to get dxyp in the parent cell
            p(i,j,1) = p(i,j,1)*dxyp(j)*xref_*yref_/dxy(jpar)
         end do
      end do
    end if

    ! fill half level boundaries:
    do l = 1, lmr+1
      phlb(:,:,l) = at(l) + bt(l)*p(:,:,1)
    end do
    
    if ( okdebug ) then
       allocate(mn(imr,jmr,lmr))
       allocate(mo(imr,jmr,lmr))
       mn = 0.0    !omitting might produce crashes on some systems (PB)
       mo = 0.0
       maxdiff = 0.0
       do l = 1,lmr
          jloop: do  j=1,jmr
             if ( yref_> 1) then
                ! not valid on edge..
                if ( j < yref_+1 .and. touch_sp(region) /= 1) cycle jloop
                if ( j > jmr-yref_ .and. touch_np(region) /= 1) cycle jloop
             endif
             iloop : do  i=1,imr
                if ( xcyc(region) /= 1 .and. (i<xref_+1)) cycle iloop
                if ( xcyc(region) /= 1 .and. (i>imr-xref_)) cycle iloop
                mnew=(phlb(i,j,l)-phlb(i,j,l+1))*dxyp(j)/grav
                mn(i,j,l) = mnew
                mo(i,j,l) = m(i,j,l)
                maxdiff = max(maxdiff,abs(mnew-m(i,j,l))/mnew)
             end do iloop
          end do jloop
       end do

       !WP! execute an allreduce over proper processors (communicator)

#ifdef MPI
       communicator=com_trac
       root_id=root_t

       call mpi_allreduce(maxdiff,maxdiff_all,1,my_real,mpi_max,communicator,ierr)

       if ( myid == root_id) &
            print*,'m2phlb: maxdIFf (mnew-m)/mnew on all pes:',maxdiff_all
       if ( maxdiff_all > 1.e-4)then
          call escape_tm('problem with m2phlb: maxdiff > 1e-4')
       endif
#else
       print*,'m2phlb: maxdIFf (mnew-m)/mnew:',maxdiff
       if ( maxdiff > 1.e-4)then
          call escape_tm('problem with m2phlb: maxdiff > 1e-4')
       endif
#endif

       deallocate(mn)
       deallocate(mo)

    end if

    nullify(p)
    nullify(m)
    nullify(phlb)
    nullify(dxyp)
    nullify(dxy)
    nullify(mpar)

  end subroutine m2phlb


  ! ********************************************************************


  subroutine advect_m(region,ntimes)
    !
    ! subroutine 
    !
    use binas,       only : grav
    use dims
    use global_data, only : region_dat, wind_dat
    use meteodata  , only : m_dat
    use meteo      , only : sp_dat
#ifdef MPI 
    use mpi_const,   only : ntracetloc,myid,root
#endif

    implicit none

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

    ! local
    real,dimension(:,:,:),pointer      :: m,am,bm,cm
    real,dimension(:,:,:),pointer      :: p
    real,dimension(:)    ,pointer      :: dxyp

    real,dimension(:,:,:),allocatable                :: mnew, msave
    real,dimension(:,:)  ,allocatable                :: pnew
    integer   :: imr
    integer   :: jmr
    integer   :: lmr
    integer   :: n,i,j,l
    real      :: ss,so

    ! start

#ifdef MPI
    !WP! only processors with nonzero ntracet allowed
    if ( ntracetloc == 0 ) return
    if(myid==root) print*, 'advect_m: called for region ', region
    if(myid==root) print*, 'advect_m: ntimes            ', ntimes
#else
    print*, 'advect_m: called for region ', region
    print*, 'advect_m: ntimes            ', ntimes
#endif

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

    allocate(mnew(imr,jmr,lmr))
    allocate(msave(imr,jmr,lmr))
    allocate(pnew(imr,jmr))

    m => m_dat(region)%data
    am => wind_dat(region)%am_t
    bm => wind_dat(region)%bm_t
    cm => wind_dat(region)%cm_t
    p  => sp_dat(region)%data
    dxyp => region_dat(region)%dxyp

    msave = m(1:imr,1:jmr,1:lmr)
    mnew  = m(1:imr,1:jmr,1:lmr)
    do i=1,ntimes
       mnew = mnew   + revert*am(0:imr-1,1:jmr  ,1:lmr  ) &   !east   !ADJ re-inserted the revert CMK
                     - revert*am(1:imr  ,1:jmr  ,1:lmr  ) &   !west
                     + revert*bm(1:imr  ,1:jmr  ,1:lmr  ) &   !south
                     - revert*bm(1:imr  ,2:jmr+1,1:lmr  ) &   !north
                     + revert*cm(1:imr  ,1:jmr  ,0:lmr-1) &   !lower
                     - revert*cm(1:imr  ,1:jmr  ,1:lmr  )     !upper
    end do
    pnew = 0.0
    do l=1,lm(region)
       do j=1,jm(region)
          do i=1,im(region)
             pnew(i,j) = pnew(i,j) + mnew(i,j,l)*grav/dxyp(j)
          end do
       end do
    end do

    p(1:imr,1:jmr,1) = pnew   !return new pressure for testing
    m(1:imr,1:jmr,1:lmr) = mnew
    nullify(m)
    nullify(am)
    nullify(bm)
    nullify(cm)
    nullify(p)
    nullify(dxyp)

    deallocate(pnew)
    deallocate(mnew)
    deallocate(msave)

  end subroutine advect_m



  subroutine coarsen_ambm(region)
   !
   ! coarsens boundary values of am,bm of all children regions
   ! in corresponding with region's am,bm
   ! written by mike botchev, march-june 1999
   ! updated by maarten krol, dec 2002
   !
   use dims
   use global_data, only: wind_dat
#ifdef MPI
   use mpi_comm
   use mpi_const
#endif

   implicit none

   ! input
   integer,intent(in) :: region

   ! local
   real,dimension(:,:,:),pointer             :: am,bm,amc,bmc
   real cell1,cellm
   integer child,i,ic,ichild,j,jc,l,lc
   integer xref_,yref_,zref_

   ! start

   if ( okdebug ) print *,'coarsen_ambm:'

   am => wind_dat(region)%am_t
   bm => wind_dat(region)%bm_t

   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)
      !WP! define these on all PE's as we will scatter them later
      !WP! amck holds am of child over levels 

      amc => wind_dat(child)%am_t

#ifdef MPI    
      if ( ntracetloc /= 0 ) then !WP! only PE's with nonzero ntracetloc 
#endif
         if ( okdebug ) print *,'coarsen_ambm: coarsening am,bm in child ',child

         ! write am: loop through the cells of x-walls of the child
         ! do this only if the region is not 360-degrees wide
         if ( im(child)/xref(child) /= im(1) ) then


            do j=jbeg(child),jend(child)
               do l=lbeg(child),lend(child)

                  jc = (j-jbeg(child))*yref_ + 1    ! jc is i_child
                  lc = (l-lbeg(child))*zref_ + 1    ! lc is l_child

                  cell1 = am(ibeg(child)-1,j,l)/(yref_*zref_)
                  cellm = am(iend(child)  ,j,l)/(yref_*zref_)

                  amc(xref_-1,  jc:jc+yref_-1, lc:lc+zref_-1) = cell1
                  amc(im(child)-xref_+1,  jc:jc+yref_-1, lc:lc+zref_-1) = cellm
               end do
            end do

         end if
#ifdef MPI    
      end if
#endif
!PB      nullify(am)
      nullify(amc)

      bmc => wind_dat(child)%bm_t
#ifdef MPI    
      if ( ntracetloc /= 0 ) then !WP! only PE's with nonzero ntracetloc 
#endif

         ! write bm: loop through the cells of y-walls of the child

         do i=ibeg(child),iend(child)
            do l=lbeg(child),lend(child)

               ic = (i-ibeg(child))*xref_ + 1    ! ic is i_child
               lc = (l-lbeg(child))*zref_ + 1    ! lc is l_child

               cell1 = bm(i,jbeg(child)  ,l)/(xref_*zref_)
               cellm = bm(i,jend(child)+1,l)/(xref_*zref_)
               ! do this only if the region does not touch the SP
               if ( touch_sp(child) /= 1 )   &
                    bmc(ic:ic+xref_-1,yref_,lc:lc+zref_-1)=cell1
               ! idem for NP... 
               if ( touch_np(child) /= 1 )  &
                    bmc(ic:ic+xref_-1,jm(child)-yref_+2,lc:lc+zref_-1)=cellm
            end do
         end do
#ifdef MPI    
      end if
#endif
!PB       nullify(bm)
      nullify(bmc)
   end do

!PBi
  nullify(am)
  nullify(bm)
!PBe

 end subroutine coarsen_ambm


end module advect_tools
