source: trunk/LMDZ.MARS/libf/phymars/simpleclouds_mod.F @ 4063

      module simpleclouds_mod
      
      implicit none
      
      contains

      subroutine simpleclouds(ngrid,nlay,ptimestep,
     &             pplay,pzlay,pt,pdt,
     &             pq,pdq,pdqcloud,pdtcloud,
     &             nq,tau,rice)
      USE updaterad, ONLY: updaterice_typ
      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, only: cpp
      use dimradmars_mod, only: naerkind
      use callkeys_mod, only: hdo, hdofrac

      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  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------------------------------------------------------------------
c     Arguments:
c     ---------
c     Inputs:
      integer,intent(in) :: ngrid ! number of atmospheric columns
      integer,intent(in) :: nlay ! number of atmospheric layers
      integer,intent(in) :: nq  ! number of tracers
      real,intent(in) :: ptimestep ! physics time step (s)
      real,intent(in) :: pplay(ngrid,nlay) ! pressure at mid-layer (Pa)
      real,intent(in) :: pzlay(ngrid,nlay) ! altitude of the layers (m)
      real,intent(in) :: pt(ngrid,nlay) ! input temperature at mid-layer (K)
      real,intent(in) :: pdt(ngrid,nlay) ! tendency on temperature from previous paramatrizations (K/s)
      real,intent(in) :: pq(ngrid,nlay,nq) ! input tracer mixing ratio (kg/kg)
      real,intent(in) :: pdq(ngrid,nlay,nq) ! tendency on tracers from previous parametrizations (kg/kg.s-1)
      real,intent(in) :: tau(ngrid,naerkind) ! Column dust optical depth in each column

c     Output:
      real,intent(out) :: rice(ngrid,nlay) ! Water ice mass mean radius (m)
                                   ! (r_c in montmessin_2004)
      real,intent(out) :: pdqcloud(ngrid,nlay,nq) ! tendencies due to H2O condensation
                                   !  and sublimation (kg/kg.s-1)
      real,intent(out) :: pdtcloud(ngrid,nlay) ! tendency on temperature due
                                   !  to latent heat (K/s)

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    Update values of water vapor and ice, and temperature.

      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     Compute saturation mixing ratio in each layer
c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

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

c     compute condensation rates (kg/kg/s-1) in each layer
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     Final tendency
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------------------------------------------------------------------
      end subroutine simpleclouds

      end module simpleclouds_mod