Index: /trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90
===================================================================
--- /trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 4147)
+++ /trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90	(revision 4148)
@@ -14,8 +14,8 @@
       use gases_h
       USE tracer_h
-      use callkeys_mod, only: global1d, szangle
       use comcstfi_mod, only: pi, mugaz, cpp
-      use callkeys_mod, only: diurnal,tracer,seashaze,corrk_recombin,   &
-                              strictboundcorrk,specOLR,diagdtau,        &
+      use callkeys_mod, only: global1d, szangle,                      &
+                              diurnal,tracer,seashaze,corrk_recombin, &
+                              strictboundcorrk,specOLR,diagdtau,      &
                               tplanckmin,tplanckmax,callclouds,Fcloudy
       use geometry_mod, only: latitude
@@ -60,5 +60,5 @@
       INTEGER,INTENT(IN) :: nq                     ! Number of tracers.
       REAL,INTENT(IN) :: qsurf(ngrid,nq)           ! Tracers on surface (kg.m-2).
-      REAL,INTENT(IN) :: zday		           ! Time elapsed since Ls=0 (sols).
+      REAL,INTENT(IN) :: zday                      ! Time elapsed since Ls=0 (sols).
       REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV)   ! Spectral Short Wavelengths Albedo. By MT2015
       REAL,INTENT(IN) :: emis(ngrid)               ! Long Wave emissivity.
@@ -198,5 +198,5 @@
 !             I.  Initialization on every call   
 !=======================================================================
- 
+
 
       ! How much light do we get ?
@@ -513,8 +513,8 @@
               dtauv_dc,tauv_dc,taucumv_dc,wbarv_dc,cosbv_dc,tauray,taugsurf,seashazefact,&
               diag_opthv,diag_opttv_dc,cdcolumn)
-         
+
          ! Mean opacity, ssa and asf :
-         where ((exp(-dtauv_cc(:,:,:)).ge.1.d-40) .and. (exp(-dtauv_dc(:,:,:)).ge.1.d-40))
-            dtauv(:,:,:) = -log((1-Fcloudy)*exp(-dtauv_cc(:,:,:)) + Fcloudy*exp(-dtauv_dc(:,:,:)))
+         where (dtauv_cc(:,:,:) .le. 100 .and. dtauv_dc(:,:,:) .le. 100.)
+            dtauv(:,:,:) = (1-Fcloudy)*dtauv_cc(:,:,:) + Fcloudy*dtauv_dc(:,:,:)
          elsewhere
             dtauv(:,:,:) = dtauv_dc(:,:,:) ! No need to average...
@@ -524,6 +524,6 @@
                taucumv(1,nw,ng) = 0.0d0
                do k = 2, L_LEVELS
-                  if ((exp(-taucumv_cc(k,nw,ng)).ge.1.d-40) .and. (exp(-taucumv_dc(k,nw,ng)).ge.1.d-40)) then
-                     taucumv(k,nw,ng) = taucumv(k-1,nw,ng) - log((1-Fcloudy)*exp(-taucumv_cc(k,nw,ng)) + Fcloudy*exp(-taucumv_dc(k,nw,ng)))
+                  if ((taucumv_cc(k,nw,ng).le.100.) .and. (taucumv_dc(k,nw,ng).le.100.)) then
+                     taucumv(k,nw,ng) = taucumv(k-1,nw,ng) + (1-Fcloudy)*taucumv_cc(k,nw,ng) + Fcloudy*taucumv_dc(k,nw,ng)
                   else
                      taucumv(k,nw,ng) = taucumv(k-1,nw,ng) + taucumv_dc(k,nw,ng) ! No need to average...
@@ -536,8 +536,14 @@
             end do
          end do
-         wbarv = (1-Fcloudy) * wbarv_cc + Fcloudy * wbarv_dc
-         cosbv = (1-Fcloudy) * cosbv_cc + Fcloudy * cosbv_dc
-        
+
+         wbarv = ((1-Fcloudy) * wbarv_cc*dtauv_cc            + Fcloudy * wbarv_dc *dtauv_dc)
+         wbarv = wbarv /((1-Fcloudy) * dtauv_cc              + Fcloudy * dtauv_dc  + 1.e-30)
+
+         cosbv = ((1-Fcloudy) * cosbv_cc * wbarv_cc*dtauv_cc + Fcloudy * cosbv_dc  * wbarv_dc *dtauv_dc)
+         cosbv = cosbv /((1-Fcloudy) * wbarv_cc*dtauv_cc     + Fcloudy * wbarv_dc *dtauv_dc + 1.e-30)
+         !------------------------------------------------------------------------------
+         
          ! Diagnostics for clouds :
+
          if (callclouds) then
             where (diag_opttv_cc(:,:,1) .lt. 1.d-30)
@@ -585,5 +591,4 @@
          end if
 
-
          ! Equivalent Albedo Calculation (for OUTPUT). MT2015
          if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight.       
@@ -620,6 +625,6 @@
 
          ! Mean opacity, ssa and asf :
-         where ((exp(-dtaui_cc(:,:,:)).ge.1.d-40) .and. (exp(-dtaui_dc(:,:,:)).ge.1.d-40))
-            dtaui(:,:,:) = -log((1-Fcloudy)*exp(-dtaui_cc(:,:,:)) + Fcloudy*exp(-dtaui_dc(:,:,:)))
+         where (dtaui_cc(:,:,:).le.100. .and. dtaui_dc(:,:,:).le.100.)
+            dtaui(:,:,:) = (1-Fcloudy)*dtaui_cc(:,:,:) + Fcloudy*dtaui_dc(:,:,:)
          elsewhere
             dtaui(:,:,:) = dtaui_dc(:,:,:) ! No need to average...
@@ -629,6 +634,6 @@
                taucumi(1,nw,ng) = 0.0d0
                do k = 2, L_LEVELS
-                  if ((exp(-taucumi_cc(k,nw,ng)).ge.1.d-40) .and. (exp(-taucumi_dc(k,nw,ng)).ge.1.d-40)) then
-                     taucumi(k,nw,ng) = taucumi(k-1,nw,ng) - log((1-Fcloudy)*exp(-taucumi_cc(k,nw,ng)) + Fcloudy*exp(-taucumi_dc(k,nw,ng)))
+                  if (taucumi_cc(k,nw,ng).le.100. .and. taucumi_dc(k,nw,ng).le.100.) then
+                     taucumi(k,nw,ng) = taucumi(k-1,nw,ng) + (1-Fcloudy)*taucumi_cc(k,nw,ng) + Fcloudy*taucumi_dc(k,nw,ng)
                   else
                      taucumi(k,nw,ng) = taucumi(k-1,nw,ng) + taucumi_dc(k,nw,ng) ! No need to average...
@@ -637,6 +642,10 @@
             end do
          end do
-         wbari = (1-Fcloudy) * wbari_cc + Fcloudy * wbari_dc
-         cosbi = (1-Fcloudy) * cosbi_cc + Fcloudy * cosbi_dc
+          wbari = ((1-Fcloudy) * wbari_cc*dtaui_cc + Fcloudy * wbari_dc *dtaui_dc)
+          wbari = wbari /((1-Fcloudy) * dtaui_cc+ Fcloudy * dtaui_dc+1.e-10)
+          cosbi = ((1-Fcloudy) * cosbi_cc * wbari_cc* dtaui_cc + Fcloudy * cosbi_dc  * wbari_dc *dtaui_dc)
+          cosbi =  cosbi /((1-Fcloudy) * wbari_cc*dtaui_cc + Fcloudy * wbari_dc *dtaui_dc+1.e-10)
+
+         !----------------------------------------------------------------------
 
          ! Diagnostics for clouds :