Index: trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F	(revision 2642)
+++ trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F	(revision 2643)
@@ -10,10 +10,10 @@
 
       SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
-     &    pq,pt,tauscaling,dust_rad_adjust,tau_pref_scenario,
-     &    tau_pref_gcm,tau,taucloudtes,aerosol,dsodust,reffrad,
-     &    QREFvis3d,QREFir3d,omegaREFir3d,
-     &    totstormfract,clearatm,dsords,dsotop,
-     &    nohmons,
-     &    clearsky,totcloudfrac)
+     &   pq,pt,tauscaling,dust_rad_adjust,IRtoVIScoef,tau_pref_scenario,
+     &   tau_pref_gcm,tau,taucloudtes,aerosol,dsodust,reffrad,
+     &   QREFvis3d,QREFir3d,omegaREFir3d,
+     &   totstormfract,clearatm,dsords,dsotop,
+     &   nohmons,
+     &   clearsky,totcloudfrac)
                                                          
       use ioipsl_getin_p_mod, only: getin_p
@@ -35,8 +35,10 @@
      &            iaer_stormdust_doubleq,
      &            iaer_topdust_doubleq
-      use dust_param_mod, only: odpref, freedust
+      use dust_param_mod, only: odpref, freedust,
+     &                          reff_driven_IRtoVIS_scenario 
       use dust_scaling_mod, only: compute_dustscaling
       use density_co2_ice_mod, only: density_co2_ice
       use surfdat_h,only: alpha_hmons,contains_mons
+      use read_dust_scenario_mod, only: read_dust_scenario
       
        IMPLICIT NONE
@@ -89,5 +91,5 @@
                           !   (for direct comparison with TES)
       REAL, INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! optical
-                           ! depth of each aerosl in each layer
+                           ! depth of each aerosol in each layer
       REAL, INTENT(OUT) ::  dsodust(ngrid,nlayer) ! density scaled opacity
                             ! of (background) dust
@@ -110,5 +112,8 @@
       REAL,INTENT(OUT) :: tauscaling(ngrid) ! Scaling factor for qdust and Ndust
       REAL,INTENT(INOUT) :: dust_rad_adjust(ngrid) ! Radiative adjustment 
-                          ! factor for dust
+                          ! factor for dust                                  
+      REAL,INTENT(INOUT) :: IRtoVIScoef(ngrid) ! conversion coefficient to apply on
+                                               ! scenario absorption IR (9.3um) CDOD
+                                               ! = tau_pref_gcm_VIS / tau_pref_gcm_IR
       REAL,INTENT(IN) :: totcloudfrac(ngrid) ! total water ice cloud fraction
       LOGICAL,INTENT(IN) :: clearsky ! true to compute RT without water ice clouds
@@ -134,4 +139,8 @@
                                !   reference wavelength
       REAL topdust0(ngrid)
+      
+      REAL aerosol_IRabs(ngrid,nlayer)
+      REAL taudust_IRabs(ngrid)
+      REAL taudust_VISext(ngrid)
 
 !      -- CO2 clouds
@@ -295,17 +304,21 @@
           tau_pref_scenario(ig) = tau_pref_scenario(1)
         end do
+!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! NB: here, IRtoVIScoef=2.6 
+! ( useful to be here only if iddist=0 (Pollack90 vertical distribution) )      
       ELSE IF ((iaervar.ge.6).and.(iaervar.le.8)) THEN
       ! clim, cold or warm synthetic scenarios
         call read_dust_scenario(ngrid,nlayer,zday,pplev,
-     &                          tau_pref_scenario)
+     &                          IRtoVIScoef,tau_pref_scenario)
       ELSE IF ((iaervar.ge.24).and.(iaervar.le.35))
      &     THEN  ! << MY... dust scenarios >>
         call read_dust_scenario(ngrid,nlayer,zday,pplev,
-     &                          tau_pref_scenario)
+     &                          IRtoVIScoef,tau_pref_scenario)
       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,
-     &                          tau_pref_scenario)
+     &                          IRtoVIScoef,tau_pref_scenario)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!     
       ELSE
         call abort_physic("aeropacity","wrong value for iaervar",1)
@@ -385,5 +398,9 @@
 c        (transport of mass and number mixing ratio)
 c==================================================================
-
+          ! Some initialisations for the IRtoVIScoef
+          aerosol_IRabs(:,:)=0.
+          taudust_IRabs(:)=0.
+          taudust_VISext(:)=0.
+      
           DO l=1,nlayer
             IF (l.LE.cstdustlevel) THEN
@@ -405,4 +422,16 @@
      &          ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
      &          pq(ig,cstdustlevel,igcm_dust_mass)
+                
+                if (reff_driven_IRtoVIS_scenario) then
+                 if ((clearatm).and.(nohmons)) then ! the IRtoVIScoef is computed only during the first call to the RT
+                 ! OPTICAL DEPTH in IR absorption to compute the IRtoVIScoef
+                  aerosol_IRabs(ig,l) =
+     &            (  0.75 * QREFir3d(ig,cstdustlevel,iaer) /
+     &            ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
+     &            ( 1. - omegaREFir3d(ig,cstdustlevel,iaer) ) *
+     &            pq(ig,cstdustlevel,igcm_dust_mass) *
+     &            ( pplev(ig,l) - pplev(ig,l+1) ) / g
+                 endif
+                endif
               ENDDO
             ELSE
@@ -422,7 +451,31 @@
      &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
      &          pq(ig,l,igcm_dust_mass)
+                
+                if (reff_driven_IRtoVIS_scenario) then
+                 if ((clearatm).and.(nohmons)) then ! the IRtoVIScoef is computed only during the first call to the RT
+                 ! OPTICAL DEPTH in IR absorption to compute the IRtoVIScoef
+                  aerosol_IRabs(ig,l) =
+     &            (  0.75 * QREFir3d(ig,l,iaer) /
+     &            ( rho_dust * reffrad(ig,l,iaer) )  ) *
+     &            ( 1. - omegaREFir3d(ig,l,iaer) ) *
+     &            pq(ig,l,igcm_dust_mass) *
+     &            ( pplev(ig,l) - pplev(ig,l+1) ) / g
+                 endif
+                endif
               ENDDO
             ENDIF
-          ENDDO
+            if (reff_driven_IRtoVIS_scenario) then
+             if ((clearatm).and.(nohmons)) then ! the IRtoVIScoef is computed only during the first call to the RT
+              taudust_VISext(:) = taudust_VISext(:) + aerosol(:,l,iaer)
+              taudust_IRabs(:) = taudust_IRabs(:) + aerosol_IRabs(:,l)
+             endif
+            endif
+          ENDDO
+          
+          if (reff_driven_IRtoVIS_scenario) then
+           if ((clearatm).and.(nohmons)) then ! the IRtoVIScoef is computed only during the first call to the RT
+            IRtoVIScoef(:) = taudust_VISext(:) / taudust_IRabs(:)
+           endif
+          endif
 
 c==================================================================
@@ -647,4 +700,22 @@
 ! 3. Specific treatments for the dust aerosols
 
+! here IRtoVIScoef has been updated, we can call again read_dust_scenario
+      if (reff_driven_IRtoVIS_scenario) then
+        IF ((iaervar.ge.6).and.(iaervar.le.8)) THEN
+        ! clim, cold or warm synthetic scenarios
+          call read_dust_scenario(ngrid,nlayer,zday,pplev,
+     &                            IRtoVIScoef,tau_pref_scenario)
+        ELSE IF ((iaervar.ge.24).and.(iaervar.le.35))
+     &       THEN  ! << MY... dust scenarios >>
+          call read_dust_scenario(ngrid,nlayer,zday,pplev,
+     &                            IRtoVIScoef,tau_pref_scenario)
+        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,
+     &                            IRtoVIScoef,tau_pref_scenario)
+        ENDIF
+      endif
+      
 #ifdef DUSTSTORM
 c -----------------------------------------------------------------
@@ -790,6 +861,6 @@
 !      coefficients and adjust aerosol() dust opacities accordingly
       call compute_dustscaling(ngrid,nlayer,naerkind,naerdust,zday,pplev
-     &                         ,tau_pref_scenario,tauscaling,
-     &                          dust_rad_adjust,aerosol)
+     &                         ,tau_pref_scenario,IRtoVIScoef,
+     &                          tauscaling,dust_rad_adjust,aerosol)
 
 ! 3.2. Recompute tau_pref_gcm, the reference dust opacity, based on dust tracer
Index: trunk/LMDZ.MARS/libf/phymars/callradite_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callradite_mod.F	(revision 2642)
+++ trunk/LMDZ.MARS/libf/phymars/callradite_mod.F	(revision 2643)
@@ -9,5 +9,5 @@
      $     dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,
      $     fluxtop_sw,tau_pref_scenario,tau_pref_gcm,
-     &     tau,aerosol,dsodust,tauscaling,dust_rad_adjust,
+     &     tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
      $     taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice,
      $     rstormdust,rtopdust,totstormfract,clearatm,dsords,dsotop,
@@ -179,4 +179,7 @@
       REAL,INTENT(INOUT) :: dust_rad_adjust(ngrid) ! Radiative adjustment 
                           ! factor for dust
+      REAL,INTENT(INOUT) :: IRtoVIScoef(ngrid) ! conversion coefficient to apply on
+                                               ! scenario absorption IR (9.3um) CDOD
+                                               ! = tau_pref_gcm_VIS / tau_pref_gcm_IR
       REAL,INTENT(IN) :: albedo(ngrid,2),emis(ngrid)
       REAL,INTENT(IN) :: ls,zday
@@ -443,6 +446,7 @@
 c     Computing aerosol optical depth in each layer:
       CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
-     &    pq,pt,tauscaling,dust_rad_adjust,tau_pref_scenario,
-     &    tau_pref_gcm,tau,taucloudtes,aerosol,dsodust,reffrad,
+     &    pq,pt,tauscaling,dust_rad_adjust,IRtoVIScoef,
+     &    tau_pref_scenario,tau_pref_gcm,tau,taucloudtes,
+     &    aerosol,dsodust,reffrad,
      &    QREFvis3d,QREFir3d,omegaREFir3d,
      &    totstormfract,clearatm,dsords,dsotop,
Index: trunk/LMDZ.MARS/libf/phymars/compute_dtau_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/compute_dtau_mod.F90	(revision 2642)
+++ trunk/LMDZ.MARS/libf/phymars/compute_dtau_mod.F90	(revision 2643)
@@ -15,5 +15,6 @@
         SUBROUTINE compute_dtau(ngrid,nlayer,                           &
                                  zday,pplev,tau_pref_gcm,               &
-                                 ptimestep,dustliftday,local_time)
+                                 ptimestep,local_time,IRtoVIScoef,      &
+                                 dustliftday)
 
         USE geometry_mod, only: longitude_deg
@@ -23,4 +24,5 @@
         USE dimradmars_mod, only: tauvis
         USE dust_param_mod, only: odpref, t_scenario_sol
+        USE read_dust_scenario_mod, only: read_dust_scenario
         
         IMPLICIT NONE
@@ -36,4 +38,7 @@
         REAL, INTENT(in) :: ptimestep 
         REAL, INTENT(in) :: local_time(ngrid)
+        REAL, INTENT(in) :: IRtoVIScoef(ngrid) ! conversion coefficient to apply on
+                                               ! scenario absorption IR (9.3um) CDOD
+                                               ! = tau_pref_gcm_VIS / tau_pref_gcm_IR
         REAL, INTENT(out) :: dustliftday(ngrid) ! Dust injection rate (s-1)
         
@@ -71,4 +76,5 @@
           zday_scenario=zday_scenario+1      ! opacity of the dust scenario is read the day after
           call read_dust_scenario(ngrid,nlayer,zday_scenario,pplev,     &
+                                         IRtoVIScoef,                   &
                                          tau_pref_target)
         endif
Index: trunk/LMDZ.MARS/libf/phymars/conf_phys.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 2642)
+++ trunk/LMDZ.MARS/libf/phymars/conf_phys.F	(revision 2643)
@@ -49,5 +49,6 @@
       use dust_param_mod, only: dustbin, doubleq, submicron, active,
      &                          lifting, freedust, callddevil,
-     &                          dustscaling_mode
+     &                          dustscaling_mode,
+     &                          reff_driven_IRtoVIS_scenario
       use aeropacity_mod, only: iddist, topdustref
       USE mod_phys_lmdz_transfert_para, ONLY: bcast
@@ -516,5 +517,5 @@
          ! OR AT LEAST NO TO FAR, TO AVOID FALLACIOUS INTERPOLATIONS.
          !
-         dustiropacity="tes" !default value - is expected to shift to mcs one day
+         dustiropacity="tes" !default value
          call getin_p("dustiropacity",dustiropacity)
          write(*,*)" dustiropacity = ",trim(dustiropacity)
@@ -530,4 +531,25 @@
      &          "invalid dustiropacity option value",1)
          end select
+! Dust scenario IR to VIS conversion
+         write(*,*)"Use an IR to VIS conversion coefficient"
+         write(*,*)"for the dust scenario, that is dependent"
+         write(*,*)"on the GCM dust effective radius,"
+         write(*,*)"instead of a fixed 2.6 coefficient ?"
+         reff_driven_IRtoVIS_scenario=.false. !default value
+         call getin_p("reff_driven_IRtoVIS_scenario",
+     &                 reff_driven_IRtoVIS_scenario)
+         write(*,*)" reff_driven_IRtoVIS_scenario = ",
+     &               reff_driven_IRtoVIS_scenario
+! Test of incompatibility:
+! if reff_driven_IRtoVIS_scenario=.true.,
+! dustrefir must be 9.3E-6 = scenarios' wavelength
+         if (reff_driven_IRtoVIS_scenario .and.
+     &      (dustrefir.ne.9.3E-6)) then
+           print*,'if reff_driven_IRtoVIS_scenario is used, then '//
+     &           'the GCM IR reference wavelength should be the one '//
+     &           'of the scenarios (dustiropacity=tes)'
+           call abort_physic(modname,
+     &        "reff_driven_IRtoVIS_scenario requires tes wavelength",1)
+         endif
 
 ! callddevil
@@ -539,5 +561,4 @@
 ! Test of incompatibility:
 ! if dustdevil is used, then dustbin should be > 0
-
          if (callddevil.and.(dustbin.lt.1)) then
            print*,'if dustdevil is used, then dustbin should > 0'
@@ -545,4 +566,5 @@
      &          "callddevil requires dustbin > 0",1)
          endif
+         
 ! sedimentation
          write(*,*) "Gravitationnal sedimentation ?"
Index: trunk/LMDZ.MARS/libf/phymars/dust_param_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dust_param_mod.F90	(revision 2642)
+++ trunk/LMDZ.MARS/libf/phymars/dust_param_mod.F90	(revision 2643)
@@ -8,4 +8,5 @@
   LOGICAL,SAVE :: freedust ! if true: no rescaling (via tauscaling) of the dust mass and number
   LOGICAL,SAVE :: callddevil ! flag to activate dust devil (dust lifing/injection) parametrization
+  LOGICAL,SAVE :: reff_driven_IRtoVIS_scenario ! use GCM dust size to convert IR scenarios to VIS
 
 !$OMP THREADPRIVATE(active, doubleq,submicron,lifting,freedust,     &
@@ -28,5 +29,5 @@
                     ! tau_pref_scenario is deemed exact
 
-!$OMP THREADPRIVATE(tauscaling, dustscaling_mode,dust_rad_adjust)
+!$OMP THREADPRIVATE(tauscaling,dustscaling_mode,dust_rad_adjust)
 
 contains
Index: trunk/LMDZ.MARS/libf/phymars/dust_rad_adjust_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dust_rad_adjust_mod.F90	(revision 2642)
+++ trunk/LMDZ.MARS/libf/phymars/dust_rad_adjust_mod.F90	(revision 2643)
@@ -13,9 +13,11 @@
 
   subroutine compute_dust_rad_adjust(ngrid,nlayer,zday,pplev, &
-                                     taudust,dust_rad_adjust)
+                                     taudust,IRtoVIScoef,     &
+                                     dust_rad_adjust)
  
   use geometry_mod, only: longitude_deg
   use time_phylmdz_mod, only: dtphys, daysec
   use dust_param_mod, only: odpref, t_scenario_sol
+  use read_dust_scenario_mod, only: read_dust_scenario
   
   implicit none
@@ -26,4 +28,7 @@
   real,intent(in) :: pplev(ngrid,nlayer+1) ! pressure (Pa) at layer boundaries
   real,intent(in) :: taudust(ngrid) ! visible dust columns opacity in the GCM
+  real,intent(in) :: IRtoVIScoef(ngrid) ! conversion coefficient to apply on
+                                        ! scenario absorption IR (9.3um) CDOD
+                                        ! = tau_pref_gcm_VIS / tau_pref_gcm_IR
   real,intent(out) :: dust_rad_adjust(ngrid) ! radiative adjustment coefficient 
                       ! for dust
@@ -74,4 +79,5 @@
 !                                       tau_pref_scenario)
   call read_dust_scenario(ngrid,nlayer,zday_scenario_next,pplev, &
+                                       IRtoVIScoef,              &
                                        tau_pref_scenario_next)
 
Index: trunk/LMDZ.MARS/libf/phymars/dust_scaling_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dust_scaling_mod.F90	(revision 2642)
+++ trunk/LMDZ.MARS/libf/phymars/dust_scaling_mod.F90	(revision 2643)
@@ -7,5 +7,6 @@
   subroutine compute_dustscaling(ngrid,nlayer,naerkind,naerdust, &
                                  zday,pplev, &
-                                 tau_pref_scenario,tauscaling, &
+                                 tau_pref_scenario,IRtoVIScoef, &
+                                 tauscaling, &
                                  dust_rad_adjust,aerosol)
     
@@ -24,4 +25,7 @@
     real,intent(in) :: tau_pref_scenario(ngrid) ! prescribed visible dust 
                        ! opacity column at odpref reference pressure
+    real,intent(in) :: IRtoVIScoef(ngrid) ! conversion coefficient to apply on
+                                          ! scenario absorption IR (9.3um) CDOD
+                                          ! = tau_pref_gcm_VIS / tau_pref_gcm_IR
     real,intent(out) :: tauscaling(ngrid) ! dust scaling factor
     real,intent(inout) :: dust_rad_adjust(ngrid) ! Radiative adjustment 
@@ -74,6 +78,7 @@
       if (zday/=zday_prev_call) then
         call compute_dust_rad_adjust(ngrid,nlayer,zday,pplev, &
-                                     taudust,dust_rad_adjust)
-       endif ! of if (zday/=zday_prev_call)
+                                     taudust,IRtoVIScoef,     &
+                                     dust_rad_adjust)
+      endif ! of if (zday/=zday_prev_call)
                                      
       ! update zday_prev_call
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2642)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2643)
@@ -65,5 +65,5 @@
      &                          tauscaling, odpref, dustbin,
      &                          dustscaling_mode, dust_rad_adjust,
-     &                          freedust
+     &                          freedust, reff_driven_IRtoVIS_scenario
       use turb_mod, only: q2, wstar, ustar, sensibFlux, 
      &                    zmax_th, hfmax_th, turb_resolved
@@ -301,4 +301,7 @@
       REAL tau_pref_scenario(ngrid) ! prescribed dust column visible opacity
                                     ! at odpref
+      REAL IRtoVIScoef(ngrid) ! conversion coefficient to apply on
+                              ! scenario absorption IR (9.3um) CDOD
+                              ! = tau_pref_gcm_VIS / tau_pref_gcm_IR
       REAL tau_pref_gcm(ngrid) ! dust column visible opacity at odpref in the GCM
 c     rocket dust storm
@@ -761,4 +764,8 @@
       dsotop(:,:)=0.
       dwatercap(:)=0
+      
+      ! Dust scenario conversion coefficient from IRabs to VISext
+      IRtoVIScoef(1:ngrid)=2.6 ! initialized with former value from Montabone et al 2015
+                               ! recomputed in aeropacity if reff_driven_IRtoVIS_scenario=.true.
 
 #ifdef DUSTSTORM
@@ -967,5 +974,5 @@
      &     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,
      &     fluxtop_sw,tau_pref_scenario,tau_pref_gcm,
-     &     tau,aerosol,dsodust,tauscaling,dust_rad_adjust,
+     &     tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
      &     taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice,
      &     rstormdust,rtopdust,totstormfract,clearatm,dsords,dsotop,
@@ -983,6 +990,6 @@
      &              dist_sol,igout,zdtlwclf,zdtswclf,fluxsurf_lwclf,
      &              fluxsurf_swclf,fluxtop_lwclf,fluxtop_swclf,
-     &              tau_pref_scenario,tau_pref_gcm,
-     &              tau,aerosol,dsodust,tauscaling,dust_rad_adjust,
+     &              tau_pref_scenario,tau_pref_gcm,tau,aerosol,
+     &              dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef,
      &              taucloudtesclf,rdust,
      &              rice,nuice,riceco2, nuiceco2,co2ice,rstormdust,
@@ -1179,4 +1186,5 @@
      &                      tsurf,igout,totstormfract,
      &                      tauscaling,dust_rad_adjust,
+     &                      IRtoVIScoef,
 c               input sub-grid scale cloud
      &                      clearsky,totcloudfrac,
@@ -1240,5 +1248,5 @@
      &                zzlay,zdtsw,zdtlw,
      &                icount,zday,zls,tsurf,igout,aerosol, 
-     &                tauscaling,dust_rad_adjust,
+     &                tauscaling,dust_rad_adjust,IRtoVIScoef,
      &                totstormfract,clearatm,
      &                clearsky,totcloudfrac,
@@ -1272,5 +1280,6 @@
              CALL compute_dtau(ngrid,nlayer,
      &        			zday,pplev,tau_pref_gcm,
-     &        			ptimestep,dustliftday,local_time)
+     &        			ptimestep,local_time,IRtoVIScoef,
+     &                          dustliftday)
            endif ! end of if (dustinjection.gt.0)
 
@@ -3247,4 +3256,10 @@
      &                 'Visible dust optical depth at 610Pa in the GCM',
      &                 'NU',2,tau_pref_gcm)
+     
+      if (reff_driven_IRtoVIS_scenario) then
+             call WRITEDIAGFI(ngrid,'IRtoVIScoef',
+     &  'conversion coeff for dust tau from abs9.3um to ext0.67um',
+     &                        '/',2,IRtoVIScoef)
+      endif
 
       if (tracer.and.(dustbin.ne.0)) then
Index: trunk/LMDZ.MARS/libf/phymars/read_dust_scenario.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/read_dust_scenario.F90	(revision 2642)
+++ 	(revision )
@@ -1,375 +1,0 @@
-subroutine read_dust_scenario(ngrid,nlayer,zday,pplev,tau_pref_scenario)
-
-! Reading of the dust scenario file
-
-use netcdf
-use geometry_mod, only: latitude, longitude ! in radians
-use datafile_mod, only: datadir
-use dust_param_mod, only: odpref, dustscaling_mode
-use planete_h, only: year_day
-USE mod_phys_lmdz_transfert_para, ONLY: bcast
-implicit none
-
-include "callkeys.h"
-
-integer, intent(in) :: ngrid,nlayer
-real, intent(in) :: zday ! date in martian days
-real, dimension(ngrid,nlayer+1), intent(in) :: pplev
-real, dimension(ngrid), intent(out) :: tau_pref_scenario ! visible dust column
-                                       ! opacity at odpref from scenario
-
-! Local variables
-
-real :: realday
-integer nid,nvarid,ierr
-integer ltloop,lsloop,iloop,jloop,varloop,ig
-real, dimension(2) :: taubuf
-real tau1(4),tau
-real alt(4)
-integer latp(4),lonp(4)
-integer yinf,ysup,xinf,xsup,tinf,tsup
-real latinf,latsup,loninf,lonsup
-real latintmp,lonintmp,latdeg,londeg
-real colat,dlat,dlon,colattmp
-logical, save :: firstcall=.true.
-logical :: timeflag
-real,save :: radeg,pi
-integer :: timedim,londim,latdim
-real, dimension(:), allocatable, save :: lat,lon,time
-real, dimension(:,:,:), allocatable, save :: tautes
-integer, save :: timelen,lonlen,latlen
-character(len=33),save :: filename
-
-!$OMP THREADPRIVATE(firstcall,radeg,pi, &
-!$OMP                lat,lon,time,tautes,  &
-!$OMP            timelen,lonlen,latlen,filename)
-
-realday=mod(zday,669.)
-
-! AS: firstcall OK absolute
-if (firstcall) then
-   firstcall=.false.
-
-   pi=acos(-1.)
-   radeg=180/pi
-   
-   ! assimilated dust file: (NB: iaervar is a common in "callkeys.h")
-   ! iaervar=4 means read dust_tes.nc file
-   ! iaervar=6 means read dust_cold.nc file
-   ! iaervar=7 means read dust_warm.nc file
-   ! iaervar=8 means read dust_clim.nc file
-   ! iaervar=24 means read dust_MY24.nc file
-   ! iaervar=25 means read dust_MY25.nc file
-   ! iaervar=26 means read dust_MY26.nc file, etc.
-   if (iaervar.eq.4) then
-   ! NB: 4: old TES assimilated MY24 dust scenarios (at 700Pa ref pressure!)
-     filename="dust_tes.nc"
-   else if (iaervar.eq.6) then
-     filename="dust_cold.nc"
-   else if (iaervar.eq.7) then
-     filename="dust_warm.nc"
-   else if (iaervar.eq.8) then
-     filename="dust_clim.nc"
-   else if ((iaervar.ge.24).and.(iaervar.le.35)) then
-     write(filename,'(a7,i2,a3)')"dust_MY",iaervar,".nc"
-   ! 124,125,126: old TES assimilated dust scenarios (at 700Pa ref pressure!)
-   else if  (iaervar.eq.124) then
-     filename="dust_tes_MY24.nc" 
-   else if (iaervar.eq.125) then
-     filename="dust_tes_MY25.nc" 
-   else if (iaervar.eq.126) then
-     filename="dust_tes_MY26.nc"
-   endif
-
-
-!$OMP MASTER
-   
-   ierr=nf90_open(trim(datadir)//"/"//trim(filename),NF90_NOWRITE,nid)
-   IF (ierr.NE.nf90_noerr) THEN
-     write(*,*)'Problem opening ',trim(filename),' (in phymars/read_dust_scenario.F90)'
-     write(*,*)'It should be in :',trim(datadir),'/'
-     write(*,*)'1) You can change this directory address in callfis.def with'
-     write(*,*)'   datadir=/path/to/datafiles'
-     write(*,*)'2) If necessary, dust*.nc (and other datafiles)'
-     write(*,*)'   can be obtained online on:'
-     write(*,*)'   http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
-     CALL abort_physic("read_dust_scenario","Cannot find file",1)
-   ENDIF
-
-   ierr=nf90_inq_dimid(nid,"Time",timedim)
-   if (ierr.ne.nf90_noerr) then
-    ! 'Time' dimension not found, look for 'time'
-    ierr=nf90_inq_dimid(nid,"time",timedim)
-    if (ierr.ne.nf90_noerr) then
-      write(*,*)"Error: read_dust_scenario <time> not found"
-    endif
-   endif
-   ierr=nf90_inquire_dimension(nid,timedim,len=timelen)
-   ierr=nf90_inq_dimid(nid,"latitude",latdim)
-   ierr=nf90_inquire_dimension(nid,latdim,len=latlen)
-   ierr=nf90_inq_dimid(nid,"longitude",londim)
-   ierr=nf90_inquire_dimension(nid,londim,len=lonlen)
-
-!$OMP END MASTER
-
-   call bcast(lonlen)
-   call bcast(latlen)
-   call bcast(timelen)
-
-   if (.not.allocated(tautes)) allocate(tautes(lonlen,latlen,timelen))
-   if (.not.allocated(lat)) allocate(lat(latlen), lon(lonlen), time(timelen))
-
-!$OMP MASTER
-
-   ! "dustop" if loading visible extinction opacity
-   ! "cdod" if loading IR absorption opacity
-   ierr=nf90_inq_varid(nid,"dustop",nvarid)
-   if (ierr.eq.nf90_noerr) then
-     ierr=nf90_get_var(nid,nvarid,tautes)
-     IF (ierr .NE. nf90_noerr) THEN
-        PRINT*, "Error: read_dust_scenario <dustop> not found"
-        write(*,*)trim(nf90_strerror(ierr))
-        call abort_physic("read_dust_scenario","dustop not found",1)
-     ENDIF
-   else
-     ! did not find "dustop" , look for "cdod"
-     ierr=nf90_inq_varid(nid,"cdod",nvarid)
-     ierr=nf90_get_var(nid,nvarid,tautes)
-     IF (ierr .NE. nf90_noerr) THEN
-        PRINT*, "Error: read_dust_scenario <cdod> not found"
-        write(*,*)trim(nf90_strerror(ierr))
-        call abort_physic("read_dust_scenario","cdod not found",1)
-     ENDIF
-     ! and multiply by 2*1.3=2.6 to convert from IR absorption
-     ! to visible extinction opacity
-     tautes(:,:,:)=2.6*tautes(:,:,:)
-   endif
-
-   ierr=nf90_inq_varid(nid,"Time",nvarid)
-   ierr=nf90_get_var(nid,nvarid,time)
-   IF (ierr .NE. nf90_noerr) THEN
-      PRINT*, "Error: read_dust_scenario <Time> not found"
-      write(*,*)trim(nf90_strerror(ierr))
-      call abort_physic("read_dust_scenario","Time not found",1)
-   ENDIF
-
-   ! Adapt time axis values if using dustinjection scheme and/or
-   ! dustscaling_mode == 2
-   IF ((dustinjection>0).or.(dustscaling_mode==2)) THEN
-     ! Sanity check: we expect to have one map per sol
-     if (timelen/=year_day) then
-       write(*,*)"read_dust_scenario error: timelen=",timelen
-       write(*,*)" whereas it was exepceted to be year_day=",year_day
-       call abort_physic("read_dust_scenario","timelen/=year_day",1)
-     endif
-     ! fill time() with intergers from 0 to timelen-1
-     write(*,*)"read_dust_scenario: adapting time axis to integer values"
-     do iloop=1,timelen
-       time(iloop)=iloop-1
-     enddo
-   ENDIF ! of IF ((dustinjection>0).or.(dustscaling_mode==2))
-
-   ierr=nf90_inq_varid(nid,"latitude",nvarid)
-   ierr=nf90_get_var(nid,nvarid,lat)
-   IF (ierr .NE. nf90_noerr) THEN
-      PRINT*, "Error: read_dust_scenario <latitude> not found"
-      write(*,*)trim(nf90_strerror(ierr))
-      call abort_physic("read_dust_scenario","latitude not found",1)
-   ENDIF
-
-   ierr=nf90_inq_varid(nid,"longitude",nvarid)
-   ierr=nf90_get_var(nid,nvarid,lon)
-   IF (ierr .NE. nf90_noerr) THEN
-      PRINT*, "Error: read_dust_scenario <longitude> not found"
-      write(*,*)trim(nf90_strerror(ierr))
-      call abort_physic("read_dust_scenario","longitude not found",1)
-   ENDIF
-
-   ierr=nf90_close(nid)
-
-!$OMP END MASTER
-
-   call bcast(tautes)
-   call bcast(time)
-   call bcast(lat)
-   call bcast(lon)
-
-endif ! of if (firstcall)
-
-do ig=1,ngrid
-
-! Find the four nearest points, arranged as follows:
-!                               1 2
-!                               3 4
-
-   latdeg=latitude(ig)*radeg  ! latitude, in degrees
-   londeg=longitude(ig)*radeg  ! longitude, in degrees east
-   colat=90-latdeg        ! colatitude, in degrees
-
-! Ehouarn: rounding effects and/or specific compiler issues
-!          sometime cause londeg to be sligthly below -180 ...
-  if (londeg.lt.-180) then
-    ! check if it is by a large amount
-    if ((londeg+180).lt.-1.e-3) then
-      write(*,*) ' ig=',ig,' londeg=',londeg
-      call abort_physic("read_dust_scenario","longitude way out of range",1)
-    else
-      londeg=-180
-    endif
-  endif
-
- ! Find encompassing latitudes
-   if (colat<(90-lat(1))) then ! between north pole and lat(1)
-      ysup=1
-      yinf=1
-   else if (colat>=90-(lat(latlen))) then ! between south pole and lat(laten)
-      ysup=latlen
-      yinf=latlen
-   else ! general case
-      do iloop=2,latlen
-         if(colat<(90-lat(iloop))) then
-           ysup=iloop-1
-           yinf=iloop
-           exit
-         endif
-      enddo
-   endif
-
-   latinf=lat(yinf)
-   latsup=lat(ysup)
-
-
- ! Find encompassing longitudes
-   ! Note: in input file, lon(1)=-180.
-   if (londeg>lon(lonlen)) then
-      xsup=1
-      xinf=lonlen
-      loninf=lon(xsup)
-      lonsup=180.0 ! since lon(1)=-180.0
-   else
-      do iloop=2,lonlen
-         if(londeg<lon(iloop)) then
-              xsup=iloop
-              xinf=iloop-1
-              exit
-         endif
-      enddo
-      loninf=lon(xinf)
-      lonsup=lon(xsup)
-   endif
-
-   if ((xsup.gt.lonlen).OR.(yinf.gt.latlen).OR.(xinf.lt.1)&
-     .OR.(ysup.lt.1)) then
-      write (*,*) "read_dust_scenario: SYSTEM ERROR on x or y in ts_gcm"
-      write (*,*) "xinf: ",xinf
-      write (*,*) "xsup: ",xsup
-      write (*,*) "yinf: ",yinf
-      write (*,*) "ysup: ",ysup
-      call abort_physic("read_dust_scenario","invalid coordinates",1)
-   endif
-
-!   loninf=lon(xinf)
-!   lonsup=lon(xsup)
-!   latinf=lat(yinf)
-!   latsup=lat(ysup)
-
-! The four neighbouring points are arranged as follows:
-!                               1 2
-!                               3 4
-
-   latp(1)=ysup
-   latp(2)=ysup
-   latp(3)=yinf
-   latp(4)=yinf
-
-   lonp(1)=xinf
-   lonp(2)=xsup
-   lonp(3)=xinf
-   lonp(4)=xsup
-
-! Linear interpolation on time, for all four neighbouring points
-
-   if ((realday<time(1)).or.(realday>time(timelen))) then
-      tinf=timelen
-      tsup=1
-      timeflag=.true.
-   else 
-      timeflag=.false.
-      do iloop=2,timelen
-         if (realday<=time(iloop)) then
-            tinf=iloop-1
-            tsup=iloop
-            exit
-         endif
-      enddo
-   endif
-
-! Bilinear interpolation on the four nearest points
-
-   if ((colat<(90-lat(1))).OR.(colat>(90-lat(latlen))).OR.(latsup==latinf)) then
-      dlat=0
-   else
-      dlat=((90-latinf)-colat)/((90-latinf)-(90-latsup))
-   endif
-
-   if (lonsup==loninf) then
-      dlon=0
-   else
-      dlon=(londeg-loninf)/(lonsup-loninf)
-   endif
-
-   do iloop=1,4
-      taubuf(1)=tautes(lonp(iloop),latp(iloop),tinf)
-      taubuf(2)=tautes(lonp(iloop),latp(iloop),tsup)
-      if (timeflag) then
-         if (realday>time(timelen)) then
-            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)+(669-time(tinf))) 
-         else
-            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*realday/(time(tsup)+(669-time(tinf)))
-         endif
-      else
-         tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)-time(tinf))
-      endif
-      if (tau1(iloop)<0) then
-          write (*,*) "read_dust_scenario: SYSTEM ERROR on tau"
-          write (*,*) "utime ",realday
-          write (*,*) "time(tinf) ",time(tinf)
-          write (*,*) "time(tsup) ",time(tsup)
-          write (*,*) "tau1 ",taubuf(1)
-          write (*,*) "tau2 ",taubuf(2)
-          write (*,*) "tau ",tau1(iloop)
-          call abort_physic("read_dust_scenario","negative tau!",1)
-      endif
-   enddo
-
-   if ((dlat>1).OR.(dlon>1) .OR. (dlat<0) .OR. (dlon<0)) then
-      write (*,*) "read_dust_scenario: SYSTEM ERROR on dlat or dlon in ts_gcm"
-      write (*,*) "dlat: ",dlat
-      write (*,*) "lat: ",latdeg
-      write (*,*) "dlon: ",dlon
-      write (*,*) "lon: ",londeg
-      call abort_physic("read_dust_scenario","negative dlat or dlon!",1)
-   endif
-
-   tau= dlat*(dlon*(tau1(2)+tau1(3)-tau1(1)-tau1(4))+tau1(1)-tau1(3)) +dlon*(tau1(4)-tau1(3))+tau1(3)
-
-   tau_pref_scenario(ig)=tau
-! 
-enddo ! of ig=1,ngrid
-
-if (filename(1:8)=="dust_tes") then
-  ! when using old TES files, correction for:
-  ! - the reference pressure was 700Pa (unlike 610Pa now)
-  ! - the 1.3 conversion factor from IR absorbtion opacity to
-  !   IR extinction opacity
-  tau_pref_scenario(:)=tau_pref_scenario(:)*1.3*(odpref/700.)
-endif
-
-if (swrtype.eq.1) then ! Fouquart (NB: swrtype is set in callkeys.h)
- ! when using old radiative transfer (like in MCD 4.x)
- ! needed to decrease opacity (*0.825) to compensate overestimation of
- ! heating rates
-  tau_pref_scenario(:)=tau_pref_scenario(:)*0.825/1.3
-endif
-
-end
Index: trunk/LMDZ.MARS/libf/phymars/read_dust_scenario_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/read_dust_scenario_mod.F90	(revision 2643)
+++ trunk/LMDZ.MARS/libf/phymars/read_dust_scenario_mod.F90	(revision 2643)
@@ -0,0 +1,397 @@
+module read_dust_scenario_mod
+
+implicit none
+
+contains
+
+subroutine read_dust_scenario(ngrid,nlayer,zday,pplev,IRtoVIScoef,tau_pref_scenario)
+
+! Reading of the dust scenario file
+
+use netcdf
+use geometry_mod, only: latitude, longitude ! in radians
+use datafile_mod, only: datadir
+use dust_param_mod, only: odpref, dustscaling_mode
+use planete_h, only: year_day
+USE mod_phys_lmdz_transfert_para, ONLY: bcast
+implicit none
+
+include "callkeys.h"
+
+integer, intent(in) :: ngrid,nlayer
+real, intent(in) :: zday ! date in martian days
+real, dimension(ngrid,nlayer+1), intent(in) :: pplev
+real, dimension(ngrid), intent(in) :: IRtoVIScoef ! conversion coefficient to apply on
+                                       ! scenario absorption IR (9.3um) CDOD
+                                       ! = tau_pref_gcm_VIS / tau_pref_gcm_IR
+real, dimension(ngrid), intent(out) :: tau_pref_scenario ! visible dust column
+                                       ! opacity at odpref from scenario
+
+! Local variables
+
+real :: realday
+integer nid,nvarid,ierr
+integer ltloop,lsloop,iloop,jloop,varloop,ig
+real, dimension(2) :: taubuf
+real tau1(4),tau
+real alt(4)
+integer latp(4),lonp(4)
+integer yinf,ysup,xinf,xsup,tinf,tsup
+real latinf,latsup,loninf,lonsup
+real latintmp,lonintmp,latdeg,londeg
+real colat,dlat,dlon,colattmp
+logical, save :: firstcall=.true.
+logical :: timeflag
+real,save :: radeg,pi
+integer :: timedim,londim,latdim
+real, dimension(:), allocatable, save :: lat,lon,time
+real, dimension(:,:,:), allocatable, save :: tautes
+integer, save :: timelen,lonlen,latlen
+character(len=33),save :: filename
+logical,save :: IRscenario ! =true if the scenarios are IR absorption opacity (variable 'cdod')
+                           ! and not VIS extinction (variable 'dustop')
+
+!$OMP THREADPRIVATE(firstcall,radeg,pi, &
+!$OMP               lat,lon,time,tautes,  &
+!$OMP               timelen,lonlen,latlen,  &
+!$OMP               filename,IRscenario)
+
+realday=mod(zday,669.)
+
+! AS: firstcall OK absolute
+if (firstcall) then
+   firstcall=.false.
+
+   pi=acos(-1.)
+   radeg=180/pi
+   
+   ! assimilated dust file: (NB: iaervar is a common in "callkeys.h")
+   ! iaervar=4 means read dust_tes.nc file
+   ! iaervar=6 means read dust_cold.nc file
+   ! iaervar=7 means read dust_warm.nc file
+   ! iaervar=8 means read dust_clim.nc file
+   ! iaervar=24 means read dust_MY24.nc file
+   ! iaervar=25 means read dust_MY25.nc file
+   ! iaervar=26 means read dust_MY26.nc file, etc.
+   if (iaervar.eq.4) then
+   ! NB: 4: old TES assimilated MY24 dust scenarios (at 700Pa ref pressure!)
+     filename="dust_tes.nc"
+   else if (iaervar.eq.6) then
+     filename="dust_cold.nc"
+   else if (iaervar.eq.7) then
+     filename="dust_warm.nc"
+   else if (iaervar.eq.8) then
+     filename="dust_clim.nc"
+   else if ((iaervar.ge.24).and.(iaervar.le.35)) then
+     write(filename,'(a7,i2,a3)')"dust_MY",iaervar,".nc"
+   ! 124,125,126: old TES assimilated dust scenarios (at 700Pa ref pressure!)
+   else if  (iaervar.eq.124) then
+     filename="dust_tes_MY24.nc" 
+   else if (iaervar.eq.125) then
+     filename="dust_tes_MY25.nc" 
+   else if (iaervar.eq.126) then
+     filename="dust_tes_MY26.nc"
+   endif
+
+
+!$OMP MASTER
+   
+   ierr=nf90_open(trim(datadir)//"/"//trim(filename),NF90_NOWRITE,nid)
+   IF (ierr.NE.nf90_noerr) THEN
+     write(*,*)'Problem opening ',trim(filename),' (in phymars/read_dust_scenario.F90)'
+     write(*,*)'It should be in :',trim(datadir),'/'
+     write(*,*)'1) You can change this directory address in callfis.def with'
+     write(*,*)'   datadir=/path/to/datafiles'
+     write(*,*)'2) If necessary, dust*.nc (and other datafiles)'
+     write(*,*)'   can be obtained online on:'
+     write(*,*)'   http://www.lmd.jussieu.fr/~lmdz/planets/mars/datadir'
+     CALL abort_physic("read_dust_scenario","Cannot find file",1)
+   ENDIF
+
+   ierr=nf90_inq_dimid(nid,"Time",timedim)
+   if (ierr.ne.nf90_noerr) then
+    ! 'Time' dimension not found, look for 'time'
+    ierr=nf90_inq_dimid(nid,"time",timedim)
+    if (ierr.ne.nf90_noerr) then
+      write(*,*)"Error: read_dust_scenario <time> not found"
+    endif
+   endif
+   ierr=nf90_inquire_dimension(nid,timedim,len=timelen)
+   ierr=nf90_inq_dimid(nid,"latitude",latdim)
+   ierr=nf90_inquire_dimension(nid,latdim,len=latlen)
+   ierr=nf90_inq_dimid(nid,"longitude",londim)
+   ierr=nf90_inquire_dimension(nid,londim,len=lonlen)
+
+!$OMP END MASTER
+
+   call bcast(lonlen)
+   call bcast(latlen)
+   call bcast(timelen)
+
+   if (.not.allocated(tautes)) allocate(tautes(lonlen,latlen,timelen))
+   if (.not.allocated(lat)) allocate(lat(latlen), lon(lonlen), time(timelen))
+
+!$OMP MASTER
+
+   ! "dustop" if loading visible extinction opacity
+   ! "cdod" if loading IR absorption opacity
+   ierr=nf90_inq_varid(nid,"dustop",nvarid)
+   if (ierr.eq.nf90_noerr) then
+     IRscenario = .false.
+     ierr=nf90_get_var(nid,nvarid,tautes)
+     IF (ierr .NE. nf90_noerr) THEN
+        PRINT*, "Error: read_dust_scenario <dustop> not found"
+        write(*,*)trim(nf90_strerror(ierr))
+        call abort_physic("read_dust_scenario","dustop not found",1)
+     ENDIF
+   else
+     ! did not find "dustop" , look for "cdod"
+     IRscenario = .true.
+     ierr=nf90_inq_varid(nid,"cdod",nvarid)
+     ierr=nf90_get_var(nid,nvarid,tautes)
+     IF (ierr .NE. nf90_noerr) THEN
+        PRINT*, "Error: read_dust_scenario <cdod> not found"
+        write(*,*)trim(nf90_strerror(ierr))
+        call abort_physic("read_dust_scenario","cdod not found",1)
+     ENDIF
+   endif
+
+   ierr=nf90_inq_varid(nid,"Time",nvarid)
+   ierr=nf90_get_var(nid,nvarid,time)
+   IF (ierr .NE. nf90_noerr) THEN
+      PRINT*, "Error: read_dust_scenario <Time> not found"
+      write(*,*)trim(nf90_strerror(ierr))
+      call abort_physic("read_dust_scenario","Time not found",1)
+   ENDIF
+
+   ! Adapt time axis values if using dustinjection scheme and/or
+   ! dustscaling_mode == 2
+   IF ((dustinjection>0).or.(dustscaling_mode==2)) THEN
+     ! Sanity check: we expect to have one map per sol
+     if (timelen/=year_day) then
+       write(*,*)"read_dust_scenario error: timelen=",timelen
+       write(*,*)" whereas it was exepceted to be year_day=",year_day
+       call abort_physic("read_dust_scenario","timelen/=year_day",1)
+     endif
+     ! fill time() with intergers from 0 to timelen-1
+     write(*,*)"read_dust_scenario: adapting time axis to integer values"
+     do iloop=1,timelen
+       time(iloop)=iloop-1
+     enddo
+   ENDIF ! of IF ((dustinjection>0).or.(dustscaling_mode==2))
+
+   ierr=nf90_inq_varid(nid,"latitude",nvarid)
+   ierr=nf90_get_var(nid,nvarid,lat)
+   IF (ierr .NE. nf90_noerr) THEN
+      PRINT*, "Error: read_dust_scenario <latitude> not found"
+      write(*,*)trim(nf90_strerror(ierr))
+      call abort_physic("read_dust_scenario","latitude not found",1)
+   ENDIF
+
+   ierr=nf90_inq_varid(nid,"longitude",nvarid)
+   ierr=nf90_get_var(nid,nvarid,lon)
+   IF (ierr .NE. nf90_noerr) THEN
+      PRINT*, "Error: read_dust_scenario <longitude> not found"
+      write(*,*)trim(nf90_strerror(ierr))
+      call abort_physic("read_dust_scenario","longitude not found",1)
+   ENDIF
+
+   ierr=nf90_close(nid)
+
+!$OMP END MASTER
+
+   call bcast(tautes)
+   call bcast(time)
+   call bcast(lat)
+   call bcast(lon)
+   call bcast(IRscenario)
+
+endif ! of if (firstcall)
+
+
+do ig=1,ngrid
+
+! Find the four nearest points, arranged as follows:
+!                               1 2
+!                               3 4
+
+   latdeg=latitude(ig)*radeg  ! latitude, in degrees
+   londeg=longitude(ig)*radeg  ! longitude, in degrees east
+   colat=90-latdeg        ! colatitude, in degrees
+
+! Ehouarn: rounding effects and/or specific compiler issues
+!          sometime cause londeg to be sligthly below -180 ...
+  if (londeg.lt.-180) then
+    ! check if it is by a large amount
+    if ((londeg+180).lt.-1.e-3) then
+      write(*,*) ' ig=',ig,' londeg=',londeg
+      call abort_physic("read_dust_scenario","longitude way out of range",1)
+    else
+      londeg=-180
+    endif
+  endif
+
+ ! Find encompassing latitudes
+   if (colat<(90-lat(1))) then ! between north pole and lat(1)
+      ysup=1
+      yinf=1
+   else if (colat>=90-(lat(latlen))) then ! between south pole and lat(laten)
+      ysup=latlen
+      yinf=latlen
+   else ! general case
+      do iloop=2,latlen
+         if(colat<(90-lat(iloop))) then
+           ysup=iloop-1
+           yinf=iloop
+           exit
+         endif
+      enddo
+   endif
+
+   latinf=lat(yinf)
+   latsup=lat(ysup)
+
+
+ ! Find encompassing longitudes
+   ! Note: in input file, lon(1)=-180.
+   if (londeg>lon(lonlen)) then
+      xsup=1
+      xinf=lonlen
+      loninf=lon(xsup)
+      lonsup=180.0 ! since lon(1)=-180.0
+   else
+      do iloop=2,lonlen
+         if(londeg<lon(iloop)) then
+              xsup=iloop
+              xinf=iloop-1
+              exit
+         endif
+      enddo
+      loninf=lon(xinf)
+      lonsup=lon(xsup)
+   endif
+
+   if ((xsup.gt.lonlen).OR.(yinf.gt.latlen).OR.(xinf.lt.1)&
+     .OR.(ysup.lt.1)) then
+      write (*,*) "read_dust_scenario: SYSTEM ERROR on x or y in ts_gcm"
+      write (*,*) "xinf: ",xinf
+      write (*,*) "xsup: ",xsup
+      write (*,*) "yinf: ",yinf
+      write (*,*) "ysup: ",ysup
+      call abort_physic("read_dust_scenario","invalid coordinates",1)
+   endif
+
+!   loninf=lon(xinf)
+!   lonsup=lon(xsup)
+!   latinf=lat(yinf)
+!   latsup=lat(ysup)
+
+! The four neighbouring points are arranged as follows:
+!                               1 2
+!                               3 4
+
+   latp(1)=ysup
+   latp(2)=ysup
+   latp(3)=yinf
+   latp(4)=yinf
+
+   lonp(1)=xinf
+   lonp(2)=xsup
+   lonp(3)=xinf
+   lonp(4)=xsup
+
+! Linear interpolation on time, for all four neighbouring points
+
+   if ((realday<time(1)).or.(realday>time(timelen))) then
+      tinf=timelen
+      tsup=1
+      timeflag=.true.
+   else 
+      timeflag=.false.
+      do iloop=2,timelen
+         if (realday<=time(iloop)) then
+            tinf=iloop-1
+            tsup=iloop
+            exit
+         endif
+      enddo
+   endif
+
+! Bilinear interpolation on the four nearest points
+
+   if ((colat<(90-lat(1))).OR.(colat>(90-lat(latlen))).OR.(latsup==latinf)) then
+      dlat=0
+   else
+      dlat=((90-latinf)-colat)/((90-latinf)-(90-latsup))
+   endif
+
+   if (lonsup==loninf) then
+      dlon=0
+   else
+      dlon=(londeg-loninf)/(lonsup-loninf)
+   endif
+
+   do iloop=1,4
+      taubuf(1)=tautes(lonp(iloop),latp(iloop),tinf)
+      taubuf(2)=tautes(lonp(iloop),latp(iloop),tsup)
+      if (timeflag) then
+         if (realday>time(timelen)) then
+            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)+(669-time(tinf))) 
+         else
+            tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*realday/(time(tsup)+(669-time(tinf)))
+         endif
+      else
+         tau1(iloop)=taubuf(1)+(taubuf(2)-taubuf(1))*(realday-time(tinf))/(time(tsup)-time(tinf))
+      endif
+      if (tau1(iloop)<0) then
+          write (*,*) "read_dust_scenario: SYSTEM ERROR on tau"
+          write (*,*) "utime ",realday
+          write (*,*) "time(tinf) ",time(tinf)
+          write (*,*) "time(tsup) ",time(tsup)
+          write (*,*) "tau1 ",taubuf(1)
+          write (*,*) "tau2 ",taubuf(2)
+          write (*,*) "tau ",tau1(iloop)
+          call abort_physic("read_dust_scenario","negative tau!",1)
+      endif
+   enddo
+
+   if ((dlat>1).OR.(dlon>1) .OR. (dlat<0) .OR. (dlon<0)) then
+      write (*,*) "read_dust_scenario: SYSTEM ERROR on dlat or dlon in ts_gcm"
+      write (*,*) "dlat: ",dlat
+      write (*,*) "lat: ",latdeg
+      write (*,*) "dlon: ",dlon
+      write (*,*) "lon: ",londeg
+      call abort_physic("read_dust_scenario","negative dlat or dlon!",1)
+   endif
+
+   tau= dlat*(dlon*(tau1(2)+tau1(3)-tau1(1)-tau1(4))+tau1(1)-tau1(3)) +dlon*(tau1(4)-tau1(3))+tau1(3)
+
+   if (IRscenario) then
+     ! if the scenarios are in infrared, multiply by IRtoVIScoef(ig)
+     ! to convert from IR absorption to visible extinction opacity
+     tau_pref_scenario(ig)=tau * IRtoVIScoef(ig)
+   else
+     tau_pref_scenario(ig)=tau
+   endif
+   
+! 
+enddo ! of ig=1,ngrid
+
+if (filename(1:8)=="dust_tes") then
+  ! when using old TES files, correction for:
+  ! - the reference pressure was 700Pa (unlike 610Pa now)
+  ! - the 1.3 conversion factor from IR absorbtion opacity to
+  !   IR extinction opacity
+  tau_pref_scenario(:)=tau_pref_scenario(:)*1.3*(odpref/700.)
+endif
+
+if (swrtype.eq.1) then ! Fouquart (NB: swrtype is set in callkeys.h)
+ ! when using old radiative transfer (like in MCD 4.x)
+ ! needed to decrease opacity (*0.825) to compensate overestimation of
+ ! heating rates
+  tau_pref_scenario(:)=tau_pref_scenario(:)*0.825/1.3
+endif
+
+end subroutine read_dust_scenario
+
+end module read_dust_scenario_mod
Index: trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90	(revision 2642)
+++ trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90	(revision 2643)
@@ -27,4 +27,5 @@
                                  tsurf,igout,totstormfract,            &
                                  tauscaling,dust_rad_adjust,           &
+                                 IRtoVIScoef,                          &
 !             input sub-grid scale cloud
                                  clearsky,totcloudfrac,                &
@@ -79,5 +80,7 @@
       REAL, INTENT(IN) :: totstormfract(ngrid)
       REAL, INTENT(INOUT) :: tauscaling(ngrid)
-      REAL,INTENT(INOUT) :: dust_rad_adjust(ngrid)
+      REAL,INTENT(INOUT) :: dust_rad_adjust(ngrid)                              
+      REAL,INTENT(INOUT) :: IRtoVIScoef(ngrid) ! NB: not modified by this call to callradite,
+                                               ! the OUT is just here because callradite needs it
       
 !     subgrid scale water ice clouds
@@ -254,5 +257,5 @@
                  zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1,     &
                  fluxtop_sw1,tau_pref_scenario,tau_pref_gcm, &
-                 tau,aerosol,dsodust,tauscaling,dust_rad_adjust,       &
+                 tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef, &
                  taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice,rstormdust,rtopdust, &
                  totstormfract,clearatm,dsords,dsotop,nohmons,&
Index: trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90	(revision 2642)
+++ trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90	(revision 2643)
@@ -22,4 +22,5 @@
                                  icount,zday,zls,tsurf,igout,aerosol,  &
                                  tauscaling,dust_rad_adjust,           &
+                                 IRtoVIScoef,                          &
 !             input sub-grid scale rocket dust storm
                                  totstormfract,clearatm,               &
@@ -75,4 +76,7 @@
       REAL, INTENT(INOUT) :: tauscaling(ngrid)
       REAL,INTENT(INOUT) :: dust_rad_adjust(ngrid)
+      REAL,INTENT(INOUT) :: IRtoVIScoef(ngrid) ! NB: not modified by this call to callradite,
+                                               ! the OUT is just here because callradite needs it
+                                               
 !     input sub-grid scale rocket dust storm
       LOGICAL, INTENT(IN) :: clearatm
@@ -284,5 +288,5 @@
                  zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1,  &
                  fluxtop_sw1,tau_pref_scenario,tau_pref_gcm, &
-                 tau,aerosol,dsodust,tauscaling,dust_rad_adjust,    &
+                 tau,aerosol,dsodust,tauscaling,dust_rad_adjust,IRtoVIScoef, &
                  taucloudtes,rdust,rice,nuice,riceco2,nuiceco2,co2ice,rstormdust,rtopdust, &
                  totstormfract,clearatm,dsords,dsotop,nohmons,&
