Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp.F90
===================================================================
--- LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp.F90	(revision 1410)
+++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp.F90	(revision 1415)
@@ -23,4 +23,5 @@
 ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
+!!#include "cosp_defs.h"
 MODULE MOD_COSP
   USE MOD_COSP_TYPES
@@ -128,12 +129,12 @@
       !     and reff_zero    == .false.  Reff use in lidar and set to 0 for radar
   endif
-  if ((gbx%use_reff) .and. (reff_zero)) then ! Inconsistent choice. Want to use Reff but not inputs passed
-        print *, '---------- COSP ERROR ------------'
-        print *, ''
-        print *, 'use_reff==.true. but Reff is always zero'
-        print *, ''
-        print *, '----------------------------------'
-        stop
-  endif
+!  if ((gbx%use_reff) .and. (reff_zero)) then ! Inconsistent choice. Want to use Reff but not inputs passed
+!        print *, '---------- COSP ERROR ------------'
+!        print *, ''
+!        print *, 'use_reff==.true. but Reff is always zero'
+!        print *, ''
+!        print *, '----------------------------------'
+!        stop
+!  endif
   if ((.not. gbx%use_reff) .and. (reff_zero)) then ! No Reff in radar. Default in lidar
         gbx%Reff = DEFAULT_LIDAR_REFF
@@ -322,4 +323,5 @@
   integer :: Niter     ! Number of calls to cosp_simulator
   integer :: i,j,k
+  integer :: I_HYDRO
   real,dimension(:,:),pointer :: column_frac_out ! Array with one column of frac_out
   integer,parameter :: scops_debug=0    !  set to non-zero value to print out inputs for debugging in SCOPS
@@ -330,4 +332,5 @@
   real,dimension(:,:),allocatable :: frac_ls,prec_ls,frac_cv,prec_cv ! Cloud/Precipitation fraction in each model level
                                                                      ! Levels are from SURFACE to TOA
+  real,dimension(:,:),allocatable :: rho ! (Npoints, Nlevels). Atmospheric dens
   type(cosp_sghydro) :: sghydro   ! Subgrid info for hydrometeors en each iteration
 
@@ -378,6 +381,6 @@
         do k=1,Nlevels,1
             do i=1,Ncolumns,1
-                if (sgx%frac_out (j,i,Nlevels+1-k) .eq. 1) frac_ls(j,k)=frac_ls(j,k)+1.
-                if (sgx%frac_out (j,i,Nlevels+1-k) .eq. 2) frac_cv(j,k)=frac_cv(j,k)+1.
+                if (sgx%frac_out (j,i,Nlevels+1-k) == I_LSC) frac_ls(j,k)=frac_ls(j,k)+1.
+                if (sgx%frac_out (j,i,Nlevels+1-k) == I_CVC) frac_cv(j,k)=frac_cv(j,k)+1.
                 if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 1) prec_ls(j,k)=prec_ls(j,k)+1.
                 if (sgx%prec_frac(j,i,Nlevels+1-k) .eq. 2) prec_cv(j,k)=prec_cv(j,k)+1.
@@ -414,5 +417,5 @@
             !--------- Mixing ratios for clouds and Reff for Clouds and precip -------
             column_frac_out => sgx%frac_out(:,k,:)
-            where (column_frac_out == 1)     !+++++++++++ LS clouds ++++++++
+            where (column_frac_out == I_LSC)     !+++++++++++ LS clouds ++++++++
                 sghydro%mr_hydro(:,k,:,I_LSCLIQ) = gbx%mr_hydro(:,:,I_LSCLIQ)
                 sghydro%mr_hydro(:,k,:,I_LSCICE) = gbx%mr_hydro(:,:,I_LSCICE)
@@ -423,5 +426,5 @@
                 sghydro%Reff(:,k,:,I_LSSNOW)     = gbx%Reff(:,:,I_LSSNOW)
                 sghydro%Reff(:,k,:,I_LSGRPL)     = gbx%Reff(:,:,I_LSGRPL)
-            elsewhere (column_frac_out == 2) !+++++++++++ CONV clouds ++++++++
+            elsewhere (column_frac_out == I_CVC) !+++++++++++ CONV clouds ++++++++
                 sghydro%mr_hydro(:,k,:,I_CVCLIQ) = gbx%mr_hydro(:,:,I_CVCLIQ) 
                 sghydro%mr_hydro(:,k,:,I_CVCICE) = gbx%mr_hydro(:,:,I_CVCICE) 
@@ -434,9 +437,9 @@
             !--------- Precip -------
             if (.not. gbx%use_precipitation_fluxes) then
-                where (column_frac_out == 1)  !+++++++++++ LS Precipitation ++++++++
+                where (column_frac_out == I_LSC)  !+++++++++++ LS Precipitation ++++++++
                     sghydro%mr_hydro(:,k,:,I_LSRAIN) = gbx%mr_hydro(:,:,I_LSRAIN)
                     sghydro%mr_hydro(:,k,:,I_LSSNOW) = gbx%mr_hydro(:,:,I_LSSNOW)
                     sghydro%mr_hydro(:,k,:,I_LSGRPL) = gbx%mr_hydro(:,:,I_LSGRPL)
-                elsewhere (column_frac_out == 2) !+++++++++++ CONV Precipitation ++++++++
+                elsewhere (column_frac_out == I_CVC) !+++++++++++ CONV Precipitation ++++++++
                     sghydro%mr_hydro(:,k,:,I_CVRAIN) = gbx%mr_hydro(:,:,I_CVRAIN) 
                     sghydro%mr_hydro(:,k,:,I_CVSNOW) = gbx%mr_hydro(:,:,I_CVSNOW) 
@@ -488,5 +491,5 @@
                         sghydro%mr_hydro(:,:,:,I_LSRAIN),sghydro%mr_hydro(:,:,:,I_LSSNOW),sghydro%mr_hydro(:,:,:,I_LSGRPL), &
                         sghydro%mr_hydro(:,:,:,I_CVRAIN),sghydro%mr_hydro(:,:,:,I_CVSNOW))
-        endif
+       endif
    !++++++++++ CRM mode ++++++++++
    else
Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_constants.F90
===================================================================
--- LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_constants.F90	(revision 1410)
+++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_constants.F90	(revision 1415)
@@ -51,4 +51,10 @@
     ! Number of possible output variables
     integer,parameter :: N_OUT_LIST = 27
+    ! Value for forward model result from a level that is under the ground
+    real,parameter :: R_GROUND = -1.0E20
+
+    ! Stratiform and convective clouds in frac_out
+    integer, parameter :: I_LSC = 1, & ! Large-scale clouds
+                          I_CVC = 2    ! Convective clouds
     
     !--- Radar constants
@@ -56,5 +62,5 @@
     integer,parameter :: DBZE_BINS     =   15   ! Number of dBZe bins in histogram (cfad)
     real,parameter    :: DBZE_MIN      = -100.0 ! Minimum value for radar reflectivity
-    real,parameter    :: DBZE_MAX      =   30.0 ! Maximum value for radar reflectivity
+    real,parameter    :: DBZE_MAX      =   80.0 ! Maximum value for radar reflectivity
     real,parameter    :: CFAD_ZE_MIN   =  -50.0 ! Lower value of the first CFAD Ze bin
     real,parameter    :: CFAD_ZE_WIDTH =    5.0 ! Bin width (dBZe)
@@ -69,6 +75,6 @@
     integer,parameter :: LIDAR_NCAT    =   4
     integer,parameter :: PARASOL_NREFL =   5 ! parasol
-!     real,parameter,dimension(PARASOL_NREFL) :: PARASOL_SZA = (/0.0, 15.0, 30.0, 45.0, 60.0/)
-    real,parameter,dimension(PARASOL_NREFL) :: PARASOL_SZA = (/1.0, 2.0, 3.0, 4.0, 5.0/)
+    real,parameter,dimension(PARASOL_NREFL) :: PARASOL_SZA = (/0.0, 20.0, 40.0, 6.0, 80.0/)
+!    real,parameter,dimension(PARASOL_NREFL) :: PARASOL_SZA = (/1.0, 2.0, 3.0, 4.0, 5.0/)
     real,parameter    :: DEFAULT_LIDAR_REFF = 30.0e-6 ! Default lidar effective radius
     
Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_isccp_simulator.F90
===================================================================
--- LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_isccp_simulator.F90	(revision 1410)
+++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_isccp_simulator.F90	(revision 1415)
@@ -85,5 +85,5 @@
      
   ! Change boxptop from hPa to Pa. This avoids using UDUNITS in CMOR
-  y%boxptop = y%boxptop*100.0
+!  y%boxptop = y%boxptop*100.0
   
   ! Check if there is any value slightly greater than 1
Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_simulator.F90
===================================================================
--- LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_simulator.F90	(revision 1410)
+++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_simulator.F90	(revision 1415)
@@ -65,20 +65,10 @@
   !+++++++++ Radar model ++++++++++  
   if (cfg%Lradar_sim) then
-    call system_clock(t0,count_rate,count_max)
     call cosp_radar(gbx,sgx,sghydro,sgradar)
-    call system_clock(t1,count_rate,count_max)
-    print *, '%%%%%%  Radar:', (t1-t0)*1.0/count_rate, ' s'
-  else 
-    print *, '%%%%%%  Radar not used'
   endif
   
   !+++++++++ Lidar model ++++++++++
   if (cfg%Llidar_sim) then
-    call system_clock(t0,count_rate,count_max)
     call cosp_lidar(gbx,sgx,sghydro,sglidar)
-    call system_clock(t1,count_rate,count_max)
-    print *, '%%%%%%  Lidar:', (t1-t0)*1.0/count_rate, ' s'
-  else 
-    print *, '%%%%%%  Lidar not used'
   endif
 
@@ -86,41 +76,20 @@
   !+++++++++ ISCCP simulator ++++++++++
   if (cfg%Lisccp_sim) then
-    call system_clock(t0,count_rate,count_max)
     call cosp_isccp_simulator(gbx,sgx,isccp)
-    call system_clock(t1,count_rate,count_max)
-    print *, '%%%%%%  ISCCP:', (t1-t0)*1.0/count_rate, ' s'
-  else 
-    print *, '%%%%%%  ISCCP not used'
   endif
   
   !+++++++++ MISR simulator ++++++++++
   if (cfg%Lmisr_sim) then
-    call system_clock(t0,count_rate,count_max)
     call cosp_misr_simulator(gbx,sgx,misr)
-    call system_clock(t1,count_rate,count_max)
-    print *, '%%%%%%  MISR:', (t1-t0)*1.0/count_rate, ' s'
-  else 
-    print *, '%%%%%%  MISR not used'
   endif
   
 
   !+++++++++++ Summary statistics +++++++++++
-!   write(*,*) 'Stats:'
-!   read(*,*) c 
   if (cfg%Lstats) then
-    call system_clock(t0,count_rate,count_max)
     call cosp_stats(gbx,sgx,cfg,sgradar,sglidar,vgrid,stradar,stlidar)
-    call system_clock(t1,count_rate,count_max)
-    print *, '%%%%%%  Stats:', (t1-t0)*1.0/count_rate, ' s'
-  endif
-  !+++++++++++ change of units after computation of statistics +++++++++++
-  if (cfg%Llidar_sim) then
-    where((sglidar%beta_tot > 0.0) .and. (sglidar%beta_tot /= R_UNDEF)) 
-        sglidar%beta_tot = log10(sglidar%beta_tot)
-    elsewhere
-        sglidar%beta_tot = R_UNDEF
-    end where
+!    print *, '%%%%%%  Stats:', (t1-t0)*1.0/count_rate, ' s'
   endif
 
+  
 END SUBROUTINE COSP_SIMULATOR
 
Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_stats.F90
===================================================================
--- LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_stats.F90	(revision 1410)
+++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_stats.F90	(revision 1415)
@@ -1,24 +1,24 @@
 ! (c) British Crown Copyright 2008, the Met Office.
 ! 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 Met Office nor the names of its contributors may be used 
-!       to endorse or promote products derived from this software without specific prior written 
+!     * Neither the name of the Met Office nor the names of its 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.
 
@@ -29,7 +29,7 @@
 ! Oct 2008 - J.-L. Dufresne   - Bug fixed. Assignment of Npoints,Nlevels,Nhydro,Ncolumns in COSP_STATS
 ! Oct 2008 - H. Chepfer       - Added PARASOL reflectance arguments
-! May 2010 - L. Fairhead      - Optimisation of COSP_CHANGE_VERTICAL_GRID routine for NEC SX8 computer
-!
-! 
+! Jun 2010 - T. Yokohata, T. Nishimura and K. Ogochi - Added NEC SXs optimisations
+!
+!
 MODULE MOD_COSP_STATS
   USE MOD_COSP_CONSTANTS
@@ -45,5 +45,5 @@
 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 SUBROUTINE COSP_STATS(gbx,sgx,cfg,sgradar,sglidar,vgrid,stradar,stlidar)
-  
+
    ! Input arguments
    type(cosp_gridbox),intent(in) :: gbx
@@ -55,7 +55,7 @@
    ! Output arguments
    type(cosp_radarstats),intent(inout) :: stradar ! Summary statistics for radar
-   type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar 
-   
-   ! Local variables 
+   type(cosp_lidarstats),intent(inout) :: stlidar ! Summary statistics for lidar
+
+   ! Local variables
    integer :: Npoints  !# of grid points
    integer :: Nlevels  !# of levels
@@ -66,5 +66,5 @@
    real,dimension(:,:,:),allocatable :: Ze_out,betatot_out,betamol_in,betamol_out,ph_in,ph_out
    real,dimension(:,:),allocatable :: ph_c,betamol_c
- 
+
    Npoints  = gbx%Npoints
    Nlevels  = gbx%Nlevels
@@ -73,5 +73,5 @@
    Nlr      = vgrid%Nlvgrid
 
-   if (cfg%Lcfad_lidarsr532) ok_lidar_cfad=.true.
+   if (cfg%Lcfad_Lidarsr532) ok_lidar_cfad=.true.
 
    if (vgrid%use_vgrid) then ! Statistics in a different vertical grid
@@ -81,5 +81,4 @@
         Ze_out = 0.0
         betatot_out  = 0.0
-        betamol_in(:,1,:) = sglidar%beta_mol(:,:)
         betamol_out= 0.0
         betamol_c  = 0.0
@@ -89,5 +88,5 @@
         !++++++++++++ Radar CFAD ++++++++++++++++
         if (cfg%Lradar_sim) then
-           call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
+            call cosp_change_vertical_grid(Npoints,Ncolumns,Nlevels,gbx%zlev,gbx%zlev_half,sgradar%Ze_tot, &
                                            Nlr,vgrid%zl,vgrid%zu,Ze_out,log_units=.true.)
             stradar%cfad_ze = cosp_cfad(Npoints,Ncolumns,Nlr,DBZE_BINS,Ze_out, &
@@ -96,4 +95,5 @@
         !++++++++++++ Lidar CFAD ++++++++++++++++
         if (cfg%Llidar_sim) then
+            betamol_in(:,1,:) = sglidar%beta_mol(:,:)
             call cosp_change_vertical_grid(Npoints,1,Nlevels,gbx%zlev,gbx%zlev_half,betamol_in, &
                                            Nlr,vgrid%zl,vgrid%zu,betamol_out)
@@ -114,5 +114,5 @@
         if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
                                     betatot_out,betamol_c,Ze_out, &
-                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)   
+                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
         ! Deallocate arrays at coarse resolution
         deallocate(Ze_out,betatot_out,betamol_in,betamol_out,betamol_c,ph_in,ph_out,ph_c)
@@ -131,11 +131,11 @@
         if (cfg%Lradar_sim.and.cfg%Llidar_sim) call cosp_lidar_only_cloud(Npoints,Ncolumns,Nlr, &
                                     sglidar%beta_tot,sglidar%beta_mol,sgradar%Ze_tot, &
-                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)   
+                                    stradar%lidar_only_freq_cloud,stradar%radar_lidar_tcc)
    endif
    ! Replace undef
-   where (stlidar%cfad_sr   == LIDAR_UNDEF) stlidar%cfad_sr   = R_UNDEF 
-   where (stlidar%lidarcld  == LIDAR_UNDEF) stlidar%lidarcld  = R_UNDEF 
-   where (stlidar%cldlayer  == LIDAR_UNDEF) stlidar%cldlayer  = R_UNDEF 
-   where (stlidar%parasolrefl == LIDAR_UNDEF) stlidar%parasolrefl = R_UNDEF 
+   where (stlidar%cfad_sr   == LIDAR_UNDEF) stlidar%cfad_sr   = R_UNDEF
+   where (stlidar%lidarcld  == LIDAR_UNDEF) stlidar%lidarcld  = R_UNDEF
+   where (stlidar%cldlayer  == LIDAR_UNDEF) stlidar%cldlayer  = R_UNDEF
+   where (stlidar%parasolrefl == LIDAR_UNDEF) stlidar%parasolrefl = R_UNDEF
 
 END SUBROUTINE COSP_STATS
@@ -162,95 +162,93 @@
 
    ! Local variables
-   integer :: i,j,k,c
+   integer :: i,j,k
    logical :: lunits
-   real :: ws(Npoints,M), ws_temp(Npoints,M)
-   real,dimension(Npoints,Nlevels) :: xl, xu ! Lower and upper boundaries of model grid
-   real,dimension(M) :: dz          ! Layer depth
-   real,dimension(Npoints,Nlevels,M) :: w   ! Weights to do the mean at each point
-   real,dimension(Npoints,Ncolumns,Nlevels) :: yp  ! Variable to be changed to a different grid.
-                                           ! Local copy at a particular point.
-                                           ! This allows for change of units.
-  
+
+   integer :: l
+   real,dimension(Npoints) :: ws,sumwyp
+   real,dimension(Npoints,Nlevels) :: xl,xu
+   real,dimension(Npoints,Nlevels) :: w
+   real,dimension(Npoints,Ncolumns,Nlevels) :: yp
 
    lunits=.false.
    if (present(log_units)) lunits=log_units
 
-   r(:,:,:) = 0.0
-   yp(:,:,:) = y(:,:,:)
-   w(:,:,:) = 0.0
-   ws(:,:) = 0.0
-   ws_temp(:,:) = 0.0
-   dz = zu - zl
-
+   r(:,:,:) = R_GROUND
+   ! Vertical grid at that point
+   xl(:,:) = zhalf(:,:)
+   xu(:,1:Nlevels-1) = xl(:,2:Nlevels)
+   xu(:,Nlevels) = zfull(:,Nlevels) +  zfull(:,Nlevels) - zhalf(:,Nlevels) ! Top level symmetric
+   yp(:,:,:) = y(:,:,:) ! Temporary variable to regrid
    ! Check for dBZ and change if necessary
    if (lunits) then
      where (y /= R_UNDEF)
-        yp = 10.0**(y/10.0)
+       yp = 10.0**(y/10.0)
+     elsewhere
+       yp = 0.0
      end where
    endif
-
-   ! Vertical grids
-   xl(:,:) = zhalf(:,:)
-   xu(:,1:Nlevels-1) = zhalf(:,2:Nlevels)
-   xu(:,Nlevels) = zfull(:,Nlevels) +  zfull(:,Nlevels) - zhalf(:,Nlevels) ! Top level symmetric
-  ! Find weights
-   do k=1, M
-     do j=1, Nlevels
-       do i=1, Npoints
+   do k=1,M
+     ! Find weights
+     w(:,:) = 0.0
+     do j=1,Nlevels
+       do i=1,Npoints
          if ((xl(i,j) < zl(k)).and.(xu(i,j) > zl(k)).and.(xu(i,j) <= zu(k))) then
            !xl(j)-----------------xu(j)
            !      zl(k)------------------------------zu(k)
-           w(i,j,k) = xu(i,j) - zl(k)
+           w(i,j) = xu(i,j) - zl(k)
          else if ((xl(i,j) >= zl(k)).and.(xu(i,j) <= zu(k))) then
            !           xl(j)-----------------xu(j)
            !      zl(k)------------------------------zu(k)
-           w(i,j,k) =  xl(i,j+1) - xl(i,j)
+           w(i,j) = xu(i,j) - xl(i,j)
          else if ((xl(i,j) >= zl(k)).and.(xl(i,j) < zu(k)).and.(xu(i,j) >= zu(k))) then
            !                           xl(j)-----------------xu(j)
            !      zl(k)------------------------------zu(k)
-           w(i,j,k) = zu(k) - xl(i,j)
+           w(i,j) = zu(k) - xl(i,j)
          else if ((xl(i,j) <= zl(k)).and.(xu(i,j) >= zu(k))) then
            !  xl(j)---------------------------xu(j)
            !        zl(k)--------------zu(k)
-           w(i,j,k) = dz(j)
+           w(i,j) = zu(k) - zl(k)
+         endif
+       enddo
+     enddo
+     ! Do the weighted mean
+     do j=1,Ncolumns
+       ws    (:) = 0.0
+       sumwyp(:) = 0.0
+       do l=1,Nlevels
+         do i=1,Npoints
+           if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
+             ws    (i) = ws    (i) + w(i,l)
+             sumwyp(i) = sumwyp(i) + w(i,l)*yp(i,j,l)
+           endif
+         enddo
+       enddo
+       do i=1,Npoints
+         if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
+           if (ws(i) > 0.0) r(i,j,k) = sumwyp(i)/ws(i)
          endif
        enddo
      enddo
    enddo
-
-   ! Do the weighted mean
-   do k=1,M
-     do i = 1, Npoints
-        ws(i,k) = sum(w(i,:,k))
-     enddo
-   enddo
-
-   ws_temp = 1.
-   where (ws(:,:) > 0.0) ws_temp(:,:)=ws(:,:)
-
-   do c=1,Ncolumns
-     do k=1,M
-       do i = 1, Npoints
-          r(i,c,k) = sum(w(i,:,k)*yp(i,c,:))/ws_temp(i,k)
-       enddo
-     enddo
-   enddo
-
-   do k=1,M
-     do i = 1, Npoints
-        if (ws(i,k) <= 0.0) r(i,:,k)=0.0
-     enddo
-   enddo
- 
    ! Check for dBZ and change if necessary
    if (lunits) then
-      where (r(:,:,:) <= 0.0)
-        r(:,:,:) = R_UNDEF
-      elsewhere
-        r(:,:,:) = 10.0*log10(r(:,:,:))
-      end where
+     do k=1,M
+       do j=1,Ncolumns
+         do i=1,Npoints
+           if (zu(k) > zhalf(i,1)) then ! Level above model bottom level
+             if (r(i,j,k) <= 0.0) then
+                 r(i,j,k) = R_UNDEF
+             else
+                 r(i,j,k) = 10.0*log10(r(i,j,k))
+             endif
+           endif
+         enddo
+       enddo
+     enddo
    endif
-   
-END SUBROUTINE COSP_CHANGE_VERTICAL_GRID 
+
+
+
+END SUBROUTINE COSP_CHANGE_VERTICAL_GRID
 
 END MODULE MOD_COSP_STATS
Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_utils.F90
===================================================================
--- LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_utils.F90	(revision 1410)
+++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/cosp_utils.F90	(revision 1415)
@@ -40,4 +40,46 @@
   END INTERFACE
 CONTAINS
+
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+!------------------- SUBROUTINE COSP_PRECIP_MXRATIO --------------
+!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+SUBROUTINE COSP_PRECIP_MXRATIO(Npoints,Nlevels,Ncolumns,p,T,prec_frac,prec_type, &
+                          n_ax,n_bx,alpha_x,c_x,d_x,g_x,a_x,b_x,gamma1,gamma2, &
+                          flux,mxratio)
+
+    ! Input arguments, (IN)
+    integer,intent(in) :: Npoints,Nlevels,Ncolumns
+    real,intent(in),dimension(Npoints,Nlevels) :: p,T,flux
+    real,intent(in),dimension(Npoints,Ncolumns,Nlevels) :: prec_frac
+    real,intent(in) :: n_ax,n_bx,alpha_x,c_x,d_x,g_x,a_x,b_x,gamma1,gamma2,prec_type
+    ! Input arguments, (OUT)
+    real,intent(out),dimension(Npoints,Ncolumns,Nlevels) :: mxratio
+    ! Local variables
+    integer :: i,j,k
+    real :: sigma,one_over_xip1,xi,rho0,rho
+    
+    mxratio = 0.0
+
+    if (n_ax >= 0.0) then ! N_ax is used to control which hydrometeors need to be computed
+        !gamma1  = gamma(alpha_x + b_x + d_x + 1.0)
+        !gamma2  = gamma(alpha_x + b_x + 1.0)
+        xi      = d_x/(alpha_x + b_x - n_bx + 1.0)
+        rho0    = 1.29
+        sigma   = (gamma2/(gamma1*c_x))*(n_ax*a_x*gamma2)**xi
+        one_over_xip1 = 1.0/(xi + 1.0)
+        
+        do k=1,Nlevels
+            do j=1,Ncolumns
+                do i=1,Npoints
+                    if ((prec_frac(i,j,k)==prec_type).or.(prec_frac(i,j,k)==3.)) then
+                        rho = p(i,k)/(287.05*T(i,k))
+                        mxratio(i,j,k)=(flux(i,k)*((rho/rho0)**g_x)*sigma)**one_over_xip1
+                        mxratio(i,j,k)=mxratio(i,j,k)/rho
+                    endif
+                enddo
+            enddo
+        enddo
+    endif
+END SUBROUTINE COSP_PRECIP_MXRATIO
 
 
Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/icarus.F
===================================================================
--- LMDZ4/branches/LMDZ4_AR5/libf/cosp/icarus.F	(revision 1410)
+++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/icarus.F	(revision 1415)
@@ -33,5 +33,5 @@
      &)
 
-!Id: icarus.f,v 4.0 2009/02/12 13:59:20 hadmw Exp $
+!$Id: icarus.f,v 4.1 2010/05/27 16:30:18 hadmw Exp $
 
 ! *****************************COPYRIGHT****************************
@@ -72,4 +72,5 @@
 ! *****************************COPYRIGHT*******************************
 ! *****************************COPYRIGHT*******************************
+! *****************************COPYRIGHT*******************************
 
       implicit none
@@ -234,5 +235,4 @@
       REAL attropmin (npoints)
       REAL atmax(npoints)
-      REAL atmin(npoints)
       REAL btcmin(npoints)
       REAL transmax(npoints)
@@ -356,5 +356,4 @@
       do j=1,npoints 
           ptrop(j)=5000.
-          atmin(j) = 400.
           attropmin(j) = 400.
           atmax(j) = 0.
@@ -373,8 +372,13 @@
                 itrop(j)=ilev
            end if
-           if (at(j,ilev) .gt. atmax(j)) atmax(j)=at(j,ilev)
-           if (at(j,ilev) .lt. atmin(j)) atmin(j)=at(j,ilev)
         enddo
 12    continue
+
+      do 13 ilev=1,nlev
+        do j=1,npoints 
+           if (at(j,ilev) .gt. atmax(j) .and.
+     &              ilev  .ge. itrop(j)) atmax(j)=at(j,ilev)
+        enddo
+13    continue
 
       end if
Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/llnl_stats.F90
===================================================================
--- LMDZ4/branches/LMDZ4_AR5/libf/cosp/llnl_stats.F90	(revision 1410)
+++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/llnl_stats.F90	(revision 1415)
@@ -24,4 +24,5 @@
 
 MODULE MOD_LLNL_STATS
+  USE MOD_COSP_CONSTANTS
   IMPLICIT NONE
 
@@ -62,6 +63,8 @@
    do j = 1, Nlevels, 1
       do k = 1, Ncolumns, 1
-         do i = 1, Npoints, 1 
-            if ((x(i,k,j) >= xmin) .and. (x(i,k,j) <= xmax)) then 
+         do i = 1, Npoints, 1
+            if (x(i,k,j) == R_GROUND) then
+               cosp_cfad(i,:,j) = R_UNDEF
+            elseif ((x(i,k,j) >= xmin) .and. (x(i,k,j) <= xmax)) then 
                ibin = ceiling((x(i,k,j) - bmin)/bwidth)
                if (ibin > Nbins) ibin = Nbins
@@ -72,5 +75,5 @@
       enddo  !k
    enddo  !j
-   cosp_cfad = cosp_cfad / Ncolumns
+   where ((cosp_cfad /= R_UNDEF).and.(cosp_cfad /= 0.0)) cosp_cfad = cosp_cfad / Ncolumns
 END FUNCTION COSP_CFAD
 
@@ -98,6 +101,6 @@
    integer :: pr,i,j
    
-!    lidar_only_freq_cloud = 0.0
-!    tcc = 0.0
+   lidar_only_freq_cloud = 0.0
+   tcc = 0.0
    do pr=1,Npoints
      do i=1,Ncolumns
Index: LMDZ4/branches/LMDZ4_AR5/libf/cosp/lmd_ipsl_stats.F90
===================================================================
--- LMDZ4/branches/LMDZ4_AR5/libf/cosp/lmd_ipsl_stats.F90	(revision 1410)
+++ LMDZ4/branches/LMDZ4_AR5/libf/cosp/lmd_ipsl_stats.F90	(revision 1415)
@@ -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
