Changeset 495 for LMDZ.3.3/branches/rel-LF/libf/phylmd/clouds_gno.F

Timestamp:

Mar 4, 2004, 4:11:16 PM (20 years ago)

Author:

lmdzadmin

Message:

Optimisation de differentes routines, IM, MAF, FH
LF

File:

• LMDZ.3.3/branches/rel-LF/libf/phylmd/clouds_gno.F (modified) (8 diffs)

LMDZ.3.3/branches/rel-LF/libf/phylmd/clouds_gno.F

@@ -48 +48 @@
       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 pi, u(klon), v(klon), erfcu(klon), erfcv(klon)
       REAL  xx1(klon), xx2(klon)
       real erf,kkk
+      real sqrtpi,sqrt2,zx1,zx2,exdel
 c lconv = true si le calcul a converge (entre autre si qsub < min_q)
        LOGICAL lconv(klon)
@@ -56 +57 @@
 
       pi = ACOS(-1.)
+      sqrtpi=sqrt(pi)
+      sqrt2=sqrt(2.)
 
       ptconv=.false.
@@ -123 +126 @@
           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.)))
+         zx1=-sqrt2*vmax(i)
+         zx2=SQRT(1.0+delta(i)/(vmax(i)**2.))
+         xx1(i)=zx1*(1.0-zx2)
+         xx2(i)=zx1*(1.0+zx2)
          xx(i) = 1.01 * xx1(i)
          if ( xx1(i) .GE. 0.0 ) xx(i) = 0.5*xx2(i)
@@ -143 +148 @@
         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.))
+          u(i) = delta(i)/(xx(i)*sqrt2) + xx(i)/(2.*sqrt2)
+          v(i) = delta(i)/(xx(i)*sqrt2) - xx(i)/(2.*sqrt2)
 
           IF ( v(i) .GT. vmax(i) ) THEN 
@@ -153 +158 @@
 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)))
+             exdel=beta(i)*EXP(delta(i))
+             aux(i) = 2.0*delta(i)*(1.-exdel)
+     :                       /(1.+exdel)
              if (aux(i).lt.0.) then
-                print*,'AUX(',i,',',k,')<0',aux(i),delta(i),beta(i)
+c                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)
+             block(i) = EXP(-v(i)*v(i)) / v(i) / sqrtpi
              dist(i) = 0.0
              fprime(i) = 1.0
@@ -169 +174 @@
 c -- erfv -> 1.0, use an asymptotic expression of erfv for v large:
 
-             erfu(i) = ERF(u(i))
+             erfcu(i) = 1.0-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.))
+             aux(i) = sqrtpi*erfcu(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)
+             block(i) = coeff(i) * EXP(-v(i)*v(i)) / v(i) / sqrtpi
              dist(i) = v(i) * aux(i) / coeff(i) - beta(i)
              fprime(i) = 2.0 / xx(i) * (v(i)**2.)
@@ -185 +190 @@
 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)
+           erfcu(i) = 1.0-ERF(u(i))
+           erfcv(i) = 1.0-ERF(v(i))
+           block(i) = erfcv(i)
+           dist(i) = erfcu(i) / erfcv(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)
+c              print*,'ATTENTION !!! xx(',i,') =', xx(i)
+c           print*,'ATTENTION !!! klon,ND,R,RS,QSUB,PTCONV,RATQSC,CLDF',
+c     .klon,ND,R(i,k),RS(i,k),QSUB(i,k),PTCONV(i,k),RATQSC(i,k),
+c     .CLDF(i,k)
+c              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)) )
+             fprime(i) = 2. /sqrtpi /xx(i) /erfcv(i)**2.
+     :           * (   erfcv(i)*v(i)*EXP(-zu2(i))
+     :               - erfcu(i)*u(i)*EXP(-zv2(i)) )
            endif
           ENDIF ! x
@@ -219 +224 @@
             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)