!copyright (C) 2001  MSC-RPN COMM  %%%RPNPHY%%%
*** S/P CLSGS2
*
#include "phy_macros_f.h"

      SUBROUTINE CLSGS2 (THL, TVE, QW, QC, FRAC, FNN, C1,  1,3
     1                  ZN, ZE, S, PS, A, B, C, N, NK)
*
#include "impnone.cdk"
*
*
      INTEGER N, NK
      REAL THL(N,NK), TVE(N,NK), QW(N,NK), QC(N,NK)
      REAL FRAC(N,NK), FNN(N,NK)
      REAL C1(N,NK), ZN(N,NK), ZE(N,NK), S(N,NK)
      REAL A(N,NK), B(N,NK), C(N,NK)
      REAL PS(N)
*
*Author
*          J. Mailhot (Jun 2002)
*
*Revision
* 001      J. Mailhot (Feb 2003) Clipping at upper levels
* 002      S. Belair  (Apr 2003) Minimum values of 50 m for ZE and ZN
*                                in calculation of sigmase.
* 002      A-M. Leduc (Jun 2003) Pass ps to blweight ---> blweight2
*
*Object
*          Calculate the boundary layer sub-grid-scale cloud properties 
*
*Arguments
*
*          - Input -
* THL      cloud water potential temperature
* TVE      virtual temperature on 'E' levels
* QW       total water content
*
*          - Output -
* QC       cloud water content
* FRAC     cloud fraction
* FNN      flux enhancement factor (fn) times cloud fraction (N)
*
*          - Input -
* C1       constant C1 in second-order moment closure
* ZN       length scale for turbulent mixing (on 'E' levels)
* ZE       length scale for turbulent dissipationa (on 'E' levels)
* S        sigma levels
* PS       surface pressure (in Pa)
* A        thermodynamic coefficient
* B        thermodynamic coefficient
* C        thermodynamic coefficient
* N        horizontal dimension
* NK       vertical dimension
*
*
*Notes
*          Implicit (i.e. subgrid-scale) cloudiness scheme for unified
*             description of stratiform and shallow, nonprecipitating
*             cumulus convection appropriate for a low-order turbulence
*             model based on Bechtold et al.:
*            - Bechtold and Siebesma 1998, JAS 55, 888-895
*            - Cuijpers and Bechtold 1995, JAS 52, 2486-2490
*            - Bechtold et al. 1995, JAS 52, 455-463
*            - Bechtold et al. 1992, JAS 49, 1723-1744
*
*
*IMPLICITS
*
#include "consphy.cdk"
*
**
*
      INTEGER J, K, ITOTAL
*
* 
      REAL EPSILON
      REAL QCMIN, QCMAX
*
*
*
***********************************************************
*     AUTOMATIC ARRAYS
**********************************************************
*
      AUTOMATIC ( DZ       , REAL    , (N,NK)  )
      AUTOMATIC ( DQWDZ    , REAL    , (N,NK)  )
      AUTOMATIC ( DTHLDZ   , REAL    , (N,NK)  )
      AUTOMATIC ( SIGMAS   , REAL    , (N,NK)  )
      AUTOMATIC ( SIGMASE  , REAL    , (N,NK)  )
      AUTOMATIC ( Q1       , REAL    , (N,NK)  )
      AUTOMATIC ( WEIGHT   , REAL    , (N,NK)  )
*
***********************************************************
*
*
*MODULES
*
      EXTERNAL DVRTDF, BLWEIGHT
*
*------------------------------------------------------------------------
*
      EPSILON = 1.0E-10
      QCMIN   = 1.0E-6
      QCMAX   = 1.0E-3
*
*
*
*       1.     Vertical derivative of THL and QW
*       ----------------------------------------
*
      DO K=1,NK-1
      DO J=1,N
        DZ(J,K) = -RGASD*TVE(J,K)*ALOG( S(J,K+1)/S(J,K) ) / GRAV
      END DO
      END DO
*
      DO J=1,N
        DZ(J,NK) = 0.0
      END DO
*
      CALL DVRTDF ( DTHLDZ, THL, DZ, N, N, N, NK)
      CALL DVRTDF ( DQWDZ, QW, DZ, N, N, N, NK)
*
*
*       2.     Standard deviation of s and normalized saturation deficit Q1
*       -------------------------------------------------------------------
*
      DO K=1,NK-1
      DO J=1,N
*                                              sigmas (cf. BCMT 1995 eq. 10)
*                                        (computation on 'E' levels stored in SIGMASE)
        SIGMASE(J,K) = 
     1      SQRT( C1(J,K)*MAX(ZN(J,K),50.)*MAX(ZE(J,K),50.) ) *  
     1           ABS( 0.5*(A(J,K)+A(J,K+1))*DQWDZ(J,K)
     1             - 0.5*(B(J,K)+B(J,K+1))*DTHLDZ(J,K) )
      END DO
      END DO
*
      DO K=2,NK-1
      DO J=1,N
*                                              (back to full levels)
        SIGMAS(J,K) = 0.5*( SIGMASE(J,K) + SIGMASE(J,K-1) )
*                                              normalized saturation deficit
        Q1(J,K) = C(J,K) / ( SIGMAS(J,K) + EPSILON )
        Q1(J,K) = MAX ( -6. , MIN ( 4. , Q1(J,K) ) )
      END DO
      END DO
*
      DO J=1,N
        SIGMAS(J,1) = 0.0  
        SIGMAS(J,NK) = 0.0  
        Q1(J,1) = 0.0
        Q1(J,NK) = 0.0
      END DO
*
*
*       3.     Cloud properties
*       -----------------------
*                                              cloud fraction, cloud water content
*                                              and flux enhancement factor
*                                              (cf. BS 1998 Appendix B)
      DO K=2,NK-1
      DO J=1,N
*
        IF( Q1(J,K) .GT. -1.2 ) THEN
          FRAC(J,K) = MAX ( 0. , MIN ( 1. ,
     1                      0.5 + 0.36*ATAN(1.55*Q1(J,K)) ) )
        ELSEIF( Q1(J,K) .GE. -6.0 ) THEN
          FRAC(J,K) = EXP ( Q1(J,K)-1.0 )
        ELSE
          FRAC(J,K) = 0.0
        ENDIF
*
        IF( Q1(J,K) .GE. 0.0 ) THEN
          QC(J,K) = EXP( -1.0 ) + 0.66*Q1(J,K) + 0.086*Q1(J,K)**2
        ELSEIF( Q1(J,K) .GE. -6.0 ) THEN
          QC(J,K) = EXP( 1.2*Q1(J,K)-1.0 )
        ELSE
          QC(J,K) = 0.0
        ENDIF
*
        QC(J,K) = MIN ( QC(J,K)*( SIGMAS(J,K) + EPSILON )
     1                  , QCMAX )
*
        FNN(J,K) = 1.0
        IF( Q1(J,K).LT.1.0 .AND. Q1(J,K).GE.-1.68 ) THEN
          FNN(J,K) = EXP( -0.3*(Q1(J,K)-1.0) )
        ELSEIF( Q1(J,K).LT.-1.68 .AND. Q1(J,K).GE.-2.5 ) THEN
          FNN(J,K) = EXP( -2.9*(Q1(J,K)+1.4) )
        ELSEIF( Q1(J,K).LT.-2.5 ) THEN
          FNN(J,K) = 23.9 + EXP( -1.6*(Q1(J,K)+2.5) )
        ENDIF
*                                              flux enhancement factor * cloud fraction
*                                              (parameterization formulation)
        FNN(J,K) = FNN(J,K)*FRAC(J,K)
        IF( Q1(J,K).LE.-2.39 .AND. Q1(J,K).GE.-4.0 ) THEN
          FNN(J,K) = 0.60
        ELSEIF( Q1(J,K).LT.-4.0 .AND. Q1(J,K).GE.-6.0 ) THEN
          FNN(J,K) = 0.30*( Q1(J,K)+6.0 )
        ELSEIF( Q1(J,K).LT.-6.0 ) THEN
          FNN(J,K) = 0.0
        ENDIF
*
*
      END DO
      END DO
*
      DO J=1,N
        FRAC(J,1) = 0.
        FRAC(J,NK) = 0.
        FNN(J,1) = 0.
        FNN(J,NK) = 0.
        QC(J,1) = 0.
        QC(J,NK) = 0.
      END DO
*
*
      CALL BLWEIGHT2 ( WEIGHT, S, PS, N, NK)
*
      DO K=1,NK
      DO J=1,N
        FRAC(J,K) = FRAC(J,K)*WEIGHT(J,K)
        FNN(J,K) = FNN(J,K)*WEIGHT(J,K)
        QC(J,K) = QC(J,K)*WEIGHT(J,K)
      END DO
      END DO
*
*
      RETURN
      END