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

      subroutine simpleclouds(ngrid,nlay,ptimestep,
     &             pplay,pzlay,pt,pdt,
     &             pq,pdq,pdqcloud,pdtcloud,
     &             nq,tau,rice)
      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 "dimensions.h"
#include "dimphys.h"
#include "comcstfi.h"
#include "callkeys.h"
#include "tracer.h"
#include "comgeomfi.h"
#include "dimradmars.h"
c------------------------------------------------------------------
c     Arguments:
c     ---------
c     Inputs:
      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 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

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:

      LOGICAL firstcall
      DATA firstcall/.true./
      SAVE firstcall

      REAL CBRT
      EXTERNAL CBRT

      INTEGER ig,l

      REAL zq(ngridmx,nlayermx,nqmx)  ! local value of tracers
      REAL zq0(ngridmx,nlayermx,nqmx) ! local initial value of tracers
      REAL zt(ngridmx,nlayermx)       ! local value of temperature
      REAL zqsat(ngridmx,nlayermx)    ! 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(ngridmx,nlayermx) ! Typical dust geom. mean radius (m)
      REAL ccntyp(ngridmx,nlayermx)
                                      ! Typical dust number density (#/kg)
c     CCN reduction factor
c      REAL, PARAMETER :: ccn_factor = 4.5  !! comme TESTS_JB // 1. avant
      
      REAL Mcon_out(ngridmx,nlayermx) ! mass to be condensed (not dMice !!)

c-----------------------------------------------------------------------
c    1. initialisation
c    -----------------

c    On "update" la valeur de q(nqmx) (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)
        enddo
      enddo

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

c         Typical dust radius profile:
          rdusttyp(ig,l)= max(.8e-6*exp(-pzlay(ig,l)/18000.),1.e-9)

c         Typical CCN profile:
c         Corrected equation, following Montmessin et al. 2004
c           (N0=2e6 m-3 has been converted into N0=1.3e8 kg-1, otherwise
c            the equation for rice is not homogeneous...)
          ccntyp(ig,l)=
     &       1.3e+8*max(tau(ig,1),0.001)/0.1*exp(-pzlay(ig,l)/10000.)
c         The previously used profile was not correct:
c         ccntyp(ig,l)=( epaisseur(ig,l)/masse(ig,l) ) *
c    &       2.e+6/0.1*max(tau(ig,1),0.001)*exp(-pzlay(ig,l)/10000.)

c         Reduce number of nuclei
!         TEMPORAIRE : decrease the number of CCNs FF 04/200
!         reduction facteur 3
!         ccntyp(ig,l) = ccntyp(ig,l) / 27.
!         reduction facteur 2
!         ccntyp(ig,l) = ccntyp(ig,l) / 8.
c -----------------------------------------------------------------
          ccntyp(ig,l) = ccntyp(ig,l) / ccn_factor

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

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

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

      call watersat(ngridmx*nlayermx,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
          
          Mcon_out(ig,l) = dzq

c         Calcul du rayon moyen des particules de glace.
c         Hypothese : Dans une couche, la glace presente se
c         repartit uniformement autour du nbre de poussieres
c         dont le rayon moyen est prescrit par rdusttyp.
c         ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
          rice(ig,l)=max( CBRT ( (zq(ig,l,igcm_h2o_ice)/rho_ice
     &      +ccntyp(ig,l)*(4./3.)*pi*rdusttyp(ig,l)**3.)
     &      /(ccntyp(ig,l)*4./3.*pi) ), rdusttyp(ig,l)) 


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

c     Tendance finale
c     ~~~~~~~~~~~~~~~
      do l=1, nlay
        do ig=1,ngridmx
          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
        end do
      end do

c------------------------------------------------------------------
c     TEST_JBM
!      IF (ngrid.eq.1) THEN
!         call WRITEDIAGFI(ngrid,"mcond","h2o condensed mass","kg",1,
!     &                    Mcon_out)
!         call WRITEDIAGFI(ngrid,"rdusttyp","rdusttyp","m",1,
!     &                    rdusttyp)
!         call WRITEDIAGFI(ngrid,"ccntyp","ccntyp","kg-1",1,
!     &                    ccntyp)
!      ENDIF
c------------------------------------------------------------------
      return
      end