ecearth3/sources/tm5mp/base/advectx__slopes.F90

!### macro's #####################################################
!
#define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr
#define IF_NOTOK_RETURN(action) if (test/=0) then; TRACEBACK; action; return; end if
#define IF_ERROR_RETURN(action) if (test> 0) then; TRACEBACK; action; return; end if
!
#include "tm5.inc"
!
!-----------------------------------------------------------------------------
!                    TM5                                                     !
!-----------------------------------------------------------------------------
!BOP
!
! !MODULE:      ADVECTX
!
! !DESCRIPTION: 
!\\
!\\
! !INTERFACE: 
!
MODULE ADVECTX
  !
  ! !USES:
  !
  USE GO,           ONLY : gol, goPr, goErr
  USE TM5_DISTGRID, ONLY : DGRID, GET_DISTGRID, REDUCE, UPDATE_HALO_JBAND
  USE PARTOOLS,     ONLY : isRoot,  npe_lon

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  PUBLIC :: advectxzoom
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=*), parameter ::  mname='advectx'
  !
  ! !REVISION HISTORY: 
  !    1 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
  !    9 Jan 2013 - P. Le Sager - works with reduced grid
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS
    
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    ADVECTXZOOM
  !
  ! !DESCRIPTION: set parameters for advectx 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE ADVECTXZOOM( REGION, TEST )
    ! !USES:
    !
    !
    use dims,          only : lm
#ifdef with_budgets    
    use budget_global, only : sum_advection
#endif
    use global_data,   only : mass_dat, wind_dat
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in) :: region
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: test
    !
    ! !REVISION HISTORY: 
    ! written by patrick berkvens and mike botchev, march-june 1999
    ! updated and modified by MK, dec 2002
    !    1 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer           :: i1, j1, i2, j2, istat, lmr
    real              :: sum_old,     sum_new

    !-------------- start

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

#ifdef with_budgets
    ! total mass tracer #1
    call REDUCE( dgrid(region), mass_dat(region)%rm(:,:,:,1), mass_dat(region)%halo, 'SUM', sum_old, istat)
#endif
    
    CALL ADVECTX_WORK(region, j1, j2, 1, lmr, test)
    IF_NOTOK_RETURN(test=1)

#ifdef with_budgets

    call REDUCE( dgrid(region), mass_dat(region)%rm(:,:,:,1), mass_dat(region)%halo, 'SUM', sum_new, istat)

    if( isRoot ) then
       sum_advection(region) = sum_advection(region) + sum_new - sum_old
    end if
#endif

    test=0

  END SUBROUTINE ADVECTXZOOM
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    ADVECTX_WORK
  !
  ! !DESCRIPTION: makes reduced grid pre-/postprocessing, call dynamu
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE ADVECTX_WORK( region, js, je, ls, le, test)
    !
    ! !USES:
    !
    use dims,        only : im, jm, lm
    use redgridZoom, only : grid_reduced, nred, uni2red, uni2red_mf, red2uni, idle_xadv
    use global_data, only : wind_dat, mass_dat
    use meteodata  , only : m_dat
    use chem_param,  only : ntracet
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in) :: region
    integer,intent(in) :: js
    integer,intent(in) :: je
    integer,intent(in) :: ls
    integer,intent(in) :: le
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: test
    !
    ! !REVISION HISTORY: 
    ! 
    !    written by mike botchev, march-june 1999
    !    1 Feb 2012 - P. Le Sager - adapted for lon-lat MPI domain decomposition
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    real, dimension(:,:,:), pointer     :: m
    real, dimension(:,:,:), pointer     :: am
    real, dimension(:,:,:), allocatable :: m_uni   ! for reduced grid...
    real, dimension(:,:,:), allocatable :: am_uni  ! for reduced grid...
    real                 :: mxval
    integer,dimension(3) :: mxloc
    integer              :: imr, jmr, lmr, i1, i2, j1, j2, halo

    ! --------- START

#ifndef slopes
    test=1
    TRACEBACK
    return
#endif

    ! Transform to the reduced grid
    if ( grid_reduced ) then 

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

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

       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(-1:0,j1:j2,:) = am(imr-1:imr,j1:j2,:)
          am(i2+1,j1:j2,:) = am(        1,j1:j2,:)
       else
          call UPDATE_HALO_JBAND(dgrid(region), am, halo, test )
          IF_NOTOK_RETURN(test=1)
       endif
       
       ! save non-reduced m and am in m_uni and am_uni:
       m_uni  = m
       am_uni = am

       ! reduce m
       call uni2red(region,i1,i2,j1,j2,test)
       IF_NOTOK_RETURN(test=1)

       ! reduce am (and update halo of reduced zonal bands)
       call uni2red_mf( region,i1,i2,j1,j2,test )
       IF_NOTOK_RETURN(test=1)
       
    end if
        
    ! transport
    if (.not.idle_xadv) then
       CALL DYNAMU (region,js,je,ls,le, test)
       IF_NOTOK_RETURN(test=1)
    endif
    
    ! Transform from reduced to uniform 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)

       ! redistribute rm,rxm,rym,rzm by using m_uni and m
       halo = m_dat(region)%halo
       call red2uni(region, i1, i2, j1, j2, halo, m_uni, test)
       
       ! 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

    test=0

  END SUBROUTINE ADVECTX_WORK
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DYNAMU
  !
  ! !DESCRIPTION: east-west tracer transport: calculate amount of tracer moved in an east-west
  !                advection substep
  !
  !   method    : slopes scheme
  !   reference : Russel and Lerner, 1979
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DYNAMU( region, js, je, ls, le, test)
    !
    ! !USES:
    !
    use dims,        only : okdebug, nregions, parent
    use dims,        only : im, jm, lm, xref, yref, zref, tref
    use dims,        only : zero, one, xi, nxi, nloop_max, limits, xcyc
    use redgridZoom, only : grid_reduced, nred, imredj
    use redgridZoom, only : mfl_alt1, m_alt1, rm_alt1, rxm_alt1, rym_alt1, rzm_alt1
    use global_data, only : wind_dat, mass_dat
    use meteodata  , only : m_dat
#ifdef with_budgets    
    use budget_global, only : budget_flux, apply_budget_global_flux
    use budget_global, only : iflux1, iflux2, jflux1, jflux2
#endif
    use chem_param,  only : ntracet, ra
    use partools,    only : par_reduce, isRowRoot
    !
    ! !INPUT PARAMETERS:
    !
    integer,intent(in) :: region
    integer,intent(in) :: js
    integer,intent(in) :: je
    integer,intent(in) :: ls
    integer,intent(in) :: le
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out) :: test
    !
    ! !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

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

    real, dimension(:,:,:,:), pointer     :: rm, rxm, rym, rzm
    real, dimension(:,:,:),   pointer     :: m, am
    real, dimension(:),       allocatable :: f, pf, fy, fz, mnew, mx
    real, dimension(:,:),     allocatable :: wrk2
    
    integer, allocatable :: ie(:)
    integer :: i, is, j, l, n, offsetn
    integer :: i1, i2, j1, j2
    logical :: WestBorder, EastBorder
    integer :: imr, jmr, lmr, tref_, xref_, yref_, zref_, ic, iloop, iemax
    integer :: nloop(jm(region),lm(region))
    real    :: min_one, max_one, min_all, max_all
    real    :: my_alpha, alpha, max_alpha, alphamax, x, rmold
    
    logical            :: cfl_ok
    integer, parameter :: max_nloop = 50
    integer            :: ns, ne, jss, jee, iie
    real               :: l_xi, g_xi, l_nl, g_nl
    character(len=2)   :: sloop

    
    !------------- Start

    ! ------------
    ! Bounds
    ! ------------
    CALL GET_DISTGRID( dgrid(region),                                     &
                       I_STRT=i1,                I_STOP=i2,               &
                       J_STRT=j1,                J_STOP=j2,               &
                       hasWestBorder=WestBorder, hasEastBorder=EastBorder )
    lmr = lm(region)
    imr = im(region)
    allocate(ie(j1:j2))

    ! ------------
    ! Work arrays
    ! ------------
    if ( ( grid_reduced ) .and. (npe_lon /= 1) .and. (nred(region)/=0) ) then
       am  => mfl_alt1
       m   =>   m_alt1
       rm  =>  rm_alt1
       rxm => rxm_alt1
       rym => rym_alt1
       rzm => rzm_alt1
       is=1
       ie=imr
       iemax=maxval(imredj(j1:j2,region))
    else
       am  => wind_dat(region)%am
       m   =>    m_dat(region)%data
       rm  => mass_dat(region)%rm
       rxm => mass_dat(region)%rxm
       rym => mass_dat(region)%rym
       rzm => mass_dat(region)%rzm
       is=i1
       ie=i2       
       iemax=i2
    endif

    allocate(   f( is-1:iemax+1 ))  ! halo=1
    allocate(  pf( is-1:iemax+1 ))  ! halo=1
    allocate(  fy( is-1:iemax+1 ))  ! halo=1
    allocate(  fz( is-1:iemax+1 ))  ! halo=1
    allocate(  mx( is-2:iemax+2 ))  ! halo=2 (was 1) 
    allocate(mnew( is-1:iemax+1 ))  ! halo=1 (was 0)

    allocate(wrk2( is-2:iemax+2, ntracet ))  ! halo=2

    ! ------------
    ! 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
          iie=ie(j)
            m( -1:0, j,:)   =   m(iie-1:iie,j,:)
           rm( -1:0, j,:,:) =  rm(iie-1:iie,j,:,:)
          rxm( -1:0, j,:,:) = rxm(iie-1:iie,j,:,:)
          rym( -1:0, j,:,:) = rym(iie-1:iie,j,:,:)
          rzm( -1:0, j,:,:) = rzm(iie-1:iie,j,:,:)
          
            m(iie+1:iie+2, j,:)   =   m(1:2,j,:)
           rm(iie+1:iie+2, j,:,:) =  rm(1:2,j,:,:)
          rxm(iie+1:iie+2, j,:,:) = rxm(1:2,j,:,:)
          rym(iie+1:iie+2, j,:,:) = rym(1:2,j,:,:)
          rzm(iie+1:iie+2, j,:,:) = rzm(1:2,j,:,:)
       end do
       ! halo of am is filled in advectx_work (thru uni2red_mf for reduced part)
    else

       CALL UPDATE_HALO_JBAND( dgrid(region),  am, wind_dat(region)%halo, test)
       CALL UPDATE_HALO_JBAND( dgrid(region),   m,    m_dat(region)%halo, test)
       CALL UPDATE_HALO_JBAND( dgrid(region),  rm, mass_dat(region)%halo, test)
       CALL UPDATE_HALO_JBAND( dgrid(region), rxm, mass_dat(region)%halo, test)
       CALL UPDATE_HALO_JBAND( dgrid(region), rym, mass_dat(region)%halo, test)
       CALL UPDATE_HALO_JBAND( dgrid(region), rzm, mass_dat(region)%halo, test)

    end if
    
    ! if requested limit zonal slopes such that no negative tracer masses
    ! should occur
    if ( limits ) then
       do n = 1,  ntracet
       do l = ls, le
       do j = js, je
       do i = is - 1, ie(j) + 1
             rxm(i,j,l,n) = max (min (rxm(i,j,l,n), rm(i,j,l,n)), -rm(i,j,l,n))
       end do
       end do
       end do
       end do
    end if

    ! init stat holders
    l_nl  = 0.0
    l_xi  = 0.0  ! max alpha

    ! ================= Find number of iterations needed
    
    LATBAND: DO l= ls, le   ! work in 1D, so final nloop can differ with latitude and levels
       DO j = js, je
          
       iie = ie(j)

       ! Determine number of loops required to avoid CFLs...   
       nloop(j,l) = 1   ! default run the loop one time!
       alpha      = 2.0
       max_alpha  = 0.0

       do while ( abs(alpha) >= one .and. nloop(j,l) < max_nloop )

          ! copy initial mass to temp array
          mx(is-2:iie+2) = m(is-2:iie+2,j,l)

          xloop: do iloop = 1, nloop(j,l)
             cfl_ok = .true.

             indloop : do i=is-1,iie   
                if (am(i,j,l)>=zero)  then
                   my_alpha=am(i,j,l)/mx(i)
                else
                   my_alpha=am(i,j,l)/mx(i+1)
                end if
                max_alpha = max(max_alpha, abs(my_alpha))
                if((abs(my_alpha)>=one)) exit indloop
             end do indloop

             ! sync
             if (grid_reduced .and. npe_lon/=1.and.nred(region)>0) then
                alpha=max_alpha
             else
                call Par_Reduce(max_alpha, 'MAX', alpha, test, all=.true., row=.true.)
             endif
             
             if (alpha>=one) then
                cfl_ok = .false.
                exit xloop
             end if

             ! cfl ok, then update mass for next iteration if any
             if(iloop/=nloop(j,l)) then
                ! is:iie updated
                ! is-1, iie+1 are BC, and also updated (assumption that
                ! xcyc(region)==1), so we do not need to call update halo
                ! every other iteration.
                mx(is-1:iie+1) = mx(is-1:iie+1) + am(is-2:iie,j,l) - am(is-1:iie+1,j,l)
                
                if (modulo(iloop,2)==0) then
                   if ( ( grid_reduced .and. npe_lon==1 ).or.(grid_reduced .and. npe_lon/=1.and.nred(region)>0)) then
                      mx(       -1:0) = mx(iie-1:iie)
                      mx(iie+1:iie+2) = mx(      1:2)
                   else
                      call UPDATE_HALO_JBAND( dgrid(region), mx, 2, test)
                   end if
                end if
             end if
          end do xloop
          

          ! CLF not ok : reduce mass flux
          if(.not.cfl_ok) then 

             am(is-2:iie+1,j,l) = am(is-2:iie+1,j,l)*nloop(j,l)/(nloop(j,l)+1)

             nloop(j,l) = nloop(j,l) + 1
             max_alpha  = 0.0

             if (nloop(j,l) == max_nloop) then
                write(gol,*)'nloop too high'; call goErr
                test=1; TRACEBACK
                return
             end if
          end if
          
       end do !while alpha>1

       l_xi = MAX (l_xi, alpha)
       
    end DO
    end DO LATBAND

    
    ! ================= UPDATE MASSES
    
    LATBAND2: DO l= ls, le
    DO j = js, je

       iie = ie(j)

       CFL: do iloop = 1,nloop(j,l)

          ! calculate new air mass distribution
          mnew(is-1:iie+1) = m(is-1:iie+1,j,l) + am(is-2:iie,j,l)-am(is-1:iie+1,j,l)

          ! loop over tracers
          TRAC: do n=1,ntracet
             
             ! calculate fluxes for rm,rxm,rym,rzm
             do i=is-1,iie
                if (am(i,j,l)>=zero)  then
                   
                   alpha = am(i,j,l)/m(i,j,l)
                   f(i)  = alpha*(rm(i,j,l,n)+(one-alpha)*rxm(i,j,l,n))
                   pf(i) = am(i,j,l)*(alpha*alpha*rxm(i,j,l,n) - 3.*f(i))
                   fy(i) = alpha*rym(i,j,l,n)
                   fz(i) = alpha*rzm(i,j,l,n)
                   
                else
                   
                   alpha = am(i,j,l)/m(i+1,j,l)
                   f(i)  = alpha*(rm(i+1,j,l,n)-(one+alpha)*rxm(i+1,j,l,n))
                   pf(i) = am(i,j,l)*(alpha*alpha*rxm(i+1,j,l,n) - 3.*f(i))
                   fy(i) = alpha*rym(i+1,j,l,n)
                   fz(i) = alpha*rzm(i+1,j,l,n)
                end if
                !xi(region,1)=max(xi(region,1),abs(alpha))

             end do


#ifdef with_budgets    
             ! 
             ! add up flux budget!
             ! FIXME: note that this woukd be wrong if grid is reduced, so set it to
             ! a quiet NaN (but possible only with F2003)
             !
             if ( apply_budget_global_flux ) then
                
                   if ( (iflux1(region)-1 >= is) .and. (iflux1(region)-1 <= iie) ) then 
                      budget_flux(region)%flux_x1(j,l,n) =  budget_flux(region)%flux_x1(j,l,n) &
                           + f(iflux1(region)-1)*1e3/ra(n)   ! moles
                   endif
                   
                   if ( (iflux2(region)-1 >= is) .and. (iflux2(region)-1 <= iie) ) then 
                      budget_flux(region)%flux_x2(j,l,n) =  budget_flux(region)%flux_x2(j,l,n) &
                           + f(iflux2(region)-1)*1e3/ra(n)
                   endif
             end if
#endif
             !   
             ! calculate new tracer mass, and tracer mass slopes
             !
             do i=is, iie
                rm(i,j,l,n)  = rm(i,j,l,n)  + (f(i-1)-f(i))
                rxm(i,j,l,n) = rxm(i,j,l,n) + (pf(i-1)-pf(i)            &
                             - (am(i-1,j,l)-am(i,j,l))*rxm(i,j,l,n)     &
                             + 3.*((am(i-1,j,l)+am(i,j,l))*rm(i,j,l,n)  &
                             - (f(i-1)+f(i))*m(i,j,l)))/mnew(i)
                !CMK: apply limits again: might be in nloop!
                if ( limits ) rxm(i,j,l,n) = &
                     max(min(rxm(i,j,l,n),rm(i,j,l,n)),-rm(i,j,l,n))
                rym(i,j,l,n) = rym(i,j,l,n) + (fy(i-1)-fy(i))
                rzm(i,j,l,n) = rzm(i,j,l,n) + (fz(i-1)-fz(i))
             end do


          end do TRAC  ! end of n-loop
          
          ! store new air mass in m array
          m(is-1:iie+1,j,l)=mnew(is-1:iie+1)


          ! update anything that changes and may be used with +/- index shift in next iloop
          if(iloop/=nloop(j,l)) then

             if ( ( grid_reduced .and. npe_lon==1 ).or.(grid_reduced .and. npe_lon/=1.and.nred(region)>0)) then
                
                  m( -1:0, j,:)   =   m(iie-1:iie,j,:)
                 rm( -1:0, j,:,:) =  rm(iie-1:iie,j,:,:)
                rxm( -1:0, j,:,:) = rxm(iie-1:iie,j,:,:)
                rym( -1:0, j,:,:) = rym(iie-1:iie,j,:,:)
                rzm( -1:0, j,:,:) = rzm(iie-1:iie,j,:,:)

                  m(iie+1:iie+2, j,:)   =   m(1:2,j,:)
                 rm(iie+1:iie+2, j,:,:) =  rm(1:2,j,:,:)
                rxm(iie+1:iie+2, j,:,:) = rxm(1:2,j,:,:)
                rym(iie+1:iie+2, j,:,:) = rym(1:2,j,:,:)
                rzm(iie+1:iie+2, j,:,:) = rzm(1:2,j,:,:)
             
             else       
                
                wrk2 = rm(:,j,l,:)
                CALL UPDATE_HALO_JBAND( dgrid(region), wrk2, mass_dat(region)%halo, test)
                rm(:,j,l,:) = wrk2

                wrk2 = rxm(:,j,l,:)
                CALL UPDATE_HALO_JBAND( dgrid(region), wrk2, mass_dat(region)%halo, test)
                rxm(:,j,l,:) = wrk2

                wrk2 = rym(:,j,l,:)
                CALL UPDATE_HALO_JBAND( dgrid(region), wrk2, mass_dat(region)%halo, test)
                rym(:,j,l,:) = wrk2

                wrk2 = rzm(:,j,l,:)
                CALL UPDATE_HALO_JBAND( dgrid(region), wrk2, mass_dat(region)%halo, test)
                rzm(:,j,l,:) = wrk2

             end if
          end if
          
       end do CFL  ! iloop

       ! restore 'old' am
       am(is-2:iie+1,j,l) = am(is-2:iie+1,j,l)*nloop(j,l)

    enddo
    enddo LATBAND2
   
    ! store algo info
    l_nl = REAL(maxval(nloop(js:je,ls:le)))

    if (grid_reduced .and. npe_lon/=1) then
       if (nred(region)==0) then
          call Par_Reduce(l_nl, 'MAX', g_nl, test, all=.false., row=.true.)
          l_nl=g_nl
       end if
       if(isRowRoot)then
          call Par_Reduce(l_nl, 'MAX', g_nl, test, all=.false., col=.true.)
       endif
    else
       call Par_Reduce(l_nl, 'MAX', g_nl, test, all=.false.)
    endif
    
    if (grid_reduced .and. npe_lon/=1) then
       if (nred(region)==0) then
          call Par_Reduce(l_xi, 'MAX', g_xi, test, all=.false., row=.true.)
          l_xi=g_xi
       end if
       if(isRowRoot)then
          call Par_Reduce(l_xi, 'MAX', g_xi, test, all=.false., col=.true.)
       endif
    else
       call Par_Reduce(l_xi, 'MAX', g_xi, test, all=.false.)
    endif
    
    deallocate(f)
    deallocate(pf)
    deallocate(fy)
    deallocate(fz)
    deallocate(mnew)
    deallocate(mx, wrk2)
    deallocate(ie)

    if ( isRoot ) then
       nloop_max(region,1) = max (nloop_max(region,1), nint (g_nl))
       xi(region,1)        = max (xi(region,1), g_xi)
    end if

    nullify(am)
    nullify(m)
    nullify(rm)
    nullify(rxm)
    nullify(rym)
    nullify(rzm)

    test=0
    
  END SUBROUTINE DYNAMU
  !EOC

END MODULE ADVECTX