Index: trunk/LMDZ.TITAN/libf/phytitan/optci.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 2133)
+++ 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
 
