Context Navigation

lmd_ipsl_stats.F90 @ 3603

Last change on this file since 3603 was 2955, checked in by musat, 7 years ago
Ajout modifs de Rodrigo Guzman et Marine Bonazolla dans la version COSP/CMIP6 IM
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: 44.6 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 &
[2955]	38	,cfad2,srbval,ncat,ntype,lidarcld,lidarcldtype,lidarcldphase,cldlayer & !OPAQ
	39	,cldtype,cldlayerphase,lidarcldtmp,parasolrefl,vgrid_z,profSR) !OPAQ !TIBO
[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)
[2955]	84	integer ntype ! nb of OPAQ products (opaque and thin clouds, z_opaque) !OPAQ
[1262]	85	integer nrefl ! nb of solar zenith angles for parasol reflectances
	86
	87	real undef ! undefined value
[1414]	88	real pnorm(npoints,ncol,llm) ! lidar ATB
[1262]	89	real pmol(npoints,llm) ! molecular ATB
[1414]	90	real land(npoints) ! Landmask [0 - Ocean, 1 - Land]
[1262]	91	real pplay(npoints,llm) ! pressure on model levels (Pa)
	92	logical ok_lidar_cfad ! true if lidar CFAD diagnostics need to be computed
	93	real refl(npoints,ncol,nrefl) ! subgrid parasol reflectance ! parasol
[2428]	94	real tmp(npoints,llm) ! temp at each levels
	95	real pnorm_perp(npoints,ncol,llm) ! lidar perpendicular ATB
[2955]	96	real vgrid_z(llm) ! mid-level altitude of the output vertical grid !OPAQ
[1262]	97
	98	! c outputs :
[1414]	99	real lidarcld(npoints,llm) ! 3D "lidar" cloud fraction
[2955]	100	real lidarcldtype(npoints,llm,ntype+1) ! 3D "lidar" OPAQ type fraction + opacity !OPAQ
[2428]	101	real sub(npoints,llm) ! 3D "lidar" indice
	102	real cldlayer(npoints,ncat) ! "lidar" cloud layer fraction (low, mid, high, total)
[2955]	103	real cldtype(npoints,ntype) ! "lidar" OPAQ type covers (opaque/thin cloud + z_opaque) !OPAQ
[2428]	104
[1414]	105	real cfad2(npoints,max_bin,llm) ! CFADs of SR
	106	real srbval(max_bin) ! SR bins in CFADs
[1262]	107	real parasolrefl(npoints,nrefl)! grid-averaged parasol reflectance
[2955]	108	! real profSR(npoints,ncol,llm) ! tableau avec les subcolumns SR !TIBO
	109	real profSR(npoints,llm,ncol) ! tableau avec les subcolumns SR !TIBO2
[1262]	110
	111	! c threshold for cloud detection :
[1414]	112	real S_clr
	113	parameter (S_clr = 1.2)
[1262]	114	real S_cld
[2428]	115	parameter (S_cld = 5.) ! Thresold for cloud detection
[1262]	116	real S_att
	117	parameter (S_att = 0.01)
[2955]	118	! parameter (S_att = 0.06) !OPAQ ! Threshold for "surface detection" equivalent
[1262]	119
	120	! c local variables :
[2428]	121	integer ic,k,i,j
[1262]	122	real x3d(npoints,ncol,llm)
	123	real x3d_c(npoints,llm),pnorm_c(npoints,llm)
	124	real xmax
[2428]	125
	126	! Output variables
	127	integer,parameter :: nphase = 6 ! nb of cloud layer phase types (ice,liquid,undefined,false ice,false liquid,Percent of ice)
	128	real lidarcldphase(npoints,llm,nphase) ! 3D "lidar" phase cloud fraction
	129	real lidarcldtmp(npoints,40,5) ! 3D "lidar" phase cloud fraction as a function of temp
	130	real cldlayerphase(npoints,ncat,nphase) ! "lidar" phase low mid high cloud fraction
	131
	132	! SR detection threshold
	133	real, parameter :: S_cld_att = 30. ! New threshold for undefine cloud phase detection
	134
	135
[1262]	136	!
	137	! c -------------------------------------------------------
	138	! c 0- Initializations
	139	! c -------------------------------------------------------
[1414]	140	!
[1262]	141	! Should be modified in future version
	142	xmax=undef-1.0
	143
	144	! c -------------------------------------------------------
	145	! c 1- Lidar scattering ratio :
	146	! c -------------------------------------------------------
[2428]	147
[1262]	148	do ic = 1, ncol
	149	pnorm_c = pnorm(:,ic,:)
[1414]	150	where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
[1262]	151	x3d_c = pnorm_c/pmol
	152	elsewhere
	153	x3d_c = undef
	154	end where
[2428]	155	x3d(:,ic,:) = x3d_c
[2955]	156	! profSR(:,ic,:) = x3d(:,ic,:) !TIBO
	157	profSR(:,:,ic) = x3d(:,ic,:) !TIBO2
[1262]	158	enddo
	159
	160	! c -------------------------------------------------------
	161	! c 2- Diagnose cloud fractions (3D, low, middle, high, total)
	162	! c from subgrid-scale lidar scattering ratios :
	163	! c -------------------------------------------------------
	164
[2428]	165	CALL COSP_CLDFRAC(npoints,ncol,llm,ncat,nphase, &
	166	tmp,x3d,pnorm,pnorm_perp,pplay, S_att,S_cld,S_cld_att,undef,lidarcld, &
	167	cldlayer,lidarcldphase,sub,cldlayerphase,lidarcldtmp)
[1262]	168
[2955]	169	CALL COSP_OPAQ(npoints,ncol,llm,ntype,x3d,S_cld,undef,lidarcldtype, & !OPAQ
	170	cldtype,vgrid_z) !OPAQ
	171
[1262]	172	! c -------------------------------------------------------
[1414]	173	! c 3- CFADs
[1262]	174	! c -------------------------------------------------------
	175	if (ok_lidar_cfad) then
	176	!
	177	! c CFADs of subgrid-scale lidar scattering ratios :
	178	! c -------------------------------------------------------
[1414]	179	CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin,undef, &
[1262]	180	x3d, &
	181	S_att,S_clr,xmax,cfad2,srbval)
	182
	183	endif ! ok_lidar_cfad
	184	! c -------------------------------------------------------
	185
	186	! c -------------------------------------------------------
	187	! c 4- Compute grid-box averaged Parasol reflectances
	188	! c -------------------------------------------------------
	189
	190	parasolrefl(:,:) = 0.0
	191
	192	do k = 1, nrefl
	193	do ic = 1, ncol
	194	parasolrefl(:,k) = parasolrefl(:,k) + refl(:,ic,k)
	195	enddo
	196	enddo
	197
	198	do k = 1, nrefl
	199	parasolrefl(:,k) = parasolrefl(:,k) / float(ncol)
	200	! if land=1 -> parasolrefl=undef
	201	! if land=0 -> parasolrefl=parasolrefl
	202	parasolrefl(:,k) = parasolrefl(:,k) * MAX(1.0-land(:),0.0) &
[1414]	203	+ (1.0 - MAX(1.0-land(:),0.0))*undef
[1262]	204	enddo
	205
	206	RETURN
	207	END SUBROUTINE diag_lidar
[1414]	208
	209
[1262]	210	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	211	!-------------------- FUNCTION COSP_CFAD_SR ------------------------
	212	! Author: Sandrine Bony (LMD/IPSL, CNRS, Paris)
	213	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[1414]	214	SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins,undef, &
[1262]	215	x,S_att,S_clr,xmax,cfad,srbval)
	216	IMPLICIT NONE
	217
	218	!--- Input arguments
	219	! Npoints: Number of horizontal points
	220	! Ncolumns: Number of subcolumns
	221	! Nlevels: Number of levels
	222	! Nbins: Number of x axis bins
	223	! xmax: maximum value allowed for x
	224	! S_att: Threshold for full attenuation
	225	! S_clr: Threshold for clear-sky layer
	226	!
	227	!--- Input-Outout arguments
	228	! x: variable to process (Npoints,Ncolumns,Nlevels), mofified where saturation occurs
	229	!
	230	! -- Output arguments
	231	! srbval : values of the histogram bins
	232	! cfad: 2D histogram on each horizontal point
	233
	234	! Input arguments
	235	integer Npoints,Ncolumns,Nlevels,Nbins
[1414]	236	real xmax,S_att,S_clr,undef
[1262]	237	! Input-output arguments
	238	real x(Npoints,Ncolumns,Nlevels)
	239	! Output :
	240	real cfad(Npoints,Nbins,Nlevels)
	241	real srbval(Nbins)
	242	! Local variables
	243	integer i, j, k, ib
[1414]	244	real srbval_ext(0:Nbins)
[1262]	245
	246	! c -------------------------------------------------------
	247	! c 0- Initializations
	248	! c -------------------------------------------------------
	249	if ( Nbins .lt. 6) return
	250
	251	srbval(1) = S_att
	252	srbval(2) = S_clr
	253	srbval(3) = 3.0
	254	srbval(4) = 5.0
	255	srbval(5) = 7.0
	256	srbval(6) = 10.0
	257	do i = 7, MIN(10,Nbins)
	258	srbval(i) = srbval(i-1) + 5.0
	259	enddo
	260	DO i = 11, MIN(13,Nbins)
	261	srbval(i) = srbval(i-1) + 10.0
	262	enddo
	263	srbval(MIN(14,Nbins)) = 80.0
	264	srbval(Nbins) = xmax
	265	cfad(:,:,:) = 0.0
	266
[1414]	267	srbval_ext(1:Nbins) = srbval
	268	srbval_ext(0) = -1.0
[1262]	269	! c -------------------------------------------------------
	270	! c c- Compute CFAD
	271	! c -------------------------------------------------------
[1414]	272	do j = 1, Nlevels
	273	do ib = 1, Nbins
	274	do k = 1, Ncolumns
	275	do i = 1, Npoints
	276	if (x(i,k,j) /= undef) then
	277	if ((x(i,k,j).gt.srbval_ext(ib-1)).and.(x(i,k,j).le.srbval_ext(ib))) &
	278	cfad(i,ib,j) = cfad(i,ib,j) + 1.0
	279	else
	280	cfad(i,ib,j) = undef
	281	endif
	282	enddo
	283	enddo
	284	enddo
	285	enddo
[1262]	286
[1414]	287	where (cfad .ne. undef) cfad = cfad / float(Ncolumns)
[1262]	288
	289	! c -------------------------------------------------------
	290	RETURN
	291	END SUBROUTINE COSP_CFAD_SR
	292
[2428]	293
[1262]	294	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	295	!-------------------- SUBROUTINE COSP_CLDFRAC -------------------
[1414]	296	! c Purpose: Cloud fraction diagnosed from lidar measurements
[1262]	297	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[2428]	298	SUBROUTINE COSP_CLDFRAC(Npoints,Ncolumns,Nlevels,Ncat,Nphase, &
	299	tmp,x,ATB,ATBperp,pplay,S_att,S_cld,S_cld_att,undef,lidarcld, &
	300	cldlayer,lidarcldphase,nsub,cldlayerphase,lidarcldtemp)
	301
	302
[1262]	303	IMPLICIT NONE
	304	! Input arguments
	305	integer Npoints,Ncolumns,Nlevels,Ncat
	306	real x(Npoints,Ncolumns,Nlevels)
[2428]	307
	308
	309	! Local parameters
	310	integer nphase ! nb of cloud layer phase types
	311	! (ice,liquid,undefined,false ice,false liquid,Percent of ice)
	312	integer,parameter :: Ntemp=40 ! indice of the temperature vector
	313	integer ip, k, iz, ic, ncol, nlev, i, itemp ! loop indice
	314	real S_cld_att ! New threshold for undefine cloud phase detection (SR=30)
	315	integer toplvlsat ! level of the first cloud with SR>30
	316	real alpha50, beta50, gamma50, delta50, epsilon50, zeta50 ! Polynomial Coef of the phase
	317	! discrimination line
	318
	319	! Input variables
	320	real tmp(Npoints,Nlevels) ! temperature
	321	real ATB(Npoints,Ncolumns,Nlevels) ! 3D Attenuated backscatter
	322	real ATBperp(Npoints,Ncolumns,Nlevels) ! 3D perpendicular attenuated backscatter
[1262]	323	real pplay(Npoints,Nlevels)
	324	real S_att,S_cld
	325	real undef
[2428]	326
	327	! Output variables
	328	real lidarcldtemp(Npoints,Ntemp,5) ! 3D Temperature 1=tot,2=ice,3=liq,4=undef,5=ice/ice+liq
	329	real tempmod(Ntemp+1) ! temperature bins
	330	real lidarcldphase(Npoints,Nlevels,Nphase) ! 3D cloud phase fraction
	331	real cldlayerphase(Npoints,Ncat,Nphase) ! low, middle, high, total cloud fractions for ice liquid and undefine phase
[1262]	332	real lidarcld(Npoints,Nlevels) ! 3D cloud fraction
	333	real cldlayer(Npoints,Ncat) ! low, middle, high, total cloud fractions
[2428]	334
[1262]	335	! Local variables
[2428]	336	real tmpi(Npoints,Ncolumns,Nlevels) ! temperature of ice cld
	337	real tmpl(Npoints,Ncolumns,Nlevels) ! temperature of liquid cld
	338	real tmpu(Npoints,Ncolumns,Nlevels) ! temperature of undef cld
	339
	340	real checktemp, ATBperp_tmp ! temporary variable
	341	real checkcldlayerphase, checkcldlayerphase2 ! temporary variable
	342	real sumlidarcldtemp(Npoints,Ntemp) ! temporary variable
	343
	344	real cldlayphase(Npoints,Ncolumns,Ncat,Nphase) ! subgrided low mid high phase cloud fraction
	345	real cldlayerphasetmp(Npoints,Ncat) ! temporary variable
	346	real cldlayerphasesum(Npoints,Ncat) ! temporary variable
	347	real lidarcldtempind(Npoints,Ntemp) ! 3D Temperature indice
	348	real lidarcldphasetmp(Npoints,Nlevels) ! 3D sum of ice and liquid cloud occurences
	349
	350
	351	! Local variables
[1262]	352	real p1
	353	real cldy(Npoints,Ncolumns,Nlevels)
	354	real srok(Npoints,Ncolumns,Nlevels)
	355	real cldlay(Npoints,Ncolumns,Ncat)
	356	real nsublay(Npoints,Ncolumns,Ncat), nsublayer(Npoints,Ncat)
	357	real nsub(Npoints,Nlevels)
	358
[2428]	359	#ifdef SYS_SX
[1414]	360	real cldlay1(Npoints,Ncolumns)
	361	real cldlay2(Npoints,Ncolumns)
	362	real cldlay3(Npoints,Ncolumns)
	363	real nsublay1(Npoints,Ncolumns)
	364	real nsublay2(Npoints,Ncolumns)
	365	real nsublay3(Npoints,Ncolumns)
[2428]	366	#endif
[1262]	367
[1414]	368
[2428]	369
	370
[1262]	371	! ---------------------------------------------------------------
[1414]	372	! 1- initialization
[1262]	373	! ---------------------------------------------------------------
	374
	375	if ( Ncat .ne. 4 ) then
	376	print *,'Error in lmd_ipsl_stats.cosp_cldfrac, Ncat must be 4, not',Ncat
	377	stop
	378	endif
	379
	380	lidarcld = 0.0
	381	nsub = 0.0
	382	cldlay = 0.0
	383	nsublay = 0.0
	384
[2428]	385	ATBperp_tmp = 0.
	386	lidarcldphase(:,:,:) = 0.
	387	cldlayphase(:,:,:,:) = 0.
	388	cldlayerphase(:,:,:) = 0.
	389	tmpi(:,:,:) = 0.
	390	tmpl(:,:,:) = 0.
	391	tmpu(:,:,:) = 0.
	392	cldlayerphasesum(:,:) = 0.
	393	lidarcldtemp(:,:,:) = 0.
	394	lidarcldtempind(:,:) = 0.
	395	sumlidarcldtemp(:,:) = 0.
	396	toplvlsat=0
	397	lidarcldphasetmp(:,:) = 0.
	398
	399	! temperature bins
	400	tempmod=(/-273.15,-90.,-87.,-84.,-81.,-78.,-75.,-72.,-69.,-66.,-63.,-60.,-57., &
	401	-54.,-51.,-48.,-45.,-42.,-39.,-36.,-33.,-30.,-27.,-24.,-21.,-18., &
	402	-15.,-12.,-9.,-6.,-3.,0.,3.,6.,9.,12.,15.,18.,21.,24.,200. /)
	403
	404	! convert C to K
	405	tempmod=tempmod+273.15
	406
	407	! Polynomial coefficient of the phase discrimination line used to separate liquid from ice
	408	! (Cesana and Chepfer, JGR, 2013)
	409	! ATBperp = ATB^5alpha50 + ATB^4beta50 + ATB^3gamma50 + ATB^2delta50 + ATB*epsilon50 + zeta50
	410	alpha50 = 9.0322e+15
	411	beta50 = -2.1358e+12
	412	gamma50 = 173.3963e06
	413	delta50 = -3.9514e03
	414	epsilon50 = 0.2559
	415	zeta50 = -9.4776e-07
	416
	417
[1262]	418	! ---------------------------------------------------------------
	419	! 2- Cloud detection
	420	! ---------------------------------------------------------------
	421
	422	do k = 1, Nlevels
	423
	424	! cloud detection at subgrid-scale:
	425	where ( (x(:,:,k).gt.S_cld) .and. (x(:,:,k).ne. undef) )
	426	cldy(:,:,k)=1.0
	427	elsewhere
	428	cldy(:,:,k)=0.0
	429	endwhere
	430
	431	! number of usefull sub-columns:
[1414]	432	where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef) )
[1262]	433	srok(:,:,k)=1.0
	434	elsewhere
	435	srok(:,:,k)=0.0
	436	endwhere
	437
	438	enddo ! k
	439
[2428]	440
[1262]	441	! ---------------------------------------------------------------
	442	! 3- grid-box 3D cloud fraction and layered cloud fractions (ISCCP pressure
	443	! categories) :
	444	! ---------------------------------------------------------------
[1414]	445	lidarcld = 0.0
	446	nsub = 0.0
[2428]	447	#ifdef SYS_SX
[1414]	448	!! XXX: Use cldlay[1-3] and nsublay[1-3] to avoid bank-conflicts.
	449	cldlay1 = 0.0
	450	cldlay2 = 0.0
	451	cldlay3 = 0.0
	452	cldlay(:,:,4) = 0.0 !! XXX: Ncat == 4
	453	nsublay1 = 0.0
	454	nsublay2 = 0.0
	455	nsublay3 = 0.0
	456	nsublay(:,:,4) = 0.0
[2428]	457
[1262]	458	do k = Nlevels, 1, -1
	459	do ic = 1, Ncolumns
	460	do ip = 1, Npoints
[2428]	461
	462	if(srok(ip,ic,k).gt.0.)then
	463	! Computation of the cloud fraction as a function of the temperature
	464	! instead of height, for ice,liquid and all clouds
	465	do itemp=1,Ntemp
	466	if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
	467	lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1.
	468	endif
	469	enddo
	470	endif
	471
	472	if (cldy(ip,ic,k).eq.1.) then
	473	do itemp=1,Ntemp
	474	if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
	475	lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1.
	476	endif
	477	enddo
	478	endif
	479
[1414]	480	p1 = pplay(ip,k)
	481
	482	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
	483	cldlay3(ip,ic) = MAX(cldlay3(ip,ic), cldy(ip,ic,k))
	484	nsublay3(ip,ic) = MAX(nsublay3(ip,ic), srok(ip,ic,k))
	485	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid clouds
	486	cldlay2(ip,ic) = MAX(cldlay2(ip,ic), cldy(ip,ic,k))
	487	nsublay2(ip,ic) = MAX(nsublay2(ip,ic), srok(ip,ic,k))
	488	else
	489	cldlay1(ip,ic) = MAX(cldlay1(ip,ic), cldy(ip,ic,k))
	490	nsublay1(ip,ic) = MAX(nsublay1(ip,ic), srok(ip,ic,k))
[1262]	491	endif
	492
[1414]	493	cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4), cldy(ip,ic,k))
[1262]	494	lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k)
	495	nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
	496	nsub(ip,k)=nsub(ip,k) + srok(ip,ic,k)
	497	enddo
	498	enddo
	499	enddo
[1414]	500	cldlay(:,:,1) = cldlay1
	501	cldlay(:,:,2) = cldlay2
	502	cldlay(:,:,3) = cldlay3
	503	nsublay(:,:,1) = nsublay1
	504	nsublay(:,:,2) = nsublay2
	505	nsublay(:,:,3) = nsublay3
[2428]	506	#else
	507	cldlay = 0.0
	508	nsublay = 0.0
	509	do k = Nlevels, 1, -1
	510	do ic = 1, Ncolumns
	511	do ip = 1, Npoints
[1262]	512
[2428]	513	! Computation of the cloud fraction as a function of the temperature
	514	! instead of height, for ice,liquid and all clouds
	515	if(srok(ip,ic,k).gt.0.)then
	516	do itemp=1,Ntemp
	517	if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
	518	lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1.
	519	endif
	520	enddo
	521	endif
	522
	523	if(cldy(ip,ic,k).eq.1.)then
	524	do itemp=1,Ntemp
	525	if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
	526	lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1.
	527	endif
	528	enddo
	529	endif
	530	!
	531
	532	iz=1
	533	p1 = pplay(ip,k)
	534	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
	535	iz=3
	536	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid clouds
	537	iz=2
	538	endif
	539
	540	cldlay(ip,ic,iz) = MAX(cldlay(ip,ic,iz),cldy(ip,ic,k))
	541	cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4),cldy(ip,ic,k))
	542	lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k)
	543
	544	nsublay(ip,ic,iz) = MAX(nsublay(ip,ic,iz),srok(ip,ic,k))
	545	nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
	546	nsub(ip,k)=nsub(ip,k) + srok(ip,ic,k)
	547
	548	enddo
	549	enddo
	550	enddo
	551	#endif
	552
	553
[1262]	554	! -- grid-box 3D cloud fraction
	555
	556	where ( nsub(:,:).gt.0.0 )
	557	lidarcld(:,:) = lidarcld(:,:)/nsub(:,:)
	558	elsewhere
	559	lidarcld(:,:) = undef
	560	endwhere
	561
	562	! -- layered cloud fractions
	563
	564	cldlayer = 0.0
	565	nsublayer = 0.0
	566
	567	do iz = 1, Ncat
	568	do ic = 1, Ncolumns
	569
[1414]	570	cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz)
	571	nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz)
[1262]	572
	573	enddo
	574	enddo
	575	where ( nsublayer(:,:).gt.0.0 )
	576	cldlayer(:,:) = cldlayer(:,:)/nsublayer(:,:)
	577	elsewhere
	578	cldlayer(:,:) = undef
	579	endwhere
	580
[2428]	581	! ---------------------------------------------------------------
	582	! 4- grid-box 3D cloud Phase :
	583	! ---------------------------------------------------------------
	584	! ---------------------------------------------------------------
	585	! 4.1 - For Cloudy pixels with 8.16km < z < 19.2km
	586	! ---------------------------------------------------------------
	587	do ncol=1,Ncolumns
	588	do i=1,Npoints
	589
	590	do nlev=Nlevels,18,-1 ! from 19.2km until 8.16km
	591	p1 = pplay(i,nlev)
	592
	593
	594	! Avoid zero values
	595	if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
	596	! Computation of the ATBperp along the phase discrimination line
	597	ATBperp_tmp = (ATB(i,ncol,nlev)*5)alpha50 + (ATB(i,ncol,nlev)*4)beta50 + &
	598	(ATB(i,ncol,nlev)*3)gamma50 + (ATB(i,ncol,nlev)*2)delta50 + &
	599	ATB(i,ncol,nlev)*epsilon50 + zeta50
	600
	601	!____________________________________________________________________________________________________
	602	!
	603	!4.1.a Ice: ATBperp above the phase discrimination line
	604	!____________________________________________________________________________________________________
	605	!
	606	if( (ATBperp(i,ncol,nlev)-ATBperp_tmp).ge.0. )then ! Ice clouds
	607	! ICE with temperature above 273,15°K = Liquid (false ice)
	608	if(tmp(i,nlev).gt.273.15)then ! Temperature above 273,15 K
	609	! Liquid: False ice corrected by the temperature to Liquid
	610	lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1. ! false ice detection ==> added to Liquid
	611	tmpl(i,ncol,nlev)=tmp(i,nlev)
	612	lidarcldphase(i,nlev,5)=lidarcldphase(i,nlev,5)+1. ! keep the information "temperature criterium used"
	613	! to classify the phase cloud
	614	cldlayphase(i,ncol,4,2) = 1. ! tot cloud
	615	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	616	cldlayphase(i,ncol,3,2) = 1.
	617	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	618	cldlayphase(i,ncol,2,2) = 1.
	619	else ! low cloud
	620	cldlayphase(i,ncol,1,2) = 1.
	621	endif
	622	cldlayphase(i,ncol,4,5) = 1. ! tot cloud
	623	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	624	cldlayphase(i,ncol,3,5) = 1.
	625	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	626	cldlayphase(i,ncol,2,5) = 1.
	627	else ! low cloud
	628	cldlayphase(i,ncol,1,5) = 1.
	629	endif
	630
	631	else
	632	! ICE with temperature below 273,15°K
	633	lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1.
	634	tmpi(i,ncol,nlev)=tmp(i,nlev)
	635	cldlayphase(i,ncol,4,1) = 1. ! tot cloud
	636	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	637	cldlayphase(i,ncol,3,1) = 1.
	638	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	639	cldlayphase(i,ncol,2,1) = 1.
	640	else ! low cloud
	641	cldlayphase(i,ncol,1,1) = 1.
	642	endif
	643
	644	endif
	645
	646	!____________________________________________________________________________________________________
	647	!
	648	! 4.1.b Liquid: ATBperp below the phase discrimination line
	649	!____________________________________________________________________________________________________
	650	!
	651	else ! Liquid clouds
	652	! Liquid with temperature above 231,15°K
	653	if(tmp(i,nlev).gt.231.15)then
	654	lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.
	655	tmpl(i,ncol,nlev)=tmp(i,nlev)
	656	cldlayphase(i,ncol,4,2) = 1. ! tot cloud
	657	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	658	cldlayphase(i,ncol,3,2) = 1.
	659	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	660	cldlayphase(i,ncol,2,2) = 1.
	661	else ! low cloud
	662	cldlayphase(i,ncol,1,2) = 1.
	663	endif
	664
	665	else
	666	! Liquid with temperature below 231,15°K = Ice (false liquid)
	667	tmpi(i,ncol,nlev)=tmp(i,nlev)
	668	lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1. ! false liquid detection ==> added to ice
	669	lidarcldphase(i,nlev,4)=lidarcldphase(i,nlev,4)+1. ! keep the information "temperature criterium used"
	670	! to classify the phase cloud
	671	cldlayphase(i,ncol,4,4) = 1. ! tot cloud
	672	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	673	cldlayphase(i,ncol,3,4) = 1.
	674	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	675	cldlayphase(i,ncol,2,4) = 1.
	676	else ! low cloud
	677	cldlayphase(i,ncol,1,4) = 1.
	678	endif
	679	cldlayphase(i,ncol,4,1) = 1. ! tot cloud
	680	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	681	cldlayphase(i,ncol,3,1) = 1.
	682	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	683	cldlayphase(i,ncol,2,1) = 1.
	684	else ! low cloud
	685	cldlayphase(i,ncol,1,1) = 1.
	686	endif
	687
	688	endif
	689
	690	endif ! end of discrimination condition
	691	endif ! end of cloud condition
	692	enddo ! end of altitude loop
	693
	694
	695
	696	! ---------------------------------------------------------------
	697	! 4.2 - For Cloudy pixels with 0km < z < 8.16km
	698	! ---------------------------------------------------------------
	699
	700	toplvlsat=0
	701	do nlev=17,1,-1 ! from 8.16km until 0km
	702	p1 = pplay(i,nlev)
	703
	704	if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
	705	! Phase discrimination line : ATBperp = ATB^5alpha50 + ATB^4beta50 + ATB^3gamma50 + ATB^2delta50
	706	! + ATB*epsilon50 + zeta50
	707	! Computation of the ATBperp of the phase discrimination line
	708	ATBperp_tmp = (ATB(i,ncol,nlev)*5)alpha50 + (ATB(i,ncol,nlev)*4)beta50 + &
	709	(ATB(i,ncol,nlev)*3)gamma50 + (ATB(i,ncol,nlev)*2)delta50 + &
	710	ATB(i,ncol,nlev)*epsilon50 + zeta50
	711	!____________________________________________________________________________________________________
	712	!
	713	! 4.2.a Ice: ATBperp above the phase discrimination line
	714	!____________________________________________________________________________________________________
	715	!
	716	! ICE with temperature above 273,15°K = Liquid (false ice)
	717	if( (ATBperp(i,ncol,nlev)-ATBperp_tmp).ge.0. )then ! Ice clouds
	718	if(tmp(i,nlev).gt.273.15)then
	719	lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1. ! false ice ==> liq
	720	tmpl(i,ncol,nlev)=tmp(i,nlev)
	721	lidarcldphase(i,nlev,5)=lidarcldphase(i,nlev,5)+1.
	722
	723	cldlayphase(i,ncol,4,2) = 1. ! tot cloud
	724	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	725	cldlayphase(i,ncol,3,2) = 1.
	726	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	727	cldlayphase(i,ncol,2,2) = 1.
	728	else ! low cloud
	729	cldlayphase(i,ncol,1,2) = 1.
	730	endif
	731
	732	cldlayphase(i,ncol,4,5) = 1. ! tot cloud
	733	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	734	cldlayphase(i,ncol,3,5) = 1.
	735	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	736	cldlayphase(i,ncol,2,5) = 1.
	737	else ! low cloud
	738	cldlayphase(i,ncol,1,5) = 1.
	739	endif
	740
	741	else
	742	! ICE with temperature below 273,15°K
	743	lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1.
	744	tmpi(i,ncol,nlev)=tmp(i,nlev)
	745
	746	cldlayphase(i,ncol,4,1) = 1. ! tot cloud
	747	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	748	cldlayphase(i,ncol,3,1) = 1.
	749	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	750	cldlayphase(i,ncol,2,1) = 1.
	751	else ! low cloud
	752	cldlayphase(i,ncol,1,1) = 1.
	753	endif
	754
	755	endif
	756
	757	!____________________________________________________________________________________________________
	758	!
	759	! 4.2.b Liquid: ATBperp below the phase discrimination line
	760	!____________________________________________________________________________________________________
	761	!
	762	else
	763	! Liquid with temperature above 231,15°K
	764	if(tmp(i,nlev).gt.231.15)then
	765	lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.
	766	tmpl(i,ncol,nlev)=tmp(i,nlev)
	767
	768	cldlayphase(i,ncol,4,2) = 1. ! tot cloud
	769	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	770	cldlayphase(i,ncol,3,2) = 1.
	771	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	772	cldlayphase(i,ncol,2,2) = 1.
	773	else ! low cloud
	774	cldlayphase(i,ncol,1,2) = 1.
	775	endif
	776
	777	else
	778	! Liquid with temperature below 231,15°K = Ice (false liquid)
	779	tmpi(i,ncol,nlev)=tmp(i,nlev)
	780	lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1. ! false liq ==> ice
	781	lidarcldphase(i,nlev,4)=lidarcldphase(i,nlev,4)+1. ! false liq ==> ice
	782
	783	cldlayphase(i,ncol,4,4) = 1. ! tot cloud
	784	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	785	cldlayphase(i,ncol,3,4) = 1.
	786	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	787	cldlayphase(i,ncol,2,4) = 1.
	788	else ! low cloud
	789	cldlayphase(i,ncol,1,4) = 1.
	790	endif
	791
	792	cldlayphase(i,ncol,4,1) = 1. ! tot cloud
	793	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	794	cldlayphase(i,ncol,3,1) = 1.
	795	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	796	cldlayphase(i,ncol,2,1) = 1.
	797	else ! low cloud
	798	cldlayphase(i,ncol,1,1) = 1.
	799	endif
	800
	801	endif
	802	endif ! end of discrimination condition
	803
	804	toplvlsat=0
	805
	806	! Find the level of the highest cloud with SR>30
	807	if(x(i,ncol,nlev).gt.S_cld_att)then ! SR > 30.
	808	toplvlsat=nlev-1
	809	goto 99
	810	endif
	811
	812	endif ! end of cloud condition
	813	enddo ! end of altitude loop
	814
	815	99 continue
	816
	817	!____________________________________________________________________________________________________
	818	!
	819	! Undefined phase: For a cloud located below another cloud with SR>30
	820	! see Cesana and Chepfer 2013 Sect.III.2
	821	!____________________________________________________________________________________________________
	822	!
	823	if(toplvlsat.ne.0)then
	824	do nlev=toplvlsat,1,-1
	825	p1 = pplay(i,nlev)
	826	if(cldy(i,ncol,nlev).eq.1.)then
	827	lidarcldphase(i,nlev,3)=lidarcldphase(i,nlev,3)+1.
	828	tmpu(i,ncol,nlev)=tmp(i,nlev)
	829
	830	cldlayphase(i,ncol,4,3) = 1. ! tot cloud
	831	if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high cloud
	832	cldlayphase(i,ncol,3,3) = 1.
	833	else if(p1.ge.(440.100.) .and. p1.lt.(680.100.)) then ! mid cloud
	834	cldlayphase(i,ncol,2,3) = 1.
	835	else ! low cloud
	836	cldlayphase(i,ncol,1,3) = 1.
	837	endif
	838
	839	endif
	840	enddo
	841	endif
	842
	843	toplvlsat=0
	844
	845	enddo
	846	enddo
	847
	848
	849
	850	!____________________________________________________________________________________________________
	851	!
	852	! Computation of final cloud phase diagnosis
	853	!____________________________________________________________________________________________________
	854	!
	855
	856	! Compute the Ice percentage in cloud = ice/(ice+liq) as a function
	857	! of the occurrences
	858	lidarcldphasetmp(:,:)=lidarcldphase(:,:,1)+lidarcldphase(:,:,2);
	859	WHERE (lidarcldphasetmp(:,:).gt. 0.)
	860	lidarcldphase(:,:,6)=lidarcldphase(:,:,1)/lidarcldphasetmp(:,:)
	861	ELSEWHERE
	862	lidarcldphase(:,:,6) = undef
	863	ENDWHERE
	864
	865	! Compute Phase 3D Cloud Fraction
	866	WHERE ( nsub(:,:).gt.0.0 )
	867	lidarcldphase(:,:,1)=lidarcldphase(:,:,1)/nsub(:,:)
	868	lidarcldphase(:,:,2)=lidarcldphase(:,:,2)/nsub(:,:)
	869	lidarcldphase(:,:,3)=lidarcldphase(:,:,3)/nsub(:,:)
	870	lidarcldphase(:,:,4)=lidarcldphase(:,:,4)/nsub(:,:)
	871	lidarcldphase(:,:,5)=lidarcldphase(:,:,5)/nsub(:,:)
	872	ELSEWHERE
	873	lidarcldphase(:,:,1) = undef
	874	lidarcldphase(:,:,2) = undef
	875	lidarcldphase(:,:,3) = undef
	876	lidarcldphase(:,:,4) = undef
	877	lidarcldphase(:,:,5) = undef
	878	ENDWHERE
	879
	880
	881	! Compute Phase low mid high cloud fractions
	882	do iz = 1, Ncat
	883	do i=1,Nphase-3
	884	do ic = 1, Ncolumns
	885	cldlayerphase(:,iz,i)=cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i)
	886	cldlayerphasesum(:,iz)=cldlayerphasesum(:,iz)+cldlayphase(:,ic,iz,i)
	887	enddo
	888	enddo
	889	enddo
	890
	891	do iz = 1, Ncat
	892	do i=4,5
	893	do ic = 1, Ncolumns
	894	cldlayerphase(:,iz,i)=cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i)
	895	enddo
	896	enddo
	897	enddo
	898
	899	! Compute the Ice percentage in cloud = ice/(ice+liq)
	900	cldlayerphasetmp(:,:)=cldlayerphase(:,:,1)+cldlayerphase(:,:,2)
	901	WHERE (cldlayerphasetmp(:,:).gt. 0.)
	902	cldlayerphase(:,:,6)=cldlayerphase(:,:,1)/cldlayerphasetmp(:,:)
	903	ELSEWHERE
	904	cldlayerphase(:,:,6) = undef
	905	ENDWHERE
	906
	907	do i=1,Nphase-1
	908	WHERE ( cldlayerphasesum(:,:).gt.0.0 )
	909	cldlayerphase(:,:,i) = (cldlayerphase(:,:,i)/cldlayerphasesum(:,:)) * cldlayer(:,:)
	910	ENDWHERE
	911	enddo
	912
	913
	914	do i=1,Npoints
	915	do iz=1,Ncat
	916	checkcldlayerphase=0.
	917	checkcldlayerphase2=0.
	918
	919	if (cldlayerphasesum(i,iz).gt.0.0 )then
	920	do ic=1,Nphase-3
	921	checkcldlayerphase=checkcldlayerphase+cldlayerphase(i,iz,ic)
	922	enddo
	923	checkcldlayerphase2=cldlayer(i,iz)-checkcldlayerphase
	924	if( (checkcldlayerphase2.gt.0.01).or.(checkcldlayerphase2.lt.-0.01) ) print *, checkcldlayerphase,cldlayer(i,iz)
	925
	926	endif
	927
	928	enddo
	929	enddo
	930
	931	do i=1,Nphase-1
	932	WHERE ( nsublayer(:,:).eq.0.0 )
	933	cldlayerphase(:,:,i) = undef
	934	ENDWHERE
	935	enddo
	936
	937
	938
	939	! Compute Phase 3D as a function of temperature
	940	do nlev=1,Nlevels
	941	do ncol=1,Ncolumns
	942	do i=1,Npoints
	943	do itemp=1,Ntemp
	944	if(tmpi(i,ncol,nlev).gt.0.)then
	945	if( (tmpi(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpi(i,ncol,nlev).lt.tempmod(itemp+1)) )then
	946	lidarcldtemp(i,itemp,2)=lidarcldtemp(i,itemp,2)+1.
	947	endif
	948	elseif(tmpl(i,ncol,nlev).gt.0.)then
	949	if( (tmpl(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpl(i,ncol,nlev).lt.tempmod(itemp+1)) )then
	950	lidarcldtemp(i,itemp,3)=lidarcldtemp(i,itemp,3)+1.
	951	endif
	952	elseif(tmpu(i,ncol,nlev).gt.0.)then
	953	if( (tmpu(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpu(i,ncol,nlev).lt.tempmod(itemp+1)) )then
	954	lidarcldtemp(i,itemp,4)=lidarcldtemp(i,itemp,4)+1.
	955	endif
	956	endif
	957	enddo
	958	enddo
	959	enddo
	960	enddo
	961
	962	! Check temperature cloud fraction
	963	do i=1,Npoints
	964	do itemp=1,Ntemp
	965	checktemp=lidarcldtemp(i,itemp,2)+lidarcldtemp(i,itemp,3)+lidarcldtemp(i,itemp,4)
	966
	967	if(checktemp.NE.lidarcldtemp(i,itemp,1))then
	968	print *, i,itemp
	969	print *, lidarcldtemp(i,itemp,1:4)
	970	endif
	971
	972	enddo
	973	enddo
	974
	975	! Compute the Ice percentage in cloud = ice/(ice+liq)
	976	! sumlidarcldtemp=sum(lidarcldtemp(:,:,2:3),3)
	977	sumlidarcldtemp(:,:)=lidarcldtemp(:,:,2)+lidarcldtemp(:,:,3)
	978
	979	WHERE(sumlidarcldtemp(:,:)>0.)
	980	lidarcldtemp(:,:,5)=lidarcldtemp(:,:,2)/sumlidarcldtemp(:,:)
	981	ELSEWHERE
	982	lidarcldtemp(:,:,5)=undef
	983	ENDWHERE
	984
	985	do i=1,4
	986	WHERE(lidarcldtempind(:,:).gt.0.)
	987	lidarcldtemp(:,:,i) = lidarcldtemp(:,:,i)/lidarcldtempind(:,:)
	988	ELSEWHERE
	989	lidarcldtemp(:,:,i) = undef
	990	ENDWHERE
	991	enddo
	992
	993	RETURN
[1262]	994	END SUBROUTINE COSP_CLDFRAC
	995	! ---------------------------------------------------------------
[1414]	996
[2955]	997	! BEGINNING OF OPAQ CHANGES
	998	! ####################################################################################
	999	! SUBROUTINE cosp_opaq
	1000	! Conventions: Ntype must be equal to 3 (opaque cloud, thin cloud, z_opaque)
	1001	! ####################################################################################
	1002	SUBROUTINE COSP_OPAQ(Npoints,Ncolumns,Nlevels,Ntype,x,S_cld,undef,lidarcldtype, &
	1003	cldtype,vgrid_z)
[2428]	1004
[2955]	1005	IMPLICIT NONE
	1006	! Input arguments
	1007	integer Npoints,Ncolumns,Nlevels,Ntype
	1008	real x(Npoints,Ncolumns,Nlevels)
	1009	real S_cld
	1010	real undef
	1011	real vgrid_z(Nlevels)
	1012	! Output :
	1013	real lidarcldtype(Npoints,Nlevels,Ntype+1) ! 3D "lidar" OPAQ type + opacity fraction
	1014	real cldtype(Npoints,Ntype) ! opaque and thin cloud covers, z_opaque
	1015	! Local variables
	1016	integer ip, k, iz, ic, zopac
	1017	real p1
	1018	real cldy(Npoints,Ncolumns,Nlevels)
	1019	real cldyopaq(Npoints,Ncolumns,Nlevels)
	1020	real srok(Npoints,Ncolumns,Nlevels)
	1021	real srokopaq(Npoints,Ncolumns,Nlevels)
	1022	real cldlay(Npoints,Ncolumns,Ntype+1) ! opaque, thin, z_opaque and all cloud cover
	1023	real nsublay(Npoints,Ncolumns,Ntype+1) ! opaque, thin, z_opaque and all cloud cover
	1024	real nsublayer(Npoints,Ntype)
	1025	real nsub(Npoints,Nlevels)
	1026	real nsubopaq(Npoints,Nlevels)
	1027	real S_att_opaq
	1028	real S_att
	1029
	1030	! ####################################################################################
	1031	! 1) Initialize
	1032	! ####################################################################################
	1033	cldtype = 0.0
	1034	lidarcldtype = 0.0
	1035	nsub = 0.0
	1036	nsubopaq = 0.0
	1037	cldlay = 0.0
	1038	nsublay = 0.0
	1039	nsublayer = 0.0
	1040	S_att_opaq = 0.06 ! Fully Attenuated threshold, from Guzman et al. 2017, JGR-A
	1041	S_att = 0.01
	1042
	1043	! ####################################################################################
	1044	! 2) Cloud detection and Fully attenuated layer detection
	1045	! ####################################################################################
	1046	do k=1,Nlevels
	1047	! Cloud detection at subgrid-scale:
	1048	where ( (x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
	1049	cldy(:,:,k)=1.0
	1050	elsewhere
	1051	cldy(:,:,k)=0.0
	1052	endwhere
	1053	! Fully attenuated layer detection at subgrid-scale:
	1054	where ( (x(:,:,k) .gt. 0.0) .and. (x(:,:,k) .lt. S_att_opaq) .and. (x(:,:,k) .ne. undef) )
	1055	cldyopaq(:,:,k)=1.0
	1056	elsewhere
	1057	cldyopaq(:,:,k)=0.0
	1058	endwhere
	1059
	1060	! Number of useful sub-column layers:
	1061	where ( (x(:,:,k) .gt. S_att) .and. (x(:,:,k) .ne. undef) )
	1062	srok(:,:,k)=1.0
	1063	elsewhere
	1064	srok(:,:,k)=0.0
	1065	endwhere
	1066	! Number of useful sub-columns layers for z_opaque 3D fraction:
	1067	where ( (x(:,:,k) .gt. 0.0) .and. (x(:,:,k) .ne. undef) )
	1068	srokopaq(:,:,k)=1.0
	1069	elsewhere
	1070	srokopaq(:,:,k)=0.0
	1071	endwhere
	1072	enddo
	1073
	1074	! ####################################################################################
	1075	! 3) Grid-box 3D OPAQ product fraction and cloud type cover (opaque/thin) + mean z_opaque
	1076	! ####################################################################################
	1077
	1078	do k= Nlevels,1,-1
	1079	do ic = 1, Ncolumns
	1080	do ip = 1, Npoints
	1081
	1082	cldlay(ip,ic,1) = MAX(cldlay(ip,ic,1),cldyopaq(ip,ic,k)) ! Opaque clouds
	1083	cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4),cldy(ip,ic,k)) ! All clouds
	1084
	1085	nsublay(ip,ic,1) = MAX(nsublay(ip,ic,1),srok(ip,ic,k))
	1086	nsublay(ip,ic,2) = MAX(nsublay(ip,ic,2),srok(ip,ic,k))
	1087	! nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
	1088	nsub(ip,k) = nsub(ip,k) + srok(ip,ic,k)
	1089	nsubopaq(ip,k) = nsubopaq(ip,k) + srokopaq(ip,ic,k)
	1090
	1091	enddo
	1092	enddo
	1093	enddo
	1094
	1095	! OPAQ variables
	1096	do ic = 1, Ncolumns
	1097	do ip = 1, Npoints
	1098
	1099	! Declaring non-opaque cloudy profiles as thin cloud profiles
	1100	if ( (cldlay(ip,ic,4) .eq. 1.0) .and. (cldlay(ip,ic,1) .eq. 0.0) ) then
	1101	cldlay(ip,ic,2) = 1.0
	1102	endif
	1103
	1104	! Filling in 3D and 2D variables
	1105
	1106	! Opaque cloud profiles
	1107	if ( cldlay(ip,ic,1) .eq. 1.0 ) then
	1108	zopac = 0.0
	1109	do k=2,Nlevels
	1110	! Declaring opaque cloud fraction and z_opaque altitude for 3D and 2D variables
	1111	if ( (cldy(ip,ic,k) .eq. 1.0) .and. (zopac .eq. 0.0) ) then
	1112	lidarcldtype(ip,k-1,3) = lidarcldtype(ip,k-1,3) + 1.0
	1113	cldlay(ip,ic,3) = vgrid_z(k-1) !z_opaque altitude
	1114	nsublay(ip,ic,3) = 1.0
	1115	zopac = 1.0
	1116	endif
	1117	if ( cldy(ip,ic,k) .eq. 1.0 ) then
	1118	lidarcldtype(ip,k,1) = lidarcldtype(ip,k,1) + 1.0
	1119	endif
	1120	enddo
	1121	endif
	1122
	1123	! Thin cloud profiles
	1124	if ( cldlay(ip,ic,2) .eq. 1.0 ) then
	1125	do k=1,Nlevels
	1126	! Declaring thin cloud fraction for 3D variable
	1127	if ( cldy(ip,ic,k) .eq. 1.0 ) then
	1128	lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1.0
	1129	endif
	1130	enddo
	1131	endif
	1132
	1133	enddo
	1134	enddo
	1135
	1136	! 3D cloud types fraction (opaque=1 and thin=2)
	1137	where ( nsub(:,:) .gt. 0.0 )
	1138	lidarcldtype(:,:,1) = lidarcldtype(:,:,1)/nsub(:,:)
	1139	lidarcldtype(:,:,2) = lidarcldtype(:,:,2)/nsub(:,:)
	1140	elsewhere
	1141	lidarcldtype(:,:,1) = undef
	1142	lidarcldtype(:,:,2) = undef
	1143	endwhere
	1144	! 3D z_opaque fraction (=3)
	1145	where ( nsubopaq(:,:) .gt. 0.0 )
	1146	lidarcldtype(:,:,3) = lidarcldtype(:,:,3)/nsubopaq(:,:)
	1147	elsewhere
	1148	lidarcldtype(:,:,3) = undef
	1149	endwhere
	1150	! 3D opacity fraction (=4) !Summing z_opaque fraction from TOA(k=Nlevels) to SFC(k=1)
	1151	lidarcldtype(:,Nlevels,4) = lidarcldtype(:,Nlevels,3)
	1152	do ip = 1, Npoints
	1153	do k = Nlevels-1, 1, -1
	1154	if ( lidarcldtype(ip,k,3) .ne. undef ) then
	1155	lidarcldtype(ip,k,4) = lidarcldtype(ip,k+1,4) + lidarcldtype(ip,k,3)
	1156	endif
	1157	enddo
	1158	enddo
	1159	where ( nsubopaq(:,:) .eq. 0.0 )
	1160	lidarcldtype(:,:,4) = undef
	1161	endwhere
	1162
	1163	! Layered cloud types (opaque, thin and z_opaque 2D variables)
	1164
	1165	do iz = 1, Ntype
	1166	do ic = 1, Ncolumns
	1167	cldtype(:,iz) = cldtype(:,iz) + cldlay(:,ic,iz)
	1168	nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz)
	1169	enddo
	1170	enddo
	1171	where (nsublayer(:,:) .gt. 0.0)
	1172	cldtype(:,:) = cldtype(:,:)/nsublayer(:,:)
	1173	elsewhere
	1174	cldtype(:,:) = undef
	1175	endwhere
	1176
	1177	END SUBROUTINE COSP_OPAQ
	1178	! END OF OPAQ CHANGES
	1179
	1180
[1262]	1181	END MODULE MOD_LMD_IPSL_STATS

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

Download in other formats:

Original Format