!copyright (C) 2001  MSC-RPN COMM  %%%RPNPHY%%%
***S/P  RIGRAD1
*

      SUBROUTINE RIGRAD1(RI,GAMA,GAMAQ,TBL,DUDZ2,T,TVE,Q,QE, 1
     %                   SIGMA, SE, WW,  N, M, NK )
#include "impnone.cdk"
      INTEGER N,M,NK
      REAL RI(N,NK),GAMA(N,NK),GAMAQ(N,NK),TBL(N)
      REAL DUDZ2(N,NK),T(M,NK),TVE(N,NK),Q(M,NK),QE(N,NK)
      REAL SIGMA(N,NK),SE(N,NK)
      REAL WW(N)
      REAL TE,VIRCOR,DZ,TVK,TVKP,FAC,DLNP
      INTEGER j,k
*
*Author
*          J. Mailhot RPN(Feb 1990)
*
*Revision
* 001      J.Mailhot RPN(Feb 1990)shallow convection (GELEYN)
* 002      G.Pellerin(August90)Adaptation to thermo functions
* 003      Y. Delage (Nov 1990)Options of shallow convection
* 004      Y. Delage  (Nov 1990)- Removal of BETA and PRI
* 005      N. Brunet  (May91)
*                New version of thermodynamic functions
*                and file of constants
* 006      C. Girard (Nov92) New parameterization of the
*          shallow convection
* 007      G. Pellerin(May95) Extract level of BLH
* 008      C. Girard (Nov95) Added calculations of GAMAQ
* 009      A. Plante (June 2003) - IBM conversion
*             - calls to vslog routine (from massvp4 library)
*
*Object
*          to calculate the gradient Richardson number
*
*Arguments
*
*          - Outputs -
* RI       gradient Richardson number
* GAMA     equilibrium gradient term for temperature
* GAMAQ    equilibrium gradient term for moisture
* TBL      Level corresponding to top of Unstable boundary layer
*
*          - Inputs -
* DUDZ2    vertical shear of the wind squared
* T        temperature
* TVE      virtual temperature at 'E' levels
* Q        specific humidity
* QE       specific humidity at 'E' levels
* SIGMA    sigma level values
* SE       intermediate sigma level values 
* WW       work field
* N        horizontal dimension
* M        1st dimension of T and Q
* NK       vertical dimension
*
*
**
*
#include "consphy.cdk"
#include "dintern.cdk"
*
#include "fintern.cdk"
*
      DO j = 1, N
         WW(j) = FOTVT( T(j,1), Q(j,1) )
      END DO
*
*
      DO k = 1, NK - 1
         DO j = 1, N
            RI(j,k)=SIGMA(j,k+1)/SIGMA(j,k)
         ENDDO
      ENDDO
      CALL VSLOG(RI,RI,N*(NK-1))

      DO k = 1, NK - 1
         DO j = 1, N
*
*           TEMPERATURES VIRTUELLES
*
            TVK = WW(j)
            WW(j) = FOTVT( T(j,k+1), Q(j,k+1) )
            TVKP = WW(j)
            TE = FOTTV( TVE(j,k), QE(j,k) )
            VIRCOR = TVE(j,k) / TE
*
            DLNP =  RI(j,k)
            DZ = RGASD * TVE(j,k) * DLNP / GRAV
*
*           RI
*
            FAC = GRAV / ( TVE(j,k) * DUDZ2(j,k) )
            RI(j,k) = FAC * ( ( TVK - TVKP ) / DZ + GRAV/CPD )
*
*           GAMA
*
            GAMA(j,k) = - GRAV / ( CPD * VIRCOR )
*
            GAMAQ(j,k) = 0.
*
*           TOP OF THE unstable BOUNDARY LAYER
*
            if( RI(j,k).gt.0. ) TBL(j) = k
*
         END DO
      END DO
*
      DO j = 1, N
*
         RI(j,NK) = RI(j,NK-1)
*
      END DO
*
      RETURN
      END