module driver 
    use const
    use comm
    use var 
    use io
    use util
    use solver
    use physics
    use filter
implicit real*8 (a-h,o-z)

contains

subroutine fsglo
    character fn1*7
    integer tt0,tt1,tt2

    dimension rhomin(k1)
    dimension scr1(i0,j0,k1)

    ! 30 days between climatological data sets
    ! we also use this to restore based on jan 1 climate
    nld=days/30.
    fnew=days/30.-nld
    fold=1.-fnew
    nld=mod(nld,12)+1
    new=nld+1
    if (new.eq.13) new=1

    call ts_nudging_for_watermass

! dog,testing

    qdot2 = 10.
    qdot = 10.

    if (lkppmix) then
      if (mod(itf,5).eq.1) call kppglo
    else 
      if (mod(itf,4).eq.1) call pp82glonew    
    endif 


    call solve_hydrostatic_equation
    ! modified pacanowski and philander vertical mixing coefficients
    ! subracting layer min rho to reduce roundoff error cannot be
    ! done before pp82glo or loop

    !       if (mod(itf,4).eq.1) call pp82glo
    !       if (mod(itf,4).eq.1) call pp82glonew

    ! smagorinsky mixing not useful or needed due to model robustness
    if (mod(itf,4).eq.1) call smagglo

    call mpi_barrier(m_cart,ierr)
    timer(9)=mpi_wtime()
    call  mpi_exch_3d_r4(u2)
    call  mpi_exch_3d_r4(v2)
    call  mpi_exch_3d_r4(t2)
    call  mpi_exch_3d_r4(s2)
    call  mpi_exch_3d_r4(u1)
    call  mpi_exch_3d_r4(v1)
    call  mpi_exch_3d_r4(t1)
    call  mpi_exch_3d_r4(s1)
    call  mpi_exch_3d_r4(ulf)
    call  mpi_exch_3d_r4(vlf)
    call  mpi_exch_3d_r4(tlf)
    call  mpi_exch_3d_r4(slf)
    call  mpi_exch_3d_r4(p)
!yc
    call  mpi_exch_3d_r4(c2)
    call  mpi_exch_3d_r4(c1)
    call  mpi_exch_3d_r4(clf)
!yc
    timer(3)=timer(3)+mpi_wtime()-timer(9)

    lb=1
    lt=2
    if (lfsrf .eq. 0) then
	do j=1,j2
	    do i=1,i2
		uz(i,j,1)=0.
		vz(i,j,1)=0.
		sz(i,j,1)=0.
		tz(i,j,1)=0.
!yc
                cfcz(i,j,1)=0.
!yc
	    enddo
	enddo
    else
	do j=1,j2
	    do i=1,i2
		uz(i,j,1)=w(i+1,j+1,1)*u2(i+1,j+1,1)
		vz(i,j,1)=w(i+1,j+1,1)*v2(i+1,j+1,1)
		sz(i,j,1)=w(i+1,j+1,1)*s2(i+1,j+1,1)
		tz(i,j,1)=w(i+1,j+1,1)*t2(i+1,j+1,1)
!yc
                cfcz(i,j,1)=w(i+1,j+1,1)*c2(i+1,j+1,1)
!yc
	    enddo
	enddo
    endif

    ! =================================================================
    !calculate all internal fluxes and substitute into conservation laws
    ! =================================================================
    ! level-by-level moving window reduces storage.
    ! note: all boundary fluxes are masked out in this loop,
    ! but are included later. this is the main computation loop.
    do k=1,k1
	l=k+1

	! variable diffusivity arrays are used for near-shore laplace filter
	do j=1,j0
	    do i=1,i1
		dhx(i,j)=dmx(i,j,k)
	    enddo
	enddo

	do j=1,j1
	    do i=1,i0
		dhy(i,j)=dmy(i,j,k)
	    enddo
	enddo

	! 4th order pressure gradient
	call compute_pressure_gradient(k)
	
	! vertical fluxes
	call compute_vertical_fluxes(k,lt)
	
	! longitudinal fluxes
	call compute_longitudinal_fluxes(k)	
	
	! latitudinal fluxes
	call compute_latitudinal_fluxes(k)
	
	! ====================
	!conservation equations
	! ====================
	! note: these are only partial updates. pressure, coriolis and boundary
	! forcings are added after statement label 500.
	do j=2,j1
	    odyj=ocs(j)*ody(j)
	    do i=2,i1
		dtin=dt*in(i,j,k)
		! longitudinal momentum
		u2(i,j,k)=u1(i,j,k)-dtin*((ux(i,j-1)-ux(i-1,j-1)+px(i,j))*odx(j)+(uy(i-1,j)-uy(i-1,j-1))*odyj+(uz(i-1,j-1,lt)-uz(i-1,j-1,lb))*odz(k))

		! latitudinal momentum
		v2(i,j,k)=v1(i,j,k)-dtin*((vx(i,j-1)-vx(i-1,j-1))*odx(j)+(vy(i-1,j)-vy(i-1,j-1))*odyj+py(i,j)*ody(j)+(vz(i-1,j-1,lt)-vz(i-1,j-1,lb))*odz(k))

		! salinity
		s2(i,j,k)=s1(i,j,k)-dtin*((sx(i,j-1)-sx(i-1,j-1))*odx(j)+(sy(i-1,j)-sy(i-1,j-1))*odyj+(sz(i-1,j-1,lt)-sz(i-1,j-1,lb))*odz(k))

		! temperature
		t2(i,j,k)=t1(i,j,k)-dtin*((tx(i,j-1)-tx(i-1,j-1))*odx(j)+(ty(i-1,j)-ty(i-1,j-1))*odyj+(tz(i-1,j-1,lt)-tz(i-1,j-1,lb))*odz(k))
!yc
                !tracer
                c2(i,j,k)=c1(i,j,k)-dtin*((cfcx(i,j-1)-cfcx(i-1,j-1))*odx(j)+(cfcy(i-1,j)-cfcy(i-1,j-1))*odyj+(cfcz(i-1,j-1,lt)-cfcz(i-1,j-1,lb))*odz(k))
!yc

	    enddo
	enddo
	i=lb
	lb=lt
	lt=i
    enddo

    ! ===============
    ! surface effect  
    ! ===============
    ! A.wind stress
    call windglo(u2,v2,days,dt,odz)
    
    ! B.heat fluxes, qdot 
    ! (longwave + latent + sensible), positive downward, no ice NOW (melth+SNOW_F+IOFF_F)
    ! call qsurfglo(t2(:,:,1))
    
    ! C.short wave radiation, qdot2 (positive downward)
    ! call qvolumeglo(t2,trans,1)
    
    ! D.evaporation and precipitation (no salt sources)
    ! call salsurfglo(S2(:,:,1))
   
    ! E.Chlorofluorocarbon (CFC) effect
    if (lcfc) call cfcsurf(c2(:,:,1),t2(:,:,1),s2(:,:,1),days,dt,odz)

    ! F.kpp non-local mixing,
    if (lkppmix) call kppsrcglo(t2,s2)

    if (lopen/=0) then 
	! ========================
	! open boundary conditions
	! ========================
	! these are all determined by "known" normal boundary veloctiy (nbv)
	! i.e. boundary normal flux is upwinded for both inflow and outflow.
	! for inflow, exterior "ghost zone" fields are fluxed inward. the ghost
	! zone fields can be set to climatology.
	! loop 506 below does all fluxes on lateral boundaries
	! loop 634 below determines nbv
	! upwind method based on nbv.
	! temporarily use abs() function.
	! later use boundary flag arrays if a vectorized function exists which
	! sets an integer array to the sign bit (i.e., 0 or 1)
	! skip on first time step because advection velocity divergence in
	! boundary zones has not been cleared out yet.

	! nbv update must be after loop 500 (to keep non-divergent bt mode)
	if (mylat .eq. ngy1) then
	    tmp=.5*dt
	    do k=1,k1
		! no open boundaries in periodic longitudinal direction
		do i=2,i1
		    tmpin=in(i,j1,k)*tmp*ody(j1)*csv(j1)*ocs(j1)
		    temp=abs(v(i,j1,k))
		    temp1=tmpin*(temp+v(i,j1,k))
		    temp2=tmpin*(temp-v(i,j1,k))
		    u2(i,j1,k)=u2(i,j1,k)-temp1*ulf(i,j1,k)+temp2*u2(i,j0,k)
		    v2(i,j1,k)=v2(i,j1,k)-temp1*vlf(i,j1,k)+temp2*v2(i,j0,k)
		    s2(i,j1,k)=s2(i,j1,k)-temp1*slf(i,j1,k)+temp2*s2(i,j0,k)
		    t2(i,j1,k)=t2(i,j1,k)-temp1*tlf(i,j1,k)+temp2*t2(i,j0,k)
!yc
                    c2(i,j1,k)=c2(i,j1,k)-temp1*clf(i,j1,k)+temp2*c2(i,j0,k)
!yc
                enddo
	    enddo

	endif

	!do j=2,j1
	!    do i=2,i1
	!        ! take out momentum with positive e-p, and add none for negative e-p
	!        ! this is for convenience only (minor effect, but not really realistic
	!        ! because this is based only on the net rather than actual fluxes)
	!        tmp=.5*(w(i,j,1)-abs(w(i,j,1)))*odz(1)
	!        u2(i,j,1)=u2(i,j,1)+dt*tmp*u2(i,j,1)
	!        v2(i,j,1)=v2(i,j,1)+dt*tmp*v2(i,j,1)
	!        t2(i,j,1)=t2(i,j,1)+dt*w(i,j,1)*t2(i,j,1)*odz(1)
	!    enddo
	!enddo
    endif
    
    do j=2,j1
	do i=2,i1
	    tmp1=dble(t_nudge(i,j))
	    tmp2=1.-tmp1
	    ! skli is no longer used (no salinity restoring)
	    ! tkli is climatological value interpolated to present time
	    tkli=fold*tclis(i,j,1,nld)+fnew*tclis(i,j,1,new)
	    ! qavg(watts/cm-cm)= (watts/deg/cc)*heat*dx(j)*dy*dz/(dx(j)*dy)/time
	    !    where qavg=heat/time in deg ! per time step
	    !======================================================================
	    ! add model-consistent climatology determined surface fluxes, which
	    ! are proporational to time averaged restoring rate for present month
	    ! add to all 3 time levels to get full and proper effect with fltw
	    ! method; otherwise "final" effect on t1 is only partial!
	    temp=in(i,j,1)*qavg(i-1,j-1,nld)
	    t2(i,j,1)=t2(i,j,1)+temp
	    t1(i,j,1)=t1(i,j,1)+temp
	    tlf(i,j,1)=tlf(i,j,1)+temp
!yc
!            c2(i,j,1)=c2(i,j,1)+temp
!            c1(i,j,1)=c1(i,j,1)+temp
!            clf(i,j,1)=clf(i,j,1)+temp
!yc
	    !======================================================================
	    ! surface restoring: nudge toward surface t climatology
	    ! to emulate eddy damping by atmospheri! heat flux.
	    temp=t2(i,j,1)
	    !      t2(i,j,1)=in(i,j,1)*(tmp*temp+tau*tkli)+(1-in(i,j,1))*tkli
	    t2(i,j,1)=in(i,j,1)*(tmp2*temp+tmp1*tkli)+(1-in(i,j,1))*tkli
!yc
!            c2(i,j,1)=in(i,j,1)*(tmp2*temp+tmp1*tkli)+(1-in(i,j,1))*tkli
!yc

	    ! accumulate nudgings for present month
	    tnudge(i-1,j-1)=tnudge(i-1,j-1)+t2(i,j,1)-temp
	enddo
    enddo

    ! ===========================================================
    ! new strategy for surface fluxes
    ! determine ensemble (climatological) monthly mean model-
    ! consistent surface heat and freshwater volume sources
    ! (update long term mean sources according to present month
    ! anomalous changes of surface temperature and salinity
    ! where anomalous changes are differences from observed
    ! climatological changes).
    ! ===========================================================

    ! procedure:
    ! a) calculate surface climate t anomaly at end of month
    ! b) add anomaly to tnudge to get qavg for present month
    ! c) blend with results from previous years
    i=30*itfday
    if (mod(itf,i).ne.i-1) go to 575
    ! nsombo=1 when qavg is initialized,
    ! and we include initial value in averaging process

    ! ww nsombo(nld)=nsombo(nld)+1
    if (itf.gt.i*12*10) nsombo(nld)=nsombo(nld)+1

    ens1=1./nsombo(nld)
    ens2=1.-ens1

    ! plot nudgings for present month
    !     do 535 j=2,j1
    !     do 535 i=2,i1
    !     call xyglo('r','present month net t nudging   ',scr(1,1,2),1,1)

    ! for long term average, temp=1/years/(30*itfday)
    !   where years=number of years included in average

    temp=1./(30*itfday)
    tmp1=temp*ens1
    tmp2=z(3)*odt*ens1

    qsum=0.
    do j=1,j2
	tmp=dx(j+1)*dy(j+1)
	do i=1,i2
	    ! for time consistencfcy with evap we also use unadjusted qavg here
	    qsum=qsum+in(i+1,j+1,1)*qavg(i,j,new)*tmp
	    ! evap is net surface outflow (positive for net positive e-p)
	    qavg(i,j,nld)=qavg(i,j,nld)+tmp1*(tnudge(i,j)+tclis(i+1,j+1,1,new)-t2(i+1,j+1,1))
	    scr(i+1,j+1,2)=daodt*30.*qavg(i,j,nld)
	    ! set monthly sum of nudgings to zero to initialize next month
	    tnudge(i,j)=0.
	enddo
    enddo
    
    if (days<=361.)then 
	! ======================================================================
	! first year only (gives reasonable advection estimate for next month;
	! also, avoids accumulation of advection-induced drift
	do j=1,j2
	    do i=1,i2
		qavg(i,j,new)=qavg(i,j,nld)
	    enddo
	enddo
    endif
    
    !convert qsum to mean watts/m-m and evap to m**3/sec
    tp=qsum
    call mpi_barrier(m_cart,ierr)
    timer(9)=mpi_wtime()
    call mpi_allreduce(tp,qsum,1,mpi_real8,mpi_sum,m_cart,ierr)
    timer(3)=timer(3)+mpi_wtime()-timer(9)
    qsum=watts*qsum/a(1)


    ! =============
    ! bottom stress
    ! =============
    575  tmp=dt*drag
    
    do j=2,j1
	do i=2,i1
	    k=kb(i,j)
	    ! don't exceed explicit limit
	    drg=min(.2,in(i,j,k)*tmp*odz(k)*sqrt(u1(i,j,k)**2+v1(i,j,k)**2))
	    u2(i,j,k)=u2(i,j,k)-drg*u1(i,j,k)
	    v2(i,j,k)=v2(i,j,k)-drg*v1(i,j,k)
	enddo
    enddo
    
    ! ====================
    ! trapezoidal coriolis
    ! ====================
    do k=1,k1
	do j=2,j1
	    do i=2,i1
		tmp=.5*(f(j-1)+ulf(i,j,k)*tanphi(j-1))*dt
		temp=1./(1.+tmp**2)
		qu=u2(i,j,k)+tmp*v1(i,j,k)
		qv=v2(i,j,k)-tmp*u1(i,j,k)
		u2(i,j,k)=temp*(qu+tmp*qv)
		v2(i,j,k)=temp*(qv-tmp*qu)
	    enddo
	enddo
    enddo

    call mpi_barrier(m_cart,ierr)
    timer(9)=mpi_wtime()

    call mpi_exch_3d_r4(u2)
    call mpi_exch_3d_r4(v2)

    timer(3)=timer(3)+mpi_wtime()-timer(9)
    ! =====================================================
    ! interpolate "a" grid quantities to "c" grid locations
    ! =====================================================
    ! 4th-order interpolations
    ! unnecessary masking has been removed here
    do k=1,k1
	do j=2,j1
	    do i=1,i1
		scr(i,j,k)=6.*(u2(i,j,k)+u2(i+1,j,k))
	    enddo
	    ! skip loop 612 for 2nd order scheme
	    do i=2,i2
		scr(i,j,k)=scr(i,j,k)+iu(i,j,k)*(-u2(i-1,j,k)+u2(i,j,k)+u2(i+1,j,k)-u2(i+2,j,k))
	    enddo
	enddo
    enddo

    call mpi_sendrecv(u2(i2,1,1),1,m_vlon3d,m_e,1,scr1(i0,1,1),1,m_vlon3d,m_w,1,m_cart,istat,ierr)
    call mpi_sendrecv(u2(3,1,1),1,m_vlon3d,m_w,1,scr1(i1,1,1),1,m_vlon3d,m_e,1,m_cart,istat,ierr)

    do k=1,k1
	do j=2,j1
	    scr(1,j,k)=scr(1,j,k)+iu(1,j,k)* (-scr1(i0,j,k)+u2(1,j,k)+u2(2,j,k)-u2(3,j,k))
	    scr(i1,j,k)=scr(i1,j,k)+iu(i1,j,k)*(-u2(i2,j,k)+u2(i1,j,k)+u2(i0,j,k)-scr1(i1,j,k))
	    do i=1,i1
		u(i,j,k)=o12*scr(i,j,k)*iu(i,j,k)
	    enddo
	enddo
    enddo


    do k=1,k1
	do j=js,jf
	    do i=2,i1
		scr(i,j,k)=6.*(v2(i,j,k)+v2(i,j+1,k))
	    enddo
	enddo
	do j=2,j2
	    do i=2,i1
		scr(i,j,k)=scr(i,j,k)-v2(i,j-1,k)+v2(i,j,k)+v2(i,j+1,k)-v2(i,j+2,k)
	    enddo
	enddo
    enddo

    call mpi_sendrecv(v2(1,j2,1),1,m_vlat3d,m_n,1,scr1(1,j0,1),1,m_vlat3d,m_s,1,m_cart,istat,ierr)
    call mpi_sendrecv(v2(1,3,1),1,m_vlat3d,m_s,1,scr1(1,j1,1),1,m_vlat3d,m_n,1,m_cart,istat,ierr)

    do k=1,k1
	! skip loop for 2nd order scheme
	if (mylat .ne. 0) then
	    do i=2,i1
		scr(i,1,k)=scr(i,1,k)-scr1(i,j0,k)+v2(i,1,k)+v2(i,2,k)-v2(i,3,k)
	    enddo
	endif

	if (mylat .ne. ngy1) then
	    do i=2,i1
		scr(i,j1,k)=scr(i,j1,k)-v2(i,j2,k)+v2(i,j1,k)+v2(i,j0,k)-scr1(i,j1,k)
	    enddo
	endif

	do j=js,jf
	    do i=2,i1
		v(i,j,k)=o12*scr(i,j,k)*iv(i,j,k)
	    enddo
	enddo
    enddo


    ! the following non-porous lid approach may be replaced by a good estimate
    ! of e-p, as in medina and north pacifi! codes, as per dietrich, et al.
    ! e-p may be significant when there is significant lateral water vapor
    ! transport and fallout by translating hurricanes and typhoons.

    ! non-porous lid approach:
    ! maintain top layer salinity by exchange with layer 2, and ignore
    ! possible vortex stretching implied by non-zero e-p (w(i,j,1)). this
    ! simplifies global incompressibility treatment in domain decomposed mpi
    ! code.
    ! do 627 j=2,j1
    ! do 627 i=2,i1
    ! 627 w(i,j,1)=0.
    ! maintain surface climatological s in the long run
    ! 90 day nudging time scale toward surface climatological s
    !conserves total salt material assuming all wet points > one layer deep

    !        tmp1=1./(90.*daodt)
    !        tmp2=tmp1*odz(2)/odz(1)
    !        do 628 j=2,j1
    !        do 628 i=2,i1
    !        tmp=in(i,j,1)*(sclis(i,j,1,new)-s1(i,j,1))
    !        s1(i,j,1)=s1(i,j,1)+tmp1*tmp
    !   628  s1(i,j,2)=s1(i,j,2)-tmp2*tmp

    ! ====================================================
    ! create non-divergent u,v,w for next step by
    ! "clearing out the divergence of the barotropic mode"
    ! ====================================================
    ! use scr and scr to store changes for incompressibility
    do k=1,4
	do j=1,j0
	    do i=1,i0
		scr(i,j,k)=0.
	    enddo
	enddo
    enddo
    ! ---------------------------------------------------
    ! iterate evp solver to get zero normal flow at shore
    ! and exactly non-divergent bt mode
    ! ---------------------------------------------------
    ! mxmask is maximum number of iterations allowed
    !check whether increasing from 8 to 9 eliminates noisy regions
    mxmask=9
    ! note: it would be slightly more efficient to separate out the bt mode
    ! (as per dan wright's suggestion) if more than one iteration is
    ! performed, but this has not yet been implemented with river sources
    ! and precipitation

    call mpi_barrier(m_cart,ierr)
    timer(9)=mpi_wtime()
    itmask=0
    652  itmask=itmask+1

    ! eventually, we should use wright's bt model approach
    ! for efficiencfcy, but the bt mode's non-zero divergence
    ! requires special treatment in that case.


    ! zero out advection velocity over land

    tmp=0.
    do j=2,j1
	do i=1,i1
	    u(i,j,1)=iu0(i,j)*u(i,j,1)
	enddo
    enddo

    do j=js,jf
	do i=2,i1
	    v(i,j,1)=iv0(i,j)*v(i,j,1)
	enddo
    enddo

    !      call write_rerun
    !      call mpi_finalize(ierr)
    !      stop 654

    if (lfsrf .eq. 0) then
	do k=1,k1
	    temp1=1./odz(k)
	    do j=2,j1
		temp=ocs(j)*ody(j)
		do i=2,i1
		    w(i,j,k+1)=w(i,j,k)-((u(i,j,k)-u(i-1,j,k))*odx(j)+(csv(j)*v(i,j,k)-csv(j-1)*v(i,j-1,k))*temp)*temp1
		end do
	    end do
	end do
	do j=2,j1
	    jj=j-1      !(j:cell, jj:size of s, easy to read)
	    do i=2,i1
		ii=i-1    !(i:cell, ii:size of s, easy to read)
		s(ii,jj)=w(i,j,kb(i,j)+1)
	    end do
	end do
    elseif (lfsrf .eq. 1) then
	do k=k1,1,-1
	    temp1=1./odz(k)
	    do j=2,j1
		temp=ocs(j)*ody(j)
		do i=2,i1
		    w(i,j,k)=w(i,j,k+1)+((u(i,j,k)-u(i-1,j,k))*odx(j)+(csv(j)*v(i,j,k)-csv(j-1)*v(i,j-1,k))*temp)*temp1
		end do
	    end do
	end do
	do j=2,j1
	    jj=j-1      !(j:cell, jj:size of s, easy to read)
	    do i=2,i1
		ii=i-1;    !(i:cell, ii:size of s, easy to read)
		s(ii,jj)=-w(i,j,1)
	    end do
	end do
    else
	if(myid.eq.0)write(*,*)'lfsrf is not set correctly'
	call mpi_finalize(ierr)
	stop
    endif


    ! ================================
    ! the rebirth of bir
    ! something for m.s. to chew on!!!
    ! ================================
    ! iter is the bir sweep number
    timer(7)=mpi_wtime()
	
    if (lsolver .eq. 1) then
	call rep(al,ab,ac,ar,at,rinv,rinv1,dum0,dum1,dum2,dumg,s,x,ie)
	!       call rep2(al,ab,ac,ar,at,rinv,rinv1,dum0,dum1,dum2,dumg,s,x,ie)
	elseif(lsolver .eq. 2) then
	call p_bicgstab_le(ab,al,ac,ar,at,s,x,cb,cl,cc,cr,ct,cgr,cgrh,cgp,cgv,cgs,cgt,cgph,cgsh)
    endif

    !    call rep(al,ab,ac,ar,at,rinv,rinv1,dum0,dum1,dum2,dumg,s,x,ie)

    timer(6)=mpi_wtime()-timer(7)
    ! rigid-lid pressure adjustment
    do j=1,j0
	do i=1,i0
	    p0(i,j)=p0(i,j)+odt*x(i,j)
	enddo
    enddo
    ! bt mode adjustment (from pressure gradient change)
    do j=2,j1
	do i=1,i1
	    scr(i,j,1)=-(x(i+1,j)-x(i,j))*odx(j)
	enddo
    enddo
    do j=js,jf
	do i=2,i1
	    scr(i,j,2)=-(x(i,j+1)-x(i,j))*odyv(j)
	enddo
    enddo

    ! top layer
    do j=2,j1
	do i=1,i1
	    scr(i,j,3)=scr(i,j,3)+scr(i,j,1)
	    u(i,j,1)=iu0(i,j)*u(i,j,1)+scr(i,j,1)
	enddo
    enddo
    do j=js,jf
	do i=2,i1
	    scr(i,j,4)=scr(i,j,4)+scr(i,j,2)
	    v(i,j,1)=iv0(i,j)*v(i,j,1)+scr(i,j,2)
	enddo
    enddo
    ! remaining layers
    do k=2,k1
	do j=2,j1
	    do i=1,i1
		u(i,j,k)=u(i,j,k)+scr(i,j,1)*iu(i,j,k)
	    enddo
	enddo
	do j=js,jf
	    do i=2,i1
		v(i,j,k)=v(i,j,k)+scr(i,j,2)*iv(i,j,k)
	    enddo
	enddo
    enddo
    !check convergence to zero flow over land
    tmp=0.
    do j=2,j1
	do i=2,i2
	    tmp=max(tmp,(1-iu0(i,j))*abs(u(i,j,1)))
	enddo
    enddo
    do j=2,j2
	do i=2,i1
	    tmp=max(tmp,(1-iv0(i,j))*abs(v(i,j,1)))
	enddo
    enddo
    tp=tmp
    call mpi_allreduce(tp,tmp,1,mpi_real8,mpi_max,m_cart,ierr)
    if(itf-it0.lt.itfday.or.mod(itf,itfday).eq.0) write(14,683) itf-it0,itmask,tmp,timer(6)
    !if (myid .eq. 0 .and. mod(itf,itfday).eq.0) write(*,683) itf-it0,itmask,tmp,timer(6)
    if (myid .eq. 0) write(*,683) itf-it0,itmask,tmp,timer(6)
    683  format('@itf-it0===',i8,',itmask=',i3,', vmx on land=',f8.4,' cm/s, solver time =', f8.4,'sec')

    if (tmp.gt..0002 .and. itmask.lt.mxmask) go to 652

    timer(4)=timer(4)+mpi_wtime()-timer(9)
    ! this completes the advanced time level advection velocity
    ! having exactly zero barotropic divergence and satisfying all
    ! kinematic and specified open boundary conditions
    if (lfsrf .eq. 0) then
	do k=1,k1                !(k:layer& top   face for w, k+1:bottom face for w)
	    tmp=1./odz(k)
	    do j=2,j1              !(j:cell & north face for v)
		temp=ocs(j)*ody(j)
		do i=2,i1            !(i:cell & east  face for u)
		    w(i,j,k+1)=iw(i,j,k+1)*(w(i,j,k)-((u(i,j,k)-u(i-1,j,k))*odx(j)+(csv(j)*v(i,j,k)-csv(j-1)*v(i,j-1,k))*temp)*tmp)
		end do
	    end do
	end do
    elseif (lfsrf .eq. 1) then
	do j=2,j1              !(j:cell & north face for v)
	    temp=ocs(j)*ody(j)
	    do i=2,i1            !(i:cell & east  face for u)
		do k=kb(i,j),1,-1
		    tmp=1./odz(k)
		    w(i,j,k)=(w(i,j,k+1)+((u(i,j,k)-u(i-1,j,k))*odx(j)+(csv(j)*v(i,j,k)-csv(j-1)*v(i,j-1,k))*temp)*tmp)
		end do
	    end do
	end do
    endif

    do k=1,k1
	do j=2,j1
	    do i=2,i1
		u2(i,j,k)=12.*(u(i-1,j,k)+u(i,j,k))
		v2(i,j,k)=12.*(v(i,j-1,k)+v(i,j,k))
	    enddo
	enddo
    enddo

    call mpi_sendrecv(u(i2,1,1),1,m_vlon3d,m_e,1,scr(i2,1,1),1,m_vlon3d,m_w,1,m_cart,istat,ierr)
    call mpi_sendrecv(u(2,1,1),1,m_vlon3d,m_w,1,scr(2,1,1),1,m_vlon3d,m_e,1,m_cart,istat,ierr)
    do k=1,k1
	do j=2,j1
	    u2(2,j,k)=o24*in(2,j,k)*(u2(2,j,k)-scr(i2,j,k)+u(1,j,k)+u(2,j,k)-u(3,j,k))
	    u2(i1,j,k)=o24*in(i1,j,k)*(u2(i1,j,k)-u(i0-3,j,k)+u(i2,j,k)+u(i1,j,k)-scr(2,j,k))
	enddo
    enddo
    do k=1,k1
	do j=2,j1
	    do i=3,i2
		u2(i,j,k)=o24*in(i,j,k)*(u2(i,j,k)-u(i-2,j,k)+u(i-1,j,k)+u(i,j,k)-u(i+1,j,k))
	    enddo
	enddo
    enddo
    call mpi_sendrecv(v(1,j2,1),1,m_vlat3d,m_n,1,scr(1,j2,1),1,m_vlat3d,m_s,1,m_cart,istat,ierr)
    call mpi_sendrecv(v(1,2,1),1,m_vlat3d,m_s,1,scr(1,2,1),1,m_vlat3d,m_n,1,m_cart,istat,ierr)
    do k=1,k1
	do i=2,i1
	    do j=3,j2
		v2(i,j,k)=v2(i,j,k)-v(i,j-2,k)+v(i,j-1,k)+v(i,j,k)-v(i,j+1,k)
	    enddo
	    if(mylat.ne.0) v2(i,2,k)=v2(i,2,k)-scr(i,j2,k)+v(i,1,k)+v(i,2,k)-v(i,3,k)
	    if(mylat.ne.ngy1) v2(i,j1,k)=v2(i,j1,k)-v(i,j0-3,k)+v(i,j2,k)+v(i,j1,k)-scr(i,2,k)
	enddo
    enddo
    do k=1,k1
	do j=js,j1
	    do i=2,i1
		v2(i,j,k)=in(i,j,k)*o24*v2(i,j,k)
	    enddo
	enddo
    enddo


    if (itf-it0.gt.itfday.and.mod(itf,itfday).ne.0) go to 712
	!check incompressibility
	tmp=0.
	err=0.
	do k=1,k1
	    do j=2,j1
		do i=2,i1
		    tmp1=(u(i,j,k)-u(i-1,j,k))*odx(j)
		    tmp2=(csv(j)*v(i,j,k)-csv(j-1)*v(i,j-1,k))*ocs(j)*ody(j)
		    tmp3=(w(i,j,k+1)-w(i,j,k))*odz(k)
		    tmp=tmp+max(abs(tmp1),abs(tmp2),abs(tmp3))*in(i,j,k)
		    err=err+abs(tmp1+tmp2+tmp3)*in(i,j,k)
		enddo
	    enddo
	enddo

	tep=err
	call mpi_allreduce(tep,err,1,mpi_real8,mpi_sum,m_cart,ierr)
	tep=tmp
	call mpi_allreduce(tep,tmp,1,mpi_real8,mpi_sum,m_cart,ierr)
	err=err/tmp
	write(14,711) err
	call flush(14)
	711  format(' *** normalized mean incompressibility error = ',1pe9.2)
    712  psm=0.

    do j=2,j1
	do i=2,i1
	    psm=psm+p0(i,j)
	enddo
    enddo
    tmp=psm
    call mpi_allreduce(tmp,psm,1,mpi_real8,mpi_sum,m_cart,ierr)

    psm=psm/dble(i2t*j2t)

    do j=2,j1
	do i=2,i1
	    p0(i,j)=p0(i,j)-psm
	enddo
    enddo


    call mpi_barrier(m_cart,ierr)  
    call mpi_exch_3d_r4(u2)
    call mpi_exch_3d_r4(v2)
    call mpi_exch_3d_r4(t2)
    call mpi_exch_3d_r4(s2)
!yc
    call mpi_exch_3d_r4(c2)
!yc

    if (.false.) then
	tt0=it0/daodt
	n=tt0/360.+.0001
	tmp=tt0/360.
	tt2=(tmp-n)*360.+.1
	tt1=tt0/360.+.0001+1
	write(fn1,'(i3.3,a1,i3.3)'),tt1,'_',tt2

	do i = 1,i0,2
	do j = 1,j0,2

	if (ydeg(j).ge.30.1.and.ydeg(j).le.40.1 .and.xdeg(i).ge.(360.-80.).and.xdeg(i).le.(360.-72.)) then
	    open(1005,file='ts/field1_'//fn1//'_'//grd,form='unformatted')
	    do k=1,10
	    write(1005) ulf(i,j,k),vlf(i,j,k),p(i,j,k),tlf(i,j,k)
	    enddo
	    if(itf.eq.it0) then
		print*,'field i,j,myid',myid,grd,i,j,ydeg(j),xdeg(i),ulf(i,j,4)
	    endif
	endif

	if (ydeg(j).ge.36.1.and.ydeg(j).le.42.1 .and.xdeg(i).ge.(360.-70.).and.xdeg(i).le.(360.-60.)) then
	    open(1006,file='ts/field2_'//fn1//'_'//grd,form='unformatted')
	    do k=1,10
	    write(1006) ulf(i,j,k),vlf(i,j,k),p(i,j,k),tlf(i,j,k)
	    enddo
	    if(itf.eq.it0) then
		print*,'field i,j,myid',myid,grd,i,j,ydeg(j),xdeg(i),ulf(i,j,4)
	    endif
	endif


	enddo
	enddo
    endif

    !   rawf=0.53
    !   fltw=0.1
    !   modified filter scheme
    call modified_filter
   
    ! ----------------------------------------------
    ! relax swamp layer and masked water to interior
    ! to reduce seepage effect  
    ! ----------------------------------------------
    call mpi_barrier(m_cart,ierr)
    call mpi_exch_2d_r4(t1(1,1,1))
    call mpi_exch_2d_r4(s1(1,1,1))
!yc
    call mpi_exch_2d_r4(c1(1,1,1))
!yc


    do j=2,j1
	do i=2,i1
	    tmp=.25*(1-in(i,j,1))
	    scr(i,j,1)=in(i,j,1)*t1(i,j,1)+tmp*(t1(i-1,j,1)+t1(i+1,j,1)+t1(i,j+1,1)+t1(i,j-1,1))
	    scr(i,j,2)=in(i,j,1)*s1(i,j,1)+tmp*(s1(i-1,j,1)+s1(i+1,j,1)+s1(i,j+1,1)+s1(i,j-1,1))
!yc
            scr(i,j,13)=in(i,j,1)*c1(i,j,1)+tmp*(c1(i-1,j,1)+c1(i+1,j,1)+c1(i,j+1,1)+c1(i,j-1,1))
!yc
	enddo
    enddo

    do j=2,j1
	do i=2,i1
	    t1(i,j,1)=scr(i,j,1)
	    s1(i,j,1)=scr(i,j,2)      
!yc
            c1(i,j,1)=scr(i,j,13)
!yc

	enddo
    enddo
end subroutine


subroutine set_windmix
!    open(8,file='prepglo/data/windmix'//grd,form='unformatted')
!    read(8) vbk,hbk
!    close(8)
     call rdnc_windmix
end subroutine


subroutine set_rundata 
!    open(9,file='prepglo/data/rundata'//grd,form='unformatted')
!    read(9) case,dscrib,dt,tlz,b,g,dm0,de0,dmz0,f,tanphi,rzmx, &
!    drag,fltw,daodt,odt,prn,orzmx,ofltw,y,yv,yvdeg,ydeg,cs,csv, &
!    ocs,ocsv,dx,dxv,odx,odxv,dy,dyv,ody,odyv,ktrm,lrstrt,mxit,isav, &
!    mxsav,lopen,lmovi,mvi,lfsrf,lsolver,xdeg,lprecon,rawf
!    close(9)
    call rdnc_rundata
   
    itfday=daodt+.0001
    !conversion of qavg (deg c per model time step) to watts/m-m
    ! use factors for deg ! per year and g/kg/yr
    watts=360.*24.*3600.*0.00134/(dt*odz(1))

    ! global salinity budget weighting factor             
end subroutine


subroutine set_zkb 
    ! --------------------------------------------------------------------
    ! metrics, z-levels(z) and bathymetry (logical depth array (kb)) data
    ! --------------------------------------------------------------------
!    open(60,file='prepglo/data/zkb'//grd,form='unformatted')
!    read(60) tlz,z,odz,odzw,kb,in,iu,iv,iw,iu0,iv0
    call rdnc_zkb
end subroutine


subroutine set_surfwinds
!    open(72,file='prepglo/data/winds'//grd,form='unformatted')
    call rdnc_winds
    do n=1,12
	! read wind stress (dynes/cm-cm), temperature (deg c), salinity(ppt),
	! q, wavg(mm/month)
!	read(72) ((taux(i,j,n),i=1,i2),j=1,j2),((tauy(i,j,n),i=1,i2),j=1,j2)

	! screen out land sources and convert from newtons/m-m to dynes/cm-cm
	!   intensifies the mixing le
	nocn=0
	do j=1,j2
	    do i=1,i2
		nocn=nocn+in(i+1,j+1,1)
		taux(i,j,n)=10.*in(i+1,j+1,1)*taux(i,j,n)
		tauy(i,j,n)=10.*in(i+1,j+1,1)*tauy(i,j,n)
	    enddo
	enddo
    enddo

    nt=nocn
    call mpi_allreduce(nt,nocn,1,mpi_integer,mpi_sum,m_cart,ierr)
end subroutine


subroutine set_area
    real*8 dtmp
    do k=1,k1
	a(k)=0.d0
	do j=2,j1
	    dtmp=dx(j)*dy(j)
	    do i=2,i1
		a(k)=a(k)+in(i,j,k)*dtmp
	    enddo
	enddo
	dtmp=a(k)
	call mpi_allreduce(dtmp,a(k),1,mpi_real8,mpi_sum,m_cart,ierr)
    enddo
end subroutine


subroutine set_ts_nudging 
    ! simple equatorial and high latitude sponge layers are needed because
    ! storm-generated penetrative convection in high lats and inappropriate
    ! hydrostati! equatorial dynamics that cannot emulate effe!   of stacked
    ! equatorial vortices on t,s near equator; the effe!   on momentum are
    ! addressed by new eddy viscosity formulation based on the time
    ! variability of vertical velocity which is big: near the equator (due to
    ! small coriolis parameter); in high lats and southern ocean due to small
    ! stratification there; and along fronts everywhere due to instabilities
    ! ======================================================================
    ! jacc,jeq,jar! refer to cell centers. assume southernmost and
    ! northernmost gridlines are at the same lats as used at 1/4 deg lon
    ! resolution, which is achieved by proper inputs, defined in
    ! prepglo/data/metgen.
    ! 1. keep the same southernmost gridline lat (y0deg=yv(1)=-69.9377 deg)
    !    as used with 1/4 deg resolution (dxmnut=15.)
    ! 2. set the new j2 value equal to the old one (792) multiplied by the
    !    ratio of the old "dxdeg" value to the new one,
    ! this puts a gridline exactly on the equator and puts the northernmost
    ! gridline, yv(j1), at 69.7654 deg lat; 1/4 deg res (dxmnut=15.), gives
    ! yv(398)= 0.0 deg (on the equator).
    ! ======================================================================
    ! ---------------------
    ! setting nudge matrix
    ! ---------------------
    scr=0.
    yscale=10.
    rate=(792./j2t)/yscale
    rate3=3.*rate
    
    ! s
    n_tmp=90.
    tmpp_s=max(1./(n_tmp*daodt),0.)
    ! t
    n_tmp=90.
    tmpp_t=max(1./(n_tmp*daodt),0.)
    ! s,t coast
    n_tmp=-1.
    tmpp_sc=max(1./(n_tmp*daodt),0.)
    n_tmp=90.
    tmpp_tc=max(1./(n_tmp*daodt),0.)
    ! s,t ac! arc
    n_tmp=90.
    tmpp_shs=max(1./(n_tmp*daodt),0.)
    n_tmp=90.
    tmpp_ths=max(1./(n_tmp*daodt),0.)


    ! hmean ac! arc
    n_tmp=90.
    tmpp_hi=max(1./(n_tmp*daodt),0.)
    ! hmean eq
    n_tmp=-1.
    tmpp_ei=max(1./(n_tmp*daodt),0.)
    ! hmean other
    n_tmp=-1.
    tmpp_o=max(1./(n_tmp*daodt),0.)
    ! hmean coast
    n_tmp=-1.
    tmpp_oc=max(1./(n_tmp*daodt),0.)


    acc_deg=-69.5
    acc_rag=acc_deg+2.

    arc_deg=83.
    arc_rag=arc_deg-2.

    do j=2,j1
	do i=2,i1
	    ! s & t
	    if(ydeg(j).ge.arc_deg.or.ydeg(j).le.acc_deg) then
		t_nudge(i,j)=tmpp_ths
		s_nudge(i,j)=tmpp_shs
		elseif((ydeg(j).ge.arc_rag)) then
		t_nudge(i,j)=(tmpp_ths-tmpp_t)*exp(-rate*abs((ydeg(j)-arc_deg)/(ydeg(j)-ydeg(j-1))))+tmpp_t
		s_nudge(i,j)=(tmpp_shs-tmpp_s)*exp(-rate*abs((ydeg(j)-arc_deg)/(ydeg(j)-ydeg(j-1))))+tmpp_s
		elseif((ydeg(j).le.acc_rag)) then
		t_nudge(i,j)=(tmpp_ths-tmpp_t)*exp(-rate*abs((ydeg(j)-acc_deg)/(ydeg(j)-ydeg(j-1))))+tmpp_t
		s_nudge(i,j)=(tmpp_shs-tmpp_s)*exp(-rate*abs((ydeg(j)-acc_deg)/(ydeg(j)-ydeg(j-1))))+tmpp_s
	    else
		t_nudge(i,j)=tmpp_t
		s_nudge(i,j)=tmpp_s
	    endif

	    if(kb(i,j).le.1.or.kb(i-1,j).le.1.or.kb(i+1,j).le.1.or.kb(i,j-1).le.1.or.kb(i,j+1).le.1) then
		t_nudge(i,j)=tmpp_tc
		s_nudge(i,j)=tmpp_sc
	    endif

	    ! hmean
	    do k=1,k0
		if(ydeg(j).ge.arc_deg.or.ydeg(j).le.acc_deg) then
		    hmean_nudge(i,j,k)=tmpp_hi
		    elseif((ydeg(j).ge.arc_rag)) then
		    hmean_nudge(i,j,k)=(tmpp_hi-tmpp_o)*exp(-rate*abs((ydeg(j)-arc_deg)/(ydeg(j)-ydeg(j-1))))+tmpp_o
		    elseif((ydeg(j).le.acc_rag)) then
		    hmean_nudge(i,j,k)=(tmpp_hi-tmpp_o)*exp(-rate*abs((ydeg(j)-acc_deg)/(ydeg(j)-ydeg(j-1))))+tmpp_o
		    elseif(ydeg(j).ge.-5..and.ydeg(j).le.5.) then
		    hmean_nudge(i,j,k)=(tmpp_ei-tmpp_o)*exp(-rate3*abs((ydeg(j)-0.)/(ydeg(j)-ydeg(j-1))))+tmpp_o
		else
		    hmean_nudge(i,j,k)=tmpp_o
		endif

		if(kb(i,j).le.1.or.kb(i-1,j).le.1.or.kb(i+1,j).le.1.or.kb(i,j-1).le.1.or.kb(i,j+1).le.1) then
		    hmean_nudge(i,j,k)=tmpp_oc
		endif
	    enddo
	enddo
    enddo
end subroutine


subroutine set_lateral_mixing 
    ! increase lateral mixing near equator and in high latitudes
    ! big northern latitude dissipation fixes initialization data problem
    it=mylon*i2
    jt=mylat*j2

    dxx=1.0
    dyy=1.0

    do k=1,k1
	do j=1,j0
	    do i=1,i0
		temp=1.e6*(exp(-.001*(2-j-jt)**2)+exp(-.001*(j1t-jt-j)**2))      ! north & south poles
		scr(i,j,k)=2.e6*exp(-.001*((ydeg(j)-0)/dyy)**2)+temp        &
		+.5e6*exp(-.01*(((xdeg(i)-109.5))**2+((ydeg(j)-15.5))**2)) &! warm bias
		+.5e6*exp(-.01*(((xdeg(i)-114.5))**2+((ydeg(j)+ 9.5))**2)) &! warm bias
		+2.e6*exp(-.01*(((xdeg(i)-300.9))**2+((ydeg(j)-45.11))**2))&! cold bias (west nova scotia)
		+2.e6*exp(-.01*(((xdeg(i)-298.7))**2+((ydeg(j)-43.68))**2))&! cold bias (sw nova scotia)
		+2.e6*exp(-.01*(((xdeg(i)-276.2))**2+((ydeg(j)-56.67))**2))&! cold bias (south hudson bay)
		+2.e6*exp(-.01*(((xdeg(i)-181.7))**2+((ydeg(j)-64.22))**2))&! cold bias (gulf of anadyr)
		+2.e6*exp(-.01*(((xdeg(i)-164.8))**2+((ydeg(j)-58.85))**2))&! cold bias (ne kamchatka)
		+1.e6*exp(-.01*(((xdeg(i)-142.3))**2+((ydeg(j)-56.67))**2))&! cold bias (sakhalin gulf)
		+1.e6*exp(-.01*(((xdeg(i)-143.4))**2+((ydeg(j)-45.11))**2))&! cold bias (aniva bay)
		+1.e6*exp(-.01*(((xdeg(i)-135.6))**2+((ydeg(j)-42.21))**2)) ! cold bias (east japan sea)
	    enddo
	enddo
    enddo

    ! note: the following may need some adjustment in the mpi code
    ! increase lateral eddy viscosity near wet/dry point interfaces
    do k=1,k1
	do i=1,i0
	    do j=1,j0
		dmx(i,j,k)=dm0+scr(i,j,k)
		dmy(i,j,k)=dmx(i,j,k)
	    enddo
	enddo
    enddo
end subroutine


subroutine set_extra_vertical_eddy_viscosity
    ! xtra vertical eddy viscosity
    do jj=2,j1
	j=min(jj,j2)
	j=max(j,3)
	do ii=2,i1
	    i=min(ii,i2)
	    i=max(i,3)
	    ! bottom slope
	    dzdx=(z(2*kb(i+1,j))-z(2*kb(i-1,j)))*odx(j) 
	    dzdy=(z(2*kb(i,j+1))-z(2*kb(i,j-1)))*ody(j)
	    zbot=z(2*kb(i,j)+1)
	    slope=sqrt(dzdx**2+dzdy**2)
	    do k=1,k2
		! regionally increased vertical eddy viscosity
		! locally increase add array to get 90-day time scale for 1 layer scale
		! tscale=90.*24.*3600.
		! after day 100 yr 3 confine xtra equatorial (j=398) mixing to top 1000m
		tmp=min(10.,2.e-4*z(2*kb(i,j)+1))
		!      evadd=.5*exp(-tmp-.01*(jt+j-j2t/2-2)**2)
		evadd=.5*exp(-tmp-.01*((ydeg(j)-0)/(ydeg(j+1)-ydeg(j)))**2)

		! eddy viscosity is big above steeply sloping bottom locations
		! near-bottom factor
		! temp=2.e-9*(zbot-z(2*k+1))**2
		! decrease vertical mixing scale at start of yr 3 to allow free
		! passage into caribbean
		temp=1.e-8*(zbot-z(2*k+1))**2
		! vertical mix time scale=o(dz*2/ev); for time scale=1 day and dz=10m,
		! ev=o(10**6/10**5)=10.

		tmp=min(5.,5.*exp(-temp)*slope)

		! biggest eddy viscosity in deep water having steeply sloping bottom
		! where breaking big amplitude internal waves are most likely to occur
		! due to low subthermocline stratification

		tmp=tmp+vbk(k)+evadd

		! explicit stability limit and ev normalization done in pp82glo
		add(ii-1,jj-1,k)=tmp
	    enddo 
	enddo 
    enddo 
end subroutine


subroutine set_seasonal_data 
    real*8 tmpp
!    open(70,file='prepglo/data/annualevitus'//grd,form='unformatted')
!    rewind 70
    ! seasonal boundary conditions
    tmpp=1./12.

    do nm=1,12
	if(myid.eq.0)write(*,15) nm
	15   format('month',i3,' in bounds')
	! read in monthly t,s boundary climatologies
!	read(70) sclim,tclim
        call rdnc_levitus(nm)

	do k=1,k1
	    do j=1,j0
		do i=1,i0
		    s2(i,j,k)=s2(i,j,k)+sclim(i,j,k)*tmpp
		    t2(i,j,k)=t2(i,j,k)+tclim(i,j,k)*tmpp
		    if(k.eq.1) sclim(i,j,k)=max(33.,sclim(i,j,k))
		    if(ini_clim.eq.13)then
			s1(i,j,k)=s1(i,j,k)+sclim(i,j,k)*tmpp
			t1(i,j,k)=t1(i,j,k)+tclim(i,j,k)*tmpp
		    else
			if(ini_clim.eq.nm)then
			    s1(i,j,k)=sclim(i,j,k)
			    t1(i,j,k)=tclim(i,j,k)
			endif
		    endif
		enddo
	    enddo
	enddo

	do j=1,j0
	    do i=1,i0
		do k=1,k1
		    sclis(i,j,k,nm)=sclim(i,j,k)
		    tclis(i,j,k,nm)=tclim(i,j,k)
		enddo
	    enddo
	enddo
    enddo 
!    close(70)
    
    do k=1,k1
	do j=1,j0
	    do i=1,i0
		sclim(i,j,k)=s2(i,j,k)
		tclim(i,j,k)=t2(i,j,k)
		s2(i,j,k)=0.
		t2(i,j,k)=0.
	    enddo 
	enddo 
    enddo 
 
    if (lrstrt.eq.0)then
	if(myid.eq.0)write(*,*)'write ncf initial'
	call write_init2nc(sclim,tclim)

	do k=1,k1
	    do j=1,j0
		do i=1,i0
		    slf(i,j,k)=s1(i,j,k)
		    tlf(i,j,k)=t1(i,j,k)
		    s2(i,j,k)=s1(i,j,k)
		    t2(i,j,k)=t1(i,j,k)
!yc
                    c1(i,j,k)=0. !t1(i,j,k) !assume c1 is same as temp
                    clf(i,j,k)=c1(i,j,k)
                    c2(i,j,k)=c1(i,j,k)
!yc
		enddo
	    enddo
	enddo

        write(14,1)case,dscrib,lrstrt,mxit,i0,j0,k0,ktrm,dt,fltw,tlz,g,dm0,de0,dmz0,rzmx,drag
        1    format('case ',a5,' (',a63,')'/'----------'/ &
        ' lrstrt   mxit   i0   j0   k0 ktrm'/2i9,4i5// &
        8x,'dt',6x,'fltw',7x,'tlz'/1p,3(1x,e9.2)// &
        7x,'g',7x,'dm0',7x,'de0',6x,'dmz0', &
        6x,'rzmx',6x,'drag'/1p,7(1x,e9.2))

	! initialize qavg
	n=30*daodt
	temp=z(3)/(30.*24.*3600.)
	do m=1,12
	    nsombo(m)=2
	    mp=m+1
	    if (m.eq.12) mp=1
	    evap=0.
	    do j=1,j2
		tmp=dx(j+1)*dy(j+1)
		do i=1,i2
		    ! this initialization is arbitrary and optional
		    qavg(i,j,m)=in(i+1,j+1,1)*(tclis(i+1,j+1,1,mp)-tclis(i+1,j+1,1,m))/n
		enddo	
	    enddo	
	enddo	

	! depth
!	read(60) ((scr(i,j,1),i=1,i0),j=1,j0)
!	close(60)
        call rdnc_zkb2
	do j=2,j1
	    do i=2,i1
		! convert to kilometers
		scr(i,j,1)=1.e-5*scr(i,j,1)
	    enddo
	enddo
    else
	! -----------------
	! read restart data
	! -----------------
	! av is for long term averages. avg is for nabcs 5-day averages.
	if (myid.eq.0) print*,'read rerun file',myid,itf
	call read_rerun

	call  mpi_exch_3d_r4(u)
	call  mpi_exch_3d_r4(v)
	call  mpi_exch_2d_r8(x)
	
	if(myid.eq.0)write(*,*)'read ncf initial'
	call read_initnc(sclim,tclim)

	call  mpi_exch_3d_r4(tclim)
	call  mpi_exch_3d_r4(sclim)

	n=days*daodt+.5
	daodt_new=daodt
	daodt_old=itf/days
	if (n/=itf) then 
	    ! tmp=(dtnew/dtold)
	    tmp=daodt_old/daodt_new
	    if (myid .eq. 0)then
		write(*,57) tmp
		write(14,57) tmp
		57   format('dtnew/dtold=',f6.3)
	    endif
	    do nm=1,12
		do j=1,j2
		    do i=1,i2
			qavg(i,j,nm)=tmp*qavg(i,j,nm)
		    enddo	
		enddo	
	    enddo	

	    if (myid .eq. 0)then
		write(*,59) itf,n
		write(14,59) itf,n
		59   format(/'time step # changed from',i7,' to',i7,' at restart to reflect dt change')
	    endif
	    itf=n
	    ! reset averaging arrays every 5 years
	endif
	
	if(myid.eq.0)write(*,61) itf,days,daodt_new,daodt_old
	61   format('itf,days,daodt_new,daodt_old',i8,3f8.1)

	n=days+.1

	if (mod(int(days+.01),1800).eq.0.) av=0.

	itf=days*daodt+0.5

	! 82 if (days.gt.460..and.mod(itf,360*itfday).ne.0) go to 85
	! wavg & qavg include advection effe!  
	! wavg is surface freshwater sources/sinks (e-p) consistent with model
	! and climatological data
	! nsombo(m)=number of ensemble-averaging years
	! assumed already included in qavg and wa[vg
	! however, we weight each new year equally to the previously completed
	! year's ensemble average, until year 11, thus given more weight to
	! recent results and thus acclerating convergence to a multi-year ensemble
	! average.
	if (n<=10*360) then 
	    ! keep old avg but speed convergence to new avg
	    do m=1,12
		nsombo(m)=2
	    enddo
	endif  
	
	if(myid .eq. 0)then
	    write(*,90) itf,days,nsombo
	    write(14,90) itf,days,nsombo
	    90   format(/'restart data read for time step ',i9,', day ',f8.2/'nsombo:',12i3)
	endif

    endif


    if (mylat .eq. 0) then
	do k=1,k1
	do i=1,i0
	t1(i,1,k)=tclim(i,1,k)
	t2(i,1,k)=tclim(i,1,k)
	tlf(i,1,k)=tclim(i,1,k)
	s1(i,1,k)=sclim(i,1,k)
	s2(i,1,k)=sclim(i,1,k)
	slf(i,1,k)=sclim(i,1,k)
	enddo
	enddo
    endif
    if (mylat .eq. ngy1) then
	do k=1,k1
	do i=1,i0
	t1(i,j0,k)=tclim(i,j0,k)
	t2(i,j0,k)=tclim(i,j0,k)
	tlf(i,j0,k)=tclim(i,j0,k)
	s1(i,j0,k)=sclim(i,j0,k)
	s2(i,j0,k)=sclim(i,j0,k)
	slf(i,j0,k)=sclim(i,j0,k)
	enddo
	enddo
    endif
end subroutine


subroutine eliminate_spurious_values
    if(lrstrt.eq.0)then
	! eliminate spurious values over land
	!      nit=2001
	nit=1
	do n=1,nit
	    do k=1,2
		tmin=1.e8
		tmax=-1.e8
		smin=1.e8
		smax=-1.e8
		do j=2,j1
		    do i=2,i1
			if (in(i,j,k)==0) then 
			    scr(i,j,1)=.25*(s1(i+1,j,k)+s1(i-1,j,k)+s1(i,j-1,k)+s1(i,j+1,k))
			    smin=min(smin,scr(i,j,1))
			    smax=max(smax,scr(i,j,1))
			    scr(i,j,2)=.25*(t1(i+1,j,k)+t1(i-1,j,k)+t1(i,j-1,k)+t1(i,j+1,k))
			    tmin=min(tmin,scr(i,j,2))
			    tmax=max(tmax,scr(i,j,2))
			endif	
		    enddo
		enddo
		do j=2,j1
		    do i=2,i1
			if (in(i,j,k)==0) then
			    s1(i,j,k)=scr(i,j,1)
			    t1(i,j,k)=scr(i,j,2)
			endif
		    enddo
		enddo

		call mpi_exch_2d_r4(t1(1,1,k))
		call mpi_exch_2d_r4(s1(1,1,k))

		if(myid.eq.0)then
		    if(mod(n,100).eq.1)write(*,2) n,k,tmin,tmax,smin,smax
		    2    format('n,k,tmin,tmax,smin,smax:',i5,i3,4f7.3)
		endif
	    enddo
	enddo
    endif

end subroutine


subroutine ts_nudging_for_watermass
    ! To emulate watermass transformation effects on THC, nudge toward
    ! eqatorial as well as hi and lo lat climatology. We do not directly
    ! address watermass tranformations closer to poles, because their most
    ! important effects are seen in the central ACC and Arctic Circle
    ! regions, but important wind forced dynamics are included closer to
    ! poles. This is valid for simulating present climate and model
    ! validation, but NOT for simulating climate change, which equires
    ! ocean coupling to full atmospheric and ice dynamics, and much more
    ! complex physics.

    ! NOTE: This gives internal effective heat and salt material "sources"
    ! not in the real ocean. However, this should give model climatology
    ! close to the real ocean with the "sources" magnitudes being a metric
    ! of model quality; they should decrease with increased resolution and
    ! better model physics and subgridscale mixing emulation. Thus, the time
    ! averaged source magnitudes are a good indicator of model applicability
    ! to climate research. Indeed, the following may eventually include time
    ! averaged "sources" to decrease damping of transients (analagous to the
    ! non-damping surface sources of Dietrich, et al., 2004), but are not
    ! simply surface boundary conditions needed for ocean-only modeling (may
    ! be derived by coupling "perfect" ocean and atmospheri! models). They
    ! also noted that the "sources" may be used to improve vertical mixing
    ! parameterizations needed to keep imperfect models on track with
    ! observed ocean climate.
    ! ======================================================================
    ! Later, check sensitivity to expanding nudging to bands of polar and
    ! eqatorial lats to invoke stronger watermass formation.
    ! Loop 95 replaces hardwired multi-region approach used by original
    ! coding before calling FSGLO.
    do k=2,k1
	do j=2,j1
	    do i=2,i1
		if(in(i,j,k)/=0) then
		    t1(i,j,k)=(1.-hmean_nudge(i,j,k))*t1(i,j,k)+hmean_nudge(i,j,k)*tclim(i,j,k)
		    s1(i,j,k)=(1.-hmean_nudge(i,j,k))*s1(i,j,k)+hmean_nudge(i,j,k)*sclim(i,j,k)
!yc
!                    c1(i,j,k)=(1.-hmean_nudge(i,j,k))*c1(i,j,k)+hmean_nudge(i,j,k)*tclim(i,j,k)
!yc
		endif
	    enddo
	enddo
    enddo

    do j=2,j1
	do i=2,i1
	    s1(i,j,1)=s1(i,j,1)+s_nudge(i,j)*in(i,j,1)*in(i,j,2)*&
			    (fold*sclis(i,j,1,nld)+fnew*sclis(i,j,1,new)-s1(i,j,1))
	enddo
    enddo
end subroutine


subroutine solve_hydrostatic_equation
    real*8 tmpd,sltd,pisd,p01,p02
    dimension rhomin(k1)
    ! ======================================================
    ! Hydrostatic equation with static balance subracted out
    ! ======================================================
    do k=1,k1
	temp=1.e8
	do j=2,j1
	    do i=2,i1
		tmpd=max(-1.8,t2(i,j,k))
		sltd=min(41.,s2(i,j,k))
		sltd=max(0.d0,sltd)
		pisd=.001*z(2*k)
		! =======================
		! dan wright full e.o.s.
		! rho(i,j,k)=r(tmpd,sltd,pisd)
		! function r(t,s,p)
		! =======================
		p02=1747.4508988+tmpd*(11.51588-0.046331033*tmpd)-sltd*(3.85429655+0.01353985*tmpd)
		p01=pisd+5884.81703666+tmpd*(39.803732+tmpd*(-0.3191477+tmpd*0.0004291133))+2.6126277*sltd
		rho(i,j,k)=p01/(p02+0.7028423*p01)
		! determine global min of rho in all layers
		if(in(i,j,k)/=0) then 
		    temp=min(dble(temp),rho(i,j,k))
		endif
	    enddo
	enddo
	tmp=temp
	call mpi_barrier(m_cart,ierr)
	timer(9)=mpi_wtime()
	call mpi_allreduce(tmp,temp,1,mpi_real8,mpi_min,m_cart,ierr)
	timer(3)=timer(3)+mpi_wtime()-timer(9)

	rhomin(k)=temp
    enddo

    ! Subtract off global min of rho (in all layers, but it is always
    ! in the top layer) to keep correct rho stratification for vertical eddy
    ! viscosity/diffusivity parameterization in PP82GLO. Later, to reduce
    ! pressure gradient errors associated with big mean density in underlying
    ! layers, we subtract the layer min rho from each layer.
    do k=1,k1
	do j=2,j1
	    do i=2,i1
		scr(i,j,k)=rho(i,j,k)-rhomin(k)
	    enddo
	enddo
    enddo

    ! Use deviation of density from layer min density to get hydrostatic
    ! pressure gradient (reduces roundoff error and does not affect
    ! horizontal pressure gradient in the primitive hydroatic equations)
    ! Non-centered top half layer approximation: a good approximation
    ! because D(RHO)/DZ is very small in mixed layer
    tmpd=g*z(2)
    do j=2,j1
	do i=2,i1
	    p(i,j,1)=p0(i,j)+tmpd*rho(i,j,1)
	enddo
    enddo
    
    do k=1,k2
	tmpd=.5*g/odzw(k+1)
	do j=2,j1
	    do i=2,i1
		p(i,j,k+1)=p(i,j,k)+tmpd*(rho(i,j,k)+rho(i,j,k+1))
	    enddo
	enddo
    enddo
end subroutine


subroutine compute_pressure_gradient(k)
    ! first pressure differences
    do j=2,j1
	! scratch array for 4th-order pressure gradient calculation
	! avoid misusing data from overlying levels
	! (put 0's where they are needed in xp)
	do i=1,i1
	    xp(i,j)=iu(i,j,k)*(p(i+1,j,k)-p(i,j,k))
	enddo
    enddo

    call mpi_sendrecv(xp(2,1),1,m_evlon2d,m_w,1,xp(i0,1),1,m_evlon2d,m_e,1,m_cart,istat,ierr)
    call mpi_sendrecv(xp(i2,1),1,m_evlon2d,m_e,1,xp(0,1),1,m_evlon2d,m_w,1,m_cart,istat,ierr)

    do j=js,jf
	do i=2,i1
	    yp(i,j)=iv(i,j,k)*(p(i,j+1,k)-p(i,j,k))
	enddo
    enddo

    call mpi_sendrecv(yp(1,2),1,m_evlat2d,m_s,1,yp(1,j0),1,m_evlat2d,m_n,1,m_cart,istat,ierr)
    call mpi_sendrecv(yp(1,j2),1,m_evlat2d,m_n,1,yp(1,0),1,m_evlat2d,m_s,1,m_cart,istat,ierr)

    if (k .eq. 1) then
	do j=2,j1
	    do i=0,i0
		xp(i,j)=iu0(i,j)*xp(i,j)
	    enddo
	enddo
	do j=js-1,jf+1
	    do i=2,i1
		yp(i,j)=iv0(i,j)*yp(i,j)
	    enddo
	enddo
    endif
    ! basi! second order
    do j=2,j1
	do i=2,i1
	    px(i,j)=6.*(xp(i,j)+xp(i-1,j))
	    py(i,j)=6.*(yp(i,j)+yp(i,j-1))
	enddo
    enddo
    ! fourth order correction
    do j=2,j1
	do i=2,i1
	    px(i,j)=px(i,j)+iu(i-1,j,k)*iu(i,j,k)*(xp(i,j)+xp(i-1,j)-xp(i+1,j)-xp(i-2,j))
	enddo
    enddo
    do j=js+1,jf
	do i=2,i1
	    py(i,j)=py(i,j)+iv(i,j-1,k)*iv(i,j,k)*(yp(i,j)+yp(i,j-1)-yp(i,j+1)-yp(i,j-2))
	enddo
    enddo
    do j=2,j1
	do i=2,i1
	    px(i,j)=o12*px(i,j)
	    py(i,j)=o12*py(i,j)
	enddo
    enddo
end subroutine


subroutine compute_vertical_fluxes(k,lt)
    ! ---------------
    ! vertical fluxes
    ! ---------------
    ! these are 4th-order-accurate in logical space ("stretched z-coord")
    ! which is logical when the stretching is near optimum.
    l=k+1
    if (k==k1) then 
	! set bottom fluxes to zero (drag law applied after label 500)
	do j=1,j2
	    do i=1,i2
		uz(i,j,lt)=0.
		vz(i,j,lt)=0.
		sz(i,j,lt)=0.
		tz(i,j,lt)=0.
!yc
                cfcz(i,j,lt)=0.
!yc
	    enddo
	enddo
    else	
	do j=2,j1
	    do i=2,i1
		scr(i,j,1)=6.*(u2(i,j,k)+u2(i,j,l))
		scr(i,j,2)=6.*(v2(i,j,k)+v2(i,j,l))
		scr(i,j,3)=6.*(s2(i,j,k)+s2(i,j,l))
		scr(i,j,4)=6.*(t2(i,j,k)+t2(i,j,l))
!yc
                scr(i,j,13)=6.*(c2(i,j,k)+c2(i,j,l))
!yc
	    enddo
	enddo

	if (k/=1.and.k/=k2) then 
	    do j=2,j1
		do i=2,i1
		    tmp=iw(i,j,k)*iw(i,j,k+2)
		    scr(i,j,1)=scr(i,j,1)+tmp*(-ulf(i,j,k-1)+u2(i,j,k)+u2(i,j,l)-u2(i,j,k+2))
		    scr(i,j,2)=scr(i,j,2)+tmp*(-vlf(i,j,k-1)+v2(i,j,k)+v2(i,j,l)-v2(i,j,k+2))
		    scr(i,j,3)=scr(i,j,3)+tmp*(-slf(i,j,k-1)+s2(i,j,k)+s2(i,j,l)-s2(i,j,k+2))
		    scr(i,j,4)=scr(i,j,4)+tmp*(-tlf(i,j,k-1)+t2(i,j,k)+t2(i,j,l)-t2(i,j,k+2))
!yc
                    scr(i,j,13)=scr(i,j,13)+tmp*(-clf(i,j,k-1)+c2(i,j,k)+c2(i,j,l)-c2(i,j,k+2))
!yc
		enddo
	    enddo
	endif

	do n=1,4
	    do j=2,j1
		do i=2,i1
		    scr(i,j,n)=o12*scr(i,j,n)
		enddo
	    enddo
	enddo
!yc
        do j=2,j1
            do i=2,i1
               scr(i,j,13)=o12*scr(i,j,13)
            enddo
        enddo
!yc


	do j=2,j1
	    do i=2,i1
		uz(i-1,j-1,lt)=w(i,j,l)*scr(i,j,1)-ev(i-1,j-1,k)*(u1(i,j,l)-u1(i,j,k))*iw(i,j,l)
		vz(i-1,j-1,lt)=w(i,j,l)*scr(i,j,2)-ev(i-1,j-1,k)*(v1(i,j,l)-v1(i,j,k))*iw(i,j,l)
		sz(i-1,j-1,lt)=w(i,j,l)*scr(i,j,3)-0.1*hv(i-1,j-1,k)*(s1(i,j,l)-s1(i,j,k))*iw(i,j,l)
		tz(i-1,j-1,lt)=w(i,j,l)*scr(i,j,4)-hv(i-1,j-1,k)*(t1(i,j,l)-t1(i,j,k))*iw(i,j,l)
!yc
                cfcz(i-1,j-1,lt)=w(i,j,l)*scr(i,j,13)-hv(i-1,j-1,k)*(c1(i,j,l)-c1(i,j,k))*iw(i,j,l)
!yc
	    enddo
	enddo
    endif

    ! set scr(3,j,5) to store u2(3,j,k) so as v2,s2 and t2
    call mpi_sendrecv(u2(3,1,k),1,m_vlon2d,m_w,1,scr(3,1,5),1,m_vlon2d,m_e,1,m_cart,istat,ierr)
    call mpi_sendrecv(v2(3,1,k),1,m_vlon2d,m_w,1,scr(3,1,6),1,m_vlon2d,m_e,1,m_cart,istat,ierr)
    call mpi_sendrecv(s2(3,1,k),1,m_vlon2d,m_w,1,scr(3,1,7),1,m_vlon2d,m_e,1,m_cart,istat,ierr)
    call mpi_sendrecv(t2(3,1,k),1,m_vlon2d,m_w,1,scr(3,1,8),1,m_vlon2d,m_e,1,m_cart,istat,ierr)
!yc
    call mpi_sendrecv(c2(3,1,k),1,m_vlon2d,m_w,1,scr(3,1,14),1,m_vlon2d,m_e,1,m_cart,istat,ierr)
!yc
    ! set scr(i2,j,5) to store u2(i2,j,k) so as v2,s2 and t2
    call mpi_sendrecv(u2(i2,1,k),1,m_vlon2d,m_e,1,scr(i2,1,5),1,m_vlon2d,m_w,1,m_cart,istat,ierr)
    call mpi_sendrecv(v2(i2,1,k),1,m_vlon2d,m_e,1,scr(i2,1,6),1,m_vlon2d,m_w,1,m_cart,istat,ierr)
    call mpi_sendrecv(s2(i2,1,k),1,m_vlon2d,m_e,1,scr(i2,1,7),1,m_vlon2d,m_w,1,m_cart,istat,ierr)
    call mpi_sendrecv(t2(i2,1,k),1,m_vlon2d,m_e,1,scr(i2,1,8),1,m_vlon2d,m_w,1,m_cart,istat,ierr)
!yc
    call mpi_sendrecv(c2(i2,1,k),1,m_vlon2d,m_e,1,scr(i2,1,14),1,m_vlon2d,m_w,1,m_cart,istat,ierr)
!yc
    ! set scr(i,3,9) to store u2(i,3,k) so as v2,s2 and t2
    call mpi_sendrecv(u2(1,3,k),1,m_vlat2d,m_s,1,scr(1,3,9),1,m_vlat2d,m_n,1,m_cart,istat,ierr)
    call mpi_sendrecv(v2(1,3,k),1,m_vlat2d,m_s,1,scr(1,3,10),1,m_vlat2d,m_n,1,m_cart,istat,ierr)
    call mpi_sendrecv(s2(1,3,k),1,m_vlat2d,m_s,1,scr(1,3,11),1,m_vlat2d,m_n,1,m_cart,istat,ierr)
    call mpi_sendrecv(t2(1,3,k),1,m_vlat2d,m_s,1,scr(1,3,12),1,m_vlat2d,m_n,1,m_cart,istat,ierr)
!yc
    call mpi_sendrecv(c2(1,3,k),1,m_vlat2d,m_s,1,scr(1,3,15),1,m_vlat2d,m_n,1,m_cart,istat,ierr)
!yc
    ! set scr(i,j2,9) to store u2(i,j2,k) so as v2,s2 and t2
    call mpi_sendrecv(u2(1,j2,k),1,m_vlat2d,m_n,1,scr(1,j2,9),1,m_vlat2d,m_s,1,m_cart,istat,ierr)
    call mpi_sendrecv(v2(1,j2,k),1,m_vlat2d,m_n,1,scr(1,j2,10),1,m_vlat2d,m_s,1,m_cart,istat,ierr)
    call mpi_sendrecv(s2(1,j2,k),1,m_vlat2d,m_n,1,scr(1,j2,11),1,m_vlat2d,m_s,1,m_cart,istat,ierr)
    call mpi_sendrecv(t2(1,j2,k),1,m_vlat2d,m_n,1,scr(1,j2,12),1,m_vlat2d,m_s,1,m_cart,istat,ierr)
!yc
    call mpi_sendrecv(c2(1,j2,k),1,m_vlat2d,m_n,1,scr(1,j2,15),1,m_vlat2d,m_s,1,m_cart,istat,ierr)
!yc
end subroutine


subroutine compute_longitudinal_fluxes(k)
    ! -------------------
    ! longitudinal fluxes
    ! -------------------
    do j=2,j1
	! higher order interpolations
	do i=1,i1
	    scr(i,j,1)=6.*(u2(i,j,k)+u2(i+1,j,k))
	    scr(i,j,2)=6.*(v2(i,j,k)+v2(i+1,j,k))
	    scr(i,j,3)=6.*(s2(i,j,k)+s2(i+1,j,k))
	    scr(i,j,4)=6.*(t2(i,j,k)+t2(i+1,j,k))
!yc
            scr(i,j,13)=6.*(c2(i,j,k)+c2(i+1,j,k))
!yc
	enddo

	tmp=iu(0,j,k)*iu(2,j,k)
	scr(1,j,1)=scr(1,j,1)+tmp*(-scr(i2,j,5)+u2(1,j,k)+u2(2,j,k)-u2(3,j,k))
	scr(1,j,2)=scr(1,j,2)+tmp*(-scr(i2,j,6)+v2(1,j,k)+v2(2,j,k)-v2(3,j,k))
	scr(1,j,3)=scr(1,j,3)+tmp*(-scr(i2,j,7)+s2(1,j,k)+s2(2,j,k)-s2(3,j,k))
	scr(1,j,4)=scr(1,j,4)+tmp* (-scr(i2,j,8)+t2(1,j,k)+t2(2,j,k)-t2(3,j,k))
!yc
        scr(1,j,13)=scr(1,j,13)+tmp*(-scr(i2,j,14)+c2(1,j,k)+c2(2,j,k)-c2(3,j,k))
!yc
	tmp=iu(i2,j,k)*iu(i0,j,k)
	scr(i1,j,1)=scr(i1,j,1)+tmp*(-u2(i2,j,k)+u2(i1,j,k)+u2(i0,j,k)-scr(3,j,5))
	scr(i1,j,2)=scr(i1,j,2)+tmp*(-v2(i2,j,k)+v2(i1,j,k)+v2(i0,j,k)-scr(3,j,6))
	scr(i1,j,3)=scr(i1,j,3)+tmp*(-s2(i2,j,k)+s2(i1,j,k)+s2(i0,j,k)-scr(3,j,7))
	scr(i1,j,4)=scr(i1,j,4)+tmp*(-t2(i2,j,k)+t2(i1,j,k)+t2(i0,j,k)-scr(3,j,8))
!yc
        scr(i1,j,13)=scr(i1,j,13)+tmp*(-c2(i2,j,k)+c2(i1,j,k)+c2(i0,j,k)-scr(3,j,14))
!yc

	do i=2,i2
	    tmp=iu(i-1,j,k)*iu(i+1,j,k)
	    scr(i,j,1)=scr(i,j,1)+tmp*(-u2(i-1,j,k)+u2(i,j,k)+u2(i+1,j,k)-u2(i+2,j,k))
	    scr(i,j,2)=scr(i,j,2)+tmp*(-v2(i-1,j,k)+v2(i,j,k)+v2(i+1,j,k)-v2(i+2,j,k))
	    scr(i,j,3)=scr(i,j,3)+tmp*(-s2(i-1,j,k)+s2(i,j,k)+s2(i+1,j,k)-s2(i+2,j,k))
	    scr(i,j,4)=scr(i,j,4)+tmp*(-t2(i-1,j,k)+t2(i,j,k)+t2(i+1,j,k)-t2(i+2,j,k))
!yc
            scr(i,j,13)=scr(i,j,13)+tmp*(-c2(i-1,j,k)+c2(i,j,k)+c2(i+1,j,k)-c2(i+2,j,k))
!yc
	enddo

	do n=1,4
	    do i=1,i1
		scr(i,j,n)=o12*scr(i,j,n)
	    enddo
	enddo
!yc
        do i=1,i1
            scr(i,j,13)=o12*scr(i,j,13)
        enddo
!yc

	do i=1,i1
	    am=dmx(i,j,k)*odx(j)
	    ah=dhx(i,j)*odx(j)
	    ! ux could be simply u(i,j,k)**2 + ...., but this does not conserve
	    ! energy and is unstable
	    ux(i,j-1)=(u(i,j,k)*scr(i,j,1)-am*(u1(i+1,j,k)-u1(i,j,k)))*iu(i,j,k)
	    vx(i,j-1)=(u(i,j,k)*scr(i,j,2)-am*(v1(i+1,j,k)-v1(i,j,k)))*iu(i,j,k)
	    sx(i,j-1)=(u(i,j,k)*scr(i,j,3)-ah*(s1(i+1,j,k)-s1(i,j,k)))*iu(i,j,k)
	    tx(i,j-1)=(u(i,j,k)*scr(i,j,4)-ah*(t1(i+1,j,k)-t1(i,j,k)))*iu(i,j,k)
!yc
            cfcx(i,j-1)=(u(i,j,k)*scr(i,j,13)-ah*(c1(i+1,j,k)-c1(i,j,k)))*iu(i,j,k)
!yc
	enddo
    enddo
end subroutine


subroutine compute_latitudinal_fluxes(k)
    ! ------------------
    ! latitudinal fluxes
    ! ------------------
    do j=1,jf
	do i=2,i1
	    scr(i,j,1)=6.*(u2(i,j,k)+u2(i,j+1,k))
	    scr(i,j,2)=6.*(v2(i,j,k)+v2(i,j+1,k))
	    scr(i,j,3)=6.*(s2(i,j,k)+s2(i,j+1,k))
	    scr(i,j,4)=6.*(t2(i,j,k)+t2(i,j+1,k))
!yc
            scr(i,j,13)=6.*(c2(i,j,k)+c2(i,j+1,k))
!yc
	enddo
    enddo

    do j=2,j2
	do i=2,i1
	    tmp=iv(i,j-1,k)*iv(i,j+1,k)
	    scr(i,j,1)=scr(i,j,1)+tmp*(-u2(i,j-1,k)+u2(i,j,k)+u2(i,j+1,k)-u2(i,j+2,k))
	    scr(i,j,2)=scr(i,j,2)+tmp*(-v2(i,j-1,k)+v2(i,j,k)+v2(i,j+1,k)-v2(i,j+2,k))
	    scr(i,j,3)=scr(i,j,3)+tmp*(-s2(i,j-1,k)+s2(i,j,k)+s2(i,j+1,k)-s2(i,j+2,k))
	    scr(i,j,4)=scr(i,j,4)+tmp*(-t2(i,j-1,k)+t2(i,j,k)+t2(i,j+1,k)-t2(i,j+2,k))
!yc
            scr(i,j,13)=scr(i,j,13)+tmp*(-c2(i,j-1,k)+c2(i,j,k)+c2(i,j+1,k)-c2(i,j+2,k))
!yc
	enddo
    enddo

    do i=2,i1
	tmp=iv(i,0,k)*iv(i,2,k)
	scr(i,1,1)=scr(i,1,1)+tmp*(-scr(i,j2,9)+u2(i,1,k)+u2(i,2,k)-u2(i,3,k))
	scr(i,1,2)=scr(i,1,2)+tmp*(-scr(i,j2,10)+v2(i,1,k)+v2(i,2,k)-v2(i,3,k))
	scr(i,1,3)=scr(i,1,3)+tmp*(-scr(i,j2,11)+s2(i,1,k)+s2(i,2,k)-s2(i,3,k))
	scr(i,1,4)=scr(i,1,4)+tmp*(-scr(i,j2,12)+t2(i,1,k)+t2(i,2,k)-t2(i,3,k))
!yc
        scr(i,1,13)=scr(i,1,13)+tmp*(-scr(i,j2,15)+c2(i,1,k)+c2(i,2,k)-c2(i,3,k))
!yc
    enddo

    do i=2,i1
	tmp=iv(i,j2,k)*iv(i,j0,k)
	scr(i,j1,1)=scr(i,j1,1)+tmp*(-u2(i,j2,k)+u2(i,j1,k)+u2(i,j0,k)-scr(i,3,9))
	scr(i,j1,2)=scr(i,j1,2)+tmp*(-v2(i,j2,k)+v2(i,j1,k)+v2(i,j0,k)-scr(i,3,10))
	scr(i,j1,3)=scr(i,j1,3)+tmp*(-s2(i,j2,k)+s2(i,j1,k)+s2(i,j0,k)-scr(i,3,11))
	scr(i,j1,4)=scr(i,j1,4)+tmp*(-t2(i,j2,k)+t2(i,j1,k)+t2(i,j0,k)-scr(i,3,12))
!yc
        scr(i,j1,13)=scr(i,j1,13)+tmp*(-c2(i,j2,k)+c2(i,j1,k)+c2(i,j0,k)-scr(i,3,15))
!yc
    enddo

    do j=js,jf
	do n=1,4
	    do i=2,i1
		scr(i,j,n)=o12*scr(i,j,n)
	    enddo
	enddo
    enddo
!yc
    do j=js,jf
        do i=2,i1
            scr(i,j,13)=o12*scr(i,j,13)
        enddo
    enddo
!yc

    do j=js,jf
	do i=2,i1
	    am=dmy(i,j,k)*odyv(j)
	    ah=dhy(i,j)*odyv(j)
	    !       am=dmx(i,j,k)*odyv(j)
	    !       ah=dhx(i,j)*odyv(j)
	    ! old free-slip
	    uy(i-1,j)=csv(j)*(v(i,j,k)*scr(i,j,1)-am*(u1(i,j+1,k)-u1(i,j,k)))*iv(i,j,k)
	    vy(i-1,j)=csv(j)*(v(i,j,k)*scr(i,j,2)-am*(v1(i,j+1,k)-v1(i,j,k)))*iv(i,j,k)
	    sy(i-1,j)=csv(j)*(v(i,j,k)*scr(i,j,3)-ah*(s1(i,j+1,k)-s1(i,j,k)))*iv(i,j,k)
	    ty(i-1,j)=csv(j)*(v(i,j,k)*scr(i,j,4)-ah*(t1(i,j+1,k)-t1(i,j,k)))*iv(i,j,k)
!yc
            cfcy(i-1,j)=csv(j)*(v(i,j,k)*scr(i,j,13)-ah*(c1(i,j+1,k)-c1(i,j,k)))*iv(i,j,k)
!yc
	enddo
    enddo

end subroutine


end module