ecearth3/sources/tm5mp/proj/domino/tm5.F90

!#################################################################
!
! TM5 as a library ...
!
!### 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 TM5

  use GO, only : gol, goPr, goErr
  use GO, only : GO_Timer_Init, GO_Timer_Done, GO_Timer_Def, GO_Timer_Start, GO_Timer_End
#ifdef with_prism
  use TM5_Prism, only : comp_name, comp_id
#endif
  
  implicit none
  private
    
  public  ::  TM5_Comm_Init,     TM5_Comm_Done, TM5_Comm_Abort
  public  ::  TM5_Messages_Init, TM5_Messages_Done
  public  ::  TM5_Model_Init,    TM5_Model_Run, TM5_Model_Done
  
#ifdef with_prism
  public  ::  comp_name, comp_id
#endif

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

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

  ! --- var ----------------------------------------
  
  integer :: itim_init, itim_done
  integer :: itim_run_init, itim_run_step, itim_run_done
  integer :: itim_run_init_cfl, itim_run_init_setmass, itim_run_init_setothers, itim_run_init_pup, itim_run_init_ssup
  integer :: itim_run_init_write_restart, itim_start, itim_output
  ! Henk Eskes
  integer :: itim_run_assim

  
contains


  ! ===================================================================
  ! === 
  ! === communication
  ! === 
  ! ===================================================================
  
  !
  ! Setup communication:
  !   o MPI_Init
  !   o fill npes and myid
  !
  
  subroutine TM5_Comm_Init( status, comm )
  
    use ParTools,     only : TM5_MPI_Init
    use OMP_ParTools, only : TM5_OMP_Init
  
    ! --- in/out ----------------------------------
    
    integer, intent(out)            ::  status
    integer, intent(in), optional   ::  comm
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/TM5_Comm_Init'
    
    ! --- begin -----------------------------------
    
    ! setup mpi stuff if necessary:
    call TM5_MPI_Init( status, comm )
    IF_NOTOK_RETURN(status=1)

    ! setup OpenMP stuff if necessary:
    call TM5_OMP_Init( status )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0
    
  end subroutine TM5_Comm_Init


  !
  ! Stop communication:
  !   o MPI_Finalize
  !
  
  subroutine TM5_Comm_Done( status, comm )
  
    use ParTools, only : TM5_MPI_Done
  
    ! --- in/out ----------------------------------
    
    integer, intent(out)            ::  status
    integer, intent(in), optional   ::  comm
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/TM5_Comm_Done'
    
    ! --- begin -----------------------------------

    ! finalize mpi stuff if necessary:
    call TM5_MPI_Done( status, comm )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0
    
  end subroutine TM5_Comm_Done


  !
  ! Abort communication:
  !   o MPI_Abort
  !
  
  subroutine TM5_Comm_Abort( errorcode, status ) 
  
    use ParTools, only : TM5_MPI_Abort
  
    ! --- in/out ----------------------------------
    
    integer, intent(in)             ::  errorcode
    integer, intent(out)            ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/TM5_Comm_Abort'
    
    ! --- begin -----------------------------------

    ! finalize mpi stuff if necessary:
    call TM5_MPI_Abort( errorcode, status )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0
    
  end subroutine TM5_Comm_Abort


  ! ===================================================================
  ! === 
  ! === arguments
  ! === 
  ! ===================================================================
  

  subroutine TM5_Arguments( status )
  
    use GO         , only : goArgCount, goGetArg
    use global_data, only : rcfile
    use partools   , only : isRoot, root, Par_Broadcast_Status, Par_Broadcast

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

    character(len=*), parameter ::  rname = mname//'/TM5_Arguments'
    
    ! --- local -----------------------------------
    
    integer               ::  narg
    integer               ::  iarg
    character(len=1024)   ::  line
    
    ! --- begin -----------------------------------
    
    ! on root only, since some mpirun version do not parse
    ! all arguments to each executable:

    ! number of arguments:
    if (isRoot) call goArgCount( narg, status )
    
    call Par_Broadcast_Status(status, root)
    IF_NOTOK_RETURN(status=1)

    call Par_Broadcast( narg, status )
    IF_NOTOK_RETURN(status=1)

    ! check ...
    if ( narg == 0 ) then
      write (gol,'("no arguments found ...")'); call goErr
      TRACEBACK; status=1; return
    end if

    ! defaults:
    rcfile = 'None'

    ! loop over arguments:
    iarg = 0
    do
      ! next:
      iarg = iarg + 1

      ! get argument:
      if (isRoot) call goGetArg( iarg, line, status )
      
      call Par_Broadcast_Status(status, root)
      IF_NOTOK_RETURN(status=1)

      call Par_Broadcast( line, status )
      IF_NOTOK_RETURN(status=1)

      ! specials ...
      select case ( trim(line) )
        ! arguments added by MPICH/mpirun :
        case ( '-p4pg', '-p4wd' )
          ! skip next argument:
          iarg = iarg + 1
        ! other ...
        case default
          ! not filled yet ?
          if ( trim(rcfile) == 'None' ) then
            rcfile = trim(line)
          else
            write (gol,'("unsupported argument : ",a)') trim(line); call goErr
            TRACEBACK; status=1; return
          end if
      end select
      ! last one is processed now ?
      if ( iarg == narg ) exit
    end do
      
    ! ok
    status = 0
    
  end subroutine TM5_Arguments
  
  
  ! ***
  
  
  subroutine TM5_Print_Usage( status )
  
    ! --- in/out ---------------------------------
    
    integer, intent(out)        ::  status
    
    ! --- begin ----------------------------------
    
    ! display usage line:
    write (*,'("Usage: tm5.x <rcfile>")')
    
    ! ok
    status = 0
    
  end subroutine TM5_Print_Usage
  


  ! ===================================================================
  ! === 
  ! === messages init/done
  ! === 
  ! ===================================================================
  

  subroutine TM5_Messages_Init( status )
  
    use GO         , only : GO_Print_Init, gol, goPr
    use GO         , only : TrcFile, Init, Done, ReadRc
    use partools   , only : npes, myid, root
    use global_data, only : rcfile

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

    character(len=*), parameter ::  rname = mname//'/TM5_Messages_Init'
    
    ! --- local -----------------------------------
    
    type(TrcFile)         ::  rcF
    logical               ::  go_print_all
    logical               ::  go_print_apply
    logical               ::  go_print_trace
    logical               ::  go_print_prompt_pe
    logical               ::  go_print_file
    character(len=256)    ::  go_print_file_base, fname

    ! --- begin -----------------------------------
    
    ! read settings:

    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    call ReadRc( rcF, 'go.print.all', go_print_all, status, default=.false. )
    IF_ERROR_RETURN(status=1)

    call ReadRc( rcF, 'go.print.prompt.pe', go_print_prompt_pe, status, default=npes>1 )
    IF_ERROR_RETURN(status=1)

    call ReadRc( rcF, 'go.print.trace', go_print_trace, status, default=.false. )
    IF_ERROR_RETURN(status=1)

    call ReadRc( rcF, 'go.print.file', go_print_file, status, default=.false. )
    IF_ERROR_RETURN(status=1)
    call ReadRc( rcF, 'go.print.file.base', go_print_file_base, status, default='go.out' )
    IF_ERROR_RETURN(status=1)

    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)

    ! standard output by root only:
    go_print_apply = go_print_all .or. (myid==root)
    
    ! write to file ?
    if ( go_print_file ) then
      if ( myid < 10 ) then
        write (fname,'(a,".",i1.1)') trim(go_print_file_base), myid
      else if ( myid < 100 ) then
        write (fname,'(a,".",i2.2)') trim(go_print_file_base), myid
      else if ( myid < 1000 ) then
        write (fname,'(a,".",i3.3)') trim(go_print_file_base), myid
      else
        write (fname,'(a,".",i6.6)') trim(go_print_file_base), myid
      end if
    else
      fname = 'stdout'
    end if

    ! setup standard output processing:
    call GO_Print_Init( status, &
                          apply=go_print_apply, trace=go_print_trace, &
                          prompt_pe=go_print_prompt_pe, pe=myid, &
                          file=go_print_file, file_name=fname )
    IF_NOTOK_RETURN(status=1)

    ! intro message ...
    write (gol,'(" ")'); call goPr
    write (gol,'("*************************************************************")'); call goPr
    write (gol,'("***                                                       ***")'); call goPr
    write (gol,'("***        Global Atmospheric Tracer Model TM5            ***")'); call goPr
    write (gol,'("***                                                       ***")'); call goPr
    write (gol,'("*************************************************************")'); call goPr
    write (gol,'(" ")'); call goPr

    ! ok
    status = 0
    
  end subroutine TM5_Messages_Init
  

  ! ***
  
    
  subroutine TM5_Messages_Done( status )
  
    use GO, only : GO_Print_Done

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

    character(len=*), parameter ::  rname = mname//'/TM5_Messages_Done'
    
    ! --- begin -----------------------------------
    
    ! final message ...
    write (gol,'(" ")'); call goPr
    write (gol,'("*************************************************************")'); call goPr
    write (gol,'("***                                                       ***")'); call goPr
    write (gol,'("***                 end message log                       ***")'); call goPr
    write (gol,'("***                                                       ***")'); call goPr
    write (gol,'("*************************************************************")'); call goPr
    write (gol,'(" ")'); call goPr

    call GO_Print_Done( status )
    IF_NOTOK_RETURN(status=1)
    
    ! ok
    status = 0
    
  end subroutine TM5_Messages_Done


  ! ====================================================================
  ! ===
  ! === Timing
  ! ===
  ! ====================================================================
  

  subroutine TM5_Timing_Init( status )

    use GO, only : GO_Timer_Init, GO_Timer_Def

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

    integer, intent(inout)              ::  status
    
    ! --- const ----------------------------------
    
    character(len=*), parameter   ::  rname = mname//'/TM5_Timing_Init'
    
    ! --- local ----------------------------------

    ! --- begin ----------------------------------
    
    call GO_Timer_Init( status )
    IF_NOTOK_RETURN(status=1)

    ! define ...
    call GO_Timer_Def( itim_init, 'init', status ) ! MODEL_INIT
    IF_NOTOK_RETURN(status=1)

    ! MODEL_RUN
    call GO_Timer_Def( itim_start, 'step start', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_run_init, 'step init', status )
    IF_NOTOK_RETURN(status=1)    
    call GO_Timer_Def( itim_run_init_cfl,           'step init check cfl', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_run_init_setmass,       'step init set mass', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_run_init_setothers,     'step init set others', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_run_init_pup,           'step init proc update', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_run_init_ssup,          'step init sources update', status )
    IF_NOTOK_RETURN(status=1) 
    call GO_Timer_Def( itim_run_init_write_restart, 'step init write restart', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_run_step, 'step run' , status )
    IF_NOTOK_RETURN(status=1)
    ! Henk Eskes, add step assim
    call GO_Timer_Def( itim_run_assim, 'step assim' , status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_run_done, 'step done', status )
    IF_NOTOK_RETURN(status=1)
    call GO_Timer_Def( itim_output,   'step output', status )
    IF_NOTOK_RETURN(status=1)
    
    call GO_Timer_Def( itim_done, 'done', status ) ! MODEL_DONE
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0
    
  end subroutine TM5_Timing_Init
  
  
  ! ***
  
  
  !--------------------------------------------------------------------------
  !                    TM5                                                  !
  !--------------------------------------------------------------------------
  !BOP
  !
  ! !IROUTINE:  TM5_Timing_Done
  !
  ! !DESCRIPTION: Interface to write profiling output. Get filename and call
  !               timer (profiler).
  !\\
  !\\
  ! !INTERFACE:
  !
  subroutine TM5_Timing_Done( status )
    !
    ! !USES:
    !
    use GO,          only : pathsep
    use GO,          only : TrcFile, Init, Done, ReadRc
    use GO,          only : GO_Timer_Done
    use Global_Data, only : rcfile
    use Partools,    only : myid
    !
    ! !INPUT/OUTPUT PARAMETERS:
    !
    integer, intent(inout)      ::  status
    !
    ! !REVISION HISTORY: 
    !   21 Sep 2010 - P. Le Sager - uses output.dir instead of outputdir
    !                 (to follow pycasso std) 
    !
    ! !REMARKS:
    !
    !EOP
    !------------------------------------------------------------------------
    !BOC

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

    character(len=*), parameter   ::  rname = mname//'/TM5_Timing_Done'
    
    ! --- local ----------------------------------

    integer               ::  l
    character(len=1024)   ::  outdir
    character(len=256)    ::  subdir
    character(len=1024)   ::  timing_file
    type(TrcFile)         ::  rcF
    logical               ::  putout

    ! --- begin ----------------------------------
    
    ! first open rcfile:
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    ! read flag; by default false to avoid problems with uncreated directories etc:
    call ReadRc( rcF, 'timing.output', putout, status, default=.false. )
    IF_ERROR_RETURN(status=1)

    ! putout ?
    if ( putout ) then

      ! output directory:
      call ReadRc( rcF, 'output.dir', outdir, status )
      IF_NOTOK_RETURN(status=1)

      ! timing subdirectory:
      call ReadRc( rcF, 'timing.output.subdir', subdir, status, default='' )
      IF_ERROR_RETURN(status=1)

      ! filename to output time profile:
      l = len_trim(rcfile)
      write (timing_file,'(5a,"_",i4.4,".prf")') &
                              trim(outdir), pathsep, trim(subdir), pathsep, &
                              rcfile(1:l-3), myid

      ! done with timers; write profile to standard output and file:
      call GO_Timer_Done( status, file=trim(timing_file) )
      IF_NOTOK_RETURN(status=1)

    end if  ! putout

    ! close:
    call Done( rcF, status )
    IF_ERROR_RETURN(status=1)

    ! ok
    status = 0
    
  end subroutine TM5_Timing_Done
  !EOC


  ! ===================================================================
  ! === 
  ! === ok file
  ! === 
  ! ===================================================================

  
  ! Write dummy file 'tm5.ok'.
  ! Existence of this file is used by the scripts to check
  ! if a run ended properly.
  ! Checking exit status would be better, but this does
  ! not trap 'stop' statements and other obscure endings.

  subroutine TM5_Write_OkFile( status )

    use GO, only : goGetFU

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

    character(len=*), parameter ::  rname = mname//'/TM5_Write_OkFile'
    
    ! --- local ----------------------------------
    
    integer           ::  fu

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

    ! get free file unit:
    call goGetFU( fu, status )
    IF_NOTOK_RETURN(status=1)
    ! open file:
    open( unit=fu, file='tm5.ok', form='formatted', status='unknown', iostat=status )
    if ( status/=0 ) then
      write (gol,'("from opening okfile")'); call goErr
    else
      ! write happy message:
      write (fu,'("Program terminated normally")',iostat=status)
      if ( status/=0 ) then
        write (gol,'("from writing to okfile")'); call goErr
      else
        ! close:
        close( fu, iostat=status )
        if ( status/=0 ) then
          write (gol,'("from closing okfile")'); call goErr
        end if
      end if
    end if

    ! ok
    status = 0
    
  end subroutine TM5_Write_OkFile


  ! ===================================================================
  ! === 
  ! === model init/done
  ! === 
  ! ===================================================================
  

  subroutine TM5_Model_Init( status )
  
    use GO,               only : TDate, NewDate
    use GO,               only : TrcFile, Init, Done, ReadRc
    use MDF,              only : MDF_Init
    use geometry,         only : calc_dxy, GeomtryH
    use dims,             only : nregions, dxy11, nlat180, okdebug
    use global_data,      only : rcfile, declare_fields, region_dat
    use tracer_data,      only : tracer_print
    use ModelIntegration, only : Proces_Init
    use Meteo,            only : Meteo_Init, Meteo_Init_Grids
    use tm5_distgrid,     only : tm5_dgrid_init
    use redgridZoom,      only : RedGrid_Init
#ifdef with_prism
    use dims            , only : nregions_all,iglbsfc!,iglbsfc_prism
    use MeteoData       , only : global_lli, levi
    use TM5_Prism       , only : TM5_Prism_Init, TM5_Prism_Init2
    use TM5_Prism       , only : ifs_nlev
#endif
#ifdef with_tendencies
    use tracer_data,      only : PLC_Init
#endif
    use initexit,         only : control_init
    use restart,          only : Restart_Init
  
    ! --- in/out ----------------------------------
    
    integer, intent(out)          ::  status
    
    ! --- const ------------------------------

    character(len=*), parameter ::  rname = mname//'/TM5_Model_Init'
    
    ! --- local ----------------------------------
    
    type(TrcFile)      ::  rcF
    integer :: region

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

    ! extract arguments:
    call TM5_Arguments( status )
    if (status/=0) then
      call TM5_Print_Usage( status )
      status=1; return
    end if

    ! setup messages
    call TM5_Messages_Init( status )
    IF_NOTOK_RETURN(status=1)

    write (gol,'(a,":   init model (read control param, init calendar/time) ...")') rname; call goPr

    ! Fisrt, read control parameters, since several are required by many inits
    CALL CONTROL_INIT( status )
    IF_NOTOK_RETURN(status=1)

    ! init parallelisation
    write (gol,'(a,":   init distributed grid ...")') rname; call goPr
    call tm5_DGRID_Init( rcfile, status )
    IF_NOTOK_RETURN(status=1)

    ! init timers:
    write (gol,'(a,":   init timers ...")') rname; call goPr
    call TM5_Timing_Init( status )
    IF_NOTOK_RETURN(status=1)
    
    ! start timing ...
    call GO_Timer_Start( itim_init, status )
    IF_NOTOK_RETURN(status=1)

    ! init MDF interface to HDF/NetCDF:
    call MDF_Init( status )
    IF_NOTOK_RETURN(status=1)

#ifdef with_prism
    ! init prism coupler: read coupling parameter from rc file
    write (gol,'(a,":   init prism ...")') rname; call goPr
    call TM5_Prism_Init( rcfile, status )
    IF_NOTOK_RETURN(status=1)
#endif
    
    ! setup restart:
    write (gol,'(a,":   init restart ...")') rname; call goPr
    call Restart_Init( status )
    IF_NOTOK_RETURN(status=1)
    
    ! setup meteo input:
    write (gol,'(a,":   init grids ...")') rname; call goPr
    call Meteo_Init_Grids( status )
    IF_NOTOK_RETURN(status=1)

#ifdef with_prism
    ! init prism coupler: grids, partition, coupled variables
    write (gol,'(a,":   init prism 2 ...")') rname; call goPr
    call TM5_Prism_Init2( nregions_all, nregions, iglbsfc, global_lli(1:nregions_all), levi, status )
    IF_NOTOK_RETURN(status=1)
#endif

#ifdef with_tendencies    
    ! init concentration, production, loss rates:
    write (gol,'(a,":   init production/loss/chemistry ...")') rname; call goPr
    call PLC_Init( rcfile, status )
    IF_NOTOK_RETURN(status=1)
#endif
    
    ! setup meteo input:
    write (gol,'(a,":   init meteo (be patient) ...")') rname; call goPr
    call Meteo_Init( status )
    IF_NOTOK_RETURN(status=1)


    ! Allocate tracers and "global data" for this run; fill tracers with 0.
    call declare_fields

    ! Fill horizontal geometry variables
    write (gol,'(a,":   horizontal geometry ...")') rname; call goPr
    call calc_dxy(dxy11, nlat180)  ! grid cell area in 1x1 grid
    do region = 1, nregions
       call GeomtryH( region )     ! grid definition (set dims:areag(1) and globa_data:region_dat%dxyp)
    end do
    region_dat(1)%zoomed = 1    ! remains from zooming capabilities
    region_dat(1)%edge = 0


    ! More geometry and extra data if needed: REDuced GRID
    write (gol,'(a,":   init reduced grid ...")') rname; call goPr
    call redgrid_init( 1, status )
    IF_NOTOK_RETURN(status=1)
    
        
    ! Now we can init processes
    write (gol,'(a,":   init processes ...")') rname; call goPr
    call Proces_Init( status )
    IF_NOTOK_RETURN(status=1)

    
    ! Note: this is the earliest that can be called since meteo fields have to be allocated
    if ( okdebug ) then
       call tracer_print(1, "process init", status)
       IF_NOTOK_RETURN(status=1)
    end if

    write (gol,'(a,":   done")') rname; call goPr
    write (gol,'(" ")') ; call goPr

    ! end timing ...
    call GO_Timer_End( itim_init, status )
    IF_NOTOK_RETURN(status=1)

    ! ok
    status = 0
    
  end subroutine TM5_Model_Init
  
  
  ! ***


  subroutine TM5_Model_Done( status )

    use MDF             , only : MDF_Done
    use ModelIntegration, only : Proces_Done
    use Meteo           , only : Meteo_Done, Meteo_Done_Grids
    use tm5_distgrid    , only : tm5_dgrid_done
#ifdef with_prism
    use TM5_Prism       , only : TM5_Prism_Done
#endif
#ifdef with_tendencies
    use tracer_data     , only : PLC_Done
#endif
    use restart         , only : Restart_Done
    use redgridZoom,      only : RedGrid_Done

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

    character(len=*), parameter ::  rname = mname//'/TM5_Model_Done'
    
    ! --- local -----------------------------------

    integer           ::  errstat
    
    ! --- begin -----------------------------------
    
    write (gol,'(a,": start")') rname; call goPr

    ! start timing ...
    call GO_Timer_Start( itim_done, status )
    IF_NOTOK_RETURN(status=1)

    ! done with restart:
    call Restart_Done( status )
    IF_NOTOK_RETURN(status=1)

#ifdef with_tendencies   
    ! done with production/loss rates
    call PLC_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif

#ifdef with_prism    
    ! done  with prism coupler
    call TM5_Prism_Done( status )
    IF_NOTOK_RETURN(status=1)
#endif

    ! do not break on error from the following routines,
    ! to rescue what could be rescued;
    ! by default, return status is ok:
    errstat = 0
    
    ! done processes
    call Proces_Done( status )
    if (status/=0) then; TRACEBACK; errstat=1; end if

    ! close meteo files etc
    call Meteo_Done( status )
    if (status/=0) then; TRACEBACK; errstat=1; end if

    ! close meteo files etc
    call Meteo_Done_Grids( status )
    if (status/=0) then; TRACEBACK; errstat=1; end if

    ! done with MDF interface to HDF/NetCDF:
    call MDF_Done( status )
    if (status/=0) then; TRACEBACK; errstat=1; end if

    ! end timing ...
    call GO_Timer_End( itim_done, status )
    if (status/=0) then; TRACEBACK; errstat=1; end if
    
    ! done with timing ...
    call TM5_Timing_Done( status )
    if (status/=0) then; TRACEBACK; errstat=1; end if

    ! done parallelisation
    call tm5_dgrid_Done( status )
    if (status/=0) then; TRACEBACK; errstat=1; end if

    ! done with standard output:
    call TM5_Messages_Done( status )
    if (status/=0) then; TRACEBACK; errstat=1; end if

    call RedGrid_Done( status )
    if (status/=0) then; TRACEBACK; errstat=1; end if

    ! write dummy file to indicate proper end:
    if ( errstat == 0 ) then
      call TM5_Write_OkFile( status )
      if (status/=0) then; TRACEBACK; errstat=1; end if
    end if

    write (gol,'(a,": end")') rname; call goPr
    write(gol,'(" ")') ; call goPr

    ! return with error status if some routines failed:
    status = errstat
    
  end subroutine TM5_Model_Done


  ! ===================================================================
  ! === 
  ! === model run
  ! === 
  ! ===================================================================
  

  subroutine TM5_Model_Run( status )

    use GO,               only : TrcFile, Init, Done, ReadRc
    use GO,               only : TDate, NewDate, IncrDate, wrtgol
    use GO,               only : rTotal, operator(+), operator(-), operator(>), operator(==)
    use dims,             only : nregions, okdebug
    use dims,             only : region_status => status
    use dims,             only : nread
    use dims,             only : idate, idatee, idatei
    use dims,             only : itau,  itaue, itaur
    use dims,             only : ndyn_max
    use dims,             only : nread, ndyn, nconv, nsrce, nchem
    use dims,             only : revert
    use dims,             only : newsrun, newmonth
    use dims,             only : nread, idate
    use global_data,      only : rcfile
    use ParTools,         only : Par_Barrier, myid
    use datetime,         only : inctime
    use Meteo,            only : Meteo_Setup_Other
    use Meteo,            only : Meteo_Setup_Mass
#ifndef without_advection
    use AdvectM_CFL,      only : Check_CFL, Setup_MassFlow
#endif
    use ModelIntegration, only : Proces_Update, Proces_Region
    use InitExit,         only : Start
    use InitExit,         only : Exitus
    use sources_sinks,    only : ss_monthly_update
    use user_output,      only : user_output_done, user_output_mean
#ifdef with_prism
    use prism_putget    , only : TM5_Prism_Puts!, TM5_Prism_Gets
#endif
#ifdef with_ecearth_optics
    use ecearth_optics     , only : ECEarth_Optics_Step
    use TM5_Prism          , only : exchange_period
#endif
#ifdef with_tendencies
    use tracer_data,       only : plc_reset_period
    use tm5_tendency_eval, only : apply_tendency, reset_tendency
#endif
    use restart,          only : Restart_Save
    use datetime,         only : tau2date,date2tau

    ! Henk Eskes: perform an NO2 retrieval + assimilation step
    use assim_interf_mod,    only : AssimNO2_interface, AssimNO2_interface_init, AssimNO2_interface_done

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

    integer, intent(out)          ::  status

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

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

    ! --- local ----------------------------------

    type(TrcFile)     ::  rcF
    type(TDate)       ::  tread1, tread2
    type(TDate)       ::  tdyn, tend, tr(2)
    logical           ::  isfirst
    integer           ::  nhalf
    logical           ::  this_is_the_end
    logical           ::  check_pressure
    integer           ::  region
    integer           ::  n
    ! CarbonTracker-specific restart quantities; not used unless
    ! the rcfile contains a "jobstep.step" key.
    integer,dimension(6) :: ct_restart_special = (/0,0,0,0,0,0/)
    integer :: ct_itau,jobstep_step
    ! --- begin -----------------------------------

    ! ~~~ rc file settings ~~~
    call GO_Timer_Start( itim_start, status )
    IF_NOTOK_RETURN(status=1)

    write (gol,'(a,": read settings ...")') rname; call goPr

    ! open rcfile:
    call Init( rcF, rcfile, status )
    IF_NOTOK_RETURN(status=1)

    ! ensure that every 'nread' seconds is at the end of a dynamic time step:
    call ReadRc( rcF, 'time.ntimestep', nread, status )
    IF_NOTOK_RETURN(status=1)

    ! a CarbonTracker-specific setting; no error if is it missing
    call ReadRc( rcF, 'jobstep.step' ,jobstep_step, status ,default=0)
    IF_ERROR_RETURN(status=1)

    ! close rcfile:
    call Done( rcF, status )
    IF_NOTOK_RETURN(status=1)

    ! ~~~~~~~~~~~~~~~~~~~~~~~~

    write (gol,'(a,": call start ..")') rname; call goPr

    ! set-up and read user input; 
    ! return time interval for which meteo was read:
    call Start( tread1, tread2, status )
    IF_NOTOK_RETURN(status=1)

    write (gol,'(a,": setup times ..")') rname; call goPr

    ! current time (begin of dynamics step)
    tdyn = NewDate( time6=idate  )
    tend = NewDate( time6=idatee )

    !synchronize time-count regions....
    itaur(:) = itau 

    write (gol,'(a,": start time loop ...")') rname; call goPr
    write (gol,'(" ")'); call goPr

    ! first step in time loop ?
    isfirst = .true.
    nhalf   = 0

    if(jobstep_step .gt. 0) then
       call date2tau(idatei,ct_itau)
       ct_itau=ct_itau + 86400*jobstep_step
       call tau2date(ct_itau,ct_restart_special)
    end if
    
    call GO_Timer_End( itim_start, status )
    IF_NOTOK_RETURN(status=1)

    ! Henk Eskes: Initialise assimilation
    call AssimNO2_interface_init ( rcFile )

    
    ! time loop over steps ndyn/2 (!)
    do 
       !*************************************************************************
       ! *** INIT STEP  ***
       !*************************************************************************
       call GO_Timer_Start( itim_run_init, status )
       IF_NOTOK_RETURN(status=1)

       ! is this the end time ?
       this_is_the_end = (revert*itau) >= (revert*itaue)

       ! next half step
       nhalf = modulo(nhalf,2) + 1

       !
       ! *** write restart ***
       !
       if ( nhalf == 1 ) then
          call GO_Timer_Start( itim_run_init_write_restart, status )
          IF_NOTOK_RETURN(status=1)

          ! eventually save extra restart file, or final save file:
          if (all((ct_restart_special - idate) .eq. 0)) then
             write (gol,'(a,": Writing restart files for CarbonTracker jobstep.step.")') rname; call goPr
             call Restart_Save( status, extra=.false., isfirst=isfirst )
          else
             ! restart file gets written when extra is .false.  
             call Restart_Save( status, extra=(.not. this_is_the_end ), isfirst=isfirst )
          endif
          IF_NOTOK_RETURN(status=1)

          call GO_Timer_End( itim_run_init_write_restart, status )
          IF_NOTOK_RETURN(status=1)
       end if
       !
       ! *** new time interval ? ***
       !
       if ( this_is_the_end ) then

          call GO_Timer_End( itim_run_init, status )
          IF_NOTOK_RETURN(status=1)
          
          exit
       endif

       !
       ! *** Perform assimilation with "ndyn_max" model time step *** (Henk Eskes)
       !
       if ( mod(itau,ndyn_max) == 0 ) then

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

          call wrtgol( '>>> NO2 retrieval-assimilation step for : ', tdyn ); call goPr

          call assimNO2_interface( rcFile, ndyn_max, status )
          IF_NOTOK_RETURN(status=1)

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

       end if

       if ( nhalf == 1 ) then
          call wrtgol( '>>> dynamics step from : ', tdyn ); call goPr
       end if

       !
       ! *** setup data ***
       !

       ! New meteo (if reached end of time interval for which meteo is valid)
       if ( tdyn == tread2 ) then

          ! reset possible reduced timestep due to CFL:
          ndyn =  ndyn_max
          nsrce = ndyn_max
          nconv = ndyn_max
          nchem = ndyn_max

          ! setup meteo data for next interval;
          ! nread is the length (in seconds) of the interval in which 
          ! surface pressure is interpolated (and mass fluxes are constant)
          tread1 = tdyn
          tread2 = tdyn + IncrDate(sec=nread)
          if ( tread2 > tend ) tread2 = tend

          ! n is the number of dynamic intervals within the
          ! time interval for which the meteo has been setup:
          n = ceiling( rTotal(tread2-tread1,'sec') / real(ndyn) )
          ndyn = nint( rTotal(tread2-tread1,'sec') / n )

          ! setup mass and mass fluxes:
          !  o skip first time; already called in 'initexit/start'
          !  o check pressure implied by advection if advection is applied
#ifdef without_advection
          check_pressure = .false.
#else
          check_pressure = .true.
#endif
          call GO_Timer_Start( itim_run_init_setmass, status )
          IF_NOTOK_RETURN(status=1)          
          call Meteo_Setup_Mass( tread1, tread2, status, check_pressure=check_pressure )
          IF_NOTOK_RETURN(status=1)
          call GO_Timer_End( itim_run_init_setmass, status )
          IF_NOTOK_RETURN(status=1)

#ifndef without_advection
          call GO_Timer_Start( itim_run_init_cfl, status )
          IF_NOTOK_RETURN(status=1)
          ! determine dynamic timestep ndyn for this interval [tread1,tread2] ;
          ! the initial number of time steps n is increased until no cfl
          ! violations occure
          call Check_CFL( tread1, tread2, n, status )
          IF_NOTOK_RETURN(status=1)
          call GO_Timer_End( itim_run_init_cfl, status )
          IF_NOTOK_RETURN(status=1)
#endif
       end if

       ! Setup meteo for dynamic step tdyn+[0,ndyn]
       if ( nhalf == 1 ) then
          ! time range of dynamic step:
          tr(1) = tdyn
          tr(2) = tdyn + IncrDate( sec=ndyn )
          
          call GO_Timer_Start( itim_run_init_setothers, status )
          IF_NOTOK_RETURN(status=1)
#ifndef without_advection
          ! convert pu/pv to am/bm/cm, eventually time interpolated
          call Setup_MassFlow( tr, status )
          IF_NOTOK_RETURN(status=1)
#endif
          ! setup (interpolate?) other meteo:
          call Meteo_Setup_Other( tr(1), tr(2), status )
          IF_NOTOK_RETURN(status=1)
          call GO_Timer_End( itim_run_init_setothers, status )
          IF_NOTOK_RETURN(status=1)

          ! Sources and sinks update (must be done after meteo setup [because of the NOx emissions
          ! vertical remapping], and before the Proces_Update [because of the getDMS call down the
          ! line]). Note that we assume that 'newmonth' is always at nhalf==1
          call GO_Timer_Start( itim_run_init_ssup, status )
          IF_NOTOK_RETURN(status=1)          
          if ( newmonth ) then
             call ss_monthly_update( status )
             IF_NOTOK_RETURN(status=1)
          end if
          call GO_Timer_End( itim_run_init_ssup, status )
          IF_NOTOK_RETURN(status=1)
   
          call GO_Timer_Start( itim_run_init_pup, status )
          IF_NOTOK_RETURN(status=1)
          ! recalculate proces dependend fields if necessary
          call Proces_Update( status )
          IF_NOTOK_RETURN(status=1)
          call GO_Timer_End( itim_run_init_pup, status )
          IF_NOTOK_RETURN(status=1)

       end if
              
#ifdef with_ecearth_optics
       if ( (modulo(tdyn%hour,exchange_period)==0) .and. all((/tdyn%min,tdyn%sec,tdyn%mili/)==0) ) then
          call wrtgol( 'Calculating aerosol radiative properties for EC-Earth at ', tdyn ); call goPr

          call ECEarth_Optics_Step ( status )
          IF_NOTOK_RETURN(status=1)

       endif
#endif

#ifdef with_prism
       ! put concentrations to IFS:
       call TM5_Prism_Puts( tdyn, status )
       IF_NOTOK_RETURN(status=1)
#endif

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

       !*************************************************************************
       ! *** RUN STEP (processes) ***
       !*************************************************************************
       call GO_Timer_Start( itim_run_step, status )
       IF_NOTOK_RETURN(status=1)

       if ( ndyn > 0 ) then

          ! reset the process status counters:
          if ( nhalf == 1 ) region_status(1:nregions) = 0

          tr(1) = tdyn
          tr(2) = tdyn + IncrDate(sec=ndyn/2)

          ! info
          if ( okdebug ) then
             if ( nhalf == 1 ) then
                call wrtgol( '--> first  half : ', tr(1), ' - ', tr(2) ); call goPr
             end if
             ! This was needed for LiNox budget computation script, when LiNOx was not written to budget file
             if ( nhalf == 2 ) then
                call wrtgol( '--> second half : ', tr(1), ' - ', tr(2) ); call goPr
             end if
          endif

          itaur(:) = itau       ! synchronize time-count regions

          call Proces_Region( 1, tr, status ) ! start recursive process for the main region = 1
          IF_NOTOK_RETURN(status=1)

       end if

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

       !*************************************************************************
       ! *** DONE STEP (next) ***
       !*************************************************************************
       !
       call GO_Timer_Start( itim_run_done, status )
       IF_NOTOK_RETURN(status=1)

       ! advance the model time with ndyn/2 seconds:
       call inctime
       tdyn = tdyn + IncrDate( sec=nint(ndyn/2.0) )

       ! update mean outputs:
       if ( mod(itau,ndyn_max) == 0) then
          call user_output_mean( status )
          IF_NOTOK_RETURN(status=1)
       end if

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

       ! end first time loop
       isfirst = .false.

    END DO   ! MAIN LOOP 

    ! Finalise assimilation (Henk Eskes)
    call AssimNO2_interface_done()

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

    ! complete user-specified output:
    call user_output_done( status )
    IF_NOTOK_RETURN(status=1)

    ! store save file etc
    call exitus( status )
    IF_NOTOK_RETURN(status=1)
    
    call GO_Timer_End( itim_output, status )
    IF_NOTOK_RETURN(status=1)
    
    write (gol,'(a,": end")') rname; call goPr
    write(gol,'(" ")') ; call goPr

    ! ok
    status = 0

  end subroutine TM5_Model_Run

end module TM5