Context Navigation

lmd_ipsl_stats.F90 @ 2906

Last change on this file since 2906 was 2428, checked in by idelkadi, 9 years ago

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

mise a jour des sources pour ISCCP, CALIPSO et PARASOL
prise en compte des changements de phases pour les nuages (Calipso)
rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp

Property copyright set to
Name of program: LMDZ
Creation date: 1984
Version: LMDZ5
License: CeCILL version 2
Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
See the license file in the root directory

File size: 37.2 KB

Rev	Line
[1262]	1	! Copyright (c) 2009, Centre National de la Recherche Scientifique
	2	! All rights reserved.
[2428]	3	! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $
	4	! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/actsim/lmd_ipsl_stats.F90 $
[1414]	5	!
	6	! Redistribution and use in source and binary forms, with or without modification, are permitted
[1262]	7	! provided that the following conditions are met:
[1414]	8	!
	9	! * Redistributions of source code must retain the above copyright notice, this list
[1262]	10	! of conditions and the following disclaimer.
	11	! * Redistributions in binary form must reproduce the above copyright notice, this list
[1414]	12	! of conditions and the following disclaimer in the documentation and/or other materials
[1262]	13	! provided with the distribution.
	14	! * Neither the name of the LMD/IPSL/CNRS/UPMC nor the names of its
[1414]	15	! contributors may be used to endorse or promote products derived from this software without
[1262]	16	! specific prior written permission.
[1414]	17	!
	18	! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
	19	! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
	20	! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
	21	! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
	22	! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
	23	! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
	24	! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
[1262]	25	! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
	26
	27
	28	!------------------------------------------------------------------------------------
	29	! Authors: Sandrine Bony and Helene Chepfer (LMD/IPSL, CNRS, UPMC, France).
	30	!------------------------------------------------------------------------------------
	31	MODULE MOD_LMD_IPSL_STATS
	32	USE MOD_LLNL_STATS
	33	IMPLICIT NONE
	34
	35	CONTAINS
	36	SUBROUTINE diag_lidar(npoints,ncol,llm,max_bin,nrefl &
[2428]	37	,tmp,pnorm,pnorm_perp,pmol,refl,land,pplay,undef,ok_lidar_cfad &
	38	,cfad2,srbval,ncat,lidarcld,lidarcldphase,cldlayer,cldlayerphase &
	39	,lidarcldtmp,parasolrefl)
[1262]	40	!
	41	! -----------------------------------------------------------------------------------
	42	! Lidar outputs :
[1414]	43	!
[2428]	44	! Diagnose cloud fraction (3D cloud fraction + low/middle/high/total cloud fraction)
	45	! and phase cloud fraction (3D, low/mid/high/total and 3D temperature)
	46	! from the lidar signals (ATB, ATBperp and molecular ATB) computed from model outputs
[1262]	47	! +
	48	! Compute CFADs of lidar scattering ratio SR and of depolarization index
[1414]	49	!
[1262]	50	! Authors: Sandrine Bony and Helene Chepfer (LMD/IPSL, CNRS, UPMC, France).
	51	!
[1414]	52	! December 2008, S. Bony, H. Chepfer and J-L. Dufresne :
[1262]	53	! - change of the cloud detection threshold S_cld from 3 to 5, for better
	54	! with both day and night observations. The optical thinest clouds are missed.
	55	! - remove of the detection of the first fully attenuated layer encountered from above.
	56	! December 2008, A. Bodas-Salcedo:
	57	! - Dimensions of pmol reduced to (npoints,llm)
[1414]	58	! August 2009, A. Bodas-Salcedo:
	59	! - Warning message regarding PARASOL being valid only over ocean deleted.
	60	! February 2010, A. Bodas-Salcedo:
	61	! - Undef passed into cosp_cfad_sr
	62	! June 2010, T. Yokohata, T. Nishimura and K. Ogochi
	63	! Optimisation of COSP_CFAD_SR
[1262]	64	!
[2428]	65	! January 2013, G. Cesana, H. Chepfer:
	66	! - Add the perpendicular component of the backscattered signal (pnorm_perp) in the arguments
	67	! - Add the temperature (tmp) in the arguments
	68	! - Add the 3D Phase cloud fraction (lidarcldphase) in the arguments
	69	! - Add the Phase low mid high cloud fraction (cldlayerphase) in the arguments
	70	! - Add the 3D Phase cloud fraction as a function of temperature (lidarcldtmp) in the arguments
	71	! - Modification of the phase diagnosis within the COSP_CLDFRAC routine to integrate the phase
	72	! diagnosis (3D, low/mid/high, 3D temperature)
	73	! Reference: Cesana G. and H. Chepfer (2013): Evaluation of the cloud water phase
	74	! in a climate model using CALIPSO-GOCCP, J. Geophys. Res., doi: 10.1002/jgrd.50376
	75	!
	76	! ------------------------------------------------------------------------------------
[1262]	77
	78	! c inputs :
	79	integer npoints
	80	integer ncol
	81	integer llm
	82	integer max_bin ! nb of bins for SR CFADs
	83	integer ncat ! nb of cloud layer types (low,mid,high,total)
	84	integer nrefl ! nb of solar zenith angles for parasol reflectances
	85
	86	real undef ! undefined value
[1414]	87	real pnorm(npoints,ncol,llm) ! lidar ATB
[1262]	88	real pmol(npoints,llm) ! molecular ATB
[1414]	89	real land(npoints) ! Landmask [0 - Ocean, 1 - Land]
[1262]	90	real pplay(npoints,llm) ! pressure on model levels (Pa)
	91	logical ok_lidar_cfad ! true if lidar CFAD diagnostics need to be computed
	92	real refl(npoints,ncol,nrefl) ! subgrid parasol reflectance ! parasol
[2428]	93	real tmp(npoints,llm) ! temp at each levels
	94	real pnorm_perp(npoints,ncol,llm) ! lidar perpendicular ATB
[1262]	95
	96	! c outputs :
[1414]	97	real lidarcld(npoints,llm) ! 3D "lidar" cloud fraction
[2428]	98	real sub(npoints,llm) ! 3D "lidar" indice
	99	real cldlayer(npoints,ncat) ! "lidar" cloud layer fraction (low, mid, high, total)
	100
[1414]	101	real cfad2(npoints,max_bin,llm) ! CFADs of SR
	102	real srbval(max_bin) ! SR bins in CFADs
[1262]	103	real parasolrefl(npoints,nrefl)! grid-averaged parasol reflectance
	104
	105	! c threshold for cloud detection :
[1414]	106	real S_clr
	107	parameter (S_clr = 1.2)
[1262]	108	real S_cld
[2428]	109	parameter (S_cld = 5.) ! Thresold for cloud detection
[1262]	110	real S_att
	111	parameter (S_att = 0.01)
	112
	113	! c local variables :
[2428]	114	integer ic,k,i,j
[1262]	115	real x3d(npoints,ncol,llm)
	116	real x3d_c(npoints,llm),pnorm_c(npoints,llm)
	117	real xmax
[2428]	118
	119	! Output variables
	120	integer,parameter :: nphase = 6 ! nb of cloud layer phase types (ice,liquid,undefined,false ice,false liquid,Percent of ice)
	121	real lidarcldphase(npoints,llm,nphase) ! 3D "lidar" phase cloud fraction
	122	real lidarcldtmp(npoints,40,5) ! 3D "lidar" phase cloud fraction as a function of temp
	123	real cldlayerphase(npoints,ncat,nphase) ! "lidar" phase low mid high cloud fraction
	124
	125	! SR detection threshold
	126	real, parameter :: S_cld_att = 30. ! New threshold for undefine cloud phase detection
	127
	128
[1262]	129	!
	130	! c -------------------------------------------------------
	131	! c 0- Initializations
	132	! c -------------------------------------------------------
[1414]	133	!
[1262]	134	! Should be modified in future version
	135	xmax=undef-1.0
	136
	137	! c -------------------------------------------------------
	138	! c 1- Lidar scattering ratio :
	139	! c -------------------------------------------------------
[2428]	140
[1262]	141	do ic = 1, ncol
	142	pnorm_c = pnorm(:,ic,:)
[1414]	143	where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
[1262]	144	x3d_c = pnorm_c/pmol
	145	elsewhere
	146	x3d_c = undef
	147	end where
[2428]	148	x3d(:,ic,:) = x3d_c
[1262]	149	enddo
	150
	151	! c -------------------------------------------------------
	152	! c 2- Diagnose cloud fractions (3D, low, middle, high, total)
	153	! c from subgrid-scale lidar scattering ratios :
	154	! c -------------------------------------------------------
	155
[2428]	156	CALL COSP_CLDFRAC(npoints,ncol,llm,ncat,nphase, &
	157	tmp,x3d,pnorm,pnorm_perp,pplay, S_att,S_cld,S_cld_att,undef,lidarcld, &
	158	cldlayer,lidarcldphase,sub,cldlayerphase,lidarcldtmp)
[1262]	159
	160	! c -------------------------------------------------------
[1414]	161	! c 3- CFADs
[1262]	162	! c -------------------------------------------------------
	163	if (ok_lidar_cfad) then
	164	!
	165	! c CFADs of subgrid-scale lidar scattering ratios :
	166	! c -------------------------------------------------------
[1414]	167	CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin,undef, &
[1262]	168	x3d, &
	169	S_att,S_clr,xmax,cfad2,srbval)
	170
	171	endif ! ok_lidar_cfad
	172	! c -------------------------------------------------------
	173
	174	! c -------------------------------------------------------
	175	! c 4- Compute grid-box averaged Parasol reflectances
	176	! c -------------------------------------------------------
	177
	178	parasolrefl(:,:) = 0.0
	179
	180	do k = 1, nrefl
	181	do ic = 1, ncol
	182	parasolrefl(:,k) = parasolrefl(:,k) + refl(:,ic,k)
	183	enddo
	184	enddo
	185
	186	do k = 1, nrefl
	187	parasolrefl(:,k) = parasolrefl(:,k) / float(ncol)
	188	! if land=1 -> parasolrefl=undef
	189	! if land=0 -> parasolrefl=parasolrefl
	190	parasolrefl(:,k) = parasolrefl(:,k) * MAX(1.0-land(:),0.0) &
[1414]	191	+ (1.0 - MAX(1.0-land(:),0.0))*undef
[1262]	192	enddo
	193
	194	RETURN
	195	END SUBROUTINE diag_lidar
[1414]	196
	197
[1262]	198	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	199	!-------------------- FUNCTION COSP_CFAD_SR ------------------------
	200	! Author: Sandrine Bony (LMD/IPSL, CNRS, Paris)
	201	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[1414]	202	SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins,undef, &
[1262]	203	x,S_att,S_clr,xmax,cfad,srbval)
	204	IMPLICIT NONE
	205
	206	!--- Input arguments
	207	! Npoints: Number of horizontal points
	208	! Ncolumns: Number of subcolumns
	209	! Nlevels: Number of levels
	210	! Nbins: Number of x axis bins
	211	! xmax: maximum value allowed for x
	212	! S_att: Threshold for full attenuation
	213	! S_clr: Threshold for clear-sky layer
	214	!
	215	!--- Input-Outout arguments
	216	! x: variable to process (Npoints,Ncolumns,Nlevels), mofified where saturation occurs
	217	!
	218	! -- Output arguments
	219	! srbval : values of the histogram bins
	220	! cfad: 2D histogram on each horizontal point
	221
	222	! Input arguments
	223	integer Npoints,Ncolumns,Nlevels,Nbins
[1414]	224	real xmax,S_att,S_clr,undef
[1262]	225	! Input-output arguments
	226	real x(Npoints,Ncolumns,Nlevels)
	227	! Output :
	228	real cfad(Npoints,Nbins,Nlevels)
	229	real srbval(Nbins)
	230	! Local variables
	231	integer i, j, k, ib
[1414]	232	real srbval_ext(0:Nbins)
[1262]	233
	234	! c -------------------------------------------------------
	235	! c 0- Initializations
	236	! c -------------------------------------------------------
	237	if ( Nbins .lt. 6) return
	238
	239	srbval(1) = S_att
	240	srbval(2) = S_clr
	241	srbval(3) = 3.0
	242	srbval(4) = 5.0
	243	srbval(5) = 7.0
	244	srbval(6) = 10.0
	245	do i = 7, MIN(10,Nbins)
	246	srbval(i) = srbval(i-1) + 5.0
	247	enddo
	248	DO i = 11, MIN(13,Nbins)
	249	srbval(i) = srbval(i-1) + 10.0
	250	enddo
	251	srbval(MIN(14,Nbins)) = 80.0
	252	srbval(Nbins) = xmax
	253	cfad(:,:,:) = 0.0
	254
[1414]	255	srbval_ext(1:Nbins) = srbval
	256	srbval_ext(0) = -1.0
[1262]	257	! c -------------------------------------------------------
	258	! c c- Compute CFAD
	259	! c -------------------------------------------------------
[1414]	260	do j = 1, Nlevels
	261	do ib = 1, Nbins
	262	do k = 1, Ncolumns
	263	do i = 1, Npoints
	264	if (x(i,k,j) /= undef) then
	265	if ((x(i,k,j).gt.srbval_ext(ib-1)).and.(x(i,k,j).le.srbval_ext(ib))) &
	266	cfad(i,ib,j) = cfad(i,ib,j) + 1.0
	267	else
	268	cfad(i,ib,j) = undef
	269	endif
	270	enddo
	271	enddo
	272	enddo
	273	enddo
[1262]	274
[1414]	275	where (cfad .ne. undef) cfad = cfad / float(Ncolumns)
[1262]	276
	277	! c -------------------------------------------------------
	278	RETURN
	279	END SUBROUTINE COSP_CFAD_SR
	280
[2428]	281
[1262]	282	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	283	!-------------------- SUBROUTINE COSP_CLDFRAC -------------------
[1414]	284	! c Purpose: Cloud fraction diagnosed from lidar measurements
[1262]	285	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[2428]	286	SUBROUTINE COSP_CLDFRAC(Npoints,Ncolumns,Nlevels,Ncat,Nphase, &
	287	tmp,x,ATB,ATBperp,pplay,S_att,S_cld,S_cld_att,undef,lidarcld, &
	288	cldlayer,lidarcldphase,nsub,cldlayerphase,lidarcldtemp)
	289
	290
[1262]	291	IMPLICIT NONE
	292	! Input arguments
	293	integer Npoints,Ncolumns,Nlevels,Ncat
	294	real x(Npoints,Ncolumns,Nlevels)
[2428]	295
	296
	297	! Local parameters
	298	integer nphase ! nb of cloud layer phase types
	299	! (ice,liquid,undefined,false ice,false liquid,Percent of ice)
	300	integer,parameter :: Ntemp=40 ! indice of the temperature vector
	301	integer ip, k, iz, ic, ncol, nlev, i, itemp ! loop indice
	302	real S_cld_att ! New threshold for undefine cloud phase detection (SR=30)
	303	integer toplvlsat ! level of the first cloud with SR>30
	304	real alpha50, beta50, gamma50, delta50, epsilon50, zeta50 ! Polynomial Coef of the phase
	305	! discrimination line
	306
	307	! Input variables
	308	real tmp(Npoints,Nlevels) ! temperature
	309	real ATB(Npoints,Ncolumns,Nlevels) ! 3D Attenuated backscatter
	310	real ATBperp(Npoints,Ncolumns,Nlevels) ! 3D perpendicular attenuated backscatter
[1262]	311	real pplay(Npoints,Nlevels)
	312	real S_att,S_cld
	313	real undef
[2428]	314
	315	! Output variables
	316	real lidarcldtemp(Npoints,Ntemp,5) ! 3D Temperature 1=tot,2=ice,3=liq,4=undef,5=ice/ice+liq
	317	real tempmod(Ntemp+1) ! temperature bins
	318	real lidarcldphase(Npoints,Nlevels,Nphase) ! 3D cloud phase fraction
	319	real cldlayerphase(Npoints,Ncat,Nphase) ! low, middle, high, total cloud fractions for ice liquid and undefine phase
[1262]	320	real lidarcld(Npoints,Nlevels) ! 3D cloud fraction
	321	real cldlayer(Npoints,Ncat) ! low, middle, high, total cloud fractions
[2428]	322
[1262]	323	! Local variables
[2428]	324	real tmpi(Npoints,Ncolumns,Nlevels) ! temperature of ice cld
	325	real tmpl(Npoints,Ncolumns,Nlevels) ! temperature of liquid cld
	326	real tmpu(Npoints,Ncolumns,Nlevels) ! temperature of undef cld
	327
	328	real checktemp, ATBperp_tmp ! temporary variable
	329	real checkcldlayerphase, checkcldlayerphase2 ! temporary variable
	330	real sumlidarcldtemp(Npoints,Ntemp) ! temporary variable
	331
	332	real cldlayphase(Npoints,Ncolumns,Ncat,Nphase) ! subgrided low mid high phase cloud fraction
	333	real cldlayerphasetmp(Npoints,Ncat) ! temporary variable
	334	real cldlayerphasesum(Npoints,Ncat) ! temporary variable
	335	real lidarcldtempind(Npoints,Ntemp) ! 3D Temperature indice
	336	real lidarcldphasetmp(Npoints,Nlevels) ! 3D sum of ice and liquid cloud occurences
	337
	338
	339	! Local variables
[1262]	340	real p1
	341	real cldy(Npoints,Ncolumns,Nlevels)
	342	real srok(Npoints,Ncolumns,Nlevels)
	343	real cldlay(Npoints,Ncolumns,Ncat)
	344	real nsublay(Npoints,Ncolumns,Ncat), nsublayer(Npoints,Ncat)
	345	real nsub(Npoints,Nlevels)
	346
[2428]	347	#ifdef SYS_SX
[1414]	348	real cldlay1(Npoints,Ncolumns)
	349	real cldlay2(Npoints,Ncolumns)
	350	real cldlay3(Npoints,Ncolumns)
	351	real nsublay1(Npoints,Ncolumns)
	352	real nsublay2(Npoints,Ncolumns)
	353	real nsublay3(Npoints,Ncolumns)
[2428]	354	#endif
[1262]	355
[1414]	356
[2428]	357
	358
[1262]	359	! ---------------------------------------------------------------
[1414]	360	! 1- initialization
[1262]	361	! ---------------------------------------------------------------
	362
	363	if ( Ncat .ne. 4 ) then
	364	print *,'Error in lmd_ipsl_stats.cosp_cldfrac, Ncat must be 4, not',Ncat
	365	stop
	366	endif
	367
	368	lidarcld = 0.0
	369	nsub = 0.0
	370	cldlay = 0.0
	371	nsublay = 0.0
	372
[2428]	373	ATBperp_tmp = 0.
	374	lidarcldphase(:,:,:) = 0.
	375	cldlayphase(:,:,:,:) = 0.
	376	cldlayerphase(:,:,:) = 0.
	377	tmpi(:,:,:) = 0.
	378	tmpl(:,:,:) = 0.
	379	tmpu(:,:,:) = 0.
	380	cldlayerphasesum(:,:) = 0.
	381	lidarcldtemp(:,:,:) = 0.
	382	lidarcldtempind(:,:) = 0.
	383	sumlidarcldtemp(:,:) = 0.
	384	toplvlsat=0
	385	lidarcldphasetmp(:,:) = 0.
	386
	387	! temperature bins
	388	tempmod=(/-273.15,-90.,-87.,-84.,-81.,-78.,-75.,-72.,-69.,-66.,-63.,-60.,-57., &
	389	-54.,-51.,-48.,-45.,-42.,-39.,-36.,-33.,-30.,-27.,-24.,-21.,-18., &
	390	-15.,-12.,-9.,-6.,-3.,0.,3.,6.,9.,12.,15.,18.,21.,24.,200. /)
	391
	392	! convert C to K
	393	tempmod=tempmod+273.15
	394
	395	! Polynomial coefficient of the phase discrimination line used to separate liquid from ice
	396	! (Cesana and Chepfer, JGR, 2013)
	397	! ATBperp = ATB^5alpha50 + ATB^4beta50 + ATB^3gamma50 + ATB^2delta50 + ATB*epsilon50 + zeta50
	398	alpha50 = 9.0322e+15
	399	beta50 = -2.1358e+12
	400	gamma50 = 173.3963e06
	401	delta50 = -3.9514e03
	402	epsilon50 = 0.2559
	403	zeta50 = -9.4776e-07
	404
	405
[1262]	406	! ---------------------------------------------------------------
	407	! 2- Cloud detection
	408	! ---------------------------------------------------------------
	409
	410	do k = 1, Nlevels
	411
	412	! cloud detection at subgrid-scale:
	413	where ( (x(:,:,k).gt.S_cld) .and. (x(:,:,k).ne. undef) )
	414	cldy(:,:,k)=1.0
	415	elsewhere
	416	cldy(:,:,k)=0.0
	417	endwhere
	418
	419	! number of usefull sub-columns:
[1414]	420	where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef) )
[1262]	421	srok(:,:,k)=1.0
	422	elsewhere
	423	srok(:,:,k)=0.0
	424	endwhere
	425
	426	enddo ! k
	427
[2428]	428
[1262]	429	! ---------------------------------------------------------------
	430	! 3- grid-box 3D cloud fraction and layered cloud fractions (ISCCP pressure
	431	! categories) :
	432	! ---------------------------------------------------------------
[1414]	433	lidarcld = 0.0
	434	nsub = 0.0
[2428]	435	#ifdef SYS_SX
[1414]	436	!! XXX: Use cldlay[1-3] and nsublay[1-3] to avoid bank-conflicts.
	437	cldlay1 = 0.0
	438	cldlay2 = 0.0
	439	cldlay3 = 0.0
	440	cldlay(:,:,4) = 0.0 !! XXX: Ncat == 4
	441	nsublay1 = 0.0
	442	nsublay2 = 0.0
	443	nsublay3 = 0.0
	444	nsublay(:,:,4) = 0.0
[2428]	445
[1262]	446	do k = Nlevels, 1, -1
	447	do ic = 1, Ncolumns
	448	do ip = 1, Npoints
[2428]	449
	450	if(srok(ip,ic,k).gt.0.)then
	451	! Computation of the cloud fraction as a function of the temperature
	452	! instead of height, for ice,liquid and all clouds
	453	do itemp=1,Ntemp
	454	if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
	455	lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1.
	456	endif
	457	enddo
	458	endif
	459
	460	if (cldy(ip,ic,k).eq.1.) then
	461	do itemp=1,Ntemp
	462	if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
	463	lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1.
	464	endif
	465	enddo
	466	endif
	467
[1414]	468	p1 = pplay(ip,k)
	469
	470	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
	471	cldlay3(ip,ic) = MAX(cldlay3(ip,ic), cldy(ip,ic,k))
	472	nsublay3(ip,ic) = MAX(nsublay3(ip,ic), srok(ip,ic,k))
	473	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid clouds
	474	cldlay2(ip,ic) = MAX(cldlay2(ip,ic), cldy(ip,ic,k))
	475	nsublay2(ip,ic) = MAX(nsublay2(ip,ic), srok(ip,ic,k))
	476	else
	477	cldlay1(ip,ic) = MAX(cldlay1(ip,ic), cldy(ip,ic,k))
	478	nsublay1(ip,ic) = MAX(nsublay1(ip,ic), srok(ip,ic,k))
[1262]	479	endif
	480
[1414]	481	cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4), cldy(ip,ic,k))
[1262]	482	lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k)
	483	nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
	484	nsub(ip,k)=nsub(ip,k) + srok(ip,ic,k)
	485	enddo
	486	enddo
	487	enddo
[1414]	488	cldlay(:,:,1) = cldlay1
	489	cldlay(:,:,2) = cldlay2
	490	cldlay(:,:,3) = cldlay3
	491	nsublay(:,:,1) = nsublay1
	492	nsublay(:,:,2) = nsublay2
	493	nsublay(:,:,3) = nsublay3
[2428]	494	#else
	495	cldlay = 0.0
	496	nsublay = 0.0
	497	do k = Nlevels, 1, -1
	498	do ic = 1, Ncolumns
	499	do ip = 1, Npoints
[1262]	500
[2428]	501	! Computation of the cloud fraction as a function of the temperature
	502	! instead of height, for ice,liquid and all clouds
	503	if(srok(ip,ic,k).gt.0.)then
	504	do itemp=1,Ntemp
	505	if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
	506	lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1.
	507	endif
	508	enddo
	509	endif
	510
	511	if(cldy(ip,ic,k).eq.1.)then
	512	do itemp=1,Ntemp
	513	if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
	514	lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1.
	515	endif
	516	enddo
	517	endif
	518	!
	519
	520	iz=1
	521	p1 = pplay(ip,k)
	522	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
	523	iz=3
	524	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid clouds
	525	iz=2
	526	endif
	527
	528	cldlay(ip,ic,iz) = MAX(cldlay(ip,ic,iz),cldy(ip,ic,k))
	529	cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4),cldy(ip,ic,k))
	530	lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k)
	531
	532	nsublay(ip,ic,iz) = MAX(nsublay(ip,ic,iz),srok(ip,ic,k))
	533	nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
	534	nsub(ip,k)=nsub(ip,k) + srok(ip,ic,k)
	535
	536	enddo
	537	enddo
	538	enddo
	539	#endif
	540
	541
[1262]	542	! -- grid-box 3D cloud fraction
	543
	544	where ( nsub(:,:).gt.0.0 )
	545	lidarcld(:,:) = lidarcld(:,:)/nsub(:,:)
	546	elsewhere
	547	lidarcld(:,:) = undef
	548	endwhere
	549
	550	! -- layered cloud fractions
	551
	552	cldlayer = 0.0
	553	nsublayer = 0.0
	554
	555	do iz = 1, Ncat
	556	do ic = 1, Ncolumns
	557
[1414]	558	cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz)
	559	nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz)
[1262]	560
	561	enddo
	562	enddo
	563	where ( nsublayer(:,:).gt.0.0 )
	564	cldlayer(:,:) = cldlayer(:,:)/nsublayer(:,:)
	565	elsewhere
	566	cldlayer(:,:) = undef
	567	endwhere
	568
[2428]	569	! ---------------------------------------------------------------
	570	! 4- grid-box 3D cloud Phase :
	571	! ---------------------------------------------------------------
	572	! ---------------------------------------------------------------
	573	! 4.1 - For Cloudy pixels with 8.16km < z < 19.2km
	574	! ---------------------------------------------------------------
	575	do ncol=1,Ncolumns
	576	do i=1,Npoints
	577
	578	do nlev=Nlevels,18,-1 ! from 19.2km until 8.16km
	579	p1 = pplay(i,nlev)
	580
	581
	582	! Avoid zero values
	583	if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
	584	! Computation of the ATBperp along the phase discrimination line
	585	ATBperp_tmp = (ATB(i,ncol,nlev)*5)alpha50 + (ATB(i,ncol,nlev)*4)beta50 + &
	586	(ATB(i,ncol,nlev)*3)gamma50 + (ATB(i,ncol,nlev)*2)delta50 + &
	587	ATB(i,ncol,nlev)*epsilon50 + zeta50
	588
	589	!____________________________________________________________________________________________________
	590	!
	591	!4.1.a Ice: ATBperp above the phase discrimination line
	592	!____________________________________________________________________________________________________
	593	!
	594	if( (ATBperp(i,ncol,nlev)-ATBperp_tmp).ge.0. )then ! Ice clouds
	595	! ICE with temperature above 273,15°K = Liquid (false ice)
	596	if(tmp(i,nlev).gt.273.15)then ! Temperature above 273,15 K
	597	! Liquid: False ice corrected by the temperature to Liquid
	598	lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1. ! false ice detection ==> added to Liquid
	599	tmpl(i,ncol,nlev)=tmp(i,nlev)
	600	lidarcldphase(i,nlev,5)=lidarcldphase(i,nlev,5)+1. ! keep the information "temperature criterium used"
	601	! to classify the phase cloud
	602	cldlayphase(i,ncol,4,2) = 1. ! tot cloud
	603	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	604	cldlayphase(i,ncol,3,2) = 1.
	605	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	606	cldlayphase(i,ncol,2,2) = 1.
	607	else ! low cloud
	608	cldlayphase(i,ncol,1,2) = 1.
	609	endif
	610	cldlayphase(i,ncol,4,5) = 1. ! tot cloud
	611	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	612	cldlayphase(i,ncol,3,5) = 1.
	613	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	614	cldlayphase(i,ncol,2,5) = 1.
	615	else ! low cloud
	616	cldlayphase(i,ncol,1,5) = 1.
	617	endif
	618
	619	else
	620	! ICE with temperature below 273,15°K
	621	lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1.
	622	tmpi(i,ncol,nlev)=tmp(i,nlev)
	623	cldlayphase(i,ncol,4,1) = 1. ! tot cloud
	624	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	625	cldlayphase(i,ncol,3,1) = 1.
	626	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	627	cldlayphase(i,ncol,2,1) = 1.
	628	else ! low cloud
	629	cldlayphase(i,ncol,1,1) = 1.
	630	endif
	631
	632	endif
	633
	634	!____________________________________________________________________________________________________
	635	!
	636	! 4.1.b Liquid: ATBperp below the phase discrimination line
	637	!____________________________________________________________________________________________________
	638	!
	639	else ! Liquid clouds
	640	! Liquid with temperature above 231,15°K
	641	if(tmp(i,nlev).gt.231.15)then
	642	lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.
	643	tmpl(i,ncol,nlev)=tmp(i,nlev)
	644	cldlayphase(i,ncol,4,2) = 1. ! tot cloud
	645	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	646	cldlayphase(i,ncol,3,2) = 1.
	647	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	648	cldlayphase(i,ncol,2,2) = 1.
	649	else ! low cloud
	650	cldlayphase(i,ncol,1,2) = 1.
	651	endif
	652
	653	else
	654	! Liquid with temperature below 231,15°K = Ice (false liquid)
	655	tmpi(i,ncol,nlev)=tmp(i,nlev)
	656	lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1. ! false liquid detection ==> added to ice
	657	lidarcldphase(i,nlev,4)=lidarcldphase(i,nlev,4)+1. ! keep the information "temperature criterium used"
	658	! to classify the phase cloud
	659	cldlayphase(i,ncol,4,4) = 1. ! tot cloud
	660	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	661	cldlayphase(i,ncol,3,4) = 1.
	662	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	663	cldlayphase(i,ncol,2,4) = 1.
	664	else ! low cloud
	665	cldlayphase(i,ncol,1,4) = 1.
	666	endif
	667	cldlayphase(i,ncol,4,1) = 1. ! tot cloud
	668	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	669	cldlayphase(i,ncol,3,1) = 1.
	670	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	671	cldlayphase(i,ncol,2,1) = 1.
	672	else ! low cloud
	673	cldlayphase(i,ncol,1,1) = 1.
	674	endif
	675
	676	endif
	677
	678	endif ! end of discrimination condition
	679	endif ! end of cloud condition
	680	enddo ! end of altitude loop
	681
	682
	683
	684	! ---------------------------------------------------------------
	685	! 4.2 - For Cloudy pixels with 0km < z < 8.16km
	686	! ---------------------------------------------------------------
	687
	688	toplvlsat=0
	689	do nlev=17,1,-1 ! from 8.16km until 0km
	690	p1 = pplay(i,nlev)
	691
	692	if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
	693	! Phase discrimination line : ATBperp = ATB^5alpha50 + ATB^4beta50 + ATB^3gamma50 + ATB^2delta50
	694	! + ATB*epsilon50 + zeta50
	695	! Computation of the ATBperp of the phase discrimination line
	696	ATBperp_tmp = (ATB(i,ncol,nlev)*5)alpha50 + (ATB(i,ncol,nlev)*4)beta50 + &
	697	(ATB(i,ncol,nlev)*3)gamma50 + (ATB(i,ncol,nlev)*2)delta50 + &
	698	ATB(i,ncol,nlev)*epsilon50 + zeta50
	699	!____________________________________________________________________________________________________
	700	!
	701	! 4.2.a Ice: ATBperp above the phase discrimination line
	702	!____________________________________________________________________________________________________
	703	!
	704	! ICE with temperature above 273,15°K = Liquid (false ice)
	705	if( (ATBperp(i,ncol,nlev)-ATBperp_tmp).ge.0. )then ! Ice clouds
	706	if(tmp(i,nlev).gt.273.15)then
	707	lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1. ! false ice ==> liq
	708	tmpl(i,ncol,nlev)=tmp(i,nlev)
	709	lidarcldphase(i,nlev,5)=lidarcldphase(i,nlev,5)+1.
	710
	711	cldlayphase(i,ncol,4,2) = 1. ! tot cloud
	712	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	713	cldlayphase(i,ncol,3,2) = 1.
	714	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	715	cldlayphase(i,ncol,2,2) = 1.
	716	else ! low cloud
	717	cldlayphase(i,ncol,1,2) = 1.
	718	endif
	719
	720	cldlayphase(i,ncol,4,5) = 1. ! tot cloud
	721	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	722	cldlayphase(i,ncol,3,5) = 1.
	723	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	724	cldlayphase(i,ncol,2,5) = 1.
	725	else ! low cloud
	726	cldlayphase(i,ncol,1,5) = 1.
	727	endif
	728
	729	else
	730	! ICE with temperature below 273,15°K
	731	lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1.
	732	tmpi(i,ncol,nlev)=tmp(i,nlev)
	733
	734	cldlayphase(i,ncol,4,1) = 1. ! tot cloud
	735	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	736	cldlayphase(i,ncol,3,1) = 1.
	737	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	738	cldlayphase(i,ncol,2,1) = 1.
	739	else ! low cloud
	740	cldlayphase(i,ncol,1,1) = 1.
	741	endif
	742
	743	endif
	744
	745	!____________________________________________________________________________________________________
	746	!
	747	! 4.2.b Liquid: ATBperp below the phase discrimination line
	748	!____________________________________________________________________________________________________
	749	!
	750	else
	751	! Liquid with temperature above 231,15°K
	752	if(tmp(i,nlev).gt.231.15)then
	753	lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.
	754	tmpl(i,ncol,nlev)=tmp(i,nlev)
	755
	756	cldlayphase(i,ncol,4,2) = 1. ! tot cloud
	757	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	758	cldlayphase(i,ncol,3,2) = 1.
	759	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	760	cldlayphase(i,ncol,2,2) = 1.
	761	else ! low cloud
	762	cldlayphase(i,ncol,1,2) = 1.
	763	endif
	764
	765	else
	766	! Liquid with temperature below 231,15°K = Ice (false liquid)
	767	tmpi(i,ncol,nlev)=tmp(i,nlev)
	768	lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1. ! false liq ==> ice
	769	lidarcldphase(i,nlev,4)=lidarcldphase(i,nlev,4)+1. ! false liq ==> ice
	770
	771	cldlayphase(i,ncol,4,4) = 1. ! tot cloud
	772	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	773	cldlayphase(i,ncol,3,4) = 1.
	774	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	775	cldlayphase(i,ncol,2,4) = 1.
	776	else ! low cloud
	777	cldlayphase(i,ncol,1,4) = 1.
	778	endif
	779
	780	cldlayphase(i,ncol,4,1) = 1. ! tot cloud
	781	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	782	cldlayphase(i,ncol,3,1) = 1.
	783	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	784	cldlayphase(i,ncol,2,1) = 1.
	785	else ! low cloud
	786	cldlayphase(i,ncol,1,1) = 1.
	787	endif
	788
	789	endif
	790	endif ! end of discrimination condition
	791
	792	toplvlsat=0
	793
	794	! Find the level of the highest cloud with SR>30
	795	if(x(i,ncol,nlev).gt.S_cld_att)then ! SR > 30.
	796	toplvlsat=nlev-1
	797	goto 99
	798	endif
	799
	800	endif ! end of cloud condition
	801	enddo ! end of altitude loop
	802
	803	99 continue
	804
	805	!____________________________________________________________________________________________________
	806	!
	807	! Undefined phase: For a cloud located below another cloud with SR>30
	808	! see Cesana and Chepfer 2013 Sect.III.2
	809	!____________________________________________________________________________________________________
	810	!
	811	if(toplvlsat.ne.0)then
	812	do nlev=toplvlsat,1,-1
	813	p1 = pplay(i,nlev)
	814	if(cldy(i,ncol,nlev).eq.1.)then
	815	lidarcldphase(i,nlev,3)=lidarcldphase(i,nlev,3)+1.
	816	tmpu(i,ncol,nlev)=tmp(i,nlev)
	817
	818	cldlayphase(i,ncol,4,3) = 1. ! tot cloud
	819	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	820	cldlayphase(i,ncol,3,3) = 1.
	821	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	822	cldlayphase(i,ncol,2,3) = 1.
	823	else ! low cloud
	824	cldlayphase(i,ncol,1,3) = 1.
	825	endif
	826
	827	endif
	828	enddo
	829	endif
	830
	831	toplvlsat=0
	832
	833	enddo
	834	enddo
	835
	836
	837
	838	!____________________________________________________________________________________________________
	839	!
	840	! Computation of final cloud phase diagnosis
	841	!____________________________________________________________________________________________________
	842	!
	843
	844	! Compute the Ice percentage in cloud = ice/(ice+liq) as a function
	845	! of the occurrences
	846	lidarcldphasetmp(:,:)=lidarcldphase(:,:,1)+lidarcldphase(:,:,2);
	847	WHERE (lidarcldphasetmp(:,:).gt. 0.)
	848	lidarcldphase(:,:,6)=lidarcldphase(:,:,1)/lidarcldphasetmp(:,:)
	849	ELSEWHERE
	850	lidarcldphase(:,:,6) = undef
	851	ENDWHERE
	852
	853	! Compute Phase 3D Cloud Fraction
	854	WHERE ( nsub(:,:).gt.0.0 )
	855	lidarcldphase(:,:,1)=lidarcldphase(:,:,1)/nsub(:,:)
	856	lidarcldphase(:,:,2)=lidarcldphase(:,:,2)/nsub(:,:)
	857	lidarcldphase(:,:,3)=lidarcldphase(:,:,3)/nsub(:,:)
	858	lidarcldphase(:,:,4)=lidarcldphase(:,:,4)/nsub(:,:)
	859	lidarcldphase(:,:,5)=lidarcldphase(:,:,5)/nsub(:,:)
	860	ELSEWHERE
	861	lidarcldphase(:,:,1) = undef
	862	lidarcldphase(:,:,2) = undef
	863	lidarcldphase(:,:,3) = undef
	864	lidarcldphase(:,:,4) = undef
	865	lidarcldphase(:,:,5) = undef
	866	ENDWHERE
	867
	868
	869	! Compute Phase low mid high cloud fractions
	870	do iz = 1, Ncat
	871	do i=1,Nphase-3
	872	do ic = 1, Ncolumns
	873	cldlayerphase(:,iz,i)=cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i)
	874	cldlayerphasesum(:,iz)=cldlayerphasesum(:,iz)+cldlayphase(:,ic,iz,i)
	875	enddo
	876	enddo
	877	enddo
	878
	879	do iz = 1, Ncat
	880	do i=4,5
	881	do ic = 1, Ncolumns
	882	cldlayerphase(:,iz,i)=cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i)
	883	enddo
	884	enddo
	885	enddo
	886
	887	! Compute the Ice percentage in cloud = ice/(ice+liq)
	888	cldlayerphasetmp(:,:)=cldlayerphase(:,:,1)+cldlayerphase(:,:,2)
	889	WHERE (cldlayerphasetmp(:,:).gt. 0.)
	890	cldlayerphase(:,:,6)=cldlayerphase(:,:,1)/cldlayerphasetmp(:,:)
	891	ELSEWHERE
	892	cldlayerphase(:,:,6) = undef
	893	ENDWHERE
	894
	895	do i=1,Nphase-1
	896	WHERE ( cldlayerphasesum(:,:).gt.0.0 )
	897	cldlayerphase(:,:,i) = (cldlayerphase(:,:,i)/cldlayerphasesum(:,:)) * cldlayer(:,:)
	898	ENDWHERE
	899	enddo
	900
	901
	902	do i=1,Npoints
	903	do iz=1,Ncat
	904	checkcldlayerphase=0.
	905	checkcldlayerphase2=0.
	906
	907	if (cldlayerphasesum(i,iz).gt.0.0 )then
	908	do ic=1,Nphase-3
	909	checkcldlayerphase=checkcldlayerphase+cldlayerphase(i,iz,ic)
	910	enddo
	911	checkcldlayerphase2=cldlayer(i,iz)-checkcldlayerphase
	912	if( (checkcldlayerphase2.gt.0.01).or.(checkcldlayerphase2.lt.-0.01) ) print *, checkcldlayerphase,cldlayer(i,iz)
	913
	914	endif
	915
	916	enddo
	917	enddo
	918
	919	do i=1,Nphase-1
	920	WHERE ( nsublayer(:,:).eq.0.0 )
	921	cldlayerphase(:,:,i) = undef
	922	ENDWHERE
	923	enddo
	924
	925
	926
	927	! Compute Phase 3D as a function of temperature
	928	do nlev=1,Nlevels
	929	do ncol=1,Ncolumns
	930	do i=1,Npoints
	931	do itemp=1,Ntemp
	932	if(tmpi(i,ncol,nlev).gt.0.)then
	933	if( (tmpi(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpi(i,ncol,nlev).lt.tempmod(itemp+1)) )then
	934	lidarcldtemp(i,itemp,2)=lidarcldtemp(i,itemp,2)+1.
	935	endif
	936	elseif(tmpl(i,ncol,nlev).gt.0.)then
	937	if( (tmpl(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpl(i,ncol,nlev).lt.tempmod(itemp+1)) )then
	938	lidarcldtemp(i,itemp,3)=lidarcldtemp(i,itemp,3)+1.
	939	endif
	940	elseif(tmpu(i,ncol,nlev).gt.0.)then
	941	if( (tmpu(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpu(i,ncol,nlev).lt.tempmod(itemp+1)) )then
	942	lidarcldtemp(i,itemp,4)=lidarcldtemp(i,itemp,4)+1.
	943	endif
	944	endif
	945	enddo
	946	enddo
	947	enddo
	948	enddo
	949
	950	! Check temperature cloud fraction
	951	do i=1,Npoints
	952	do itemp=1,Ntemp
	953	checktemp=lidarcldtemp(i,itemp,2)+lidarcldtemp(i,itemp,3)+lidarcldtemp(i,itemp,4)
	954
	955	if(checktemp.NE.lidarcldtemp(i,itemp,1))then
	956	print *, i,itemp
	957	print *, lidarcldtemp(i,itemp,1:4)
	958	endif
	959
	960	enddo
	961	enddo
	962
	963	! Compute the Ice percentage in cloud = ice/(ice+liq)
	964	! sumlidarcldtemp=sum(lidarcldtemp(:,:,2:3),3)
	965	sumlidarcldtemp(:,:)=lidarcldtemp(:,:,2)+lidarcldtemp(:,:,3)
	966
	967	WHERE(sumlidarcldtemp(:,:)>0.)
	968	lidarcldtemp(:,:,5)=lidarcldtemp(:,:,2)/sumlidarcldtemp(:,:)
	969	ELSEWHERE
	970	lidarcldtemp(:,:,5)=undef
	971	ENDWHERE
	972
	973	do i=1,4
	974	WHERE(lidarcldtempind(:,:).gt.0.)
	975	lidarcldtemp(:,:,i) = lidarcldtemp(:,:,i)/lidarcldtempind(:,:)
	976	ELSEWHERE
	977	lidarcldtemp(:,:,i) = undef
	978	ENDWHERE
	979	enddo
	980
	981	RETURN
[1262]	982	END SUBROUTINE COSP_CLDFRAC
	983	! ---------------------------------------------------------------
[1414]	984
[2428]	985
[1262]	986	END MODULE MOD_LMD_IPSL_STATS

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format