Index: LMDZ5/trunk/libf/phylmd/cosp/lmd_ipsl_stats.F90
===================================================================
--- LMDZ5/trunk/libf/phylmd/cosp/lmd_ipsl_stats.F90	(revision 2428)
+++ LMDZ5/trunk/libf/phylmd/cosp/lmd_ipsl_stats.F90	(revision 2955)
@@ -36,6 +36,6 @@
       SUBROUTINE diag_lidar(npoints,ncol,llm,max_bin,nrefl &
                   ,tmp,pnorm,pnorm_perp,pmol,refl,land,pplay,undef,ok_lidar_cfad &
-                  ,cfad2,srbval,ncat,lidarcld,lidarcldphase,cldlayer,cldlayerphase &
-                  ,lidarcldtmp,parasolrefl)
+                  ,cfad2,srbval,ncat,ntype,lidarcld,lidarcldtype,lidarcldphase,cldlayer & !OPAQ
+                  ,cldtype,cldlayerphase,lidarcldtmp,parasolrefl,vgrid_z,profSR)          !OPAQ !TIBO
 !
 ! -----------------------------------------------------------------------------------
@@ -82,4 +82,5 @@
       integer max_bin               ! nb of bins for SR CFADs
       integer ncat                  ! nb of cloud layer types (low,mid,high,total)
+      integer ntype                 ! nb of OPAQ products (opaque and thin clouds, z_opaque) !OPAQ
       integer nrefl                 ! nb of solar zenith angles for parasol reflectances
 
@@ -93,13 +94,18 @@
       real tmp(npoints,llm)         ! temp at each levels
       real pnorm_perp(npoints,ncol,llm)  ! lidar perpendicular ATB
+      real vgrid_z(llm)             ! mid-level altitude of the output vertical grid         !OPAQ
 
 ! c outputs :
       real lidarcld(npoints,llm)     ! 3D "lidar" cloud fraction
+      real lidarcldtype(npoints,llm,ntype+1)   ! 3D "lidar" OPAQ type fraction + opacity     !OPAQ
       real sub(npoints,llm)     ! 3D "lidar" indice
       real cldlayer(npoints,ncat)    ! "lidar" cloud layer fraction (low, mid, high, total)
+      real cldtype(npoints,ntype)  ! "lidar" OPAQ type covers (opaque/thin cloud + z_opaque) !OPAQ
 
       real cfad2(npoints,max_bin,llm) ! CFADs of SR
       real srbval(max_bin)           ! SR bins in CFADs
       real parasolrefl(npoints,nrefl)! grid-averaged parasol reflectance
+!     real profSR(npoints,ncol,llm)  ! tableau avec les subcolumns SR !TIBO
+      real profSR(npoints,llm,ncol)  ! tableau avec les subcolumns SR !TIBO2
 
 ! c threshold for cloud detection :
@@ -110,4 +116,5 @@
       real S_att
       parameter (S_att = 0.01)
+!      parameter (S_att = 0.06)  !OPAQ ! Threshold for "surface detection" equivalent
 
 ! c local variables :
@@ -147,4 +154,6 @@
         end where
          x3d(:,ic,:) = x3d_c
+!	profSR(:,ic,:) = x3d(:,ic,:) !TIBO
+	profSR(:,:,ic) = x3d(:,ic,:) !TIBO2
       enddo
 
@@ -157,4 +166,7 @@
               tmp,x3d,pnorm,pnorm_perp,pplay, S_att,S_cld,S_cld_att,undef,lidarcld, &
               cldlayer,lidarcldphase,sub,cldlayerphase,lidarcldtmp)
+
+    CALL COSP_OPAQ(npoints,ncol,llm,ntype,x3d,S_cld,undef,lidarcldtype,            & !OPAQ
+                   cldtype,vgrid_z)                                                  !OPAQ
 
 ! c -------------------------------------------------------
@@ -983,4 +995,187 @@
 ! ---------------------------------------------------------------
 
+! BEGINNING OF OPAQ CHANGES
+    ! ####################################################################################
+    ! SUBROUTINE cosp_opaq
+    ! Conventions: Ntype must be equal to 3 (opaque cloud, thin cloud, z_opaque)
+    ! ####################################################################################
+    SUBROUTINE COSP_OPAQ(Npoints,Ncolumns,Nlevels,Ntype,x,S_cld,undef,lidarcldtype,   &
+                         cldtype,vgrid_z)
+
+      IMPLICIT NONE
+! Input arguments
+      integer Npoints,Ncolumns,Nlevels,Ntype
+      real x(Npoints,Ncolumns,Nlevels)
+      real S_cld
+      real undef
+      real vgrid_z(Nlevels)
+! Output :
+      real lidarcldtype(Npoints,Nlevels,Ntype+1) ! 3D "lidar" OPAQ type + opacity fraction
+      real cldtype(Npoints,Ntype)              ! opaque and thin cloud covers, z_opaque
+! Local variables
+      integer ip, k, iz, ic, zopac
+      real p1
+      real cldy(Npoints,Ncolumns,Nlevels)
+      real cldyopaq(Npoints,Ncolumns,Nlevels)
+      real srok(Npoints,Ncolumns,Nlevels)
+      real srokopaq(Npoints,Ncolumns,Nlevels)
+      real cldlay(Npoints,Ncolumns,Ntype+1)  ! opaque, thin, z_opaque and all cloud cover
+      real nsublay(Npoints,Ncolumns,Ntype+1) ! opaque, thin, z_opaque and all cloud cover
+      real nsublayer(Npoints,Ntype)
+      real nsub(Npoints,Nlevels)
+      real nsubopaq(Npoints,Nlevels)
+      real S_att_opaq
+      real S_att
+  
+    ! ####################################################################################
+	! 1) Initialize    
+    ! ####################################################################################
+    cldtype               = 0.0
+    lidarcldtype          = 0.0
+    nsub                  = 0.0
+    nsubopaq              = 0.0
+    cldlay                = 0.0
+    nsublay               = 0.0
+    nsublayer             = 0.0
+    S_att_opaq            = 0.06 ! Fully Attenuated threshold, from Guzman et al. 2017, JGR-A
+    S_att                 = 0.01 
+
+    ! ####################################################################################
+    ! 2) Cloud detection and Fully attenuated layer detection
+    ! ####################################################################################
+    do k=1,Nlevels
+       ! Cloud detection at subgrid-scale:
+       where ( (x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
+          cldy(:,:,k)=1.0
+       elsewhere
+          cldy(:,:,k)=0.0
+       endwhere
+       ! Fully attenuated layer detection at subgrid-scale:
+       where ( (x(:,:,k) .gt. 0.0) .and. (x(:,:,k) .lt. S_att_opaq) .and. (x(:,:,k) .ne. undef) )
+          cldyopaq(:,:,k)=1.0
+       elsewhere
+          cldyopaq(:,:,k)=0.0
+       endwhere
+
+       ! Number of useful sub-column layers:
+       where ( (x(:,:,k) .gt. S_att) .and. (x(:,:,k) .ne. undef) )
+          srok(:,:,k)=1.0
+       elsewhere
+          srok(:,:,k)=0.0
+       endwhere
+       ! Number of useful sub-columns layers for z_opaque 3D fraction:
+       where ( (x(:,:,k) .gt. 0.0) .and. (x(:,:,k) .ne. undef) )
+          srokopaq(:,:,k)=1.0
+       elsewhere
+          srokopaq(:,:,k)=0.0
+       endwhere
+    enddo
+
+    ! ####################################################################################
+    ! 3) Grid-box 3D OPAQ product fraction and cloud type cover (opaque/thin) + mean z_opaque
+    ! ####################################################################################
+
+    do k= Nlevels,1,-1
+       do ic = 1, Ncolumns
+          do ip = 1, Npoints
+
+             cldlay(ip,ic,1)   = MAX(cldlay(ip,ic,1),cldyopaq(ip,ic,k)) ! Opaque clouds
+             cldlay(ip,ic,4)   = MAX(cldlay(ip,ic,4),cldy(ip,ic,k))     ! All clouds
+
+             nsublay(ip,ic,1)  = MAX(nsublay(ip,ic,1),srok(ip,ic,k))
+             nsublay(ip,ic,2)  = MAX(nsublay(ip,ic,2),srok(ip,ic,k))
+!             nsublay(ip,ic,4)  = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
+             nsub(ip,k)        = nsub(ip,k) + srok(ip,ic,k)
+             nsubopaq(ip,k)    = nsubopaq(ip,k) + srokopaq(ip,ic,k)
+
+          enddo
+       enddo
+    enddo   
+
+! OPAQ variables
+     do ic = 1, Ncolumns
+        do ip = 1, Npoints
+
+     ! Declaring non-opaque cloudy profiles as thin cloud profiles
+	   if ( (cldlay(ip,ic,4) .eq. 1.0) .and. (cldlay(ip,ic,1) .eq. 0.0) ) then
+	      cldlay(ip,ic,2)  =  1.0
+ 	   endif
+
+     ! Filling in 3D and 2D variables
+
+     ! Opaque cloud profiles
+	   if ( cldlay(ip,ic,1) .eq. 1.0 ) then
+	      zopac = 0.0
+	      do k=2,Nlevels
+     ! Declaring opaque cloud fraction and z_opaque altitude for 3D and 2D variables
+	         if ( (cldy(ip,ic,k) .eq. 1.0) .and. (zopac .eq. 0.0) ) then
+		    lidarcldtype(ip,k-1,3) = lidarcldtype(ip,k-1,3) + 1.0
+		    cldlay(ip,ic,3)        = vgrid_z(k-1) !z_opaque altitude
+		    nsublay(ip,ic,3)       = 1.0
+		    zopac = 1.0
+		 endif
+	         if ( cldy(ip,ic,k) .eq. 1.0 ) then
+		    lidarcldtype(ip,k,1)   = lidarcldtype(ip,k,1) + 1.0
+                 endif
+	      enddo
+	   endif
+
+     ! Thin cloud profiles
+	   if ( cldlay(ip,ic,2) .eq. 1.0 ) then
+	      do k=1,Nlevels
+     ! Declaring thin cloud fraction for 3D variable
+                 if ( cldy(ip,ic,k) .eq. 1.0 ) then
+                    lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1.0
+                 endif
+	      enddo
+           endif
+
+       enddo
+    enddo   
+
+    ! 3D cloud types fraction (opaque=1 and thin=2)
+    where ( nsub(:,:) .gt. 0.0 )
+       lidarcldtype(:,:,1) = lidarcldtype(:,:,1)/nsub(:,:)
+       lidarcldtype(:,:,2) = lidarcldtype(:,:,2)/nsub(:,:)
+    elsewhere
+       lidarcldtype(:,:,1) = undef
+       lidarcldtype(:,:,2) = undef
+    endwhere
+    ! 3D z_opaque fraction (=3)
+    where ( nsubopaq(:,:) .gt. 0.0 )
+       lidarcldtype(:,:,3) = lidarcldtype(:,:,3)/nsubopaq(:,:)
+    elsewhere
+       lidarcldtype(:,:,3) = undef
+    endwhere
+    ! 3D opacity fraction (=4) !Summing z_opaque fraction from TOA(k=Nlevels) to SFC(k=1)
+       lidarcldtype(:,Nlevels,4) = lidarcldtype(:,Nlevels,3)
+    do ip = 1, Npoints
+     	do k = Nlevels-1, 1, -1
+           if ( lidarcldtype(ip,k,3) .ne. undef ) then
+	      lidarcldtype(ip,k,4) = lidarcldtype(ip,k+1,4) + lidarcldtype(ip,k,3)
+           endif
+	enddo
+    enddo
+    where ( nsubopaq(:,:) .eq. 0.0 )
+       lidarcldtype(:,:,4) = undef
+    endwhere
+
+    ! Layered cloud types (opaque, thin and z_opaque 2D variables)
+
+    do iz = 1, Ntype
+       do ic = 1, Ncolumns
+          cldtype(:,iz)   = cldtype(:,iz)   + cldlay(:,ic,iz)
+          nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz)
+       enddo
+    enddo
+    where (nsublayer(:,:) .gt. 0.0)
+       cldtype(:,:) = cldtype(:,:)/nsublayer(:,:)
+    elsewhere
+       cldtype(:,:) = undef
+    endwhere
+
+  END SUBROUTINE COSP_OPAQ
+! END OF OPAQ CHANGES
+
 
 END MODULE MOD_LMD_IPSL_STATS
