Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

lmd_ipsl_stats.F90 @ 3331

Last change on this file since 3331 was 3331, checked in by acozic, 6 years ago
Add modification for isotopes
Property svn:executable set to ``*
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