!#################################################################
!
! routine to perform main processing and      
! handles time ordering of model processes    
!
!### 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 modelIntegration

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

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

  private

  public :: Proces_Init, Proces_Done, Proces_Update
  public :: Proces_Region
  
  ! local ..

  character(len=1)        ::  output_after_step
  
  ! --- const ----------------------------------
  
  ! module name
  character(len=*), parameter  ::  mname = 'ModelIntegration'
  

  ! --- var ----------------------------------
  
  ! timer handles:
  integer ::  itim_advectx, itim_advecty, itim_advectz
  integer ::  itim_vertical, itim_chemistry, itim_source_sink



contains

  
  ! =====================================================
  
  
  subroutine Proces_Init( status )
  
    use GO,                  only : GO_Timer_Def
    use GO,                  only : TrcFile, Init, Done, ReadRc
    use global_data,         only : rcfile

    use Binas      , only : p_global
    use dims       , only : nregions, nregions_all, iglbsfc
    use global_data, only : rcfile
    use Meteo      , only : sp_dat, phlb_dat, m_dat
    use Meteo      , only : sp1_dat, sp2_dat, tsp_dat, spm_dat
    use Meteo      , only : mfu_dat, mfv_dat
    use Meteo      , only : mfw_dat
    use Meteo      , only : temper_dat, humid_dat
    use Meteo      , only : entu_dat, entd_dat, detu_dat, detd_dat
    use Meteo      , only : lwc_dat, iwc_dat, cc_dat, cco_dat, ccu_dat
    use Meteo      , only : gph_dat, omega_dat
    use Meteo      , only : oro_dat, lsmask_dat
    use Meteo      , only : Set
    use Meteo      , only : Meteo_Alloc
#ifndef without_advection
    use Advect        , only : Advect_Init
#endif
#ifndef without_convection
    use Convection    , only : Convection_Init
#ifndef without_diffusion
    use Diffusion     , only : Diffusion_Init
#endif
#endif
#ifndef without_chemistry
    use Chemistry     , only : Chemistry_Init
#endif
    use Sources_Sinks , only : Sources_Sinks_Init
#ifndef without_dry_deposition
    use Dry_Deposition, only : Dry_Deposition_Init
#endif
#ifndef without_wet_deposition
    use Wet_Deposition, only : Wet_Deposition_Init
#endif
    
    ! --- in/out --------------------------------
    
    integer, intent(out)    ::  status
    
    ! --- const -------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/Proces_Init'
    
    ! --- local ---------------------------------
    
    type(TrcFile)     ::  rcF
    integer           ::  n
    
    ! --- begin ---------------------------------

    call goLabel(rname)
  
    !
    ! select meteo to be used
    !

    ! loop over all (zoom) regions:    
    do n = 1, nregions

      ! current (surface) pressure and mass:  
      call Set(   sp_dat(n), status, used=.true. )   
      call Set(    m_dat(n), status, used=.true. )
      call Set( phlb_dat(n), status, used=.true. )  

      ! surface pressures, tendency, and average (mid)
      call Set( sp1_dat(n), status, used=.true. )
      call Set( sp2_dat(n), status, used=.true. )
      call Set( tsp_dat(n), status, used=.true. )
      call Set( spm_dat(n), status, used=.true. )  

    end do
    

    !
    ! check fpp macro's
    !

    ! macro 'with_zoom' should be defined if nregions>1
    if ( nregions > 1 ) then
#ifndef with_zoom
      write (gol,'("macro with_zoom not defined while nregons = ",i2)') nregions; call goErr
      TRACEBACK; status=1; return
#endif
    end if

    ! chance to have faster code ...
    if ( nregions == 1 ) then
#ifdef with_zoom
      write (gol,'("WARNING - macro with_zoom defined while nregions=1; slower run")'); call goPr
#endif
    end if


    ! 
    ! init processes
    !
    
    ! NOTE: these inits should select the appropriate meteo 
    
#ifndef without_advection
    ! init advection ?
    call Advect_Init( status )
    IF_NOTOK_RETURN(status=1)
#endif
    
#ifndef without_convection
    ! init convec:
    call Convection_Init( status )
    IF_NOTOK_RETURN(status=1)
#ifndef without_diffusion
    ! init diffusion:
    call Diffusion_Init( status )
    IF_NOTOK_RETURN(status=1)
#endif
#endif

#ifndef without_chemistry
    ! init chemistry:
    call Chemistry_Init( status )
    IF_NOTOK_RETURN(status=1)
#endif
    
    ! init sources and sinks:
    call Sources_Sinks_Init( status )
    IF_NOTOK_RETURN(status=1)

#ifndef without_dry_deposition
    ! setup dry depos:
    call Dry_Deposition_Init( status )
    IF_NOTOK_RETURN(status=1)
#endif
    
#ifndef without_wet_deposition
    ! setup wet depos:
    call Wet_Deposition_Init( status )
    IF_NOTOK_RETURN(status=1)
#endif
    
    
    !
    ! allocate meteo
    !
    
    ! some meteo implies other ...
    do n = 1, nregions_all
      ! convec requires gph and omega:
      if ( entu_dat(n)%used ) then
        call Set(    gph_dat(n), status, used=.true. )
        call Set(  omega_dat(n), status, used=.true. )
      end if
      ! omega (Pa/s) is computed form mfw (kg/s):
      if ( omega_dat(n)%used ) then
        call Set(    mfw_dat(n), status, used=.true. )
      end if
      ! gph is computed from oro/sp/temper/humid
      if ( gph_dat(n)%used ) then
        call Set(    oro_dat(n), status, used=.true. )
        call Set(     sp_dat(n), status, used=.true. )
        call Set( temper_dat(n), status, used=.true. )
        call Set(  humid_dat(n), status, used=.true. )
      end if
      ! sp is interpolated in time between sp1 and sp2:
      if ( sp_dat(n)%used ) then
        call Set(    sp1_dat(n), status, used=.true. )
        call Set(    sp2_dat(n), status, used=.true. )
      end if
      ! cloud covers and overhead/underfeet cloud covers
      if ( cc_dat(n)%used ) then
        call Set( cco_dat(n), status, used=.true. )
        call Set( ccu_dat(n), status, used=.true. )
      end if
    end do

    !>>> Please use obscure meteo fields only if you really need them ...
    !    If too many fields are used by default, even a simple meteo production run
    !    of an extra surface field will request fields that are not used anyway ...
    !
    write (gol,'("WARNING - ")'); call goPr
    write (gol,'("WARNING - global orography not in use by default anymore;")'); call goPr
    write (gol,'("WARNING - see the comment in ",a)') rname; call goPr
    write (gol,'("WARNING - ")'); call goPr
    !
    ! always store also the orography on 1x1:
    ! call Set(    oro_dat(iglbsfc), status, used=.true. )
    !
    !<<<
    
    ! allocate used meteo
    call Meteo_Alloc( status )
    IF_NOTOK_RETURN(status=1)
    

    !
    ! default values
    !

    ! Fill default values for surface pressure;
    ! might be required for mass fields used in init routines
    do n = 1, nregions_all
      if (  sp_dat(n)%used )  sp_dat(n)%data = p_global
      if ( spm_dat(n)%used ) spm_dat(n)%data = p_global
    end do
    
    !
    ! output sampling
    !
    
    ! read how to sample the output:
    ! x,y,z,c,v,e   : always after this process
    ! a             : after all steps (testing, not recommended)
    ! d             : the 'old' way (default)
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    call ReadRc( rcF, 'output.after.step', output_after_step, status, default = 'd' )
    IF_ERROR_RETURN(status=1)
    if (status == 0) then
       write(gol, *) '****************************************' ; call goPr
       write(gol, *) 'Output will be collected after step :  ',output_after_step ; call goPr
       write(gol, *) '****************************************' ; call goPr
    elseif ( status < 0) then
       write(gol, *) '**********************************************' ; call goPr
       write(gol, *) 'Output will be collected in the OLD way'  ; call goPr
       write(gol, *) 'consider to include e.g. output.after.step : v' ; call goPr
       write(gol, *) '**********************************************' ; call goPr
    endif
    ! close rcfile:
    call Done( rcF, status  )
    IF_NOTOK_RETURN(status=1)
    
    !
    ! timing
    !

    ! define ...
    call GO_Timer_Def( itim_advectx    , 'advectx'     , status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_advecty    , 'advecty'     , status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_advectz    , 'advectz'     , status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_vertical   , 'vertical'    , status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_chemistry  , 'chemistry'   , status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_source_sink, 'source_sink' , status )
    IF_NOTOK_RETURN(status=1)
    
    !
    ! done
    !

    ! end
    call goLabel()
    
    ! ok
    status = 0

  end subroutine Proces_Init
  
  
  ! ***
  
     
  subroutine Proces_Done( status )
  
#ifndef without_chemistry
    use Chemistry     , only : Chemistry_Done
#endif
#ifndef without_advection
    use Advect        , only : Advect_Done
#endif
#ifndef without_convection
    use Convection    , only : Convection_Done
#ifndef without_diffusion
    use Diffusion     , only : Diffusion_Done
#endif
#endif
    use Sources_Sinks , only : Sources_Sinks_Done
#ifndef without_dry_deposition
    use Dry_Deposition, only : Dry_Deposition_Done
#endif
#ifndef without_wet_deposition
    use Wet_Deposition, only : Wet_Deposition_Done
#endif
    
    ! --- in/out --------------------------------
    
    integer, intent(out)    ::  status
    
    ! --- const -------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/Proces_Done'
    
    ! --- begin ---------------------------------
    
#ifndef without_advection
    call Advect_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif
    
#ifndef without_convection
    ! done with convec:
     call Convection_Done( status )
    IF_NOTOK_RETURN(status=1)
#ifndef without_diffusion
    ! done with diffusion:
    call Diffusion_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif
#endif

#ifndef without_chemistry
    ! done with chemistry:
    call Chemistry_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif

    ! done with sources and sinks:    
    call Sources_Sinks_Done( status )
    IF_NOTOK_RETURN(status=1)
    
#ifndef without_dry_deposition
    ! done with dry depos
    call Dry_Deposition_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif
    
#ifndef without_wet_deposition
    ! done with wet depos
    call Wet_Deposition_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif
        
    ! ok
    status = 0

  end subroutine Proces_Done
  
  
  ! ***
  
  
  !
  ! update process fields, for example because meteo has been changed
  ! 
     
  subroutine Proces_Update( status )
  
    use dims          , only : nregions
#ifndef without_diffusion
    use Diffusion     , only : Calc_Kzz, diffusion_write, read_diffusion, write_diffusion
    use toolbox       , only : escape_tm
#endif
    use sources_sinks , only : trace_after_read
#ifndef without_dry_deposition
    use dry_deposition, only : dd_surface_fields
#endif
#ifndef without_wet_deposition
    use wet_deposition, only : calculate_lsp_scav
#endif
  
    ! --- in/out --------------------------------
    
    integer, intent(out)    ::  status
    
    ! --- const -------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/Proces_Update'
    
    ! --- local --------------------------------
    
    integer         ::  region
    
    ! --- begin ---------------------------------
    
#ifndef without_diffusion
    ! recalculate diffusion coeffiecents if necessary

    ! If diffusion_write, then read existing DKG file or
    ! write it if it does not already exist.

    if (diffusion_write) then
       call read_diffusion( status )
       if ( status > 0 ) then
          ! error reading                                                                                                                                                                                                                       
          call escape_tm('Tracer/calc_kzz: error reading diffusion')
       else if ( status < 0 ) then
          ! dkg file could not be found --> calculate                                                                                                                                                                                           
          call Calc_Kzz
          ! write to file                                                                                                                                                                                                                       
          call write_diffusion( status )
          if ( status /= 0 ) call escape_tm('Tracer: error writing diffusion')
       end if
    else
       call Calc_Kzz
    endif
#endif

#ifndef without_dry_deposition
    ! calculate the dry_deposition fields from stuff in fsurf file (will coarsen...)
    call dd_surface_fields( status )
    IF_NOTOK_RETURN(status=1)
#endif

#ifndef without_wet_deposition
    ! pre-calculate the lsp scavenging rloss1, rloss2, rloss3 
    do region = 1, nregions
       call calculate_lsp_scav(region)
    end do
#endif

    ! recalculate chemistry fields if necessary
    call trace_after_read( status )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0

  end subroutine Proces_Update
  
  
  ! ***
  
  
  !
  ! for a given region, performs tref(region)/tref(parent(region)) 
  ! advection steps for each direction x,y,z
  ! including all region's children, granchildren, etc. 
  ! written by patrick berkvens and mike botchev, march-june 1999
  !WP! Also chemistry and sources_sinks are called during advection, 
  !WP! wouter peters, july 2000
  !
  !example (cmk)
  !   XYZ                              ECCEZYX
  !      xyz          eccezyx                  cezyx          xyz          ec
  !         xyzeccezyx       ecxyzzyxec             cezyxxyzec   zyxeccexyz
  ! in this case the operations e and c should be executed in the 
  ! interface region
  ! but not in the core of the zoom region (otherwise double counting)
  ! This results in the most simple algorithm, because you should leave
  ! the DO_STEPS routine AFTER 
  ! (1) finishing all the steps...
  ! (2) XYZ
  ! MK 2000...
  !

  recursive subroutine Proces_Region( region, tr, status )

    use GO          , only : TDate, NewDate
    use dims        , only : tref, revert, ndyn, itaur
    use dims        , only : parent, children
    use dims        , only : statusr => status
    use dims        , only : n_operators
    use datetime    , only : tau2date
    use user_output , only : user_output_step
    use zoom_tools  , only : update_parent
    use advect_tools, only : m2phlb, m2phlb1
    use tracer_data , only : Par_Check_Domain

    ! --- in/out -------------------------------
    
    integer, intent(in)   ::  region
    type(TDate)           ::  tr(2)
    integer, intent(out)  ::  status

    ! --- const -------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/proces_region'
    
    ! --- local ---------------------------------
    
    integer :: child, i, ichild, tref_, n_children, dtime

    integer       ::  time6(6)
    type(TDate)   ::  tri(2)

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

    !main time step for all processes. note that all timesteps are default/2
    dtime = revert*ndyn/(2*tref(region))
    !e.g. xyzvec  (firts call) and cevzyx (second call)

    do i = 1, tref_

       ! time step: itaur(region)+[0,dtime]
       call tau2date( itaur(region), time6 )
       tri(1) = NewDate( time6=time6 )
       call tau2date( itaur(region)+dtime, time6 )
       tri(2) = NewDate( time6=time6 )

       call do_steps( region, tri, status )
       IF_NOTOK_RETURN(status=1)

       ! CALL advect_region for all the children (IF there are any)

       !
       ! children(region,0                   )  : number of children
       ! children(region,1:children(region,0))  : list with region numbers
       !
       ! example for glb6x4/eur3x2/eur1x1 :
       !   children(1,:)= 1, 2 0 0
       !   children(2,:)= 1, 3 0 0
       !   children(3,:)= 0, 0 0 0
       !
       ichild = 0
       do while ( ichild < children(region,0) )
          ichild = ichild + 1
          child = children(region,ichild)
          call proces_region( child, tri, status )
          IF_NOTOK_RETURN(status=1)
       end do
       ! do the remaining steps if necessary...!
       if ( mod(statusr(region),n_operators) /= 0 ) then
          call do_steps( region, tri, status )
          IF_NOTOK_RETURN(status=1)
       end if
       
       ! update region time.....
       itaur(region) = itaur(region) + dtime

       !>>> moved back here (MK)
       !>>> also moved to de_steps for specific sampling.....
       !! Accumulate output, write info at stations
       !! Note: the edges in the interface may not have been updated: 
       !!   this is the task of the parent.
       if (output_after_step == 'd') then    !old default
          call Par_Check_domain(region,'n','tracer') 
          call user_output_step( region, status )   
          IF_NOTOK_RETURN(status=1)
       endif
       !<<<

    end do

    ! update parent using child data ?
    ! compute pressure grid from (changed) air mass:
    if ( region /= 1 ) then
       call update_parent( region )
       call m2phlb( region )
    else
       call m2phlb1( 1 )
    end if

    !if(okdebug) call Par_Check_mass(region,'after proces_region')
    
    ! ok
    status = 0

  end subroutine proces_region


  ! ***
  

  !-------------------------------------------------------------
  ! performs next three (x-, y- and z-) advection steps
  ! written by patrick berkvens and mike botchev, march-june 1999
  !wp! added the call to chemistry. subroutine now does four 
  !wp! steps (c-x-y-z) wp, july 2000
  !mk! changed the order for proper processing of the interfaces.
  !mk! feb-2001
  !-------------------------------------------------------------

  subroutine do_steps( region, tr, status )

    use GO            , only : TDate
    use GO            , only : GO_Timer_Start, GO_Timer_End
    use dims          , only : statusr => status
    use dims          , only : okdebug
    use dims          , only : splitorderzoom, n_operators
#ifndef without_chemistry
    use chemistry     , only : chemie
    use sources_sinks , only : source2
#endif
    use sources_sinks , only : source1
#ifdef with_budgets
    use budget_global , only : budget_transportg
#endif
    use toolbox       , only : escape_tm
#ifdef MPI 
    use mpi_const     , only : previous_par
#endif
    use ParTools      , only : myid, root
    use ParTools      , only : ntracetloc, ntracet_ar
    use tracer_data   , only : Par_Check_Domain
    use global_data   , only : mass_dat
    use chem_param    , only : names, ntracet
#ifndef without_convection
    use convection    , only : convec
#endif
#ifndef without_advection
    use advectx       , only : advectxzoom, put_xedges
    use advecty       , only : advectyzoom, put_yedges
    use advectz       , only : advectzzoom
#endif
#ifndef without_wet_deposition
    use wet_deposition, only : lspscav
#endif
#ifdef with_tendencies
    use tracer_data, only : PLC_Set, plc_ntr, plc_ipr_tchem
#endif
    use user_output , only : user_output_step

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

    integer, intent(in)   ::  region
    type(TDate)           ::  tr(2)
    integer, intent(out)  ::  status

    ! --- const -------------------------------
    
    character(len=*), parameter  ::  rname = mname//'/do_steps'
    
    ! --- local -----------------------------------
    
    integer          :: child, i123, ichild, tref_child, reg, rgi
    character        :: tobedone
    character(len=1) :: next_step, prev_step   !cmk
    integer          :: start_x, stop_x, start_y, stop_y, n !CMKTEMP
#ifdef with_tendencies
    integer          ::  itr
#endif

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

    if ( okdebug .and. myid==root ) print *,'do_steps: region ',region

    !example (cmk)
    !   XYZ          EC CEX
    !      xyzeccezyx      cexyzzyxec
    ! in this case the operations e and c should be executed in 
    ! the interface region but not in the core of the zoom region 
    ! (otherwise double counting)
    ! This results in the most simple algorithm, because you should leave
    ! the DO_STEPS routine AFTER a Z OR after finishing all steps...

    ! note that the parent is responsible for the interface with the children
    ! this means that the interface may not be 'update' in case of (e.g.)
    ! xyz        a
    !    xyzaazyx
    ! in this case, the edges of the child are not updated with the proces a.
    ! this has consequences for the mmix and save output....
    ! THE INTERFACE IS PART OF THE PARENT, NOT THE CHILD.....

    do i123=1,n_operators  !WP! changed 3 to n_operators

      next_step = splitorderzoom(region,statusr(region)+1)
      if ( statusr(region) /= 0 ) then
        prev_step = splitorderzoom(region,statusr(region))
      else
        prev_step = ' '
      end if

      if ( okdebug .and. myid == root ) 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)
          print *, 'do_steps: ',reg,': ', &
               splitorderzoom(reg,1:statusr(reg)),tobedone
        end do

      end if

      select case(next_step)

        case ( 'v' )
        
          call GO_Timer_Start( itim_vertical, status )
          IF_NOTOK_RETURN(status=1)

#ifndef without_convection

          !WP! data must be on tracers for convection 
          call Par_Check_domain(region,'n','tracer') 
          !if ( okdebug ) call Par_Check_mass(region,'bef_convec')

#ifdef with_budgets
          call budget_transportg(region,0,'conv',prev_step)
#endif

          !if ( okdebug ) print *, 'do_steps: Call convection'
          call convec(region, status )
          IF_NOTOK_RETURN(status=1)

#ifdef with_budgets
          call budget_transportg(region,1,'conv',prev_step)
#endif

          !if ( okdebug ) call Par_Check_mass(region,'aft_convec')

#endif

          call GO_Timer_End( itim_vertical, status )
          IF_NOTOK_RETURN(status=1)

        case ( 'c' )

          call GO_Timer_Start( itim_chemistry, status )
          IF_NOTOK_RETURN(status=1)

#ifndef without_chemistry

          !WP! data must be on levels for chemistry 
          call Par_Check_domain(region,'n','levels')
          !if ( okdebug ) call Par_Check_mass(region,'bef_chem')

#ifdef with_tendencies    
          ! store concentrations before chemistry in plc process '0'
          do itr = 1, plc_ntr
            call PLC_Set( 'fill-tracer', itr, 0, status )
            IF_NOTOK_RETURN(status=1)
          end do
#endif

          if ( okdebug ) print*, 'do_steps: Call chemie'
          call Chemie( region, tr, status )
          IF_NOTOK_RETURN(status=1)

          !if ( okdebug ) call Par_Check_mass(region,'aft_chem')

#ifdef with_tendencies    
          ! add change due to chemistry to 'tchem'
          do itr = 1, plc_ntr
            call PLC_Set( 'add-change', itr, plc_ipr_tchem, status, ipr2=0, t=tr(2) )
            IF_NOTOK_RETURN(status=1)
          end do
#endif

#endif

          call GO_Timer_End( itim_chemistry, status )
          IF_NOTOK_RETURN(status=1)

        case( 's' ) 

          call GO_Timer_Start( itim_source_sink, status )
          IF_NOTOK_RETURN(status=1)

          !WP! data must be on tracers for convection 
          call Par_Check_domain(region,'n','tracer')
          !if ( okdebug ) call Par_Check_mass(region,'bef_source')

          if ( okdebug ) print*, 'do_steps: Call source'

          call source1( region, tr, status )
          IF_NOTOK_RETURN(status=1)

#ifndef without_chemistry
          call source2( region, tr, status ) 
          IF_NOTOK_RETURN(status=1)
#endif

#ifndef without_wet_deposition
          ! wet removal by lsp; note that dry deposition is applied in chemistry
          call lspscav( region )
#endif

          !if ( okdebug ) call Par_Check_mass(region,'aft_source')

          call GO_Timer_End( itim_source_sink, status )
          IF_NOTOK_RETURN(status=1)

        case( 'x' )

          call GO_Timer_Start( itim_advectx, status )
          IF_NOTOK_RETURN(status=1)

#ifndef without_advection

          !WP! data must be on tracers for X
          call Par_Check_domain(region,'c','tracer')
          !if ( okdebug ) call Par_Check_mass(region,'befX_tracers')

#ifdef with_budgets
          call budget_transportg(region,0,'advx',prev_step)
#endif

          call advectxzoom(region)

#ifdef with_budgets
          call budget_transportg(region,1,'advx',prev_step)
#endif

          !if ( okdebug ) call Par_Check_mass(region,'aftX_tracers')

#endif

          call GO_Timer_End( itim_advectx, status )
          IF_NOTOK_RETURN(status=1)

        case( 'y' )

          call GO_Timer_Start( itim_advecty, status )
          IF_NOTOK_RETURN(status=1)

#ifndef without_advection

          !WP! data must be on tracers for Y
          call Par_Check_domain(region,'c','tracer')
          !if ( okdebug ) call Par_Check_mass(region,'befY_tracers')

#ifdef with_budgets
          call budget_transportg(region,0,'advy',prev_step)
#endif

          call advectyzoom(region)

#ifdef with_budgets
          call budget_transportg(region,1,'advy',prev_step)
#endif

          !if ( okdebug ) call Par_Check_mass(region,'aftY_tracers')

#endif

          call GO_Timer_End( itim_advecty, status )
          IF_NOTOK_RETURN(status=1)

        case( 'z' )

          call GO_Timer_Start( itim_advectz, status )
          IF_NOTOK_RETURN(status=1)

#ifndef without_advection

          !WP! data must be on tracers for Z
          call Par_Check_domain(region,'n','tracer')
          !if ( okdebug ) call Par_Check_mass(region,'befZ_tracers')

#ifdef with_budgets
          call budget_transportg(region,0,'advz',prev_step)
#endif

          call advectzzoom(region)

#ifdef with_budgets
          call budget_transportg(region,1,'advz',prev_step)
#endif

          !if ( okdebug ) call Par_Check_mass(region,'aftZ_tracers')

#endif

          call GO_Timer_End( itim_advectz, status )
          IF_NOTOK_RETURN(status=1)

        case default

          print *,'do_steps:  strange value in splitorderzoom: ', next_step
          print *,'do_steps:  (must be c, e, v, x, y or z)'
          call escape_tm('do_steps: Error') 

      end select
      
      statusr(region) = statusr(region)+1

      ! e.g. xyzvec...........cevzyx
      ! exit after xyz or vec of cevzyx....
      ! if the exit is after vec...the interface of child should be updated...
      ! 
      if ( mod(statusr(region),n_operators) == 0 ) then
        ! after xyz....vec...update child interface CMK removed feb 2003...
        ! BUT replaced in march 2003 by MK: budgets were corrupted when
        ! the edges are NOT updated!!!  still to find out why???
#ifndef without_advection
        if ( next_step /= 'x' ) then 
          call put_xedges(region)
          call put_yedges(region)
        end if    ! if not .... zyxcev
#endif
        exit ! e.g. vceecvzyx 
      end if

      ! flexible sample scheme ('a' means after all steps), d = old 'default',
      ! output_after_step = v,c,e,x,y,z,d,a
      if ((output_after_step == next_step) .or. (output_after_step == 'a')) then
         call Par_Check_domain(region,'n','tracer') 
         call user_output_step( region, status )   
         IF_NOTOK_RETURN(status=1)
      endif

      !exit after 'yz' sequence to proces xyz children
      if ( next_step == 'z' ) then
         prev_step = splitorderzoom(region,statusr(region)-1)
         if ( prev_step == 'y' ) exit
      end if

    end do
    
    ! ok
    status = 0

  end subroutine do_steps

 
  ! ***
  

  character function upper(xyz)

    implicit none
    
    character(1),intent(in) :: xyz

    if (xyz=='x') then
       upper = 'X'
    else if (xyz=='y') then
       upper = 'Y'
    else if (xyz=='z') then
       upper = 'Z'
    else if (xyz=='c') then
       upper = 'C'
    else if (xyz=='s') then
       upper = 'S'
    else if (xyz=='v') then
       upper = 'V'
    else
       upper = '_'
    end if

  end function upper
 

end module modelIntegration
