Index: trunk/LMDZ.TITAN/libf/phytitan/bilinearbig.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/bilinearbig.F90	(revision 1681)
+++ trunk/LMDZ.TITAN/libf/phytitan/bilinearbig.F90	(revision 1722)
@@ -56,11 +56,18 @@
 !     ... and for y within the temperature range
       if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
-         write(*,*) 'Warning from bilinearbig routine:'
-         write(*,*) 'Outside continuum temperature range!'
+         print*,y_arr(1),y_arr(nY)
+         !write(*,*) 'Warning from bilinearbig routine:'
+         !write(*,*) 'Outside continuum temperature range!'
          if(y.lt.y_arr(1))then
             y=y_arr(1)+0.01
+            b=1
+            y1=y_arr(b)
+            y2=y_arr(b+1)
          endif
          if(y.gt.y_arr(nY))then
             y=y_arr(nY)-0.01
+            b=nY-1
+            y1=y_arr(b)
+            y2=y_arr(b+1)
          endif
       else
Index: trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 1681)
+++ trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 1722)
@@ -8,4 +8,5 @@
           tau_col,firstcall,lastcall)
 
+      use mod_phys_lmdz_para, only : is_master
       use radinc_h
       use radcommon_h
@@ -17,5 +18,5 @@
       USE tracer_h
       use comcstfi_mod, only: pi, mugaz, cpp
-      use callkeys_mod, only: diurnal,tracer,nosurf,        &
+      use callkeys_mod, only: diurnal,tracer,        &
                               strictboundcorrk,specOLR
 
@@ -75,5 +76,5 @@
       REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)        ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
       REAL,INTENT(OUT) :: tau_col(ngrid)                 ! Diagnostic from aeropacity.
-      REAL,INTENT(OUT) :: albedo_equivalent(ngrid)       ! Spectrally Integrated Albedo. For Diagnostic. By MT2015  
+      REAL,INTENT(OUT) :: albedo_equivalent(ngrid)       ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
       
       
@@ -113,5 +114,5 @@
       REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
 
-      REAL*8 tauaero(L_LEVELS+1,naerkind)
+      REAL*8 tauaero(L_LEVELS,naerkind)
       REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
       REAL*8 nfluxoutv_nu(L_NSPECTV)                 ! Outgoing band-resolved VI flux at TOA (W/m2).
@@ -167,10 +168,10 @@
 
         ! test on allocated necessary because of CLFvarying (two calls to callcorrk in physiq)
-        if(.not.allocated(QXVAER)) allocate(QXVAER(L_LEVELS+1,L_NSPECTV,naerkind))
-        if(.not.allocated(QSVAER)) allocate(QSVAER(L_LEVELS+1,L_NSPECTV,naerkind))
-        if(.not.allocated(GVAER)) allocate(GVAER(L_LEVELS+1,L_NSPECTV,naerkind))
-        if(.not.allocated(QXIAER)) allocate(QXIAER(L_LEVELS+1,L_NSPECTI,naerkind))
-        if(.not.allocated(QSIAER)) allocate(QSIAER(L_LEVELS+1,L_NSPECTI,naerkind))
-        if(.not.allocated(GIAER)) allocate(GIAER(L_LEVELS+1,L_NSPECTI,naerkind))
+        if(.not.allocated(QXVAER)) allocate(QXVAER(L_LEVELS,L_NSPECTV,naerkind))
+        if(.not.allocated(QSVAER)) allocate(QSVAER(L_LEVELS,L_NSPECTV,naerkind))
+        if(.not.allocated(GVAER)) allocate(GVAER(L_LEVELS,L_NSPECTV,naerkind))
+        if(.not.allocated(QXIAER)) allocate(QXIAER(L_LEVELS,L_NSPECTI,naerkind))
+        if(.not.allocated(QSIAER)) allocate(QSIAER(L_LEVELS,L_NSPECTI,naerkind))
+        if(.not.allocated(GIAER)) allocate(GIAER(L_LEVELS,L_NSPECTI,naerkind))
 
          !!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...)
@@ -376,5 +377,5 @@
             ! Test / Correct for freaky s. s. albedo values.
             do iaer=1,naerkind
-               do k=1,L_LEVELS+1
+               do k=1,L_LEVELS
 
                   do nw=1,L_NSPECTV
@@ -419,7 +420,5 @@
             ! boundary conditions
             tauaero(1,iaer)          = tauaero(2,iaer)
-            tauaero(L_LEVELS+1,iaer) = tauaero(L_LEVELS,iaer)
             !tauaero(1,iaer)          = 0.
-            !tauaero(L_LEVELS+1,iaer) = 0.
             
          end do ! naerkind
@@ -430,14 +429,4 @@
             albv(nw)=albedo(ig,nw)
          ENDDO
-
-         if (nosurf) then ! Case with no surface.
-            DO nw=1,L_NSPECTV
-               if(albv(nw).gt.0.0) then
-                  print*,'For open lower boundary in callcorrk must'
-                  print*,'have spectral surface band albedos all set to zero!'
-                  call abort
-               endif
-            ENDDO         
-         endif
 
       if ((ngrid.eq.1).and.(global1d)) then ! Fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight.
Index: trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 1681)
+++ trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90	(revision 1722)
@@ -71,4 +71,7 @@
       real,save :: MassPlanet
 !$OMP THREADPRIVATE(flatten,Rmean,J2,MassPlanet)
+      real,save :: surfalbedo
+      real,save :: surfemis
+!$OMP THREADPRIVATE(surfalbedo,surfemis)
 
       logical,save :: iscallphys=.false.!existence of callphys.def
Index: trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 1681)
+++ trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90	(revision 1722)
@@ -10,5 +10,4 @@
 
   use radinc_h, only: ini_radinc_h, naerkind
-  use radcommon_h, only: ini_radcommon_h
   use datafile_mod, only: datadir
   use comdiurn_h, only: sinlat, coslat, sinlon, coslon
@@ -543,9 +542,7 @@
   ENDDO
 
+  ! initialize variables in radinc_h
   call ini_radinc_h(nlayer)
  
-  ! allocate "radcommon_h" arrays
-  call ini_radcommon_h()
-
   ! allocate "comsoil_h" arrays
   call ini_comsoil_h(ngrid)
Index: trunk/LMDZ.TITAN/libf/phytitan/optci.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 1681)
+++ trunk/LMDZ.TITAN/libf/phytitan/optci.F90	(revision 1722)
@@ -32,5 +32,5 @@
 
   real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
-  real*8 DTAUKI(L_LEVELS+1,L_NSPECTI,L_NGAUSS)
+  real*8 DTAUKI(L_LEVELS,L_NSPECTI,L_NGAUSS)
   real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
   real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
@@ -42,9 +42,9 @@
 
   ! for aerosols
-  real*8  QXIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
-  real*8  QSIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
-  real*8  GIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
-  real*8  TAUAERO(L_LEVELS+1,NAERKIND)
-  real*8  TAUAEROLK(L_LEVELS+1,L_NSPECTI,NAERKIND)
+  real*8  QXIAER(L_LEVELS,L_NSPECTI,NAERKIND)
+  real*8  QSIAER(L_LEVELS,L_NSPECTI,NAERKIND)
+  real*8  GIAER(L_LEVELS,L_NSPECTI,NAERKIND)
+  real*8  TAUAERO(L_LEVELS,NAERKIND)
+  real*8  TAUAEROLK(L_LEVELS,L_NSPECTI,NAERKIND)
   real*8  TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
 
@@ -52,8 +52,8 @@
   ! J. Vatant d'Ollone (2016)
   real*8 GWEIGHT(L_NGAUSS)
-  real*8 DHAZE_T(L_LEVELS+1,L_NSPECTI)
-  real*8 DHAZES_T(L_LEVELS+1,L_NSPECTI)
-  real*8 SSA_T(L_LEVELS+1,L_NSPECTI)
-  real*8 ASF_T(L_LEVELS+1,L_NSPECTI)
+  real*8 DHAZE_T(L_LEVELS,L_NSPECTI)
+  real*8 DHAZES_T(L_LEVELS,L_NSPECTI)
+  real*8 SSA_T(L_LEVELS,L_NSPECTI)
+  real*8 ASF_T(L_LEVELS,L_NSPECTI)
   real*8 INT_DTAU(L_NLAYRAD,L_NSPECTI)
   real*8 K_HAZE(L_NLAYRAD,L_NSPECTI)
@@ -115,5 +115,5 @@
   end do                    ! levels
 
-
+  ! Spectral dependance of aerosol absorption
   do iaer=1,naerkind
      DO NW=1,L_NSPECTI
@@ -127,9 +127,14 @@
 
      do K=2,L_LEVELS
-
+     
         ilay = k / 2 ! int. arithmetic => gives the gcm layer index
         
-! continuum absorption
-        DCONT = 0.0d0 
+        DAERO=SUM(TAEROS(K,NW,1:naerkind)) ! aerosol absorption
+
+        !================= Titan customisation ========================================
+        call disr_haze(dz(k),plev(k),wnoi(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
+        ! =============================================================================
+
+        DCONT = 0.0d0 ! continuum absorption
 
         if(continuum.and.(.not.graybody))then
@@ -207,12 +212,4 @@
         endif
 
-! aerosol absorption
-	DAERO=SUM(TAEROS(K,NW,1:naerkind))
-
-!================= Titan customisation ========================================
-        call disr_haze(dz(k),plev(k),wnoi(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
-! =============================================================================
-
-
         do ng=1,L_NGAUSS-1
 
@@ -251,7 +248,4 @@
   end do
 
-  DTAUKI(L_LEVELS+1,1:L_NSPECTI,1:L_NGAUSS)=0.d0
-  
-  
   !=======================================================================
   !     Now the full treatment for the layers, where besides the opacity
@@ -261,6 +255,6 @@
   do iaer=1,naerkind
     DO NW=1,L_NSPECTI
-     DO K=2,L_LEVELS+1
-           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
+     DO K=2,L_LEVELS
+           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER) ! effect of scattering albedo
      ENDDO
     ENDDO
@@ -269,5 +263,5 @@
   ! Haze scattering
   DO NW=1,L_NSPECTI
-    DO K=2,L_LEVELS+1
+    DO K=2,L_LEVELS
       DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
     ENDDO
@@ -281,13 +275,11 @@
      END DO ! L vertical loop
      
-     !last level
-     L              = L_NLAYRAD
-     K              = 2*L+1
-     
+     ! Last level
+     L           = L_NLAYRAD
+     K           = 2*L+1
      btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + DHAZES_T(K,NW)
      
   END DO                    ! NW spectral loop
   
-! ======================================================================
 
   DO NW=1,L_NSPECTI
@@ -296,5 +288,4 @@
 
         K              = 2*L+1
-
         DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
 
@@ -306,5 +297,8 @@
                    GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
            end do
-           atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
+           atemp = atemp +                   &
+                ASF_T(K,NW)*DHAZES_T(K,NW) + &
+                ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
+
            WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
         else
@@ -321,33 +315,32 @@
      END DO ! L vertical loop
      
-     !     No vertical averaging on bottom layer
-
-     L = L_NLAYRAD
-     K = 2*L+1
-     DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)
+     ! Last level
      
-        atemp = 0.
-        if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
-           do iaer=1,naerkind
-              atemp = atemp +                                     &
-                   GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER)
-           end do
-           atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
-           WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
-        else
-           WBARI(L,nw,ng) = 0.0D0
-           DTAUI(L,NW,NG) = 1.0D-9
-        endif
-
-        if(btemp(L,nw) .GT. 0.0d0) then
-           cosbi(L,NW,NG) = atemp/btemp(L,nw)
-        else
-           cosbi(L,NW,NG) = 0.0D0
-        end if
+     L              = L_NLAYRAD
+     K              = 2*L+1
+     DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) ! + 1.e-50
+
+     atemp = 0.
+     if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
+        do iaer=1,naerkind
+           atemp = atemp + GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER)
+        end do
+        atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
+        WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
+     else
+        WBARI(L,nw,ng) = 0.0D0
+        DTAUI(L,NW,NG) = 1.0D-9
+     endif
+
+     if(btemp(L,nw) .GT. 0.0d0) then
+        cosbi(L,NW,NG) = atemp/btemp(L,nw)
+     else
+        cosbi(L,NW,NG) = 0.0D0
+     end if
+
 
      ! Now the other Gauss points, if needed.
 
      DO NG=1,L_NGAUSS-1
-     
         IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
 
@@ -357,5 +350,7 @@
 
               if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
-                WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
+
+                 WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
+
               else
                  WBARI(L,nw,ng) = 0.0D0
@@ -366,16 +361,19 @@
            END DO ! L vertical loop
            
-            !     No vertical averaging on bottom layer
-
-            L = L_NLAYRAD
-            K = 2*L+1
-            DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)
-            if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
+           ! Last level 
+           L              = L_NLAYRAD
+           K              = 2*L+1
+           DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)! + 1.e-50
+
+           if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
+
               WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
-            else
-               WBARI(L,nw,ng) = 0.0D0
-               DTAUI(L,NW,NG) = 1.0D-9
-            endif
-            cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
+
+           else
+              WBARI(L,nw,ng) = 0.0D0
+              DTAUI(L,NW,NG) = 1.0D-9
+           endif
+
+           cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
            
         END IF
Index: trunk/LMDZ.TITAN/libf/phytitan/optcv.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 1681)
+++ trunk/LMDZ.TITAN/libf/phytitan/optcv.F90	(revision 1722)
@@ -24,5 +24,5 @@
   !     
   !     THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL  
-  !     IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VISUAL
+  !     IT CALCULATES FOR EACH LAYER, FOR EACH SPECTRAL INTERVAL IN THE VISUAL
   !     LAYER: WBAR, DTAU, COSBAR
   !     LEVEL: TAU
@@ -39,5 +39,5 @@
 
   real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
-  real*8 DTAUKV(L_LEVELS+1,L_NSPECTV,L_NGAUSS)
+  real*8 DTAUKV(L_LEVELS,L_NSPECTV,L_NGAUSS)
   real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
   real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
@@ -48,9 +48,9 @@
 
   ! for aerosols
-  real*8  QXVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
-  real*8  QSVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
-  real*8  GVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
-  real*8  TAUAERO(L_LEVELS+1,NAERKIND)
-  real*8  TAUAEROLK(L_LEVELS+1,L_NSPECTV,NAERKIND)
+  real*8  QXVAER(L_LEVELS,L_NSPECTV,NAERKIND)
+  real*8  QSVAER(L_LEVELS,L_NSPECTV,NAERKIND)
+  real*8  GVAER(L_LEVELS,L_NSPECTV,NAERKIND)
+  real*8  TAUAERO(L_LEVELS,NAERKIND)
+  real*8  TAUAEROLK(L_LEVELS,L_NSPECTV,NAERKIND)
   real*8  TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
 
@@ -58,8 +58,8 @@
   ! J. Vatant d'Ollone (2016)
   real*8 GWEIGHT(L_NGAUSS)
-  real*8 DHAZE_T(L_LEVELS+1,L_NSPECTI)
-  real*8 DHAZES_T(L_LEVELS+1,L_NSPECTI)
-  real*8 SSA_T(L_LEVELS+1,L_NSPECTI)
-  real*8 ASF_T(L_LEVELS+1,L_NSPECTI)
+  real*8 DHAZE_T(L_LEVELS,L_NSPECTI)
+  real*8 DHAZES_T(L_LEVELS,L_NSPECTI)
+  real*8 SSA_T(L_LEVELS,L_NSPECTI)
+  real*8 ASF_T(L_LEVELS,L_NSPECTI)
   real*8 INT_DTAU(L_NLAYRAD,L_NSPECTI)
   real*8 K_HAZE(L_NLAYRAD,L_NSPECTI)
@@ -73,5 +73,4 @@
   real*8  TAURAY(L_NSPECTV)
   real*8  TRAY(L_LEVELS,L_NSPECTV)
-  real*8  TRAYAER
   real*8  DPR(L_LEVELS), U(L_LEVELS)
   real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
@@ -79,4 +78,5 @@
   real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
   real*8 DCONT,DAERO
+  real*8 DRAYAER
   double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
   double precision p_cross
@@ -127,5 +127,5 @@
   end do                    ! levels
 
-
+  ! Spectral dependance of aerosol absorption
   do iaer=1,naerkind
      do NW=1,L_NSPECTV
@@ -135,4 +135,6 @@
      end do
   end do
+  
+  ! Rayleigh scattering
   do NW=1,L_NSPECTV
      do K=2,L_LEVELS
@@ -140,5 +142,5 @@
      end do                    ! levels
   end do
-
+  
   !     we ignore K=1...
   do K=2,L_LEVELS
@@ -148,9 +150,15 @@
      do NW=1,L_NSPECTV
 
-        TRAYAER = TRAY(K,NW)
-        !     TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
+        !================= Titan customisation ========================================
+ 	call disr_haze(dz(k),plev(k),wnov(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
+        ! =============================================================================
+
+        DRAYAER = TRAY(K,NW)
+        !     DRAYAER is Tau RAYleigh scattering, plus AERosol opacity
         do iaer=1,naerkind
-           TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
+           DRAYAER = DRAYAER + TAEROS(K,NW,IAER)
         end do
+
+        DRAYAER = DRAYAER + DHAZE_T(K,NW) ! Titan's aerosol
 
         DCONT = 0.0 ! continuum absorption
@@ -221,8 +229,4 @@
         endif
 
-!================= Titan customisation ========================================
- 	call disr_haze(dz(k),plev(k),wnov(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
-! =============================================================================
-
         do ng=1,L_NGAUSS-1
 
@@ -239,11 +243,11 @@
                 LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
 
+
            TAUGAS  = U(k)*ANS
 
            TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
            DTAUKV(K,nw,ng) = TAUGAS & 
-                             + TRAYAER & ! TRAYAER includes all scattering contributions
-                             + DCONT    & ! For parameterized continuum aborption
-			     + DHAZE_T(K,NW)  ! For Titan haze
+                             + DRAYAER & ! DRAYAER includes all scattering contributions
+                             + DCONT ! For parameterized continuum aborption
 
         end do
@@ -253,10 +257,5 @@
 
         NG              = L_NGAUSS
-        DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT  & ! For parameterized continuum absorption
-			      + DHAZE_T(K,NW)  ! For Titan haze
-
-        do iaer=1,naerkind
-           DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) +  TAEROS(K,NW,IAER)
-        end do ! a bug was here!
+        DTAUKV(K,nw,ng) = DRAYAER + DCONT ! Scattering + parameterized continuum absorption, including Titan's haze
 
      end do
@@ -267,10 +266,9 @@
   !     Now the full treatment for the layers, where besides the opacity
   !     we need to calculate the scattering albedo and asymmetry factors
-  ! ======================================================================
-  
+
   do iaer=1,naerkind
     DO NW=1,L_NSPECTV
-      DO K=2,L_LEVELS   ! AS: shouldn't this be L_LEVELS+1 ? (see optci)
-           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
+      DO K=2,L_LEVELS
+           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER) ! effect of scattering albedo
       ENDDO
     ENDDO
@@ -279,6 +277,6 @@
   ! Haze scattering
   DO NW=1,L_NSPECTV
-    DO K=2,L_LEVELS+1
-      DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
+    DO K=2,L_LEVELS
+      DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW) ! effect of scattering albedo on haze
     ENDDO
   ENDDO
@@ -288,26 +286,21 @@
      DO L=1,L_NLAYRAD-1
         K              = 2*L+1
-        
-	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind)) &
-                    + SUM(GVAER(K+1,NW,1:naerkind) * TAUAEROLK(K+1,NW,1:naerkind)) &
-                    + ASF_T(K,NW)*DHAZES_T(K,NW)
-                    
+	atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))+SUM(GVAER(K+1,NW,1:naerkind) * TAUAEROLK(K+1,NW,1:naerkind)) &
+                    + ASF_T(K,NW)*DHAZES_T(K,NW) + ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
         btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind)) &
-                    + DHAZES_T(K,NW)
-        
-	ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
+                    + DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
+	ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ?
 	btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
 	COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
      END DO ! L vertical loop
      
-     !last level
-     L              = L_NLAYRAD
-     K              = 2*L+1
-     
+     ! Last level
+     L           = L_NLAYRAD
+     K           = 2*L+1
      atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind)) &
                  + ASF_T(K,NW)*DHAZES_T(K,NW)
      btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) &
                  + DHAZES_T(K,NW)
-     ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW)
+     ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ?
      btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
      COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
@@ -316,6 +309,4 @@
   END DO                    ! NW spectral loop
 
-! ===========================================================================================
-
   DO NG=1,L_NGAUSS
     DO NW=1,L_NSPECTV
@@ -328,5 +319,5 @@
       END DO ! L vertical loop
 
-        !     No vertical averaging on bottom layer
+        ! Last level
 
         L              = L_NLAYRAD
@@ -356,5 +347,5 @@
 
 
-!  Titan's outputs (J.V.O, 2016)===============================================
+!  Titan's outputs (JVO, 2016)===============================================
 !      do l=1,L_NLAYRAD
 !         do nw=1,L_NSPECTV
Index: trunk/LMDZ.TITAN/libf/phytitan/phyetat0_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/phyetat0_mod.F90	(revision 1681)
+++ trunk/LMDZ.TITAN/libf/phytitan/phyetat0_mod.F90	(revision 1722)
@@ -10,4 +10,6 @@
                      emis,q2,qsurf)
 
+! to use  'getin_p'
+      use ioipsl_getin_p_mod, only: getin_p
 
   use tabfi_mod, only: tabfi
@@ -17,4 +19,5 @@
                      get_field, get_var, inquire_field, &
                      inquire_dimension, inquire_dimension_length
+  use callkeys_mod, only: surfalbedo,surfemis
 
   implicit none
@@ -106,5 +109,5 @@
   endif
 else
-  phisfi(:)=0
+  phisfi(:)=0.
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: surface geopotential <phisfi> range:", &
@@ -118,5 +121,9 @@
   endif
 else
-  albedodat(:)=0.5 ! would be better to read value from def file...
+  ! If no startfi file, use parameter surfalbedo in def file
+  surfalbedo=0.5
+  call getin_p("surfalbedo",surfalbedo)
+  print*,"surfalbedo",surfalbedo
+  albedodat(:)=surfalbedo
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: Bare ground albedo <albedodat> range:", &
@@ -130,5 +137,5 @@
   endif
 else
-  zmea(:)=0
+  zmea(:)=0.
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: <ZMEA> range:", &
@@ -142,5 +149,5 @@
   endif
 else
-  zstd(:)=0
+  zstd(:)=0.
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: <ZSTD> range:", &
@@ -154,5 +161,5 @@
   endif
 else
-  zsig(:)=0
+  zsig(:)=0.
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: <ZSIG> range:", &
@@ -166,5 +173,5 @@
   endif
 else
-  zgam(:)=0
+  zgam(:)=0.
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: <ZGAM> range:", &
@@ -178,5 +185,5 @@
   endif
 else
-  zthe(:)=0
+  zthe(:)=0.
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: <ZTHE> range:", &
@@ -190,5 +197,5 @@
   endif
 else
-  tsurf(:)=0 ! will be updated afterwards in physiq !
+  tsurf(:)=0. ! will be updated afterwards in physiq !
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: Surface temperature <tsurf> range:", &
@@ -202,5 +209,9 @@
   endif
 else
-  emis(:)=1 ! would be better to read value from def file...
+  ! If no startfi file, use parameter surfemis in def file
+  surfemis=1.0
+  call getin_p("surfemis",surfemis)
+  print*,"surfemis",surfemis
+  emis(:)=surfemis
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: Surface emissivity <emis> range:", &
@@ -214,5 +225,5 @@
   endif
 else
-  q2(:,:)=0
+  q2(:,:)=0.
 endif ! of if (startphy_file)
 write(*,*) "phyetat0: PBL wind variance <q2> range:", &
@@ -231,5 +242,5 @@
       endif
     else
-      qsurf(:,iq)=0
+      qsurf(:,iq)=0.
     endif ! of if (startphy_file)
     write(*,*) "phyetat0: Surface tracer <",trim(txt),"> range:", &
Index: trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 1681)
+++ trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90	(revision 1722)
@@ -1186,9 +1186,11 @@
       if(callsoil)then
          TsS = SUM(cell_area(:)*tsoil(:,nsoilmx))/totarea        ! mean temperature at bottom soil layer
-         print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]     ave[Tdeep]'
-         print*,Ts1,Ts2,Ts3,TsS
+         if (is_master) then
+            print*,'          ave[Tsurf]             min[Tsurf]             max[Tsurf]             ave[Tdeep]'
+            print*,Ts1,Ts2,Ts3,TsS
+         end if
       else
-               if (is_master) then
-            print*,'  ave[Tsurf]     min[Tsurf]     max[Tsurf]'
+         if (is_master) then
+            print*,'          ave[Tsurf]             min[Tsurf]             max[Tsurf]'
             print*,Ts1,Ts2,Ts3
          endif
@@ -1279,5 +1281,5 @@
 
 
-      print*,'--> Ls =',zls*180./pi
+      if (is_master) print*,'--> Ls =',zls*180./pi
       
       
Index: trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90	(revision 1681)
+++ trunk/LMDZ.TITAN/libf/phytitan/radcommon_h.F90	(revision 1722)
@@ -7,5 +7,5 @@
 !
 !                             radcommon.h
-!v
+!
 !----------------------------------------------------------------------C
 !
@@ -130,5 +130,5 @@
       real*8,save :: gweight(L_NGAUSS)
 !$OMP THREADPRIVATE(QREFvis,QREFir,omegaREFvis,omegaREFir,& 	! gweight read by master in sugas_corrk
-		!$OMP tstellar,planckir,PTOP,TAUREF)
+		!$OMP tstellar,planckir,PTOP)
 
 !     If the gas optical depth (top to the surface) is less than
@@ -154,13 +154,3 @@
 !$OMP THREADPRIVATE(glat,eclipse)
 
-contains
-
-      subroutine ini_radcommon_h
-      use radinc_h, only: L_LEVELS
-      implicit none
-      
-        allocate(TAUREF(L_LEVELS+1))
-      
-      end subroutine ini_radcommon_h
-
 end module radcommon_h
Index: trunk/LMDZ.TITAN/libf/phytitan/wstats.F90
===================================================================
--- trunk/LMDZ.TITAN/libf/phytitan/wstats.F90	(revision 1681)
+++ trunk/LMDZ.TITAN/libf/phytitan/wstats.F90	(revision 1722)
@@ -219,5 +219,5 @@
       start=(/1,1,indx,0/)
       if (klon_glo>1) then !general case
-        sizes=(/nbp_lon+1,nbp_lev,1,0/)
+        sizes=(/nbp_lon+1,nbp_lat,1,0/)
       else
         sizes=(/1,1,1,0/)
