subroutine largescale(ptimestep, pplev, pplay, pt, pq,   &
                        pdt, pdq, pdtlsc, pdqvaplsc, pdqliqlsc, rneb)

      use watercommon_h, only : RLVTT, RCPD, RVTMP2,  &
          T_h2O_ice_clouds,T_h2O_ice_liq,Psat_water,Lcpdqsat_water
      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 ptimestep                 ! intervalle du temps (s)
      REAL pplev(ngridmx,nlayermx+1) ! pression a inter-couche
      REAL pplay(ngridmx,nlayermx)   ! pression au milieu de couche
      REAL pt(ngridmx,nlayermx)      ! temperature (K)
      real pq(ngridmx,nlayermx,nqmx) ! tracer mixing ratio (kg/kg)
      REAL pdt(ngridmx,nlayermx)     ! physical temperature tenedency (K/s)
      REAL pdq(ngridmx,nlayermx,nqmx)! physical tracer tenedency (K/s)
      REAL pdtlsc(ngridmx,nlayermx)  ! incrementation de la temperature (K)
      REAL pdqvaplsc(ngridmx,nlayermx) ! incrementation de la vapeur d'eau
      REAL pdqliqlsc(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 , nn
      INTEGER,PARAMETER :: niter=4
      REAL zt(ngridmx), zq(ngridmx)
      REAL zcond(ngridmx),zcond_iter
      REAL zdelq(ngridmx)
      REAL zqs(ngridmx), zdqs(ngridmx)
      REAL psat_tmp
      
! evaporation calculations
      REAL dqevap(ngridmx,nlayermx),dtevap(ngridmx,nlayermx)     
      REAL qevap(ngridmx,nlayermx,nqmx)
      REAL tevap(ngridmx,nlayermx)

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

!     GCM -----> subroutine variables, initialisation of outputs

      pdtlsc(1:ngridmx,1:nlayermx)  = 0.0
      pdqvaplsc(1:ngridmx,1:nlayermx)  = 0.0
      pdqliqlsc(1:ngridmx,1:nlayermx) = 0.0
      rneb(1:ngridmx,1:nlayermx) = 0.0


      ! Evaporate cloud water/ice
      call evap(ptimestep,pt,pq,pdq,pdt,dqevap,dtevap,qevap,tevap)
      ! note: we use qevap but not tevap in largescale/moistadj
            ! otherwise is a big mess


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

      zt(1:ngridmx)=pt(1:ngridmx,k)+(pdt(1:ngridmx,k)+dtevap(1:ngridmx,k))*ptimestep
      zq(1:ngridmx)=qevap(1:ngridmx,k,igcm_h2o_vap) !liquid water is included in qevap

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

         if(zt(i).le.15.) zt(i)=15.   ! check too low temperatures
!         call watersat(zt(i),pplay(i,k),zqs(i))
         call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
 
         zdelq(i) = ratqs * zq(i)
	 if(zq(i)+zdelq(i).lt.0.999) then
	    rneb(i,k) = (zq(i)+zdelq(i)-zqs(i)) / (2.0*zdelq(i))
	    zx_q(i) = (zq(i)+zdelq(i)+zqs(i))/2.0 !water vapor in cloudy sky
            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)))
	 else
	    if(zq(i).gt.zqs(i)) then
	       rneb(i,k)=1.
	       zx_q(i)=zq(i)
	    else
	       rneb(i,k)=0.
	       zx_q(i)=zqs(i)  !no condensation needed
	    Endif
	 Endif

!        iterative process to stabilize the scheme when large water amounts JL12
         zcond(i) = 0.0
         Do nn=1,niter  
!            call watersat_grad(zt(i),zqs(i),zdqs(i))
	    call Lcpdqsat_water(zt(i),pplay(i,k),psat_tmp,zqs(i),zdqs(i))
            zcond_iter = MAX(0.0,(zx_q(i)-zqs(i))*rneb(i,k)/(1.+zdqs(i)))	   
               !zcond always postive! cannot evaporate clouds!
               !this is why we must reevaporate before largescale
            zx_q(i) = zx_q(i) - zcond_iter
	    zcond(i) = zcond(i) + zcond_iter
	    zt(i) = zt(i) + zcond_iter*RLVTT/RCPD
            call Psat_water(zt(i),pplay(i,k),psat_tmp,zqs(i))
	 End do ! niter


         zcond(i) = zcond(i)/ptimestep ! added by RDW
      ENDDO

!     Tendances de t et q
         pdqvaplsc(1:ngridmx,k)  = dqevap(1:ngridmx,k) - zcond(1:ngridmx)
         pdqliqlsc(1:ngridmx,k) = - pdqvaplsc(1:ngridmx,k)
         pdtlsc(1:ngridmx,k)  = pdqliqlsc(1:ngridmx,k)*RLVTT/RCPD

   Enddo ! k= nlayermx, 1, -1
   
      !print*,'qsat=',zqs
      !print*,'q=',q
      !print*,'dq=',pdqvaplsc*ptimestep
      !print*,'dT in LS=',pdtlsc*ptimestep

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

      return
      end