c-----------------------------------------------------------------------
      subroutine kvmix (kpp_type,npt,nz,ntrac,h,nzi,fcor,u,v,t,s,
     *                   taux,tauy,wnsp,sswt,hflx,sflx,trflx,
     *                   sourct,sourctr,vvc,vdc,hblt,zmld)

c  input : 
c       npt     number of grid points on top level
c       nz      total number of layers
c       ntrac   total number of passive tracers
c       nzi     number of layers at grid point ipt
c       h       layer thicknesses
c       u       zonal velocity
c       v       meridional velocity
c       t       potential temperature
c       s       salinity
c       taux    surface zonal momentum  flux
c       tauy    surface meridional mom. flux 
c       wnsp    surface windspeed
c       hflx    surface heat flux
c       sflx    surface salinity flux
c       trflx   surface tracer flux
c       sswt    solar sw surface heat flux 
c       fcor    coriolis parameter
c  output:
c       vvc     vertical viscosity
c       vdc     vertical diffusivity for (temperature,salinity)
c       sourct  "source" terms for (temperature,salinity,tracers)
c       hblt    boundary layer depth
c       

# include "comm_kpp.h"
       parameter (GRAVTY = 9.8) 
       parameter (km   = 100)
       parameter (kmp1 = km + 1)
       parameter (kmm1 = km - 1)
       parameter (CD = 1.3e-3)
       parameter (CD2_over_rho = 1.123e-3)  ! sqrt(CD/rho(ocean))
       parameter (CD2o4 = 4.492e-3)  ! 4*sqrt(CD/rho(ocean))

       dimension nzi(npt),u(npt,nz),v(npt,nz),h(npt,nz),taux(npt),tauy(npt),
     *    wnsp(npt),sswt(npt),fcor(npt),t(npt,nz),s(npt,nz),
     *    hflx(npt),sflx(npt)
     *   ,trflx(npt,ntrac),zmld(npt)
       dimension vvc(npt,nz),vdc(npt,nz,2),hblt(npt),
     *    sourct(npt,nz,2),sourctr(npt,nz,ntrac)


c=======================================================================
c
c     written  by  : jan morzel, august  11, 1994
c     modified by  : jan morzel, january 25, 1995 : "dVsq" and 1d code
c
c     compute vertical mixing coefficients based on the k-profile
c     and oceanic planetary boundary layer scheme by large & mcwilliams.
c
c     summary:
c     - compute interior mixing everywhere:
c       interior mixing gets computed at all interfaces due to constant
c       internal wave background activity ("fkpm" and "fkph"), which
c       is enhanced in places of static instability (local richardson
c       number < 0).
c       Additionally, mixing can be enhanced by adding contribution due
c       to shear instability which is a function of the local richardson
c       number 
c     - k-profile parameterization (kpp) scheme in the boundary layer:
c 
c       a.boundary layer depth:
c         at every gridpoint the depth of the oceanic boundary layer 
c         ("hbl") gets computed by evaluating bulk richardson numbers.
c       b.boundary layer mixing:
c         within the boundary layer, above hbl, vertical mixing is 
c         determined by turbulent surface fluxes, and interior mixing at
c         the lower boundary, i.e. at hbl.
c     
c     this subroutine provides the interface between mom and the 
c     subroutine "kppmix", where boundary layer depth, vertical 
c     viscosity, vertical diffusivity, and counter gradient term (ghat)
c     are computed slabwise.
c     note: subroutine "kppmix" uses m-k-s units.
c
c     time level:
c     input tracer and velocity profiles are evaluated at time level 
c     tau, surface fluxes come from tau or tau-1.
c
c     grid option:
c     in this "1-grid" implementation, diffusivity and viscosity
c     profiles are computed on the "t-grid" (by using velocity shear
c     profiles averaged from the "u-grid" onto the "t-grid"; note, that
c     the averaging includes zero values on coastal and seafloor grid 
c     points).  viscosity on the "u-grid" is computed by averaging the 
c     "t-grid" viscosity values onto the "u-grid".
c
c     vertical grid:
c     mixing coefficients get evaluated at the bottom of the lowest 
c     layer, i.e. at depth zw(km).  these values are only useful when 
c     the model ocean domain does not include the entire ocean down to
c     the seafloor ("upperocean" setup) and allows flux through the
c     bottom of the domain.  for full-depth runs, these mixing 
c     cefficients are being zeroed out before leaving this subroutine.
c             upperocean disabled -nhn
c
c
c     local arrays:
c
c     zgrid (0:kmp1)= vert.levels of tracers on jp (<=0)         (m)
c     hwide (0:kmp1)= layer thicknesses on jp      (>=0)         (m)
c
c     ustar        = sfc friction velocity                    (m/s)
c     Bo           = sfc turbulent buoyancy forcing         (m2/s3)
c     Bosol        = sfc radiative buoyancy forcing         (m2/s3)
c
c     Shsq  (km)    = (local velocity shear)**2
c                                 at interfaces , for "ri_iwmix" (m2/s2)
c     dVsq  (km)    = (velocity shear re sfc)**2
c                                 at grid levels, for "bldepth"  (m2/s2)
c     rho          = density of surface layer               (kg/m3)
c     talpha(kmp1)  = d(rho)/ d(potential temp.)           (kg/m3/C)
c     sbeta (kmp1)  = d(rho)/ d(salinity)                (kg/m3/PSU)
c
c     dbloc (km)    = local delta buoyancy at interfaces,
c                         for "ri_iwmix" & "bldepth"              (m/s2)
c     dbsfc (km)    = delta buoyancy w/ respect to
c                         sfc, at grid levels, for Ritop          (m/s2)
c     Ritop (km)    = numerator of bulk richardson 
c                         number  at grid levels, for "bldepth"    
c
c     vvcc  (0:kmp1)= vertical viscosity on current row       (m2/s)

c     vvcp  (0:kmp1)= vertical viscosity on next row          (m2/s)
c     vdctp (0:kmp1)= vert. diff on next row for temperature  (m2/s)
c     vdcsp (0:kmp1)= vert. diff on next row for salt&tracers (m2/s)
c     ghatp (km)    = nonlocal transport coeff. on next row   (s/m2)
c
c     hbltp        = bdry lyr depth on "t-grid" on next row     (m)
c
c=======================================================================
c
c-----------------------------------------------------------------------
c     local arrays
c-----------------------------------------------------------------------
c
      real zgrid(0:kmp1),hwide(0:kmp1),
     $     ustar,bo,bosol,temp(km),salt(km),
     $     Shsq(km),dVsq(km),
     $     rho,talpha(kmp1),sbeta(kmp1),
     $     dbloc(km),dbsfc(km),Ritop(km)
 
         
      common /kvmixc/
     $     vvcc(0:kmp1),vvcp(0:kmp1),
     $     vdctp(0:kmp1),vdcsp(0:kmp1),
     $     ghatp(km),hbltp
 
 
c
c-----------------------------------------------------------------------
c     set jerlov water type (needed for solar penetration, otherwise
c     dummy array)
c-----------------------------------------------------------------------
c
      jwtype = jwatertype
c
      do ipt = 1, npt
         nzl = nzi(ipt)
         nzp = nzl + 1
         nzm = nzl - 1

c-----------------------------------------------------------------------
c     prepare input arrays for subroutine "kppmix" to compute 
c     viscosity and diffusivity and ghatp on next row.
c     all input arrays need to be in m-k-s units.
c
c     note: for the computation of the bulk richardson number in the
c     "bldepth" subroutine, gradients of velocity and buoyancy are 
c     required at every depth. in the case of very fine vertical grids
c     (thickness of top layer < 2m), the surface reference depth must
c     be set to zref=epsilon/2*zgrid(i,k), and the reference value
c     of velocity and buoyancy must be computed as vertical average 
c     between the surface and 2*zref.  in the case of coarse vertical 
c     grids zref is zgrid(i,1)/2., and the surface reference value is
c     simply the surface value at zgrid(i,1).
c-----------------------------------------------------------------------
c
c     vertical grid
c     -------------	
c
         zgrid(0)    = epsln
         hwide(0)    = epsln
         
         do 1100 k=1,nzl
               hwide(k) = h(ipt,k)
               zgrid(k) = zgrid(k-1) - 0.5*(hwide(k)+hwide(k-1))
               temp(k)  = t(ipt,k)
               salt(k)  = s(ipt,k)
 1100    continue

         k = nzp
         hwide(k) = epsln
         zgrid(k) = zgrid(k-1) - 0.5*hwide(k-1)
         
c
c     velocity shear
c     --------------
c     get velocity shear squared, averaged from "u-grid" 
c     onto "w-grid" (in (m/s)**2):

c     dVsq(k)=(Uref-U(k))**2+(Vref-V(k))**2      at grid levels
c     Shsq(k)=(U(k)-U(k+1))**2+(V(k)-V(k+1))**2  at interfaces  
c
c
         do 1200 k=1,nzm
            kp1 = k+1
            dVsq(k) =   (u(ipt,1)-u(ipt,k))**2
     $                + (v(ipt,1)-v(ipt,k))**2
            Shsq(k) =   (u(ipt,k)-u(ipt,kp1))**2
     $                + (v(ipt,k)-v(ipt,kp1))**2
 1200    continue
         
         k=nzl
         
         dVsq(k) =   (u(ipt,1)-u(ipt,k))**2
     $             + (v(ipt,1)-v(ipt,k))**2
         Shsq(k) = 0.0
         
c
c
c-----------------------------------------------------------------------
c
c     density related quantities
c     --------------------------
c     based on mom equation of state, which computes normalized
c     density  drho{k,kr} for layer k with respect to reference layer
c     kr, i.e. drho{T(k),S(k), tr(kr), sr(kr), zt(kr)}, as the 
c     difference of drho=rho{k,kr}-rhor(kr):
c 
c     density of surface layer                                   (kg/m3)
c           rho   = drho{T(1),S(1),zt(1)} + 1000. + rhor(zt(1))
c     local buoyancy gradient at km interfaces:                   (m/s2)
c           dbloc = g/rho{k+1,k+1} * [ drho{k,k+1}-drho{k+1,k+1} ]
c     buoyancy difference with respect to "zref", i.e. the surface(m/s2)
c           dbsfc = g * [ drho{1,k}/rho{1,k} - drho{k,k}/rho{k,k} ]
c     thermal expansion coefficient without 1/rho factor       (kg/m3/C)
c           talpha= d(rho{k,k})/d(T(k))
c     salt expansion coefficient without 1/rho factor        (kg/m3/PSU)
c           sbeta = d(rho{k,k})/d(S(k))
c-----------------------------------------------------------------------
c
         
         call statekpp (nzl, temp, salt, zgrid,
     $                  rho,  dbloc,  dbsfc, talpha, sbeta  )
         
c
c     numerator of bulk richardson number on grid levels
c     --------------------------------------------------
c     note: land and ocean bottom values need to be set to zero
c     so that the subroutine "bldepth" works correctly
c
         do 1420 k=1,nzl
            Ritop(k) = (zgrid(1)-zgrid(k)) * dbsfc(k) 
 1420    continue   

c
c     kinematic surface fluxes on the "t-grid" in m-k-s units:
c     --------------------------------------------------------
c
c     wusfc = kinematic zonal velocity sfc flux on jp            (m2/s2)
c           = -taux(i,jp1)
c     wvsfc = kinematic merid.velocity ..                        (m2/s2)
c           = -tauy(i,jp1)
c     wtsfc = kinematic temperature    ..                        (C*m/s)
c           = -hflx(i,jp1)
c     wssfc = kinematic salinity       ..                      (PSU*m/s)
c           = -sflx(i,jp1)
c     wtsol = kinematic solar temp.    ..                        (C*m/s)
c           = -sswt(i,jp1)*??
c
c     friction velocity, turbulent and radiative sfc buoyancy forcing
c     ---------------------------------------------------------------
c
c     (note: g is gravity in m-k-s units)
c
c     ustar(i) = sqrt(sqrt(wusfc(i)**2 + wvsfc(i)**2))             (m/s)
c     Bo(i)    = -g*(talpha(i,1)*wtsfc(i)-sbeta(i,1)*wssfc(i))/rho(i)
c     Bosol(i) = -g*(talpha(i,1)*wtsol(i)                    )/rho(i)
c                                                                (m2/s3)
c


         if (kpp_type.eq.0) then
            ustar = CD2_over_rho*wnsp(ipt)
         elseif (kpp_type.eq.1) then
            ustar = max(CD2o4,sqrt(sqrt( taux(ipt)**2 + tauy(ipt)**2 )))
         elseif (kpp_type.eq.2) then
            ustar = CD2o4
         else
            ustar = CD2o4*1e-8
         endif
         
c
         Bo    =  -GRAVTY*( talpha(1) * hflx(ipt) 
     $        - sbeta (1) * sflx(ipt) )
     $        / rho
c
         Bosol =  -GRAVTY*  talpha(1) * sswt(ipt) 
     $        / rho

c
c     coriolis parameter, fcor is defined on the "u-grid"
c
         Coriol = fcor(ipt)
c
c-----------------------------------------------------------------------
c     solve for viscosity, diffusivity, and ghat 
c-----------------------------------------------------------------------
c
         call kppmix(nz, zgrid , hwide , nzl , Shsq  , dVsq ,
     $        ustar , Bo    , Bosol , 
     $        dbloc , Ritop , Coriol, jwtype ,
     $        vvcp  , vdcsp , vdctp , ghatp  , hbltp ) 
         
         
         do 1700 k=1,nzl
            vvc(ipt,k)   = vvcp(k)
            vdc(ipt,k,1) = vdctp(k)
            vdc(ipt,k,2) = vdcsp(k)
 1700    continue
         
         hblt(ipt) = hbltp
         zmld(ipt) = max(zmld(ipt),hbltp)
         
c-----------------------------------------------------------------------
c     add ghatp term to "source" term
c-----------------------------------------------------------------------
         k =1
         sourct(ipt,k,1) = hflx(ipt)/h(ipt,k)*(-vdc(ipt,k,1)*ghatp(k) )
         sourct(ipt,k,2) = sflx(ipt)/h(ipt,k)*(-vdc(ipt,k,2)*ghatp(k) )
         
         do k=2,nzl
            sourct(ipt,k,1)= hflx(ipt)/h(ipt,k)*( vdc(ipt,k-1,1)*ghatp(k-1) 
     $                                           -vdc(ipt,k  ,1)*ghatp(k  ) )
            sourct(ipt,k,2)= sflx(ipt)/h(ipt,k)*( vdc(ipt,k-1,2)*ghatp(k-1) 
     $                                           -vdc(ipt,k  ,2)*ghatp(k  ) )
         enddo
 
          do 500 m=1,ntrac
            k =1
            sourctr(ipt,k,m) = trflx(ipt,m)/h(ipt,k)*(-vdc(ipt,k,2)*ghatp(k))
            
            do 480 k=2,nzl
               sourctr(ipt,k,m) = trflx(ipt,m)/h(ipt,k)
     $             *( vdc(ipt,k-1,2)*ghatp(k-1) - vdc(ipt,k  ,2)*ghatp(k  ) )
  480       continue  
  500    continue


      enddo
c
c
c
      return
      end
c
c

c-----------------------------------------------------------------------
      subroutine kppmix(nz, zgrid , hwide , kmtj  , Shsq   , dVsq  ,
     $                  ustar , Bo    , Bosol , 
     $                  dbloc , Ritop , Coriol, jwtype ,
     $                  visc  , difs  , dift  , ghats  , hbl 
     $                                                      ) 

c
c-----------------------------------------------------------------------
c
c     Main driver subroutine for kpp vertical mixing scheme and 
c     interface to greater ocean model
c
c     written  by: bill large,    june  6, 1994
c     modified by: jan morzel,    june 30, 1994
c                  bill large,  august 11, 1994
c                  bill large, january 25, 1995 : "dVsq" and 1d code
c
c-----------------------------------------------------------------------
c     
      parameter (mdiff = 3)  ! number of diffusivities for local arrays
c
c input
      real zgrid(0:nz+1)     ! vertical grid (<= 0)            (m)
      real hwide(0:nz+1)     ! layer thicknesses               (m)
      integer kmtj           ! number of vertical layers on this row
      real Shsq(nz)          ! (local velocity shear)^2        (m/s)^2
      real dVsq(nz)          ! (velocity shear re sfc)^2       (m/s)^2
      real ustar             ! surface friction velocity       (m/s)
      real Bo                ! surface turbulent buoy. forcing (m^2/s^3)
      real Bosol             ! radiative buoyancy forcing      (m^2/s^3)
      real dbloc(nz)         ! local delta buoyancy            (m/s^2)
c                              across interfaces
      real Ritop(nz)     ! numerator of bulk Richardson 
c                              Number:                         (m/s)^2
c          Ritop = (-z - -zref)* delta buoyancy w/ respect to sfc
      real Coriol            ! Coriolis parameter              (1/s)
      integer jwtype         ! Jerlov water type               (1 to 5)
c                              (only used ifdef "sdflxsw")
c
c output
      real visc(0:nz+1)  ! vertical viscosity coefficient  (m^2/s)
      real difs(0:nz+1)  ! vertical scalar diffusivity     (m^2/s)
      real dift(0:nz+1)  ! vertical temperature diffusivity(m^2/s)
      real ghats(nz)     ! nonlocal transport              (s/m^2)
      real hbl           ! boundary layer depth (m)
c
c local
      real bfsfc        ! surface buoyancy forcing        (m^2/s^3)
      real ws           ! momentum velocity scale
      real wm           ! scalar   velocity scale 
      real caseA        ! = 1 in case A; =0 in case B
      real stable       ! = 1 in stable forcing; =0 in unstable
      real dkm1(mdiff)   ! boundary layer difs at kbl-1 level
      real gat1(mdiff)   ! shape function at sigma=1
      real dat1(mdiff)   ! derivative of shape function at sigma=1
      real blmc(mdiff,100)! boundary layer mixing coefficients
      real sigma        ! normalized depth (d / hbl)
      real Rib(2)        ! bulk Richardson number
c1d                            output for 1d code:
      integer kbl       ! index of first grid level below hbl
 
c-----------------------------------------------------------------------
c compute interior mixing coefficients everywhere, due to constant 
c internal wave activity, static instability, and local shear 
c instability.
c Zero surface values
c-----------------------------------------------------------------------
c

      kmtjp = kmtj+1

      call ri_iwmix (kmtj, kmtjp, Shsq, dbloc, zgrid , visc, difs , dift  )
c
c
c-----------------------------------------------------------------------
c Zero values at seafloor and below for for blmix
c-----------------------------------------------------------------------
c
      do 210 ki=1,nz
         ktemp = max  (kmtj,1)
         wait  = 1.0 - MIN( ki/ktemp, 1 )
         visc(ki) = wait * visc(ki) 
         difs(ki) = wait * difs(ki)
         dift(ki) = wait * dift(ki)
 210  continue

c 
c-----------------------------------------------------------------------
c compute boundary layer mixing coefficients:
c
c diagnose the new boundary layer depth
c-----------------------------------------------------------------------
c
      call bldepth (kmtj, kmtjp, zgrid, hwide, dVsq  ,
     $              dbloc, Ritop, ustar , Bo   , Bosol, Coriol, jwtype,
     $              hbl  , bfsfc, stable, caseA, kbl  ,
     $              Rib  , sigma, wm    , ws   )
c 
c-----------------------------------------------------------------------
c compute boundary layer diffusivities
c-----------------------------------------------------------------------
c
      call blmix    (kmtj   , kmtjp , mdiff , zgrid, hwide , 
     $               ustar, bfsfc, hbl  , stable, caseA, 
     $               visc , difs , dift , kbl   , 
     $               gat1 , dat1 , dkm1 , blmc  , ghats, 
     $               sigma, wm   , ws   )
 
c-----------------------------------------------------------------------
c enhance diffusivity at interface kbl - 1
c-----------------------------------------------------------------------
c
      call enhance  (kmtj   , kmtjp , mdiff , dkm1 , visc  ,
     $               difs , dift , hbl  , kbl   , zgrid, caseA ,
     $               blmc ,ghats )
c 
c-----------------------------------------------------------------------
c combine interior and boundary layer coefficients and nonlocal term
c-----------------------------------------------------------------------
c

      do 300 ki= 1,kmtj
         if(ki.lt.kbl) then
            visc(ki)=blmc(1,ki)
            difs(ki)=blmc(2,ki)
            dift(ki)=blmc(3,ki)
         else
            ghats(ki)=0.
         endif
  300 continue
      

      return
      end

c *********************************************************************

      subroutine bldepth (
     $               km   , kmp1 , zgrid, hwide, dVsq  ,
     $               dbloc, Ritop, ustar , Bo   , Bosol, Coriol, jwtype,
     $               hbl  , bfsfc, stable, caseA, kbl  ,
     $               Rib  , sigma, wm    , ws   )

 
c
c     the oceanic planetray boundary layer depth, hbl, is determined as
c     the shallowest depth where the bulk richardson number is
c     equal to the critical value, Ricr.
c
c     bulk richardson numbers are evaluated by computing velocity and
c     buoyancy differences between values at zgrid(kl) < 0 and surface
c     reference values.
c     in this configuration, the reference values are equal to the
c     values in the surface layer.  
c     when using a very fine vertical grid, these values should be 
c     computed as the vertical average of velocity and buoyancy from 
c     the surface down to epsilon*zgrid(kl).
c
c     when the bulk richardson number at k exceeds Ricr, hbl is
c     linearly interpolated between grid levels zgrid(k) and zgrid(k-1).
c
c     The water column and the surface forcing are diagnosed for 
c     stable/ustable forcing conditions, and where hbl is relative 
c     to grid points (caseA), so that conditional branches can be 
c     avoided in later subroutines.
c
# include "comm_kpp.h"

c  model  
      integer km,kmp1      ! number of vertical levels
      real zgrid(0:kmp1)   ! vertical grid (<= 0)            (m)
      real hwide(0:kmp1)   ! layer thicknesses               (m)
c
c  input
      real dVsq(km)        ! (velocity shear re sfc)^2         (m/s)^2
      real dbloc(km)       ! local delta buoyancy              (m/s^2)
      real Ritop(km)       ! numerator of bulk Richardson Number
c                            =(z-zref)*dbsfc, where dbsfc=delta 
c                            buoyancy with respect to surface  (m/s)^2
      real ustar           ! surface friction velocity         (m/s)
      real Bo              ! surface turbulent buoyancy forcing(m^2/s^3)
      real Bosol           ! radiative buoyancy forcing        (m^2/s^3)
      real Coriol          ! Coriolis parameter                (1/s)
      integer jwtype       ! Jerlov water type                 (1 to 5)
c
c  output
      real hbl             ! boundary layer depth              (m)
      real bfsfc           ! Bo+radiation absorbed to d=hbf*hbl(m^2/s^3)
      real stable          ! =1 in stable forcing; =0 unstable
      real caseA           ! =1 in case A, =0 in case B 
      integer kbl          ! index of first grid level below hbl 
c 
c  local
      real Rib(2)          ! Bulk Richardson number
      real sigma           ! normalized depth (d/hbl)
      real wm,ws           ! turbulent velocity scales         (m/s)
c
c
c find bulk Richardson number at every grid level until > Ricr
c
c note: the reference depth is -epsilon/2.*zgrid(i,k), but the reference
c       u,v,t,s values are simply the surface layer values,
c       and not the averaged values from 0 to 2*ref.depth,
c       which is necessary for very fine grids(top layer < 2m thickness)
c note: max values when Ricr never satisfied are
c       kbl(i)=kmtj(i) and hbl(i)=-zgrid(i,kmtj(i))
c
c
c     indices for array Rib(i,k), the bulk Richardson number.
c
      ka = 1
      ku = 2
c
c     initialize hbl and kbl to bottomed out values
c
      Rib(ka) = 0.0
      kbl    = MAX(km,1)
      hbl    = -zgrid(kbl)
      
      do 200 kl = 2,km
c
c        compute bfsfc = sw fraction at hbf * zgrid
c
# ifdef sflxsw
         call swfrac(imt,hbf,zgrid(kl),jwtype,bfsfc)
# endif
         
c
c           use caseA as temporary array
c
         caseA  = -zgrid(kl)
c
c           compute bfsfc= Bo + radiative contribution down to hbf * hbl
c
         bfsfc  = Bo
# ifdef sflxsw
     $        + Bosol * (1. - bfsfc)
# endif
         stable = 0.5 + SIGN( 0.5, bfsfc )
         sigma  = stable * 1. + (1.-stable) * epsilon
         
         
c-----------------------------------------------------------------------
c        compute velocity scales at sigma, for hbl= caseA = -zgrid(i,kl)
c-----------------------------------------------------------------------
         
         call wscale(sigma, caseA, ustar, bfsfc,   wm, ws)
         
c-----------------------------------------------------------------------
c           compute the turbulent shear contribution to Rib
c-----------------------------------------------------------------------
c
         bvsq =0.5*
     $        ( dbloc(kl-1) / (zgrid(kl-1)-zgrid(kl  ))+
     $          dbloc(kl  ) / (zgrid(kl  )-zgrid(kl+1)) )
         Vtsq = - zgrid(kl) * ws * sqrt(abs(bvsq)) * Vtc
c 
c           compute bulk Richardson number at new level, dunder
c           note: Ritop needs to be zero on land and ocean bottom
c           points so that the following if statement gets triggered
c           correctly. otherwise, hbl might get set to (big) negative
c           values, that might exceed the limit for the "exp" function
c           in "swfrac"
c
         Rib(ku) = Ritop(kl) / (dVsq(kl)+Vtsq+epsl)
         if((kbl.eq.km).and.(Rib(ku).gt.Ricr)) then
            
            linearly interpolate to find hbl where Rib = Ricr
            
            hbl = -zgrid(kl-1) + (zgrid(kl-1)-zgrid(kl)) *
     $           (Ricr - Rib(ka)) / (Rib(ku)-Rib(ka))
            kbl = kl
            
         endif
         
         ksave = ka
         ka    = ku
         ku    = ksave
         
  200 continue
      
c -----------------------------------------------------------------------
c find stability and buoyancy forcing for boundary layer
c-----------------------------------------------------------------------
c
# ifdef sflxsw
      call swfrac(imt,-1.0,hbl,jwtype,bfsfc)
# endif
c
# ifdef sflxsw
      bfsfc  = Bo 
     $     + Bosol * (1. - bfsfc)
# endif
      stable = 0.5 + SIGN( 0.5, bfsfc )
      bfsfc  = bfsfc + stable * epsl !ensures bfsfc never=0
c
# ifdef kmixcheckekmo
c
c-----------------------------------------------------------------------
c check hbl limits for hekman or hmonob
c-----------------------------------------------------------------------
c
      
      if(bfsfc.gt.0.0) then
         hekman = cekman * ustar / (abs(Coriol)+epsl)
         hmonob = cmonob * ustar*ustar*ustar 
     $        /vonk / (bfsfc+epsl) 
         hlimit = stable     * AMIN1(hekman,hmonob) +
     $        (stable-1.) * zgrid(km)
         hbl = AMIN1(hbl,hlimit)
         hbl = AMAX1(hbl,-zgrid(1))
      endif
c
c-----------------------------------------------------------------------
c      find new kbl
c-----------------------------------------------------------------------
c
      kbl = km
      do 405 kl=2,km
         if((kbl.eq.km).and.(-zgrid(kl).gt.hbl)) then
            kbl = kl
         endif
  405 continue
c
c-----------------------------------------------------------------------
c      find stability and buoyancy forcing for final hbl values
c-----------------------------------------------------------------------
c
#  ifdef sflxsw
      call swfrac(imt,-1.0,hbl,jwtype,bfsfc)
#  endif
      
#  ifdef sflxsw
      bfsfc  = Bo 
     $     + Bosol * (1. - bfsfc)
#  endif
      stable = 0.5 + SIGN( 0.5, bfsfc )
      bfsfc  = bfsfc + stable * epsl 
      
# endif
      
c
c-----------------------------------------------------------------------
c determine caseA and caseB
c-----------------------------------------------------------------------
c
      caseA  = 0.5 + 
     $     SIGN( 0.5,-zgrid(kbl) - 0.5*hwide(kbl) - hbl)
      
      
      return
      end
      
c *********************************************************************
      
      subroutine wscale(sigma, hbl, ustar, bfsfc, wm , ws   )
      
c     compute turbulent velocity scales.
c     use a 2D-lookup table for wm and ws as functions of ustar and
c     zetahat (=vonk*sigma*hbl*bfsfc).
c
c     note: the lookup table is only used for unstable conditions
c     (zehat.le.0), in the stable domain wm (=ws) gets computed
c     directly.
c
# include "comm_kpp.h"

c  input
      real sigma      ! normalized depth (d/hbl)
      real hbl        ! boundary layer depth (m)
      real ustar      ! surface friction velocity         (m/s)
      real bfsfc      ! total surface buoyancy flux       (m^2/s^3)

c  output
      real wm,ws      ! turbulent velocity scales at sigma

c local
      real zehat           ! = zeta *  ustar**3

c
c-----------------------------------------------------------------------
c use lookup table for zehat < zmax only; otherwise use
c stable formulae
c-----------------------------------------------------------------------
c
      zehat = vonk * sigma * hbl * bfsfc
      
      if (zehat.le.zmax) then
         
         zdiff  = zehat-zmin
         iz = int( zdiff/deltaz )
         iz = min( iz , nni )
         iz = max( iz , 0  )
         izp1=iz+1
         
         udiff  = ustar-umin
         ju = int( udiff/deltau)
         ju = min( ju , nnj )
         ju = max( ju , 0  )
         jup1=ju+1
         
         zfrac = zdiff/deltaz - float(iz)
         ufrac = udiff/deltau - float(ju)
         
         fzfrac= 1.-zfrac
         wam   = (fzfrac)  * wmt(iz,jup1) + zfrac*wmt(izp1,jup1)
         wbm   = (fzfrac)  * wmt(iz,ju  ) + zfrac*wmt(izp1,ju  )
         wm = (1.-ufrac)* wbm          + ufrac*wam
         
         was   = (fzfrac)  * wst(iz,jup1) + zfrac*wst(izp1,jup1)
         wbs   = (fzfrac)  * wst(iz,ju  ) + zfrac*wst(izp1,ju  )
         ws = (1.-ufrac)* wbs          + ufrac*was
         
      else
         
         u3    = ustar*ustar*ustar
         wm = vonk * ustar * u3 / ( u3 + conc1*zehat )
         ws = wm
         
      endif
      
      return
      end
      
c **********************************************************************
      subroutine ri_iwmix (km, kmp1, Shsq, dbloc, zgrid, visc, difs, dift)
      
      
c     compute interior viscosity diffusivity coefficients due 
c     to shear instability (dependent on a local richardson number),
c     to background internal wave activity, and 
c     to static instability (local richardson number < 0).
c
# include "comm_kpp.h"
      
c  input
      real Shsq(km)          ! (local velocity shear)^2          (m/s)^2
      real dbloc(km)         ! local delta buoyancy              (m/s^2)
      real zgrid(0:kmp1)     ! vertical grid (<= 0)              (m)
      integer km,kmp1        ! number of vertical levels
 
c output
      real visc(0:kmp1)      ! vertical viscosivity coefficient  (m^2/s)
      real difs(0:kmp1)      ! vertical scalar diffusivity       (m^2/s)
      real dift(0:kmp1)      ! vertical temperature diffusivity  (m^2/s)
 
c local variables
      real Rigg              ! local richardson number
      real fri,fcon          ! function of Rig

c-----------------------------------------------------------------------
c     compute interior gradient Ri at all interfaces ki=1,km, (not surface)
c       use visc(0      ) as temporary storage for local Ri 
c       use visc(ki=1,km) as temporary storage to be smoothed
c       use dift(ki=1,km) as temporary storage of N^2 = BVSQ    
c       set values at bottom and below to nearest value above bottom
c-----------------------------------------------------------------------
      
      do 100 ki = 1, km
         
         visc(0)  = dbloc(ki) * (zgrid(ki)-zgrid(ki+1)) /
     $        ( Shsq(ki) + epsl)
         kterm =  MAX (km,1)
         wait  = AMIN1(float( ki / kterm)  , c1 )
         visc(ki) = (c1-wait) * visc(0)  + wait * visc(ki-1)
         dift(ki) = dbloc(ki) / (zgrid(ki)-zgrid(ki+1))

         
  100 continue

c-----------------------------------------------------------------------
c     vertically smooth Ri num_v_smooth_Ri times
c-----------------------------------------------------------------------
      
      do mr = 1,num_v_smooth_Ri
         call z121(km+1,visc)
      enddo
      
c-----------------------------------------------------------------------
c                           after smoothing loop
c-----------------------------------------------------------------------
      do 210 ki = 1, km
c-----------------------------------------------------------------------
c  evaluate f of Vaisala squared    for convection        store in fcon
c  evaluate f of   smooth Ri (fri) for shear instability store in fri
c-----------------------------------------------------------------------
         
         Rigg  = AMAX1( dift(ki) , BVSQcon )
         ratio = AMIN1( (BVSQcon-Rigg)/BVSQcon , c1 )
         fcon  = (c1 - ratio*ratio)
         fcon  = fcon * fcon * fcon
         
         Rigg  = AMAX1( visc(ki) , c0 )
         ratio = AMIN1( Rigg/Riinfty , c1 )
         fri   = (c1 - ratio*ratio)
         fri   = fri  * fri   * fri

c ----------------------------------------------------------------------
c            evaluate diffusivities and viscosity
c    mixing due to internal waves, and shear and static instability
c-----------------------------------------------------------------------
         
         visc(ki) = (difmiw + fcon * difmcon + fri * difm0)
         difs(ki) = (difsiw + fcon * difscon + fri * difs0)
         dift(ki) = difs(ki)

         
  210 continue
      
      
c-----------------------------------------------------------------------
c         set surface values to 0.0
c-----------------------------------------------------------------------
      
      visc(0)    = c0
      dift(0)    = c0
      difs(0)    = c0
      
      return
      end
      
      
      subroutine z121 (kmp1,v)
      
c-----------------------------------------------------------------------
c    Apply 121 smoothing in k to 2-d array V(i,k=1,km)
c     top (0) value is used as a dummy
c     bottom (kmp1) value is set to input value from above.
c-----------------------------------------------------------------------
 
c  input
       real V(0:kmp1)           ! 2-D array to be smoothed in kmp1 direction
       
       km  = kmp1 - 1
       p25 = 0.25
       p5  = 0.5
       
       V(0)    =  p25 * V(1)
       V(kmp1) =        V(km)
       
       do k=1,km
          tmp      =  V(k)
          V(k)   =  V(0)  + p5 * V(k) + p25 * V(k+1)
          V(0)   =  p25 * tmp
       enddo
       
       return
       end

c *********************************************************************

      subroutine blmix 
     $             (km   , kmp1 , mdiff , zgrid, hwide , 
     $              ustar, bfsfc, hbl , stable, caseA, 
     $              visc , difs , dift, kbl   , 
     $              gat1 , dat1 , dkm1, blmc  , ghats, 
     $              sigma, wm   , ws  )

 
c
c     mixing coefficients within boundary layer depend on surface
c     forcing and the magnitude and gradient of interior mixing below
c     the boundary layer ("matching").
c
c     caution: if mixing bottoms out at hbl = -zgrid(km) then
c     fictitious layer at kmp1 is needed with small but finite width 
c     hwide(kmp1) (eg. epsl = 1.e-20).
c
# include "comm_kpp.h"

c model
      integer km,kmp1        ! number of vertical levels
      integer mdiff          ! number of viscosities + diffusivities
      real zgrid(0:kmp1) ! vertical grid (<=0)               (m)
      real hwide(0:kmp1) ! layer thicknesses                 (m)
c
c input
      real ustar        ! surface friction velocity         (m/s)
      real bfsfc        ! surface buoyancy forcing        (m^2/s^3)
      real hbl          ! boundary layer depth              (m)
      real stable       ! = 1 in stable forcing
      real caseA        ! = 1 in case A
      real visc(0:kmp1)  ! vertical viscosity coefficient    (m^2/s)
      real difs(0:kmp1)  ! vertical scalar diffusivity       (m^2/s)
      real dift(0:kmp1)  ! vertical temperature diffusivity  (m^2/s)
      integer kbl       ! index of first grid level below hbl
c
c output
      real gat1(mdiff)
      real dat1(mdiff)
      real dkm1(mdiff)   ! boundary layer difs at kbl-1 level
      real blmc(mdiff,km)    ! boundary layer mixing coefficients(m^2/s)
      real ghats(km)     ! nonlocal scalar transport
c
c  local
      real sigma        ! normalized depth (d / hbl)
      real ws, wm       ! turbulent velocity scales         (m/s)
 
c
c-----------------------------------------------------------------------
c compute velocity scales at hbl
c-----------------------------------------------------------------------
c
      sigma = stable * 1.0 + (1.-stable) * epsilon
c
      call wscale(sigma, hbl, ustar, bfsfc,   wm, ws)
c 
      kn    = ifix(caseA+epsl) *(kbl -1) +
     $     (1-ifix(caseA+epsl)) * kbl
c
c-----------------------------------------------------------------------
c find the interior viscosities and derivatives at hbl(i) 
c-----------------------------------------------------------------------
c
      delhat = 0.5*hwide(kn) - zgrid(kn) - hbl
      R      = 1.0 - delhat / hwide(kn)
      dvdzup = (visc(kn-1) - visc(kn)) / hwide(kn) 
      dvdzdn = (visc(kn)   - visc(kn+1)) / hwide(kn+1)
      viscp  = 0.5 * ( (1.-R) * (dvdzup + abs(dvdzup))+
     $     R  * (dvdzdn + abs(dvdzdn)) )
c
      dvdzup = (difs(kn-1) - difs(kn)) / hwide(kn) 
      dvdzdn = (difs(kn)   - difs(kn+1)) / hwide(kn+1)
      difsp  = 0.5 * ( (1.-R) * (dvdzup + abs(dvdzup))+
     $     R  * (dvdzdn + abs(dvdzdn)) )
c
      dvdzup = (dift(kn-1) - dift(kn)) / hwide(kn) 
      dvdzdn = (dift(kn)   - dift(kn+1)) / hwide(kn+1)
      diftp  = 0.5 * ( (1.-R) * (dvdzup + abs(dvdzup))+
     $     R  * (dvdzdn + abs(dvdzdn)) )
c
      visch  = visc(kn) + viscp * delhat
      difsh  = difs(kn) + difsp * delhat
      difth  = dift(kn) + diftp * delhat
c
      f1 = stable * conc1 * bfsfc / (ustar**4+epsl) 
      
      gat1(1) = visch / hbl / (wm+epsl)
      dat1(1) = -viscp / (wm+epsl) + f1 * visch
      dat1(1) = min(dat1(1),0.) 
      
      gat1(2) = difsh  / hbl / (ws+epsl)
      dat1(2) = -difsp / (ws+epsl) + f1 * difsh 
      dat1(2) = min(dat1(2),0.) 
      
      gat1(3) = difth /  hbl / (ws+epsl)
      dat1(3) = -diftp / (ws+epsl) + f1 * difth 
      dat1(3) = min(dat1(3),0.) 
      
c
      do 300 ki = 1,km       
c
c-----------------------------------------------------------------------
c     compute turbulent velocity scales on the interfaces
c-----------------------------------------------------------------------
c
         sig     = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl
         sigma= stable*sig + (1.-stable)*AMIN1(sig,epsilon)
         call wscale(sigma, hbl, ustar, bfsfc,   wm,  ws)
c
c-----------------------------------------------------------------------
c     compute the dimensionless shape functions at the interfaces
c-----------------------------------------------------------------------
c
         sig = (-zgrid(ki) + 0.5 * hwide(ki)) / hbl
         a1 = sig - 2.
         a2 = 3.-2.*sig
         a3 = sig - 1.
c
         Gm = a1 + a2 * gat1(1) + a3 * dat1(1) 
         Gs = a1 + a2 * gat1(2) + a3 * dat1(2)
         Gt = a1 + a2 * gat1(3) + a3 * dat1(3)
c
c-----------------------------------------------------------------------
c     compute boundary layer diffusivities at the interfaces
c-----------------------------------------------------------------------
c
         blmc(1,ki) = hbl * wm * sig * (1. + sig * Gm)
         blmc(2,ki) = hbl * ws * sig * (1. + sig * Gs)
         blmc(3,ki) = hbl * ws * sig * (1. + sig * Gt)
c
c-----------------------------------------------------------------------
c     nonlocal transport term = ghats * <ws>o
c-----------------------------------------------------------------------
c
         ghats(ki) = (1.-stable) * cg / (ws * hbl +epsl)
         
         
  300 continue
c
c-----------------------------------------------------------------------
c find diffusivities at kbl-1 grid level 
c-----------------------------------------------------------------------
c
      sig      =  -zgrid(kbl-1)  / hbl
      sigma =stable * sig + (1.-stable) * MIN(sig,epsilon)
c
      call wscale(sigma, hbl, ustar, bfsfc,   wm, ws)
c
      sig = -zgrid(kbl-1) / hbl
      a1= sig - 2.
      a2 = 3.-2.*sig
      a3 = sig - 1.
      Gm = a1 + a2 * gat1(1) + a3 * dat1(1)
      Gs = a1 + a2 * gat1(2) + a3 * dat1(2)
      Gt = a1 + a2 * gat1(3) + a3 * dat1(3)
      dkm1(1) = hbl * wm * sig * (1. + sig * Gm)
      dkm1(2) = hbl * ws * sig * (1. + sig * Gs)
      dkm1(3) = hbl * ws * sig * (1. + sig * Gt)
      
      
      return
      end
      
c ******************************************************************

      subroutine enhance (km   , kmp1  , mdiff , dkm1  , visc ,
     &                    difs , dift  , hbl   , kbl   , zgrid , caseA,
     &                    blmc , ghats )

c
c enhance the diffusivity at the kbl-.5 interface
c
c input
      integer km,kmp1        ! number of vertical levels
      integer mdiff          ! number of viscosities + diffusivities
      integer kbl       ! grid above hbl
      real hbl          ! boundary layer depth             (m)
      real dkm1(mdiff)   ! bl diffusivity at kbl-1 grid level
      real zgrid(0:kmp1) ! vertical grid (<= 0)             (m)
      real visc(0:kmp1)  ! enhanced viscosity               (m^2/s) 
      real difs(0:kmp1)  ! enhanced thermal diffusivity     (m^2/s)
      real dift(0:kmp1)  ! enhanced scalar  diffusivity     (m^2/s)
      real caseA        ! = 1 in caseA, = 0 in case B
c
c input/output
      real ghats(km)     ! nonlocal transport               (s/m**2)
c                              modified ghats at kbl(i)-1 interface
c output
      real blmc(mdiff,km)! enhanced bound. layer mixing coeff.
c
c local
      real delta             ! fraction hbl lies beteen zgrid neighbors
         
      ki  = kbl - 1
            
      delta = (hbl+zgrid(ki)) / (zgrid(ki)-zgrid(ki+1))
      
      dkmp5 = caseA * visc(ki) + (1.-caseA) * blmc(1,ki)
      dstar = (1.-delta)**2 * dkm1(1) + delta**2 * dkmp5      
      blmc(1,ki) = (1.-delta) * visc(ki) + delta * dstar
      
      dkmp5 = caseA * difs(ki) + (1.-caseA) * blmc(2,ki)
      dstar = (1.-delta)**2 * dkm1(2) + delta**2 * dkmp5    
      blmc(2,ki) = (1.-delta) * difs(ki) + delta * dstar
      
      dkmp5 = caseA * dift(ki) + (1.-caseA) * blmc(3,ki)
      dstar = (1.-delta)**2 * dkm1(3) + delta**2 * dkmp5     
      blmc(3,ki) = (1.-delta) * dift(ki) + delta * dstar
      
      ghats(ki) = (1.-caseA) * ghats(ki)
      
      return
      end
      
c***********************************************************************
      
      subroutine kmixinit
c
c     initialize some constants for kmix subroutines, and initialize
c     for kmix subroutine "wscale" the 2D-lookup table for wm and ws 
c     as functions of ustar and zetahat (=vonk*sigma*hbl*bfsfc).
c
 
# include "comm_kpp.h"
# include "kpp.h"

c local
      real zehat                        ! = zeta *  ustar**3
      real zeta                         ! = stability parameter d/L

      epsilon = epsilon_def
      vonk  = vonk_def
      conc1 = conc1_def
      conam =  conam_def
      concm  = concm_def
      conc2  = conc2_def
      zetam  = zetam _def
      conas  = conas_def
      concs  = concs_def
      conc3  = conc3_def
      zetas  = zetas_def

      Ricr =  Ricr_def
      cekman =  cekman_def
      cmonob =  cmonob_def
      concv =  concv _def
      hbf =  hbf_def

      zmin =  zmin_def
      zmax =  zmax_def
      umin =  umin_def
      umax =  umax_def

      Riinfty =  Riinfty_def
      BVSQcon =  BVSQcon_def
      num_v_smooth_Ri =  num_v_smooth_Ri_def
 
      fkph =  fkph_def
      fkpm =  fkpm_def
      vvclim =  vvclim_def
      vdclim =  vdclim_def
      vvcric =  vvcric_def
      vdcric =  vdcric_def

      cstar =  cstar_def
      jwatertype =  jwatertype_def
c
c-----------------------------------------------------------------------
c define some non-dimensional constants and 
c the vertical mixing coefficients in m-k-s units
c-----------------------------------------------------------------------
c
      epsl    = epsln
      Vtc     = concv * sqrt(0.2/concs/epsilon) / vonk**2 / Ricr
      cg      = cstar * vonk * (concs * vonk * epsilon)**(1./3.)
c
      difm0   = vvcric / c10000 
      difs0   = vdcric / c10000
      difmiw  = fkpm   / c10000
      difsiw  = fkph   / c10000
      difmcon = vvclim / c10000
      difscon = vdclim / c10000
c
c-----------------------------------------------------------------------
c construct the wm and ws lookup tables
c-----------------------------------------------------------------------
c
      deltaz = (zmax-zmin)/(nni+1) 
      deltau = (umax-umin)/(nnj+1)
      
      do 100 i=0,nni+1
         zehat = deltaz*(i) + zmin
         do 90 j=0,nnj+1
            usta = deltau*(j) + umin
            zeta = zehat/(usta**3+epsl)

            if(zehat.ge.0.) then
               wmt(i,j) = vonk*usta/(1.+conc1*zeta)
               wst(i,j) = wmt(i,j)
            else
               if(zeta.gt.zetam) then
                  wmt(i,j) = vonk* usta * (1.-conc2*zeta)**(1./4.)
               else
                  wmt(i,j) = vonk* (conam*usta**3-concm*zehat)**(1./3.)
               endif
               if(zeta.gt.zetas) then
                  wst(i,j) = vonk* usta * (1.-conc3*zeta)**(1./2.)
               else
                  wst(i,j) = vonk* (conas*usta**3-concs*zehat)**(1./3.)
               endif
            endif   
 90      continue   
 100  continue
c
 
      return
      end


      subroutine statekpp (nz, temp, salt, zgrid,
     $                  rho,  dbloc,  dbsfc, talpha, sbeta  )

c-------------------------------------------
c in situ density (kg/m^3) as a function of
c potential temperature, salinity & pressure.

      dimension temp(nz), salt(nz), zgrid(nz), 
     $          dbloc(nz), dbsfc(nz), talpha(nz), sbeta(nz)

      parameter (GRAVTY = 9.8) 

      data a0  /  9.9984495075397717E+02 /
      data a1  /  4.7242360502556559E-02 /
      data a2  / -2.9355135788249624E-03 /
      data a3  / -7.4801638387318654E-03 /
      data a4  /  3.3806262972062488E-05 /
      data a5  /  3.6995743312330793E-07 /
      data a6  /  3.5589087697669534E-05 /
      data a7  /  8.0747615263133631E-01 /
      data a8  / -8.0022984647632515E-06 /
      data a9  /  5.0377429128311419E-03 /
      data a10 / -2.9566247219662127E-05 /
      data a11 / -5.1678728359375112E-08 /
      data a12 /  3.7388128423099962E-10 /


c   polynomial approximation function for density
      pdens12(tq,sq,pq) = a0 +   
     *     tq*(a1 + sq*a2  + tq*(a3  + sq*a4 + tq*a6)) +
     *     sq*(a7 + pq*a8) + 
     *     pq*(a9 + tq*(a10 + tq*a5) + pq*(a11 + tq*a12)) 


c   d(pdens)/d(temp)
      alpha12(tq,sq,pq) = a1 + sq*a2  + pq*(a10 + pq*a12) + 
     *     tq*(2.*(a3 + sq*a4 + pq*a5) + tq*3.*tq*a6)

c   d(pdens)/d(salt)
      beta12(tq,sq,pq) =  a7 + tq*(a2 + tq*a4) + pq*a8

      k = 1
      presk      = 0.
      rho = pdens12(temp(k),  salt(k),  presk)
      talpha(k)  = alpha12(temp(k),  salt(k),  presk)
      sbeta(k)   = beta12 (temp(k),  salt(k),  presk)
      dbsfc(k)   = 0.

      do k = 2, nz
         presk      = -zgrid(k)
         rhok       = pdens12(temp(k),  salt(k),  presk)
         rhokm1     = pdens12(temp(k-1),salt(k-1),presk)
         rho1       = pdens12(temp(1),  salt(1),  presk)
         talpha(k)  = alpha12(temp(k),  salt(k),  presk)
         sbeta(k)   = beta12 (temp(k),  salt(k),  presk)
         dbloc(k-1) = GRAVTY * (1. - rhokm1/rhok)
         dbsfc(k)   = GRAVTY * (1. - rho1/rhok)
      enddo

      k = nz
      talpha(k+1)  = talpha(k)
      sbeta(k+1)   = sbeta (k)
      dbloc(k)     = 0.

      return
      end
