ecearth3/sources/tm5mp/base/modelIntegration.F90

!
#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:      MODELINTEGRATION
!
! !DESCRIPTION: handles time ordering of model processes
!\\
!\\
! !INTERFACE: 
!
MODULE MODELINTEGRATION
  !
  ! !USES:
  !
  use GO, only : gol, goPr, goErr, goLabel

  IMPLICIT NONE
  PRIVATE
  !
  ! !PUBLIC MEMBER FUNCTIONS:
  !
  public :: Proces_Init, Proces_Done, Proces_Update
  public :: Proces_Region
  !
  ! !PRIVATE DATA MEMBERS:
  !
  character(len=1)             ::  output_after_step
  character(len=*), parameter  ::  mname = 'ModelIntegration'    ! module name
  !
  ! timer handles
  integer ::  itim_advectx,  itim_advecty,   itim_advectz
  integer ::  itim_vertical, itim_chemistry, itim_source_sink
  !
  !
  ! !REMARKS:
  !
  !EOP
  !------------------------------------------------------------------------

CONTAINS

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PROCES_INIT
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PROCES_INIT( status )
    !
    ! !USES:
    !
    use GO,          only : GO_Timer_Def
    use GO,          only : TrcFile, Init, Done, ReadRc
    use user_output, only : User_output_Init    
    use Binas,       only : p_global
    use dims,        only : nregions, nregions_all, iglbsfc
    use global_data, only : rcfile
    use Meteodata,   only : sp_dat,  phlb_dat, m_dat
    use Meteodata,   only : sp1_dat, sp2_dat,  tsp_dat, spm_dat
    use Meteodata,   only : mfu_dat, mfv_dat,  mfw_dat
    use Meteodata,   only : temper_dat, humid_dat
    use Meteodata,   only : entu_dat, entd_dat, detu_dat, detd_dat
    use Meteodata,   only : lwc_dat,  iwc_dat,  cc_dat,   cco_dat, ccu_dat
    use Meteodata,   only : gph_dat,  omega_dat
    use Meteodata,   only : oro_dat,  lsmask_dat
    use Meteodata,   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
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)    ::  status
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    type(TrcFile) ::  rcF
    integer       ::  n

    ! --- BEGIN ---------------------------------

    call goLabel(rname)  
    !
    ! (1) meteo always used
    !
    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
    !
    ! (2) Init processes (these inits should select the appropriate meteo to use)
    !
#ifndef without_advection
    call Advect_Init( status )
    IF_NOTOK_RETURN(status=1)
#endif

#ifndef without_convection
    call Convection_Init( status )
    IF_NOTOK_RETURN(status=1)
#ifndef without_diffusion
    call Diffusion_Init( status )
    IF_NOTOK_RETURN(status=1)
#endif
#endif

#ifndef without_chemistry
    call Chemistry_Init( status )
    IF_NOTOK_RETURN(status=1)
#endif

    call Sources_Sinks_Init( status )
    IF_NOTOK_RETURN(status=1)

#ifndef without_dry_deposition
    call Dry_Deposition_Init( status )
    IF_NOTOK_RETURN(status=1)
#endif

#ifndef without_wet_deposition
    call Wet_Deposition_Init( status )
    IF_NOTOK_RETURN(status=1)
#endif

    call user_output_init( status )
    IF_NOTOK_RETURN(status=1)    
    !
    ! (3) Check for (and set-to-use) required meteo
    !
    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: not here - set it in process init that really need it
    !    ! call Set(    oro_dat(iglbsfc), status, used=.true. )
    !    !
    !    !<<<

    !
    ! (4) Allocate used meteo
    ! 
    call Meteo_Alloc( status )
    IF_NOTOK_RETURN(status=1)
    !
    ! (5) 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
    !
    ! (6) When do we sample/accumulate output?
    !    
    !      - x,y,z,c,v,e   : always after this process
    !      - a             : after all steps (testing, not recommended)
    !      - d             : the 'old' way (after all, outside do_steps)
    !      - Recommended value: 'v' (default)
    ! 
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    call ReadRc( rcF, 'output.after.step', output_after_step, status, default = 'v' )
    IF_ERROR_RETURN(status=1)
    if (output_after_step == 'd') 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
    else
       write(gol, *) ' ****************************************'                  ; call goPr
       write(gol, *) ' Output will be collected after step :  ',output_after_step ; call goPr
       write(gol, *) ' ****************************************'                  ; call goPr
    endif

    call Done( rcF, status  )
    IF_NOTOK_RETURN(status=1)
    !
    ! (7) timing
    !
    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)
    !
    ! (8) done
    !
    call goLabel()
    status = 0

  END SUBROUTINE PROCES_INIT
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PROCES_DONE
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PROCES_DONE( status )
    !
    ! !USES:
    !
#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
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)    ::  status
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    ! --- BEGIN ---------------------------------

#ifndef without_advection
    call Advect_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif

#ifndef without_convection
    call Convection_Done( status )
    IF_NOTOK_RETURN(status=1)
#ifndef without_diffusion
    call Diffusion_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif
#endif

#ifndef without_chemistry
    call Chemistry_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif

    call Sources_Sinks_Done( status )
    IF_NOTOK_RETURN(status=1)

#ifndef without_dry_deposition
    call Dry_Deposition_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif

#ifndef without_wet_deposition
    call Wet_Deposition_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif

    status = 0

  END SUBROUTINE PROCES_DONE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PROCES_UPDATE
  !
  ! !DESCRIPTION: update process fields, for example because meteo has been changed
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PROCES_UPDATE( status )
    !
    ! !USES:
    !
    use dims          , only : nregions
#ifndef without_diffusion
    use Diffusion     , only : Calc_Kzz, diffusion_write, read_diffusion, write_diffusion
#endif
    use sources_sinks , only : ss_after_read_meteo_update
#ifndef without_dry_deposition
    use dry_deposition, only : dd_surface_fields
#endif
#ifndef without_wet_deposition
    use wet_deposition, only : calculate_lsp_scav
#endif
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)    ::  status
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

    character(len=*), parameter  ::  rname = mname//'/Proces_Update'
    integer ::  region

    ! --- BEGIN ---------------------------------

#ifndef without_diffusion
    ! Read or recalculate diffusion coefficients 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_ERROR_RETURN(status=1)  ! if error reading

       ! else if dkg file could not be found --> calculate and write
       if (status < 0) then

          call Calc_Kzz(status)
          IF_NOTOK_RETURN(status=1)

          call write_diffusion( status )
          IF_NOTOK_RETURN(status=1)

       end if
    else
       call Calc_Kzz(status)
       IF_NOTOK_RETURN(status=1)
    end if
#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

    ! update sources and sinks that depends on meteo
    call ss_after_read_meteo_update( status )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0

  END SUBROUTINE PROCES_UPDATE
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    PROCES_REGION
  !
  ! !DESCRIPTION: for a given region, performs tref(region)/tref(parent(region)) steps.
  !
  !
  !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
  !  
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE PROCES_REGION( region, tr, status )
    !
    ! !USES:
    !
    use GO,           only : TDate, NewDate
    use dims,         only : tref, revert, ndyn, itaur, parent
    use dims,         only : statusr => status
    use dims,         only : n_operators
    use datetime,     only : tau2date
    use user_output,  only : user_output_step
    use advect_tools, only : m2phlb1
    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)   ::  region
    type(TDate)           ::  tr(2)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)  ::  status
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer     :: i, tref_, 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))

    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)

       ! 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

       ! Accumulate or Sample various output (note: used to be the default output_after_step)
       if (output_after_step == 'd') then
          call user_output_step( region, status )   
          IF_NOTOK_RETURN(status=1)
       endif

    end do

    ! compute pressure grid from (changed) air mass:
    call m2phlb1

    ! ok
    status = 0

  END SUBROUTINE PROCES_REGION
  !EOC

  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:    DO_STEPS
  !
  ! !DESCRIPTION: 
  !\\
  !\\
  ! !INTERFACE:
  !
  SUBROUTINE DO_STEPS( region, tr, status )
    !
    ! !USES:
    !
    use GO,             only : TDate
    use GO,             only : GO_Timer_Start, GO_Timer_End
    use dims,           only : statusr => status
    use dims,           only : okdebug, splitorderzoom, n_operators
#ifndef without_chemistry
    use chemistry,      only : chemie
#endif
    use sources_sinks,  only : sources_sinks_apply
#ifdef with_budgets
    use budget_global,  only : budget_transportg
#endif
    use chem_param,     only : names, ntracet
#ifndef without_convection
    use convection,     only : convec
#endif
#ifndef without_advection
    use advectx,        only : advectxzoom
    use advecty,        only : advectyzoom
    use advectz,        only : dynamw
#endif
#ifndef without_wet_deposition
    use wet_deposition, only : lspscav
#endif
    use user_output,    only : user_output_step

!DBG     ! DEBUG STATEMENT
!DBG     use ParTools,       only : isRoot
!DBG     use tracer_data,    only : Par_Check_mass,  tracer_print ! debug
!DBG !    use restart,        only : Restart_Write                 ! debug

    !
    ! !INPUT PARAMETERS:
    !
    integer, intent(in)   ::  region
    type(TDate)           ::  tr(2)
    !
    ! !OUTPUT PARAMETERS:
    !
    integer, intent(out)  ::  status
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    integer          :: i123, reg, rgi
    character        :: tobedone
    character(len=1) :: next_step, prev_step   !cmk
    integer          :: start_x, stop_x, start_y, stop_y, n !CMKTEMP

    ! --- BEGIN -----------------------------------

    !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

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

       !DBG      ! DEBUG STATEMENT
       !DBG      if ( okdebug ) then
       !DBG         if (isRoot) then
       !DBG            ! print current step with capital letter (X,Y or Z)
       !DBG            do reg=1,region
       !DBG               tobedone = ' '
       !DBG               if ( reg == region ) tobedone = upper(next_step)
       !DBG               write(gol,*) ' do_steps: ',reg,': ', splitorderzoom(reg,1:statusr(reg)),tobedone
       !DBG               call goPr
       !DBG            end do
       !DBG         end if
       !DBG 
       !DBG         call tracer_print(region, "do_steps befor "//next_step, status)
       !DBG         IF_NOTOK_RETURN(status=1)
       !DBG         call Par_Check_mass(region, "bef_"//next_step, debug=.true.)
       !DBG      end if

       SELECT CASE(next_step)

       CASE ( 'v' )

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

#ifndef without_convection
#ifdef with_budgets
          call budget_transportg(region,0,'conv',prev_step)
#endif
          call convec(region, status )
          IF_NOTOK_RETURN(status=1)
#ifdef with_budgets
          call budget_transportg(region,1,'conv',prev_step)
#endif
#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
          CALL CHEMIE( region, tr, status )
          IF_NOTOK_RETURN(status=1)
#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)

          call sources_sinks_apply( region, tr, status )
          IF_NOTOK_RETURN(status=1)
 
#ifndef without_wet_deposition
          ! wet removal by lsp (dry deposition is applied in chemistry)
          call lspscav( region )
#endif
          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
#ifdef with_budgets
          call budget_transportg(region,0,'advx',prev_step)
#endif
          CALL ADVECTXZOOM(region, status)
          IF_NOTOK_RETURN(status=1)

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

          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
#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
#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
#ifdef with_budgets
          call budget_transportg(region,0,'advz',prev_step)
#endif
          CALL DYNAMW

#ifdef with_budgets
          call budget_transportg(region,1,'advz',prev_step)
#endif
#endif
          call GO_Timer_End( itim_advectz, status )
          IF_NOTOK_RETURN(status=1)


       CASE default

          write(gol,*)'do_steps:  strange value in splitorderzoom: ', next_step ; call goErr 
          write(gol,*)'do_steps:  (must be c, e, v, x, y or z)'                 ; call goErr 
          TRACEBACK
          status=1; return

       END SELECT

       !DBG       ! DEBUG STATEMENT
       !DBG       if (okdebug) then 
       !DBG          call tracer_print(region, "do_steps after "//next_step, status)
       !DBG          IF_NOTOK_RETURN(status=1)
       !DBG          call Par_Check_mass(region, "aft_"//next_step, debug=.true.)
       !DBG !         call Restart_Write( status, key="aft_"//next_step )
       !DBG !         IF_NOTOK_RETURN(status=1)
       !DBG       end if

       statusr(region) = statusr(region)+1

       ! e.g. xyzvec...........cevzyx
       ! exit after xyz or vec of cevzyx....
       if ( mod(statusr(region),n_operators) == 0 ) then
          exit
       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 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
  !EOC

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

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

END MODULE MODELINTEGRATION