Index: trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90	(revision 2241)
+++ trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90	(revision 2242)
@@ -92,5 +92,5 @@
 
   ! Initialization of zdq here since intent=out and no action performed on every tracers
-  zdq(:,:,:) = 0.0
+  zdq(:,:,:) = 0.D0
 
   ! Initialize tracers updated with former processes from physics
@@ -102,8 +102,24 @@
     ! We suppose a given order of tracers !
     int2ext(ilon,:) = ( plev(ilon,1:nlay)-plev(ilon,2:nlay+1) ) / g3d(ilon,1:nlay)
-    m0as(:) = zq(ilon,:,1) * int2ext(ilon,:)
-    m3as(:) = zq(ilon,:,2) * int2ext(ilon,:)
-    m0af(:) = zq(ilon,:,3) * int2ext(ilon,:)
-    m3af(:) = zq(ilon,:,4) * int2ext(ilon,:)
+    
+    ! Check because of the threshold of small tracers values in the dynamics
+    ! WARNING : With this patch it enables to handles the small values required
+    ! by YAMMS, but still it might still leads to some unphysical values of
+    ! radii inside YAMMS, harmless for now, but who might be a problem the day
+    ! you'll want to compute optics from the radii.
+    WHERE (pq(ilon,:,1) > 2.D-200 .AND. pq(ilon,:,2) > 2.D-200) ! Test on both moments to avoid divergences if one hit the threshold but not the other
+      m0as(:) = zq(ilon,:,1) * int2ext(ilon,:)                  ! It can still be a pb if both m0 and m3 has been set to epsilon at the beginning of dynamics
+      m3as(:) = zq(ilon,:,2) * int2ext(ilon,:)                  ! then mixed, even though both are above the threshold, their ratio can be nonsense
+    ELSEWHERE
+      m0as(:)=0.D0
+      m3as(:)=0.D0
+    ENDWHERE
+    WHERE (pq(ilon,:,3) > 2.D-200 .AND. pq(ilon,:,4) > 2.D-200)
+      m0af(:) = zq(ilon,:,3) * int2ext(ilon,:)
+      m3af(:) = zq(ilon,:,4) * int2ext(ilon,:)
+    ELSEWHERE
+      m0af(:)=0.D0
+      m3af(:)=0.D0
+    ENDWHERE
     
     if (callclouds) then ! if call clouds
@@ -131,6 +147,6 @@
     !dm0n(:) = 0._mm_wp ; dm3n(:) = 0._mm_wp ; dm3i(:,:) = 0._mm_wp ; dgazs(:,:) = 0._mm_wp
 
-    dm0as(:) = 0.0 ; dm3as(:) = 0.0 ; dm0af(:) = 0.0 ; dm3af(:) = 0.0
-    dm0n(:) = 0.0 ; dm3n(:) = 0.0 ; dm3i(:,:) = 0.0 ; dgazs(:,:) = 0.0
+    dm0as(:) = 0.D0 ; dm3as(:) = 0.D0 ; dm0af(:) = 0.D0 ; dm3af(:) = 0.D0
+    dm0n(:) = 0.D0 ; dm3n(:) = 0.D0 ; dm3i(:,:) = 0.D0 ; dgazs(:,:) = 0.D0
     ! call microphysics
 
@@ -165,8 +181,5 @@
       enddo
     endif
-
-    ! Sanity check ( way safer to be done here rather than within YAMMS )
-    WHERE (zq+zdq < 0.0) ; zdq = -zq ; END WHERE
-
+    
   END DO ! loop on ilon
 
Index: trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90	(revision 2241)
+++ trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90	(revision 2242)
@@ -15,5 +15,5 @@
       
       ! Default directories for correlated-k data
-      ! Set in inifis_mod
+      ! Read in inifis_mod, set in physiq_mod
       character(LEN=100),save :: corrkdir = 'datagcm/corrk_data/50X50'
       character(LEN=100),save :: banddir  = 'datagcm/corrk_data/50X50/50x50'
@@ -25,4 +25,9 @@
 !$OMP THREADPRIVATE(config_mufi)
 
+      ! Default file for coupled microphysics optical properties
+      ! Set in physiq_mod
+      character(LEN=100),save :: haze_opt_file ='datagcm/HAZE_OPT_23x23.DAT'
+!$OMP THREADPRIVATE(haze_opt_file)
+
       end module datafile_mod
 !-----------------------------------------------------------------------
Index: trunk/LMDZ.TITAN/libf/phytitan/optci.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 2241)
+++ trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 2242)
@@ -6,9 +6,9 @@
                          tgasref,pfgasref,wnoi,scalep,indi
   use gases_h
+  use datafile_mod, only: haze_opt_file
   use comcstfi_mod, only: r
   use callkeys_mod, only: continuum,graybody,corrk_recombin,               &
                           callclouds,callmufi,seashaze,uncoupl_optic_haze
   use tracer_h, only : nmicro,nice
-  use MMP_OPTICS
 
   implicit none
@@ -83,5 +83,4 @@
   ! variables for k in units m^-1
   real*8 dz(L_LEVELS)
-  !real*8 rho !! see test below
 
   integer igas, jgas, ilay
@@ -89,9 +88,19 @@
   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
+  ! Variables for haze optics
+  character(len=200) file_path
+  logical file_ok
+  integer dumch
+  real*8  dumwvl
+
+  real*8 m3as,m3af
+  real*8 dtauaer_s,dtauaer_f
+  real*8,save :: rhoaer_s(L_NSPECTI),ssa_s(L_NSPECTI),asf_s(L_NSPECTI)
+  real*8,save :: rhoaer_f(L_NSPECTI),ssa_f(L_NSPECTI),asf_f(L_NSPECTI)
+!$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,asf_s,asf_f)
+
   logical,save :: firstcall=.true.
-  !$OMP THREADPRIVATE(firstcall)
+!$OMP THREADPRIVATE(firstcall)
+
 
   !! AS: to save time in computing continuum (see bilinearbig)
@@ -101,10 +110,20 @@
   ENDIF
 
-  ! Some initialisation beacause there's a pb with disr_haze at the limits (nw=1)
+  ! Some initialisation because there's a pb with disr_haze at the limits (nw=1)
   ! I should check this - For now we set vars to zero : better than nans - JVO 2017
-  
-  dhaze_t(:,:) = 0.
-  ssa_t(:,:) = 0.
-  asf_t(:,:) = 0.
+  DHAZE_T(:,:) = 0.0
+  SSA_T(:,:)   = 0.0
+  ASF_T(:,:)   = 0.0
+
+  ! Load tabulated haze optical properties if needed.
+  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  IF (firstcall .AND. callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
+     OPEN(12,file=TRIM(haze_opt_file),form='formatted') ! The file has been inquired in physiq_mod firstcall
+     READ(12,*) ! dummy header
+     DO NW=1,L_NSPECTI
+       READ(12,*) dumch, dumwvl, rhoaer_f(nw), ssa_f(nw), asf_f(nw), rhoaer_s(nw), ssa_s(nw), asf_s(nw)
+     ENDDO
+     CLOSE(12)
+  ENDIF
 
   !=======================================================================
@@ -118,5 +137,5 @@
   do K=2,L_LEVELS
   
-     ilay = k / 2 ! int. arithmetic => gives the gcm layer index
+     ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
      
      DPR(k) = PLEV(K)-PLEV(K-1)
@@ -137,21 +156,29 @@
      do K=2,L_LEVELS
      
-        ilay = k / 2 ! int. arithmetic => gives the gcm layer index
+        ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
         
-        ! 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 = pqmo(ilay,1)
-          m3as = pqmo(ilay,2)
-          m0af = pqmo(ilay,3)
-          m3af = pqmo(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)
+          m3as = pqmo(ilay,2) / 2.0
+          m3af = pqmo(ilay,4) / 2.0
+          
+          IF ( ilay .lt. 18 ) THEN
+            m3as = pqmo(18,2) / 2.0 *dz(k)/dz(18)
+            m3af = pqmo(18,4) / 2.0 *dz(k)/dz(18)
+          ENDIF
+
+          dtauaer_s     = m3as*rhoaer_s(nw)
+          dtauaer_f     = m3af*rhoaer_f(nw)
+          DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
+          
+          IF ( dtauaer_s + dtauaer_f .GT. 1.D-30 ) THEN
+            SSA_T(k,nw)   = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
+            ASF_T(k,nw)   = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) )  &
+                            / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f )
+          ELSE
+             DHAZE_T(k,nw) = 0.D0
+             SSA_T(k,nw)   = 1.0
+             ASF_T(k,nw)   = 1.0
+          ENDIF
+          
           IF (callclouds.and.firstcall) & 
             WRITE(*,*) 'WARNING: In optci, optical properties &
@@ -159,6 +186,6 @@
         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))
-          if (seashaze) dhaze_T(k,nw) = dhaze_T(k,nw)*seashazefact(k)
+          call disr_haze(dz(k),plev(k),wnoi(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
+          if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
         ENDIF
 
Index: trunk/LMDZ.TITAN/libf/phytitan/optcv.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 2241)
+++ trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 2242)
@@ -6,4 +6,5 @@
                          tgasref,pfgasref,wnov,scalep,indv
   use gases_h
+  use datafile_mod, only: haze_opt_file
   use comcstfi_mod, only: r
   use callkeys_mod, only: continuum,graybody,callgasvis,corrk_recombin,     &
@@ -99,9 +100,18 @@
   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
+  ! Variables for haze optics
+  character(len=200) file_path
+  logical file_ok
+  integer dumch
+  real*8  dumwvl
+
+  real*8 m3as,m3af
+  real*8 dtauaer_s,dtauaer_f
+  real*8,save :: rhoaer_s(L_NSPECTV),ssa_s(L_NSPECTV),asf_s(L_NSPECTV)
+  real*8,save :: rhoaer_f(L_NSPECTV),ssa_f(L_NSPECTV),asf_f(L_NSPECTV)
+!$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,asf_s,asf_f)
+  
   logical,save :: firstcall=.true.
-  !$OMP THREADPRIVATE(firstcall)
+!$OMP THREADPRIVATE(firstcall)
 
 
@@ -114,9 +124,21 @@
   ! Some initialisation beacause there's a pb with disr_haze at the limits (nw=1)
   ! I should check this - For now we set vars to zero : better than nans - JVO 2017
-  
-  dhaze_t(:,:) = 0.
-  ssa_t(:,:) = 0.
-  asf_t(:,:) = 0.
-
+  DHAZE_T(:,:) = 0.0
+  SSA_T(:,:)   = 0.0
+  ASF_T(:,:)   = 0.0
+  
+  ! Load tabulated haze optical properties if needed.
+  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  IF (firstcall .AND. callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
+     OPEN(12,file=TRIM(haze_opt_file),form='formatted') ! The file has been inquired in physiq_mod firstcall
+     READ(12,*) ! dummy header
+     DO NW=1,L_NSPECTI
+       READ(12,*) ! there's IR 1st
+     ENDDO
+     DO NW=1,L_NSPECTV
+       READ(12,*) dumch, dumwvl, rhoaer_f(nw), ssa_f(nw), asf_f(nw), rhoaer_s(nw), ssa_s(nw), asf_s(nw)
+     ENDDO
+     CLOSE(12)
+  ENDIF
 
   !=======================================================================
@@ -132,5 +154,5 @@
   do K=2,L_LEVELS
   
-     ilay = k / 2 ! int. arithmetic => gives the gcm layer index
+     ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
   
      DPR(k) = PLEV(K)-PLEV(K-1)
@@ -150,5 +172,5 @@
   ! Rayleigh scattering
   do NW=1,L_NSPECTV
-     TRAY(1:4,NW)   = 1d-30
+     TRAY(1:4,NW)   = 1.d-30
      do K=5,L_LEVELS
         TRAY(K,NW)   = TAURAY(NW) * DPR(K)
@@ -159,23 +181,31 @@
   do K=2,L_LEVELS
   
-     ilay = k / 2 ! int. arithmetic => gives the gcm layer index
+     ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
 
      do NW=1,L_NSPECTV
      
-        ! 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 = pqmo(ilay,1)
-          m3as = pqmo(ilay,2)
-          m0af = pqmo(ilay,3)
-          m3af = pqmo(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)
+          m3as = pqmo(ilay,2) / 2.0
+          m3af = pqmo(ilay,4) / 2.0
+          
+          IF ( ilay .lt. 18 ) THEN
+            m3as = pqmo(18,2) / 2.0 * dz(k) / dz(18)
+            m3af = pqmo(18,4) / 2.0 * dz(k) / dz(18)
+          ENDIF
+
+          dtauaer_s     = m3as*rhoaer_s(nw)
+          dtauaer_f     = m3af*rhoaer_f(nw)
+          DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
+
+          IF ( dtauaer_s + dtauaer_f .GT. 1.D-30 ) THEN
+            SSA_T(k,nw)   = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
+            ASF_T(k,nw)   = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) )  &
+                            / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f )
+          ELSE
+             DHAZE_T(k,nw) = 0.D0
+             SSA_T(k,nw)   = 1.0
+             ASF_T(k,nw)   = 1.0
+          ENDIF
+          
           IF (callclouds.and.firstcall) & 
             WRITE(*,*) 'WARNING: In optcv, optical properties &
@@ -183,14 +213,14 @@
         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))
-          if (seashaze) dhaze_T(k,nw) = dhaze_T(k,nw)*seashazefact(k)
+          call disr_haze(dz(k),plev(k),wnov(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
+          if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
         ENDIF
         
         !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
-	!   but visible does not handle very well diffusion in first layer.
-	!   The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
-	!   in the 4 first semilayers in optcv, but not optci.
-	!   This solves random variations of the sw heating at the model top. 
-        if (k<5)  dhaze_T(K,:) = 0.0
+        !   but visible does not handle very well diffusion in first layer.
+        !   The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
+        !   in the 4 first semilayers in optcv, but not optci.
+        !   This solves random variations of the sw heating at the model top. 
+        if (k<5)  DHAZE_T(K,:) = 0.0
          
         DRAYAER = TRAY(K,NW)
@@ -316,9 +346,9 @@
   ! Haze scattering
             !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
-	    !   but not in the visible
-	    !   The dhaze_s is thus set to 0 in the 4 first semilayers in optcv, but not optci.
-	    !   This solves random variations of the sw heating at the model top. 
+            !   but not in the visible
+            !   The dhaze_s is thus set to 0 in the 4 first semilayers in optcv, but not optci.
+            !   This solves random variations of the sw heating at the model top. 
   DO NW=1,L_NSPECTV
-      DHAZES_T(1:4,NW) = 0.d0
+    DHAZES_T(1:4,NW) = 0.d0
     DO K=5,L_LEVELS
       DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW) ! effect of scattering albedo on haze
@@ -386,5 +416,4 @@
   END DO                 ! end full gauss loop
 
-
   if(firstcall) firstcall = .false.
 
Index: trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 2241)
+++ trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 2242)
@@ -21,5 +21,5 @@
       use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
       use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
-      use datafile_mod, only: datadir, corrkdir, banddir
+      use datafile_mod, only: datadir, corrkdir, banddir, haze_opt_file
       use geometry_mod, only: latitude, longitude, cell_area
       USE comgeomfi_h, only: totarea, totarea_planet
@@ -402,4 +402,5 @@
       real,dimension(ngrid,nlayer,nq) :: dtpq ! (temporary in 01/18)
 
+      logical file_ok
 
 !-----------------------------------------------------------------------------
@@ -417,7 +418,7 @@
         REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: g3d   !! Latitude-Altitude depending gravitational acceleration (m.s-2).
         REAL(kind=8), DIMENSION(:,:), INTENT(IN) :: temp  !! Temperature at the center of each layer (K).  
-        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq    !! Tracers (\(kg.kg^{-1}}\)).
-        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: zdqfi !! Tendency from former processes for tracers (\(kg.kg^{-1}}\)).
-        REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq   !! Microphysical tendency for tracers (\(kg.kg^{-1}}\)).
+        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: pq    !! Tracers (\(X.kg^{-1}}\)).
+        REAL(kind=8), DIMENSION(:,:,:), INTENT(IN)  :: zdqfi !! Tendency from former processes for tracers (\(X.kg^{-1}}\)).
+        REAL(kind=8), DIMENSION(:,:,:), INTENT(OUT) :: zdq   !! Microphysical tendency for tracers (\(X.kg^{-1}}\)).
       END SUBROUTINE calmufi
     END INTERFACE
@@ -484,12 +485,12 @@
 !        Variables set to 0
 !        ~~~~~~~~~~~~~~~~~~
-         dtrad(:,:) = 0.0
-         fluxrad(:) = 0.0
-         zdtsw(:,:) = 0.0
-         zdtlw(:,:) = 0.0
+         dtrad(:,:) = 0.D0
+         fluxrad(:) = 0.D0
+         zdtsw(:,:) = 0.D0
+         zdtlw(:,:) = 0.D0
 
 !        Initialize setup for correlated-k radiative transfer
-!        JVO 17 : Was in callcorrk firstcall, but we need spectral intervals for microphysics
-!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!        JVO 17 : Was in callcorrk firstcall, but we need spectral intervals for microphysics.
+!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
          if (corrk) then
@@ -512,4 +513,16 @@
            int_dtauv(:,:,:) = 0.D0
            
+           IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
+             haze_opt_file=trim(datadir)//'/HAZE_OPTIC_'//trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))//'.DAT'
+             inquire(file=trim(haze_opt_file),exist=file_ok)
+             if(.not.file_ok) then
+               write(*,*) 'The file ',TRIM(haze_opt_file),' with the haze optical properties'
+               write(*,*) 'was not found by optci.F90 ! Check in ', TRIM(datadir)
+               write(*,*) 'that you have the one corresponding to the given spectral resolution !!'
+               write(*,*) 'Meanwhile I abort ...'
+               call abort
+              endif
+           ENDIF
+           
          endif
 
@@ -536,10 +549,6 @@
 
          IF ( callmufi ) THEN
-           ! WARNING JB (27/03/2018): inimufi AND mmp_initialize_optics should be wrapped in an OMP SINGLE statement.
+           ! WARNING JB (27/03/2018): inimufi should be wrapped in an OMP SINGLE statement.
            call inimufi(ptimestep)
-           ! Optical coupling of YAMMS is plugged but inactivated for now
-           ! as long as the microphysics only isn't fully debugged -- JVO 01/18
-           ! NOTE JB (03/18): mmp_optic_file is initialized in inimufi => mmp_initialize_optics must be called AFTER inimufi
-           IF (.NOT.uncoupl_optic_haze) call mmp_initialize_optics(mmp_optic_file)
 
            ! initialize microphysics diagnostics arrays.
@@ -547,5 +556,4 @@
 
          ENDIF
-
 
 #ifdef CPP_XIOS
@@ -582,6 +590,6 @@
 !        Initialize albedo calculation. 
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-         albedo(:,:)=0.0
-          albedo_bareground(:)=0.0
+         albedo(:,:)=0.D0
+          albedo_bareground(:)=0.D0
          call surfini(ngrid,nq,qsurf,albedo,albedo_bareground) 
          
@@ -660,14 +668,14 @@
       ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
-      pdt(:,:)    = 0.0     
-      zdtsurf(:)  = 0.0
-      pdq(:,:,:)  = 0.0
-      dqsurf(:,:) = 0.0
-      pdu(:,:)    = 0.0
-      pdv(:,:)    = 0.0
-      pdpsrf(:)   = 0.0
-      zflubid(:)  = 0.0      
-      taux(:)     = 0.0
-      tauy(:)     = 0.0
+      pdt(:,:)    = 0.D0     
+      zdtsurf(:)  = 0.D0
+      pdq(:,:,:)  = 0.D0
+      dqsurf(:,:) = 0.D0
+      pdu(:,:)    = 0.D0
+      pdv(:,:)    = 0.D0
+      pdpsrf(:)   = 0.D0
+      zflubid(:)  = 0.D0      
+      taux(:)     = 0.D0
+      tauy(:)     = 0.D0
 
       zday=pday+ptime ! Compute time, in sols (and fraction thereof).
@@ -906,5 +914,5 @@
                fluxtop_lw(:)     = emis(:)*sigma*tsurf(:)**4
 
-               dtrad(:,:)=0.0 ! no atmospheric radiative heating
+               dtrad(:,:)=0.D0 ! no atmospheric radiative heating
 
             endif ! end of corrk
@@ -1044,7 +1052,7 @@
 
          zdh(:,:) = pdt(:,:)/zpopsk(:,:)
-         zduadj(:,:)=0.0
-         zdvadj(:,:)=0.0
-         zdhadj(:,:)=0.0
+         zduadj(:,:)=0.D0
+         zdvadj(:,:)=0.D0
+         zdhadj(:,:)=0.D0
 
 
@@ -1088,12 +1096,32 @@
  
          if (callmufi) then
-            zzlev(:,nlayer+1)=zzlay(:,nlayer)+(zzlay(:,nlayer)-zzlev(:,nlayer)) ! JVO 19 : We assume zzlev isn't reused later on
+
+            zzlev(:,nlayer+1)=zzlay(:,nlayer)+(zzlay(:,nlayer)-zzlev(:,nlayer)) ! JVO 19 : We assume zzlev isn't reused later on (could be done cleaner)
+
 #ifdef USE_QTEST
-            dtpq(:,:,:) = 0.0 ! we want tpq to go only through mufi
+            dtpq(:,:,:) = 0.D0 ! we want tpq to go only through mufi
             call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,tpq,dtpq,zdqmufi)
             tpq(:,:,:) = tpq(:,:,:) + zdqmufi(:,:,:)*ptimestep ! only manipulation of tpq->*ptimestep here
 #else
+            
             call calmufi(ptimestep,pplev,zzlev,pplay,zzlay,gzlat,pt,pq,pdq,zdqmufi)
+
             pdq(:,:,:) = pdq(:,:,:) + zdqmufi(:,:,:)
+
+            ! Sanity check ( way safer to be done here rather than within YAMMS )
+            ! Important : the sanity check intentionally include the former processes tendency !
+            ! NB : Despite this sanity check there might be still some unphysical values going through :
+            !        - Negatives, but harmless as it will be only for the output files
+            !          just remove them in post-proc.
+            !        - Weird unphysical ratio of m0 and m3, ok for now, but take care of them if
+            !          you want to compute optics from radii.
+            WHERE ( (pq(:,:,1)+pdq(:,:,1)*ptimestep < 0.D0) .OR. (pq(:,:,2)+pdq(:,:,2)*ptimestep < 0.D0) )
+                    pdq(:,:,1) = (epsilon(1.0)-1.D0)*pq(:,:,1)/ptimestep 
+                    pdq(:,:,2) = (epsilon(1.0)-1.D0)*pq(:,:,2)/ptimestep
+            ENDWHERE 
+            WHERE ( (pq(:,:,3)+pdq(:,:,3)*ptimestep < 0.D0) .OR. (pq(:,:,4)+pdq(:,:,4)*ptimestep < 0.D0) )
+                    pdq(:,:,3) = (epsilon(1.0)-1.D0)*pq(:,:,3)/ptimestep 
+                    pdq(:,:,4) = (epsilon(1.0)-1.D0)*pq(:,:,4)/ptimestep
+            ENDWHERE
 #endif
 
@@ -1103,4 +1131,5 @@
               call calmufi(ptimestep,zplevbar,zzlevbar,zplaybar,zzlaybar, &
                            gzlat,ztfibar,zqfibar,zdqfibar,zdqmufibar)
+              ! TODO : Add a sanity check here !
             endif
 
@@ -1215,5 +1244,5 @@
             ! v. Updates and positivity check
             ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-            zdqchi(:,:,:)  = 0.D0 ! -> for the microphysical tracers we need to update to 0 each time ;-)
+            zdqchi(:,:,:)  = 0.D0 ! -> dycchi is saved but for the nmicro tracers we must update to 0 at each step
 
             do iq=1,nkim
@@ -1315,6 +1344,6 @@
 
       if(firstcall)then
-         zdtdyn(:,:)=0.0
-         zdudyn(:,:)=0.0
+         zdtdyn(:,:)=0.D0
+         zdudyn(:,:)=0.D0
       endif
 
@@ -1674,15 +1703,15 @@
         CALL send_xios_field("dqcond_C3H6",dyccond(:,:,17+nmicro))
         CALL send_xios_field("dqcond_C4H4",dyccond(:,:,21+nmicro))
-        CALL send_xios_field("dqcond_CH3CCH",dyccond(:,:,24+nmicro))
-        CALL send_xios_field("dqcond_C3H8",dyccond(:,:,25+nmicro))
-        CALL send_xios_field("dqcond_C4H2",dyccond(:,:,26+nmicro))
-        CALL send_xios_field("dqcond_C4H6",dyccond(:,:,27+nmicro))
-        CALL send_xios_field("dqcond_C4H10",dyccond(:,:,28+nmicro))
-        CALL send_xios_field("dqcond_AC6H6",dyccond(:,:,29+nmicro))
-        CALL send_xios_field("dqcond_HCN",dyccond(:,:,36+nmicro))
-        CALL send_xios_field("dqcond_CH3CN",dyccond(:,:,40+nmicro))
-        CALL send_xios_field("dqcond_HC3N",dyccond(:,:,42+nmicro))
-        CALL send_xios_field("dqcond_NCCN",dyccond(:,:,43+nmicro))
-        CALL send_xios_field("dqcond_C4N2",dyccond(:,:,44+nmicro))
+        CALL send_xios_field("dqcond_CH3CCH",dyccond(:,:,23+nmicro))
+        CALL send_xios_field("dqcond_C3H8",dyccond(:,:,24+nmicro))
+        CALL send_xios_field("dqcond_C4H2",dyccond(:,:,25+nmicro))
+        CALL send_xios_field("dqcond_C4H6",dyccond(:,:,26+nmicro))
+        CALL send_xios_field("dqcond_C4H10",dyccond(:,:,27+nmicro))
+        CALL send_xios_field("dqcond_AC6H6",dyccond(:,:,28+nmicro))
+        CALL send_xios_field("dqcond_HCN",dyccond(:,:,35+nmicro))
+        CALL send_xios_field("dqcond_CH3CN",dyccond(:,:,39+nmicro))
+        CALL send_xios_field("dqcond_HC3N",dyccond(:,:,41+nmicro))
+        CALL send_xios_field("dqcond_NCCN",dyccond(:,:,42+nmicro))
+        CALL send_xios_field("dqcond_C4N2",dyccond(:,:,43+nmicro))
 
         ! Pseudo-evaporation flux (mol/mol/s)
