Changeset 1816

Timestamp:

Nov 2, 2017, 4:40:47 PM (7 years ago)

Author:

jaudouard

Message:

Commit for CO2 clouds microphysics.

Location:

trunk/LMDZ.MARS

Files:

• README (modified) (1 diff)
• libf/phymars/co2cloud.F (modified) (31 diffs)
• libf/phymars/improvedCO2clouds.F (modified) (28 diffs)
• libf/phymars/massflowrateCO2.F (modified) (34 diffs)
• libf/phymars/microphys.h (modified) (2 diffs)
• libf/phymars/nucleaCO2.F (modified) (1 diff)
• libf/phymars/physiq_mod.F (modified) (8 diffs)
• libf/phymars/tcondco2.F90 (added)
• libf/phymars/updaterad.F90 (modified) (4 diffs)

trunk/LMDZ.MARS/README

@@ -2470 +2470 @@
 - Some cleaning of the CO2 ice clouds routine has been done. Not perfect yet!
  
+== 02/10/2017 == JA
+- Update for CO2 clouds scheme.
+- Several bugs have been corrected, further cleaning performed.
+- The main routine of co2 clouds is co2cloud.F. Please read its comments if interested.
+
+ 

trunk/LMDZ.MARS/libf/phymars/co2cloud.F

@@ -4 +4 @@
      &                nq,tau,tauscaling,rdust,rice,riceco2,nuice,
      &                rsedcloudco2,rhocloudco2,
-     &                rsedcloud,rhocloud,pzlev,pdqs_sedco2)
+     &                rsedcloud,rhocloud,pzlev,pdqs_sedco2,
+     &                pdu,pu)
 ! to use  'getin'
       use dimradmars_mod, only: naerkind
@@ -10 +11 @@
       USE ioipsl_getincom
       USE updaterad
-      use conc_mod, only: mmean
+      use conc_mod, only: mmean,rnew
       use tracer_mod, only: nqmx, igcm_co2, igcm_co2_ice,
      &     igcm_dust_mass, igcm_dust_number,igcm_h2o_ice,
@@ -18 +19 @@
      &     rho_ice_co2,r3n_q,rho_ice,nuice_sed
      
-
       IMPLICIT NONE
-
 
 #include "datafile.h"
 #include "callkeys.h"
 #include "microphys.h"
-
 
 c=======================================================================
@@ -33 +31 @@
 c  due to timescales smaller than the GCM integration timestep.
 c  microphysics subroutine is improvedCO2clouds.F
-c  
+c  the microphysics time step is a fraction of the physical one
+c  the integer imicroco2 must be set in callphys.def  
+c
 c  The co2 clouds tracers (co2_ice, ccn mass and concentration) are 
 c  sedimented at each microtimestep. pdqs_sedco2 keeps track of the 
@@ -43 +43 @@
 c
 c 07/2017 J.Audouard
-c Several logicals can be set to .true. in callphys.def
+c Several logicals and integer must be set to .true. in callphys.def 
+c uf not, default values are .false.
 c co2clouds=.true. call this routine
 c co2useh2o=.true. allow the use of water ice particles as CCN for CO2
@@ -49 +50 @@
 c CLFvaryingCO2=.true. allows a subgrid temperature distribution
 c of amplitude spantCO2(=integer in callphys.def)
+c imicroco2=50 
 c
+c The subgrid Temperature distribution is modulated (0 or 1) by Spiga et
+c al. (GRL 2012) Saturation Index to account for GW propagation or 
+c dissipation upwards.
+c
+c 4D and column opacities are computed using Qext values at 1µm. 
 c=======================================================================
 
@@ -69 +76 @@
       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)
@@ -75 +81 @@
       real 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)
+      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
@@ -86 +93 @@
       REAL pdtcloudco2(ngrid,nlay)    ! tendence temperature due
                                    ! a la chaleur latente
-
       DOUBLE PRECISION riceco2(ngrid,nlay)    ! Ice mass mean radius (m)
                                ! (r_c in montmessin_2004)
@@ -94 +100 @@
       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 pdqs_sedco2(ngrid) ! CO2 flux at the surface
+
 c   local:
 c   ------
@@ -101 +108 @@
       real rhocloud(ngrid,nlay)  ! Cloud density (kg.m-3)
       ! for ice radius computation
-      
-      REAl ccntyp
-      
+        
       ! for time loop
       INTEGER microstep  ! time subsampling step variable
-      INTEGER imicro     ! time subsampling for coupled water microphysics & sedimentation
-      SAVE imicro
+      INTEGER imicroco2     ! time subsampling for coupled water microphysics & sedimentation
+      SAVE imicroco2
       REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
       SAVE microtimestep
@@ -125 +130 @@
       INTEGER iq,ig,l,i
       LOGICAL,SAVE :: firstcall=.true.
-      DOUBLE PRECISION Nccnco2, Niceco2,Nco2,mdustJA,ndustJA
-      DOUBLE PRECISION Qccnco2
-      real :: beta
+      DOUBLE PRECISION Nccnco2, Niceco2,Nco2,Qccnco2
+      real :: beta ! for sedimentation
 
       real epaisseur (ngrid,nlay) ! Layer thickness (m)
       real masse (ngrid,nlay) ! Layer mass (kg.m-2)
-      double precision diff,diff0
-      real tempo_traceur_t(ngrid,nlay)
+      real tempo_traceur_t(ngrid,nlay) ! tracers with real-time value in microtimeloop
       real tempo_traceurs(ngrid,nlay,nq)
-      real sav_trac(ngrid,nlay,nq)
+      real sav_trac(ngrid,nlay,nq) !For sedimentation tendancy
       real pdqsed(ngrid,nlay,nq)
-      REAL lw                   !Latent heat of sublimation (J.kg-1) 
-      REAL,save :: l0,l1,l2,l3,l4 !Coeffs for lw
-
-      DOUBLE PRECISION,allocatable,save :: memdMMccn(:,:) ! Nb particules intégré
-      DOUBLE PRECISION,allocatable,save :: memdMMh2o(:,:)
-      DOUBLE PRECISION,allocatable,save :: memdNNccn(:,:)
-      DOUBLE PRECISION :: myT
-!     What we need for Qext reading and tau computation
+
+      DOUBLE PRECISION,allocatable,save :: memdMMccn(:,:) !memory of h2o particles 
+      DOUBLE PRECISION,allocatable,save :: memdMMh2o(:,:) !only if co2useh2o=.true.
+      DOUBLE PRECISION,allocatable,save :: memdNNccn(:,:) !Nb particules H2O intégré 
+     
+!     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, PARAMETER :: rmin_cld = 1.e-11 ! Minimum radius (m)
-      DOUBLE PRECISION, PARAMETER :: rmax_cld = 3.e-6 ! Maximum radius (m)
-      DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-12
-! Minimum boundary radius (m)
-      DOUBLE PRECISION, PARAMETER :: rbmax_cld = 5.e-6 ! Maximum boundary radius (m)
-
+      DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9 ! Minimum radius (m)
+      DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6 ! Maximum radius (m)
+      DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-10! Minimum boundary radius (m)
+      DOUBLE PRECISION, PARAMETER :: rbmax_cld = 2.e-4 ! Maximum boundary radius (m)
       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
-
-      logical :: file_ok
+      logical :: file_ok !Qext file reading
       double precision :: radv(10000),Qextv1mic(10000)
       double precision :: Qext1bins(100),Qtemp
@@ -163 +161 @@
       integer :: nelem,lebon1,lebon2,uQext
       save Qext1bins
-      character(len=100) scanline  
       DOUBLE PRECISION n_aer(nbinco2_cld),Rn,No,n_derf,dev2
       DOUBLE PRECISION Qext1bins2(ngrid,nlay)   
@@ -173 +170 @@
       REAL :: zq(ngrid, nlay,nq)
 
-      REAL :: tcond(ngrid,nlay)
+      DOUBLE PRECISION :: tcond(ngrid,nlay) !CO2 condensation temperature
       REAL ::  zqvap(ngrid,nlay)
-      REAL :: zqice(ngrid,nlay)
+      REAL ::  zqice(ngrid,nlay)
       REAL ::  spant,zdelt ! delta T for the temperature distribution
-      REAL :: zteff(ngrid, nlay)! effective temperature in the cloud,neb
-      REAL :: pqeff(ngrid, nlay, nq)! effective tracers quantities in the cloud
-      REAL :: cloudfrac(ngrid,nlay) ! cloud fraction
-      REAL :: mincloud ! min cloud frac
+      REAL ::  zteff(ngrid, nlay)! effective temperature in the cloud,neb
+      REAL ::  pqeff(ngrid, nlay, nq)! effective tracers quantities in the cloud
+      REAL ::  cloudfrac(ngrid,nlay) ! cloud fraction
+      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)
+      
 c      logical :: CLFvaryingCO2
 c     ** un petit test de coherence
@@ -197 +198 @@
         write(*,*) "time subsampling for microphysic ?"
 #ifdef MESOSCALE
-        imicro = 2
+        imicroco2 = 2
 #else
-        imicro = 30
+        imicroco2 = 30
 #endif
-        call getin("imicro",imicro)
-        write(*,*)"imicro = ",imicro
+        call getin("imicroco2",imicroco2)
+        write(*,*)"imicroco2 = ",imicroco2
         
-        microtimestep = ptimestep/real(imicro)
+        microtimestep = ptimestep/real(imicroco2)
         write(*,*)"Physical timestep is",ptimestep 
         write(*,*)"CO2 Microphysics timestep is",microtimestep 
       
-        ! Values for latent heat computation
-        l0=595594d0      
-        l1=903.111d0    
-        l2=-11.5959d0   
-        l3=0.0528288d0
-        l4=-0.000103183d0 
-c$$$        
-c$$$        if (meteo_flux_number .ne. 0) then 
-c$$$
-c$$$        write(*,*) "Warning ! Meteoritic flux of dust is turned on"
-c$$$        write(*,*) "Dust mass = ",meteo_flux_mass 
-c$$$        write(*,*) "Dust number = ",meteo_flux_number
-c$$$        write(*,*) "Are added at the z-level (km)",meteo_alt
-c$$$        write(*,*) "Every timestep in co2cloud.F"
-c$$$         endif
-         
+   
         if (.not. allocated(memdMMccn)) allocate(memdMMccn(ngrid,nlay))
         if (.not. allocated(memdNNccn)) allocate(memdNNccn(ngrid,nlay))
@@ -289 +275 @@
            ltemp2=abs(radv(:)-rb_cldco2(i+1))
            lebon1=minloc(ltemp1,DIM=1)
-           lebon2=minloc(ltemp2,DIM=1)
+           lebon2=min(minloc(ltemp2,DIM=1),10000)
            nelem=lebon2-lebon1+1.
            Qtemp=0d0
            do l=0,nelem
-              Qtemp=Qtemp+Qextv1mic(lebon1+l) !mean value in the interval
+              Qtemp=Qtemp+Qextv1mic(min(lebon1+l,10000)) !mean value in the interval
            enddo
            Qtemp=Qtemp/nelem
@@ -320 +306 @@
           write(*,*) "for the CO2 microphysics "
         endif
+
         firstcall=.false.
       ENDIF                     ! of IF (firstcall)
@@ -330 +317 @@
       subpdqcloudco2(1:ngrid,1:nlay,1:nq) = 0
       subpdtcloudco2(1:ngrid,1:nlay)      = 0
-c$$$
-c$$$       call WRITEDIAGFI(ngrid,"pzlay","pzlay","km",3,
-c$$$     &        pzlay)
-c$$$       call WRITEDIAGFI(ngrid,"pplay","pplay","Pa",3,
-c$$$     &        pplay)
 
       wq(:,:)=0
@@ -353 +335 @@
       enddo
 
-      !CLFvaryingCO2=.true.
-      
 c-------------------------------------------------------------------
 c   0.  Representation of sub-grid water ice clouds
-c------------------
+c-------------------------------------------------------------------
       IF (CLFvaryingCO2) THEN
+
          spant=spantCO2         ! delta T for the temprature distribution
          mincloud=0.1           ! min cloudfrac when there is ice  
-        ! IF (flagcloudco2) THEN
-        !    write(*,*) "CLFCO2 Delta T", spant
-        !    write(*,*) "CLFCO2 mincloud", mincloud
-        !    flagcloudco2=.false.
-        ! END IF
+         zteff(:,:)=pt(:,:) 
+         cloudfrac(:,:)=mincloud
          
 c-----Tendencies
@@ -382 +360 @@
          zqvap=zq(:,:,igcm_co2)
          zqice=zq(:,:,igcm_co2_ice)
-!! AS : this routine is not present? 
-!         CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond)
+
+
+         call WRITEDIAGFI(ngrid,"pzlev","pzlev","km",3,
+     &        pzlev)
+         call WRITEDIAGFI(ngrid,"pzlay","pzlay","km",3,
+     &        pzlay)
+         call WRITEDIAGFI(ngrid,"pplay","pplay","Pa",3,
+     &        pplay)
+
+         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
+            !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)
+
+!Modulate the DeltaT by GW propagation index : 
+         ! Saturation index S in Spiga 2012 paper 
+         !Assuming like in the paper, 
+         !GW phase speed (stationary waves) c=0 m.s-1
+         !lambdaH =150 km
+         !Fo=7.5e-7 J.m-3
+        
+         CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond)
 ! A tester:            CALL tcondco2(ngrid,nlay,pplay,zqvap,tcond)
          zdelt=spant  
-         DO l=1,nlay
-            DO ig=1,ngrid
-               IF (tcond(ig,l) .ge. (zt(ig,l)+zdelt)) THEN !Toute la fraction est saturée
+           DO ig=1,ngrid
+              
+            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
                   zteff(ig,l)=zt(ig,l)
                   cloudfrac(ig,l)=1.
@@ -397 +417 @@
                   cloudfrac(ig,l)=(tcond(ig,l)-zt(ig,l)+zdelt)/
      &                 (2.0*zdelt)
-                  zteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2. !Ja Temperature moyenne de la fraction nuageuse 
-                  
-               END IF
+                  zteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2. !Temperature moyenne de la fraction nuageuse 
+               END IF           !ig if (tcond(ig,l) ...
                zteff(ig,l)=zteff(ig,l)-pdt(ig,l)*ptimestep
                IF (cloudfrac(ig,l).le. mincloud) THEN
                   cloudfrac(ig,l)=mincloud
-               ELSE IF (cloudfrac(ig,l).ge. 1) THEN
+               ELSE IF (cloudfrac(ig,l).gt. 1) THEN
                   cloudfrac(ig,l)=1.
                END IF
             ENDDO
-         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
 !     Totalcloud frac of the column missing here
 c-----------------------
@@ -417 +441 @@
             END DO
          END DO
-      END IF                    ! end if (CLFvaryingco2)--------------------------------------------
-c   1.  Tendencies: 
-c------------------
-  
-c------ Effective tracers quantities in the cloud fraction
-        IF (CLFvaryingCO2) THEN     
-           pqeff(:,:,:)=pq(:,:,:) ! prevent from buggs (A. Pottier)
-   
-c           pqeff(:,:,igcm_ccn_mass) =pq(:,:,igcm_ccn_mass)/
-c     &                                cloudfrac(:,:)
-c           pqeff(:,:,igcm_ccn_number)= 
-c     &     pq(:,:,igcm_ccn_number)/cloudfrac(:,:)
-c           pqeff(:,:,igcm_h2o_ice)= pq(:,:,igcm_h2o_ice)/
-c     &                               cloudfrac(:,:)
-           pqeff(:,:,igcm_ccnco2_mass) =pq(:,:,igcm_ccnco2_mass)/
-     &                                cloudfrac(:,:)
-           pqeff(:,:,igcm_ccnco2_number)= 
-     &     pq(:,:,igcm_ccnco2_number)/cloudfrac(:,:)
-           pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice)/
-     &                               cloudfrac(:,:)
-
-        ELSE
-           pqeff(:,:,:)=pq(:,:,:)
-c           pqeff(:,:,igcm_ccn_mass)= pq(:,:,igcm_ccn_mass)
-c           pqeff(:,:,igcm_ccn_number)= pq(:,:,igcm_ccn_number)
-c           pqeff(:,:,igcm_h2o_ice)= pq(:,:,igcm_h2o_ice)
-c           pqeff(:,:,igcm_ccnco2_mass)= pq(:,:,igcm_ccnco2_mass)
-c           pqeff(:,:,igcm_ccnco2_number)= pq(:,:,igcm_ccnco2_number)
-c           pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice)
-        END IF      
-        tempo_traceur_t(1:ngrid,1:nlay)=zteff(1:ngrid,1:nlay)
-      tempo_traceurs(1:ngrid,1:nlay,1:nq)=pqeff(1:ngrid,1:nlay,1:nq)
-
+      END IF                    ! end if (CLFvaryingco2)
 c------------------------------------------------------------------
-c Time subsampling for microphysics 
+c microtimestep timeloop for microphysics:
+c 0.Stepped entry for tendancies
+c 1.Compute sedimentation and update tendancies
+c 2.Call co2clouds microphysics
+c 3.Update tendancies 
 c------------------------------------------------------------------
-      DO microstep=1,imicro 
+      DO microstep=1,imicroco2
 c------ Temperature tendency subpdt
-        ! Each microtimestep we give the cloud scheme a stepped entry subpdt instead of pdt
         ! If imicro=1 subpdt is the same as pdt
-        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
-             subpdq(ig,l,igcm_dust_mass) = 
-     &            subpdq(ig,l,igcm_dust_mass)
-     &            + pdq(ig,l,igcm_dust_mass)
-             subpdq(ig,l,igcm_dust_number) = 
-     &            subpdq(ig,l,igcm_dust_number)
-     &            + pdq(ig,l,igcm_dust_number)
-             subpdq(ig,l,igcm_ccnco2_mass) = 
-     &            subpdq(ig,l,igcm_ccnco2_mass)
-     &            + pdq(ig,l,igcm_ccnco2_mass)
-             subpdq(ig,l,igcm_ccnco2_number) = 
-     &            subpdq(ig,l,igcm_ccnco2_number)
-     &            + pdq(ig,l,igcm_ccnco2_number)
-             subpdq(ig,l,igcm_co2_ice) = 
-     &            subpdq(ig,l,igcm_co2_ice)
-     &            + pdq(ig,l,igcm_co2_ice)
-             subpdq(ig,l,igcm_co2) = 
-     &            subpdq(ig,l,igcm_co2)
-     &            + pdq(ig,l,igcm_co2)
-             subpdq(ig,l,igcm_h2o_ice) = 
-     &            subpdq(ig,l,igcm_h2o_ice)
-     &            + pdq(ig,l,igcm_h2o_ice)
-             subpdq(ig,l,igcm_ccn_mass) = 
-     &            subpdq(ig,l,igcm_ccn_mass)
-     &            + pdq(ig,l,igcm_ccn_mass)
-             subpdq(ig,l,igcm_ccn_number) = 
-     &            subpdq(ig,l,igcm_ccn_number)
-     &            + pdq(ig,l,igcm_ccn_number)
-          ENDDO
-       ENDDO
-
-c add meteoritic flux of dust (old version)
-cNew version (John Plane values) are added in improvedCO2clouds
-    !convert meteo_alt (in km) to z-level 
-        !pzlay altitudes of the layers
-c$$$!zlev altitudes (nlayl+1) of the boundaries
-c$$$       if (meteo_flux_number .ge. 0) then 
-c$$$          do ig=1, ngrid
-c$$$             l=1
-c$$$             do  while ( (((meteo_alt .ge. pplev(ig,l)) .and.
-c$$$     &            (meteo_alt .le. pplev(ig,l+1))) .eq. .false.)
-c$$$     &            .and. (l .lt. nlay) )
-c$$$                l=l+1
-c$$$             enddo
-c$$$             meteo_lvl=l
-c$$$             subpdq(ig,meteo_lvl,igcm_dust_mass)=
-c$$$     &            subpdq(ig,meteo_lvl,igcm_dust_mass)
-c$$$     &            +meteo_flux_mass
-c$$$             subpdq(ig,meteo_lvl,igcm_dust_number)=
-c$$$     &            subpdq(ig,meteo_lvl,igcm_dust_number)
-c$$$     &            +meteo_flux_number
-c$$$          enddo
-c$$$       endif
-c-------------------------------------------------------------------
-c   2.  Main call to the cloud schemes:
-c------------------------------------------------
-      CALL improvedCO2clouds(ngrid,nlay,microtimestep,
-     &     pplay,pzlev,zteff,subpdt,
-     &     pqeff,subpdq,subpdqcloudco2,subpdtcloudco2,
-     &     nq,tauscaling,memdMMccn,memdMMh2o,memdNNccn)
-c-------------------------------------------------------------------
-c   3.  Updating tendencies after cloud scheme:
-c-----------------------------------------------
-          DO l=1,nlay
+         DO l=1,nlay
             DO ig=1,ngrid
-               subpdt(ig,l) =
-     &              subpdt(ig,l) + subpdtcloudco2(ig,l)
-               subpdq(ig,l,igcm_dust_mass) =
+               subpdt(ig,l) = subpdt(ig,l)
+     &              + pdt(ig,l) ! At each micro timestep we add pdt in order to have a stepped entry
+               subpdq(ig,l,igcm_dust_mass) = 
      &              subpdq(ig,l,igcm_dust_mass)
-     &              + subpdqcloudco2(ig,l,igcm_dust_mass)
-               subpdq(ig,l,igcm_dust_number) =
+     &              + pdq(ig,l,igcm_dust_mass)
+               subpdq(ig,l,igcm_dust_number) = 
      &              subpdq(ig,l,igcm_dust_number)
-     &              + subpdqcloudco2(ig,l,igcm_dust_number)
-               subpdq(ig,l,igcm_ccnco2_mass) =
+     &              + pdq(ig,l,igcm_dust_number)
+
+               subpdq(ig,l,igcm_ccnco2_mass) = 
      &              subpdq(ig,l,igcm_ccnco2_mass)
-     &              + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
-               subpdq(ig,l,igcm_ccnco2_number) =
+     &              + pdq(ig,l,igcm_ccnco2_mass)
+               subpdq(ig,l,igcm_ccnco2_number) = 
      &              subpdq(ig,l,igcm_ccnco2_number)
-     &              + subpdqcloudco2(ig,l,igcm_ccnco2_number)
-               subpdq(ig,l,igcm_co2_ice) =
+     &              + pdq(ig,l,igcm_ccnco2_number)
+
+               subpdq(ig,l,igcm_co2_ice) = 
      &              subpdq(ig,l,igcm_co2_ice)
-     &              + subpdqcloudco2(ig,l,igcm_co2_ice)
-               subpdq(ig,l,igcm_co2) =
+     &              + pdq(ig,l,igcm_co2_ice)
+               subpdq(ig,l,igcm_co2) = 
      &              subpdq(ig,l,igcm_co2)
-     &              + subpdqcloudco2(ig,l,igcm_co2)
-               subpdq(ig,l,igcm_h2o_ice) =
+     &              + pdq(ig,l,igcm_co2)
+
+               subpdq(ig,l,igcm_h2o_ice) = 
      &              subpdq(ig,l,igcm_h2o_ice)
-     &              + subpdqcloudco2(ig,l,igcm_h2o_ice)
-               subpdq(ig,l,igcm_ccn_mass) =
+     &              + pdq(ig,l,igcm_h2o_ice)
+               subpdq(ig,l,igcm_ccn_mass) = 
      &              subpdq(ig,l,igcm_ccn_mass)
-     &              + subpdqcloudco2(ig,l,igcm_ccn_mass)
-               subpdq(ig,l,igcm_ccn_number) =
+     &              + pdq(ig,l,igcm_ccn_mass)
+               subpdq(ig,l,igcm_ccn_number) = 
      &              subpdq(ig,l,igcm_ccn_number)
-     &              + subpdqcloudco2(ig,l,igcm_ccn_number)
+     &              + pdq(ig,l,igcm_ccn_number)
             ENDDO
          ENDDO
-        
-!Sedimentation
-!Update traceurs, compute rsed 
+c- Effective tracers quantities in the cloud fraction
+         IF (CLFvaryingCO2) THEN     
+            pqeff(:,:,:)=pq(:,:,:) ! prevent from buggs (A. Pottier)
+            pqeff(:,:,igcm_ccnco2_mass) =pq(:,:,igcm_ccnco2_mass)/
+     &           cloudfrac(:,:)
+            pqeff(:,:,igcm_ccnco2_number)= 
+     &           pq(:,:,igcm_ccnco2_number)/cloudfrac(:,:)
+            pqeff(:,:,igcm_co2_ice)= pq(:,:,igcm_co2_ice)/
+     &           cloudfrac(:,:)
+         ELSE
+            pqeff(:,:,:)=pq(:,:,:)
+         END IF  
        
+c------------------------------------------------------
+c 1.SEDIMENTATION : update tracers, compute parameters,
+c   call to sedimentation routine, update tendancies
+c------------------------------------------------------
        DO l=1, nlay
           DO ig=1,ngrid             
@@ -567 +511 @@
              tempo_traceurs(ig,l,:)=pqeff(ig,l,:)
      &            +subpdq(ig,l,:)*microtimestep
-
              rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4*
      &            tempo_traceur_t(ig,l)-2.87e-6*
@@ -578 +521 @@
              Qccnco2=max(tempo_traceurs(ig,l,igcm_ccnco2_mass),
      &            1.e-30)
-             mdustJA= tempo_traceurs(ig,l,igcm_dust_mass) 
-             ndustJA=tempo_traceurs(ig,l,igcm_dust_number)
-             if ((Nccnco2 .lt. tauscaling(ig)) .or. (Qccnco2 .lt. 
-     &            1.e-30 *tauscaling(ig))) then 
-                rdust(ig,l)=1.e-10
-             else
-                rdust(ig,l)=(3./4./pi/2500.*Qccnco2/Nccnco2)**(1./3.)
-                rdust(ig,l)=max(rdust(ig,l),1.e-10) 
-c     rdust(ig,l)=min(rdust(ig,l),5.e-6)   
-             end if
-             rhocloudco2t(ig,l) = (Niceco2 *rho_ice_co2 
-     &            + Qccnco2*tauscaling(ig)*rho_dust)
-     &            / (Niceco2 + Qccnco2*tauscaling(ig))
-             
+             call updaterice_microco2(Niceco2,
+     &            Qccnco2,Nccnco2,
+     &            tauscaling(ig),riceco2(ig,l),rhocloudco2t(ig,l))
+             if (Niceco2 .le. 1.e-25 
+     &            .or. Nccnco2*tauscaling(ig) .le. 1) THEN
+                riceco2(ig,l)=1.e-9
+             endif
              rhocloudco2t(ig,l)=min(max(rhocloudco2t(ig,l)
      &            ,rho_ice_co2),rho_dust)
-             riceco2(ig,l)=(Niceco2*3.0/
-     &            (4.0*rho_ice_co2*pi*Nccnco2
-     &            *tauscaling(ig)) +rdust(ig,l)*rdust(ig,l)
-     &            *rdust(ig,l))**(1.0/3.0)
-             riceco2(ig,l)=max(1.e-10,riceco2(ig,l))
-             
              rsedcloudco2(ig,l)=max(riceco2(ig,l)*
      &            (1.+nuiceco2_sed)*(1.+nuiceco2_sed)*(1.+nuiceco2_sed),
-     &            rdust(ig,l))
-             
+     &            riceco2(ig,l))
           ENDDO
        ENDDO
-       
-!     Gravitational sedimentation
-       
+!     Gravitational sedimentation       
        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)=
-     &      tempo_traceurs(:,:,igcm_ccnco2_number)
-       
-      call newsedim(ngrid,nlay,ngrid*nlay,ngrid*nlay,
+     &      tempo_traceurs(:,:,igcm_ccnco2_number) 
+       !We save actualized tracer values to compute sedimentation tendancies
+       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,
      &     microtimestep,pplev,masse,epaisseur,tempo_traceur_t,
      &     rsedcloudco2,rhocloudco2t,
      &     tempo_traceurs(:,:,igcm_ccnco2_mass),wq,beta) 
-      
       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 l = 1, nlay            !Compute tendencies
          DO ig=1,ngrid
             pdqsed(ig,l,igcm_ccnco2_mass)=
@@ -648 +571 @@
          ENDDO
       ENDDO
-            !pdqsed est la tendance due a la sedimentation
-     
-      DO l=1,nlay
+      !update subtimestep tendencies with sedimentation input
+        DO l=1,nlay
          DO ig=1,ngrid
             subpdq(ig,l,igcm_ccnco2_mass) =
@@ -663 +585 @@
          ENDDO
       ENDDO   
- 
+c------------------------------------------------------
+c      2.  Main call to the cloud schemes:
+c------------------------------------------------------
+      CALL improvedCO2clouds(ngrid,nlay,microtimestep,
+     &     pplay,pplev,zteff,subpdt,
+     &     pqeff,subpdq,subpdqcloudco2,subpdtcloudco2,
+     &     nq,tauscaling,memdMMccn,memdMMh2o,memdNNccn)
+c-----------------------------------------------------
+c      3.  Updating tendencies after cloud scheme:
+c-----------------------------------------------------
+          DO l=1,nlay
+            DO ig=1,ngrid
+               subpdt(ig,l) =
+     &              subpdt(ig,l) + subpdtcloudco2(ig,l)
+
+               subpdq(ig,l,igcm_dust_mass) =
+     &              subpdq(ig,l,igcm_dust_mass)
+     &              + subpdqcloudco2(ig,l,igcm_dust_mass)
+               subpdq(ig,l,igcm_dust_number) =
+     &              subpdq(ig,l,igcm_dust_number)
+     &              + subpdqcloudco2(ig,l,igcm_dust_number)
+
+               subpdq(ig,l,igcm_ccnco2_mass) =
+     &              subpdq(ig,l,igcm_ccnco2_mass)
+     &              + subpdqcloudco2(ig,l,igcm_ccnco2_mass)
+               subpdq(ig,l,igcm_ccnco2_number) =
+     &              subpdq(ig,l,igcm_ccnco2_number)
+     &              + subpdqcloudco2(ig,l,igcm_ccnco2_number)
+
+               subpdq(ig,l,igcm_co2_ice) =
+     &              subpdq(ig,l,igcm_co2_ice)
+     &              + subpdqcloudco2(ig,l,igcm_co2_ice)
+               subpdq(ig,l,igcm_co2) =
+     &              subpdq(ig,l,igcm_co2)
+     &              + subpdqcloudco2(ig,l,igcm_co2)
+
+               subpdq(ig,l,igcm_h2o_ice) =
+     &              subpdq(ig,l,igcm_h2o_ice)
+     &              + subpdqcloudco2(ig,l,igcm_h2o_ice)
+               subpdq(ig,l,igcm_ccn_mass) =
+     &              subpdq(ig,l,igcm_ccn_mass)
+     &              + subpdqcloudco2(ig,l,igcm_ccn_mass)
+               subpdq(ig,l,igcm_ccn_number) =
+     &              subpdq(ig,l,igcm_ccn_number)
+     &              + subpdqcloudco2(ig,l,igcm_ccn_number)
+            ENDDO
+         ENDDO
       ENDDO                     ! of DO microstep=1,imicro
       
-c-------------------------------------------------------------------
-c   6.  Compute final tendencies after time loop:
+c------------------------------------------------
+c   Compute final tendencies after time loop:
 c------------------------------------------------
 c CO2 flux at surface (kg.m-2.s-1)
       do ig=1,ngrid 
-         pdqs_sedco2(ig)=pdqs_sedco2(ig)/real(imicro)
+         pdqs_sedco2(ig)=pdqs_sedco2(ig)/real(imicroco2)
       enddo
-
 c------ Temperature tendency pdtcloud
        DO l=1,nlay
          DO ig=1,ngrid
              pdtcloudco2(ig,l) =
-     &         subpdt(ig,l)/real(imicro)-pdt(ig,l)
+     &         subpdt(ig,l)/real(imicroco2)-pdt(ig,l)
           ENDDO
        ENDDO
-
 c------ Tracers tendencies pdqcloud
        DO l=1,nlay
-         DO ig=1,ngrid
-         
-            pdqcloudco2(ig,l,igcm_co2_ice) = 
-     &        subpdq(ig,l,igcm_co2_ice)/real(imicro)
-     &       - pdq(ig,l,igcm_co2_ice)
-            pdqcloudco2(ig,l,igcm_co2) = 
-     &        subpdq(ig,l,igcm_co2)/real(imicro)
-     &       - pdq(ig,l,igcm_co2)
-            pdqcloudco2(ig,l,igcm_h2o_ice) = 
-     &        subpdq(ig,l,igcm_h2o_ice)/real(imicro)
-     &       - pdq(ig,l,igcm_h2o_ice)
-         ENDDO
+          DO ig=1,ngrid        
+             pdqcloudco2(ig,l,igcm_co2_ice) = 
+     &            subpdq(ig,l,igcm_co2_ice)/real(imicroco2)
+     &            - pdq(ig,l,igcm_co2_ice)
+             pdqcloudco2(ig,l,igcm_co2) = 
+     &            subpdq(ig,l,igcm_co2)/real(imicroco2)
+     &            - pdq(ig,l,igcm_co2)
+             pdqcloudco2(ig,l,igcm_h2o_ice) = 
+     &            subpdq(ig,l,igcm_h2o_ice)/real(imicroco2)
+     &            - pdq(ig,l,igcm_h2o_ice)
+          ENDDO
        ENDDO
-
-        DO l=1,nlay
-         DO ig=1,ngrid
-            pdqcloudco2(ig,l,igcm_ccnco2_mass) = 
-     &        subpdq(ig,l,igcm_ccnco2_mass)/real(imicro)
-     &       - pdq(ig,l,igcm_ccnco2_mass)
-            pdqcloudco2(ig,l,igcm_ccnco2_number) = 
-     &        subpdq(ig,l,igcm_ccnco2_number)/real(imicro)
-     &       - pdq(ig,l,igcm_ccnco2_number)
-            pdqcloudco2(ig,l,igcm_ccn_mass) = 
-     &        subpdq(ig,l,igcm_ccn_mass)/real(imicro)
-     &       - pdq(ig,l,igcm_ccn_mass)
-            pdqcloudco2(ig,l,igcm_ccn_number) = 
-     &        subpdq(ig,l,igcm_ccn_number)/real(imicro)
-     &       - pdq(ig,l,igcm_ccn_number)
+       DO l=1,nlay
+          DO ig=1,ngrid
+             pdqcloudco2(ig,l,igcm_ccnco2_mass) = 
+     &            subpdq(ig,l,igcm_ccnco2_mass)/real(imicroco2)
+     &            - pdq(ig,l,igcm_ccnco2_mass)
+             pdqcloudco2(ig,l,igcm_ccnco2_number) = 
+     &            subpdq(ig,l,igcm_ccnco2_number)/real(imicroco2)
+     &            - pdq(ig,l,igcm_ccnco2_number)
+             pdqcloudco2(ig,l,igcm_ccn_mass) = 
+     &            subpdq(ig,l,igcm_ccn_mass)/real(imicroco2)
+     &            - pdq(ig,l,igcm_ccn_mass)
+             pdqcloudco2(ig,l,igcm_ccn_number) = 
+     &            subpdq(ig,l,igcm_ccn_number)/real(imicroco2)
+     &            - pdq(ig,l,igcm_ccn_number)
+          ENDDO
+       ENDDO
+       DO l=1,nlay
+          DO ig=1,ngrid
+             pdqcloudco2(ig,l,igcm_dust_mass) = 
+     &            subpdq(ig,l,igcm_dust_mass)/real(imicroco2)
+     &            - pdq(ig,l,igcm_dust_mass)
+             pdqcloudco2(ig,l,igcm_dust_number) = 
+     &            subpdq(ig,l,igcm_dust_number)/real(imicroco2)
+     &            - pdq(ig,l,igcm_dust_number)
+          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
+             IF ((pqeff(ig,l,igcm_ccnco2_number) + 
+     &            ptimestep* (pdq(ig,l,igcm_ccnco2_number) + 
+     &            pdqcloudco2(ig,l,igcm_ccnco2_number))
+     &            .lt. 1.)
+     &            .or. (pqeff(ig,l,igcm_ccnco2_mass) + 
+     &            ptimestep* (pdq(ig,l,igcm_ccnco2_mass) + 
+     &            pdqcloudco2(ig,l,igcm_ccnco2_mass))
+     &            .lt. 1.e-20)) THEN
+                pdqcloudco2(ig,l,igcm_ccnco2_number) =
+     &               - pqeff(ig,l,igcm_ccnco2_number)/ptimestep 
+     &               - pdq(ig,l,igcm_ccnco2_number)+1.
+                pdqcloudco2(ig,l,igcm_dust_number) =  
+     &               -pdqcloudco2(ig,l,igcm_ccnco2_number)
+                pdqcloudco2(ig,l,igcm_ccnco2_mass) =
+     &               - pqeff(ig,l,igcm_ccnco2_mass)/ptimestep
+     &               - pdq(ig,l,igcm_ccnco2_mass)+1.e-20
+                pdqcloudco2(ig,l,igcm_dust_mass) = 
+     &               -pdqcloudco2(ig,l,igcm_ccnco2_mass)
+             ENDIF
+          ENDDO
+       ENDDO
+       DO l=1,nlay
+          DO ig=1,ngrid
+             IF ( (pqeff(ig,l,igcm_dust_number) + 
+     &            ptimestep* (pdq(ig,l,igcm_dust_number) + 
+     &            pdqcloudco2(ig,l,igcm_dust_number)) .le. 1.)
+     &            .or. (pqeff(ig,l,igcm_dust_mass)+ 
+     &            ptimestep* (pdq(ig,l,igcm_dust_mass) + 
+     &            pdqcloudco2(ig,l,igcm_dust_mass))
+     &            .le. 1.e-20)) then                  
+                pdqcloudco2(ig,l,igcm_dust_number) =
+     &               - pqeff(ig,l,igcm_dust_number)/ptimestep 
+     &               - pdq(ig,l,igcm_dust_number)+1.
+                pdqcloudco2(ig,l,igcm_ccnco2_number) =  
+     &               -pdqcloudco2(ig,l,igcm_dust_number)
+                pdqcloudco2(ig,l,igcm_dust_mass) =
+     &               - pqeff(ig,l,igcm_dust_mass)/ptimestep
+     &               - pdq(ig,l,igcm_dust_mass) +1.e-20
+                pdqcloudco2(ig,l,igcm_ccnco2_mass) = 
+     &               -pdqcloudco2(ig,l,igcm_dust_mass)
+             ENDIF
+          ENDDO
+       ENDDO
+        !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq      
+       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)) 
+     &       .lt. 1.e-15) THEN
+           pdqcloudco2(ig,l,igcm_co2_ice) = 
+     &     - pqeff(ig,l,igcm_co2_ice)/ptimestep-pdq(ig,l,igcm_co2_ice)
+           pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice)
+          ENDIF    
+          IF (pqeff(ig,l,igcm_co2) + ptimestep*
+     &       (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2)) 
+     &       .lt. 0.1) THEN
+           pdqcloudco2(ig,l,igcm_co2) = 
+     &     - pqeff(ig,l,igcm_co2)/ptimestep - pdq(ig,l,igcm_co2)
+           pdqcloudco2(ig,l,igcm_co2_ice)= -pdqcloudco2(ig,l,igcm_co2)
+          ENDIF
          ENDDO
         ENDDO
 
-       
-        DO l=1,nlay
-         DO ig=1,ngrid
-            pdqcloudco2(ig,l,igcm_dust_mass) = 
-     &        subpdq(ig,l,igcm_dust_mass)/real(imicro)
-     &       - pdq(ig,l,igcm_dust_mass)
-            pdqcloudco2(ig,l,igcm_dust_number) = 
-     &        subpdq(ig,l,igcm_dust_number)/real(imicro)
-     &       - pdq(ig,l,igcm_dust_number)
-         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
-              IF ((pqeff(ig,l,igcm_ccnco2_number) + 
-     &             ptimestep* (pdq(ig,l,igcm_ccnco2_number) + 
-     &             pdqcloudco2(ig,l,igcm_ccnco2_number))
-     &             .lt. 1.e-20)
-     &             .or. (pqeff(ig,l,igcm_ccnco2_mass) + 
-     &             ptimestep* (pdq(ig,l,igcm_ccnco2_mass) + 
-     &             pdqcloudco2(ig,l,igcm_ccnco2_mass))
-     &             .lt. 1.e-30)) THEN
-                 
-                 pdqcloudco2(ig,l,igcm_ccnco2_number) =
-     &                - pqeff(ig,l,igcm_ccnco2_number)/ptimestep 
-     &                - pdq(ig,l,igcm_ccnco2_number)+1.e-20
-                 pdqcloudco2(ig,l,igcm_dust_number) =  
-     &                -pdqcloudco2(ig,l,igcm_ccnco2_number)
-                 pdqcloudco2(ig,l,igcm_ccnco2_mass) =
-     &                - pqeff(ig,l,igcm_ccnco2_mass)/ptimestep
-     &                - pdq(ig,l,igcm_ccnco2_mass)+1.e-30
-                 pdqcloudco2(ig,l,igcm_dust_mass) = 
-     &                -pdqcloudco2(ig,l,igcm_ccnco2_mass)
-
-              ENDIF
-           ENDDO
-        ENDDO
-        
-
-      
-        DO l=1,nlay
-           DO ig=1,ngrid
-              IF ( (pqeff(ig,l,igcm_dust_number) + 
-     &             ptimestep* (pdq(ig,l,igcm_dust_number) + 
-     &             pdqcloudco2(ig,l,igcm_dust_number)) .le. 1.e-30)
-     &             .or. (pqeff(ig,l,igcm_dust_mass)+ 
-     &             ptimestep* (pdq(ig,l,igcm_dust_mass) + 
-     &             pdqcloudco2(ig,l,igcm_dust_mass))
-     &             .le. 1.e-30)) then 
-                 
-                 pdqcloudco2(ig,l,igcm_dust_number) =
-     &                - pqeff(ig,l,igcm_dust_number)/ptimestep 
-     &                - pdq(ig,l,igcm_dust_number)+1.e-30
-                 pdqcloudco2(ig,l,igcm_ccnco2_number) =  
-     &                  -pdqcloudco2(ig,l,igcm_dust_number)
-                 pdqcloudco2(ig,l,igcm_dust_mass) =
-     &                  - pqeff(ig,l,igcm_dust_mass)/ptimestep
-     &                  - pdq(ig,l,igcm_dust_mass) +1.e-30
-                 pdqcloudco2(ig,l,igcm_ccnco2_mass) = 
-     &                -pdqcloudco2(ig,l,igcm_dust_mass)
-
-              ENDIF
-           ENDDO
-        ENDDO
-        !pq+ptime*(pdq+pdqc)=1 ! pdqc=1-pq/ptime-pdq
-c$$$
-c$$$      
-c$$$        DO l=1,nlay
-c$$$         DO ig=1,ngrid
-c$$$          IF (pq(ig,l,igcm_co2_ice) + ptimestep*
-c$$$     &       (pdq(ig,l,igcm_co2_ice) + pdqcloudco2(ig,l,igcm_co2_ice)) 
-c$$$     &       .lt. 1.e-30) THEN
-c$$$           pdqcloudco2(ig,l,igcm_co2_ice) = 
-c$$$     &     - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice)
-c$$$           pdqcloudco2(ig,l,igcm_co2) = -pdqcloudco2(ig,l,igcm_co2_ice)
-c$$$           !write(*,*) "WARNING CO2 ICE in co2cloud.F"
-c$$$
-c$$$          ENDIF    
-c$$$          IF (pq(ig,l,igcm_co2) + ptimestep*
-c$$$     &       (pdq(ig,l,igcm_co2) + pdqcloudco2(ig,l,igcm_co2)) 
-c$$$     &       .lt. 0.5) THEN
-c$$$           pdqcloudco2(ig,l,igcm_co2) = 
-c$$$     &     - pdq(ig,l,igcm_co2_ice) !- pdq(ig,l,igcm_co2)
-c$$$c           pdqcloudco2(ig,l,igcm_co2_ice) = -pdqcloudco2(ig,l,igcm_co2)
-c$$$           pdqcloudco2(ig,l,igcm_co2_ice) = 
-c$$$     &     - pq(ig,l,igcm_co2_ice)/ptimestep - pdq(ig,l,igcm_co2_ice)      
-c$$$          ! write(*,*) "WARNING CO2 in co2cloud.F"
-c$$$          ENDIF
-c$$$         ENDDO
-c$$$        ENDDO
-c$$$  
-        
+c Update clouds parameters values in the cloud fraction (for output)
         DO l=1, nlay
            DO ig=1,ngrid
 
-c     call updaterice_microco2(
-c     &    pq(ig,l,igcm_co2_ice) +                    ! ice mass
-c     &   (pdq(ig,l,igcm_co2_ice) +                   ! ice mass
-c     &    pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep,    ! ice mass
-c     &    pq(ig,l,igcm_ccnco2_mass) +                   ! ccn mass
-c     &   (pdq(ig,l,igcm_ccnco2_mass) +                  ! ccn mass
-c     &    pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep,   ! ccn mass
-c     &    pq(ig,l,igcm_ccnco2_number) +                 ! ccn number
-c     &   (pdq(ig,l,igcm_ccnco2_number) +                ! ccn number
-c     &    pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep, ! ccn number
-c     &    tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
-c        write(*,*) "in co2clouds, riceco2(ig,l)= ",riceco2(ig,l)
               Niceco2=pqeff(ig,l,igcm_co2_ice) +                   
      &             (pdq(ig,l,igcm_co2_ice) + 
@@ -835 +767 @@
               Nccnco2=max((pqeff(ig,l,igcm_ccnco2_number) +                 
      &             (pdq(ig,l,igcm_ccnco2_number) +               
-     &             pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep)*
-     &             tauscaling(ig),1.e-30)
+     &             pdqcloudco2(ig,l,igcm_ccnco2_number))*ptimestep)
+     &             ,1.e-30)
               Qccnco2=max((pqeff(ig,l,igcm_ccnco2_mass) +                 
      &             (pdq(ig,l,igcm_ccnco2_mass) +               
-     &             pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep)*
-     &             tauscaling(ig),1.e-30)
+     &             pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep)
+     &             ,1.e-30)
               
-              if (Nccnco2 .lt. 0.1) then 
-                 rdust(ig,l)=1.e-10
-              else
-        
-                 rdust(ig,l)= Qccnco2
-     &                *0.75/pi/rho_dust
-     &                / Nccnco2
-                 rdust(ig,l)= rdust(ig,l)**(1./3.)            
-                 rdust(ig,l)=max(1.e-10,rdust(ig,l))
-c     rdust(ig,l)=min(5.e-6,rdust(ig,l))
-              endif
               myT=zteff(ig,l)+(pdt(ig,l)+pdtcloudco2(ig,l))*ptimestep
               rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4* 
      &             myT-2.87e-6* myT* myT)
               rho_ice_co2=rho_ice_co2T(ig,l)
-
-              lw = l0 + l1 * myT + l2 *myT**2 + 
-     &             l3 * myT**3 + l4 * myT**4 !J.kg-1
-
-              riceco2(ig,l)=(Niceco2*3.0/
-     &             (4.0*rho_ice_co2*pi*Nccnco2)
-     &             +rdust(ig,l)*rdust(ig,l)
-     &             *rdust(ig,l) )**(1.0/3.0)
-c    &      .or. (riceco2(ig,l) .le. rdust(ig,l))
-              if ( (Niceco2 .le. 1.e-25).or. (Nccnco2 .le. 0.1) )THEN !anciennement 200
-c           riceco2(ig,l)=0.
-
-c     &       .or. (Nccnco2 .le. 1.)
-c        endif
-!Flux chaleur latente <0 quand sublimation  
-
-           pdtcloudco2(ig,l)= pdtcloudco2(ig,l)-Niceco2*lw/cpp/ptimestep
-c$$$  !NO CLOUD : RESET TRACERS AND CONSERVE MASS
-c           if (pq(ig,l,igcm_co2_ice)+(pdq(ig,l,igcm_co2_ice)+
-c     &          pdqcloudco2(ig,l,igcm_co2_ice))*ptimestep .le. 0.) then
-c              pdqcloudco2(ig,l,igcm_co2)=0.
-c              pdqcloudco2(ig,l,igcm_co2_ice)=0.
-c           else
-             pdqcloudco2(ig,l,igcm_co2_ice)=-1.*pqeff(ig,l,igcm_co2_ice)
-     &             /ptimestep-pdq(ig,l,igcm_co2_ice)
-              pdqcloudco2(ig,l,igcm_co2)=-1.*
-     &             pdqcloudco2(ig,l,igcm_co2_ice)
-c           endif
-!     Reverse h2o ccn and ices into h2o tracers
-           if (memdMMccn(ig,l) .gt. 0) then
-              pdqcloudco2(ig,l,igcm_ccn_mass)=memdMMccn(ig,l)/ptimestep
-           else
-              memdMMccn(ig,l)=0.
-              pdqcloudco2(ig,l,igcm_ccn_mass)=0.
-           endif 
-           if (memdNNccn(ig,l) .gt. 0) then
-             pdqcloudco2(ig,l,igcm_ccn_number)=memdNNccn(ig,l)/ptimestep
-           else 
-              memdNNccn(ig,l)=0.
-              pdqcloudco2(ig,l,igcm_ccn_number)=0.
-           endif
-           if (memdMMh2o(ig,l) .gt. 0) then
-              pdqcloudco2(ig,l,igcm_h2o_ice)=memdMMh2o(ig,l)/ptimestep
-           else
-              memdMMh2o(ig,l)=0.
-              pdqcloudco2(ig,l,igcm_h2o_ice)=0.
-           endif
-c           if (pq(ig,l,igcm_ccnco2_mass)+(pdq(ig,l,igcm_ccnco2_mass)+
-c     &          pdqcloudco2(ig,l,igcm_ccnco2_mass))*ptimestep 
-c     &          .le. 1.e-30) then
-c              pdqcloudco2(ig,l,igcm_ccnco2_mass)=0.
-c              pdqcloudco2(ig,l,igcm_ccnco2_number)=0.
-c              pdqcloudco2(ig,l,igcm_co2)=0.
-c              pdqcloudco2(ig,l,igcm_co2_ice)=0.
-c           else
-              pdqcloudco2(ig,l,igcm_ccnco2_mass)=
-     &             -1.*pqeff(ig,l,igcm_ccnco2_mass)
-     &             /ptimestep-pdq(ig,l,igcm_ccnco2_mass)
-c           endif
-c           if (pq(ig,l,igcm_ccnco2_number)+
-c     &          (pdq(ig,l,igcm_ccnco2_number)+
-c     &          pdqcloudco2(ig,l,igcm_ccnco2_number))
-c     &          *ptimestep.le. 1.e-30) 
-c     &          then
-c              pdqcloudco2(ig,l,igcm_ccnco2_mass)=0.
-c              pdqcloudco2(ig,l,igcm_ccnco2_number)=0.
-c              pdqcloudco2(ig,l,igcm_co2)=0.
-c              pdqcloudco2(ig,l,igcm_co2_ice)=0.
-c           else
-              pdqcloudco2(ig,l,igcm_ccnco2_number)=
-     &             -1.*pqeff(ig,l,igcm_ccnco2_number)
-     &             /ptimestep-pdq(ig,l,igcm_ccnco2_number)
-c           endif
-c           if (pq(ig,l,igcm_dust_number)+
-c     &          (pdq(ig,l,igcm_dust_number)+
-c     &          pdqcloudco2(ig,l,igcm_dust_number))
-c     &          *ptimestep.le. 1.e-30) 
-c     &          then
-c              pdqcloudco2(ig,l,igcm_dust_number)=0.
-c              pdqcloudco2(ig,l,igcm_dust_mass)=0.
-c           else
-              pdqcloudco2(ig,l,igcm_dust_number)=
-     &             pqeff(ig,l,igcm_ccnco2_number)
-     &             /ptimestep+pdq(ig,l,igcm_ccnco2_number)
-     &             -memdNNccn(ig,l)/ptimestep
-c           endif
-c           if (pq(ig,l,igcm_dust_mass)+
-c     &          (pdq(ig,l,igcm_dust_mass)+
-c     &          pdqcloudco2(ig,l,igcm_dust_mass))
-c     &          *ptimestep .le. 1.e-30) 
-c     &          then
-c              pdqcloudco2(ig,l,igcm_dust_number)=0.
-c              pdqcloudco2(ig,l,igcm_dust_mass)=0.
-c           else
-              pdqcloudco2(ig,l,igcm_dust_mass)=
-     &             pqeff(ig,l,igcm_ccnco2_mass)
-     &             /ptimestep+pdq(ig,l,igcm_ccnco2_mass)
-     &             -(memdMMh2o(ig,l)+memdMMccn(ig,l))/ptimestep
-c           endif
-           memdMMccn(ig,l)=0.
-           memdMMh2o(ig,l)=0.
-           memdNNccn(ig,l)=0.
-           riceco2(ig,l)=0. 
-        endif
+c             rho_ice_co2 is shared by tracer_mod and used in updaterice
+c     Compute particle size
+              call updaterice_microco2(Niceco2,
+     &             Qccnco2,Nccnco2,
+     &             tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
+              
+              if ( (Niceco2 .le. 1.e-25 .or. 
+     &             Nccnco2*tauscaling(ig) .le. 1.) )THEN
+                 riceco2(ig,l)=0.
+              endif
 c     Compute opacities
-      No=Nccnco2
-      Rn=-log(riceco2(ig,l))
-      n_derf = erf( (rb_cldco2(1)+Rn) *dev2)
-      Qext1bins2(ig,l)=0.
-      do i = 1, nbinco2_cld !Qext below 50 is negligible
-         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
-      !l'opacité de la case ig est la somme sur l de Qext1bins2
-   
-   
+           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(
@@ -988 +813 @@
      &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
           
-
         call updaterdust(
      &    pqeff(ig,l,igcm_dust_mass) +                   ! dust mass
@@ -999 +823 @@
 
          ENDDO
-       ENDDO
-       tau1mic(:)=0.
-       do l=1,nlay
-          do ig=1,ngrid 
-             tau1mic(ig)=tau1mic(ig)+Qext1bins2(ig,l)
-          enddo
-      enddo
-
-c$$$
-c$$$      if (riceco2(725,22) .ge. 1.e-10) then
-c$$$        
-c$$$         write(*,*) 'DIAGJA co2 ',pqeff(725,22,igcm_co2) +                   
-c$$$     &        (pdq(725,22,igcm_co2) + 
-c$$$     &        pdqcloudco2(725,22,igcm_co2))*ptimestep
-c$$$         write(*,*) 'DIAGJA co2_ice',pqeff(725,22,igcm_co2_ice) +
-c$$$     &        (pdq(725,22,igcm_co2_ice) + 
-c$$$     &        pdqcloudco2(725,22,igcm_co2_ice))*ptimestep
-c$$$         
-c$$$         write(*,*) 'DIAGJA riceco2',riceco2(725,22)
-c$$$         write(*,*) 'DIAGJA T',zteff(725,22) +
-c$$$     &        (pdt(725,22) + pdtcloudco2(725,22))*ptimestep
-c$$$         write(*,*) 'DIAG pdtcloud',pdtcloudco2(725,22)*ptimestep
-c$$$         write(*,*) 'DIAGJA ccnNco2',pqeff(725,22,igcm_ccnco2_number)+
-c$$$     &        (pdq(725,22,igcm_ccnco2_number) + 
-c$$$     &        pdqcloudco2(725,22,igcm_ccnco2_number))*ptimestep
-c$$$         
-c$$$         write(*,*) 'DIAGJA dustN',pqeff(725,22,igcm_dust_number) +
-c$$$     &        (pdq(725,22,igcm_dust_number) + 
-c$$$     &        pdqcloudco2(725,22,igcm_dust_number))*ptimestep
-c$$$      ENDIF
-c$$$      
-
-
+       ENDDO 
      
 c     A correction if a lot of subliming CO2 fills the 1st layer FF04/2005
 c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 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
-          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)
+       
+       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
       end do
        
-       
-       DO l=1,nlay
+      DO l=1,nlay
          DO ig=1,ngrid
            rsedcloud(ig,l)=max(rice(ig,l)*
@@ -1071 +863 @@
      &            pdqcloudco2(ig,l,igcm_co2))*ptimestep)*
      &            (mmean(ig,l)/44.01)*pplay(ig,l)/zqsatco2(ig,l)
-             !write(*,*) "In CO2 pt,sat ",pt(ig,l),satuco2(ig,l)
             enddo
          enddo
@@ -1097 +888 @@
      &             pdqcloudco2(ig,l,igcm_co2)*cloudfrac(ig,l)
                   pdtcloudco2(ig,l)=pdtcloudco2(ig,l)*cloudfrac(ig,l)
+                  Qext1bins2(ig,l)=Qext1bins2(ig,l)*cloudfrac(ig,l)
                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
+!Outputs:
+         call WRITEDIAGFI(ngrid,"SatIndex","SatIndex"," ",3,
+     &        SatIndex)
          call WRITEDIAGFI(ngrid,"satuco2","vap in satu","kg/kg",3,
      &        satuco2)
          call WRITEdiagfi(ngrid,"riceco2","ice radius","m"
-     &        ,3,riceco2)
-!     or output in diagfi.nc (for testphys1d)
-c         call WRITEDIAGFI(ngrid,'ps','Surface pressure','Pa',0,ps)
-c         call WRITEDIAGFI(ngrid,'temp','Temperature ',
-c     &                       'K JA',1,pt)
+     &        ,3,riceco2)         
+         call WRITEdiagfi(ngrid,"cloudfrac","co2 cloud fraction"
+     &        ," ",3,cloudfrac)
+         call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2"
+     &        ,"m",3,rsedcloudco2)
+         call WRITEdiagfi(ngrid,"Tau3D1mic"," co2 ice opacities"
+     &        ," ",3,Qext1bins2)
+         call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron"
+     &        ," ",2,tau1mic)
          
-      call WRITEdiagfi(ngrid,"rsedcloudco2","rsed co2","m",3,
-     &   rsedcloudco2)
-      
-      call WRITEdiagfi(ngrid,"tau1mic","co2 ice opacity 1 micron"," ",2,
-     &   tau1mic)
-      call WRITEdiagfi(ngrid,"cloudfrac","co2 cloud fraction"," ",3,
-     &   cloudfrac)
-! used for rad. transfer calculations
-! nuice is constant because a lognormal distribution is prescribed
-c      nuice(1:ngrid,1:nlay)=nuice_ref 
-
-
-
-c=======================================================================
-
-      END
- 
+         END

trunk/LMDZ.MARS/libf/phymars/improvedCO2clouds.F

@@ -1 +1 @@
       subroutine improvedCO2clouds(ngrid,nlay,ptimestep,
-     &             pplay,pzlev,pt,pdt,
+     &             pplay,pplev,pt,pdt,
      &             pq,pdq,pdqcloudco2,pdtcloudco2,
      &             nq,tauscaling,
      &             memdMMccn,memdMMh2o,memdNNccn)
-! to use  'getin'
       USE comcstfi_h
       USE ioipsl_getincom
       USE updaterad
       use tracer_mod
-!, only: rho_ice_co2, nuiceco2_sed, igcm_co2,
-!     &                      rho_ice,igcm_h2o_ice, igcm_ccn_number,
-!     &                      igcm_co2_ice, igcm_dust_mass,
-!     &                      igcm_dust_number, igcm_ccnco2_mass,
-!    &                      igcm_ccnco2_number
       use conc_mod, only: mmean
 
@@ -45 +39 @@
 c Adaptation for CO2 clouds by Joachim Audouard (09/16), based on the work
 c of Constantino Listowski 
+c Warning:
+c If meteoritic particles are activated and turn into co2 ice particles,
+c then they will be reversed in the dust tracers if the cloud sublimates
+ 
 c------------------------------------------------------------------
-!#include "dimensions.h"
-!#include "dimphys.h"
 #include "callkeys.h"
-!#include "tracer.h"
-!#include "comgeomfi.h"
-!#include "dimradmars.h"
 #include "microphys.h"
 #include "datafile.h"
-!#include "microphysCO2.h"
-!#include "conc.h"
 c------------------------------------------------------------------
 c     Inputs:
@@ -63 +54 @@
       REAL ptimestep             ! pas de temps physique (s)
       REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
-      REAL pzlev(ngrid,nlay)     ! altitude au milieu des couches ()
-
+      REAL pplev(ngrid,nlay+1)   ! pression inter couches (Pa)
+      real masse (ngrid,nlay) ! Layer mass (kg.m-2)
       REAL pt(ngrid,nlay)        ! temperature at the middle of the
                                  !   layers (K)
@@ -73 +64 @@
                                  !   (kg/kg.s-1)
       REAL tauscaling(ngrid)     ! Convertion factor for qdust and Ndust
-
       REAL rice(ngrid,nlay)    ! Water Ice mass mean radius (m)
                                 ! used for nucleation of CO2 on ice-coated ccns
       REAL rccnh2o(ngrid,nlay)    ! Water Ice mass mean radius (m)
-
 c     Outputs:
       REAL pdqcloudco2(ngrid,nlay,nq) ! tendance de la condensation
@@ -87 +76 @@
 c------------------------------------------------------------------
 c     Local variables:
-
       LOGICAL firstcall
       DATA firstcall/.true./
       SAVE firstcall
-
       REAL*8   derf ! Error function
-      !external derf
-
-      !REAL*8   massflowrateCO2 
-      !external massflowrateCO2
-     
       INTEGER ig,l,i,flag_pourri
       
@@ -104 +86 @@
       REAL zt(ngrid,nlay)       ! local value of temperature
       REAL zqsat(ngrid,nlay)    ! saturation vapor pressure for CO2
+      real tcond(ngrid,nlay)
+      real zqco2(ngrid,nlay)
       REAL lw                         !Latent heat of sublimation (J.kg-1) 
       REAL,save :: l0,l1,l2,l3,l4
-      REAL cste
       DOUBLE PRECISION dMice           ! mass of condensed ice
       DOUBLE PRECISION sumcheck
-      DOUBLE PRECISION pco2     ! Co2 vapor partial pressure (Pa)
+      DOUBLE PRECISION facteurmax!for energy limit on mass growth
+      DOUBLE PRECISION pco2,psat  ! Co2 vapor partial pressure (Pa)
       DOUBLE PRECISION satu ! Co2 vapor saturation ratio over ice
       DOUBLE PRECISION Mo,No
@@ -116 +100 @@
       DOUBLE PRECISION memdMMh2o(ngrid,nlay)
       DOUBLE PRECISION memdNNccn(ngrid,nlay)
-
+      
 !     Radius used by the microphysical scheme (m)
       DOUBLE PRECISION n_aer(nbinco2_cld) ! number concentration volume-1 of particle/each size bin
@@ -149 +133 @@
 c     Parameters of the size discretization
 c       used by the microphysical scheme
-      DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-11 ! Minimum radius (m)
-      DOUBLE PRECISION, PARAMETER :: rmax_cld = 3.e-6 ! Maximum radius (m)
-      DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-12
-                                           ! Minimum bounary radius (m)
-      DOUBLE PRECISION, PARAMETER :: rbmax_cld = 5.e-6 ! Maximum boundary radius (m)
+      DOUBLE PRECISION, PARAMETER :: rmin_cld = 1.e-9 ! Minimum radius (m)
+      DOUBLE PRECISION, PARAMETER :: rmax_cld = 5.e-6 ! Maximum radius (m)
+      DOUBLE PRECISION, PARAMETER :: rbmin_cld =1.e-10! Minimum bounary radius (m)
+      DOUBLE PRECISION, PARAMETER :: rbmax_cld = 2.e-4 ! Maximum boundary radius (m)
       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 dr_cld(nbinco2_cld)   ! width of each rad_cldco2 bin (m)
+      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)
 
       DOUBLE PRECISION Proba,Masse_atm,drsurdt,reff,Probah2o
-      REAL sigma_iceco2 ! Variance of the ice and CCN distributions
-      SAVE sigma_iceco2
-      
-
-      REAL sigma_ice ! Variance of the ice and CCN distributions
+      REAL sigma_iceco2 ! Variance of the co2 ice and CCN distributions
+      SAVE sigma_iceco2 
+      REAL sigma_ice ! Variance of the h2o ice and CCN distributions
       SAVE sigma_ice
-
       DOUBLE PRECISION Niceco2,Qccnco2,Nccnco2
-      
-      integer coeffh2o
-!Variables for the meteoritic flux
-
-        double precision meteo_ccn(ngrid,nlay,100) !100=nbinco2_cld !!! 
-        double precision,save :: meteo(130,100)
-        double precision mtemp(100)
-        logical file_ok
-        integer altitudes_meteo(130),nelem,lebon1,lebon2
-        double precision :: ltemp1(130),ltemp2(130)
-        integer ibin,uMeteo,j
-c----------------------------------      
-c TESTS
-
-
-
-      REAL satubf(ngrid,nlay),satuaf(ngrid,nlay) 
-      REAL res_out(ngrid,nlay)
- 
-
-c------------------------------------------------------------------
+      integer,save :: coeffh2o  !will be =0 is co2useh2o=.false.
+!     Variables for the meteoritic flux:
+      double precision meteo_ccn(ngrid,nlay,100) !100=nbinco2_cld !!! 
+      double precision,save :: meteo(130,100)
+      double precision mtemp(100),pression_meteo(130)
+      logical file_ok
+      integer nelem,lebon1,lebon2
+      double precision :: ltemp1(130),ltemp2(130)
+      integer ibin,uMeteo,j
 
       IF (firstcall) THEN
@@ -204 +171 @@
 c       dr_cld is the width of each rad_cldco2 bin.
 
-c       Volume ratio between two adjacent bins
-   !     vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
-   !     vrat_cld = exp(vrat_cld)
-        vrat_cld = log(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
-        vrat_cld = exp(vrat_cld)
-c        write(*,*) "vrat_cld", vrat_cld
-
+        vrat_cld = dlog(rmax_cld/rmin_cld) / float(nbinco2_cld-1) *3.
+        vrat_cld = dexp(vrat_cld)
         rb_cldco2(1)  = rbmin_cld
         rad_cldco2(1) = rmin_cld
         vol_cld(1) = 4./3. * dble(pi) * rmin_cld*rmin_cld*rmin_cld
-   !     vol_cld(1) = 4./3. * pi * rmin_cld*rmin_cld*rmin_cld
-
         do i=1,nbinco2_cld-1
           rad_cldco2(i+1)  = rad_cldco2(i) * vrat_cld**(1./3.)
           vol_cld(i+1)  = vol_cld(i) * vrat_cld
-        enddo
-        
+        enddo        
         do i=1,nbinco2_cld
           rb_cldco2(i+1)= ( (2.*vrat_cld) / (vrat_cld+1.) )**(1./3.) *
@@ -229 +188 @@
         dr_cld(nbinco2_cld)   = rb_cldco2(nbinco2_cld+1) -
      &       rb_cldco2(nbinco2_cld)
-
         print*, ' '
         print*,'Microphysics co2: size bin information:'
@@ -240 +198 @@
         write(*,'(i3,3x,e12.6)') nbinco2_cld+1,rb_cldco2(nbinco2_cld+1)
         print*,'-----------------------------------'
-
         do i=1,nbinco2_cld+1
-            rb_cldco2(i) = log(rb_cldco2(i))  !! we save that so that it is not computed
+            rb_cldco2(i) = dlog(rb_cldco2(i))  !! we save that so that it is not computed
                                          !! at each timestep and gridpoint
         enddo
@@ -250 +207 @@
 c        mtetaco2 = 0.952
         write(*,*) 'co2_param contact parameter:', mtetaco2
-
 c       Volume of a co2 molecule (m3)
         vo1 = m0co2 / dble(rho_ice_co2) ! m0co2 et non mco2
+c       below, vo1 will be recomputed using real rho ice co2 (T-dependant)
         vo1co2=vo1 ! AJOUT JA 
 c       Variance of the ice and CCN distributions
@@ -258 +215 @@
         sigma_ice = sqrt(log(1.+nuice_sed))
 
- 
-c        write(*,*) 'Variance of ice & CCN distribs :', sigma_iceco2
-c        write(*,*) 'nuice for sedimentation:', nuiceco2_sed
-c        write(*,*) 'Volume of a co2 molecule:', vo1co2
   
         write(*,*) 'Variance of ice & CCN CO2 distribs :', sigma_iceco2
@@ -275 +228 @@
            coeffh2o=1
         endif
-        meteo_ccn(:,:,:)=0.
-
         if (meteo_flux) then
            write(*,*)
@@ -282 +233 @@
            write(*,*) "meteoritic dust particles are available"
            write(*,*) "for co2 ice nucleation! "
-           write(*,*) "Flux given by J. Plane (altitude,size bins)"
+           write(*,*) "Flux given by J. Plane (pressions,size bins)"
            ! Initialisation of the flux: it is constant and is it saved
-           !We must interpolate the table to the GCM altitudes
+           !We must interpolate the table to the GCM pressures
            INQUIRE(FILE=datafile(1:LEN_TRIM(datafile))//
      &       '/Meteo_flux_Plane.dat', EXIST=file_ok)
@@ -296 +247 @@
      &          '/Meteo_flux_Plane.dat'
      &          ,FORM='formatted')
-!13000 records (130 altitudes x 100 bin sizes)
+!13000 records (130 pressions x 100 bin sizes)
+           read(uMeteo,*) !skip 1 line
+           do i=1,130 
+              read(uMeteo,*) pression_meteo(i)
+              write(*,*) pression_meteo(i)
+           enddo
+           read(uMeteo,*)       !skip 1 line
            do i=1,130 
               do j=1,100        ! les mêmes 100 bins size que la distri nuclea: on touche pas
                  read(uMeteo,'(F12.6)') meteo(i,j)
               enddo
-              altitudes_meteo(i)=i ! On doit maintenant réinterpoler le tableau(130,100) sur les altitudes du GCM (nlay,100)
+!On doit maintenant réinterpoler le tableau(130,100) sur les pressions du GCM (nlay,100)
            enddo
            close(uMeteo)
-        endif !of if meteo_flux
-
+        write(*,*) "File meteo_flux read, end of firstcall in co2 micro"
+        endif                     !of if meteo_flux
+      
+        ! Parameter values for Latent heat computation
         l0=595594d0      
         l1=903.111d0    
@@ -312 +271 @@
         l4=-0.000103183d0
         firstcall=.false.
-        
-        
       END IF
-c      write(*,*) "max memdNN =",maxval(memdNNccn)
 !=============================================================
 ! 1. Initialisation
 !=============================================================
-      !cste = 4*pi*rho_ice*ptimestep !not used for co2
       flag_pourri=0
-
-      res_out(:,:) = 0
+      meteo_ccn(:,:,:)=0.
       rice(:,:) = 1.e-8
       riceco2(:,:) = 1.e-11
@@ -335 +289 @@
      &      pt(1:ngrid,1:nlay) + 
      &      pdt(1:ngrid,1:nlay) * ptimestep 
-c      call WRITEDIAGFI(ngrid,"Ztclouds",
-c     &     "Ztclouds",'K',3,zt)
-c       call WRITEDIAGFI(ngrid,"pdtclouds",
-c     &     "pdtclouds",'K',3,pdt*ptimestep)
-      
       zq(1:ngrid,1:nlay,1:nq) = 
      &      pq(1:ngrid,1:nlay,1:nq) + 
      &      pdq(1:ngrid,1:nlay,1:nq) * ptimestep
-      
-      
       WHERE( zq(1:ngrid,1:nlay,1:nq) < 1.e-30 )
-     &       zq(1:ngrid,1:nlay,1:nq) = 1.e-30
-
-      zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
-      
+     &     zq(1:ngrid,1:nlay,1:nq) = 1.e-30
+         
+         zq0(1:ngrid,1:nlay,1:nq) = zq(1:ngrid,1:nlay,1:nq)
 !=============================================================
 ! 2. Compute saturation
 !=============================================================
-   
-          
-
       dev2 = 1. / ( sqrt(2.) * sigma_iceco2 )
       dev3 = 1. / ( sqrt(2.) * sigma_ice )
-
       call co2sat(ngrid*nlay,zt,pplay,zqsat) !zqsat is psat(co2)
-        
-
-!=============================================================
-!   Bonus: additional meteoritic particles for nucleation
+      zqco2=zq(:,:,igcm_co2)+zq(:,:,igcm_co2_ice)
+      CALL tcondco2(ngrid,nlay,pplay,zqco2,tcond)
+!=============================================================
+! 3. Bonus: additional meteoritic particles for nucleation
 !=============================================================
       if (meteo_flux) then
-         !altitude_meteo(130)
-         !pzlev(ngrid,nlay)
+         !pression_meteo(130)
+         !pplev(ngrid,nlay+1)
          !meteo(130,100)
          !resultat: meteo_ccn(ngrid,nlay,100)
          do l=1,nlay
             do ig=1,ngrid
-               ltemp1=abs(altitudes_meteo(:)-pzlev(ig,l))
-               ltemp2=abs(altitudes_meteo(:)-pzlev(ig,l+1))
+               masse(ig,l)=(pplev(ig,l) - pplev(ig,l+1)) /g 
+               ltemp1=abs(pression_meteo(:)-pplev(ig,l))
+               ltemp2=abs(pression_meteo(:)-pplev(ig,l+1))
                lebon1=minloc(ltemp1,DIM=1) 
                lebon2=minloc(ltemp2,DIM=1)
@@ -381 +324 @@
                   mtemp(ibin)=sum(meteo(lebon1:lebon2,ibin))
                enddo
-               meteo_ccn(ig,l,:)=mtemp(:)/nelem
+               meteo_ccn(ig,l,:)=mtemp(:)/nelem/masse(ig,l) !Par kg air
+csi par m carre, x epaisseur/masse pour par kg/air 
+               !write(*,*) "masse air ig l=",masse(ig,l)
+               !check original unit with J. Plane
             enddo
          enddo
       endif
-
-c     Main loop over the GCM's grid
+c ------------------------------------------------------------------------
+c ---------  Actual microphysics : Main loop over the GCM's grid ---------
+c ------------------------------------------------------------------------
        DO l=1,nlay
          DO ig=1,ngrid
 c       Get the partial pressure of co2 vapor and its saturation ratio
            pco2 = zq(ig,l,igcm_co2) * (mmean(ig,l)/44.01) * pplay(ig,l)
-c        satu = zq(ig,l,igcm_co2) / zqsat(ig,l)
            satu = pco2 / zqsat(ig,l)
-!=============================================================
-! 3. Nucleation
-!=============================================================
+
            rho_ice_co2T(ig,l)=1000.*(1.72391-2.53e-4*zt(ig,l)
-     &          -2.87e-6*zt(ig,l)*zt(ig,l))
+     &          -2.87e-6*zt(ig,l)*zt(ig,l)) !T-dependant CO2 ice density
            vo2co2 = m0co2 / dble(rho_ice_co2T(ig,l))
            rho_ice_co2=rho_ice_co2T(ig,l)
 
-           IF ( satu .ge. 1d0 ) THEN ! if there is condensation
-
+!=============================================================
+!4. Nucleation
+!=============================================================
+           IF ( satu .ge. 1 ) THEN ! if there is condensation
+
+c     call updaterccn(zq(ig,l,igcm_dust_mass),
+c     &             zq(ig,l,igcm_dust_number),rdust(ig,l),tauscaling(ig))
+
+              !We do rdust computation "by hand" because we don't want to
+              ! change the mininumum rdust value in updaterad... It would
+              ! mess up various part of the GCM !
+              
               rdust(ig,l)= zq(ig,l,igcm_dust_mass)
      &             *0.75/pi/rho_dust
      &             / zq(ig,l,igcm_dust_number)
               rdust(ig,l)= rdust(ig,l)**(1./3.)
-c             write(*,*) "Improved2, l,Rdust = ",l,rdust(ig,l)
-            rdust(ig,l)=max(1.e-9,rdust(ig,l))
-c            rdust(ig,l)=min(1.e-5,rdust(ig,l))
-             ! write(*,*) "Improved3,Rdust = ",rdust(ig,l)
+              if (zq(ig,l,igcm_dust_mass)*tauscaling(ig) .le. 1.e-20 
+     &             .or.  zq(ig,l,igcm_dust_number)*tauscaling(ig) .le.1
+     &             .or. rdust(ig,l) .le. 1.e-9) then
+                 rdust(ig,l)=1.e-9
+              endif
+              rdust(ig,l)=min(5.e-4,rdust(ig,l))
       
 c       Expand the dust moments into a binned distribution
@@ -416 +372 @@
               No = zq(ig,l,igcm_dust_number)* tauscaling(ig)+1.e-30
               Rn = rdust(ig,l)
-              Rn = -log(Rn) 
+              Rn = -dlog(Rn) 
               Rm = Rn - 3. * sigma_iceco2*sigma_iceco2  
-              n_derf = erf( (rb_cldco2(1)+Rn) *dev2)
-              m_derf = erf( (rb_cldco2(1)+Rm) *dev2)
+              n_derf = derf( (rb_cldco2(1)+Rn) *dev2)
+              m_derf = derf( (rb_cldco2(1)+Rm) *dev2)
       
               do i = 1, nbinco2_cld
@@ -427 +383 @@
                  m_derf = derf((rb_cldco2(i+1)+Rm) *dev2)
                  n_aer(i) = n_aer(i) + 0.5 * No * n_derf +
-     &                meteo_ccn(ig,l,i) !Ajout meteo_ccn
+     &                meteo_ccn(ig,l,i) !Ajout meteo_ccn particles
                  m_aer(i) = m_aer(i) + 0.5 * Mo * m_derf
                  m_aer2(i)=4./3.*pi*rho_dust
-     &                *n_aer(i)*rad_cldco2(i)*rad_cldco2(i)
+     &                *meteo_ccn(ig,l,i)*rad_cldco2(i)*rad_cldco2(i)
      &                *rad_cldco2(i)
-c                 write(*,*) "diff =",rad_cldco2(i),m_aer(i),m_aer2(i)
               enddo
-c              write(*,*) "Bilan =",sum(m_aer),sum(m_aer2)
+              if (meteo_flux) m_aer(:)=m_aer(:)+m_aer2(:)
               
-c              sumcheck = 0
-c              do i = 1, nbinco2_cld
-c                 sumcheck = sumcheck + n_aer(i)
-c              enddo
-c              sumcheck = abs(sumcheck/No - 1)
-c              if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
-c                 print*, "WARNING, No sumcheck PROBLEM"
-c                 print*, "sumcheck, No, rdust",sumcheck, No,rdust(ig,l)
-c                 print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
-c                 print*, "Dust binned distribution", n_aer
-c                 STOP
-c              endif
-              
-c              sumcheck = 0
-c              do i = 1, nbinco2_cld
-c                 sumcheck = sumcheck + m_aer(i)
-c              enddo
-c              sumcheck = abs(sumcheck/Mo - 1)
-c              if ((sumcheck .gt. 1e-5) .and.  (1./Rn .gt. rmin_cld))
-c     &             then
-c                 print*, "WARNING, Mo sumcheck PROBLEM"
-c                 print*, "sumcheck, Mo",sumcheck, Mo
-c                 print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig,l
-c                 print*, "Dust binned distribution", m_aer
-c                 STOP
-c              endif
-              m_aer(:)=m_aer2(:)
-              
-              
-              rccnh2o(ig,l)= zq(ig,l,igcm_ccn_mass)
-     &             *0.75/pi/rho_dust
-     &             / zq(ig,l,igcm_ccn_number)
-              
-              rice(ig,l)=( zq(ig,l,igcm_h2o_ice)*3.0/
-     &             (4.0*rho_ice_co2*zq(ig,l,igcm_ccn_number)
-     &             *pi*tauscaling(ig)) +rccnh2o(ig,l)*rccnh2o(ig,l)
-     &             *rccnh2o(ig,l) )**(1.0/3.0)
-              rhocloud(ig,l)=( zq(ig,l,igcm_h2o_ice)*rho_ice
-     &            +zq(ig,l,igcm_ccn_mass) *tauscaling(ig)*rho_dust)
-     &            / (zq(ig,l,igcm_h2o_ice)+zq(ig,l,igcm_ccn_mass)
-     &             *tauscaling(ig))
-c              call updaterice_micro(
-c     &             zq(ig,l,igcm_h2o_ice), ! ice mass
-c     &             zq(ig,l,igcm_ccn_mass), ! ccn mass
-c     &             zq(ig,l,igcm_ccn_number), ! ccn number
-c     &             tauscaling(ig),rice(ig,l),rhocloud(ig,l))
-              ! rice radius of CCN + H20 crystal 
-c              write(*,*) "Improved1 Rice=",rice(ig,l)
-              rice(ig,l)=max(1.e-10,rice(ig,l))
-c              rice(ig,l)=min(1.e-5,rice(ig,l))
-              !write(*,*) "Improved2 Rice=",rice(ig,l)
+              !Same but with h2o particles as CCN
+              call updaterice_micro(zq(ig,l,igcm_h2o_ice),
+     &             zq(ig,l,igcm_ccn_mass),zq(ig,l,igcm_ccn_number),
+     &             tauscaling(ig),rice(ig,l),rhocloud(ig,l))
               Mo = zq(ig,l,igcm_h2o_ice) +
      &             zq(ig,l,igcm_ccn_mass)*tauscaling(ig)+1.e-30!*tauscaling(ig)
@@ -492 +400 @@
               No = zq(ig,l,igcm_ccn_number)* tauscaling(ig) + 1.e-30
               Rn = rice(ig,l)
-              Rn = -log(Rn) 
+              Rn = -dlog(Rn) 
               Rm = Rn - 3. * sigma_ice*sigma_ice  
-              n_derf = erf( (rb_cldco2(1)+Rn) *dev3)
-              m_derf = erf( (rb_cldco2(1)+Rm) *dev3)
+              n_derf = derf( (rb_cldco2(1)+Rn) *dev3)
+              m_derf = derf( (rb_cldco2(1)+Rm) *dev3)
               do i = 1, nbinco2_cld
                  n_aer_h2oice(i) = -0.5 * No * n_derf 
@@ -504 +412 @@
                  m_aer_h2oice(i) = m_aer_h2oice(i) + 0.5 * Mo * m_derf 
                  rad_h2oice(i) = rad_cldco2(i)
-                 m_aer_h2oice2(i)=4./3.*pi*rhocloud(ig,l)
-     &                *n_aer_h2oice(i)*rad_h2oice(i)*rad_h2oice(i)
-     &                *rad_h2oice(i)
-                 
-
-c                 write(*,*) "before nuc, i,rad_h2o(i)= ",i,rad_cldco2(i)
-c     &                ,m_aer_h2oice(i),n_aer_h2oice(i)
               enddo
-c              sumcheck = 0
-c              do i = 1, nbinco2_cld
-c                 sumcheck = sumcheck + n_aer_h2oice(i)
-c              enddo
-c              sumcheck = abs(sumcheck/No - 1)
-c              if ((sumcheck .gt. 1e-5).and. (1./Rn .gt. rmin_cld)) then
-c                 print*, "WARNING, Noh2o sumcheck PROBLEM"
-c                 print*, "sumcheck, No, rice",sumcheck, No,rice(ig,l)
-c                 print*, "min radius, Rn, ig, l", rmin_cld, 1./Rn, ig, l
-c                 print*, "Dust binned distribution", n_aer_h2oice
-c                 STOP
-c              endif
-              
-c            sumcheck = 0
-c              do i = 1, nbinco2_cld
-c                 sumcheck = sumcheck + m_aer_h2oice(i)
-c              enddo
-c              sumcheck = abs(sumcheck/Mo - 1)
-c              if ((sumcheck .gt. 1e-5) .and.  (1./Rn .gt. rmin_cld))
-c     &             then
-c                 print*, "WARNING, Moh2o sumcheck PROBLEM"
-c                 print*, "sumcheck, Mo",sumcheck, Mo
-c                 print*, "min radius, Rm, ig, l", rmin_cld, 1./Rm, ig,l
-c                 print*, "Dust binned distribution", m_aer_h2oice
-c                 STOP
-c              endif
-c       Get the rates of nucleation
-
+              !Call to nucleation routine
               call nucleaCO2(dble(pco2),zt(ig,l),dble(satu)
      &             ,n_aer,rate,n_aer_h2oice
      &             ,rad_h2oice,rateh2o,vo2co2)
-              m_aer_h2oice(:)=m_aer_h2oice2(:)
               dN = 0.
               dM = 0.
@@ -550 +423 @@
               do i = 1, nbinco2_cld
                  Proba=1.0-dexp(-1.*ptimestep*rate(i))
-                 Probah2o=coeffh2o*(1.0-dexp(-1.*ptimestep*rateh2o(i)))
+                 Probah2o=coeffh2o*(1.0-dexp(-1.*ptimestep*rateh2o(i))) !if co2useh2o=.false., this is =0
                  dNh2o    = dNh2o + n_aer_h2oice(i) * Probah2o
                  dMh2o    = dMh2o + m_aer_h2oice(i) * Probah2o
                  dN       = dN + n_aer(i) * Proba
-                 dM       = dM + m_aer(i) * Proba
-             
+                 dM       = dM + m_aer(i) * Proba             
               enddo
         ! dM  masse activée (kg) et dN nb particules par  kg d'air
-            
+        ! Now increment CCN tracers and update dust tracers
               dNN= dN/tauscaling(ig)
               dMM= dM/tauscaling(ig)
               dNN=min(dNN,zq(ig,l,igcm_dust_number))
               dMM=min(dMM,zq(ig,l,igcm_dust_mass))
-           !   if (dNN .gt. 0) then
-           !      write(*,*) "Nuclea dNN crees=",dNN
-           !      write(*,*) "Nuclea dMM crees=",dMM
-           !   endif
-              dNNh2o=dNh2o/tauscaling(ig)
-              dNNh2o=min(dNNh2o,zq(ig,l,igcm_ccn_number))
-
-              ratioh2o_ccn=1./(zq(ig,l,igcm_h2o_ice)
-     &             +zq(ig,l,igcm_ccn_mass)*tauscaling(ig))          
-    
-              dMh2o_ice=dMh2o*zq(ig,l,igcm_h2o_ice)*ratioh2o_ccn
-              dMh2o_ccn=dMh2o*zq(ig,l,igcm_ccn_mass)*
-     &             tauscaling(ig)*ratioh2o_ccn
-  
-              dMh2o_ccn=dMh2o_ccn/tauscaling(ig)
-              dMh2o_ccn=min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
-              dMh2o_ice=min(dMh2o_ice,zq(ig,l,igcm_h2o_ice))
-
               zq(ig,l,igcm_ccnco2_mass)   = 
      &             zq(ig,l,igcm_ccnco2_mass)   + dMM
               zq(ig,l,igcm_ccnco2_number) =
      &             zq(ig,l,igcm_ccnco2_number) + dNN
-
               zq(ig,l,igcm_dust_mass)= zq(ig,l,igcm_dust_mass)-dMM 
               zq(ig,l,igcm_dust_number)=zq(ig,l,igcm_dust_number)-dNN
@@ -591 +444 @@
 c     Update CCN for CO2 nucleating on H2O CCN :
               ! Warning: must keep memory of it 
-              zq(ig,l,igcm_ccnco2_mass)   = 
-     &             zq(ig,l,igcm_ccnco2_mass)  + dMh2o_ice+dMh2o_ccn
-              zq(ig,l,igcm_ccnco2_number) = 
-     &             zq(ig,l,igcm_ccnco2_number) + dNNh2o
-
-              zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)-dNNh2o 
-              zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)-dMh2o_ice
-              zq(ig,l,igcm_ccn_mass)= zq(ig,l,igcm_ccn_mass)-dMh2o_ccn
-          
-
-              memdMMh2o(ig,l)=memdMMh2o(ig,l)+dMh2o_ice
-              memdMMccn(ig,l)=memdMMccn(ig,l)+dMh2o_ccn
-              memdNNccn(ig,l)=memdNNccn(ig,l)+dNNh2o
+              if (co2useh2o) then
+                 dNNh2o=dNh2o/tauscaling(ig)
+                 dNNh2o=min(dNNh2o,zq(ig,l,igcm_ccn_number))
+                 ratioh2o_ccn=1./(zq(ig,l,igcm_h2o_ice)
+     &                +zq(ig,l,igcm_ccn_mass)*tauscaling(ig))  
+                 dMh2o_ice=dMh2o*zq(ig,l,igcm_h2o_ice)*ratioh2o_ccn
+                 dMh2o_ccn=dMh2o*zq(ig,l,igcm_ccn_mass)*
+     &                tauscaling(ig)*ratioh2o_ccn
+                 dMh2o_ccn=dMh2o_ccn/tauscaling(ig)
+                 dMh2o_ccn=min(dMh2o_ccn,zq(ig,l,igcm_ccn_mass))
+                 dMh2o_ice=min(dMh2o_ice,zq(ig,l,igcm_h2o_ice))
+                 zq(ig,l,igcm_ccnco2_mass)   = 
+     &                zq(ig,l,igcm_ccnco2_mass)  + dMh2o_ice+dMh2o_ccn
+                 zq(ig,l,igcm_ccnco2_number) = 
+     &                zq(ig,l,igcm_ccnco2_number) + dNNh2o
+                zq(ig,l,igcm_ccn_number)=zq(ig,l,igcm_ccn_number)-dNNh2o 
+                zq(ig,l,igcm_h2o_ice) = zq(ig,l,igcm_h2o_ice)-dMh2o_ice
+                zq(ig,l,igcm_ccn_mass)= zq(ig,l,igcm_ccn_mass)-dMh2o_ccn
+                memdMMh2o(ig,l)=memdMMh2o(ig,l)+dMh2o_ice
+                memdMMccn(ig,l)=memdMMccn(ig,l)+dMh2o_ccn
+                memdNNccn(ig,l)=memdNNccn(ig,l)+dNNh2o
+             endif !of if co2useh2o
            ENDIF                ! of is satu >1
 !=============================================================
@@ -612 +474 @@
 c Indeed, if we are supersaturated and still don't have at least one nuclei, we should better wait
 c to avoid unrealistic value for nuclei radius and so on for cases that remain negligible.
-           IF (zq(ig,l,igcm_ccnco2_number)*tauscaling(ig) .ge. 1)THEN 
-           ! we trigger crystal growth
-c
-c              write(*,*) "ccn number mass=",zq(ig,l,igcm_ccnco2_number),
-c     &             zq(ig,l,igcm_ccnco2_mass)
-
-c              Niceco2 = zq(ig,l,igcm_co2_ice)
-c              Qccnco2 = zq(ig,l,igcm_ccnco2_mass)
-c              Nccnco2 = zq(ig,l,igcm_ccnco2_number)
-c              call updaterice_microco2(Niceco2,Qccnco2,Nccnco2,
-c     &             tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) 
-c              write(*,*) "updater rice=",riceco2(ig,l)
-              
-              rdust(ig,l)= zq(ig,l,igcm_ccnco2_mass)
-     &             *0.75/pi/rho_dust
-     &             / zq(ig,l,igcm_ccnco2_number)
-              rdust(ig,l)= rdust(ig,l)**(1./3.)            
-              rdust(ig,l)=max(1.e-10,rdust(ig,l))
-c              rdust(ig,l)=min(5.e-6,rdust(ig,l))
-
-              riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/
-     &             (4.0*rho_ice_co2*zq(ig,l,igcm_ccnco2_number)
-     &             *pi*tauscaling(ig)) +rdust(ig,l)*rdust(ig,l)
-     &             *rdust(ig,l) )**(1.0/3.0)
-
-c              riceco2(ig,l)=max(1.e-10,riceco2(ig,l))
-c              riceco2(ig,l)=min(1.e-5,riceco2(ig,l))
-          ! WATCH OUT: CO2 nuclei is supposed to be dust
-          ! only when deriving rhocloud (otherwise would need to keep info on  water embedded in co2) - listo
-c              write(*,*) "Rdust before growth = ",rdust(ig,l)
-c              write(*,*) "Riceco2 before growth = ",riceco2(ig,l)
-
-              !! Niceco2,Qccnco2,Nccnco2
-c              Niceco2 = zq(ig,l,igcm_co2_ice)
-c              Qccnco2 = zq(ig,l,igcm_ccnco2_mass)
-c              Nccnco2 = zq(ig,l,igcm_ccnco2_number)
-c              call updaterice_microco2(Niceco2,Qccnco2,Nccnco2,
-c     &             tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) 
-c              write(*,*) "Riceco2 before growth = ",riceco2(ig,l)
-c              write(*,*) "rdust before growth = ",rdust(ig,l)
-c               write(*,*) "co2ice before growth=",zq(ig,l,igcm_co2_ice)
-c               write(*,*) "co2 before growth=",zq(ig,l,igcm_co2)
-c              write(*,*) "pplay before growth=",pplay(ig,l)
-c              write(*,*) "zt before growth =",zt(ig,l)
-c              write(*,*) "Satu before growth=",satu
-c              riceco2(ig,l)=max(riceco2(ig,l),rdust(ig,l))
-              No   = zq(ig,l,igcm_ccnco2_number)* tauscaling(ig)+1.e-30
-! No nb de particules de poussieres mis à l'échelle pour donner une opacité optique
-
-c       saturation at equilibrium
-c       rice should not be too small, otherwise seq value is not valid
-c              seq  = exp(2.*sigco2*mco2 / (rho_ice_co2*rgp*zt(ig,l)*
-c     &             max(riceco2(ig,l),1.e-7))) !Exponant sans unité OK
-
-ccccccc  Scheme of microphys. mass growth for CO2
-
+           No   = zq(ig,l,igcm_ccnco2_number)* tauscaling(ig)+1.e-30
+           IF (No .ge. 1)THEN   ! we trigger crystal growth
+              call updaterice_microco2(zq(ig,l,igcm_co2_ice),
+     &            zq(ig,l,igcm_ccnco2_mass),zq(ig,l,igcm_ccnco2_number),
+     &            tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l))
+
+              Ic_rice=0.
+              flag_pourri=0
+              lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 + 
+     &             l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4 !J.kg-1
+              facteurmax=abs(Tcond(ig,l)-zt(ig,l))*cpp/lw
+              !specific heat of co2 ice = 1000 J.kg-1.K-1 
+              !specific heat of atm cpp = 744.5 J.kg-1.K-1
+
+ccccccc  call scheme of microphys. mass growth for CO2
               call massflowrateCO2(pplay(ig,l),zt(ig,l),
-     &             satu,riceco2(ig,l),mmean(ig,l),Ic_rice) ! Mass transfer rate (kg/s) for a rice particle
-         ! Ic_rice mass flux kg.s-1 <0 si croissance ! 
-              if (isnan(Ic_rice)) then 
+     &             satu,riceco2(ig,l),mmean(ig,l),Ic_rice) 
+c     Ic_rice Mass transfer rate (kg/s) for a rice particle >0 si croissance ! 
+             
+              if (isnan(Ic_rice) .or. Ic_rice .eq. 0.) then 
                  Ic_rice=0.
                  flag_pourri=1
+                 pdtcloudco2(ig,l)=-pdt(ig,l)
+                 dMice=0
+                 
+              else
+                 dMice=zq(ig,l,igcm_ccnco2_number)*Ic_rice*ptimestep
+     &                *tauscaling(ig) ! Kg par kg d'air, >0 si croissance !
+                 !kg.s-1 par particule * nb particule par kg air*s
+                 ! = kg par kg air 
+              
+              dMice = max(dMice,max(-facteurmax,-zq(ig,l,igcm_co2_ice)))
+              dMice = min(dMice,min(facteurmax,zq(ig,l,igcm_co2)))
+!facteurmax maximum quantity of CO2 that can sublime/condense according to available thermal energy
+! latent heat release       >0 if growth i.e. if dMice >0
+              pdtcloudco2(ig,l)=dMice*lw/cpp/ptimestep
+! kgco2/kgair* J/kgco2 * 1/(J.kgair-1.K-1)/s= K par seconde
+              !Now update tracers 
+              zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice)+dMice
+              zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)-dMice
               endif
-c     drsurdt=-1.0/(4.0*pi*riceco2(ig,l)*
-c     &             riceco2(ig,l)*rho_ice_co2)*Ic_rice
-              dMice =  No * Ic_rice*ptimestep ! Kg par kg d'air, <0 si croissance !           
-c              write(*,*) "dMicev0 in improved = " , dMice
-
-             if (dMice .ge. 0d0) then
-                dMice = min(dMice,abs(zq(ig,l,igcm_co2_ice)))
-             else
-                dMice =-1.* min(abs(dMice),abs(zq(ig,l,igcm_co2)))
-             endif
-c             riceco2(ig,l)=riceco2(ig,l)+drsurdt*ptimestep
-c              write(*,*) "riceco2+dr/dt = ", riceco2(ig,l)
-              !write(*,*) "dMice in improved = " , dMice             
-              
-             zq(ig,l,igcm_co2_ice) = zq(ig,l,igcm_co2_ice)-dMice
-             zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2)+dMice
-
-! latent heat release       >0 if growth i.e. if dMice <0
-
-
-              lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 + 
-     &             l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4 !J.kg-1
-c              write(*,*) "CPP= ",cpp    ! = 744.5
-              !Peut etre un probleme de signe ici!
-              pdtcloudco2(ig,l)= -1.*dMice*lw/cpp/ptimestep ! kg par couche * J par kg /J par K / s = K par seconde 
-              
-              !write(*,*) "pdtcloudco2 after growth = ",pdtcloudco2(ig,l)
-              
-              !write(*,*) "co2 after growth = ",zq(ig,l,igcm_co2)
-              !write(*,*) "co2_ice after growth = ",zq(ig,l,igcm_co2_ice)
-
-          !deltaT par condens/subli. qui remplace le dT du CO2 du newcondens pré-constantino
-          !PDT should be in K/s. 
-!=============================================================
-! 5. Dust cores released, tendancies, latent heat, etc ...
-!=============================================================
-
-c              if (abs(dMice) .ge. 1.e-10) then
-              
-c                 write(*,*) "DIAG zq pdt",(zq(ig,l,igcm_co2_ice)- 
-c     &         zq0(ig,l,igcm_co2_ice))/ptimestep,pdtcloudco2(ig,l)
-c           endif
-           ENDIF                ! of if NCCN > 1
-
-              rdust(ig,l)= zq(ig,l,igcm_ccnco2_mass)
-     &             *0.75/pi/rho_dust
-     &             / zq(ig,l,igcm_ccnco2_number)
-              rdust(ig,l)= rdust(ig,l)**(1./3.)            
-              rdust(ig,l)=max(1.e-9,rdust(ig,l))
-c              rdust(ig,l)=min(5.e-6,rdust(ig,l))
-    
-              riceco2(ig,l)=( zq(ig,l,igcm_co2_ice)*3.0/
-     &             (4.0*rho_ice_co2*zq(ig,l,igcm_ccnco2_number)
-     &             *pi*tauscaling(ig)) +rdust(ig,l)*rdust(ig,l)
-     &             *rdust(ig,l) )**(1.0/3.0)
-              !Niceco2 = zq(ig,l,igcm_co2_ice)
-              !Qccnco2 = zq(ig,l,igcm_ccnco2_mass)
-              !Nccnco2 = zq(ig,l,igcm_ccnco2_number)
-c              
-c              call updaterice_microCO2(Niceco2,Qccnco2,Nccnco2,
-c     &             tauscaling(ig),riceco2(ig,l),rhocloudco2(ig,l)) 
-
-c         If there is no more co2 ice, all the ice particles sublimate, all the condensation nuclei are released:
-
-              if ((zq(ig,l,igcm_co2_ice).le. 1.e-25)) then 
-c  &             .or. flag_pourri .eq. 1
-c     &             .or.(riceco2(ig,l) .le. rdust(ig,l))
-c     &             .or. (l .ge.1 .and. l .le. 5)
-c     &             .or. (zq(ig,l,igcm_co2_ice) .ge. 0.1)
-                 lw = l0 + l1 * zt(ig,l) + l2 * zt(ig,l)**2 + 
-     &                l3 * zt(ig,l)**3 + l4 * zt(ig,l)**4 !J.kg-1
-c              write(*,*) "CPP= ",cpp    ! = 744.5
-              
-              pdtcloudco2(ig,l)=pdtcloudco2(ig,l)-1.
-     &               *zq(ig,l,igcm_co2_ice)*lw/cpp/ptimestep ! 
- !On sublime tout !
-c        write(*,*) "Riceco2 improved before reset=",riceco2(ig,l)
-c        write(*,*) "Niceco2 improved before reset=",
-c     &       zq(ig,l,igcm_co2_ice)
-c        write(*,*) "Rdust improved before reset=",rdust(ig,l)
-                
+
+!=============================================================
+! 5. Dust cores releasing if no more co2 ice :
+!=============================================================
+
+              if (zq(ig,l,igcm_co2_ice).le. 1.e-25)THEN
+ !On sublime tout
+                 if (co2useh2o) then
                  if (memdMMccn(ig,l) .gt. 0) then
                     zq(ig,l,igcm_ccn_mass)=zq(ig,l,igcm_ccn_mass)
@@ -769 +536 @@
      &                   +memdNNccn(ig,l)
                  endif
-                 
-c                 if (zq(ig,l,igcm_ccnco2_mass) .gt. 1.e-30) then
+              endif
                     zq(ig,l,igcm_dust_mass) = 
      &                   zq(ig,l,igcm_dust_mass)
      &                   + zq(ig,l,igcm_ccnco2_mass)-
      &                   (memdMMh2o(ig,l)+memdMMccn(ig,l))
-c                 endif
-c                 if (zq(ig,l,igcm_ccnco2_number) .gt. 1.e-30) then
                     zq(ig,l,igcm_dust_number) = 
      &                   zq(ig,l,igcm_dust_number)
      &                   + zq(ig,l,igcm_ccnco2_number)-memdNNccn(ig,l)
-c                 endif
                  
-c                 if (zq(ig,l,igcm_co2_ice) .gt. 1.e-30) then
                     zq(ig,l,igcm_co2) = zq(ig,l,igcm_co2) 
      &                   + zq(ig,l,igcm_co2_ice)
-c                 endif
                  
                  zq(ig,l,igcm_ccnco2_mass)=0.
@@ -795 +556 @@
                  riceco2(ig,l)=0.
 
-              endif
+              endif !of if co2_ice <1e-25
+
+              ENDIF                ! of if NCCN > 1
           ENDDO                ! of ig loop
         ENDDO                   ! of nlayer loop  
-            ! Get cloud tendencies
-      pdqcloudco2(1:ngrid,1:nlay,igcm_co2) =
-     &     (zq(1:ngrid,1:nlay,igcm_co2) - 
-     &     zq0(1:ngrid,1:nlay,igcm_co2))/ptimestep
-      pdqcloudco2(1:ngrid,1:nlay,igcm_co2_ice) =
-     &     (zq(1:ngrid,1:nlay,igcm_co2_ice) -
-     &     zq0(1:ngrid,1:nlay,igcm_co2_ice))/ptimestep
-      pdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) =
-     &     (zq(1:ngrid,1:nlay,igcm_h2o_ice) - 
-     &     zq0(1:ngrid,1:nlay,igcm_h2o_ice))/ptimestep
-      pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) =
-     &     (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
-     &     zq0(1:ngrid,1:nlay,igcm_ccn_mass))/ptimestep
-      pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) =
-     &     (zq(1:ngrid,1:nlay,igcm_ccn_number) -
-     &     zq0(1:ngrid,1:nlay,igcm_ccn_number))/ptimestep
-      pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) =
-     &     (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) -
-     &     zq0(1:ngrid,1:nlay,igcm_ccnco2_mass))/ptimestep
-      pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_number) =
-     &     (zq(1:ngrid,1:nlay,igcm_ccnco2_number) -
-     &     zq0(1:ngrid,1:nlay,igcm_ccnco2_number))/ptimestep
-      pdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) =
-     &     (zq(1:ngrid,1:nlay,igcm_dust_mass) -
-     &     zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep 
-      pdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) =
-     &     (zq(1:ngrid,1:nlay,igcm_dust_number) -
-     &     zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep
-      return
-      end
-      
-      
-      
+
+!=============================================================
+! 6. END: get cloud tendencies 
+!=============================================================
+
+          ! Get cloud tendencies
+        pdqcloudco2(1:ngrid,1:nlay,igcm_co2) =
+     &       (zq(1:ngrid,1:nlay,igcm_co2) - 
+     &       zq0(1:ngrid,1:nlay,igcm_co2))/ptimestep
+        pdqcloudco2(1:ngrid,1:nlay,igcm_co2_ice) =
+     &       (zq(1:ngrid,1:nlay,igcm_co2_ice) -
+     &       zq0(1:ngrid,1:nlay,igcm_co2_ice))/ptimestep
+        pdqcloudco2(1:ngrid,1:nlay,igcm_h2o_ice) =
+     &       (zq(1:ngrid,1:nlay,igcm_h2o_ice) - 
+     &       zq0(1:ngrid,1:nlay,igcm_h2o_ice))/ptimestep
+        pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_mass) =
+     &       (zq(1:ngrid,1:nlay,igcm_ccn_mass) -
+     &       zq0(1:ngrid,1:nlay,igcm_ccn_mass))/ptimestep
+        pdqcloudco2(1:ngrid,1:nlay,igcm_ccn_number) =
+     &       (zq(1:ngrid,1:nlay,igcm_ccn_number) -
+     &       zq0(1:ngrid,1:nlay,igcm_ccn_number))/ptimestep
+        pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_mass) =
+     &       (zq(1:ngrid,1:nlay,igcm_ccnco2_mass) -
+     &       zq0(1:ngrid,1:nlay,igcm_ccnco2_mass))/ptimestep
+        pdqcloudco2(1:ngrid,1:nlay,igcm_ccnco2_number) =
+     &       (zq(1:ngrid,1:nlay,igcm_ccnco2_number) -
+     &       zq0(1:ngrid,1:nlay,igcm_ccnco2_number))/ptimestep
+        pdqcloudco2(1:ngrid,1:nlay,igcm_dust_mass) =
+     &       (zq(1:ngrid,1:nlay,igcm_dust_mass) -
+     &       zq0(1:ngrid,1:nlay,igcm_dust_mass))/ptimestep 
+        pdqcloudco2(1:ngrid,1:nlay,igcm_dust_number) =
+     &       (zq(1:ngrid,1:nlay,igcm_dust_number) -
+     &       zq0(1:ngrid,1:nlay,igcm_dust_number))/ptimestep
+        return
+        end
+      
+      
+      

trunk/LMDZ.MARS/libf/phymars/massflowrateCO2.F

@@ -1 @@
-
-
 c=======================================================================
       subroutine massflowrateCO2(P,T,Sat,Radius,Matm,Ic)
 c
-c     Determination of the mass transfer rate
+c     Determination of the mass transfer rate for CO2 condensation & 
+c     sublimation
 c
 c     newton-raphson method
-
 c     CLASSICAL  (no SF etc.)
-      
-c     AUTOMATIC SETTING OF RANGES FOR NEWTON-RAPHSON FOR THE PAPER
-
-c     MASS FLUX Ic 
-
+c 
+c   inputs: Pressure (P), Temperature (T), saturation ratio (Sat),
+c           particle radius (Radius), molecular mass of the atm (Matm)
+c   output:  MASS FLUX Ic
+c
+c   Authors: C. Listowski (2014) then J. Audouard (2016-2017)
 c=======================================================================
       USE comcstfi_h
@@ -20 +19 @@
 
       include "microphys.h"
-c      include "microphysCO2.h"
 
 
@@ -26 +24 @@
 c   arguments: INPUT
 c   ----------
-
       REAL T,Matm
       REAL*8 SAT
@@ -33 +30 @@
 c   arguments: OUTPUT
 c   ----------
-
       DOUBLE PRECISION   Ic
-
 c   Local Variables
 c   ----------
-
       DOUBLE PRECISION   Tcm
       DOUBLE PRECISION   T_inf, T_sup, T_dT
@@ -45 +39 @@
       DOUBLE PRECISION   rtsafe
       DOUBLE PRECISION   left, fval, dfval
-
 c    function for newton-raphson iterative method
 c    --------------------------
-
       EXTERNAL classical
-
          
-      Tcm = 80!dble(T)    ! initialize pourquoi 0 et pas t(i)
-
+      Tcm = 110!dble(T)    ! initialize pourquoi 0 et pas t(i)
       T_inf = 0d0
       T_sup = 800d0
-
-      T_dT = 0.01  ! precision - mettre petit et limiter nb iteration? 
+      T_dT = 0.001  ! precision - mettre petit et limiter nb iteration? 
       
 666   CONTINUE
 
-c      print*, 'Radius ', Radius
-c      print*, 'SAT = ', Sat
       call  coefffunc(P,T,Sat,Radius,Matm,kmix,Lsub,C0,C1,C2)
 
-      if (isnan(C0) .eqv. .true.)  C0=0d0
-
-   
+      if (isnan(C0) .eqv. .true.)  C0=0d0   
 c     FIND SURFACE TEMPERATURE (Tc) : iteration sur t 
-
       cond = 4.*pi*Radius*kmix
- 
       Tcm = rtsafe(classical,T_inf,T_sup,T_dT,Radius,C0,C1,C2)
-
-
-      if (Tcm.LE.0d0) then ! unsignificant cases where S<<<Seq and Ncores <<1e-10
-
-            Tcm = 0d0
-
+      if (Tcm.LE.0d0 .or. isnan(Tcm)) then ! unsignificant cases where S<<<Seq and Ncores <<1e-10
+         Tcm = T
       endif
-
-
 c     THEN COMPUTE MASS FLUX Ic from FINAL Tsurface (Tc)
-
       Ic = (Tcm-T)
-
-      Ic = cond*Ic/(-Lsub)
+      Ic = cond*Ic/Lsub
 c regarder de combien varie la solution Ic entre Tcm et Tcm+T_dT
-    
+  
       RETURN
 
@@ -95 +70 @@
 
 c****************************************************************
-
       FUNCTION rtsafe(funcd,x1,x2,xacc,Radius,C0,C1,C2)
-*
 *
 *     Newton Raphsen routine (Numerical Recipe)
@@ -114 +87 @@
       EXTERNAL funcd
 
-      PARAMETER (MAXIT=50000)   !Maximum allowed number of iterations. Using a combination of Newton-Raphson and bisection, 
+      PARAMETER (MAXIT=10000)   !Maximum allowed number of iterations. Using a combination of Newton-Raphson and bisection, 
                                 !find the root of a function bracketed between x1 and x2. The root, returned as the function value rtsafe,
                                 !will be refined until its accuracy is known within !!±xacc. funcd is a user-supplied subroutine which 
@@ -128 +101 @@
 
       if ((fl.gt.0..and.fh.gt.0.).or.(fl.lt.0..and.fh.lt.0.) ) then
-         
-         
          x1=0d0
-         x2=500d0
-         
+         x2=1200d0         
          call funcd(x1,fl,df,C0,C1,C2)
          call funcd(x2,fh,df,C0,C1,C2)
-         
-         write(*,*) 'root must be bracketed in rtsafe'
+c         write(*,*) 'root must be bracketed in rtsafe'
       endif
 
 
       if (fl.eq.0.) then
-
          rtsafe=x1
          return 
-
       else if (fh.eq.0.) then
-
          rtsafe=x2
          return 
-
-      else if (fl.lt.0.) then   !Orient the search so that f(xl) < 0.
-                                                                                                
+      else if (fl.lt.0.) then   !Orient the search so that f(xl) < 0.                                                                   
          xl=x1
          xh=x2 
-
       else
-
         xh=x1
         xl=x2 
-
       endif
-
       rtsafe = .5*(x1+x2)       !Initialize the guess for root,
       dxold  = abs(x2-x1)       !the stepsize before last,
       dx     = dxold            ! and the last step.
-
-
       call funcd(rtsafe,f,df,C0,C1,C2)
-
       DO 11 j=1,MAXIT           !Loop over allowed iterations. 
-
-         
-         !print*, 'iteration:', j
-         !print*, rtsafe
-
 
          if (((rtsafe-xh)*df-f)*((rtsafe-xl)*df-f).gt.0.   ! Bisect if Newton out of range
@@ -182 +134 @@
              dx=0.5*(xh-xl) 
              rtsafe=xl+dx 
-
              if (xl.eq.rtsafe) return                   !Change in root is negligible. Newton step acceptable. Take it.
-
          else 
-
             dxold=dx
             dx=f/df 
             temp=rtsafe
-
             rtsafe=rtsafe-dx
-
             if(temp.eq.rtsafe) return 
-
          endif
 
-
         if(abs(dx).lt.xacc) return !Convergence criterion. The one new function evaluation per iteration. Maintain the bracket on the root.
-
          call funcd(rtsafe,f,df,C0,C1,C2)
-
         if(f.lt.0.) then
          xl=rtsafe 
@@ -207 +150 @@
          xh=rtsafe 
         endif
-
-
 11    ENDDO
 
-      write(*,*) 'rtsafe exceeding maximum iterations'
-
+      rtsafe=0d0
+      write(*,*) 'rtsafe exceeding maximum iterations,Tcm=',rtsafe
       return 
 
@@ -219 +160 @@
 
 c********************************************************************************
-
       subroutine classical(x,f,df,C0,C1,C2)
-             
 c     Function to give as input to RTSAFE (NEWTON-RAPHOEN)
-
-
 c********************************************************************************
 
@@ -231 +168 @@
       DOUBLE PRECISION   x
       DOUBLE PRECISION  C0,C1,C2
-
       DOUBLE PRECISION f
       DOUBLE PRECISION  df
-      
-
-     
 
       f   = x + C0*exp(C1*x) - C2   ! start f
@@ -242 +175 @@
   
       return
-
       END  
 
@@ -250 +182 @@
 
 c********************************************************************************
-c defini la fonction eq 6 papier 2014 
+c defini la fonction eq 6 papier 2014  (Listowski et al., 2014)
       use tracer_mod, only: rho_ice_co2
       USE comcstfi_h
 
       implicit none
-
       include "microphys.h"
-c      include "microphysCO2.h"
-    
 
 c   arguments: INPUT
@@ -270 +199 @@
 c   local:
 c   ------
-
-
       DOUBLE PRECISION Cpatm,Cpn2,Cpco2
       DOUBLE PRECISION psat, xinf, pco2
@@ -298 +225 @@
 
 c     Equilibirum pressure over a flat surface
-
       psat = 1.382 * 1.00e12 * exp(-3182.48/dble(T))  ! (Pa)
-
 c     Compute transport coefficient
-
       pco2 = psat * dble(S)
-
 c     Latent heat of sublimation if CO2  co2 (J.kg-1)
 c     version Azreg_Ainou (J/kg) :
-
       l0=595594.      
       l1=903.111     
@@ -313 +235 @@
       l3=0.0528288 
       l4=-0.000103183
-  
       Lsub = l0 + l1 * dble(T) + l2 * dble(T)**2 + l3 * 
      &     dble(T)**3 + l4 * dble(T)**4 ! J/kg
-
 c     atmospheric density
-
       rhoatm = dble(P*Matm)/(rgp*dble(T))   ! g.m-3
       rhoatm = rhoatm * 1.00e-3 !kg.m-3
-
       call  KthMixNEW(kmix,T,pco2/dble(P),rhoatm) ! compute thermal cond of mixture co2/N2
       call  Diffcoeff(P, T, Dv)  
@@ -328 +246 @@
 
 c     ----- FS correction for Diff
-
       vthco2  = sqrt(8d0*kbz*dble(T)/(dble(pi) * mco2/nav)) ! units OK: m.s-1
-
       knudsen = 3*Dv / (vthco2 * rc) 
-
       lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black)
-      
       Dv      = Dv / (1. + lambda * knudsen)
-        
 c     ----- FS correction for Kth 
-
       vthatm = sqrt(8d0*kbz*dble(T)/(pi * 1.00e-3*dble(Matm)/nav)) ! Matm/nav = mass of "air molecule" in G , *1e-3 --> kg 
-
       Cpatm = Cpco2 * pco2/dble(P) + Cpn2 * (1d0 - pco2/dble(P)) !J.kg-1.K-1
-
       lpmt = 3 * kmix / (rhoatm * vthatm * (Cpatm - 0.5*rgp/
      &     (dble(Matm)*1.00e-3))) ! mean free path related to heat transfer
-
       knudsen = lpmt / rc
-
       lambda  = (1.333+0.71/knudsen) / (1.+1./knudsen) ! pas adaptée, Dahneke 1983? en fait si (Monschick&Black)
-
       kmix    = kmix /  (1. + lambda * knudsen)
-
 c     --------------------- ASSIGN coeff values for FUNCTION
-
       xinf = dble(S) * psat / dble(P)
-
       Ak = exp(2d0*sigco2*mco2/(rgp* dble(rho_ice_co2*T* rc) ))
-  
       C0 = mco2*Dv*psat*Lsub/(rgp*dble(T)*kmix)*Ak*exp(-Lsub*mco2/
      &     (rgp*dble(T)))
- 
       C1 = Lsub*mco2/(rgp*dble(T)**2)
-
       C2 = dble(T) + dble(P)*mco2*Dv*Lsub*xinf/(kmix*rgp*dble(T))
-     
       RETURN
       END
@@ -370 +270 @@
 
 c======================================================================
-
       subroutine Diffcoeff(P, T, Diff)
-
 c     Compute diffusion coefficient CO2/N2
 c     cited in Ilona's lecture - from Reid et al. 1987
 c======================================================================
-
-
        IMPLICIT NONE
 
        include "microphys.h"
-c       include "microphysCO2.h"
 
 c      arguments
@@ -425 +320 @@
 c======================================================================
 
-
        implicit none
        
        include "microphys.h"
-c       include "microphysCO2.h"
-
 c      arguments
 c     -----------
@@ -457 +349 @@
          DOUBLE PRECISION  kco2, kn2
 
-
       x1 = x
       x2 = 1d0 - x
 
-
       M1 = mco2
       M2 = mn2
 
-
       Tc1 =  304.1282 !(Scalabrin et al. 2006)
       Tc2 =  126.192  ! (Lemmon & Jacobsen 2003)
@@ -471 +360 @@
       Pc1 =  73.773   ! (bars)
       Pc2 =  33.958   ! (bars)
-  
-  
+    
       Gamma1 = 210.*(Tc1*M1**(3.)/Pc1**(4.))**(1./6.)
       Gamma2 = 210.*(Tc2*M2**(3.)/Pc2**(4.))**(1./6.)
 
-
 c Translational conductivities
-
-
 
       lambda_trans1 = ( exp(0.0464 * T/Tc1) - exp(-0.2412 * T/Tc1) )
@@ -487 +372 @@
      &                                                          /Gamma2
       
-         
 c     Coefficient of Mason and Saxena
-
-
       epsilon = 1.
-
 
       A11 = 1.
@@ -498 +379 @@
       A22 = 1.
 
-
       A12 = epsilon * (1. + sqrt(lambda_trans1/lambda_trans2)*
      &                    (M1/M2)**(1./4.))**(2.) / sqrt(8*(1.+ M1/M2))
@@ -505 +385 @@
      &                    (M2/M1)**(1./4.))**(2.) / sqrt(8*(1.+ M2/M1))
 
-
 c     INDIVIDUAL COND.
 
@@ -511 +390 @@
          call KthN2LemJac(kn2,T,rho)
 
-
 c     MIXTURE COND.
-
         Kthmix = kco2*x1 /(x1*A11 + x2*A12) + kn2*x2 /(x1*A21 + x2*A22)
         Kthmix = Kthmix*1e-3   ! from mW.m-1.K-1 to  W.m-1.K-1
-
         RETURN
 
         END
 
-
-c======================================================================
-
+c======================================================================
          subroutine KthN2LemJac(kthn2,T,rho)
-
 c        Compute thermal cond of N2 (Lemmon and Jacobsen, 2003)
 cWITH viscosity
@@ -565 +438 @@
 
         DOUBLE PRECISION k1, k2
-
 
          N1 = 1.511d0
@@ -612 +484 @@
          gamma9 = 1.
 
-
-
-
-c----------------------------------------------------------------------  
-         
+c----------------------------------------------------------------------           
          call viscoN2(T,visco)  !! v given in microPa.s
-
-       
+      
          Tc   = 126.192d0
          rhoc = 11.1839  * 1000 * mn2   !!!from mol.dm-3 to kg.m-3
@@ -626 +493 @@
          delta = rho/rhoc 
 
-
-
-         k1 =  N1 * visco + N2 * tau**t2 + N3 * tau**t3  !!! mW m-1 K-1
-      
-
+         k1 =  N1 * visco + N2 * tau**t2 + N3 * tau**t3  !!! mW m-1 K-1     
 c--------- residual thermal conductivity
-
-
 
          k2 = N4 * tau**t4 * delta**d4 * exp(-gamma4*delta**l4)         
@@ -642 +503 @@
      &  +     N9 * tau**t9 * delta**d9 * exp(-gamma9*delta**l9) 
 
-
          kthn2 = k1 + k2
-
 
          RETURN
@@ -744 +603 @@
       DOUBLE PRECISION h1,h2,h3,h4,h5,h6,h7,h8,h9,h10
       DOUBLE PRECISION n1,n2,n3,n4,n5,n6,n7,n8,n9,n10
-
-      
 
       Tc   = 304.1282   !(K) 
@@ -763 +620 @@
       g10 = 5.5
 
-
       h1 = 1.
       h2 = 5.
@@ -786 +642 @@
       n10 = 0.00433043347
 
-
-
       Tr   = T/Tc
       rhor = rho/rhoc
-
-
 
       k1 = n1*Tr**(g1)*rhor**(h1) + n2*Tr**(g2)*rhor**(h2) 
@@ -802 +654 @@
     
       k2  = exp(-5.*rhor**(2.)) * k2
-        
-        
+                
       kthco2 = (k1 + k2) *  Lambdac   ! mW
 
-
       RETURN
 

trunk/LMDZ.MARS/libf/phymars/microphys.h

@@ -2 +2 @@
 ! INCLUDE 'microphys.h'
 ! Parameters and physical constants used by the microphysal scheme;
+! Parameters for CO2 microphysics are also in this file
 !-----------------------------------------------------------------------
 
@@ -75 +76 @@
 !       (initialized in improvedCO2clouds.F)
 !    bachnar 2016 value :0.78 
-!old value 0.952
-      REAL, parameter :: mtetaco2 = 0.952
+!old value 0.95
+      REAL, parameter :: mtetaco2 = 0.95
 !     Volume of a co2 molecule (m3)
        DOUBLE PRECISION :: vo1co2

trunk/LMDZ.MARS/libf/phymars/nucleaCO2.F

@@ -15 +15 @@
 *     Adapted for the LMD/GCM by J.-B. Madeleine      *
 *     (October 2011)                                  *
-*     Optimisation by A. Spiga (February 2012)        *  
+*     Optimisation by A. Spiga (February 2012)        * 
+*   CO2 nucleation routine dev. by Constantino        *
+*     Listowski and Joachim Audouard (2016-2017),     *
+*     adapted from the water ice nucleation           *
 *******************************************************
  ! nucrate = output
-      ! nucrate_h2o en sortie aussi : 
+ ! nucrate_h2o en sortie aussi : 
 !nucleation sur dust et h2o separement ici
-!#include "tracer.h"
 #include "microphys.h"
-c#include "microphysCO2.h"
 
 c     Inputs
       DOUBLE PRECISION pco2,sat,vo2co2
       DOUBLE PRECISION n_ccn(nbinco2_cld), n_ccn_h2oice(nbinco2_cld)
-      REAL temp
-
+      REAL temp !temperature
 c     Output
-   !   DOUBLE PRECISION nucrate(nbinco2_cld)
       DOUBLE PRECISION nucrate(nbinco2_cld)
       DOUBLE PRECISION nucrate_h2oice(nbinco2_cld) ! h2o as substrate
-
       double precision rad_h2oice(nbinco2_cld) ! h2o ice grid (as substrate)
 
 c     Local variables
       DOUBLE PRECISION nco2
-c      DOUBLE PRECISION sigco2      ! Water-ice/air surface tension  (N.m)
-c      external sigco2
       DOUBLE PRECISION rstar    ! Radius of the critical germ (m)
       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 fshapeco2   ! function defined at the end of the file
-      DOUBLE PRECISION deltaf
-
-c     Ratio rstar/radius of the nucleating dust particle
-c     double precision xratio
-      
+      DOUBLE PRECISION deltaf      
       double precision mtetalocal,mtetalocalh ! local mteta in double precision
-
       double precision fshapeco2simple,zefshapeco2
-
-
       integer i
-      
-      LOGICAL firstcall
-      DATA firstcall/.true./
-      SAVE firstcall
-
 c     *************************************************
 

trunk/LMDZ.MARS/libf/phymars/physiq_mod.F

@@ -338 +338 @@
       REAL Mccntot(ngrid)       ! Total mass of ccn (kg/m2)
       REAL rave(ngrid)          ! Mean water ice effective radius (m)
-      REAL raveco2(ngrid)       ! Mean co2 ice effective radius (m)
       REAL opTES(ngrid,nlayer)  ! abs optical depth at 825 cm-1
       REAL tauTES(ngrid)        ! column optical depth at 825 cm-1
@@ -1218 +1217 @@
 
          
-         IF (co2clouds) THEN
+         IF (co2clouds ) THEN
             
             call co2cloud(ngrid,nlayer,ptimestep,
@@ -1225 +1224 @@
      &           nq,tau,tauscaling,rdust,rice,riceco2,nuice,
      &           rsedcloudco2,rhocloudco2,
-     &           rsedcloud,rhocloud,zzlev,zdqssed_co2)
+     &           rsedcloud,rhocloud,zzlev,zdqssed_co2,
+     &           pdu,pu)
                 
 
 c Temperature variation due to latent heat release
-            if (activice) then  !Maybe create activice_co2 ?
+c            if (activice) then  !Maybe create activice_co2 ?
                pdt(1:ngrid,1:nlayer) = 
      &              pdt(1:ngrid,1:nlayer) + 
-     &              zdtcloudco2(1:ngrid,1:nlayer)
-            endif
+     &              zdtcloudco2(1:ngrid,1:nlayer)! --> in newcondens
+c            endif
             
 
@@ -1240 +1240 @@
 ! This is due to single precision rounding problems
               
-              pdq(1:ngrid,1:nlayer,igcm_co2) = 
-     &             pdq(1:ngrid,1:nlayer,igcm_co2) + 
-     &             zdqcloudco2(1:ngrid,1:nlayer,igcm_co2)
-              pdq(1:ngrid,1:nlayer,igcm_co2_ice) = 
-     &             pdq(1:ngrid,1:nlayer,igcm_co2_ice) + 
-     &             zdqcloudco2(1:ngrid,1:nlayer,igcm_co2_ice)
-              pdq(1:ngrid,1:nlayer,igcm_ccnco2_mass) = 
-     &             pdq(1:ngrid,1:nlayer,igcm_ccnco2_mass) + 
-     &             zdqcloudco2(1:ngrid,1:nlayer,igcm_ccnco2_mass)
-             pdq(1:ngrid,1:nlayer,igcm_ccnco2_number) = 
-     &             pdq(1:ngrid,1:nlayer,igcm_ccnco2_number) + 
-     &             zdqcloudco2(1:ngrid,1:nlayer,igcm_ccnco2_number)
-             !Update water ice clouds values as well
-             pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = 
-     &            pdq(1:ngrid,1:nlayer,igcm_h2o_ice) + 
-     &            zdqcloudco2(1:ngrid,1:nlayer,igcm_h2o_ice)
-             
-! increment dust and ccn masses and numbers
-! We need to check that we have Nccn & Ndust > 0
-! This is due to single precision rounding problems
-             pdq(1:ngrid,1:nlayer,igcm_ccn_mass) = 
-     &            pdq(1:ngrid,1:nlayer,igcm_ccn_mass) + 
-     &            zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_mass)
-             pdq(1:ngrid,1:nlayer,igcm_ccn_number) = 
-     &            pdq(1:ngrid,1:nlayer,igcm_ccn_number) + 
-     &            zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
-             where (pq(:,:,igcm_ccn_mass) + 
-     &            ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
-             pdq(:,:,igcm_ccn_mass) = 
-     &            - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
-             pdq(:,:,igcm_ccn_number) = 
-     &            - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
+! increment dust tracers tendancies
+               pdq(1:ngrid,1:nlayer,igcm_dust_mass) = 
+     &              pdq(1:ngrid,1:nlayer,igcm_dust_mass) + 
+     &              zdqcloudco2(1:ngrid,1:nlayer,igcm_dust_mass)          
+               pdq(1:ngrid,1:nlayer,igcm_dust_number) = 
+     &              pdq(1:ngrid,1:nlayer,igcm_dust_number) + 
+     &              zdqcloudco2(1:ngrid,1:nlayer,igcm_dust_number)
+               pdq(1:ngrid,1:nlayer,igcm_co2) = 
+     &              pdq(1:ngrid,1:nlayer,igcm_co2) + 
+     &              zdqcloudco2(1:ngrid,1:nlayer,igcm_co2)
+               pdq(1:ngrid,1:nlayer,igcm_co2_ice) = 
+     &              pdq(1:ngrid,1:nlayer,igcm_co2_ice) + 
+     &              zdqcloudco2(1:ngrid,1:nlayer,igcm_co2_ice)
+               pdq(1:ngrid,1:nlayer,igcm_ccnco2_mass) = 
+     &              pdq(1:ngrid,1:nlayer,igcm_ccnco2_mass) + 
+     &              zdqcloudco2(1:ngrid,1:nlayer,igcm_ccnco2_mass)
+               pdq(1:ngrid,1:nlayer,igcm_ccnco2_number) = 
+     &              pdq(1:ngrid,1:nlayer,igcm_ccnco2_number) + 
+     &              zdqcloudco2(1:ngrid,1:nlayer,igcm_ccnco2_number)
+!Update water ice clouds values as well
+             if (co2useh2o) then
+                pdq(1:ngrid,1:nlayer,igcm_h2o_ice) = 
+     &               pdq(1:ngrid,1:nlayer,igcm_h2o_ice) + 
+     &               zdqcloudco2(1:ngrid,1:nlayer,igcm_h2o_ice)
+                pdq(1:ngrid,1:nlayer,igcm_ccn_mass) = 
+     &               pdq(1:ngrid,1:nlayer,igcm_ccn_mass) + 
+     &               zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_mass)
+                pdq(1:ngrid,1:nlayer,igcm_ccn_number) = 
+     &               pdq(1:ngrid,1:nlayer,igcm_ccn_number) + 
+     &               zdqcloudco2(1:ngrid,1:nlayer,igcm_ccn_number)
+                where (pq(:,:,igcm_ccn_mass) + 
+     &               ptimestep*pdq(:,:,igcm_ccn_mass) < 0.)
+                pdq(:,:,igcm_ccn_mass) = 
+     &               - pq(:,:,igcm_ccn_mass)/ptimestep + 1.e-30
+                pdq(:,:,igcm_ccn_number) = 
+     &               - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
              end where
              where (pq(:,:,igcm_ccn_number) + 
@@ -1279 +1283 @@
              pdq(:,:,igcm_ccn_number) = 
      &            - pq(:,:,igcm_ccn_number)/ptimestep + 1.e-30
-             end where
+          end where
+       endif
 c     Negative values?
              where (pq(:,:,igcm_ccnco2_mass) + 
@@ -1296 +1301 @@
           end where
        
-! increment dust tracers tendancies
-          pdq(1:ngrid,1:nlayer,igcm_dust_mass) = 
-     &         pdq(1:ngrid,1:nlayer,igcm_dust_mass) + 
-     &         zdqcloudco2(1:ngrid,1:nlayer,igcm_dust_mass)          
-          pdq(1:ngrid,1:nlayer,igcm_dust_number) = 
-     &         pdq(1:ngrid,1:nlayer,igcm_dust_number) + 
-     &         zdqcloudco2(1:ngrid,1:nlayer,igcm_dust_number)
 c     Negative values?
           where (pq(:,:,igcm_dust_mass) + 
@@ -1537 +1535 @@
      $              co2ice,albedo,emis,
      $              zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc,
-     $              fluxsurf_sw,zls) 
+     $              fluxsurf_sw,zls)!,
+c     &              zzlev,zdqssed_co2,zdqcloudco2,
+c     &              zdtcloudco2)
 
          DO l=1,nlayer
@@ -2260 +2260 @@
       !state  real  RICE      ikj   misc  1  -  h  "RICE"      "ICE RADIUS"                      "m"
       comm_RICE(1:ngrid,1:nlayer) = rice(1:ngrid,1:nlayer)
-      if (co2clouds) then
-        comm_RICE(1:ngrid,1:nlayer) = riceco2(1:ngrid,1:nlayer)
-      endif
    
       !! calculate sensible heat flux in W/m2 for outputs

trunk/LMDZ.MARS/libf/phymars/updaterad.F90

@@ -10 +10 @@
 ! So, subroutines are designed for scalar values instead of tables
 
-
 ! T. Navarro, June 2012
 ! CO2 clouds added 09/16 by J. Audouard
@@ -19 +18 @@
 real, parameter :: r3icemin = 1.e-30   ! ie ricemin  = 0.0001 microns
 real, parameter :: ricemin  = 1.e-10
-
 real, parameter :: r3icemax = 125.e-12 ! ie ricemax  = 500 microns
 real, parameter :: ricemax  = 500.e-6 
@@ -25 +23 @@
 double precision, parameter :: r3iceco2min = 1.e-30
 double precision, parameter :: riceco2min  = 1.e-10
-
 double precision, parameter :: r3iceco2max = 125.e-12
 double precision, parameter :: riceco2max  = 500.e-6
-
 
 real, parameter :: qice_threshold  = 1.e-15 ! 1.e-10
@@ -108 +104 @@
 USE comcstfi_h, only:  pi
 implicit none
-!By CL and JA 09/16
+!CO2 clouds parameter update by CL and JA 09/16
 
 DOUBLE PRECISION, intent(in)  :: qice,qccn,nccn