SUBROUTINE watercloud(ngrid,nlay, ptimestep, 
     &                pplev,pplay,pdpsrf,pzlev,pzlay,pt,pdt,
     &                pq,pdq,pdqcloud,pdqscloud,pdtcloud,
     &                nq,tau,tauscaling,rdust,rice,nuice,
     &                rsedcloud,rhocloud)
      IMPLICIT NONE

c=======================================================================
c  Water-ice cloud formation
c  
c  Includes two different schemes:
c    - A simplified scheme (see simpleclouds.F)
c    - An improved microphysical scheme (see improvedclouds.F)
c
c  Authors: Franck Montmessin, Francois Forget, Ehouarn Millour, 
c           J.-B. Madeleine
c
c  2004 - Oct. 2011
c=======================================================================

c-----------------------------------------------------------------------
c   declarations:
c   -------------

#include "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"
#include "tracer.h"
#include "comgeomfi.h"
#include "dimradmars.h"

c   Inputs:
c   ------

      INTEGER ngrid,nlay
      integer nq                 ! nombre de traceurs 
      REAL ptimestep             ! pas de temps physique (s)
      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
      REAL pdpsrf(ngrid)         ! tendance surf pressure
      REAL pzlev(ngrid,nlay+1)   ! altitude at layer boundaries
      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
      REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
      REAL pdt(ngrid,nlay)       ! tendance temperature des autres param.

      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
      real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)

      REAL tau(ngridmx,naerkind)   ! Column dust optical depth at each point
      REAL tauscaling(ngridmx)     ! Convertion factor for dust amount
      real rdust(ngridmx,nlayermx) ! Dust geometric mean radius (m)

c   Outputs:
c   -------

      real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation H2O(kg/kg.s-1)
      real pdqscloud(ngrid,nq)     ! flux en surface (kg.m-2.s-1)
      REAL pdtcloud(ngrid,nlay)    ! tendance temperature due
                                   !   a la chaleur latente

      REAL rice(ngrid,nlay)    ! Ice mass mean radius (m)
                               ! (r_c in montmessin_2004)
      REAL nuice(ngrid,nlay)   ! Estimated effective variance
                               !   of the size distribution
      real rsedcloud(ngridmx,nlayermx) ! Cloud sedimentation radius
      real rhocloud(ngridmx,nlayermx)  ! Cloud density (kg.m-3)

c   local:
c   ------

      INTEGER ig,l
      LOGICAL,SAVE :: firstcall=.true.

c    ** un petit test de coherence
c       --------------------------

      IF (firstcall) THEN
        IF(ngrid.NE.ngridmx) THEN
            PRINT*,'STOP dans watercloud'
            PRINT*,'probleme de dimensions :'
            PRINT*,'ngrid  =',ngrid
            PRINT*,'ngridmx  =',ngridmx
            STOP
        ENDIF
         
        if (nq.gt.nqmx) then
           write(*,*) 'stop in watercloud (nq.gt.nqmx)!'
           write(*,*) 'nq=',nq,' nqmx=',nqmx
           stop
        endif
         
        write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
        write(*,*) "            igcm_h2o_ice=",igcm_h2o_ice

        firstcall=.false.
      ENDIF ! of IF (firstcall)

c     Main call to the different cloud schemes:
      IF (microphys) THEN
        CALL improvedclouds(ngrid,nlay,ptimestep,
     &             pplev,pplay,pt,pdt,
     &             pq,pdq,pdqcloud,pdqscloud,pdtcloud,
     &             nq,tauscaling,rdust,rice,nuice,
     &             rsedcloud,rhocloud)
      ELSE
        CALL simpleclouds(ngrid,nlay,ptimestep,
     &             pplev,pplay,pzlev,pzlay,pt,pdt,
     &             pq,pdq,pdqcloud,pdqscloud,pdtcloud,
     &             nq,tau,rice,nuice,rsedcloud)
      ENDIF

c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
c     Then that should not affect the ice particle radius
      do ig=1,ngridmx
        if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
          if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
     &    rice(ig,2)=rice(ig,3) 
          rice(ig,1)=rice(ig,2)
        end if
      end do

c=======================================================================

!!!!!!!!!! FOR PHOTOCHEMISTRY, REIMPLEMENT output of surfdust/surfice
!!      if (photochem) then
!!c        computation of dust and ice surface area (micron2/cm3)
!!c        for heterogeneous chemistry
!!
!!         do l = 1,nlay
!!            do ig = 1,ngrid
!!c                                                                      
!!c              npart: number density of ccn in #/cm3    
!!c                                                     
!!               npart(ig,l) = 1.e-6*ccn(ig,l)  
!!     $                       *masse(ig,l)/epaisseur(ig,l) 
!!c
!!c              dust and ice surface area
!!c                                                                  
!!               surfdust(ig,l) = npart(ig,l)*4.*pi*1.e12*rdust(ig,l)**2
!!c                                                                  
!!               if (rice(ig,l) .ge. rdust(ig,l)) then                   
!!                  surfice(ig,l)  = npart(ig,l)*4.*pi*1.e12*rice(ig,l)**2
!!                  surfdust(ig,l) = 0.
!!               else                                                    
!!                  surfice(ig,l) = 0.                                 
!!               end if                                             
!!            end do      ! of do ig=1,ngrid
!!         end do         ! of do l=1,nlay
!!      end if            ! of photochem


      RETURN
      END