Index: trunk/LMDZ.MARS/README
===================================================================
--- trunk/LMDZ.MARS/README	(revision 2414)
+++ trunk/LMDZ.MARS/README	(revision 2415)
@@ -3159,2 +3159,7 @@
 Correction of a bad unit for qdusttotal0 and qdusttotal1 (diagnostic outputs of the stormdust scheme)
 + some corrections of comments about aerosol and tauref in aeropacity
+
+== 13/10/2020 == EM
+Code cleanup: remove "tauref" and replace it with two distinct variables,
+"tau_pref_gcm" and "tau_pref_scenario" which are repectively the visible
+dust opacity column (at 610 Pa) in the GCM and prescribed by a dust scenario.
Index: trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F	(revision 2414)
+++ trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F	(revision 2415)
@@ -10,5 +10,6 @@
 
       SUBROUTINE aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
-     &    pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
+     &    pq,tauscaling,tau_pref_scenario,tau_pref_gcm,
+     &    tau,taucloudtes,aerosol,dsodust,reffrad,
      &    QREFvis3d,QREFir3d,omegaREFir3d,
      &    totstormfract,clearatm,dsords,dsotop,
@@ -53,26 +54,4 @@
 c         - reference pressure is now set to 610Pa (not 700Pa)
 c   
-c   input:
-c   ----- 
-c   ngrid             Number of gridpoint of horizontal grid
-c   nlayer            Number of layer
-c   nq                Number of tracer
-c   zday                  Date (time since Ls=0, in martian days)
-c   ls                Solar longitude (Ls) , radian
-c   pplay,pplev       pressure (Pa) in the middle and boundary of each layer
-c   pq                Dust mixing ratio (used if tracer =T and active=T).
-c   reffrad(ngrid,nlayer,naerkind)  Aerosol effective radius
-c   QREFvis3d(ngrid,nlayer,naerkind) \ 3d extinction coefficients
-c   QREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths;
-c   omegaREFir3d(ngrid,nlayer,naerkind)  / at reference wavelengths;
-c
-c   output:
-c   -------
-c   tauref       Prescribed mean column optical depth at 610 Pa 
-c   tau          Column total visible dust optical depth at each point
-c   aerosol      aerosol(ig,l,1) is the dust optical
-c                depth in layer l, grid point ig
-c   taualldust   CW17 total opacity for all dust scatterer stormdust included
-c
 c=======================================================================
       include "callkeys.h"
@@ -85,30 +64,49 @@
 c    Input/Output
 c    ------------
-      INTEGER, INTENT(IN) ::  ngrid,nlayer,nq
-      REAL, INTENT(IN) :: ls,zday
-      REAL, INTENT(IN) ::  pplev(ngrid,nlayer+1),pplay(ngrid,nlayer)
-      REAL, INTENT(IN) ::  pq(ngrid,nlayer,nq)
-      REAL, INTENT(OUT) :: tauref(ngrid) ! prescribed or computed
-                                         ! column dust opacity
-      REAL, INTENT(OUT) :: tau(ngrid,naerkind)
-      REAL,INTENT(OUT) :: taucloudtes(ngrid)! Cloud opacity at infrared
-                               !   reference wavelength using
-                               !   Qabs instead of Qext
-                               !   (direct comparison with TES)
-      REAL, INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind)
-      REAL, INTENT(OUT) ::  dsodust(ngrid,nlayer)
+      INTEGER,INTENT(IN) :: ngrid ! number of atmospheric columns
+      INTEGER,INTENT(IN) :: nlayer ! number of atmospheric layers
+      INTEGER,INTENT(IN) :: nq ! number of tracers
+      REAL,INTENT(IN) :: ls ! Solar Longitude (rad)
+      REAL,INTENT(IN) :: zday ! date (in martian sols) since Ls=0
+      REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! pressure (Pa) in the middle of
+                         ! each atmospheric layer
+      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! pressure (Pa) at the boundaries
+                         ! of the atmospheric layers
+      REAL,INTENT(IN) ::  pq(ngrid,nlayer,nq) ! tracers
+      REAL,INTENT(OUT) :: tau_pref_scenario(ngrid) ! prescribed dust column
+                          ! visible opacity at odpref from scenario
+      REAL,INTENT(OUT) :: tau_pref_gcm(ngrid) ! computed dust column
+                          ! visible opacity at odpref in the GCM
+      REAL,INTENT(OUT) :: tau(ngrid,naerkind) ! column total visible
+                          ! optical depth of each aerosol
+      REAL,INTENT(OUT) :: taucloudtes(ngrid)! Water ice cloud opacity at 
+                          !   infrared reference wavelength using
+                          !   Qabs instead of Qext
+                          !   (for direct comparison with TES)
+      REAL, INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! optical
+                           ! depth of each aerosl in each layer
+      REAL, INTENT(OUT) ::  dsodust(ngrid,nlayer) ! density scaled opacity
+                            ! of (background) dust
       REAL, INTENT(OUT) ::  dsords(ngrid,nlayer) !dso of stormdust
       REAL, INTENT(OUT) ::  dsotop(ngrid,nlayer) !dso of topdust  
-      REAL, INTENT(INOUT) :: reffrad(ngrid,nlayer,naerkind)
-      REAL, INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind)
-      REAL, INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind)
-      REAL, INTENT(IN) ::  omegaREFir3d(ngrid,nlayer,naerkind)
-      LOGICAL, INTENT(IN) :: clearatm
-      REAL, INTENT(IN) :: totstormfract(ngrid)
-      LOGICAL, INTENT(IN) :: nohmons
+      REAL, INTENT(INOUT) :: reffrad(ngrid,nlayer,naerkind) ! effective radius
+                             ! of the aerosols in the grid boxes
+      REAL, INTENT(IN) :: QREFvis3d(ngrid,nlayer,naerkind) ! 3D extinction
+                          ! coefficients (in the visible) of aerosols
+      REAL, INTENT(IN) :: QREFir3d(ngrid,nlayer,naerkind) ! 3D extinction
+                          ! coefficients (in the infra-red) of aerosols
+      REAL, INTENT(IN) :: omegaREFir3d(ngrid,nlayer,naerkind) ! at the
+                          ! reference wavelengths
+      LOGICAL, INTENT(IN) :: clearatm ! true to compute RT without stormdust
+                             ! and false to compute RT in rocket dust storms
+      REAL, INTENT(IN) :: totstormfract(ngrid) ! mesh fraction with a rocket
+                          ! dust storm
+      LOGICAL, INTENT(IN) :: nohmons ! true to compute RT without slope wind
+                             ! topdust, false to compute RT in the topdust
       REAL, INTENT(IN) :: alpha_hmons(ngrid)
       REAL, INTENT(OUT) ::  tauscaling(ngrid) ! Scaling factor for qdust and Ndust
-      REAL,INTENT(IN) :: totcloudfrac(ngrid) ! total cloud fraction
-      LOGICAL,INTENT(IN) :: clearsky ! true for part without clouds,false for part with clouds (total or sub-grid clouds)
+      REAL,INTENT(IN) :: totcloudfrac(ngrid) ! total water ice cloud fraction
+      LOGICAL,INTENT(IN) :: clearsky ! true to compute RT without water ice clouds
+      ! false to compute RT with clouds (total or sub-grid clouds)
 c
 c    Local variables :
@@ -237,19 +235,17 @@
       END IF ! end of if firstcall
 
-! 1. Get prescribed tauref, Dust column optical depth at "odpref" Pa 
-!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-      IF(freedust) THEN
-         tauref(:) = 0. ! tauref is computed after, instead of being forced
-
-      ELSE IF(iaervar.eq.1) THEN 
+! 1. Get prescribed tau_pref_scenario, Dust column optical depth at "odpref" Pa 
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+      IF(iaervar.eq.1) THEN 
          do ig=1, ngrid
-          tauref(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
+          tau_pref_scenario(ig)=max(tauvis,1.e-9) ! tauvis=cste (set in callphys.def
                                        ! or read in starfi
         end do
       ELSE IF (iaervar.eq.2) THEN   ! << "Viking" Scenario>>
 
-        tauref(1) = 0.7+.3*cos(ls+80.*pi/180.) ! like seen by VL1
+        tau_pref_scenario(1) = 0.7+.3*cos(ls+80.*pi/180.) ! like seen by VL1
         do ig=2,ngrid
-          tauref(ig) = tauref(1)
+          tau_pref_scenario(ig) = tau_pref_scenario(1)
         end do
 
@@ -262,30 +258,33 @@
           if (latitude(ig).ge.0) then
           ! Northern hemisphere
-            tauref(ig)= tauN +
+            tau_pref_scenario(ig)= tauN +
      &      (taueq-tauN)*0.5*(1+tanh((45-latitude(ig)*180./pi)*6/60))
           else
           ! Southern hemisphere
-            tauref(ig)= tauS +
+            tau_pref_scenario(ig)= tauS +
      &      (taueq-tauS)*0.5*(1+tanh((45+latitude(ig)*180./pi)*6/60))
           endif
         enddo ! of do ig=1,ngrid
       ELSE IF (iaervar.eq.5) THEN   ! << Escalier Scenario>>
-        tauref(1) = 2.5
+        tau_pref_scenario(1) = 2.5
         if ((ls.ge.30.*pi/180.).and.(ls.le.150.*pi/180.))
-     &                              tauref(1) = .2
+     &                              tau_pref_scenario(1) = .2
 
         do ig=2,ngrid
-          tauref(ig) = tauref(1)
+          tau_pref_scenario(ig) = tau_pref_scenario(1)
         end do
       ELSE IF ((iaervar.ge.6).and.(iaervar.le.8)) THEN
       ! clim, cold or warm synthetic scenarios
-        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
+        call read_dust_scenario(ngrid,nlayer,zday,pplev,
+     &                          tau_pref_scenario)
       ELSE IF ((iaervar.ge.24).and.(iaervar.le.34))
      &     THEN  ! << MY... dust scenarios >>
-        call read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
+        call read_dust_scenario(ngrid,nlayer,zday,pplev,
+     &                          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,tauref)
+        call read_dust_scenario(ngrid,nlayer,zday,pplev,
+     &                          tau_pref_scenario)
       ELSE
         call abort_physic("aeropacity","wrong value for iaervar",1)
@@ -354,5 +353,5 @@
              DO ig=1,ngrid
                 zp=odpref/pplay(ig,l)
-                aerosol(ig,l,1)= tauref(ig)/odpref *
+                aerosol(ig,l,1)= tau_pref_scenario(ig)/odpref *
      s           (pplev(ig,l)-pplev(ig,l+1))
      s           *max( exp(.03*(1.-max(zp,1.))) , 1.E-3 )
@@ -589,10 +588,10 @@
           DO l=1,nlayer
             DO ig=1,ngrid
-              tauref(ig) = tauref(ig) +
+              tau_pref_gcm(ig) = tau_pref_gcm(ig) +
      &                    aerosol(ig,l,iaerdust(iaer))
             ENDDO
           ENDDO
         ENDDO
-        tauref(:) = tauref(:) * odpref / pplev(:,1)
+        tau_pref_gcm(:) = tau_pref_gcm(:) * odpref / pplev(:,1)
 
 c--------------------------------------------------
@@ -656,5 +655,5 @@
       !! -- the additional reference opacity will
       !!      thus be taulocref*odpref/pplev
-      tauuser(ig)=max( tauref(ig) * pplev(ig,1) /odpref , 
+      tauuser(ig)=max( tau_pref_gcm(ig) * pplev(ig,1) /odpref , 
      &          taulocref * yeah )
 
@@ -696,9 +695,9 @@
           !! Mass mixing ratio perturbation due to local dust storm in each layer
           more_dust(ig,l,1)=
-     &                     (tauuser(ig)-(tauref(ig)
+     &                     (tauuser(ig)-(tau_pref_gcm(ig)
      &                      * pplev(ig,1) /odpref)) /
      &                      int_factor(ig)
           more_dust(ig,l,2)=
-     &                     (tauuser(ig)-(tauref(ig) *
+     &                     (tauuser(ig)-(tau_pref_gcm(ig) *
      &                     pplev(ig,1) /odpref))
      &                      / int_factor(ig) *
@@ -716,5 +715,5 @@
       ENDDO
       ENDIF !! IF (localstorm)
-      tauref(:)=0.
+      tau_pref_gcm(:)=0.
       ENDIF !! IF (firstcall)
 #endif
@@ -723,12 +722,13 @@
 ! 3.1. Compute "tauscaling", the dust rescaling coefficient and adjust
 !      aerosol() dust opacities accordingly
-      call compute_tauscaling(ngrid,nlayer,naerkind,naerdust,
-     &                        zday,pplev,tauref,tauscaling,aerosol)
-
-! 3.2. Recompute tauref, the reference dust opacity, based on dust tracer
+      call compute_tauscaling(ngrid,nlayer,naerkind,naerdust,zday,pplev,
+     &                        tau_pref_scenario,tauscaling,aerosol)
+
+! 3.2. Recompute tau_pref_gcm, the reference dust opacity, based on dust tracer
 !      mixing ratios and their optical properties
 
       IF (freedust) THEN
-        ! tauref has been initialized to 0 before.
+        ! Initialisation :
+        tau_pref_gcm(:)=0
         DO iaer=1,naerdust
           DO l=1,nlayer
@@ -744,10 +744,10 @@
       ENDIF
 #endif
-c      MV19: tauref must ALWAYS contain the opacity of all dust tracers
-       ! GCM DUST OPTICAL DEPTH tauref is to be compared
-       ! with the observation CDOD tauref_scenario
+c      MV19: tau_pref_gcm must ALWAYS contain the opacity of all dust tracers
+       ! GCM DUST OPTICAL DEPTH tau_pref_gcm is to be compared
+       ! with the observation CDOD tau_pref_scenario
        ! => visible wavelength
        IF (name_iaer(iaerdust(iaer)).eq."dust_doubleq") THEN
-              tauref(ig) = tauref(ig) +
+              tau_pref_gcm(ig) = tau_pref_gcm(ig) +
      &  (  0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
      &  ( rho_dust * reffrad(ig,l,iaerdust(iaer)) )  ) *
@@ -755,5 +755,5 @@
      &  ( pplev(ig,l) - pplev(ig,l+1) ) / g
        ELSE IF (name_iaer(iaerdust(iaer)).eq."stormdust_doubleq") THEN
-              tauref(ig) = tauref(ig) +
+              tau_pref_gcm(ig) = tau_pref_gcm(ig) +
      &  (  0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
      &  ( rho_dust * reffrad(ig,l,iaerdust(iaer)) )  ) *
@@ -761,5 +761,5 @@
      &  ( pplev(ig,l) - pplev(ig,l+1) ) / g 
        ELSE IF (name_iaer(iaerdust(iaer)).eq."topdust_doubleq") THEN
-              tauref(ig) = tauref(ig) +
+              tau_pref_gcm(ig) = tau_pref_gcm(ig) +
      &  (  0.75 * QREFvis3d(ig,l,iaerdust(iaer)) /
      &  ( rho_dust * reffrad(ig,l,iaerdust(iaer)) )  ) *
@@ -771,5 +771,8 @@
           ENDDO
         ENDDO
-        tauref(:) = tauref(:) * odpref / pplev(:,1)
+        tau_pref_gcm(:) = tau_pref_gcm(:) * odpref / pplev(:,1)
+      ELSE
+        ! dust opacity strictly follows what is imposed by the dust scenario
+        tau_pref_gcm(:)=tau_pref_scenario(:)
       ENDIF ! of IF (freedust)
 
Index: trunk/LMDZ.MARS/libf/phymars/callradite_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/callradite_mod.F	(revision 2414)
+++ trunk/LMDZ.MARS/libf/phymars/callradite_mod.F	(revision 2415)
@@ -8,5 +8,6 @@
      $     emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,
      $     dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,
-     $     fluxtop_sw,tauref,tau,aerosol,dsodust,tauscaling,
+     $     fluxtop_sw,tau_pref_scenario,tau_pref_gcm,
+     &     tau,aerosol,dsodust,tauscaling,
      $     taucloudtes,rdust,rice,nuice,co2ice,rstormdust,rtopdust,
      $     totstormfract,clearatm,dsords,dsotop,alpha_hmons,nohmons,
@@ -153,5 +154,4 @@
 c fluxtop_sw(ngrid,2)       outgoing upward flux SW for solar band#2 (W.m-2)
 
-c   tauref       Prescribed mean column optical depth at 610 Pa 
 c   tau          Column total visible dust optical depth at each point
 c   aerosol(ngrid,nlayer,naerkind)    aerosol extinction optical depth
@@ -187,6 +187,9 @@
       REAL,INTENT(OUT) :: fluxsurf_lw(ngrid), fluxtop_lw(ngrid)
       REAL,INTENT(OUT) :: fluxsurf_sw(ngrid,2), fluxtop_sw(ngrid,2)
-
-      REAL,INTENT(OUT) :: tauref(ngrid), tau(ngrid,naerkind)
+      REAL,INTENT(OUT) :: tau_pref_scenario(ngrid) ! prescribed dust column
+                          ! visible opacity at odpref from scenario
+      REAL,INTENT(OUT) :: tau_pref_gcm(ngrid) ! computed dust column
+                          ! visible opacity at odpref in the GCM
+      REAL,INTENT(OUT) :: tau(ngrid,naerkind)
       REAL,INTENT(OUT) :: taucloudtes(ngrid)! Cloud opacity at infrared
                                !   reference wavelength using
@@ -204,5 +207,5 @@
       REAL,INTENT(IN) :: totstormfract(ngrid) ! dust storm mesh fraction
       REAL,INTENT(OUT) :: rstormdust(ngrid,nlayer)  ! Storm dust geometric mean radius (m)
-      REAL,INTENT(INOUT) :: dsords(ngrid,nlayer) ! density scaled opacity for rocket dust storm dust
+      REAL,INTENT(OUT) :: dsords(ngrid,nlayer) ! density scaled opacity for rocket dust storm dust
       
 c     entrainment by slope wind
@@ -210,5 +213,5 @@
       REAL, INTENT(IN) :: alpha_hmons(ngrid) ! sub-grid scale topography mesh fraction
       REAL,INTENT(OUT) :: rtopdust(ngrid,nlayer)  ! Topdust geometric mean radius (m)
-      REAL,INTENT(INOUT) :: dsotop(ngrid,nlayer) ! density scaled opacity for topmons dust
+      REAL,INTENT(OUT) :: dsotop(ngrid,nlayer) ! density scaled opacity for topmons dust
       
 c     sub-grid scale water ice clouds
@@ -423,5 +426,6 @@
 c     Computing aerosol optical depth in each layer:
       CALL aeropacity(ngrid,nlayer,nq,zday,pplay,pplev,ls,
-     &    pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
+     &    pq,tauscaling,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 2414)
+++ trunk/LMDZ.MARS/libf/phymars/compute_dtau_mod.F90	(revision 2415)
@@ -5,5 +5,6 @@
         REAL,SAVE :: ti_injection_sol ! time of beginning injection
         REAL,SAVE :: tf_injection_sol ! time of end injection
-        REAL,PARAMETER :: t_scenario_sol=14/24.   ! time of tauref_scenario
+        REAL,PARAMETER :: t_scenario_sol=14/24. ! time of day at which
+                          ! tau_pref_scenario is deemed exact
 
         REAL,SAVE,ALLOCATABLE :: dtau(:) ! Dust opacity difference (at 610Pa)
@@ -13,5 +14,5 @@
 
         SUBROUTINE compute_dtau(ngrid,nlayer,                           &
-                                 zday,pplev,tauref,                     &
+                                 zday,pplev,tau_pref_gcm,               &
                                  ptimestep,dustliftday,local_time)
 
@@ -31,5 +32,6 @@
         REAL, INTENT(in) :: zday ! date at lon=0, in fraction of sols
         REAL, INTENT(in) :: pplev(ngrid,nlayer+1) ! pressure (Pa)
-        REAL, INTENT(in) :: tauref(ngrid) ! Computed dust opacity at 610Pa
+        REAL, INTENT(in) :: tau_pref_gcm(ngrid) ! Visible dust opacity column
+                            ! at 610Pa as computed in the GCM
         REAL, INTENT(in) :: ptimestep 
         REAL, INTENT(in) :: local_time(ngrid)
@@ -38,5 +40,6 @@
         INTEGER :: ig, l
         INTEGER, SAVE :: nb_daystep ! nomber of step a day
-        REAL :: tauref_scenario(ngrid) ! from dust scenario
+        REAL :: tau_pref_target(ngrid) ! dust opacity column at odpref=610 Pa
+                ! as extracted from dust scenario
         REAL :: zday_scenario
         REAL,ALLOCATABLE,SAVE :: local_time_prev(:)
@@ -59,16 +62,17 @@
         ENDIF
         
-        ! 1. Obtain tauref_scenario from dust scenario at zday+1
+        ! 1. Obtain tau_pref_target from dust scenario at zday+1
         if (iaervar.eq.1) then
-          tauref_scenario = tauvis
+          tau_pref_target = tauvis
         else
           zday_scenario=zday-modulo(zday,1.) ! integer value of the day: the scenario opacity is measured at 14:00
           zday_scenario=zday_scenario+1      ! opacity of the dust scenario is read the day after
           call read_dust_scenario(ngrid,nlayer,zday_scenario,pplev,     &
-                                         tauref_scenario)
+                                         tau_pref_target)
         endif
        ! for diagnostics
-        call WRITEDIAGFI(ngrid,"tauref_scenario","tauref_scenario",     &
-                          "",2,tauref_scenario)
+        call WRITEDIAGFI(ngrid,"tau_pref_target", &
+                          "target visible dust opacity column at 610Pa", &
+                          "",2,tau_pref_target)
 
         ! 2. Compute dtau() and dustliftday()
@@ -76,8 +80,8 @@
          IF ((local_time(ig).ge.t_scenario_sol).and.                    &
                  (local_time_prev(ig).lt.(t_scenario_sol)))THEN
-                 dtau(ig)=tauref_scenario(ig)-tauref(ig)
+                 dtau(ig)=tau_pref_target(ig)-tau_pref_gcm(ig)
          ENDIF
 
-        ! Use dtau (when positiove) to compute dustliftday
+        ! Use dtau (when positive) to compute dustliftday
          IF (dtau(ig).LT.0) THEN
              dustliftday(ig)=0.
Index: trunk/LMDZ.MARS/libf/phymars/dust_param_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dust_param_mod.F90	(revision 2414)
+++ trunk/LMDZ.MARS/libf/phymars/dust_param_mod.F90	(revision 2415)
@@ -12,5 +12,5 @@
 
   REAL,PARAMETER :: odpref = 610. ! Reference pressure (Pa) of
-                     ! DOD (Dust optical Depth) tauref
+                     ! DOD (Dust optical Depth) tau_pref_*
   
   REAL,SAVE,ALLOCATABLE :: tauscaling(:)   ! Convertion factor for qdust and Ndust
Index: trunk/LMDZ.MARS/libf/phymars/dust_scaling_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/dust_scaling_mod.F90	(revision 2414)
+++ trunk/LMDZ.MARS/libf/phymars/dust_scaling_mod.F90	(revision 2415)
@@ -7,5 +7,5 @@
   subroutine compute_tauscaling(ngrid,nlayer,naerkind,naerdust, &
                                 zday,pplev, &
-                                tauref,tauscaling,aerosol)
+                                tau_pref_scenario,tauscaling,aerosol)
     
     use dust_param_mod, only: tauscaling_mode, odpref
@@ -20,5 +20,6 @@
     real,intent(in) :: zday
     real,intent(in) :: pplev(ngrid,nlayer+1) ! inter-layer pressure (Pa)
-    real,intent(in) :: tauref(ngrid) ! prescribed dust opacity
+    real,intent(in) :: tau_pref_scenario(ngrid) ! prescribed visible dust 
+                       ! opacity column at odpref reference pressure
     real,intent(out) :: tauscaling(ngrid) ! dust scaling factor
     real,intent(inout) :: aerosol(ngrid,nlayer,naerkind) ! opacities
@@ -36,5 +37,5 @@
     if (tauscaling_mode == 1) then
       ! GCM v5.3 style: tauscaling is computed so that
-      ! aerosol() opacities correspond to the prescribed tauref()
+      ! aerosol() opacities correspond to the prescribed tau_pref_scenario()
       
       ! 1. compute dust column opacity using aerosol() dusts
@@ -49,5 +50,5 @@
       
       ! 2. compute the scaling factor
-      tauscaling(:)=tauref(:)*pplev(:,1)/odpref/taudust(:)
+      tauscaling(:)=tau_pref_scenario(:)*pplev(:,1)/odpref/taudust(:)
     endif ! of if (tauscaling_mode == 1)
     
Index: trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2414)
+++ trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2415)
@@ -281,5 +281,7 @@
       REAL fluxtop_lw(ngrid)       !Outgoing LW (IR) flux to space (W.m-2)
       REAL fluxtop_sw(ngrid,2)     !Outgoing SW (solar) flux to space (W.m-2)
-      REAL tauref(ngrid)           ! Reference column optical depth at odpref
+      REAL tau_pref_scenario(ngrid) ! prescribed dust column visible opacity
+                                    ! at odpref
+      REAL tau_pref_gcm(ngrid) ! dust column visible opacity at odpref in the GCM
 c     rocket dust storm
       REAL totstormfract(ngrid)     ! fraction of the mesh where the dust storm is contained
@@ -928,5 +930,6 @@
      &     emis,mu0,zplev,zplay,pt,tsurf,fract,dist_sol,igout,
      &     zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,
-     &     fluxtop_sw,tauref,tau,aerosol,dsodust,tauscaling,
+     &     fluxtop_sw,tau_pref_scenario,tau_pref_gcm,
+     &     tau,aerosol,dsodust,tauscaling,
      &     taucloudtes,rdust,rice,nuice,co2ice,rstormdust,rtopdust,
      &     totstormfract,clearatm,dsords,dsotop,alpha_hmons,nohmons,
@@ -943,5 +946,6 @@
      &              albedo,emis,mu0,zplev,zplay,pt,tsurf,fract,
      &              dist_sol,igout,zdtlwclf,zdtswclf,fluxsurf_lwclf,
-     &              fluxsurf_swclf,fluxtop_lwclf,fluxtop_swclf,tauref,
+     &              fluxsurf_swclf,fluxtop_lwclf,fluxtop_swclf,
+     &              tau_pref_scenario,tau_pref_gcm,
      &              tau,aerosol,dsodust,tauscaling,taucloudtesclf,rdust,
      &              rice,nuice,co2ice,rstormdust,rtopdust,totstormfract,
@@ -999,7 +1003,7 @@
      &                       longitude(igout)*180/pi
 
-           write(*,'(" tauref(",f4.0," Pa) =",f9.6,
+           write(*,'(" tau_pref_gcm(",f4.0," Pa) =",f9.6,
      &             " tau(",f4.0," Pa) =",f9.6)')
-     &            odpref,tauref(igout),
+     &            odpref,tau_pref_gcm(igout),
      &            odpref,tau(igout,1)*odpref/zplev(igout,1)
 c          ---------------------------------------------------------
@@ -1140,5 +1144,5 @@
 c               output
      &                      pdqrds,wspeed,dsodust,dsords,dsotop,
-     &                      tauref)
+     &                      tau_pref_scenario,tau_pref_gcm)
 
 c      update the tendencies of both dust after vertical transport
@@ -1204,5 +1208,5 @@
      &                nohmons,hsummit,
      &                pdqtop,wtop,dsodust,dsords,dsotop,
-     &                tauref)
+     &                tau_pref_scenario,tau_pref_gcm)
      
                    
@@ -1231,5 +1235,5 @@
 
              CALL compute_dtau(ngrid,nlayer,
-     &        			zday,pplev,tauref,
+     &        			zday,pplev,tau_pref_gcm,
      &        			ptimestep,dustliftday,local_time)
            endif ! end of if (dustinjection.gt.0)
@@ -2495,6 +2499,10 @@
         call wstats(ngrid,"watercap","H2O ice cover",
      &                "kg.m-2",2,watercap)
-        call wstats(ngrid,"tauref","reference dod at 610 Pa","NU",
-     &                2,tauref)
+        call wstats(ngrid,"tau_pref_scenario",
+     &              "prescribed visible dod at 610 Pa","NU",
+     &                2,tau_pref_scenario)
+        call wstats(ngrid,"tau_pref_gcm",
+     &              "visible dod at 610 Pa in the GCM","NU",
+     &                2,tau_pref_gcm)
         call wstats(ngrid,"fluxsurf_lw",
      &                "Thermal IR radiative flux to surface","W.m-2",2,
@@ -2752,5 +2760,5 @@
       comm_SWDOWNZ(1:ngrid) = fluxsurf_sw_tot(1:ngrid)
       !state  real  TAU_DUST   ij   misc  1  -  h  "TAU_DUST"  "REFERENCE VISIBLE DUST OPACITY"  ""
-      comm_TAU_DUST(1:ngrid) = tauref(1:ngrid)
+      comm_TAU_DUST(1:ngrid) = tau_pref_gcm(1:ngrid)
       !state  real  RDUST     ikj   misc  1  -  h  "RDUST"     "DUST RADIUS"                     "m"
       comm_RDUST(1:ngrid,1:nlayer) = rdust(1:ngrid,1:nlayer)
@@ -3086,6 +3094,11 @@
 c        ----------------------------------------------------------
 
-      call WRITEDIAGFI(ngrid,'tauref',
-     &                 'Dust ref opt depth','NU',2,tauref)
+      call WRITEDIAGFI(ngrid,'tau_pref_scenario',
+     &                 'Prescribed visible dust optical depth at 610Pa',
+     &                 'NU',2,tau_pref_scenario)
+
+      call WRITEDIAGFI(ngrid,'tau_pref_gcm',
+     &                 'Visible dust optical depth at 610Pa in the GCM',
+     &                 'NU',2,tau_pref_gcm)
 
       if (tracer.and.(dustbin.ne.0)) then
@@ -3382,6 +3395,7 @@
 
 #ifndef MESOSCALE
-         write(*,'("Ls =",f11.6," tauref(",f4.0," Pa) =",f9.6)')
-     &    zls*180./pi,odpref,tauref
+         write(*,
+     &    '("Ls =",f11.6," tau_pref_scenario(",f4.0," Pa) =",f9.6)')
+     &    zls*180./pi,odpref,tau_pref_scenario
 c      ----------------------------------------------------------------------
 c      Output in grads file "g1d" (ONLY when using testphys1d)
@@ -3491,6 +3505,10 @@
      &                        'total mass of dust only',
      &                        'kg.m-2',0,mdusttot)
-             call WRITEDIAGFI(ngrid,'tauref',
-     &                    'Dust ref opt depth','NU',0,tauref)
+             call WRITEDIAGFI(ngrid,'tau_pref_scenario',
+     &                        'Prescribed dust ref opt depth at 610 Pa',
+     &                        'NU',0,tau_pref_scenario)
+             call WRITEDIAGFI(ngrid,'tau_pref_gcm',
+     &                        'Dust ref opt depth at 610 Pa in the GCM',
+     &                        'NU',0,tau_pref_gcm)
              call WRITEDIAGFI(ngrid,'rdsdqsdust',
      &                        'deposited surface stormdust mass',
Index: trunk/LMDZ.MARS/libf/phymars/read_dust_scenario.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/read_dust_scenario.F90	(revision 2414)
+++ trunk/LMDZ.MARS/libf/phymars/read_dust_scenario.F90	(revision 2415)
@@ -1,3 +1,3 @@
-subroutine read_dust_scenario(ngrid,nlayer,zday,pplev,tauref)
+subroutine read_dust_scenario(ngrid,nlayer,zday,pplev,tau_pref_scenario)
 
 ! Reading of the dust scenario file
@@ -6,4 +6,5 @@
 use geometry_mod, only: latitude, longitude ! in radians
 use datafile_mod, only: datadir
+use dust_param_mod, only: odpref
 implicit none
 
@@ -13,5 +14,6 @@
 real, intent(in) :: zday ! date in martian days
 real, dimension(ngrid,nlayer+1), intent(in) :: pplev
-real, dimension(ngrid), intent(out) :: tauref
+real, dimension(ngrid), intent(out) :: tau_pref_scenario ! visible dust column
+                                       ! opacity at odpref from scenario
 
 ! Local variables
@@ -333,5 +335,5 @@
    tau= dlat*(dlon*(tau1(2)+tau1(3)-tau1(1)-tau1(4))+tau1(1)-tau1(3)) +dlon*(tau1(4)-tau1(3))+tau1(3)
 
-   tauref(ig)=tau
+   tau_pref_scenario(ig)=tau
 ! 
 enddo ! of ig=1,ngrid
@@ -342,5 +344,5 @@
   ! - the 1.3 conversion factor from IR absorbtion opacity to
   !   IR extinction opacity
-  tauref(:)=tauref(:)*1.3*(610./700.)
+  tau_pref_scenario(:)=tau_pref_scenario(:)*1.3*(odpref/700.)
 endif
 
@@ -349,5 +351,5 @@
  ! needed to decrease opacity (*0.825) to compensate overestimation of
  ! heating rates
-  tauref(:)=tauref(:)*0.825/1.3
+  tau_pref_scenario(:)=tau_pref_scenario(:)*0.825/1.3
 endif
 
Index: trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90	(revision 2414)
+++ trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90	(revision 2415)
@@ -31,5 +31,5 @@
 !             output
                                  pdqrds,wrad,dsodust,dsords,dsotop,    &
-                                 tauref)
+                                 tau_pref_scenario,tau_pref_gcm)
 
       USE tracer_mod, only: igcm_stormdust_mass,igcm_stormdust_number &
@@ -92,9 +92,11 @@
       REAL, INTENT(OUT) :: pdqrds(ngrid,nlayer,nq) ! tendancy field for dust when detraining
       REAL, INTENT(OUT) :: wrad(ngrid,nlayer+1)   ! vertical speed within the rocket dust storm
-      REAL, INTENT(INOUT) :: dsodust(ngrid,nlayer) ! density scaled opacity of env. dust
-      REAL, INTENT(INOUT) :: dsords(ngrid,nlayer) ! density scaled opacity of storm dust
-      REAL, INTENT(INOUT) :: dsotop(ngrid,nlayer) ! density scaled opacity of topmons dust
-      REAL, INTENT(OUT) :: tauref(ngrid)
-
+      REAL, INTENT(OUT) :: dsodust(ngrid,nlayer) ! density scaled opacity of env. dust
+      REAL, INTENT(OUT) :: dsords(ngrid,nlayer) ! density scaled opacity of storm dust
+      REAL, INTENT(OUT) :: dsotop(ngrid,nlayer) ! density scaled opacity of topmons dust
+      REAL,INTENT(OUT) :: tau_pref_scenario(ngrid) ! prescribed dust column 
+                                               ! visible opacity at odpref
+      REAL,INTENT(OUT) :: tau_pref_gcm(ngrid) ! dust column visible opacity at
+                                              ! odpref in the GCM
 !--------------------------------------------------------
 ! Local variables 
@@ -246,5 +248,6 @@
                  emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,      &
                  zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1,     &
-                 fluxtop_sw1,tauref,tau,aerosol,dsodust,tauscaling,       &
+                 fluxtop_sw1,tau_pref_scenario,tau_pref_gcm, &
+                 tau,aerosol,dsodust,tauscaling,       &
                  taucloudtes,rdust,rice,nuice,co2ice,rstormdust,rtopdust, &
                  totstormfract,clearatm,dsords,dsotop,alpha_hmons,nohmons,&
Index: trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90
===================================================================
--- trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90	(revision 2414)
+++ trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90	(revision 2415)
@@ -33,5 +33,5 @@
 !             output
                                  pdqtop,wfin,dsodust,dsords,dsotop,    &
-                                 tauref)
+                                 tau_pref_scenario,tau_pref_gcm)
 
       USE tracer_mod, only: igcm_topdust_mass,igcm_topdust_number &
@@ -96,8 +96,11 @@
 !     output for second radiative transfer 
       REAL, INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind)
-      REAL, INTENT(INOUT) :: dsodust(ngrid,nlayer)
-      REAL, INTENT(INOUT) :: dsords(ngrid,nlayer)
-      REAL, INTENT(INOUT) :: dsotop(ngrid,nlayer)
-      REAL, INTENT(OUT) :: tauref(ngrid)
+      REAL, INTENT(OUT) :: dsodust(ngrid,nlayer)
+      REAL, INTENT(OUT) :: dsords(ngrid,nlayer)
+      REAL, INTENT(OUT) :: dsotop(ngrid,nlayer)
+      REAL,INTENT(OUT) :: tau_pref_scenario(ngrid) ! prescribed dust column 
+                                               ! visible opacity at odpref
+      REAL,INTENT(OUT) :: tau_pref_gcm(ngrid) ! dust column visible opacity at
+                                              ! odpref in the GCM
 
 !--------------------------------------------------------
@@ -271,5 +274,6 @@
                  emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout,   &
                  zdtlw1,zdtsw1,fluxsurf_lw1,fluxsurf_sw1,fluxtop_lw1,  &
-                 fluxtop_sw1,tauref,tau,aerosol,dsodust,tauscaling,    &
+                 fluxtop_sw1,tau_pref_scenario,tau_pref_gcm, &
+                 tau,aerosol,dsodust,tauscaling,    &
                  taucloudtes,rdust,rice,nuice,co2ice,rstormdust,rtopdust, &
                  totstormfract,clearatm,dsords,dsotop,alpha_hmons,nohmons,&
