Index: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/newmicro.F90
===================================================================
--- LMDZ6/branches/DYNAMICO-conv/libf/phylmd/newmicro.F90	(revision 3024)
+++ LMDZ6/branches/DYNAMICO-conv/libf/phylmd/newmicro.F90	(revision 3356)
@@ -1,7 +1,5 @@
 ! $Id$
 
-
-
-SUBROUTINE newmicro(ok_cdnc, bl95_b0, bl95_b1, paprs, pplay, t, pqlwp, pclc, &
+SUBROUTINE newmicro(flag_aerosol, ok_cdnc, bl95_b0, bl95_b1, paprs, pplay, t, pqlwp, pclc, &
     pcltau, pclemi, pch, pcl, pcm, pct, pctlwp, xflwp, xfiwp, xflwc, xfiwc, &
     mass_solu_aero, mass_solu_aero_pi, pcldtaupi, re, fl, reliq, reice, &
@@ -10,7 +8,12 @@
   USE dimphy
   USE phys_local_var_mod, ONLY: scdnc, cldncl, reffclwtop, lcc, reffclws, &
-    reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra
+      reffclwc, cldnvi, lcc3d, lcc3dcon, lcc3dstra, icc3dcon, icc3dstra,  & 
+      zfice, dNovrN
   USE phys_state_var_mod, ONLY: rnebcon, clwcon
   USE icefrac_lsc_mod ! computes ice fraction (JBM 3/14)
+  USE ioipsl_getin_p_mod, ONLY : getin_p
+  USE print_control_mod, ONLY: lunout
+
+
   IMPLICIT NONE
   ! ======================================================================
@@ -139,4 +142,5 @@
   ! within the grid cell)
 
+  INTEGER flag_aerosol
   LOGICAL ok_cdnc
   REAL bl95_b0, bl95_b1 ! Parameter in B&L 95-Formula
@@ -152,5 +156,4 @@
   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 rad_chaud_pi(klon, klev) !--rayon pour les nuages chauds pre-industriels
@@ -162,4 +165,11 @@
   REAL reliq_pi(klon, klev), reice_pi(klon, klev)
 
+  REAL,SAVE :: cdnc_min=-1.
+  REAL,SAVE :: cdnc_min_m3
+  !$OMP THREADPRIVATE(cdnc_min,cdnc_min_m3)
+  REAL,SAVE :: cdnc_max=-1.
+  REAL,SAVE :: cdnc_max_m3
+  !$OMP THREADPRIVATE(cdnc_max,cdnc_max_m3)
+
   ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! FH : 2011/05/24
@@ -173,4 +183,15 @@
   ! Pour retrouver les résultats numériques de la version d'origine,
   ! on impose 0.71 quand on est proche de 0.71
+
+  if (first) THEN
+      call getin_p('cdnc_min',cdnc_min)
+      cdnc_min_m3=cdnc_min*1.E6
+      IF (cdnc_min_m3<0.) cdnc_min_m3=20.E6 ! astuce pour retrocompatibilite
+      write(lunout,*)'cdnc_min=', cdnc_min_m3/1.E6
+      call getin_p('cdnc_max',cdnc_max)
+      cdnc_max_m3=cdnc_max*1.E6
+      IF (cdnc_max_m3<0.) cdnc_max_m3=1000.E6 ! astuce pour retrocompatibilite
+      write(lunout,*)'cdnc_max=', cdnc_max_m3/1.E6
+  ENDIF
 
   d_rei_dt = (rei_max-rei_min)/81.4
@@ -204,6 +225,6 @@
         xflwc(i, k) = (1.-zfice(i,k))*pqlwp(i, k)
         xfiwc(i, k) = zfice(i, k)*pqlwp(i, k)
-      END DO
-    END DO
+      ENDDO
+    ENDDO
   ELSE ! of IF (iflag_t_glace.EQ.0)
     DO k = 1, klev
@@ -222,6 +243,6 @@
         xflwc(i, k) = (1.-zfice(i,k))*pqlwp(i, k)
         xfiwc(i, k) = zfice(i, k)*pqlwp(i, k)
-      END DO
-    END DO
+      ENDDO
+    ENDDO
   ENDIF
 
@@ -232,18 +253,50 @@
     DO k = 1, klev
       DO i = 1, klon
-
         ! Formula "D" of Boucher and Lohmann, Tellus, 1995
         ! Cloud droplet number concentration (CDNC) is restricted
         ! to be within [20, 1000 cm^3]
 
-        ! --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
-        cdnc(i, k) = min(1000.E6, max(20.E6,cdnc(i,k)))
-
         ! --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)))
+        cdnc_pi(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc_pi(i,k)))
+
+      ENDDO
+    ENDDO
+
+    !--flag_aerosol=7 => MACv2SP climatology 
+    !--in this case there is an enhancement factor
+    IF (flag_aerosol .EQ. 7) THEN 
+
+      !--present-day
+      DO k = 1, klev
+        DO i = 1, klon
+          cdnc(i, k) = cdnc_pi(i,k)*dNovrN(i)
+        ENDDO
+      ENDDO
+
+    !--standard case
+    ELSE
+
+      DO k = 1, klev
+        DO i = 1, klon
+
+          ! Formula "D" of Boucher and Lohmann, Tellus, 1995
+          ! Cloud droplet number concentration (CDNC) is restricted
+          ! to be within [20, 1000 cm^3]
+
+          ! --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
+          cdnc(i, k) = min(cdnc_max_m3, max(cdnc_min_m3,cdnc(i,k)))
+
+        ENDDO
+      ENDDO
+
+    ENDIF !--flag_aerosol
+
+    !--computing cloud droplet size
+    DO k = 1, klev
+      DO i = 1, klon
 
         ! --present-day case
@@ -280,8 +333,8 @@
             zfiwp_var*(3.448E-03+2.431/rei)
 
-        END IF
-
-      END DO
-    END DO
+        ENDIF
+
+      ENDDO
+    ENDDO
 
   ELSE !--not ok_cdnc
@@ -293,14 +346,14 @@
         rad_chaud(i, k) = rad_chau2
         rad_chaud_pi(i, k) = rad_chau2
-      END DO
-    END DO
+      ENDDO
+    ENDDO
     DO k = min(3, klev) + 1, klev
       DO i = 1, klon
         rad_chaud(i, k) = rad_chau1
         rad_chaud_pi(i, k) = rad_chau1
-      END DO
-    END DO
-
-  END IF !--ok_cdnc
+      ENDDO
+    ENDDO
+
+  ENDIF !--ok_cdnc
 
   ! --computation of cloud optical depth and emissivity
@@ -377,5 +430,5 @@
         pclemi(i, k) = 1.0 - exp(-coef_chau*zflwp_var-df*k_ice*zfiwp_var)
 
-      END IF
+      ENDIF
 
       reice(i, k) = rei
@@ -384,6 +437,6 @@
       xfiwp(i) = xfiwp(i) + xfiwc(i, k)*rhodz(i, k)
 
-    END DO
-  END DO
+    ENDDO
+  ENDDO
 
   ! --if cloud droplet radius is fixed, then pcldtaupi=pcltau
@@ -394,7 +447,7 @@
         pcldtaupi(i, k) = pcltau(i, k)
         reice_pi(i, k) = reice(i, k)
-      END DO
-    END DO
-  END IF
+      ENDDO
+    ENDDO
+  ENDIF
 
   DO k = 1, klev
@@ -403,6 +456,6 @@
       reliq_pi(i, k) = rad_chaud_pi(i, k)
       reice_pi(i, k) = reice(i, k)
-    END DO
-  END DO
+    ENDDO
+  ENDDO
 
   ! COMPUTE CLOUD LIQUID PATH AND TOTAL CLOUDINESS
@@ -420,5 +473,5 @@
     pcl(i) = 1.0
     pctlwp(i) = 0.0
-  END DO
+  ENDDO
 
   ! --calculation of liquid water path
@@ -427,6 +480,6 @@
     DO i = 1, klon
       pctlwp(i) = pctlwp(i) + pqlwp(i, k)*rhodz(i, k)
-    END DO
-  END DO
+    ENDDO
+  ENDDO
 
   ! --calculation of cloud properties with cloud overlap
@@ -450,8 +503,8 @@
             (i),kind=8),1.-zepsec))
           zcloudl(i) = pclc(i, k)
-        END IF
+        ENDIF
         zcloud(i) = pclc(i, k)
-      END DO
-    END DO
+      ENDDO
+    ENDDO
   ELSE IF (novlp==2) THEN
     DO k = klev, 1, -1
@@ -465,7 +518,7 @@
         ELSE IF (paprs(i,k)>=prlmc) THEN
           pcl(i) = min(pclc(i,k), pcl(i))
-        END IF
-      END DO
-    END DO
+        ENDIF
+      ENDDO
+    ENDDO
   ELSE IF (novlp==3) THEN
     DO k = klev, 1, -1
@@ -479,8 +532,8 @@
         ELSE IF (paprs(i,k)>=prlmc) THEN
           pcl(i) = pcl(i)*(1.0-pclc(i,k))
-        END IF
-      END DO
-    END DO
-  END IF
+        ENDIF
+      ENDDO
+    ENDDO
+  ENDIF
 
   DO i = 1, klon
@@ -488,5 +541,5 @@
     pcm(i) = 1. - pcm(i)
     pcl(i) = 1. - pcl(i)
-  END DO
+  ENDDO
 
   ! ========================================================
@@ -509,8 +562,8 @@
         ELSE
           lcc3d(i, k) = pclc(i, k)*phase3d(i, k)
-        END IF
+        ENDIF
         scdnc(i, k) = lcc3d(i, k)*cdnc(i, k) ! m-3
-      END DO
-    END DO
+      ENDDO
+    ENDDO
 
     DO i = 1, klon
@@ -520,5 +573,5 @@
       IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1.
       IF (novlp.EQ.2) tcc(i) = 0.
-    END DO
+    ENDDO
 
     DO i = 1, klon
@@ -534,8 +587,8 @@
               WRITE (*, *) 'Hypothese de recouvrement: MAXIMUM'
               first = .FALSE.
-            END IF
+            ENDIF
             flag_max = -1.
             ftmp(i) = max(tcc(i), pclc(i,k))
-          END IF
+          ENDIF
 
           IF (novlp.EQ.3) THEN
@@ -543,8 +596,8 @@
               WRITE (*, *) 'Hypothese de recouvrement: RANDOM'
               first = .FALSE.
-            END IF
+            ENDIF
             flag_max = 1.
             ftmp(i) = tcc(i)*(1-pclc(i,k))
-          END IF
+          ENDIF
 
           IF (novlp.EQ.1) THEN
@@ -554,9 +607,9 @@
                 &                                          RANDOM'
               first = .FALSE.
-            END IF
+            ENDIF
             flag_max = 1.
             ftmp(i) = tcc(i)*(1.-max(pclc(i,k),pclc(i,k+1)))/(1.-min(pclc(i, &
               k+1),1.-thres_neb))
-          END IF
+          ENDIF
           ! Effective radius of cloud droplet at top of cloud (m)
           reffclwtop(i) = reffclwtop(i) + rad_chaud(i, k)*1.0E-06*phase3d(i, &
@@ -570,10 +623,10 @@
           tcc(i) = ftmp(i)
 
-        END IF ! is there a visible, not-too-small cloud?
-      END DO ! loop over k
+        ENDIF ! is there a visible, not-too-small cloud?
+      ENDDO ! loop over k
 
       IF (novlp.EQ.3 .OR. novlp.EQ.1) tcc(i) = 1. - tcc(i)
 
-    END DO ! loop over i
+    ENDDO ! loop over i
 
     ! ! Convective and Stratiform Cloud Droplet Effective Radius (REFFCLWC
@@ -586,4 +639,11 @@
         lcc3dstra(i, k) = lcc3dstra(i, k) - lcc3dcon(i, k) ! eau liquide stratiforme
         lcc3dstra(i, k) = max(lcc3dstra(i,k), 0.0)
+        !FC pour la glace (CAUSES)
+        icc3dcon(i, k) = rnebcon(i, k)*(1-phase3d(i, k))*clwcon(i, k) !  glace convective
+        icc3dstra(i, k)= pclc(i, k)*pqlwp(i, k)*(1-phase3d(i, k))
+        icc3dstra(i, k) = icc3dstra(i, k) - icc3dcon(i, k) ! glace stratiforme
+        icc3dstra(i, k) = max( icc3dstra(i, k), 0.0)
+        !FC (CAUSES)
+
         ! 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.* &
@@ -596,6 +656,6 @@
         reffclws(i, k) = radius
         reffclws(i, k) = reffclws(i, k)*lcc3dstra(i, k)
-      END DO !klev
-    END DO !klon
+      ENDDO !klev
+    ENDDO !klon
 
     ! Column Integrated Cloud Droplet Number (CLDNVI) : variable 2D
@@ -609,5 +669,5 @@
         lcc_integrat(i) = lcc_integrat(i) + lcc3d(i, k)*dh(i, k)
         height(i) = height(i) + dh(i, k)
-      END DO ! klev
+      ENDDO ! klev
       lcc_integrat(i) = lcc_integrat(i)/height(i)
       IF (lcc_integrat(i)<=1.0E-03) THEN
@@ -615,6 +675,6 @@
       ELSE
         cldnvi(i) = cldnvi(i)*lcc(i)/lcc_integrat(i)
-      END IF
-    END DO ! klon
+      ENDIF
+    ENDDO ! klon
 
     DO i = 1, klon
@@ -626,12 +686,18 @@
         IF (lcc3dcon(i,k)<=0.0) lcc3dcon(i, k) = 0.0
         IF (lcc3dstra(i,k)<=0.0) lcc3dstra(i, k) = 0.0
-      END DO
+!FC (CAUSES)
+        IF (icc3dcon(i,k)<=0.0) icc3dcon(i, k) = 0.0
+        IF (icc3dstra(i,k)<=0.0) icc3dstra(i, k) = 0.0
+!FC (CAUSES)
+      ENDDO
       IF (reffclwtop(i)<=0.0) reffclwtop(i) = 0.0
       IF (cldncl(i)<=0.0) cldncl(i) = 0.0
       IF (cldnvi(i)<=0.0) cldnvi(i) = 0.0
       IF (lcc(i)<=0.0) lcc(i) = 0.0
-    END DO
-
-  END IF !ok_cdnc
+    ENDDO
+
+  ENDIF !ok_cdnc
+
+  first=.false. !to be sure
 
   RETURN
