function minval (q,im)
      implicit none
#include "CPP_EEOPTIONS.h"
      integer  im, i
      real q(im), minval
      minval = 1.e15
      do i=1,im
      if( q(i).lt.minval ) minval = q(i)
      enddo
      return
      end
      FUNCTION ERRF (ARG)
C***********************************************************************
C  FUNCTION ERRF
C  PURPOSE
C     COMPUTES ERROR FUNCTION OF ARGUMENT
C  USAGE
C     CALLED BY TRBFLX
C  DESCRIPTION OF PARAMETERS
C     ARG   -  INPUTED ARGUMENT
C  REMARKS:
C        USED TO COMPUTE FRACTIONAL CLOUD COVER AND LIQUID WATER CONTENT
C          FROM TURBULENCE STATISTICS
C **********************************************************************
      implicit none
      real arg,errf
 
      real aa1,aa2,aa3,aa4,aa5,pp,x2,x3,x4,x5
      PARAMETER ( AA1  =   0.254829592 )
      PARAMETER ( AA2  =  -0.284496736 )
      PARAMETER ( AA3  =   1.421413741 )
      PARAMETER ( AA4  =  -1.453152027 )
      PARAMETER ( AA5  =   1.061405429 )
      PARAMETER ( PP   =   0.3275911   )
      PARAMETER ( X2   =   AA2 / AA1   )
      PARAMETER ( X3   =   AA3 / AA2   )
      PARAMETER ( X4   =   AA5 / AA3   )
      PARAMETER ( X5   =   AA5 / AA4   )
 
      real aarg,tt
      
      ERRF = 1.
      AARG=ABS(ARG)
 
      IF ( AARG .LT. 4.0 )  THEN
       TT = 1./(1.+PP*AARG)
       ERRF = 1. -
     1  (AA1*TT*(1.+X2*TT*(1.+X3*TT*(1.+X4*TT*(1.+X5*TT)))))
     2   * EXP(-AARG*AARG)
      ENDIF
 
      IF ( ARG .LT. 0.0 )  ERRF = -ERRF
 
      RETURN
      END

      SUBROUTINE STRIP(A,B,IA,IB,L,K)                                           
      implicit none
#include "CPP_EEOPTIONS.h"
      integer ia,ib,L,K
      real A(IA,L), B(IB,L)                                                
      
      INTEGER OFFSET,Len,i,j
                                                                                
      OFFSET = IB*(K-1)                                                         
      LEN    = MIN(IB,IA-OFFSET)                                                
      OFFSET = OFFSET+1                                                         
                                                                                
      IF(LEN.EQ.IB) THEN                                                        
        DO 100 J=1,L                                                            
        DO 100 I=1,LEN                                                          
        B(I,J) = A(I+OFFSET-1,J)                                                
100     CONTINUE                                                                
      ELSE                                                                      
        DO 200 J=1,L                                                            
        DO 300 I=1,LEN                                                          
        B(I,J) = A(I+OFFSET-1,J)                                                
300     CONTINUE                                                                
        DO 400 I=1,IB-LEN                                                       
        B(LEN+I,J) = A(LEN+OFFSET-1,J)                                            
400     CONTINUE                                                                
200     CONTINUE                                                                
      ENDIF                                                                     
                                                                                
      RETURN                                                                    
      END                                                                       
      SUBROUTINE PASTE(B,A,IB,IA,L,K)                                           
      implicit none
#include "CPP_EEOPTIONS.h"
      integer ia,ib,L,K
      real A(IA,L), B(IB,L)                                                

      INTEGER OFFSET,LEN,i,j
                                                                                
      OFFSET = IB*(K-1)                                                         
      LEN    = MIN(IB,IA-OFFSET)                                                
      OFFSET = OFFSET+1                                                         
                                                                                
      DO 100 J=1,L                                                              
      DO 100 I=1,LEN                                                            
      A(I+OFFSET-1,J) = B(I,J)                                                  
100   CONTINUE                                                                  
                                                                                
      RETURN                                                                    
      END                                                                       
      SUBROUTINE PSTBMP(B,A,IB,IA,L,K)                                          
      implicit none
#include "CPP_EEOPTIONS.h"
      integer ia,ib,L,K
      real A(IA,L), B(IB,L)                                                

      INTEGER OFFSET,LEN,i,j
                                                                                
      OFFSET = IB*(K-1)                                                         
      LEN    = MIN(IB,IA-OFFSET)                                                
      OFFSET = OFFSET+1                                                         
                                                                                
      DO 100 J=1,L                                                              
      DO 100 I=1,LEN                                                            
      A(I+OFFSET-1,J) = A(I+OFFSET-1,J) + B(I,J)                                
100   CONTINUE                                                                  
C                                                                               
      RETURN                                                                    
      END                                                                       
      SUBROUTINE STRINT(A,B,IA,IB,L,K)
      implicit none
      integer ia,ib,L,K
      INTEGER A(IA,L), B(IB,L)

      INTEGER OFFSET,LEN,i,j

      OFFSET = IB*(K-1)
      LEN = MIN(IB,IA-OFFSET)
      OFFSET = OFFSET+1

      IF(LEN.EQ.IB) THEN
        DO 100 J=1,L
        DO 100 I=1,LEN
        B(I,J) = A(I+OFFSET-1,J)
100     CONTINUE
      ELSE
        DO 200 J=1,L
        DO 300 I=1,LEN
        B(I,J) = A(I+OFFSET-1,J)
300     CONTINUE
        DO 400 I=1,IB-LEN
        B(LEN+I,J) = A(LEN+OFFSET-1,J)
400     CONTINUE
200     CONTINUE
      ENDIF

      RETURN
      END
      SUBROUTINE QSAT (TT,P,Q,DQDT,LDQDT)
C***********************************************************************
C
C  PURPOSE:
C  ========
C    Compute Saturation Specific Humidity
C
C  INPUT:
C  ======
C    TT ......... Temperature (Kelvin)
C    P .......... Pressure (mb)
C    LDQDT ...... Logical Flag to compute QSAT Derivative
C
C  OUTPUT:
C  =======
C    Q .......... Saturation Specific Humidity
C    DQDT ....... Saturation Specific Humidity Derivative wrt Temperature
C
C
C***********************************************************************

      IMPLICIT NONE
#include "CPP_EEOPTIONS.h"
      real TT, P, Q, DQDT
      LOGICAL LDQDT

      real AIRMW, H2OMW
      
      PARAMETER ( AIRMW  = 28.97      )                                         
      PARAMETER ( H2OMW  = 18.01      )                                         

      real ESFAC, ERFAC
      PARAMETER ( ESFAC = H2OMW/AIRMW       )
      PARAMETER ( ERFAC = (1.0-ESFAC)/ESFAC )

      real aw0, aw1, aw2, aw3, aw4, aw5, aw6
      real bw0, bw1, bw2, bw3, bw4, bw5, bw6
      real ai0, ai1, ai2, ai3, ai4, ai5, ai6
      real bi0, bi1, bi2, bi3, bi4, bi5, bi6

      real d0, d1, d2, d3, d4, d5, d6
      real e0, e1, e2, e3, e4, e5, e6
      real f0, f1, f2, f3, f4, f5, f6
      real g0, g1, g2, g3, g4, g5, g6

c ********************************************************
c ***  Polynomial Coefficients WRT Water (Lowe, 1977) ****
c ***              (Valid +50 C to -50 C)             ****
c ********************************************************

      parameter ( aw0 =  6.107799961e+00 * esfac )
      parameter ( aw1 =  4.436518521e-01 * esfac )
      parameter ( aw2 =  1.428945805e-02 * esfac )
      parameter ( aw3 =  2.650648471e-04 * esfac )
      parameter ( aw4 =  3.031240396e-06 * esfac )
      parameter ( aw5 =  2.034080948e-08 * esfac )
      parameter ( aw6 =  6.136820929e-11 * esfac )

      parameter ( bw0 = +4.438099984e-01 * esfac )
      parameter ( bw1 = +2.857002636e-02 * esfac )
      parameter ( bw2 = +7.938054040e-04 * esfac )
      parameter ( bw3 = +1.215215065e-05 * esfac )
      parameter ( bw4 = +1.036561403e-07 * esfac )
      parameter ( bw5 = +3.532421810e-10 * esfac )
      parameter ( bw6 = -7.090244804e-13 * esfac )


c ********************************************************
c ***   Polynomial Coefficients WRT Ice  (Lowe, 1977) ****
c ***              (Valid  +0 C to -50 C)             ****
c ********************************************************

      parameter ( ai0 = +6.109177956e+00 * esfac )
      parameter ( ai1 = +5.034698970e-01 * esfac )
      parameter ( ai2 = +1.886013408e-02 * esfac )
      parameter ( ai3 = +4.176223716e-04 * esfac )
      parameter ( ai4 = +5.824720280e-06 * esfac )
      parameter ( ai5 = +4.838803174e-08 * esfac )
      parameter ( ai6 = +1.838826904e-10 * esfac )

      parameter ( bi0 = +5.030305237e-01 * esfac )
      parameter ( bi1 = +3.773255020e-02 * esfac )
      parameter ( bi2 = +1.267995369e-03 * esfac )
      parameter ( bi3 = +2.477563108e-05 * esfac )
      parameter ( bi4 = +3.005693132e-07 * esfac )
      parameter ( bi5 = +2.158542548e-09 * esfac )
      parameter ( bi6 = +7.131097725e-12 * esfac )


c ********************************************************
c ***         Polynomial Coefficients WRT Ice         ****
c ***   Starr and Cox (1985) (Valid -40 C to -70 C)   ****
c ********************************************************


      parameter ( d0 = 0.535098336e+01 * esfac )
      parameter ( d1 = 0.401390832e+00 * esfac )
      parameter ( d2 = 0.129690326e-01 * esfac )
      parameter ( d3 = 0.230325039e-03 * esfac )
      parameter ( d4 = 0.236279781e-05 * esfac )
      parameter ( d5 = 0.132243858e-07 * esfac )
      parameter ( d6 = 0.314296723e-10 * esfac )

      parameter ( e0 = 0.469290530e+00 * esfac )
      parameter ( e1 = 0.333092511e-01 * esfac )
      parameter ( e2 = 0.102164528e-02 * esfac )
      parameter ( e3 = 0.172979242e-04 * esfac )
      parameter ( e4 = 0.170017544e-06 * esfac )
      parameter ( e5 = 0.916466531e-09 * esfac )
      parameter ( e6 = 0.210844486e-11 * esfac )


c ********************************************************
c ***         Polynomial Coefficients WRT Ice         ****
c ***   Starr and Cox (1985) (Valid -65 C to -95 C)   ****
c ********************************************************

      parameter ( f0 = 0.298152339e+01 * esfac )
      parameter ( f1 = 0.191372282e+00 * esfac )
      parameter ( f2 = 0.517609116e-02 * esfac )
      parameter ( f3 = 0.754129933e-04 * esfac )
      parameter ( f4 = 0.623439266e-06 * esfac )
      parameter ( f5 = 0.276961083e-08 * esfac )
      parameter ( f6 = 0.516000335e-11 * esfac )

      parameter ( g0 = 0.312654072e+00 * esfac )
      parameter ( g1 = 0.195789002e-01 * esfac )
      parameter ( g2 = 0.517837908e-03 * esfac )
      parameter ( g3 = 0.739410547e-05 * esfac )
      parameter ( g4 = 0.600331350e-07 * esfac )
      parameter ( g5 = 0.262430726e-09 * esfac )
      parameter ( g6 = 0.481960676e-12 * esfac )

      real TMAX, TICE
      PARAMETER ( TMAX=323.15, TICE=273.16)
      
      real T, D, W, QX, DQX
      T = MIN(TT,TMAX) - TICE
      DQX = 0.
      QX  = 0.

c Fitting for temperatures above 0 degrees centigrade
c ---------------------------------------------------
      if(t.gt.0.) then
       qx = aw0+T*(aw1+T*(aw2+T*(aw3+T*(aw4+T*(aw5+T*aw6)))))
      if (ldqdt)  then
      dqx = bw0+T*(bw1+T*(bw2+T*(bw3+T*(bw4+T*(bw5+T*bw6)))))
      endif
      endif

c Fitting for temperatures between 0 and -40
c ------------------------------------------
      if( t.le.0. .and. t.gt.-40.0 ) then
        w = (40.0 + t)/40.0
       qx =     w *(aw0+T*(aw1+T*(aw2+T*(aw3+T*(aw4+T*(aw5+T*aw6))))))
     .    + (1.-w)*(ai0+T*(ai1+T*(ai2+T*(ai3+T*(ai4+T*(ai5+T*ai6))))))
      if (ldqdt)  then
      dqx =     w *(bw0+T*(bw1+T*(bw2+T*(bw3+T*(bw4+T*(bw5+T*bw6))))))
     .    + (1.-w)*(bi0+T*(bi1+T*(bi2+T*(bi3+T*(bi4+T*(bi5+T*bi6))))))
      endif
      endif

c Fitting for temperatures between -40 and -70
c --------------------------------------------
      if( t.le.-40.0 .and. t.ge.-70.0 ) then
       qx = d0+T*(d1+T*(d2+T*(d3+T*(d4+T*(d5+T*d6)))))
      if (ldqdt) then
      dqx = e0+T*(e1+T*(e2+T*(e3+T*(e4+T*(e5+T*e6)))))
      endif
      endif

c Fitting for temperatures less than -70
c --------------------------------------
      if(t.lt.-70.0) then
       qx = f0+t*(f1+t*(f2+t*(f3+t*(f4+t*(f5+t*f6)))))
      if (ldqdt) then
      dqx = g0+t*(g1+t*(g2+t*(g3+t*(g4+t*(g5+t*g6)))))
      endif
      endif

c Compute Saturation Specific Humidity
c ------------------------------------
      D = (P-ERFAC*QX)
      IF(D.LT.0.) THEN
       Q = 1.0
       IF (LDQDT)  DQDT = 0.
      ELSE
       D = 1.0 / D
       Q = MIN(QX * D,1.0)
       IF (LDQDT)  DQDT = (1.0 + ERFAC*Q) * D * DQX
      ENDIF
      RETURN
      END
      subroutine vqsat (tt,p,q,dqdt,ldqdt,n)
      implicit none
#include "CPP_EEOPTIONS.h"
      integer  i,n
      logical  ldqdt
      real tt(n), p(n), q(n), dqdt(n)
#if CRAY
#if f77
cfpp$ expand (qsat)
#endif
#endif
      do i=1,n
      call qsat ( tt(i),p(i),q(i),dqdt(i),ldqdt )
      enddo
      return
      end

      subroutine stripit(a,b,irun,ia,ib,l,k)
      implicit none
#include "CPP_EEOPTIONS.h"
      integer ia,ib,irun,l,k
      real a(ia,l), b(ib,l)
      integer   i,j,len,offset

      offset = ib*(k-1)
      len    = min(ib,irun-offset)
      offset = offset+1

      if(len.eq.ib) then
        do 100 j=1,l
        do 100 i=1,len
        b(i,j) = a(i+offset-1,j)
100     continue
      else
        do 200 j=1,l
        do 300 i=1,len
        b(i,j) = a(i+offset-1,j)
300     continue
        do 400 i=1,ib-len
        b(len+i,j) = a(len+offset-1,j)
400     continue
200     continue
      endif
      return
      end

      subroutine stripitint(a,b,irun,ia,ib,l,k)
      implicit none
#include "CPP_EEOPTIONS.h"
      integer ia,ib,irun,l,k,a(ia,l),b(ib,l)
      integer i,j,len,offset

      offset = ib*(k-1)
      len = min(ib,irun-offset)
      offset = offset+1

      if(len.eq.ib) then
        do 100 j=1,l
        do 100 i=1,len
        b(i,j) = a(i+offset-1,j)
100     continue
      else
        do 200 j=1,l
        do 300 i=1,len
        b(i,j) = a(i+offset-1,j)
300     continue
        do 400 i=1,ib-len
        b(len+i,j) = a(len+offset-1,j)
400     continue
200     continue
      endif
      return
      end

      subroutine pastit(b,a,ib,ia,irun,L,k)
      implicit none
#include "CPP_EEOPTIONS.h"
      integer ib,ia,L,k,irun,len,offset
      integer i,j
      real a(ia,l), b(ib,l)

      offset = ib*(k-1)
      len    = min(ib,irun-offset)
      offset = offset+1

      do 100 j=1,L
      do 100 i=1,len
      a(i+offset-1,j) = b(i,j)
100   continue
      return
      end

      subroutine pstbitint(b,a,ib,ia,irun,l,k)
      implicit none
#include "CPP_EEOPTIONS.h"
      integer ib,ia,L,k,irun,len,offset
      real a(ia,l)
      integer b(ib,l)
      integer i,j

      offset = ib*(k-1)
      len    = min(ib,irun-offset)
      offset = offset+1
               
      do 100 j=1,L
      do 100 i=1,len
      a(i+offset-1,j) = a(i+offset-1,j) + float(b(i,j))
100   continue
      return  
      end        


      subroutine pstbmpit(b,a,ib,ia,irun,l,k)
      implicit none
#include "CPP_EEOPTIONS.h"
      integer ib,ia,L,k,irun,len,offset
      real a(ia,l), b(ib,l)
      integer i,j

      offset = ib*(k-1)
      len    = min(ib,irun-offset)
      offset = offset+1
               
      do 100 j=1,L
      do 100 i=1,len
      a(i+offset-1,j) = a(i+offset-1,j) + b(i,j)
100   continue
      return  
      end        

      subroutine strip2tile(a,index,b,irun,ia,ib,levs,npeice)
c-----------------------------------------------------------------------
c subroutine strip2tile - extract one processors worth of grid points
c                         from a grid space array to a stripped tile
c                         space array
c 
c  input:
c      a      - array to be stripped FROM [ia,levs]
c      index  - array of horizontal indeces of grid points to convert to
c               tile space
c      irun   - number of points in array a that need to be stripped
c      ia     - inner of dimension of source array
c      ib     - inner dimension of target array AND the number of points
c               in the target array to be filled
c      levs   - number of vertical levels AND outer dimension of arrays
c      npeice - the current strip number to be filled
c output:
c      b      - array to be filled, ie, one processors field [ib,levs]
c-----------------------------------------------------------------------
      implicit none
#include "CPP_EEOPTIONS.h"
      integer ia,ib,irun,levs,npeice
      real a(ia,levs), b(ib,levs)
      integer index(irun)
      integer i,k,len,offset

      offset = ib*(npeice-1)
      len    = min(ib,irun-offset)
      offset = offset+1

      if(len.eq.ib) then
        do 100 k=1,levs
        do 100 i=1,len
        b(i,k) = a(index(i+offset-1),k)
100     continue
      else
        do 200 k=1,levs
        do 300 i=1,len
        b(i,k) = a(index(i+offset-1),k)
300     continue
        do 400 i=1,ib-len
        b(len+i,k) = a(index(len+offset-1),k)
400     continue
200     continue
      endif
      return
      end

      subroutine paste2grd_old(b,index,chfr,ib,numpts,a,ia,levs,npeice)
c-----------------------------------------------------------------------
c subroutine paste2grd - paste one processors worth of grid points
c                        from a stripped tile array to a grid
c                        space array
c 
c  input:
c      b      - array to be pasted back into grid space [ib,levs]
c      index  - array of horizontal indeces of grid points to convert to
c               tile space[numpts]
c      chfr   - fractional area covered by the tile [ib]
c      ib     - inner dimension of source array AND number of points in 
c               array a that need to be pasted
c      numpts - total number of points which were stripped
c      ia     - inner of dimension of target array
c      levs   - number of vertical levels AND outer dimension of arrays
c      npeice - the current strip number to be filled
c output:
c      a      - grid space array to be filled [ia,levs]
c
c IMPORTANT NOTE:
c
c This routine will result in roundoff differences if called from
c within a parallel region.
c-----------------------------------------------------------------------

      implicit none
#include "CPP_EEOPTIONS.h"
      integer ia,ib,levs,numpts,npeice
      integer index(numpts)
      real a(ia,levs), b(ib,levs), chfr(ib)

      integer i,L,offset,len

      offset = ib*(npeice-1)
      len    = min(ib,numpts-offset)
      offset = offset+1

      do L = 1,levs
      do i=1,len
       a(index(i+offset-1),L) = a(index(i+offset-1),L) + b(i,L)*chfr(i)
      enddo
      enddo
      return
      end
      subroutine paste2grd (b,index,chfr,ib,numpts,a,ia,levs,npeice,
     .                                                            check)
c-----------------------------------------------------------------------
c subroutine paste2grd - paste one processors worth of grid points
c                        from a stripped tile array to a grid
c                        space array
c 
c  input:
c      b      - array to be pasted back into grid space [ib,levs]
c      index  - array of horizontal indeces of grid points to convert to
c               tile space[numpts]
c      chfr   - fractional area covered by the tile [ib]
c      ib     - inner dimension of source array AND number of points in 
c               array a that need to be pasted
c      numpts - total number of points which were stripped
c      ia     - inner of dimension of target array
c      levs   - number of vertical levels AND outer dimension of arrays
c      npeice - the current strip number to be filled
c      check  - logical to check for undefined values
c output:
c      a      - grid space array to be filled [ia,levs]
c
c IMPORTANT NOTE:
c
c This routine will result in roundoff differences if called from
c within a parallel region.
c-----------------------------------------------------------------------

      implicit none
#include "CPP_EEOPTIONS.h"
      integer ia,ib,levs,numpts,npeice
      integer index(numpts)
      real a(ia,levs), b(ib,levs), chfr(ib)
      logical check

      integer i,L,offset,len
      real    undef,getcon

      offset = ib*(npeice-1)
      len    = min(ib,numpts-offset)
      offset = offset+1

      if( check ) then
      undef = getcon('UNDEF')
       do L= 1,levs
       do i= 1,len
        if( a(index(i+offset-1),L).eq.undef .or. b(i,L).eq.undef ) then
         a(index(i+offset-1),L) = undef
        else
         a(index(i+offset-1),L)=a(index(i+offset-1),L) + b(i,L)*chfr(i)
        endif
       enddo
       enddo
      else
       do L= 1,levs
       do i= 1,len
        a(index(i+offset-1),L)=a(index(i+offset-1),L) + b(i,L)*chfr(i)
       enddo
       enddo
      endif

      return
      end
      SUBROUTINE GRD2MSC(A,IM,JM,IGRD,B,MXCHPS,NCHP)

      implicit none
#include "CPP_EEOPTIONS.h"
      integer im,jm,mxchps,nchp
      integer igrd(mxchps)
      real A(IM,JM), B(MXCHPS)

      integer i

      IF(NCHP.GE.0) THEN
        DO I=1,NCHP
         B(I) = A(IGRD(I),1)
        ENDDO
      ELSE
        PRINT *, 'ERROR IN GRD2MSC'
      ENDIF

      RETURN
      END

      SUBROUTINE MSC2GRD(IGRD,CHFR,B,MXCHPS,NCHP,FRACG,A,IM,JM)

      implicit none
#include "CPP_EEOPTIONS.h"
      real zero,one
      parameter ( one = 1.)
      parameter (zero = 0.)
      integer im,jm,mxchps,nchp
      integer igrd(mxchps)
      real A(IM,JM), B(MXCHPS), CHFR(MXCHPS), FRACG(IM,JM)

      real VT1(IM,JM)
      integer i

      IF(NCHP.GE.0) THEN
        DO I=1,IM*JM
         VT1(I,1) = ZERO
        ENDDO

        DO I=1,NCHP
         VT1(IGRD(I),1) = VT1(IGRD(I),1) + B(I)*CHFR(I)
        ENDDO

        DO I=1,IM*JM
         A(I,1) = A(I,1)*(ONE-FRACG(I,1)) + VT1(I,1)
        ENDDO
      ELSE
        PRINT *, 'ERROR IN MSC2GRD'
      ENDIF

      RETURN
      END

      subroutine chpprm(nymd,nhms,mxchps,nchp,chlt,ityp,alai,
     1       agrn,zoch,z2ch,cdrc,cdsc,sqsc,ufac,rsl1,rsl2,rdcs)

       implicit none
#include "CPP_EEOPTIONS.h"

       integer nymd,nhms,nchp,mxchps,ityp(mxchps)
       real chlt(mxchps)
       real alai(mxchps),agrn(mxchps)
       real zoch(mxchps), z2ch(mxchps), cdrc(mxchps), cdsc(mxchps)
       real sqsc(mxchps), ufac(mxchps), rsl1(mxchps), rsl2(mxchps)
       real rdcs(mxchps)

C*********************************************************************
C********************* SUBROUTINE CHPPRM  ****************************
C**********************  14 JUNE 1991   ******************************
C*********************************************************************

      integer ntyps
      parameter (ntyps=10)

      real pblzet
      parameter (pblzet = 50.)
      integer k1,k2,nymd1,nhms1,nymd2,nhms2,i
      real getcon,vkrm,rootl,vroot,dum1,dum2,alphaf
      real facm,facp
      real scat,d

      real
     &     vgdd(12, ntyps),    vgz0(12, ntyps),
     &     vgrd(12, ntyps),    vgrt(12, ntyps),

     &     vgrf11(ntyps),         vgrf12(ntyps),
     &     vgtr11(ntyps),         vgtr12(ntyps),
     &     vgroca(ntyps),         vgrotd(ntyps),
     &     vgrdrs(ntyps),         vgz2  (ntyps)


      data vgz0  /
     1   2.6530, 2.6530, 2.6530, 2.6530, 2.6530, 2.6530, 2.6530, 
     1   2.6530, 2.6530, 2.6530, 2.6530, 2.6530,  
     2   0.5200, 0.5200, 0.6660, 0.9100, 1.0310, 1.0440, 1.0420,
     2   1.0370, 1.0360, 0.9170, 0.6660, 0.5200, 
     3   1.1120, 1.1030, 1.0880, 1.0820, 1.0760, 1.0680, 1.0730,
     3   1.0790, 1.0820, 1.0880, 1.1030, 1.1120, 
     4   0.0777, 0.0778, 0.0778, 0.0779, 0.0778, 0.0771, 0.0759,
     4   0.0766, 0.0778, 0.0779, 0.0778, 0.0778, 
     5   0.2450, 0.2450, 0.2270, 0.2000, 0.2000, 0.2000, 0.2000,
     5   0.267,  0.292,  0.280,  0.258,  0.2450,  
     6   0.0752, 0.0752, 0.0752, 0.0752, 0.0752, 0.0757, 0.0777,
     6   0.0778, 0.0774, 0.0752, 0.0752, 0.0752, 
     7   0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 
     7   0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 
     8   0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 
     8   0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 
     9   0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 
     9   0.0112, 0.0112, 0.0112, 0.0112, 0.0112,
     1   0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 0.0112, 
     1   0.0112, 0.0112, 0.0112, 0.0112, 0.0112
     &  /


      data vgrd    /
     1   285.87, 285.87, 285.87, 285.87, 285.87, 285.87, 285.87,
     1   285.87, 285.87, 285.87, 285.87, 285.87, 
     2   211.32, 211.32, 218.78, 243.40, 294.87, 345.90, 355.18,
     2   341.84, 307.22, 244.84, 218.78, 211.32, 
     3   565.41, 587.05, 623.46, 638.13, 652.86, 675.04, 660.24, 
     3   645.49, 638.13, 623.46, 587.05, 565.41, 
     4    24.43,  24.63,  24.80,  24.96,  25.72,  27.74,  30.06,
     4    28.86,  25.90,  25.11,  24.80,  24.63, 
     5   103.60, 103.60, 102.35, 100.72, 100.72, 100.72, 100.72,
     5   105.30, 107.94, 106.59, 104.49, 103.60, 
     6    22.86,  22.86,  22.86,  22.86,  22.86,  23.01,  24.36,
     6    24.69,  24.04,  22.86,  22.86,  22.86, 
     7    23.76,  23.76,  23.76,  23.76,  23.76,  23.76,  23.76,
     7    23.76,  23.76,  23.76,  23.76,  23.76, 
     8    23.76,  23.76,  23.76,  23.76,  23.76,  23.76,  23.76,
     8    23.76,  23.76,  23.76,  23.76,  23.76, 
     9    23.76,  23.76,  23.76,  23.76,  23.76,  23.76,  23.76,
     9    23.76,  23.76,  23.76,  23.76,  23.76, 
     1    23.76,  23.76,  23.76,  23.76,  23.76,  23.76,  23.76,
     1    23.76,  23.76,  23.76,  23.76,  23.76
     &  /

      data vgrt  /
     1   19737.8, 19737.8, 19737.8, 19737.8, 19737.8, 19737.8, 19737.8,
     1   19737.8, 19737.8, 19737.8, 19737.8, 19737.8,
     2    5010.0,  5010.0,  5270.0,  6200.0,  8000.0,  9700.0,  9500.0,
     2    8400.0,  6250.0,  5270.0,  5010.0,  5010.0,
     3    9000.0,  9200.0,  9533.3,  9666.7,  9800.0,  9866.7,  9733.3, 
     3    9666.7,  9533.3,  9200.0,  9000.0,  9000.0,
     4    5500.0,  5625.0,  5750.0,  5875.0,  6625.0,  8750.0,  9375.0,
     4    6875.0,  6000.0,  5750.0,  5625.0,  5500.0,
     5    6500.0,  6000.0,  5500.0,  5500.0,  5500.0,  5500.0,  5500.0,
     5    7500.0,  8500.0,  7000.0,  6500.0,  6500.0,
     6   10625.0, 10625.0, 10625.0, 10625.0, 10625.0, 11250.0, 18750.0,
     6   17500.0, 10625.0, 10625.0, 10625.0, 10625.0,
     7       1.0,     1.0,     1.0,     1.0,     1.0,     1.0,     1.0,
     7       1.0,     1.0,     1.0,     1.0,     1.0,
     8       1.0,     1.0,     1.0,     1.0,     1.0,     1.0,     1.0,
     8       1.0,     1.0,     1.0,     1.0,     1.0,
     9       1.0,     1.0,     1.0,     1.0,     1.0,     1.0,     1.0,
     9       1.0,     1.0,     1.0,     1.0,     1.0,
     1       1.0,     1.0,     1.0,     1.0,     1.0,     1.0,     1.0,
     1       1.0,     1.0,     1.0,     1.0,     1.0
     &  /


      data vgdd  /
     1   27.37,   27.37,   27.37,   27.37,   27.37,   27.37,  27.37,
     1   27.37,   27.37,   27.37,   27.37,   27.37, 
     2   13.66,   13.66,   14.62,   15.70,   16.33,   16.62,   16.66,
     2   16.60,   16.41,   15.73,   14.62,   13.66, 
     3   13.76,   13.80,   13.86,   13.88,   13.90,   13.93,   13.91,
     3   13.89,   13.88,   13.86,   13.80,   13.76, 
     4    0.218,   0.227,   0.233,   0.239,   0.260,   0.299,   0.325,
     4    0.313,   0.265,   0.244,   0.233,   0.227,  
     5    2.813,   2.813,   2.662,   2.391,   2.391,   2.391,   2.391,
     5    2.975,   3.138,   3.062,   2.907,   2.813,  
     6    0.10629, 0.10629, 0.10629, 0.10629, 0.10629, 0.12299, 0.21521,
     6    0.22897, 0.19961, 0.10629, 0.10629, 0.10629,  
     7    0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,
     7    0.0001,  0.0001,  0.0001,  0.0001,  0.0001, 
     8    0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,
     8    0.0001,  0.0001,  0.0001,  0.0001,  0.0001, 
     9    0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,
     9    0.0001,  0.0001,  0.0001,  0.0001,  0.0001,
     1    0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,  0.0001,
     1    0.0001,  0.0001,  0.0001,  0.0001,  0.0001 
     &  /


      data vgrf11 /0.10,0.10,0.07,0.105,0.10,0.10,.001,.001,.001,.001/

      data vgrf12 /0.16,0.16,0.16,0.360,0.16,0.16,.001,.001,.001,.001/

      data vgtr11 /0.05,0.05,0.05,0.070,0.05,0.05,.001,.001,.001,.001/

      data vgtr12 /.001,.001,.001, .220,.001,.001,.001,.001,.001,.001/

      data vgroca /
     &  0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6, 0.384E-6,
     &   .1E-6, .1E-6, .1E-6, .1E-6  / 

      data vgrotd /1.00,1.00,0.50,0.50,0.50,0.20,0.10,0.10,0.10,0.10/ 

      data vgrdrs  /
     &   0.75E13, 0.75E13, 0.75E13, 0.40E13, 0.75E13, 0.75E13,
     &   0.10E13, 0.10E13, 0.10E13, 0.10E13  /

      data vgz2 /35.0, 20.0, 17.0, 0.6, 5.0, 0.6, 0.1, 0.1, 0.1, 0.1/

      vkrm = GETCON('VON KARMAN')

      call time_bound  ( nymd,nhms, nymd1,nhms1, nymd2,nhms2, k1,k2 )
      call interp_time ( nymd,nhms, nymd1,nhms1, nymd2,nhms2, facm,facp)

      do i=1,nchp

      zoch(i) = vgz0(k2,ityp(i))*facp + vgz0(k1,ityp(i))*facm
      rdcs(i) = vgrd(k2,ityp(i))*facp + vgrd(k1,ityp(i))*facm

      rootl   = vgrt(k2,ityp(i))*facp + vgrt(k1,ityp(i))*facm

      vroot  = rootl * vgroca(ityp (i))
      dum1   = log (vroot / (1. - vroot))
      dum2   = 1. / (8. * 3.14159 * rootl)
      alphaf = dum2 * (vroot - 3. -2. * dum1)

      rsl1(i) = vgrdrs (ityp (i)) / (rootl * vgrotd (ityp (i)))
      rsl2(i) = alphaf / vgrotd (ityp (i))

      scat =       agrn(i) *(vgtr11(ityp(i)) + vgrf11(ityp(i)))
     &     + (1. - agrn(i))*(vgtr12(ityp(i)) + vgrf12(ityp(i)))
      sqsc(i) = sqrt (1. - scat)

      d       = vgdd(k2,ityp(i))*facp + vgdd(k1,ityp(i))*facm
      ufac(i) = log( (vgz2(ityp(i)) - d) / zoch(i) ) 
     *        / log(  pblzet             / zoch(i) )

      z2ch(i) = vgz2(ityp (i))

      cdsc(i) = pblzet/zoch(i)+1.
      cdrc(i) = vkrm/log(cdsc(i))
      cdrc(i) = cdrc(i)*cdrc(i)
      cdsc(i) = sqrt(cdsc(i))
      cdsc(i) = cdrc(i)*cdsc(i)

      enddo

      return
      end

      subroutine pkappa (ple,pkle,im,jm,lm,pkz)
C***********************************************************************
C  Purpose
C     Calculate Phillips P**Kappa
C
C  Arguments
C     PLE .... edge-level pressure
C     PKLE ... edge-level pressure**kappa
C     IM ..... longitude  dimension
C     JM ..... latitude   dimension
C     LM ..... vertical   dimension
C     PKZ .... mid-level pressure**kappa
C***********************************************************************
      implicit  none

      integer   im,jm,lm
      real      ple (im,jm,lm+1)
      real      pkle(im,jm,lm+1)
      real      pkz (im,jm,lm)

      real akap1,getcon
      integer i,j,L
 
      akap1 = 1.0 + getcon('KAPPA')
                  
      do L = 1,lm
      do j = 1,jm
      do i = 1,im
       pkz(i,j,L) = ( ple (i,j,l+1)*pkle(i,j,l+1)
     .             - ple (i,j,l)*pkle(i,j,l) )
     .           /  ( akap1* (ple (i,j,l+1)-ple (i,j,l)) )
      enddo
      enddo
      enddo

      return
      end