#include "options.inc"



      subroutine setup (domains_,sub_domain_,
     &                      sub_domain_pe0_,my_comm_)
c
c---------------------------------------------------------------------
c     set up everything for the ocean model
c     c.eden
c---------------------------------------------------------------------
c
      use spflame_module
      implicit none
      integer domains_,sub_domain_,sub_domain_pe0_(domains_),my_comm_

      domains=domains_
      sub_domain=sub_domain_
      sub_domain_pe0(1:domains)=sub_domain_pe0_(1:domains)
      my_comm=my_comm_
c
c-----------------------------------------------------------------------
c     set up main module, the grid, the equations of state,
c     the topography, the poisson solver or free surface, 
c     initial conditions and the diagnostic output.
c-----------------------------------------------------------------------
c
      call sub_flush(6)
      call init_spflame_module
      call sub_flush(6)
      call grids
      call init_dens_module(zt,km,my_pe)
      call topog
      call init_congrad
      if (enable_freesurf) call init_freesurf
      call initial_cond
c
c-----------------------------------------------------------------------
c     set up the time manager
c-----------------------------------------------------------------------
c
      time_manager_verbose=.false.               ! be quiet
      if (my_pe==0) time_manager_verbose=.true.  ! be verbose
      call init_time_manager(time_step,runlen,snap_int,
     &                       restart_stamp,itt_restart )
c
c-----------------------------------------------------------------------
c     set up various modules
c-----------------------------------------------------------------------
c
      call init_diag
      call sub_flush(6)
      if (enable_cgh_vert_mixing)      call init_cghmix 
      if (enable_tkemix)               call tkemix_init
      if (enable_bbl)                  call bbl_init
      if (enable_diffusion_isoneutral) call isoneutral_init
      if (enable_diffusion_isopycnic)  call isopycnic_init
      if (enable_blue)                 call blue_init
      if (enable_blue_mean)            call blue_mean_init
      if (enable_quicker_advection)    call init_quicker
      if (enable_upstream_advection)   call init_adv_flux_upstream
      if (enable_4th_advection)        call init_adv_flux_4th
      if (enable_fct_advection)        call init_adv_flux_fct
      if (enable_flux_delimiter)       call flux_delimit_init
                                       call implicitvmix_init()
                                       call passive_tracer_init()
      call sub_flush(6)
c
c-----------------------------------------------------------------------
c     set up the io-procedure for SBC/OBC and sponge layers
c-----------------------------------------------------------------------
c
      if (mother< 0) then
#ifdef couple
        call sbc_couple_init()
#else
       if (.not.enable_simple_sbc) call sbc_initialize
#endif
       if (.not.enable_simple_obc) call obc_initialize
      else
       ! child domains get their data for SBC and OBS from mother domain
      endif

      if (enable_sponge .and. .not. enable_simple_spg) then
        call sponge_initialize
      endif

      call sub_flush(6)

c
c-----------------------------------------------------------------------
c      some output about compiler directives
c-----------------------------------------------------------------------
c
      if (my_pe==0) then
       print*,' '
       print*,' you have switched on the following compiler directives:'
#ifdef partial_cell
       print*,'  -- partial_cell --'
       print*,'      enables partial filled bottom cells'
#endif
#ifdef relax_traditional_approx
       print*,'  -- relax_traditional_approx --'
       print*,'      relax traditional approximation'
#endif
       print*,' end of list'
       print*,' '
       print*,' the following compiler directives were derived:'
#ifdef vector_host
       print*,'  -- vector_host --'
       print*,'      enables vectorized code fragments '
#endif
       print*,' end of list'
       print*,' '
      endif
c
c-----------------------------------------------------------------------
c     some checks
c-----------------------------------------------------------------------
c
      call checks()

      if (my_pe==0) then
       print*,''
       print*,' ok, ready to launch the integration ...'
       print*,''
      endif
      call sub_flush(6)

      end subroutine setup



      subroutine initial_cond()
c
c-----------------------------------------------------------------------
c      set up initial conditions or read restart
c-----------------------------------------------------------------------
c
      use spflame_module
      use freesurf_module
      implicit none
      integer :: i,j,k,n, is,ie,js,je,ioT,ioS,kk
      real :: x
      character(len=32) :: s
      real, allocatable    :: buf(:,:)
      character (len=40) :: tr_name

      is=max(is_pe,2); ie=min(ie_pe,imt-1);
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      u=0.; t=0.; psi=0.

      if (mother<0) then ! mother domain reads first initial conditions

       if (read_restart) then
        call read_restart_file
       else

        restart_stamp ='m/d/y= 1/ 1/1900, h:m:s= 0: 0: 0' 
        itt_restart=1
        u=0.
        if (enable_ktmix) dml(is:ie,js:je)=0.
        if (enable_tkemix) eke=0.

        if (enable_simple_initial_cond) then

         call init_cond_template()

        else 
c
c-----------------------------------------------------------------------
c      read in T and S from files
c      which have to be prepared beforehand. Here we are using
c      the FLAME format for initial conditions.
c-----------------------------------------------------------------------
c
         allocate(buf(imt,km))
         do n=1,nt
          if (my_pe==0) then
           if (n==1) tr_name = 'temp.mom.ic'
           if (n==2) tr_name = 'salt.mom.ic'
           if (n>2) then
            write(tr_name, '("tracer_",i2,".mom.ic")') n
            do i=1,len_trim(tr_name)
              if (tr_name(i:i) == ' ') tr_name(i:i)='0'
            enddo
           endif
           print*,' '
           print*,' reading initial conditions from file ',
     &            tr_name(1:len_trim(tr_name))
           call getunit (ioT,tr_name, 'u s r ieee')
           read(ioT)  
          endif
          do j=1,jmt
           if (my_pe==0) then
            read(ioT)  
            read (ioT) s,x,kk,kk,kk,x,kk,(x,i=1,imt),(x,k=1,km),buf
           endif
           call bcast_real(buf,imt*km,0)
           if (js_pe  <= j .and. je_pe >= j ) then
            do i=is_pe,ie_pe
             t(i,:,j,n,0) = buf(i,:)*tmask(i,:,j)
             t(i,:,j,n,1) = buf(i,:)*tmask(i,:,j)
             t(i,:,j,n,2) = buf(i,:)*tmask(i,:,j)
            enddo
           endif
          enddo ! j loop
          if (my_pe==0) close(ioT)
         enddo ! n loop
         deallocate(buf)

         if (my_pe==0) then
           print*,' done '
           print*,' '
         endif

        endif ! enable_simple_initial_cond

       endif  ! read_restart

      endif  ! mother <0
c
c     exchange initial conditions and whatelse necessary
c     over the sub domains
c
      call initial_exchg
c
c     exchange inner boundaries here
c
      do n=1,nt
c       call border_exchg(t(:,:,:,n,0),km,2)
c       call border_exchg(t(:,:,:,n,1),km,2)
c       call border_exchg(t(:,:,:,n,2),km,2)
       call border_exchg(t(is_pe-2,1,js_pe-2,n,0),km,2)
       call border_exchg(t(is_pe-2,1,js_pe-2,n,1),km,2)
       call border_exchg(t(is_pe-2,1,js_pe-2,n,2),km,2)
      enddo
      do n=1,2
c       call border_exchg(u(:,:,:,n,0),km,2)
c       call border_exchg(u(:,:,:,n,1),km,2)
c       call border_exchg(u(:,:,:,n,2),km,2)
       call border_exchg(u(is_pe-2,1,js_pe-2,n,0),km,2)
       call border_exchg(u(is_pe-2,1,js_pe-2,n,1),km,2)
       call border_exchg(u(is_pe-2,1,js_pe-2,n,2),km,2)
      enddo

      if (enable_freesurf) then
       do n=1,2
        call border_exchg(ubar(:,:,n),1,1)
        call border_exchg(ps(:,:,n),1,1)
       enddo
       do n=1,3
        call border_exchg(etat(:,:,n),1,1)
       enddo
      else
       call border_exchg(ptd(:,:),1,1)
       call border_exchg(guess(:,:),1,1)
       call border_exchg(psi(:,:,1),1,1)
       call border_exchg(psi(:,:,2),1,1)
      endif
      is=max(is_pe,2); ie=min(ie_pe,imt-2) ! inconsistency with clinic !
      js=max(2,js_pe); je=min(je_pe,jmt-2)
      do j=js-1,je+1
       do k=1,km
        call model_dens(t(is-1,k,j,1,tau), t(is-1,k,j,2,tau),
     &                  rho(is-1,k,j),k,ie-is+3
#ifdef partial_cell
     &                       ,ztp(is-1,k,j)
#endif
     &                       )
       enddo
      enddo
      end subroutine initial_cond



      subroutine write_restart_file
c
c-----------------------------------------------------------------------
c      write a restart file
c-----------------------------------------------------------------------
c
      use spflame_module
      use freesurf_module
      implicit none
      integer i,j,k,io,n, is,ie,js,je
      character (len=80) :: filename
      real buf(imt,jmt)

      is=is_pe; ie=ie_pe; js=js_pe; je=je_pe

      write(filename, '("restart_",i3,".dta")') sub_domain
      do i=1,len_trim(filename)
         if (filename(i:i) == ' ') filename(i:i)='0'
      enddo

      if (my_pe==0) then
       print*,' writing to restart file ',filename(1:len_trim(filename))
      endif

      call get_free_iounit(io)
      open(io,file=filename,form='unformatted',status='unknown')

      if (my_pe==0) write(io) itt,current_stamp,nt

      if (enable_freesurf) then

       do n=1,2
        buf(is:ie,js:je)=ubar(is:ie,js:je,n)
        call pe0_recv_2D(buf)
        if (my_pe==0) write(io) buf

        buf(is:ie,js:je)=ps(is:ie,js:je,n)
        call pe0_recv_2D(buf)
        if (my_pe==0) write(io) buf
       enddo

       do n=1,3
        buf(is:ie,js:je)=etat(is:ie,js:je,n)
        call pe0_recv_2D(buf)
        if (my_pe==0) write(io) buf
       enddo

      else

       buf(is:ie,js:je)=ptd(is:ie,js:je)
       call pe0_recv_2D(buf)
       if (my_pe==0) write(io) buf

       buf(is:ie,js:je)=guess(is:ie,js:je)
       call pe0_recv_2D(buf)
       if (my_pe==0) write(io) buf

       buf(is:ie,js:je)=psi(is:ie,js:je,1)
       call pe0_recv_2D(buf)
       if (my_pe==0) write(io) buf

       buf(is:ie,js:je)=psi(is:ie,js:je,2)
       call pe0_recv_2D(buf)
       if (my_pe==0) write(io) buf
      endif

      do k=1,km
       do n=1,nt
        buf(is:ie,js:je)=t(is:ie,k,js:je,n,taum1)
        call pe0_recv_2D(buf)
        if (my_pe==0) write(io) buf
        buf(is:ie,js:je)=t(is:ie,k,js:je,n,tau)
        call pe0_recv_2D(buf)
        if (my_pe==0) write(io) buf
       enddo
       do n=1,2
        buf(is:ie,js:je)=u(is:ie,k,js:je,n,taum1)
        call pe0_recv_2D(buf)
        if (my_pe==0) write(io) buf
        buf(is:ie,js:je)=u(is:ie,k,js:je,n,tau)
        call pe0_recv_2D(buf)
        if (my_pe==0) write(io) buf
       enddo
      enddo

      if (enable_ktmix) then
       buf(is:ie,js:je)=dml(is:ie,js:je)
       call pe0_recv_2D(buf)
       if (my_pe==0) write(io) buf
      endif

      if (enable_tkemix) then
       do k=1,km
        buf(is:ie,js:je)=eke(is:ie,k,js:je,taum1)
        call pe0_recv_2D(buf)
        if (my_pe==0) write(io) buf
        buf(is:ie,js:je)=eke(is:ie,k,js:je,tau)
        call pe0_recv_2D(buf)
        if (my_pe==0) write(io) buf
       enddo
      endif

      close(io)

      if (my_pe==0) then
c
c      write a little file indicating the itt of the last
c      time step. This file will only be written for a successful
c      end of the integration.
c
       write(filename, '("ritt_",i3)') sub_domain
       do i=1,len_trim(filename)
         if (filename(i:i) == ' ') filename(i:i)='0'
       enddo
       call get_free_iounit(io)
       open(io,file=filename,form='formatted',status='unknown')
       write(io,*) itt
       close(io)

       print*,' all done '
      endif

      end subroutine write_restart_file



      subroutine read_restart_file
c
c-----------------------------------------------------------------------
c      read a restart file
c-----------------------------------------------------------------------
c
      use spflame_module
      use freesurf_module
      implicit none
      integer i,j,k,io,n,is,ie,js,je,nt_in
      character (len=80) :: filename
      real buf(imt,jmt)
      logical :: okay=.true.

      is=is_pe; ie=ie_pe; js=js_pe; je=je_pe

      write(filename, '("restart_",i3,".dta")') sub_domain
      do i=1,len_trim(filename)
         if (filename(i:i) == ' ') filename(i:i)='0'
      enddo

      if (my_pe==0) then
       print*,' reading from restart file ',
     &          filename(1:len_trim(filename))
      endif

      call get_free_iounit(io)
      open(io,file=filename,form='unformatted',status='old')

      if (my_pe==0) read(io) itt_restart,restart_stamp,nt_in
      call bcast_integer(itt_restart,1,0)
      call bcast_char(restart_stamp,32,0)

      if (my_pe==0) then 
       print*,' time stamp = ',restart_stamp
       print*,' itt        = ',itt_restart
       print*,' nt         = ',nt_in
       print*,''
       if (nt_in<nt) then
        print*,''
        print*,' WARNING WARNING WARNING WARNING WARNING'
        print*,' found only ',nt_in,' tracer in file'
        print*,' but number of tracer was set to ',nt
        print*,' this means that the rest is set to zero'
        print*,' WARNING WARNING WARNING WARNING WARNING'
        print*,''
       endif
       if (nt_in>nt) then
        print*,''
        print*,' WARNING WARNING WARNING WARNING WARNING'
        print*,' found ',nt_in,' tracer in file'
        print*,' but number of tracer was set to ',nt
        print*,' this means that the rest is not read'
        print*,' WARNING WARNING WARNING WARNING WARNING'
        print*,''
        nt_in = nt      
       endif
      endif

      call bcast_integer(nt_in,1,0)

      if (enable_freesurf) then

       do n=1,2
        if (my_pe==0) read(io) buf
        call pe0_send_2D(buf)
        ubar(is:ie,js:je,n)=buf(is:ie,js:je)

        if (my_pe==0) read(io) buf
        call pe0_send_2D(buf)
        ps(is:ie,js:je,n)=buf(is:ie,js:je)
       enddo

       do n=1,3
        if (my_pe==0) read(io) buf
        call pe0_send_2D(buf)
        etat(is:ie,js:je,n)=buf(is:ie,js:je)
       enddo

      else
       if (my_pe==0) read(io) buf
       call pe0_send_2D(buf)
       ptd(is:ie,js:je)=buf(is:ie,js:je)

       if (my_pe==0) read(io) buf
       call pe0_send_2D(buf)
       guess(is:ie,js:je)=buf(is:ie,js:je)

       if (my_pe==0) read(io) buf
       call pe0_send_2D(buf)
       psi(is:ie,js:je,1)=buf(is:ie,js:je)

       if (my_pe==0) read(io) buf
       call pe0_send_2D(buf)
       psi(is:ie,js:je,2)=buf(is:ie,js:je)
      endif

      do k=1,km
       do n=1,nt_in ! read only nt_in tracers
        if (my_pe==0) read(io) buf
        call pe0_send_2D(buf)
        t(is:ie,k,js:je,n,taum1)=buf(is:ie,js:je)
        if (my_pe==0) read(io) buf
        call pe0_send_2D(buf)
        t(is:ie,k,js:je,n,tau)=buf(is:ie,js:je)
       enddo
c      set remaining tracers to zero
       do n=nt_in+1,nt
        t(is:ie,k,js:je,n,:)=0.
       enddo
c    skip remaining tracers in restart
       do n=nt+1,nt_in 
        if (my_pe==0) read(io) 
       enddo
       do n=1,2
        if (my_pe==0) read(io) buf
        call pe0_send_2D(buf)
        u(is:ie,k,js:je,n,taum1)=buf(is:ie,js:je)
        if (my_pe==0) read(io) buf
        call pe0_send_2D(buf)
        u(is:ie,k,js:je,n,tau)=buf(is:ie,js:je)
       enddo
      enddo

      if (enable_ktmix) then
       if (my_pe==0) read(io,end=2001) buf
       goto 2002
 2001  okay=.false. 
 2002  call bcast_logical(okay,1,0)
       if (.not. okay) goto 2099
       call pe0_send_2D(buf)
       dml(is:ie,js:je)=buf(is:ie,js:je)
      endif

      if (enable_tkemix) then
       do k=1,km
        if (my_pe==0) read(io,end=1001) buf
        goto 1002
 1001   okay=.false. 
 1002   call bcast_logical(okay,1,0)
        if (.not. okay) goto 1099
        call pe0_send_2D(buf)
        eke(is:ie,k,js:je,taum1)=buf(is:ie,js:je)

        if (my_pe==0) read(io,end=1011) buf
        goto 1012
 1011   okay = .false.
 1012   call bcast_logical(okay,1,0)
        if (.not. okay) goto 1099
        call pe0_send_2D(buf)
        eke(is:ie,k,js:je,tau)=buf(is:ie,js:je)
       enddo
      endif

      close(io)
      if (my_pe==0) then
       print*,' done'
      endif
      return

 1099 if (my_pe==0) then
       print*,''
       print*,'---------------------------------------------'
       print*,' WARNING : cannot read TKE from file '
       print*,' ', filename(1:len_trim(filename))
       print*,' setting TKE = 0 '
       print*,'---------------------------------------------'
       print*,''
      endif
      close(io)
      eke=0.
      return

 2099 if (my_pe==0) then
       print*,''
       print*,'---------------------------------------------'
       print*,' WARNING : cannot read mixed layer from file '
       print*,' ', filename(1:len_trim(filename))
       print*,' setting dml = 0 '
       print*,'---------------------------------------------'
       print*,''
      endif
      close(io)
      dml=0.
      return


      end subroutine read_restart_file



