SUBROUTINE largescale2(dtime, paprs, pplay, t, q,
     s                   d_t, d_q, d_ql, rneb)

      use watercommon_h, only : RLVTT

      IMPLICIT none

!==================================================================
!     
!     Purpose
!     -------
!     Calculates large-scale (stratiform) H2O condensation.
!     
!     Authors
!     -------
!     Adapted from the LMDTERRE code by R. Wordsworth (2009)
!     Original author Z. X. Li (1993)
!     
!==================================================================

#include "dimensions.h" 
#include "dimphys.h"
#include "comcstfi.h"

#include "fisice.h"
#include "callkeys.h"
#include "tracer.h"


!     Arguments
      REAL dtime ! intervalle du temps (s)
      REAL paprs(ngridmx,nlayermx+1) ! pression a inter-couche
      REAL pplay(ngridmx,nlayermx) ! pression au milieu de couche
      REAL t(ngridmx,nlayermx) ! temperature (K)
      REAL q(ngridmx,nlayermx) ! humidite specifique (kg/kg)
      REAL d_t(ngridmx,nlayermx) ! incrementation de la temperature (K)
      REAL d_q(ngridmx,nlayermx) ! incrementation de la vapeur d'eau
      REAL d_ql(ngridmx,nlayermx) ! incrementation de l'eau liquide
      REAL rneb(ngridmx,nlayermx) ! fraction nuageuse

!     Options du programme
      REAL ratqs   ! determine largeur de la distribution de vapeur
      PARAMETER (ratqs=0.2)

!     Variables locales
      REAL CBRT
      EXTERNAL CBRT
      INTEGER i, k
      REAL zt(ngridmx), zq(ngridmx)
      REAL zcond(ngridmx)
      REAL zdelq(ngridmx)
      REAL zqs(ngridmx), zdqs(ngridmx)

      REAL zcor(ngridmx), zdelta(ngridmx), zcvm5(ngridmx)
      REAL zx_q(ngridmx)

!     Initialisation des sorties
      DO k = 1, nlayermx
      DO i = 1, ngridmx
         d_t(i,k)=0.
         d_q(i,k)=0.
         d_ql(i,k)=0.
         rneb(i,k) = 0.0
      ENDDO
      ENDDO

!     Boucle verticale (du haut vers le bas)
      DO 9999 k = nlayermx, 1, -1

      DO i = 1, ngridmx
         zt(i)=t(i,k)
         zq(i)=q(i,k)
      ENDDO

!     Calculer la vapeur d'eau saturante et 
!     determiner la condensation partielle
      DO i = 1, ngridmx

         call watersat_2(zt(i),pplay(i,k),zqs(i))
         call watersat_grad(zt(i),zqs(i),zdqs(i))

         !IF (zt(i).LT.t_coup) THEN
         !   zqs(i) = qsats(zt(i))/pplay(i,k)
         !   zdqs(i) = dqsats(zt(i),zqs(i))
         !ELSE
         !   zqs(i) = qsatl(zt(i))/pplay(i,k)
         !   zdqs(i) = dqsatl(zt(i),zqs(i))
         !ENDIF

         zdelq(i) = ratqs * zq(i)
         rneb(i,k) = (zq(i)+zdelq(i)-zqs(i)) / (2.0*zdelq(i))
         zcond(i) = 0.0
         zx_q(i) = (zq(i)+zdelq(i)+zqs(i))/2.0
         if (rneb(i,k) .LE. 0.0) zx_q(i) = 0.0
         if (rneb(i,k) .GE. 1.0) zx_q(i) = zq(i)
         rneb(i,k) = MAX(0.0,MIN(1.0,rneb(i,k)))
         zcond(i) = MAX(0.0,zx_q(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i))
         zcond(i) = zcond(i)/dtime ! added by RDW

!     for varying particle size in rad tran and (possibly) sedimentation
!     to be dealt with in next version...

!         rice(i,k) = CBRT( 3*zcond(i)/( 4*Nmix_h2o*pi*rho_ice))
!         rice(i,k) = max(rice(i,k),1.e-16)

      ENDDO

!     Tendances de t et q
      DO i = 1, ngridmx
         d_q(i,k) = - zcond(i)
         d_ql(i,k) = zcond(i)
         d_t(i,k) = zcond(i)*RLVTT/cpp
      ENDDO

 9999 CONTINUE

      !print*,'rice=',rice
      !print*,'rneb=',rneb

      RETURN
      END