#include "options.inc"


      subroutine tracer ()
c
c=======================================================================
c     compute tracers at "tau+1"
c     original author:   r.c.pacanowski       e-mail rcp@gfdl.gov
c     SPFLAME version:   c.eden
c=======================================================================
c
      use spflame_module
      use bbl_module
      implicit none
c-----------------------------------------------------------------------
c     local arrays
c-----------------------------------------------------------------------
      integer is,ie,js,je,is1,ie1,js1,je1,ip,jp
      integer i,j,k,n
      real :: tend(is_pe:ie_pe), fxa
      real T_i,T_j
#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,taum1) - dhwe(i,k-1,j)
     &            *(t(i+ip,max(1,k-1),j,n,taum1) - t(i+ip,k,j,n,taum1))
     &               /dhwt(i+ip,k-1,j)
      T_j(i,k,j,n,jp) = t(i,max(1,k-1),j+jp,n,taum1) - dhwn(i,k-1,j)
     &            *(t(i,max(1,k-1),j+jp,n,taum1) - t(i,k,j+jp,n,taum1))
     &               /dhwt(i,k-1,j+jp)
#else
      T_i(i,k,j,n,ip) = t(i+ip,k,j,n,taum1)
      T_j(i,k,j,n,jp) = t(i,k,j+jp,n,taum1)
#endif
c
c-----------------------------------------------------------------------
c     advective terms
c-----------------------------------------------------------------------
c
      real adv_fe(is_pe-1:ie_pe,km,js_pe:je_pe)
      real adv_fn(is_pe:ie_pe,km,js_pe-1:je_pe)
      real adv_fb(is_pe:ie_pe,0:km,js_pe:je_pe)
      real adv_tx
      ADV_Tx(i,k,j) = (adv_fe(i,k,j) - adv_fe(i-1,k,j))*cstdxt2r(i,j)
#ifdef partial_cell
     &                             /dht(i,k,j)
#endif
      real adv_ty
      ADV_Ty(i,k,j) = (adv_fn(i,k,j) - adv_fn(i,k,j-1))*cstdyt2r(j)
# ifdef partial_cell
     &                             /dht(i,k,j)
# endif
      real adv_tz
      ADV_Tz(i,k,j) = (adv_fb(i,k-1,j) - adv_fb(i,k,j))
#ifdef partial_cell
     &               *.5/dht(i,k,j)
#else
     &               *dzt2r(k)
#endif
c
c-----------------------------------------------------------------------
c     diffusive terms
c-----------------------------------------------------------------------
c
      real :: diff_fe(is_pe-1:ie_pe,km,js_pe:je_pe) 
      real :: diff_fn(is_pe:ie_pe,km,js_pe-1:je_pe) 
      real :: diff_fb(is_pe:ie_pe,0:km,js_pe:je_pe)
      real diff_tx
      DIFF_Tx(i,k,j) = (diff_fe(i,  k,j)*tmask(i+1,k,j)
     &            - diff_fe(i-1,k,j)*tmask(i-1,k,j))*cstdxtr(i,j)
#ifdef partial_cell
     &                             /dht(i,k,j)
#endif
      real diff_ty
      DIFF_Ty(i,k,j) = (diff_fn(i,k,j  )*tmask(i,k,j+1)
     &             - diff_fn(i,k,j-1)*tmask(i,k,j-1))*cstdytr(j)
#ifdef partial_cell
     &                             /dht(i,k,j)
#endif
      real DIFF_Tz
      DIFF_Tz(i,k,j) = (diff_fb(i,k-1,j) - diff_fb(i,k,j))
     &                 *(1.-aidif_tracer)
#ifdef partial_cell
     &               /dht(i,k,j)
#else
     &               *dztr(k)
#endif
c
c     part of vertical diffusive flux due to isopyncal mixing
c     which is solved explicitly, other part is contained
c     in the vertical diffusivity diff_cbt and solved implicitly
c
      real, allocatable :: diff_fbiso(:,:,:)
      real DIFF_Tz_iso
      DIFF_Tz_iso(i,k,j) = (diff_fbiso(i,k-1,j) - diff_fbiso(i,k,j))
#ifdef partial_cell
     &               /dht(i,k,j)
#else
     &               *dztr(k)
#endif
c

#ifdef detailed_timing
      call tic('tracer main')
#endif
      if (enable_diffusion_isoneutral.or. 
     &    enable_diffusion_isopycnic) then
       allocate( diff_fbiso(is_pe:ie_pe,0:km,js_pe:je_pe) )
       diff_fbiso=0.
      endif

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

#ifdef detailed_timing
      call tic('tracer obc')
#endif

c
c     care about TS values for open boundaries, read also psi
c
      if (.not.enable_simple_obc) call obc_read()
c
      js1=js
      if (my_blk_j == 1 .and. enable_obc_south) then
       call cobc(2); js1=3
      endif

      je1=je       
      if (my_blk_j == n_pes_j .and. enable_obc_north) then
       je1=jmt-2
      endif

#ifdef detailed_timing
      call toc('tracer obc')
#endif

      do n=1,nt
c
c-----------------------------------------------------------------------
c       calculate 2*advective flux across 
c       east/north/bottom faces of "T" cells.
c-----------------------------------------------------------------------
c
#ifdef detailed_timing
        call tic('tracer adv')
#endif

        if (enable_quicker_advection) then
c
c        quicker scheme
c
         call adv_east_flux_quicker(n,adv_fe)
         call adv_north_flux_quicker(n,adv_fn)
         call adv_bottom_flux_quicker(n,adv_fb)

        elseif (enable_upstream_advection) then
c
c        upstream advection scheme
c
         call adv_flux_upstream(n,adv_fe,adv_fn,adv_fb)

        elseif (enable_4th_advection) then
c
c        4 th order advection scheme
c
         call adv_flux_4th(n,adv_fe,adv_fn,adv_fb)

        elseif (enable_fct_advection) then
c
c         FCT advection scheme
c
         call adv_flux_fct(n,adv_fe,adv_fn,adv_fb)

        else
c
c        centered differences
c
         do j=js,je
#ifdef vector_host
          adv_fe(is-1:ie,1:km,j) = adv_vet(is-1:ie,1:km,j)*
     &           (t(is-1:ie,1:km,j,n,tau)+t(is:ie+1,1:km,j,n,tau))
#else
          do k=1,km
            do i=is-1,ie
	      adv_fe(i,k,j) = adv_vet(i,k,j)*(t(i,  k,j,n,tau) + 
     &                                        t(i+1,k,j,n,tau))
            enddo
          enddo
#endif
         enddo

         do j=js-1,je
#ifdef vector_host
	  adv_fn(is:ie,1:km,j) = adv_vnt(is:ie,1:km,j)*
     &      (t(is:ie,1:km,j,n,tau) + t(is:ie,1:km,j+1,n,tau))
#else
          do k=1,km
           do i=is,ie
	    adv_fn(i,k,j) = adv_vnt(i,k,j)*(t(i,k,j  ,n,tau) + 
     &                                      t(i,k,j+1,n,tau))
           enddo
          enddo
#endif
         enddo

         do j=js,je
#ifdef vector_host
	  adv_fb(is:ie,1:km-1,j) = adv_vbt(is:ie,1:km-1,j)*
     &          (t(is:ie,1:km-1,j,n,tau) + t(is:ie,2:km,j,n,tau))
#else
          do k=1,km-1
            do i=is,ie
	      adv_fb(i,k,j)  = adv_vbt(i,k,j)*(t(i,k,  j,n,tau) +
     &                                         t(i,k+1,j,n,tau))
            enddo
          enddo
#endif
         enddo
         
        endif

        if (enable_diffusion_isopycnic .and. n>2) then
          call isopycnic_add_gm(n,adv_fe,adv_fn,adv_fb)
        endif

#ifdef detailed_timing
        call toc('tracer adv')
#endif
c
c-----------------------------------------------------------------------
c       calculate diffusive flux across eastern and northern faces 
c       of "T" cells due to various parameterizations for diffusion.
c-----------------------------------------------------------------------
c
#ifdef detailed_timing
        call tic('tracer diff')
#endif
        if (enable_diffusion_harmonic) then
c
c        diffusive flux on eastern face of "T" cells
c
         do j=js,je
          do k=1,km
            do i=is-1,ie              
              diff_fe(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
         enddo
c
c        north face 
c
          do j=js-1,je
            fxa=ah*csu_dyur(j)
            do k=1,km
              do i=is,ie
                diff_fn(i,k,j) =fxa*
     &            (T_j(i,k,j,n,1)-T_j(i,k,j,n,0))*ah_scale(i,j)
#ifdef partial_cell
     &                           *dhtn(i,k,j)
#endif
              enddo
            enddo
          enddo
c 
        else
         diff_fe=0.
         diff_fn=0.
        endif
c
c-----------------------------------------------------------------------
c       calculate diffusive flux across bottom face of "T" cells
c-----------------------------------------------------------------------
c
        do j=js,je
          do k=1,km-1
            do i=is,ie
              diff_fb(i,k,j) = diff_cbt(i,k,j)*
     &                         (t(i,k,j,n,taum1) - t(i,k+1,j,n,taum1))
#ifdef partial_cell
     &                   /dhwt(i,k,j)
#else
     &                   *dzwr(k)
#endif
            enddo
          enddo
        enddo

#ifdef detailed_timing
        call toc('tracer diff')
#endif
c
c-----------------------------------------------------------------------
c       set boundary conditions on "t" cells for vertical diffusion
c       and vertical advection
c-----------------------------------------------------------------------
c
        diff_fb(is:ie,km,js:je)=0.
        do j=js,je
          do i=is,ie
            diff_fb(i,0,j)  = stf(i,j,n)
            diff_fb(i,kmt(i,j),j) = btf(i,j,n)
	    adv_fb(i,0,j)   = adv_vbt(i,0,j)*(t(i,1,j,n,tau) +
     &                                       t(i,1,j,n,tau)) 
	    adv_fb(i,km,j)  = 0.
          enddo
        enddo
c
c-----------------------------------------------------------------------
c       add diffusive fluxes along isopycnals/neutral surfaces
c-----------------------------------------------------------------------
c
        if (enable_diffusion_isoneutral) then
          call isoneutral_flux(n,diff_fe,diff_fn,diff_fbiso)
        endif

        if (enable_diffusion_isopycnic) then
         call isopycnic_flux(n,diff_fe,diff_fn,diff_fbiso)
        endif
c
c-----------------------------------------------------------------------
c       add biharmonic diffusive fluxes
c-----------------------------------------------------------------------
c
        if (enable_diffusion_biharmonic) then
          call delsqt(n,diff_fe,diff_fn)
        endif
c
c-----------------------------------------------------------------------
c       limit advective fluxes for passive tracer
c       (works for positive tracers only)
c-----------------------------------------------------------------------
c
        if (n>2 .and. enable_flux_delimiter) then
         call delimit_adv_flux(n,adv_fe,adv_fn,adv_fb)
         if (enable_diffusion_isoneutral.or.
     &    enable_diffusion_isopycnic) then
          call delimit_all_fluxes(n,adv_fe,adv_fn,adv_fb,
     &                                     diff_fe,diff_fn,diff_fbiso)
         endif
        endif
c
c-----------------------------------------------------------------------
c       set source term for "T" cells
c-----------------------------------------------------------------------
c
        do j=js,je
          do k=1,km
            do i=is,ie
	      source(i,k,j) = 0.
            enddo
          enddo
        enddo

        if (enable_sponge) then
         if (enable_simple_spg) then
          call sponge_template(n)
         else
          call sponge_read(n)
         endif
        endif

        if (enable_shortwave_sbc .and. n==1 ) then
c
c=======================================================================
c       incorporate short wave penetration via the "source" 
c       term. note that the divergence of shortwave for the first
c       level "divpen(1)" is compensating for the effect of having
c       the shortwave component already included in the total
c       surface tracer flux "stf(i,j,1)"
c
c       incorporating the penetrative shortwave radiative flux into
c       the vertical diffusive flux term directly is not correct when
c       vertical diffusion is solved implicitly. The tri-diagonal
c       solver takes diffusion coefficients at tops and bottoms of
c       cells but not the fluxes.
c=======================================================================
c
         do j=js,je
	  do k=1,km
	    do i=is,ie
              source(i,k,j) = source(i,k,j)
     &                       + qsol(i,j)*divpen_shortwave(k) 
            enddo
          enddo
         enddo

        endif
c
c-----------------------------------------------------------------------
c       solve for "tau+1" tracer at center of "T" cells using statement
c       functions to represent each component of the calculation 
c       choose start and end indices for "full" calculation
c       or calculate zonal open boundary conditions
c-----------------------------------------------------------------------
c
#ifdef detailed_timing
        call tic('tracer obc')
#endif

        is1=is
        if ( my_blk_i == 1 .and. enable_obc_west) then
           call cobc2 (is,ie,js,je,2,n)
           is1=max(is,3)
        endif

        ie1=ie
        if ( my_blk_i == n_pes_i .and. enable_obc_east) then
          call cobc2 (is,ie,js,je,imt-1,n)
          ie1=min(ie,imt-2)
        endif

#ifdef detailed_timing
        call toc('tracer obc')
c
        call tic('tracer tendency')
#endif

        do j=js1,je1
         do k=1,km
          if (enable_bbl.and.enable_bbl_advection) then
           do i=is1,ie1
	    tend(i) = DIFF_Tz(i,k,j)+ source(i,k,j) +
     &         (- ADV_Tx(i,k,j) - ADV_Ty(i,k,j) - ADV_Tz(i,k,j))
     &                 *(1.0-sigma_weight(i,k,j))
     & 	           + DIFF_Tx(i,k,j)+ DIFF_Ty(i,k,j)
           enddo
          else
           do i=is1,ie1
	    tend(i) = DIFF_Tz(i,k,j)+ source(i,k,j)
     &             - ADV_Tx(i,k,j)-ADV_Ty(i,k,j)-ADV_Tz(i,k,j)
     & 	           + DIFF_Tx(i,k,j)+ DIFF_Ty(i,k,j)
           enddo
          endif
          if (enable_diffusion_isoneutral.or.
     &        enable_diffusion_isopycnic   ) then
           do i=is1,ie1
            tend(i)=tend(i)+DIFF_Tz_iso(i,k,j)
           enddo
          endif
          do i=is1,ie1
	   t(i,k,j,n,taup1) = t(i,k,j,n,taum1) + 
     &                 c2dt*tend(i)*tmask(i,k,j)
          enddo
         enddo
        enddo

#ifdef detailed_timing
        call toc('tracer tendency')
#endif
c
c-----------------------------------------------------------------------
c       add in dT/dt component due to explicit bolus advection
c       (for isoneutral this part is absorbed in diffusive (skew) flux.
c-----------------------------------------------------------------------
c
        if (enable_diffusion_isopycnic) then
         if (n>2) then
c  advection velocities already added to advective fluxes for passive tracers
         else
          call isopycnic_add_div_gm(n,is1,ie1,js1,je1)
         endif
        endif
c
c-----------------------------------------------------------------------
c       add in dT/dt component due to Bottom boundary layer
c       of Beckmann + Doescher.
c-----------------------------------------------------------------------
c
        if (enable_bbl) then
         call add_bbl(n,is1,ie1,js1,je1,adv_fe,adv_fb)
        endif

#ifdef detailed_timing
        call tic('tracer implicit')
#endif
c
c-----------------------------------------------------------------------
c       add in dT/dt component due to implicit vertical diffusion
c       (open boundaries use always explicit vertical diffusion with kappa_h)
c-----------------------------------------------------------------------
c
        if (enable_implicit_vert_diff)
     &   call implicit_vert_diff(n,is1,ie1,js1,je1)

#ifdef detailed_timing
        call toc('tracer implicit')
#endif
c
c-----------------------------------------------------------------------
c       Diagnostics which depends on the advective/ diffusive
c       fluxes for tracer n.
c-----------------------------------------------------------------------
c
        if (enable_heat_tr .and. snapshot_time_step) 
     &         call diag_heat_tr(n,adv_fn,diff_fn)

      enddo  ! n loop

      if (my_blk_j == n_pes_j.and. enable_obc_north) call cobc(jmt-1)
c
c-----------------------------------------------------------------------
c     explicit convection: adjust column if gravitationally unstable
c-----------------------------------------------------------------------
c
#ifdef detailed_timing
      call tic('tracer convect')
#endif
      if (enable_expl_convection) call convect (taup1)
#ifdef detailed_timing
      call toc('tracer convect')
#endif
c
c-----------------------------------------------------------------------
c     Kraus-Turner-Mixing
c-----------------------------------------------------------------------
c
#ifdef detailed_timing
      call tic('tracer ktmix')
#endif
      if (enable_ktmix) call Ktmix
#ifdef detailed_timing
      call toc('tracer ktmix')
#endif
c
c-----------------------------------------------------------------------
c      Add tendencies due to "unphysical" sources :)
c-----------------------------------------------------------------------
c
       call tic('passive tracer')
       call passive_tracer_add()
       call toc('passive tracer')
c
c-----------------------------------------------------------------------
c      set boundary conditions at each PE
c-----------------------------------------------------------------------
c
#ifdef detailed_timing
      call tic('tracer mpp')
#endif
      do n=1,nt
c        call set_cyclic(t(:,:,:,n,taup1),km,2)
        call set_cyclic(t(is_pe-2,1,js_pe-2,n,taup1),km,2)
        if (enable_quicker_advection.or.enable_4th_advection
     &                             .or. enable_fct_advection) 
     &      call set_cyclic_4th_order(t(:,:,:,n,taup1),km,2)
        if (my_blk_j == 1 .and. enable_obc_south) 
     &        t(:,:,1,n,taup1) = t(:,:,2,n,taup1)
        if (my_blk_j == n_pes_j .and. enable_obc_north)
     &        t(:,:,jmt,n,taup1) = t(:,:,jmt-1,n,taup1)
        if (my_blk_i == 1 .and. enable_obc_west) 
     &        t(1,:,:,n,taup1) = t(2,:,:,n,taup1)
        if (my_blk_i == n_pes_i .and. enable_obc_east) 
     &        t(imt,:,:,n,taup1) = t(imt-1,:,:,n,taup1)
c        call border_exchg(t(:,:,:,n,taup1),km,2)
        call border_exchg(t(is_pe-2,1,js_pe-2,n,taup1),km,2)
      enddo
#ifdef detailed_timing
      call toc('tracer mpp')
#endif
c
      if (enable_diffusion_isoneutral.or.
     &    enable_diffusion_isopycnic) 
     &    deallocate( diff_fbiso )

#ifdef detailed_timing
      call toc('tracer main')
#endif
      end subroutine tracer





      subroutine delsqt (n,diff_fe,diff_fn)
      use spflame_module
      implicit none
      integer is,ie,js,je,i,j,k,n,ip,jp,ism1,iep1,jsm1,jep1
      real :: diff_fe(is_pe-1:ie_pe,km,js_pe:je_pe) 
      real :: diff_fn(is_pe:ie_pe,km,js_pe-1:je_pe) 
      real T_i,T_j
#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,taum1) - dhwe(i,k-1,j)
     &            *(t(i+ip,max(1,k-1),j,n,taum1) - t(i+ip,k,j,n,taum1))
     &               /dhwt(i+ip,k-1,j)
      T_j(i,k,j,n,jp) = t(i,max(1,k-1),j+jp,n,taum1) - dhwn(i,k-1,j)
     &            *(t(i,max(1,k-1),j+jp,n,taum1) - t(i,k,j+jp,n,taum1))
     &               /dhwt(i,k-1,j+jp)
#else
      T_i(i,k,j,n,ip) = t(i+ip,k,j,n,taum1)
      T_j(i,k,j,n,jp) = t(i,k,j+jp,n,taum1)
#endif

      real ::  diff_fe_biha(is_pe-2:ie_pe+1,km,js_pe-1:je_pe+1) 
      real ::  diff_fn_biha(is_pe-1:ie_pe+1,km,js_pe-2:je_pe+1) 
      real ::  del2(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1) 
      real diff_tx,diff_ty
      DIFF_Tx(i,k,j) = (diff_fe_biha(i,  k,j)*tmask(i+1,k,j)
     &            - diff_fe_biha(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) = (diff_fn_biha(i,k,j  )*tmask(i,k,j+1)
     &                - diff_fn_biha(i,k,j-1)*tmask(i,k,j-1))*cstdytr(j)
#ifdef partial_cell
     &                             /dht(i,k,j)
#endif

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
      ism1=max(2,is-1); iep1=min(imt-1,ie+1)
      jsm1=max(2,js-1); jep1=min(jmt-1,je+1)
c
c=======================================================================
c     compute del**2 of prognostic variables on MW rows "js" ... "je"
c     and gradients of del**2 quantities for variable "n" in the
c     tracer equation
c=======================================================================

c
c-----------------------------------------------------------------------
c       diffusive flux across eastern face of "T" cell
c       diffusive flux across northern face of "T" cell
c-----------------------------------------------------------------------
c
c       fe ism1-1:iep1 jsm1:jep1
        do j=js-1,je+1
          do k=1,km
            do i=ism1-1,iep1
              diff_fe_biha(i,k,j) = diff_cet*cstr(j)*dxur(i)*
     &                   (t_i(i,k,j,n,1) - t_i(i,k,j,n,0))
     &                        *ahbi_scale(i,j)
#ifdef partial_cell
     &                        *dhte(i,k,j)
#endif
            enddo
          enddo
        enddo

c       fn jsm1:iep1 jsm1-1:jep1
        do j=jsm1-1,jep1
          do k=1,km
            do i=is-1,ie+1
              diff_fn_biha(i,k,j) = diff_cnt*csu(j)*dyur(j)* 
     &                     (t_j(i,k,j,n,1) - t_j(i,k,j,n,0))
     &                         *ahbi_scale(i,j)
#ifdef partial_cell
     &                         *dhtn(i,k,j)
#endif
            enddo
          enddo
        enddo
c
c       compute -ah*del**2 of tracer
c
c       del2 is-1:ie+1 js-1:je+1
        do j=jsm1,jep1
          do k=1,km
            do i=ism1,iep1
	      del2(i,k,j) = -DIFF_Tx(i,k,j) - DIFF_Ty(i,k,j)
            enddo
          enddo
          if (my_blk_i == 1 .and. enable_obc_west) then
             del2(1,:,:) = del2(2,:,:)
          elseif (my_blk_i == 1)   then
             del2(1,:,:) = 0.
          endif
          if (my_blk_i == n_pes_i .and. enable_obc_east) then
            del2(imt,:,:) = del2(imt-1,:,:)
          elseif (my_blk_i == n_pes_i) then
            del2(imt,:,:) = 0.
          endif
        enddo

        call set_cyclic(del2(:,:,:),km,1)
c
c       set -del**2 = 0 on southern wall
c
        if (my_blk_j == 1)  del2(is-1:ie+1,:,1) = 0.
c
c       set -del**2 = 0 on northern wall
c
        if (my_blk_j == n_pes_j)  del2(is-1:ie+1,:,jmt) = 0.

c-----------------------------------------------------------------------
c       diffusive flux across eastern face of "T" cell
c       diffusive flux across northern face of "T" cell
c-----------------------------------------------------------------------
c      fe  is-1:ie js:je
        do j=js,je
          do k=1,km
            do i=is-1,ie
              diff_fe(i,k,j) = diff_fe(i,k,j)+
     &                         diff_cet*cstr(j)*dxur(i)*
     &                         (del2(i+1,k,j)-del2(i,k,j))        
     &                        *ahbi_scale(i,j)
#ifdef partial_cell
     &                        *dhte(i,k,j)
#endif
            enddo
          enddo
        enddo
c      fn  is:ie js-1:je
        do j=js-1,je
          do k=1,km
            do i=is,ie
              diff_fn(i,k,j) = diff_fn(i,k,j) +
     &                         diff_cnt*csu(j)*dyur(j)*
     &                        (del2(i,k,j+1) - del2(i,k,j))
     &                         *ahbi_scale(i,j)
#ifdef partial_cell
     &                         *dhtn(i,k,j)
#endif
            enddo
          enddo
        enddo
	end subroutine delsqt

