Index: trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 1896)
+++ trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 1897)
@@ -45,5 +45,5 @@
       INTEGER,INTENT(IN) :: ngrid                  ! Number of atmospheric columns.
       INTEGER,INTENT(IN) :: nlayer                 ! Number of atmospheric layers.
-      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)       ! Tracers (kg/kg_of_air).
+      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)       ! Tracers (X/m2).
       INTEGER,INTENT(IN) :: nq                     ! Number of tracers.
       REAL,INTENT(IN) :: qsurf(ngrid,nq)           ! Tracers on surface (kg.m-2).
@@ -298,5 +298,5 @@
             endif
             
-            call optcv(pq,plevrad,tmid,pmid,                       &
+            call optcv(pq(ig,:,1:nmicro),nlayer,plevrad,tmid,pmid,        &
                  dtauv,tauv,taucumv,wbarv,cosbv,tauray,taugsurf)
 
@@ -340,5 +340,5 @@
 !-----------------------------------------------------------------------
 
-         call optci(pq,plevrad,tlevrad,tmid,pmid,                   &
+         call optci(pq(ig,:,1:nmicro),nlayer,plevrad,tlevrad,tmid,pmid,    &
               dtaui,taucumi,cosbi,wbari,taugsurfi)
 
Index: trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 1896)
+++ trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 1897)
@@ -12,6 +12,9 @@
       logical,save :: callgasvis,continuum,graybody
 !$OMP THREADPRIVATE(callgasvis,continuum,graybody)
-      logical,save :: strictboundcorrk                                     
+      logical,save :: strictboundcorrk
 !$OMP THREADPRIVATE(strictboundcorrk)
+      logical,save :: uncoupl_optic_haze
+!$OMP THREADPRIVATE(uncoupl_optic_haze)
+
       logical,save :: callchim, callmufi, callclouds
 !$OMP THREADPRIVATE(callchim,callmufi,callclouds)
@@ -45,9 +48,8 @@
 !$OMP THREADPRIVATE(iddist,iradia,startype)
       
-      real,save :: df_mufi, rm_mufi, rho_aer_mufi
       real,save :: p_prod, tx_prod, rc_prod
       real,save :: air_rad
-!$OMP THREADPRIVATE(df_mufi, rm_mufi, rho_aer_mufi,p_prod,tx_prod,rc_prod,air_rad)
-
+!$OMP THREADPRIVATE(p_prod,tx_prod,rc_prod,air_rad)
+      
       real,save :: szangle
 !$OMP THREADPRIVATE(szangle)
Index: trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90	(revision 1896)
+++ trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90	(revision 1897)
@@ -1,6 +1,8 @@
+
+
 SUBROUTINE calmufi(plev, zlev, play, zlay, temp, pq, zdq)
   !! Interface subroutine to YAMMS model for Titan LMDZ GCM.
   !!
-  !! The subroutine computes the microphysics processes for a single vertical column.
+  !! The subroutine computes the microphysics processes for a si:le vertical column.
   !!
   !! - All input vectors are assumed to be defined from GROUND to TOP of the atmosphere.
@@ -52,5 +54,4 @@
 
   INTEGER :: ilon, i,nices
-
   INTEGER :: nlon,nlay
 
@@ -89,6 +90,5 @@
     ! Convert tracers to extensive ( except for gazs where we work with molar mass ratio )
     ! We suppose a given order of tracers !
-    int2ext(:) = ( plev(ilon,1:nlay) - plev(ilon,2:nlay+1) ) / g
-
+    int2ext(:) = ( plev(ilon,1:nlay)-plev(ilon,2:nlay+1) ) / g
     m0as(:) = pq(ilon,:,1) * int2ext(:)
     m3as(:) = pq(ilon,:,2) * int2ext(:)
@@ -97,13 +97,12 @@
     
     if (callclouds) then ! if call clouds
-      dm0n(:) = pq(ilon,:,5) * int2ext(:)
-      dm3n(:) = pq(ilon,:,6) * int2ext(:)
+      m0n(:) = pq(ilon,:,5) * int2ext(:)
+      m3n(:) = pq(ilon,:,6) * int2ext(:)
       do i=1,nices
-        dm3i(:,nices) = pq(ilon,:,6+i) * int2ext(:)
-        dgazs(:,i)    = pq(ilon,:,ices_indx(i)) * rat_mmol(ices_indx(i)) ! For gazs we work on the full tracer array !!
+        m3i(:,nices) = pq(ilon,:,6+i) * int2ext(:)
+        gazs(:,i)    = pq(ilon,:,ices_indx(i)) * rat_mmol(ices_indx(i)) ! For gazs we work on the full tracer array !!
         ! We use the molar mass ratio from GCM in case there is discrepancy with the mm one
       enddo
     endif
-
 
     ! Initialize YAMMS atmospheric column
@@ -135,5 +134,5 @@
     ! Convert tracers back to intensives ( except for gazs where we work with molar mass ratio )
     ! We suppose a given order of tracers !
-
+ 
     zdq(ilon,:,1) = dm0as(:) / int2ext(:)
     zdq(ilon,:,2) = dm3as(:) / int2ext(:)
@@ -150,5 +149,4 @@
       enddo
     endif
-
   END DO ! loop on ilon
 
Index: trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 1896)
+++ trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 1897)
@@ -377,18 +377,9 @@
      endif
 
-     write(*,*) "Fractal dimension ?"
-     df_mufi=2.0 ! default value
-     call getin_p("df_mufi",df_mufi)
-     write(*,*)" df_mufi = ",df_mufi
-
-     write(*,*) "Monomer radius (m) ?"
-     rm_mufi=6.66e-08 ! default value
-     call getin_p("rm_mufi",rm_mufi)
-     write(*,*)" rm_mufi = ",rm_mufi
-
-     write(*,*) "Aerosol density (kg.m-3)?"
-     rho_aer_mufi=1.e3 ! default value
-     call getin_p("rho_aer_mufi",rho_aer_mufi)
-     write(*,*)" rho_aer_mufi = ",rho_aer_mufi
+     write(*,*) "Disable the coupling of microphysics within rad. transf. ?"
+     write(*,*) "If disabled we will assume a planetwide vert. profile of extinction ..."
+     uncoupl_optic_haze=.true. ! default value - true as long as the microphysics is bugged
+     call getin_p("uncoupl_optic_haze",uncoupl_optic_haze)
+     write(*,*)" uncoupl_optic_haze = ",uncoupl_optic_haze
 
      write(*,*) "Pressure level of aer. production (Pa) ?"
Index: trunk/LMDZ.TITAN/libf/phytitan/inimufi.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/inimufi.F90	(revision 1896)
+++ trunk/LMDZ.TITAN/libf/phytitan/inimufi.F90	(revision 1897)
@@ -2,6 +2,5 @@
 
   use mmp_gcm
-  use callkeys_mod, only : callclouds, df_mufi, &
-        rm_mufi, rho_aer_mufi, p_prod, tx_prod, rc_prod, air_rad
+  use callkeys_mod, only : callclouds, p_prod, tx_prod, rc_prod, air_rad
   use tracer_h
   use comcstfi_mod, only : g, rad, mugaz
@@ -42,4 +41,7 @@
   logical :: err
 
+  ! PATCH : YAMMS now allows to enable/disable effective g computations:
+  !  In the library it defaults to .true., in the GCM we do not want it !
+  mm_use_effg = .false.
 
   !----------------------------------------------------
@@ -47,5 +49,5 @@
   ! ---------------------------------------------------
 
-   call mmp_initialize(ptimestep,df_mufi,rm_mufi,rho_aer_mufi,p_prod,tx_prod,rc_prod, &
+   call mmp_initialize(ptimestep,p_prod,tx_prod,rc_prod, &
         rad,g,air_rad,mugaz,callclouds,config_mufi)
 
Index: trunk/LMDZ.TITAN/libf/phytitan/optci.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 1896)
+++ trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 1897)
@@ -1,3 +1,3 @@
-subroutine optci(PQ,PLEV,TLEV,TMID,PMID,      &
+subroutine optci(PQO,NLAY,PLEV,TLEV,TMID,PMID,      &
      DTAUI,TAUCUMI,COSBI,WBARI,TAUGSURF)
 
@@ -6,5 +6,8 @@
   use gases_h
   use comcstfi_mod, only: g, r
-  use callkeys_mod, only: continuum,graybody
+  use callkeys_mod, only: continuum,graybody,callclouds,callmufi, uncoupl_optic_haze
+  use tracer_h, only : nmicro,nice
+  use MMP_OPTICS
+
   implicit none
 
@@ -34,5 +37,6 @@
   ! Input/Output
   !==========================================================
-  REAL*8, INTENT(IN)  :: PQ ! Tracers (kg/kg_of_air).
+  REAL*8, INTENT(IN)  :: PQO(nlay,nmicro) ! Tracers (X/m2).
+  INTEGER, INTENT(IN) :: NLAY             ! Number of pressure layers (for pqo)
   REAL*8, INTENT(IN)  :: PLEV(L_LEVELS), TLEV(L_LEVELS)
   REAL*8, INTENT(IN)  :: TMID(L_LEVELS), PMID(L_LEVELS)
@@ -86,4 +90,10 @@
   integer interm
 
+  real*8 m0as,m3as,m0af,m3af
+  real*8 ext_s,sca_s,ssa_s,asf_s
+  real*8 ext_f,sca_f,ssa_f,asf_f
+  logical,save :: firstcall=.true.
+  !$OMP THREADPRIVATE(firstcall)
+
   !! AS: to save time in computing continuum (see bilinearbig)
   IF (.not.ALLOCATED(indi)) THEN
@@ -127,7 +137,26 @@
         ilay = k / 2 ! int. arithmetic => gives the gcm layer index
         
-        !================= Titan customisation ========================================
-        call disr_haze(dz(k),plev(k),wnoi(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
-        ! =============================================================================
+        ! Optical coupling of YAMMS is plugged but inactivated for now
+        ! as long as the microphysics only isn't fully debugged -- JVO 01/18
+        IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
+          m0as = pqo(ilay,1)
+          m3as = pqo(ilay,2)
+          m0af = pqo(ilay,3)
+          m3af = pqo(ilay,4)
+
+          IF (.NOT.mmp_sph_optics_ir(m0as,m3as,nw,ext_s,sca_s,ssa_s,asf_s)) &
+          CALL abort_gcm("optcv", "Fatal error in mmp_sph_optics_ir", 12)
+          IF (.NOT.mmp_fra_optics_ir(m0af,m3af,nw,ext_f,sca_f,ssa_f,asf_f)) &
+          CALL abort_gcm("optcv", "Fatal error in mmp_fra_optics_ir", 12)
+          dhaze_T(k,nw) = ext_s+ext_f
+          SSA_T(k,nw)   = (sca_s+sca_f)/dhaze_T(k,nw)
+          ASF_T(k,nw)   = (asf_s*sca_s + asf_f*sca_f) /(sca_s+sca_f)
+          IF (callclouds.and.firstcall) & 
+            WRITE(*,*) 'WARNING: In optci, optical properties &
+                       &calculations are not implemented yet'
+        ELSE
+          ! Call fixed vertical haze profile of extinction - same for all columns
+          call disr_haze(dz(k),plev(k),wnoi(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
+        ENDIF
 
         DCONT = 0.0d0 ! continuum absorption
@@ -401,4 +430,6 @@
 ! ==============================================================================  
 
+  if(firstcall) firstcall = .false.
+
   return
 
Index: trunk/LMDZ.TITAN/libf/phytitan/optcv.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 1896)
+++ trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 1897)
@@ -1,3 +1,3 @@
-SUBROUTINE OPTCV(PQ,PLEV,TMID,PMID,  &
+SUBROUTINE OPTCV(PQO,NLAY,PLEV,TMID,PMID,  &
      DTAUV,TAUV,TAUCUMV,WBARV,COSBV,TAURAY,TAUGSURF)
 
@@ -6,5 +6,7 @@
   use gases_h
   use comcstfi_mod, only: g, r
-  use callkeys_mod, only: continuum,graybody,callgasvis
+  use callkeys_mod, only: continuum,graybody,callgasvis,callclouds,callmufi,uncoupl_optic_haze
+  use tracer_h, only: nmicro,nice
+  use MMP_OPTICS
 
   implicit none
@@ -41,5 +43,6 @@
   ! Input/Output
   !==========================================================
-  REAL*8, INTENT(IN)  :: PQ ! Tracers (kg/kg_of_air).
+  REAL*8, INTENT(IN)  :: PQO(nlay,nmicro) ! Tracers (X/m2).
+  INTEGER, INTENT(IN) :: NLAY             ! Number of pressure layers (for pqo)
   REAL*8, INTENT(IN)  :: PLEV(L_LEVELS)
   REAL*8, INTENT(IN)  :: TMID(L_LEVELS), PMID(L_LEVELS)
@@ -97,4 +100,11 @@
   integer interm
 
+  real*8 m0as,m3as,m0af,m3af
+  real*8 ext_s,sca_s,ssa_s,asf_s
+  real*8 ext_f,sca_f,ssa_f,asf_f
+  logical,save :: firstcall=.true.
+  !$OMP THREADPRIVATE(firstcall)
+
+
   !! AS: to save time in computing continuum (see bilinearbig)
   IF (.not.ALLOCATED(indv)) THEN
@@ -149,8 +159,27 @@
 
      do NW=1,L_NSPECTV
-
-        !================= Titan customisation ========================================
-        call disr_haze(dz(k),plev(k),wnov(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
-        ! =============================================================================
+     
+        ! Optical coupling of YAMMS is plugged but inactivated (if false) for now
+        ! as long as the microphysics only isn't fully debugged -- JVO 01/18
+        IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
+          m0as = pqo(ilay,1)
+          m3as = pqo(ilay,2)
+          m0af = pqo(ilay,3)
+          m3af = pqo(ilay,4)
+
+          IF (.NOT.mmp_sph_optics_vis(m0as,m3as,nw,ext_s,sca_s,ssa_s,asf_s)) &
+          CALL abort_gcm("optcv", "Fatal error in mmp_sph_optics_vis", 12)
+          IF (.NOT.mmp_fra_optics_vis(m0af,m3af,nw,ext_f,sca_f,ssa_f,asf_f)) &
+          CALL abort_gcm("optcv", "Fatal error in mmp_fra_optics_vis", 12)
+          dhaze_T(k,nw) = ext_s+ext_f
+          SSA_T(k,nw)   = (sca_s+sca_f)/dhaze_T(k,nw)
+          ASF_T(k,nw)   = (asf_s*sca_s + asf_f*sca_f) /(sca_s+sca_f)
+          IF (callclouds.and.firstcall) & 
+            WRITE(*,*) 'WARNING: In optcv, optical properties &
+                       &calculations are not implemented yet'
+        ELSE
+          ! Call fixed vertical haze profile of extinction - same for all columns
+          call disr_haze(dz(k),plev(k),wnov(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
+        ENDIF
          
         DRAYAER = TRAY(K,NW)
@@ -355,4 +384,5 @@
 ! ==============================================================================  
 
+  if(firstcall) firstcall = .false.
 
   return
Index: trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 1896)
+++ trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 1897)
@@ -1,2 +1,4 @@
+#define USE_QTEST
+
       module physiq_mod
       
@@ -46,4 +48,5 @@
       use wxios, only: wxios_context_init, xios_context_finalize
 #endif
+      use MMP_OPTICS
       implicit none
 
@@ -150,4 +153,5 @@
 !                + clean of all too-generic (ocean, water, co2 ...) routines
 !                + Titan's chemistry
+!           Microphysical moment model - J.Burgalat / J.Vatant d'Ollone (2017-2018)
 !============================================================================================
 
@@ -383,8 +387,18 @@
 !$OMP THREADPRIVATE(tankCH4)
 
+      ! -----******----- FOR MUPHYS OPTICS -----******----- 
+      integer :: i,j
+      real :: pqo(ngrid,nlayer,nq)   ! Tracers for the optics (X/m2).
+      real :: i2e(nlayer)            ! int 2 ext factor
+      ! -----******----- END FOR MUPHYS OPTICS -----******----- 
+
+      real,save,dimension(:,:,:), allocatable :: tpq ! Tracers for decoupled microphysical tests ( temporary in 01/18 )
+!$OMP THREADPRIVATE(tpq)
+
+
 !-----------------------------------------------------------------------------
     ! Interface to calmufi 
     !   --> needed in order to pass assumed-shape arrays. Otherwise we must put calmufi in a module 
-    !       (to have an explicit generated by the compiler).
+    !       (to have an explicit interface generated by the compiler).
     !   Or one can put calmufi in MMP_GCM module (in muphytitan).
     INTERFACE
@@ -410,4 +424,6 @@
 ! --------------------------------
       if (firstcall) then
+        allocate(tpq(ngrid,nlayer,nq))
+        tpq(:,:,:) = pq(:,:,:)
 
         ! Initialisation of nmicro as well as tracers names, indexes ...
@@ -518,4 +534,8 @@
 
            call inimufi(nq,ptimestep)
+           
+           ! Optical coupling of YAMMS is plugged but inactivated for now
+           ! as long as the microphysics only isn't fully debugged -- JVO 01/18
+           IF (.NOT.uncoupl_optic_haze) call mmp_initialize_optics("/path/to/mmp_optic_table.nc")
 
          ENDIF
@@ -820,6 +840,18 @@
                call call_profilgases(nlayer)
 
+               ! Convert (microphysical) tracers for optics: X.kg-1 --> X.m-2_
+               ! NOTE: it should be moved somewhere else: calmufi performs the same kind of
+               ! computations... waste of time...
+               DO i = 1, ngrid
+                 i2e(:) = ( pplev(i,1:nlayer)-pplev(i,2:nlayer+1) ) / g 
+                 pqo(i,:,:) = 0.0
+                 DO j=1,nmicro-nice
+                   pqo(i,:,j) = pq(i,:,j)*i2e(:)
+                 ENDDO
+               ENDDO
+
+
                ! standard callcorrk
-               call callcorrk(ngrid,nlayer,pq,nq,qsurf,                           &
+               call callcorrk(ngrid,nlayer,pqo,nq,qsurf,                          &
                               albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
                               tsurf,fract,dist_star,                              &
@@ -1108,5 +1140,12 @@
  
          if (callmufi) then
-
+#ifdef USE_QTEST
+               if (ngrid.eq.1) then ! We obviously don't have access to (and don't need) zonal means in 1D
+                  call calmufi(pplev,zzlev,pplay,zzlay,pt,tpq,zdqmufi)
+               else
+                  call calmufi(zplevbar,zzlevbar,zplaybar,zzlaybar,ztfibar,tpq,zdqmufi)
+               endif
+               tpq(:,:,:) = tpq(:,:,:) + zdqmufi(1:ngrid,1:nlayer,1:nq)*ptimestep
+#else
             ! Inside this routine we will split 2D->1D, intensive->extensive and separate different types of tracers
             ! Should be put in phytrac
@@ -1119,5 +1158,5 @@
 
             pdq(1:ngrid,1:nlayer,1:nq) = pdq(1:ngrid,1:nlayer,1:nq) + zdqmufi(1:ngrid,1:nlayer,1:nq)   
-
+#endif
          endif ! end of 'callmufi'
       
@@ -1340,66 +1379,4 @@
 
 
-!-----------------------------------
-!        Saving statistics :
-!-----------------------------------
-
-!    Note :("stats" stores and accumulates 8 key variables in file "stats.nc"
-!           which can later be used to make the statistic files of the run:
-!           "stats")          only possible in 3D runs !!!
-
-         
-      if (callstats) then
-
-         call wstats(ngrid,"ps","Surface pressure","Pa",2,ps)
-         call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf)
-         call wstats(ngrid,"fluxsurf_lw",                               &
-                     "Thermal IR radiative flux to surface","W.m-2",2,  &
-                     fluxsurf_lw)
-         call wstats(ngrid,"fluxtop_lw",                                &
-                     "Thermal IR radiative flux to space","W.m-2",2,    &
-                     fluxtop_lw)
-                     
-!            call wstats(ngrid,"fluxsurf_sw",                               &
-!                        "Solar radiative flux to surface","W.m-2",2,       &
-!                         fluxsurf_sw_tot)                     
-!            call wstats(ngrid,"fluxtop_sw",                                &
-!                        "Solar radiative flux to space","W.m-2",2,         &
-!                        fluxtop_sw_tot)
-
-
-         call wstats(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
-         call wstats(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
-         call wstats(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
-         !call wstats(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
-         !call wstats(ngrid,"ALB_1st","First Band Surface albedo"," ",2,albedo(:,1))
-         call wstats(ngrid,"p","Pressure","Pa",3,pplay)
-         call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt)
-         call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu)
-         call wstats(ngrid,"v","Meridional (North-South) wind","m.s-1",3,zv)
-         call wstats(ngrid,"w","Vertical (down-up) wind","m.s-1",3,pw)
-         call wstats(ngrid,"q2","Boundary layer eddy kinetic energy","m2.s-2",3,q2)
-
-         if (tracer) then
-            do iq=1,nq
-               call wstats(ngrid,noms(iq),noms(iq),'kg/kg',3,zq(1,1,iq))
-               call wstats(ngrid,trim(noms(iq))//'_surf',trim(noms(iq))//'_surf',  & 
-                           'kg m^-2',2,qsurf(1,iq) )
-                          
-!              call wstats(ngrid,trim(noms(iq))//'_reff',                          & 
-!                          trim(noms(iq))//'_reff',                                   & 
-!                          'm',3,reffrad(1,1,iq))
-
-            end do
-
-         endif ! end of 'tracer' 
-
-         if(lastcall) then
-            write (*,*) "Writing stats..."
-            call mkstats(ierr)
-         endif
-         
-      endif ! end of 'callstats'
-
-
 !-----------------------------------------------------------------------------------------------------
 !           OUTPUT in netcdf file "DIAGFI.NC", containing any variable for diagnostic
@@ -1475,5 +1452,5 @@
         ! Temporary inclusions for winds diagnostics.
         call writediagfi(ngrid,"zdudif","Turbdiff tend. zon. wind","m s-2",3,zdudif)
-	call writediagfi(ngrid,"zdudyn","Dyn. tend. zon. wind","m s-2",3,zdudyn)
+        call writediagfi(ngrid,"zdudyn","Dyn. tend. zon. wind","m s-2",3,zdudyn)
 
         ! Temporary inclusions for heating diagnostics.
@@ -1491,8 +1468,21 @@
 
          if (callmufi) then ! For now we assume an given order for tracers !
+#ifdef USE_QTEST
+            ! Microphysical tracers passed through dyn+phys(except mufi)
+            call writediagfi(ngrid,"mu_m0as_dp","Dynphys only spherical mode 0th order moment",'kg/kg',3,zq(:,:,1))
+            call writediagfi(ngrid,"mu_m3as_dp","Dynphys only spherical mode 3rd order moment",'kg/kg',3,zq(:,:,2))
+            call writediagfi(ngrid,"mu_m0af_dp","Dynphys only fractal mode 0th order moment",'kg/kg',3,zq(:,:,3))
+            call writediagfi(ngrid,"mu_m3af_dp","Dynphys only fractal mode 3rd order moment",'kg/kg',3,zq(:,:,4))
+            ! Microphysical tracers passed through mufi only
+            call writediagfi(ngrid,"mu_m0as_mo","Mufi only spherical mode 0th order moment",'kg/kg',3,tpq(:,:,1))
+            call writediagfi(ngrid,"mu_m3as_mo","Mufi only spherical mode 3rd order moment",'kg/kg',3,tpq(:,:,2))
+            call writediagfi(ngrid,"mu_m0af_mo","Mufi only fractal mode 0th order moment",'kg/kg',3,tpq(:,:,3))
+            call writediagfi(ngrid,"mu_m3af_mo","Mufi only fractal mode 3rd order moment",'kg/kg',3,tpq(:,:,4))
+#else
             call writediagfi(ngrid,"mu_m0as","Spherical mode 0th order moment",'kg/kg',3,zq(:,:,1))
             call writediagfi(ngrid,"mu_m3as","Spherical mode 3rd order moment",'kg/kg',3,zq(:,:,2))
             call writediagfi(ngrid,"mu_m0af","Fractal mode 0th order moment",'kg/kg',3,zq(:,:,3))
             call writediagfi(ngrid,"mu_m3af","Fractal mode 3rd order moment",'kg/kg',3,zq(:,:,4))
+#endif
          endif ! end of 'callmufi'
 
Index: trunk/LMDZ.TITAN/libf/phytitan/setspi.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/setspi.F90	(revision 1896)
+++ trunk/LMDZ.TITAN/libf/phytitan/setspi.F90	(revision 1897)
@@ -138,4 +138,8 @@
       end do
 !     note M=L_NSPECTI+1 after loop due to Fortran bizarreness
+!        --> No fortran bizarreness here... incrementation is performed at the end of the loop...
+!        --> then when M reached L_NSPECTI, we initialiaze the last element of each array and
+!            ... increment one last time M... tadaaaa, mystery solved !
+!            The same logic is applied on for loop in C !
 
 !=======================================================================
