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

  use GO, only : gol, goPr, goErr

  implicit none

  ! ------------ interface ---------------------

  private

  public  ::  Advect_Init, Advect_Done
  public  ::  dynam0

  ! --- const --------------------------------------

  character(len=*), parameter  ::  mname = 'advect'


contains


  ! ================================================================
  
  
  subroutine Advect_Init( status )
  
    use dims , only : nregions
    use Meteo, only : Set
    use Meteo, only : sp1_dat, sp2_dat, sp_dat
#ifdef with_prism
    use Meteo, only : tsp_dat
#endif
    use Meteo, only : mfu_dat, mfv_dat, mfw_dat
    use Meteo, only :  pu_dat,  pv_dat,  pw_dat

    ! --- in/out --------------------------------
    
    integer, intent(out)           ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Advect_Init'
    
    ! --- local --------------------------------
    
    integer          ::  n
    
    ! --- begin --------------------------------
    
    ! loop over all (zoom) regions:
    do n = 1, nregions
    
      ! surface pressures and mass fluxes for advection
      call Set(  sp_dat(n), status, used=.true. )
      call Set( sp1_dat(n), status, used=.true. )
      call Set( sp2_dat(n), status, used=.true. )
#ifdef with_prism
      call Set( tsp_dat(n), status, used=.true. )
#endif
      call Set( mfu_dat(n), status, used=.true. )
      call Set( mfv_dat(n), status, used=.true. )
      call Set( mfw_dat(n), status, used=.true. )
      call Set(  pu_dat(n), status, used=.true. )
      call Set(  pv_dat(n), status, used=.true. )
      call Set(  pw_dat(n), status, used=.true. )
    
    end do  ! regions
    
    ! ok
    status = 0
    
  end subroutine Advect_Init
  
  
  ! ***
  
  
  subroutine Advect_Done( status )
  
    ! --- in/out --------------------------------
    
    integer, intent(out)           ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/Advect_Done'

    ! --- begin --------------------------------
    
    ! nothing to be done

    ! ok
    status = 0
    
  end subroutine Advect_Done
  
  
  
  ! ================================================================
  

  !-----------------------------------------------------------------------
  !
  !****  dynam0          - calculates vertical massfluxes        v 9.1
  !
  !      programmed by           mh      mpi HH          1-oct-1991
  !
  !      purpose
  !      -------
  !      This sr is CALLed each time after new massfluxes are read in.
  !      Calculate the new air-masses in each grid box from the
  !      surface pressure.
  !      Calculate the amount of air moved in each advection substep,
  !      also in the vertical direction.
  !
  !      interface
  !      ---------
  !      CALL dynam0
  !
  !      method
  !      ------
  !      integrate air conservation equation in the vertical direction
  !      in order to obtain the vertical massfluxes.
  !
  !      externals
  !      ---------
  !      none
  !
  !      reference
  !      ---------       
  !      see manual
  !
  ! modified by mk (1999) for zoom version.
  !-----------------------------------------------------------------------

  subroutine dynam0(region)

    use dims,        only : im, jm, lm, ndyn, tref, zero, bt, xcyc
    use global_data, only : wind_dat
    use meteo      , only : pu_dat, pv_dat
#ifdef MPI 
    use mpi_const,   only : ntracetloc
    use mpi_comm,    only : barrier_t
#endif

    implicit none

    ! input
    integer,intent(in)   :: region 

    ! local
    real,dimension(:,:,:),pointer           :: pu,pv,am,bm,cm
    real,dimension(im(region),jm(region))   :: pit
    real,dimension(im(region),jm(region),lm(region)-1)   :: sd
    real,dimension(im(region),jm(region),lm(region))     :: conv_adv

    real    :: dtu, dtv, dtw
    integer :: i, j, l, lmr

    ! start

#ifdef MPI
    if ( ntracetloc == 0 ) return !WP! only processors with nonzero ntracet
#endif

    ! leave if fluxes are not in use:
    if ( .not. pu_dat(region)%used ) return
    if ( .not. pv_dat(region)%used ) return
    ! set pointers:
    pu => pu_dat(region)%data
    pv => pv_dat(region)%data
    am => wind_dat(region)%am_t
    bm => wind_dat(region)%bm_t
    cm => wind_dat(region)%cm_t

    !---- length of advection substeps (in seconds) in the three directions
    dtu=ndyn/(2.*tref(region))  ! again removed the revert! cmk
    dtv=ndyn/(2.*tref(region))
    dtw=ndyn/(2.*tref(region))

    ! number of levels in this region:
    lmr = lm(region)
    
    ! init to zero:
    am = zero
    bm = zero
    cm = zero

    !
    !     compute conv_adv, the horizontal mass convergence
    !      the arrays pu and pv contain the horizontal air mass fluxes crossing 
    !      the grid box boundaries. pu(i,j,l) is the eastward mass flux in kg/sec
    !      at the eastern edge of box i,j,l; pv(i,j,l) is the northward
    !      mass flux at the southern edge of box i,j,l
    !      
    do l=1,lmr
       do  j=1,jm(region)
          do  i=1,im(region)
             conv_adv(i,j,l)=pu(i-1,j,l)-pu(i,j,l)+pv(i,j,l)-pv(i,j+1,l)
          end do
       end do
    end do
    !
    !     compute pit, the vertiCALLy integrated conv_adv
    !
    do j=1,jm(region)
       do  i=1,im(region)
          pit(i,j)=conv_adv(i,j,1)
          do l=2,lmr
             pit(i,j)=pit(i,j)+conv_adv(i,j,l)
          end do
          !
          !     compute the vertical massflux on the box boundaries (in kg/s)
          !mk  in the hybrid system, the tendency in ps is given by bt. Bt is
          !mk  already normalized. Now distribute pit such that the new mass 
          !mk  distribution after advection exactly matches the 
          !mk  exactly coincides with the distribution that is obtained 
          !mk  from the new surface pressure.
          !
          sd(i,j,lmr-1)=conv_adv(i,j,lmr) - (bt(lmr)-bt(lmr+1))*pit(i,j)
          do  l=lmr-2,1,-1
             sd(i,j,l)=sd(i,j,l+1)+conv_adv(i,j,l+1)-(bt(l+1)-bt(l+2))*pit(i,j)
          end do
       end do
    end do
    !
    !     compute amount of air moved each advection substep, am, bm or cm (kg)
    !
    !      am(i,j,l) is the amount of air moved eastward at the eastern
    !      boundary of grid box i,j,l
    !     compute cm
    !      cm(i,j,l) is the amount of air moved in upward direction at the
    !      upper boundary of grid box i,j,l
    !     compute bm
    !      bm(i,j,l) is the amount of air moved northward at the southern
    !      boundary of grid box i,j,l
    !
    do l=1,lmr
       do j=1,jm(region)
          do i=0,im(region)
             am(i,j,l)=dtu*pu(i,j,l)
          end do
       end do
    end do

    do l=1,lmr
       do j=1,jm(region)+1
          do i=1,im(region)
             bm(i,j,l)=dtv*pv(i,j,l)
          end do
       end do
    end do

    do l=1,lmr-1
       do j=1,jm(region)
          do i=1,im(region)
             cm(i,j,l)=-dtw*sd(i,j,l)
          end do
       end do
    end do 

    if ( xcyc(region) == 1 ) then
       am(0,:,:) = am(im(region),:,:)
       am(im(region)+1,:,:) = am(1,:,:)
    end if

    nullify(pu)
    nullify(pv)
    nullify(am)
    nullify(bm)
    nullify(cm)

#ifdef MPI 
    call barrier_t
#endif

  end subroutine dynam0


end module advect
