Index: /trunk/LMDZ.MARS/README
===================================================================
--- /trunk/LMDZ.MARS/README	(revision 2245)
+++ /trunk/LMDZ.MARS/README	(revision 2246)
@@ -2861,2 +2861,6 @@
 == 10/02/2020 == CM
 Initialization of mu0_int and ztim1 in 'nirco2abs.F' to prevent crash during debug mode simulations.
+
+== 27/02/2020 == AB
+Resolved ticket #32 : 1) dsodust is now calculated only once in the InfraRed by aeropacity_mod (used to be wrongly calculated twice, such as dsodust=IR_part+Visible_part) ; 2) dsords is now calculated in the IR by aeropacity_mod (used to be calculated in the Visible) ; 3) dsotop is added and calculated in the IR in aeropacity_mod
++ some cleaning and commenting of the code
Index: /trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F	(revision 2245)
+++ /trunk/LMDZ.MARS/libf/phymars/aeropacity_mod.F	(revision 2246)
@@ -8,5 +8,5 @@
      &    pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
      &    QREFvis3d,QREFir3d,omegaREFir3d,
-     &    totstormfract,clearatm,dsords,
+     &    totstormfract,clearatm,dsords,dsotop,
      &    alpha_hmons,nohmons,
      &    clearsky,totcloudfrac)
@@ -90,5 +90,6 @@
       REAL, INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind)
       REAL, INTENT(OUT) ::  dsodust(ngrid,nlayer)
-      REAL, INTENT(OUT) ::  dsords(ngrid,nlayer) !dso of stormdust  
+      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)
@@ -171,4 +172,5 @@
       dsords(:,:)=0. !CW17: initialize dsords
       dsodust(:,:)=0.
+      dsotop(:,:)=0.
 
 ! identify tracers
@@ -376,5 +378,5 @@
 
 c==================================================================
-        CASE("dust_doubleq") aerkind! Two-moment scheme for dust
+        CASE("dust_doubleq") aerkind! Two-moment scheme for background dust
 c        (transport of mass and number mixing ratio)
 c==================================================================
@@ -385,4 +387,7 @@
 c             avoid unrealistic values due to constant lifting:
               DO ig=1,ngrid
+	      ! OPTICAL DEPTH for the computation of tauref, 
+	      ! which is to be compared with tauref_scenario
+	      ! => visible wavelength
                 aerosol(ig,l,iaer) = 
      &          (  0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
@@ -390,5 +395,7 @@
      &          pq(ig,cstdustlevel,igcm_dust_mass) *
      &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
-                ! DENSITY SCALED OPACITY IN THE INFRARED:
+              ! DENSITY SCALED OPACITY :
+	      ! GCM output to be compared with observations
+	      ! => infrared wavelength
                 dsodust(ig,l) =
      &          (  0.75 * QREFir3d(ig,cstdustlevel,iaer) /
@@ -398,10 +405,15 @@
             ELSE
               DO ig=1,ngrid
-                aerosol(ig,l,iaer) =
+              ! OPTICAL DEPTH for the computation of tauref, 
+	      ! which is to be compared with tauref_scenario
+	      ! => visible wavelength 
+	        aerosol(ig,l,iaer) =
      &          (  0.75 * QREFvis3d(ig,l,iaer) /
      &          ( rho_dust * reffrad(ig,l,iaer) )  ) *
      &          pq(ig,l,igcm_dust_mass) *
      &          ( pplev(ig,l) - pplev(ig,l+1) ) / g
-                ! DENSITY SCALED OPACITY IN THE INFRARED:
+              ! DENSITY SCALED OPACITY :
+	      ! GCM output to be compared with observations
+	      ! => infrared wavelength
                 dsodust(ig,l) =
      &          (  0.75 * QREFir3d(ig,l,iaer) /
@@ -505,12 +517,25 @@
 c           avoid unrealistic values due to constant lifting:
                DO ig=1,ngrid
-                  aerosol(ig,l,iaer) = 
+	       ! OPTICAL DEPTH for the computation of tauref, 
+	       ! which is to be compared with tauref_scenario
+	       ! => visible wavelength
+                 aerosol(ig,l,iaer) = 
      &           (  0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
      &           ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
      &           pq(ig,cstdustlevel,igcm_stormdust_mass) *
      &           ( pplev(ig,l) - pplev(ig,l+1) ) / g
+               ! DENSITY SCALED OPACITY :
+	       ! GCM output to be compared with observations
+	       ! => infrared wavelength
+                 dsords(ig,l) =
+     &           (  0.75 * QREFir3d(ig,cstdustlevel,iaer) /
+     &           ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
+     &           pq(ig,cstdustlevel,igcm_stormdust_mass)
                ENDDO
              ELSE
-              DO ig=1,ngrid
+               DO ig=1,ngrid
+	       ! OPTICAL DEPTH for the computation of tauref, 
+	       ! which is to be compared with tauref_scenario
+	       ! => visible wavelength
                  aerosol(ig,l,iaer) =
      &           (  0.75 * QREFvis3d(ig,l,iaer) /
@@ -518,5 +543,12 @@
      &           pq(ig,l,igcm_stormdust_mass) *
      &           ( pplev(ig,l) - pplev(ig,l+1) ) / g
-              ENDDO
+               ! DENSITY SCALED OPACITY :
+	       ! GCM output to be compared with observations
+	       ! => infrared wavelength
+                 dsords(ig,l) =
+     &           (  0.75 * QREFir3d(ig,l,iaer) /
+     &           ( rho_dust * reffrad(ig,l,iaer) )  ) *
+     &           pq(ig,l,igcm_stormdust_mass)
+               ENDDO
              ENDIF
           ENDDO
@@ -533,24 +565,40 @@
         ELSE  ! part of the mesh with concentred dust storm
           DO l=1,nlayer
-!             IF (l.LE.cstdustlevel) THEN
-!c          Opacity in the first levels is held constant to 
-!c           avoid unrealistic values due to constant lifting:
-!               DO ig=1,ngrid
-!                  aerosol(ig,l,iaer) = 
-!     &           (  0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
-!     &           ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
-!     &           pq(ig,cstdustlevel,igcm_topdust_mass) *
-!     &           ( pplev(ig,l) - pplev(ig,l+1) ) / g
-!               ENDDO
-!             ELSE
-              DO ig=1,ngrid
-                 aerosol(ig,l,iaer) =
+             IF (l.LE.cstdustlevel) THEN
+c          Opacity in the first levels is held constant to 
+c           avoid unrealistic values due to constant lifting:
+               DO ig=1,ngrid
+                  aerosol(ig,l,iaer) = 
+     &           (  0.75 * QREFvis3d(ig,cstdustlevel,iaer) /
+     &           ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
+     &           pq(ig,cstdustlevel,igcm_topdust_mass) *
+     &           ( pplev(ig,l) - pplev(ig,l+1) ) / g
+               ! DENSITY SCALED OPACITY :
+	       ! GCM output to be compared with observations
+	       ! => infrared wavelength
+                 dsotop(ig,l) =
+     &           (  0.75 * QREFir3d(ig,cstdustlevel,iaer) /
+     &           ( rho_dust * reffrad(ig,cstdustlevel,iaer) )  ) *
+     &           pq(ig,cstdustlevel,igcm_topdust_mass)
+               ENDDO
+             ELSE
+               DO ig=1,ngrid
+               ! OPTICAL DEPTH for the computation of tauref, 
+	       ! which is to be compared with tauref_scenario
+	       ! => visible wavelength
+	         aerosol(ig,l,iaer) =
      &           (  0.75 * QREFvis3d(ig,l,iaer) /
      &           ( rho_dust * reffrad(ig,l,iaer) )  ) *
      &           pq(ig,l,igcm_topdust_mass) *
      &           ( pplev(ig,l) - pplev(ig,l+1) ) / g
-              ENDDO
-!             ENDIF
-
+               ! DENSITY SCALED OPACITY :
+	       ! GCM output to be compared with observations
+	       ! => infrared wavelength
+                 dsotop(ig,l) =
+     &           (  0.75 * QREFir3d(ig,l,iaer) /
+     &           ( rho_dust * reffrad(ig,l,iaer) )  ) *
+     &           pq(ig,l,igcm_topdust_mass)
+               ENDDO
+             ENDIF
           ENDDO
         ENDIF
@@ -852,31 +900,4 @@
 #endif
 
-c -----------------------------------------------------------------
-c Density scaled opacity and column opacity output
-c -----------------------------------------------------------------
-          IF (rdstorm) then
-            DO l=1,nlayer
-              IF (l.LE.cstdustlevel) THEN
-                DO ig=1,ngrid 
-                  dsodust(ig,l)=dsodust(ig,l) +
-     &                      aerosol(ig,l,iaer_dust_doubleq) * g /
-     &                      (pplev(ig,l) - pplev(ig,l+1))
-
-                  dsords(ig,l) = dsords(ig,l) +
-     &              aerosol(ig,l,iaer_stormdust_doubleq)* g/
-     &              (pplev(ig,l) - pplev(ig,l+1))
-                ENDDO
-              ELSE
-                DO ig=1,ngrid
-                  dsodust(ig,l) =dsodust(ig,l) +
-     &                      aerosol(ig,l,iaer_dust_doubleq) * g /
-     &                      (pplev(ig,l) - pplev(ig,l+1))
-                   dsords(ig,l) = dsords(ig,l) +
-     &                        aerosol(ig,l,iaer_stormdust_doubleq)* g/
-     &                        (pplev(ig,l) - pplev(ig,l+1))
-                ENDDO
-              ENDIF
-            ENDDO
-          ENDIF
 
 c -----------------------------------------------------------------
Index: /trunk/LMDZ.MARS/libf/phymars/callradite_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/callradite_mod.F	(revision 2245)
+++ /trunk/LMDZ.MARS/libf/phymars/callradite_mod.F	(revision 2246)
@@ -10,5 +10,5 @@
      $     fluxtop_sw,tauref,tau,aerosol,dsodust,tauscaling,
      $     taucloudtes,rdust,rice,nuice,co2ice,rstormdust,rtopdust,
-     $     totstormfract,clearatm,dsords,alpha_hmons,nohmons,
+     $     totstormfract,clearatm,dsords,dsotop,alpha_hmons,nohmons,
      $     clearsky,totcloudfrac)
 
@@ -35,6 +35,6 @@
 c
 c   These calculations are only valid on the part of the atmosphere
-c   where Local Thermal Equilibrium (NLTE) is verified. In practice
-c   The calculations are only performed for the first "nlaylte"
+c   where Local Thermal Equilibrium (LTE) is verified. In practice
+c   the calculations are only performed for the first "nlaylte"
 c   parameters (nlaylte is calculated by subroutine "nlthermeq"
 c   and stored in module "yomlw_h").
@@ -203,11 +203,12 @@
       REAL,INTENT(IN) :: totstormfract(ngrid) ! dust storm mesh fraction
       REAL,INTENT(OUT) :: rstormdust(ngrid,nlayer)  ! Storm dust geometric mean radius (m)
-      REAL 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
       LOGICAL, INTENT(IN) :: nohmons ! true for background dust 
       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(OUT) :: dsotop(ngrid,nlayer) ! density scaled opacity for topmons dust
+      
 c     sub-grid scale water ice clouds
       LOGICAL,INTENT(IN) :: clearsky
@@ -269,5 +270,7 @@
 
       REAL :: QREFvis3d(ngrid,nlayer,naerkind)
+      ! QREFvis3d : Extinction efficiency at the VISible reference wavelength
       REAL :: QREFir3d(ngrid,nlayer,naerkind)
+      ! QREFir3d : Extinction efficiency at the InfraRed reference wavelength
 
       REAL :: omegaREFvis3d(ngrid,nlayer,naerkind)
@@ -421,5 +424,5 @@
      &    pq,tauscaling,tauref,tau,taucloudtes,aerosol,dsodust,reffrad,
      &    QREFvis3d,QREFir3d,omegaREFir3d,
-     &    totstormfract,clearatm,dsords,
+     &    totstormfract,clearatm,dsords,dsotop,
      &    alpha_hmons,nohmons,
      &    clearsky,totcloudfrac)
Index: /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2245)
+++ /trunk/LMDZ.MARS/libf/phymars/physiq_mod.F	(revision 2246)
@@ -407,4 +407,5 @@
       REAL dsodust(ngrid,nlayer) 
       REAL dsords(ngrid,nlayer)
+      REAL dsotop(ngrid,nlayer)
       REAL wtop(ngrid,nlayer+1) ! vertical velocity topdust tracer
 
@@ -855,5 +856,5 @@
      &     fluxtop_sw,tauref,tau,aerosol,dsodust,tauscaling,
      &     taucloudtes,rdust,rice,nuice,co2ice,rstormdust,rtopdust,
-     &     totstormfract,clearatm,dsords,alpha_hmons,nohmons,
+     &     totstormfract,clearatm,dsords,dsotop,alpha_hmons,nohmons,
      &     clearsky,totcloudfrac)
 
@@ -871,5 +872,5 @@
      &              tau,aerosol,dsodust,tauscaling,taucloudtesclf,rdust,
      &              rice,nuice,co2ice,rstormdust,rtopdust,totstormfract,
-     &              clearatm,dsords,alpha_hmons,nohmons,
+     &              clearatm,dsords,dsotop,alpha_hmons,nohmons,
      &              clearsky,totcloudfrac)
                clearsky = .false.  ! just in case.
@@ -1069,5 +1070,5 @@
      &                      nohmons,alpha_hmons,
 c               output
-     &                      pdqrds,wspeed,dsodust,dsords,
+     &                      pdqrds,wspeed,dsodust,dsords,dsotop,
      &                      tauref)
                    
@@ -1129,8 +1130,8 @@
      &                zzlay,zdtsw,zdtlw,
      &                icount,zday,zls,tsurf,igout,aerosol, 
-     &                tauscaling,totstormfract,clearatm,dsords,
+     &                tauscaling,totstormfract,clearatm,
      &                clearsky,totcloudfrac,
      &                nohmons,hsummit,
-     &                pdqtop,wtop,dsodust,
+     &                pdqtop,wtop,dsodust,dsords,dsotop,
      &                tauref)
      
@@ -2924,10 +2925,11 @@
      &                        'part/kg',3,ndust)
              call WRITEDIAGFI(ngrid,'dsodust',
-     &                       'density scaled optical depth in the IR',
-     &                       'm2.kg-1',3,dsodust)
+     &                        'density scaled optical depth in the IR',
+     &                        'm2.kg-1',3,dsodust)
              call WRITEDIAGFI(ngrid,'dso',
-     &                       'density scaled optical depth in the IR',
-     &                       'm2.kg-1',3,dsodust+dsords)
-           else
+     &                        'density scaled optical depth in the IR',
+     &                        'm2.kg-1',3,dsodust+dsords+dsotop)
+	     
+           else ! (doubleq=.false.)
              do iq=1,dustbin
                write(str2(1:2),'(i2.2)') iq
@@ -2990,4 +2992,10 @@
            endif ! (rdstorm)
 
+           if (slpwind) then
+	     call WRITEDIAGFI(ngrid,'dsotop',
+     &                       'density scaled opacity of topdust',
+     &                       'm2.kg-1',3,dsotop)
+           endif ! (slpwind)
+	   
            if (scavenging) then
              call WRITEDIAGFI(ngrid,'ccnq','CCN mass mr',
Index: /trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90	(revision 2245)
+++ /trunk/LMDZ.MARS/libf/phymars/rocketduststorm_mod.F90	(revision 2246)
@@ -30,5 +30,5 @@
                                  nohmons,alpha_hmons,                  &
 !             output
-                                 pdqrds,wrad,dsodust,dsords,           &
+                                 pdqrds,wrad,dsodust,dsords,dsotop,    &
                                  tauref)
 
@@ -94,4 +94,5 @@
       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) :: tauref(ngrid)
 
@@ -247,5 +248,5 @@
                  fluxtop_sw1,tauref,tau,aerosol,dsodust,tauscaling,       &
                  taucloudtes,rdust,rice,nuice,co2ice,rstormdust,rtopdust, &
-                 totstormfract,clearatm,dsords,alpha_hmons,nohmons,       &
+                 totstormfract,clearatm,dsords,dsotop,alpha_hmons,nohmons,&
                  clearsky,totcloudfrac)
 
Index: /trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90
===================================================================
--- /trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90	(revision 2245)
+++ /trunk/LMDZ.MARS/libf/phymars/topmons_mod.F90	(revision 2246)
@@ -26,5 +26,5 @@
                                  tauscaling,                           &
 !             input sub-grid scale rocket dust storm
-                                 totstormfract,clearatm,dsords,        &
+                                 totstormfract,clearatm,               &
 !             input sub-grid scale cloud
                                  clearsky,totcloudfrac,                &
@@ -32,5 +32,5 @@
                                  nohmons,hsummit,                      &
 !             output
-                                 pdqtop,wfin,dsodust,                   &
+                                 pdqtop,wfin,dsodust,dsords,dsotop,    &
                                  tauref)
 
@@ -98,4 +98,5 @@
       REAL, INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind)
       REAL, INTENT(OUT) :: dsords(ngrid,nlayer)
+      REAL, INTENT(OUT) :: dsotop(ngrid,nlayer)
       REAL, INTENT(OUT) :: tauref(ngrid)
 
@@ -211,4 +212,5 @@
       dsodust(:,:)=0.
       dsords(:,:)=0.
+      dsotop(:,:)=0.
       pdqtop(:,:,:) = 0.
       dqvl_topdust_mass(:,:)=0.
@@ -274,5 +276,5 @@
                  fluxtop_sw1,tauref,tau,aerosol,dsodust,tauscaling,    &
                  taucloudtes,rdust,rice,nuice,co2ice,rstormdust,rtopdust, &
-                 totstormfract,clearatm,dsords,alpha_hmons,nohmons,    &
+                 totstormfract,clearatm,dsords,dsotop,alpha_hmons,nohmons,&
                  clearsky,totcloudfrac)
        ! **********************************************************************
