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

      SUBROUTINE CONDS(TE,QE,SRR,SSR,FN, 1
     *               TP1,QP1,PSP1,KBL,
     +               SIGMA, TAU, N, NI, NK, DBGCOND, SATUCO)
#include "impnone.cdk"
*
      LOGICAL LO,DBGCOND
      LOGICAL SATUCO
      INTEGER IERR
      REAL TAU
*
*
      INTEGER N,NI,NK
      REAL TEMP1, TEMP2, TEMP3
C
      REAL TP1(N,NK),QP1(N,NK),
     *          PSP1(N),KBL(N),SIGMA(NI,NK)
      REAL TE(NI,NK),QE(NI,NK),SRR(NI),SSR(NI)
      REAL FN(NI,NK)
      INTEGER NIP,NIKP,NKP1,NKM1,JL,JK,IS
*
*Author
*          J.Mailhot 11/03/85.  Adapted from E.C.M.W.F.
*
*Revision
* 001      J. Mailhot (Mar 1987) base of stable condensation
* 002      G.Pellerin(Nov87)Adaptation to revised code
*                     (Mar88)Standard documentations
* 003      J. Mailhot (Mar 1988) threshold of evaporation
* 004      J P Toviessi(May1990)Conversion in CFT77
* 005      G.Pellerin(August90)Standardization of thermo functions
* 006      N. Brunet  (May91)
*                 New version of thermodynamic functions
*                 and file of constants
* 007      B. Bilodeau  (August 1991)- Adaptation to UNIX
*
* 008      J. Mailhot (Dec.1992) - Newton Method Bug Correction
*            (ref. Revision 005)
* 009      C. Girard (Nov.1992) - Small clean-up, more correction
*          to the thickness of the 1st layer and "implicit"
*          evaporation of the precipitation
* 010      A. Methot (Aug 93) - L/Cp added in calculation of evap.
* 011      G. Lemay (Oct 93) - Dynamic memory allocation with stkmemw
* 012      R. Benoit (Aug 93) - Local Sigma
* 013      B. Bilodeau (June 94) - New physics interface
* 014      B. Bilodeau (Jan  01) - Automatic arrays
* 015      M. Lepine  (March 2003) -  CVMG... Replacements
* 016      L. Spacek (June 2003) - IBM conversion
*                                - boucle 331 defectuosite
*
*Object
*          to calculate the T and Q tendencies due to large scale
*          precipitation
*
*Arguments
*
*          - Output -
* TE       temperature tendency due to stratiform processes
* QE       specific humidity tendency due to stratiform processes
*
*          - Output -
* SRR      rate of liquid precipitation
* SSR      rate of solid precipitation
* FN       cloud fraction
*
*          - Input -
* TP1      temperature
* QP1      specific humidity
* PSP1     surface pressure
* KBL      index of first level in boundary layer
* SIGMA    sigma levels
* TAU      FACTDT * timestep (see common bloc OPTIONS)
* N        dimension of some arrays
* NI       1st horizontal dimension
* NK       vertical dimension
* DBGCOND  .TRUE. to have debugging for condensation
*          .FALSE. to have no debugging
* SATUCO   .TRUE. to have water/ice phase for saturation
*          .FALSE. to have water phase only for saturation
*
*Notes
*          During the process, the variables T and Q at (T+DT) are
*          adjusted instantly. There is no storage of water or snow
*          in the cloud. In super-saturated layers, Q is restored
*          to Q* and the difference between Q and Q* is added to
*          the fluxes of rain and snow from the top to the bottom
*          of the layer. The evaporation, condensation or melting
*          can affect the divergence of the precipitation fluxes.
*          Reference in ECMWF Research Manual (Vol.3)
*          Physical Parameterization Chapter 5)
*
**
*
************************************************************************
*     AUTOMATIC ARRAYS
************************************************************************
*
      AUTOMATIC ( IQCD   , INTEGER , (NI   ) )
*
      AUTOMATIC ( LO1    , LOGICAL , (NI   ) )
*
      AUTOMATIC ( ZPP1   , REAL    , (NI,NK) )
      AUTOMATIC ( ZDSG   , REAL    , (NI,NK) )
      AUTOMATIC ( ZDPP1  , REAL    , (NI,NK) )
      AUTOMATIC ( ZTP1   , REAL    , (NI,NK) )
      AUTOMATIC ( ZQP1   , REAL    , (NI,NK) )
      AUTOMATIC ( ZQSATE , REAL    , (NI,NK) )
      AUTOMATIC ( ZTC    , REAL    , (NI,NK) )
      AUTOMATIC ( ZQC    , REAL    , (NI,NK) )
      AUTOMATIC ( ZRFL   , REAL    , (NI   ) )
      AUTOMATIC ( ZSFL   , REAL    , (NI   ) )
      AUTOMATIC ( ZRFLN  , REAL    , (NI   ) )
      AUTOMATIC ( ZSFLN  , REAL    , (NI   ) )
      AUTOMATIC ( ZFLN   , REAL    , (NI   ) )
*
************************************************************************
*
      REAL ZSTPRO,ZDIP,ZEVAP,ZMELT,ZSQFLN,ZNIMP
      REAL ZEPFLM,ZEPFLS,ZCONS1,ZCONS2
      REAL ZTMST,ZLDCPE,ZQCD
      REAL ZQSATC,ZCOR,ZRITO,ZLVDCP,ZLSDCP,ZRIT
*
C
C*    PHYSICAL CONSTANTS.
C     -------- ----------
#include "comphy.cdk"
*
      REAL CHLS
#include "consphy.cdk"
#include "dintern.cdk"
#include "fintern.cdk"
*
C     -------------
C
C     *ZEVAP*IS A CONSTANT FOR THE EVAPORATION OF
C     TOTAL PRECIPITATION, *ZMELT* IS THE CONSTANT OF THE FORMULA FOR
C     THE RATE OF CHANGE OF THE LIQUID WATER/ICE COMPOSITION OF THESE
C     PRECIPITATIONS.
C
      NKP1=NK+1
      NKM1=NK-1
      ZEVAP=CEVAP
      ZMELT=CMELT
C
C*    SECURITY PARAMETERS.
C     --------------------
C
C         *ZEPFLM* IS A MINIMUM FLUX TO AVOID DIVIDING BY ZERO IN THE IC
C     PROPORTION CALCULATIONS.
C
      ZEPFLM=1.E-24
      ZEPFLS=SQRT(ZEPFLM)
C
C*    COMPUTATIONAL CONSTANTS.
C     ------------- ----------
C
      ZTMST= TAU
      CHLS = CHLC + CHLF
C
      ZCONS1=CPV - CPD
      ZCONS2 = 1./(ZTMST*GRAV)
C
C
C     ------------------------------------------------------------------
C
C
C*         1.     PRELIMINARY COMPUTATIONS.
C                 ----------- -------------
C
  200 CONTINUE
C
C*         1.1     INITIAL VALUES FOR ACCUMULATION
C
  210 CONTINUE
C
      DO 211 JL=1,NI
      ZRFL(JL)=0.
      ZSFL(JL)=0.
      ZFLN(JL)=0.
  211 CONTINUE
C
C
C
C     ------------------------------------------------------------------
C
C*         2.     CLOUD VARIABLES, RAIN/SNOW FLUXES.
C                 ----- ---------- --------- -------
C
  300 CONTINUE
C
C*         2.1     T+1 T,Q VARIABLES AND SATURATION MIXING RATIO.
C
  310 CONTINUE
      DO 3150 JL=1,NI
      ZDSG(JL,1)=0.5*(SIGMA(JL,2)-SIGMA(JL,1))
      DO 315 JK=2,NKM1
      ZDSG(JL,JK)=0.5*(SIGMA(JL,JK+1)-SIGMA(JL,JK-1))
  315 CONTINUE
      ZDSG(JL,NK)=0.5*(1.-SIGMA(JL,NKM1))+0.5*(1.-SIGMA(JL,NK))
 3150 continue
C
      DO 311 JK=1,NK
      DO 311 JL=1,NI
      ZPP1(JL,JK)=SIGMA(JL,JK)*PSP1(JL)
      ZTP1(JL,JK)=TP1(JL,JK)
      ZQP1(JL,JK)=QP1(JL,JK)
      ZDPP1(JL,JK)=ZDSG(JL,JK)*PSP1(JL)
  311 CONTINUE
      IF(SATUCO)THEN
      DO 312 JK=1,NK
      DO 312 JL=1,NI
      TEMP1 = ZTP1(JL,JK)
      TEMP2 = ZPP1(JL,JK)
  312 ZQSATE(JL,JK)=FOQST(TEMP1,TEMP2)
      ELSE
      DO 314 JK=1,NK
      DO 314 JL=1,NI
      TEMP1= ZTP1(JL,JK)
      TEMP2= ZPP1(JL,JK)
  314 ZQSATE(JL,JK)=FOQSA(TEMP1,TEMP2)
      ENDIF
C
C
C*         2.2     CALCULATE TC AND QC IN SUPERSATURATED LAYERS. THE
C*                 CONDENSATION CALCULATIONS ARE DONE WITH TWO ITERATION
C
  320 CONTINUE
C
C***
      DO 329 JK=1,NK
C***
      DO 322 JL=1,NI
      ZTC(JL,JK)=ZTP1(JL,JK)
      ZQC(JL,JK)=ZQP1(JL,JK)
  322 CONTINUE
      IF(SATUCO)THEN
      DO 323 JL=1,NI
      ZQSATC=ZQSATE(JL,JK)
!      ZLDCPE = CVMGT(CHLC,CHLS,ZTC(JL,JK)-TCDK .GT. 0.)
!     *      /(CPD+ZCONS1*ZQC(JL,JK))
      if (ZTC(JL,JK)-TCDK .GT. 0.) then
         ZLDCPE = CHLC /(CPD+ZCONS1*ZQC(JL,JK))
      else
         ZLDCPE = CHLS /(CPD+ZCONS1*ZQC(JL,JK))
      endif
      TEMP1 = ZTC(JL,JK)
      ZCOR=ZLDCPE*FODQS(ZQSATC,TEMP1)
      ZQCD=AMAX1(0.,(ZQC(JL,JK)-ZQSATC)/(1.+ZCOR))
      LO=ZQCD.EQ.0.
!      IQCD(JL) = CVMGT(0.,1.,LO)
      if (ZQCD.EQ.0.) then
         IQCD(JL) = 0.
      else
         IQCD(JL) = 1.
      endif
      ZQC(JL,JK)=ZQC(JL,JK)-ZQCD
      ZTC(JL,JK)=ZTC(JL,JK)+ZQCD*ZLDCPE
  323 CONTINUE
      ELSE
      DO 327 JL=1,NI
      ZQSATC=ZQSATE(JL,JK)
!      ZLDCPE= CVMGT(CHLC, CHLS,  ZTC(JL,JK)-TCDK .GT. 0.)
!     *      /(CPD+ZCONS1*ZQC(JL,JK))
      if (ZTC(JL,JK)-TCDK .GT. 0.) then
         ZLDCPE = CHLC /(CPD+ZCONS1*ZQC(JL,JK))
      else
         ZLDCPE = CHLS /(CPD+ZCONS1*ZQC(JL,JK))
      endif
      TEMP1 = ZTC(JL,JK)
      ZCOR=ZLDCPE*FODQA(ZQSATC,TEMP1)
      ZQCD=AMAX1(0.,(ZQC(JL,JK)-ZQSATC)/(1.+ZCOR))
      LO=ZQCD.EQ.0.
!      IQCD(JL) = CVMGT(0.,1.,LO)
      if (ZQCD.EQ.0.) then
        IQCD(JL) = 0.
      else
        IQCD(JL) = 1.
      endif
      ZQC(JL,JK)=ZQC(JL,JK)-ZQCD
      ZTC(JL,JK)=ZTC(JL,JK)+ZQCD*ZLDCPE
  327 CONTINUE
      ENDIF
      IS=0
      DO 324 JL=1,NI
      IS=IS+IQCD(JL)
  324 CONTINUE
      IF (IS.NE.0) THEN
      IF(SATUCO)THEN
      DO 325 JL=1,NI
      TEMP1 = ZTC (JL,JK)
      TEMP2 = ZPP1(JL,JK)
      ZQSATC= FOQST(TEMP1,TEMP2)
!      ZLDCPE = CVMGT(CHLC,CHLS,ZTC(JL,JK)-TCDK .GT. 0.)
!     *      /(CPD+ZCONS1*ZQC(JL,JK))
      if (ZTC(JL,JK)-TCDK .GT. 0.) then
         ZLDCPE = CHLC /(CPD+ZCONS1*ZQC(JL,JK))
      else
         ZLDCPE = CHLS /(CPD+ZCONS1*ZQC(JL,JK))
      endif
      ZCOR=ZLDCPE*FODQS(ZQSATC,TEMP1)
      ZQCD=(ZQC(JL,JK)-ZQSATC)/(1.+ZCOR)
      LO1(JL)=IQCD(JL).NE.0
!      ZQCD = CVMGT(ZQCD,0.,LO1(JL))
      if (IQCD(JL) .eq. 0) ZQCD = 0.
      ZQC(JL,JK)=ZQC(JL,JK)-ZQCD
      ZTC(JL,JK)=ZTC(JL,JK)+ZQCD*ZLDCPE
  325 CONTINUE
      ELSE
      DO 328 JL=1,NI
      TEMP1 = ZTC (JL,JK)
      TEMP2 = ZPP1(JL,JK)
      ZQSATC= FOQSA(TEMP1,TEMP2)
!      ZLDCPE= CVMGT(CHLC, CHLS,  ZTC(JL,JK)-TCDK .GT. 0.)
!     *      /(CPD+ZCONS1*ZQC(JL,JK))
      if (ZTC(JL,JK)-TCDK .GT. 0.) then
         ZLDCPE = CHLC /(CPD+ZCONS1*ZQC(JL,JK))
      else
         ZLDCPE = CHLS /(CPD+ZCONS1*ZQC(JL,JK))
      endif
      ZCOR=ZLDCPE*FODQA(ZQSATC,TEMP1)
      ZQCD=(ZQC(JL,JK)-ZQSATC)/(1.+ZCOR)
      LO1(JL)=IQCD(JL).NE.0
!      ZQCD = CVMGT(ZQCD,0.,LO1(JL))
      if (IQCD(JL) .eq. 0) ZQCD = 0.
      ZQC(JL,JK)=ZQC(JL,JK)-ZQCD
      ZTC(JL,JK)=ZTC(JL,JK)+ZQCD*ZLDCPE
  328 CONTINUE
      ENDIF
      ENDIF
      DO 326 JL=1,NI
      LO1(JL)=ZQP1(JL,JK).LE.ZQSATE(JL,JK)
!      ZTC(JL,JK) = CVMGT(ZTP1(JL,JK),ZTC(JL,JK),LO1(JL))
      if (LO1(JL)) ZTC(JL,JK) = ZTP1(JL,JK)
!      ZQC(JL,JK) = CVMGT(ZQP1(JL,JK),ZQC(JL,JK),LO1(JL))
      if (LO1(JL)) ZQC(JL,JK) = ZQP1(JL,JK)
  326 CONTINUE
C
C
  329 CONTINUE
C
C
C
C***
      DO 645 JK=1,NK
C
C
C*         3.3     CALCULATE RAIN/SNOW FLUX IN SUPERSATURATED LAYERS.
C
  330 CONTINUE
C
C***
      DO 331 JL=1,NI
      LO = ZTC(JL,JK) .GT. TGL
      ZSTPRO    =AMAX1((ZQP1(JL,JK)-ZQC(JL,JK)),0.)
      TEMP1=ZSTPRO*ZDPP1(JL,JK)*ZCONS2
!      ZRFLN(JL)=ZRFL(JL)+CVMGT(ZSTPRO    *ZDPP1(JL,JK)*ZCONS2,0.,LO)
!      ZSFLN(JL)=ZSFL(JL)+CVMGT(0.,ZSTPRO    *ZDPP1(JL,JK)*ZCONS2,LO)
      ZRFLN(JL)=ZRFL(JL)
      ZSFLN(JL)=ZSFL(JL)
      IF (LO) THEN
         ZRFLN(JL) = ZRFLN(JL) + TEMP1
      ELSE
         ZSFLN(JL) = ZSFLN(JL) + TEMP1
      ENDIF
  331 CONTINUE
C
C
C***
      IF (JK.GT.1) THEN
C
C     ------------------------------------------------------------------
C
C*         3.     EVAPORATION OF PRECIPITATIONS.
C                 ----------- -- ---------------
C
  400 CONTINUE
C
      DO 521 JL=1,NI
C***
C
C*         3.2     EVAPORATION OF PRECIPITATIONS.
C

  420 CONTINUE
      ZSQFLN = SQRT( ZRFLN(JL)+ZSFLN(JL) )
      TEMP1 = ZQSATE(JL,JK)
      TEMP2 = ZTP1(JL,JK)
!      ZLDCPE= CVMGT(CHLC, CHLS, TEMP2-TCDK .GT. 0.)
!     *      /(CPD+ZCONS1*ZQC(JL,JK))
      if (TEMP2-TCDK .GT. 0.) then
         ZLDCPE = CHLC /(CPD+ZCONS1*ZQC(JL,JK))
      else
         ZLDCPE = CHLS /(CPD+ZCONS1*ZQC(JL,JK))
      endif

C
      ZNIMP = 1. + 2.*(1.+ ZLDCPE*FODQS(TEMP1,TEMP2))
     *              *ZEVAP*ZSQFLN/ZCONS2
C
      ZFLN(JL) = (AMAX1(0.,ZSQFLN-ZEVAP*ZDPP1(JL,JK)
     *           *AMAX1(0.,ZQSATE(JL,JK)-ZQP1(JL,JK))/ZNIMP ))**2
C
C
C     ------------------------------------------------------------------
C
C*         4.     MELTING/FREEZING, OUTGOING RAIN/SNOW FLUXES.
C                 ----------------- -------- --------- ------
C
  500 CONTINUE
C
C
C
C*         5.1     MELTING/FREEZING OF PRECIPITATIONS.
C*         5.2     OUTGOING FLUXES AT THE BOTTOM OF THE LAYER.
C
  520 CONTINUE
      ZDIP    =ZDPP1(JL,JK)/ZPP1(JL,JK)**2
      ZRITO=(ZSFLN(JL)/AMAX1(ZSFLN(JL)+ZRFLN(JL),ZEPFLM))
      ZRIT=ZRITO-ZMELT*ZDIP    *(ZTC(JL,JK)-TGL)/AMAX1(ZEPFLS,0.5
     *     *(SQRT(ZFLN(JL))+SQRT(ZRFL(JL)+ZSFL(JL))))
      ZRIT=AMIN1(1.,AMAX1(0.,ZRIT))
      ZSFLN(JL)=ZRIT*ZFLN(JL)
      ZRFLN(JL)=ZFLN(JL)-ZSFLN(JL)
  521 CONTINUE
C
C
C
C     ------------------------------------------------------------------
C
C*         5.     TENDENCIES DUE TO CONDENSATION, SURFACE FLUXES.
C                 ---------- --- -- ------------- ------- -------
C
  600 CONTINUE
C
C
C*         5.2     UPDATE FLAG IN ACTIVE LAYERS.
C
  620 CONTINUE
C***
      ENDIF
C***
      DO 621 JL=1,NI
      ZLVDCP=CHLC/(CPD+ZCONS1*ZQC(JL,JK))
      ZLSDCP=CHLS/(CPD+ZCONS1*ZQC(JL,JK))
      QE(JL,JK)= -((ZRFLN(JL)-ZRFL(JL))+(ZSFLN(JL)
     +           - ZSFL(JL)))*(GRAV/ZDPP1(JL,JK))
      TE(JL,JK)=(ZLVDCP*(ZRFLN(JL)-ZRFL(JL))+ZLSDCP*(ZSFLN(JL)
     +          - ZSFL(JL)))*(GRAV/ZDPP1(JL,JK))
  621 CONTINUE
C
C*         5.3     DO ZONAL MEAN AND BOX DIAGNOSTICS.
C
  630 CONTINUE
C
C*         5.4     SWAP OF FLUXES, END OF VERTICAL LOOP AND STABLE
C*                 RAIN AND SNOW RATES.
C
  640 CONTINUE
      DO 641 JL=1,NI
      ZRFL(JL)=ZRFLN(JL)
      ZSFL(JL)=ZSFLN(JL)
  641 CONTINUE
C***
  645 CONTINUE
C***
      DO 647 JL=1,NI
      SRR(JL)=ZRFL(JL)
      SSR(JL)=ZSFL(JL)
  647 CONTINUE
C
C
*
      RETURN
      END