#include "options.inc"

c
c
c
c     changes in the main code:
c
c       zu and ubar in clinic are computed differently
c       in adv_vel call of convU and changed adv_vbt(:,0,:)
c       in spflame call of freesurf instead of tropic
c          (also changes for timing analysis)
c           and call of calc_thickness
c       in setup call of init_freesurf
c       in diag, write eta
c
c       to do : restart IO, diagnostics, domain exchg
c
c
      module freesurf_module
      implicit none

!----------------------------------------------------------------------------
c
c     Compute the barotropic fields:
c     vertically averaged velocity ubar and free surface height etat  
c
c     Two main options
c      
c     A. linear
c     Keeps the ocean volume constant.  Does not conserve tracers. 
c     Similar to the equations solved with Killworth etal
c     1991 and Dukowicz and Smith 1993.
c     Is enabled using not option partial_cell
c
c     B. nonlinear
c     Top ocean cell undulates in time using partial cell methods. 
c     Total tracers constant if no boundary fluxes and dtts=dtuv.
c     Is enabled using option partial_cell
c
c     Basics of the algorithm:
c
c     1. The baroclinic equations have been previously solved for
c     time taup1 = T0+dtuv = itt*dtuv. The vertically integrated
c     forcing terms zu, centered at baroclinic time tau=T0, are known.  
c
c     2. Using zu as forcing terms, held fixed at time tau=T0,
c     the barotropic equations are integrated for N=2*dtuv/dtsf 
c     barotropic timesteps t, each of length dtsf. 
c     The barotropic integration goes from baroclinic 
c     time T0 to T0+2*dtuv (tau to tau+2).
c
c     3. To update the ocean surface height onto the baroclinic
c     time step taup1=T0+dtuv, a leap-frog is used.  However, the 
c     eta(taum1) value is taken from the time averaged eta(taum1),
c     instead of using a Robert time filter.  This approach
c     leads to a smoother and more stable solution.  If the
c     nonlinear free surface is used, then eta(taup1) is used 
c     to update the top cell thicknesses dht(i,1,j) & dhu(i,1,j).
c
c     4. For the nonlinear free surface, the updated full velocity is
c     computed in loadmw.  The corrected updated barotropic velocity
c     is diagnosed from the updated full velocity and eta(taup1).
c     This corrected ubar is used to initialize the next barotropic
c     time step.  
c
c     For spin-ups to equlibrium in coarse models, long tracer 
c     time steps are available, so far as the scheme remains
c     numerically stable.  However, conservation is manifest
c     only with the nonlinear free surface and with dtts=dtuv.  
c
c     authors:
c             R. C. Pacanowski e-mail=> rcp@gfdl.gov    
c             S. M. Griffies   e-mail=> smg@gfdl.gov 
c
c     main variables
c
!     ps    = surface pressure (,,1) is for tau, (,,2) is for tau-1 
!     etat  = surface height on T-cell. 
!             (,,1) is for leap-frog tau and taup1,
!             (,,2) is for time mean tau and tau-1, 
!             (,,3) is for time mean tau and taup1
!     etau  = surface height on U-cells at tau.  diagnosed from etat
!     depthu = time varying depth of a U column; updated on baroclinic times.
!     depthur= 1/depthu
!     thicktau  =top level thickness at tau; from leap-frog etat 
!     emsk  = surface land/sea mask on T-cells
!     umsk  = surface land/sea mask on U-cells
!     ubar  = barotropic velocity defined on U point for "tau" 
!             (,,1) is zonal and (,,2) is meridional velocity
!     convU = convergence of vertically integrated flow (cm/sec)
!     rhosurf = density in surface level.
!----------------------------------------------------------------------------

      real :: volume_tot, salt_tot

      real, allocatable :: umsk(:,:),emsk(:,:)
      real, allocatable :: depthu(:,:),depthur(:,:)
      real, allocatable :: convU(:,:)
      real, allocatable :: rhosurf(:,:)
      real, allocatable :: etau(:,:)
      real, allocatable :: thicktau(:,:)

      real, allocatable :: ubar(:,:,:)
      real, allocatable :: etat(:,:,:)
      real, allocatable :: ps(:,:,:)

      real :: dt_fs 
      real, parameter :: thickmin = 10.0

      real :: obc_etanull = 0.            ! reference water level at open boundaries
c      real :: obc_rest    = 1./(86400.*5) ! restoring time scale for that reference
      real :: obc_rest    = 0.5 ! restoring time scale for that reference

      end module freesurf_module




      subroutine init_freesurf
      use spflame_module
      use freesurf_module
      implicit none
      integer i,j
      real :: dt_min,dt_loc,c_max,c_loc
c
c=======================================================================
c     Initialize some terms for free surface.
c=======================================================================
c
      if (my_pe==0) then
        print*,''
        print*,' Initialising free surface formulation '
        print*,''
      endif

      dt_min =  1e20; c_max  = -1e20
      do j=1,jmt
        do i=1,imt
         if (kmt_big(i,j) >0) then
          c_loc = sqrt(grav*zw(kmt_big(i,j) ))
          dt_loc=min(dxt(i)*cst(j),dyt(j)) / (2*c_loc)
          dt_min = min( dt_loc, dt_min)
          c_max  = max( c_loc , c_max)
         endif
        enddo
      enddo
      dt_fs = time_step / nint( time_step/dt_min ) 
      if (my_pe==0) then
        print*,' maxmimal gravity phase speed ',c_max/100.,' m/s'
        print*,' time step restriction        ',dt_min,' s'
        print*,' time step of barotropic mode ',dt_fs,' s'
        print*,' time step of baroclinic mode ',time_step,' s'
      endif

#ifdef partial_cell
      if (my_pe==0) then
       print*,' Using the nonlinear formulation: '
       print*,'  Top ocean cell undulates in time using partial cells.'
       print*,'  Total tracers should be constant now '
       print*,'  if no there are no boundary fluxes. '
      endif
#else
      if (my_pe==0) then
       print*,' Using the linear formulation: '
       print*,'  The ocean volume should be constant,'
       print*,'  however,  tracers are not conserved'
       print*,'  Similar to the equations solved with Killworth etal'
       print*,'  1991 and Dukowicz and Smith 1993.'
      endif
#endif

      volume_tot = 0.; salt_tot   = 0.

      allocate( umsk    (is_pe-1:ie_pe+1,js_pe-1:je_pe+1) )
      allocate( emsk    (is_pe-1:ie_pe+1,js_pe-1:je_pe+1) )
      allocate( depthu  (is_pe-1:ie_pe+1,js_pe-1:je_pe+1) )
      allocate( depthur (is_pe-1:ie_pe+1,js_pe-1:je_pe+1) )
      allocate( convU   (is_pe-1:ie_pe+1,js_pe-1:je_pe+1) )
      allocate( rhosurf (is_pe-1:ie_pe+1,js_pe-1:je_pe+1) )
      allocate( etau    (is_pe-1:ie_pe+1,js_pe-1:je_pe+1) )
      allocate( thicktau(is_pe-1:ie_pe+1,js_pe-1:je_pe+1) )

      allocate( ubar  (is_pe-1:ie_pe+1,js_pe-1:je_pe+1,2) )
      allocate( etat  (is_pe-1:ie_pe+1,js_pe-1:je_pe+1,3) )
      allocate( ps    (is_pe-1:ie_pe+1,js_pe-1:je_pe+1,2) )

      convU      = 0.
      etau       = 0.
      thicktau   = dzt(1)
      umsk=0.
      emsk=0.
      depthu=0.
      depthur=0.

      do j=js_pe,je_pe
        do i=is_pe,ie_pe
          umsk(i,j)       = min(kmu(i,j),1)
          emsk(i,j)       = min(kmt(i,j),1)
          depthur(i,j)    = hr(i,j)
        enddo
      enddo
      do j=js_pe,je_pe
        do i=is_pe,ie_pe
         if (kmu(i,j)>0) then
#ifdef partial_cell
          depthu(i,j) = hup(i,j)
#else
          depthu(i,j) = zw(kmu(i,j))
#endif
         else
          depthu(i,j) = 0.
         endif
        enddo
      enddo

      call border_exchg(umsk,1,1)
      call set_cyclic  (umask,1,1)
      call border_exchg(emsk,1,1)
      call set_cyclic  (emsk,1,1)
      call border_exchg(depthu,1,1)
      call set_cyclic  (depthu,1,1)
      call border_exchg(depthur,1,1)
      call set_cyclic  (depthur,1,1)
c
c     this is weird and taken from MOM3, no barotropic vel. at obc?
c     do not use that
c
      if (enable_obc_west.and.my_blk_i==1) then  
c       umsk(1,:)   = 0.; umsk(2,:)   = 0.
      endif
      if (enable_obc_east.and.my_blk_i==n_pes_i) then  
c       umsk(imt,:)   = 0.; umsk(imt-1,:) = 0.
      endif
      if (enable_obc_south.and.my_blk_j==1) then  
c       umsk(:,1)   = 0.; umsk(:,2)   = 0.
      endif
      if (enable_obc_north.and.my_blk_j==n_pes_j) then  
c       umsk(:,jmt)     = 0.; umsk(:,jmt-1)   = 0.
      endif

      rhosurf    = rho0*emsk
c
c     the following variables should be read from the restart
c     or set by as initial conditions
c
      ubar(:,:,:)   = 0.
      etat(:,:,:)   = 0.
      ps(:,:,:)     = 0.
c
      if (my_pe==0) then
        print*,''
        print*,' Free surface initialization done '
        print*,''
      endif
      end subroutine init_freesurf




      subroutine calc_convU 
      use spflame_module
      use freesurf_module
      implicit none
      integer :: is,ie,js,je
      integer :: i,j,jm1
      real :: factor
c
c=======================================================================
c     Compute surface vertical velocity centered on bottom of T-cell
c=======================================================================
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      do j=js,je
        jm1=max(1,j-1)
        do i=is,ie
          factor = cstr(j)*dxt2r(i)*dytr(j)
          convU(i,j) = -factor
     &         *(
     &         dyu(j)*(depthu(i,j)*ubar(i,j,1)
     &         -depthu(i-1,j)*ubar(i-1,j,1))
     &         + dyu(jm1)*(depthu(i,jm1)*ubar(i,jm1,1)
     &         -depthu(i-1,jm1)*ubar(i-1,jm1,1))
     &         + dxu(i)*(csu(j)*depthu(i,j)*ubar(i,j,2)
     &         -csu(jm1)*depthu(i,jm1)*ubar(i,jm1,2))
     &         + dxu(i-1)*(csu(j)*depthu(i-1,j)*ubar(i-1,j,2)
     &         -csu(jm1)*depthu(i-1,jm1)*ubar(i-1,jm1,2))
     &         )       
        enddo
      enddo
      call border_exchg(convU,1,1)
      call set_cyclic(  convU,1,1)
      end subroutine calc_convU




      subroutine freesurf
      use spflame_module
      use freesurf_module
      implicit none
c
c=======================================================================
c     time stepping of the explicit free surface
c=======================================================================
c
c     eta    = sea surface elevation in cm at barotropic time t
      real ::  eta   (is_pe-1:ie_pe+1,js_pe-1:je_pe+1)
c     etap   = sea surface elevation in cm at barotropic time t+1
      real ::  etap  (is_pe-1:ie_pe+1,js_pe-1:je_pe+1)
c     hubar  = depth over u-cells * barotropic velocity
      real ::  hubar(is_pe-1:ie_pe+1,js_pe-1:je_pe+1,2)
c
      real :: etarel(is_pe-1:ie_pe+1,js_pe-1:je_pe+1)
      real :: conv_hu(is_pe:ie_pe,js_pe:je_pe),conv_bar
      real :: utempm(is_pe:ie_pe,2)
      real :: gravrho0, gravrho0r  
      real :: fc,fx,fy
      integer :: is,ie,js,je,i,j,k,n,nb,itbt
      real :: de,cp,cmax,etamin,etamax,etashift

      gravrho0   = grav*rho0  
      gravrho0r  = 1./gravrho0

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

      if (eulerback2) then
c
c       just add external mode to u at tau+1
c
        goto 20
      endif

      if (eulerback1 .or. eulerfore) then
c
c       just go ahead
c
      endif

      call border_exchg(zu(:,:,1),1,1)
      call border_exchg(zu(:,:,2),1,1)
      call set_cyclic(zu(:,:,1),1,1)
      call set_cyclic(zu(:,:,2),1,1)
c
c-----------------------------------------------------------------------
c     initialization for barotropic equations
c-----------------------------------------------------------------------
c
c     calculate surface density
c
      do j=js-1,je+1
       call model_dens(t(is-1,1,j,1,tau), t(is-1,1,j,2,tau),
     &                 rhosurf(is-1,j),1,ie-is+3
#ifdef partial_cell
     &                       ,ztp(is-1,1,j)
#endif
     &                       )
       rhosurf(is-1:ie+1,j) = (rhosurf(is-1:ie+1,j) + rho0)
     &                         *tmask(is-1:ie+1,1,j)      
      enddo
c
c     shift (tau) into (tau-1) for next baroclinic time step
c
      ps(:,:,2)     = ps(:,:,1) 
c
c     save leap-frog etat(tau) in thicktau
c     save time mean etat(tau) temporarily in etat(..,1)
c
      thicktau      = max(thickmin,dzt(1)+etat(:,:,1))
      etat(:,:,1)   = etat(:,:,3)
c
c     initialize fields used in the barotropic time stepping 
c
      eta           = etat(:,:,3)
      etap          = 0.
      hubar(:,:,1)  = depthu(:,:)*ubar(:,:,1)
      hubar(:,:,2)  = depthu(:,:)*ubar(:,:,2)
      etarel        = 0.
c
c     initialize ubar for temporary accumulation of hubar
c
      ubar = hubar 
c
c=======================================================================
c     integrate barotropic equations for "nb" barotropic time steps
c     using an Euler forward-backward scheme from baroclinic time T0
c     to baroclinic time T0+2*dtuv.  
c=======================================================================
c
      nb = 2*nint(dt/dt_fs)

      do itbt=1,nb
c
c-----------------------------------------------------------------------
c       time step free surface height using a forward barotropic step
c-----------------------------------------------------------------------
c
        do j=js,je
	  do i=is,ie
	    conv_hu(i,j)  = -dxt2r(i)*dytr(j)*cstr(j)
     &                    *( hubar(i  ,j,1  )*dyu(j)
     &                      +hubar(i,j-1,1  )*dyu(j-1)
     &                    -( hubar(i-1,j,1  )*dyu(j)
     &                      +hubar(i-1,j-1,1)*dyu(j-1))
     &                    +( hubar(i,j,2)*dxu(i) + 
     &                       hubar(i-1,j,2)*dxu(i-1))*csu(j)
     &                    -( hubar(i,j-1,2)*dxu(i) + 
     &                       hubar(i-1,j-1,2)*dxu(i-1))*csu(j-1))
	    etap(i,j) = emsk(i,j) * (eta(i,j) + dt_fs*conv_hu(i,j) )
	  enddo
        enddo
c
c-----------------------------------------------------------------------
c       radiation condition for free surface at open boundaries
c-----------------------------------------------------------------------
c
        if (enable_obc_south.and. my_blk_j==1) then  
         j = 2
         do i=is,ie
           de= eta(i,j+2) - eta(i,j+1)
#ifdef partial_cell
           cmax = -2.*sqrt(grav*htp(i,j) )    
#else
           cmax = -2.*sqrt(grav*zw(max(1,kmt(i,j))))    
#endif
           if (de == 0.) then
             cp = cmax
           else
            cp =-dyu(j+1)/dt_fs*(etap(i,j+1)-eta(i,j+1))/de
            if (cp >= 0.)  cp = 0.
            if (cp < cmax) cp = cmax
           endif
           etap(i,j)=eta(i,j)-cp*dt_fs*dyur(j)*(eta(i,j+1)-eta(i,j))
         enddo
        endif

        if (enable_obc_north.and.my_blk_j==n_pes_j) then
         j = jmt-1
         do i=is,ie
          de = eta(i,j-1) - eta(i,j-2)
#ifdef partial_cell
          cmax = 2.*sqrt(grav*htp(i,j))
#else
          cmax = 2.*sqrt(grav*zw(max(1,kmt(i,j))))
#endif
          if (de == 0.) then
           cp = cmax
          else
           cp = -dyu(j-1)/dt_fs*(etap(i,j-1)-eta(i,j-1))/de
           if (cp <= 0.)  cp = 0.
           if (cp > cmax) cp = cmax
          endif
          etap(i,j)=eta(i,j)-cp*dt_fs*dyur(j)*(eta(i,j)-eta(i,j-1))
         enddo
        endif

        if (enable_obc_west.and.my_blk_i==1) then
         i = 2
         do j=js,je
          de = etap(i+2,j) - etap(i+1,j)
#ifdef partial_cell
          cmax = -2.*sqrt(grav*htp(i,j))
#else
          cmax = -2.*sqrt(grav*zw(max(1,kmt(i,j))))
#endif
          if (de == 0.) then
           cp = cmax
          else
           cp=-dxu(i+1)*csu(j)/dt_fs*(etap(i+1,j)-eta(i+1,j))/de
           if (cp >= 0.)  cp = 0.
           if (cp < cmax) cp = cmax
          endif
          etap(i,j)=eta(i,j) -
     &              cp*dt_fs*dxur(i)*csur(j)*(eta(i+1,j)-eta(i,j))
         enddo
        endif


        if (enable_obc_east.and.my_blk_i==n_pes_i) then
         i = imt-1
         do j=js,je
          de= etap(i-1,j) - etap(i-2,j)
#ifdef partial_cell
          cmax = 2.*sqrt(grav*htp(i,j))
#else
          cmax = 2.*sqrt(grav*zw(max(1,kmt(i,j))))
#endif
          if (de == 0.) then
           cp = cmax
          else
           cp=-dxu(i-1)*csu(j)/dt_fs*(etap(i-1,j)-eta(i-1,j))/de
           if (cp <= 0.)  cp = 0.
           if (cp > cmax) cp = cmax
          endif
          etap(i,j)  = eta(i,j) -
     &      cp*dt_fs*dxur(i)*csur(j)*(eta(i,j)-eta(i-1,j))
         enddo
        endif
c
c-----------------------------------------------------------------------
c       boundary exchange for free surface
c-----------------------------------------------------------------------
c
        call border_exchg(etap,1,1)
        call set_cyclic(  etap,1,1)


#ifdef notdef
c
c        convergence of hubar over entire domain
c
         conv_bar=0.
         do j=js,je
          do i=is,ie
            fx=  dxt(i)*cst(j)*dyt(j)*tmask(i,1,j)
            conv_bar=conv_bar+ conv_hu(i,j)*fx
          enddo
         enddo
         call global_sum(conv_bar)
         conv_bar=conv_bar/tcella(1)
         if (my_pe==0) print*,' conv  = ',conv_bar

         conv_bar=0.
         do j=max(3,js),je
          do i=is,ie
            conv_bar=conv_bar+ (etap(i,j)-eta(i,j) )/dt_fs*emsk(i,j)
     &         *dxt(i)*cst(j)*dyt(j)*tmask(i,1,j)
          enddo
         enddo
         call global_sum(conv_bar)
         conv_bar=conv_bar/tcella(1)
         if (my_pe==0) print*,' eta_t = ',conv_bar
#endif

c
c-----------------------------------------------------------------------
c       diagnose surface pressure and absorb rho0
c-----------------------------------------------------------------------
c
        ps(:,:,1) = etap(:,:)*grav*rhosurf(:,:)/rho0


c#ifdef notdef
c
c-----------------------------------------------------------------------
c       calculate reference water level at open boundaries
c-----------------------------------------------------------------------
c
        if (enable_obc_south.and..not.prescribe_psi_obc_south) then

         if (my_blk_j == 1) then
          j  = 2 
          etamin = minval( eta(is:ie,j+1) )
          etamax = maxval( eta(is:ie,j+1) )
         else
          etamin = 1e20; etamax = -1e12
         endif

         call global_max(etamax)
         call global_min(etamin)

         if ( my_blk_j == 1) then
          etashift = etamin
          if (ABS(etamin) > ABS(etamax)) etashift = etamax
          etarel(:,j) = eta(:,j+1)-(etashift-obc_etanull)
         endif

        endif

        if (enable_obc_south.and.prescribe_psi_obc_south
     &      .and. my_blk_j == 1) then
         etarel(is:ie,2)    = psi_wall_south(is:ie)
        endif


        if (enable_obc_north.and..not.prescribe_psi_obc_north ) then
         if ( my_blk_j == n_pes_j) then
          j  = jmt-1
          etamin = minval( eta(is:ie,j-1) )
          etamax = maxval( eta(is:ie,j-1) )
         else
          etamin = 1e20; etamax = -1e12
         endif

         call global_max(etamax)
         call global_min(etamin)

         if ( my_blk_j == n_pes_j) then
          etashift = etamin
          if (ABS(etamin) > ABS(etamax)) etashift = etamax
          etarel(:,j) = eta(:,j-1)-(etashift-obc_etanull)
         endif
        endif

        if (enable_obc_north.and.prescribe_psi_obc_north
     &      .and. my_blk_j == n_pes_j) then
         etarel(is:ie,jmt-1) = psi_wall_north(is:ie)
        endif

        if (enable_obc_west.and..not.prescribe_psi_obc_west) then
         if ( my_blk_i == 1) then
          i  = 2
          etamin = minval( eta(i+1,js:je) )
          etamax = maxval( eta(i+1,js:je) )
         else
          etamin = 1e20; etamax = -1e12
         endif

         call global_max(etamax)
         call global_min(etamin)

         if (my_blk_i == 1) then
          etashift = etamin
          if (ABS(etamin) > ABS(etamax)) etashift = etamax
          etarel(i,:) = eta(i+1,:)-(etashift-obc_etanull)
         endif
        endif

        if (enable_obc_west.and.prescribe_psi_obc_west
     &      .and. my_blk_i == 1) then
          etarel(2,js:je) = psi_wall_west(js:je)
        endif

        if (enable_obc_east.and..not.prescribe_psi_obc_east) then
         if (my_blk_i == n_pes_i) then
          i  = imt-1
          etamin = minval( eta(i-1,js:je) )
          etamax = maxval( eta(i-1,js:je) )
         else
          etamin = 1e20; etamax = -1e12
         endif

         call global_max(etamax)
         call global_min(etamin)

         if ( my_blk_i == n_pes_i) then
          etashift = etamin
          if (ABS(etamin) > ABS(etamax)) etashift = etamax
          etarel(i,:) = eta(i-1,:)-(etashift-obc_etanull)
         endif
        endif

        if (enable_obc_east.and.prescribe_psi_obc_east
     &      .and. my_blk_i == n_pes_i) then
          etarel(imt-1,js:je) = psi_wall_east(js:je)
        endif

c#endif

c
c-----------------------------------------------------------------------
c       time step hubar using updated ps.
c       this is the "backward" step in the forward-backward scheme.
c-----------------------------------------------------------------------
c
        do j=js,je
	  do i=is,ie
	    fc = cori(i,j,1)*0.5
c
c           surface pressure gradients + forcing + explicit coriolis
c
	    utempm(i,1) = hubar(i,j,1) + dt_fs*(
     &                  - depthu(i,j)*csur(j)*dxu2r(i)
     &                  *( ps(i+1,j+1,1) + ps(i+1,j,1)
     &                   - ps(i  ,j+1,1) - ps(i  ,j,1) )
     &                  + zu(i,j,1)+fc*hubar(i,j,2)
     &                  )
	    utempm(i,2) = hubar(i,j,2) + dt_fs*(
     &                  - depthu(i,j)*dyu2r(j)
     &                  *( ps(i,j+1,1) + ps(i+1,j+1,1)
     &                   - ps(i,j,1) - ps(i+1,j,1) )
     &                  + zu(i,j,2)-fc*hubar(i,j,1)
     &                  )
	  enddo

c#ifdef notdef
c
c         restore to reference water level at open boundaries
c
          if (enable_obc_south .and. j == 2 ) then
	   do i=is,ie
             utempm(i,2) = utempm(i,2) - depthu(i,j)*grav*dyur(j)
     &              *dt_fs*(eta(i,j+1) - etarel(i,j)) * obc_rest
           enddo
          endif

          if (enable_obc_north .and. j == jmt-1) then
	   do i=is,ie
             utempm(i,2) = utempm(i,2) - depthu(i,j)*grav*dyur(j)
     &              *dt_fs*(etarel(i,j) - eta(i,j-1)) * obc_rest
           enddo
          endif
          if (enable_obc_west .and. my_blk_i==1) then
           i=2
           utempm(i,1) = utempm(i,1) - depthu(i,j)*grav*dxur(i)
     &              *csur(j)*dt_fs*(eta(i+1,j) - etarel(i,j)) 
     &              * obc_rest
          endif
          if (enable_obc_east .and. my_blk_i==n_pes_i) then
           i=imt-1
           utempm(i,1) = utempm(i,1)  - depthu(i,j)*grav*dxur(i)
     &              *csur(j)*dt_fs*(etarel(i,j) - eta(i-1,j))
     &              * obc_rest
          endif
c#endif
c
c         solve for hubar using implicit coriolis piece
c
	  do i=is,ie
	    fc = cori(i,j,1)*0.5
	    fx = fc*dt_fs
	    fy = 1./(1.+fx**2)
	    hubar(i,j,1) = umsk(i,j)*(utempm(i,1) + fx*utempm(i,2))*fy
	    hubar(i,j,2) = umsk(i,j)*(utempm(i,2) - fx*utempm(i,1))*fy
	  enddo
        enddo
c
c       boundary exchange
c
        call border_exchg(hubar(:,:,1),1,1)
        call set_cyclic(  hubar(:,:,1),1,1)
        call border_exchg(hubar(:,:,2),1,1)
        call set_cyclic(  hubar(:,:,2),1,1)

        if (enable_obc_south .and. my_blk_j==1) then
         hubar(:,1,:) = hubar(:,2,:)
        endif
        if (my_blk_j == n_pes_j .and. enable_obc_north) then
         hubar(:,jmt-1,:) = hubar(:,jmt-2,:)
         hubar(:,jmt  ,:) = hubar(:,jmt-2,:)
        endif
        if (my_blk_i == 1 .and. enable_obc_west) then
         hubar(1,:,:) = hubar(2,:,:)
        endif
        if (my_blk_i == n_pes_i .and. enable_obc_east) then
         hubar(imt  ,:,:) = hubar(imt-2,:,:)
         hubar(imt-1,:,:) = hubar(imt-2,:,:)
        endif
c
c-----------------------------------------------------------------------
c       accumulate hubar using ubar as a temporary field 
c       accumulate etat in etat(,,3)
c       shift etap = eta(t+1) into eta(t) for next barotropic step
c-----------------------------------------------------------------------
c
        ubar        = ubar        + hubar
        etat(:,:,3) = etat(:,:,3) + etap(:,:)
        eta         = etap           
c
      enddo  !end of barotropic time loop
c
c-----------------------------------------------------------------------
c     Compute time averaged ubar and eta centered at baroclinic 
c     step taup1 = T0+dtuv. 
c-----------------------------------------------------------------------
c    
      ubar(:,:,1) = ubar(:,:,1)*depthur(:,:) / (float(nb)+1.0)
      ubar(:,:,2) = ubar(:,:,2)*depthur(:,:) / (float(nb)+1.0)
      etat(:,:,3) = etat(:,:,3)              / (float(nb)+1.0)
c
c-----------------------------------------------------------------------
c     Update etat(,,1) to baroclinic time T0+dtuv using a leap-frog.  
c     convU, sff, and source_eta are at baroclinic time tau=T0.  
c     Also update etat(,,2) to be the time mean etat(tau).
c     No need to Robert time filter since the etat(taum1) value is
c     a time average.
c-----------------------------------------------------------------------
c
      eta(:,:)    = etat(:,:,1)
      etat(:,:,1) = etat(:,:,2) + c2dt*(convU(:,:))*emsk(:,:)
      etat(:,:,2) = eta(:,:)
      call border_exchg(etat(:,:,1),1,1)
      call set_cyclic(  etat(:,:,1),1,1)
c
c-----------------------------------------------------------------------
c     update surface pressure to new baroclinic time step taup1
c-----------------------------------------------------------------------
c
      ps(:,:,1) = etat(:,:,1)*grav*rhosurf(:,:)
c
c-----------------------------------------------------------------------
c     update internal velocities
c-----------------------------------------------------------------------
c
 20   continue   ! jump from above for eulerback2

      do n=1,2
       do j=js,je
        do k=1,km
         do i=is,ie
          u(i,k,j,n,taup1) = (u(i,k,j,n,taup1) + ubar(i,j,n))
     &                             *umask(i,k,j)
         enddo
        enddo
       enddo
      enddo

      do n=1,2
       call set_cyclic(u(:,:,:,n,taup1),km,2)
       if (my_blk_j == 1 .and. enable_obc_south) 
     &        u(:,:,1,n,taup1) = u(:,:,2,n,taup1)
       if (my_blk_j == n_pes_j .and. enable_obc_north) then
              u(:,:,jmt-1,n,taup1) = u(:,:,jmt-2,n,taup1)
              u(:,:,jmt  ,n,taup1) = u(:,:,jmt-2,n,taup1)
       endif
       if (my_blk_i == 1 .and. enable_obc_west) 
     &        u(1,:,:,n,taup1) = u(2,:,:,n,taup1)
       if (my_blk_i == n_pes_i .and. enable_obc_east) then
              u(imt  ,:,:,n,taup1) = u(imt-2,:,:,n,taup1)
              u(imt-1,:,:,n,taup1) = u(imt-2,:,:,n,taup1)
       endif
       call border_exchg(u(:,:,:,n,taup1),km,2)
      enddo
      end subroutine freesurf


    
      subroutine calc_thickness 
#if defined partial_cell 
      use spflame_module
      use freesurf_module
      use isoneutral_module
      implicit none
      integer :: is,ie,js,je,i,j,k,n
c
c=======================================================================
c     Update top grid cell thickness as well as ocean depth and volume
c     k=1 throughout this routine.  Uses partial cells.
c=======================================================================
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c     thickness of upper t-cell at baroclinic time tau
c
      do j=js,je
       do i=is,ie
        dht(i,1,j) = max(thickmin,dzt(1)+etat(i,j,1))
       enddo
      enddo
      call border_exchg(dht(:,1,:),1,2)
      call set_cyclic(  dht(:,1,:),1,2)
c       
c     thickness of upper u-cell, isoneutral geometric factor,
c     and depth of u-cell ocean.
c
      if (enable_diffusion_isoneutral) then
       do j=js,je
        do i=is,ie
          delqc(i,1,j,0) = 2.*fracdz(1,0)*dht(i,1,j)
          delqc(i,1,j,1) = 2.*fracdz(1,1)*dht(i,1,j)
        enddo
       enddo
       call border_exchg(delqc(:,1,:,0),1,1)
       call border_exchg(delqc(:,1,:,1),1,1)
       call set_cyclic(  delqc(:,1,:,0),1,1)
       call set_cyclic(  delqc(:,1,:,1),1,1)
      endif   

      do j=js,je
       do i=is,ie
        dhu(i,1,j)     = min(dht(i,1,j), dht(i+1,1,j), dht(i,1,j+1),
     &                       dht(i+1,1,j+1))
       enddo
       if (my_blk_i==1)       dhu(1,1,j)   = dht(1,1,j)
       if (my_blk_i==n_pes_i) dhu(imt,1,j) = dht(imt,1,j)
      enddo
      call border_exchg(dhu(:,1,:),1,2)
      call set_cyclic(  dhu(:,1,:),1,2)

      do j=js,je
        do i=is-1,ie+1
          etau(i,j)   = (dhu(i,1,j)-dzt(1))*umsk(i,j)  
          depthu(i,j) = hup(i,j) + etau(i,j)
          if(umsk(i,j) == 1.) then  
            depthur(i,j) = 1./depthu(i,j)
          else
            depthur(i,j) = 0.
          endif
        enddo
      enddo
      call border_exchg(etau,1,1)
      call set_cyclic(  etau,1,1)
      call border_exchg(depthu,1,1)
      call set_cyclic(  depthu,1,1)
      call border_exchg(depthur,1,1)
      call set_cyclic(  depthur,1,1)
c
c     time dependent volume of the ocean for T and U cells
c
c
c     ignore this at the moment
c

c
c     diagnose ubar from full velocity
c
      do n=1,2
       do i=is,ie
        ubar(i,j,n) = 0.
	do k=1,km
         ubar(i,j,n) = ubar(i,j,n)+ u(i,k,j,n,tau)*dhu(i,k,j)
        enddo
        ubar(i,j,n)   = ubar(i,j,n)*depthur(i,j)
       enddo
       call border_exchg(ubar(:,:,n),1,1)
       call set_cyclic(  ubar(:,:,n),1,1)
      enddo

#endif
      end subroutine  calc_thickness




