ecearth3/sources/tm5mp/base/advectm_cfl.F90

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

!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      ADVECTM_CFL
!
! !DESCRIPTION: module to check for CFL violation.
!
!  Exactly the same operations are performed on m as in the 'real' advection
! routines. However, rm, rxm, rym, rzm are not updated. The 'old' m is saved
! using store_masses and restored with 'restore_masses'.
!\\
!\\
! !INTERFACE: 
!
MODULE ADVECTM_CFL
  !
  ! !USES:
  !
  use GO,           only : gol, goErr, goPr
  use dims,         only : nregions, maxref
  use dims,         only : revert
  use dims,         only : okdebug
  use tm5_distgrid, only : dgrid, Get_DistGrid, update_halo, update_halo_jband

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public        ::  Init_CFL, Done_CFL
  public        ::  Check_CFL         
  public        ::  Setup_MassFlow    
  !
  ! !PRIVATE TYPES:
  !
  type oldmass_data
     real, pointer                    :: m(:,:,:)
     real, dimension(:,:,:),  pointer :: am
     real, dimension(:,:,:),  pointer :: bm
     real, dimension(:,:,:),  pointer :: cm
  end type oldmass_data
  !
  ! !PRIVATE DATA MEMBERS:
  !
  integer, parameter      ::  max_global_iteration    = 32
  integer, parameter      ::  n_operators             = 3  ! xyz
  integer, parameter      ::  nsplitsteps             = 6  ! xyzzyx
  character, parameter    ::  splitorder(nsplitsteps) = (/'x','y','z','z','y','x'/) 

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

  integer                 ::  global_iteration
  integer                 ::  ndyn_save
  integer                 ::  regionm_status(nregions)

  integer                 ::  mloop_max(nregions,3) = 0     ! local CFL counters
  real                    ::  xim      (nregions,3) = 0.0   !

  logical                 ::  cfl_ok

  type(oldmass_data)      ::  oldmass_dat(nregions)
  
  ! alternate mass array (used on row_roots, to deal with reduced grid). With and w/o halo.
  real, dimension(:,:,:), allocatable, target :: m_alt1
  real, dimension(:,:,:), allocatable :: m_alt2

  character               ::  splitorderzoom(nregions,maxref*nsplitsteps)
  !
  ! !REVISION HISTORY: 
  !   10 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !   
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    INIT_CFL
  !
  ! !DESCRIPTION: initialize data for CFL
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE INIT_CFL
    !
    ! !USES:
    !
    use dims,        only : lm, im
    use global_data, only : wind_dat
    use meteodata,   only : m_dat
    use partools,    only : isRowRoot, npe_lon
    use redgridZoom, only : nred
    !
    ! !REVISION HISTORY: 
    !   10 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer :: region, imr, lmr, halo
    integer :: i1, i2, j1, j2

    ! --- begin -------------------------------

    DO region = 1, nregions

       CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

       lmr = lm(region)

       halo=2
       ALLOCATE (oldmass_dat(region)%m(i1-halo:i2+halo, j1-halo:j2+halo, lmr))

       halo = wind_dat(region)%halo
       
       ALLOCATE (oldmass_dat(region)%am(i1-halo:i2+halo, j1-halo:j2+halo, 0:lmr+1))
       ALLOCATE (oldmass_dat(region)%bm(i1-halo:i2+halo, j1-halo:j2+halo, 0:lmr+1))
       ALLOCATE (oldmass_dat(region)%cm(i1-halo:i2+halo, j1-halo:j2+halo, 0:lmr+1))

    END DO

    ! array to work on row_root (only if really needed)
    if ( nred(1)/=0 .and. (npe_lon/=1)) then
       if(isRowRoot)then
          imr = im(1)
          lmr = lm(1)
          halo=m_dat(1)%halo
          allocate(m_alt1( 1-halo:imr+halo, j1:j2, lmr))
          allocate(m_alt2( imr, j1:j2, lmr))
       else
          allocate(m_alt2(1,1,1))
       endif
    endif
    
    ! set splitorder
    splitorderzoom = ' '
    splitorderzoom(1,1:nsplitsteps) = splitorder

  END SUBROUTINE INIT_CFL
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DONE_CFL
  !
  ! !DESCRIPTION: deallocate data for CFL
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DONE_CFL
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer :: region

    ! info ...
    write (gol,'(" ")'); call goPr
    write (gol,'("CFL info from check_cfl:")'); call goPr
    write (gol,'("    x: mloop_max, xim",i4,f10.4)') mloop_max(1,1), xim(1,1); call goPr
    write (gol,'("    y: mloop_max, xim",i4,f10.4)') mloop_max(1,2), xim(1,2); call goPr
    write (gol,'("    z: mloop_max, xim",i4,f10.4)') mloop_max(1,3), xim(1,3); call goPr
    
    ! clear    
    do region=1, nregions
       deallocate(oldmass_dat(region)%m )
       deallocate(oldmass_dat(region)%am)
       deallocate(oldmass_dat(region)%bm)
       deallocate(oldmass_dat(region)%cm)
    end do

    if (allocated(m_alt1)) deallocate(m_alt1)
    if (allocated(m_alt2)) deallocate(m_alt2)

  END SUBROUTINE DONE_CFL
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    CHECK_CFL
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE CHECK_CFL( t1, t2, n, status )
    !
    ! !USES:
    !
    use GO,          only : TDate, IncrDate, operator(+)
    use global_data, only : wind_dat    
    use dims,        only : ndyn, nsrce, nconv, nchem
    use dims,        only : nread    ! <-- 'ntimestep' from rcfile, read in 'start'
    use datetime,    only : new_valid_timestep
    !
    ! !INPUT PARAMETERS:
    !
    type(TDate), intent(in)  ::  t1, t2
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer, intent(inout)   ::  n  !  number of time steps: n is increased until no cfl
    ! 
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)     ::  status
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    type(TDate)              ::  tr(2)
    integer                  ::  i, ndyn_old
    integer                  ::  region
    real                     ::  fraction

    ! --- begin -------------------------------
    if ( okdebug ) then
       write (gol,'(" checking CFL : initial timestep of  ",i5," sec with ", i3," timesteps")') ndyn, n; call goPr
    endif
    
    ! increase number of time steps until cfl ok or maximum exceeded
    global_iteration = 0
    do 

       ! increase number of time steps:
       global_iteration = global_iteration + 1

       ! check ...
       if ( global_iteration > max_global_iteration ) then
          write (gol,'("exceeded maximum number of time steps")'); call goErr
          write (gol,'("in ",a)') rname; status=1; return; call goErr
       end if

       ! save current air masses in 'oldmass_dat'
       call store_masses( n )   

       ! loop over number of ndyn timesteps;
       ! apply advection over complete (large) time interval
       do i = 1, n

          ! small time interval
          tr(1) = t1 + IncrDate(sec=(i-1)*ndyn)
          tr(2) = t1 + IncrDate(sec=   i *ndyn)

          ! fill am/bm/cm with balanced mass flows, eventually time interpolated
          call Setup_MassFlow( tr, status )
          IF_ERROR_RETURN(status=1)

          ! first half step of advection:
          call determine_cfl_iter(1, status)     ! xyz ...
          IF_NOTOK_RETURN(status=1)
          
          ! second half step if first is ok:
          if ( cfl_ok ) then
             call determine_cfl_iter(1, status)  ! ... zyx
             IF_NOTOK_RETURN(status=1)
          end if
          
          if ( okdebug ) then
             write (gol,'(" checking CFL : time step ",i2," of ",i2," ; cfl_ok : ",l1)') i, n, cfl_ok; call goPr
          endif

          ! cfl not ok ? then leave loop and decrease time step:
          if (.not.cfl_ok) exit

          ! flag ...
          regionm_status(:) = 0

       end do

       ! restore massa at begin of time interval
       call restore_masses

       if ( .not. cfl_ok) then
          ndyn_old = ndyn
          ! new ndyn should be a denominator of e.g. 3 hour 
          call new_valid_timestep( ndyn, nread, nread ) 
          if ( okdebug ) then
             write (gol,'(" checking CFL : reducing timestep to ",i5," sec")') ndyn; call goPr
          endif
          nsrce = ndyn
          nconv = ndyn
          nchem = ndyn
          n = nint( real(n*ndyn_old)/real(ndyn) ) ! should work!
          fraction = real(ndyn)/real(ndyn_old)
          do region = 1, nregions
             wind_dat(region)%am = wind_dat(region)%am*fraction
             wind_dat(region)%bm = wind_dat(region)%bm*fraction
             wind_dat(region)%cm = wind_dat(region)%cm*fraction
          enddo
       else 
          exit
       end if

    end do   ! global_iteration

    ! ok
    status = 0

  end subroutine Check_CFL
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    SETUP_MASSFLOW
  !
  ! !DESCRIPTION: Fill am/bm/cm :
  ! 
  !  o interpolate balanced mass fluxes pu/pv in time (if necessary)
  !  o multiply with dt; store in am/bm/cm
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine Setup_MassFlow( tr, status )
    !
    ! !USES:
    !
    use GO,           only : TDate
    use meteodata,    only : pu_dat, pv_dat
    use meteo,        only : TimeInterpolation
    use advect,       only : dynam0
    !
    ! !INPUT PARAMETERS:
    !
    type(TDate), intent(in)    ::  tr(2)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)       ::  status
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter ::  rname = mname//'/Setup_MassFlow'
    integer                     ::  n

    ! --- begin ----------------------------------
    !
    ! time interpolation
    !
    do n = 1, nregions
       !
       call TimeInterpolation( pu_dat(n), tr, status )
       IF_ERROR_RETURN(status=1)
       !
       call TimeInterpolation( pv_dat(n), tr, status )
       IF_ERROR_RETURN(status=1)
       !
    end do
    !
    ! fill am/bm/cm
    !
    do n = 1, nregions
       ! pu/pv to am/bm/cm
       call dynam0( n, status ) 
       IF_NOTOK_RETURN(status=1)
    end do

    status = 0    

  end subroutine Setup_MassFlow
  !EOC

  
  ! ****************************************************************
  ! ***  
  ! *** mass data
  ! ***  
  ! ****************************************************************
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    STORE_MASSES
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE STORE_MASSES(n)
    !
    ! !USES:
    !
    use dims,          only : nregions, tref, lm
    use global_data,   only : wind_dat
    use meteodata,     only : m_dat
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: n    ! number of ndym timesteps...
    !
    ! !REVISION HISTORY: 
    !   10 Nov 2011 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real, dimension(:,:,:), pointer :: m, am, bm, cm
    integer   :: region, lmr, halo
    integer   :: i1, i2, j1, j2
    integer   :: ntimes

    ! --- begin --------------------------

    do region = 1, nregions

       CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )

       lmr = lm(region)
       
       oldmass_dat(region)%m  = m_dat(region)%data
       oldmass_dat(region)%am = wind_dat(region)%am
       oldmass_dat(region)%bm = wind_dat(region)%bm
       oldmass_dat(region)%cm = wind_dat(region)%cm
       regionm_status(region) = 0

       if (revert == -1) then
          m => m_dat(region)%data
          am => wind_dat(region)%am
          bm => wind_dat(region)%bm
          cm => wind_dat(region)%cm

          ! backward advection (two times since xyz..zyx)

          ntimes = 2*n*tref(region)   ! fluxes for ndyn/2 s
          m(i1:i2,j1:j2,1:lmr) =  m(   i1:i2  ,   j1:j2  , 1:lmr  ) &
                        - ntimes*am( i1-1:i2-1,   j1:j2  , 1:lmr  ) &   !east
                        + ntimes*am(   i1:i2  ,   j1:j2  , 1:lmr  ) &   !west
                        - ntimes*bm(   i1:i2  ,   j1:j2  , 1:lmr  ) &   !south
                        + ntimes*bm(   i1:i2  , j1+1:j2+1, 1:lmr  ) &   !north
                        - ntimes*cm(   i1:i2  ,   j1:j2  , 0:lmr-1) &   !lower
                        + ntimes*cm(   i1:i2  ,   j1:j2  , 1:lmr  )     !upper
          nullify(am)
          nullify(bm)
          nullify(cm)
          nullify(m)
       endif
    end do
    
  END SUBROUTINE STORE_MASSES
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  RESTORE_MASSES
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE RESTORE_MASSES
    !
    ! !USES:
    !
    use dims,        only : nregions
    use global_data, only : wind_dat
    use meteodata,   only : m_dat
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    integer :: region

    do region = 1, nregions
       m_dat(region)%data  = oldmass_dat(region)%m
       wind_dat(region)%am = oldmass_dat(region)%am
       wind_dat(region)%bm = oldmass_dat(region)%bm
       wind_dat(region)%cm = oldmass_dat(region)%cm
    enddo

  END SUBROUTINE RESTORE_MASSES
  !EOC

  
  ! ****************************************************************
  ! ***  
  ! *** mass advection
  ! ***  
  ! ****************************************************************

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DETERMINE_CFL_ITER
  !
  ! !DESCRIPTION: 
  !
  ! Determine a 'global' iteration that is needed to prevent CFL violations.
  ! Normally the operator splitting sequence is:  
  ! xyz        vvzyx
  !    xyzvvzyx     vzyxxyzv  
  ! If a CFL violation will occur anywhere, the global timestep is reduced and
  ! the whole sequence is called two (or more) times. The advantage is that
  ! 'masses' are largely restored by advection from neighbouring cells
  !
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DETERMINE_CFL_ITER( REGION, STATUS )
    !
    ! !USES:
    !
    use dims,   only: parent, tref, revert, ndyn
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)        :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)       :: status
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer                  :: i, tref_
    integer, dimension(3)    :: ll,uu

    ! set flag to signal cfl violations
    cfl_ok = .true.

    ! determine refinement factor with respect to the parent
    tref_ = tref(region)/tref(parent(region))

    do i=1,tref_

       call do_steps(region, status)
       IF_NOTOK_RETURN(status=1)
       if (.not. cfl_ok) return
       
       !do the remaining steps if necessary...!
       if (mod(regionm_status(region),n_operators) /= 0 ) then
          call do_steps(region, status)
          IF_NOTOK_RETURN(status=1)
          if(.not. cfl_ok) return   
       end if
       
    end do

  END SUBROUTINE DETERMINE_CFL_ITER
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DO_STEPS
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine do_steps(region, status)
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer          :: i123, reg, rgi, j
    character        :: tobedone
    character(len=1) :: next_step, prev_step   

    ! --- begin ---
    status=0
    
    if ( okdebug ) then
       write(gol,*) 'do_steps: region ',region
       call goPr
    end if

    do i123=1,n_operators  

       next_step = splitorderzoom(region,regionm_status(region)+1)

       if ( regionm_status(region) /= 0 ) then
          prev_step = splitorderzoom(region,regionm_status(region))
       else
          prev_step = ' '
       end if

       if (okdebug) then
          ! it's only to make work of the code visible -- step to be done 
          ! will be printed with capital letter (X,Y or Z)
          do reg=1,region
             tobedone = ' '
             if ( reg == region ) tobedone = upper(next_step)
             write(gol,*) 'do_steps: ',reg,': ', &
                  splitorderzoom(reg,1:regionm_status(reg)),tobedone
             call goPr
          enddo
       endif

       tobedone = upper(next_step)

       select case(next_step)
       case('x')
          call advectmx
       case('y')
          call dynamvm
       case('z')
          call dynamwm
       case default
          write(gol,*)'do_steps:  strange value in splitorderzoom: ',  &
               splitorderzoom(region,regionm_status(region)) ; call goErr
          write(gol,*)'do_steps:  (must be x, y or z)'; call goErr
          status=1; TRACEBACK
          return
       end select

       if (.not. cfl_ok) return
       regionm_status(region) = regionm_status(region)+1
       ! these statements are needed when
       ! more processes are involved that change rm
       ! for instance emissions.
       ! for m-advection only just leave when ready
       if ( mod(regionm_status(region),n_operators) == 0 ) then
          exit ! e.g.after zyx or xyz 
       end if
    end do

  END SUBROUTINE DO_STEPS
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !FUNCTION:   UPPER
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  character function upper(xyz)
    !
    ! !INPUT PARAMETERS:
    !
    character(1), intent(in) :: xyz
    !
    ! !REVISION HISTORY: 
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    if (xyz=='x') then
       upper = 'X'
    else if (xyz=='y') then
       upper = 'Y'
    else if (xyz=='z') then
       upper = 'Z'
    else
       upper = '_'
    end if

  end function upper
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:     ADVECTMX
  !
  ! !DESCRIPTION:  does reduced grid pre-/postprocessing and calls dynamu
  !
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE ADVECTMX
    !
    ! !USES:
    !
    use dims,         only : im, jm, lm
    use redgridZoom,  only : grid_reduced, nred, uni2red_mf, idle_xadv
    use global_data,  only : wind_dat
    use meteodata,    only : m_dat
    use partools,     only : npe_lon, myid, par_broadcast
    !
    ! !REVISION HISTORY: 
    !    written by mike botchev, march-june 1999
    !   30 Jan 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !   20 Dec 2012 - P. Le Sager - fixed for reduced grid
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/advectmx'
    real, dimension(:,:,:), pointer     :: m
    real, dimension(:,:,:), pointer     :: am
    real, dimension(:,:,:), allocatable :: m_uni   ! m backup
    real, dimension(:,:,:), allocatable :: am_uni  ! am backup
    integer :: region, imr,jmr,lmr, i1, i2, j1, j2, halo, status

    ! Block bounds
    region=1
    CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2 )
    lmr = lm(region)
    
    ! Transform to the reduced grid
    if ( grid_reduced ) then    

       halo = m_dat(region)%halo
       allocate( m_uni( i1-halo:i2+halo, j1-halo:j2+halo, lmr    ))

       halo = wind_dat(region)%halo
       allocate(am_uni( i1-halo:i2+halo, j1-halo:j2+halo, 0:lmr+1))
       
       am => wind_dat(region)%am
       m => m_dat(region)%data

       ! take care of BC for uniform am, so that am_uni is correct when
       !  updating m_uni after the call to dynamun
       if (npe_lon==1) then
          am(0,j1:j2,:) = am(im(region),j1:j2,:)
       else
          call UPDATE_HALO_JBAND(dgrid(region), am, halo, status )
          IF_NOTOK_RETURN(status=1)
       endif
       
       ! save non-reduced m and am in m_uni and am_uni:
       m_uni  = m
       am_uni = am
       
       ! reduce m (simplified local routine)
       call uni2red(region,i1,i2,j1,j2,status)
       IF_NOTOK_RETURN(status=1)
       
       ! reduce am (and update halo of reduced zonal bands)
       call uni2red_mf(region,i1,i2,j1,j2, status)
       IF_NOTOK_RETURN(status=1)
       
    end if
    
    if (.not.idle_xadv) then
       CALL DYNAMUM (region,j1,j2,1,lmr)
    endif
    
    ! idle cores need to know about cfl_ok
    if (nred(region)/=0 .and. (npe_lon/=1)) then ! if on a row of many cores that has reduced grid
       call par_broadcast(cfl_ok, status, row=.true.)
    endif
    
    ! Transform from the reduced grid
    if ( grid_reduced ) then
       
       ! advection on uniform grid
       m_uni(i1:i2,j1:j2,1:lmr) = m_uni(i1:i2,j1:j2,1:lmr) + &
            am_uni(i1-1:i2-1,j1:j2,1:lmr) - am_uni(i1:i2,j1:j2,1:lmr)

       ! refill m with m_uni, but only boxes that have been reduced
       halo = m_dat(region)%halo
       call red2uni(region, i1, i2, j1, j2, halo, m_uni, status)
       IF_NOTOK_RETURN(status=1)

       ! recreate rm/rxm/rym/rzm by using equal masses (not m_uni) 
       !call red2uni_em(region)
              
       ! recreate am:
       am = am_uni
       
       nullify(am)
       nullify(m)
       deallocate(am_uni)
       deallocate(m_uni)
    endif
    
  END SUBROUTINE ADVECTMX
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    UNI2RED
  !
  ! !DESCRIPTION: transforms data (air mass only) from uniform grid to reduced grid
  !
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine uni2red(region,i1,i2,j1,j2,status)
    !
    ! !USES:
    !
    use dims,         only : im,jm,lm
    use redgridZoom,  only : nred, clustsize, jred, imredj
    use meteodata,    only : m_dat
    use partools,     only : isRowRoot, npe_lon
    use tm5_distgrid, only : dgrid, gather_j_band
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)  :: region,i1,i2,j1,j2
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)  :: status
    !
    ! !REVISION HISTORY: 
    ! written by mike botchev, march-june 1999 ; modified by Maarten Krol, dec 2002
    ! 
    !   19 Dec 2012 - P. Le Sager - modified for TM5-MP (decomposition along latitudes only)
    !
    ! !REMARKS:
    !  - simplified version of redgridZoom/uni2red : applied to air mass only, not to the tracers.
    !     
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/uni2red'
    
    real, dimension(:,:,:), pointer :: m, m_alt
    integer :: halo,i,ie,is,j,l,lrg,redfact,lmr,imr
    real    :: summ

    ! start
    m => m_dat(region)%data
    lmr=lm(region)
    imr=im(region)

    
    ! Two scenarios
    if (npe_lon /= 1) then
       !
       ! decomposition along longitudes => work on row_root
       !
       if (nred(region)/=0) then  ! some reduce grid on this tile
       
          ! gather m (need to remove i-halo)
          allocate(m_alt( i1:i2, j1:j2, lmr))
          m_alt = m(i1:i2,j1:j2,:)
       
          CALL GATHER_J_BAND(dgrid(region), m_alt, m_alt2, status, jref=j1, rowcom=.true.)
          IF_NOTOK_RETURN(status=1)

          deallocate(m_alt)
          
          if (isRowRoot) then

             m_alt1( 1:imr,:,:) = m_alt2
             m_alt1(     0,:,:) = m_alt1(  imr,:,:)
             m_alt1(    -1,:,:) = m_alt1(imr-1,:,:)

             do lrg=1,nred(region)

                j = jred(lrg,region)
                redfact=clustsize(j,region)

                do l=1,lmr 
                do i = 1,imredj(j,region)
                      ! the is:ie  array section will be reduced to i
                      is = (i-1)*redfact + 1  
                      ie = i*redfact
                      summ = sum(m_alt1(is:ie,j,l))
                      m_alt1(i,j,l) = summ
                end do  !i
                end do  !l
             enddo
          endif
          
       endif
       
    else
       !
       ! Reduced grid without longitudinal decomposition
       !
       do lrg=1,nred(region)

          j = jred(lrg,region)
          redfact=clustsize(j,region)

          do l=1,lmr 
             do i = 1,imredj(j,region)
                ! the is:ie  array section will be reduced to i
                is = (i-1)*redfact + 1  
                ie = i*redfact
                summ = sum(m(is:ie,j,l))
                m(i,j,l) = summ
             end do  !i
          end do  !l

       enddo
       
    endif
    
    nullify(m)
    status=0
    
  END SUBROUTINE UNI2RED
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RED2UNI
  !
  ! !DESCRIPTION: transforms data from reduced grid back to uniform grid
  !
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine red2uni(region, i1, i2, j1, j2, halo, m_uni, status)
    !
    ! !USES:
    !
    use dims,         only : im, lm
    use redgridZoom,  only : nred, clustsize, jred, imredj
    use meteodata,    only : m_dat
    use partools,     only : isRowRoot, npe_lon
    use tm5_distgrid, only : dgrid, scatter_j_band
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: region, i1, i2, j1, j2, halo
    real,    intent(in) :: m_uni(i1-halo:,j1-halo:,1:)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: status
    !
    ! !REVISION HISTORY: 
    ! written by mike botchev, march-june 1999
    !   20 Dec 2012 - P. Le Sager - TM5-MP version
    !
    ! !REMARKS:
    !  - simplified version of redgridZoom/red2uni : applied to air mass only, not to the tracers
    !     
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'_MOD_red2uni'
    
    real,dimension(:,:,:),  pointer    :: m, m_
    integer :: i, j, l, ie, is, lrg, n, redfact, lmr, imr

    ! start
    m => m_dat(region)%data
    lmr=lm(region)
    imr=im(region)

    ! Two scenarios
    if ((npe_lon /= 1).and.(nred(region)>0)) then
       !
       ! decomposition along longitudes => work on row_root
       !
       
       ! (1) use m_alt1 to fill m in bands that are not reduced - only if necessary
       if (nred(region)< (j2-j1+1)) then
          
          allocate(m_( i1:i2, j1:j2, lmr))
          if (isRowRoot) m_alt2 = m_alt1( 1:imr,:,:)
          
          CALL SCATTER_J_BAND(dgrid(region), m_, m_alt2, status, jref=j1, rowcom=.true.)
          IF_NOTOK_RETURN(status=1)

          m(i1:i2,j1:j2,:) = m_

          deallocate(m_)
       endif
          
       ! (2) bands with reduced grid, just use m_uni
       if (nred(region)/=0) then
          do lrg=1,nred(region)
             
             j = jred(lrg,region)
             m(:,j,:)= m_uni(:,j,:)
                
          end do
       endif
       
       ! update halo
       call UPDATE_HALO_JBAND(dgrid(region), m, m_dat(region)%halo, status )

    else
       !
       ! Reduced grid without longitudinal decomposition
       !
       do lrg=1,nred(region)
    
          j = jred(lrg,region)
          redfact=clustsize(j,region)
          
          m(i1:i2,j,:)= m_uni(i1:i2,j,:)
          m(0,j,:) = m(imr,j,:)
       
       end do
    endif
 
    nullify(m)
    status=0

  end subroutine red2uni
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    RED2UNI_EM
  !
  ! !DESCRIPTION: transforms data from reduced grid back to uniform grid
  !
  !    ****************************************************************
  !    ***** COMMENTED - NOT PORTED TO TM5MP WITH NPE_LON /=1 !! ******
  !    ****************************************************************
  !
  !
  !\\
  !\\
  ! !INTERFACE:
  !
!PLS  subroutine red2uni_em(region, i1, i2, j1, j2)
!PLS    !
!PLS    ! !USES:
!PLS    !
!PLS    use dims, only        : im,jm,lm
!PLS    use redgridZoom, only : nred, clustsize, jred, imredj
!PLS    use meteodata,   only : m_dat
!PLS    !
!PLS    ! !INPUT PARAMETERS:
!PLS    !
!PLS    integer, intent(in) :: region, i1, i2, j1, j2
!PLS    !
!PLS    ! !REVISION HISTORY: 
!PLS    ! written by mike botchev, march-june 1999
!PLS    !   20 Dec 2012 - P. Le Sager - TM5-MP version
!PLS    ! 
!PLS    ! !REMARKS:
!PLS    !  - simplified version of redgridZoom/red2uni_em : applied to air mass only,
!PLS    !     no the tracers
!PLS    !
!PLS    !EOP
!PLS    !------------------------------------------------------------------------
!PLS    !BOC
!PLS 
!PLS    ! local
!PLS    real,dimension(:,:,:),  pointer    :: m
!PLS    integer i,ie,is,j,l,lrg,redfact,lmr
!PLS    real mass
!PLS    ! start
!PLS    lmr=lm(region)
!PLS    m => m_dat(region)%data
!PLS    do lrg=1,nred(region)
!PLS       j = jred(lrg,region)
!PLS       redfact=clustsize(j,region)
!PLS       do l=1,lmr 
!PLS          do i =  imredj(j,region),1,-1
!PLS             is = (i-1)*redfact + 1  
!PLS             ie = i*redfact
!PLS             mass=m(i,j,l); m(is:ie,j,l)= mass/(ie-is+1)
!PLS          enddo
!PLS          m(0,j,l) = m(im(region),j,l)
!PLS       end do
!PLS    end do
!PLS    nullify(m)
!PLS  end subroutine red2uni_em
  !EOC


  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DYNAMUM
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine dynamum(region,js,je,ls,le)
    !
    ! !USES:
    !
    use dims,        only : parent, im, lm, zero, one
    use redgridZoom, only : grid_reduced, nred, imredj, mfl_alt1
    use global_data, only : wind_dat
    use meteodata,   only : m_dat
    use partools,    only : Par_Reduce, npe_lon, par_broadcast, isRowRoot
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in) :: region
    integer,intent(in) :: js
    integer,intent(in) :: je
    integer,intent(in) :: ls
    integer,intent(in) :: le
    !
    ! !REMARKS:
    !   - no J-loop anymore to work on lat-lon domain decomposition [we could keep
    !     the L-loop if faster]
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real, dimension(:,:,:),  pointer    :: m, am
    real, dimension(:,:,:), allocatable :: mnew, mx
    integer :: i, is, j, l, n, i1, i2, j1, j2, status
    integer :: lmr, nloop, iloop, iemax
    integer, allocatable :: ie(:)
    real :: max_alpha
    real :: my_alpha, alpha, alpha2, x, rmold
    logical            :: cfl_okl
    integer, parameter :: max_nloop = 10

    ! ------------
    ! Bounds
    ! ------------
    CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2)
    lmr = lm(region)
    allocate(ie(j1:j2))

    ! ------------
    ! Work arrays
    ! ------------    
    if ( ( grid_reduced ) .and. (npe_lon /= 1) .and. (nred(region)/=0) ) then
       m  => m_alt1
       am => mfl_alt1
       is=1
       ie=im(region)
       iemax=maxval(imredj(j1:j2,region))
    else
       m  => m_dat(region)%data  ! halo=2, 1:lmr
       am => wind_dat(region)%am ! halo=1; 0:lmr+1
       is=i1
       ie=i2       
       iemax=i2
    endif
    
    allocate(mnew( is:iemax,     j1:j2,   ls:le))  ! halo=0; ls:le = 1:lmr    
    allocate( mx(is-1:iemax+1, j1-1:j2+1, ls:le))  ! halo=1, temp arr to store initial mass
        
    ! ------------
    ! Reduced grid & halo
    ! --------------
    if ( ( grid_reduced .and. npe_lon==1 ).or.(grid_reduced .and. npe_lon/=1.and.nred(region)>0)) then
       ie = imredj(j1:j2,region)
       do j=j1,j2
          m(      0, j,:) = m(ie(j),j,:)
          m(ie(j)+1, j,:) = m(    1,j,:)
       end do
       ! halo of am is filled in advectmx (thru uni2red_mf for reduced part)
    else
       CALL UPDATE_HALO_JBAND( dgrid(region), am, wind_dat(region)%halo, status)
       ! IF_NOTOK_RETURN(status=1)
       CALL UPDATE_HALO_JBAND( dgrid(region), m, m_dat(region)%halo, status)
       ! IF_NOTOK_RETURN(status=1)
    end if

    ! ------------
    ! Determine number of loops required to avoid CFLs
    ! -----------------
    nloop = 1   ! default run the loop one time!
    alpha = 2.0
    max_alpha = 0.0
    
    !CMK added oct2003: nloop limit..
    do while ( abs(alpha) >= one .and. nloop < max_nloop )

       ! copy intial mass to temp array
       do j=js,je
          mx(is-1:ie(j)+1,j,:) = m(is-1:ie(j)+1,j,ls:le)
       end do
       
       xloop: do iloop = 1, nloop
          cfl_okl = .true.

          lloop: do l=ls,le
          do j=js,je
          do i=is-1,ie(j)

             if (am(i,j,l)>=zero)  then
                my_alpha=am(i,j,l)/mx(i,j,l)
             else
                my_alpha=am(i,j,l)/mx(i+1,j,l)
             end if

             max_alpha = max(max_alpha, abs(my_alpha))

             if((abs(my_alpha)>=one)) then
                exit lloop
             end if

          end do
          end do
          end do lloop

          if (grid_reduced .and. npe_lon/=1) then
             if (nred(region)==0) then
                call Par_Reduce(max_alpha, 'MAX', alpha2, status, all=.true., row=.true.)
                max_alpha=alpha2
             end if
             if(isRowRoot)then
                call Par_Reduce(max_alpha, 'MAX', alpha, status, all=.true., col=.true.)
             endif
             if (nred(region)==0) then
                call par_broadcast(alpha, status, row=.true.)
             endif
          else
             call Par_Reduce(max_alpha, 'MAX', alpha, status, all=.true.)
          endif
          
          if(status==1) then
             write(gol,*) "error par_all_reduce";call goPr
          end if
                    
          if (alpha>=one) then
             cfl_okl = .false.
             exit xloop
          end if
                    
          ! update mass for next iteration
          if (iloop/=nloop) then
             ! [is:ie]
             do j=js,je
                mx(is:ie(j),j,:) = mx(is:ie(j),j,:) + am(is-1:ie(j)-1,j,ls:le) - am(is:ie(j),j,ls:le)
             end do

             ! bc
             if ( ( grid_reduced .and. npe_lon==1 ).or.(grid_reduced .and. npe_lon/=1.and.nred(region)>0)) then
                do j=js,je
                   mx(      0, j,:) = mx(ie(j),j,ls:le)
                   mx(ie(j)+1, j,:) = mx(    1,j,ls:le)
                end do
             else                
                call UPDATE_HALO_JBAND( dgrid(region), mx, 1, status)
                ! IF_NOTOK_RETURN(status=1)
             end if
          end if
          
       end do  xloop


       if(.not.cfl_okl) then 

          do j=js,je
             ! reduce mass flux
             am(is-1:ie(j),j,ls:le) = am(is-1:ie(j),j,ls:le)*nloop/(nloop+1)
             !FASTER:  am(is-1:ie,:,:) = am(is-1:ie,:,:)*(nloop/(nloop+1))
          end do

          if ( okdebug ) then
             write(gol,*) 'dynamu: nloop increased to', nloop + 1; call goPr
             write(gol,*) 'dynamu: alpha is', alpha              ; call goPr
          end if
          
          nloop = nloop + 1
          max_alpha = 0.0
          
          if (nloop == max_nloop) then
             ! not good solution found: set flag and return
             ! note that am is restored in the calling routines
             ! and m will be restored also!
             cfl_ok = .false.
             deallocate(mx, ie)
             deallocate(mnew)
             nullify(am)
             nullify(m)
             return
          endif

       end if
    end do !while alpha>1

    mloop_max(region,1) = max(mloop_max(region,1),nloop)
    xim(region,1) = max(xim(region,1),alpha)

    ! CFL loop
    !    do iloop = 1,nloop
    !       
    !       ! calculate new air mass distribution
    !       mnew(is:ie,js:je,ls:le) = m(is:ie,js:je,ls:le) + am(is-1:ie-1,js:je,ls:le)-am(is:ie,js:je,ls:le)
    !       m(is:ie,js:je,ls:le) = mnew(is:ie,js:je,ls:le)
    !       
    !    end do  ! iloop : SHOULD BE THE SAME AS
    do j=js,je
       ! calculate new air mass distribution
       m(is:ie(j),j,ls:le) = m(is:ie(j),j,ls:le) + nloop*( am(is-1:ie(j)-1,j,ls:le)-am(is:ie(j),j,ls:le) )
    
       ! restore 'old' am
       am(is-1:ie(j),j,ls:le) = am(is-1:ie(j),j,ls:le)*nloop
    end do

    !Done
    deallocate(mnew, mx, ie)
    nullify(am)
    nullify(m)
    
  end subroutine dynamum
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DYNAMVM
  !
  ! !DESCRIPTION: south-north tracer transport: calculate amount of tracer moved in a
  !                south-north advection substep
  !
  !       method    : slopes scheme
  !       reference : Russel and Lerner, 1979
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DYNAMVM
    !
    ! !USES:
    !
    use dims,        only : yref, zref, jm, lm
    use dims,        only : parent, nregions 
    use dims,        only : zero, one
    use global_data, only : wind_dat
    use meteodata,   only : m_dat
    use chem_param,  only : ntracet
    use partools,    only : Par_Reduce
    !
    ! !REVISION HISTORY: 
    !       programmed by           mh      mpi HH          23-feb-1994
    !       fixed bug in calculating maximum courant number
    !                               mh, Thu, Feb 24, 1994 12:49:27
    !     included code for limits of slopes to prevent negative tracer 
    !     masses                   mh, 20-jun-1994
    !
    !     zoom version written by mike botchev, march-june 1999
    !   30 Jan 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real, dimension(:,:,:),  pointer     :: m,bm
    real, dimension(:,:),    allocatable :: mx
    real, dimension(:,:,:),  allocatable :: mnew
    integer  :: i,j,je,js,l,n,iee,iss
    integer  :: lmr
    integer  :: yref_
    real     :: sfs,sfzs,sfn,sfzn
    real     :: my_beta, beta, max_beta
    integer  :: iloop, nloop
    logical  :: cfl_okl, SouthPole, NorthPole, SouthBorder, NorthBorder
    integer  :: lrg, redfact, ixe, ixs
    real     :: summ
    integer  :: i1, i2, j1, j2, jss, jee, status
    integer  :: is,ie,ls,le
    
    integer, parameter :: region=1
    integer, parameter :: max_nloop = 1

    !--- start ---
        
    ! compute refinement factors with respect to the parent
    yref_ = yref(region)/yref(parent(region))

    ! Indices
    CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2, J_STRT=j1, J_STOP=j2,&
                       hasSouthPole=SouthPole,     hasNorthPole=NorthPole,       & 
                       hasSouthBorder=SouthBorder, hasNorthBorder=NorthBorder    )
    lmr = lm(region)
    js = j1
    je = j2
    
    is=i1
    ie=i2

    ls=1
    le=lmr
    
    ! Work arrays
    allocate(mnew( i1:i2, j1:j2, ls:le)) ! halo=0; ls:le = 1:lmr
    m  => m_dat(region)%data             ! halo=2, 1:lmr
    bm => wind_dat(region)%bm            ! halo=1; 0:lmr+1
    
    call UPDATE_HALO( dgrid(region),  m,    m_dat(region)%halo, status)
    call UPDATE_HALO( dgrid(region), bm, wind_dat(region)%halo, status)

    ! Reduce work domain if border of zoom region, and take care of poles    
    if (NorthBorder) then ! j2=jmr by design
       je = j2 - yref_ + 1
       if (NorthPole) je=j2-1
    endif
          
    if (SouthBorder) then ! j1=1 by design
       js = j1 + yref_ - 1
       if (SouthPole) js=j1+1
    endif

    ! indices for mass update
    jss=js
    jee=je
    if (SouthPole) jss=js-1  ! =j1=1
    if (NorthPole) jee=je+1  ! =j2=jmr

    
    ! loop over vertical layers
    LEVELS: do l=ls,le

       ! compute all inner fluxes
       ! f(*,j,*) is flux entering j-th cell from below?cmk
       nloop = 1   ! determine nloop per layer!
       beta = 2.0
       max_beta = 0.0
       ! allocate(mx(is:ie,js-1:je+1))
       allocate(mx( i1-1:i2+1, j1-1:j2+1) )
       
       do while(abs(beta) >= one .and. nloop <= max_nloop)

          mx(is:ie,js-1:je+1) = m(is:ie,js-1:je+1,l)

          xloop: do iloop = 1, nloop
             cfl_okl = .true.

             uloop:do j=js, je+1
             do i=is, ie
                   if (bm(i,j,l)>=zero) then
                      my_beta=bm(i,j,l)/mx(i,j-1)
                   else
                      my_beta=bm(i,j,l)/mx(i,j)
                   end if
                   max_beta = max(max_beta, abs(my_beta))
                   if (abs(my_beta)>=one) then
                      exit uloop
                   end if
             end do !i
             end do uloop !j

             call Par_Reduce(max_beta, 'MAX', beta, status, all=.true.)
             if(status==1) then
                write(gol,*) "error par_all_reduce";call goErr
             end if

             xim(region,2) = max(xim(region,2), beta)       

             if (beta>=one) then
                cfl_okl = .false.
                exit xloop
             end if
             
             if (iloop/=nloop) then
                mx(is:ie,jss:jee) = mx(is:ie,jss:jee) + bm(is:ie,jss:jee,l) - bm(is:ie,jss+1:jee+1,l)
                CALL UPDATE_HALO( dgrid(region), mx, 1, status)
             end if
             
          end do xloop
          
          if ( .not. cfl_okl ) then 
             ! allow NO iterations in y-dir
             cfl_ok = .false.
             if ( okdebug) print *,'dynamv: ', i,j,l, beta
             deallocate(mx)
             deallocate(mnew)
             nullify(m)
             nullify(bm)
             return
          end if

       end do !while 

       deallocate(mx)
       mloop_max(region,2) = max(mloop_max(region,2),nloop) 

       ! calculate new air mass distribution, and store it in m array
       do iloop = 1,nloop
          
          mnew(is:ie,jss:jee,l) = m(is:ie,jss:jee,l) + bm(is:ie,jss:jee,l) - bm(is:ie,jss+1:jee+1,l)      
          m(is:ie,jss:jee,l) = mnew(is:ie,jss:jee,l)
          
       end do

       ! restore old bm's
       bm(:,:,l) = bm(:,:,l)*nloop

    end do LEVELS ! end of l-loop over vertical layers....

    deallocate(mnew)
    nullify(m)
    nullify(bm)

  END SUBROUTINE DYNAMVM
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DYNAMWM
  !
  ! !DESCRIPTION: vertical tracer transport: calculate amount of tracer moved in a
  !                vertical advection substep
  !
  !       method    : slopes scheme
  !       reference : Russel and Lerner, 1979  
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DYNAMWM
    !
    ! !USES:
    !
    use dims,        only : lm, zref
    use dims,        only : zbeg, zend, epsz, zero, one
    use dims,        only : parent
    use global_data, only : wind_dat
    use meteodata,   only : m_dat
    use partools,    only : Par_Reduce
    !
    ! !REVISION HISTORY:
    ! 
    ! - programmed by mh, mpi HH, 1994-1995
    ! - zoom version written by mike botchev, march-june 1999
    ! - 1 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    real, dimension(:,:,:), pointer     :: m, cm
    real, dimension(:,:,:), allocatable :: mnew, cmold, mold
    real, dimension(:,:),   allocatable :: gamma
    integer  :: i,j,l,le,lee,ls,lss,n,lmr,i1,i2,j1,j2
    integer  :: zref_, halo
    real     :: summ
    real     :: mxval, gam, beta, max_beta
    integer  :: n_advz, iter, lrg, ixs, ixe, redfact, status
    integer  :: is,ie,js,je

    integer, parameter :: iter_limit=1 ! no iterations!
    integer, parameter :: region=1

    !----------- START

    ! Indices & work arrays
    CALL GET_DISTGRID( dgrid(region), I_STRT=i1, I_STOP=i2,J_STRT=j1, J_STOP=j2)
    lmr=lm(region)

    allocate( mnew( i1-1:i2+1, j1-1:j2+1, lmr) )     ! halo=1
    allocate(gamma( i1-1:i2+1, j1-1:j2+1     ) )     ! halo=1

    halo = wind_dat(region)%halo
    allocate(cmold( i1-halo:i2+halo, j1-halo:j2+halo, 0:lmr+1) ) ! same halo as wind_dat)

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

    m  =>    m_dat(region)%data
    cm => wind_dat(region)%cm
    
    is=i1 ; ie=i2 ; js=j1 ; je=j2
    

    ! compute refinement factors with respect to the parent
    zref_ = zref(region)/zref(parent(region))

    ! 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
    
    !    
    ! --- calculate new air mass distribution
    !  
    lss = ls; lee = le    
    cmold=cm
    mold=m
    n_advz=1

    
    cm = cmold/float(n_advz)
    m  = mold

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

       if (abs(zbeg(region)-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)
          lss = 2
       end if

       if (abs(zend(region)-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)
          lee = 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)

       ! reset gamma
       max_beta = 0.

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

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

          max_beta = max(max_beta, abs(gamma(i,j)))

          if ( abs(gamma(i,j)) >= one ) then
             if ( okdebug ) then
                write(gol,*) 'dynamwm: CFL violation'; call goErr
                write(gol,*) ' gamma, m, cm=', max_beta, m(i,j,l), m(i,j,l+1), cm(i,j,l), &
                     '    in region,i,j,l: ',region,i,j,l; call goErr
             endif
             exit cloop
          end if

       end do
       end do
       end do cloop

       call Par_Reduce(max_beta, 'MAX', beta, status, all=.true.)
       xim(region,3)=max(xim(region,3), beta)

       if ( beta >= one ) exit advz

       ! do not allow iterations in Z
       ! no CFLs upto now.....possibly in next?
       m(is:ie,js:je,ls:le)=mnew(is:ie,js:je,ls:le)

    enddo 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 ( beta >= one ) then
       cfl_ok = .false.
       deallocate(mnew)
       deallocate(gamma)
       deallocate(cmold)
       deallocate(mold)
       nullify(m)
       nullify(cm)
       return
    end if


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

    ! do the iteration (Could get rid of mnew)
    mloop_max(region,3) = max(mloop_max(region,3), n_advz)
    do iter=1,n_advz
       if (abs(zbeg(region)-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)
          lss = 2
       end if
       if (abs(zend(region)-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)
          lee = 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)

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

    enddo !iter

    ! finally :
    ! mix according to reduced grid (to allow for the adjoint implementation!)
    !cmk  do l=1,lm(region)
    !cmk     do lrg=1,nred(region)
    !cmk        redfact=clustsize(lrg,region)
    !cmk        j = jred(lrg,region)
    !cmk        do i =  1,imredj(lrg,region)
    !cmk           ixs = (i-1)*redfact + 1
    !cmk           ixe = i*redfact 
    !cmk           summ = sum(m(ixs:ixe,j,l))
    !cmk           m(ixs:ixe,j,l) = summ/(ixe-ixs+1)
    !cmk        end do ! I
    !cmk     end do ! lrg
    !cmk  end do ! l

    cm=cmold    ! reset cm

    deallocate(mnew)
    deallocate(gamma)
    deallocate(cmold)
    deallocate(mold)

    nullify(m)
    nullify(cm)

  END SUBROUTINE DYNAMWM
  !EOC

END MODULE ADVECTM_CFL