Index: LMDZ4/trunk/libf/cosp/lmd_ipsl_stats.F90
===================================================================
--- LMDZ4/trunk/libf/cosp/lmd_ipsl_stats.F90	(revision 1279)
+++ LMDZ4/trunk/libf/cosp/lmd_ipsl_stats.F90	(revision 1414)
@@ -1,24 +1,24 @@
 ! Copyright (c) 2009, Centre National de la Recherche Scientifique
 ! All rights reserved.
-! 
-! Redistribution and use in source and binary forms, with or without modification, are permitted 
+!
+! Redistribution and use in source and binary forms, with or without modification, are permitted
 ! provided that the following conditions are met:
-! 
-!     * Redistributions of source code must retain the above copyright notice, this list 
+!
+!     * Redistributions of source code must retain the above copyright notice, this list
 !       of conditions and the following disclaimer.
 !     * Redistributions in binary form must reproduce the above copyright notice, this list
-!       of conditions and the following disclaimer in the documentation and/or other materials 
+!       of conditions and the following disclaimer in the documentation and/or other materials
 !       provided with the distribution.
 !     * Neither the name of the LMD/IPSL/CNRS/UPMC nor the names of its
-!       contributors may be used to endorse or promote products derived from this software without 
+!       contributors may be used to endorse or promote products derived from this software without
 !       specific prior written permission.
-! 
-! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR 
-! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND 
-! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR 
-! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL 
-! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 
-! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER 
-! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT 
+!
+! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
+! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
+! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
+! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
+! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
 ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
@@ -39,13 +39,13 @@
 ! -----------------------------------------------------------------------------------
 ! Lidar outputs :
-! 
+!
 ! Diagnose cloud fraction (3D cloud fraction + low/middle/high/total cloud fraction
 ! from the lidar signals (ATB and molecular ATB) computed from model outputs
 !      +
 ! Compute CFADs of lidar scattering ratio SR and of depolarization index
-! 
+!
 ! Authors: Sandrine Bony and Helene Chepfer (LMD/IPSL, CNRS, UPMC, France).
 !
-! December 2008, S. Bony,  H. Chepfer and J-L. Dufresne : 
+! December 2008, S. Bony,  H. Chepfer and J-L. Dufresne :
 ! - change of the cloud detection threshold S_cld from 3 to 5, for better
 ! with both day and night observations. The optical thinest clouds are missed.
@@ -53,4 +53,10 @@
 ! December 2008, A. Bodas-Salcedo:
 ! - Dimensions of pmol reduced to (npoints,llm)
+! August 2009, A. Bodas-Salcedo:
+! - Warning message regarding PARASOL being valid only over ocean deleted.
+! February 2010, A. Bodas-Salcedo:
+! - Undef passed into cosp_cfad_sr
+! June 2010, T. Yokohata, T. Nishimura and K. Ogochi
+! Optimisation of COSP_CFAD_SR
 !
 ! Version 1.0 (June 2007)
@@ -70,7 +76,7 @@
 
       real undef                    ! undefined value
-      real pnorm(npoints,ncol,llm)  ! lidar ATB 
+      real pnorm(npoints,ncol,llm)  ! lidar ATB
       real pmol(npoints,llm)        ! molecular ATB
-      real land(npoints)            ! Landmask [0 - Ocean, 1 - Land]    
+      real land(npoints)            ! Landmask [0 - Ocean, 1 - Land]
       real pplay(npoints,llm)       ! pressure on model levels (Pa)
       logical ok_lidar_cfad         ! true if lidar CFAD diagnostics need to be computed
@@ -78,13 +84,13 @@
 
 ! c outputs :
-      real lidarcld(npoints,llm)     ! 3D "lidar" cloud fraction 
+      real lidarcld(npoints,llm)     ! 3D "lidar" cloud fraction
       real cldlayer(npoints,ncat)    ! "lidar" cloud fraction (low, mid, high, total)
-      real cfad2(npoints,max_bin,llm) ! CFADs of SR  
-      real srbval(max_bin)           ! SR bins in CFADs  
+      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
 
 ! c threshold for cloud detection :
-      real S_clr 
-      parameter (S_clr = 1.2) 
+      real S_clr
+      parameter (S_clr = 1.2)
       real S_cld
 !      parameter (S_cld = 3.0)  ! Previous thresold for cloud detection
@@ -102,10 +108,5 @@
 ! c 0- Initializations
 ! c -------------------------------------------------------
-! Parasol reflectance algorithm is not valid over land. Write
-! a warning if there is no land. Landmask [0 - Ocean, 1 - Land] 
-      IF ( MAXVAL(land(:)) .EQ. 0.0) THEN
-          WRITE (*,*) 'WARNING. PARASOL reflectance is not valid over land' &
-             & ,' and there is only land'
-      END IF
+!
 
 !  Should be modified in future version
@@ -116,5 +117,5 @@
 ! c -------------------------------------------------------
 !
-!       where ((pnorm.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 )) 
+!       where ((pnorm.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
 !          x3d = pnorm/pmol
 !       elsewhere
@@ -124,5 +125,5 @@
       do ic = 1, ncol
         pnorm_c = pnorm(:,ic,:)
-        where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 )) 
+        where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
             x3d_c = pnorm_c/pmol
         elsewhere
@@ -142,5 +143,5 @@
 
 ! c -------------------------------------------------------
-! c 3- CFADs 
+! c 3- CFADs
 ! c -------------------------------------------------------
       if (ok_lidar_cfad) then
@@ -148,5 +149,5 @@
 ! c CFADs of subgrid-scale lidar scattering ratios :
 ! c -------------------------------------------------------
-      CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin, &
+      CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin,undef, &
                  x3d, &
                  S_att,S_clr,xmax,cfad2,srbval)
@@ -172,16 +173,16 @@
 ! if land=0 -> parasolrefl=parasolrefl
         parasolrefl(:,k) = parasolrefl(:,k) * MAX(1.0-land(:),0.0) &
-                           + (1.0 - MAX(1.0-land(:),0.0))*undef 
+                           + (1.0 - MAX(1.0-land(:),0.0))*undef
       enddo
 
       RETURN
       END SUBROUTINE diag_lidar
-	  
-	  
+
+
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 !-------------------- FUNCTION COSP_CFAD_SR ------------------------
 ! Author: Sandrine Bony (LMD/IPSL, CNRS, Paris)
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-      SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins, &
+      SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins,undef, &
                       x,S_att,S_clr,xmax,cfad,srbval)
       IMPLICIT NONE
@@ -205,5 +206,5 @@
 ! Input arguments
       integer Npoints,Ncolumns,Nlevels,Nbins
-      real xmax,S_att,S_clr
+      real xmax,S_att,S_clr,undef
 ! Input-output arguments
       real x(Npoints,Ncolumns,Nlevels)
@@ -213,4 +214,5 @@
 ! Local variables
       integer i, j, k, ib
+      real srbval_ext(0:Nbins)
 
 ! c -------------------------------------------------------
@@ -235,26 +237,26 @@
       cfad(:,:,:) = 0.0
 
-
+      srbval_ext(1:Nbins) = srbval
+      srbval_ext(0) = -1.0
 ! c -------------------------------------------------------
 ! c c- Compute CFAD
 ! c -------------------------------------------------------
 
-        do j = Nlevels, 1, -1 
-          do k = 1, Ncolumns
-              where ( x(:,k,j).le.srbval(1) ) &
-                        cfad(:,1,j) = cfad(:,1,j) + 1.0
-          enddo  !k
-        enddo  !j
-
-      do ib = 2, Nbins
-        do j = Nlevels, 1, -1 
-          do k = 1, Ncolumns
-              where ( ( x(:,k,j).gt.srbval(ib-1) ) .and. ( x(:,k,j).le.srbval(ib) ) ) &
-                        cfad(:,ib,j) = cfad(:,ib,j) + 1.0
-          enddo  !k
-        enddo  !j
-      enddo  !ib
-
-      cfad(:,:,:) = cfad(:,:,:) / float(Ncolumns)
+      do j = 1, Nlevels
+         do ib = 1, Nbins
+            do k = 1, Ncolumns
+               do i = 1, Npoints
+                  if (x(i,k,j) /= undef) then
+                     if ((x(i,k,j).gt.srbval_ext(ib-1)).and.(x(i,k,j).le.srbval_ext(ib))) &
+                          cfad(i,ib,j) = cfad(i,ib,j) + 1.0
+                  else 
+                     cfad(i,ib,j) = undef
+                  endif
+               enddo
+            enddo
+         enddo
+      enddo
+
+      where (cfad .ne. undef)  cfad = cfad / float(Ncolumns)
 
 ! c -------------------------------------------------------
@@ -264,5 +266,5 @@
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 !-------------------- SUBROUTINE COSP_CLDFRAC -------------------
-! c Purpose: Cloud fraction diagnosed from lidar measurements 
+! c Purpose: Cloud fraction diagnosed from lidar measurements
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
       SUBROUTINE COSP_CLDFRAC(Npoints,Ncolumns,Nlevels,Ncat, &
@@ -288,7 +290,14 @@
       real nsub(Npoints,Nlevels)
 
-
-! ---------------------------------------------------------------
-! 1- initialization 
+      real cldlay1(Npoints,Ncolumns)
+      real cldlay2(Npoints,Ncolumns)
+      real cldlay3(Npoints,Ncolumns)
+      real nsublay1(Npoints,Ncolumns)
+      real nsublay2(Npoints,Ncolumns)
+      real nsublay3(Npoints,Ncolumns)
+
+
+! ---------------------------------------------------------------
+! 1- initialization
 ! ---------------------------------------------------------------
 
@@ -317,5 +326,5 @@
 
 ! number of usefull sub-columns:
-         where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef)  ) 
+         where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef)  )
            srok(:,:,k)=1.0
          elsewhere
@@ -329,27 +338,45 @@
 ! categories) :
 ! ---------------------------------------------------------------
-! Test abderr
+      lidarcld = 0.0
+      nsub = 0.0
+
+!! XXX: Use cldlay[1-3] and nsublay[1-3] to avoid bank-conflicts.
+      cldlay1 = 0.0
+      cldlay2 = 0.0
+      cldlay3 = 0.0
+      cldlay(:,:,4) = 0.0 !! XXX: Ncat == 4
+      nsublay1 = 0.0
+      nsublay2 = 0.0
+      nsublay3 = 0.0
+      nsublay(:,:,4) = 0.0
       do k = Nlevels, 1, -1
        do ic = 1, Ncolumns
         do ip = 1, Npoints
-          iz=1
-          p1 = pplay(ip,k)
-          if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
-            iz=3
-          else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid clouds
-            iz=2
+         p1 = pplay(ip,k)
+
+         if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
+            cldlay3(ip,ic) = MAX(cldlay3(ip,ic), cldy(ip,ic,k))
+            nsublay3(ip,ic) = MAX(nsublay3(ip,ic), srok(ip,ic,k))
+         else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid clouds
+            cldlay2(ip,ic) = MAX(cldlay2(ip,ic), cldy(ip,ic,k))
+            nsublay2(ip,ic) = MAX(nsublay2(ip,ic), srok(ip,ic,k))
+         else
+            cldlay1(ip,ic) = MAX(cldlay1(ip,ic), cldy(ip,ic,k))
+            nsublay1(ip,ic) = MAX(nsublay1(ip,ic), srok(ip,ic,k))
          endif
 
-         cldlay(ip,ic,iz) = MAX(cldlay(ip,ic,iz),cldy(ip,ic,k))
-         cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4),cldy(ip,ic,k))
+         cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4), cldy(ip,ic,k))
          lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k)
-
-         nsublay(ip,ic,iz) = MAX(nsublay(ip,ic,iz),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)
-
         enddo
        enddo
       enddo
+      cldlay(:,:,1) = cldlay1
+      cldlay(:,:,2) = cldlay2
+      cldlay(:,:,3) = cldlay3
+      nsublay(:,:,1) = nsublay1
+      nsublay(:,:,2) = nsublay2
+      nsublay(:,:,3) = nsublay3
 
 ! -- grid-box 3D cloud fraction
@@ -369,6 +396,6 @@
        do ic = 1, Ncolumns
 
-          cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz)    
-          nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz) 
+          cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz)
+          nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz)
 
        enddo
@@ -383,4 +410,4 @@
       END SUBROUTINE COSP_CLDFRAC
 ! ---------------------------------------------------------------
-	  
+
 END MODULE MOD_LMD_IPSL_STATS
