Index: /LMDZ5/trunk/libf/phylmd/conf_phys_m.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/conf_phys_m.F90	(revision 1711)
+++ /LMDZ5/trunk/libf/phylmd/conf_phys_m.F90	(revision 1712)
@@ -18,5 +18,5 @@
                        iflag_cldcon, &
                        iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs, &
-  		       ok_ade, ok_aie, aerosol_couple, &
+  		       ok_ade, ok_aie, ok_cdnc, aerosol_couple, &
                        flag_aerosol, new_aod, &
                        bl95_b0, bl95_b1,&
@@ -60,4 +60,5 @@
 ! ok_instan:  sorties instantanees
 ! ok_ade, ok_aie: apply or not aerosol direct and indirect effects
+! ok_cdnc, ok cloud droplet number concentration
 ! bl95_b*: parameters in the formula to link CDNC to aerosol mass conc 
 !
@@ -70,5 +71,5 @@
   logical              :: ok_LES
   LOGICAL              :: callstats
-  LOGICAL              :: ok_ade, ok_aie, aerosol_couple
+  LOGICAL              :: ok_ade, ok_aie, ok_cdnc, aerosol_couple
   INTEGER              :: flag_aerosol
   LOGICAL              :: new_aod
@@ -84,5 +85,5 @@
   logical,SAVE        :: ok_LES_omp   
   LOGICAL,SAVE        :: callstats_omp
-  LOGICAL,SAVE        :: ok_ade_omp, ok_aie_omp, aerosol_couple_omp
+  LOGICAL,SAVE        :: ok_ade_omp, ok_aie_omp, ok_cdnc_omp, aerosol_couple_omp
   INTEGER, SAVE       :: flag_aerosol_omp
   LOGICAL, SAVE       :: new_aod_omp
@@ -272,4 +273,12 @@
   call getin('ok_aie', ok_aie_omp)
 
+!
+!Config Key  = ok_cdnc
+!Config Desc = ok cloud droplet number concentration
+!Config Def  = .false.
+!Config Help = Used in newmicro.F
+!
+  ok_cdnc_omp = .false.
+  call getin('ok_cdnc', ok_cdnc_omp)
 !
 !Config Key  = aerosol_couple
@@ -1677,4 +1686,5 @@
     ok_ade = ok_ade_omp
     ok_aie = ok_aie_omp
+    ok_cdnc = ok_cdnc_omp
     aerosol_couple = aerosol_couple_omp
     flag_aerosol=flag_aerosol_omp
@@ -1774,4 +1784,9 @@
        END IF
     END IF
+
+! ok_cdnc must be set to y if ok_aie is activated
+    IF (ok_aie .AND. .NOT. ok_cdnc) THEN
+       CALL abort_gcm('conf_phys', 'ok_cdnc must be set to y if ok_aie is activated',1)
+    ENDIF 
 
 !$OMP MASTER
Index: /LMDZ5/trunk/libf/phylmd/etat0_netcdf.F90
===================================================================
--- /LMDZ5/trunk/libf/phylmd/etat0_netcdf.F90	(revision 1711)
+++ /LMDZ5/trunk/libf/phylmd/etat0_netcdf.F90	(revision 1712)
@@ -99,5 +99,5 @@
   REAL    :: dummy
   LOGICAL :: ok_newmicro, ok_journe, ok_mensuel, ok_instan, ok_hf
-  LOGICAL :: ok_LES, ok_ade, ok_aie, aerosol_couple, new_aod, callstats
+  LOGICAL :: ok_LES, ok_ade, ok_aie, ok_cdnc, aerosol_couple, new_aod, callstats
   INTEGER :: iflag_radia, flag_aerosol
   REAL    :: bl95_b0, bl95_b1, fact_cldcon, facttemps, ratqsbas, ratqshaut
@@ -136,5 +136,5 @@
                    iflag_cldcon,                                        &
                    iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,            &
-                   ok_ade, ok_aie, aerosol_couple,                      &
+                   ok_ade, ok_aie, ok_cdnc, aerosol_couple,             &
                    flag_aerosol, new_aod,                               &
                    bl95_b0, bl95_b1,                                    &
Index: /LMDZ5/trunk/libf/phylmd/newmicro.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/newmicro.F	(revision 1711)
+++ /LMDZ5/trunk/libf/phylmd/newmicro.F	(revision 1712)
@@ -2,15 +2,13 @@
 
 
-
 !     
-      SUBROUTINE newmicro (paprs, pplay,ok_newmicro,
+      SUBROUTINE newmicro (ok_cdnc, bl95_b0, bl95_b1, 
+     .                  paprs, pplay,
      .                  t, pqlwp, pclc, pcltau, pclemi,
      .                  pch, pcl, pcm, pct, pctlwp,
-     s                  xflwp, xfiwp, xflwc, xfiwc,
-     e                  ok_aie, 
-     e                  mass_solu_aero, mass_solu_aero_pi, 
-     e                  bl95_b0, bl95_b1,
-     s                  cldtaupi, re, fl, reliq, reice)
-
+     .                  xflwp, xfiwp, xflwc, xfiwc,
+     .                  mass_solu_aero, mass_solu_aero_pi, 
+     .                  pcldtaupi, re, fl, reliq, reice)
+c
       USE dimphy
       USE phys_local_var_mod, only: scdnc,cldncl,reffclwtop,lcc,
@@ -21,61 +19,68 @@
 c======================================================================
 c Auteur(s): Z.X. Li (LMD/CNRS) date: 19930910
+c            O.   Boucher (LMD/CNRS) mise a jour en 201212
 c Objet: Calculer epaisseur optique et emmissivite des nuages
 c======================================================================
 c Arguments:
+c ok_cdnc-input-L-flag pour calculer les rayons a partir des aerosols
+c
 c t-------input-R-temperature
-c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere (kg/kg)
+c pqlwp---input-R-eau liquide nuageuse dans l'atmosphere dans la partie nuageuse (kg/kg)
 c pclc----input-R-couverture nuageuse pour le rayonnement (0 a 1)
-c 
-c ok_aie--input-L-apply aerosol indirect effect or not
 c mass_solu_aero-----input-R-total mass concentration for all soluble aerosols[ug/m^3]
-c mass_solu_aero_pi--input-R-dito, pre-industrial value
-c bl95_b0-input-R-a parameter, may be varied for tests (s-sea, l-land)
-c bl95_b1-input-R-a parameter, may be varied for tests (    -"-      )
+c mass_solu_aero_pi--input-R-ditto, pre-industrial value
+c
+c bl95_b0-input-R-a PARAMETER, may be varied for tests (s-sea, l-land)
+c bl95_b1-input-R-a PARAMETER, may be varied for tests (    -"-      )
 c      
-c cldtaupi-output-R-pre-industrial value of cloud optical thickness, 
-c                   needed for the diagnostics of the aerosol indirect 
-c                   radiative forcing (see radlwsw)
 c re------output-R-Cloud droplet effective radius multiplied by fl [um]
 c fl------output-R-Denominator to re, introduced to avoid problems in
 c                  the averaging of the output. fl is the fraction of liquid
 c                  water clouds within a grid cell           
+c
 c pcltau--output-R-epaisseur optique des nuages
 c pclemi--output-R-emissivite des nuages (0 a 1)
+c pcldtaupi-output-R-pre-industrial value of cloud optical thickness, 
+c
+c pcl-output-R-2D low-level cloud cover
+c pcm-output-R-2D mid-level cloud cover
+c pch-output-R-2D high-level cloud cover
+c pct-output-R-2D total cloud cover
 c======================================================================
 C
 #include "YOMCST.h"
-c
-cym#include "dimensions.h"
-cym#include "dimphy.h"
 #include "nuage.h"
-cIM cf. CR: include pour NOVLP et ZEPSEC
 #include "radepsi.h"
 #include "radopt.h"
+
 c choix de l'hypothese de recouvrememnt nuageuse
-      LOGICAL RANDOM,MAXIMUM_RANDOM,MAXIMUM
-      parameter (RANDOM=.FALSE., MAXIMUM_RANDOM=.TRUE., MAXIMUM=.FALSE.)
+      LOGICAL RANDOM, MAXIMUM_RANDOM, MAXIMUM
+      PARAMETER (RANDOM=.FALSE., MAXIMUM_RANDOM=.TRUE., MAXIMUM=.FALSE.)
+c
       LOGICAL, SAVE :: FIRST=.TRUE.
 !$OMP THREADPRIVATE(FIRST)
-c Hypoyhese de recouvrement : MAXIMUM_RANDOM
       INTEGER flag_max
-      REAL phase3d(klon, klev),dh(klon, klev),pdel(klon, klev),
-     .     zrho(klon, klev)
-      REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon)
+c
+c threshold PARAMETERs
       REAL thres_tau,thres_neb
       PARAMETER (thres_tau=0.3, thres_neb=0.001)
-      REAL t_tmp
-      REAL gravit
-      PARAMETER (gravit=9.80616)  !m/s2
-      REAL pqlwpcon(klon, klev), pqlwpstra(klon, klev) 
-c
-      REAL paprs(klon,klev+1), pplay(klon,klev)
+c
+      REAL phase3d(klon, klev)
+      REAL tcc(klon), ftmp(klon), lcc_integrat(klon), height(klon)
+c
+      REAL paprs(klon,klev+1)
+      REAL pplay(klon,klev)
       REAL t(klon,klev)
-c
       REAL pclc(klon,klev)
       REAL pqlwp(klon,klev)
-      REAL pcltau(klon,klev), pclemi(klon,klev)
-c
-      REAL pct(klon), pctlwp(klon), pch(klon), pcl(klon), pcm(klon)
+      REAL pcltau(klon,klev)
+      REAL pclemi(klon,klev)
+      REAL pcldtaupi(klon, klev)
+c
+      REAL pct(klon)
+      REAL pcl(klon)
+      REAL pcm(klon)
+      REAL pch(klon)
+      REAL pctlwp(klon)
 c
       LOGICAL lo
@@ -85,45 +90,34 @@
 !      PARAMETER (cetahb = 0.45, cetamb = 0.80)
 ! Remplacer
-!cetahb*paprs(i,1) par  prmhc
-!cetamb*paprs(i,1) par  prlmc 
-      REAL prmhc    ! Pressure between medium and high level cloud
-      REAL prlmc    ! Pressure between low and medium level cloud
+! cetahb*paprs(i,1) par  prmhc
+! cetamb*paprs(i,1) par  prlmc 
+      REAL prmhc    ! Pressure between medium and high level cloud in Pa
+      REAL prlmc    ! Pressure between low and medium level cloud in Pa
       PARAMETER (prmhc = 440.*100., prlmc = 680.*100.)
-
 C
       INTEGER i, k
-cIM: 091003   REAL zflwp, zradef, zfice, zmsac
-      REAL zflwp(klon), zradef, zfice, zmsac
-cIM: 091003 rajout
       REAL xflwp(klon), xfiwp(klon)
       REAL xflwc(klon,klev), xfiwc(klon,klev)
 c
-      REAL radius, rad_chaud
-cc      PARAMETER (rad_chau1=13.0, rad_chau2=9.0, rad_froid=35.0)
-ccc      PARAMETER (rad_chaud=15.0, rad_froid=35.0)
-c sintex initial      PARAMETER (rad_chaud=10.0, rad_froid=30.0)
-      REAL coef, coef_froi, coef_chau
+      REAL radius
+c
+      REAL coef_froi, coef_chau
       PARAMETER (coef_chau=0.13, coef_froi=0.09)
+c
       REAL seuil_neb
       PARAMETER (seuil_neb=0.001)
+c
       INTEGER nexpo ! exponentiel pour glace/eau
       PARAMETER (nexpo=6)
-ccc      PARAMETER (nexpo=1)
-
-c -- sb:
-      logical ok_newmicro
-c     parameter (ok_newmicro=.FALSE.)
-cIM: 091003   real rel, tc, rei, zfiwp
-      real rel, tc, rei, zfiwp(klon)
-      real k_liq, k_ice0, k_ice, DF
-      parameter (k_liq=0.0903, k_ice0=0.005) ! units=m2/g
-      parameter (DF=1.66) ! diffusivity factor
-c sb --
+c      PARAMETER (nexpo=1)
+
+      REAL rel, tc, rei
+      REAL k_ice0, k_ice, DF
+      PARAMETER (k_ice0=0.005) ! units=m2/g
+      PARAMETER (DF=1.66) ! diffusivity factor
+c
 cjq for the aerosol indirect effect
 cjq introduced by Johannes Quaas (quaas@lmd.jussieu.fr), 27/11/2003
 cjq      
-      LOGICAL ok_aie            ! Apply AIE or not?
-      LOGICAL ok_a1lwpdep       ! a1 LWP dependent?
-      
       REAL mass_solu_aero(klon, klev)    ! total mass concentration for all soluble aerosols [ug m-3]
       REAL mass_solu_aero_pi(klon, klev) ! - " - (pre-industrial value)
@@ -135,7 +129,7 @@
       REAL fl(klon, klev)       ! xliq * rneb (denominator to re; fraction of liquid water clouds within the grid cell)
       
+      LOGICAL ok_cdnc
       REAL bl95_b0, bl95_b1     ! Parameter in B&L 95-Formula
       
-      REAL cldtaupi(klon, klev) ! pre-industrial cloud opt thickness for diag
 cjq-end    
 cIM cf. CR:parametres supplementaires
@@ -145,14 +139,10 @@
       REAL zcloudm(klon) 
       REAL zcloudl(klon) 
-
-
-c **************************
-c *                        *
-c * DEBUT PARTIE OPTIMISEE *
-c *                        *
-c **************************
-
-      REAL diff_paprs(klon, klev), zfice1, zfice2(klon, klev)
-      REAL rad_chaud_tab(klon, klev), zflwp_var, zfiwp_var
+      REAL rhodz(klon, klev)  !--rho*dz pour la couche
+      REAL zrho(klon, klev)   !--rho pour la couche
+      REAL dh(klon, klev)     !--dz pour la couche
+      REAL zfice(klon, klev)
+      REAL rad_chaud(klon, klev) !--rayon pour les nuages chauds
+      REAL zflwp_var, zfiwp_var
       REAL d_rei_dt
 
@@ -171,286 +161,213 @@
 ! Pour retrouver les résultats numériques de la version d'origine,
 ! on impose 0.71 quand on est proche de 0.71
-
+c
       d_rei_dt=(rei_max-rei_min)/81.4
       if (abs(d_rei_dt-0.71)<1.e-4) d_rei_dt=0.71
-!      print*,'d_rei_dT ',d_rei_dt,rei_min,rei_max
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 c
 c Calculer l'epaisseur optique et l'emmissivite des nuages
-c
 c     IM inversion des DO
+c
       xflwp = 0.d0
       xfiwp = 0.d0
       xflwc = 0.d0
       xfiwc = 0.d0
-
-! Initialisation 
+c
       reliq=0.
       reice=0.
-
+c
       DO k = 1, klev
-         DO i = 1, klon
-            diff_paprs(i,k) = (paprs(i,k)-paprs(i,k+1))/RG
+        DO i = 1, klon
+c-layer calculation
+          rhodz(i,k) = (paprs(i,k)-paprs(i,k+1))/RG  ! kg/m2
+          zrho(i,k)=pplay(i,k)/t(i,k)/RD             ! kg/m3
+          dh(i,k)=rhodz(i,k)/zrho(i,k)               ! m
+c-Fraction of ice in cloud using a linear transition
+          zfice(i,k) = 1.0 - (t(i,k)-t_glace_min) / 
+     &                       (t_glace_max-t_glace_min)
+          zfice(i,k) = MIN(MAX(zfice(i,k),0.0),1.0)
+c-IM Total Liquid/Ice water content                                    
+          xflwc(i,k) = (1.-zfice(i,k))*pqlwp(i,k)
+          xfiwc(i,k) = zfice(i,k)*pqlwp(i,k)
          ENDDO
       ENDDO
 
-      IF (ok_newmicro) THEN
-
-
-         DO k = 1, klev
-            DO i = 1, klon
-c               zfice2(i,k) = 1.0 - (t(i,k)-t_glace) / (273.13-t_glace)
-               zfice2(i,k) = 1.0 - (t(i,k)-t_glace_min) / 
-     &                             (t_glace_max-t_glace_min)
-               zfice2(i,k) = MIN(MAX(zfice2(i,k),0.0),1.0)
-c     IM Total Liquid/Ice water content                                    
-               xflwc(i,k) = (1.-zfice2(i,k))*pqlwp(i,k)
-               xfiwc(i,k) = zfice2(i,k)*pqlwp(i,k)
-c     IM In-Cloud Liquid/Ice water content
-c     xflwc(i,k) = xflwc(i,k)+(1.-zfice)*pqlwp(i,k)/pclc(i,k)
-c     xfiwc(i,k) = xfiwc(i,k)+zfice*pqlwp(i,k)/pclc(i,k)
-            ENDDO
-         ENDDO
-
-         IF (ok_aie) THEN
-            DO k = 1, klev
-               DO i = 1, klon
-                                ! Formula "D" of Boucher and Lohmann, Tellus, 1995
-                                !             
-                  cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
+      IF (ok_cdnc) THEN
+c
+c--we compute cloud properties as a function of the aerosol load
+c
+        DO k = 1, klev
+            DO i = 1, klon
+c
+c Formula "D" of Boucher and Lohmann, Tellus, 1995
+c Cloud droplet number concentration (CDNC) is restricted
+c to be within [20, 1000 cm^3]
+c
+c--present-day case
+                cdnc(i,k) = 10.**(bl95_b0+bl95_b1*
      &               log(MAX(mass_solu_aero(i,k),1.e-4))/log(10.))*1.e6 !-m-3
-                                ! Cloud droplet number concentration (CDNC) is restricted
-                                ! to be within [20, 1000 cm^3]
-                                ! 
-                  cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
-                                !
-                                !
-                  cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
+                cdnc(i,k)=MIN(1000.e6,MAX(20.e6,cdnc(i,k)))
+c
+c--pre-industrial case
+                cdnc_pi(i,k) = 10.**(bl95_b0+bl95_b1*
      &               log(MAX(mass_solu_aero_pi(i,k),1.e-4))/log(10.))
      &               *1.e6 !-m-3
-                  cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
-               ENDDO
-            ENDDO
-            DO k = 1, klev
-               DO i = 1, klon
-!                  rad_chaud_tab(i,k) = 
-!     &                 MAX(1.1e6 
-!     &                 *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k)))  
-!     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.),5.)
-                  rad_chaud_tab(i,k) = 
-     &                 1.1
-     &                 *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k)))  
-     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.)
-                  rad_chaud_tab(i,k) = MAX(rad_chaud_tab(i,k) * 1e6, 5.) 
-               ENDDO            
-            ENDDO
-         ELSE
-            DO k = 1, MIN(3,klev)
-               DO i = 1, klon
-                  rad_chaud_tab(i,k) = rad_chau2
-               ENDDO            
-            ENDDO
-            DO k = MIN(3,klev)+1, klev
-               DO i = 1, klon
-                  rad_chaud_tab(i,k) = rad_chau1
-               ENDDO            
-            ENDDO
-
-         ENDIF
-         
-         DO k = 1, klev
-!            IF(.not.ok_aie) THEN
-            rad_chaud = rad_chau1
-            IF (k.LE.3) rad_chaud = rad_chau2
-!            ENDIF
-            DO i = 1, klon
-               IF (pclc(i,k) .LE. seuil_neb) THEN
-               
-c     -- effective cloud droplet radius (microns):
-               
-c     for liquid water clouds: 
-                                ! For output diagnostics
-                                !
-                                ! Cloud droplet effective radius [um]
-                                !
-                                ! we multiply here with f * xl (fraction of liquid water
-                                ! clouds in the grid cell) to avoid problems in the
-                                ! averaging of the output.
-                                ! In the output of IOIPSL, derive the real cloud droplet 
-                                ! effective radius as re/fl
-                                !
-                                    
-                  fl(i,k) = seuil_neb*(1.-zfice2(i,k))            
-                  re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
-                  
-                  rel = 0.
-                  rei = 0.
-                  pclc(i,k) = 0.0
-                  pcltau(i,k) = 0.0
-                  pclemi(i,k) = 0.0
-                  cldtaupi(i,k) = 0.0                  
-               ELSE
-
+                cdnc_pi(i,k)=MIN(1000.e6,MAX(20.e6,cdnc_pi(i,k)))
+c
+c--present-day case
+                rad_chaud(i,k) = 
+     &               1.1*((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k)))  
+     &               /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.)
+                rad_chaud(i,k) = MAX(rad_chaud(i,k) * 1.e6, 5.) 
+c
+c--pre-industrial case
+                radius = 
+     &               1.1*((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k)))  
+     &               /(4./3.*RPI*1000.*cdnc_pi(i,k)))**(1./3.)
+                radius = MAX(radius*1.e6, 5.) 
+c 
+c--pre-industrial case
+c--liquid/ice cloud water paths:
+                IF (pclc(i,k) .LE. seuil_neb) THEN
+c
+                pcldtaupi(i,k) = 0.0                  
+c
+                ELSE
+c                  
+                zflwp_var= 1000.*(1.-zfice(i,k))*pqlwp(i,k)/pclc(i,k)
+     &               *rhodz(i,k)
+                zfiwp_var= 1000.*zfice(i,k)*pqlwp(i,k)/pclc(i,k)
+     &               *rhodz(i,k)
+                tc = t(i,k)-273.15
+                rei = d_rei_dt*tc + rei_max
+                if (tc.le.-81.4) rei = rei_min
+c
+c-- cloud optical thickness :
+c [for liquid clouds, traditional formula, 
+c for ice clouds, Ebert & Curry (1992)] 
+c                  
+                if (zflwp_var.eq.0.) radius = 1. 
+                if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1. 
+                pcldtaupi(i,k) = 3.0/2.0 * zflwp_var / radius
+     &                 + zfiwp_var * (3.448e-03  + 2.431/rei)
+c
+                ENDIF
+c
+          ENDDO            
+        ENDDO
+c
+      ELSE !--not ok_cdnc
+c
+c-prescribed cloud droplet radius
+c
+         DO k = 1, MIN(3,klev)
+           DO i = 1, klon
+               rad_chaud(i,k) = rad_chau2
+           ENDDO            
+         ENDDO
+         DO k = MIN(3,klev)+1, klev
+           DO i = 1, klon
+               rad_chaud(i,k) = rad_chau1
+           ENDDO            
+         ENDDO
+
+      ENDIF !--ok_cdnc
+c
+c--computation of cloud optical depth and emissivity  
+c--in the general case
+c
+       DO k = 1, klev
+          DO i = 1, klon
+c
+             IF (pclc(i,k) .LE. seuil_neb) THEN
+c               
+c effective cloud droplet radius (microns) for liquid water clouds: 
+c For output diagnostics cloud droplet effective radius [um]
+c we multiply here with f * xl (fraction of liquid water
+c clouds in the grid cell) to avoid problems in the averaging of the output.
+c In the output of IOIPSL, derive the REAL cloud droplet 
+c effective radius as re/fl
+c
+                fl(i,k) = seuil_neb*(1.-zfice(i,k))            
+                re(i,k) = rad_chaud(i,k)*fl(i,k)
+                pclc(i,k)   = 0.0
+                pcltau(i,k) = 0.0
+                pclemi(i,k) = 0.0
+c
+             ELSE
+c
 c     -- liquid/ice cloud water paths:
                   
-                  zflwp_var= 1000.*(1.-zfice2(i,k))*pqlwp(i,k)/pclc(i,k)
-     &                 *diff_paprs(i,k)
-                  zfiwp_var= 1000.*zfice2(i,k)*pqlwp(i,k)/pclc(i,k)
-     &                 *diff_paprs(i,k)
-                  
-c     -- effective cloud droplet radius (microns):
-               
-c     for liquid water clouds: 
-                                    
-                  IF (ok_aie) THEN
-                     radius = 
-     &                    1.1
-     &                    *((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k)))  
-     &                    /(4./3.*RPI*1000.*cdnc_pi(i,k)))**(1./3.)
-                     radius = MAX(radius*1e6, 5.) 
-                  
-                     tc = t(i,k)-273.15
-                     rei = d_rei_dt*tc + rei_max
-                     if (tc.le.-81.4) rei = rei_min
-                     if (zflwp_var.eq.0.) radius = 1. 
-                     if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1. 
-                     cldtaupi(i,k) = 3.0/2.0 * zflwp_var / radius
-     &                    + zfiwp_var * (3.448e-03  + 2.431/rei)
-
-                  ENDIF         ! ok_aie
-                                ! For output diagnostics
-                                !
-                                ! Cloud droplet effective radius [um]
-                                !
-                                ! we multiply here with f * xl (fraction of liquid water
-                                ! clouds in the grid cell) to avoid problems in the
-                                ! averaging of the output.
-                                ! In the output of IOIPSL, derive the real cloud droplet 
-                                ! effective radius as re/fl
-                                !
- 
-                  fl(i,k) = pclc(i,k)*(1.-zfice2(i,k))            
-                  re(i,k) = rad_chaud_tab(i,k)*fl(i,k)
-                  
-                  rel = rad_chaud_tab(i,k)
-c     for ice clouds: as a function of the ambiant temperature
-c     [formula used by Iacobellis and Somerville (2000), with an 
-c     asymptotical value of 3.5 microns at T<-81.4 C added to be 
-c     consistent with observations of Heymsfield et al. 1986]:
-c  2011/05/24 : rei_min = 3.5 becomes a free parameter as well as rei_max=61.29
-                  tc = t(i,k)-273.15
-                  rei = d_rei_dt*tc + rei_max
-                  if (tc.le.-81.4) rei = rei_min
-c     -- cloud optical thickness :
-               
-c     [for liquid clouds, traditional formula, 
-c     for ice clouds, Ebert & Curry (1992)] 
-                  
+                zflwp_var= 1000.*(1.-zfice(i,k))*pqlwp(i,k)/pclc(i,k)
+     &               *rhodz(i,k)
+                zfiwp_var= 1000.*zfice(i,k)*pqlwp(i,k)/pclc(i,k)
+     &               *rhodz(i,k)
+c               
+c effective cloud droplet radius (microns) for liquid water clouds: 
+c For output diagnostics cloud droplet effective radius [um]
+c we multiply here with f * xl (fraction of liquid water
+c clouds in the grid cell) to avoid problems in the averaging of the output.
+c In the output of IOIPSL, derive the REAL cloud droplet 
+c effective radius as re/fl
+c 
+                fl(i,k) = pclc(i,k)*(1.-zfice(i,k))            
+                re(i,k) = rad_chaud(i,k)*fl(i,k)
+c 
+                rel = rad_chaud(i,k)
+c
+c for ice clouds: as a function of the ambiant temperature
+c [formula used by Iacobellis and Somerville (2000), with an 
+c asymptotical value of 3.5 microns at T<-81.4 C added to be 
+c consistent with observations of Heymsfield et al. 1986]:
+c 2011/05/24 : rei_min = 3.5 becomes a free PARAMETER as well as rei_max=61.29
+c
+                tc = t(i,k)-273.15
+                rei = d_rei_dt*tc + rei_max
+                if (tc.le.-81.4) rei = rei_min
+c
+c-- cloud optical thickness :
+c [for liquid clouds, traditional formula, 
+c for ice clouds, Ebert & Curry (1992)] 
+c                  
                  if (zflwp_var.eq.0.) rel = 1. 
                  if (zfiwp_var.eq.0. .or. rei.le.0.) rei = 1. 
                  pcltau(i,k) = 3.0/2.0 * ( zflwp_var/rel )
      &                 + zfiwp_var * (3.448e-03  + 2.431/rei)
+c
 c     -- cloud infrared emissivity:
-               
-c     [the broadband infrared absorption coefficient is parameterized
+c     [the broadband infrared absorption coefficient is PARAMETERized
 c     as a function of the effective cld droplet radius]
-               
 c     Ebert and Curry (1992) formula as used by Kiehl & Zender (1995):
+c
                   k_ice = k_ice0 + 1.0/rei
-                  
+c                 
                   pclemi(i,k) = 1.0
      &                 - EXP( -coef_chau*zflwp_var - DF*k_ice*zfiwp_var)
-
-               ENDIF
-              reliq(i,k)=rel
-              reice(i,k)=rei 
-!              if (i.eq.1) then
-!              print*,'Dans newmicro rel, rei :',rel, rei
-!              print*,'Dans newmicro reliq, reice :',
-!     $             reliq(i,k),reice(i,k)
-!              endif
-
-            ENDDO
-         ENDDO
-
+c
+             ENDIF
+c
+             reliq(i,k)=rel
+             reice(i,k)=rei 
+c
+             xflwp(i) = xflwp(i)+ xflwc(i,k) * rhodz(i,k)
+             xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * rhodz(i,k)
+c
+          ENDDO
+       ENDDO
+c
+c--if cloud droplet radius is fixed, then pcldtaupi=pcltau
+c
+      IF (.NOT.ok_cdnc) THEN
          DO k = 1, klev
             DO i = 1, klon
-               xflwp(i) = xflwp(i)+ xflwc(i,k) * diff_paprs(i,k)
-               xfiwp(i) = xfiwp(i)+ xfiwc(i,k) * diff_paprs(i,k)
-            ENDDO
-         ENDDO
-
-      ELSE
-         DO k = 1, klev
-            rad_chaud = rad_chau1
-            IF (k.LE.3) rad_chaud = rad_chau2
-            DO i = 1, klon
-                              
-               IF (pclc(i,k) .LE. seuil_neb) THEN
-
-                  pclc(i,k) = 0.0
-                  pcltau(i,k) = 0.0
-                  pclemi(i,k) = 0.0
-                  cldtaupi(i,k) = 0.0
-
-               ELSE
-
-                  zflwp_var = 1000.*pqlwp(i,k)*diff_paprs(i,k)
-     &                 /pclc(i,k)
-                  
-                  zfice1 = MIN(
-     &                 MAX( 1.0 - (t(i,k)-t_glace_min) / 
-     &                    (t_glace_max-t_glace_min),0.0),1.0)**nexpo
-                  
-                  radius = rad_chaud * (1.-zfice1) + rad_froid * zfice1
-                  coef   = coef_chau * (1.-zfice1) + coef_froi * zfice1
-
-                  pcltau(i,k) = 3.0 * zflwp_var / (2.0 * radius)
-                  pclemi(i,k) = 1.0 - EXP( - coef * zflwp_var)
-
-               ENDIF
-                              
-            ENDDO
-         ENDDO
-      ENDIF
-      
-      IF (.NOT.ok_aie) THEN
-         DO k = 1, klev
-            DO i = 1, klon
-               cldtaupi(i,k)=pcltau(i,k)
+               pcldtaupi(i,k)=pcltau(i,k)
             ENDDO
          ENDDO               
       ENDIF
-
-ccc   DO k = 1, klev
-ccc   DO i = 1, klon
-ccc   t(i,k) = t(i,k)
-ccc   pclc(i,k) = MAX( 1.e-5 , pclc(i,k) )
-ccc   lo = pclc(i,k) .GT. (2.*1.e-5)
-ccc   zflwp = pqlwp(i,k)*1000.*(paprs(i,k)-paprs(i,k+1))
-ccc   .          /(rg*pclc(i,k))
-ccc   zradef = 10.0 + (1.-sigs(k))*45.0
-ccc   pcltau(i,k) = 1.5 * zflwp / zradef
-ccc   zfice=1.0-MIN(MAX((t(i,k)-263.)/(273.-263.),0.0),1.0)
-ccc   zmsac = 0.13*(1.0-zfice) + 0.08*zfice
-ccc   pclemi(i,k) = 1.-EXP(-zmsac*zflwp)
-ccc   if (.NOT.lo) pclc(i,k) = 0.0
-ccc   if (.NOT.lo) pcltau(i,k) = 0.0
-ccc   if (.NOT.lo) pclemi(i,k) = 0.0
-ccc   ENDDO
-ccc   ENDDO
-ccccc print*, 'pas de nuage dans le rayonnement'
-ccccc DO k = 1, klev
-ccccc DO i = 1, klon
-ccccc pclc(i,k) = 0.0
-ccccc pcltau(i,k) = 0.0
-ccccc pclemi(i,k) = 0.0
-ccccc ENDDO
-ccccc ENDDO
 C     
 C     COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
-C     
 c     IM cf. CR:test: calcul prenant ou non en compte le recouvrement
 c     initialisations
+c
       DO i=1,klon
          zclear(i)=1.
@@ -465,30 +382,32 @@
       ENDDO
 C
-cIM cf CR DO k=1,klev
+c--calculation of liquid water path
+c
       DO k = klev, 1, -1
          DO i = 1, klon
-            pctlwp(i) = pctlwp(i) 
-     &           + pqlwp(i,k)*diff_paprs(i,k)
+            pctlwp(i) = pctlwp(i)+ pqlwp(i,k)*rhodz(i,k)
          ENDDO
       ENDDO
-c     IM cf. CR
+c     
+c--calculation of cloud properties with cloud overlap
+c
       IF (NOVLP.EQ.1) THEN
          DO k = klev, 1, -1
             DO i = 1, klon
                zclear(i)=zclear(i)*(1.-MAX(pclc(i,k),zcloud(i)))
-     &              /(1.-MIN(real(zcloud(i), kind=8),1.-ZEPSEC))
+     &              /(1.-MIN(REAL(zcloud(i), kind=8),1.-ZEPSEC))
                pct(i)=1.-zclear(i) 
                IF (paprs(i,k).LT.prmhc) THEN
                   pch(i) = pch(i)*(1.-MAX(pclc(i,k),zcloudh(i)))
-     &                 /(1.-MIN(real(zcloudh(i), kind=8),1.-ZEPSEC))
+     &                 /(1.-MIN(REAL(zcloudh(i), kind=8),1.-ZEPSEC))
                   zcloudh(i)=pclc(i,k)
                ELSE IF (paprs(i,k).GE.prmhc .AND.
      &                 paprs(i,k).LT.prlmc) THEN
                   pcm(i) = pcm(i)*(1.-MAX(pclc(i,k),zcloudm(i)))
-     &                 /(1.-MIN(real(zcloudm(i), kind=8),1.-ZEPSEC))
+     &                 /(1.-MIN(REAL(zcloudm(i), kind=8),1.-ZEPSEC))
                   zcloudm(i)=pclc(i,k)
                ELSE IF (paprs(i,k).GE.prlmc) THEN 
                   pcl(i) = pcl(i)*(1.-MAX(pclc(i,k),zcloudl(i)))
-     &                 /(1.-MIN(real(zcloudl(i), kind=8),1.-ZEPSEC))
+     &                 /(1.-MIN(REAL(zcloudl(i), kind=8),1.-ZEPSEC))
                   zcloudl(i)=pclc(i,k)
                endif
@@ -527,27 +446,26 @@
          ENDDO
       ENDIF
-      
 C     
       DO i = 1, klon
-c     IM cf. CR          pct(i)=1.-pct(i)
          pch(i)=1.-pch(i)
          pcm(i)=1.-pcm(i)
          pcl(i)=1.-pcl(i)
       ENDDO
-
+c
 c ========================================================
-! DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
+c DIAGNOSTICS CALCULATION FOR CMIP5 PROTOCOL
 c ========================================================
-!! change by Nicolas Yan (LSCE) 
-!! Cloud Droplet Number Concentration (CDNC) : 3D variable
-!! Fractionnal cover by liquid water cloud (LCC3D) : 3D variable
-!! Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable
-!! Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable
-!! Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable
-      IF (ok_newmicro) THEN
-         IF (ok_aie) THEN
+c change by Nicolas Yan (LSCE) 
+c Cloud Droplet Number Concentration (CDNC) : 3D variable
+c Fractionnal cover by liquid water cloud (LCC3D) : 3D variable
+c Cloud Droplet Number Concentration at top of cloud (CLDNCL) : 2D variable
+c Droplet effective radius at top of cloud (REFFCLWTOP) : 2D variable
+c Fractionnal cover by liquid water at top of clouds (LCC) : 2D variable
+c
+         IF (ok_cdnc) THEN
+c
             DO k = 1, klev
                DO i = 1, klon
-                  phase3d(i,k)=1-zfice2(i,k)
+                  phase3d(i,k)=1-zfice(i,k)
                   IF (pclc(i,k) .LE. seuil_neb) THEN
                      lcc3d(i,k)=seuil_neb*phase3d(i,k)
@@ -558,5 +476,5 @@
                ENDDO
             ENDDO
-
+c
             DO i=1,klon
                lcc(i)=0.
@@ -566,23 +484,16 @@
                IF(MAXIMUM) tcc(i) = 0.
             ENDDO
-     
-
+c
             DO i=1,klon
                DO k=klev-1,1,-1 !From TOA down
-
-
+c
             ! Test, if the cloud optical depth exceeds the necessary
             ! threshold:
 
-             IF (pcltau(i,k).GT.thres_tau .AND. pclc(i,k).GT.thres_neb)
-     .                                                             THEN
-               ! To calculate the right Temperature at cloud top,
-               ! interpolate it between layers:
-                  t_tmp = t(i,k) +
-     .              (paprs(i,k+1)-pplay(i,k))/(pplay(i,k+1)-pplay(i,k))
-     .              * ( t(i,k+1) - t(i,k) )
-
-                  IF(MAXIMUM) THEN
-                    IF(FIRST) THEN
+             IF (pcltau(i,k).GT.thres_tau 
+     .            .AND. pclc(i,k).GT.thres_neb) THEN
+
+                  IF (MAXIMUM) THEN
+                    IF (FIRST) THEN
                        write(*,*)'Hypothese de recouvrement: MAXIMUM'
                        FIRST=.FALSE.
@@ -592,6 +503,6 @@
                   ENDIF
 
-                  IF(RANDOM) THEN
-                    IF(FIRST) THEN
+                  IF (RANDOM) THEN
+                    IF (FIRST) THEN
                        write(*,*)'Hypothese de recouvrement: RANDOM'
                        FIRST=.FALSE.
@@ -601,6 +512,6 @@
                   ENDIF
 
-                  IF(MAXIMUM_RANDOM) THEN
-                    IF(FIRST) THEN
+                  IF (MAXIMUM_RANDOM) THEN
+                    IF (FIRST) THEN
                        write(*,*)'Hypothese de recouvrement: MAXIMUM_
      .                         RANDOM'
@@ -613,5 +524,5 @@
                   ENDIF
 c Effective radius of cloud droplet at top of cloud (m)
-                  reffclwtop(i) = reffclwtop(i) + rad_chaud_tab(i,k) * 
+                  reffclwtop(i) = reffclwtop(i) + rad_chaud(i,k) * 
      .           1.0E-06 * phase3d(i,k) * ( tcc(i) - ftmp(i))*flag_max
 c CDNC at top of cloud (m-3)
@@ -626,6 +537,7 @@
           ENDIF ! is there a visible, not-too-small cloud?  
           ENDDO ! loop over k
-
-          IF(RANDOM .OR. MAXIMUM_RANDOM) tcc(i)=1.-tcc(i)
+c
+          IF (RANDOM .OR. MAXIMUM_RANDOM) tcc(i)=1.-tcc(i)
+c
          ENDDO ! loop over i
 
@@ -633,48 +545,23 @@
             DO i = 1, klon
                DO k = 1, klev
-                  pqlwpcon(i,k)=rnebcon(i,k)*clwcon(i,k) ! fraction eau liquide convective
-                  pqlwpstra(i,k)=pclc(i,k)*phase3d(i,k)-pqlwpcon(i,k) ! fraction eau liquide stratiforme
-                  IF (pqlwpstra(i,k) .LE. 0.0) pqlwpstra(i,k)=0.0
+! Weight to be used for outputs: eau_liquide*couverture nuageuse
+                  lcc3dcon(i,k) =rnebcon(i,k)*phase3d(i,k)*clwcon(i,k)  ! eau liquide convective
+                  lcc3dstra(i,k)=pclc(i,k)*pqlwp(i,k)*phase3d(i,k)
+                  lcc3dstra(i,k)=lcc3dstra(i,k)-lcc3dcon(i,k)           ! eau liquide stratiforme
+                  lcc3dstra(i,k)=MAX(lcc3dstra(i,k),0.0)
+! Compute cloud droplet radius as above in meter
+                  radius=1.1*((pqlwp(i,k)*pplay(i,k)/(RD * T(i,k)))  
+     &               /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.)
+                  radius=MAX(radius, 5.e-6) 
 ! Convective Cloud Droplet Effective Radius (REFFCLWC) : variable 3D
-                  reffclwc(i,k)=1.1
-     &                 *((pqlwpcon(i,k)*pplay(i,k)/(RD * T(i,k)))
-     &                 /(4./3*RPI*1000.*cdnc(i,k)) )**(1./3.)
-                  reffclwc(i,k) = MAX(reffclwc(i,k) * 1e6, 5.)
-
+                  reffclwc(i,k)=radius
+                  reffclwc(i,k)=reffclwc(i,k)*lcc3dcon(i,k)
 ! Stratiform Cloud Droplet Effective Radius (REFFCLWS) : variable 3D
-                  IF ((pclc(i,k)-rnebcon(i,k)) .LE. seuil_neb) THEN ! tout sous la forme convective
-                     reffclws(i,k)=0.0
-                     lcc3dstra(i,k)= 0.0
-                  ELSE
-                     reffclws(i,k) = (pclc(i,k)*phase3d(i,k)*
-     &                               rad_chaud_tab(i,k)-
-     &                            pqlwpcon(i,k)*reffclwc(i,k))
-                     IF(reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0
-                     lcc3dstra(i,k)=pqlwpstra(i,k)
-                 ENDIF
-!Convertion from um to m
-                  IF(rnebcon(i,k). LE. seuil_neb) THEN
-                    reffclwc(i,k) = reffclwc(i,k)*seuil_neb*clwcon(i,k)
-     &                              *1.0E-06
-                    lcc3dcon(i,k)= seuil_neb*clwcon(i,k)
-                  ELSE
-                    reffclwc(i,k) = reffclwc(i,k)*pqlwpcon(i,k)
-     &                              *1.0E-06
-                    lcc3dcon(i,k) = pqlwpcon(i,k)
-                  ENDIF
-
-                  reffclws(i,k) = reffclws(i,k)*1.0E-06
-
+                  reffclws(i,k)=radius
+                  reffclws(i,k)=reffclws(i,k)*lcc3dstra(i,k)
                ENDDO !klev
             ENDDO !klon 
-
-!! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
-            DO k = 1, klev
-               DO i = 1, klon
-                   pdel(i,k) = paprs(i,k)-paprs(i,k+1)
-                   zrho(i,k)=pplay(i,k)/t(i,k)/RD                  ! kg/m3
-                   dh(i,k)=pdel(i,k)/(gravit*zrho(i,k)) ! hauteur de chaque boite (m)
-               ENDDO
-            ENDDO
+c
+c Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
 c
             DO i = 1, klon
@@ -697,21 +584,20 @@
             DO i = 1, klon
                DO k = 1, klev
-                IF (scdnc(i,k) .LE. 0.0) scdnc(i,k)=0.0
-                IF (reffclws(i,k) .LE. 0.0) reffclws(i,k)=0.0
-                IF (reffclwc(i,k) .LE. 0.0) reffclwc(i,k)=0.0
-                IF (lcc3d(i,k) .LE. 0.0) lcc3d(i,k)=0.0
-                IF (lcc3dcon(i,k) .LE. 0.0) lcc3dcon(i,k)=0.0
+                IF (scdnc(i,k) .LE. 0.0)     scdnc(i,k)=0.0
+                IF (reffclws(i,k) .LE. 0.0)  reffclws(i,k)=0.0
+                IF (reffclwc(i,k) .LE. 0.0)  reffclwc(i,k)=0.0
+                IF (lcc3d(i,k) .LE. 0.0)     lcc3d(i,k)=0.0
+                IF (lcc3dcon(i,k) .LE. 0.0)  lcc3dcon(i,k)=0.0
                 IF (lcc3dstra(i,k) .LE. 0.0) lcc3dstra(i,k)=0.0
                ENDDO
-               IF (reffclwtop(i) .LE. 0.0) reffclwtop(i)=0.0
-               IF (cldncl(i) .LE. 0.0) cldncl(i)=0.0
-               IF (cldnvi(i) .LE. 0.0) cldnvi(i)=0.0
-               IF (lcc(i) .LE. 0.0) lcc(i)=0.0
+               IF (reffclwtop(i) .LE. 0.0)  reffclwtop(i)=0.0
+               IF (cldncl(i) .LE. 0.0)      cldncl(i)=0.0
+               IF (cldnvi(i) .LE. 0.0)      cldnvi(i)=0.0
+               IF (lcc(i) .LE. 0.0)         lcc(i)=0.0
             ENDDO
-
-         ENDIF !ok_aie
-      ENDIF !ok newmicro
-c
-C
+c
+      ENDIF !ok_cdnc
+c
       RETURN
+c
       END
Index: /LMDZ5/trunk/libf/phylmd/physiq.F
===================================================================
--- /LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1711)
+++ /LMDZ5/trunk/libf/phylmd/physiq.F	(revision 1712)
@@ -1094,7 +1094,8 @@
       ! Parameters
       LOGICAL ok_ade, ok_aie    ! Apply aerosol (in)direct effects or not
+      LOGICAL ok_cdnc          ! ok cloud droplet number concentration (O. Boucher 01-2013)
       REAL bl95_b0, bl95_b1   ! Parameter in Boucher and Lohmann (1995)
-      SAVE ok_ade, ok_aie, bl95_b0, bl95_b1
-c$OMP THREADPRIVATE(ok_ade, ok_aie, bl95_b0, bl95_b1)
+      SAVE ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1
+c$OMP THREADPRIVATE(ok_ade, ok_aie, ok_cdnc, bl95_b0, bl95_b1)
       LOGICAL, SAVE :: aerosol_couple ! true  : calcul des aerosols dans INCA
                                       ! false : lecture des aerosol dans un fichier
@@ -1251,5 +1252,5 @@
      .     fact_cldcon, facttemps,ok_newmicro,iflag_radia,
      .     iflag_cldcon,iflag_ratqs,ratqsbas,ratqshaut,tau_ratqs,
-     .     ok_ade, ok_aie, aerosol_couple, 
+     .     ok_ade, ok_aie, ok_cdnc, aerosol_couple, 
      .     flag_aerosol, new_aod,
      .     bl95_b0, bl95_b1,
@@ -3145,12 +3146,10 @@
 
       if (ok_newmicro) then
-      CALL newmicro (paprs, pplay,ok_newmicro,
-     .            t_seri, cldliq, cldfra, cldtau, cldemi,
-     .            cldh, cldl, cldm, cldt, cldq,
-     .            flwp, fiwp, flwc, fiwc,
-     e            ok_aie,
-     e            mass_solu_aero, mass_solu_aero_pi,
-     e            bl95_b0, bl95_b1,
-     s            cldtaupi, re, fl, ref_liq, ref_ice)
+      CALL newmicro (ok_cdnc, bl95_b0, bl95_b1,
+     .              paprs, pplay, t_seri, cldliq, cldfra,
+     .              cldtau, cldemi, cldh, cldl, cldm, cldt, cldq,
+     e              flwp, fiwp, flwc, fiwc,
+     e              mass_solu_aero, mass_solu_aero_pi,
+     s              cldtaupi, re, fl, ref_liq, ref_ice) 
       else
       CALL nuage (paprs, pplay,
@@ -3161,5 +3160,4 @@
      e            bl95_b0, bl95_b1,
      s            cldtaupi, re, fl)
-      
       endif
 c
