Index: /trunk/LMDZ.GENERIC/README
===================================================================
--- /trunk/LMDZ.GENERIC/README	(revision 725)
+++ /trunk/LMDZ.GENERIC/README	(revision 726)
@@ -688,2 +688,19 @@
 - Some cleanup in start2archive.F and ini_archive.F to get them to work for
   the generic model (removed some "Martian" specificities).
+  
+== 17/07/2012 == JL for LK
+- Generalization of aerosol scheme:
+    - any number of aerosols can be used and id numbers are determined consistently by the code. Aerosol order
+    not important anymore.
+    - addition of a module with the id numbers for aerosols (aerosol_mod.F90).
+    - initialization of aerosols id numbers in iniaerosol.F90 
+    - compile with -s x where x *must* be equal to the number of aerosols turned on in callphys.def (either by a
+    flag or by dusttau>0 for dust).
+      => may have to erase object files when compiling with s option for the first time.
+- For no aerosols, run with aeroco2=.true. and aerofixco2=.true (the default distribution for fixed co2
+aerosols is 1.e-9; can be changed in aeropacity).
+- If starting from an old start file, recreate start file with the q=0 option in newstart.e.
+- update callphys.def with aeroXXX and aerofixXXX options (only XXX=co2,h2o supported for
+now). Dust is activated by setting dusttau>0. See the early mars case in deftank. 
+- To add other aerosols, see Laura Kerber.
+
Index: /trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 725)
+++ /trunk/LMDZ.GENERIC/libf/phystd/aeropacity.F90	(revision 726)
@@ -2,5 +2,6 @@
          aerosol,reffrad,QREFvis3d,QREFir3d,tau_col,cloudfrac,totcloudfrac,clearsky)
 
-       use radinc_h, only : naerkind, L_TAUMAX
+       use radinc_h, only : L_TAUMAX,naerkind
+       use aerosol_mod
                   
        implicit none
@@ -43,4 +44,5 @@
 #include "comgeomfi.h"
 #include "tracer.h"
+#include "comvert.h"
 
       INTEGER ngrid,nlayer,nq
@@ -76,5 +78,5 @@
       real topdust, expfactor, zp
       REAL taudusttmp(ngridmx) ! Temporary dust opacity used before scaling
-
+      REAL tauh2so4tmp(ngridmx) ! Temporary h2so4 opacity used before scaling
       ! BENJAMIN MODIFS
       real CLFtot,CLF1,CLF2
@@ -100,30 +102,37 @@
         endif
 
+        write(*,*) "Tracers found in aeropacity:"
         do iq=1,nqmx
           tracername=noms(iq)
           if (tracername.eq."co2_ice") then
             i_co2ice=iq
+          write(*,*) "i_co2ice=",i_co2ice
+
           endif
           if (tracername.eq."h2o_ice") then
             i_h2oice=iq
+            write(*,*) "i_h2oice=",i_h2oice
           endif
         enddo
-        
-        write(*,*) "aeropacity: i_co2ice=",i_co2ice
-        write(*,*) "            i_h2oice=",i_h2oice
-        
-      if(watercond.and.(.not.aerofixed))then
-         if(naerkind.lt.2)then
-            print*,'Cannot have active H2O clouds with naerkind less than 2!'
-            call abort
-         endif
-      endif
-
-      if(dusttau.gt.0.0)then
-         if(naerkind.lt.3)then
-            print*,'Cannot have active dust with naerkind less than 3!'
-            call abort
-         endif
-      endif
+
+        if (naerkind.ne.0) then
+          print*, "If you would like to use aerosols, make sure any old"
+          print*, "start files are updated in newstart using the option"
+          print*, "q=0"
+          write(*,*) "Aerosols found in aeropacity:"
+        endif
+
+        if (iaero_co2.ne.0) then
+          print*, 'iaero_co2=  ',iaero_co2
+        endif
+        if (iaero_h2o.ne.0) then
+          print*,'iaero_h2o=  ',iaero_h2o    
+        endif
+        if (iaero_dust.ne.0) then
+          print*,'iaero_dust= ',iaero_dust
+        endif
+        if (iaero_h2so4.ne.0) then
+          print*,'iaero_h2so4= ',iaero_h2so4
+        endif
 
         firstcall=.false.
@@ -131,16 +140,15 @@
 
 
-      DO iaer = 1, naerkind ! Loop on aerosol kind (NOT tracer)
 !     ---------------------------------------------------------
-         aerkind: SELECT CASE (iaer)
-!==================================================================
-         CASE(1) aerkind                         ! CO2 ice (iaer=1)
-!==================================================================
-
+!==================================================================
+!    CO2 ice aerosols
+!==================================================================
+
+        if (iaero_co2.ne.0) then
+           iaer=iaero_co2
 !       1. Initialization
             aerosol(:,:,iaer)=0.0
-
 !       2. Opacity calculation
-            if (aerofixed.or.(i_co2ice.eq.0)) then !  CO2 ice cloud prescribed
+            if (aerofixco2.or.(i_co2ice.eq.0)) then !  CO2 ice cloud prescribed
                do ig=1, ngrid
                   do l=1,nlayer
@@ -162,5 +170,5 @@
                      aerosol(ig,l,iaer) = aerosol0
 !                     aerosol(ig,l,iaer) = 0.0
-
+!                     print*, aerosol(ig,l,iaer)
 !        using cloud fraction
 !                     aerosol(ig,l,iaer) = -log(1 - CLF + CLF*exp(-aerosol0/CLF))
@@ -170,15 +178,16 @@
                   ENDDO
                ENDDO
-            end if
-            
-!==================================================================
-         CASE(2) aerkind              ! Water ice / liquid (iaer=2)
-!==================================================================
-
+            end if ! if fixed or varying
+        end if ! if CO2 aerosols   
+!==================================================================
+!     Water ice / liquid 
+!==================================================================
+
+        if (iaero_h2o.ne.0) then 
+           iaer=iaero_h2o
 !       1. Initialization
             aerosol(:,:,iaer)=0.0
-
 !       2. Opacity calculation
-            if (aerofixed.or.(i_h2oice.eq.0).or.clearsky) then
+            if (aerofixh2o.or.(i_h2oice.eq.0).or.clearsky) then
                do ig=1,ngrid ! to stop the rad tran bug
                   do l=1,nlayer
@@ -252,33 +261,121 @@
                enddo
             end if
-
-
-!==================================================================
-         CASE(3) aerkind                            ! Dust (iaer=3)
-!==================================================================
+        end if ! End if h2o aerosol
+
+!==================================================================
+!             Dust 
+!==================================================================
+        if (iaero_dust.ne.0) then
+          iaer=iaero_dust
+
+!         1. Initialization 
+            aerosol(:,:,iaer)=0.0
+          
+            topdust=30.0 ! km  (used to be 10.0 km) LK
+
+!       2. Opacity calculation
+
+!           expfactor=0.
+           DO l=1,nlayer-1
+             DO ig=1,ngrid
+!             Typical mixing ratio profile
+
+                 zp=(preff/pplay(ig,l))**(70./topdust)
+                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
+
+!             Vertical scaling function
+              aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
+               *expfactor
+
+
+             ENDDO
+           ENDDO
+
+!          Rescaling each layer to reproduce the choosen (or assimilated)
+!          dust extinction opacity at visible reference wavelength, which
+!          is scaled to the "preff" reference surface pressure available in comvert.h
+!          and stored in startfi.nc
+
+            taudusttmp(1:ngrid)=0.
+              DO l=1,nlayer
+                DO ig=1,ngrid
+                   taudusttmp(ig) = taudusttmp(ig) &
+                          +  aerosol(ig,l,iaer)
+                ENDDO
+              ENDDO
+            DO l=1,nlayer-1
+               DO ig=1,ngrid
+                  aerosol(ig,l,iaer) = max(1E-20, &
+                          dusttau &
+                       *  pplev(ig,1) / preff &
+                       *  aerosol(ig,l,iaer) &
+                       /  taudusttmp(ig))
+
+              ENDDO
+            ENDDO
+        end if ! If dust aerosol   
+
+!==================================================================
+!           H2SO4 
+!==================================================================
+! added by LK
+      if (iaero_h2so4.ne.0) then
+          iaer=iaero_h2so4
 
 !       1. Initialization
             aerosol(:,:,iaer)=0.0
 
-            topdust=10.0 ! km
-
-            print*,'WARNING, dust is experimental for Early Mars'
 
 !       2. Opacity calculation
 
-            do ig=1,ngrid
-               do l=1,nlayer
-                  zp=700./pplay(ig,l)
-                  aerosol(ig,l,iaer)=(dusttau/700.)*(pplev(ig,l)-pplev(ig,l+1)) &
-                       *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 ) &
-                       *QREFvis3d(ig,l,iaer) / QREFvis3d(ig,1,iaer) 
-                  ! takes into account particle variation
-              enddo
-            enddo
-
-!==================================================================
-         END SELECT aerkind
-      ENDDO ! iaer (loop on aerosol kind)
-
+!           expfactor=0.
+           DO l=1,nlayer-1
+             DO ig=1,ngrid
+!             Typical mixing ratio profile
+
+                 zp=(preff/pplay(ig,l))**(70./30) !emulating topdust
+                 expfactor=max(exp(0.007*(1.-max(zp,1.))),1.e-3)
+
+!             Vertical scaling function
+              aerosol(ig,l,iaer)= (pplev(ig,l)-pplev(ig,l+1)) &
+               *expfactor
+
+
+             ENDDO
+           ENDDO
+            tauh2so4tmp(1:ngrid)=0.
+              DO l=1,nlayer
+                DO ig=1,ngrid
+                   tauh2so4tmp(ig) = tauh2so4tmp(ig) &
+                          +  aerosol(ig,l,iaer)
+                ENDDO
+              ENDDO
+            DO l=1,nlayer-1
+               DO ig=1,ngrid
+                  aerosol(ig,l,iaer) = max(1E-20, &
+                          1 &
+                       *  pplev(ig,1) / preff &
+                       *  aerosol(ig,l,iaer) &
+                       /  tauh2so4tmp(ig))
+
+              ENDDO
+            ENDDO
+
+! 1/700. is assuming a "sulfurtau" of 1
+! Sulfur aerosol routine to be improved.
+!                     aerosol0 =                         &
+!                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
+!                          ( rho_h2so4 * reffrad(ig,l,iaer) )  ) *   &
+!                          ( pq(ig,l,i_h2so4) + 1.E-9 ) *         &
+!                          ( pplev(ig,l) - pplev(ig,l+1) ) / g
+!                     aerosol0           = max(aerosol0,1.e-9)
+!                     aerosol0           = min(aerosol0,L_TAUMAX)
+!                     aerosol(ig,l,iaer) = aerosol0
+
+!                  ENDDO
+!               ENDDO
+      end if
+ 
+           
 
 ! --------------------------------------------------------------------------
@@ -316,5 +413,4 @@
          endif
       end do
-
       return
     end subroutine aeropacity
Index: /trunk/LMDZ.GENERIC/libf/phystd/aeroptproperties.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/aeroptproperties.F90	(revision 725)
+++ /trunk/LMDZ.GENERIC/libf/phystd/aeroptproperties.F90	(revision 726)
@@ -5,5 +5,5 @@
 !                                  omegaREFvis3d,omegaREFir3d)
 
-      use radinc_h,    only: L_NSPECTI,L_NSPECTV,naerkind,nsizemax
+      use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind
       use radcommon_h, only: QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir
       use radcommon_h, only: qrefvis,qrefir,omegarefvis,omegarefir
@@ -208,5 +208,5 @@
 
           if (firstcall) then
-             print*,'Optical properties of the aerosol are homogenous for:'
+             print*,'Optical prop. of the aerosol are homogenous for:'
              print*,'iaer = ',iaer
           endif
Index: /trunk/LMDZ.GENERIC/libf/phystd/aerosol_mod.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/aerosol_mod.F90	(revision 726)
+++ /trunk/LMDZ.GENERIC/libf/phystd/aerosol_mod.F90	(revision 726)
@@ -0,0 +1,14 @@
+!==================================================================
+module aerosol_mod
+!==================================================================
+
+!  aerosol indexes: these are initialized in iniaerosol.F and should be 0 if the
+!                 corresponding aerosol was not activated in callphys.def
+      integer :: iaero_co2 
+      integer :: iaero_h2o
+      integer :: iaero_dust
+      integer :: iaero_h2so4
+
+!==================================================================
+end module aerosol_mod
+!==================================================================
Index: /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 725)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90	(revision 726)
@@ -14,4 +14,5 @@
       use ioipsl_getincom 
       use gases_h
+      use aerosol_mod
 
       implicit none
@@ -189,5 +190,6 @@
       REAL*8 muvarrad(L_LEVELS)
 
-!=======================================================================
+         
+!===============================================================
 !     Initialization on first call
 
@@ -211,45 +213,44 @@
 !     Effective radius and variance of the aerosols
 
-         do iaer=1,naerkind
+
 !     these values will change once the microphysics gets to work
 !     UNLESS tracer=.false., in which case we should be working with
 !     a fixed aerosol layer, and be able to define reffrad in a 
 !     .def file. To be improved!
-
-            if(iaer.eq.1)then ! CO2 ice
-               do l=1,nlayer
-                  do ig=1,ngrid
+!     Generalization of aerosols added by LK
+            do iaer = 1,naerkind
+              if (iaer.eq.iaero_co2) then ! active CO2 aerosols  
+                do l=1,nlayer
+                   do ig=1,ngrid
                      reffrad(ig,l,iaer)  = 1.e-4
                      nueffrad(ig,l,iaer) = 0.1 
-                  enddo
-               enddo
-            endif
-
-            if(iaer.eq.2)then ! H2O ice
-               do l=1,nlayer
-                  do ig=1,ngrid
+                   enddo
+                enddo
+              endif
+              if (iaer.eq.iaero_h2o) then ! active H2O aerosols
+                do l=1,nlayer
+                   do ig=1,ngrid
                      reffrad(ig,l,iaer)  = 1.e-5
                      nueffrad(ig,l,iaer) = 0.1 
-                  enddo
-               enddo
-            endif
-
-            if(iaer.eq.3)then ! dust
-               do l=1,nlayer
-                  do ig=1,ngrid
+                   enddo
+                enddo
+              endif
+              if(iaer.eq.iaero_dust)then ! active dust
+                do l=1,nlayer
+                   do ig=1,ngrid
                      reffrad(ig,l,iaer)  = 1.e-5
                      nueffrad(ig,l,iaer) = 0.1 
-                  enddo
-               enddo
-            endif
- 
-            if(iaer.gt.3)then
-               print*,'Error in callcorrk, naerkind is too high.'
-               print*,'The code still needs generalisation to arbitrary'
-               print*,'aerosol kinds and number.'
-               call abort
-            endif
-
-         enddo
+                   enddo
+                enddo
+              endif
+              if(iaer.eq.iaero_h2so4)then ! Active h2so4 aerosols
+                do l=1,nlayer
+                   do ig=1,ngrid
+                     reffrad(ig,l,iaer)  = 1.e-6
+                     nueffrad(ig,l,iaer) = 0.1
+                   enddo
+                enddo
+              endif
+            enddo
 
          print*, "callcorrk: Correlated-k data base folder:",trim(datadir)
@@ -313,101 +314,103 @@
 
 
-
-      if(radfixed)then
-         do l=1,nlayer
+      do iaer=1,naerkind
+
+       if(iaer.eq.iaero_co2)then
+         if(radfixed)then
+           do l=1,nlayer
             do ig=1,ngrid
-               reffrad(ig,l,1)    = 5.e-5 ! CO2 ice
+               reffrad(ig,l,iaer)    = 5.e-5 ! CO2 ice
             enddo
-         enddo
-         print*,'CO2 ice particle size = ',reffrad(1,1,1)/1.e-6,' um'
-         if(naerkind.ge.2)then
-            do l=1,nlayer
-               do ig=1,ngrid
-                  !reffrad(ig,l,2) = 2.e-5 ! H2O ice
-                   reffrad(ig,l,2) = rad_chaud ! H2O ice
-
-                   zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
-                   zfice = MIN(MAX(zfice,0.0),1.0)
-                   reffrad(ig,l,2)= rad_chaud * (1.-zfice) + rad_froid * zfice
-                   nueffrad(ig,l,2) = coef_chaud * (1.-zfice) + coef_froid * zfice
-               enddo
-            enddo
-            print*,'H2O ice particle size = ',reffrad(1,1,2)/1.e-6,' um'
-         endif
-         if(naerkind.eq.3)then
-            do l=1,nlayer
-               do ig=1,ngrid
-                  reffrad(ig,l,naerkind) = 2.e-6 ! dust
-               enddo
-            enddo
-            print*,'Dust particle size = ',reffrad(1,1,naerkind)/1.e-6,' um'
-         endif
-         if(naerkind.gt.3)then
-            print*,'Code not general enough to deal with naerkind > 3 yet.'
-            call abort
-         endif
-
-      else
-
-         maxrad=0.0
-         !averad=0.0
-         minrad=1.0
-         do l=1,nlayer
-
-            !masse = (pplev(ig,l) - pplev(ig,l+1))/g
+           enddo
+           print*,'CO2 ice particle size =',reffrad(1,1,iaer)/1.e-6,'um'
+         else 
+           maxrad=0.0
+           !averad=0.0
+           minrad=1.0
+           do l=1,nlayer
+
+             !masse = (pplev(ig,l) - pplev(ig,l+1))/g
 
             do ig=1,ngrid
-               !if(tracer)then
                if(tracer.and.igcm_co2_ice.gt.0)then
 
                   if(igcm_co2_ice.lt.1)then
-                     print*,'Tracers but no CO2 ice still seems to be a problem...'
+                     print*,'Tracers but no CO2 ice still'
+                     print*,'seems to be a problem...'
                      print*,'Aborting in callcorrk.'
                      stop
                   endif
-
-                  reffrad(ig,l,1) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ &
+                  reffrad(ig,l,iaer) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ &
                        (4*Nmix_co2*pi*rho_co2) )
                endif
-               reffrad(ig,l,1) = max(reffrad(ig,l,1),1.e-6)
-               reffrad(ig,l,1) = min(reffrad(ig,l,1),500.e-6)
-
-               !averad = averad + reffrad(ig,l,1)*area(ig)
-               maxrad = max(reffrad(ig,l,1),maxrad)
-               minrad = min(reffrad(ig,l,1),minrad)
+               reffrad(ig,l,iaer) = max(reffrad(ig,l,iaer),1.e-6)
+               reffrad(ig,l,iaer) = min(reffrad(ig,l,iaer),500.e-6)
+               !averad = averad + reffrad(ig,l,iaer)*area(ig)
+               maxrad = max(reffrad(ig,l,iaer),maxrad)
+               minrad = min(reffrad(ig,l,iaer),minrad)
             enddo
-         enddo
-         if(igcm_co2_ice.gt.0)then
-            print*,'Max. CO2 ice particle size = ',maxrad/1.e-6,' um'
-            print*,'Min. CO2 ice particle size = ',minrad/1.e-6,' um'
+           enddo
+           if(igcm_co2_ice.gt.0)then
+             print*,'Max. CO2 ice particle size = ',maxrad/1.e-6,' um'
+             print*,'Min. CO2 ice particle size = ',minrad/1.e-6,' um'
+           endif
          endif
-
-         if((naerkind.ge.2).and.water)then
-            maxrad=0.0
-            minrad=1.0
+       endif
+
+
+       if(iaer.eq.iaero_h2o)then
+          if(radfixed) then
             do l=1,nlayer
                do ig=1,ngrid
-                  reffrad(ig,l,2) = CBRT( 3*pq(ig,l,igcm_h2o_ice)/ &
-                       (4*Nmix_h2o*pi*rho_ice) )
-                  reffrad(ig,l,2) = max(reffrad(ig,l,2),1.e-6)
-                  reffrad(ig,l,2) = min(reffrad(ig,l,2),100.e-6)
-
-                  maxrad = max(reffrad(ig,l,2),maxrad)
-                  minrad = min(reffrad(ig,l,2),minrad)
+                  !reffrad(ig,l,iaer) = 2.e-5 ! H2O ice
+                  reffrad(ig,l,iaer) = rad_chaud ! H2O ice
+
+                  zfice = 1.0 - (pt(ig,l)-T_h2O_ice_clouds) / (T_h2O_ice_liq-T_h2O_ice_clouds)
+                  zfice = MIN(MAX(zfice,0.0),1.0)
+                  reffrad(ig,l,2)= rad_chaud * (1.-zfice) + rad_froid * zfice
+                  nueffrad(ig,l,2) = coef_chaud * (1.-zfice) + coef_froid * zfice
                enddo
             enddo
-            print*,'Max. H2O ice particle size = ',maxrad/1.e-6,' um'
-            print*,'Min. H2O ice particle size = ',minrad/1.e-6,' um'
+            print*,'H2O ice particle size=',reffrad(1,1,iaer)/1.e-6,'um'
+          else
+             if(water)then
+               maxrad=0.0
+               minrad=1.0
+               do l=1,nlayer
+                do ig=1,ngrid
+                  reffrad(ig,l,iaer) = CBRT( 3*pq(ig,l,igcm_h2o_ice)/ &
+                       (4*Nmix_h2o*pi*rho_ice) )
+                  reffrad(ig,l,iaer) = max(reffrad(ig,l,iaer),1.e-6)
+                  reffrad(ig,l,iaer) = min(reffrad(ig,l,iaer),100.e-6)
+
+                  maxrad = max(reffrad(ig,l,iaer),maxrad)
+                  minrad = min(reffrad(ig,l,iaer),minrad)
+                enddo
+               enddo
+               print*,'Max. H2O ice particle size = ',maxrad/1.e-6,' um'
+               print*,'Min. H2O ice particle size = ',minrad/1.e-6,' um'
+             endif
          endif
-
-         if(naerkind.eq.3)then
+       endif
+
+         if(iaer.eq.iaero_dust)then
             do l=1,nlayer
                do ig=1,ngrid
-                  reffrad(ig,l,naerkind) = 2.e-6 ! dust
+                  reffrad(ig,l,iaer) = 2.e-6 ! dust
                enddo
             enddo
+            print*,'Dust particle size = ',reffrad(1,1,iaer)/1.e-6,' um'
          endif
 
-      endif
+         if(iaer.eq.iaero_h2so4)then
+            do l=1,nlayer
+               do ig=1,ngrid
+                  reffrad(ig,l,iaer) = 1.e-6 ! h2so4
+               enddo
+            enddo
+            print*,'H2SO4 particle size =',reffrad(1,1,iaer)/1.e-6,' um'
+         endif
+      enddo ! do iaer=1,naerkind
+         ! added by LK
 
 
@@ -523,5 +526,6 @@
                   do nw=1,L_NSPECTV
                      if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then
-                        print*,'Serious problems with qsvaer values in callcorrk'
+                        print*,'Serious problems with qsvaer values' 
+                        print*,'in callcorrk'
                         call abort
                      endif
@@ -533,5 +537,6 @@
                   do nw=1,L_NSPECTI 
                      if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then
-                        print*,'Serious problems with qsiaer values in callcorrk'
+                        print*,'Serious problems with qsiaer values'
+                        print*,'in callcorrk'
                         call abort
                      endif
Index: /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 725)
+++ /trunk/LMDZ.GENERIC/libf/phystd/callkeys.h	(revision 726)
@@ -34,5 +34,6 @@
      &   , sedimentation,water,watercond,waterrain                      &
      &   , rainthreshold                                                &
-     &   , aerofixed                                                    &
+     &   , aeroco2,aeroh2o,aeroh2so4                                    &
+     &   , aerofixco2,aerofixh2o                                        &
      &   , hydrology                                                    &
      &   , sourceevol                                                   &
@@ -72,5 +73,6 @@
       logical sedimentation
       logical water,watercond,waterrain
-      logical aerofixed
+      logical aeroco2,aeroh2o,aeroh2so4 
+      logical aerofixco2,aerofixh2o
       logical hydrology
       logical sourceevol
Index: /trunk/LMDZ.GENERIC/libf/phystd/iniaerosol.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/iniaerosol.F	(revision 726)
+++ /trunk/LMDZ.GENERIC/libf/phystd/iniaerosol.F	(revision 726)
@@ -0,0 +1,70 @@
+      SUBROUTINE iniaerosol()
+
+
+      use radinc_h, only: naerkind
+      use aerosol_mod
+
+      IMPLICIT NONE
+c=======================================================================
+c   subject:
+c   --------
+c   Initialization related to aerosols 
+c   (CO2 aerosols, dust, water, chemical species, ice...)   
+c
+c   author: Laura Kerber
+c   ------
+c        
+c=======================================================================
+
+
+#include "callkeys.h"
+      integer ia
+
+      ! Initialize aerosol indexes to zero:
+      ! NB: iaero_* indexes are in module aerosol_mod
+
+      iaero_co2   =0
+      iaero_h2o   =0
+      iaero_dust  =0
+      iaero_h2so4 =0
+      
+      ia=0
+      if (aeroco2) then
+      ia=ia+1
+      iaero_co2=ia
+      write(*,*) 'CO2 aerosol = ', iaero_co2
+      endif
+ 
+      if (aeroh2o) then
+      ia=ia+1
+      iaero_h2o=ia
+      write(*,*) 'H2O aerosol = ', iaero_h2o
+      endif
+
+      if (dusttau.gt.0) then
+      ia=ia+1
+      iaero_dust=ia
+      write(*,*) 'Dust aerosol = ', iaero_dust
+      endif
+
+      if (aeroh2so4) then
+      ia=ia+1
+      iaero_h2so4=ia
+      write(*,*) 'H2SO4 aerosol = ', iaero_h2so4
+      endif
+      
+      write(*,*) 'Number of aerosols= ', ia
+      
+      if (ia.eq.0) then  !For the zero aerosol case
+      ia = 1
+      endif
+
+      if (ia.ne.naerkind) then
+          print*, 'Aerosols counted not equal to naerkind'
+          print*, 'Compile with tag -s',ia,'to run'
+          print*, 'or change options in callphys.def'
+          print*, 'Abort in iniaerosol.F'
+      call abort
+      endif
+
+      end
Index: /trunk/LMDZ.GENERIC/libf/phystd/inifis.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 725)
+++ /trunk/LMDZ.GENERIC/libf/phystd/inifis.F	(revision 726)
@@ -367,9 +367,4 @@
 ! TRACERS:
 
-         write(*,*)"Fixed / inactive aerosol distribution?"
-         aerofixed=.true.       ! default value
-         call getin("aerofixed",aerofixed)
-         write(*,*)" aerofixed = ",aerofixed
-
          write(*,*)"Varying H2O cloud fraction?"
          CLFvarying=.false.     ! default value
@@ -392,4 +387,9 @@
          write(*,*)" Nmix_h2o = ",Nmix_h2o
 
+!         write(*,*)"Number of radiatively active aerosols:"
+!         naerkind=0. ! default value
+!         call getin("naerkind",naerkind)
+!         write(*,*)" naerkind = ",naerkind
+
          write(*,*)"Opacity of dust (if used):"
          dusttau=0. ! default value
@@ -397,10 +397,28 @@
          write(*,*)" dusttau = ",dusttau
 
-! Test of incompatibility:
-! if less than 3 aerosols, then dusttau should = 0
-         if ((naerkind.lt.3).and.dusttau.gt.0.) then
-           print*,'Need naer>2 for radiatively active dust!'
-           stop
-         endif
+         write(*,*)"Radiatively active CO2 aerosols?"
+         aeroco2=.false.     ! default value
+         call getin("aeroco2",aeroco2)
+         write(*,*)" aeroco2 = ",aeroco2
+
+         write(*,*)"Fixed CO2 aerosol distribution?"
+         aerofixco2=.false.     ! default value
+         call getin("aerofixco2",aerofixco2)
+         write(*,*)" aerofixco2 = ",aerofixco2
+
+         write(*,*)"Radiatively active water ice?"
+         aeroh2o=.false.     ! default value
+         call getin("aeroh2o",aeroh2o)
+         write(*,*)" aeroh2o = ",aeroh2o
+
+         write(*,*)"Fixed H2O aerosol distribution?"
+         aerofixh2o=.false.     ! default value
+         call getin("aerofixh2o",aerofixh2o)
+         write(*,*)" aerofixh2o = ",aerofixh2o
+
+         write(*,*)"Radiatively active sulfuric acid aersols?"
+         aeroh2so4=.false.     ! default value
+         call getin("aeroh2so4",aeroh2so4)
+         write(*,*)" aeroh2so4 = ",aeroh2so4
 
          write(*,*)"Cloud pressure level (with kastprof only):"
@@ -426,9 +444,18 @@
 
 ! Test of incompatibility:
-! if no tracers, then aerofixed should be true
-         if ((.not.tracer).and.(.not.aerofixed)) then
-           print*,'if tracers are off, aerofixed must be ON!'
+! if no tracers, then the variable aerosols should be turned off
+         if ((.not.tracer).and.(.not.aerofixco2)) then
+           print*,'if tracers are off, you must turn aerofixco2 on.'
+           print*, '(i.e., you must fix the aerosol distribution)'
+           print*, 'This option was formerly known as aerofixed'
            stop
          endif
+         if ((.not.tracer).and.(.not.aerofixh2o)) then
+           print*,'if tracers are off, you must turn aerofixh2o on.'
+           print*, '(i.e., you must fix the aerosol distribution)'
+           print*, 'This option was formerly known as aerofixed'
+           stop
+         endif
+
 
 ! Test of incompatibility:
Index: /trunk/LMDZ.GENERIC/libf/phystd/largescale.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/largescale.F	(revision 725)
+++ /trunk/LMDZ.GENERIC/libf/phystd/largescale.F	(revision 726)
@@ -111,5 +111,6 @@
 
 !     for varying particle size in rad tran and (possibly) sedimentation
-         if(.not.aerofixed)then
+         if(aeroh2o.and.(.not.aerofixh2o)) then
+
             reffH2O(i,k) = CBRT( 3*zcond(i)/( 4*Nmix_h2o*pi*rho_ice))
             reffH2O(i,k) = max(reffH2O(i,k),1.e-8)
Index: /trunk/LMDZ.GENERIC/libf/phystd/physdem1.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/physdem1.F	(revision 725)
+++ /trunk/LMDZ.GENERIC/libf/phystd/physdem1.F	(revision 726)
@@ -5,5 +5,4 @@
      .                   cloudfrac,totcloudfrac,hice)
 
-      use radcommon_h, only: tauvis
  
       implicit none
@@ -215,7 +214,4 @@
       tab_cntrl(33) = dtemisice(1) ! time scale for snow metamorphism (north)
       tab_cntrl(34) = dtemisice(2) ! time scale for snow metamorphism (south)
-
-c dust aerosol properties
-      tab_cntrl(27) = tauvis      ! mean visible optical depth
 
       tab_cntrl(28) = 0. 
Index: /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 725)
+++ /trunk/LMDZ.GENERIC/libf/phystd/physiq.F90	(revision 726)
@@ -7,8 +7,9 @@
                   pdu,pdv,pdt,pdq,pdpsrf,tracerdyn)
  
-      use radinc_h, only : naerkind,L_NSPECTI,L_NSPECTV
+      use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind
       use watercommon_h, only : RLVTT
       use gases_h
       use radcommon_h, only: sigma
+      use aerosol_mod
       implicit none
 
@@ -389,5 +390,5 @@
 !     included by RW for variable H2O particle sizes
       real reffH2O(ngridmx,nlayermx)
-      real reffcol(ngridmx,2)
+      real reffcol(ngridmx,naerkind)
 
 !     included by RW for sourceevol
@@ -420,4 +421,9 @@
          zdtsw(:,:) = 0.0
          zdtlw(:,:) = 0.0
+
+!        initialize aerosol indexes 
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+            call iniaerosol()
+
       
 !        initialize tracer names, indexes and properties
@@ -427,4 +433,6 @@
             call initracer()
          endif ! end tracer
+
+! 
 
 !        read startfi (initial parameters)
@@ -811,5 +819,5 @@
             dEtotLW  = cpp*SUM(massarea(:,:)*zdtlw(:,:))/totarea
             dEtotsSW = SUM(fluxsurf_sw(:)*(1.-albedo(:))*area(:))/totarea
-            dEtotsLW = SUM(fluxsurf_lw(:)*emis(:)*area(:)-zplanck(:))/totarea
+            dEtotsLW = SUM((fluxsurf_lw(:)*emis(:)-zplanck(:))*area(:))/totarea
 	    dEzRadsw(:,:)=cpp*mass(:,:)*zdtsw(:,:)
 	    dEzRadlw(:,:)=cpp*mass(:,:)*zdtlw(:,:)
@@ -1045,5 +1053,5 @@
             call abort
          endif
-
+     print*, 'we are in co2cond!!!!!'
          call condense_cloud(ngrid,nlayer,nq,ptimestep,   &
               capcal,pplay,pplev,tsurf,pt,                  &
@@ -1161,7 +1169,7 @@
                      pdq(ig,l,igcm_h2o_ice) = pdq(ig,l,igcm_h2o_ice)+dqcldlscale(ig,l)
                      pdt(ig,l)              = pdt(ig,l)+dtlscale(ig,l)
-
-                     if(.not.aerofixed)then
-                        reffrad(ig,l,2)=reffH2O(ig,l)
+                     
+                     if(aeroh2o.and.(.not.aerofixh2o)) then
+                        reffrad(ig,l,iaero_h2o)=reffH2O(ig,l)
                      endif
 
@@ -1541,5 +1549,6 @@
          enddo
          close(38)
-         print*,'As testradtimes enabled, exiting physics on first call'
+         print*,'As testradtimes enabled,'
+         print*,'exiting physics on first call'
          call abort
       endif
@@ -1559,16 +1568,16 @@
          enddo
 
-         ! not generalised for arbitrary aerosols yet!!!
+         ! Generalised for arbitrary aerosols now. (LK)
          reffcol(:,:)=0.0
          do ig=1,ngrid
             do l=1,nlayer
-               if(co2cond)then
-                  reffcol(ig,1) = reffcol(ig,1) + zq(ig,l,igcm_co2_ice) * &
-                                  reffrad(ig,l,1) *    &
+               if(co2cond.and.(iaero_co2.ne.0)) then
+                  reffcol(ig,iaero_co2) = reffcol(ig,iaero_co2) + zq(ig,l,igcm_co2_ice) * &
+                                  reffrad(ig,l,iaero_co2) *    &
                                   (pplev(ig,l) - pplev(ig,l+1)) / g
                endif
-               if(water)then
-                  reffcol(ig,2) = reffcol(ig,2) + zq(ig,l,igcm_h2o_ice) * &
-                                  reffrad(ig,l,2) *    &
+               if(water.and.(iaero_h2o.ne.0)) then
+                  reffcol(ig,iaero_h2o) = reffcol(ig,iaero_h2o) + zq(ig,l,igcm_h2o_ice) * &
+                                  reffrad(ig,l,iaero_h2o) *    &
                                   (pplev(ig,l) - pplev(ig,l+1)) / g
                endif
@@ -1732,7 +1741,7 @@
                 call wstats(ngridmx,trim(noms(iq))//'_col',trim(noms(iq))//'_col',    & 
                           'kg m^-2',2,qcol(1,iq) )
-                call wstats(ngridmx,trim(noms(iq))//'_reff',                          & 
-                          trim(noms(iq))//'_reff',                                   & 
-                          'm',3,reffrad(1,1,iq))
+!                call wstats(ngridmx,trim(noms(iq))//'_reff',                          & 
+!                          trim(noms(iq))//'_reff',                                   & 
+!                          'm',3,reffrad(1,1,iq))
               end do
             if (water) then
@@ -1809,8 +1818,16 @@
 
 !     Output aerosols
-        if (igcm_co2_ice.ne.0) call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3,reffrad(1,1,1))
-        if (igcm_h2o_ice.ne.0) call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3,reffrad(1,1,2))
-        if (igcm_co2_ice.ne.0) call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol','um kg m^-2',2,reffcol(1,1))
-        if (igcm_h2o_ice.ne.0) call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol','um kg m^-2',2,reffcol(1,2))
+        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0)   &
+           call writediagfi(ngridmx,'CO2ice_reff','CO2ice_reff','m',3, &
+           reffrad(1,1,iaero_co2))
+        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0)   &
+           call writediagfi(ngridmx,'H2Oice_reff','H2Oice_reff','m',3, &
+           reffrad(1,1,iaero_h2o))
+        if (igcm_co2_ice.ne.0.and.iaero_co2.ne.0)   &
+           call writediagfi(ngridmx,'CO2ice_reffcol','CO2ice_reffcol', &
+           'um kg m^-2',2,reffcol(1,iaero_co2))
+        if (igcm_h2o_ice.ne.0.and.iaero_h2o.ne.0)   &
+           call writediagfi(ngridmx,'H2Oice_reffcol','H2Oice_reffcol', &
+           'um kg m^-2',2,reffcol(1,iaero_h2o))
 
 !     Output tracers
Index: /trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90	(revision 725)
+++ /trunk/LMDZ.GENERIC/libf/phystd/radcommon_h.F90	(revision 726)
@@ -46,7 +46,4 @@
 !   ~~~~~~~~~
 ! 
-! tauvis: dust optical depth at reference wavelength  ("longrefvis" set
-! in dimradmars.h : typically longrefvis = 0.67E-6 m, as measured by Viking )
-!
 ! For the "naerkind" kind of aerosol radiative properties : 
 ! QVISsQREF  :  Qext / Qext("longrefvis")
@@ -61,5 +58,4 @@
 ! omegaIR  :  mean single scattering albedo
 ! gIR      :  mean assymetry factor
-
 
       REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI)
@@ -78,5 +74,4 @@
       real*8 tgasmin, tgasmax
 
-      real tauvis
       real QVISsQREF(L_NSPECTV,naerkind,nsizemax)
       real omegavis(L_NSPECTV,naerkind,nsizemax)
Index: /trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90	(revision 725)
+++ /trunk/LMDZ.GENERIC/libf/phystd/radinc_h.F90	(revision 726)
@@ -5,5 +5,5 @@
 #include "dimensions.h"
 #include "bands.h"
-#include"scatterers.h"
+#include "scatterers.h"
 
 !======================================================================
Index: /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 725)
+++ /trunk/LMDZ.GENERIC/libf/phystd/rcm1d.F	(revision 726)
@@ -1,5 +1,4 @@
       program rcm1d
 
-      use radcommon_h, only: tauvis
 ! to use  'getin'
       use ioipsl_getincom 
@@ -469,12 +468,4 @@
       write(*,*)" specOLR = ",specOLR
 
-!!! AS12: what shall we do with this?
-c Proprietes des poussiere aerosol
-c --------------------------------
-      tauvis=0.2 ! default value for tauvis
-      print *,'Reference dust opacity at 700 Pa ?'
-      call getin("tauvis",tauvis)
-      write(*,*) " tauvis = ",tauvis
- 
 c
 c  pour le schema d'ondes de gravite
@@ -709,10 +700,8 @@
 !         write(103,*) 'Ls=',zls*180./pi
 !         write(103,*) 'Lat=', lati(1)*180./pi
-!         write(103,*) 'Tau=', tauvis/700*psurf
 !         write(103,*) 'RunEnd - Atmos. Temp. File'
 !         write(103,*) 'RunEnd - Atmos. Temp. File'
 !         write(104,*) 'Ls=',zls*180./pi
 !         write(104,*) 'Lat=', lati(1)
-!         write(104,*) 'Tau=', tauvis/700*psurf
 !         write(104,*) 'RunEnd - Atmos. Temp. File'
         ENDIF
Index: /trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90	(revision 725)
+++ /trunk/LMDZ.GENERIC/libf/phystd/suaer_corrk.F90	(revision 726)
@@ -2,5 +2,5 @@
 
       ! inputs
-      use radinc_h,    only: L_NSPECTI,L_NSPECTV,naerkind,nsizemax,iim,jjm
+      use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,iim,jjm,naerkind
       use radcommon_h, only: blamv,blami,lamrefir,lamrefvis
       use datafile_mod, only: datadir
@@ -10,4 +10,5 @@
       use radcommon_h, only: radiustab,nsize,tstellar
       use radcommon_h, only: qrefvis,qrefir,omegarefvis,omegarefir
+      use aerosol_mod
 
       implicit none
@@ -44,5 +45,4 @@
 
 #include "callkeys.h"
-
 !     Optical properties (read in external ASCII files)
       INTEGER          :: nwvl  ! Number of wavelengths in
@@ -101,55 +101,68 @@
 !     Please also choose the reference wavelengths of each aerosol
 
-      forgetit=.true.
-
-!     naerkind=1
+      do iaer=1,naerkind
+       if (iaer.eq.iaero_co2)then
+        forgetit=.true.
+        print*, 'naerkind= co2', iaer
+
 !     visible
-      if(forgetit)then
-         file_id(1,1) = 'optprop_co2_vis_ff.dat' ! Francois' values
-      else
-         file_id(1,1) = 'optprop_co2ice_vis_n50.dat'
-      endif
-      lamrefvis(1)=1.5E-6   ! 1.5um OMEGA/MEx ???
+        if(forgetit)then
+           file_id(iaer,1) = 'optprop_co2_vis_ff.dat' ! Francois' values
+        else
+           file_id(iaer,1) = 'optprop_co2ice_vis_n50.dat'
+        endif
+        lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx ???
 
 !     infrared
-      if(forgetit)then
-         file_id(1,2) = 'optprop_co2_ir_ff.dat' ! Francois' values
-      else
-         file_id(1,2) = 'optprop_co2ice_ir_n50.dat'
-      endif
-      lamrefir(1)=12.1E-6   ! 825cm-1 TES/MGS ???
-         
+        if(forgetit)then
+           file_id(iaer,2) = 'optprop_co2_ir_ff.dat' ! Francois' values
+        else
+           file_id(iaer,2) = 'optprop_co2ice_ir_n50.dat'
+        endif
+        lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS ???
+       endif ! CO2 aerosols  
 ! NOTE: these lamref's are currently for dust, not CO2 ice.
 ! JB thinks this shouldn't matter too much, but it needs to be 
 ! fixed / decided for the final version.
 
-      IF (naerkind .gt. 1) THEN
-!     naerkind=2
+       if (iaer.eq.iaero_h2o) then
+        print*, 'naerkind= h2o', iaer
+
 !     visible
-         file_id(2,1) = 'optprop_icevis_n50.dat'
-         lamrefvis(2)=1.5E-6   ! 1.5um OMEGA/MEx
+         file_id(iaer,1) = 'optprop_icevis_n50.dat'
+         lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx
 !     infrared
-         file_id(2,2) = 'optprop_iceir_n50.dat'
-         lamrefir(2)=12.1E-6   ! 825cm-1 TES/MGS
-      ENDIF
-
-      IF (naerkind .eq. 3) THEN
-!     naerkind=3
+         file_id(iaer,2) = 'optprop_iceir_n50.dat'
+         lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
+       endif
+
+       if (iaer.eq.iaero_dust) then
+        print*, 'naerkind= dust', iaer
+
 !     visible
-         file_id(naerkind,1) = 'optprop_dustvis_n50.dat'
-         !lamrefvis(3)=1.5E-6   ! 1.5um OMEGA/MEx
-         lamrefvis(naerkind)=0.67e-6
+         file_id(iaer,1) = 'optprop_dustvis_n50.dat'
+         !lamrefvis(iaer)=1.5E-6   ! 1.5um OMEGA/MEx
+         lamrefvis(iaer)=0.67e-6
 !     infrared
-         file_id(naerkind,2) = 'optprop_dustir_n50.dat'
-         !lamrefir(3)=12.1E-6   ! 825cm-1 TES/MGS
-         lamrefir(naerkind)=9.3E-6
-
-      ENDIF
-
-      IF (naerkind .gt. 3) THEN
-         print*,'naerkind = ',naerkind
-         print*,'but we only have data for 3 types, exiting.'
-         call abort
-      ENDIF
+         file_id(iaer,2) = 'optprop_dustir_n50.dat'
+         !lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
+         lamrefir(iaer)=9.3E-6
+       endif 
+
+       if (iaer.eq.iaero_h2so4) then
+         print*, 'naerkind= h2so4', iaer
+
+!     visible
+         file_id(iaer,1) = 'optprop_h2so4vis_n20.dat'
+         !lamrefir(iaer)=  doesn't exist?
+         lamrefvis(iaer)=1.5E-6   ! no idea, must find
+!     infrared
+         file_id(iaer,2) = 'optprop_h2so4ir_n20.dat'
+         !lamrefir(iaer)=12.1E-6   ! 825cm-1 TES/MGS
+         lamrefir(iaer)=9.3E-6 ! no idea, must find
+! added by LK
+       endif
+       
+      enddo
 
 !------------------------------------------------------------------
@@ -282,5 +295,5 @@
 
       jfile = jfile+1
-      IF ((idomain.NE.1).OR.(iaer.NE.1)) THEN
+      IF ((idomain.NE.iaero_co2).OR.(iaer.NE.iaero_co2)) THEN
          endwhile = .true.
       ENDIF
@@ -300,7 +313,7 @@
 
 !     1.5 If Francois' data, convert wvl to metres
-      if(forgetit)then
+       if(iaer.eq.iaero_co2)then
+         if(forgetit)then
          !   print*,'please sort out forgetit for naerkind>1'
-         if(iaer.eq.1)then
             do iwvl=1,nwvl
                wvl(iwvl)=wvl(iwvl)*1.e-6
Index: /trunk/LMDZ.GENERIC/libf/phystd/tabfi.F
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/tabfi.F	(revision 725)
+++ /trunk/LMDZ.GENERIC/libf/phystd/tabfi.F	(revision 726)
@@ -45,5 +45,4 @@
 ! to use  'getin'
       use ioipsl_getincom 
-      use radcommon_h, only: tauvis
 
       implicit none
@@ -138,5 +137,4 @@
       emisice(2) = tab_cntrl(tab0+25)
       emissiv    = tab_cntrl(tab0+26)
-      tauvis     = tab_cntrl(tab0+27)  ! dust opt. depth vis.
       iceradius(1)= tab_cntrl(tab0+31) ! mean scat radius of CO2 snow (north)
       iceradius(2)= tab_cntrl(tab0+32) ! mean scat radius of CO2 snow (south)
@@ -200,6 +198,4 @@
       write(*,5) '(33)   dtemisice(1)',tab_cntrl(tab0+33),dtemisice(1)
       write(*,5) '(34)   dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2)
-
-      write(*,5) '(27)         tauvis',tab_cntrl(tab0+27),tauvis
 
       write(*,5) '(35)        volcapa',tab_cntrl(tab0+35),volcapa
@@ -228,8 +224,6 @@
       write(*,*) '(31 et 32) iceradius : mean scat radius of CO2 snow'
       write(*,*) '(33 et 34) dtemisice : time scale for snow',
-     &           '                       metamorphism'
-      write(*,*) '(27)        tauvis : mean dust vis. reference ',
-     &           '                     opacity'
-      write(*,*) '(35)     volcapa  : soil volumetric heat capacity'
+     &           'metamorphism'
+      write(*,*) '(35)      volcapa : soil volumetric heat capacity'
       write(*,*) '(18)     obliquit : planet obliquity (deg)'
       write(*,*) '(17)     peri_day : periastron date (sols since Ls=0)'
@@ -365,12 +359,4 @@
           write(*,*) ' dtemisice(2) (new value):',dtemisice(2)
 
-        else if (modif(1:len_trim(modif)) .eq. 'tauvis') then
-          write(*,*) 'current value:',tauvis
-          write(*,*) 'enter new value:'
- 114      read(*,*,iostat=ierr) tauvis
-          if(ierr.ne.0) goto 114
-          write(*,*) 
-          write(*,*) ' tauvis (new value):',tauvis
-
         else if (modif(1:len_trim(modif)) .eq. 'obliquit') then
           write(*,*) 'current value:',obliquit
@@ -535,6 +521,4 @@
       write(*,5) '(34)   dtemisice(2)',tab_cntrl(tab0+34),dtemisice(2)
  
-      write(*,5) '(27)         tauvis',tab_cntrl(tab0+27),tauvis
-
       write(*,5) '(35)        volcapa',tab_cntrl(tab0+35),volcapa
 
Index: /trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90
===================================================================
--- /trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90	(revision 725)
+++ /trunk/LMDZ.GENERIC/libf/phystd/turbdiff.F90	(revision 726)
@@ -137,9 +137,14 @@
 
          if(water)then
-                ivap=igcm_h2o_vap
-                iliq=igcm_h2o_ice
-		iliq_surf=igcm_h2o_vap
-		iice_surf=igcm_h2o_ice ! simply to make the code legible               
-                                  ! to be generalised later
+             ivap=igcm_h2o_vap
+             iliq=igcm_h2o_ice
+             iliq_surf=igcm_h2o_vap
+	     iice_surf=igcm_h2o_ice ! simply to make the code legible               
+                                  ! to be generalised
+	 else
+             ivap=0
+             iliq=0
+             iliq_surf=0
+	     iice_surf=0 ! simply to make the code legible               	  
          endif
          sensibFlux(:)=0.
