Context Navigation

lmd_ipsl_stats.F90 @ 5360

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

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,lidarcld,lidarcldphase,cldlayer,cldlayerphase &
39	,lidarcldtmp,parasolrefl)
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 nrefl ! nb of solar zenith angles for parasol reflectances
85
86	real undef ! undefined value
87	real pnorm(npoints,ncol,llm) ! lidar ATB
88	real pmol(npoints,llm) ! molecular ATB
89	real land(npoints) ! Landmask [0 - Ocean, 1 - Land]
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
93	real tmp(npoints,llm) ! temp at each levels
94	real pnorm_perp(npoints,ncol,llm) ! lidar perpendicular ATB
95
96	! c outputs :
97	real lidarcld(npoints,llm) ! 3D "lidar" cloud fraction
98	real sub(npoints,llm) ! 3D "lidar" indice
99	real cldlayer(npoints,ncat) ! "lidar" cloud layer fraction (low, mid, high, total)
100
101	real cfad2(npoints,max_bin,llm) ! CFADs of SR
102	real srbval(max_bin) ! SR bins in CFADs
103	real parasolrefl(npoints,nrefl)! grid-averaged parasol reflectance
104
105	! c threshold for cloud detection :
106	real S_clr
107	parameter (S_clr = 1.2)
108	real S_cld
109	parameter (S_cld = 5.) ! Thresold for cloud detection
110	real S_att
111	parameter (S_att = 0.01)
112
113	! c local variables :
114	integer ic,k,i,j
115	real x3d(npoints,ncol,llm)
116	real x3d_c(npoints,llm),pnorm_c(npoints,llm)
117	real xmax
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
129	!
130	! c -------------------------------------------------------
131	! c 0- Initializations
132	! c -------------------------------------------------------
133	!
134	! Should be modified in future version
135	xmax=undef-1.0
136
137	! c -------------------------------------------------------
138	! c 1- Lidar scattering ratio :
139	! c -------------------------------------------------------
140
141	do ic = 1, ncol
142	pnorm_c = pnorm(:,ic,:)
143	where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
144	x3d_c = pnorm_c/pmol
145	elsewhere
146	x3d_c = undef
147	end where
148	x3d(:,ic,:) = x3d_c
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
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)
159
160	! c -------------------------------------------------------
161	! c 3- CFADs
162	! c -------------------------------------------------------
163	if (ok_lidar_cfad) then
164	!
165	! c CFADs of subgrid-scale lidar scattering ratios :
166	! c -------------------------------------------------------
167	CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin,undef, &
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) &
191	+ (1.0 - MAX(1.0-land(:),0.0))*undef
192	enddo
193
194	RETURN
195	END SUBROUTINE diag_lidar
196
197
198	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
199	!-------------------- FUNCTION COSP_CFAD_SR ------------------------
200	! Author: Sandrine Bony (LMD/IPSL, CNRS, Paris)
201	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
202	SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins,undef, &
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
224	real xmax,S_att,S_clr,undef
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
232	real srbval_ext(0:Nbins)
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
255	srbval_ext(1:Nbins) = srbval
256	srbval_ext(0) = -1.0
257	! c -------------------------------------------------------
258	! c c- Compute CFAD
259	! c -------------------------------------------------------
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
274
275	where (cfad .ne. undef) cfad = cfad / float(Ncolumns)
276
277	! c -------------------------------------------------------
278	RETURN
279	END SUBROUTINE COSP_CFAD_SR
280
281
282	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
283	!-------------------- SUBROUTINE COSP_CLDFRAC -------------------
284	! c Purpose: Cloud fraction diagnosed from lidar measurements
285	!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
291	IMPLICIT NONE
292	! Input arguments
293	integer Npoints,Ncolumns,Nlevels,Ncat
294	real x(Npoints,Ncolumns,Nlevels)
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
311	real pplay(Npoints,Nlevels)
312	real S_att,S_cld
313	real undef
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
320	real lidarcld(Npoints,Nlevels) ! 3D cloud fraction
321	real cldlayer(Npoints,Ncat) ! low, middle, high, total cloud fractions
322
323	! Local variables
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
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
347	#ifdef SYS_SX
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)
354	#endif
355
356
357
358
359	! ---------------------------------------------------------------
360	! 1- initialization
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
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
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:
420	where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef) )
421	srok(:,:,k)=1.0
422	elsewhere
423	srok(:,:,k)=0.0
424	endwhere
425
426	enddo ! k
427
428
429	! ---------------------------------------------------------------
430	! 3- grid-box 3D cloud fraction and layered cloud fractions (ISCCP pressure
431	! categories) :
432	! ---------------------------------------------------------------
433	lidarcld = 0.0
434	nsub = 0.0
435	#ifdef SYS_SX
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
445
446	do k = Nlevels, 1, -1
447	do ic = 1, Ncolumns
448	do ip = 1, Npoints
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
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))
479	endif
480
481	cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4), cldy(ip,ic,k))
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
488	cldlay(:,:,1) = cldlay1
489	cldlay(:,:,2) = cldlay2
490	cldlay(:,:,3) = cldlay3
491	nsublay(:,:,1) = nsublay1
492	nsublay(:,:,2) = nsublay2
493	nsublay(:,:,3) = nsublay3
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
500
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
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
558	cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz)
559	nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz)
560
561	enddo
562	enddo
563	where ( nsublayer(:,:).gt.0.0 )
564	cldlayer(:,:) = cldlayer(:,:)/nsublayer(:,:)
565	elsewhere
566	cldlayer(:,:) = undef
567	endwhere
568
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
982	END SUBROUTINE COSP_CLDFRAC
983	! ---------------------------------------------------------------
984
985
986	END MODULE MOD_LMD_IPSL_STATS

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

Download in other formats:

Original Format