#include "options.inc"


      program driver
c
c--------------------------------------------------------------
c     drive SPFLAME from here
c     c.eden
c--------------------------------------------------------------
c
      use timing_module
      implicit none
      integer ierr,domains,my_comm,sub_domain,n
      integer,allocatable :: n_pes_domains(:)
      integer,allocatable :: sub_domain_pe0(:)
      logical :: eorun = .false.

      call mpi_init(ierr)
      call mpi_test()

      call tic('setup')
c
c     read the command line
c
      call get_commandline_arg_int(1,domains,ierr)
      if (ierr/=0) goto 999
      allocate(n_pes_domains(domains),sub_domain_pe0(domains))
      do n=1,domains
       call get_commandline_arg_int(n+1,n_pes_domains(n),ierr)
       if (ierr/=0) goto 999
      enddo
c
c     orchestrating all the processors
c
      call split_pes(domains,n_pes_domains,
     &                   my_comm,sub_domain,sub_domain_pe0)
c
c     redirect stdout to files for sub domains
c
      if (sub_domain >0) call redirect_stdout(sub_domain)
c
c     set up the model and do first domain exchange
c
      call setup(domains,sub_domain,sub_domain_pe0,my_comm)
      call setvbc
      call domain_exchg

      call toc('setup')
c
c     here is the main loop
c
  20  continue
      call spflame(eorun)
#ifdef couple
      if (.not.eorun) call sbc_couple
#endif
      if (.not. eorun ) goto 20
c
c     clean up, write restart
c
      call tic('cleanup')
#ifdef couple
c
c     write also restart for OASIS
c
      call sbc_couple_write_restart
#endif
      call write_restart_file
      call bcast_world_logical(eorun,1,0)  ! a barrier
      call toc('cleanup')

      call show_timing()
c
c     and stop
c
      call mpi_finalize(ierr)
      stop
 999  print*,' Usage:  ./program_name "Number of domains"',
     &       ' "Number of PEs for 1. domain"',
     &       ' "Number of PEs for 2. domain" ...'
      print*,' '
      call halt_stop('error in command line arguments')
      end program driver



      subroutine spflame(eorun)
c
c--------------------------------------------------------------
c     main loop for SPFLAME, measure time for the main routines
c--------------------------------------------------------------
c
      use spflame_module

      implicit none
      logical :: eorun 
      integer :: nmax,n
      character (len=32) :: stamp

      call tic('domexchg')
      call nmax_exchg(nmax)
      call toc('domexchg')

      call tic('main')
      do n=1,nmax
 20    call time_step_begin()

       if (enable_freesurf) call calc_thickness

       call tic('setvbc')
       call setvbc
       call toc('setvbc')

       call tic('isopyc')
       if (enable_diffusion_isoneutral) call isoneutral
       if (enable_diffusion_isopycnic)  call isopycnic
       call toc('isopyc')

       call tic('vmixc')
       call vmixc
       call toc('vmixc')
     
       call tic('adv_vel')
       call adv_vel
       call toc('adv_vel')

       call tic('clinic')
       call clinic
       call toc('clinic')

       call tic('tracer')
       call tracer
       call toc('tracer')

       call tic('tropic')
       if (enable_freesurf) then
        call freesurf
       else
        call rigid_lid
       endif
       call toc('tropic')


       call tic('diag')
       call diag
       call toc('diag')


       call time_step_end()
       if (eulerback2) goto 20    ! do the second part of the euler backw. 

      enddo
      call toc('main')

      call tic('domexchg')
      call domain_exchg
      call toc('domexchg')

      eorun=.false.;if (end_of_run) eorun=.true.
c
c     make sure that all pes quit their loops
c     PE 0 of sub domain 0 is the reference
c
      call bcast_world_logical(eorun,1,0) 
      call sub_flush(6)
      end subroutine spflame





      subroutine show_timing()
c
c--------------------------------------------------------------
c     show timing results here
c--------------------------------------------------------------
c
      use spflame_module
      implicit none
      real :: x1,x2,x3,x4
      integer :: npe,iret


#ifdef timing
      do npe = 0,n_pes
       call barrier

      if (my_pe == npe) then
      x1= timing_secs('main')/100.

      print*,''
      print*,'Timing summary for PE #',my_pe, 
     &  ' on sub_domain #',sub_domain

      print*,' MOM loop time summary    = ',
     &   timing_secs('main') ,' s'

      print*,' setup time summary       = ',
     &   timing_secs('setup'),' s'

      print*,' cleanup time summary     = ',
     &   timing_secs('cleanup'),' s'

      print*,' Domain exchange time     = ',
     &   timing_secs('domexchg'),' s',
     &   '=',timing_secs('domexchg')/x1, '%'

      print*,' adv_vel time summary     = ',
     &   timing_secs('adv_vel'),' s',' =',
     &   timing_secs('adv_vel')/x1,'%'

      print*,' clinic time summary      = ',
     &   timing_secs('clinic'),' s',' =',
     &   timing_secs('clinic')/x1,'%'

      print*,' tracer time summary      = ',
     &   timing_secs('tracer'),' s',' =',
     &   timing_secs('tracer')/x1,'%'

      print*,' passive tr. time summary = ',
     &   timing_secs('passive tracer'),' s',' =',
     &   timing_secs('passive tracer')/x1,'%'

      if (enable_freesurf) then
       print*,' free surf time summary  = ',
     &   timing_secs('tropic'),' s',' =',
     &   timing_secs('tropic')/x1,'%'
      else
       print*,' tropic time summary      = ',
     &   timing_secs('tropic'),' s',' =',
     &   timing_secs('tropic')/x1,'%'
      endif

      print*,' isopyc time summary      = ',
     &   timing_secs('isopyc'),' s',' =',
     &   timing_secs('isopyc')/x1,'%'

      print*,' setvbc time summary      = ',
     &   timing_secs('setvbc'),' s',' =',
     &   timing_secs('setvbc')/x1,'%'

      print*,' diag time summary        = ',
     &   timing_secs('diag'),' s',' =',
     &   timing_secs('diag')/x1,'%'

      print*,' vmixc time summary       = ',
     &   timing_secs('vmixc'),' s',' =',
     &   timing_secs('vmixc')/x1,'%'

      print*,' costs for measuring      = ',
     &   timing_secs('tictoc'),' s',' =',
     &   timing_secs('tictoc')/x1,'%'


#ifdef detailed_timing

        print*,''
        print*,'More details for PE #',my_pe
        print*,''
        print*,' for tracer :'
        x1=timing_secs('tracer main')/100.
        print*,' tend  = ',timing_secs('tracer tendency')/x1 
        print*,' convt = ',timing_secs('tracer convect')/x1 
        print*,' ktmix = ',timing_secs('tracer ktmix')/x1 
        print*,' adv   = ',timing_secs('tracer adv')/x1 
        print*,' diff  = ',timing_secs('tracer diff')/x1 
        print*,' mpp   = ',timing_secs('tracer mpp')/x1 
        print*,' obc   = ',timing_secs('tracer obc')/x1 
        print*,' impli = ',timing_secs('tracer implicit')/x1 
c        print*,' c rho = ',timing_secs('convect rho')/x1 
c        print*,' c mix = ',timing_secs('convect mix')/x1 

        x2=0.
        x2=x2+timing_secs('tracer tendency')/x1 
        x2=x2+timing_secs('tracer convect')/x1 
        x2=x2+timing_secs('tracer ktmix')/x1 
        x2=x2+timing_secs('tracer adv')/x1 
        x2=x2+timing_secs('tracer diff')/x1 
        x2=x2+timing_secs('tracer mpp')/x1 
        x2=x2+timing_secs('tracer obc')/x1 
        x2=x2+timing_secs('tracer implicit')/x1 
        print*,' sum = ',x2
        print*,' for clinic :'
        x1=timing_secs('clinic main')/100.
        print*,' dens = ',timing_secs('clinic density')/x1 
        print*,' tend = ',timing_secs('clinic tend')/x1 
        print*,' biha = ',timing_secs('clinic biha')/x1 
        print*,' adv  = ',timing_secs('clinic adv')/x1 
        print*,' diff = ',timing_secs('clinic diff')/x1 
        print*,' obc  = ',timing_secs('clinic obc')/x1 
        print*,' impli= ',timing_secs('clinic impli')/x1 
        if (enable_freesurf) then
         print*,' for free surface :'
        else
         print*,' for tropic :'
         x1=timing_secs('tropic')/100.
         print*,' zu exchg = ',timing_secs('tropic zu exchg')/x1 
         print*,' forcing  = ',timing_secs('tropic forc')/x1 
         print*,' congrad  = ',timing_secs('tropic congr')/x1 
         print*,' phase vel= ',timing_secs('tropic phase')/x1 
         print*,' obc      = ',timing_secs('tropic obc')/x1 
         print*,' add ext  = ',timing_secs('tropic addext')/x1 
        endif
#endif

       endif
       call sub_flush(6)
       call barrier
      enddo

#endif

      end subroutine show_timing

