Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 529)
+++ /trunk/LMDZ.MARS/README	(revision 530)
@@ -1399,3 +1399,10 @@
 >> Update in deftank: callphys.def.outliers,run.def.1d; added traceur.def.scavenging
 
-
+== 15/02/12 == AS
+>> Optimization of watercycle new capabilities, e.g.:
+   ... in nuclea, removed 'if' from function fshape, modified **, decrease number of overall variables
+   ... in growthrate and newsedim, modified ** 
+   ... in improvedclouds, modified **, improved recursive calls to erf (now loop w/ 2 erf + 1 log instead of 4 erf + 4 log)
+>> Checked consistency with previous formulations
+>> Code 35% faster with a Valles Marineris mesoscale run in springtime (set imicro=25)
+   ... probably even better at cloudier seasons
Index: /trunk/LMDZ.MARS/libf/phymars/growthrate.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/growthrate.F	(revision 529)
+++ /trunk/LMDZ.MARS/libf/phymars/growthrate.F	(revision 530)
@@ -58,9 +58,11 @@
       k  = (0.17913 * temp - 13.9789) * 4.184e-4
 c     - Latent heat of h2o (J.kg-1)
-      Lv = (2834.3 - 0.28 * (temp-To) - 0.004 * (temp-To)**2 ) * 1.e+3
+      Lv = (2834.3 
+     &        - 0.28  * (temp-To) 
+     &        - 0.004 * (temp-To) * (temp-To) ) * 1.e+3
 
 c     - Constant to compute gas mean free path
 c     l= (T/P)*a, with a = (  0.707*8.31/(4*pi*molrad**2 * avogadro))
-      afactor = 0.707*rgp/(4 * pi* molco2**2  * nav)
+      afactor = 0.707*rgp/(4 * pi * molco2 * molco2 * nav)
 
 c     - Compute Dv, water vapor diffusion coefficient
@@ -69,5 +71,6 @@
 
       Dv = 1./3. * sqrt( 8*kbz*temp/(pi*mh2o/nav) )* kbz * temp / 
-     &   ( pi * pmid * (molco2+molh2o)**2 * sqrt(1.+mh2o/mco2) )
+     &   ( pi * pmid * (molco2+molh2o)*(molco2+molh2o) 
+     &        * sqrt(1.+mh2o/mco2) )
 
       knudsen = temp / pmid * afactor / rcrystal
@@ -76,5 +79,5 @@
 
 c     - Compute Rk
-      Rk = Lv**2 * rho_ice * mh2o / (k*rgp*temp**2.)
+      Rk = Lv**2 * rho_ice * mh2o / (k*rgp*temp*temp)
 c     - Compute Rd
       Rd = rgp * temp *rho_ice / (Dv*psat*mh2o)
Index: /trunk/LMDZ.MARS/libf/phymars/improvedclouds.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/improvedclouds.F	(revision 529)
+++ /trunk/LMDZ.MARS/libf/phymars/improvedclouds.F	(revision 530)
@@ -25,4 +25,5 @@
 c           (October 2011)
 c           T. Navarro, debug & correction (October-December 2011)
+c           A. Spiga, optimization (February 2012)
 c------------------------------------------------------------------
 #include "dimensions.h"
@@ -94,5 +95,5 @@
       REAL*8 Mo,No
       REAL*8 dN,dM,newvap
-      REAL*8 Rn, Rm, dev2
+      REAL*8 Rn, Rm, dev2, yeah, yeahn, yeahm
       REAL*8 n_aer(nbin_cld) ! number conc. of particle/each size bin
       REAL*8 m_aer(nbin_cld) ! mass mixing ratio of particle/each size bin
@@ -166,13 +167,13 @@
         rb_cld(1)  = rbmin_cld
         rad_cld(1) = rmin_cld
-        vol_cld(1) = 4./3. * dble(pi) * rmin_cld**3.
+        vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
 
         do i=1,nbin_cld-1
-          rad_cld(i+1)  = rad_cld(i) * vrat_cld**(1./3.)
+          rad_cld(i+1)  = rad_cld(i) * CBRT(vrat_cld) !!**(1./3.)
           vol_cld(i+1)  = vol_cld(i) * vrat_cld
         enddo
         
         do i=1,nbin_cld
-          rb_cld(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) *
+          rb_cld(i+1)= CBRT( (2.*vrat_cld) / (vrat_cld+1.) ) *
      &      rad_cld(i)
           dr_cld(i)  = rb_cld(i+1) - rb_cld(i)
@@ -285,16 +286,33 @@
         Mo = zq(ig,l,igcm_dust_mass)* tauscaling(ig)
         No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+ 1.e-30
+        dev2 = 1. / ( sqrt(2.) * sigma_ice )
         Rn = rdust(ig,l)
-        Rm = Rn * exp( 3. * sigma_ice**2. )
-        Rn = 1. / Rn
-        Rm = 1. / Rm
-        dev2 = 1. / ( sqrt(2.) * sigma_ice )
+        Rn = -dlog(Rn) 
+        Rm = Rn - 3. * sigma_ice*sigma_ice  
+        yeah = dlog(rb_cld(1))
+        yeahn = derf( (yeah+Rn) *dev2)
+        yeahm = derf( (yeah+Rm) *dev2)
         do i = 1, nbin_cld
-          n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)*Rn) * dev2 )
-     &                           -derf( dlog(rb_cld(i) * Rn) * dev2 ) )
-          m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)*Rm) * dev2 )
-     &                           -derf( dlog(rb_cld(i) * Rm) * dev2 ) )
+          n_aer(i) = -0.5 * No * yeahn !! this ith previously computed
+          m_aer(i) = -0.5 * Mo * yeahm !! this ith previously computed
+          yeah = dlog(rb_cld(i+1))     !! this (i+1)th now computed
+          yeahn = derf( (yeah+Rn) *dev2)
+          yeahm = derf( (yeah+Rm) *dev2)
+          n_aer(i) = n_aer(i) + 0.5 * No * yeahn
+          m_aer(i) = m_aer(i) + 0.5 * Mo * yeahm
         enddo
-                
+!!! MORE EFFICIENT COMPUTATIONALLY THAN
+!        Rn = rdust(ig,l)
+!        Rm = Rn * exp( 3. * sigma_ice**2. )
+!        Rn = 1. / Rn
+!        Rm = 1. / Rm
+!        do i = 1, nbin_cld
+!          n_aer(i) = 0.5 * No * ( derf( dlog(rb_cld(i+1)*Rn) * dev2 )
+!     &                           -derf( dlog(rb_cld(i) * Rn) * dev2 ) )
+!          m_aer(i) = 0.5 * Mo * ( derf( dlog(rb_cld(i+1)*Rm) * dev2 )
+!     &                           -derf( dlog(rb_cld(i) * Rm) * dev2 ) )
+!        enddo
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+        
 !        sumcheck = 0
 !        do i = 1, nbin_cld
@@ -346,5 +364,5 @@
           rhocloud(ig,l) = min(max(rhocloud(ig,l),rho_ice),rho_dust)
           rice(ig,l) =
-     &      ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
+     &      CBRT( Mo / No * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
 c          nuice(ig,l)=nuice_ref ! used for rad. transfer calculations
           if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
@@ -466,5 +484,5 @@
           rice_out(ig,l)=rice(ig,l)
           rice(ig,l) =
-     &      ( Mo / No * 0.75 / pi / rhocloud(ig,l) ) **(1./3.)
+     &      CBRT( Mo / No * 0.75 / pi / rhocloud(ig,l) ) !**(1./3.)
           if ((Mo.lt.1.e-20) .or. (No.le.1)) rice(ig,l) = 1.e-8
           rice_out(ig,l)=rice(ig,l)-rice_out(ig,l)
@@ -478,5 +496,6 @@
           rdust(ig,l)=min(max(rdust(ig,l),1.e-10),500.e-6)
           
-          rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)**3.,
+          rsedcloud(ig,l)=max( rice(ig,l)*(1.+nuice_sed)*
+     &                         (1.+nuice_sed)*(1.+nuice_sed),
      &                         rdust(ig,l) )
           rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
@@ -540,5 +559,4 @@
       
       drice_col(ig) = drice_col(ig)/icetot(ig)
-
 
       IF (ngrid.ne.1) THEN ! 3D
Index: /trunk/LMDZ.MARS/libf/phymars/newsedim.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/newsedim.F	(revision 529)
+++ /trunk/LMDZ.MARS/libf/phymars/newsedim.F	(revision 530)
@@ -91,5 +91,5 @@
 c       - Constant  to compute gas mean free path
 c        l= (T/P)*a, avec a = (  0.707*8.31/(4*pi*molrad**2 * avogadro))
-         a = 0.707*8.31/(4*3.1416* molrad**2  * 6.023e23)
+         a = 0.707*8.31/(4*3.1416* molrad*molrad * 6.023e23)
 
 c       - Correction to account for non-spherical shape (Murphy et al.  1990)
@@ -124,5 +124,5 @@
             endif  
 c           vstokes
-            vstokes(ig,l) = b * rhofall * rfall**2 *
+            vstokes(ig,l) = b * rhofall * rfall*rfall *
      &      (1 + 1.333* ( a*pt(ig,l)/pplev(ig,l) )/rfall)
 
Index: /trunk/LMDZ.MARS/libf/phymars/nuclea.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/nuclea.F	(revision 529)
+++ /trunk/LMDZ.MARS/libf/phymars/nuclea.F	(revision 530)
@@ -11,4 +11,5 @@
 *     Adapted for the LMD/GCM by J.-B. Madeleine      *
 *     (October 2011)                                  *
+*     Optimisation by A. Spiga (February 2012)        *  
 *******************************************************
 
@@ -34,12 +35,15 @@
       DOUBLE PRECISION gstar    ! # of molecules forming a critical embryo
       DOUBLE PRECISION fistar   ! Activation energy required to form a critical embryo (J)
-      DOUBLE PRECISION zeldov   ! Zeldovitch factor (no dim)
+!      DOUBLE PRECISION zeldov   ! Zeldovitch factor (no dim)
       DOUBLE PRECISION fshape   ! function defined at the end of the file
       DOUBLE PRECISION deltaf
 
 c     Ratio rstar/radius of the nucleating dust particle
-      double precision xratio
+c     double precision xratio
       
       double precision mtetalocal ! local mteta in double precision
+
+      double precision fshapesimple,zefshape
+
 
       integer i
@@ -51,5 +55,5 @@
 c     *************************************************
 
-      mtetalocal = mteta
+      mtetalocal = mteta  !! use mtetalocal for better performance
 
 cccccccccccccccccccccccccccccccccccccccccccccccccc
@@ -86,5 +90,8 @@
         nh2o   = ph2o / kbz / temp
         rstar  = 2. * sig(temp) * vo1 / (rgp*temp*dlog(sat))
-        gstar  = 4. * nav * pi * (rstar**3) / (3.*vo1)
+        gstar  = 4. * nav * pi * (rstar * rstar * rstar) / (3.*vo1)
+        
+        fshapesimple = (2.+mtetalocal)*(1.-mtetalocal)*(1.-mtetalocal)
+     &                   / 4.
 
 c       Loop over size bins
@@ -97,7 +104,12 @@
           endif
 
-          xratio   = rad_cld(i) / rstar
-          fistar = (4./3.*pi) * sig(temp) * (rstar**2.) *
-     &             fshape(mtetalocal,xratio)
+          if (rad_cld(i).gt.3000.*rstar) then
+            zefshape = fshapesimple
+          else
+            zefshape = fshape(mtetalocal,rad_cld(i)/rstar)
+          endif
+
+          fistar = (4./3.*pi) * sig(temp) * (rstar * rstar) * 
+     &             zefshape
           deltaf = (2.*desorp-surfdif-fistar)/
      &             (kbz*temp)
@@ -107,10 +119,11 @@
             nucrate(i) = 0.
           else
-            zeldov = sqrt ( fistar /
-     &               (3.*pi*kbz*temp*(gstar**2.)) )
-            nucrate(i)= zeldov * kbz * temp * rstar
+            nucrate(i)= sqrt ( fistar /
+     &               (3.*pi*kbz*temp*(gstar*gstar)) )
+     &                  * kbz * temp * rstar
      &                  * rstar * 4. * pi
-     &                  * ( nh2o*rad_cld(i) )**2.
-     &                  / ( fshape(mtetalocal,xratio) * nus * m0 )
+     &                  * ( nh2o*rad_cld(i) )
+     &                  * ( nh2o*rad_cld(i) )
+     &                  / ( zefshape * nus * m0 )
      &                  * dexp (deltaf)
           endif
@@ -137,16 +150,20 @@
 
       double precision cost,rap
-      double precision phi
-      double precision a,b,c
+      double precision yeah
 
-      phi = sqrt( 1. - 2.*cost*rap + rap**2 )
-      a = 1. + ( (1.-cost*rap)/phi )**3
-      b = (rap**3) * (2.-3.*(rap-cost)/phi+((rap-cost)/phi)**3)
-      c = 3. * cost * (rap**2) * ((rap-cost)/phi-1.)
+          !! PHI
+          yeah = sqrt( 1. - 2.*cost*rap + rap*rap )
+          !! FSHAPE = TERM A
+          fshape = (1.-cost*rap) / yeah
+          fshape = fshape * fshape * fshape
+          fshape = 1. + fshape
+          !! ... + TERM B
+          yeah = (rap-cost)/yeah
+          fshape = fshape + rap*rap*rap*(2.-3.*yeah+yeah*yeah*yeah)
+          !! ... + TERM C 
+          fshape = fshape + 3. * cost * rap * rap * (yeah-1.)
+          !! FACTOR 1/2
+          fshape = 0.5*fshape
 
-      fshape = 0.5*(a+b+c)
-
-      if (rap.gt.3000.) fshape = ((2.+cost)*(1.-cost)**2)/4.
-
-      return
+      return 
       end
