Context Navigation

lmd_ipsl_stats.F90 @ 3377

Last change on this file since 3377 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

Line
1	! Copyright (c) 2009, Centre National de la Recherche Scientifique
2	! All rights reserved.
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 $
5	!
6	! Redistribution and use in source and binary forms, with or without modification, are permitted
7	! provided that the following conditions are met:
8	!
9	! * Redistributions of source code must retain the above copyright notice, this list
10	! of conditions and the following disclaimer.
11	! * Redistributions in binary form must reproduce the above copyright notice, this list
12	! of conditions and the following disclaimer in the documentation and/or other materials
13	! provided with the distribution.
14	! * Neither the name of the LMD/IPSL/CNRS/UPMC nor the names of its
15	! contributors may be used to endorse or promote products derived from this software without
16	! specific prior written permission.
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
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 &
37	,tmp,pnorm,pnorm_perp,pmol,refl,land,pplay,undef,ok_lidar_cfad &
38	,cfad2,srbval,ncat,ntype,lidarcld,lidarcldtype,lidarcldphase,cldlayer & !OPAQ
39	,cldtype,cldlayerphase,lidarcldtmp,parasolrefl,vgrid_z,profSR) !OPAQ !TIBO
40	!
41	! -----------------------------------------------------------------------------------
42	! Lidar outputs :
43	!
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
47	! +
48	! Compute CFADs of lidar scattering ratio SR and of depolarization index
49	!
50	! Authors: Sandrine Bony and Helene Chepfer (LMD/IPSL, CNRS, UPMC, France).
51	!
52	! December 2008, S. Bony, H. Chepfer and J-L. Dufresne :
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)
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
64	!
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	! ------------------------------------------------------------------------------------
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 ntype ! nb of OPAQ products (opaque and thin clouds, z_opaque) !OPAQ
85	integer nrefl ! nb of solar zenith angles for parasol reflectances
86
87	real undef ! undefined value
88	real pnorm(npoints,ncol,llm) ! lidar ATB
89	real pmol(npoints,llm) ! molecular ATB
90	real land(npoints) ! Landmask [0 - Ocean, 1 - Land]
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
94	real tmp(npoints,llm) ! temp at each levels
95	real pnorm_perp(npoints,ncol,llm) ! lidar perpendicular ATB
96	real vgrid_z(llm) ! mid-level altitude of the output vertical grid !OPAQ
97
98	! c outputs :
99	real lidarcld(npoints,llm) ! 3D "lidar" cloud fraction
100	real lidarcldtype(npoints,llm,ntype+1) ! 3D "lidar" OPAQ type fraction + opacity !OPAQ
101	real sub(npoints,llm) ! 3D "lidar" indice
102	real cldlayer(npoints,ncat) ! "lidar" cloud layer fraction (low, mid, high, total)
103	real cldtype(npoints,ntype) ! "lidar" OPAQ type covers (opaque/thin cloud + z_opaque) !OPAQ
104
105	real cfad2(npoints,max_bin,llm) ! CFADs of SR
106	real srbval(max_bin) ! SR bins in CFADs
107	real parasolrefl(npoints,nrefl)! grid-averaged parasol reflectance
108	! real profSR(npoints,ncol,llm) ! tableau avec les subcolumns SR !TIBO
109	real profSR(npoints,llm,ncol) ! tableau avec les subcolumns SR !TIBO2
110
111	! c threshold for cloud detection :
112	real S_clr
113	parameter (S_clr = 1.2)
114	real S_cld
115	parameter (S_cld = 5.) ! Thresold for cloud detection
116	real S_att
117	parameter (S_att = 0.01)
118	! parameter (S_att = 0.06) !OPAQ ! Threshold for "surface detection" equivalent
119
120	! c local variables :
121	integer ic,k,i,j
122	real x3d(npoints,ncol,llm)
123	real x3d_c(npoints,llm),pnorm_c(npoints,llm)
124	real xmax
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
136	!
137	! c -------------------------------------------------------
138	! c 0- Initializations
139	! c -------------------------------------------------------
140	!
141	! Should be modified in future version
142	xmax=undef-1.0
143
144	! c -------------------------------------------------------
145	! c 1- Lidar scattering ratio :
146	! c -------------------------------------------------------
147
148	do ic = 1, ncol
149	pnorm_c = pnorm(:,ic,:)
150	where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
151	x3d_c = pnorm_c/pmol
152	elsewhere
153	x3d_c = undef
154	end where
155	x3d(:,ic,:) = x3d_c
156	! profSR(:,ic,:) = x3d(:,ic,:) !TIBO
157	profSR(:,:,ic) = x3d(:,ic,:) !TIBO2
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
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)
168
169	CALL COSP_OPAQ(npoints,ncol,llm,ntype,x3d,S_cld,undef,lidarcldtype, & !OPAQ
170	cldtype,vgrid_z) !OPAQ
171
172	! c -------------------------------------------------------
173	! c 3- CFADs
174	! c -------------------------------------------------------
175	if (ok_lidar_cfad) then
176	!
177	! c CFADs of subgrid-scale lidar scattering ratios :
178	! c -------------------------------------------------------
179	CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin,undef, &
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) &
203	+ (1.0 - MAX(1.0-land(:),0.0))*undef
204	enddo
205
206	RETURN
207	END SUBROUTINE diag_lidar
208
209
210	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
211	!-------------------- FUNCTION COSP_CFAD_SR ------------------------
212	! Author: Sandrine Bony (LMD/IPSL, CNRS, Paris)
213	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
214	SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins,undef, &
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
236	real xmax,S_att,S_clr,undef
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
244	real srbval_ext(0:Nbins)
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
267	srbval_ext(1:Nbins) = srbval
268	srbval_ext(0) = -1.0
269	! c -------------------------------------------------------
270	! c c- Compute CFAD
271	! c -------------------------------------------------------
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
286
287	where (cfad .ne. undef) cfad = cfad / float(Ncolumns)
288
289	! c -------------------------------------------------------
290	RETURN
291	END SUBROUTINE COSP_CFAD_SR
292
293
294	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
295	!-------------------- SUBROUTINE COSP_CLDFRAC -------------------
296	! c Purpose: Cloud fraction diagnosed from lidar measurements
297	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
303	IMPLICIT NONE
304	! Input arguments
305	integer Npoints,Ncolumns,Nlevels,Ncat
306	real x(Npoints,Ncolumns,Nlevels)
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
323	real pplay(Npoints,Nlevels)
324	real S_att,S_cld
325	real undef
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
332	real lidarcld(Npoints,Nlevels) ! 3D cloud fraction
333	real cldlayer(Npoints,Ncat) ! low, middle, high, total cloud fractions
334
335	! Local variables
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
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
359	#ifdef SYS_SX
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)
366	#endif
367
368
369
370
371	! ---------------------------------------------------------------
372	! 1- initialization
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
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
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:
432	where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef) )
433	srok(:,:,k)=1.0
434	elsewhere
435	srok(:,:,k)=0.0
436	endwhere
437
438	enddo ! k
439
440
441	! ---------------------------------------------------------------
442	! 3- grid-box 3D cloud fraction and layered cloud fractions (ISCCP pressure
443	! categories) :
444	! ---------------------------------------------------------------
445	lidarcld = 0.0
446	nsub = 0.0
447	#ifdef SYS_SX
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
457
458	do k = Nlevels, 1, -1
459	do ic = 1, Ncolumns
460	do ip = 1, Npoints
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
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))
491	endif
492
493	cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4), cldy(ip,ic,k))
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
500	cldlay(:,:,1) = cldlay1
501	cldlay(:,:,2) = cldlay2
502	cldlay(:,:,3) = cldlay3
503	nsublay(:,:,1) = nsublay1
504	nsublay(:,:,2) = nsublay2
505	nsublay(:,:,3) = nsublay3
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
512
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
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
570	cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz)
571	nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz)
572
573	enddo
574	enddo
575	where ( nsublayer(:,:).gt.0.0 )
576	cldlayer(:,:) = cldlayer(:,:)/nsublayer(:,:)
577	elsewhere
578	cldlayer(:,:) = undef
579	endwhere
580
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
994	END SUBROUTINE COSP_CLDFRAC
995	! ---------------------------------------------------------------
996
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)
1004
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
1181	END MODULE MOD_LMD_IPSL_STATS

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

Download in other formats:

Original Format