Index: trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/lect_start_archive.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/lect_start_archive.F	(revision 1708)
+++ trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/lect_start_archive.F	(revision 1711)
@@ -2,5 +2,5 @@
      &     date,tsurf,tsoil,emis,q2,
      &     t,ucov,vcov,ps,co2ice,h,phisold_newgrid,
-     &     q,qsurf,tauscaling,surfith,nid)
+     &     q,qsurf,tauscaling,totcloudfrac,surfith,nid)
 c=======================================================================
 c
@@ -49,4 +49,5 @@
       REAL,INTENT(OUT) :: q2(ngrid,nlayer+1),qsurf(ngrid,nqtot)
       REAL,INTENT(OUT) :: tauscaling(ngrid) ! dust conversion factor
+      REAL,INTENT(OUT) :: totcloudfrac(ngrid) ! sub grid cloud fraction
       REAL,INTENT(OUT) :: phisold_newgrid(iip1,jjp1)
       REAL,INTENT(OUT) :: t(iip1,jjp1,llm)
@@ -130,4 +131,5 @@
       REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
       real tauscalingS(iip1,jjp1)
+      real totcloudfracS(iip1,jjp1)
 
       real ptotal, co2icetotal
@@ -158,4 +160,5 @@
       real, dimension(:,:,:,:), allocatable :: qold
       real, dimension(:,:), allocatable :: tauscalingold
+      real, dimension(:,:), allocatable :: totcloudfracold
 
       real tab_cntrl(100)
@@ -358,4 +361,5 @@
       allocate(qsurfold(imold+1,jmold+1,nqtot))
       allocate(tauscalingold(imold+1,jmold+1))
+      allocate(totcloudfracold(imold+1,jmold+1))
 
       allocate(var (imold+1,jmold+1,llm))
@@ -683,5 +687,5 @@
 
 c-----------------------------------------------------------------------
-c 5.1 Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], q2surf, tauscaling)
+c 5.1 Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], q2surf, tauscaling,totcloudfrac)
 c-----------------------------------------------------------------------
  
@@ -777,4 +781,23 @@
         IF (ierr .NE. NF_NOERR) THEN
            PRINT*, "lect_start_archive: Failed loading <tauscaling>"
+           PRINT*, NF_STRERROR(ierr)
+           CALL abort
+        ENDIF
+      ENDIF
+      
+      ierr = NF_INQ_VARID (nid, "totcloudfrac", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: <totcloudfrac> not in file"
+         totcloudfracold(:,:) = 1
+      ELSE
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
+     &                            totcloudfracold)
+#else
+        ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
+     &                          totcloudfracold)
+#endif
+        IF (ierr .NE. NF_NOERR) THEN
+           PRINT*, "lect_start_archive: Failed loading <totcloudfrac>"
            PRINT*, NF_STRERROR(ierr)
            CALL abort
@@ -1063,4 +1086,10 @@
      &                   1,rlonuold,rlatvold,rlonu,rlatv)
       call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tauscalings,tauscaling)
+
+c Sub grid cloud fraction
+      call interp_horiz (totcloudfracold,totcloudfracs,imold,jmold,iim,
+     &                   jjm,1,rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,totcloudfracs,totcloudfrac)
+      
 c     write(46,*) 'emis',emis
 c-----------------------------------------------------------------------
@@ -1422,4 +1451,5 @@
       deallocate(qsurfold)
       deallocate(tauscalingold)
+      deallocate(totcloudfracold)
       deallocate(var,varp1)
 
Index: trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/newstart.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/newstart.F	(revision 1708)
+++ trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/newstart.F	(revision 1711)
@@ -41,5 +41,5 @@
       USE temps_mod, ONLY: day_ini,hour_ini
       USE ener_mod, ONLY: etot0,ptot0,ztot0,stot0,ang0
-
+      
       implicit none
 
@@ -168,4 +168,6 @@
       real,allocatable :: coefvmr(:)  ! Correction coefficient when changing composition
       integer :: iloc(1), iqmax
+! sub-grid cloud fraction
+      real,allocatable :: totcloudfrac(ngridmx)
 
 c sortie visu pour les champs dynamiques
@@ -419,5 +421,5 @@
      &   date,tsurf,tsoil,emis,q2,
      &   t,ucov,vcov,ps,co2ice,teta,phisold_newgrid,q,qsurf,
-     &   tauscaling,surfith,nid)
+     &   tauscaling,totcloudfrac,surfith,nid)
         write(*,*) "OK, read start_archive file"
 	! copy soil thermal inertia
@@ -440,5 +442,6 @@
         CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,ngridmx,llm,nqtot,
      .        day_ini,time,
-     .        tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling)
+     .        tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling,
+     .        totcloudfrac)
         
         ! copy albedo and soil thermal inertia
@@ -1657,5 +1660,6 @@
       call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
      .              dtphys,hour_ini,
-     .              tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling)
+     .              tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling,
+     .              totcloudfrac)
 
 c=======================================================================
Index: trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/start2archive.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/start2archive.F	(revision 1708)
+++ trunk/LMDZ.MARS/libf/dynlonlat_phylonlat/phymars/start2archive.F	(revision 1711)
@@ -63,4 +63,5 @@
       REAL co2ice(ngridmx)	! CO2 ice layer
       REAL tauscaling(ngridmx) ! dust conversion factor
+      REAL totcloudfrac(ngridmx) ! sub-grid cloud fraction
       REAL q2(ngridmx,llm+1)
       REAL emis(ngridmx)
@@ -78,4 +79,5 @@
       REAL co2iceS(ip1jmp1)
       REAL tauscalingS(ip1jmp1)
+      REAL totcloudfracS(ip1jmp1)
       REAL q2S(ip1jmp1,llm+1)
       REAL,ALLOCATABLE :: qsurfS(:,:)
@@ -145,5 +147,5 @@
       CALL phyetat0 (fichnom,0,Lmodif,nsoilmx,ngridmx,llm,nqtot,
      &      day_ini_fi,timefi,tsurf,tsoil,emis,q2,qsurf,co2ice,
-     &      tauscaling)
+     &      tauscaling,totcloudfrac)
 
        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
@@ -241,4 +243,5 @@
 c qsurf --> qsurfS
 c tauscaling --> tauscalingS
+c totcloudfrac --> totcloudfracS
 c
 c-----------------------------------------------------------------------
@@ -253,5 +256,6 @@
       call gr_fi_dyn(nqtot,ngridmx,iip1,jjp1,qsurf,qsurfS)
       call gr_fi_dyn(1,ngridmx,iip1,jjp1,tauscaling,tauscalingS)
-
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,totcloudfrac,totcloudfracS)
+      
 c=======================================================================
 c Info pour controler
@@ -343,4 +347,6 @@
       call write_archive(nid,ntime,'tauscaling',
      &  'dust conversion factor',' ',2,tauscalingS)
+      call write_archive(nid,ntime,'totcloudfrac',
+     &  'sub grid cloud fraction',' ',2,totcloudfracS)
       call write_archive(nid,ntime,'emis','grd emis',' ',2,emisS)
       call write_archive(nid,ntime,'ps','Psurf','Pa',2,ps)
Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/lect_start_archive.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/lect_start_archive.F	(revision 1708)
+++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/lect_start_archive.F	(revision 1711)
@@ -2,5 +2,5 @@
      &     date,tsurf,tsoil,emis,q2,
      &     t,ucov,vcov,ps,co2ice,h,phisold_newgrid,
-     &     q,qsurf,tauscaling,surfith,nid)
+     &     q,qsurf,tauscaling,totcloudfrac,surfith,nid)
 c=======================================================================
 c
@@ -49,4 +49,5 @@
       REAL,INTENT(OUT) :: q2(ngrid,nlayer+1),qsurf(ngrid,nqtot)
       REAL,INTENT(OUT) :: tauscaling(ngrid) ! dust conversion factor
+      REAL,INTENT(OUT) :: totcloudfrac(ngrid) ! sub grid cloud fraction
       REAL,INTENT(OUT) :: phisold_newgrid(iip1,jjp1)
       REAL,INTENT(OUT) :: t(iip1,jjp1,llm)
@@ -130,4 +131,5 @@
       REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
       real tauscalingS(iip1,jjp1)
+      real totcloudfracS(iip1,jjp1)
 
       real ptotal, co2icetotal
@@ -158,4 +160,5 @@
       real, dimension(:,:,:,:), allocatable :: qold
       real, dimension(:,:), allocatable :: tauscalingold
+      real, dimension(:,:), allocatable :: totcloudfracold
 
       real tab_cntrl(100)
@@ -358,4 +361,5 @@
       allocate(qsurfold(imold+1,jmold+1,nqtot))
       allocate(tauscalingold(imold+1,jmold+1))
+      allocate(totcloudfracold(imold+1,jmold+1))
 
       allocate(var (imold+1,jmold+1,llm))
@@ -683,5 +687,5 @@
 
 c-----------------------------------------------------------------------
-c 5.1 Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], q2surf, tauscaling)
+c 5.1 Lecture des champs 2D (co2ice, emis,ps,tsurf,Tg[10], q2surf, tauscaling,totcloudfrac)
 c-----------------------------------------------------------------------
  
@@ -777,4 +781,23 @@
         IF (ierr .NE. NF_NOERR) THEN
            PRINT*, "lect_start_archive: Failed loading <tauscaling>"
+           PRINT*, NF_STRERROR(ierr)
+           CALL abort
+        ENDIF
+      ENDIF
+
+      ierr = NF_INQ_VARID (nid, "totcloudfrac", nvarid)
+      IF (ierr .NE. NF_NOERR) THEN
+         PRINT*, "lect_start_archive: <totcloudfrac> not in file"
+         totcloudfracold(:,:) = 1
+      ELSE
+#ifdef NC_DOUBLE
+        ierr = NF_GET_VARA_DOUBLE(nid, nvarid,start,count,
+     &                            totcloudfracold)
+#else
+        ierr = NF_GET_VARA_REAL(nid, nvarid,start,count,
+     &                          totcloudfracold)
+#endif
+        IF (ierr .NE. NF_NOERR) THEN
+           PRINT*, "lect_start_archive: Failed loading <totcloudfrac>"
            PRINT*, NF_STRERROR(ierr)
            CALL abort
@@ -1063,4 +1086,10 @@
      &                   1,rlonuold,rlatvold,rlonu,rlatv)
       call gr_dyn_fi (1,iim+1,jjm+1,ngrid,tauscalings,tauscaling)
+      
+c Sub grid cloud fraction
+      call interp_horiz (totcloudfracold,totcloudfracs,imold,jmold,iim,
+     &                   jjm,1,rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (1,iim+1,jjm+1,ngrid,totcloudfracs,totcloudfrac)
+      
 c     write(46,*) 'emis',emis
 c-----------------------------------------------------------------------
@@ -1422,4 +1451,5 @@
       deallocate(qsurfold)
       deallocate(tauscalingold)
+      deallocate(totcloudfracold)
       deallocate(var,varp1)
 
Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F	(revision 1708)
+++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/newstart.F	(revision 1711)
@@ -170,4 +170,6 @@
       real,allocatable :: coefvmr(:)  ! Correction coefficient when changing composition
       integer :: iloc(1), iqmax
+! sub-grid cloud fraction
+      real :: totcloudfrac(ngridmx)
 
 c sortie visu pour les champs dynamiques
@@ -402,5 +404,5 @@
      &   date,tsurf,tsoil,emis,q2,
      &   t,ucov,vcov,ps,co2ice,teta,phisold_newgrid,q,qsurf,
-     &   tauscaling,surfith,nid)
+     &   tauscaling,totcloudfrac,surfith,nid)
         write(*,*) "OK, read start_archive file"
 	! copy soil thermal inertia
@@ -423,5 +425,5 @@
         CALL phyetat0 (fichnom,tab0,Lmodif,nsoilmx,ngridmx,llm,nqtot,
      .        day_ini,time,
-     .        tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling)
+     .        tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling,totcloudfrac)
         
         ! copy albedo and soil thermal inertia
@@ -1640,5 +1642,6 @@
       call physdem1("restartfi.nc",nsoilmx,ngridmx,llm,nqtot,
      .              dtphys,hour_ini,
-     .              tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling)
+     .              tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling,
+     .              totcloudfrac)
 
 c=======================================================================
Index: trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/start2archive.F
===================================================================
--- trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/start2archive.F	(revision 1708)
+++ trunk/LMDZ.MARS/libf/dynphy_lonlat/phymars/start2archive.F	(revision 1711)
@@ -64,4 +64,5 @@
       REAL co2ice(ngridmx)        ! CO2 ice layer
       REAL tauscaling(ngridmx) ! dust conversion factor
+      REAL totcloudfrac(ngridmx) ! sub-grid cloud fraction
       REAL q2(ngridmx,llm+1)
       REAL emis(ngridmx)
@@ -79,4 +80,5 @@
       REAL co2iceS(ip1jmp1)
       REAL tauscalingS(ip1jmp1)
+      REAL totcloudfracS(ip1jmp1)
       REAL q2S(ip1jmp1,llm+1)
       REAL,ALLOCATABLE :: qsurfS(:,:)
@@ -156,5 +158,5 @@
       CALL phyetat0 (fichnom,0,Lmodif,nsoilmx,ngridmx,llm,nqtot,
      &      day_ini_fi,timefi,tsurf,tsoil,emis,q2,qsurf,co2ice,
-     &      tauscaling)
+     &      tauscaling,totcloudfrac)
 
        ierr = NF_OPEN (fichnom, NF_NOWRITE,nid1)
@@ -244,4 +246,5 @@
 c qsurf --> qsurfS
 c tauscaling --> tauscalingS
+c totcloudfrac --> totcloudfracS
 c
 c-----------------------------------------------------------------------
@@ -256,4 +259,5 @@
       call gr_fi_dyn(nqtot,ngridmx,iip1,jjp1,qsurf,qsurfS)
       call gr_fi_dyn(1,ngridmx,iip1,jjp1,tauscaling,tauscalingS)
+      call gr_fi_dyn(1,ngridmx,iip1,jjp1,totcloudfrac,totcloudfracS)
 
 c=======================================================================
@@ -346,4 +350,6 @@
       call write_archive(nid,ntime,'tauscaling',
      &  'dust conversion factor',' ',2,tauscalingS)
+      call write_archive(nid,ntime,'totcloudfrac',
+     &  'sub grid cloud fraction',' ',2,totcloudfracS)
       call write_archive(nid,ntime,'emis','grd emis',' ',2,emisS)
       call write_archive(nid,ntime,'ps','Psurf','Pa',2,ps)
Index: trunk/LMDZ.MARS/libf/phymars/aeropacity.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/aeropacity.F	(revision 1708)
+++ 	(revision )
@@ -1,729 +1,0 @@
-      SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
-     &    pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
-     &    nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d)
-                                                   
-! to use  'getin'
-      USE ioipsl_getincom, only: getin
-      use tracer_mod, only: noms, igcm_h2o_ice, igcm_dust_mass,
-     &                      igcm_dust_submicron, rho_dust, rho_ice,
-     &                      nqdust
-      use geometry_mod, only: latitude ! grid point latitudes (rad)
-      use comgeomfi_h, only: sinlat ! sines of grid point latitudes
-#ifdef DUSTSTORM
-      use geometry_mod, only: longitude
-      use tracer_mod, only: r3n_q, ref_r0, igcm_dust_number
-#endif
-      use planete_h
-      USE comcstfi_h
-      use dimradmars_mod, only: naerkind, name_iaer, 
-     &            iaerdust,tauvis,
-     &            iaer_dust_conrath,iaer_dust_doubleq,
-     &            iaer_dust_submicron,iaer_h2o_ice
-       IMPLICIT NONE
-c=======================================================================
-c   subject:
-c   --------
-c   Computing aerosol optical depth in each gridbox.
-c
-c   author: F.Forget 
-c   ------
-c   update F. Montmessin (water ice scheme) 
-c      and S. Lebonnois (12/06/2003) compatibility dust/ice/chemistry
-c   update J.-B. Madeleine 2008-2009:
-c       - added 3D scattering by aerosols;
-c       - dustopacity transferred from physiq.F to callradite.F,
-c           and renamed into aeropacity.F;
-c   update E. Millour, march 2012:
-c         - reference pressure is now set to 610Pa (not 700Pa)
-c   
-c   input:
-c   ----- 
-c   ngrid             Number of gridpoint of horizontal grid
-c   nlayer            Number of layer
-c   nq                Number of tracer
-c   zday                  Date (time since Ls=0, in martian days)
-c   ls                Solar longitude (Ls) , radian
-c   pplay,pplev       pressure (Pa) in the middle and boundary of each layer
-c   pq                Dust mixing ratio (used if tracer =T and active=T).
-c   reffrad(ngrid,nlayer,naerkind)  Aerosol effective radius
-c   QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
-c   QREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths;
-c   omegaREFvis3d(ngrid,nlayer,naerkind) \ 3d single scat. albedo
-c   omegaREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths;
-c
-c   output:
-c   -------
-c   tauref       Prescribed mean column optical depth at 610 Pa 
-c   tau          Column total visible dust optical depth at each point
-c   aerosol      aerosol(ig,l,1) is the dust optical
-c                depth in layer l, grid point ig
-
-c
-c=======================================================================
-#include "callkeys.h"
-
-c-----------------------------------------------------------------------
-c
-c    Declarations :
-c    --------------
-c
-c    Input/Output
-c    ------------
-      INTEGER ngrid,nlayer,nq
-
-      REAL ls,zday,expfactor    
-      REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
-      REAL pq(ngrid,nlayer,nq)
-      REAL tauref(ngrid), tau(ngrid,naerkind)
-      REAL aerosol(ngrid,nlayer,naerkind)
-      REAL dsodust(ngrid,nlayer)
-      REAL reffrad(ngrid,nlayer,naerkind)
-      REAL nueffrad(ngrid,nlayer,naerkind)
-      REAL QREFvis3d(ngrid,nlayer,naerkind)
-      REAL QREFir3d(ngrid,nlayer,naerkind)
-      REAL omegaREFvis3d(ngrid,nlayer,naerkind)
-      REAL omegaREFir3d(ngrid,nlayer,naerkind)
-c
-c    Local variables :
-c    -----------------
-      INTEGER l,ig,iq,i,j
-      INTEGER iaer           ! Aerosol index
-      real topdust(ngrid)
-      real zlsconst, zp
-      real taueq,tauS,tauN
-c     Mean Qext(vis)/Qext(ir) profile
-      real msolsir(nlayer,naerkind)
-c     Mean Qext(ir)/Qabs(ir) profile
-      real mqextsqabs(nlayer,naerkind)
-c     Variables used when multiple particle sizes are used
-c       for dust or water ice particles in the radiative transfer
-c       (see callradite.F for more information).
-      REAL taudusttmp(ngrid)! Temporary dust opacity
-                               !   used before scaling
-      REAL tauscaling(ngrid) ! Scaling factor for qdust and Ndust
-      REAL taudustvis(ngrid) ! Dust opacity after scaling
-      REAL taudusttes(ngrid) ! Dust opacity at IR ref. wav. as
-                               !   "seen" by the GCM.
-      REAL taucloudvis(ngrid)! Cloud opacity at visible
-                               !   reference wavelength
-      REAL taucloudtes(ngrid)! Cloud opacity at infrared
-                               !   reference wavelength using
-                               !   Qabs instead of Qext
-                               !   (direct comparison with TES)
-      REAL topdust0(ngrid)
-
-#ifdef DUSTSTORM
-!! Local dust storms
-      logical localstorm        ! =true to create a local dust storm
-      real taulocref,ztoploc,radloc,lonloc,latloc  ! local dust storm parameters
-      real reffstorm, yeah
-      REAL ray(ngrid) ! distance from dust storm center
-      REAL tauuser(ngrid) ! opacity perturbation due to dust storm
-      REAL more_dust(ngrid,nlayer,2) ! Mass mixing ratio perturbation due to the dust storm
-      REAL int_factor(ngrid) ! useful factor to compute mmr perturbation
-      real l_top ! layer of the storm's top
-      REAL zalt(ngrid, nlayer) ! useful factor to compute l_top 
-#endif
-
-c   local saved variables
-c   ---------------------
-
-c     Level under which the dust mixing ratio is held constant
-c       when computing the dust opacity in each layer
-c       (this applies when doubleq and active are true)
-      INTEGER, PARAMETER :: cstdustlevel0 = 7
-      INTEGER, SAVE      :: cstdustlevel
-
-      LOGICAL,SAVE :: firstcall=.true.
-
-! indexes of water ice and dust tracers:
-      INTEGER,SAVE :: i_ice=0  ! water ice
-      real,parameter :: odpref=610. ! DOD reference pressure (Pa)
-      CHARACTER(LEN=20) :: txt ! to temporarly store text
-      CHARACTER(LEN=1) :: txt2 ! to temporarly store text
-! indexes of dust scatterers:
-      INTEGER,SAVE :: naerdust ! number of dust scatterers
-
-      tau(1:ngrid,1:naerkind)=0
-
-! identify tracers
-
-      IF (firstcall) THEN
-        ! identify scatterers that are dust
-        naerdust=0
-        DO iaer=1,naerkind
-          txt=name_iaer(iaer)
-          IF (txt(1:4).eq."dust") THEN
-            naerdust=naerdust+1
-            iaerdust(naerdust)=iaer
-          ENDIF
-        ENDDO
-        ! identify tracers which are dust
-        i=0
-        DO iq=1,nq
-          txt=noms(iq)
-          IF (txt(1:4).eq."dust") THEN
-          i=i+1
-          nqdust(i)=iq
-          ENDIF
-        ENDDO
-
-        IF (water.AND.activice) THEN
-          i_ice=igcm_h2o_ice
-          write(*,*) "aeropacity: i_ice=",i_ice
-        ENDIF
-
-c       typical profile of solsir and (1-w)^(-1):
-        msolsir(1:nlayer,1:naerkind)=0
-        mqextsqabs(1:nlayer,1:naerkind)=0
-        WRITE(*,*) "Typical profiles of Qext(vis)/Qext(IR)"
-        WRITE(*,*) "  and Qext(IR)/Qabs(IR):"
-        DO iaer = 1, naerkind ! Loop on aerosol kind
-          WRITE(*,*) "Aerosol # ",iaer
-          DO l=1,nlayer
-            DO ig=1,ngrid
-              msolsir(l,iaer)=msolsir(l,iaer)+
-     &              QREFvis3d(ig,l,iaer)/
-     &              QREFir3d(ig,l,iaer)
-              mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
-     &              (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
-            ENDDO
-            msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngrid)
-            mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngrid)
-          ENDDO
-          WRITE(*,*) "solsir: ",msolsir(:,iaer)
-          WRITE(*,*) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
-        ENDDO
-
-!       load value of tauvis from callphys.def (if given there,
-!       otherwise default value read from starfi.nc file will be used)
-        call getin("tauvis",tauvis)
-
-        IF (freedust) THEN
-          cstdustlevel = 1
-        ELSE
-          cstdustlevel = cstdustlevel0
-        ENDIF
-
-
-#ifndef DUSTSTORM
-        firstcall=.false.
-#endif
-
-      END IF
-
-c     Vertical column optical depth at "odpref" Pa 
-c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      IF(freedust) THEN
-         tauref(:) = 0. ! tauref is computed after, instead of being forced
-
-      ELSE IF(iaervar.eq.1) THEN 
-         do ig=1, ngrid
-          tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
-                                       ! or read in starfi
-        end do
-      ELSE IF (iaervar.eq.2) THEN   ! << "Viking" Scenario>>
-
-        tauref(1) = 0.7+.3*cos(ls+80.*pi/180.) ! like seen by VL1
-        do ig=2,ngrid
-          tauref(ig) = tauref(1)
-        end do
-
-      ELSE IF (iaervar.eq.3) THEN  ! << "MGS" scenario >>
-
-        taueq= 0.2 +(0.5-0.2) *(cos(0.5*(ls-4.363)))**14
-        tauS= 0.1 +(0.5-0.1)  *(cos(0.5*(ls-4.363)))**14
-        tauN = 0.1
-c	   if (peri_day.eq.150) then
-c	     tauS=0.1
-c	     tauN=0.1 +(0.5-0.1)  *(cos(0.5*(ls+pi-4.363)))**14
-c	     taueq= 0.2 +(0.5-0.2) *(cos(0.5*(ls+pi-4.363)))**14
-c           endif
-        do ig=1,ngrid
-          if (latitude(ig).ge.0) then
-          ! Northern hemisphere
-            tauref(ig)= tauN +
-     &      (taueq-tauN)*0.5*(1+tanh((45-latitude(ig)*180./pi)*6/60))
-          else
-          ! Southern hemisphere
-            tauref(ig)= tauS +
-     &      (taueq-tauS)*0.5*(1+tanh((45+latitude(ig)*180./pi)*6/60))
-          endif
-        enddo ! of do ig=1,ngrid
-      ELSE IF (iaervar.eq.5) THEN   ! << Escalier Scenario>>
-c         tauref(1) = 0.2
-c         if ((ls.ge.210.*pi/180.).and.(ls.le.330.*pi/180.))
-c    &                              tauref(1) = 2.5
-        tauref(1) = 2.5
-        if ((ls.ge.30.*pi/180.).and.(ls.le.150.*pi/180.))
-     &                              tauref(1) = .2
-
-        do ig=2,ngrid
-          tauref(ig) = tauref(1)
-        end do
-      ELSE IF ((iaervar.ge.6).and.(iaervar.le.8)) THEN
-      ! clim, cold or warm synthetic scenarios
-        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
-      ELSE IF ((iaervar.ge.24).and.(iaervar.le.32))
-     &     THEN  ! << MY... dust scenarios >>
-        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
-      ELSE IF ((iaervar.eq.4).or.
-     &         ((iaervar.ge.124).and.(iaervar.le.126))) THEN
-       ! "old" TES assimation dust scenario (values at 700Pa in files!)
-        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
-      ELSE
-        stop 'problem with iaervar in aeropacity.F'
-      ENDIF
-
-c -----------------------------------------------------------------
-c Computing the opacity in each layer
-c -----------------------------------------------------------------
-
-      DO iaer = 1, naerkind ! Loop on aerosol kind
-c     --------------------------------------------
-        aerkind: SELECT CASE (name_iaer(iaer))
-c==================================================================
-        CASE("dust_conrath") aerkind      ! Typical dust profile
-c==================================================================
-
-c       Altitude of the top of the dust layer
-c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        zlsconst=SIN(ls-2.76)
-        if (iddist.eq.1) then
-          do ig=1,ngrid
-             topdust(ig)=topdustref         ! constant dust layer top
-          end do
-
-        else if (iddist.eq.2) then          ! "Viking" scenario
-          do ig=1,ngrid
-            ! altitude of the top of the aerosol layer (km) at Ls=2.76rad:
-            ! in the Viking year scenario
-            topdust0(ig)=60. -22.*sinlat(ig)**2
-            topdust(ig)=topdust0(ig)+18.*zlsconst
-          end do
-
-        else if(iddist.eq.3) then         !"MGS" scenario
-          do ig=1,ngrid
-            topdust(ig)=60.+18.*zlsconst
-     &                -(32+18*zlsconst)*sin(latitude(ig))**4
-     &                 - 8*zlsconst*(sin(latitude(ig)))**5
-          end do
-        endif
-
-c       Optical depth in each layer :
-c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-        if(iddist.ge.1) then
-
-          expfactor=0.
-          DO l=1,nlayer
-            DO ig=1,ngrid
-c             Typical mixing ratio profile 
-              if(pplay(ig,l).gt.odpref
-     $                        /(988.**(topdust(ig)/70.))) then
-                zp=(odpref/pplay(ig,l))**(70./topdust(ig))
-                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
-              else    
-                expfactor=1.e-3
-              endif
-c             Vertical scaling function
-              aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) *
-     &          expfactor *
-     &          QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
-            ENDDO
-          ENDDO
-
-        else if(iddist.eq.0) then   
-c         old dust vertical distribution function (pollack90)
-          DO l=1,nlayer
-             DO ig=1,ngrid
-                zp=odpref/pplay(ig,l)
-                aerosol(ig,l,1)= tauref(ig)/odpref *
-     s           (pplev(ig,l)-pplev(ig,l+1))
-     s           *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 )
-             ENDDO
-          ENDDO
-        end if
-
-c==================================================================
-        CASE("dust_doubleq") aerkind! Two-moment scheme for dust
-c        (transport of mass and number mixing ratio)
-c==================================================================
-
-          DO l=1,nlayer
-            IF (l.LE.cstdustlevel) THEN
-c           Opacity in the first levels is held constant to 
-c             avoid unrealistic values due to constant lifting:
-              DO ig=1,ngrid
-                aerosol(ig,l,iaer) = 
-     &          (  0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
-     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
-     &          pq(ig,cstdustlevel,igcm_dust_mass) *
-     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
-                ! DENSITY SCALED OPACITY IN INFRARED:
-                dsodust(ig,l) =
-     &          (  0.75 * QREFir3d(ig,cstdustlevel,iaer) /
-     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
-     &          pq(ig,cstdustlevel,igcm_dust_mass)
-              ENDDO
-            ELSE
-              DO ig=1,ngrid
-                aerosol(ig,l,iaer) =
-     &          (  0.75 * QREFvis3d(ig,l,iaer) /
-     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
-     &          pq(ig,l,igcm_dust_mass) *
-     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
-                ! DENSITY SCALED OPACITY IN INFRARED:
-                dsodust(ig,l) =
-     &          (  0.75 * QREFir3d(ig,l,iaer) /
-     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
-     &          pq(ig,l,igcm_dust_mass)
-              ENDDO
-            ENDIF
-          ENDDO
-
-c==================================================================
-        CASE("dust_submicron") aerkind   ! Small dust population
-c==================================================================
-
-          DO l=1,nlayer
-            IF (l.LE.cstdustlevel) THEN
-c           Opacity in the first levels is held constant to 
-c             avoid unrealistic values due to constant lifting:
-              DO ig=1,ngrid
-                aerosol(ig,l,iaer) = 
-     &          (  0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
-     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
-     &          pq(ig,cstdustlevel,igcm_dust_submicron) *
-     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
-              ENDDO
-            ELSE
-              DO ig=1,ngrid
-                aerosol(ig,l,iaer) = 
-     &          (  0.75 * QREFvis3d(ig,l,iaer) /
-     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
-     &          pq(ig,l,igcm_dust_submicron) *
-     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
-              ENDDO
-            ENDIF
-          ENDDO
-
-c==================================================================
-        CASE("h2o_ice") aerkind             ! Water ice crystals
-c==================================================================
-
-c       1. Initialization
-        aerosol(1:ngrid,1:nlayer,iaer) = 0.
-        taucloudvis(1:ngrid) = 0.
-        taucloudtes(1:ngrid) = 0.
-c       2. Opacity calculation
-        DO ig=1, ngrid
-          DO l=1,nlayer
-            aerosol(ig,l,iaer) = max(1E-20,
-     &        (  0.75 * QREFvis3d(ig,l,iaer) /
-     &        ( rho_ice * reffrad(ig,l,iaer) )  ) *
-     &        pq(ig,l,i_ice) *
-     &        ( pplev(ig,l) - pplev(ig,l+1) ) / g
-     &                              )
-            taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
-            taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
-     &        QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
-     &        ( 1.E0 - omegaREFir3d(ig,l,iaer) )
-          ENDDO
-        ENDDO
-c       3. Outputs -- Now done in physiq.F
-!        IF (ngrid.NE.1) THEN
-!          CALL WRITEDIAGFI(ngrid,'tauVIS','tauext VIS refwvl',
-!     &      ' ',2,taucloudvis)
-!          CALL WRITEDIAGFI(ngrid,'tauTES','tauabs IR refwvl',
-!     &      ' ',2,taucloudtes)
-!          IF (callstats) THEN
-!            CALL wstats(ngrid,'tauVIS','tauext VIS refwvl',
-!     &        ' ',2,taucloudvis)
-!            CALL wstats(ngrid,'tauTES','tauabs IR refwvl',
-!     &        ' ',2,taucloudtes)
-!          ENDIF
-!        ELSE
-!         CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
-!        ENDIF
-c==================================================================
-        END SELECT aerkind
-c     -----------------------------------
-      ENDDO ! iaer (loop on aerosol kind)
-
-c -----------------------------------------------------------------
-c Rescaling each layer to reproduce the choosen (or assimilated)
-c   dust extinction opacity at visible reference wavelength, which
-c   is originally scaled to an equivalent odpref Pa pressure surface.
-c -----------------------------------------------------------------
-
-#ifdef DUSTSTORM
-c -----------------------------------------------------------------
-! Calculate reference opacity without perturbation
-c -----------------------------------------------------------------
-      IF (firstcall) THEN
-        DO iaer=1,naerdust
-          DO l=1,nlayer
-            DO ig=1,ngrid
-              tauref(ig) = tauref(ig) +
-     &                    aerosol(ig,l,iaerdust(iaer))
-            ENDDO
-          ENDDO
-        ENDDO
-        tauref(:) = tauref(:) * odpref / pplev(:,1)
-
-c--------------------------------------------------
-c Get parameters of the opacity perturbation
-c--------------------------------------------------
-        iaer=1  ! just change dust
-
-        write(*,*) "Add a local storm ?"
-        localstorm=.true. ! default value
-        call getin("localstorm",localstorm)
-        write(*,*) " localstorm = ",localstorm
-
-        IF (localstorm) THEN
-          WRITE(*,*) "********************"
-          WRITE(*,*) "ADDING A LOCAL STORM"
-          WRITE(*,*) "********************"
-
-          write(*,*) "ref opacity of local dust storm"
-              taulocref = 4.25 ! default value
-              call getin("taulocref",taulocref)
-              write(*,*) " taulocref = ",taulocref
-
-          write(*,*) "target altitude of local storm (km)"
-              ztoploc = 10.0 ! default value
-              call getin("ztoploc",ztoploc)
-              write(*,*) " ztoploc = ",ztoploc
-
-          write(*,*) "radius of dust storm (degree)"
-              radloc = 0.5 ! default value
-              call getin("radloc",radloc)
-              write(*,*) " radloc = ",radloc
-
-          write(*,*) "center longitude of storm (deg)"
-              lonloc = 25.0 ! default value
-              call getin("lonloc",lonloc)
-              write(*,*) " lonloc = ",lonloc
-
-          write(*,*) "center latitude of storm (deg)"
-              latloc = -2.5 ! default value
-              call getin("latloc",latloc)
-              write(*,*) " latloc = ",latloc
-        
-          write(*,*) "reff storm (mic) 0. for background"
-              reffstorm = 0.0 ! default value
-              call getin("reffstorm",reffstorm)
-              write(*,*) " reffstorm = ",reffstorm
-
-!! LOOP: modify opacity
-      DO ig=1,ngrid
-
-      !! distance to the center:
-      ray(ig)=SQRT((latitude(ig)*180./pi-latloc)**2 +
-     &          (longitude(ig)*180./pi -lonloc)**2)
-
-      !! transition factor for storm
-      !! factor is hardcoded -- increase it to steepen
-      yeah = (TANH(2.+(radloc-ray(ig))*10.)+1.)/2.
-
-      !! new opacity field
-      !! -- add an opacity set to taulocref
-      !! -- the additional reference opacity will
-      !!      thus be taulocref*odpref/pplev
-      tauuser(ig)=max( tauref(ig) * pplev(ig,1) /odpref , 
-     &          taulocref * yeah )
-
-      !! compute l_top
-          DO l=1,nlayer
-            zalt(ig,l) = LOG( pplev(ig,1)/pplev(ig,l) )
-     &                      / g / 44.01
-     &                    * 8.31 * 210.
-                IF (     (ztoploc .lt. zalt(ig,l)  )
-     &          .and. (ztoploc .gt. zalt(ig,l-1)) ) l_top=l-1
-          ENDDO
-
-     !! change reffrad if ever needed
-      IF (reffstorm .gt. 0.) THEN
-          DO l=1,nlayer
-             IF (l .lt. l_top+1) THEN
-                reffrad(ig,l,iaer) = max( reffrad(ig,l,iaer), reffstorm
-     &          * 1.e-6 * yeah )
-             ENDIF
-          ENDDO
-      ENDIF
-
-      ENDDO
-!! END LOOP
-
-      !! compute perturbation in each layer (equation 8 in Spiga et al. JGR 2013)
-      DO ig=1,ngrid
-          int_factor(ig)=0.
-          DO l=1,nlayer
-             IF (l .lt. l_top+1) THEN
-                      int_factor(ig) =
-     &                int_factor(ig) +
-     &          (  0.75 * QREFvis3d(ig,l,iaer) /
-     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
-     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
-             ENDIF
-          ENDDO
-          DO l=1, nlayer
-          !! Mass mixing ratio perturbation due to local dust storm in each layer
-          more_dust(ig,l,1)=
-     &                     (tauuser(ig)-(tauref(ig)
-     &                      * pplev(ig,1) /odpref)) /
-     &                      int_factor(ig)
-          more_dust(ig,l,2)=
-     &                     (tauuser(ig)-(tauref(ig) *
-     &                     pplev(ig,1) /odpref))
-     &                      / int_factor(ig) *
-     &                     ((ref_r0/reffrad(ig,l,iaer))**3)
-     &                      * r3n_q 
-          ENDDO
-      ENDDO
-
-      !! quantity of dust for each layer with the addition of the perturbation
-      DO l=1, l_top
-          pq(:,l,igcm_dust_mass)= pq(:,l,igcm_dust_mass)
-     .            + more_dust(:,l,1)
-          pq(:,l,igcm_dust_number)= pq(:,l,igcm_dust_number)
-     .            + more_dust(:,l,2)
-      ENDDO
-      ENDIF !! IF (localstorm)
-      tauref(:)=0.
-      ENDIF !! IF (firstcall)
-#endif
-
-      IF (freedust) THEN
-          tauscaling(:) = 1.
-
-      ELSE
-c       Temporary scaling factor
-        taudusttmp(1:ngrid)=0.
-        DO iaer=1,naerdust
-          DO l=1,nlayer
-            DO ig=1,ngrid
-c             Scaling factor
-              taudusttmp(ig) = taudusttmp(ig) + 
-     &                         aerosol(ig,l,iaerdust(iaer))
-            ENDDO
-          ENDDO
-        ENDDO
-
-c       Saved scaling factor
-        DO ig=1,ngrid
-            tauscaling(ig) = tauref(ig) * 
-     &                       pplev(ig,1) / odpref / taudusttmp(ig)
-        ENDDO
-
-      ENDIF ! IF (freedust)
-
-c     Opacity computation
-      DO iaer=1,naerdust
-        DO l=1,nlayer
-          DO ig=1,ngrid
-            aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
-     &                aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
-          ENDDO
-        ENDDO
-      ENDDO
-
-      DO l=1,nlayer
-        DO ig=1,ngrid
-          dsodust(ig,l) = max(1E-20,dsodust(ig,l)* tauscaling(ig))
-        ENDDO
-      ENDDO
-
-      IF (freedust) THEN
-        ! tauref has been initialized to 0 before.
-        DO iaer=1,naerdust
-          DO l=1,nlayer
-            DO ig=1,ngrid
-#ifdef DUSTSTORM
-      !! recalculate opacity because storm perturbation has been added
-      IF (firstcall) THEN
-              aerosol(ig,l,iaer) =
-     &          (  0.75 * QREFvis3d(ig,l,iaer) /
-     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
-     &          pq(ig,l,igcm_dust_mass) *
-     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
-      ENDIF
-#endif
-              tauref(ig) = tauref(ig) + 
-     &                    aerosol(ig,l,iaerdust(iaer))
-            ENDDO
-          ENDDO
-        ENDDO
-        tauref(:) = tauref(:) * odpref / pplev(:,1)
-      ENDIF
-      
-c output for debug
-c        IF (ngrid.NE.1) THEN
-c             CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
-c     &      '#',2,taudusttmp)
-c             CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
-c     &      '#',2,tauscaling)
-c        ELSE
-c             CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
-c     &      '#',0,taudusttmp)
-c             CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
-c     &      '#',0,tauscaling)
-c        ENDIF
-c -----------------------------------------------------------------
-c Column integrated visible optical depth in each point
-c -----------------------------------------------------------------
-      DO iaer=1,naerkind
-        do l=1,nlayer
-           do ig=1,ngrid
-             tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
-           end do
-        end do
-      ENDDO
-
-#ifdef DUSTSTORM
-      IF (firstcall) THEN
-        firstcall=.false.
-      ENDIF
-#endif
-
-c -----------------------------------------------------------------
-c Density scaled opacity and column opacity output
-c -----------------------------------------------------------------
-c      dsodust(1:ngrid,1:nlayer) = 0.
-c      DO iaer=1,naerdust
-c        DO l=1,nlayer
-c          DO ig=1,ngrid
-c            dsodust(ig,l) = dsodust(ig,l) +
-c     &                      aerosol(ig,l,iaerdust(iaer)) * g /
-c     &                      (pplev(ig,l) - pplev(ig,l+1))
-c          ENDDO
-c        ENDDO
-c        IF (ngrid.NE.1) THEN
-c          write(txt2,'(i1.1)') iaer
-c          call WRITEDIAGFI(ngrid,'taudust'//txt2,
-c     &                    'Dust col opacity',
-c     &                    ' ',2,tau(1,iaerdust(iaer)))
-c          IF (callstats) THEN
-c            CALL wstats(ngrid,'taudust'//txt2,
-c     &                 'Dust col opacity',
-c     &                 ' ',2,tau(1,iaerdust(iaer)))
-c          ENDIF
-c        ENDIF
-c      ENDDO
-
-c      IF (ngrid.NE.1) THEN
-c       CALL WRITEDIAGFI(ngrid,'dsodust','tau*g/dp',
-c    &                    'm2.kg-1',3,dsodust)
-c        IF (callstats) THEN
-c          CALL wstats(ngrid,'dsodust',
-c     &               'tau*g/dp',
-c     &               'm2.kg-1',3,dsodust)
-c        ENDIF
-c      ELSE
-c        CALL WRITEDIAGFI(ngrid,"dsodust","dsodust","m2.kg-1",1,
-c     &                   dsodust)
-c      ENDIF ! of IF (ngrid.NE.1)
-c -----------------------------------------------------------------
-      return
-      end 
Index: trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F	(revision 1711)
+++ trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F	(revision 1711)
@@ -0,0 +1,771 @@
+      MODULE aeropacity_mod
+
+      IMPLICIT NONE
+
+      CONTAINS
+
+      SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
+     &    pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
+     &    nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d,
+     &    clearsky,totcloudfrac)
+                                                   
+! to use  'getin'
+      USE ioipsl_getincom, only: getin
+      use tracer_mod, only: noms, igcm_h2o_ice, igcm_dust_mass,
+     &                      igcm_dust_submicron, rho_dust, rho_ice,
+     &                      nqdust
+      use geometry_mod, only: latitude ! grid point latitudes (rad)
+      use comgeomfi_h, only: sinlat ! sines of grid point latitudes
+#ifdef DUSTSTORM
+      use geometry_mod, only: longitude
+      use tracer_mod, only: r3n_q, ref_r0, igcm_dust_number
+#endif
+      use planete_h
+      USE comcstfi_h
+      use dimradmars_mod, only: naerkind, name_iaer, 
+     &            iaerdust,tauvis,
+     &            iaer_dust_conrath,iaer_dust_doubleq,
+     &            iaer_dust_submicron,iaer_h2o_ice
+       IMPLICIT NONE
+c=======================================================================
+c   subject:
+c   --------
+c   Computing aerosol optical depth in each gridbox.
+c
+c   author: F.Forget 
+c   ------
+c   update F. Montmessin (water ice scheme) 
+c      and S. Lebonnois (12/06/2003) compatibility dust/ice/chemistry
+c   update J.-B. Madeleine 2008-2009:
+c       - added 3D scattering by aerosols;
+c       - dustopacity transferred from physiq.F to callradite.F,
+c           and renamed into aeropacity.F;
+c   update E. Millour, march 2012:
+c         - reference pressure is now set to 610Pa (not 700Pa)
+c   
+c   input:
+c   ----- 
+c   ngrid             Number of gridpoint of horizontal grid
+c   nlayer            Number of layer
+c   nq                Number of tracer
+c   zday                  Date (time since Ls=0, in martian days)
+c   ls                Solar longitude (Ls) , radian
+c   pplay,pplev       pressure (Pa) in the middle and boundary of each layer
+c   pq                Dust mixing ratio (used if tracer =T and active=T).
+c   reffrad(ngrid,nlayer,naerkind)  Aerosol effective radius
+c   QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
+c   QREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths;
+c   omegaREFvis3d(ngrid,nlayer,naerkind) \ 3d single scat. albedo
+c   omegaREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths;
+c
+c   output:
+c   -------
+c   tauref       Prescribed mean column optical depth at 610 Pa 
+c   tau          Column total visible dust optical depth at each point
+c   aerosol      aerosol(ig,l,1) is the dust optical
+c                depth in layer l, grid point ig
+
+c
+c=======================================================================
+#include "callkeys.h"
+
+c-----------------------------------------------------------------------
+c
+c    Declarations :
+c    --------------
+c
+c    Input/Output
+c    ------------
+      INTEGER ngrid,nlayer,nq
+
+      REAL ls,zday,expfactor    
+      REAL pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
+      REAL pq(ngrid,nlayer,nq)
+      REAL tauref(ngrid), tau(ngrid,naerkind)
+      REAL aerosol(ngrid,nlayer,naerkind)
+      REAL dsodust(ngrid,nlayer)
+      REAL reffrad(ngrid,nlayer,naerkind)
+      REAL nueffrad(ngrid,nlayer,naerkind)
+      REAL QREFvis3d(ngrid,nlayer,naerkind)
+      REAL QREFir3d(ngrid,nlayer,naerkind)
+      REAL omegaREFvis3d(ngrid,nlayer,naerkind)
+      REAL omegaREFir3d(ngrid,nlayer,naerkind)
+      real,intent(in) :: totcloudfrac(ngrid) ! total cloud fraction
+      logical,intent(in) :: clearsky ! true for part without clouds, 
+                                     ! false for part with clouds (total or sub-grid clouds)
+      real :: CLFtot ! total cloud fraction
+c
+c    Local variables :
+c    -----------------
+      INTEGER l,ig,iq,i,j
+      INTEGER iaer           ! Aerosol index
+      real topdust(ngrid)
+      real zlsconst, zp
+      real taueq,tauS,tauN
+c     Mean Qext(vis)/Qext(ir) profile
+      real msolsir(nlayer,naerkind)
+c     Mean Qext(ir)/Qabs(ir) profile
+      real mqextsqabs(nlayer,naerkind)
+c     Variables used when multiple particle sizes are used
+c       for dust or water ice particles in the radiative transfer
+c       (see callradite.F for more information).
+      REAL taudusttmp(ngrid)! Temporary dust opacity
+                               !   used before scaling
+      REAL tauscaling(ngrid) ! Scaling factor for qdust and Ndust
+      REAL taudustvis(ngrid) ! Dust opacity after scaling
+      REAL taudusttes(ngrid) ! Dust opacity at IR ref. wav. as
+                               !   "seen" by the GCM.
+      REAL taucloudvis(ngrid)! Cloud opacity at visible
+                               !   reference wavelength
+      REAL taucloudtes(ngrid)! Cloud opacity at infrared
+                               !   reference wavelength using
+                               !   Qabs instead of Qext
+                               !   (direct comparison with TES)
+      REAL topdust0(ngrid)
+
+#ifdef DUSTSTORM
+!! Local dust storms
+      logical localstorm        ! =true to create a local dust storm
+      real taulocref,ztoploc,radloc,lonloc,latloc  ! local dust storm parameters
+      real reffstorm, yeah
+      REAL ray(ngrid) ! distance from dust storm center
+      REAL tauuser(ngrid) ! opacity perturbation due to dust storm
+      REAL more_dust(ngrid,nlayer,2) ! Mass mixing ratio perturbation due to the dust storm
+      REAL int_factor(ngrid) ! useful factor to compute mmr perturbation
+      real l_top ! layer of the storm's top
+      REAL zalt(ngrid, nlayer) ! useful factor to compute l_top 
+#endif
+
+c   local saved variables
+c   ---------------------
+
+c     Level under which the dust mixing ratio is held constant
+c       when computing the dust opacity in each layer
+c       (this applies when doubleq and active are true)
+      INTEGER, PARAMETER :: cstdustlevel0 = 7
+      INTEGER, SAVE      :: cstdustlevel
+
+      LOGICAL,SAVE :: firstcall=.true.
+
+! indexes of water ice and dust tracers:
+      INTEGER,SAVE :: i_ice=0  ! water ice
+      real,parameter :: odpref=610. ! DOD reference pressure (Pa)
+      CHARACTER(LEN=20) :: txt ! to temporarly store text
+      CHARACTER(LEN=1) :: txt2 ! to temporarly store text
+! indexes of dust scatterers:
+      INTEGER,SAVE :: naerdust ! number of dust scatterers
+
+      tau(1:ngrid,1:naerkind)=0
+
+! identify tracers
+
+      IF (firstcall) THEN
+        ! identify scatterers that are dust
+        naerdust=0
+        DO iaer=1,naerkind
+          txt=name_iaer(iaer)
+          IF (txt(1:4).eq."dust") THEN
+            naerdust=naerdust+1
+            iaerdust(naerdust)=iaer
+          ENDIF
+        ENDDO
+        ! identify tracers which are dust
+        i=0
+        DO iq=1,nq
+          txt=noms(iq)
+          IF (txt(1:4).eq."dust") THEN
+          i=i+1
+          nqdust(i)=iq
+          ENDIF
+        ENDDO
+
+        IF (water.AND.activice) THEN
+          i_ice=igcm_h2o_ice
+          write(*,*) "aeropacity: i_ice=",i_ice
+        ENDIF
+
+c       typical profile of solsir and (1-w)^(-1):
+        msolsir(1:nlayer,1:naerkind)=0
+        mqextsqabs(1:nlayer,1:naerkind)=0
+        WRITE(*,*) "Typical profiles of Qext(vis)/Qext(IR)"
+        WRITE(*,*) "  and Qext(IR)/Qabs(IR):"
+        DO iaer = 1, naerkind ! Loop on aerosol kind
+          WRITE(*,*) "Aerosol # ",iaer
+          DO l=1,nlayer
+            DO ig=1,ngrid
+              msolsir(l,iaer)=msolsir(l,iaer)+
+     &              QREFvis3d(ig,l,iaer)/
+     &              QREFir3d(ig,l,iaer)
+              mqextsqabs(l,iaer)=mqextsqabs(l,iaer)+
+     &              (1.E0-omegaREFir3d(ig,l,iaer))**(-1)
+            ENDDO
+            msolsir(l,iaer)=msolsir(l,iaer)/REAL(ngrid)
+            mqextsqabs(l,iaer)=mqextsqabs(l,iaer)/REAL(ngrid)
+          ENDDO
+          WRITE(*,*) "solsir: ",msolsir(:,iaer)
+          WRITE(*,*) "Qext/Qabs(IR): ",mqextsqabs(:,iaer)
+        ENDDO
+
+!       load value of tauvis from callphys.def (if given there,
+!       otherwise default value read from starfi.nc file will be used)
+        call getin("tauvis",tauvis)
+
+        IF (freedust) THEN
+          cstdustlevel = 1
+        ELSE
+          cstdustlevel = cstdustlevel0
+        ENDIF
+
+
+#ifndef DUSTSTORM
+        firstcall=.false.
+#endif
+
+      END IF
+
+c     Vertical column optical depth at "odpref" Pa 
+c     ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+      IF(freedust) THEN
+         tauref(:) = 0. ! tauref is computed after, instead of being forced
+
+      ELSE IF(iaervar.eq.1) THEN 
+         do ig=1, ngrid
+          tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
+                                       ! or read in starfi
+        end do
+      ELSE IF (iaervar.eq.2) THEN   ! << "Viking" Scenario>>
+
+        tauref(1) = 0.7+.3*cos(ls+80.*pi/180.) ! like seen by VL1
+        do ig=2,ngrid
+          tauref(ig) = tauref(1)
+        end do
+
+      ELSE IF (iaervar.eq.3) THEN  ! << "MGS" scenario >>
+
+        taueq= 0.2 +(0.5-0.2) *(cos(0.5*(ls-4.363)))**14
+        tauS= 0.1 +(0.5-0.1)  *(cos(0.5*(ls-4.363)))**14
+        tauN = 0.1
+c	   if (peri_day.eq.150) then
+c	     tauS=0.1
+c	     tauN=0.1 +(0.5-0.1)  *(cos(0.5*(ls+pi-4.363)))**14
+c	     taueq= 0.2 +(0.5-0.2) *(cos(0.5*(ls+pi-4.363)))**14
+c           endif
+        do ig=1,ngrid
+          if (latitude(ig).ge.0) then
+          ! Northern hemisphere
+            tauref(ig)= tauN +
+     &      (taueq-tauN)*0.5*(1+tanh((45-latitude(ig)*180./pi)*6/60))
+          else
+          ! Southern hemisphere
+            tauref(ig)= tauS +
+     &      (taueq-tauS)*0.5*(1+tanh((45+latitude(ig)*180./pi)*6/60))
+          endif
+        enddo ! of do ig=1,ngrid
+      ELSE IF (iaervar.eq.5) THEN   ! << Escalier Scenario>>
+c         tauref(1) = 0.2
+c         if ((ls.ge.210.*pi/180.).and.(ls.le.330.*pi/180.))
+c    &                              tauref(1) = 2.5
+        tauref(1) = 2.5
+        if ((ls.ge.30.*pi/180.).and.(ls.le.150.*pi/180.))
+     &                              tauref(1) = .2
+
+        do ig=2,ngrid
+          tauref(ig) = tauref(1)
+        end do
+      ELSE IF ((iaervar.ge.6).and.(iaervar.le.8)) THEN
+      ! clim, cold or warm synthetic scenarios
+        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
+      ELSE IF ((iaervar.ge.24).and.(iaervar.le.32))
+     &     THEN  ! << MY... dust scenarios >>
+        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
+      ELSE IF ((iaervar.eq.4).or.
+     &         ((iaervar.ge.124).and.(iaervar.le.126))) THEN
+       ! "old" TES assimation dust scenario (values at 700Pa in files!)
+        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
+      ELSE
+        stop 'problem with iaervar in aeropacity.F'
+      ENDIF
+
+c -----------------------------------------------------------------
+c Computing the opacity in each layer
+c -----------------------------------------------------------------
+
+      DO iaer = 1, naerkind ! Loop on aerosol kind
+c     --------------------------------------------
+        aerkind: SELECT CASE (name_iaer(iaer))
+c==================================================================
+        CASE("dust_conrath") aerkind      ! Typical dust profile
+c==================================================================
+
+c       Altitude of the top of the dust layer
+c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        zlsconst=SIN(ls-2.76)
+        if (iddist.eq.1) then
+          do ig=1,ngrid
+             topdust(ig)=topdustref         ! constant dust layer top
+          end do
+
+        else if (iddist.eq.2) then          ! "Viking" scenario
+          do ig=1,ngrid
+            ! altitude of the top of the aerosol layer (km) at Ls=2.76rad:
+            ! in the Viking year scenario
+            topdust0(ig)=60. -22.*sinlat(ig)**2
+            topdust(ig)=topdust0(ig)+18.*zlsconst
+          end do
+
+        else if(iddist.eq.3) then         !"MGS" scenario
+          do ig=1,ngrid
+            topdust(ig)=60.+18.*zlsconst
+     &                -(32+18*zlsconst)*sin(latitude(ig))**4
+     &                 - 8*zlsconst*(sin(latitude(ig)))**5
+          end do
+        endif
+
+c       Optical depth in each layer :
+c       ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+        if(iddist.ge.1) then
+
+          expfactor=0.
+          DO l=1,nlayer
+            DO ig=1,ngrid
+c             Typical mixing ratio profile 
+              if(pplay(ig,l).gt.odpref
+     $                        /(988.**(topdust(ig)/70.))) then
+                zp=(odpref/pplay(ig,l))**(70./topdust(ig))
+                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
+              else    
+                expfactor=1.e-3
+              endif
+c             Vertical scaling function
+              aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) *
+     &          expfactor *
+     &          QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer)
+            ENDDO
+          ENDDO
+
+        else if(iddist.eq.0) then   
+c         old dust vertical distribution function (pollack90)
+          DO l=1,nlayer
+             DO ig=1,ngrid
+                zp=odpref/pplay(ig,l)
+                aerosol(ig,l,1)= tauref(ig)/odpref *
+     s           (pplev(ig,l)-pplev(ig,l+1))
+     s           *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 )
+             ENDDO
+          ENDDO
+        end if
+
+c==================================================================
+        CASE("dust_doubleq") aerkind! Two-moment scheme for dust
+c        (transport of mass and number mixing ratio)
+c==================================================================
+
+          DO l=1,nlayer
+            IF (l.LE.cstdustlevel) THEN
+c           Opacity in the first levels is held constant to 
+c             avoid unrealistic values due to constant lifting:
+              DO ig=1,ngrid
+                aerosol(ig,l,iaer) = 
+     &          (  0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
+     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
+     &          pq(ig,cstdustlevel,igcm_dust_mass) *
+     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
+                ! DENSITY SCALED OPACITY IN INFRARED:
+                dsodust(ig,l) =
+     &          (  0.75 * QREFir3d(ig,cstdustlevel,iaer) /
+     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
+     &          pq(ig,cstdustlevel,igcm_dust_mass)
+              ENDDO
+            ELSE
+              DO ig=1,ngrid
+                aerosol(ig,l,iaer) =
+     &          (  0.75 * QREFvis3d(ig,l,iaer) /
+     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
+     &          pq(ig,l,igcm_dust_mass) *
+     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
+                ! DENSITY SCALED OPACITY IN INFRARED:
+                dsodust(ig,l) =
+     &          (  0.75 * QREFir3d(ig,l,iaer) /
+     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
+     &          pq(ig,l,igcm_dust_mass)
+              ENDDO
+            ENDIF
+          ENDDO
+
+c==================================================================
+        CASE("dust_submicron") aerkind   ! Small dust population
+c==================================================================
+
+          DO l=1,nlayer
+            IF (l.LE.cstdustlevel) THEN
+c           Opacity in the first levels is held constant to 
+c             avoid unrealistic values due to constant lifting:
+              DO ig=1,ngrid
+                aerosol(ig,l,iaer) = 
+     &          (  0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
+     &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
+     &          pq(ig,cstdustlevel,igcm_dust_submicron) *
+     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
+              ENDDO
+            ELSE
+              DO ig=1,ngrid
+                aerosol(ig,l,iaer) = 
+     &          (  0.75 * QREFvis3d(ig,l,iaer) /
+     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
+     &          pq(ig,l,igcm_dust_submicron) *
+     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
+              ENDDO
+            ENDIF
+          ENDDO
+
+c==================================================================
+        CASE("h2o_ice") aerkind             ! Water ice crystals
+c==================================================================
+
+c       1. Initialization
+        aerosol(1:ngrid,1:nlayer,iaer) = 0.
+        taucloudvis(1:ngrid) = 0.
+        taucloudtes(1:ngrid) = 0.
+c       2. Opacity calculation
+        ! NO CLOUDS
+        IF (clearsky) THEN
+            aerosol(1:ngrid,1:nlayer,iaer) =1.e-9
+        ! CLOUDSs
+        ELSE ! else (clearsky)
+          DO ig=1, ngrid
+            DO l=1,nlayer
+              aerosol(ig,l,iaer) = max(1E-20,
+     &          (  0.75 * QREFvis3d(ig,l,iaer) /
+     &          ( rho_ice * reffrad(ig,l,iaer) )  ) *
+     &          pq(ig,l,i_ice) *
+     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
+     &                              )
+              taucloudvis(ig) = taucloudvis(ig) + aerosol(ig,l,iaer)
+              taucloudtes(ig) = taucloudtes(ig) + aerosol(ig,l,iaer)*
+     &          QREFir3d(ig,l,iaer) / QREFvis3d(ig,l,iaer) *
+     &          ( 1.E0 - omegaREFir3d(ig,l,iaer) )
+            ENDDO
+          ENDDO
+          ! SUB-GRID SCALE CLOUDS
+          IF (CLFvarying) THEN
+             DO ig=1, ngrid
+                DO l=1,nlayer-1
+                   CLFtot  = max(totcloudfrac(ig),0.01) 
+                   aerosol(ig,l,iaer)=
+     &                    aerosol(ig,l,iaer)/CLFtot
+                   aerosol(ig,l,iaer) = 
+     &                    max(aerosol(ig,l,iaer),1.e-9)
+                ENDDO
+             ENDDO
+!          ELSE ! else (CLFvarying)
+!             DO ig=1, ngrid
+!                DO l=1,nlayer-1 ! to stop the rad tran bug
+!                   CLFtot  = CLFfixval
+!                   aerosol(ig,l,iaer)=
+!     &                    aerosol(ig,l,iaer)/CLFtot
+!                   aerosol(ig,l,iaer) = 
+!     &                    max(aerosol(ig,l,iaer),1.e-9)
+!                ENDDO
+!             ENDDO
+          ENDIF ! end (CLFvarying)             
+        ENDIF ! end (clearsky)
+
+c       3. Outputs -- Now done in physiq.F
+!        IF (ngrid.NE.1) THEN
+!          CALL WRITEDIAGFI(ngrid,'tauVIS','tauext VIS refwvl',
+!     &      ' ',2,taucloudvis)
+!          CALL WRITEDIAGFI(ngrid,'tauTES','tauabs IR refwvl',
+!     &      ' ',2,taucloudtes)
+!          IF (callstats) THEN
+!            CALL wstats(ngrid,'tauVIS','tauext VIS refwvl',
+!     &        ' ',2,taucloudvis)
+!            CALL wstats(ngrid,'tauTES','tauabs IR refwvl',
+!     &        ' ',2,taucloudtes)
+!          ENDIF
+!        ELSE
+!         CALL writeg1d(ngrid,1,taucloudtes,'tautes','NU')
+!        ENDIF
+c==================================================================
+        END SELECT aerkind
+c     -----------------------------------
+      ENDDO ! iaer (loop on aerosol kind)
+
+c -----------------------------------------------------------------
+c Rescaling each layer to reproduce the choosen (or assimilated)
+c   dust extinction opacity at visible reference wavelength, which
+c   is originally scaled to an equivalent odpref Pa pressure surface.
+c -----------------------------------------------------------------
+
+#ifdef DUSTSTORM
+c -----------------------------------------------------------------
+! Calculate reference opacity without perturbation
+c -----------------------------------------------------------------
+      IF (firstcall) THEN
+        DO iaer=1,naerdust
+          DO l=1,nlayer
+            DO ig=1,ngrid
+              tauref(ig) = tauref(ig) +
+     &                    aerosol(ig,l,iaerdust(iaer))
+            ENDDO
+          ENDDO
+        ENDDO
+        tauref(:) = tauref(:) * odpref / pplev(:,1)
+
+c--------------------------------------------------
+c Get parameters of the opacity perturbation
+c--------------------------------------------------
+        iaer=1  ! just change dust
+
+        write(*,*) "Add a local storm ?"
+        localstorm=.true. ! default value
+        call getin("localstorm",localstorm)
+        write(*,*) " localstorm = ",localstorm
+
+        IF (localstorm) THEN
+          WRITE(*,*) "********************"
+          WRITE(*,*) "ADDING A LOCAL STORM"
+          WRITE(*,*) "********************"
+
+          write(*,*) "ref opacity of local dust storm"
+              taulocref = 4.25 ! default value
+              call getin("taulocref",taulocref)
+              write(*,*) " taulocref = ",taulocref
+
+          write(*,*) "target altitude of local storm (km)"
+              ztoploc = 10.0 ! default value
+              call getin("ztoploc",ztoploc)
+              write(*,*) " ztoploc = ",ztoploc
+
+          write(*,*) "radius of dust storm (degree)"
+              radloc = 0.5 ! default value
+              call getin("radloc",radloc)
+              write(*,*) " radloc = ",radloc
+
+          write(*,*) "center longitude of storm (deg)"
+              lonloc = 25.0 ! default value
+              call getin("lonloc",lonloc)
+              write(*,*) " lonloc = ",lonloc
+
+          write(*,*) "center latitude of storm (deg)"
+              latloc = -2.5 ! default value
+              call getin("latloc",latloc)
+              write(*,*) " latloc = ",latloc
+        
+          write(*,*) "reff storm (mic) 0. for background"
+              reffstorm = 0.0 ! default value
+              call getin("reffstorm",reffstorm)
+              write(*,*) " reffstorm = ",reffstorm
+
+!! LOOP: modify opacity
+      DO ig=1,ngrid
+
+      !! distance to the center:
+      ray(ig)=SQRT((latitude(ig)*180./pi-latloc)**2 +
+     &          (longitude(ig)*180./pi -lonloc)**2)
+
+      !! transition factor for storm
+      !! factor is hardcoded -- increase it to steepen
+      yeah = (TANH(2.+(radloc-ray(ig))*10.)+1.)/2.
+
+      !! new opacity field
+      !! -- add an opacity set to taulocref
+      !! -- the additional reference opacity will
+      !!      thus be taulocref*odpref/pplev
+      tauuser(ig)=max( tauref(ig) * pplev(ig,1) /odpref , 
+     &          taulocref * yeah )
+
+      !! compute l_top
+          DO l=1,nlayer
+            zalt(ig,l) = LOG( pplev(ig,1)/pplev(ig,l) )
+     &                      / g / 44.01
+     &                    * 8.31 * 210.
+                IF (     (ztoploc .lt. zalt(ig,l)  )
+     &          .and. (ztoploc .gt. zalt(ig,l-1)) ) l_top=l-1
+          ENDDO
+
+     !! change reffrad if ever needed
+      IF (reffstorm .gt. 0.) THEN
+          DO l=1,nlayer
+             IF (l .lt. l_top+1) THEN
+                reffrad(ig,l,iaer) = max( reffrad(ig,l,iaer), reffstorm
+     &          * 1.e-6 * yeah )
+             ENDIF
+          ENDDO
+      ENDIF
+
+      ENDDO
+!! END LOOP
+
+      !! compute perturbation in each layer (equation 8 in Spiga et al. JGR 2013)
+      DO ig=1,ngrid
+          int_factor(ig)=0.
+          DO l=1,nlayer
+             IF (l .lt. l_top+1) THEN
+                      int_factor(ig) =
+     &                int_factor(ig) +
+     &          (  0.75 * QREFvis3d(ig,l,iaer) /
+     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
+     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
+             ENDIF
+          ENDDO
+          DO l=1, nlayer
+          !! Mass mixing ratio perturbation due to local dust storm in each layer
+          more_dust(ig,l,1)=
+     &                     (tauuser(ig)-(tauref(ig)
+     &                      * pplev(ig,1) /odpref)) /
+     &                      int_factor(ig)
+          more_dust(ig,l,2)=
+     &                     (tauuser(ig)-(tauref(ig) *
+     &                     pplev(ig,1) /odpref))
+     &                      / int_factor(ig) *
+     &                     ((ref_r0/reffrad(ig,l,iaer))**3)
+     &                      * r3n_q 
+          ENDDO
+      ENDDO
+
+      !! quantity of dust for each layer with the addition of the perturbation
+      DO l=1, l_top
+          pq(:,l,igcm_dust_mass)= pq(:,l,igcm_dust_mass)
+     .            + more_dust(:,l,1)
+          pq(:,l,igcm_dust_number)= pq(:,l,igcm_dust_number)
+     .            + more_dust(:,l,2)
+      ENDDO
+      ENDIF !! IF (localstorm)
+      tauref(:)=0.
+      ENDIF !! IF (firstcall)
+#endif
+
+      IF (freedust) THEN
+          tauscaling(:) = 1.
+
+      ELSE
+c       Temporary scaling factor
+        taudusttmp(1:ngrid)=0.
+        DO iaer=1,naerdust
+          DO l=1,nlayer
+            DO ig=1,ngrid
+c             Scaling factor
+              taudusttmp(ig) = taudusttmp(ig) + 
+     &                         aerosol(ig,l,iaerdust(iaer))
+            ENDDO
+          ENDDO
+        ENDDO
+
+c       Saved scaling factor
+        DO ig=1,ngrid
+            tauscaling(ig) = tauref(ig) * 
+     &                       pplev(ig,1) / odpref / taudusttmp(ig)
+        ENDDO
+
+      ENDIF ! IF (freedust)
+
+c     Opacity computation
+      DO iaer=1,naerdust
+        DO l=1,nlayer
+          DO ig=1,ngrid
+            aerosol(ig,l,iaerdust(iaer)) = max(1E-20,
+     &                aerosol(ig,l,iaerdust(iaer))* tauscaling(ig))
+          ENDDO
+        ENDDO
+      ENDDO
+
+      DO l=1,nlayer
+        DO ig=1,ngrid
+          dsodust(ig,l) = max(1E-20,dsodust(ig,l)* tauscaling(ig))
+        ENDDO
+      ENDDO
+
+      IF (freedust) THEN
+        ! tauref has been initialized to 0 before.
+        DO iaer=1,naerdust
+          DO l=1,nlayer
+            DO ig=1,ngrid
+#ifdef DUSTSTORM
+      !! recalculate opacity because storm perturbation has been added
+      IF (firstcall) THEN
+              aerosol(ig,l,iaer) =
+     &          (  0.75 * QREFvis3d(ig,l,iaer) /
+     &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
+     &          pq(ig,l,igcm_dust_mass) *
+     &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
+      ENDIF
+#endif
+              tauref(ig) = tauref(ig) + 
+     &                    aerosol(ig,l,iaerdust(iaer))
+            ENDDO
+          ENDDO
+        ENDDO
+        tauref(:) = tauref(:) * odpref / pplev(:,1)
+      ENDIF
+      
+c output for debug
+c        IF (ngrid.NE.1) THEN
+c             CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
+c     &      '#',2,taudusttmp)
+c             CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
+c     &      '#',2,tauscaling)
+c        ELSE
+c             CALL WRITEDIAGFI(ngrid,'taudusttmp','virtual tau dust',
+c     &      '#',0,taudusttmp)
+c             CALL WRITEDIAGFI(ngrid,'tausca','tauscaling',
+c     &      '#',0,tauscaling)
+c        ENDIF
+c -----------------------------------------------------------------
+c Column integrated visible optical depth in each point
+c -----------------------------------------------------------------
+      DO iaer=1,naerkind
+        do l=1,nlayer
+           do ig=1,ngrid
+             tau(ig,iaer) = tau(ig,iaer) + aerosol(ig,l,iaer)
+           end do
+        end do
+      ENDDO
+
+#ifdef DUSTSTORM
+      IF (firstcall) THEN
+        firstcall=.false.
+      ENDIF
+#endif
+
+c -----------------------------------------------------------------
+c Density scaled opacity and column opacity output
+c -----------------------------------------------------------------
+c      dsodust(1:ngrid,1:nlayer) = 0.
+c      DO iaer=1,naerdust
+c        DO l=1,nlayer
+c          DO ig=1,ngrid
+c            dsodust(ig,l) = dsodust(ig,l) +
+c     &                      aerosol(ig,l,iaerdust(iaer)) * g /
+c     &                      (pplev(ig,l) - pplev(ig,l+1))
+c          ENDDO
+c        ENDDO
+c        IF (ngrid.NE.1) THEN
+c          write(txt2,'(i1.1)') iaer
+c          call WRITEDIAGFI(ngrid,'taudust'//txt2,
+c     &                    'Dust col opacity',
+c     &                    ' ',2,tau(1,iaerdust(iaer)))
+c          IF (callstats) THEN
+c            CALL wstats(ngrid,'taudust'//txt2,
+c     &                 'Dust col opacity',
+c     &                 ' ',2,tau(1,iaerdust(iaer)))
+c          ENDIF
+c        ENDIF
+c      ENDDO
+
+c      IF (ngrid.NE.1) THEN
+c       CALL WRITEDIAGFI(ngrid,'dsodust','tau*g/dp',
+c    &                    'm2.kg-1',3,dsodust)
+c        IF (callstats) THEN
+c          CALL wstats(ngrid,'dsodust',
+c     &               'tau*g/dp',
+c     &               'm2.kg-1',3,dsodust)
+c        ENDIF
+c      ELSE
+c        CALL WRITEDIAGFI(ngrid,"dsodust","dsodust","m2.kg-1",1,
+c     &                   dsodust)
+c      ENDIF ! of IF (ngrid.NE.1)
+c -----------------------------------------------------------------
+
+      END SUBROUTINE aeropacity
+      
+      END MODULE aeropacity_mod
Index: trunk/LMDZ.MARS/libf/phymars/callkeys.h
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 1708)
+++ trunk/LMDZ.MARS/libf/phymars/callkeys.h	(revision 1711)
@@ -14,5 +14,5 @@
      &   ,activice,water,tifeedback,microphys,supersat,caps,photochem   &
      &   ,calltherm,callrichsl,callslope,tituscap,callyamada4,co2clouds,&
-     &   microphysco2
+     &   microphysco2,CLFvarying
      
       COMMON/callkeys_i/iradia,iaervar,iddist,ilwd,ilwb,ilwn,ncouche    &
@@ -20,5 +20,5 @@
      
       COMMON/callkeys_r/topdustref,semi,alphan,euveff,                  &
-     &   tke_heat_flux,dustrefir,fixed_euv_value
+     &   tke_heat_flux,dustrefir,fixed_euv_value,CLFfixval
      
       LOGICAL callrad,calldifv,calladj,callcond,callsoil,               &
@@ -41,4 +41,5 @@
       real euveff
       real tke_heat_flux
+      real CLFfixval
 
       integer iddist
@@ -59,4 +60,5 @@
       logical sedimentation
       logical activice,tifeedback,supersat,caps
+      logical CLFvarying
       logical co2clouds
       logical water
Index: trunk/LMDZ.MARS/libf/phymars/callradite.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callradite.F	(revision 1708)
+++ 	(revision )
@@ -1,548 +1,0 @@
-      SUBROUTINE callradite(icount,ngrid,nlayer,nq,zday,ls,pq,albedo,
-     $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
-     $     dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
-     &     tauref,tau,aerosol,dsodust,tauscaling,taucloudtes,rdust,rice,
-     &     nuice,co2ice)
-
-      use dimradmars_mod, only: ndomainsz, nflev, nsun, nir
-      use dimradmars_mod, only: naerkind, name_iaer,
-     &            iaer_dust_conrath,iaer_dust_doubleq,
-     &            iaer_dust_submicron,iaer_h2o_ice
-      use yomlw_h, only: gcp, nlaylte
-      use comcstfi_h, only: g,cpp
-      use time_phylmdz_mod, only: daysec
-      IMPLICIT NONE
-c=======================================================================
-c   subject:
-c   --------
-c   Subroutine designed to call the main canonic
-c   radiative transfer subroutine "lwmain" et "swmain"
-c   to compute radiative heating and cooling rate and
-c   radiative fluxes to the surface.
-c
-c   These calculations are only valid on the part of the atmosphere
-c   where Local Thermal Equilibrium (NLTE) is verified. In practice
-c   The calculations are only performed for the first "nlaylte"
-c   parameters (nlaylte is calculated by subroutine "nlthermeq"
-c   and stored in module "yomlw_h").
-c
-c   The purpose of this subroutine is to:
-c      1) Make some initial calculation at first call
-c      2) Split the calculation in several sub-grid
-c        ("sub-domain") to save memory and
-c        be able run on a workstation at high resolution
-c        The sub-grid size is defined in dimradmars_mod
-c      3) Compute the 3D scattering parameters depending on the
-c        size distribution of the different tracers (added by JBM)
-c      4) call "lwmain" and "swmain"
-c
-c
-c   authors:   
-c   ------
-c   Francois Forget / Christophe Hourdin / J.-B. Madeleine (2009)
-c
-c
-c   3D scattering scheme user's guide (J.-B. Madeleine)
-c   ---------------------------------
-c
-c   This routine has been modified to take into account 3D, time
-c   dependent scattering properties of the aerosols.
-c---- The look-up tables that contain the scattering parameters
-c   of a given tracer, for different sizes, are read by SUAER.F90.
-c   The names of the corresponding ASCII files have to be set in
-c   this subroutine (file_id variable), and files must be in the
-c   directory specified in datafile.h. Please make sure that the
-c   ASCII files are correctly written, and that the range
-c   of particle sizes is consistent with what you would expect.
-c---- SUAER.F90 is in charge of reading the ASCII files and averaging
-c   the scattering parameters in each GCM channel, using the three last
-c   equations of Forget et al. 1998 (GRL 25, No.7, p.1105-1108).
-c---- These look-up tables, loaded during the firstcall, are then
-c   constantly used by the subroutine "aeroptproperties.F" to compute,
-c   online, the 3D scattering parameters, based on the size distribution
-c   (reffrad and nueffrad) of the different tracers, in each grid box.
-c   These 3D size distributions are loaded by the "updatereffrad.F"
-c   subroutine. A log-normal distribution is then assumed in
-c   "aeroptproperties.F", along with a Gauss-Legendre integration.
-c---- The optical depth at the visible reference wavelength (set in
-c   SUAER.F90, after the file_id variable) is then computed by
-c   the subroutine "aeropacity.F", by using the size and spatial
-c   distribution of the corresponding tracer. This connection has to
-c   be implemented in "aeropacity.F" when adding a new tracer. To do so,
-c   one can use equation 2 of Forget et al. 1998 (Icarus 131, p.302-316).
-c---- The resulting variables "aerosol", "QVISsQREF3d", "omegaVIS3d" and
-c   "gVIS3d" (same in the infrared) are finally used by lwmain.F and 
-c   swmain.F to solve the radiative transfer equation.
-c
-c   changes:
-c   -------
-c
-c   > SRL 7/2000
-c   
-c   This version has been modified to only calculate radiative tendencies
-c   over layers 1..NFLEV (set in dimradmars_mod).  Returns zero for higher
-c   layers, if any.
-c   In other routines, nlayer -> nflev.
-c   Routines affected: lwflux, lwi, lwmain, lwxb, lwxd, lwxn.
-c
-c   > J.-B. Madeleine 10W12
-c   This version uses the variable's splitting, which can be usefull
-c     when performing very high resolution simulation like LES.
-c
-c   ----------
-c   Here, solar band#1 is spectral interval between "long1vis" and "long2vis"
-c   set in dimradmars_mod 
-c   Here, solar band#2 is spectral interval between "long2vis" and "long3vis"
-c   set in dimradmars_mod
-c
-c   input:
-c   ----- 
-c   icount                counter of call to subroutine physic by gcm
-c   ngrid                 number of gridpoint of horizontal grid
-c   nlayer                Number of layer
-c   nq                    Number of tracer
-c   ls                    Solar longitude (Ls) , radian
-c   zday                  Date (time since Ls=0, in martian days)
-c   pq(ngrid,nlayer,nq)   Advected fields
-c
-c   albedo (ngrid,2)      hemispheric surface albedo
-c                         albedo (i,1) : mean albedo for solar band#1 
-c                                        (see below)
-c                         albedo (i,2) : mean albedo for solar band#2
-c                                        (see below)
-c   emis                  Thermal IR surface emissivity (no unit)
-c   mu0(ngrid)           cos of solar zenith angle
-c                           (=1 when sun at zenith)
-c   pplay(ngrid,nlayer)    pressure (Pa) in the middle of each layer
-c   pplev(ngrid,nlayer+1)  pressure (Pa) at boundaries of each layer
-c   pt(ngrid,nlayer)       atmospheric temperature in each layer (K)
-c   tsurf(ngrid)           surface temperature (K)
-c   fract(ngrid)         day fraction of the time interval 
-c                          =1 during the full day ; =0 during the night
-c   declin                 latitude of subsolar point
-c   dist_sol               sun-Mars distance (AU)
-c   igout                  coordinate of analysed point for debugging
-c   reffrad(ngrid,nlayer,naerkind)  Aerosol effective radius
-c   nueffrad(ngrid,nlayer,naerkind) Aerosol effective variance
-
-c
-c  output:
-c  -------
-c dtlw (ngrid,nlayer)       longwave (IR) heating rate (K/s)
-c dtsw(ngrid,nlayer)        shortwave (Solar) heating rate (K/s)
-c fluxsurf_lw(ngrid)        surface downward flux tota LW (thermal IR) (W.m-2)
-c fluxsurf_sw(ngrid,1)      surface downward flux SW for solar band#1 (W.m-2)
-c fluxsurf_sw(ngrid,2)      surface downward flux SW for solar band#2 (W.m-2)
-c
-c fluxtop_lw(ngrid)         outgoing upward flux tota LW (thermal IR) (W.m-2)
-c fluxtop_sw(ngrid,1)       outgoing upward flux SW for solar band#1 (W.m-2)
-c fluxtop_sw(ngrid,2)       outgoing upward flux SW for solar band#2 (W.m-2)
-
-c   tauref       Prescribed mean column optical depth at 610 Pa 
-c   tau          Column total visible dust optical depth at each point
-c   aerosol(ngrid,nlayer,naerkind)    aerosol extinction optical depth
-c                         at reference wavelength "longrefvis" set
-c                         in dimradmars_h , in each layer, for one of
-c                         the "naerkind" kind of aerosol optical
-c                         properties.
-
-c=======================================================================
-c
-c    Declarations :
-c    -------------
-c
-#include "callkeys.h"
-
-c-----------------------------------------------------------------------
-c    Input/Output
-c    ------------
-      INTEGER,INTENT(IN) :: icount        
-      INTEGER,INTENT(IN) :: ngrid,nlayer,nq 
-      INTEGER,INTENT(IN) :: igout
-
-      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
-      REAL,INTENT(IN) :: tauscaling(ngrid) ! Conversion factor for 
-                               ! qdust and Ndust
-      REAL,INTENT(IN) :: albedo(ngrid,2),emis(ngrid)
-      REAL,INTENT(IN) :: ls,zday
-
-      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
-      REAL,INTENT(IN) :: pt(ngrid,nlayer)
-      REAL,INTENT(IN) :: tsurf(ngrid)
-      REAL,INTENT(IN) :: dist_sol,mu0(ngrid),fract(ngrid)
-      REAL,INTENT(OUT) :: dtlw(ngrid,nlayer),dtsw(ngrid,nlayer)
-      REAL,INTENT(OUT) :: fluxsurf_lw(ngrid), fluxtop_lw(ngrid)
-      REAL,INTENT(OUT) :: fluxsurf_sw(ngrid,2), fluxtop_sw(ngrid,2)
-
-      REAL,INTENT(OUT) :: tauref(ngrid), tau(ngrid,naerkind)
-      REAL,INTENT(OUT) :: taucloudtes(ngrid)! Cloud opacity at infrared
-                               !   reference wavelength using
-                               !   Qabs instead of Qext
-                               !   (direct comparison with TES)
-      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind)
-      REAL,INTENT(OUT) :: dsodust(ngrid,nlayer)
-      REAL,INTENT(OUT) :: rdust(ngrid,nlayer)  ! Dust geometric mean radius (m)
-      REAL,INTENT(OUT) :: rice(ngrid,nlayer)   ! Ice geometric mean radius (m)
-      REAL,INTENT(OUT) :: nuice(ngrid,nlayer)  ! Estimated effective variance
-      REAL,INTENT(IN) :: co2ice(ngrid)           ! co2 ice surface layer (kg.m-2)
-
-c
-c    Local variables :
-c    -----------------
-
-      INTEGER j,l,ig,n,ich
-      INTEGER aer_count,iaer
-      INTEGER jd,ig0,nd
-
-      real  cste_mars ! solar constant on Mars (Wm-2)
-      REAL ptlev(ngrid,nlayer+1)
-
-      INTEGER,SAVE :: ndomain
-
-c     Thermal IR net radiative budget (W m-2)
-
-      real znetrad(ndomainsz,nflev)
-
-      real zfluxd_sw(ndomainsz,nflev+1,2)
-      real zfluxu_sw(ndomainsz,nflev+1,2)
-
-      REAL zplev(ndomainsz,nflev+1)
-      REAL zztlev(ndomainsz,nflev+1)
-      REAL zplay(ndomainsz,nflev)
-      REAL zt(ndomainsz,nflev)
-      REAL zaerosol(ndomainsz,nflev,naerkind)
-      REAL zalbedo(ndomainsz,2)
-      REAL zdp(ndomainsz,nflev)
-      REAL zdt0(ndomainsz)
-
-      REAL zzdtlw(ndomainsz,nflev)
-      REAL zzdtsw(ndomainsz,nflev)
-      REAL zzflux(ndomainsz,6)
-      real zrmuz
-
-      REAL :: zQVISsQREF3d(ndomainsz,nflev,nsun,naerkind)
-      REAL :: zomegaVIS3d(ndomainsz,nflev,nsun,naerkind)
-      REAL :: zgVIS3d(ndomainsz,nflev,nsun,naerkind)
-
-      REAL :: zQIRsQREF3d(ndomainsz,nflev,nir,naerkind)
-      REAL :: zomegaIR3d(ndomainsz,nflev,nir,naerkind)
-      REAL :: zgIR3d(ndomainsz,nflev,nir,naerkind)
-
-c     Aerosol size distribution
-      REAL :: reffrad(ngrid,nlayer,naerkind)
-      REAL :: nueffrad(ngrid,nlayer,naerkind)
-c     Aerosol optical properties
-      REAL :: QVISsQREF3d(ngrid,nlayer,nsun,naerkind)
-      REAL :: omegaVIS3d(ngrid,nlayer,nsun,naerkind)
-      REAL :: gVIS3d(ngrid,nlayer,nsun,naerkind)
-
-      REAL :: QIRsQREF3d(ngrid,nlayer,nir,naerkind)
-      REAL :: omegaIR3d(ngrid,nlayer,nir,naerkind)
-      REAL :: gIR3d(ngrid,nlayer,nir,naerkind)
-
-      REAL :: QREFvis3d(ngrid,nlayer,naerkind)
-      REAL :: QREFir3d(ngrid,nlayer,naerkind)
-
-      REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
-      REAL :: omegaREFir3d(ngrid,nlayer,naerkind)
-
-c   local saved variables
-c   ---------------------
-
-      real,save,allocatable :: pview(:)
-      
-      real zco2   ! volume fraction of CO2 in Mars atmosphere
-      DATA zco2/0.95/
-      SAVE zco2
-
-      LOGICAL firstcall
-      DATA firstcall/.true./
-      SAVE firstcall
-
-c----------------------------------------------------------------------
-
-c     Initialisation
-c     --------------
-
-      IF (firstcall) THEN
-        ! compute ndomain and allocate local saved arrays
-        ndomain= (ngrid-1) / ndomainsz + 1
-        allocate(pview(ngrid))
-
-c        Assign a number to the different scatterers
-c        -------------------------------------------
-
-         iaer_dust_conrath=0
-         iaer_dust_doubleq=0
-         iaer_dust_submicron=0
-         iaer_h2o_ice=0
-
-         aer_count=0
-         if (.NOT.active) then
-           do iaer=1,naerkind
-             if (name_iaer(iaer).eq."dust_conrath") then
-               iaer_dust_conrath = iaer
-               aer_count = aer_count + 1
-             endif
-           enddo
-         endif
-         if (doubleq.AND.active) then
-           do iaer=1,naerkind
-             if (name_iaer(iaer).eq."dust_doubleq") then
-               iaer_dust_doubleq = iaer
-               aer_count = aer_count + 1
-             endif
-           enddo
-         endif
-         if (submicron.AND.active) then
-           do iaer=1,naerkind
-             if (name_iaer(iaer).eq."dust_submicron") then
-               iaer_dust_submicron = iaer
-               aer_count = aer_count + 1
-             endif
-           enddo
-         endif
-         if (water.AND.activice) then
-           do iaer=1,naerkind
-             if (name_iaer(iaer).eq."h2o_ice") then
-               iaer_h2o_ice = iaer
-               aer_count = aer_count + 1
-             endif
-           enddo
-         endif
-
-c        Check that we identified all tracers:
-         if (aer_count.ne.naerkind) then
-           write(*,*) "callradite: found only ",aer_count," scatterers"
-           write(*,*) "               expected ",naerkind
-           write(*,*) "please make sure that the number of"
-           write(*,*) "scatterers in scatterers.h, the names"
-           write(*,*) "in callradite.F, and the flags in"
-           write(*,*) "callphys.def are all consistent!"
-           do iaer=1,naerkind
-             write(*,*)'      ',iaer,' ',trim(name_iaer(iaer))
-           enddo
-           stop
-         else
-           write(*,*) "callradite: found all scatterers, namely:"
-           do iaer=1,naerkind
-             write(*,*)'      ',iaer,' ',trim(name_iaer(iaer))
-           enddo
-         endif
-c        -------------------------------------------
-
-         DO ig=1,ngrid
-            pview(ig)=1.66     ! cosecant of viewing angle
-         ENDDO
-         gcp = g/cpp
-
-c        Loading the optical properties in external look-up tables:
-         CALL SUAER
-!         CALL SULW ! this step is now done in ini_yomlw_h
-
-         write(*,*) 'Splitting radiative calculations: ',
-     $              ' ngrid,ndomainsz,ndomain',
-     $                ngrid,ndomainsz,ndomain
-         if (ngrid .EQ. 1) then
-           if (ndomainsz .NE. 1) then
-             print*
-             print*,'ATTENTION !!!'
-             print*,'pour tourner en 1D, '
-             print*,'fixer ndomainsz=1 dans phymars/dimradmars_h'
-             print*
-             call exit(1)
-           endif
-         endif
-         firstcall=.false.
-      END IF
-
-c     Computing aerosol optical properties and opacity
-c     ------------------------------------------------
-
-c     Updating aerosol size distributions:
-      CALL updatereffrad(ngrid,nlayer,
-     &                rdust,rice,nuice,
-     &                reffrad,nueffrad,
-     &                pq,tauscaling,tau,pplay)
-
-c     Computing 3D scattering parameters:
-      CALL aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
-     &                      QVISsQREF3d,omegaVIS3d,gVIS3d,
-     &                      QIRsQREF3d,omegaIR3d,gIR3d,
-     &                      QREFvis3d,QREFir3d,
-     &                      omegaREFvis3d,omegaREFir3d)
-
-c     Computing aerosol optical depth in each layer:
-      CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
-     &     pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
-     &     nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d)
-
-c     Starting loop on sub-domain
-c     ----------------------------
-
-      DO jd=1,ndomain
-        ig0=(jd-1)*ndomainsz
-        if (jd.eq.ndomain) then
-         nd=ngrid-ig0
-        else
-         nd=ndomainsz
-        endif
-
-c       Spliting input variable in sub-domain input variables
-c       ---------------------------------------------------
-
-        do l=1,nlaylte
-         do ig = 1,nd
-           do iaer = 1, naerkind
-             do ich = 1, nsun
-               zQVISsQREF3d(ig,l,ich,iaer) = 
-     &                           QVISsQREF3d(ig0+ig,l,ich,iaer)
-               zomegaVIS3d(ig,l,ich,iaer) = 
-     &                           omegaVIS3d(ig0+ig,l,ich,iaer)
-               zgVIS3d(ig,l,ich,iaer) = 
-     &                           gVIS3d(ig0+ig,l,ich,iaer)
-             enddo
-             do ich = 1, nir
-               zQIRsQREF3d(ig,l,ich,iaer) = 
-     &                           QIRsQREF3d(ig0+ig,l,ich,iaer)
-               zomegaIR3d(ig,l,ich,iaer) = 
-     &                           omegaIR3d(ig0+ig,l,ich,iaer)
-               zgIR3d(ig,l,ich,iaer) = 
-     &                           gIR3d(ig0+ig,l,ich,iaer)
-             enddo
-           enddo
-         enddo
-        enddo
-
-        do l=1,nlaylte+1
-         do ig = 1,nd
-          zplev(ig,l) = pplev(ig0+ig,l)
-         enddo
-        enddo
-
-        do l=1,nlaylte
-         do ig = 1,nd
-          zplay(ig,l) = pplay(ig0+ig,l)
-          zt(ig,l) = pt(ig0+ig,l)
-c         Thickness of each layer (Pa) :
-          zdp(ig,l)= pplev(ig0+ig,l) - pplev(ig0+ig,l+1)
-         enddo
-        enddo
-
-        do n=1,naerkind
-          do l=1,nlaylte
-            do ig=1,nd
-              zaerosol(ig,l,n) = aerosol(ig0+ig,l,n)
-            enddo
-          enddo
-        enddo
-
-        do j=1,2
-          do ig = 1,nd
-           zalbedo(ig,j) = albedo(ig0+ig,j)
-          enddo
-        enddo
-
-c       Intermediate  levels: (computing tlev)
-c       ---------------------------------------
-c       Extrapolation for the air temperature above the surface
-        DO ig=1,nd
-              zztlev(ig,1)=zt(ig,1)+
-     s        (zplev(ig,1)-zplay(ig,1))*
-     s        (zt(ig,1)-zt(ig,2))/(zplay(ig,1)-zplay(ig,2))
-
-              zdt0(ig) = tsurf(ig0+ig) - zztlev(ig,1)
-        ENDDO
-
-        DO l=2,nlaylte
-         DO ig=1, nd
-               zztlev(ig,l)=0.5*(zt(ig,l-1)+zt(ig,l)) 
-         ENDDO
-        ENDDO
-
-        DO ig=1, nd
-           zztlev(ig,nlaylte+1)=zt(ig,nlaylte)
-        ENDDO
-
-
-c       Longwave ("lw") radiative transfer (= thermal infrared)
-c       -------------------------------------------------------
-        call lwmain (ig0,icount,nd,nflev
-     .        ,zdp,zdt0,emis(ig0+1),zplev,zztlev,zt
-     .        ,zaerosol,zzdtlw
-     .        ,fluxsurf_lw(ig0+1),fluxtop_lw(ig0+1)
-     .        ,znetrad
-     &        ,zQIRsQREF3d,zomegaIR3d,zgIR3d
-     &        ,co2ice(ig0+1))
-
-c       Shortwave ("sw") radiative transfer (= solar radiation)
-c       -------------------------------------------------------
-c          Mars solar constant (W m-2)
-c          1370 W.m-2 is the solar constant at 1 AU.
-           cste_mars=1370./(dist_sol*dist_sol)
-
-           call swmain ( nd, nflev,
-     S     cste_mars, zalbedo,
-     S     mu0(ig0+1), zdp, zplev, zaerosol, fract(ig0+1),
-     S     zzdtsw, zfluxd_sw, zfluxu_sw,
-     &     zQVISsQREF3d,zomegaVIS3d,zgVIS3d)
-
-c       ------------------------------------------------------------
-c       Un-spliting output variable from sub-domain input variables
-c       ------------------------------------------------------------
-
-        do l=1,nlaylte
-         do ig = 1,nd
-          dtlw(ig0+ig,l) = zzdtlw(ig,l)
-          dtsw(ig0+ig,l) = zzdtsw(ig,l)
-         enddo
-        enddo
-
-        do l=1,nlaylte+1
-         do ig = 1,nd
-          ptlev(ig0+ig,l) = zztlev(ig,l)
-         enddo
-        enddo
-
-        do ig = 1,nd
-          fluxsurf_sw(ig0+ig,1) = zfluxd_sw(ig,1,1)
-          fluxsurf_sw(ig0+ig,2) = zfluxd_sw(ig,1,2)
-          fluxtop_sw(ig0+ig,1) = zfluxu_sw(ig,nlaylte+1,1)
-          fluxtop_sw(ig0+ig,2) = zfluxu_sw(ig,nlaylte+1,2)
-        enddo
-
-      ENDDO         !   (boucle jd=1, ndomain)
-
-c     Zero tendencies for any remaining layers between nlaylte and nlayer
-      if (nlayer.gt.nlaylte) then
-         do l = nlaylte+1, nlayer
-            do ig = 1, ngrid
-               dtlw(ig, l) = 0.
-               dtsw(ig, l) = 0.
-            enddo
-         enddo
-      endif
-
-c     Output for debugging if lwrite=T
-c     --------------------------------
-c     Write all nlayer layers, even though only nlaylte layers may have
-c     non-zero tendencies.
-
-         IF(lwrite) THEN
-            PRINT*,'Diagnotique for the radiation'
-            PRINT*,'albedo, emissiv, mu0,fract,fluxsurf_lw,fluxsurf_sw'
-            PRINT*,albedo(igout,1),emis(igout),mu0(igout),
-     s           fract(igout), fluxsurf_lw(igout),
-     $     fluxsurf_sw(igout,1)+fluxsurf_sw(igout,2)
-            PRINT*,'Tlay Tlev Play Plev dT/dt SW dT/dt LW (K/s)'
-            PRINT*,'daysec',daysec
-            DO l=1,nlayer
-               PRINT*,pt(igout,l),ptlev(igout,l),
-     s         pplay(igout,l),pplev(igout,l),
-     s         dtsw(igout,l),dtlw(igout,l)
-            ENDDO
-         ENDIF
-
-
-      return
-      end 
Index: trunk/LMDZ.MARS/libf/phymars/callradite_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callradite_mod.F	(revision 1711)
+++ trunk/LMDZ.MARS/libf/phymars/callradite_mod.F	(revision 1711)
@@ -0,0 +1,560 @@
+      MODULE callradite_mod
+
+      IMPLICIT NONE
+
+      CONTAINS
+
+      SUBROUTINE callradite(icount,ngrid,nlayer,nq,zday,ls,pq,albedo,
+     $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
+     $     dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
+     &     tauref,tau,aerosol,dsodust,tauscaling,taucloudtes,rdust,rice,
+     &     nuice,co2ice,clearsky,totcloudfrac)
+
+      use aeropacity_mod
+      use dimradmars_mod, only: ndomainsz, nflev, nsun, nir
+      use dimradmars_mod, only: naerkind, name_iaer,
+     &            iaer_dust_conrath,iaer_dust_doubleq,
+     &            iaer_dust_submicron,iaer_h2o_ice
+      use yomlw_h, only: gcp, nlaylte
+      use comcstfi_h, only: g,cpp
+      use time_phylmdz_mod, only: daysec
+      IMPLICIT NONE
+c=======================================================================
+c   subject:
+c   --------
+c   Subroutine designed to call the main canonic
+c   radiative transfer subroutine "lwmain" et "swmain"
+c   to compute radiative heating and cooling rate and
+c   radiative fluxes to the surface.
+c
+c   These calculations are only valid on the part of the atmosphere
+c   where Local Thermal Equilibrium (NLTE) is verified. In practice
+c   The calculations are only performed for the first "nlaylte"
+c   parameters (nlaylte is calculated by subroutine "nlthermeq"
+c   and stored in module "yomlw_h").
+c
+c   The purpose of this subroutine is to:
+c      1) Make some initial calculation at first call
+c      2) Split the calculation in several sub-grid
+c        ("sub-domain") to save memory and
+c        be able run on a workstation at high resolution
+c        The sub-grid size is defined in dimradmars_mod
+c      3) Compute the 3D scattering parameters depending on the
+c        size distribution of the different tracers (added by JBM)
+c      4) call "lwmain" and "swmain"
+c
+c
+c   authors:   
+c   ------
+c   Francois Forget / Christophe Hourdin / J.-B. Madeleine (2009)
+c
+c
+c   3D scattering scheme user's guide (J.-B. Madeleine)
+c   ---------------------------------
+c
+c   This routine has been modified to take into account 3D, time
+c   dependent scattering properties of the aerosols.
+c---- The look-up tables that contain the scattering parameters
+c   of a given tracer, for different sizes, are read by SUAER.F90.
+c   The names of the corresponding ASCII files have to be set in
+c   this subroutine (file_id variable), and files must be in the
+c   directory specified in datafile.h. Please make sure that the
+c   ASCII files are correctly written, and that the range
+c   of particle sizes is consistent with what you would expect.
+c---- SUAER.F90 is in charge of reading the ASCII files and averaging
+c   the scattering parameters in each GCM channel, using the three last
+c   equations of Forget et al. 1998 (GRL 25, No.7, p.1105-1108).
+c---- These look-up tables, loaded during the firstcall, are then
+c   constantly used by the subroutine "aeroptproperties.F" to compute,
+c   online, the 3D scattering parameters, based on the size distribution
+c   (reffrad and nueffrad) of the different tracers, in each grid box.
+c   These 3D size distributions are loaded by the "updatereffrad.F"
+c   subroutine. A log-normal distribution is then assumed in
+c   "aeroptproperties.F", along with a Gauss-Legendre integration.
+c---- The optical depth at the visible reference wavelength (set in
+c   SUAER.F90, after the file_id variable) is then computed by
+c   the subroutine "aeropacity.F", by using the size and spatial
+c   distribution of the corresponding tracer. This connection has to
+c   be implemented in "aeropacity.F" when adding a new tracer. To do so,
+c   one can use equation 2 of Forget et al. 1998 (Icarus 131, p.302-316).
+c---- The resulting variables "aerosol", "QVISsQREF3d", "omegaVIS3d" and
+c   "gVIS3d" (same in the infrared) are finally used by lwmain.F and 
+c   swmain.F to solve the radiative transfer equation.
+c
+c   changes:
+c   -------
+c
+c   > SRL 7/2000
+c   
+c   This version has been modified to only calculate radiative tendencies
+c   over layers 1..NFLEV (set in dimradmars_mod).  Returns zero for higher
+c   layers, if any.
+c   In other routines, nlayer -> nflev.
+c   Routines affected: lwflux, lwi, lwmain, lwxb, lwxd, lwxn.
+c
+c   > J.-B. Madeleine 10W12
+c   This version uses the variable's splitting, which can be usefull
+c     when performing very high resolution simulation like LES.
+c
+c   ----------
+c   Here, solar band#1 is spectral interval between "long1vis" and "long2vis"
+c   set in dimradmars_mod 
+c   Here, solar band#2 is spectral interval between "long2vis" and "long3vis"
+c   set in dimradmars_mod
+c
+c   input:
+c   ----- 
+c   icount                counter of call to subroutine physic by gcm
+c   ngrid                 number of gridpoint of horizontal grid
+c   nlayer                Number of layer
+c   nq                    Number of tracer
+c   ls                    Solar longitude (Ls) , radian
+c   zday                  Date (time since Ls=0, in martian days)
+c   pq(ngrid,nlayer,nq)   Advected fields
+c
+c   albedo (ngrid,2)      hemispheric surface albedo
+c                         albedo (i,1) : mean albedo for solar band#1 
+c                                        (see below)
+c                         albedo (i,2) : mean albedo for solar band#2
+c                                        (see below)
+c   emis                  Thermal IR surface emissivity (no unit)
+c   mu0(ngrid)           cos of solar zenith angle
+c                           (=1 when sun at zenith)
+c   pplay(ngrid,nlayer)    pressure (Pa) in the middle of each layer
+c   pplev(ngrid,nlayer+1)  pressure (Pa) at boundaries of each layer
+c   pt(ngrid,nlayer)       atmospheric temperature in each layer (K)
+c   tsurf(ngrid)           surface temperature (K)
+c   fract(ngrid)         day fraction of the time interval 
+c                          =1 during the full day ; =0 during the night
+c   declin                 latitude of subsolar point
+c   dist_sol               sun-Mars distance (AU)
+c   igout                  coordinate of analysed point for debugging
+c   reffrad(ngrid,nlayer,naerkind)  Aerosol effective radius
+c   nueffrad(ngrid,nlayer,naerkind) Aerosol effective variance
+
+c
+c  output:
+c  -------
+c dtlw (ngrid,nlayer)       longwave (IR) heating rate (K/s)
+c dtsw(ngrid,nlayer)        shortwave (Solar) heating rate (K/s)
+c fluxsurf_lw(ngrid)        surface downward flux tota LW (thermal IR) (W.m-2)
+c fluxsurf_sw(ngrid,1)      surface downward flux SW for solar band#1 (W.m-2)
+c fluxsurf_sw(ngrid,2)      surface downward flux SW for solar band#2 (W.m-2)
+c
+c fluxtop_lw(ngrid)         outgoing upward flux tota LW (thermal IR) (W.m-2)
+c fluxtop_sw(ngrid,1)       outgoing upward flux SW for solar band#1 (W.m-2)
+c fluxtop_sw(ngrid,2)       outgoing upward flux SW for solar band#2 (W.m-2)
+
+c   tauref       Prescribed mean column optical depth at 610 Pa 
+c   tau          Column total visible dust optical depth at each point
+c   aerosol(ngrid,nlayer,naerkind)    aerosol extinction optical depth
+c                         at reference wavelength "longrefvis" set
+c                         in dimradmars_h , in each layer, for one of
+c                         the "naerkind" kind of aerosol optical
+c                         properties.
+
+c=======================================================================
+c
+c    Declarations :
+c    -------------
+c
+#include "callkeys.h"
+
+c-----------------------------------------------------------------------
+c    Input/Output
+c    ------------
+      INTEGER,INTENT(IN) :: icount        
+      INTEGER,INTENT(IN) :: ngrid,nlayer,nq 
+      INTEGER,INTENT(IN) :: igout
+
+      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)
+      REAL,INTENT(IN) :: tauscaling(ngrid) ! Conversion factor for 
+                               ! qdust and Ndust
+      REAL,INTENT(IN) :: albedo(ngrid,2),emis(ngrid)
+      REAL,INTENT(IN) :: ls,zday
+
+      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
+      REAL,INTENT(IN) :: pt(ngrid,nlayer)
+      REAL,INTENT(IN) :: tsurf(ngrid)
+      REAL,INTENT(IN) :: dist_sol,mu0(ngrid),fract(ngrid)
+      REAL,INTENT(OUT) :: dtlw(ngrid,nlayer),dtsw(ngrid,nlayer)
+      REAL,INTENT(OUT) :: fluxsurf_lw(ngrid), fluxtop_lw(ngrid)
+      REAL,INTENT(OUT) :: fluxsurf_sw(ngrid,2), fluxtop_sw(ngrid,2)
+
+      REAL,INTENT(OUT) :: tauref(ngrid), tau(ngrid,naerkind)
+      REAL,INTENT(OUT) :: taucloudtes(ngrid)! Cloud opacity at infrared
+                               !   reference wavelength using
+                               !   Qabs instead of Qext
+                               !   (direct comparison with TES)
+      REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind)
+      REAL,INTENT(OUT) :: dsodust(ngrid,nlayer)
+      REAL,INTENT(OUT) :: rdust(ngrid,nlayer)  ! Dust geometric mean radius (m)
+      REAL,INTENT(OUT) :: rice(ngrid,nlayer)   ! Ice geometric mean radius (m)
+      REAL,INTENT(OUT) :: nuice(ngrid,nlayer)  ! Estimated effective variance
+      REAL,INTENT(IN) :: co2ice(ngrid)           ! co2 ice surface layer (kg.m-2)
+c     sub-grid scale water ice clouds
+      real,intent(out) :: totcloudfrac(ngrid)
+      logical,intent(in) :: clearsky
+
+c
+c    Local variables :
+c    -----------------
+
+      INTEGER j,l,ig,n,ich
+      INTEGER aer_count,iaer
+      INTEGER jd,ig0,nd
+
+      real  cste_mars ! solar constant on Mars (Wm-2)
+      REAL ptlev(ngrid,nlayer+1)
+
+      INTEGER,SAVE :: ndomain
+
+c     Thermal IR net radiative budget (W m-2)
+
+      real znetrad(ndomainsz,nflev)
+
+      real zfluxd_sw(ndomainsz,nflev+1,2)
+      real zfluxu_sw(ndomainsz,nflev+1,2)
+
+      REAL zplev(ndomainsz,nflev+1)
+      REAL zztlev(ndomainsz,nflev+1)
+      REAL zplay(ndomainsz,nflev)
+      REAL zt(ndomainsz,nflev)
+      REAL zaerosol(ndomainsz,nflev,naerkind)
+      REAL zalbedo(ndomainsz,2)
+      REAL zdp(ndomainsz,nflev)
+      REAL zdt0(ndomainsz)
+
+      REAL zzdtlw(ndomainsz,nflev)
+      REAL zzdtsw(ndomainsz,nflev)
+      REAL zzflux(ndomainsz,6)
+      real zrmuz
+
+      REAL :: zQVISsQREF3d(ndomainsz,nflev,nsun,naerkind)
+      REAL :: zomegaVIS3d(ndomainsz,nflev,nsun,naerkind)
+      REAL :: zgVIS3d(ndomainsz,nflev,nsun,naerkind)
+
+      REAL :: zQIRsQREF3d(ndomainsz,nflev,nir,naerkind)
+      REAL :: zomegaIR3d(ndomainsz,nflev,nir,naerkind)
+      REAL :: zgIR3d(ndomainsz,nflev,nir,naerkind)
+
+c     Aerosol size distribution
+      REAL :: reffrad(ngrid,nlayer,naerkind)
+      REAL :: nueffrad(ngrid,nlayer,naerkind)
+c     Aerosol optical properties
+      REAL :: QVISsQREF3d(ngrid,nlayer,nsun,naerkind)
+      REAL :: omegaVIS3d(ngrid,nlayer,nsun,naerkind)
+      REAL :: gVIS3d(ngrid,nlayer,nsun,naerkind)
+
+      REAL :: QIRsQREF3d(ngrid,nlayer,nir,naerkind)
+      REAL :: omegaIR3d(ngrid,nlayer,nir,naerkind)
+      REAL :: gIR3d(ngrid,nlayer,nir,naerkind)
+
+      REAL :: QREFvis3d(ngrid,nlayer,naerkind)
+      REAL :: QREFir3d(ngrid,nlayer,naerkind)
+
+      REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
+      REAL :: omegaREFir3d(ngrid,nlayer,naerkind)
+
+c   local saved variables
+c   ---------------------
+
+      real,save,allocatable :: pview(:)
+      
+      real zco2   ! volume fraction of CO2 in Mars atmosphere
+      DATA zco2/0.95/
+      SAVE zco2
+
+      LOGICAL firstcall
+      DATA firstcall/.true./
+      SAVE firstcall
+
+c----------------------------------------------------------------------
+
+c     Initialisation
+c     --------------
+
+      IF (firstcall) THEN
+        ! compute ndomain and allocate local saved arrays
+        ndomain= (ngrid-1) / ndomainsz + 1
+        allocate(pview(ngrid))
+
+c        Assign a number to the different scatterers
+c        -------------------------------------------
+
+         iaer_dust_conrath=0
+         iaer_dust_doubleq=0
+         iaer_dust_submicron=0
+         iaer_h2o_ice=0
+
+         aer_count=0
+         if (.NOT.active) then
+           do iaer=1,naerkind
+             if (name_iaer(iaer).eq."dust_conrath") then
+               iaer_dust_conrath = iaer
+               aer_count = aer_count + 1
+             endif
+           enddo
+         endif
+         if (doubleq.AND.active) then
+           do iaer=1,naerkind
+             if (name_iaer(iaer).eq."dust_doubleq") then
+               iaer_dust_doubleq = iaer
+               aer_count = aer_count + 1
+             endif
+           enddo
+         endif
+         if (submicron.AND.active) then
+           do iaer=1,naerkind
+             if (name_iaer(iaer).eq."dust_submicron") then
+               iaer_dust_submicron = iaer
+               aer_count = aer_count + 1
+             endif
+           enddo
+         endif
+         if (water.AND.activice) then
+           do iaer=1,naerkind
+             if (name_iaer(iaer).eq."h2o_ice") then
+               iaer_h2o_ice = iaer
+               aer_count = aer_count + 1
+             endif
+           enddo
+         endif
+
+c        Check that we identified all tracers:
+         if (aer_count.ne.naerkind) then
+           write(*,*) "callradite: found only ",aer_count," scatterers"
+           write(*,*) "               expected ",naerkind
+           write(*,*) "please make sure that the number of"
+           write(*,*) "scatterers in scatterers.h, the names"
+           write(*,*) "in callradite.F, and the flags in"
+           write(*,*) "callphys.def are all consistent!"
+           do iaer=1,naerkind
+             write(*,*)'      ',iaer,' ',trim(name_iaer(iaer))
+           enddo
+           stop
+         else
+           write(*,*) "callradite: found all scatterers, namely:"
+           do iaer=1,naerkind
+             write(*,*)'      ',iaer,' ',trim(name_iaer(iaer))
+           enddo
+         endif
+c        -------------------------------------------
+
+         DO ig=1,ngrid
+            pview(ig)=1.66     ! cosecant of viewing angle
+         ENDDO
+         gcp = g/cpp
+
+c        Loading the optical properties in external look-up tables:
+         CALL SUAER
+!         CALL SULW ! this step is now done in ini_yomlw_h
+
+         write(*,*) 'Splitting radiative calculations: ',
+     $              ' ngrid,ndomainsz,ndomain',
+     $                ngrid,ndomainsz,ndomain
+         if (ngrid .EQ. 1) then
+           if (ndomainsz .NE. 1) then
+             print*
+             print*,'ATTENTION !!!'
+             print*,'pour tourner en 1D, '
+             print*,'fixer ndomainsz=1 dans phymars/dimradmars_h'
+             print*
+             call exit(1)
+           endif
+         endif
+         firstcall=.false.
+      END IF
+
+c     Computing aerosol optical properties and opacity
+c     ------------------------------------------------
+
+c     Updating aerosol size distributions:
+      CALL updatereffrad(ngrid,nlayer,
+     &                rdust,rice,nuice,
+     &                reffrad,nueffrad,
+     &                pq,tauscaling,tau,pplay)
+
+c     Computing 3D scattering parameters:
+      CALL aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
+     &                      QVISsQREF3d,omegaVIS3d,gVIS3d,
+     &                      QIRsQREF3d,omegaIR3d,gIR3d,
+     &                      QREFvis3d,QREFir3d,
+     &                      omegaREFvis3d,omegaREFir3d)
+
+c     Computing aerosol optical depth in each layer:
+      CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
+     &     pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
+     &     nueffrad,QREFvis3d,QREFir3d,omegaREFvis3d,omegaREFir3d,
+     &     clearsky,totcloudfrac)
+
+c     Starting loop on sub-domain
+c     ----------------------------
+
+      DO jd=1,ndomain
+        ig0=(jd-1)*ndomainsz
+        if (jd.eq.ndomain) then
+         nd=ngrid-ig0
+        else
+         nd=ndomainsz
+        endif
+
+c       Spliting input variable in sub-domain input variables
+c       ---------------------------------------------------
+
+        do l=1,nlaylte
+         do ig = 1,nd
+           do iaer = 1, naerkind
+             do ich = 1, nsun
+               zQVISsQREF3d(ig,l,ich,iaer) = 
+     &                           QVISsQREF3d(ig0+ig,l,ich,iaer)
+               zomegaVIS3d(ig,l,ich,iaer) = 
+     &                           omegaVIS3d(ig0+ig,l,ich,iaer)
+               zgVIS3d(ig,l,ich,iaer) = 
+     &                           gVIS3d(ig0+ig,l,ich,iaer)
+             enddo
+             do ich = 1, nir
+               zQIRsQREF3d(ig,l,ich,iaer) = 
+     &                           QIRsQREF3d(ig0+ig,l,ich,iaer)
+               zomegaIR3d(ig,l,ich,iaer) = 
+     &                           omegaIR3d(ig0+ig,l,ich,iaer)
+               zgIR3d(ig,l,ich,iaer) = 
+     &                           gIR3d(ig0+ig,l,ich,iaer)
+             enddo
+           enddo
+         enddo
+        enddo
+
+        do l=1,nlaylte+1
+         do ig = 1,nd
+          zplev(ig,l) = pplev(ig0+ig,l)
+         enddo
+        enddo
+
+        do l=1,nlaylte
+         do ig = 1,nd
+          zplay(ig,l) = pplay(ig0+ig,l)
+          zt(ig,l) = pt(ig0+ig,l)
+c         Thickness of each layer (Pa) :
+          zdp(ig,l)= pplev(ig0+ig,l) - pplev(ig0+ig,l+1)
+         enddo
+        enddo
+
+        do n=1,naerkind
+          do l=1,nlaylte
+            do ig=1,nd
+              zaerosol(ig,l,n) = aerosol(ig0+ig,l,n)
+            enddo
+          enddo
+        enddo
+
+        do j=1,2
+          do ig = 1,nd
+           zalbedo(ig,j) = albedo(ig0+ig,j)
+          enddo
+        enddo
+
+c       Intermediate  levels: (computing tlev)
+c       ---------------------------------------
+c       Extrapolation for the air temperature above the surface
+        DO ig=1,nd
+              zztlev(ig,1)=zt(ig,1)+
+     s        (zplev(ig,1)-zplay(ig,1))*
+     s        (zt(ig,1)-zt(ig,2))/(zplay(ig,1)-zplay(ig,2))
+
+              zdt0(ig) = tsurf(ig0+ig) - zztlev(ig,1)
+        ENDDO
+
+        DO l=2,nlaylte
+         DO ig=1, nd
+               zztlev(ig,l)=0.5*(zt(ig,l-1)+zt(ig,l)) 
+         ENDDO
+        ENDDO
+
+        DO ig=1, nd
+           zztlev(ig,nlaylte+1)=zt(ig,nlaylte)
+        ENDDO
+
+
+c       Longwave ("lw") radiative transfer (= thermal infrared)
+c       -------------------------------------------------------
+        call lwmain (ig0,icount,nd,nflev
+     .        ,zdp,zdt0,emis(ig0+1),zplev,zztlev,zt
+     .        ,zaerosol,zzdtlw
+     .        ,fluxsurf_lw(ig0+1),fluxtop_lw(ig0+1)
+     .        ,znetrad
+     &        ,zQIRsQREF3d,zomegaIR3d,zgIR3d
+     &        ,co2ice(ig0+1))
+
+c       Shortwave ("sw") radiative transfer (= solar radiation)
+c       -------------------------------------------------------
+c          Mars solar constant (W m-2)
+c          1370 W.m-2 is the solar constant at 1 AU.
+           cste_mars=1370./(dist_sol*dist_sol)
+
+           call swmain ( nd, nflev,
+     S     cste_mars, zalbedo,
+     S     mu0(ig0+1), zdp, zplev, zaerosol, fract(ig0+1),
+     S     zzdtsw, zfluxd_sw, zfluxu_sw,
+     &     zQVISsQREF3d,zomegaVIS3d,zgVIS3d)
+
+c       ------------------------------------------------------------
+c       Un-spliting output variable from sub-domain input variables
+c       ------------------------------------------------------------
+
+        do l=1,nlaylte
+         do ig = 1,nd
+          dtlw(ig0+ig,l) = zzdtlw(ig,l)
+          dtsw(ig0+ig,l) = zzdtsw(ig,l)
+         enddo
+        enddo
+
+        do l=1,nlaylte+1
+         do ig = 1,nd
+          ptlev(ig0+ig,l) = zztlev(ig,l)
+         enddo
+        enddo
+
+        do ig = 1,nd
+          fluxsurf_sw(ig0+ig,1) = zfluxd_sw(ig,1,1)
+          fluxsurf_sw(ig0+ig,2) = zfluxd_sw(ig,1,2)
+          fluxtop_sw(ig0+ig,1) = zfluxu_sw(ig,nlaylte+1,1)
+          fluxtop_sw(ig0+ig,2) = zfluxu_sw(ig,nlaylte+1,2)
+        enddo
+
+      ENDDO         !   (boucle jd=1, ndomain)
+
+c     Zero tendencies for any remaining layers between nlaylte and nlayer
+      if (nlayer.gt.nlaylte) then
+         do l = nlaylte+1, nlayer
+            do ig = 1, ngrid
+               dtlw(ig, l) = 0.
+               dtsw(ig, l) = 0.
+            enddo
+         enddo
+      endif
+
+c     Output for debugging if lwrite=T
+c     --------------------------------
+c     Write all nlayer layers, even though only nlaylte layers may have
+c     non-zero tendencies.
+
+         IF(lwrite) THEN
+            PRINT*,'Diagnotique for the radiation'
+            PRINT*,'albedo, emissiv, mu0,fract,fluxsurf_lw,fluxsurf_sw'
+            PRINT*,albedo(igout,1),emis(igout),mu0(igout),
+     s           fract(igout), fluxsurf_lw(igout),
+     $     fluxsurf_sw(igout,1)+fluxsurf_sw(igout,2)
+            PRINT*,'Tlay Tlev Play Plev dT/dt SW dT/dt LW (K/s)'
+            PRINT*,'daysec',daysec
+            DO l=1,nlayer
+               PRINT*,pt(igout,l),ptlev(igout,l),
+     s         pplay(igout,l),pplev(igout,l),
+     s         dtsw(igout,l),dtlw(igout,l)
+            ENDDO
+         ENDIF
+
+
+      END SUBROUTINE callradite
+
+      END MODULE callradite_mod 
Index: trunk/LMDZ.MARS/libf/phymars/conf_phys.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 1708)
+++ trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 1711)
@@ -438,4 +438,15 @@
          write(*,*) " water = ",water
 
+! sub-grid cloud fraction: fixed
+         write(*,*) "Fixed cloud fraction?"
+         CLFfixval=1.0 ! default value
+         call getin("CLFfixval",CLFfixval)
+         write(*,*) "CLFfixval=",CLFfixval
+! sub-grid cloud fraction: varying
+         write(*,*) "Use partial nebulosity?"
+         CLFvarying=.false. ! default value
+         call getin("CLFvarying",CLFvarying)
+         write(*,*)"CLFvarying=",CLFvarying
+
 !CO2 clouds scheme?
          write(*,*) "Compute CO2 clouds ?"
Index: trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 1708)
+++ trunk/LMDZ.MARS/libf/phymars/dyn1d/testphys1d.F	(revision 1711)
@@ -118,4 +118,8 @@
       character(len=44) :: txt
 
+!MV17 nuages alizee
+!     included by AP14 for cloud fraction computations- see BC generic
+      !real,allocatable,dimension(:,:),save :: cloudfrac
+      real,allocatable,dimension(:),save :: totcloudfrac
 c=======================================================================
 
@@ -225,4 +229,7 @@
         allocate(dqdyn(nlayer,nq))
         allocate(mqtot(nq))
+!MV17 nuages alizee
+        !allocate(cloudfrac(ngrid,nlayermx))!essai allocate AP14
+        allocate(totcloudfrac(ngrid))
         
         ! read tracer names from file traceur.def
@@ -711,5 +718,6 @@
       call physdem1("startfi.nc",nsoilmx,ngrid,llm,nq,
      .              dtphys,time,
-     .              tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling)
+     .              tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling,
+     .              totcloudfrac) !MV17 nuages alizee
 
 c=======================================================================
Index: trunk/LMDZ.MARS/libf/phymars/phyetat0.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyetat0.F90	(revision 1708)
+++ trunk/LMDZ.MARS/libf/phymars/phyetat0.F90	(revision 1711)
@@ -1,5 +1,5 @@
 subroutine phyetat0 (fichnom,tab0,Lmodif,nsoil,ngrid,nlay,nq, &
                      day_ini,time0,tsurf,tsoil,emis,q2,qsurf,co2ice, &
-                     tauscaling)
+                     tauscaling,totcloudfrac)
 !  use netcdf
   use tracer_mod, only: noms ! tracer names
@@ -46,5 +46,5 @@
   real,intent(out) :: co2ice(ngrid) ! co2 ice cover
   real,intent(out) :: tauscaling(ngrid) ! dust conversion factor
-
+  real,intent(out) :: totcloudfrac(ngrid) ! total cloud fraction
 !======================================================================
 !  Local variables:
@@ -241,4 +241,13 @@
 endif
 
+! Sub-grid cloud fraction
+call get_field("totcloudfrac",totcloudfrac,found,indextime)
+if (.not.found) then
+  write(*,*) "phyetat0: <totcloudfrac> not in file WARNING put to 1"
+  totcloudfrac(:) = 1.0 !valeur par defaut (CLFfixval par defaut)
+else
+  write(*,*) "phyetat0: total cloud fraction <totcloudfrac> range:", &
+             minval(totcloudfrac), maxval(totcloudfrac)
+endif
 
 ! Surface temperature :
Index: trunk/LMDZ.MARS/libf/phymars/phyredem.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 1708)
+++ trunk/LMDZ.MARS/libf/phymars/phyredem.F90	(revision 1711)
@@ -145,5 +145,5 @@
 subroutine physdem1(filename,nsoil,ngrid,nlay,nq, &
                     phystep,time,tsurf,tsoil,co2ice,emis,q2,qsurf,&
-                    tauscaling)
+                    tauscaling,totcloudfrac)
   ! write time-dependent variable to restart file
   use iostart, only : open_restartphy, close_restartphy, & 
@@ -165,4 +165,5 @@
   real,intent(in) :: qsurf(ngrid,nq)
   real,intent(in) :: tauscaling(ngrid)
+  real, intent(in) :: totcloudfrac(ngrid)
   
   integer :: iq
@@ -194,5 +195,8 @@
   ! Planetary Boundary Layer
   call put_field("q2","pbl wind variance",q2,time)
-  
+
+  ! Sub-grid cloud fraction
+  call put_field("totcloudfrac","Total cloud fraction",totcloudfrac,time) 
+ 
   ! Dust conversion factor
   ! Only to be read by newstart to convert to actual dust values
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 1708)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 1711)
@@ -14,4 +14,7 @@
      $            ,pdu,pdv,pdt,pdq,pdpsrf)
 
+      use watercloud_mod
+      use aeropacity_mod
+      use callradite_mod 
       use tracer_mod, only: noms, mmol, igcm_co2, igcm_n2, igcm_co2_ice,
      &                      igcm_co, igcm_o, igcm_h2o_vap, igcm_h2o_ice,
@@ -386,4 +389,19 @@
 !      REAL zu2(ngrid)
 
+c sub-grid scale water ice clouds (A. Pottier 2013)
+      logical clearsky
+      ! flux for the part without clouds
+      real zdtsw1(ngrid,nlayer)
+      real zdtlw1(ngrid,nlayer)
+      real fluxsurf_lw1(ngrid)     
+      real fluxsurf_sw1(ngrid,2)  
+      real fluxtop_lw1(ngrid)
+      real fluxtop_sw1(ngrid,2) 
+      REAL taucloudtes1(ngrid)
+      ! tf: fraction of clouds, ntf: fraction without clouds
+      real tf, ntf
+      real,allocatable,dimension(:),save :: totcloudfrac ! total cloud fraction over the column
+      REAL rave2(ngrid), totrave2(ngrid) ! Mean water ice mean radius (m)
+
 c=======================================================================
 
@@ -394,4 +412,9 @@
 c  ---------------------------------------
       IF (firstcall) THEN
+        ALLOCATE(totcloudfrac(ngrid)) !A. Pottier: souci avec ca a priori, tentative
+                                      ! de les allouer dans dimradmars_mod.F90 
+        !ouverture de nebuldata.out et nettoyage du precedent
+         open(144,file="nebuldata.out")
+         close(144,status="delete")
 
 c        variables set to 0
@@ -412,5 +435,6 @@
      &         nsoilmx,ngrid,nlayer,nq,
      &         day_ini,time_phys,
-     &         tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling)
+     &         tsurf,tsoil,emis,q2,qsurf,co2ice,tauscaling,
+     &         totcloudfrac)
 
          if (pday.ne.day_ini) then
@@ -552,4 +576,8 @@
          call solarlong(float(day_ini),zls)
       end if
+      !A. Pottier: checking nebuldata.out
+      open(14,file="nebuldata.out",position="append")
+      write(14,*) "LsEtJour ",zls*180/pi,zday
+      close(14)
 
 c     Initialize pressure levels
@@ -699,11 +727,51 @@
 c          Radiative transfer
 c          ------------------
-
+           ! callradite for the part with clouds
+           clearsky=.false. ! part with clouds for both cases CLFvarying true/false 
            CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo,
      $     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
      $     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw,
      $     tauref,tau,aerosol,dsodust,tauscaling,taucloudtes,rdust,rice,
-     $     nuice,co2ice)
-
+     $     nuice,co2ice,clearsky,totcloudfrac)
+           ! case of sub-grid water ice clouds: callradite for the clear case
+            IF (CLFvarying) THEN
+               ! ---> PROBLEMS WITH ALLOCATED ARRAYS 
+               ! (temporary solution in callcorrk: do not deallocate
+               ! if
+               ! CLFvarying ...) ?? AP ??
+               clearsky=.true. ! 
+               CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,
+     &              albedo,emis,mu0,zplev,zplay,pt,tsurf,fract,
+     &              dist_sol,igout,zdtlw1,zdtsw1,fluxsurf_lw1,
+     &              fluxsurf_sw1,fluxtop_lw1,fluxtop_sw1,tauref,tau,
+     &              aerosol,dsodust,tauscaling,taucloudtes1,rdust,
+     &              rice,nuice,co2ice, clearsky, totcloudfrac)
+               clearsky = .false.  ! just in case.
+               ! Sum the fluxes and heating rates from cloudy/clear
+               ! cases
+               DO ig=1,ngrid
+                  tf=totcloudfrac(ig)
+                  ntf=1.-tf
+                  fluxsurf_lw(ig) = ntf*fluxsurf_lw1(ig) 
+     &                                          + tf*fluxsurf_lw(ig)
+                  fluxsurf_sw(ig,1) = ntf*fluxsurf_sw1(ig,1) 
+     &                                          + tf*fluxsurf_sw(ig,1)
+                  fluxsurf_sw(ig,2) = ntf*fluxsurf_sw1(ig,2) 
+     &                                          + tf*fluxsurf_sw(ig,2)
+                  fluxtop_lw(ig)  = ntf*fluxtop_lw1(ig) 
+     &                                          + tf*fluxtop_lw(ig)
+                  fluxtop_sw(ig,1)  = ntf*fluxtop_sw1(ig,1)  
+     &                                          + tf*fluxtop_sw(ig,1)
+                  fluxtop_sw(ig,2)  = ntf*fluxtop_sw1(ig,2)  
+     &                                          + tf*fluxtop_sw(ig,2)
+                  taucloudtes(ig) = ntf*taucloudtes1(ig) 
+     &                                          + tf*taucloudtes(ig)
+                  zdtlw(ig,1:nlayer) = ntf*zdtlw1(ig,1:nlayer) 
+     &                                         + tf*zdtlw(ig,1:nlayer)
+                  zdtsw(ig,1:nlayer) = ntf*zdtsw1(ig,1:nlayer) 
+     &                                         + tf*zdtsw(ig,1:nlayer)
+               ENDDO
+
+            ENDIF ! (CLFvarying)
 
 #ifdef DUSTSTORM
@@ -1082,5 +1150,5 @@
      &                pq,pdq,zdqcloud,zdtcloud,
      &                nq,tau,tauscaling,rdust,rice,nuice,
-     &                rsedcloud,rhocloud)
+     &                rsedcloud,rhocloud,totcloudfrac)
      
 c Temperature variation due to latent heat release
@@ -1740,5 +1808,6 @@
           call physdem1("restartfi.nc",nsoilmx,ngrid,nlayer,nq,
      .                ptimestep,ztime_fin,
-     .                tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling)
+     .                tsurf,tsoil,co2ice,emis,q2,qsurf,tauscaling,
+     .               totcloudfrac)
           
          ENDIF
@@ -1844,4 +1913,15 @@
                enddo
              endif ! of if (scavenging)
+
+           !Alternative A. Pottier weighting
+           rave2(:) = 0.
+           totrave2(:) = 0.
+           do ig=1,ngrid
+              do l=1,nlayer
+              rave2(ig) =rave2(ig)+ zq(ig,l,igcm_h2o_ice)*rice(ig,l)
+              totrave2(ig) = totrave2(ig) + zq(ig,l,igcm_h2o_ice)
+              end do
+              rave2(ig)=max(rave2(ig)/max(totrave2(ig),1.e-30),1.e-30)
+           end do
 
            endif ! of if (water)
@@ -2388,4 +2468,8 @@
      &                       'Mean reff',
      &                       'm',2,rave)
+!A. Pottier
+             CALL WRITEDIAGFI(ngrid,'rmoym',
+     &                      'alternative reffice',
+     &                      'm',2,rave2)
             call WRITEDIAGFI(ngrid,'saturation',
      &           'h2o vap saturation ratio','dimless',3,satu)
@@ -2414,4 +2498,11 @@
      &                       'J.s-1/2.m-2.K-1',3,inertiesoil)
               endif
+!A. Pottier
+              if (CLFvarying) then !AP14 nebulosity
+                 call WRITEDIAGFI(ngrid,'totcloudfrac',
+     &                'Total cloud fraction',
+     &                       ' ',2,totcloudfrac)
+              endif !clf varying
+
            endif !(water)
 
@@ -2885,4 +2976,16 @@
      &                      'm',0,rave)
 
+          !Alternative A. Pottier weighting
+           rave2 = 0.
+           totrave2 = 0.
+           do l=1,nlayer
+              rave2 =rave2+ zq(1,l,igcm_h2o_ice)*rice(1,l)
+              totrave2 = totrave2 + zq(1,l,igcm_h2o_ice)
+           end do
+           rave2=max(rave2/max(totrave2,1.e-30),1.e-30)
+          CALL WRITEDIAGFI(ngrid,'rmoym',
+     &                      'reffice',
+     &                      'm',0,rave2)
+
            do iq=1,nq
              call WRITEDIAGFI(ngrid,trim(noms(iq))//'_s',
@@ -2898,6 +3001,13 @@
      &                          +zdqcloud(1,:,igcm_h2o_vap))
      
-          call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
+        call WRITEDIAGFI(ngrid,"rice","ice radius","m",1,
      &                    rice)
+              
+              if (CLFvarying) then
+                 call WRITEDIAGFI(ngrid,'totcloudfrac',
+     &                'Total cloud fraction',
+     &                       ' ',0,totcloudfrac)
+              endif !clfvarying
+
          ENDIF ! of IF (water)
          
Index: trunk/LMDZ.MARS/libf/phymars/tcondwater.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/tcondwater.F90	(revision 1711)
+++ trunk/LMDZ.MARS/libf/phymars/tcondwater.F90	(revision 1711)
@@ -0,0 +1,71 @@
+SUBROUTINE tcondwater(nbpts,p,q,tcond)
+IMPLICIT NONE 
+
+!---------------------------------------------------
+! Condensation temperature for water ice; based on 
+! the saturation in watersat.F AP14
+!--------------------------------------------------i
+
+integer, intent(in) :: nbpts
+real, intent(in), dimension(nbpts):: p,q
+real, intent(out), dimension(nbpts):: tcond
+
+real, dimension(nbpts):: res
+real:: x
+integer:: i
+
+do i=1,nbpts
+   !write(*,*) "i en cours", i, "sur nbpts=", nbpts
+   !write(*,*) "q(i)",q(i),"p(i)",p(i)
+   x=alog10(MAX(q(i),1e-16)*p(i)/(100.*0.41)) ! max pour erreur q<=0
+   ! attention change le 0.41 de place le 10 juin 2014 car priorites
+   ! fortran dans watersat.F
+   !write(*,*) "x tcondwater AP14! :) :) :) :P", x
+   !res(i) = 2.52826991e+02+ 2.39287870e+01*x+ 2.27275932e+00*x**2
+   !        + 2.21832905e-01*x**3+ 2.23453930e-02*x**4+2.26075106e-03*x**5
+   !        + 2.12411064e-04*x**6+1.64642075e-05*x**7+9.22615632e-07*x**8
+   !        + 3.18958825e-08*x**9+5.00656720e-10*x**10 degre 10: trop!
+   res(i) = 2.52846556e+02+ 2.39229653e+01*x+ 2.21333897e+00*x**2  &
+                + 1.79977992e-01*x**3+ 1.00068175e-02*x**4+2.55145012e-04*x**5
+   !write(*,*) "rex(x) tcondwater AP14! :) :) :) :P", res(i)
+enddo
+
+tcond=res
+
+return 
+end
+!polynome de degre 5 pas 0.0001
+!polynomial coefs [  2.52846556e+02   2.39229653e+01   2.21333897e+00
+!1.79977992e-01
+!   1.00068175e-02   2.55145012e-04]
+!maximum des abs(difference) 0.0604390646026
+
+! polynome de degre 12
+!polynomial coefs [  2.52826992e+02   2.39287716e+01   2.27274870e+00
+!2.21863471e-01
+ !  2.23765903e-02   2.26605393e-03   2.07841624e-04   1.37374700e-05
+ !  2.45106231e-07  -6.16151111e-08  -6.96026651e-09  -3.22690558e-10
+!  -5.86804217e-12]   maximum des abs(difference) 2.73827428146e-05
+! polynome de degre 6
+!polynomial coefs [  2.52831053e+02   2.39333049e+01   2.26006967e+00
+!2.06350715e-01
+!   1.56882616e-02   7.83034223e-04   1.77637297e-05] 
+!maximum des abs(difference) 0.013723768413
+! polynome de degre 8
+! polynomial coefs [  2.52827042e+02   2.39294477e+01   2.27281607e+00
+!2.20637577e-01
+!   2.12491384e-02   1.82810482e-03   1.19638702e-04   4.91167244e-06
+!   9.12272144e-08]
+!maximum des abs(difference) 0.000725035990172
+!polynome de degre 9
+!polynomial coefs [  2.52826985e+02   2.39289186e+01   2.27286626e+00
+!2.21623591e-01
+!   2.20538443e-02   2.11335665e-03   1.73307328e-04   1.04920165e-05
+!   3.94201385e-07   6.70874574e-09]
+!maximum des abs(difference) 0.000168806463876
+!polynome de degre 10
+!polynomial coefs [  2.52826991e+02   2.39287870e+01   2.27275932e+00
+!2.21832905e-01
+!   2.23453930e-02   2.26075106e-03   2.12411064e-04   1.64642075e-05
+!   9.22615632e-07   3.18958825e-08   5.00656720e-10]
+!maximum des abs(difference) 3.96286844477e-05
Index: trunk/LMDZ.MARS/libf/phymars/watercloud.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/watercloud.F	(revision 1708)
+++ 	(revision )
@@ -1,486 +1,0 @@
-       SUBROUTINE watercloud(ngrid,nlay,ptimestep, 
-     &                pplev,pplay,pdpsrf,pzlay,pt,pdt,
-     &                pq,pdq,pdqcloud,pdtcloud,
-     &                nq,tau,tauscaling,rdust,rice,nuice,
-     &                rsedcloud,rhocloud)
-! to use  'getin'
-      USE ioipsl_getincom
-      USE updaterad
-      USE comcstfi_h
-      use tracer_mod, only: nqmx, igcm_h2o_vap, igcm_h2o_ice,
-     &                      igcm_dust_mass, igcm_dust_number,
-     &                      igcm_ccn_mass, igcm_ccn_number,
-     &                      rho_dust, nuice_sed, nuice_ref
-      use dimradmars_mod, only: naerkind
-      IMPLICIT NONE
-
-
-c=======================================================================
-c  Water-ice cloud formation
-c  
-c  Includes two different schemes:
-c    - A simplified scheme (see simpleclouds.F)
-c    - An improved microphysical scheme (see improvedclouds.F)
-c
-c  There is a time loop specific to cloud formation 
-c  due to timescales smaller than the GCM integration timestep.
-c
-c  Authors: Franck Montmessin, Francois Forget, Ehouarn Millour, 
-c           J.-B. Madeleine, Thomas Navarro
-c
-c  2004 - 2012
-c=======================================================================
-
-c-----------------------------------------------------------------------
-c   declarations:
-c   -------------
-
-#include "callkeys.h"
-
-c   Inputs:
-c   ------
-
-      INTEGER ngrid,nlay
-      INTEGER nq                 ! nombre de traceurs 
-      REAL ptimestep             ! pas de temps physique (s)
-      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
-      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
-      REAL pdpsrf(ngrid)         ! tendence surf pressure
-      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
-      REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
-      REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
-
-      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
-      real pdq(ngrid,nlay,nq)    ! tendence avant condensation  (kg/kg.s-1)
-
-      REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
-      REAL tauscaling(ngrid)   ! Convertion factor for dust amount
-      real rdust(ngrid,nlay)   ! Dust geometric mean radius (m)
-
-c   Outputs:
-c   -------
-
-      real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
-      REAL pdtcloud(ngrid,nlay)    ! tendence temperature due
-                                   ! a la chaleur latente
-
-      REAL rice(ngrid,nlay)    ! Ice mass mean radius (m)
-                               ! (r_c in montmessin_2004)
-      REAL nuice(ngrid,nlay)   ! Estimated effective variance
-                               !   of the size distribution
-      real rsedcloud(ngrid,nlay) ! Cloud sedimentation radius
-      real rhocloud(ngrid,nlay)  ! Cloud density (kg.m-3)
-
-c   local:
-c   ------
-      
-      ! for ice radius computation
-      REAL Mo,No
-      REAl ccntyp
-      
-      ! for time loop
-      INTEGER microstep  ! time subsampling step variable
-      INTEGER imicro     ! time subsampling for coupled water microphysics & sedimentation
-      SAVE imicro
-      REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
-      SAVE microtimestep
-      
-      ! tendency given by clouds (inside the micro loop)
-      REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud
-      REAL subpdtcloud(ngrid,nlay)    ! cf. pdtcloud
-
-      ! global tendency (clouds+physics)
-      REAL subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
-      REAL subpdt(ngrid,nlay)         ! cf. pdtcloud
-
-      ! no supersaturation when option supersat is false
-      REAL zt(ngrid,nlay)       ! local value of temperature
-      REAL zqsat(ngrid,nlay)    ! saturation
-
-      INTEGER iq,ig,l
-      LOGICAL,SAVE :: firstcall=.true.
-
-c    ** un petit test de coherence
-c       --------------------------
-
-      IF (firstcall) THEN
-         
-        if (nq.gt.nqmx) then
-           write(*,*) 'stop in watercloud (nq.gt.nqmx)!'
-           write(*,*) 'nq=',nq,' nqmx=',nqmx
-           stop
-        endif
-         
-        write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
-        write(*,*) "            igcm_h2o_ice=",igcm_h2o_ice
-                
-        write(*,*) "time subsampling for microphysic ?"
-#ifdef MESOSCALE
-        imicro = 2
-#else
-        imicro = 30
-#endif
-        call getin("imicro",imicro)
-        write(*,*)"imicro = ",imicro
-        
-        microtimestep = ptimestep/real(imicro)
-        write(*,*)"Physical timestep is",ptimestep 
-        write(*,*)"Microphysics timestep is",microtimestep 
-
-        firstcall=.false.
-      ENDIF ! of IF (firstcall)
-      
-c-----Initialization
-      subpdq(1:ngrid,1:nlay,1:nq) = 0
-      subpdt(1:ngrid,1:nlay)      = 0
-      
-      ! default value if no ice
-      rhocloud(1:ngrid,1:nlay) = rho_dust
-
-
-
-c------------------------------------------------------------------
-c Time subsampling for microphysics 
-c------------------------------------------------------------------
-      DO microstep=1,imicro 
-      
-c-------------------------------------------------------------------
-c   1.  Tendencies: 
-c------------------
-
-
-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
-          ENDDO
-        ENDDO
-c------ Tracers tendencies subpdq
-c------ At each micro timestep we add pdq in order to have a stepped entry
-        IF (microphys) THEN
-          DO l=1,nlay
-            DO ig=1,ngrid
-              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_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
-        ENDIF
-        DO l=1,nlay
-          DO ig=1,ngrid
-            subpdq(ig,l,igcm_h2o_ice) = 
-     &          subpdq(ig,l,igcm_h2o_ice)
-     &        + pdq(ig,l,igcm_h2o_ice)
-            subpdq(ig,l,igcm_h2o_vap) = 
-     &          subpdq(ig,l,igcm_h2o_vap)
-     &        + pdq(ig,l,igcm_h2o_vap)
-          ENDDO
-        ENDDO
-        
-        
-c-------------------------------------------------------------------
-c   2.  Main call to the different cloud schemes:
-c------------------------------------------------
-        IF (microphys) THEN
-           CALL improvedclouds(ngrid,nlay,microtimestep,
-     &             pplay,pt,subpdt,
-     &             pq,subpdq,subpdqcloud,subpdtcloud,
-     &             nq,tauscaling)
-
-        ELSE
-           CALL simpleclouds(ngrid,nlay,microtimestep,
-     &             pplay,pzlay,pt,subpdt,
-     &             pq,subpdq,subpdqcloud,subpdtcloud,
-     &             nq,tau,rice)
-        ENDIF
-
-
-c-------------------------------------------------------------------
-c   3.  Updating tendencies after cloud scheme:
-c-----------------------------------------------
-
-        IF (microphys) THEN
-          DO l=1,nlay
-            DO ig=1,ngrid
-              subpdq(ig,l,igcm_dust_mass) =
-     &            subpdq(ig,l,igcm_dust_mass)
-     &          + subpdqcloud(ig,l,igcm_dust_mass)
-              subpdq(ig,l,igcm_dust_number) =
-     &            subpdq(ig,l,igcm_dust_number)
-     &          + subpdqcloud(ig,l,igcm_dust_number)
-              subpdq(ig,l,igcm_ccn_mass) =
-     &            subpdq(ig,l,igcm_ccn_mass)
-     &          + subpdqcloud(ig,l,igcm_ccn_mass)
-              subpdq(ig,l,igcm_ccn_number) =
-     &            subpdq(ig,l,igcm_ccn_number)
-     &          + subpdqcloud(ig,l,igcm_ccn_number)
-            ENDDO
-          ENDDO
-        ENDIF
-        DO l=1,nlay
-          DO ig=1,ngrid
-            subpdq(ig,l,igcm_h2o_ice) =
-     &          subpdq(ig,l,igcm_h2o_ice)
-     &        + subpdqcloud(ig,l,igcm_h2o_ice)
-            subpdq(ig,l,igcm_h2o_vap) =
-     &          subpdq(ig,l,igcm_h2o_vap)
-     &        + subpdqcloud(ig,l,igcm_h2o_vap)
-          ENDDO
-        ENDDO
-        
-        IF (activice) THEN
-          DO l=1,nlay
-            DO ig=1,ngrid
-              subpdt(ig,l) =
-     &            subpdt(ig,l) + subpdtcloud(ig,l)
-            ENDDO
-          ENDDO
-        ENDIF
-     
- 
-      ENDDO ! of DO microstep=1,imicro
-      
-c-------------------------------------------------------------------
-c   6.  Compute final tendencies after time loop:
-c------------------------------------------------
-
-c------ Temperature tendency pdtcloud
-       DO l=1,nlay
-         DO ig=1,ngrid
-             pdtcloud(ig,l) =
-     &         subpdt(ig,l)/real(imicro)-pdt(ig,l)
-          ENDDO
-       ENDDO
-       
-c------ Tracers tendencies pdqcloud
-       DO l=1,nlay
-         DO ig=1,ngrid
-            pdqcloud(ig,l,igcm_h2o_ice) = 
-     &        subpdq(ig,l,igcm_h2o_ice)/real(imicro) 
-     &       - pdq(ig,l,igcm_h2o_ice)
-            pdqcloud(ig,l,igcm_h2o_vap) = 
-     &        subpdq(ig,l,igcm_h2o_vap)/real(imicro) 
-     &       - pdq(ig,l,igcm_h2o_vap)
-         ENDDO
-       ENDDO
-       
-       IF(microphys) THEN
-        DO l=1,nlay
-         DO ig=1,ngrid
-            pdqcloud(ig,l,igcm_ccn_mass) = 
-     &        subpdq(ig,l,igcm_ccn_mass)/real(imicro) 
-     &       - pdq(ig,l,igcm_ccn_mass)
-            pdqcloud(ig,l,igcm_ccn_number) = 
-     &        subpdq(ig,l,igcm_ccn_number)/real(imicro) 
-     &       - pdq(ig,l,igcm_ccn_number)
-         ENDDO
-        ENDDO
-       ENDIF
-       
-       IF(scavenging) THEN
-        DO l=1,nlay
-         DO ig=1,ngrid
-            pdqcloud(ig,l,igcm_dust_mass) = 
-     &        subpdq(ig,l,igcm_dust_mass)/real(imicro) 
-     &       - pdq(ig,l,igcm_dust_mass)
-            pdqcloud(ig,l,igcm_dust_number) = 
-     &        subpdq(ig,l,igcm_dust_number)/real(imicro)
-     &       - pdq(ig,l,igcm_dust_number)
-         ENDDO
-        ENDDO
-       ENDIF
-
-c------- Due to stepped entry, other processes tendencies can add up to negative values
-c------- Therefore, enforce positive values and conserve mass
-
-
-       IF(microphys) THEN
-        DO l=1,nlay
-         DO ig=1,ngrid
-          IF ((pq(ig,l,igcm_ccn_number) + 
-     &      ptimestep* (pdq(ig,l,igcm_ccn_number) + 
-     &        pdqcloud(ig,l,igcm_ccn_number)) .le. 1.)
-     &   .or. (pq(ig,l,igcm_ccn_mass) + 
-     &      ptimestep* (pdq(ig,l,igcm_ccn_mass) + 
-     &        pdqcloud(ig,l,igcm_ccn_mass)) .le. 1.e-20)) THEN
-         pdqcloud(ig,l,igcm_ccn_number) =
-     &     - pq(ig,l,igcm_ccn_number)/ptimestep 
-     &     - pdq(ig,l,igcm_ccn_number) + 1.
-         pdqcloud(ig,l,igcm_dust_number) =  
-     &     -pdqcloud(ig,l,igcm_ccn_number)
-         pdqcloud(ig,l,igcm_ccn_mass) =
-     &     - pq(ig,l,igcm_ccn_mass)/ptimestep
-     &     - pdq(ig,l,igcm_ccn_mass) + 1.e-20
-         pdqcloud(ig,l,igcm_dust_mass) = 
-     &     -pdqcloud(ig,l,igcm_ccn_mass)
-          ENDIF
-         ENDDO
-        ENDDO
-       ENDIF
-
-       IF(scavenging) THEN
-        DO l=1,nlay
-         DO ig=1,ngrid
-          IF ((pq(ig,l,igcm_dust_number) + 
-     &      ptimestep* (pdq(ig,l,igcm_dust_number) + 
-     &        pdqcloud(ig,l,igcm_dust_number)) .le. 1.)
-     &   .or. (pq(ig,l,igcm_dust_mass) + 
-     &      ptimestep* (pdq(ig,l,igcm_dust_mass) + 
-     &        pdqcloud(ig,l,igcm_dust_mass)) .le. 1.e-20)) THEN
-         pdqcloud(ig,l,igcm_dust_number) =
-     &     - pq(ig,l,igcm_dust_number)/ptimestep 
-     &     - pdq(ig,l,igcm_dust_number) + 1.
-         pdqcloud(ig,l,igcm_ccn_number) =  
-     &     -pdqcloud(ig,l,igcm_dust_number)
-         pdqcloud(ig,l,igcm_dust_mass) =
-     &     - pq(ig,l,igcm_dust_mass)/ptimestep
-     &     - pdq(ig,l,igcm_dust_mass) + 1.e-20
-         pdqcloud(ig,l,igcm_ccn_mass) = 
-     &     -pdqcloud(ig,l,igcm_dust_mass)
-          ENDIF
-         ENDDO
-        ENDDO
-       ENDIF
-
-        DO l=1,nlay
-         DO ig=1,ngrid
-          IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
-     &       (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice)) 
-     &       .le. 1.e-8) THEN
-           pdqcloud(ig,l,igcm_h2o_ice) = 
-     &     - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
-           pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
-          ENDIF 
-          IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
-     &       (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap)) 
-     &       .le. 1.e-8) THEN
-           pdqcloud(ig,l,igcm_h2o_vap) = 
-     &     - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
-           pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
-          ENDIF
-         ENDDO
-        ENDDO
-
-
-c------Update the ice and dust particle size "rice" for output or photochemistry
-c------Only rsedcloud is used for the water cycle
-
-      IF(scavenging) THEN 
-        DO l=1, nlay
-         DO ig=1,ngrid
-
-        call updaterdust(
-     &    pq(ig,l,igcm_dust_mass) +                   ! dust mass
-     &   (pdq(ig,l,igcm_dust_mass) +                  ! dust mass
-     &    pdqcloud(ig,l,igcm_dust_mass))*ptimestep,   ! dust mass
-     &    pq(ig,l,igcm_dust_number) +                 ! dust number
-     &   (pdq(ig,l,igcm_dust_number) +                ! dust number
-     &    pdqcloud(ig,l,igcm_dust_number))*ptimestep, ! dust number
-     &    rdust(ig,l))
-
-         ENDDO
-        ENDDO
-      ENDIF
-
-      IF(microphys) THEN
-
-       ! In case one does not want to allow supersatured water when using microphysics.
-       ! Not done by default.
-       IF(.not.supersat) THEN     
-        zt  = pt + (pdt+pdtcloud)*ptimestep
-        call watersat(ngrid*nlay,zt,pplay,zqsat)
-        DO l=1, nlay
-         DO ig=1,ngrid
-          IF (pq(ig,l,igcm_h2o_vap)
-     &      + (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
-     &      * ptimestep .ge. zqsat(ig,l)) THEN
-             pdqcloud(ig,l,igcm_h2o_vap) = 
-     &         (zqsat(ig,l) - pq(ig,l,igcm_h2o_vap))/ptimestep
-     &        - pdq(ig,l,igcm_h2o_vap)
-             pdqcloud(ig,l,igcm_h2o_ice) = 
-     &         -pdqcloud(ig,l,igcm_h2o_vap)
-             ! no need to correct ccn_number, updaterad can handle this properly.
-          ENDIF
-         ENDDO
-        ENDDO        
-       ENDIF
-       
-       DO l=1, nlay
-         DO ig=1,ngrid
-
-        call updaterice_micro(
-     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
-     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
-     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
-     &    pq(ig,l,igcm_ccn_mass) +                   ! ccn mass
-     &   (pdq(ig,l,igcm_ccn_mass) +                  ! ccn mass
-     &    pdqcloud(ig,l,igcm_ccn_mass))*ptimestep,   ! ccn mass
-     &    pq(ig,l,igcm_ccn_number) +                 ! ccn number
-     &   (pdq(ig,l,igcm_ccn_number) +                ! ccn number
-     &    pdqcloud(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
-     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
-          
-         ENDDO
-       ENDDO
-        
-      ELSE ! no microphys
-        
-        DO l=1,nlay
-          DO ig=1,ngrid
-          
-        call updaterice_typ(
-     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
-     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
-     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
-     &    tau(ig,1),pzlay(ig,l),rice(ig,l)) 
-
-          ENDDO
-         ENDDO
-       
-       ENDIF ! of IF(microphys)
-      
-      
-      
-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)))
-     &    rice(ig,2)=rice(ig,3) 
-          rice(ig,1)=rice(ig,2)
-        end if
-      end do
-       
-       
-       DO l=1,nlay
-         DO ig=1,ngrid
-           rsedcloud(ig,l)=max(rice(ig,l)*
-     &                 (1.+nuice_sed)*(1.+nuice_sed)*(1.+nuice_sed),
-     &                    rdust(ig,l))
-!          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
-         ENDDO
-       ENDDO
-       
-! used for rad. transfer calculations
-! nuice is constant because a lognormal distribution is prescribed
-      nuice(1:ngrid,1:nlay)=nuice_ref 
-
-
-
-c=======================================================================
-
-      END
- 
Index: trunk/LMDZ.MARS/libf/phymars/watercloud_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/watercloud_mod.F	(revision 1711)
+++ trunk/LMDZ.MARS/libf/phymars/watercloud_mod.F	(revision 1711)
@@ -0,0 +1,629 @@
+       MODULE watercloud_mod
+
+       IMPLICIT NONE
+
+       CONTAINS
+
+       SUBROUTINE watercloud(ngrid,nlay,ptimestep, 
+     &                pplev,pplay,pdpsrf,pzlay,pt,pdt,
+     &                pq,pdq,pdqcloud,pdtcloud,
+     &                nq,tau,tauscaling,rdust,rice,nuice,
+     &                rsedcloud,rhocloud,totcloudfrac)
+! to use  'getin'
+      USE ioipsl_getincom
+      USE updaterad
+      USE comcstfi_h
+      use tracer_mod, only: nqmx, igcm_h2o_vap, igcm_h2o_ice,
+     &                      igcm_dust_mass, igcm_dust_number,
+     &                      igcm_ccn_mass, igcm_ccn_number,
+     &                      rho_dust, nuice_sed, nuice_ref
+      use dimradmars_mod, only: naerkind
+      IMPLICIT NONE
+
+
+c=======================================================================
+c  Water-ice cloud formation
+c  
+c  Includes two different schemes:
+c    - A simplified scheme (see simpleclouds.F)
+c    - An improved microphysical scheme (see improvedclouds.F)
+c
+c  There is a time loop specific to cloud formation 
+c  due to timescales smaller than the GCM integration timestep.
+c
+c  Authors: Franck Montmessin, Francois Forget, Ehouarn Millour, 
+c           J.-B. Madeleine, Thomas Navarro
+c
+c  2004 - 2012
+c=======================================================================
+
+c-----------------------------------------------------------------------
+c   declarations:
+c   -------------
+
+#include "callkeys.h"
+
+c   Inputs:
+c   ------
+
+      INTEGER ngrid,nlay
+      INTEGER nq                 ! nombre de traceurs 
+      REAL ptimestep             ! pas de temps physique (s)
+      REAL pplev(ngrid,nlay+1)   ! pression aux inter-couches (Pa)
+      REAL pplay(ngrid,nlay)     ! pression au milieu des couches (Pa)
+      REAL pdpsrf(ngrid)         ! tendence surf pressure
+      REAL pzlay(ngrid,nlay)     ! altitude at the middle of the layers
+      REAL pt(ngrid,nlay)        ! temperature at the middle of the layers (K)
+      REAL pdt(ngrid,nlay)       ! tendence temperature des autres param.
+
+      real pq(ngrid,nlay,nq)     ! traceur (kg/kg)
+      real pdq(ngrid,nlay,nq)    ! tendence avant condensation  (kg/kg.s-1)
+
+      REAL tau(ngrid,naerkind) ! Column dust optical depth at each point
+      REAL tauscaling(ngrid)   ! Convertion factor for dust amount
+      real rdust(ngrid,nlay)   ! Dust geometric mean radius (m)
+
+c   Outputs:
+c   -------
+
+      real pdqcloud(ngrid,nlay,nq) ! tendence de la condensation H2O(kg/kg.s-1)
+      REAL pdtcloud(ngrid,nlay)    ! tendence temperature due
+                                   ! a la chaleur latente
+
+      REAL rice(ngrid,nlay)    ! Ice mass mean radius (m)
+                               ! (r_c in montmessin_2004)
+      REAL nuice(ngrid,nlay)   ! Estimated effective variance
+                               !   of the size distribution
+      real rsedcloud(ngrid,nlay) ! Cloud sedimentation radius
+      real rhocloud(ngrid,nlay)  ! Cloud density (kg.m-3)
+
+      REAL, INTENT(INOUT):: totcloudfrac(ngrid) ! Cloud fraction (A. Pottier 2013)
+c   local:
+c   ------
+      
+      ! for ice radius computation
+      REAL Mo,No
+      REAl ccntyp
+      
+      ! for time loop
+      INTEGER microstep  ! time subsampling step variable
+      INTEGER imicro     ! time subsampling for coupled water microphysics & sedimentation
+      SAVE imicro
+      REAL microtimestep ! integration timestep for coupled water microphysics & sedimentation
+      SAVE microtimestep
+      
+      ! tendency given by clouds (inside the micro loop)
+      REAL subpdqcloud(ngrid,nlay,nq) ! cf. pdqcloud
+      REAL subpdtcloud(ngrid,nlay)    ! cf. pdtcloud
+
+      ! global tendency (clouds+physics)
+      REAL subpdq(ngrid,nlay,nq)      ! cf. pdqcloud
+      REAL subpdt(ngrid,nlay)         ! cf. pdtcloud
+
+      ! no supersaturation when option supersat is false
+      REAL zt(ngrid,nlay)       ! local value of temperature
+      REAL zqsat(ngrid,nlay)    ! saturation
+
+      INTEGER iq,ig,l
+      LOGICAL,SAVE :: firstcall=.true.
+
+! Representation of sub-grid water ice clouds A. Pottier 2013
+!      REAL :: zt(ngrid, nlay)
+      REAL :: zq(ngrid, nlay,nq)
+      REAL :: zdelt  
+      REAL :: norm
+      REAL :: ponder
+      REAL :: tcond(ngrid,nlay)
+      REAL ::  zqvap(ngrid,nlay)
+      REAL :: zqice(ngrid,nlay)
+      REAL ::  spant ! delta T for the temperature distribution
+!      REAL :: zqsat(ngrid,nlay) ! saturation
+      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
+      LOGICAL, save :: flagcloud=.true. 
+c    ** un petit test de coherence
+c       --------------------------
+
+      IF (firstcall) THEN
+         
+        if (nq.gt.nqmx) then
+           write(*,*) 'stop in watercloud (nq.gt.nqmx)!'
+           write(*,*) 'nq=',nq,' nqmx=',nqmx
+           stop
+        endif
+         
+        write(*,*) "watercloud: igcm_h2o_vap=",igcm_h2o_vap
+        write(*,*) "            igcm_h2o_ice=",igcm_h2o_ice
+                
+        write(*,*) "time subsampling for microphysic ?"
+#ifdef MESOSCALE
+        imicro = 2
+#else
+        imicro = 30
+#endif
+        call getin("imicro",imicro)
+        write(*,*)"imicro = ",imicro
+        
+        microtimestep = ptimestep/real(imicro)
+        write(*,*)"Physical timestep is",ptimestep 
+        write(*,*)"Microphysics timestep is",microtimestep 
+
+        firstcall=.false.
+      ENDIF ! of IF (firstcall)
+      
+c-----Initialization
+      subpdq(1:ngrid,1:nlay,1:nq) = 0
+      subpdt(1:ngrid,1:nlay)      = 0
+      
+      ! default value if no ice
+      rhocloud(1:ngrid,1:nlay) = rho_dust
+
+c-------------------------------------------------------------------
+c   0.  Representation of sub-grid water ice clouds
+c------------------
+c-----Tendencies
+      DO l=1,nlay
+       DO ig=1,ngrid
+          zt(ig,l)=pt(ig,l)+ pdt(ig,l)*ptimestep
+       ENDDO
+      ENDDO
+      DO l=1,nlay
+        DO ig=1,ngrid
+          DO iq=1,nq
+             zq(ig,l,iq)=pq(ig,l,iq)+pdq(ig,l,iq)*ptimestep 
+          ENDDO
+        ENDDO
+      ENDDO
+c-----Effective temperature calculation
+      IF (CLFvarying) THEN
+         spant=3.0 ! delta T for the temprature distribution
+         mincloud=0.1 ! min cloudfrac when there is ice
+         IF (flagcloud) THEN
+             write(*,*) "Delta T", spant
+             write(*,*) "mincloud", mincloud
+             flagcloud=.false.
+         END IF
+         CALL watersat(ngrid*nlay,zt,pplay,zqsat)
+         zqvap=zq(:,:,igcm_h2o_vap)
+         zqice=zq(:,:,igcm_h2o_ice)
+         CALL tcondwater(ngrid*nlay,pplay,zqvap+zqice,tcond)
+         DO l=1,nlay
+           DO ig=1,ngrid
+              zdelt=spant !MAX(spant*zt(ig,l),1.e-12), now totally in K. Fixed width
+              IF (tcond(ig,l) .ge. (zt(ig,l)+zdelt)) THEN
+                 zteff(ig,l)=zt(ig,l)
+                 cloudfrac(ig,l)=1.
+              ELSE IF (tcond(ig,l) .le. (zt(ig,l)-zdelt)) THEN
+                 zteff(ig,l)=zt(ig,l)-zdelt
+                 cloudfrac(ig,l)=mincloud
+              ELSE
+                 cloudfrac(ig,l)=(tcond(ig,l)-zt(ig,l)+zdelt)/
+     &                           (2.0*zdelt)
+                 zteff(ig,l)=(tcond(ig,l)+zt(ig,l)-zdelt)/2.
+              END IF
+              zteff(ig,l)=zteff(ig,l)-pdt(ig,l)*ptimestep
+              IF (cloudfrac(ig,l).le.0) THEN
+                 cloudfrac(ig,l)=mincloud
+              ELSE IF (cloudfrac(ig,l).gt.1) THEN
+                 cloudfrac(ig,l)=1.
+              END IF
+           ENDDO
+         ENDDO
+c-----Calculation of the total cloud coverage of the column
+         DO ig=1,ngrid
+            totcloudfrac(ig) = 0.
+            norm=0.
+            DO l=1,nlay
+               ponder=zqice(ig,l)*(pplev(ig,l) - pplev(ig,l+1))
+               totcloudfrac(ig) = totcloudfrac(ig) 
+     &                   + cloudfrac(ig,l)*ponder
+               norm=norm+ponder
+            ENDDO
+            totcloudfrac(ig)=MAX(totcloudfrac(ig)/norm,1.e-12) ! min value if NaNs
+         ENDDO
+c-----No sub-grid cloud representation (CLFvarying=false)
+      ELSE
+         DO l=1,nlay
+            DO ig=1,ngrid
+               zteff(ig,l)=pt(ig,l)
+            END DO
+         END DO
+      END IF ! end if (CLFvarying)
+
+c------------------------------------------------------------------
+c Time subsampling for microphysics 
+      rhocloud(1:ngrid,1:nlay) = rho_dust
+      DO microstep=1,imicro 
+      
+c-------------------------------------------------------------------
+c   1.  Tendencies: 
+c------------------
+
+
+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
+          ENDDO
+        ENDDO
+c------ Tracers tendencies subpdq
+c------ At each micro timestep we add pdq in order to have a stepped entry
+        IF (microphys) THEN
+          DO l=1,nlay
+            DO ig=1,ngrid
+              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_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
+        ENDIF
+        DO l=1,nlay
+          DO ig=1,ngrid
+            subpdq(ig,l,igcm_h2o_ice) = 
+     &          subpdq(ig,l,igcm_h2o_ice)
+     &        + pdq(ig,l,igcm_h2o_ice)
+            subpdq(ig,l,igcm_h2o_vap) = 
+     &          subpdq(ig,l,igcm_h2o_vap)
+     &        + pdq(ig,l,igcm_h2o_vap)
+          ENDDO
+        ENDDO
+c------ Effective tracers quantities in the cloud fraction
+        IF (CLFvarying) THEN     
+           pqeff(:,:,:)=pq(:,:,:) ! prevent from buggs (A. Pottier)
+           pqeff(:,:,igcm_ccn_mass) =pq(:,:,igcm_ccn_mass)/
+     &                                cloudfrac(:,:)
+           pqeff(:,:,igcm_ccn_number)= 
+     &     pq(:,:,igcm_ccn_number)/cloudfrac(:,:)
+           pqeff(:,:,igcm_h2o_ice)= pq(:,:,igcm_h2o_ice)/
+     &                               cloudfrac(:,:)
+        ELSE
+           pqeff(:,:,:)=pq(:,:,:)
+           pqeff(:,:,igcm_ccn_mass)= pq(:,:,igcm_ccn_mass)
+           pqeff(:,:,igcm_ccn_number)= pq(:,:,igcm_ccn_number)
+           pqeff(:,:,igcm_h2o_ice)= pq(:,:,igcm_h2o_ice)
+        END IF      
+        
+c-------------------------------------------------------------------
+c   2.  Main call to the different cloud schemes:
+c------------------------------------------------
+        IF (microphys) THEN
+           CALL improvedclouds(ngrid,nlay,microtimestep,
+     &             pplay,zteff,subpdt, 
+     &             pqeff,subpdq,subpdqcloud,subpdtcloud,
+     &             nq,tauscaling)
+
+        ELSE
+           CALL simpleclouds(ngrid,nlay,microtimestep,
+     &             pplay,pzlay,zteff,subpdt,
+     &             pqeff,subpdq,subpdqcloud,subpdtcloud,
+     &             nq,tau,rice)
+        ENDIF
+
+c-------------------------------------------------------------------
+c   3.  Updating tendencies after cloud scheme:
+c-----------------------------------------------
+
+        IF (microphys) THEN
+          DO l=1,nlay
+            DO ig=1,ngrid
+              subpdq(ig,l,igcm_dust_mass) =
+     &            subpdq(ig,l,igcm_dust_mass)
+     &          + subpdqcloud(ig,l,igcm_dust_mass)
+              subpdq(ig,l,igcm_dust_number) =
+     &            subpdq(ig,l,igcm_dust_number)
+     &          + subpdqcloud(ig,l,igcm_dust_number)
+              subpdq(ig,l,igcm_ccn_mass) =
+     &            subpdq(ig,l,igcm_ccn_mass)
+     &          + subpdqcloud(ig,l,igcm_ccn_mass)
+              subpdq(ig,l,igcm_ccn_number) =
+     &            subpdq(ig,l,igcm_ccn_number)
+     &          + subpdqcloud(ig,l,igcm_ccn_number)
+            ENDDO
+          ENDDO
+        ENDIF
+        DO l=1,nlay
+          DO ig=1,ngrid
+            subpdq(ig,l,igcm_h2o_ice) =
+     &          subpdq(ig,l,igcm_h2o_ice)
+     &        + subpdqcloud(ig,l,igcm_h2o_ice)
+            subpdq(ig,l,igcm_h2o_vap) =
+     &          subpdq(ig,l,igcm_h2o_vap)
+     &        + subpdqcloud(ig,l,igcm_h2o_vap)
+          ENDDO
+        ENDDO
+        
+        IF (activice) THEN
+          DO l=1,nlay
+            DO ig=1,ngrid
+              subpdt(ig,l) =
+     &            subpdt(ig,l) + subpdtcloud(ig,l)
+            ENDDO
+          ENDDO
+        ENDIF
+     
+ 
+      ENDDO ! of DO microstep=1,imicro
+      
+c-------------------------------------------------------------------
+c   6.  Compute final tendencies after time loop:
+c------------------------------------------------
+
+c------ Temperature tendency pdtcloud
+       DO l=1,nlay
+         DO ig=1,ngrid
+             pdtcloud(ig,l) =
+     &         subpdt(ig,l)/real(imicro)-pdt(ig,l)
+          ENDDO
+       ENDDO
+       
+c------ Tracers tendencies pdqcloud
+       DO l=1,nlay
+         DO ig=1,ngrid
+            pdqcloud(ig,l,igcm_h2o_ice) = 
+     &        subpdq(ig,l,igcm_h2o_ice)/real(imicro) 
+     &       - pdq(ig,l,igcm_h2o_ice)
+            pdqcloud(ig,l,igcm_h2o_vap) = 
+     &        subpdq(ig,l,igcm_h2o_vap)/real(imicro) 
+     &       - pdq(ig,l,igcm_h2o_vap)
+         ENDDO
+       ENDDO
+       
+       IF(microphys) THEN
+        DO l=1,nlay
+         DO ig=1,ngrid
+            pdqcloud(ig,l,igcm_ccn_mass) = 
+     &        subpdq(ig,l,igcm_ccn_mass)/real(imicro) 
+     &       - pdq(ig,l,igcm_ccn_mass)
+            pdqcloud(ig,l,igcm_ccn_number) = 
+     &        subpdq(ig,l,igcm_ccn_number)/real(imicro) 
+     &       - pdq(ig,l,igcm_ccn_number)
+         ENDDO
+        ENDDO
+       ENDIF
+       
+       IF(scavenging) THEN
+        DO l=1,nlay
+         DO ig=1,ngrid
+            pdqcloud(ig,l,igcm_dust_mass) = 
+     &        subpdq(ig,l,igcm_dust_mass)/real(imicro) 
+     &       - pdq(ig,l,igcm_dust_mass)
+            pdqcloud(ig,l,igcm_dust_number) = 
+     &        subpdq(ig,l,igcm_dust_number)/real(imicro)
+     &       - pdq(ig,l,igcm_dust_number)
+         ENDDO
+        ENDDO
+       ENDIF
+
+c------- Due to stepped entry, other processes tendencies can add up to negative values
+c------- Therefore, enforce positive values and conserve mass
+       IF(microphys) THEN
+        DO l=1,nlay
+         DO ig=1,ngrid
+          IF ((pq(ig,l,igcm_ccn_number) + 
+     &      ptimestep* (pdq(ig,l,igcm_ccn_number) + 
+     &        pdqcloud(ig,l,igcm_ccn_number)) .le. 1.)
+     &   .or. (pq(ig,l,igcm_ccn_mass) + 
+     &      ptimestep* (pdq(ig,l,igcm_ccn_mass) + 
+     &        pdqcloud(ig,l,igcm_ccn_mass)) .le. 1.e-20)) THEN
+         pdqcloud(ig,l,igcm_ccn_number) =
+     &     - pq(ig,l,igcm_ccn_number)/ptimestep 
+     &     - pdq(ig,l,igcm_ccn_number) + 1.
+         pdqcloud(ig,l,igcm_dust_number) =  
+     &     -pdqcloud(ig,l,igcm_ccn_number)
+         pdqcloud(ig,l,igcm_ccn_mass) =
+     &     - pq(ig,l,igcm_ccn_mass)/ptimestep
+     &     - pdq(ig,l,igcm_ccn_mass) + 1.e-20
+         pdqcloud(ig,l,igcm_dust_mass) = 
+     &     -pdqcloud(ig,l,igcm_ccn_mass)
+          ENDIF
+         ENDDO
+        ENDDO
+       ENDIF
+
+       IF(scavenging) THEN
+        DO l=1,nlay
+         DO ig=1,ngrid
+          IF ((pq(ig,l,igcm_dust_number) + 
+     &      ptimestep* (pdq(ig,l,igcm_dust_number) + 
+     &        pdqcloud(ig,l,igcm_dust_number)) .le. 1.)
+     &   .or. (pq(ig,l,igcm_dust_mass) + 
+     &      ptimestep* (pdq(ig,l,igcm_dust_mass) + 
+     &        pdqcloud(ig,l,igcm_dust_mass)) .le. 1.e-20)) THEN
+         pdqcloud(ig,l,igcm_dust_number) =
+     &     - pq(ig,l,igcm_dust_number)/ptimestep 
+     &     - pdq(ig,l,igcm_dust_number) + 1.
+         pdqcloud(ig,l,igcm_ccn_number) =  
+     &     -pdqcloud(ig,l,igcm_dust_number)
+         pdqcloud(ig,l,igcm_dust_mass) =
+     &     - pq(ig,l,igcm_dust_mass)/ptimestep
+     &     - pdq(ig,l,igcm_dust_mass) + 1.e-20
+         pdqcloud(ig,l,igcm_ccn_mass) = 
+     &     -pdqcloud(ig,l,igcm_dust_mass)
+          ENDIF
+         ENDDO
+        ENDDO
+       ENDIF
+
+        DO l=1,nlay
+         DO ig=1,ngrid
+          IF (pq(ig,l,igcm_h2o_ice) + ptimestep*
+     &       (pdq(ig,l,igcm_h2o_ice) + pdqcloud(ig,l,igcm_h2o_ice)) 
+     &       .le. 1.e-8) THEN
+           pdqcloud(ig,l,igcm_h2o_ice) = 
+     &     - pq(ig,l,igcm_h2o_ice)/ptimestep - pdq(ig,l,igcm_h2o_ice)
+           pdqcloud(ig,l,igcm_h2o_vap) = -pdqcloud(ig,l,igcm_h2o_ice)
+          ENDIF 
+          IF (pq(ig,l,igcm_h2o_vap) + ptimestep*
+     &       (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap)) 
+     &       .le. 1.e-8) THEN
+           pdqcloud(ig,l,igcm_h2o_vap) = 
+     &     - pq(ig,l,igcm_h2o_vap)/ptimestep - pdq(ig,l,igcm_h2o_vap)
+           pdqcloud(ig,l,igcm_h2o_ice) = -pdqcloud(ig,l,igcm_h2o_vap)
+          ENDIF
+         ENDDO
+        ENDDO
+
+
+c------Update the ice and dust particle size "rice" for output or photochemistry
+c------Only rsedcloud is used for the water cycle
+
+      IF(scavenging) THEN 
+        DO l=1, nlay
+         DO ig=1,ngrid
+
+        call updaterdust(
+     &    pq(ig,l,igcm_dust_mass) +                   ! dust mass
+     &   (pdq(ig,l,igcm_dust_mass) +                  ! dust mass
+     &    pdqcloud(ig,l,igcm_dust_mass))*ptimestep,   ! dust mass
+     &    pq(ig,l,igcm_dust_number) +                 ! dust number
+     &   (pdq(ig,l,igcm_dust_number) +                ! dust number
+     &    pdqcloud(ig,l,igcm_dust_number))*ptimestep, ! dust number
+     &    rdust(ig,l))
+
+         ENDDO
+        ENDDO
+      ENDIF
+
+      IF(microphys) THEN
+
+       ! In case one does not want to allow supersatured water when using microphysics.
+       ! Not done by default.
+       IF(.not.supersat) THEN     
+        zt  = pt + (pdt+pdtcloud)*ptimestep
+        call watersat(ngrid*nlay,zt,pplay,zqsat)
+        DO l=1, nlay
+         DO ig=1,ngrid
+          IF (pq(ig,l,igcm_h2o_vap)
+     &      + (pdq(ig,l,igcm_h2o_vap) + pdqcloud(ig,l,igcm_h2o_vap))
+     &      * ptimestep .ge. zqsat(ig,l)) THEN
+             pdqcloud(ig,l,igcm_h2o_vap) = 
+     &         (zqsat(ig,l) - pq(ig,l,igcm_h2o_vap))/ptimestep
+     &        - pdq(ig,l,igcm_h2o_vap)
+             pdqcloud(ig,l,igcm_h2o_ice) = 
+     &         -pdqcloud(ig,l,igcm_h2o_vap)
+             ! no need to correct ccn_number, updaterad can handle this properly.
+          ENDIF
+         ENDDO
+        ENDDO        
+       ENDIF
+       
+       DO l=1, nlay
+         DO ig=1,ngrid
+
+        call updaterice_micro(
+     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
+     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
+     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
+     &    pq(ig,l,igcm_ccn_mass) +                   ! ccn mass
+     &   (pdq(ig,l,igcm_ccn_mass) +                  ! ccn mass
+     &    pdqcloud(ig,l,igcm_ccn_mass))*ptimestep,   ! ccn mass
+     &    pq(ig,l,igcm_ccn_number) +                 ! ccn number
+     &   (pdq(ig,l,igcm_ccn_number) +                ! ccn number
+     &    pdqcloud(ig,l,igcm_ccn_number))*ptimestep, ! ccn number
+     &    tauscaling(ig),rice(ig,l),rhocloud(ig,l))
+          
+         ENDDO
+       ENDDO
+        
+      ELSE ! no microphys
+        
+        DO l=1,nlay
+          DO ig=1,ngrid
+          
+        call updaterice_typ(
+     &    pq(ig,l,igcm_h2o_ice) +                    ! ice mass
+     &   (pdq(ig,l,igcm_h2o_ice) +                   ! ice mass
+     &    pdqcloud(ig,l,igcm_h2o_ice))*ptimestep,    ! ice mass
+     &    tau(ig,1),pzlay(ig,l),rice(ig,l)) 
+
+          ENDDO
+         ENDDO
+       
+       ENDIF ! of IF(microphys)
+      
+      
+      
+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)))
+     &    rice(ig,2)=rice(ig,3) 
+          rice(ig,1)=rice(ig,2)
+        end if
+      end do
+       
+       
+       DO l=1,nlay
+         DO ig=1,ngrid
+           rsedcloud(ig,l)=max(rice(ig,l)*
+     &                 (1.+nuice_sed)*(1.+nuice_sed)*(1.+nuice_sed),
+     &                    rdust(ig,l))
+!          rsedcloud(ig,l)=min(rsedcloud(ig,l),1.e-4)
+         ENDDO
+       ENDDO
+       
+! used for rad. transfer calculations
+! nuice is constant because a lognormal distribution is prescribed
+      nuice(1:ngrid,1:nlay)=nuice_ref 
+
+c------Update tendencies for sub-grid water ice clouds
+      IF (CLFvarying) THEN
+        DO ig=1,ngrid
+          DO l=1,nlay
+            pdqcloud(ig,l,igcm_dust_mass)=pdqcloud(ig,l,igcm_dust_mass)
+     &          *cloudfrac(ig,l)
+            pdqcloud(ig,l,igcm_ccn_mass)=pdqcloud(ig,l,igcm_ccn_mass)
+     &          *cloudfrac(ig,l)
+            pdqcloud(ig,l,igcm_dust_number)=pdqcloud(ig,l,
+     &           igcm_dust_number) *cloudfrac(ig,l)
+            pdqcloud(ig,l,igcm_ccn_number)=pdqcloud(ig,l,
+     &           igcm_ccn_number) *cloudfrac(ig,l)
+            pdqcloud(ig,l,igcm_h2o_vap)=pdqcloud(ig,l,
+     &           igcm_h2o_vap) *cloudfrac(ig,l)
+            pdqcloud(ig,l,igcm_h2o_ice)=pdqcloud(ig,l,
+     &           igcm_h2o_ice) *cloudfrac(ig,l)
+          ENDDO
+        ENDDO   
+        pdtcloud(:,:)=pdtcloud(:,:)*cloudfrac(:,:)
+      ENDIF
+c=======================================================================
+      call WRITEDIAGFI(ngrid,"pdqice2","pdqcloudice apres microphysique"
+     &      ,"kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,igcm_h2o_ice))
+      call WRITEDIAGFI(ngrid,"pdqvap2","pdqcloudvap apres microphysique"
+     &      ,"kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
+     &      igcm_h2o_vap))
+      call WRITEDIAGFI(ngrid,"pdqccn2","pdqcloudccn apres microphysique"
+     &      ,"kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
+     &      igcm_ccn_mass)) 
+      call WRITEDIAGFI(ngrid,"pdqccnN2","pdqcloudccnN apres 
+     &      microphysique","nb/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
+     &      igcm_ccn_number))
+      call WRITEDIAGFI(ngrid,"pdqdust2", "pdqclouddust apres
+     &      microphysique","kg/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
+     &      igcm_dust_mass))
+      call WRITEDIAGFI(ngrid,"pdqdustN2", "pdqclouddustN apres
+     &      microphysique","nb/kg.s-1",3,pdqcloud(1:ngrid,1:nlay,
+     &      igcm_dust_number)) 
+
+
+c=======================================================================
+
+      END SUBROUTINE watercloud
+      
+      END MODULE watercloud_mod 
