C **********************************************************************
C *                                                                    *
C *                                SOFTWARE LICENSING                  *
C *                                                                    *
C *                This program is free software; you can redistribute *
C *                it and/or modify it under the terms of the GNU      *
C *                General Public License as published by the Free     *
C *                Software Foundation, either Version 2 of the        *
C *                license, or (at your option) any later version.     *
C *                                                                    *
C *                This program is distributed in the hope that it     *
C *                will be useful, but without any warranty; without   *
C *                even the implied warranty of merchantability or     *
C *                fitness for a particular purpose. See the GNU       *
C *                General Public License for more details.            *
C *                                                                    *
C *                A copy of the GNU General Public License is         *
C *                available at http://www.gnu.org/copyleft/gpl.html   *
C *                or by writing to the Free Software Foundation, Inc.,*
C *                59 Temple Place - Suite 330, Boston, MA 02111, USA. *
C *                                                                    *
C **********************************************************************

      SUBROUTINE PROFE
C     VERSION(03/15/03)
C          Mixing formulation of PRITCHARD(19??) and OFFICER (1976)

c     ======>  constants to be read in by user:  prit(im,jm)

      INCLUDE 'comdeck'
      SAVE
Co
      REAL Khomo(IM,JM,KB)
      DIMENSION ZDIST(IM,JM,KB),SHEAR(IM,JM,KB)
      DIMENSION RI(IM,JM,KB),buoy(IM,JM,KB),STAB(IM,JM,KB)
      DIMENSION DH(IM,JM),PRIT(IM,JM)
c
      gee=9.806
      small=1.e-10
C
      DO 10 J=2,JMM1
      DO 10 I=2,IMM1
      PRIT(I,J)=VPRNU    ! HARDWIRED
 10   DH(I,J)=H(I,J)+ETF(I,J)
C
      DO 105 J=2,JMM1
      DO 105 K=1,KB
      DO 105 I=2,IMM1
 105  ZDIST(I,J,K)=-DH(I,J)*DH(I,J)*Z(K)*(1.+Z(K))**2
c105  ZDIST(I,J,K)=((DH(I,J)*Z(K)+ETF(I,J))+(1.+Z(K)))
C
      DO 115 J=2,JMM1
      DO 115 K=2,KBM1
      DO 115 I=2,IMM1
      Buoy(I,J,K)=-GEE*(RHO(I,J,K-1)-RHO(I,J,K))/(DZZ(K-1)*DH(I,J))
c     next line is a crude assumption, we should 
c         mix the unstable region via "convective adjustment"
      if(buoy(i,j,k).lt.0.) buoy(i,j,k)=0.0
 115  CONTINUE
C
      DO 120 J=2,JMM1
      DO 120 K=2,KBM1
      DO 120 I=2,IMM1
      SHEAR(I,J,K)=.25*
     .       ( (U(I,J,K)-U(I,J,K-1)+U(I+1,J,K)-U(I+1,J,K-1))**2
     .         +(V(I,J,K)-V(I,J,K-1)+V(I,J+1,K)-V(I,J+1,K-1))**2)
     .              /(DZZ(K-1)*DH(I,J))**2
 120  RI(I,J,K)=Buoy(I,J,K)/(SHEAR(I,J,K)+small)
C
      DO 110 K=2,KBM1
      DO 110 J=2,JMM1
      DO 110 I=2,IMM1
 110  KHOMO(I,J,K)=0.41*0.41*ZDIST(I,J,K)*sqrt(SHEAR(I,J,K))
c110  KHOMO(I,J,K)=prit(I,J)*ZDIST(I,J,K)**2*sqrt(SHEAR(I,J,K))
C
      DO 280 J=2,JMM1
      DO 280 K=2,KBM1
      DO 280 I=2,IMM1
      STAB(I,J,K)=(1.+RI(I,J,K))
      KM(I,J,K)=PRIT(I,J)*KHOMO(I,J,K)/STAB(I,J,K)
      KH(I,J,K)=PRIT(I,J)*KHOMO(I,J,K)/STAB(I,J,K)**2
 280  CONTINUE
C
      RETURN
      END