#include "options.inc"


      subroutine vmixc()
c
c=======================================================================
c     set the vertical diffusivity and viscosity
c     by choosing one of several possible parameterisations
c=======================================================================
c
      use spflame_module
      implicit none

      if (enable_cgh_vert_mixing) then
        call cghmix()
      elseif (enable_const_vert_mixing) then
        visc_cbu=kappa_m
        diff_cbt=kappa_h
        visc_cbu(:,km,:)=0.
        diff_cbt(:,km,:)=0.
        if (enable_impl_convection) then
          call implicit_convection()
        endif
      elseif (enable_tkemix) then
        call tkemix()
      endif

      if (enable_diffusion_isoneutral) call isoneutral_add_K33  
      if (enable_diffusion_isopycnic)  call isopycnic_add_K33  
      end subroutine vmixc


      subroutine setvbc
c
c=======================================================================
c     set the vertical boundary conditions
c=======================================================================
c
      use spflame_module
      implicit none
      real            :: uvmag,tshift
      integer         :: i,j,is,ie,js,je,n,kz
      real, parameter :: rinf = 1.0e+32
      real :: frozen(is_pe:ie_pe),cd

      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
c
c     the surface fluxes
c
      if (mother<0) then

#ifdef couple
#else
       if (enable_simple_sbc) then

         call sbc_template()

       else
c
c       read in SBC conditions
c
        call sbc_read

       endif ! enable_simple_sbc
#endif

      else
c
c      sub domains get taux/tauy, ustar, stf_clim and stf_rest
c      from mother domain
c
      endif ! mother <0


#ifndef couple
c
c     apply restoring boundary conditions for T and S
c
      do j=js,je
       do i=is,ie
        stf(i,j,1)=stf_rest(i,j,1)*(stf_clim(i,j,1)-t(i,1,j,1,taum1))
        stf(i,j,2)=stf_rest(i,j,2)*(stf_clim(i,j,2)-t(i,1,j,2,taum1))
       enddo
      enddo

      if (nt>2) then
c
c      apply surface boundary conditions for other tracers
c
       call passive_tracer_sflx()

      endif

#endif


c
c     the bottom fluxes are zero
c
      btf=0.; bmf=0.
c
c     set bottom drag 
c
      do j=js,je
       do i=is,ie
        kz = kmu(i,j)
        if (kz .ne. 0) then
         uvmag    = sqrt(u(i,kz,j,1,taum1)**2 + 
     &                   u(i,kz,j,2,taum1)**2 + tidaloff**2)
         bmf(i,j,1) = cdbot*u(i,kz,j,1,taum1)*uvmag
         bmf(i,j,2) = cdbot*u(i,kz,j,2,taum1)*uvmag
        endif
       enddo
      enddo

      if (enable_windstress_feedback) then
       do j=js,je
        do i=is,ie
         kz = kmu(i,j)
         if (kz .ne. 0) then
          uvmag = sqrt(u(i,kz,j,1,taum1)**2 + 
     &                 u(i,kz,j,2,taum1)**2 + 1.e-8)/100. ! in m/s
          cd = (0.52+1.56/uvmag)/1000. *1.35 *10./100. ! C_d * rho *10/100 (to convert )
          smf(i,j,1) = smf(i,j,1)+cd*u(i,1,j,1,taum1)*uvmag
          smf(i,j,2) = smf(i,j,2)+cd*u(i,1,j,2,taum1)*uvmag
         endif
        enddo
       enddo
      endif


      if (enable_icemask) then
c
c      apply ice mask, momentum and salt fluxes are not affected
c
       do j=js,je
        call freezing_point(t(is_pe,1,j,2,tau),frozen,ie_pe-is_pe+1)
        do i=is,ie
c        set tshift  = 0    when T <= frozen
c            tshift  = huge when T  > frozen 
         tshift = rinf * amax1((t(i,1,j,1,tau)-frozen(i)), 0.) 
c        set icemask = 0 when tshift = 0 
c            icemask = 1 when tshift = huge
         icemask(i,j) = amin1(tshift,1.)* tmask(i,1,j)
        enddo
       enddo
c
c     apply icemask to heat fluxes only if there are cooling
c
       where (stf(is:ie,js:je,1) < 0.)
     &        stf(is:ie,js:je,1) =
     &        stf(is:ie,js:je,1) * icemask(is:ie,js:je)
c
c      apply icemask to other tracer fluxes
c
       do n=2,nt
        stf(is:ie,js:je,n) =
     &  stf(is:ie,js:je,n) * icemask(is:ie,js:je)
       enddo
c
c     apply icemask to wind stirring
c
       if (enable_ktmix.or.enable_tkemix) then
        ustar(is:ie,js:je) = ustar(is:ie,js:je)*icemask(is:ie,js:je)
       endif

      endif  ! enable icemask

      end subroutine setvbc



      subroutine init_cghmix 
      use spflame_module
      implicit none
c
c=======================================================================
c     Initialization for the Cummins/Holloway/Gargett vertical mixing
c     scheme
c     Cummins, Holloway and Gargett (JPO, 20, 1990, 817-830).
c    
c     Christian Dieterich Tue Nov 11 17:25:58 MET 1997
c     adapted from "ppmixi" by
c     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
c     SPFLAME version : c.eden
c=======================================================================
c
      if (my_pe==0) then
       print*,''
       print*,' Initialization of Cummins et al vertical mixing scheme' 
       print*,''
       print*,' proportionality factors to the inverse of the'
       print*,' square root of Brunt-Vaisaelae frequency ='
       print*,' (for diffusivity/viscosity) ',cgh_vdcfac,' '
       print*,' min value for mixing at surface to simulate high freq'
       print*,' wind mixing (if absent in forcing)= ',wndmix,' cm^2/s'
       print*,' diffusivity is bounded by ',
     &          diff_cbt_back,diff_cbt_cut,' cm^2/s'
       print*,' viscosity is bounded by ',
     &          visc_cbu_back,visc_cbu_cut,' cm^2/s'
       if (enable_cgh_vert_momentum_mixing) then
         print*,' using Cummins et al also for momentum'
       endif
       if (enable_cgh_impl_convection) then
         print*,' using high vertical diffusivities in case of'
         print*,' an unstable stratification '
       endif
       print*,''
       print*,' done' 
       print*,''
      endif
      end subroutine init_cghmix 


      subroutine cghmix
      use spflame_module
      implicit none
c
c=======================================================================
c     Compute vertical mixing coefficients based on...
c     Cummins, Holloway and Gargett (JPO, 20, 1990, 817-830).
c
c     Note: By default only vertical diffusivities are computed as a
c     function of Brunt-Vaisaelae frequency. Vertical viscosity is set
c     to the constant value "kappa_m".
c     If "cghvmix_momentum" is enabled the vertical viscosities also
c     dependend on stability.  
c
c     inputs:
c
c      cgh_vdcfac     = proportionality factor to the inverse of the
c                       square root of Brunt-Vaisaelae frequency
c      diff_cbt_cut   = maximum "diff_cbt" (cm**2/sec)
c      diff_cbt_back  = background "diff_cbt" (cm**2/sec)
c if enable_cgh_vert_momentum_mixing
c      wndmix         = min value for mixing at surface to simulate high freq
c                       wind mixing (if absent in forcing). (cm**2/sec)
c      visc_cbu_cut   = maximum "visc_cbu" (cm**2/sec)
c      visc_cbu_back  = background "visc_cbu" (cm**2/sec)
c else
c      kappa_m        = constant "visc_cbu" (cm**2/sec)
c endif
c
c     outputs:
c
c      visc_cbu = viscosity coefficient at bottom of "u" cells (cm**2/s)
c      diff_cbt = diffusion coefficient at bottom of "t" cells (cm**2/s)
c
c     Christian Dieterich Tue Nov 11 17:02:23 MET 1997
c     adapted from "ppmix" by
c     author:      r. c. pacanowski      e-mail=> rcp@gfdl.gov
c=======================================================================
c
      integer :: is,ie,js,je,i,j,k
      real    :: ro1(is_pe:ie_pe+1),ro2(is_pe:ie_pe+1)
      real    :: bvfrqs(is_pe:ie_pe+1,km,js_pe:je_pe+1)
      real    :: fx,diff_loc, visc_loc
      real    :: bvfrqs_cbu
      real    :: diff_cbt_limit

      diff_cbt_limit = diff_cbt_cut
      if (enable_cgh_impl_convection) diff_cbt_limit = 1.0e6
         
      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     compute squared Brunt-Vaisaelae frequency across bottom of "t"
c     cells at tau-1 for rows js through je in the MW. Set squared
c     Brunt-Vaisaelae frequency = zero across bottom and in land areas 
c-----------------------------------------------------------------------
c
      do j=js,je+1
       do k=1,km-1
        call model_dens(t(is,k  ,j,1,taum1),t(is,k  ,j,2,taum1),
     &                  ro1(is),k,ie-is+2
#ifdef partial_cell
     &                       ,ztp(is,k,j)
#endif
     &                       )
        call model_dens(t(is,k+1,j,1,taum1),t(is,k+1,j,2,taum1),
     &                  ro2(is),k,ie-is+2
#ifdef partial_cell
     &                       ,ztp(is,k,j)
#endif
     &                       )
        do i=is,ie+1
#ifdef partial_cell
         fx = -grav/(rho0*dhwt(i,k,j))
#else
         fx = -grav/rho0*dzwr(k)
#endif
         bvfrqs(i,k,j) =fx*(ro1(i)-ro2(i))*tmask(i,k+1,j)
        enddo
       enddo
      enddo
c
c-----------------------------------------------------------------------
c     set no-flux condition on density difference across bottom level
c-----------------------------------------------------------------------
c
      bvfrqs(:,km,:) = 0.
c
c-----------------------------------------------------------------------
c     compute vertical diffusivity coeff on "t" cell bottoms
c     in regions of gravitational instability, reset the vertical
c     mixing coefficient to the limit
c     Cummins, Holloway and Gargett (JPO, 20, 1990, 817-830). Eq: 10
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km-1
          do i=is,ie
            diff_cbt(i,k,j) = diff_cbt_cut
            if (bvfrqs(i,k,j) > 0.) then
              diff_loc = cgh_vdcfac(1)/sqrt(bvfrqs(i,k,j))
              diff_loc = max(diff_loc, diff_cbt_back)
              diff_loc = min(diff_loc, diff_cbt_cut)
              diff_cbt(i,k,j) = diff_loc
            elseif (bvfrqs(i,k,j) < 0.) then !????? ----- <<<<<
c             set kappa to a high value for unstable stratification
c             (just in case of enable_cgh_impl_convection = .true.)
c             (else diff_cbt_cut )
              diff_cbt(i,k,j) = diff_cbt_limit
            endif
          enddo
        enddo
      enddo

      if (enable_cgh_vert_momentum_mixing) then
c
c-----------------------------------------------------------------------
c     compute vertical viscosity coeff on "u" cell bottoms
c     in regions of gravitational instability, reset the vertical
c     mixing coefficient to the limit
c     Cummins, Holloway and Gargett (JPO, 20, 1990, 817-830). Eq: 10
c-----------------------------------------------------------------------
c
       do j=js,je
        do k=1,km-1
          do i=is,ie
            visc_cbu(i,k,j) = visc_cbu_cut
            bvfrqs_cbu = .25*(bvfrqs(i,k,j)+
     &           bvfrqs(i+1,k,j)+bvfrqs(i,k,j+1)+bvfrqs(i+1,k,j+1))
            if (bvfrqs_cbu .gt. 0.) then
              visc_loc = cgh_vdcfac(2)/sqrt(bvfrqs_cbu)
              visc_loc = max(visc_loc, visc_cbu_back)
              visc_loc = min(visc_loc, visc_cbu_cut)
              visc_cbu(i,k,j) = visc_loc
            else
              visc_cbu(i,k,j) = visc_cbu_cut
            endif
          enddo
        enddo
       enddo
      else
c
c-----------------------------------------------------------------------
c     set vertical viscosity coeff on "u" cell bottoms
c-----------------------------------------------------------------------
c
       visc_cbu=kappa_m
      endif
c
c-----------------------------------------------------------------------
c     if enable_cgh_vert_momentum_mixing
c      approximation for high freq wind mixing near the surface
c     endif
c     set no flux through bottom of bottom level "km"
c-----------------------------------------------------------------------
c
      if (enable_cgh_vert_momentum_mixing) then
       do j=js,je
        do i=is,ie
          if (visc_cbu(i,1,j) .lt. wndmix) visc_cbu(i,1,j) = wndmix
c   account for deeper ekman layer in sub_domain with higher
c   vertical resolution
          if (mother >= 0 .and. km>km_mother) then
           do k=1,zoom_fac_k
            if (visc_cbu(i,k,j) .lt. wndmix) visc_cbu(i,k,j) = wndmix
           enddo
          endif
        enddo
       enddo
      endif

      do j=js,je
        do i=is,ie
	  diff_cbt(i,km,j) = 0.
	  visc_cbu(i,km,j) = 0.
        enddo
      enddo
c
c-----------------------------------------------------------------------
c       set diffusion and viscosity coeffs to zero on land box bottoms
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km-1
          do i=is,ie
            visc_cbu(i,k,j) = umask(i,k+1,j)*visc_cbu(i,k,j)
            diff_cbt(i,k,j) = tmask(i,k+1,j)*diff_cbt(i,k,j)
          enddo
        enddo
      enddo
      end subroutine cghmix

      
      subroutine implicit_convection
      use spflame_module
      implicit none
c
c=======================================================================
c     set kappa to diff_cbt_cut for unstable stratification
c=======================================================================
c
      integer :: is,ie,js,je,i,j,k
      real    :: ro1(is_pe:ie_pe+1),ro2(is_pe:ie_pe+1)
      real    :: bvfrqs(is_pe:ie_pe+1,km,js_pe:je_pe+1)
      real    :: fx

         
      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     compute squared Brunt-Vaisaelae frequency across bottom of "t"
c     cells at tau-1 for rows js through je in the MW. Set squared
c     Brunt-Vaisaelae frequency = zero across bottom and in land areas 
c-----------------------------------------------------------------------
c
      do j=js,je+1
       do k=1,km-1
        call model_dens(t(is,k  ,j,1,taum1),t(is,k  ,j,2,taum1),
     &                  ro1(is),k,ie-is+2
#ifdef partial_cell
     &                       ,ztp(is,k,j)
#endif
     &                       )
        call model_dens(t(is,k+1,j,1,taum1),t(is,k+1,j,2,taum1),
     &                  ro2(is),k,ie-is+2
#ifdef partial_cell
     &                       ,ztp(is,k,j)
#endif
     &                       )
        do i=is,ie+1
#ifdef partial_cell
         fx = -grav/(rho0*dhwt(i,k,j))
#else
         fx = -grav/rho0*dzwr(k)
#endif
         bvfrqs(i,k,j) =fx*(ro1(i)-ro2(i))*tmask(i,k+1,j)
        enddo
       enddo
      enddo
c
c-----------------------------------------------------------------------
c     set no-flux condition on density difference across bottom level
c-----------------------------------------------------------------------
c
      bvfrqs(:,km,:) = 0.
c
c-----------------------------------------------------------------------
c     set kappa to a high value for unstable stratification
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km-1
          do i=is,ie
            if (bvfrqs(i,k,j) <=0.) then
              diff_cbt(i,k,j) = diff_cbt_cut
            endif
          enddo
        enddo
      enddo

      do j=js,je
        do i=is,ie
	  diff_cbt(i,km,j) = 0.
        enddo
      enddo
c
c-----------------------------------------------------------------------
c       set diffusion coeffs to zero on land box bottoms
c-----------------------------------------------------------------------
c
      do j=js,je
        do k=1,km-1
          do i=is,ie
            diff_cbt(i,k,j) = tmask(i,k+1,j)*diff_cbt(i,k,j)
          enddo
        enddo
      enddo
      end subroutine implicit_convection


#if !defined vector_host || !defined flame_density || defined partial_cell

      subroutine convect (tauxx)
c
c-----------------------------------------------------------------------
c     The following convection scheme is an alternative to the
c     standard scheme. In contrast to the standard scheme,
c     it totally removes all gravitational instability in the
c     water column. It does that in one pass, so the parameter
c     ncon becomes irrelevant if this option is selected.
c     The scheme is equivalent to those used by Rahmstorf 
c     (jgr 96,6951-6963) and by Marotzke (jpo 21,903-907).
c     It is discussed in a note to Ocean Modelling (101). It uses
c     as much cpu time as 1-3 passes of the standard scheme, 
c     depending on the amount of static instability found in the
c     model, and is much faster than using "implicitvmix".
c
c     Written by Stefan Rahmstorf, Institut fuer Meereskunde,
c     Kiel, Germany.            Comments welcome:
c                         srahmstorf@meereskunde.uni-kiel.d400.de
c
c     definition of internal variables:
c     kcon = maximum number of levels at this location
c     lcon = counts levels down
c     lcona = upper layer of a convective part of water column
c     lconb = lower layer of a convective part of water column
c     rhoup = density referenced to same level
c     rholo = density referenced to level below
c                (note that densities are not absolute!)
c     dztsum = sum of layer thicknesses
c     trasum = sum of layer tracer values
c     tramix = mixed tracer value after convection
ce    the following two arrays are now in spflame_module
c     lctot = total of number of levels involved in convection
c     lcven = number of levels ventilated (convection to surface)
c     note: lctot can in rare cases count some levels twice, if they
c           get involved in two originally separate, but then
c           overlapping convection areas in the water column! It
c           is a control parameter; the sensible parameter to plot
c           is lcven. Lcven is 0 on land, 1 on ocean points with no
c           convection, and anything up to km on convecting points. 
c
c     Written by Stefan Rahmstorf, Institut fuer Meereskunde,
c                         srahmstorf@meereskunde.uni-kiel.d400.de
c     SPFLAME version : c.eden
c-----------------------------------------------------------------------
c
      use spflame_module
      implicit none
      integer :: is,ie,js,je,tauxx
      real :: rhoup(is_pe:ie_pe,km), rholo(is_pe:ie_pe,km)
      real :: trasum(nt), tramix,dztsum
      integer :: n,i,j,k,l,l1,kcon,lcon,lcona,lconb,lmix
c
c     check each row column by column; note that 'goto 1310'
c     finishes a particular column and moves to the next one.
c
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)

      do j=js,je
c
c       find density of entire row for stability determination
c
        do l=1,km-1
          l1=l+1
          call model_dens(t(is,l1,j,1,tauxx), t(is,l1,j,2,tauxx),
     &                    rhoup(is,l1),l1,ie-is+1
#ifdef partial_cell
     &                       ,ztp(is,l1,j)
#endif
     &                       )
          call model_dens(t(is,l ,j,1,tauxx), t(is,l ,j,2,tauxx),
     &                    rholo(is,l ),l1,ie-is+1
#ifdef partial_cell
     &                       ,ztp(is,l1,j)
#endif
     &                       )
	enddo

        do i=is,ie
          kcon     = kmt(i,j)
          lctot(i,j) = 0
          lcven(i,j) = 0
	  if (kcon == 0) goto 1310
          lcven(i,j) = 1
          lcon     = 0
c
c         1. initial search for uppermost unstable pair; if none is
c            found, move on to next column
c
          do k=kcon-1,1,-1
            if (rholo(i,k) .gt. rhoup(i,k+1)) lcon = k
          enddo
c
          if (lcon .eq. 0) goto 1310
c
1319      lcona = lcon
          lconb = lcon + 1
c
c         2. mix the first two unstable layers
c
#ifdef partial_cell
          dztsum = dht(i,lcona,j) + dht(i,lconb,j)
#else
          dztsum = dzt(lcona) + dzt(lconb)
#endif
          do n=1,nt
#ifdef partial_cell
            trasum(n)        = t(i,lcona,j,n,tauxx)*dht(i,lcona,j) + 
     &                         t(i,lconb,j,n,tauxx)*dht(i,lconb,j)
#else
            trasum(n)        = t(i,lcona,j,n,tauxx)*dzt(lcona) + 
     &                         t(i,lconb,j,n,tauxx)*dzt(lconb)
#endif
            tramix           = trasum(n) / dztsum
            t(i,lcona,j,n,tauxx) = tramix
            t(i,lconb,j,n,tauxx) = tramix
          enddo
c
c         3. test layer below lconb
c
1306      continue
          if (lconb .eq. kcon) goto 1308
c
          l1 = lconb + 1
          rholo(i,lconb) = model_dens_scalar (t(i,lconb,j,1,tauxx)
     &,                              t(i,lconb,j,2,tauxx), l1
#ifdef partial_cell
     &                       ,ztp(i,l1,j)
#endif
     &                       )
c
          if (rholo(i,lconb) .gt. rhoup(i,l1)) then
            lconb = lconb+1
#ifdef partial_cell
            dztsum = dztsum + dht(i,lconb,j)
#else
            dztsum = dztsum + dzt(lconb)
#endif
            do n=1,nt
#ifdef partial_cell
          trasum(n) = trasum(n) + t(i,lconb,j,n,tauxx)*dht(i,lconb,j)
#else
              trasum(n) = trasum(n) + t(i,lconb,j,n,tauxx)*dzt(lconb)
#endif
              tramix = trasum(n) / dztsum
              do lmix=lcona,lconb
                t(i,lmix,j,n,tauxx) = tramix
              enddo
            enddo
            goto 1306
          end if
c
c         4. test layer above lcona
c
1308      continue
          if (lcona .gt. 1) then
            l1 = lcona-1
            rholo(i,l1) = model_dens_scalar(t(i,l1,j,1,tauxx) 
     &,                            t(i,l1,j,2,tauxx),lcona
#ifdef partial_cell
     &                       ,ztp(i,lcona,j)
#endif
     &                       )
            rhoup(i,lcona) = model_dens_scalar(t(i,lcona,j,1,tauxx)
     &,                               t(i,lcona,j,2,tauxx),lcona
#ifdef partial_cell
     &                       ,ztp(i,lcona,j)
#endif
     &                       )
            if (rholo(i,lcona-1) .gt. rhoup(i,lcona)) then
              lcona = lcona-1
#ifdef partial_cell
              dztsum = dztsum + dht(i,lcona,j)
#else
              dztsum = dztsum + dzt(lcona)
#endif
              do n=1,nt
#ifdef partial_cell
            trasum(n) = trasum(n) + t(i,lcona,j,n,tauxx)*dht(i,lcona,j)
#else
                trasum(n) = trasum(n) + t(i,lcona,j,n,tauxx)*dzt(lcona)
#endif
                tramix = trasum(n) / dztsum 
                do lmix=lcona,lconb
                  t(i,lmix,j,n,tauxx) = tramix
                enddo
              enddo
              goto 1306
            end if
          end if
c
c         5. remember the total number of levels mixed by convection
c            in this water column, as well as the ventilated column
c
          lctot(i,j) = lctot(i,j) + lconb - lcona + 1
          if (lcona .eq. 1) lcven(i,j) = lconb - lcona + 1
c
c         6. resume search if step 3. and 4. have been passed and this
c            unstable part of the water column has thus been removed,
c            i.e. find further unstable areas further down the column
c
          if (lconb .eq. kcon) goto 1310
          lcon = lconb
c
1302      continue
          lcon = lcon + 1
c
          if (lcon .eq. kcon) goto 1310
c
          if (rholo(i,lcon) .le. rhoup(i,lcon+1)) goto 1302
c
          goto 1319

1310      continue
        enddo
      enddo
      end subroutine convect

#else

c !!! partly vectorized version !!!

      subroutine convect(tauxx)

!-----------------------------------------------------------------------
!   the following code acts only for the first kmax levels
!   therefore, another scheme is necessary below 
!   intention is to use it with ktmix scheme (see below) only
!
!   Changes made at CSIRO Marine Research:
!   The main loop of the original routine cannot be vectorized. However,
!   it is possible to vectorize part of the algorithm. If convective mixing
!   occurs it is nearly always due to surface cooling rather than an internal 
!   instability. We need only search in downwards in this case. As soon as
!   a column is stabilised the search is discontinued for that column. Once
!   all columns have been 'ventilated' we then proceed to the original
!   algorithm. Since we have removed the surface instabilities it is highly
!   unlikely that we enter the main loop. If no convection at all occurs then
!   there is a small overhead equivalent to calculating a row of densities.
!   It may be worthwhile escaping the vectorized ventilation section when the
!   number of unstable columns is below some threshold and enter an optimised
!   scalar section of code.
!
!   Note: Some Fortran90 constructs have been used for clarity. It is
!   straightforward to recode in f77.
!   
!   icasc  = logical variable indicating whether a column is still unstable
!   lcven1 = levels ventilated in in new 'cascade' sweep. Plot the maximum
!            of lcven and lcven1 now.
!
!   author: Russ Fiedler fiedler@marine.csiro.au
!
c     SPFLAME version : c.eden
!-----------------------------------------------------------------------

      use spflame_module
      implicit none
      integer :: tauxx
      integer,parameter :: kmax=15
      integer :: is,ie,js,je,i,j,k,k2,n,lev,l1
      integer :: kcon,lcon(is_pe:ie_pe),lcona,lconb,lmix
      integer :: lcven1(is_pe:ie_pe)
      real :: rhoup(is_pe:ie_pe,km), rholo(is_pe:ie_pe,km)
      real :: trasum(nt), tramix,dztsum
      real :: dzsum, dzsum1, ratio1, ratio2
      logical :: icasc(is_pe:ie_pe)
      real :: tq,sq,dens
      dens(tq,sq,k)  = (dn_coef(k,1) + 
     &    (dn_coef(k,4)+dn_coef(k,7)*sq)*sq +
     &    (dn_coef(k,3)+dn_coef(k,8)*sq+dn_coef(k,6)*tq)*tq)*tq +
     &    (dn_coef(k,2)+(dn_coef(k,5)+ dn_coef(k,9)*sq)*sq)*sq

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

      if (kmax >km ) then
        if(my_pe==0) print*,' parameter kmax greater than km  '
        call halt_stop(' in convect (vector version)')
      endif

#ifdef partial_cell
        call halt_stop(' in convect (vector version) -- 4')
#endif
#ifndef flame_density
        call halt_stop(' in convect (vector version) -- 5')
#endif

      do j=js,je

! 0. Pre search for upper levels. And allow cascade downwards.
!    Most convection will occur here. if we can stabilize
!    then the next search will skip the scalar loop

        lcven1(is:ie)=0.0
        do i=is,ie
          icasc(i)=.false.
          if ( kmt(i,j) /= 0 ) icasc(i)=.true.
        enddo

        dzsum1=dzt(1)
        do lev=1,kmax-1
          dzsum=dzsum1
          dzsum1=dzsum1+dzt(lev+1)
          ratio1=dzsum/dzsum1
          ratio2=dzt(lev+1)/dzsum1
          do i=is,ie
! Test if cascade is possible
            if( lev >= kmt(i,j) ) icasc(i) = .false.
            if (icasc(i)) then
          rhoup(i,lev+1)=dens(t(i,lev+1,j,1,tauxx)-dn_to(lev+1)
     &,                       t(i,lev+1,j,2,tauxx)-dn_so(lev+1),lev+1)
          rholo(i,lev)  =dens(t(i,lev  ,j,1,tauxx)-dn_to(lev+1)
     &,                       t(i,lev  ,j,2,tauxx)-dn_so(lev+1),lev+1)
! Test whether to mix else mark end of cascade
! NB only need to adjust lower value. Will do last part of mixing later.
              if( rholo(i,lev) >= rhoup(i,lev+1) ) then
                t(i,lev+1,j,1,tauxx)   = t(i,lev  ,j,1,tauxx)*ratio1 + 
     &                                   t(i,lev+1,j,1,tauxx)*ratio2
                t(i,lev+1,j,2,tauxx)   = t(i,lev  ,j,2,tauxx)*ratio1 + 
     &                                   t(i,lev+1,j,2,tauxx)*ratio2
                lcven1(i)=lev+1      
              else
                icasc(i)=.false.
              endif
            endif
          enddo ! i
          if( .not. any(icasc(is:ie)) ) exit
        enddo ! lev

!     the other tracers
       do n = 3, nt
        dzsum1=dzt(1)
        do lev=1,kmax-1
          dzsum=dzsum1
          dzsum1=dzsum1+dzt(lev+1)
          ratio1=dzsum/dzsum1
          ratio2=dzt(lev+1)/dzsum1
          do i=is,ie
            if(lev <= lcven1(i)-1) then
                t(i,lev+1,j,n,tauxx)   = t(i,lev  ,j,n,tauxx)*ratio1 + 
     &                                   t(i,lev+1,j,n,tauxx)*ratio2
            endif
          enddo
        enddo
       enddo

! Adjust values up the column 

       do n = 1, nt
        do lev=1,maxval(lcven1(is:ie))-1
          do i=is,ie
            if(lev.le.lcven1(i)-1) then
               t(i,lev,j,n,tauxx) = t(i,lcven1(i),j,n,tauxx)
            endif
          enddo
        enddo
       enddo

! Redo all the densities 

        do k=1,kmax-1
         do i=is,ie
          k2=k+1
          rhoup(i,k2) = dens(t(i,k2,j,1,tauxx)-dn_to(k2)
     &,                      t(i,k2,j,2,tauxx)-dn_so(k2), k2)
          rholo(i,k)  = dens(t(i,k ,j,1,tauxx)-dn_to(k2)
     &,                      t(i,k ,j,2,tauxx)-dn_so(k2), k2)
	 enddo
	enddo

!     check each row column by column; note that 'goto 1310'
!     finishes a particular column and moves to the next one.

!         1. initial search for uppermost unstable pair;

        lcon(is:ie)=0
        lctot(is:ie,j) = 0
        lcven(is:ie,j) = 0

        do k=kmax-1,1,-1
          do i=is,ie
            if ( rholo(i,k) > rhoup(i,k+1)) lcon(i) = k
          enddo
        enddo
! Correct for land and get flag to bail out of 
        do i=is,ie
          if( lcon(i) >= kmt(i,j) .or. lcon(i) >= kmax ) lcon(i)=0
        enddo

       if( maxval(lcon(is:ie)) /= 0 ) then 

        do i=is,ie
          if ( lcon(i) == 0 ) goto 1310

          kcon     = min(kmax,kmt(i,j))
          lcven(i,j) = 1

1319      lcona = lcon(i)
          lconb = lcon(i) + 1

!         2. mix the first two unstable layers

          dztsum = dzt(lcona) + dzt(lconb)
          do n=1,nt
            trasum(n)        = t(i,lcona,j,n,tauxx)*dzt(lcona) + 
     &                         t(i,lconb,j,n,tauxx)*dzt(lconb)
            tramix           = trasum(n) / dztsum
            t(i,lcona,j,n,tauxx) = tramix
            t(i,lconb,j,n,tauxx) = tramix
          enddo

!         3. test layer below lconb

1306      continue
          if ( lconb == kcon ) goto 1308

          l1 = lconb + 1
          rholo(i,lconb) = dens(t(i,lconb,j,1,tauxx)-dn_to(l1)
     &,                      t(i,lconb,j,2,tauxx)-dn_so(l1), l1)
          if ( rholo(i,lconb) > rhoup(i,l1) ) then
            lconb = lconb+1
            dztsum = dztsum + dzt(lconb)
            do n=1,nt
              trasum(n) = trasum(n) + t(i,lconb,j,n,tauxx)*dzt(lconb)
              tramix = trasum(n) / dztsum
              do lmix=lcona,lconb
                t(i,lmix,j,n,tauxx) = tramix
              enddo
            enddo
            goto 1306
          end if

!         4. test layer above lcona

1308      continue
          if ( lcona > 1 ) then
            l1 = lcona-1
            rholo(i,l1) = dens(t(i,l1,j,1,tauxx) -dn_to(lcona)
     &,                       t(i,l1,j,2,tauxx)-dn_so(lcona),lcona)
            rhoup(i,lcona) = dens(t(i,lcona,j,1,tauxx)-dn_to(lcona)
     &,                       t(i,lcona,j,2,tauxx)-dn_so(lcona),lcona)
            if ( rholo(i,lcona-1) > rhoup(i,lcona) ) then
              lcona = lcona-1
              dztsum = dztsum + dzt(lcona)
              do n=1,nt
                trasum(n) = trasum(n) + t(i,lcona,j,n,tauxx)*dzt(lcona)
                tramix = trasum(n) / dztsum 
                do lmix=lcona,lconb
                  t(i,lmix,j,n,tauxx) = tramix
                enddo
              enddo
              goto 1306
            end if
          end if

!         5. remember the total number of levels mixed by convection
!            in this water column, as well as the ventilated column

          lctot(i,j) = lctot(i,j) + lconb - lcona + 1
          if ( lcona == 1 ) lcven(i,j) = lconb

!         6. resume search if step 3. and 4. have been passed and this
!            unstable part of the water column has thus been removed,
!            i.e. find further unstable areas further down the column

          if ( lconb == kcon ) goto 1310
          lcon(i) = lconb

1302      continue
          lcon(i) = lcon(i) + 1

          if (lcon(i) .eq. kcon) goto 1310

          if ( rholo(i,lcon(i)) <= rhoup(i,lcon(i)+1) ) goto 1302

          goto 1319

1310      continue
        enddo

       endif

c       do i=is,ie
c         lctot(i,j)=zt(max(1,min(km,lcven1(i))))/100.  +lcven(i,j)
c       enddo
      enddo ! jrow
      end subroutine convect
#endif




      subroutine Ktmix
 
!------------------------------------------------------------------------------
! 
!     Kraus-Turner mixing, wind-driven part only, 
!     Sterl & Kattenberg numerical scheme.
! 
!     Subroutine for GFDL-MOM, V. 2.1
! 
!     Doc. of equations used:           LaTeX-file "Ktmix.tex".
!     Doc. of changes in code:          ASCII-file "Changes_Ktmix.doc"
! 
!     original version:           M. Cox, 1986
!                   Joachim Dengg, jdengg@awi-bremerhaven.de
!     SPFLAME version : c.eden
!
!------------------------------------------------------------------------------
!
      use spflame_module
      implicit none

!     -------------------------------------------------------------------------
!     Specify some necessary parameters:
!     -------------------------------------------------------------------------

!     - Maximum level of penetration of wind mixing
!
!       (to avoid penetrative mixing to too large depths due to asymptotic
!        decay of TKE-profile, choose the depth where TKE has decayed to 
!        -say - 1 percent of its surface value:  zt(kmax) = z_scale*ln(0.01) )
 
      integer,parameter :: kmax = 15       ! 250m: 1% level for z_scale=50m

!     - Small number to prevent floating overflows in layers
!       that are already mixed

      real,parameter :: small1 = 1.e-20
!     -------------------------------------------------------------------------
!     Define model variables
!     -------------------------------------------------------------------------

      real ::       pe(is_pe:ie_pe,kmax)   ! potential energy of water column 
!                                    !                  down to level k
      real ::      dPE(is_pe:ie_pe,kmax)   ! difference of PE before and after
!                                    !                            mixing
      real ::      TKE(is_pe:ie_pe)        ! turbulent kinetic energy
!                                    !      available from wind
      real ::     dTKE(is_pe:ie_pe,kmax)   ! difference between available
!                                    ! TKE and energy necessary to
!                                    ! mix next level completely

      real ::   rho_kt(is_pe:ie_pe,kmax)   ! pot. density [g/cm**3], 
!                                    !       referenced to z=0
      real ::     mass(is_pe:ie_pe,kmax)   ! mass in level k
      real :: sum_mass(is_pe:ie_pe,kmax)   ! total mass down to level k

      real ::    delta(is_pe:ie_pe)        ! thickness of partial mixing [cm]
      real ::     frac(is_pe:ie_pe)        ! fract. of level thickness [%] 
      real ::  c1mfrac(is_pe:ie_pe)        ! (1.-frac)

      real ::    cmask(is_pe:ie_pe,kmax)   ! mask for completely mixed levels
      real ::    pmask(is_pe:ie_pe,kmax)   ! mask for partially  mixed levels
      real ::   cpmask(is_pe:ie_pe,kmax)   ! mask for completely or 
!                                    !          partially  mixed levels
!                                    ! (0.: do not mix, 1.: mix)

      real ::     ztot(is_pe:ie_pe)        ! total depth of completely 
!                                    !                   mixed layers
      real ::      zcp(is_pe:ie_pe)        ! total depth of completely and 
!                                    !         partially mixed layers
      real ::     rzcp(is_pe:ie_pe)        ! (1./zcp)
      real ::    rztot(is_pe:ie_pe)        ! (1./ztot)

      real ::  dml_new(is_pe:ie_pe)        ! new depth of mixed layer [cm]

      real ::      t_c(is_pe:ie_pe,nt)     ! tracer on completely mixed levels
      real ::     t_cp(is_pe:ie_pe,nt)     ! tracer on comp. or part. 
!                                    !                      mixed levels

      real :: fxa, fxb               ! auxiliary variables

      logical, save :: first_time = .true.

!     - depth of center of gravity of vertical T-boxes [cm]
!           (In MOM2, T-points do not necessarily coincide with
!            midpoints of layers.)
      real, allocatable, save :: zkt(:)
!
!     ustar(:,:)  =  ustar**3 [(cm/s)**3]
!     effwind     =  m1, coefficient of effectiveness of wind forcing
!     z_scale     =  vertical scale length of TKE dissipation [cm]
!     c2dtts      =  2*timestep for MOM leapfrog, 1*timestep for Euler
!     rho0        =  reference density: physically almost irrelevant, 
!                    but numerically useful for greater precision.
!                    (physical importance: computation of TKE)
!     dml         =  old depth of mixed layer (only used to calculate 
!                    dissipation of TKE with exponential length scale) [cm]
!
!     tmask(i,k)  =  MOM`s land mask; prevents mixing into topography
!                 =  (0., 1.) over "t" grid (land, ocean) points
!     t(k,nt)     =  tracers from MOM (nt=1: T, nt=2: S)

      integer is,ie,js,je,i,j,k,m

#ifdef partial_cell
c      would have to define 3D array here ... skip this
c      anyway   zw(k) -> zw(k-1)+dhwt(i,k,j)
c      anyway   zkt(k) = zw(k-1)+dhwt(i,k,j)-0.5*dht(i,k,j)
#else
      if (first_time) then
c       if (my_pe == 0) print*, 'Initialization for Mixed Layer ktmix'
       first_time = .false.
       allocate( zkt(km))
!      define new vertical axis for centre of gravity of each T-box
       do k=1, km
       enddo
       do k=1, km
        zkt(k) = zw(k)-0.5*dzt(k)
       enddo
      endif
#endif

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

!     -------------------------------------------------------------------------
!     Apply Ktmix to the whole memory window
!     -------------------------------------------------------------------------

      j_loop: do j=js, je

!     -------------------------------------------------------------------------
!     Compute TKE provided by the wind 
!     -------------------------------------------------------------------------

      fxa = effwind * rho0 * c2dt / grav  
      fxb = -1./z_scale

#ifdef testktmix
      do i=is, ie  
        if (ustar(i,j) == -1.e10 .and. tmask(i,1,j) == 1.)
     &       write (6,*) 'ALARM: i,j,tmask: ',
     &                                i,j,tmask(i,1,j)
      end do
#endif
      TKE=0.
      TKE(is:ie) = fxa * exp(fxb*dml(is:ie,j)) * ustar(is:ie,j)

!     -------------------------------------------------------------------------
!     Compute potential densities referenced to the surface
!     -------------------------------------------------------------------------

      rho_kt=0.
      do k=1, kmax
       call model_dens(t(is,k,j,1,taup1),t(is,k,j,2,taup1),
     &                 rho_kt(is,k),1,ie-is+1 
#ifdef partial_cell
     &                       ,ztp(is,1,j)
#endif
     &                       )
      end do

!     -------------------------------------------------------------------------
!     Compute cumulative mass and PE from surface to maximum depth of ML
!     -------------------------------------------------------------------------

!     ... to determine how many layers can be mixed completely.

!     (Note: all sums are calculated down to level kmax to
!            allow vectorization.)

      mass=0.
      do k=1, kmax
#ifdef partial_cell
        mass(:,k) = rho_kt(:,k)*dht(is_pe:ie_pe,k,j)
#else
        mass(:,k) = rho_kt(:,k)*dzt(k)             ! rho(k)*dz(k)
#endif
      end do
 
      sum_mass=0.;pe=0.
      sum_mass(:,1) = mass(:,1)                    ! rho(1)     *dz(1)
#ifdef partial_cell
      pe(:,1) = mass(:,1)*
     &       (dhwt(is_pe:ie_pe,1,j)-0.5*dht(is_pe:ie_pe,1,j)) 
#else
      pe(:,1) = mass(:,1)*zkt(1)             ! rho(1)*z(1)*dz(1)
#endif

!                                                  ! Sum(rho(k)*dz(k))
!                                                  ! Sum(rho(k)*z(k)*dz(k))
      do k=2, kmax
        sum_mass(:,k) = sum_mass(:,k-1) + mass(:,k)  
#ifdef partial_cell
              pe(:,k) =       pe(:,k-1) + mass(:,k)*
     &       (zw(k-1)+dhwt(is_pe:ie_pe,k,j)-0.5*dht(is_pe:ie_pe,k,j)) 
#else
              pe(:,k) =       pe(:,k-1) + mass(:,k)*zkt(k)  
#endif
      end do

!     -------------------------------------------------------------------------
!     Compute difference dTKE between available TKE 
!                                 and PE necessary to mix next level
!     -------------------------------------------------------------------------
!
!     del TKE = TKE - del PE(k) 

      dPe=0.
      do k=1, kmax
#ifdef partial_cell
        fxa = zw(max(1,k-1)); if (k==1) fxa=0.
        dPE(:,k) = pe(:,k) 
     &      -0.5*(fxa+dhwt(is_pe:ie_pe,k,j))*sum_mass(:,k)
#else
        dPE(:,k) = pe(:,k) - 0.5*zw(k)*sum_mass(:,k)
#endif
      end do

      dPE = max(dPE, 0.)                        ! limit dPE to positive 

      dTKE=0.
      do k=1, kmax
        dTKE(:,k) = TKE(:) - dPE(:,k)
      end do

!     -------------------------------------------------------------------------
!     Compute mask for layers which are to be completely mixed
!     -------------------------------------------------------------------------
!
!     Stop mixing at the depth where del TKE is smaller than 0!
!
!     (Note: even if the density is already homogenized, del PE(k)
!            may not be exactly equal to zero due to round-off errors 
!            in the computation of PE. Patch: Check density as well.)
!
!      Also: the computation starts at k=2, because even without any 
!            exterior forcing, the mixing mask for level k=1 will always 
!            be 1 and the partial mixing for level 2 will be turned 
!            on, because at level 1, dPE=0. This does not have any  
!            serious consequences, though, because the partial mixing 
!            depth of level 2 still depends on the energy actually  
!            available, which in this case is negligible.)
       
      pmask = 0.; cmask=0.
      cmask(is:ie,1:kmax) = tmask(is:ie,1:kmax,j)        ! defaults

      where (dTKE(:,2:kmax) < 0. .and.           ! eliminate spurious cases
     &       rho_kt(:,2:kmax)-rho_kt(:,1:kmax-1) == 0.)
         dTKE(:,2:kmax) = 0.
      end where

      where (dTKE(:,2:kmax) < 0.) cmask(:,2:kmax) = 0.

      do k=2, kmax                               ! fill spurious gaps
        cmask(:,k) = cmask(:,k-1)*cmask(:,k)
      end do

!     -------------------------------------------------------------------------
!     Compute mask for layers which are to be partially mixed
!     -------------------------------------------------------------------------

      do k=2, kmax
        pmask(is:ie,k) = (cmask(is:ie,k-1)-cmask(is:ie,k)) 
     &                   * tmask(is:ie,k,j)
      end do

      if (minval(pmask) < -tiny(pmask)) then              ! security check
        write (6,*) 'Ktmix probably received an unstable ',
     &                   'density profile at:'
        do k=2, kmax
          do i=is,ie
            if (pmask(i,k)==-1.)
     &           write (6,*) 'i, j, k:', i, j, k
          end do
        end do

#ifdef testktmix
        do i=is, ie                                  ! test
          if (minval(pmask(i,:)) < -tiny(pmask)) then
            do k=1, kmax
              write (6,*) 'i,k:  ', i, k
              write (6,*) 'rho:  ', rho_kt(i,k)
              write (6,*) 'PE:   ', pe(i,k), .5*zw(k)*sum_mass(i,k)
              write (6,*) 'dPE:  ', dpe(i,k)
              write (6,*) 'dTKE: ', dTKE(i,k)
              write (6,*) 'mask: ', cmask(i,k), pmask(i,k)
            end do
          end if
        end do
#endif

        write (6,*) 'make sure that convct2 is called!'
        write (6,*) 

        call halt_stop('ktmix ==> partial mixing mask negative!')

      end if

!     -------------------------------------------------------------------------
!     Compute mask for either complete or partial mixing
!     -------------------------------------------------------------------------

      cpmask = cmask + pmask

!     -------------------------------------------------------------------------
!     Find depth of penetration of mixing
!
!      ... using the scheme described in Sterl and Kattenberg, 1994,
!                                                JGR, 99, 1439-14157
!     -------------------------------------------------------------------------
!
!     frac = fraction =  thickness of sublayer to be partially mixed relative
!                        to thickness of corresponding model level [%]
!          = (left-over energy del TKE after mixing last level completely)
!          / (energy needed to completely mix the partially mixed level 
!             at the base of the mixed layer)


         frac = 0.
        delta = 0.
#ifdef partial_cell
      dml_new = dht(is_pe:ie_pe,1,j)
          zcp = dht(is_pe:ie_pe,1,j)
         ztot = dht(is_pe:ie_pe,1,j)
#else
      dml_new = dzt(1)
          zcp = dzt(1)
         ztot = dzt(1)
#endif

      do k=2, kmax
!                                      ! frac=0 for completely mixed levels
          frac =    frac +  pmask(:,k)*dTKE(:,k-1)/     
     &                                (dTKE(:,k-1)-dTKE(:,k)+small1)
#ifdef partial_cell
         delta =   delta +  pmask(:,k)*frac*dht(is_pe:ie_pe,k,j)
       dml_new = dml_new +  cmask(:,k)*dht(is_pe:ie_pe,k,j) 
     &                   + pmask(:,k)*delta
          ztot =    ztot +  cmask(:,k)*dht(is_pe:ie_pe,k,j)
           zcp =     zcp + cpmask(:,k)*dht(is_pe:ie_pe,k,j)
#else
         delta =   delta +  pmask(:,k)*frac*dzt(k)
       dml_new = dml_new +  cmask(:,k)*dzt(k) + pmask(:,k)*delta
          ztot =    ztot +  cmask(:,k)*dzt(k)
           zcp =     zcp + cpmask(:,k)*dzt(k)
#endif
      end do

!     -------------------------------------------------------------------------
!     Assign new mixed-layer depth 
!     -------------------------------------------------------------------------

      dml(is:ie,j) = dml_new(is:ie)

!     -------------------------------------------------------------------------
!     Mix tracers
!
!      (Note: this scheme ensures that the potential energy created
!             exactly equals the TKE put in)
!     -------------------------------------------------------------------------

      c1mfrac = 1.-frac                 ! Speed up calculation by
        rztot = 1./ztot                 ! pre-computing constants
         rzcp = 1./zcp

      t_cp=0.; t_c=0.
      m_loop: do m=1, nt


!       Calculate mean values of tracers:
!         for mixing of completely mixed layers only:          t_c
!         for mixing of completely and partially mixed layers: t_cp

#ifdef partial_cell
        t_cp(is:ie,m) = t(is:ie,1,j,m,taup1)*dht(is:ie,1,j)
         t_c(is:ie,m) = t(is:ie,1,j,m,taup1)*dht(is:ie,1,j)
#else
        t_cp(is:ie,m) = t(is:ie,1,j,m,taup1)*dzt(1)
         t_c(is:ie,m) = t(is:ie,1,j,m,taup1)*dzt(1)
#endif

        do k=2, kmax
#ifdef partial_cell
          t_cp(is:ie,m) = t_cp(is:ie,m)
     &                    + cpmask(is:ie,k)*t(is:ie,k,j,m,taup1)
     &                     *dht(is:ie,k,j)
           t_c(is:ie,m) = t_c(is:ie,m)
     &                    +  cmask(is:ie,k)*t(is:ie,k,j,m,taup1)
     &                     *dht(is:ie,k,j)
#else
          t_cp(is:ie,m) = t_cp(is:ie,m)
     &                    + cpmask(is:ie,k)*t(is:ie,k,j,m,taup1)*dzt(k)
           t_c(is:ie,m) = t_c(is:ie,m)
     &                    +  cmask(is:ie,k)*t(is:ie,k,j,m,taup1)*dzt(k)
#endif
        end do

        t_cp(:,m) = t_cp(:,m)*rzcp
         t_c(:,m) =  t_c(:,m)*rztot

!       Assign new tracer values to layers by 
!
!         (a) leaving unaffected layers unchanged,
!         (b) interpolating between t_c and t_cp 
!                for completely mixed layers,
!         (c) interpolating between the old tracer value t of this level 
!             and t_cp for the partially mixed layer.

        do k=1, kmax
          t(is:ie,k,j,m,taup1) =
     &      (1.-cpmask(is:ie,k))*t(is:ie,k,j,m,taup1)      ! (a)
     &    + cmask(is:ie,k)                                 ! (b)
     &            *(c1mfrac(is:ie)*t_c(is:ie,m)            ! (b)
     &             +frac(is:ie)*t_cp(is:ie,m))             ! (b)
     &    + pmask(is:ie,k)                                 ! (c)
     &            *(c1mfrac(is:ie)*t(is:ie,k,j,m,taup1)    ! (c)
     &             +frac(is:ie)*t_cp(is:ie,m))             ! (c)
        end do

!       (rely on calling routine to apply horizontal boundary 
!        conditions for tracers)

      end do m_loop
      end do j_loop
      end subroutine Ktmix



      subroutine tkemix_init()
      use spflame_module
      implicit none
      !
      !----------------------------------------------------------------
      ! Document parameters for TKE scheme in output
      !
      ! wrapping of parameters to the old MOM1/OPA code of A.Oschlies:
      !
      ! C_eps        = ediss
      ! C_tke        = efave
      ! C_kappa      = ediff
      ! tke_surf_fac = ebb
      ! tke_min      = emin
      ! tke_min_surf = emin0
      !----------------------------------------------------------------
      !
      if (my_pe==0) then
       print*,''
       print*,' Initialization of the Gaspar et al.',
     &          ' vertical mixing scheme' 
       print*,' Coefficient for the Kolmogoroff dissipation',C_eps
       print*,' Coefficient for vertical eddy viscosity', C_kappa
       print*,' Coefficient for vertical TKE viscosity', C_tke
       print*,' Minimum value of TKE', tke_min,' cm^2/s^2'
       print*,' Minimum value of TKE at the surface', tke_min_surf,
     &              ' cm^2/s^2'
       print*,' Coefficient for input of TKE at surface ',tke_surf_fac
       if (nmxl==0) then
        print*,' mixing length is bound by distance to surface/bottom'
       elseif (nmxl == 1) then
        print*,' mixing length is bound by vertical scale factor'
       elseif (nmxl == 2) then
        print*,' abs( d/dz(mixing length) ) is bounded by dzt'
       endif
       print*,' diffusivity is bounded by ',
     &          diff_cbt_back,diff_cbt_cut,' cm^2/s'
       print*,' viscosity is bounded by ',
     &          visc_cbu_back,visc_cbu_cut,' cm^2/s'
       print*,''
       print*,' done' 
       print*,''
      endif
      end subroutine tkemix_init



      subroutine tkemix()

!-----------------------------------------------------------------------
!     Compute the vertical eddy viscosity and diffusivity coefficients
!     using a 1.5 turbulent closure scheme after
!     Gaspar et al., jgr, 95, 1990 and Blanke and Delecluse, jpo, 1991
!
!     The time evolution of the Turbulent Kinetic Energy (TKE) is
!     computed from the following equation:
!
!     d(tke)/dt = C_tke d( K_m d(tke)/dz )/dz + K_m (d(U)/dz)**2
!               +g/rho_0 /P_rt K_m d(RHO)/dz   - C_eps tke**(2/3) / mxl
!
!           K_m =  C_kappa mxl tke**(1/2) 
!           mxl = 1/(sqrt(tke)*N)
!     with the boundary conditions:
!     tke = max( tke_min_surf,tke_surf_fac (TAUX**2 + TAUY**2)**(1/2) 
!                                                     at the surface
!     tke = tke_min                                   at the bottom
!   
!     The time differencing is implicit for the vertical diffusion
!     term, implicit and linearized for the Kolmogoroff dissipation
!     term, and explicit forward for both buoyancy and dynamic
!     production terms.
!
!     The mixing length is the usual bouyancy length scale
!               mxl = 1/(sqrt(tke)*N)
!     where N is the Brunt-Vaisala frequency.
!     Three different physical limitations can be imposed on the
!     mixing length depending on NMXL:
!           NMXL=0 : mxl bounded by the distance to surface/bottom
!           NMXL=1 : mxl bounded by the vertical scale factor
!           NMXL=2 : vertical derivative of mxl bounded by 1.
!
!     The eddy diffusion coefficients are deduced from TKE using
!     the following equations:
!           visc_cbu = max( visc_cbu_back, C_kappa mxl tke**(1/2) )
!           diff_cbt = max( diff_cbt_back, visc_cbu/P_rt )
!           K_m = max( diff_cbu_back, visc_cbu)
!
!     diff_cbt and visc_cbu are horizontally averaged to avoid numerical
!     instabilities, K_m is not averaged.
!
!     Code was taken from OPA and/or MOM1
!     original author B. Blanke, G. Madec, ... A. Oschlies 
!     adapted for SPFLAME:   C. Eden
!-----------------------------------------------------------------------
      use spflame_module
      implicit none

#if defined  IRIX_host || defined ALPHA_host || defined LINUX_host ||defined SUN_host
!     single precision calculation
      real,parameter :: eps   =1.e-15
      real,parameter :: eps2  =1.e-15
#else
!     double precision calculation
      real,parameter :: eps   =1.e-23
      real,parameter :: eps2  =1.e-20
#endif
      real :: K_m(is_pe-1:ie_pe+1,0:km,js_pe-1:je_pe+1) 
      real :: mxl(is_pe:ie_pe,km,js_pe:je_pe)  
      real :: shear(is_pe:ie_pe,km,js_pe:je_pe)
      real :: Nsqr(is_pe:ie_pe,km,js_pe:je_pe)  
      real :: wkx(is_pe:ie_pe,km), wky(is_pe:ie_pe,km) 
      real, dimension(is_pe:ie_pe,km) :: wi,ws,wd,wy,wx,wt
      real, dimension(is_pe:ie_pe)    :: ro1,ro2
      real :: w_surf(is_pe:ie_pe)   
      real,parameter :: zmlmin=1.e-8*100.
      real :: fxa, P_rt, Ri
      integer :: is,ie,js,je,i,j,k

!     -------------------------------------------------------------------------
!     limit the longitude indices
!     -------------------------------------------------------------------------
      is=max(is_pe,2); ie=min(ie_pe,imt-1)
      js=max(2,js_pe); je=min(je_pe,jmt-1)
!
!     1. Mixing length
!     ----------------
!     1.1 Brunt-Vaisala frequency
!
      Nsqr=0.
      do j=js,je
       do k=1,km-1
        call model_dens(t(is,k  ,j,1,tau),t(is,k  ,j,2,tau),
     &                  ro1(is),k, ie-is+1
#ifdef partial_cell
     &                       ,ztp(is,k,j)
#endif
     &                       )
        call model_dens(t(is,k+1,j,1,tau),t(is,k+1,j,2,tau),
     &                  ro2(is),k, ie-is+1
#ifdef partial_cell
     &                       ,ztp(is,k,j)
#endif
     &                       )
        do i=is,ie
#ifndef partial_cell
         fxa = -grav/rho0*dzwr(k)
#endif
#ifdef partial_cell
         fxa = -grav/(rho0*dhwt(i,k,j))
#endif
         Nsqr(i,k,j) =fxa*(ro1(i)-ro2(i))*tmask(i,k+1,j)
        enddo
       enddo
       Nsqr(is:ie,km,j)=0.
      enddo

      do j=js,je
!
!     1.2 bouyancy length scale
!
       do k=1,km-1
        do i=is,ie
          mxl(i,k,j) = sqrt( 2.*eke(i,k,j,tau)/max(Nsqr(i,k,j),eps))
     &                                       *tmask(i,k+1,j)
          mxl(i,k,j) = max(mxl(i,k,j),zmlmin)
        enddo
       enddo
       mxl(is:ie,km,j) = zmlmin
!
!      1.3 physical limits for mixing length
!
       if (nmxl==0) then
!
!       bounded by the distance to surface/bottom
!
        do k=1,km-1
          do i=is,ie
#ifdef partial_cell
            fxa = zw(max(1,k-1)); if (k==1) fxa=0.
            fxa = fxa + dhwt(i,k,j)
            mxl(i,k,j) = min(fxa,mxl(i,k,j),htp(i,j)-fxa)
#else
            mxl(i,k,j) = min(zw(k),mxl(i,k,j),zw(kmt(i,j))-zw(k))
#endif
            mxl(i,k,j) = max(mxl(i,k,j),zmlmin)
          enddo
        enddo

       elseif (nmxl == 1) then
!
!       bounded by the vertical scale factor
!
        do k=1,km-1
          do i=is,ie
#ifdef partial_cell
            mxl(i,k,j) = min(mxl(i,k,j),dhwt(i,k,j))
#else
            mxl(i,k,j) = min(mxl(i,k,j),dzw(k))
#endif
          enddo
        enddo

       elseif (nmxl == 2) then
!
!      abs( d/dz(EXML) ) bounded by dzt
!
!      top to bottom
!
        do i=is,ie
#ifdef partial_cell
          mxl(i,1,j) = min(dht(i,1,j),mxl(i,1,j))
#else
          mxl(i,1,j) = min(zw(1),mxl(i,1,j))
#endif
        enddo
        do k=2,km-1
          do i=is,ie
#ifdef partial_cell
            mxl(i,k,j) = min(mxl(i,k-1,j)+dht(i,k,j),mxl(i,k,j))
#else
            mxl(i,k,j) = min(mxl(i,k-1,j)+dzt(k),mxl(i,k,j))
#endif
          enddo
        enddo
!
!       bottom to top
!
        do k=km-1,1,-1
          do i=is,ie
#ifdef partial_cell
            mxl(i,k,j) = min(mxl(i,k+1,j)+dht(i,k+1,j),mxl(i,k,j))
#else
            mxl(i,k,j) = min(mxl(i,k+1,j)+dzt(k+1),mxl(i,k,j))
#endif
          enddo
        enddo

       endif

      enddo ! j
!
!     2. vertical eddy diffusivity and viscosity
!     ----------------------------------------
!
      K_m=0.
      do j=js,je
       do k=1,km-1
        do i=is,ie
          fxa = C_kappa*mxl(i,k,j)*sqrt(eke(i,k,j,tau))
          K_m(i,k,j) = max(fxa,diff_cbt_back)*tmask(i,k+1,j)
        enddo
       enddo
!
!      set minimum mixing length at surface
!
       K_m(is:ie,0,j) = K_m(is:ie,1,j)
      enddo

      call border_exchg(K_m,km+1,1)
      call set_cyclic(K_m,km+1,1)
!
!     2.2 eddy viscosity: horizontal average
!
      do j=js,je
       do k=1,km-1
        do i=is,ie ! imt-1
          visc_cbu(i,k,j) = (K_m(i,k,j)+K_m(i+1,k,j)
     &                      +K_m(i,k,j+1)+K_m(i+1,k,j+1))/
     &  max(1.,tmask(i,k+1,j)+tmask(i,k+1,j+1)+
     &         tmask(i+1,k+1,j)+tmask(i+1,k+1,j+1) )
          visc_cbu(i,k,j) = max(visc_cbu_back,visc_cbu(i,k,j))
          visc_cbu(i,k,j) = min(visc_cbu(i,k,j),visc_cbu_cut)
     &                      *umask(i,k+1,j)
        enddo
       enddo
       visc_cbu(is:ie,km,j)= 0.
      enddo
!
!     2.3 eddy diffusivity: horizontal average
!
      do j=js,je
       do k=1,km-1
        do i=is,ie 
          diff_cbt(i,k,j) = ( K_m(i-1,k,j)+K_m(i,k,j)+K_m(i+1,k,j)
     &             + K_m(i-1,k,j+1)+K_m(i,k,j+1)+K_m(i+1,k,j+1)
     &             + K_m(i-1,k,j-1)+K_m(i,k,j-1)+K_m(i+1,k,j-1) )/
     &  max(1.,tmask(i-1,k+1,j)+tmask(i,k+1,j)+tmask(i+1,k+1,j)
     &     +tmask(i-1,k+1,j+1)+tmask(i,k+1,j+1)+tmask(i+1,k+1,j+1)
     &     +tmask(i-1,k+1,j-1)+tmask(i,k+1,j-1)+tmask(i+1,k+1,j-1) )
          diff_cbt(i,k,j) = max(diff_cbt_back,diff_cbt(i,k,j)) 
          diff_cbt(i,k,j) = min(diff_cbt(i,k,j),diff_cbt_cut)
     &                     *tmask(i,k+1,j)
        enddo
       enddo
       diff_cbt(is:ie,km,j)= 0.
      enddo
!
!     3. square of vertical shear of horizontal velocity at T-points
!     --------------------------------------------------------------
!
      do j=js,je
       do k=1,km
        do i=is,ie
          wkx(i,k) = 0.25*(u(i-1,k,j-1,1,tau) + u(i,k,j-1,1,tau)
     &                   + u(i-1,k,j  ,1,tau) + u(i,k,j  ,1,tau))
          wky(i,k) = 0.25*(u(i-1,k,j-1,2,tau) + u(i,k,j-1,2,tau)
     &                   + u(i-1,k,j  ,2,tau) + u(i,k,j  ,2,tau))
        enddo
       enddo
       do k=1,km-1
        do i=is,ie
#ifdef partial_cell
          shear(i,k,j) = ( (wkx(i,k)-wkx(i,k+1))/dhwt(i,k,j) )**2
     &                  +( (wky(i,k)-wky(i,k+1))/dhwt(i,k,j) )**2
#else
          shear(i,k,j) = ( (wkx(i,k)-wkx(i,k+1))*dzwr(k) )**2
     &                  +( (wky(i,k)-wky(i,k+1))*dzwr(k) )**2
#endif
     &                 *tmask(i,k+1,j)
        enddo
       enddo
      enddo
!
!     5.  multiplying Nsqr and VVC by the inverse of the Prandtl number.
!     ------------------------------------------------------------------
!     The Prandtl number is given by an empirical function of the
!     Richardson number and is bounded by 10.
!
      do j=js,je
       do k=1,km-1
        do i=is,ie
          Ri  = max(Nsqr(i,k,j),eps2)/(shear(i,k,j)+eps2)
          P_rt = 1./max(1.,5.*Ri)
          P_rt = max(0.1,P_rt)
          Nsqr(i,k,j)      = P_rt*Nsqr(i,k,j)
          diff_cbt(i,k,j)  = max( diff_cbt_back,P_rt*diff_cbt(i,k,j) )
     &                      *tmask(i,k+1,j)
        enddo
       enddo
      enddo
!
!     6. surface boundary condition on TKE
!     ------------------------------------
!
      do j=js,je 

       do i=is,ie
        fxa = ustar(i,j) *tke_surf_fac*tmask(i,k,j)
#ifdef partial_cell
        w_surf(i) = eke(i,1,j,taum1)+dht(i,1,j)
#else
        w_surf(i) = eke(i,1,j,taum1)+dzt(1)
#endif
     &     /(C_tke*max(diff_cbt_back ,K_m(i,1,j)) )*fxa
        w_surf(i) = max(w_surf(i),tke_min_surf)
       enddo
!
!      7. implicit solution for TKE
!      ----------------------------
!      computation from level 1 to km-1 (tke(0) is given above,
!      and tke(km)=0)
!      The coefficients of the tridiagonal matrix are in WI, WS,
!      and WD. The right hand side (explicit part) is in WY.
!      The result is given in WX.
!
       do k=1,km-1
        do i=is,ie 
         fxa = -0.5*c2dt*C_tke
        wi(i,k) =fxa*( K_m(i,k  ,j)+K_m(i,k-1,j)) 
#ifdef partial_cell
     &                /(dht(i,k,j)*dhwt(i,k,j)) *tmask(i,k+1,j)
#else
     &                /(dzt(k  )*dzw(k)) *tmask(i,k+1,j)
#endif
        ws(i,k) =fxa*( K_m(i,k+1,j)+K_m(i,k  ,j) )
#ifdef partial_cell
     &                /(dht(i,k+1,j)*dhwt(i,k,j)) *tmask(i,k+1,j)
#else
     &                /(dzt(k+1)*dzw(k)) *tmask(i,k+1,j)
#endif
        fxa=sqrt(eke(i,k,j,tau))/mxl(i,k,j)
        wd(i,k) = 1. - wi(i,k) - ws(i,k) +1.5*c2dt*C_eps*fxa
        wy(i,k) = eke(i,k,j,taum1) + 0.5*c2dt*C_eps*fxa*eke(i,k,j,tau)
     &                + c2dt*K_m(i,k,j)*( shear(i,k,j)-Nsqr(i,k,j) )
        enddo
       enddo
!
!      7.1 surface boundary conditions on TKE
!
       wy(is:ie,1) = wy(is:ie,1) - wi(is:ie,1)*w_surf(is:ie)
!
!      7.2 matrix inversion
!
!      Solve M*X = Y , where M is a tri-diagonal K*K matrix (K=km-1)
!
!      ( wd1 ws1  0   0   0  )( wx1 ) ( wy1 )
!      ( wi2 wd2 ws2  0   0  )( wx2 ) ( wy2 )
!      (  0  wi3 wd3 ws3  0  )( wx3 )=( wy3 )
!      (       ...           )( ... ) ( ... )
!      (  0   0   0  wiK wdK )( wxK ) ( wyK )
!
!      M is decomposed into the product of an upper and a lower
!      diagonal matrix.
!
       wt(is:ie,1) = wd(is:ie,1)
       do k=2,km-1
        wt(is:ie,k)= wd(is:ie,k)-wi(is:ie,k)*ws(is:ie,k-1)/wt(is:ie,k-1)
       enddo
       wd(is:ie,1) = wy(is:ie,1)
       do k=2,km-1
        wd(is:ie,k)= wy(is:ie,k)-wi(is:ie,k)/wt(is:ie,k-1)*wd(is:ie,k-1)
       enddo
       wx(is:ie,km-1) = wd(is:ie,km-1)/wt(is:ie,km-1)
       do k= km-2,1,-1
        wx(is:ie,k)= (wd(is:ie,k)-ws(is:ie,k)*wx(is:ie,k+1))/wt(is:ie,k)
       enddo
!
!      7.3 minimum value of TKE : tke_min
!
       do k=1,km-1
        do i=is,ie
          eke(i,k,j,taup1) = max(wx(i,k),tke_min)
        enddo
       enddo
       eke(is:ie,km,j,taup1) = 0.

      enddo ! j

      end subroutine tkemix
