Index: /LMDZ.3.3/trunk/libf/phylmd/clouds_gno.F
===================================================================
--- /LMDZ.3.3/trunk/libf/phylmd/clouds_gno.F	(revision 381)
+++ /LMDZ.3.3/trunk/libf/phylmd/clouds_gno.F	(revision 381)
@@ -0,0 +1,250 @@
+C
+C================================================================================
+C
+      SUBROUTINE CLOUDS_GNO(klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF)
+      IMPLICIT NONE
+C     
+C--------------------------------------------------------------------------------
+C
+C Inputs:
+C
+C  ND----------: Number of vertical levels
+C  R--------ND-: Domain-averaged mixing ratio of total water 
+C  RS-------ND-: Mean saturation humidity mixing ratio within the gridbox
+C  QSUB-----ND-: Mixing ratio of condensed water within clouds associated
+C                with SUBGRID-SCALE condensation processes (here, it is
+C                predicted by the convection scheme)
+C Outputs:
+C
+C  PTCONV-----ND-: Point convectif = TRUE
+C  RATQSC-----ND-: Largeur normalisee de la distribution
+C  CLDF-----ND-: Fraction nuageuse
+C
+C--------------------------------------------------------------------------------
+
+
+      INTEGER klon,ND
+      REAL R(klon,ND),  RS(klon,ND), QSUB(klon,ND)
+      LOGICAL PTCONV(klon,ND)
+      REAL RATQSC(klon,ND)
+      REAL CLDF(klon,ND)
+
+c -- parameters controlling the iteration:
+c --    nmax    : maximum nb of iterations (hopefully never reached)
+c --    epsilon : accuracy of the numerical resolution 
+c --    vmax    : v-value above which we use an asymptotic expression for ERF(v)
+
+      INTEGER nmax
+      PARAMETER ( nmax = 10) 
+      REAL epsilon, vmax0, vmax(klon)
+      PARAMETER ( epsilon = 0.02, vmax0 = 2.0 ) 
+
+      REAL min_mu, min_Q
+      PARAMETER ( min_mu =  1.e-12, min_Q=1.e-12 )
+     
+      INTEGER i,K, n, m
+      REAL mu(klon), qsat(klon), delta(klon), beta(klon) 
+      real zu2(klon),zv2(klon)
+      REAL xx(klon), aux(klon), coeff(klon), block(klon)
+      REAL  dist(klon), fprime(klon), det(klon)
+      REAL pi, u(klon), v(klon), erfu(klon), erfv(klon)
+      REAL  xx1(klon), xx2(klon)
+      real erf,kkk
+c lconv = true si le calcul a converge (entre autre si qsub < min_q)
+       LOGICAL lconv(klon)
+
+
+      pi = ACOS(-1.)
+
+      ptconv=.false.
+      ratqsc=0.
+
+
+      DO 500 K = 1, ND
+
+                                    do i=1,klon ! vector
+      mu(i) = R(i,K)
+      mu(i) = MAX(mu(i),min_mu)
+      qsat(i) = RS(i,K) 
+      qsat(i) = MAX(qsat(i),min_mu)
+      delta(i) = log(mu(i)/qsat(i))
+                                    enddo ! vector
+
+C
+C ***          There is no subgrid-scale condensation;        ***
+C ***   the scheme becomes equivalent to an "all-or-nothing"  *** 
+C ***             large-scale condensation scheme.            ***
+C
+
+C
+C ***     Some condensation is produced at the subgrid-scale       ***
+C ***                                                              ***
+C ***       PDF = generalized log-normal distribution (GNO)        ***
+C ***   (k<0 because a lower bound is considered for the PDF)      ***
+C ***                                                              ***
+C ***  -> Determine x (the parameter k of the GNO PDF) such        ***
+C ***  that the contribution of subgrid-scale processes to         ***
+C ***  the in-cloud water content is equal to QSUB(K)              ***
+C ***  (equations (13), (14), (15) + Appendix B of the paper)      ***
+C ***                                                              ***
+C ***    Here, an iterative method is used for this purpose        ***
+C ***    (other numerical methods might be more efficient)         ***
+C ***                                                              ***
+C ***          NB: the "error function" is called ERF              ***
+C ***                 (ERF in double precision)                   ***
+C
+
+c  On commence par eliminer les cas pour lesquels on n'a pas
+c  suffisamment d'eau nuageuse.
+
+                                    do i=1,klon ! vector
+
+      IF ( QSUB(i,K) .lt. min_Q ) THEN
+        ptconv(i,k)=.false.
+        ratqsc(i,k)=0.
+        lconv(i)  = .true.
+
+c   Rien on a deja initialise
+
+      ELSE 
+
+        lconv(i)  = .FALSE. 
+        vmax(i) = vmax0
+
+        beta(i) = QSUB(i,K)/mu(i) + EXP( -MIN(0.0,delta(i)) )
+
+c --  roots of equation v > vmax:
+
+        det(i) = delta(i) + vmax(i)**2.
+        if (det(i).LE.0.0) vmax(i) = vmax0 + 1.0
+        det(i) = delta(i) + vmax(i)**2.
+
+        if (det(i).LE.0.) then
+          xx(i) = -0.0001
+        else 
+         xx1(i)=-SQRT(2.)*vmax(i)*(1.0-SQRT(1.0+delta(i)/(vmax(i)**2.)))
+         xx2(i)=-SQRT(2.)*vmax(i)*(1.0+SQRT(1.0+delta(i)/(vmax(i)**2.)))
+         xx(i) = 1.01 * xx1(i)
+         if ( xx1(i) .GE. 0.0 ) xx(i) = 0.5*xx2(i)
+        endif
+        if (delta(i).LT.0.) xx(i) = -0.5*SQRT(log(2.)) 
+
+      ENDIF
+
+                                    enddo       ! vector
+
+c----------------------------------------------------------------------
+c   Debut des nmax iterations pour trouver la solution.
+c----------------------------------------------------------------------
+
+      DO n = 1, nmax 
+
+                                    do i=1,klon ! vector
+        if (.not.lconv(i)) then
+
+          u(i) = delta(i)/(xx(i)*sqrt(2.)) + xx(i)/(2.*sqrt(2.))
+          v(i) = delta(i)/(xx(i)*sqrt(2.)) - xx(i)/(2.*sqrt(2.))
+
+          IF ( v(i) .GT. vmax(i) ) THEN 
+
+            IF (     ABS(u(i))  .GT. vmax(i) 
+     :          .AND.  delta(i) .LT. 0. ) THEN
+
+c -- use asymptotic expression of erf for u and v large:
+c ( -> analytic solution for xx )
+
+             aux(i) = 2.0*delta(i)*(1.-beta(i)*EXP(delta(i)))
+     :                       /(1.+beta(i)*EXP(delta(i)))
+             if (aux(i).lt.0.) then
+                print*,'AUX(',i,',',k,')<0',aux(i),delta(i),beta(i)
+                aux(i)=0.
+             endif
+             xx(i) = -SQRT(aux(i))
+             block(i) = EXP(-v(i)*v(i)) / v(i) / SQRT(pi)
+             dist(i) = 0.0
+             fprime(i) = 1.0
+
+            ELSE
+
+c -- erfv -> 1.0, use an asymptotic expression of erfv for v large:
+
+             erfu(i) = ERF(u(i))
+c  !!! ATTENTION : rajout d'un seuil pour l'exponentiel
+             aux(i) = SQRT(pi)*(1.0-erfu(i))*EXP(min(v(i)*v(i),100.))
+             coeff(i) = 1.0 - 1./2./(v(i)**2.) + 3./4./(v(i)**4.)
+             block(i) = coeff(i) * EXP(-v(i)*v(i)) / v(i) / SQRT(pi)
+             dist(i) = v(i) * aux(i) / coeff(i) - beta(i)
+             fprime(i) = 2.0 / xx(i) * (v(i)**2.)
+     :           * ( coeff(i)*EXP(-delta(i)) - u(i) * aux(i) )
+     :           / coeff(i) / coeff(i)
+            
+            ENDIF ! ABS(u)
+
+          ELSE
+
+c -- general case:
+
+           erfu(i) = ERF(u(i))
+           erfv(i) = ERF(v(i))
+           block(i) = 1.0-erfv(i)
+           dist(i) = (1.0 - erfu(i)) / (1.0 - erfv(i)) - beta(i)
+           zu2(i)=u(i)*u(i)
+           zv2(i)=v(i)*v(i)
+           if(zu2(i).gt.20..or. zv2(i).gt.20.) then
+              print*,'ATTENTION !!! xx(',i,') =', xx(i)
+           print*,'ATTENTION !!! klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF',
+     .klon,ND,R(i,k),RS(i,k),QSUB(i,k),PTCONV(i,k),RATQSC(i,k),
+     .CLDF(i,k)
+              print*,'ATTENTION !!! zu2 zv2 =',zu2(i),zv2(i)
+              zu2(i)=20.
+              zv2(i)=20.
+             fprime(i) = 0.
+           else
+             fprime(i) = 2. /SQRT(pi) /xx(i) /(1.0-erfv(i))**2.
+     :           * (   (1.0-erfv(i))*v(i)*EXP(-zu2(i))
+     :               - (1.0-erfu(i))*u(i)*EXP(-zv2(i)) )
+           endif
+          ENDIF ! x
+
+c -- test numerical convergence:
+
+c         print*,'avant test ',i,k,lconv(i),u(i),v(i)
+          if ( ABS(dist(i)/beta(i)) .LT. epsilon ) then 
+c           print*,'v-u **2',(v(i)-u(i))**2
+c           print*,'exp v-u **2',exp((v(i)-u(i))**2)
+            ptconv(i,K) = .TRUE. 
+            lconv(i)=.true.
+c  borne pour l'exponentielle
+            ratqsc(i,k)=min(2.*(v(i)-u(i))**2,20.)
+            ratqsc(i,k)=sqrt(exp(ratqsc(i,k))-1.)
+            CLDF(i,K) = 0.5 * block(i)
+cccccccccccccccccccccccc
+c           kkk=-sqrt(log(1.+ratqsc(i,k)**2))
+c           u(i)=delta(i)/(kkk*sqrt(2.))-kkk/(2.*sqrt(2.))
+c           v(i)=delta(i)/(kkk*sqrt(2.))+kkk/(2.*sqrt(2.))
+c           erfu(i)=erf(u(i))
+c           erfv(i)=erf(v(i))
+c           print*,'SIG ',k,qsub(i,k)
+c    s      ,mu(i)*((1.-erfv(i))/(1.-erfu(i))-qsat(i)/mu(i))
+c    s      ,0.5*erfu(i)
+cccccccccccccccccccccccc
+          else
+            xx(i) = xx(i) - dist(i)/fprime(i)
+          endif
+c         print*,'apres test ',i,k,lconv(i)
+
+        endif ! lconv
+                                    enddo       ! vector
+
+c----------------------------------------------------------------------
+c   Fin des nmax iterations pour trouver la solution.
+        ENDDO ! n
+c----------------------------------------------------------------------
+
+500   CONTINUE  ! K
+
+       RETURN
+       END
+
+
+
