#include "options.inc"


      subroutine clinic 
c
c=======================================================================
c     compute internal mode velocity components for "tau+1"
c     original author:   r.c.pacanowski       e-mail rcp@gfdl.gov
c     SPFLAME version:   c.eden
c=======================================================================
c
      use spflame_module
      use freesurf_module
      implicit none
      integer :: is,ie,js,je,i,j,k,n
      integer :: is1,ie1,js1,je1
      real    :: baru(is_pe:ie_pe),fxa,fxb
      real    :: grad_p(is_pe:ie_pe,km,js_pe:je_pe,2) 
#ifdef partial_cell
      real dhue,dhun
      dhue(i,k,j)   = min(dhu(i+1,k,j),dhu(i,k,j))
      dhun(i,k,j)   = min(dhu(i,k,j+1),dhu(i,k,j))
#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_ux,adv_uy,adv_uz,adv_metric

      ADV_Ux(i,k,j)   = (adv_fe(i,k,j) - adv_fe(i-1,k,j))*csudxu2r(i,j)
#ifdef partial_cell
     &                             /dhu(i,k,j)
#endif

      ADV_Uy(i,k,j) =(adv_fn(i,k,j) - adv_fn(i,k,j-1))*csudyu2r(j)
#ifdef partial_cell
     &                             /dhu(i,k,j)
#endif

      ADV_Uz(i,k,j) = (adv_fb(i,k-1,j) - adv_fb(i,k,j))
#ifdef partial_cell
     &                *0.5/dhu(i,k,j)
#else
     &                *dzt2r(k)
#endif
      ADV_metric(i,k,j,n)=advmet(j,n)*u(i,k,j,1,tau)*u(i,k,j,3-n,tau)
c
c-----------------------------------------------------------------------
c     viscosity
c-----------------------------------------------------------------------
c
      real :: DIFF_ux,diff_fe(is_pe-1:ie_pe,  km,js_pe  :je_pe) 
      real :: DIFF_Uy,diff_fn(is_pe  :ie_pe,  km,js_pe-1:je_pe) 
      real :: DIFF_uz,diff_fb(is_pe  :ie_pe,0:km,js_pe  :je_pe)
      real :: diff_metric(is_pe:ie_pe,km,js_pe:je_pe)
      real :: del2(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1,2)

      DIFF_Ux(i,k,j) = (diff_fe(i,k,j) - diff_fe(i-1,k,j))
     &                      *rho0csudxur(i,j)
#ifdef partial_cell
     &                            /dhu(i,k,j)
#endif

      DIFF_Uy(i,k,j) = (diff_fn(i,k,j)-diff_fn(i,k,j-1))
     &            *rho0csudyur(j)
#if defined partial_cell
     &            /dhu(i,k,j)
#endif

      DIFF_Uz(i,k,j) = (diff_fb(i,k-1,j) - diff_fb(i,k,j))
     &                 *(1.-aidif_momentum)
#if defined partial_cell
     &                 /(rho0*dhu(i,k,j))
#else
     &                 *rho0dztr(k)
#endif


#ifdef detailed_timing
      call tic('clinic main')
#endif

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     calculate hydrostatic pressure gradient
c-----------------------------------------------------------------------
c
      call pressure_gradient(grad_p)
c
c-----------------------------------------------------------------------
c     calculate del**2 u for biharmonic friction first to insure
c     correct treatment of biharmonic metric terms.
c-----------------------------------------------------------------------
c
#ifdef detailed_timing
      call tic('clinic biha')
#endif

      if (enable_friction_biharmonic) call delsqu(del2)

#ifdef detailed_timing
      call toc('clinic biha')
#endif
c
c-----------------------------------------------------------------------
c     solve for one component of velocity at a time
c     n = 1 => zonal component
c     n = 2 => meridional component
c-----------------------------------------------------------------------
c
      do n=1,2
c
c-----------------------------------------------------------------------
c       calculate 2*advective flux (for speed) across east face of
c       "u" cells.
c-----------------------------------------------------------------------
c
#ifdef detailed_timing
      call tic('clinic adv')
#endif
        do j=js,je
          do k=1,km
            do i=is-1,ie
	      adv_fe(i,k,j) = adv_veu(i,k,j)*(u(i,  k,j,n,tau) + 
     &                                        u(i+1,k,j,n,tau))
            enddo
          enddo
        enddo
c
c-----------------------------------------------------------------------
c       2*advective flux across northern face of "u" cells 
c-----------------------------------------------------------------------
c
        do j=js-1,je
         do k=1,km
          do i=is,ie
	    adv_fn(i,k,j) = adv_vnu(i,k,j)*(u(i,k,j  ,n,tau) + 
     &                                      u(i,k,j+1,n,tau))
          enddo
         enddo
        enddo
#ifdef detailed_timing
        call toc('clinic adv')
#endif
c
c-----------------------------------------------------------------------
c       diffusive flux across east face of "u" cell
c       diffusive flux across north face of "u" cell
c-----------------------------------------------------------------------

#ifdef detailed_timing
        call tic('clinic diff')
#endif
        if (enable_friction_harmonic) then
c
c       calculate diffusive flux on northern faces of "u" cells 
c       with laplacian friction
c
         do j=js-1,je
          fxa=am*cst(j+1)*dytr(j+1)
          do k=1,km
           do i=is,ie
              diff_fn(i,k,j) = fxa*am_scale(i,j+1)*
     &                          (u(i,k,j+1,n,taum1) - u(i,k,j,n,taum1))
#if defined partial_cell
     &                           *dhun(i,k,j)          
#endif
           enddo
          enddo
         enddo
c
c       calculate diffusive flux on eastern faces of "u" cells 
c       with laplacian friction
c
         do j=js,je
          do k=1,km
            do i=is-1,ie
              diff_fe(i,k,j) = am_csudxtr(i,j)*
     &                          (u(i+1,k,j,n,taum1) - u(i,k,j,n,taum1))
     &                           *am_scale(i,j)
#if defined partial_cell
     &                           *dhue(i,k,j)          
#endif
            enddo
          enddo
         enddo
c
c        calculate metric term
c
         do j=js,je
          do k=1,km
           do i=is,ie
            diff_metric(i,k,j)=am3_ha(j)*u(i,k,j,n,taum1)*am_scale(i,j)
     &                       + am4_ha(j,n)*dxmetr(i)*am_scale(i,j)
     &                    *(u(i+1,k,j,3-n,taum1)-u(i-1,k,j,3-n,taum1))
#ifdef partial_cell
     &                        +pc_sink_ha(i,k,j)*u(i,k,j,n,taum1)
#endif
            enddo
          enddo
         enddo

        else
         diff_fn=0.
         diff_fe=0.
         diff_metric=0.
        endif

#ifdef detailed_timing
        call toc('clinic diff')
#endif
c
c-----------------------------------------------------------------------
c       calculate diffusive flux on eastern and northern faces of
c      "u" cells with biharmonic friction
c-----------------------------------------------------------------------
c
#ifdef detailed_timing
        call tic('clinic biha')
#endif
        if (enable_friction_biharmonic) then

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)+
     &                         visc_ceu*csur(j)*dxtr(i+1)*
     &                        (del2(i+1,k,j,n)-del2(i,k,j,n))
     &                        *ambi_scale(i+1,j)
#if defined partial_cell
     &                        *dhue(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)+
     &                         visc_cnu*cst(j+1)*dytr(j+1)*
     &                        (del2(i,k,j+1,n) - del2(i,k,j,n))
     &                        *ambi_scale(i,j+1)
#if defined partial_cell
     &                        *dhun(i,k,j)          
#endif
            enddo
          enddo
         enddo
c
c        calculate metric term
c
         do j=js,je
          do k=1,km
           do i=is,ie
            diff_metric(i,k,j) = diff_metric(i,k,j)+
     &         am3_biha(j)*del2(i,k,j,n)*ambi_scale(i,j)
     & +am4_biha(j,n)*(del2(i+1,k,j,3-n) - del2(i-1,k,j,3-n))*dxmetr(i)
     &                            *ambi_scale(i,j)
#if defined partial_cell
     &                        +pc_sink_biha(i,k,j)*del2(i,k,j,n)
#endif
           enddo
          enddo
         enddo

        endif

#ifdef detailed_timing
        call toc('clinic biha')
#endif
c
c-----------------------------------------------------------------------
c       calculate 2*advective flux (for speed) on bottom face of
c       "u" cell. also diffusive flux on bottom face of "u" cell
c-----------------------------------------------------------------------
c
        do j=js,je
          do k=1,km-1
            do i=is,ie
	      adv_fb(i,k,j) = adv_vbu(i,k,j)*(u(i,k,  j,n,tau) +
     &                                        u(i,k+1,j,n,tau))
              diff_fb(i,k,j) = visc_cbu(i,k,j)*
     &                         (u(i,k,j,n,taum1) - u(i,k+1,j,n,taum1))
# ifdef partial_cell
     &                      /min(dhwt(i,k,j),   dhwt(i+1,k,j)
     &,                           dhwt(i,k,j+1), dhwt(i+1,k,j+1))
# else
     &                        *dzwr(k)
# endif

            enddo
          enddo
        enddo
        diff_fb(is:ie,km,js:je)=0.
c
c-----------------------------------------------------------------------
c       set surface and bottom vert b.c. on "u" cells for mixing
c       set surface and bottom vert b.c. on "u" cells for advection
c-----------------------------------------------------------------------
c
        do j=js,je
          do i=is,ie
            diff_fb(i,0,j)  = smf(i,j,n)
            diff_fb(i,kmu(i,j),j) = bmf(i,j,n)
	    adv_fb(i,0,j)   = adv_vbu(i,0,j)*(u(i,1,j,n,tau) +
     &                                        u(i,1,j,n,tau)) 
	    adv_fb(i,km,j)  = 0.
          enddo
        enddo

c
c-----------------------------------------------------------------------
c       solve for the internal mode part of du/dt at center of 
c       "u" cells by neglecting the surface pressure gradients. use
c       statement functions to represent each component of the 
c       calculation. 
c       Open boundaries need special treatment.
c-----------------------------------------------------------------------
c
#ifdef detailed_timing
        call tic('clinic obc')
#endif

        if (enable_obc_south.and. my_blk_j==1) then
         js1=3; j=2
         do k=1,km
            do i=is,ie
	      u(i,k,j,n,taup1) = (DIFF_Uz(i,k,j) - grad_p(i,k,j,n) 
     &           + DIFF_Ux(i,k,j) + DIFF_Uy(i,k,j) +diff_metric(i,k,j)
     &           + cori(i,j,n)*u(i,k,j,3-n,tau) 
     &              )*umask(i,k,j)
            enddo
         enddo
	 u(:,:,1,n,taup1) = u(:,:,2,n,taup1)
        else
         js1=js
        endif

        if (enable_obc_north.and. my_blk_j==n_pes_j) then
         je1=jmt-3; j=jmt-2
         do k=1,km
            do i=is,ie
	      u(i,k,j,n,taup1) = (DIFF_Uz(i,k,j) - grad_p(i,k,j,n) 
     &           + DIFF_Ux(i,k,j) + DIFF_Uy(i,k,j) +diff_metric(i,k,j)
     &           + cori(i,j,n)*u(i,k,j,3-n,tau) 
     &              )*umask(i,k,j)
            enddo
         enddo
	 u(:,:,jmt  ,n,taup1) = u(:,:,jmt-2,n,taup1)
	 u(:,:,jmt-1,n,taup1) = u(:,:,jmt-2,n,taup1)
        else
         je1=je
        endif

        if (enable_obc_west.and. my_blk_i==1) then
         is1=3; i=2
         do j=js,je
           do k=1,km
	      u(i,k,j,n,taup1) = (DIFF_Uz(i,k,j)- grad_p(i,k,j,n) 
     &           + DIFF_Ux(i,k,j) + DIFF_Uy(i,k,j) +diff_metric(i,k,j)
     &           + cori(i,j,n)*u(i,k,j,3-n,tau) 
     &              )*umask(i,k,j)
            enddo
         enddo
	 u(i-1,:,:,n,taup1) = u(i,:,:,n,taup1)
        else
         is1=is
        endif

        if (enable_obc_east.and. my_blk_i==n_pes_i) then
         ie1=imt-3; i=imt-2
         do j=js,je
           do k=1,km
	      u(i,k,j,n,taup1) = (DIFF_Uz(i,k,j) - grad_p(i,k,j,n) 
     &           + DIFF_Ux(i,k,j) + DIFF_Uy(i,k,j) +diff_metric(i,k,j)
     &           + cori(i,j,n)*u(i,k,j,3-n,tau) 
     &              )*umask(i,k,j)
           enddo
         enddo
         u(i+1,:,:,n,taup1) = u(i,:,:,n,taup1)
         u(i+2,:,:,n,taup1) = u(i,:,:,n,taup1)
        else
         ie1=ie
        endif
#ifdef detailed_timing
        call toc('clinic obc')
#endif
c
c       and now the interior
c
#ifdef detailed_timing
        call tic('clinic tend')
#endif
        do j=js1,je1
          do k=1,km
            do i=is1,ie1
	      u(i,k,j,n,taup1) = (
     &           DIFF_Uz(i,k,j) 
     &           -  ADV_Ux(i,k,j) - ADV_Uy(i,k,j) - ADV_Uz(i,k,j) 
     &           + ADV_metric(i,k,j,n) 
     &           + DIFF_Ux(i,k,j) + DIFF_Uy(i,k,j) +diff_metric(i,k,j)
     &          - grad_p(i,k,j,n) + cori(i,j,n)*u(i,k,j,3-n,tau) 
     &              )*umask(i,k,j)
            enddo
          enddo
        enddo

#ifdef detailed_timing
        call toc('clinic tend')
#endif

#ifdef detailed_timing
        call tic('clinic impli')
#endif

#ifdef relax_traditional_approx
        if (n==1) then
         do j=js,je
          do k=1,km-1
            do i=is,ie
	      u(i,k,j,n,taup1) = u(i,k,j,n,taup1)+ 
     &                  2.0*omega*csu(j)*adv_vbu(i,k,j)
            enddo
          enddo
         enddo
        endif
#endif
c
c-----------------------------------------------------------------------
c       add in du/dt component due to implicit vertical diffusion
c-----------------------------------------------------------------------
c
        if (enable_implicit_vert_fric)
     &    call implicit_vert_fric(n,is,ie,js,je)

#ifdef detailed_timing
        call toc('clinic impli')
#endif
c
c-----------------------------------------------------------------------
c       construct the vertical average of du/dt for forcing
c       the barotropic equation
c-----------------------------------------------------------------------
c
        if (enable_freesurf) then

         do j=js,je
          do i=is,ie
           zu(i,j,n) =  .5*adv_fb(i,0,j)*umask(i,1,j)
          enddo
          do k=1,km
           do i=is,ie
	    zu(i,j,n) = zu(i,j,n) +
     &       (  u(i,k,j,n,taup1)-cori(i,j,n)*u(i,k,j,3-n,tau) )
#if defined partial_cell
     &                       *dhu(i,k,j) 
#else
     &                       *dzt(k) 
#endif

           enddo
          enddo
         enddo

        else

         zu(:,:,n) = 0.
         do j=js,je
          do k=1,km
            do i=is,ie
              zu(i,j,n) = zu(i,j,n) + u(i,k,j,n,taup1)
#ifdef partial_cell
     &                       *dhu(i,k,j)
#else
     &                       *dzt(k)
#endif
            enddo
          enddo
         enddo
         do j=js,je
          do i=is,ie
            zu(i,j,n) = zu(i,j,n)*hr(i,j)
          enddo
         enddo

        endif
c
c-----------------------------------------------------------------------
c       end of velocity component "n" loop
c-----------------------------------------------------------------------
c
      enddo ! n
c
c-----------------------------------------------------------------------
c     compute "tau+1" velocities accounting for implicit part of the
c     coriolis term if treated implicitly. velocities are in error by an
c     arbitrary constant related to neglecting the unknown surface
c     pressure gradients
c-----------------------------------------------------------------------
c
      do j=js,je
       do k=1,km
        do i=is,ie
         u(i,k,j,1,taup1) = u(i,k,j,1,taum1) + c2dt*u(i,k,j,1,taup1)
         u(i,k,j,2,taup1) = u(i,k,j,2,taum1) + c2dt*u(i,k,j,2,taup1)
        enddo
       enddo
      enddo
c
c-----------------------------------------------------------------------
c     subtract incorrect vertical means (related to ignoring horizontal
c     gradients of the surface pressure) to get pure internal modes.
c-----------------------------------------------------------------------
c
      do n=1,2
       do j=js,je
        baru=0.
        do k=1,km
         do i=is,ie
          baru(i) = baru(i) + u(i,k,j,n,taup1)*umask(i,k,j)
#ifdef partial_cell
     &                                    *dhu(i,k,j)
#else
     &                                    *dzt(k)
#endif
         enddo
        enddo
        if (enable_freesurf) then
c
c        account for changing total depth
c
         do k=1,km
          do i=is,ie
           u(i,k,j,n,taup1) = u(i,k,j,n,taup1)
     &                      - umask(i,k,j)*baru(i)*depthur(i,j)
          enddo
         enddo
        else
         do k=1,km
          do i=is,ie
           u(i,k,j,n,taup1) = u(i,k,j,n,taup1)
     &                      - umask(i,k,j)*baru(i)*hr(i,j)
          enddo
         enddo
        endif
       enddo

c       boundary exchange for taup1 is done in tropic or freesurf
c        call border_exchg(u(:,:,:,n,taup1),km,2)

      enddo

#ifdef detailed_timing
      call toc('clinic main')
#endif

      end subroutine clinic 





      subroutine delsqu (del2)
      use spflame_module
      implicit none
c
c=======================================================================
c     compute del**2 of momentum
c=======================================================================
c
      integer is,ie,js,je,i,j,k,n
      integer ism1,iep1,jsm1,jep1

      real :: del2(is_pe-1:ie_pe+1,km,js_pe-1:je_pe+1,2)
      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 diff_ux,diff_uy
      DIFF_Ux(i,k,j) = (diff_fe_biha(i,k,j) - diff_fe_biha(i-1,k,j))
     &                 *rho0csudxur(i,j)
#ifdef partial_cell
     &                 /dhu(i,k,j)
#endif
      DIFF_Uy(i,k,j) = (diff_fn_biha(i,k,j) - diff_fn_biha(i,k,j-1))
     &                 *rho0csudyur(j)
#ifdef partial_cell
     &                 /dhu(i,k,j)
#endif
#ifdef partial_cell
      real dhue,dhun
      dhue(i,k,j)   = min(dhu(i+1,k,j),dhu(i,k,j))
      dhun(i,k,j)   = min(dhu(i,k,j+1),dhu(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)
      
      do n=1,2
c
c-----------------------------------------------------------------------
c       diffusive flux across east face of "u" cell
c       diffusive flux across north face of "u" 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) = visc_ceu*csur(j)*dxtr(i+1)*
     &                         (u(i+1,k,j,n,taum1) - u(i,k,j,n,taum1))
     &                         *ambi_scale(i,j)
#ifdef partial_cell
     &                         *dhue(i,k,j)
#endif
            enddo
          enddo
       enddo

c      fn ism1: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) =visc_cnu*cst(j+1)*dytr(j+1)* 
     &                        (u(i,k,j+1,n,taum1) - u(i,k,j,n,taum1))
     &                         *ambi_scale(i,j+1)
#ifdef partial_cell
     &                         *dhun(i,k,j)
#endif
            enddo
          enddo
       enddo
c
c       compute -am*del**2
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,n) = (-DIFF_Ux(i,k,j) - DIFF_Uy(i,k,j) 
     &            - am3_biha(j)*u(i,k,j,n,taum1) 
     &                         *ambi_scale(i,j)
     &            - am4_biha(j,n)*
     &                         ambi_scale(i,j)*
     &         (u(i+1,k,j,3-n,taum1) - u(i-1,k,j,3-n,taum1))*dxmetr(i)
#if defined partial_cell
     &                        -pc_sink_biha(i,k,j)*u(i,k,j,n,taum1)
#endif
     &                        )*umask(i,k,j)
            enddo
          enddo

          if (my_blk_i == 1) then
            if (enable_obc_west) then
             do k=1,km
               del2(1,k,j,n) = del2(2,k,j,n) + DIFF_Ux(2,k,j)
             enddo
            else
             del2(1,:,j,n) = 0.
            endif
          endif
          if (my_blk_i == n_pes_i) then
            if (enable_obc_east) then
             do k=1,km
               del2(imt,k,j,n) = del2(imt-1,k,j,n) + DIFF_Ux(imt-1,k,j)
             enddo
            else
             del2(imt,:,j,n) = 0.
            endif
          endif
       enddo

       call set_cyclic(del2(:,:,:,n),km,1)

       if (my_blk_j == 1)  then
         if (enable_obc_south) then
          do k=1,km
           do i=is-1,ie+1
            del2(i,k,1,n) = del2(i,k,2,n) + DIFF_Uy(i,k,2) 
           enddo
          enddo
         else
          del2(is-1:ie+1,:,1,n) = 0.
         endif
       endif

       if (my_blk_j == n_pes_j)  then
         if (enable_obc_north) then
          do k=1,km
            do i=is-1,ie+1
	      del2(i,k,jmt,n) = del2(i,k,jmt-1,n) + DIFF_Uy(i,k,jmt-1) 
            enddo
          enddo
         else
	  del2(is-1:ie+1,:,jmt,n) = 0.
         endif
       endif

      enddo ! n

      end subroutine delsqu 






      subroutine pressure_gradient(grad_p)
      use spflame_module
      implicit none
      real :: grad_p(is_pe:ie_pe,km,js_pe:je_pe,2) 
      integer :: i,j,k,is,ie,js,je,kb
      real    :: fxa,fxb,t1,t2,t3,t4

#ifdef partial_cell
      real :: rhoij,rhoip1,rhojp1,rhoip1jp1
      real :: pressure(is_pe:ie_pe+1,km,js_pe:je_pe+1)
#endif
      integer :: ii,ik,ij
      real :: tempik
      tempik(ii,ik,ij) = rho(ii,ik-1,ij) + rho(ii,ik,ij)  !kk
c

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c-----------------------------------------------------------------------
c     construct the hydrostatic pressure gradients: 1 = dp/dx; 2 = dp/dy
c-----------------------------------------------------------------------
c
#ifdef detailed_timing
      call tic('clinic density')
#endif

      do j=js-1,je+1
        do k=1,km
         call model_dens(t(is-1,k,j,1,tau), t(is-1,k,j,2,tau),
     &                   rho(is-1,k,j),k,ie-is+3
#ifdef partial_cell
     &                       ,ztp(is-1,k,j)
#endif
     &                       )
        enddo
      enddo

      if (enable_blue)      call add_blue_to_rho()
      if (enable_blue_mean) call add_blue_mean_to_rho()

#ifdef detailed_timing
      call toc('clinic density')
#endif

#ifndef partial_cell
c
c     compute horizontal pressure gradient at the first level
c
      do j=js,je
        fxa  = grav*rho0r*dzw(0)*csur(j)
        fxb  = grav*rho0r*dzw(0)*dyu2r(j)
        do i=is,ie
          t1              = rho(i+1,1,j+1) - rho(i  ,1,j)
          t2              = rho(i  ,1,j+1) - rho(i+1,1,j)
          grad_p(i,1,j,1) = (t1-t2)*fxa*dxu2r(i)
          grad_p(i,1,j,2) = (t1+t2)*fxb
        enddo
      enddo

c
c     compute the change in pressure gradient between levels
c     (faster this way)
c
      do j=js,je
        fxa = grav*rho0r*csur(j)*0.5
        fxb = grav*rho0r*dyu4r(j)
        do k=2,km
          do i=is,ie
            t1              = tempik(i+1,k,j+1) - tempik(i  ,k,j)
            t2              = tempik(i  ,k,j+1) - tempik(i+1,k,j)
            t3              = fxa*(t1-t2)*dxu2r(i)*dzw(k-1)
            t4              = fxb*(t1+t2)   *dzw(k-1)
            grad_p(i,k,j,1) = grad_p(i,k-1,j,1) + t3
            grad_p(i,k,j,2) = grad_p(i,k-1,j,2) + t4
          enddo
        enddo
      enddo

#else
c
c     compute hydrostatic pressure anomaly
c
      do j=js,je+1
	do i=is,ie+1
         pressure(i,1,j) = rho(i,1,j)*grav*dzw(0)
	 kb = max(kmt_big(i,j),1)
         do k=2,kb
          pressure(i,k,j) = pressure(i,k-1,j) 
     &              + (rho(i,k-1,j)+rho(i,k,j)) 
     &                  *0.5*grav*dhwt(i,k-1,j)
         enddo
	 fxa = pressure(i,kb,j)
         do k=kb+1,km
            pressure(i,k,j) = fxa
         enddo
       enddo
      enddo
c
c     compute pressure gradient
c
      do j=js,je
        fxa  = rho0r*csur(j)
        fxb  = rho0r*dyu2r(j)
	do k=1,km
          do i=is,ie
            t1              = pressure(i+1,k,j+1) - pressure(i  ,k,j)
            t2              = pressure(i  ,k,j+1) - pressure(i+1,k,j)
            grad_p(i,k,j,1) = (t1-t2)*fxa*dxu2r(i)*umask(i,k,j)
            grad_p(i,k,j,2) = (t1+t2)*fxb*umask(i,k,j)
          enddo
        enddo
      enddo
c
c     add in correction due to variable thickness partial cell
c
      do j=js,je
        fxa = grav*rho0r*csur(j)
        fxb = grav*rho0r*dyu4r(j)
        do i=is,ie
          k = kmu(i,j)
	  if (k .gt. 1) then
	    rhoij     = rho(i,k,j)
	    rhoip1    = rho(i+1,k,j)
	    rhojp1    = rho(i,k,j+1)
	    rhoip1jp1 = rho(i+1,k,j+1)
            grad_p(i,k,j,1) = grad_p(i,k,j,1) - umask(i,k,j)*fxa
     &         *dxu4r(i)*((rhoip1 + rhoij)*(ztp(i+1,k,j) - ztp(i,k,j))
     &           +(rhoip1jp1 + rhojp1)*(ztp(i+1,k,j+1) - ztp(i,k,j+1)))
            grad_p(i,k,j,2) = grad_p(i,k,j,2) - umask(i,k,j)*fxb*(  
     &                    (rhojp1 + rhoij)*(ztp(i,k,j+1) - ztp(i,k,j))
     &           +(rhoip1jp1 + rhoip1)*(ztp(i+1,k,j+1) - ztp(i+1,k,j)))
          endif
	enddo
      enddo
#endif
c

      end subroutine pressure_gradient
