c
c-----------------------------------------------------------------------
c     file contains main driver and restart file IO
c
c-----------------------------------------------------------------------
c
      program barbi
c=======================================================================
c     The main driver for BARBI
c
c                      C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c=======================================================================
      use barbi_module
      implicit none
      integer :: n
c     setup the model
      call setup
c     loop over time steps; use subr. of time manager for time stepping
 10     call time_step_begin()
c       calculate barotropic advection velocities
        call adv_vel
        do n=1,n_order,2
c        time step baroclinic velocity moments
         call clinic(n)
c        time step density moments
         call moments(n)
        enddo
c       time step barotropic mode
        call tropic()
c       write diagnostics
        call diag()
c       shift time stepping indicees, do evetually a mixing time step
        call time_step_end()
        if (.not. end_of_run .or. eulerback2 ) goto 10
c     end of integration, save restart and stop
      if (write_restart) call write_rest
      print*,' end of integration...'
      stop
      end program barbi


      subroutine read_rest(itt_in)
c=======================================================================
c     read NetCDF restart file for BARBI
c
c                      C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c=======================================================================
      use barbi_module
      implicit none
      include "netcdf.inc"
      integer iret,ncid,n,itt_in
      integer ipsiid,iuid,iubid,ivid,ivbid,igid,itid
      integer dims(4), corner(4), edges(4)
      real (kind=4) var(imt,jmt)
      print*,''
      print*,' reading from restart file ',restart_file
      iret=nf_open(restart_file,NF_nowrite,ncid)

      iret= nf_get_att_int(ncid,nf_global,'itt',itt_in)

      iret=nf_inq_varid(ncid,'psi',ipsiid)
      iret=nf_inq_varid(ncid,'u',iuid)
      iret=nf_inq_varid(ncid,'v',ivid)
      iret=nf_inq_varid(ncid,'ub',iubid)
      iret=nf_inq_varid(ncid,'vb',ivbid)
      iret=nf_inq_varid(ncid,'g',igid)

      edges=1
      edges(1)=imt
      edges(2)=jmt
      corner=1
      iret= nf_get_vara_real (ncid,ipsiid,corner,edges,var)
      psi(:,:,2)=var
      iret= nf_get_vara_real (ncid,iuid,corner,edges,var)
      u(:,:,1,0)=var
      iret= nf_get_vara_real (ncid,ivid,corner,edges,var)
      u(:,:,2,0)=var

      corner(3)=2
      iret= nf_get_vara_real (ncid,ipsiid,corner,edges,var)
      psi(:,:,1)=var
      iret= nf_get_vara_real (ncid,iuid,corner,edges,var)
      u(:,:,1,1)=var
      iret= nf_get_vara_real (ncid,ivid,corner,edges,var)
      u(:,:,2,1)=var

      do n=0,n_order
       corner(3)=n+1
       corner(4)=1
       iret= nf_get_vara_real (ncid,igid,corner,edges,var)
       g(:,:,n,0)=var
       iret= nf_get_vara_real (ncid,iubid,corner,edges,var)
       ub(:,:,1,0,n)=var
       iret= nf_get_vara_real (ncid,ivbid,corner,edges,var)
       ub(:,:,2,0,n)=var

       corner(4)=2
       iret= nf_get_vara_real (ncid,igid,corner,edges,var)
       g(:,:,n,1) =var
       iret= nf_get_vara_real (ncid,iubid,corner,edges,var)
       ub(:,:,1,1,n) =var
       iret= nf_get_vara_real (ncid,ivbid,corner,edges,var)
       ub(:,:,2,1,n) =var
      enddo

      call ncclos (ncid, iret)

      print*,'done'
      print*,''
      end subroutine read_rest



      subroutine write_rest()
c=======================================================================
c     write NetCDF restart file for BARBI
c
c                      C.Eden Feb. 2002 (mailto ceden@phys.ocean.dal.ca)
c=======================================================================
      use barbi_module
      implicit none
      include "netcdf.inc"
      integer iret,ncid
      integer ipsiid,iuid,iubid,ivid,ivbid,igid
      integer dims(4), corner(4), edges(4)
      integer lon_dim,lat_dim
      integer order_dim,timedim,n
      real (kind=4) :: var(imt,jmt)
      print*,''
      print*,'  writing to restart file ',restart_file
      ncid = nccre (restart_file, NCCLOB, iret)

      lon_dim  = ncddef(ncid, 'x', imt, iret)
      Lat_dim  = ncddef(ncid, 'y',  jmt, iret)
      order_dim = ncddef(ncid, 'Order',  n_order+1, iret)
      Timedim  = ncddef(ncid, 'Time',2, iret)

      dims(1) = Lon_dim
      dims(2) = Lat_dim
      dims(3) = Timedim
      ipsiid  = ncvdef (ncid,'psi', NCFLOAT,3,dims,iret)
      iuid  = ncvdef (ncid,'u', NCFLOAT,3,dims,iret)
      ivid  = ncvdef (ncid,'v', NCFLOAT,3,dims,iret)
      dims(3) = order_dim
      dims(4) = Timedim
      igid  = ncvdef (ncid,'g', NCFLOAT,4,dims,iret)
      iubid  = ncvdef (ncid,'ub', NCFLOAT,4,dims,iret)
      ivbid  = ncvdef (ncid,'vb', NCFLOAT,4,dims,iret)
      iret= nf_put_att_int(ncid,nf_global,'itt',nf_int,1,itt)

      call ncendf(ncid, iret)
      corner=1
      edges=1
      edges(1)=imt
      edges(2)=jmt
      var= psi(:,:,2)
      iret= nf_put_vara_real (ncid,ipsiid,corner,edges,var)
      var= u(:,:,1,taum1)
      iret= nf_put_vara_real (ncid,iuid,corner,edges,var)
      var= u(:,:,2,taum1)
      iret= nf_put_vara_real (ncid,ivid,corner,edges,var)

      corner(3)=2
      var= psi(:,:,1)
      iret= nf_put_vara_real (ncid,ipsiid,corner,edges,var)
      var= u(:,:,1,tau)
      iret= nf_put_vara_real (ncid,iuid,corner,edges,var)
      var= u(:,:,2,tau)
      iret= nf_put_vara_real (ncid,ivid,corner,edges,var)


      do n=0,n_order
       corner(3)=n+1
       corner(4)=1
       var= g(:,:,n,taum1)
       iret= nf_put_vara_real (ncid,igid,corner,edges,var)
       var= ub(:,:,1,taum1,n)
       iret= nf_put_vara_real (ncid,iubid,corner,edges,var)
       var= ub(:,:,2,taum1,n)
       iret= nf_put_vara_real (ncid,ivbid,corner,edges,var)

       corner(4)=2
       var= g(:,:,n,tau)
       iret= nf_put_vara_real (ncid,igid,corner,edges,var)
       var= ub(:,:,1,tau,n)
       iret= nf_put_vara_real (ncid,iubid,corner,edges,var)
       var= ub(:,:,2,tau,n)
       iret= nf_put_vara_real (ncid,ivbid,corner,edges,var)
      enddo

      call ncclos (ncid, iret)

      print*,'done'
      print*,''
      end subroutine write_rest


