Index: trunk/LMDZ.MARS/libf/phymars/co2cloud.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/co2cloud.F	(revision 1818)
+++ trunk/LMDZ.MARS/libf/phymars/co2cloud.F	(revision 1820)
@@ -6,9 +6,9 @@
      &                rsedcloud,rhocloud,pzlev,pdqs_sedco2,
      &                pdu,pu)
-! to use  'getin'
+      USE ioipsl_getincom, only: getin
       use dimradmars_mod, only: naerkind
-      USE comcstfi_h
-      USE ioipsl_getincom
-      USE updaterad
+      USE comcstfi_h, only: pi, g, cpp
+      USE updaterad, only: updaterice_microco2, updaterice_micro,
+     &                     updaterdust
       use conc_mod, only: mmean,rnew
       use tracer_mod, only: nqmx, igcm_co2, igcm_co2_ice,
@@ -21,7 +21,7 @@
       IMPLICIT NONE
 
-#include "datafile.h"
-#include "callkeys.h"
-#include "microphys.h"
+      include "datafile.h"
+      include "callkeys.h"
+      include "microphys.h"
 
 c=======================================================================
@@ -68,59 +68,44 @@
 
 c-----------------------------------------------------------------------
-c   declarations:
+c   arguments:
 c   -------------
 
-c   Inputs:
-c   ------
-
-      INTEGER ngrid,nlay
-      INTEGER nq                 ! nombre de traceurs 
-      REAL ptimestep            ! pas de temps physique (s)
-      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
-      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
-      REAL pdpsrf(ngrid)         ! tendence surf pressure
-      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
-      REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
-      REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
-      real,intent(in) :: pzlev(ngrid,nlay+1) ! altitude at the boundaries of the layers
-      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
-      real pdq(ngrid,nlay,nq)    ! tendance avant condensation  (kg/kg.s-1)
-
-      real rice(ngrid,nlay)    ! Water Ice mass mean radius (m)
+      INTEGER, INTENT(IN) :: ngrid,nlay
+      REAL, INTENT(IN) :: ptimestep            ! pas de temps physique (s)
+      REAL, INTENT(IN) :: pplev(ngrid,nlay+1)   ! Inter-layer pressures (Pa)
+      REAL, INTENT(IN) :: pplay(ngrid,nlay)     ! mid-layer pressures (Pa)
+      REAL, INTENT(IN) :: pdpsrf(ngrid)         ! tendency on surface pressure
+      REAL, INTENT(IN) :: pzlay(ngrid,nlay)     ! altitude at the middle of the layers
+      REAL, INTENT(IN) :: pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
+      REAL, INTENT(IN) :: pdt(ngrid,nlay)       ! tendency on temperature from other parametrizations
+      real, INTENT(IN) :: pq(ngrid,nlay,nq)     ! tracers (kg/kg)
+      real, INTENT(IN) :: pdq(ngrid,nlay,nq)    ! tendencies before condensation  (kg/kg.s-1)
+      real, intent(out) :: pdqcloudco2(ngrid,nlay,nq) ! tendency due to CO2 condensation (kg/kg.s-1)
+      real, intent(out) :: pdtcloudco2(ngrid,nlay)    ! tendency on temperature due to latent heat
+      INTEGER, INTENT(IN) :: nq                 ! number of tracers
+      REAL, INTENT(IN) :: tau(ngrid,naerkind) ! Column dust optical depth at each point
+      REAL, INTENT(IN) :: tauscaling(ngrid)   ! Convertion factor for dust amount
+      REAL, INTENT(OUT) :: rdust(ngrid,nlay)   ! Dust geometric mean radius (m)
+      real, intent(OUT) :: rice(ngrid,nlay)    ! Water Ice mass mean radius (m)
                                 ! used for nucleation of CO2 on ice-coated ccns
-      DOUBLE PRECISION rho_ice_co2T(ngrid,nlay) !T-dependant CO2 ice density
-      DOUBLE PRECISION :: myT   ! temperature scalar for co2 density computation
-      REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
-      REAL tauscaling(ngrid)   ! Convertion factor for dust amount
-      real rdust(ngrid,nlay)   ! Dust geometric mean radius (m)
-
-c   Outputs:
-c   -------
-
-      real pdqcloudco2(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
-      REAL pdtcloudco2(ngrid,nlay)    ! tendence temperature due
-                                   ! a la chaleur latente
-      DOUBLE PRECISION riceco2(ngrid,nlay)    ! Ice mass mean radius (m)
+      DOUBLE PRECISION, INTENT(out) :: riceco2(ngrid,nlay)    ! Ice mass mean radius (m)
                                ! (r_c in montmessin_2004)
-      REAL nuice(ngrid,nlay)   ! Estimated effective variance
+      REAL, INTENT(in) :: nuice(ngrid,nlay)   ! Estimated effective variance
                                !   of the size distribution
-      real rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius
-      real rhocloudco2(ngrid,nlay)  ! Cloud density (kg.m-3)
-      real rhocloudco2t(ngrid,nlay)  ! Cloud density (kg.m-3)
-      real pdqs_sedco2(ngrid) ! CO2 flux at the surface
+      real, intent(out) :: rsedcloudco2(ngrid,nlay) ! Cloud sedimentation radius
+      real, intent(out) :: rhocloudco2(ngrid,nlay)  ! Cloud density (kg.m-3)
+      real, intent(out) :: rsedcloud(ngrid,nlay) ! Water Cloud sedimentation radius
+      real, intent(out) :: rhocloud(ngrid,nlay)  ! Water Cloud density (kg.m-3)
+      real, intent(in) :: pzlev(ngrid,nlay+1) ! altitude at the boundaries of the layers
+      real, intent(out) :: pdqs_sedco2(ngrid) ! CO2 flux at the surface
+      REAL, INTENT(IN) :: pdu(ngrid,nlay),pu(ngrid,nlay) !Zonal Wind: zu=pu+pdu*ptimestep
 
 c   local:
 c   ------
-      !water
-      real rsedcloud(ngrid,nlay) ! Cloud sedimentation radius
-      real rhocloud(ngrid,nlay)  ! Cloud density (kg.m-3)
-      ! for ice radius computation
         
       ! for time loop
       INTEGER microstep  ! time subsampling step variable
-      INTEGER imicroco2     ! time subsampling for coupled water microphysics & sedimentation
-      SAVE imicroco2
-      REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
-      SAVE microtimestep
+      INTEGER, SAVE :: imicroco2     ! time subsampling for coupled water microphysics & sedimentation
+      REAL, SAVE :: microtimestep ! integration timestep for coupled water microphysics & sedimentation
       
       ! tendency given by clouds (inside the micro loop)
@@ -135,4 +120,7 @@
       REAL satuco2(ngrid,nlay)  ! co2 satu ratio for output
       REAL zqsatco2(ngrid,nlay) ! saturation co2
+
+      DOUBLE PRECISION rho_ice_co2T(ngrid,nlay) !T-dependant CO2 ice density
+      DOUBLE PRECISION :: myT   ! temperature scalar for co2 density computation
 
       INTEGER iq,ig,l,i
@@ -154,6 +142,5 @@
 !     What we need for Qext reading and tau computation : size distribution
       DOUBLE PRECISION vrat_cld ! Volume ratio
-      DOUBLE PRECISION rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
-      SAVE rb_cldco2
+      DOUBLE PRECISION, SAVE :: rb_cldco2(nbinco2_cld+1) ! boundary values of each rad_cldco2 bin (m)
       DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9 ! Minimum radius (m)
       DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6 ! Maximum radius (m)
@@ -162,11 +149,11 @@
       DOUBLE PRECISION dr_cld(nbinco2_cld) ! width of each rad_cldco2 bin (m)
       DOUBLE PRECISION vol_cld(nbinco2_cld) ! particle volume for each bin (m3)
-      REAL sigma_iceco2 ! Variance of the ice and CCN distributions
+      REAL, SAVE :: sigma_iceco2 ! Variance of the ice and CCN distributions
       logical :: file_ok !Qext file reading
       double precision :: radv(10000),Qextv1mic(10000)
-      double precision :: Qext1bins(100),Qtemp
+      double precision, save :: Qext1bins(100)
+      double precision :: Qtemp
       double precision :: ltemp1(10000),ltemp2(10000)
       integer :: nelem,lebon1,lebon2,uQext
-      save Qext1bins
       DOUBLE PRECISION n_aer(nbinco2_cld),Rn,No,n_derf,dev2
       DOUBLE PRECISION Qext1bins2(ngrid,nlay)   
@@ -177,4 +164,6 @@
       REAL zt(ngrid,nlay)       ! local value of temperature
       REAL :: zq(ngrid, nlay,nq)
+
+      real :: rhocloudco2t(ngrid,nlay)  ! Cloud density (kg.m-3)
 
       DOUBLE PRECISION :: tcond(ngrid,nlay) !CO2 condensation temperature
@@ -187,5 +176,4 @@
       REAL ::  mincloud ! min cloud frac
       DOUBLE PRECISION:: rho,zu,NN,gradT !For Saturation Index computation
-      REAL :: pdu(ngrid,nlay),pu(ngrid,nlay) !Wind field zu=pu+pdu*ptimestep
       DOUBLE PRECISION :: SatIndex(ngrid,nlay),SatIndexmap(ngrid)
       
@@ -378,31 +366,32 @@
 
          if (satindexco2) then !logical in callphys.def
-         DO l=12,26
-            ! layers 12 --> 26 ~ 12->85 km
-            DO ig=1,ngrid
-         ! Calcul de N^2 static stability
-            gradT=(zt(ig,l+1)-zt(ig,l))/(pzlev(ig,l+1)-pzlev(ig,l))
-            NN=sqrt(g/zt(iq,l)*(max(gradT,-g/cpp)+g/cpp))
-            !calcul of wind field 
-            zu=pu(ig,l)  + pdu(ig,l)*ptimestep
-!     calcul of background density
-            rho=pplay(ig,l)/(rnew(ig,l)*zt(ig,l))
-            !saturation index
-            SatIndex(ig,l)=sqrt(7.5e-7*150.e3/(2.*pi)*NN/(rho*zu*zu*zu))
-         enddo
-         enddo
+           DO l=12,26
+             ! layers 12 --> 26 ~ 12->85 km
+             DO ig=1,ngrid
+             ! compute N^2 static stability
+             gradT=(zt(ig,l+1)-zt(ig,l))/(pzlev(ig,l+1)-pzlev(ig,l))
+             NN=sqrt(g/zt(iq,l)*(max(gradT,-g/cpp)+g/cpp))
+             ! compute absolute value of zonal wind field 
+             zu=abs(pu(ig,l)  + pdu(ig,l)*ptimestep)
+             ! compute background density
+             rho=pplay(ig,l)/(rnew(ig,l)*zt(ig,l))
+             !saturation index
+             SatIndex(ig,l)=sqrt(7.5e-7*150.e3/(2.*pi)*NN/
+     &                                               (rho*zu*zu*zu))
+             ENDDO
+           ENDDO
             !Then compute Satindex map 
             ! layers 12 --> 26 ~ 12->85 km
-         DO ig=1,ngrid
-            SatIndexmap(ig)=maxval(SatIndex(ig,12:26))
-         enddo
-
-         call WRITEDIAGFI(ngrid,"SatIndexmap","SatIndexmap","km",2,
-     &        SatIndexmap)
+           DO ig=1,ngrid
+             SatIndexmap(ig)=maxval(SatIndex(ig,12:26))
+           ENDDO
+
+           call WRITEDIAGFI(ngrid,"SatIndexmap","SatIndexmap","km",2,
+     &         SatIndexmap)
          else
-            DO ig=1,ngrid
-               SatIndexmap(ig)=0.05 !maxval(SatIndex(ig,12:26))
-            enddo
-         endif
+           do ig=1,ngrid
+             SatIndexmap(ig)=0.05 !maxval(SatIndex(ig,12:26))
+           enddo
+         endif ! of if (satindexco2)
 
 !Modulate the DeltaT by GW propagation index : 
@@ -414,16 +403,15 @@
         
          CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond)
-! A tester:            CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond)
          zdelt=spant  
-           DO ig=1,ngrid
+         DO ig=1,ngrid
               
-            if (SatIndexmap(ig) .le. 0.1) THEN
-              DO l=1,nlay-1
+           IF (SatIndexmap(ig) .le. 0.1) THEN
+             DO l=1,nlay-1
          
                IF (tcond(ig,l) .ge. (zt(ig,l)+zdelt) 
-     &             .or. tcond(ig,l) .le. 0 ) THEN !Toute la fraction est saturée
+     &             .or. tcond(ig,l) .le. 0 ) THEN !The entire fraction is saturated
                   zteff(ig,l)=zt(ig,l)
                   cloudfrac(ig,l)=1.
-               ELSE IF (tcond(ig,l) .le. (zt(ig,l)-zdelt)) THEN !Rien n'est saturé
+               ELSE IF (tcond(ig,l) .le. (zt(ig,l)-zdelt)) THEN ! No saturation at all
                   zteff(ig,l)=zt(ig,l)-zdelt
                   cloudfrac(ig,l)=mincloud
@@ -431,5 +419,5 @@
                   cloudfrac(ig,l)=(tcond(ig,l)-zt(ig,l)+zdelt)/
      &                 (2.0*zdelt)
-                  zteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2. !Temperature moyenne de la fraction nuageuse 
+                  zteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2. !Mean temperature of the cloud fraction
                END IF           !ig if (tcond(ig,l) ...
                zteff(ig,l)=zteff(ig,l)-pdt(ig,l)*ptimestep
@@ -439,11 +427,11 @@
                   cloudfrac(ig,l)=1.
                END IF
-            ENDDO
-         ELSE
+             ENDDO
+           ELSE
 !SatIndex not favorable for GW : leave pt untouched
-            zteff(ig,l)=pt(ig,l)
-            cloudfrac(ig,l)=mincloud
-         END IF                 ! of if(SatIndexmap...
-      ENDDO
+             zteff(ig,l)=pt(ig,l)
+             cloudfrac(ig,l)=mincloud
+           ENDIF                 ! of if(SatIndexmap...
+         ENDDO ! of DO ig=1,ngrid
 !     Totalcloud frac of the column missing here
 c-----------------------
@@ -466,6 +454,6 @@
 c------ Temperature tendency subpdt
         ! If imicro=1 subpdt is the same as pdt
-         DO l=1,nlay
-            DO ig=1,ngrid
+        DO l=1,nlay
+          DO ig=1,ngrid
                subpdt(ig,l) = subpdt(ig,l)
      &              + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
@@ -500,8 +488,8 @@
      &              subpdq(ig,l,igcm_ccn_number)
      &              + pdq(ig,l,igcm_ccn_number)
-            ENDDO
-         ENDDO
+          ENDDO
+        ENDDO
 c- Effective tracers quantities in the cloud fraction
-         IF (CLFvaryingCO2) THEN     
+        IF (CLFvaryingCO2) THEN     
             pqeff(:,:,:)=pq(:,:,:) ! prevent from buggs (A. Pottier)
             pqeff(:,:,igcm_ccnco2_mass) =pq(:,:,igcm_ccnco2_mass)/
@@ -511,7 +499,7 @@
             pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice)/
      &           cloudfrac(:,:)
-         ELSE
+        ELSE
             pqeff(:,:,:)=pq(:,:,:)
-         END IF  
+        END IF  
        
 c------------------------------------------------------
@@ -519,5 +507,5 @@
 c   call to sedimentation routine, update tendancies
 c------------------------------------------------------
-       DO l=1, nlay
+        DO l=1, nlay
           DO ig=1,ngrid             
              tempo_traceur_t(ig,l)=zteff(ig,l)+subpdt(ig,l)
@@ -548,30 +536,30 @@
      &            riceco2(ig,l))
           ENDDO
-       ENDDO
+        ENDDO
 !     Gravitational sedimentation       
-       sav_trac(:,:,igcm_co2_ice)=tempo_traceurs(:,:,igcm_co2_ice)
-       sav_trac(:,:,igcm_ccnco2_mass)=
+        sav_trac(:,:,igcm_co2_ice)=tempo_traceurs(:,:,igcm_co2_ice)
+        sav_trac(:,:,igcm_ccnco2_mass)=
      &      tempo_traceurs(:,:,igcm_ccnco2_mass)
-       sav_trac(:,:,igcm_ccnco2_number)=
+        sav_trac(:,:,igcm_ccnco2_number)=
      &      tempo_traceurs(:,:,igcm_ccnco2_number) 
        !We save actualized tracer values to compute sedimentation tendancies
-       call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
+        call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
      &     microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
      &     rsedcloudco2,rhocloudco2t,
      &     tempo_traceurs(:,:,igcm_co2_ice),wq,beta) !  3 traceurs
 !     sedim at the surface of co2 ice : keep track of it for physiq_mod
-      do ig=1,ngrid 
-         pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1)/microtimestep
-      end do
-      call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
+        do ig=1,ngrid 
+          pdqs_sedco2(ig)=pdqs_sedco2(ig)+ wq(ig,1)/microtimestep
+        end do
+        call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
      &     microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
      &     rsedcloudco2,rhocloudco2t,
      &     tempo_traceurs(:,:,igcm_ccnco2_mass),wq,beta) 
-      call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
+        call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
      &     microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
      &     rsedcloudco2,rhocloudco2t,
      &     tempo_traceurs(:,:,igcm_ccnco2_number),wq,beta) 
-      DO l = 1, nlay            !Compute tendencies
-         DO ig=1,ngrid
+        DO l = 1, nlay            !Compute tendencies
+          DO ig=1,ngrid
             pdqsed(ig,l,igcm_ccnco2_mass)=
      &           (tempo_traceurs(ig,l,igcm_ccnco2_mass)-
@@ -583,6 +571,6 @@
      &           (tempo_traceurs(ig,l,igcm_co2_ice)-
      &           sav_trac(ig,l,igcm_co2_ice))/microtimestep
-         ENDDO
-      ENDDO
+          ENDDO
+        ENDDO
       !update subtimestep tendencies with sedimentation input
         DO l=1,nlay
@@ -598,9 +586,9 @@
      &           +pdqsed(ig,l,igcm_co2_ice)
          ENDDO
-      ENDDO   
+        ENDDO   
 c------------------------------------------------------
 c      2.  Main call to the cloud schemes:
 c------------------------------------------------------
-      CALL improvedCO2clouds(ngrid,nlay,microtimestep,
+        CALL improvedCO2clouds(ngrid,nlay,microtimestep,
      &     pplay,pplev,zteff,subpdt,
      &     pqeff,subpdq,subpdqcloudco2,subpdtcloudco2,
@@ -609,6 +597,6 @@
 c      3.  Updating tendencies after cloud scheme:
 c-----------------------------------------------------
-          DO l=1,nlay
-            DO ig=1,ngrid
+        DO l=1,nlay
+          DO ig=1,ngrid
                subpdt(ig,l) =
      &              subpdt(ig,l) + subpdtcloudco2(ig,l)
@@ -644,6 +632,6 @@
      &              subpdq(ig,l,igcm_ccn_number)
      &              + subpdqcloudco2(ig,l,igcm_ccn_number)
-            ENDDO
-         ENDDO
+          ENDDO
+        ENDDO
       ENDDO                     ! of DO microstep=1,imicro
       
@@ -656,13 +644,13 @@
       enddo
 c------ Temperature tendency pdtcloud
-       DO l=1,nlay
-         DO ig=1,ngrid
+      DO l=1,nlay
+        DO ig=1,ngrid
              pdtcloudco2(ig,l) =
      &         subpdt(ig,l)/real(imicroco2)-pdt(ig,l)
-          ENDDO
-       ENDDO
+        ENDDO
+      ENDDO
 c------ Tracers tendencies pdqcloud
-       DO l=1,nlay
-          DO ig=1,ngrid        
+      DO l=1,nlay
+        DO ig=1,ngrid        
              pdqcloudco2(ig,l,igcm_co2_ice) = 
      &            subpdq(ig,l,igcm_co2_ice)/real(imicroco2)
@@ -674,8 +662,8 @@
      &            subpdq(ig,l,igcm_h2o_ice)/real(imicroco2)
      &            - pdq(ig,l,igcm_h2o_ice)
-          ENDDO
-       ENDDO
-       DO l=1,nlay
-          DO ig=1,ngrid
+        ENDDO
+      ENDDO
+      DO l=1,nlay
+        DO ig=1,ngrid
              pdqcloudco2(ig,l,igcm_ccnco2_mass) = 
      &            subpdq(ig,l,igcm_ccnco2_mass)/real(imicroco2)
@@ -690,8 +678,8 @@
      &            subpdq(ig,l,igcm_ccn_number)/real(imicroco2)
      &            - pdq(ig,l,igcm_ccn_number)
-          ENDDO
-       ENDDO
-       DO l=1,nlay
-          DO ig=1,ngrid
+        ENDDO
+      ENDDO
+      DO l=1,nlay
+        DO ig=1,ngrid
              pdqcloudco2(ig,l,igcm_dust_mass) = 
      &            subpdq(ig,l,igcm_dust_mass)/real(imicroco2)
@@ -700,10 +688,10 @@
      &            subpdq(ig,l,igcm_dust_number)/real(imicroco2)
      &            - pdq(ig,l,igcm_dust_number)
-          ENDDO
-       ENDDO
+        ENDDO
+      ENDDO
 c-------Due to stepped entry, other processes tendencies can add up to negative values
 c-------Therefore, enforce positive values and conserve mass
-       DO l=1,nlay
-          DO ig=1,ngrid
+      DO l=1,nlay
+        DO ig=1,ngrid
              IF ((pqeff(ig,l,igcm_ccnco2_number) + 
      &            ptimestep* (pdq(ig,l,igcm_ccnco2_number) + 
@@ -725,8 +713,8 @@
      &               -pdqcloudco2(ig,l,igcm_ccnco2_mass)
              ENDIF
-          ENDDO
-       ENDDO
-       DO l=1,nlay
-          DO ig=1,ngrid
+        ENDDO
+      ENDDO
+      DO l=1,nlay
+        DO ig=1,ngrid
              IF ( (pqeff(ig,l,igcm_dust_number) + 
      &            ptimestep* (pdq(ig,l,igcm_dust_number) + 
@@ -747,9 +735,9 @@
      &               -pdqcloudco2(ig,l,igcm_dust_mass)
              ENDIF
-          ENDDO
-       ENDDO
+        ENDDO
+      ENDDO
         !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq      
-       DO l=1,nlay
-          DO ig=1,ngrid
+      DO l=1,nlay
+        DO ig=1,ngrid
              IF (pqeff(ig,l,igcm_co2_ice) + ptimestep*
      &       (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice)) 
@@ -766,10 +754,10 @@
            pdqcloudco2(ig,l,igcm_co2_ice)= -pdqcloudco2(ig,l,igcm_co2)
           ENDIF
-         ENDDO
-        ENDDO
+        ENDDO
+      ENDDO
 
 c Update clouds parameters values in the cloud fraction (for output)
-        DO l=1, nlay
-           DO ig=1,ngrid
+      DO l=1, nlay
+        DO ig=1,ngrid
 
               Niceco2=pqeff(ig,l,igcm_co2_ice) +                   
@@ -801,19 +789,21 @@
      &             Nccnco2*tauscaling(ig) .le. 1.) )THEN
                  riceco2(ig,l)=0.
+                 Qext1bins2(ig,l)=0.
+              else
+c     Compute opacities
+                No=Nccnco2*tauscaling(ig)
+                Rn=-dlog(riceco2(ig,l))
+                n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
+                Qext1bins2(ig,l)=0.
+                do i = 1, nbinco2_cld 
+                 n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
+                 n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
+                 n_aer(i) = n_aer(i) + 0.5 * No * n_derf
+                 Qext1bins2(ig,l)=Qext1bins2(ig,l)+Qext1bins(i)*n_aer(i)
+                enddo
               endif
-c     Compute opacities
-           No=Nccnco2*tauscaling(ig)
-           Rn=-dlog(riceco2(ig,l))
-           n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
-           Qext1bins2(ig,l)=0.
-           do i = 1, nbinco2_cld 
-              n_aer(i) = -0.5 * No * n_derf !! this ith previously computed
-              n_derf = derf((rb_cldco2(i+1)+Rn) *dev2)
-              n_aer(i) = n_aer(i) + 0.5 * No * n_derf
-              Qext1bins2(ig,l)=Qext1bins2(ig,l)+Qext1bins(i)*n_aer(i)
-           enddo
      
       !update rice water
-        call updaterice_micro(
+          call updaterice_micro(
      &    pqeff(ig,l,igcm_h2o_ice) +                    ! ice mass
      &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
@@ -827,5 +817,5 @@
      &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
           
-        call updaterdust(
+          call updaterdust(
      &    pqeff(ig,l,igcm_dust_mass) +                   ! dust mass
      &   (pdq(ig,l,igcm_dust_mass) +                  ! dust mass
@@ -836,6 +826,6 @@
      &    rdust(ig,l))
 
-         ENDDO
-       ENDDO 
+        ENDDO
+      ENDDO 
      
 c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
@@ -843,10 +833,10 @@
 c     Then that should not affect the ice particle radius
        
-       do ig=1,ngrid
-          if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
+      do ig=1,ngrid
+        if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,2)))then
              if(pdpsrf(ig)*ptimestep.gt.0.9*(pplev(ig,1)-pplev(ig,3)))
      &       riceco2(ig,2)=riceco2(ig,3) 
              riceco2(ig,1)=riceco2(ig,2)
-        end if
+        endif
       end do
        
@@ -860,5 +850,5 @@
       ENDDO
        
-       DO l=1,nlay
+      DO l=1,nlay
          DO ig=1,ngrid
            rsedcloudco2(ig,l)=max(riceco2(ig,l)*
@@ -867,18 +857,18 @@
 c          rsedcloudco2(ig,l)=min(rsedcloudco2(ig,l),1.e-5)
          ENDDO
-       ENDDO
+      ENDDO
        
-       call co2sat(ngrid*nlay,zteff+(pdt+pdtcloudco2)*ptimestep
+      call co2sat(ngrid*nlay,zteff+(pdt+pdtcloudco2)*ptimestep
      &      ,pplay,zqsatco2)
-       do l=1,nlay
-          do ig=1,ngrid 
+      do l=1,nlay
+        do ig=1,ngrid 
              satuco2(ig,l) = (pqeff(ig,l,igcm_co2)  + 
      &            (pdq(ig,l,igcm_co2) +
      &            pdqcloudco2(ig,l,igcm_co2))*ptimestep)*
      &            (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l)
-            enddo
-         enddo
-!Tout ce qui est modifié par la microphysique de CO2 doit être rapporté a cloudfrac 
-         IF (CLFvaryingCO2) THEN
+        enddo
+      enddo
+!Everything modified by CO2 microphysics must be wrt cloudfrac
+      IF (CLFvaryingCO2) THEN
             DO l=1,nlay
                DO ig=1,ngrid
@@ -905,27 +895,27 @@
                ENDDO
             ENDDO   
-         ENDIF
-!l'opacité de la case ig est la somme sur l de Qext1bins2: Est-ce-vrai?
-         tau1mic(:)=0.
-         do l=1,nlay
-            do ig=1,ngrid 
-               tau1mic(ig)=tau1mic(ig)+Qext1bins2(ig,l)
-            enddo
-         enddo
+      ENDIF
+! opacity in mesh ig is the sum over l of Qext1bins2: Is this true ?
+      tau1mic(:)=0.
+      do l=1,nlay
+        do ig=1,ngrid 
+          tau1mic(ig)=tau1mic(ig)+Qext1bins2(ig,l)
+        enddo
+      enddo
 !Outputs:
-         call WRITEDIAGFI(ngrid,"SatIndex","SatIndex"," ",3,
+      call WRITEDIAGFI(ngrid,"SatIndex","SatIndex"," ",3,
      &        SatIndex)
-         call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3,
+      call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3,
      &        satuco2)
-         call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
+      call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
      &        ,3,riceco2)         
-         call WRITEdiagfi(ngrid,"cloudfrac","co2 cloud fraction"
+      call WRITEdiagfi(ngrid,"cloudfrac","co2 cloud fraction"
      &        ," ",3,cloudfrac)
-         call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2"
+      call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2"
      &        ,"m",3,rsedcloudco2)
-         call WRITEdiagfi(ngrid,"Tau3D1mic"," co2 ice opacities"
+      call WRITEdiagfi(ngrid,"Tau3D1mic"," co2 ice opacities"
      &        ," ",3,Qext1bins2)
-         call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron"
+      call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron"
      &        ," ",2,tau1mic)
          
-         END
+      END
