#include "options.inc"


      subroutine cobc (j)
      use spflame_module
      implicit none
      integer, intent(in) :: j
      integer i,k,n,is,ie
      real var,var1
      real diff_fb(is_pe:ie_pe,0:km)
      real ADV_Ty(is_pe:ie_pe,km)
      real vad(is_pe:ie_pe,km),phase_vel(is_pe:ie_pe,km,nt)
c
c-----------------------------------------------------------------------
c     at open boundaries always laplacian formulation is used (at tau)
c     (same as for interior but for tau !!)
c
c     notice: I have changed the standard code here and have
c             used now scaling for the diffusivity
c
c     SPFLAME version : c.eden
c-----------------------------------------------------------------------
c
      real diff_fe_ha(is_pe-1:ie_pe,km,j)
      real diff_tx,diff_ty,diff_tz
      real T_i,T_j
      integer ip,jp
#ifdef partial_cell
      real dhte,dhtn,dhwe,dhwn
      dhte(i,k,j)   = min(dht(i+1,k,j),dht(i,k,j))
      dhtn(i,k,j)   = min(dht(i,k,j+1),dht(i,k,j))
      dhwe(i,k,j)   = min(dhwt(i+1,k,j),dhwt(i,k,j))
      dhwn(i,k,j)   = min(dhwt(i,k,j+1),dhwt(i,k,j))

      T_i(i,k,j,n,ip) = t(i+ip,max(1,k-1),j,n,tau) - dhwe(i,k-1,j)
     &            *(t(i+ip,max(1,k-1),j,n,tau) - t(i+ip,k,j,n,tau))
     &               /dhwt(i+ip,k-1,j)
      T_j(i,k,j,n,jp) = t(i,max(1,k-1),j+jp,n,tau) - dhwn(i,k-1,j)
     &            *(t(i,max(1,k-1),j+jp,n,tau) - t(i,k,j+jp,n,tau))
     &               /dhwt(i,k-1,j+jp)
#else
      T_i(i,k,j,n,ip) = t(i+ip,k,j,n,tau)
      T_j(i,k,j,n,jp) = t(i,k,j+jp,n,tau)
#endif
      DIFF_Tx(i,k,j) = (diff_fe_ha(i,  k,j)*tmask(i+1,k,j)
     &                - diff_fe_ha(i-1,k,j)*tmask(i-1,k,j))*cstdxtr(i,j)
#ifdef partial_cell
     &                             /dht(i,k,j)
#endif
      DIFF_Ty(i,k,j) = ahc_north_ha(j)*tmask(i,k,j+1)*
     &                     (t_j(i,k,j,n,1) - t_j(i,k,j,n,0))
     &                         *ah_scale(i,j)
#ifdef partial_cell
     &                         /dht(i,k,j)*dhtn(i,k,j)
#endif
     &                   - ahc_south_ha(j)*tmask(i,k,j-1)*
     &                     (t_j(i,k,j,n,0) - t_j(i,k,j,n,-1))
     &                         *ah_scale(i,j)
#ifdef partial_cell
     &                         /dht(i,k,j)*dhtn(i,k,j-1)
#endif
      DIFF_Tz(i,k) = (diff_fb(i,k-1) - diff_fb(i,k))
#ifdef partial_cell
     &               /dht(i,k,j)
#else
     &               *dztr(k)
#endif
      real obcs,obcn,data

      obcs=1./(86400.*5.); obcn=1./(86400.*5.) ! hardwired 5 day time scale
 
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
c
      if (j .eq. 2) then
c
c       my_blk_j was already tested
c
c-----------------------------------------------------------------------
c       1) compute the advective velocity "vad" 
c          at the north face of the "t" grid box
c
c       2) compute phase velocity at the southern boundary: c1s
c-----------------------------------------------------------------------
c
        var = -dyu(j+1)/dt 
        do k=1,km
          do i=is,ie
             vad(i,k)=(u(i,k,j,2,tau)*dxt(i)+u(i-1,k,j,2,tau)
     &            *dxt(i-1))/(dxt(i)+dxt(i-1))*tmask(i,k,j)
             if (vad(i,k).gt.0) vad(i,k)=0.
c	     
	     do n=1,nt
	      var1=t(i,k,j+2,n,taum1)-t(i,k,j+1,n,taum1)
	      if (var1.eq.0.) then
		phase_vel(i,k,n)=var
	      else
		phase_vel(i,k,n)=var*(t(i,k,j+1,n,tau)-t(i,k,j+1,n,taum1))
     &                            /var1*tmask(i,k,j)
                if (phase_vel(i,k,n).gt. 0.) phase_vel(i,k,n)=0.
		if (phase_vel(i,k,n).lt.var) phase_vel(i,k,n)=var
	      endif
	     enddo
          enddo
        enddo
      endif
c
      if (j .eq. jmt-1) then
c
c-----------------------------------------------------------------------
c       1) compute the advective velocity "vad" 
c          at the south face of the "t" grid box
c
c       2) compute phase velocity at the northern boundary: c1n
c-----------------------------------------------------------------------
c
        var=dyu(j-1)/dt
        do k=1,km
          do i=is,ie
           vad(i,k)=(u(i,k,j-1,2,tau)*dxt(i)+u(i-1,k,j-1,2,tau)
     &             *dxt(i-1))/(dxt(i)+dxt(i-1))*tmask(i,k,j)
           if (vad(i,k).lt.0) vad(i,k)=0.
           do n=1,nt
            var1=t(i,k,j-1,n,taum1)-t(i,k,j-2,n,taum1)
	    if (var1.eq.0.) then
 	      phase_vel(i,k,n)=var
	    else
 	      phase_vel(i,k,n)=-var*(t(i,k,j-1,n,tau)-t(i,k,j-1,n,taum1))
     &                           /var1*tmask(i,k,j)
              if (phase_vel(i,k,n).lt.  0) phase_vel(i,k,n)=0.
	      if (phase_vel(i,k,n).gt.var) phase_vel(i,k,n)=var
	    endif
	   enddo
          enddo
        enddo
      endif
c
      do n=1,nt
c
c-----------------------------------------------------------------------
c           diffusive flux across eastern face of "T" cells
c           (use constant horizontal diffusion)
c-----------------------------------------------------------------------
c
        do k=1,km
         do i=is-1,ie
          diff_fe_ha(i,k,j)  = ah_cstdxur(i,j)*
     &                     (t_i(i,k,j,n,1) - t_i(i,k,j,n,0))
     &                         *ah_scale(i,j)
#ifdef partial_cell
     &                         *dhte(i,k,j)
#endif
         enddo
        enddo
c
c-----------------------------------------------------------------------
c           diffusive flux across bottom face of "T" cells
c           use constant vertical diffusion
c-----------------------------------------------------------------------
c
        do k=1,km-1
          do i=is,ie
             diff_fb(i,k) =  kappa_h*
     &                         (t(i,k,j,n,tau) - t(i,k+1,j,n,tau))
#ifdef partial_cell
     &                   /dhwt(i,k,j)
#else
     &                   *dzwr(k)
#endif
          enddo
        enddo

        diff_fb(is:ie,km)=0.
        do i=is,ie
          diff_fb(i,0)        = stf(i,j,n)
          diff_fb(i,kmt(i,j)) = btf(i,j,n)
        enddo
c
c-----------------------------------------------------------------------
c    calculate values for meridional diffusion and advection of tracers
c-----------------------------------------------------------------------
c
        if (j .eq. 2) then
c
c    radiation condition at the southern wall
c    -------------------------------------------
c
         do k=1,km
           do i=is,ie
              ADV_Ty(i,k)=-(phase_vel(i,k,n)+vad(i,k))
     &                  *(t(i,k,j+1,n,tau)-t(i,k,j,n,tau))/dyu(j)
           enddo
         enddo
        end if
c
        if (j .eq. jmt-1) then
c
c    radiation condition at the northern wall
c    -------------------------------------------
c
         do k=1,km
           do i=is,ie
              ADV_Ty(i,k)=-(phase_vel(i,k,n)+vad(i,k))
     &                  *(t(i,k,j,n,tau)-t(i,k,j-1,n,tau))/dyu(j-1)
           enddo
         enddo
        end if
c
c-----------------------------------------------------------------------
c           construct tracer source terms
c           here: force n/s wall to observed values
c-----------------------------------------------------------------------
c
        do k=1,km
         do i=is,ie
          source(i,k,j) = 0.
         enddo
        enddo
c
	if (j == 2 .and. restore_TS_obc_south.and. n<=2 ) then
	    do k=1,km
	      do i=is,ie
                   data=ts_obc_south(i,k,n)
                   source(i,k,j) = source(i,k,j) - 
     &                               obcs*(t(i,k,j,n,tau) - data)
              enddo
	    enddo
	endif
c
	if (j .eq. jmt-1 .and. restore_TS_obc_north.and. n<=2 ) then
	    do k=1,km
	      do i=is,ie
                   data=ts_obc_north(i,k,n)
                   source(i,k,j) = source(i,k,j) - 
     &                               obcn*(t(i,k,j,n,tau) - data)
              enddo
	    enddo
	endif
c
c-----------------------------------------------------------------------
c      calculate the new tracer quantities 
c-----------------------------------------------------------------------
c
        do k=1,km
          do i=is,ie
                t(i,k,j,n,taup1) = t(i,k,j,n,tau) + dt*(
     &                  DIFF_Tx(i,k,j) + DIFF_Tz(i,k)
     &                + DIFF_Ty(i,k,j)+ ADV_Ty(i,k) 
     &                + source(i,k,j) )*tmask(i,k,j)
          enddo
        enddo

      enddo  ! tracer loop

      end subroutine cobc






      subroutine cobc2 (is,ie, js, je, i, n)
      use spflame_module
      implicit none
      integer, intent(in) :: is,ie,js,je,i,n
      integer ii,j,k
      real diff_fb(0:km,js:je)
      real ADV_Tx(km,js:je)
      real uad(km,js:je)
      real phase_vel(km,js:je)
      real var,var2
      real T_i,T_j
      integer ip,jp
#ifdef partial_cell
      real dhte,dhtn,dhwe,dhwn
      dhte(i,k,j)   = min(dht(i+1,k,j),dht(i,k,j))
      dhtn(i,k,j)   = min(dht(i,k,j+1),dht(i,k,j))
      dhwe(i,k,j)   = min(dhwt(i+1,k,j),dhwt(i,k,j))
      dhwn(i,k,j)   = min(dhwt(i,k,j+1),dhwt(i,k,j))

      T_i(i,k,j,n,ip) = t(i+ip,max(1,k-1),j,n,tau) - dhwe(i,k-1,j)
     &            *(t(i+ip,max(1,k-1),j,n,tau) - t(i+ip,k,j,n,tau))
     &               /dhwt(i+ip,k-1,j)
      T_j(i,k,j,n,jp) = t(i,max(1,k-1),j+jp,n,tau) - dhwn(i,k-1,j)
     &            *(t(i,max(1,k-1),j+jp,n,tau) - t(i,k,j+jp,n,tau))
     &               /dhwt(i,k-1,j+jp)
#else
      T_i(i,k,j,n,ip) = t(i+ip,k,j,n,tau)
      T_j(i,k,j,n,jp) = t(i,k,j+jp,n,tau)
#endif
c
c-----------------------------------------------------------------------
c     at open boundaries always laplacian formulation is used (at tau)
c-----------------------------------------------------------------------
c
      real diff_fe_ha(i-1:i,km,js_pe:je_pe)
      real diff_ty,diff_tx,diff_tz

      DIFF_Tx(i,k,j) = (diff_fe_ha(i,  k,j)*tmask(i+1,k,j)
     &                - diff_fe_ha(i-1,k,j)*tmask(i-1,k,j))*cstdxtr(i,j)
#ifdef partial_cell
     &                             /dht(i,k,j)
#endif

      DIFF_Ty(i,k,j) = ahc_north_ha(j)*tmask(i,k,j+1)*
     &                     (t_j(i,k,j,n,1) - t_j(i,k,j,n,0))
     &                         *ah_scale(i,j)
#ifdef partial_cell
     &                         /dht(i,k,j)*dhtn(i,k,j)
#endif
     &                   - ahc_south_ha(j)*tmask(i,k,j-1)*
     &                     (t_j(i,k,j,n,0) - t_j(i,k,j,n,-1))
     &                         *ah_scale(i,j)
#ifdef partial_cell
     &                         /dht(i,k,j)*dhtn(i,k,j-1)
#endif

      DIFF_Tz(i,k,j) = (diff_fb(k-1,j) - diff_fb(k,j))
#ifdef partial_cell
     &               /dht(i,k,j)
#else
     &               *dztr(k)
#endif

      real obce,obcw,data

      obce=1./(86400.*5.); obcw=1./(86400.*5.)
c
      if (i .eq. 2  ) then
        do k=1,km
          do j=js,je
             var = -dxu(i+1)*csu(j)/dt 
             uad(k,j)=(u(i,k,j,1,tau)*dyt(j)+u(i,k,j-1,1,tau)
     &            *dyt(j-1))/(dyt(j)+dyt(j-1))*tmask(i,k,j)
             if (uad(k,j).gt.0) uad(k,j)=0.
	     var2=t(i+2,k,j,n,taum1)-t(i+1,k,j,n,taum1)
	     if (var2.eq.0.) then
		phase_vel(k,j)=var
             else
		phase_vel(k,j)=var*(t(i+1,k,j,n,tau)-t(i+1,k,j,n,taum1))
     &                            /var2*tmask(i,k,j)
                if (phase_vel(k,j).gt. 0.) phase_vel(k,j)=0.
		if (phase_vel(k,j).lt.var) phase_vel(k,j)=var
	     endif
             if (j .eq. 2)     uad(k,j)=0.
             if (j .eq. jmt-1) uad(k,j)=0.
          enddo
        enddo
      endif
c
      if (i .eq. imt-1  ) then
        do k=1,km
          do j=js,je
             var = dxu(i-2)*csu(j)/dt 
             uad(k,j)=(u(i-1,k,j,1,tau)*dyt(j)+u(i-1,k,j-1,1,tau)
     &            *dyt(j-1))/(dyt(j)+dyt(j-1))*tmask(i,k,j)
             if (uad(k,j).lt.0) uad(k,j)=0.
	     var2=t(i-1,k,j,n,taum1)-t(i-2,k,j,n,taum1)
	     if (var2.eq.0.) then
		phase_vel(k,j)=var
	     else
		phase_vel(k,j)=-var*(t(i-1,k,j,n,tau)
     &                -t(i-1,k,j,n,taum1))/var2*tmask(i,k,j)
                if (phase_vel(k,j).lt. 0.) phase_vel(k,j)=0.
		if (phase_vel(k,j).gt.var) phase_vel(k,j)=var
	     endif
             if (j .eq. 2)     uad(k,j)=0.
             if (j .eq. jmt-1) uad(k,j)=0.
          enddo
        enddo
      endif
c
c-----------------------------------------------------------------------
c           diffusive flux across eastern face of "T" cells
c           (use constant horizontal diffusion)
c-----------------------------------------------------------------------
c
      do j=js,je
         do k=1,km
           do ii=i-1,i
             diff_fe_ha(ii,k,j)  =ah_cstdxur(ii,j)*
     &                         (t_i(ii,k,j,n,1) - t_i(ii,k,j,n,0))
     &                         *ah_scale(i,j)
#ifdef partial_cell
     &                         *dhte(i,k,j)
#endif
           enddo
         enddo
      enddo
c
c-----------------------------------------------------------------------
c           diffusive flux across bottom face of "T" cells
c           use constant vertical diffusion
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km-1
	        diff_fb(k,j) = kappa_h*
     &                         (t(i,k,j,n,tau) - t(i,k+1,j,n,tau))
#ifdef partial_cell
     &                   /dhwt(i,k,j)
#else
     &                   *dzwr(k)
#endif
        enddo
      enddo
      diff_fb(km,js:je)=0.
      do j=js,je
         diff_fb(0,j)  = stf(i,j,n)
         diff_fb(kmt(i,j),j) = btf(i,j,n)
      enddo
c
c-----------------------------------------------------------------------
c    radiation condition for lateral advection of tracers
c-----------------------------------------------------------------------
c
      if (i .eq.  2) then
c
c         western wall:
c         -------------
c
        do j=js,je
           do k=1,km
              ADV_Tx(k,j)=-(phase_vel(k,j)+uad(k,j))
     &         *(t(i+1,k,j,n,tau)-t(i,k,j,n,tau))*dxur(i)*csur(j)
           enddo
        enddo
      endif
c
      if (i .eq. imt-1  ) then
c
c         eastern wall:
c         -------------
c
        do j=js,je
           do k=1,km
              ADV_Tx(k,j)=-(phase_vel(k,j)+uad(k,j))
     &         *(t(i,k,j,n,tau)-t(i-1,k,j,n,tau))*dxur(i-1)*csur(j)
           enddo
        enddo
      endif
c
c-----------------------------------------------------------------------
c           construct tracer source terms
c           here: force w/e wall to observed values
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km
          source(i,k,j) = 0.
        enddo
      enddo
c
      if (i .eq. 2 .and. restore_TS_obc_west .and. n<=2) then
	    do j=js,je
	      do k=1,km
                   data=ts_obc_west(j,k,n)
                   source(i,k,j) = source(i,k,j) - 
     &                               obcw*(t(i,k,j,n,tau) - data)
              enddo
	    enddo
      endif
c
      if (i .eq.  imt-1 .and. restore_TS_obc_east .and. n<=2) then
	    do j=js,je
	      do k=1,km
                   data=ts_obc_east(j,k,n)
                   source(i,k,j) = source(i,k,j) - 
     &                               obce*(t(i,k,j,n,tau) - data)
              enddo
	    enddo
      endif
c
c-----------------------------------------------------------------------
c     calculate the new tracer quantities 
c-----------------------------------------------------------------------
c
      do j=js,je
         do k=1,km
                t(i,k,j,n,taup1) = t(i,k,j,n,tau) + dt*(
     &                  DIFF_Tx(i,k,j) 
     &                + DIFF_Tz(i,k,j) 
     &                + DIFF_Ty(i,k,j) 
     &                + ADV_Tx(k,j) 
     &                + source(i,k,j)  )*tmask(i,k,j)
         enddo
      enddo
c
      end subroutine cobc2
