#include "options.inc"



      program cpflame
c=======================================================================
c      Top level driver for cpflame
c=======================================================================
      use cpflame_module
      implicit none
      integer :: ierr
c---------------------------------------------------------------------------------
c       Initialize MPI setting
c---------------------------------------------------------------------------------
      call mpi_init(ierr)
      call my_mpi_test()
      call get_mpi_comm_world(my_comm)
      call mpi_comm_rank(my_comm,my_pe,ierr)
      call mpi_comm_size(my_comm,n_pes,ierr)
c---------------------------------------------------------------------------------
c      Initialize model
c---------------------------------------------------------------------------------
      call tic('setup')
      call setup()
#ifdef enable_ekecalc
      call init_ekecalc
#endif
#ifdef enable_ekecalc_3D
      call eke3D_init
#endif
#ifdef enable_linear_closure
      call init_linear_closure
#endif
#ifdef enable_bolus_velocity
      call init_bolus
#endif
#ifdef enable_vert_friction_trm
      call init_vert_friction_trm
#endif
#ifdef enable_tke_closure
      call init_tke_closure
#endif
#ifdef enable_diag_particle
      call init_particle
#endif
#ifdef enable_diag_tracer
      call init_tracer
#endif
#ifdef enable_diag_timeave
      call init_time_averages
#endif
#ifdef enable_diag_zonalave
      call init_zonal_averages
#endif
#ifdef enable_diag_spectral
      call init_spectral
#endif
#ifdef enable_diag_numbers
      call init_numbers
#endif
#ifdef enable_diag_over
      call init_overturning
#endif
#ifdef enable_back_stratification
      call init_back_stratification
#endif

      if (my_pe>0) time_manager_verbose = .false.
      call init_time_manager_raw(dt_in,runlen,snap_int,
     &                       restart_time_in,restart_itt )
      call disable_mixing_time_steps
c       mixing_type = 'euler_foreward' 
c       gamma=0

#ifdef enable_press_perturb
      call init_press_perturb
#endif

      call toc('setup')

c---------------------------------------------------------------------------------
c      Begin main model loop
c---------------------------------------------------------------------------------
 10   call time_step_begin()
      call tic('main loop')
      call boundary_conditions
      if (enable_hydrostatic) call vertical_velocity

      call tic('buoyancy')
      call integrate_buoyancy
#ifdef enable_back_stratification
      call back_stratification
#endif
#ifdef enable_diag_tracer
      call integrate_tracer
#endif
      call toc('buoyancy')

      call tic('momentum')
      call momentum_tendency
#ifdef enable_ekecalc
      call ekecalc
#endif
#ifdef enable_ekecalc_3D
      call eke3D
#endif
#ifdef enable_linear_closure
      call linear_closure
#endif
#ifdef enable_bolus_velocity
      call bolus_velocity
#endif
#ifdef enable_tke_closure
      call tke_closure
#endif
#ifdef enable_vert_friction_trm
      call vert_friction_trm
#endif
#ifdef enable_back_stratification
      call back_strat_momentum
#endif
      call toc('momentum')

      call tic('pressure')
      call solve_pressure
#ifdef enable_press_perturb
      call press_perturb
#endif
      call toc('pressure')

      call tic('integrate')
      call integrate
      call toc('integrate')

      call tic('diagnose')
      call diagnose
#ifdef enable_diag_particle
      call integrate_particle
#endif
#ifdef enable_diag_timeave
      call time_averages
#endif
#ifdef enable_diag_zonalave
      call zonal_averages
#endif
#ifdef enable_diag_spectral
      call spectral
#endif
#ifdef enable_diag_numbers
      call diag_numbers
#endif
#ifdef enable_diag_over
      call overturning
#endif
      call toc('diagnose')
      call toc('main loop')
      call tic('cleanup')
      if (last_time_step) call write_restart
      call toc('cleanup')
      call time_step_end()
      if (.not. end_of_run ) goto 10
c---------------------------------------------------------------------------------
c       End main model loop
c---------------------------------------------------------------------------------

      call show_timing()

      call mpi_finalize(ierr)
      end program cpflame



      subroutine read_restart
c ----------------------------------
c       read the restart file
c ----------------------------------
      use cpflame_module
      implicit none
      integer :: j,io,n,js,je

      if (my_pe==0) print*,' reading from restart file '
      call get_free_iounit(io)
      open(io,file='restart.dta',form='unformatted',status='old')
      read(io) restart_itt,restart_time_in
c        restart_itt=0 !!!!
c        call set_calendar_type(no_leap)
c        restart_time_in = set_date(1900,1,1,0,0,0.0)
      if (my_pe==0) then 
        print*,' days since initial time ',restart_time_in%days
        print*,' secs since initial time ',restart_time_in%seconds
        print*,' itt        = ',restart_itt
      endif
      do j=1,jmt
        read(io) u(:,j,:,:,tau), u(:,j,:,:,taum1)
        read(io) b(:,j,:,tau),   b(:,j,:,taum1)
        read(io) eta(:,j,tau),   eta(:,j,taum1), psi(:,j,:)
        read(io) p_full(:,j,:,tau), p_full(:,j,:,taum1) ! this is new
        if (enable_expl_free_surf) then
         read(io) bu(:,j,tau),bv(:,j,tau)
         read(io) bu(:,j,taum1),bv(:,j,taum1)
        endif
      enddo
      close(io)
      if (my_pe==0) print*,' done reading restart'
      end subroutine read_restart


      subroutine write_restart
c ----------------------------------
c      write a restart file
c ----------------------------------
      use cpflame_module
      implicit none
      integer :: j,io,n,js,je

      if (my_pe==0) then
       print*,' writing restart file '
       call get_free_iounit(io)
       open(io,file='restart.dta',form='unformatted',status='unknown')
       write(io) itt,current_time
       do j=js_pe,je_pe
        write(io) u(:,j,:,:,tau), u(:,j,:,:,taum1)
        write(io) b(:,j,:,tau),   b(:,j,:,taum1)
        write(io) eta(:,j,tau),   eta(:,j,taum1), psi(:,j,:)
        write(io) p_full(:,j,:,tau), p_full(:,j,:,taum1) ! this is new
        if (enable_expl_free_surf) then
         write(io) bu(:,j,tau),  bv(:,j,tau)
         write(io) bu(:,j,taum1),  bv(:,j,taum1)
        endif
       enddo
      endif

      do n=1,n_pes-1
        if (my_pe==n) then
         call send_integer(js_pe,1,0,0)
         call send_integer(je_pe,1,0,0)
         do j=js_pe,je_pe
          call send_real(u(:,j,:,:,tau  ),imt*km*3,0,0)
          call send_real(u(:,j,:,:,taum1),imt*km*3,0,0)
          call send_real(b(:,j,:,tau  )  ,imt*km  ,0,0)
          call send_real(b(:,j,:,taum1)  ,imt*km  ,0,0)
          call send_real(eta(:,j,tau  )  ,imt     ,0,0)
          call send_real(eta(:,j,taum1)  ,imt     ,0,0)
          call send_real(psi(:,j,:)      ,imt*km  ,0,0)
          call send_real(p_full(:,j,:,tau),imt*km ,0,0)
          call send_real(p_full(:,j,:,taum1),imt*km ,0,0)
          if (enable_expl_free_surf) then
           call send_real(bu(:,j,tau  )  ,imt     ,0,0)
           call send_real(bu(:,j,taum1)  ,imt     ,0,0)
           call send_real(bv(:,j,tau  )  ,imt     ,0,0)
           call send_real(bv(:,j,taum1)  ,imt     ,0,0)
          endif
         enddo
        elseif (my_pe==0) then
         call recv_integer(js,1,n,0)
         call recv_integer(je,1,n,0)
         do j=js,je
          call recv_real(u(:,j,:,:,tau  ),imt*km*3,n,0)
          call recv_real(u(:,j,:,:,taum1),imt*km*3,n,0)
          call recv_real(b(:,j,:,tau  )  ,imt*km  ,n,0)
          call recv_real(b(:,j,:,taum1)  ,imt*km  ,n,0)
          call recv_real(eta(:,j,tau  )  ,imt     ,n,0)
          call recv_real(eta(:,j,taum1)  ,imt     ,n,0)
          call recv_real(psi(:,j,:)      ,imt*km  ,n,0)
          call recv_real(p_full(:,j,:,tau),imt*km  ,n,0)
          call recv_real(p_full(:,j,:,taum1),imt*km,n,0)
          if (enable_expl_free_surf) then
           call recv_real(bu(:,j,tau  )  ,imt     ,n,0)
           call recv_real(bu(:,j,taum1)  ,imt     ,n,0)
           call recv_real(bv(:,j,tau  )  ,imt     ,n,0)
           call recv_real(bv(:,j,taum1)  ,imt     ,n,0)
          endif
          write(io) u(:,j,:,:,tau), u(:,j,:,:,taum1)
          write(io) b(:,j,:,tau),   b(:,j,:,taum1)
          write(io) eta(:,j,tau),   eta(:,j,taum1), psi(:,j,:)
          write(io) p_full(:,j,:,tau), p_full(:,j,:,taum1) ! this is new
          if (enable_expl_free_surf) then
           write(io) bu(:,j,tau),  bv(:,j,tau)
           write(io) bu(:,j,taum1),  bv(:,j,taum1)
          endif
         enddo
        endif
        call barrier
      enddo

      if (my_pe==0) close(io)

      if (my_pe==0) then
       call get_free_iounit(io)
       open(io,file='ritt',form='formatted',status='unknown')
       write(io,*) itt
       close(io)
       print*,' bye'
       print*,''
      endif
      end subroutine write_restart







      subroutine show_timing()
c
c--------------------------------------------------------------
c     show timing results here
c--------------------------------------------------------------
c
      use cpflame_module
      implicit none
      real :: fxa,fxb
      integer :: npe,iret
      character (len=128) :: s

      do npe = 0,n_pes
       call barrier

      if (my_pe == npe) then

       print*,''
       print*,'Timing summary for PE #',my_pe 


       print*,' setup time summary       = ',
     &   timing_secs('setup'),' s'
       print*,' cleanup time summary     = ',
     &   timing_secs('cleanup'),' s'
       print*,' costs for measuring      = ',
     &   timing_secs('tictoc'),' s'
       print*,' main loop time summary    = ',
     &   timing_secs('main loop') ,' s'
       fxa= timing_secs('main loop') 

       print*,' '
       print*,' in main loop: '
       print*,' '

       s='buoyancy';fxb=timing_secs(s)
       print*,' ',s(1:len_trim(s)),' time summary     = ',
     &   ' s',fxb, ' = ',fxb/fxa*100,'% of main loop'

       s='momentum';fxb=timing_secs(s)
       print*,' ',s(1:len_trim(s)),' time summary     = ',
     &   ' s',fxb, ' = ',fxb/fxa*100,'% of main loop'

       s='pressure';fxb=timing_secs(s)
       print*,' ',s(1:len_trim(s)),' time summary     = ',
     &   ' s',fxb, ' = ',fxb/fxa*100,'% of main loop'

       s='integrate';fxb=timing_secs(s)
       print*,' ',s(1:len_trim(s)),' time summary     = ',
     &   ' s',fxb, ' = ',fxb/fxa*100,'% of main loop'

       s='diagnose';fxb=timing_secs(s)
       print*,' ',s(1:len_trim(s)),' time summary     = ',
     &   ' s',fxb, ' = ',fxb/fxa*100,'% of main loop'

       endif
       call sub_flush(6)
       call barrier
      enddo
      end subroutine show_timing


