 
      Subroutine PUTVEC(CDFID,NAMVAR,RECN,VALUES,NK)
      Implicit None
      Integer CDFID, RECN, NK, ID, RCODE, NCVID
      Integer START(2), COUNT(2)
      Real VALUES(NK)
      Character*(*) NAMVAR
C***********************************************************************
C
C     PURPOSE: TO WRITE A VECTOR OF DIMENSION NK FROM  A PREVIOUSLY
C     DEFINED 'NETCDF' FILE .
C
C     ARGUMENTS :
C                -INPUT :
C                        - CDFID  : NETCDF FILE NAME 
C                        - NAMVAR :  'NETCDF' NAME OF VARIABLE TO INSERT
C                                     VALUES       
C                        -   RECN :  'NETCDF' RECORD NUMBER OF VARIABLE
C                        - VALUES :   VALUES TO INSERT      
C                        - NK     :  DIMENSION OF VARIABLE      
C
C                -OUTPUT: ---NONE---
C
C
C     NOTE: NETCDF LIBRARY IS NEEDED FOR EXECUTION OF THIS ROUTINE.
C           THIS ROUTINE WILL WORK WITHIN THE CONTEXT OF NETCDF FILES
C           DESIGNED FOR THE PBL EVALUATION PROJECT.
C
C     LANGUAGE  : FORTRAN
C     PROGRAMMER: P. KOCLAS , FEBRUARY 1991
C
C***********************************************************************
      ID = NCVID(CDFID,NAMVAR,RCODE)
      START(1) = 1
      START(2) = RECN
      COUNT(1) = NK
      COUNT(2) = 1
      Call NCVPT(CDFID,ID,START,COUNT,VALUES,RCODE)
      Return
      End


      Subroutine GETVEC(CDFID,NAMVAR,RECN,VALUES,NK)
      Implicit None
      Integer CDFID, RECN, NK, ID, RCODE, NCVID
      Integer START(2), COUNT(2)
      Real VALUES(NK)
      Character*(*) NAMVAR
C***********************************************************************
C
C     PURPOSE: TO READ A VECTOR OF DIMENSION NK FROM  A PREVIOUSLY
C     DEFINED 'NETCDF' FILE .
C
C     ARGUMENTS :
C                -INPUT :
C                        - CDFID  : NETCDF FILE NAME 
C                        - NAMVAR :  'NETCDF' NAME OF VARIABLE TO READ 
C                        -   RECN :  'NETCDF' RECORD NUMBER OF VARIABLE 
C                        - VALUES :   VALUES TO READ 
C                        - NK     :  DIMENSION OF VARIABLE      
C
C                        - VALUES :   VALUES TO READ
C
C
C     NOTE: NETCDF LIBRARY IS NEEDED FOR EXECUTION OF THIS ROUTINE.
C           THIS ROUTINE WILL WORK WITHIN THE CONTEXT OF NETCDF FILES
C           DESIGNED FOR THE PBL EVALUATION PROJECT.
C
C     LANGUAGE  : FORTRAN
C     PROGRAMMER: P. KOCLAS , FEBRUARY 1991
C
C***********************************************************************
      ID = NCVID(CDFID,NAMVAR,RCODE)
      START(1) = 1
      START(2) = RECN
      COUNT(1) = NK
      COUNT(2) = 1
      Call NCVGT(CDFID,ID,START,COUNT,VALUES,RCODE)
      Return
      End


      Subroutine PUTIVEC(CDFID,NAMVAR,VALUES,NVALS)
      Implicit None
      Integer CDFID, NVALS, ID, RCODE, NCVID, START, COUNT
      Integer VALUES(NVALS)
      Character*(*) NAMVAR
C***********************************************************************
C
C     PURPOSE: TO WRITE A VECTOR OF DIMENSION NVALS FROM A PREVIOUSLY
C     DEFINED 'NETCDF' FILE .
C
C
C     ARGUMENTS :
C                -INPUT :
C                        - CDFID  : NETCDF FILE NAME 
C                        - NAMVAR :  'NETCDF' NAME OF VARIABLE TO READ
C                        - VALUES :   VALUES TO READ 
C                        - NVALS  :  DIMENSION OF VARIABLE      
C
C                -OUTPUT: ---NONE---
C
C
C     NOTE: NETCDF LIBRARY IS NEEDED FOR EXECUTION OF THIS ROUTINE.
C           THIS ROUTINE WILL WORK WITHIN THE CONTEXT OF NETCDF FILES
C           DESIGNED FOR THE PBL EVALUATION PROJECT.
C
C     LANGUAGE  : FORTRAN
C     PROGRAMMER: P. KOCLAS , FEBRUARY 1991
C
C***********************************************************************
      ID = NCVID(CDFID,NAMVAR,RCODE)
      START = 1
      COUNT = NVALS
      Call NCVPT(CDFID,ID,START,COUNT,VALUES,RCODE)
      Return
      End


      Subroutine GETIVEC(CDFID,NAMVAR,VALUES,NVALS)
      Implicit None
      Integer CDFID, NVALS, ID, RCODE, NCVID, START, COUNT
      Integer VALUES(NVALS)
      Character*(*) NAMVAR
C***********************************************************************
C
C     PURPOSE: TO READ A VECTOR OF DIMENSION NVALS FROM A PREVIOUSLY
C     DEFINED 'NETCDF' FILE .
C
C     ARGUMENTS :
C                -INPUT :
C                        - CDFID  : NETCDF FILE NAME
C                        - NAMVAR :  'NETCDF' NAME OF VARIABLE TO READ
C                        - VALUES :   VALUES TO READ
C                        - NVALS  :  DIMENSION OF VARIABLE
C
C                        - VALUES :   VALUES TO READ
C
C
C     NOTE: NETCDF LIBRARY IS NEEDED FOR EXECUTION OF THIS ROUTINE.
C           THIS ROUTINE WILL WORK WITHIN THE CONTEXT OF NETCDF FILES
C           DESIGNED FOR THE PBL EVALUATION PROJECT.
C
C     LANGUAGE  : FORTRAN
C     PROGRAMMER: P. KOCLAS , FEBRUARY 1991
C
C
C***********************************************************************
      ID = NCVID(CDFID,NAMVAR,RCODE)
      START = 1
      COUNT = NVALS
      Call NCVGT(CDFID,ID,START,COUNT,VALUES,RCODE)
      Return
      End


      Subroutine PUTSCA(CDFID,NAMVAR,RECN,VALUE)
      Implicit None
      Integer CDFID, RECN, ID, RCODE, NCVID
      Real VALUE
      Character*(*) NAMVAR
C***********************************************************************
C
C     PURPOSE: TO INSERT A SINGLE VALUE INTO A PREVIOUSLY
C     DEFINED 'NETCDF' FILE .
C
C     ARGUMENTS :
C                -INPUT :
C                        -  CDFID : NETCDF FILE NAME
C                        - NAMVAR :  'NETCDF' NAME OF VARIABLE TO INSERT
C                                     VALUE
C                        -   RECN :  'NETCDF' RECORD NUMBER OF VARIABLE
C                        -  VALUE :   VALUE TO INSERT
C
C                -OUTPUT: ---NONE---
C
C
C     NOTE: NETCDF LIBRARY IS NEEDED FOR EXECUTION OF THIS ROUTINE.
C           THIS ROUTINE WILL WORK WITHIN THE CONTEXT OF NETCDF FILES
C           DESIGNED FOR THE PBL EVALUATION PROJECT.
C
C     LANGUAGE  : FORTRAN
C     PROGRAMMER: P. KOCLAS , FEBRUARY 1991
C
C***********************************************************************
      ID = NCVID(CDFID,NAMVAR,RCODE)
      Call NCVPT1(CDFID,ID,RECN,VALUE,RCODE)
      Return
      End


      Subroutine GETSCA(CDFID,NAMVAR,RECN,VALUE)
      Implicit None
      Integer CDFID, RECN, ID, RCODE, NCVID
      Real VALUE
      Character*(*) NAMVAR
C***********************************************************************
C
C     PURPOSE: TO READ A SINGLE VALUE INTO A PREVIOUSLY
C     DEFINED 'NETCDF' FILE .
C
C     ARGUMENTS :
C                -INPUT :
C                        - CDFID  : NETCDF FILE NAME
C                        - NAMVAR :  'NETCDF' NAME OF VARIABLE TO READ
C                        -   RECN :  'NETCDF' RECORD NUMBER OF VARIABLE
C                        - VALUE  :   VALUE TO READ
C
C                -OUTPUT:
C                        - NAMVAR :  'NETCDF' NAME OF VARIABLE TO READ
C
C
C     NOTE: NETCDF LIBRARY IS NEEDED FOR EXECUTION OF THIS ROUTINE.
C           THIS ROUTINE WILL WORK WITHIN THE CONTEXT OF NETCDF FILES
C           DESIGNED FOR THE PBL EVALUATION PROJECT.
C
C     LANGUAGE  : FORTRAN
C     PROGRAMMER: P. KOCLAS , FEBRUARY 1991
C
C***********************************************************************
      ID = NCVID(CDFID,NAMVAR,RCODE)
      Call NCVGT1(CDFID,ID,RECN,VALUE,RCODE)
      Return
      End


      Subroutine GETREC(RECN,STEP,IREC,TIME,NSTEPS,OK,TOL)
      Implicit None
      Integer RECN, STEP, NSTEPS, I
      Integer IREC(NSTEPS), TIME(NSTEPS)
      Real TOL
      Logical OK
C***********************************************************************
C
C     PURPOSE: TO FIND THE (RECN) NUMBER CORRESPONDING TO TIMESTEP
C     (STEP) OF A PBL MODEL INTEGRATION.
C
C     ARGUMENTS :
C                -INPUT :
C                        -   RECN :  'NETCDF' RECORD NUMBER OF VARIABLE
C                        -   IREC :  ARRAY CONTAINING RECORD NUMBERS
C                        -   TIME :  ARRAY CONTAINING TIMESTEPS
C                        -  NSTEPS:  DIMENSION GIVING NUMBER OF
C                                    TIMESTEPS
C                        -     OK : LOGICAL VARIABLE =TRUE IF RECORD
C                                   NUMBER WAS FOUND =FALSE OTHERWISE.
C
C                -OUTPUT:
C                        -   RECN :  RECORD NUMBER
C
C
C     LANGUAGE  : FORTRAN
C     PROGRAMMER: P. KOCLAS , FEBRUARY 1991
C
C***********************************************************************
      OK = .TRUE.
      Do 10 I = 1, NSTEPS
        Write (*,*) TIME(I), STEP, I, NSTEPS
        If (ABS(TIME(I)-STEP).LE.TOL) Then
          RECN = IREC(I)
          Return
        End If
   10 Continue
      Print *, ' record # not found in file'
      OK = .FALSE.
      Return
      End

      Subroutine GETNAM(CDFID,VARNAM,DIMTYP,NAME,DIMSIZ,OK)
      Implicit None
      Integer ID, RCODE, VARTYP, NVDIMS, NVATTS, K, K1, K2
      Integer NCVID, CDFID, DIMSIZ, MAXDIMS
      Parameter (MAXDIMS = 2)
      Integer VDIMS(MAXDIMS), SIZE(MAXDIMS), LENGTH
      Character*(*) VARNAM, DIMTYP, NAME
      Character*20 DUMNAM, NAM(MAXDIMS)
      Logical OK
C***********************************************************************
C
C     PURPOSE : OBTAIN THE NAME AND SIZE OF A VARIABLE'S HEIGHT OR TIME
C               COORDINATE VIA IT'S NETCDF DIMENSIONS.
C     ARGUMENTS :
C                  INPUT:
C                        -CDFID : NETCDF FILE ID
C                        -VARNAM: NETCDF VARIABLE NAME
C                        -DIMTYP: TYPE OF DIMENSION ( I.E. 'Z' FOR HEIGHT
C                                  'TIME' FOR TIME
C                 OUTPUT:
C                        -NAME  : NAME OF COORDINATE
C                        -DIMSIZ: SIZE OF COORDINATE
C                        -    OK: = .TRUE. IF COORDINATE IS FOUND,
C                                 = .FALSE. IF COORDINATE IS NOT FOUND.
C
C     PROGRAMMER: P. KOCLAS (NCAR MMM) MAY 1991.
C
C     NOTE: NETCDF LIBRARY IS NEEDED FOR EXECUTION OF THIS ROUTINE.
C      ---> THIS ROUTINE WILL WORK WITHIN THE CONTEXT OF NETCDF FILES
C           DESIGNED FOR THE PBL EVALUATION PROJECT. <---
C
C           IT IS ASSUMED THAT VARIABLES ARE SAVED AS VAR=VAR(Z,TIME)
C           A LOGICAL FLAG IS RETURNED IF THE COORDINATE NAME IS
C           NOT FOUND.
C
C***********************************************************************

C     CHECK IF WE ARE ASKING FOR THE RIGHT TYPE OF COORDINATE
      OK = .TRUE.
      If (DIMTYP.NE.'Z'.AND.DIMTYP.NE.'TIME') Then
        OK = .FALSE.
        Return
      End If

C***********************************************************************
C     GET THE DIMENSION NAMES
C***********************************************************************
      ID = NCVID(CDFID,VARNAM,RCODE)
      Call NCVINQ(CDFID,ID,DUMNAM,VARTYP,NVDIMS,VDIMS,NVATTS,RCODE)

      Do 10 K = 1, NVDIMS
        Call NCDINQ(CDFID,VDIMS(K),NAM(K),SIZE(K),RCODE)
   10 Continue

      If (NVDIMS.GT.MAXDIMS) Then
        Print *, ' # of dimensions >', MAXDIMS, ' !!'
        OK = .FALSE.
      End If

C***********************************************************************
C     SET THE PROPER NAME AND SIZE OF COORDINATE WANTED
C***********************************************************************
      LENGTH = INDEX(NAM(1),' ') - 1
      If (NAM(1)(1:LENGTH).EQ.'TIME') Then
        K1 = 1
        K2 = 2
      Else
        K1 = 2
        K2 = 1
      End If

      If (DIMTYP.EQ.'TIME') Then
        NAME = NAM(K1)
        DIMSIZ = SIZE(K1)
      Else
        NAME = NAM(K2)
        DIMSIZ = SIZE(K2)
      End If

      Return
      End


      Subroutine GETGLOB(CDFID,ATTVAL,LENGST,ATTNAM,NCGLOBAL)
      Implicit None

      Integer CDFID, LENGST, NCGLOBAL, ID, RCODE, IDUM
      Character*(*) ATTNAM, ATTVAL
C***********************************************************************
C
C     PURPOSE : OBTAIN THE VALUE OF A GLOBAL ATTRIBUTE GIVEN IT'S NAME.
C
C     ARGUMENTS :
C                  INPUT:
C                        -CDFID : NETCDF FILE ID
C                        -ATTNAM: NAME OF ATTRIBUTE.
C                        -NCGLOBAL: NETCDF ID FOR GLOBAL ATTRIBUTES
C                 OUTPUT:
C                        -ATTVAL: VALUE OF GLOBAL ATTRIBUTE.
C                        -LENGST: LENGTH OF CHARACTER STRING OF ATTVAL.
C
C     PROGRAMMER: P. KOCLAS (NCAR MMM) 1991.
C
C     NOTE: NETCDF LIBRARY IS NEEDED FOR EXECUTION OF THIS ROUTINE.
C           THIS ROUTINE WILL WORK WITHIN THE CONTEXT OF NETCDF FILES
C           DESIGNED FOR THE PBL EVALUATION PROJECT.
C
C
C***********************************************************************
      ID = NCGLOBAL
      Call NCAINQ(CDFID,ID,ATTNAM,IDUM,LENGST,RCODE)
      Call NCAGTC(CDFID,ID,ATTNAM,ATTVAL,LENGST,RCODE)
      Return
      End


      Subroutine MAKENAM(FILNAM,LENGTH,CASE,MODEL)
      Implicit None

      Character*(*) FILNAM, CASE, MODEL
      Character*20 DIR_NAME, EXTEN
      Integer LENGTH, L1, L2, L3, L4
C***********************************************************************
C
C     PURPOSE : TO GENERATE A FILE NAME FOR FILES USED IN THE 'PBL
C               EVALUATION PROJECT'.
C     ARGUMENTS :
C
C             INPUT:
C                   -MODEL : NAME OF THE MODEL DATA CONTAINED IN THE FILE
C                   -CASE  : NAME OF THE PBL SIMULATION
C
C             OUTPUT:
C                   -FILNAME: THE FILE NAME
C                   - LENGTH: THE NUMBER OF CHARACTERS OF THE FILE NAME
C
C     PROGRAMMER: P. KOCLAS (NCAR MMM) 1991.
C
C***********************************************************************
      DIR_NAME = '../data/'
      EXTEN = '.cdf'

      L1 = INDEX(MODEL,' ') - 1
      L2 = INDEX(CASE,' ') - 1
      L3 = INDEX(DIR_NAME,' ') - 1
      L4 = INDEX(EXTEN,' ') - 1

      FILNAM = DIR_NAME(1:L3) // MODEL(1:L1) // CASE(1:L2) // 
     *    EXTEN(1:L4)
      LENGTH = L1 + L2 + L3 + L4 + 1

      Return
      End


      Subroutine SURFACE(U1,V1,T1,Q1,B1,C1,TSURF,QSURF,BSURF,CSURF,
     *    USTAR,TSTAR,QSTAR,BSTAR,CSTAR,Z1,ZA,Z0M,Z0H,Z0,UA,VA,TA,QA,BA,
     *    CA,ISURF,ISURFB,ISURFC,WUSFC,WVSFC,WTSFC,WQSFC,WBSFC,WCSFC,
     *    GRAV,TREF,VK,DEL_Z,ZI)

      Implicit None

C***********************************************************************
C
C     CALCULATES SURFACE BOUNDARY CONDITIONS USING M-O SIMILARITY
C     AND THE PROFILE RELATIONS OF BUSINGER ET AL. (1971)
C
C     INPUTS:  1) MEAN QUANTITES AT 1ST GRIDPOINT.
C              2) ROUGHNESS LENGTH.
C              3) SURFACE FLUXES OR SURFACE VALUES OF TEMPERATURE,
C                 HUMIDITY, AND SCALARS (DEPENDING ON SURFACE BOUNDARY
C                 CONDITION - CONSTANT FLUX OR CONSTANT VALUE AT SURFACE.
C              4) INTEGERS WHICH SPECIFY WHAT BOUNDARY CONDITION TO
C                 USE (ISURF FOR TEMPERATURE AND HUMIDITY, ISURFB FOR
C                 SCALAR B, AND ISURFC FOR SCALAR C.
C              5) GRAVITY, REFERENCE TEMPERATURE, VON KARMAN'S
C                 CONSTANT
C
C     OUTPUTS: 1) SURFACE LAYER SCALES (USTAR, TSTAR, QSTAR, BSTAR,
C                 CSTAR)
C              2) SURFACE FLUXES OR SURFACE VALUES (DEPENDING ON
C                 BOUNDARY CONDITION SPECIFIED)
C              3) ANEMOMETER HEIGHT VALUES (NEEDED FOR KPROF SCHEME)
C
C**********************************************************************

C     MEAN QUANTITIES AT THE LOWEST GRIDPOINT
      Real U1, V1, T1, Q1, B1, C1

C     MEAN QUANTITIES AT THE SURFACE
      Real TSURF, QSURF, BSURF, CSURF

C     SURFACE LAYER SCALES
      Real USTAR, TSTAR, QSTAR, BSTAR, CSTAR, U_FREE, SMAL

C     HEIGHT OF LOWEST GRIDPOINT, ANEMOMETER HEIGHT, AND ROUGHNESS
C     HEIGHTS
      Real Z1, ZA, Z0, Z0M, Z0H, ZI, DEL_Z

C     MEAN QUANTITIES AT THE ANEMOMETER HEIGHT
      Real UA, VA, TA, QA, BA, CA, WSPDA

C     SURFACE FLUXES
      Real WUSFC, WVSFC, WTSFC, WQSFC, WBSFC, WCSFC

C     VALUES AT THE LOWEST GRIDPOINT MINUS VALUES AT THE
C     SURFACE
      Real TMEAN, QMEAN, BMEAN, CMEAN

      Real GRAV, TREF, VK, WIND, XM1, ZETA, PSI_M, PSI_H, RATIO
      Real XM1OLD, THTA, TEP, BUSPSI_M1, BUSPSI_H2, PROFILE
      Real FPSI_M, FPSI_H
      Real USTAR1, USTAR2

C     PBL SIMILARITY VARIABLES
      Real CUN, CTHETAN, CU, CTHETA

C     NUMERICAL RECIPES SECANT ROOT FINDING ROUTINE AND RESIDUAL FUNCTION
      Real ZBRENT, USTAR_RES, ZBRAK
      External USTAR_RES, USTARPBL_RES

      Integer ISURF, ISURFB, ISURFC, ITEMP, ITEMPB
      Integer ITEMPC, ITER

      WIND = SQRT((U1**2)+(V1**2))
      Z0M = Z0
      Z0H = Z0
      SMAL = 0.00001

C     SET LOWER LIMIT ON wind TO GIVE u_free FOR THE FREE
C     CONVECTION CASE
      u_free = 0.0
      if(abs(wtsfc) .gt. smal) then
         U_FREE = 0.07 * (GRAV/TREF*WTSFC*DEL_Z) ** (1.0/3.0)
      endif
      WIND = AMAX1(WIND,U_FREE)

C     USE THE NEUTRAL VALUE OF ustar AS THE FIRST GUESS
      USTAR = VK * WIND / (ALOG((Z1+Z0M)/Z0M))

C//////////////////////////////////////////////////////////////////////
C     IF FIRST GRID POINT IS NEAR ENOUGH TO THE SURFACE (50 M OR LESS),
C     USE M-O SIMILARITY THEORY
C//////////////////////////////////////////////////////////////////////

      If (Z1.LE.50.0) Then
      
C        CHECK TO SEE IF CONDITIONS ARE NEUTRAL/NEARLY NEUTRAL      
        ZETA = -Z1 * VK * GRAV * WTSFC / (TREF*USTAR**3)
	 
C        IF ZETA IS VERY SMALL, DON'T DO NON-LINEAR ITERATION 
C        NOTE: 1.0e-4 AND 0.01 ARE ARBITRARY LIMITSS
        If (ABS(ZETA).GT.1.0E-4.AND.USTAR.GT.1.0E-2) Then
      
C           PUT A SIMPLE LOWER LIMIT ON ustar FOR FREE CONVECTION CASE
          If (USTAR.LT.0.01) USTAR = 0.01
	    
C           FIRST BRACKET THE ROOT
          USTAR1 = 0.0
          USTAR2 = 0.0
          Call ZBRAK(USTAR_RES,USTAR*0.5,USTAR*1.5,100,USTAR1,USTAR2,1,
     *        USTAR,WIND,TREF,GRAV,WTSFC,Z1,Z0M,VK,ZI)
     
C           CHECK FOR BRACKETED ROOT
          If (USTAR1.LT.1.0E-4.AND.USTAR2.LT.1.0E-4) Then
            USTAR1 = USTAR * 0.5
            USTAR2 = USTAR * 1.5
          End If
     
C           FIND ROOT (ustar)
          USTAR = ZBRENT(USTAR_RES,USTAR1,USTAR2,1.0E-6,USTAR,WIND,TREF,
     *        GRAV,WTSFC,Z1,Z0M,VK,ZI)

C           CALCULATE NEW VALUE OF ZETA (z/L) FOR NEW ustar
          ZETA = -Z1 * VK * GRAV * WTSFC / (TREF*USTAR**3)
        End If
	 
C        CALCULATE M-O PARAMETERS 
CTEMP         
C         psi_m = buspsi_m1(zeta)
        PSI_M = FPSI_M(ZETA)
CTEMP

C         psi_h = buspsi_h2(zeta)
        PSI_H = FPSI_H(ZETA)
        RATIO = PROFILE(VK,Z1,Z0H,PSI_H)

C        CALCULATE T* AND Q* OR (T1-T0) AND (Q1-Q0), DEPENDING ON ISURF
        If (ISURF.EQ.0) Then
          TSTAR = -WTSFC / USTAR
          QSTAR = -WQSFC / USTAR
          TMEAN = RATIO * TSTAR
          QMEAN = RATIO * QSTAR
          TSURF = T1 - TMEAN
          QSURF = Q1 - QMEAN
        Else
          TMEAN = T1 - TSURF
          QMEAN = Q1 - QSURF
          TSTAR = TMEAN / RATIO
          QSTAR = QMEAN / RATIO
          WTSFC = -USTAR * TSTAR
          WQSFC = -USTAR * QSTAR
        End If

C        CALCUATE B* AND C* OR (B1-B0) AND (C1-C0), DEPENDING ON
C        ISURFB AND ISURFC
        If (ISURFB.EQ.0) Then
          BSTAR = -WBSFC / USTAR
          BMEAN = RATIO * BSTAR
          BSURF = B1 - BMEAN
        Else
          BMEAN = B1 - BSURF
          BSTAR = BMEAN / RATIO
          WBSFC = -USTAR * BSTAR
        End If

        If (ISURFC.EQ.0) Then
          CSTAR = -WCSFC / USTAR
          CMEAN = RATIO * CSTAR
          CSURF = C1 - CMEAN
        Else
          CMEAN = C1 - CSURF
          CSTAR = CMEAN / RATIO
          WCSFC = -USTAR * CSTAR
        End If

      Else

C//////////////////////////////////////////////////////////////////////
C     IF FIRST MODEL GRID POINT IS ABOVE 50 M THEN USE PBL SIMILARITY
C     THEORY AS FOR THE MIXED LAYER MODEL
C//////////////////////////////////////////////////////////////////////

C        CHECK TO SEE IF CONDITIONS ARE NEUTRAL/NEARLY NEUTRAL      
        ZETA = -Z1 * VK * GRAV * WTSFC / (TREF*USTAR**3)
	 
C        IF ZETA IS VERY SMALL, DON'T DO NON-LINEAR ITERATION 
C        NOTE: 1.0e-4 IS ARBITRARY        
        If (ABS(ZETA).GT.1.0E-4) Then
      
C           PUT A SIMPLE LOWER LIMIT ON ustar FOR FREE CONVECTION CASE
          If (USTAR.LT.0.01) USTAR = 0.01

C           FIRST BRACKET THE ROOT 
          USTAR1 = 0.0
          USTAR2 = 0.0
          Call ZBRAK(USTARPBL_RES,USTAR*0.5,USTAR*1.5,100,USTAR1,USTAR2,
     *        1,USTAR,WIND,TREF,GRAV,WTSFC,Z1,Z0M,VK,ZI)
     
C           CHECK FOR BRACKETED ROOT. 
          If (USTAR1.LT.1.0E-4.AND.USTAR2.LT.1.0E-4) Then
            USTAR1 = USTAR * 0.5
            USTAR2 = USTAR * 1.5
          End If
     
C           FIND ROOT (ustar)
          USTAR = ZBRENT(USTARPBL_RES,USTAR1,USTAR2,1.0E-6,USTAR,WIND,
     *        TREF,GRAV,WTSFC,Z1,Z0M,VK,ZI)


        End If
	 
C        PUT LOWER LIMIT ON ustar
C         ustar = amax1(ustar, vk * u_free / (alog((z1+z0m)/z0m)))


C        CALCULATE NEUTRAL EXCHANGE COEFFICIENTS
        CUN = 1.0 / (ALOG(0.025*ZI/Z0M)/VK+8.4)
        CTHETAN = 1.0 / (0.74/VK*ALOG(0.025*ZI/Z0M)+7.3)
        CU = USTAR / WIND
        CTHETA = 1.0 / (1.0/CTHETAN+1.0/CU-1.0/CUN)

C        CALCULATE T* AND Q* OR (T1-T0) AND (Q1-Q0), DEPENDING ON ISURF
        If (ISURF.EQ.0) Then
          TSTAR = -WTSFC / USTAR
          QSTAR = -WQSFC / USTAR
          TMEAN = TSTAR / CTHETA
          QMEAN = QSTAR / CTHETA
          TSURF = T1 - TMEAN
          QSURF = Q1 - QMEAN
        Else
          TMEAN = T1 - TSURF
          QMEAN = Q1 - QSURF
          TSTAR = TMEAN * CTHETA
          QSTAR = QMEAN * CTHETA
          WTSFC = -USTAR * TSTAR
          WQSFC = -USTAR * QSTAR
        End If

C        CALCUATE B* AND C* OR (B1-B0) AND (C1-C0), DEPENDING ON
C        ISURFB AND ISURFC
        If (ISURFB.EQ.0) Then
          BSTAR = -WBSFC / USTAR
          BMEAN = BSTAR / CTHETA
          BSURF = B1 - BMEAN
        Else
          BMEAN = B1 - BSURF
          BSTAR = BMEAN * CTHETA
          WBSFC = -USTAR * BSTAR
        End If

        If (ISURFC.EQ.0) Then
          CSTAR = -WCSFC / USTAR
          CMEAN = RATIO / CTHETA
          CSURF = C1 - CMEAN
        Else
          CMEAN = C1 - CSURF
          CSTAR = CMEAN * CTHETA
          WCSFC = -USTAR * CSTAR
        End If

      End If

C     CALCULATE ANEMOMETER VALUES OF U, V, Q, T, B, C
      ZETA = -ZA * VK * GRAV * WTSFC / (TREF*USTAR**3)
CTEMP     
C      psi_m = buspsi_m1(zeta)
      PSI_M = FPSI_M(ZETA)
CTEMP
C      psi_h = buspsi_h2(zeta)
      PSI_H = FPSI_H(ZETA)
      RATIO = PROFILE(VK,ZA,Z0H,PSI_H)

      WSPDA = USTAR * (ALOG(ZA/Z0M)-PSI_M) / VK
      TA = TSTAR * RATIO + TSURF
      QA = QSTAR * RATIO + QSURF
      BA = BSTAR * RATIO + BSURF
      CA = CSTAR * RATIO + CSURF

      TEP = U1 / WIND
      TEP = AMIN1(TEP,1.)
      TEP = AMAX1(TEP,-1.)
      THTA = ACOS(TEP)
      WUSFC = -(USTAR**2) * COS(THTA)
      WVSFC = -(USTAR**2) * SIN(THTA) * SIGN(1.,V1)
      UA = WSPDA * COS(THTA)
      VA = WSPDA * SIN(THTA) * SIGN(1.,V1)

      Return
      End

      Subroutine MIX_SURFACE(U1,V1,T1,Q1,B1,C1,TSURF,QSURF,BSURF,CSURF,
     *    USTAR,TSTAR,QSTAR,BSTAR,CSTAR,Z1,ZA,Z0M,Z0H,Z0,UA,VA,TA,QA,BA,
     *    CA,ISURF,ISURFB,ISURFC,WUSFC,WVSFC,WTSFC,WQSFC,WBSFC,WCSFC,
     *    GRAV,TREF,VK,DEL_Z,ZI,Z1_FLAG)

      Implicit None

C***********************************************************************
C
C     CALCULATES SURFACE BOUNDARY CONDITIONS USING M-O SIMILARITY
C     AND THE PROFILE RELATIONS OF BUSINGER ET AL. (1971)
C
C     INPUTS:  1) MEAN QUANTITES AT 1ST GRIDPOINT.
C              2) ROUGHNESS LENGTH.
C              3) SURFACE FLUXES OR SURFACE VALUES OF TEMPERATURE,
C                 HUMIDITY, AND SCALARS (DEPENDING ON SURFACE BOUNDARY
C                 CONDITION - CONSTANT FLUX OR CONSTANT VALUE AT SURFACE.
C              4) INTEGERS WHICH SPECIFY WHAT BOUNDARY CONDITION TO
C                 USE (ISURF FOR TEMPERATURE AND HUMIDITY, ISURFB FOR
C                 SCALAR B, AND ISURFC FOR SCALAR C.
C              5) GRAVITY, REFERENCE TEMPERATURE, VON KARMAN'S
C                 CONSTANT
C
C     OUTPUTS: 1) SURFACE LAYER SCALES (USTAR, TSTAR, QSTAR, BSTAR,
C                 CSTAR)
C              2) SURFACE FLUXES OR SURFACE VALUES (DEPENDING ON
C                 BOUNDARY CONDITION SPECIFIED)
C              3) ANEMOMETER HEIGHT VALUES (NEEDED FOR KPROF SCHEME)
C
C**********************************************************************

C     MEAN QUANTITIES AT THE LOWEST GRIDPOINT
      Real U1, V1, T1, Q1, B1, C1

C     MEAN QUANTITIES AT THE SURFACE
      Real TSURF, QSURF, BSURF, CSURF

C     SURFACE LAYER SCALES
      Real USTAR, TSTAR, QSTAR, BSTAR, CSTAR, U_FREE, SMAL

C     HEIGHT OF LOWEST GRIDPOINT, ANEMOMETER HEIGHT, AND ROUGHNESS
C     HEIGHTS
      Real Z1, ZA, Z0, Z0M, Z0H, ZI, DEL_Z

C     MEAN QUANTITIES AT THE ANEMOMETER HEIGHT
      Real UA, VA, TA, QA, BA, CA, WSPDA

C     SURFACE FLUXES
      Real WUSFC, WVSFC, WTSFC, WQSFC, WBSFC, WCSFC

C     VALUES AT THE LOWEST GRIDPOINT MINUS VALUES AT THE
C     SURFACE
      Real TMEAN, QMEAN, BMEAN, CMEAN

C     NUMERICAL RECIPES SECANT ROOT FINDING ROUTINE AND RESIDUAL FUNCTION
      Real ZBRENT, USTARPBL_RES, ZBRAK, Z1_RES
      External USTARPBL_RES, Z1_RES

      Real GRAV, TREF, VK, WIND, XM1, ZETA, PSI_M, PSI_H, RATIO
      Real XM1OLD, THTA, TEP, BUSPSI_M1, BUSPSI_H2, PROFILE
      Real FPSI_M, FPSI_H
      Real ZM1, ZM2, Z
      Real USTAR1, USTAR2

C     EXCHANGE COEFFICIENTS
      Real CTHETA, CU, CUN, CTHETAN

C     FLAGS
      Integer ISURF, ISURFB, ISURFC, ITEMP, ITEMPB, Z1_FLAG
      Integer ITEMPC, ITER

      WIND = SQRT((U1**2)+(V1**2))
      Z0M = Z0
      Z0H = Z0
      SMAL = 0.00001

C     SET LOWER LIMIT ON wind TO GIVE u_free FOR THE FREE
C     CONVECTION CASE
      u_free = 0.0
      if(abs(wtsfc) .gt. smal) then
         U_FREE = 0.07 * (GRAV/TREF*WTSFC*DEL_Z) ** (1.0/3.0)
      endif
      WIND = AMAX1(WIND,U_FREE)
      write(6,6001) wind,u_free,ustar
6001  format('6001, wind = ',e15.6,/,
     +       '      u_free = ',e15.6,/,
     +       '      ustar = ',e15.6)

C     FIRST BRACKET THE ROOT
c     Call ZBRAK(USTARPBL_RES,0.75*USTAR,USTAR*3.0,100,USTAR1,USTAR2,1,
      Call ZBRAK(USTARPBL_RES,USTAR/2.0,USTAR*3.0,100,USTAR1,USTAR2,1,
     *    USTAR,WIND,TREF,GRAV,WTSFC,Z1,Z0M,VK,ZI)
      write(6,6002) wind,u_free,ustar,ustar1,ustar2
6002  format('6002, wind = ',e15.6,/,
     +       '      u_free = ',e15.6,/,
     +       '      ustar = ',e15.6,/,
     +       '      ustar1 = ',e15.6,/,
     +       '      ustar2 = ',e15.6)
C     FIND ROOT (ustar)
      USTAR = ZBRENT(USTARPBL_RES,USTAR1,USTAR2,1.0E-6,USTAR,WIND,TREF,
     *    GRAV,WTSFC,Z1,Z0M,VK,ZI)
      write(6,44001) ustar
44001 format('44001, ustar = ',e15.6)

C     CALCULATE NEUTRAL EXCHANGE COEFFICIENTS
      CUN = 1.0 / (ALOG(0.025*ZI/Z0M)/VK+8.4)
      CTHETAN = 1.0 / (0.74/VK*ALOG(0.025*ZI/Z0M)+7.3)
      CU = USTAR / WIND
      CTHETA = 1.0 / (1.0/CTHETAN+1.0/CU-1.0/CUN)
     
C     CALCULATE T* AND Q* OR (T1-T0) AND (Q1-Q0), DEPENDING ON ISURF
      If (ISURF.EQ.0) Then
        TSTAR = -WTSFC / USTAR
        QSTAR = -WQSFC / USTAR
        TMEAN = TSTAR / CTHETA
        QMEAN = QSTAR / CTHETA
        TSURF = T1 - TMEAN
        QSURF = Q1 - QMEAN
      Else
        TMEAN = T1 - TSURF
        QMEAN = Q1 - QSURF
        TSTAR = TMEAN * CTHETA
        QSTAR = QMEAN * CTHETA
        WTSFC = -USTAR * TSTAR
        WQSFC = -USTAR * QSTAR
      End If

C     CALCUATE B* AND C* OR (B1-B0) AND (C1-C0), DEPENDING ON
C     ISURFB AND ISURFC
      If (ISURFB.EQ.0) Then
        BSTAR = -WBSFC / USTAR
        BMEAN = BSTAR / CTHETA
        BSURF = B1 - BMEAN
      Else
        BMEAN = B1 - BSURF
        BSTAR = BMEAN * CTHETA
        WBSFC = -USTAR * BSTAR
      End If

      If (ISURFC.EQ.0) Then
        CSTAR = -WCSFC / USTAR
        CMEAN = RATIO / CTHETA
        CSURF = C1 - CMEAN
      Else
        CMEAN = C1 - CSURF
        CSTAR = CMEAN * CTHETA
        WCSFC = -USTAR * CSTAR
      End If

C     CALCULATE ANEMOMETER VALUES OF U, V, Q, T, B, C
      ZETA = -ZA * VK * GRAV * WTSFC / (TREF*USTAR**3)
CTEMP      
C      psi_m = buspsi_m1(zeta)
      PSI_M = FPSI_M(ZETA)
CTEMP
C      psi_h = buspsi_h2(zeta)
      PSI_H = FPSI_H(ZETA)
      RATIO = PROFILE(VK,ZA,Z0H,PSI_H)

      WSPDA = USTAR * (ALOG(ZA/Z0M)-PSI_M) / VK
      TA = TSTAR * RATIO + TSURF
      QA = QSTAR * RATIO + QSURF
      BA = BSTAR * RATIO + BSURF
      CA = CSTAR * RATIO + CSURF

      TEP = U1 / WIND
      TEP = AMIN1(TEP,1.)
      TEP = AMAX1(TEP,-1.)
      THTA = ACOS(TEP)
      WUSFC = -(USTAR**2) * COS(THTA)
      WVSFC = -(USTAR**2) * SIN(THTA) * SIGN(1.,V1)
      UA = WSPDA * COS(THTA)
      VA = WSPDA * SIN(THTA) * SIGN(1.,V1)

      Return
      End

      Function PROFILE(VK,Z,Z0H,PSI_H)
      Implicit None

C*********************************************************************
C     CALCULATES M-O SIMILARITY PROFILE RELATION FOR TEMPERATURE
C     AND SCALARS
C**********************************************************************
      Real Z, Z0H, PSI_H, VK, PROFILE
      PROFILE = 0.74 * (ALOG(Z/Z0H)-PSI_H) / VK
      Return
      End

      Subroutine INTERP(ZM,ZT,DELZ,NKLES,NK,U,V,THETA,Q,B,C,UG,VG,ZLES,
     *    ULES,VLES,THETAL,QLES,BLES,CLES,UGLES,VGLES,WU0,WV0,WT0,WQ0,
     *    WB0,WC0,GRAV,TREF,VK,Z0,CASE,UGSURF,VGSURF,MODEL,USTAR,TSURF,
     *    BSURF,CSURF)

      Implicit None

C***********************************************************************
C
C     CALCULATES MEAN AND TURBULENT GRIDS BASED ON SPECIFIED DELTA-Z
C     AND INTERPOLATES MEAN QUANTITIES TO THIS GRID
C
C     INPUTS:  1) SPECIFIED DELTA-Z
C              2) LES GRID
C              3) MEAN QUANTITIES ON LES GRID
C              4) GEOSTROPHIC WIND ON LES GRID
C              5) INITIAL SURFACE FLUXES FROM LES
C              6) ROUGHNESS LENGTH, GRAVITY, REFERENCE TEMPERATURE,
C                 VON KARMAN'S CONSTANT
C              7) PBL MODEL AND CASE NAME
C              8) CONVECTIVE VELOCITY
C
C     OUTPUT: 1) MEAN AND TURBULENT PBL GRIDS
C             2) VALUES OF MEAN QUANTITIES AND GEOSTROPHIC WIND 
C                INTERPOLATED TO THE NEW GRID
C
C
C NOTE:  SINCE THE STABILITY BOUNDS SCHEME IS A MIXED-LAYER MODEL,
C        THE FIRST GRID POINT MUST NOT BE IN THE SURFACE LAYER. THUS,
C        THE FIRST GRIDPOINT IS SET NO LOWER THAN 60 M FOR THE STABILITY
C        BOUNDS SCHEME.
C        
C************************************************************************

      Integer K, NK, NKLES, J, JJ, ITOPM, ITOPT, ISTART

      Parameter (ITOPM = 250)
      Parameter (ITOPT = 249)

C     MEAN AND TURBULENT PBL GRID, DELTA-Z, ROUGHNESS LENGTH
      Real ZM(ITOPM), ZT(ITOPT), DELZ, Z0

C     MEAN PBL VARIABLES
      Real U(ITOPM), V(ITOPM), THETA(ITOPM), Q(ITOPM), B(ITOPM), 
     *    C(ITOPM)

C     GEOSTROPHIC WIND FOR PBL MODELS
      Real UG(ITOPM), VG(ITOPM)
      
C     LES GRID AND  LES DELTA-Z
      Real ZLES(NKLES), DZLES

C     LES MEAN VARIABLES
      Real ULES(NKLES), VLES(NKLES), THETAL(NKLES), QLES(NKLES), 
     *    BLES(NKLES), CLES(NKLES)

C     LES GEOSTROPHIC WIND
      Real UGLES(NKLES), VGLES(NKLES)

C     SURFACE VALUES
      Real TSURF, QSURF, BSURF, CSURF, UGSURF, VGSURF, ZSURF

C     LES SURFACE FLUXES
      Real WU0, WV0, WT0, WQ0, WC0, WB0

C     SCALING VARIABLES
      Real USTAR, TSTAR, QSTAR, BSTAR, CSTAR

C     FUNCTIONS CALLED
      Real LOGINT, XSURF, PSI_M, PSI_H, RATIO, PROFILE, BUSPSI_M1, 
     *    BUSPSI_H2, LININT
      Real FPSI_M, FPSI_H

C     GRAVITY, REFERENCE TEMPERATURE, VON KARMAN'S CONSTANT
      Real GRAV, TREF, VK
      
C     Z/MONIN LENGTH
      Real ZETA

C     MISC VARIABLES 
      Real WIND, TEP, THTA
      Character*80 CASE, MODEL

C     AVERAGING VARIABLES
      Integer K_POINTER, K_BEGIN, K_END
      Real Z_BEGIN, Z_END, W1
      Real U_M1, V_M1, THETA_M1, Q_M1, B_M1, C_M1

C     SET FIRST GRIDPOINT
      If (MODEL.NE.'STABIL'.OR.DELZ.GE.60.0) Then
        ZM(1) = DELZ / 2.0
      Else
        ZM(1) = 60.0
      End If

C     CALCULATE THE MEAN GRID
      Do 10 K = 2, ITOPM
        ZM(K) = ZM(K-1) + DELZ
        If ((ZM(K)+DELZ/2.0).GT.ZLES(NKLES)) Go To 20
        If (K.EQ.ITOPM) Then
          Write (*,*) 'please use a larger delta-z.  not enough space   
     *        allocated for arrays'
          Stop
        End If
   10 Continue

C     REDEFINE TOP GRID POINT AS THE NEAREST LEVEL BELOW OR EQUAL TO THE
C     THE TOP LES GRID
   20 NK = K - 1

C     CALCULATE THE TURBULENT GRID
      Do 30 J = 1, NK - 1
        ZT(J) = .5 * (ZM(J+1)+ZM(J))
   30 Continue

C     INTERPOLATE QUANTITIES TO THE MEAN GRID
      DZLES = ZLES(2) - ZLES(1)
      ISTART = 1

C     CALCULATE QUANTITES AT LEVELS ABOVE THE FIRST LES GRIDPOINT   
      QSURF = 0.0
      Do 50 J = ISTART, NK
        K = 1
   40   If (ZLES(K).GT.ZM(J)) Then

C           GEOSTROPHIC WIND LINEARLY INTERPOLATED
          UG(J) = LININT(ZM(J),ZLES(K),ZLES(K-1),UGLES(K),UGLES(K-1))
          VG(J) = LININT(ZM(J),ZLES(K),ZLES(K-1),VGLES(K),VGLES(K-1))
        Else
          K = K + 1
          Go To 40
        End If
   50 Continue

C     MAIN AVERAGING LOOP
      K_BEGIN = 0
      K_END = 0
      Z_BEGIN = ZM(1) / 2.0
      Z_END = Z_BEGIN
      K_POINTER = 1

C     FIND THE LOWER GRID POINT LIMIT ON THE FIRST PASS
   60 If (.NOT.((ZLES(K_POINTER).GT.Z_BEGIN).OR.(K_POINTER.EQ.NKLES))) 
     *    Then
        K_POINTER = K_POINTER + 1
        Go To 60
      End If

C     SET k_begin TO BE INCREMENTED BY 1 BEFORE USE       
      K_END = K_POINTER - 2

C     RESET k_pointer
      K_POINTER = 0

      Do 90 J = ISTART, NK
        K_BEGIN = K_END + 1
        Z_BEGIN = Z_END
        Z_END = (ZM(J)+ZM(J+1)) / 2.0
C        FIND THE GRID POINT LIMITS OF POINTS TO INTERPOLATE
   70   If (.NOT.((ZLES(K_POINTER).GT.Z_END).OR.(K_POINTER.EQ.NKLES))) 
     *      Then
          K_POINTER = K_POINTER + 1
          Go To 70
        End If
        K_END = K_POINTER - 1

C        SET UP POINT BELOW LOWER LES POINT IN INTEGRATION
        U_M1 = ULES(K_BEGIN)
        V_M1 = VLES(K_BEGIN)
        THETA_M1 = THETAL(K_BEGIN)
        Q_M1 = QLES(K_BEGIN)
        B_M1 = BLES(K_BEGIN)
        C_M1 = CLES(K_BEGIN)
	    
C        k_begin AND k_end NOW POINT TO THE BEGINING AND THE END OF
C        THE POINTS WHICH ARE TO BE INCLUDED IN THE AVERAGED MEAN VALUES.
C        MEAN VALUE IS NOW CALCULATED USING WEIGHTING

C        CALCULATE WEIGHTED CONTRIBUTION FROM FIRST POINT
        W1 = (ZLES(K_BEGIN)-Z_BEGIN) / (Z_END-Z_BEGIN)
        U(J) = W1 * (ULES(K_BEGIN)+U_M1) / 2.0
        V(J) = W1 * (VLES(K_BEGIN)+V_M1) / 2.0
        THETA(J) = W1 * (THETAL(K_BEGIN)+THETA_M1) / 2.0
        Q(J) = W1 * (QLES(K_BEGIN)+Q_M1) / 2.0
        B(J) = W1 * (BLES(K_BEGIN)+B_M1) / 2.0
        C(J) = W1 * (CLES(K_BEGIN)+C_M1) / 2.0

C        CALCULATE CONTRIBUTION FROM INTERIOR POINTS
        Do 80 K = K_BEGIN + 1, K_END

C           SET UP POINT BELOW LOWER LES POINT IN INTEGRATION
          U_M1 = ULES(K-1)
          V_M1 = VLES(K-1)
          THETA_M1 = THETAL(K-1)
          Q_M1 = QLES(K-1)
          B_M1 = BLES(K-1)
          C_M1 = CLES(K-1)

          W1 = (ZLES(K)-ZLES(K-1)) / (Z_END-Z_BEGIN)
          U(J) = U(J) + W1 * (ULES(K)+U_M1) / 2.0
          V(J) = V(J) + W1 * (VLES(K)+V_M1) / 2.0
          THETA(J) = THETA(J) + W1 * (THETAL(K)+THETA_M1) / 2.0
          Q(J) = Q(J) + W1 * (QLES(K)+Q_M1) / 2.0
          B(J) = B(J) + W1 * (BLES(K)+B_M1) / 2.0
          C(J) = C(J) + W1 * (CLES(K)+C_M1) / 2.0
   80   Continue

C        CALCULATE THE WEIGHTED CONTRIBUTION FROM THE LAST POINT
        W1 = (Z_END-ZLES(K_END)) / (Z_END-Z_BEGIN)
        U(J) = U(J) + W1 * (ULES(K_END)+ULES(K_END+1)) / 2.0
        V(J) = V(J) + W1 * (VLES(K_END)+VLES(K_END+1)) / 2.0
        THETA(J) = THETA(J) + W1 * (THETAL(K_END)+THETAL(K_END+1)) / 2.0
        Q(J) = Q(J) + W1 * (QLES(K_END)+QLES(K_END+1)) / 2.0
        B(J) = B(J) + W1 * (BLES(K_END)+BLES(K_END+1)) / 2.0
        C(J) = C(J) + W1 * (CLES(K_END)+CLES(K_END+1)) / 2.0

   90 Continue

      Return
      End


      Function LOGINT(ZMID,ZTOP,ZBOT,VTOP,VBOT)
      Implicit None
C*********************************************************************
C     LOGARITHMICALLY INTERPOLATES QUANTITIIES TO A HEIGHT
C     `ZMID' GIVEN THE VALUE `VTOP' AT A HEIGHT `ZTOP' AND
C      A VALUES `VBOT' AT A HEIGHT `ZBOT'
C*********************************************************************
      Real LOGINT, ZMID, ZBOT, ZTOP, VTOP, VBOT
      LOGINT = ((ALOG(ZMID)-ALOG(ZBOT))/(ALOG(ZTOP)-ALOG(ZBOT))) * (VTOP
     *    -VBOT) + VBOT
      Return
      End


      Function LININT(ZMID,ZTOP,ZBOT,VTOP,VBOT)
      Implicit None
C*********************************************************************
C     LINEARLY INTERPOLATES QUANTITIIES TO A HEIGHT `ZMID'
C     GIVEN THE VALUES `VTOP' AT A HEIGHT `ZTOP' AND `VBOT' 
C     AT A HEIGHT `ZBOT'
C*********************************************************************
      Real LININT, ZMID, ZBOT, ZTOP, VTOP, VBOT
      LININT = ((ZMID-ZBOT)/(ZTOP-ZBOT)) * (VTOP-VBOT) + VBOT
      Return
      End

      Function XSURF(WU0,WV0,WT0,WX0,G,TREF,Z1LES,VK,Z0,X1,USTAR,XSTAR,
     *    ZETA)

      Implicit None
C********************************************************************
C     CALCULATES SURFACE VALUES FROM M-O SIMILARITY AND RELATIONS
C     OF BUSINGER ET AL (1971) GIVEN SURFACE FLUXES AND VALUES AT
C     FIRST GRIDPOINT
C********************************************************************

C     LES INITIAL SURFACE FLUXES
      Real WU0, WV0, WT0, WX0, USTAR, XSTAR, G, TREF
      Real ZETA, Z1LES, VK, PROFILE, PSI_H, RATIO, BUSPSI_H2
      Real FPSI_H
      Real XSURF, X1, Z0

C     NOTE: ustar IS CALCULATED ELSEWHERE AND PASSED INTO xsurf

      XSTAR = -WX0 / USTAR
      ZETA = -G / TREF * WT0 * Z1LES * VK / USTAR ** 3
      ZETA = AMAX1(-2.0,ZETA)

C      psi_h=buspsi_h2(zeta)
      PSI_H = FPSI_H(ZETA)
      
      RATIO = PROFILE(VK,Z1LES,Z0,PSI_H)
      XSURF = X1 - RATIO * XSTAR
      Return
      End


      Subroutine MIXINIT(ZM,ZT,NKLES,NK,U,V,THETA,Q,B,C,UG,VG,ZLES,ULES,
     *    VLES,THETAL,QLES,BLES,CLES,UGLES,VGLES,PBLHGT,LESZI,USLOPE,
     *    VSLOPE,TSLOPE,QSLOPE,BSLOPE,CSLOPE,UB,VB,TB,QB,BB,CB)
     
      Implicit None
C***********************************************************************
C
C     INITIALIZES QUANTITIES NEEEDED BY MIXED LAYER MODEL
C
C     INPUT:  1) LES GRID
C             2) VALUES OF LES MEAN QUANTITES
C             3) LES PBL HEIGHT
C
C     OUTPUT: 1) PBL HEIGHT
C             2) MIXED-LAYER AVERAGES OF MEAN QUANTITES
C             3) SLOPE AND Y-INTERCEPT OF MEAN QUANTITES ABOVE INVERSION
C                (QUANTITIES ASSUMED LINEAR ABOVE INVERSION)
C             4) PLOTTING GRID FOR PBL QUANTITIES (NOTE - THIS GRID USED
C                ONLY FOR PLOTTING MEAN QUANTITIES IN PBLPLOT; GRID NOT
C                USED IN MIXED LAYER SCHEME EXCEPT TO CREATE INITIAL
C                MIXED-LAYER AVERAGES)
C             5) VALUES OF MEAN QUANTITES ON THE PLOTTING GRID
C
C*************************************************************************
      Integer ITOP, NKLES, NK, J, K
      Parameter (ITOP = 250)
      
C     MEAN AND TURBULENT PLOTTING GRIDS, DELTA-Z
      Real ZM(ITOP), ZT(ITOP), DELZ

C     MEAN PBL VARIABLES
      Real U(ITOP), V(ITOP), THETA(ITOP), Q(ITOP), B(ITOP), C(ITOP)

C     GEOSTROPHIC WIND FOR PBL MODELS
      Real UG(ITOP), VG(ITOP)
      
C     LES GRID 
      Real ZLES(NKLES)

C     LES MEAN VARIABLES
      Real ULES(NKLES), VLES(NKLES), THETAL(NKLES), QLES(NKLES), 
     *    BLES(NKLES), CLES(NKLES)
      
C     LES GEOSTROPHIC WIND
      Real UGLES(NKLES), VGLES(NKLES)

C     FUNCTIONS CALLED
      Real LOGINT, LININT

C     LES PBL HEIGHT, PBL HEIGHT
      Real LESZI, PBLHGT

C     FIRST LEVEL BELOW PBL HEIGHT
      Integer IPBL

C     MIXED-LAYER AVERAGED QUANTITIES
      Real UMEAN, VMEAN, TMEAN, QMEAN, BMEAN, CMEAN

C     SLOPE OF MEAN VARIABLES ABOVE INVERSION
      Real USLOPE, VSLOPE, TSLOPE, QSLOPE, BSLOPE, CSLOPE

C     Y-INTERCEPT OF MEAN VARIABLES ABOVE INVERSION
      Real UB, VB, TB, QB, BB, CB

C     USED IN CALCULATION OF MIXED-LAYER AVERAGES
      Real USUM, VSUM, TSUM, QSUM, BSUM, CSUM

C     POTENTIAL TEMPERATURE JUST ABOVE INVERSION, TEMPERATURE
C     JUMP AT INVERSION
      Real TTOP, DTHETA

C     SET PBLHEIGHT TO LES VALUES
      PBLHGT = LESZI

c     CALCULATE THE PBL GRIDS
      DELZ = 50.0
      ZM(1) = 50.0
      ZT(1) = ZM(1)
      Do 10 K = 2, ITOP
        ZM(K) = ZM(K-1) + DELZ
        ZT(K) = ZM(K)
        If (ZM(K).GT.ZLES(NKLES)) Go To 20
   10 Continue

c     REDEFINE TOP GRID AS THE NEAREST LEVEL BELOW OR EQUAL TO THE
c     THE TOP LES GRID
   20 NK = K - 1

 
C     INTERPOLATE PBL QUANTITIES TO THE NEW GRID
      Do 40 J = 1, NK
        K = 2
   30   If (ZLES(K).GT.ZM(J)) Then
          U(J) = LOGINT(ZM(J),ZLES(K),ZLES(K-1),ULES(K),ULES(K-1))
          V(J) = LOGINT(ZM(J),ZLES(K),ZLES(K-1),VLES(K),VLES(K-1))
          THETA(J) = LOGINT(ZM(J),ZLES(K),ZLES(K-1),THETAL(K),
     *        THETAL(K-1))
          Q(J) = LOGINT(ZM(J),ZLES(K),ZLES(K-1),QLES(K),QLES(K-1))
          B(J) = LOGINT(ZM(J),ZLES(K),ZLES(K-1),BLES(K),BLES(K-1))
          C(J) = LOGINT(ZM(J),ZLES(K),ZLES(K-1),CLES(K),CLES(K-1))
          UG(J) = LININT(ZM(J),ZLES(K),ZLES(K-1),UGLES(K),UGLES(K-1))
          VG(J) = LININT(ZM(J),ZLES(K),ZLES(K-1),VGLES(K),VGLES(K-1))
        Else
          K = K + 1
          Go To 30
        End If
   40 Continue

C     FIND LEVEL CORRESPONDING TO PBL HEIGHT
      Do 50 J = 1, NK
        If (ZM(J).GT.LESZI) Then
          IPBL = (J-1)
          Go To 60
        End If
   50 Continue

   60 Continue
C     CALCULATE MIXED-LAYER AVERAGES
      USUM = 0.0
      VSUM = 0.0
      TSUM = 0.0
      QSUM = 0.0
      BSUM = 0.0
      CSUM = 0.0
      Do 70 J = 1, IPBL - 1
        USUM = USUM + U(J) * DELZ
        VSUM = VSUM + V(J) * DELZ
        TSUM = TSUM + THETA(J) * DELZ
        QSUM = QSUM + Q(J) * DELZ
        BSUM = BSUM + B(J) * DELZ
        CSUM = CSUM + C(J) * DELZ
   70 Continue
      USUM = USUM + U(IPBL) * (PBLHGT-ZM(IPBL))
      VSUM = VSUM + V(IPBL) * (PBLHGT-ZM(IPBL))
      TSUM = TSUM + THETA(IPBL) * (PBLHGT-ZM(IPBL))
      QSUM = QSUM + Q(IPBL) * (PBLHGT-ZM(IPBL))
      BSUM = BSUM + B(IPBL) * (PBLHGT-ZM(IPBL))
      CSUM = CSUM + C(IPBL) * (PBLHGT-ZM(IPBL))

      UMEAN = USUM / (PBLHGT-ZM(1))
      VMEAN = VSUM / (PBLHGT-ZM(1))
      TMEAN = TSUM / (PBLHGT-ZM(1))
      QMEAN = QSUM / (PBLHGT-ZM(1))
      BMEAN = BSUM / (PBLHGT-ZM(1))
      CMEAN = CSUM / (PBLHGT-ZM(1))

C     CALCULATE SLOPE AND INTERCEPT OF MEAN VARIABLES ABOVE INVERSION
      USLOPE = (U(NK)-U(NK-5)) / (ZM(NK)-ZM(NK-5))
      UB = U(NK-2) - USLOPE * ZM(NK-2)
      VSLOPE = (V(NK)-V(NK-5)) / (ZM(NK)-ZM(NK-5))
      VB = V(NK-2) - VSLOPE * ZM(NK-2)
      TSLOPE = (THETA(NK)-THETA(NK-5)) / (ZM(NK)-ZM(NK-5))
      TB = THETA(NK-2) - TSLOPE * ZM(NK-2)
      QSLOPE = (Q(NK)-Q(NK-5)) / (ZM(NK)-ZM(NK-5))
      QB = Q(NK-2) - QSLOPE * ZM(NK-2)
      BSLOPE = (B(NK)-B(NK-5)) / (ZM(NK)-ZM(NK-5))
      BB = B(NK-2) - BSLOPE * ZM(NK-2)
      CSLOPE = (C(NK)-C(NK-5)) / (ZM(NK)-ZM(NK-5))
      CB = C(NK-2) - CSLOPE * ZM(NK-2)
     
C     CALCULATE JUMP IN POTENTIAL TEMP. AT INVERSION
      TTOP = TSLOPE * PBLHGT + TB
      DTHETA = TTOP - TMEAN

C     MAKE SURE DTHETA DOES NOT START OUT NEGATIVE
      If (DTHETA.LT.0) Then
        PBLHGT = (TMEAN-TB) / TSLOPE + 5.0

C        FIND LEVEL CORRESPONDING TO PBL HEIGHT
        Do 80 J = 1, NK
          If (ZM(J).GT.PBLHGT) Then
            IPBL = (J-1)
            Go To 90
          End If
   80   Continue
   90   Continue
      End If

C     PUT QUANTITIES INTO PBL GRID FOR PLOTTING PURPOSES
C     VALUES IN MIXED LAYER ARE EQUAL TO THEIR AVERAGE
      Do 100 J = 1, IPBL
        U(J) = UMEAN
        V(J) = VMEAN
        THETA(J) = TMEAN
        Q(J) = QMEAN
        B(J) = BMEAN
        C(J) = CMEAN
  100 Continue

C     VALUES ABOVE INVERSION ARE CALCULATED FROM THEIR SLOPE AND
C     INTERCEPT
      Do 110 J = IPBL + 1, NK
        U(J) = USLOPE * ZM(J) + UB
        V(J) = VSLOPE * ZM(J) + VB
        THETA(J) = TSLOPE * ZM(J) + TB
        Q(J) = QSLOPE * ZM(J) + QB
        B(J) = BSLOPE * ZM(J) + BB
        C(J) = CSLOPE * ZM(J) + CB
  110 Continue

      Return
      End




C***********************************************************************
C
C     CALCULATION OF STABILITY FUNCTION psi_m
C     REFERENCE: ARYA P. 167
C
C**********************************************************************
      Function FPSI_M(ZETA)
      Implicit None
      Real ZETA, FPSI_M, PI, X

      Data PI /3.1415926/
      If (ZETA.LT.0.) Then
        X = (1.0-15.0*ZETA) ** 0.2
        FPSI_M = ALOG(((1.0+X*X)/2.0)*((1.0+X)/2.0)**2) - ATAN(X) + PI /
     *      2.0
      Else
        FPSI_M = -4.7 * ZETA
      End If
      
      Return
      End
      
C***********************************************************************
C
C     CALCULATION OF STABILITY FUNCTION psi_h
C     REFERENCE: ARYA P. 167
C
C**********************************************************************
      Function FPSI_H(ZETA)
      Implicit None
      Real ZETA, FPSI_H, X

      If (ZETA.LT.0.) Then
        X = SQRT(1.0-15.0*ZETA)
        FPSI_H = 2.0 * ALOG((1.0+X**2)/2.0)
      Else
        FPSI_H = -4.7 * ZETA
      End If
      Return
      End








      Function BUSPSI_M1(ZETA)
      Implicit None
      Real ZETA, BUSPSI_M1, PII, X
C***********************************************************************
C
C     PURPOSE: TO CALCULATE STABILITY CORRECTION FUNCTION IN THE 
C               SURFACE LAYER WIND PROFILE.
C
C     ARGUMENTS:
C
C          INPUT:
C                 -ETA: STABILITY PARAMETER (Z/L)
C
C         OUTPUT:
C                 -BUSFIM1 :CORRECTION FUNCTION FOR WIND.
C
C     PROGRAMMER P. KOCLAS NCAR MMM (1991)
C
C***********************************************************************
      Data PII /1.5707963/
      If (ZETA.LT.0.) Then
        X = (1.-15.*ZETA) ** 0.2
        BUSPSI_M1 = 2. * ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2. * 
     *      ATAN(X) + PII
C        if(buspsi_m1 .gt. 2.0) buspsi_m1=2.0
      Else
        BUSPSI_M1 = -4.7 * ZETA
      End If
      Return
      End








      Function BUSPSI_H2(ZETA)
      Implicit None
      Real ZETA, BUSPSI_H2, Y
C***********************************************************************
C
C     PURPOSE: TO CALCULATE STABILITY CORRECTION FUNCTION IN THE 
C               SURFACE LAYER TEMPERATURE PROFILE.
C
C     ARGUMENTS:
C
C          INPUT:
C                 -ETA: STABILITY PARAMETER (Z/L)
C
C         OUTPUT:
C                 -BUSFIH2 :CORRECTION FUNCTION FOR TEMPERATURE.
C
C     PROGRAMMER P. KOCLAS NCAR MMM (1991)
C
C***********************************************************************
      If (ZETA.LT.0.) Then
        Y = SQRT(1.-9.*ZETA)
        BUSPSI_H2 = ALOG((1.+Y)/2.) * 2.
      Else
        BUSPSI_H2 = -4.7 * ZETA
      End If
      Return
      End








      Subroutine BLACKA(SCAL,LAMBDA,VK,ZT,NBOT,NTOP)
      Implicit None
      Integer NBOT, NTOP, K
      Real LAMBDA, VK, KZ
      Real ZT(NBOT:NTOP), SCAL(NBOT:NTOP)
C***********************************************************************
C
C     PURPOSE: TO CALCULATE LENGTH SCALES VIA THE BLACKADAR
C              FORMULATION.
C
C     ARGUMENTS:
C
C      INPUT:
C            -LAMBDA: SCALING CONSTANT
C            -    VK: VON KARMAN CONSTANT
C            -    ZT:  VERTICAL GRID FOR TURBULENCE VARIABLES.
C            -  NBOT: INDEX OF CLOSEST POINT TO THE SURFACE
C            -  NTOP: INDEX OF CLOSEST POINT TO THE MODEL TOP
C     OUTPUT:
C            -  SCAL: LENGTH SCALE.
C
C     PROGRAMMER P. KOCLAS NCAR MMM (1991)
C
C***********************************************************************
      Do 10 K = NBOT, NTOP - 1
        KZ = VK * ZT(K)
        SCAL(K) = KZ / (1.+KZ/LAMBDA)
   10 Continue
      Return
      End









      Subroutine GRADRI(RI,DUDZ,DVDZ,BUOYZ,NBOT,NTOP)
      Implicit None
      Integer NBOT, NTOP, K
      Real SHEAR, BACKGR
      Real DUDZ(NTOP), DVDZ(NTOP), BUOYZ(NTOP), RI(NTOP)
C***********************************************************************
C
C     PURPOSE: TO CALCULATE BULK RICHARSON NUMBERS.
C
C     ARGUMENTS:
C
C     INPUT:
C            -  DUDZ: VERTICAL GRADIENT OF 'X' COMPONENT OF WIND.     (1/S)
C            -  DVDZ: VERTICAL GRADIENT OF 'Y' COMPONENT OF WIND.     (1/S)
C            - BUOYZ: BUOYANCY.                                     (1/S*S)
C            -  NBOT: INDEX OF CLOSEST POINT TO THE SURFACE
C            -  NTOP: INDEX OF CLOSEST POINT TO THE MODEL TOP
C     OUTPUT:
C            -    RI: GRADIENT RICHARDSON NUMBER.
C
C     PROGRAMMER P. KOCLAS NCAR MMM (1991)
C
C      -NOTE- : A BACKGROUND SHEAR IS ADDED TO AVOID DIVISION BY ZERO.
C
C
C***********************************************************************
      Data BACKGR /1.E-9/
      Do 10 K = NBOT, NTOP - 1
        SHEAR = (DUDZ(K)**2+DVDZ(K)**2) + BACKGR
        RI(K) = BUOYZ(K) / SHEAR
   10 Continue
      Return
      End


      Subroutine ENTRFLX(FLUX,PROFILE,Z,FLX_TYPE,NBOT,NTOP,GUESS)
      Implicit None
      Integer NBOT, NTOP, FLX_TYPE
      Integer SEARCH_TOP, SEARCH_BOT, FOUND
      Real PROFILE(NBOT:NTOP), Z(NBOT:NTOP)
      Real FLUX, GUESS

      Integer I

      Real MAXMIN
      
C     FIND THE GRID POINT NEAREST THE PBL HEIGHT
      FOUND = -1
      I = NBOT
      Do While (FOUND.LT.0)
        If (Z(I).GT.GUESS.OR.I.EQ.NTOP) Then
          FOUND = I
        End If
        I = I + 1
      End Do
       
      
C     DEFINE THE SEARCH RANGE FOR THE MAX/MIN FLUX 
C     (7 IS AN ARBITRARY NUMBER)    
      SEARCH_TOP = FOUND + 7
      SEARCH_BOT = FOUND - 7
      
C     DO A CHECK ON THE BOUNDS
      If (SEARCH_BOT.LT.NBOT) SEARCH_BOT = NBOT
      If (SEARCH_TOP.GT.NTOP) SEARCH_TOP = NTOP
      
C     LOOK FOR TOP DOWN TYPE SCALAR MINIMUM     
      If (FLX_TYPE.EQ.0) Then
        MAXMIN = 1.0E10
        Do 10 I = SEARCH_BOT, SEARCH_TOP
          If (PROFILE(I).LE.MAXMIN) MAXMIN = PROFILE(I)
   10   Continue

      Else
C     LOOK FOR BOTTOM UP TYPE SCALAR MAXIMUM
        MAXMIN = -1.0E10
        Do 20 I = SEARCH_BOT, SEARCH_TOP
          If (PROFILE(I).GE.MAXMIN) MAXMIN = PROFILE(I)
   20   Continue

      End If


C     SET ENTRAINMENT FLUX TO MAXIMUM/MINIMUM IN PROFILE
      FLUX = MAXMIN

      Return
      End
      

C***********************************************************************
C     CALCULATION OF ustar
C**********************************************************************
      Subroutine CALC_USTAR(U,V,U_OLD,V_OLD,DEL_T,UG,VG,Z,NTOP,F_COR,
     *    USTAR)
      Implicit None
      Integer NTOP
      Real U(NTOP), V(NTOP), UG(NTOP), VG(NTOP), Z(NTOP), USTAR, F_COR
      Real U_OLD(NTOP), V_OLD(NTOP), DEL_T
      Integer I
      Real SUM_UW, SUM_VW
      Real UG_SFC, VG_SFC
      Real DU_DT1, DV_DT1, DU_DT2, DV_DT2

C     EXTRAPOLATE ug AND vg TO SURFACE
      UG_SFC = UG(1) - Z(1) * (UG(2)-UG(1)) / (Z(2)-Z(1))
      VG_SFC = VG(1) - Z(1) * (VG(2)-VG(1)) / (Z(2)-Z(1))

C     ADD LOWEST LAYER TO SUMS 
      DU_DT2 = (U(1)-U_OLD(1)) / DEL_T
      DV_DT2 = (V(1)-V_OLD(1)) / DEL_T

      SUM_UW = (F_COR*VG_SFC+DU_DT2-F_COR*(V(1)-VG(1))) * Z(1) / 2.0
      SUM_VW = (-F_COR*UG_SFC+DV_DT2-F_COR*(UG(1)-U(1))) * Z(1) / 2.0

C     COMPLETE INTEGRAL TO TOP OF MODEL
      Do 10 I = 2, NTOP
        DU_DT1 = (U(I-1)-U_OLD(I-1)) / DEL_T
        DV_DT1 = (V(I-1)-V_OLD(I-1)) / DEL_T
        DU_DT2 = (U(I)-U_OLD(I)) / DEL_T
        DV_DT2 = (V(I)-V_OLD(I)) / DEL_T
	
        SUM_UW = SUM_UW + (DU_DT1-F_COR*(V(I-1)-VG(I-1))+DU_DT2-F_COR*(
     *      V(I)-VG(I))) * (Z(I)-Z(I-1)) / 2.0

        SUM_VW = SUM_VW + (DV_DT1-F_COR*(UG(I-1)-U(I-1))+DV_DT2-F_COR*(
     *      UG(I)-U(I))) * (Z(I)-Z(I-1)) / 2.0

   10 Continue

C     CALCULATE ustar
      USTAR = (SUM_UW**2+SUM_VW**2) ** 0.25

      Return
      End


C***********************************************************************
C     SECANT ROOT FINDING ROUTINE FROM NUMERICAL RECIPES THIS
C     ROUTINE IS USED TO FIND A NON-NEUTRAL VALUE OF ustar
C***********************************************************************
      Function ZBRENT2(FUNC,X1,X2,TOL,USTAR,WIND,TREF,GRAV,WTSFC,Z1,Z0M,
     *    VK,ZI)
c     Parameter (ITMAX = 100,EPS = 3.E-8)
      Parameter (ITMAX = 100,EPS = 3.E-6)
      Real WIND, TREF, GRAV, WTSFC, Z1, Z0M, VK, ZI

      Real FUNC
      External FUNC

      A = X1
      B = X2
      FA = FUNC(A,USTAR,WIND,TREF,GRAV,WTSFC,Z1,Z0M,VK,ZI)
      FB = FUNC(B,USTAR,WIND,TREF,GRAV,WTSFC,Z1,Z0M,VK,ZI)
      If (FB*FA.GT.0.) Pause 'Root must be bracketed for ZBRENT2.'
      FC = FB
      Do 10 ITER = 1, ITMAX
        If (FB*FC.GT.0.) Then
          C = A
          FC = FA
          D = B - A
          E = D
        End If
        If (ABS(FC).LT.ABS(FB)) Then
          A = B
          B = C
          C = A
          FA = FB
          FB = FC
          FC = FA
        End If
        TOL1 = 2. * EPS * ABS(B) + 0.5 * TOL
        XM = .5 * (C-B)
        If (ABS(XM).LE.TOL1.OR.FB.EQ.0.) Then
          ZBRENT2 = B
          Return
        End If
        If (ABS(E).GE.TOL1.AND.ABS(FA).GT.ABS(FB)) Then
          S = FB / FA
          If (A.EQ.C) Then
            P = 2. * XM * S
            Q = 1. - S
          Else
            Q = FA / FC
            R = FB / FC
            P = S * (2.*XM*Q*(Q-R)-(B-A)*(R-1.))
            Q = (Q-1.) * (R-1.) * (S-1.)
          End If
          If (P.GT.0.) Q = -Q
          P = ABS(P)
          If (2.*P.LT.MIN(3.*XM*Q-ABS(TOL1*Q),ABS(E*Q))) Then
            E = D
            D = P / Q
          Else
            D = XM
            E = D
          End If
        Else
          D = XM
          E = D
        End If
        A = B
        FA = FB
        If (ABS(D).GT.TOL1) Then
          B = B + D
        Else
          B = B + SIGN(TOL1,XM)
        End If
        FB = FUNC(B,USTAR,WIND,TREF,GRAV,WTSFC,Z1,Z0M,VK)
   10 Continue
      Pause 'ZBRENT2 exceeding maximum iterations.'
      ZBRENT2 = B
c
      End
      Function zbrent(FUNC,X1,X2,TOL,USTAR,WIND,TREF,GRAV,WTSFC,Z1,Z0M,
     *    VK,ZI)
      Parameter (ITMAX = 100,EPS = 1.E-4)
      Real FUNC
      External FUNC
      a = x1
      b = x2
      ta = func(a,ustar,wind,tref,grav,wtsfc,z1,z0m,vk,zi)
      tb = func(b,ustar,wind,tref,grav,wtsfc,z1,z0m,vk,zi)
      if(ta .gt. 0) then
         fa = ta
         fb = tb
         a  = x1
         b  = x2
      else
         fa = tb
         fb = ta
         a  = x2
         b  = x1
      endif
c     write(6,7009) fa,fb
c7009 format(' 7009 zbrent, fa = ',e15.6,' fb = ',e15.6)
      if (fb*fa.gt.0.)  then
         write(6,8000) a,b,fa,fb
 8000    format(' 8000, Trouble ZBRENT first guess roots',/,
     +          '        must be of opposite sign',/,
     +          ' a = ',e15.6,5x,' fa = ',e15.6,/,
     +          ' b = ',e15.6,5x,' fb = ',e15.6)
         stop
      endif
      iter = 0
   10 continue
        iter = iter + 1
        c = 0.5*(a+b)
        fc = func(c,ustar,wind,tref,grav,wtsfc,z1,z0m,vk)
        diff = abs(c-b)
c          write(6,9000) a,b,c,fa,fb,fc,diff
c9000      format(' a = ',e15.6,5x,' b = ',e15.6,5x,' c = ',e15.6,/, 
c    +            ' fa = ',e15.6,5x,' fb = ',e15.6,5x,' fc = ',e15.6,/,
c    +            ' diff = ',e15.6)
        if(diff .lt. eps*abs(b)) go to 99
        if(iter .gt. itmax) then
           write(6,6000) a,b,c,diff,fa,fb,fc,itmax
 6000      format('6000 Trouble in ZBRENT',/,
     +            ' a = ',e15.6,/,
     +            ' b = ',e15.6,/,
     +            ' c = ',e15.6,/,
     +            ' fa = ',e15.6,/,
     +            ' fb = ',e15.6,/,
     +            ' fc = ',e15.6,/,
     +            ' itmax = ',i5)
           stop
        endif
        if(fc .gt. 0.0) then
           a = c
        else
           b = c
        endif
        go to 10
c
   99 continue
      zbrent = c
c
      end


C***********************************************************************
C     BRACKETING ROUTINE FOR BRACKETING THE ROOT(S) OF ustar
C***********************************************************************
      Subroutine ZBRAK(FX,X1,X2,N,XB1,XB2,NB,USTAR,WIND,TREF,GRAV,WTSFC,
     *    Z1,Z0M,VK,ZI)
      Real WIND, TREF, GRAV, WTSFC, Z1, Z0M, VK
      Integer N, NB
      Real X1, X2, XB1, XB2, FX
      External FX
      Integer I, NBB
      Real DX, FC, FP, X

      NBB = 0
      X = X1
      DX = (X2-X1) / N
      FP = FX(X,USTAR,WIND,TREF,GRAV,WTSFC,Z1,Z0M,VK,ZI)
      Do 10 I = 1, N
        X = X + DX
        FC = FX(X,USTAR,WIND,TREF,GRAV,WTSFC,Z1,Z0M,VK,ZI)
        If (FC*FP.LT.0.) Then
          NBB = NBB + 1
          XB1 = X - DX
          XB2 = X
        End If
        FP = FC
        If (NBB.EQ.NB) Return
   10 Continue
      Return
      End

C***********************************************************************
C     MINIMIZATION FUNCTION CONTAINING ROOT OF ustar USING SURFACE
C     LAYER SIMILARITY FUNCTIONS
C***********************************************************************
      Function USTAR_RES(A,USTAR,WIND,TREF,GRAV,WTSFC,Z1,Z0M,VK,ZI)
      Implicit None

      Real A, USTAR, WIND, TREF, GRAV, WTSFC, Z1, Z0M, VK
      Real USTAR_RES, BUSPSI_M1
      Real FPSI_M
      Real ZETA, USTR, ZI

C     ITERATION VARIABLE IS PASSED IN a ARGUMENT AND IS COPIED
c     FOR USE IN RESIDUAL CALCULATION
      USTR = A

C     CALCULATE zeta (z/L)
      ZETA = -Z1 * VK * GRAV * WTSFC / (TREF*USTR**3)
      ZETA = AMAX1(-2.0,ZETA)

      USTAR_RES = WIND - USTR / VK * (ALOG((Z1+Z0M)/Z0M)-FPSI_M(ZETA))
C     1            (alog((z1+z0m)/z0m) - buspsi_m1(zeta))

      Return
      End


C***********************************************************************
C     MINIMIZATION FUNCTION CONTAINING ROOT OF ustar USING PBL 
C     SIMILARITY FUNCTIONS
C***********************************************************************
      Function USTARPBL_RES(A,USTAR,WIND,TREF,GRAV,WTSFC,Z1,Z0M,VK,ZI)
      Implicit None

      Real A, USTAR, WIND, TREF, GRAV, WTSFC, Z1, Z0M, VK, ZI
      Real USTARPBL_RES
      Real ZETA, USTR
      Real L, XI, ETA, SMAL, zi_l
      Real F, G, RIB

C     ITERATION VARIABLE IS PASSED IN a ARGUMENT AND IS COPIED
c     FOR USE IN RESIDUAL CALCULATION
      USTR = A

      SMAL = 0.00001
      if(abs(wtsfc) .lt. smal) then
         L = -9999.0
         zi_l = 0.0
      else
         L = -USTR ** 3 * TREF / (VK*GRAV*WTSFC)
         zi_l = zi/l
      endif
      zi_l  = zi/l
      write(6,5001) wtsfc, l, zi_l
 5001 format(' 5001 wtsfc = ',e15.6,' L = ',e15.6,
     +       ' zi/l = ',e15.6)
      XI = (1.0-0.025*15.0*ZI_L) ** 0.25
      ETA = (1.0-0.025*9.0*ZI_L) ** 0.25
      write(6,9001) xi,eta
 9001 format(' 9001 xi = ',e15.6,' eta = ',e15.6)

C     CALCULATE F AND G
      F = 1.0 / VK * (ALOG(0.025*ZI/Z0M)-ALOG((1.0+XI**2)/2.0)-2.0*
     *    ALOG((1.0+XI)/2.0)+2.0*ATAN(XI)-1.5707963) + 8.4 * (1.0-50.0*
     *    ZI_L) ** (-0.16)

      G = 0.74 / VK * (ALOG(0.025*ZI/Z0M)-2.0*ALOG((1.0+ETA**2)/2.0)) +
     *    7.3 * (1.0-5.8*ZI_L) ** (-0.47)
      write(6,9002) f,g   
 9002 format(' 9002 f = ',e15.6,' g = ',e15.6)

C     CALCULATE THE BULK RICHARDSON NUMBER (Rib)
      RIB = ZI / (VK*L) * G / (F**2)
c     RIB = ZI_L/VK * G / (F**2)

C     CALCULATE zeta
      ZETA = LOG(ABS(RIB)) - 3.5

      USTARPBL_RES = USTR / WIND - 1.0 / ((ALOG(0.025*ZI/Z0M)/VK+8.4)-
     *    25.0*EXP(0.26*ZETA-0.03*ZETA**2))
c
      write(6,7001) USTARPBL_RES
7001  format(' 7001 USTARPBL_RES = ',e15.6)

      Return
      End


C***********************************************************************
C     MINIMIZATION FUNCTION CONTAINING ROOT OF z1
C***********************************************************************
      Function Z1_RES(A,USTAR,WIND,TREF,GRAV,WTSFC,Z1,Z0M,VK)
      Implicit None

      Real A, USTAR, WIND, TREF, GRAV, WTSFC, Z1, Z0M, VK
      Real Z1_RES, BUSPSI_M1
      Real FPSI_M
      Real ZETA, Z

C     ITERATION VARIABLE IS PASSED IN a ARGUMENT AND IS COPIED
C     FOR USE IN RESIDUAL CALCULATION
      Z = A

C     CALCULATE zeta (z/L)
      ZETA = -Z * VK * GRAV * WTSFC / (TREF*USTAR**3)
      ZETA = AMAX1(-2.0,ZETA)

      Z1_RES = WIND - USTAR / VK * (ALOG((Z+Z0M)/Z0M)-FPSI_M(ZETA))
C     1            (alog((z+z0m)/z0m) - buspsi_m1(zeta))
C      z1_res = wind - ustar / vk * (alog((z+z0m)/z0m))


      Return
      End


C***********************************************************************
C     SMOOTHING ROUTINE FOR STABILITY BOUNDS MODEL
C**********************************************************************
      Subroutine MIXSTAB(NK,Z,U,V,T,Q,B,C,ZI)
      Implicit None
      
      Integer NK
      Real Z(NK), U(NK), V(NK), T(NK), Q(NK), B(NK), C(NK), ZI
      Real UTOTAL, VTOTAL, TTOTAL, QTOTAL, BTOTAL, CTOTAL
      Integer I, IZI

      Do 10 I = 1, NK
        If (Z(I).GT.ZI) Then
          IZI = I - 1
          Go To 20
        End If
   10 Continue

C     CALCULATE AVERAGES
   20 UTOTAL = 0.0
      VTOTAL = 0.0
      TTOTAL = 0.0
      QTOTAL = 0.0
      BTOTAL = 0.0
      CTOTAL = 0.0
      Do 30 I = 1, IZI
        UTOTAL = UTOTAL + U(I) / IZI
        VTOTAL = VTOTAL + V(I) / IZI
        TTOTAL = TTOTAL + T(I) / IZI
        QTOTAL = QTOTAL + Q(I) / IZI
        BTOTAL = BTOTAL + B(I) / IZI
        CTOTAL = CTOTAL + C(I) / IZI
   30 Continue

C     ASSIGN AVERAGES TO PROFILE BELOW zi
      Do 40 I = 1, IZI
        U(I) = UTOTAL
        V(I) = VTOTAL
        T(I) = TTOTAL
        Q(I) = QTOTAL
        B(I) = BTOTAL
        C(I) = CTOTAL
   40 Continue

      Return
      End
      
C***********************************************************************
C     MEAN PROFILE OUTPUT ROUTINE
C**********************************************************************
      Subroutine OUTMEAN(NK,Z,U,V,T,B,C,TIME,OUTUNIT)
      Implicit None
      
      Integer NK, TIME, OUTUNIT, I
      Real Z(NK), U(NK), V(NK), T(NK), B(NK), C(NK)
      
C     OUTPUT TIME AND VARIABLE LABLES 
      Write (OUTUNIT,5000) TIME

      Do 10, I = 1, NK
        Write (OUTUNIT,5100) Z(I), U(I), V(I), T(I), B(I), C(I)
   10 Continue
	       
      Return
 5000 Format (I10,' zm/u/v/theta/B/C')
 5100 Format (6E12.5)
      End
