source: trunk/LMDZ.MARS/libf/phymars/simpleclouds.F @ 3026

      subroutine simpleclouds(ngrid,nlay,ptimestep,
     &             pplay,pzlay,pt,pdt,
     &             pq,pdq,pdqcloud,pdtcloud,
     &             nq,tau,rice)
      USE updaterad
      USE watersat_mod, ONLY: watersat
      use tracer_mod, only: igcm_h2o_vap, igcm_h2o_ice,
     &                      igcm_hdo_vap, igcm_hdo_ice,
     &                      qperemin
      USE comcstfi_h
      use dimradmars_mod, only: naerkind

      implicit none
c------------------------------------------------------------------
c  This routine is used to form clouds when a parcel of the GCM is
c    saturated. It is a simplified scheme, and there is almost no
c    microphysics involved. When the air is saturated, water-ice
c    clouds form on a fraction of the dust particles, specified by
c    the constant called "ccn_factor". There is no supersaturation,
c    and no nucleation rates computed. A more accurate scheme can
c    be found in the routine called "improvedclouds.F".

c  Modif de zq si saturation dans l'atmosphere
c  si zq(ig,l)> zqsat(ig,l) ->    zq(ig,l)=zqsat(ig,l)
c  Le test est effectue de bas en haut. L'eau condensee
c    (si saturation) est remise dans la couche en dessous.
c  L'eau condensee dans la couche du bas est deposee a la surface

c  Authors: Franck Montmessin (water ice scheme)
c           Francois Forget (changed nuclei density & outputs)
c           Ehouarn Millour (sept.2008, tracers are now handled
c                                   by name and not fixed index)
c           J.-B. Madeleine (developed a single routine called
c                            simpleclouds.F, and corrected calculations
c                            of the typical CCN profile, Oct. 2011)
c------------------------------------------------------------------
#include "callkeys.h"

c------------------------------------------------------------------
c     Arguments:
c     ---------
c     Inputs:
      INTEGER ngrid,nlay
      integer nq                 ! nombre de traceurs
      REAL ptimestep             ! pas de temps physique (s)
      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
      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(ngrid,naerkind)   ! Column dust optical depth at each point

c     Output:
      REAL rice(ngrid,nlay)      ! Ice mass mean radius (m)
                                 ! (r_c in montmessin_2004)
      real pdqcloud(ngrid,nlay,nq) ! tendance de la condensation
                                   !   H2O(kg/kg.s-1)
      REAL pdtcloud(ngrid,nlay)    ! tendance temperature due
                                   !   a la chaleur latente

c------------------------------------------------------------------
c     Local variables:

      REAL rhocloud(ngrid,nlay)  ! Cloud density (kg.m-3)

      INTEGER ig,l

      REAL zq(ngrid,nlay,nq)    ! local value of tracers
      REAL zq0(ngrid,nlay,nq)   ! local initial value of tracers
      REAL zt(ngrid,nlay)       ! local value of temperature
      REAL zqsat(ngrid,nlay)    ! saturation
      REAL*8 dzq                      ! masse de glace echangee (kg/kg)
      REAL lw                         !Latent heat of sublimation (J.kg-1) 
      REAL,PARAMETER :: To=273.15     ! reference temperature, T=273.15 K
      real rdusttyp(ngrid,nlay) ! Typical dust geom. mean radius (m)
      REAL ccntyp(ngrid,nlay)
                                      ! Typical dust number density (#/kg)
      REAL alpha_c(ngrid,nlay) !!MARGAUX: alpha_c as a spatial variable

c     CCN reduction factor
c      REAL, PARAMETER :: ccn_factor = 4.5  !! comme TESTS_JB // 1. avant
      REAL DoH_vap(ngrid,nlay)
      REAL DoH_ice(ngrid,nlay)
c-----------------------------------------------------------------------
c    1. initialisation
c    -----------------

c    On "update" la valeur de q(nq) (water vapor) et temperature.
c    On effectue qqes calculs preliminaires sur les couches : 

      do l=1,nlay
        do ig=1,ngrid
          zq(ig,l,igcm_h2o_vap)=
     &      pq(ig,l,igcm_h2o_vap)+pdq(ig,l,igcm_h2o_vap)*ptimestep
          zq(ig,l,igcm_h2o_vap)=max(zq(ig,l,igcm_h2o_vap),1.E-30) ! FF 12/2004 
          zq0(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)
          zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep

          zq(ig,l,igcm_h2o_ice)=
     &      pq(ig,l,igcm_h2o_ice)+pdq(ig,l,igcm_h2o_ice)*ptimestep
          zq(ig,l,igcm_h2o_ice)=max(zq(ig,l,igcm_h2o_ice),0.) ! FF 12/2004 
          zq0(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)

          if (hdo) then
          zq(ig,l,igcm_hdo_vap)=
     &      pq(ig,l,igcm_hdo_vap)+pdq(ig,l,igcm_hdo_vap)*ptimestep
          zq(ig,l,igcm_hdo_vap)=max(zq(ig,l,igcm_hdo_vap),1e-30) ! FF 12/2004 
          zq0(ig,l,igcm_hdo_vap)=zq(ig,l,igcm_hdo_vap)

          zq(ig,l,igcm_hdo_ice)=
     &      pq(ig,l,igcm_hdo_ice)+pdq(ig,l,igcm_hdo_ice)*ptimestep
          zq(ig,l,igcm_hdo_ice)=max(zq(ig,l,igcm_hdo_ice),1e-30) ! FF 12/2004 
          zq0(ig,l,igcm_hdo_ice)=zq(ig,l,igcm_hdo_ice)

          endif !hdo
        enddo
      enddo

      pdqcloud(1:ngrid,1:nlay,1:nq)=0
      pdtcloud(1:ngrid,1:nlay)=0

      alpha_c(1:ngrid,1:nlay)=0.

c     ----------------------------------------------
c
c     Rapport de melange a saturation dans la couche l : -------
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      call watersat(ngrid*nlay,zt,pplay,zqsat)

c     taux de condensation (kg/kg/s-1) dans les differentes couches
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

      do l=1,nlay
        do ig=1,ngrid

          if (zq(ig,l,igcm_h2o_vap).ge.zqsat(ig,l))then  !  Condensation
            dzq=zq(ig,l,igcm_h2o_vap)-zqsat(ig,l)

          elseif(zq(ig,l,igcm_h2o_vap).lt.zqsat(ig,l))then  ! Sublimation
            dzq=-min(zqsat(ig,l)-zq(ig,l,igcm_h2o_vap),
     &               zq(ig,l,igcm_h2o_ice))
          endif
            
c         Water Mass change
c         ~~~~~~~~~~~~~~~~~
          zq(ig,l,igcm_h2o_ice)=zq(ig,l,igcm_h2o_ice)+dzq
          zq(ig,l,igcm_h2o_vap)=zq(ig,l,igcm_h2o_vap)-dzq
        
        enddo ! of do ig=1,ngrid
      enddo ! of do l=1,nlay


c     Tendance finale
c     ~~~~~~~~~~~~~~~
      do l=1, nlay
        do ig=1,ngrid
          pdqcloud(ig,l,igcm_h2o_vap)=(zq(ig,l,igcm_h2o_vap)
     &                            -zq0(ig,l,igcm_h2o_vap))/ptimestep
          pdqcloud(ig,l,igcm_h2o_ice) =
     &      (zq(ig,l,igcm_h2o_ice) - zq0(ig,l,igcm_h2o_ice))/ptimestep

          lw=(2834.3-0.28*(zt(ig,l)-To)-0.004*(zt(ig,l)-To)**2)*1.e+3
          pdtcloud(ig,l)=-pdqcloud(ig,l,igcm_h2o_vap)*lw/cpp

          if (hdo) then
            if (pdqcloud(ig,l,igcm_h2o_ice).gt.0.0) then !condens

                if (hdofrac) then ! do we use fractionation?
c               alpha_c(ig,l) = exp(16288./zt(ig,l)**2.-9.34d-2)
                alpha_c(ig,l) = exp(13525./zt(ig,l)**2.-5.59d-2)  !Lamb
                else
                alpha_c(ig,l) = 1.d0
                endif
               
                if (zq0(ig,l,igcm_h2o_vap).gt.qperemin) then
                  pdqcloud(ig,l,igcm_hdo_ice)= 
     &              pdqcloud(ig,l,igcm_h2o_ice)*alpha_c(ig,l)*
     &         ( zq0(ig,l,igcm_hdo_vap)
     &                 /zq0(ig,l,igcm_h2o_vap) )
                else
                   pdqcloud(ig,l,igcm_hdo_ice)= 0.0
                endif

                pdqcloud(ig,l,igcm_hdo_ice) =
     &            min(pdqcloud(ig,l,igcm_hdo_ice),
     &              zq0(ig,l,igcm_hdo_vap)/ptimestep)

                pdqcloud(ig,l,igcm_hdo_vap)= 
     &               -pdqcloud(ig,l,igcm_hdo_ice)       

          else  ! sublimation

             if (zq0(ig,l,igcm_h2o_ice).gt.qperemin) then 
                pdqcloud(ig,l,igcm_hdo_ice)= 
     &               pdqcloud(ig,l,igcm_h2o_ice)*
     &      ( zq0(ig,l,igcm_hdo_ice)
     &              /zq0(ig,l,igcm_h2o_ice) )
             else
                pdqcloud(ig,l,igcm_hdo_ice)= 0.
             endif

              pdqcloud(ig,l,igcm_hdo_ice) =
     &          max(pdqcloud(ig,l,igcm_hdo_ice),
     &            -zq0(ig,l,igcm_hdo_ice)/ptimestep)

              pdqcloud(ig,l,igcm_hdo_vap)= 
     &             -pdqcloud(ig,l,igcm_hdo_ice)        

            endif ! condensation/sublimation

          endif ! hdo

        enddo ! of do ig=1,ngrid
      enddo ! of do l=1,nlay

c     ice crystal radius
      do l=1, nlay
        do ig=1,ngrid
          call updaterice_typ(zq(ig,l,igcm_h2o_ice),
     &       tau(ig,1),pzlay(ig,l),rice(ig,l))
        end do
      end do

c     if (hdo) then
c           CALL WRITEDIAGFI(ngrid,'alpha_c',
c    &                       'alpha_c',
c    &                       ' ',3,alpha_c) 
c     endif !hdo
c------------------------------------------------------------------
      return
      end