My Project
 All Classes Files Functions Variables Macros
lmd_ipsl_stats.F90
Go to the documentation of this file.
1 ! Copyright (c) 2009, Centre National de la Recherche Scientifique
2 ! All rights reserved.
3 !
4 ! Redistribution and use in source and binary forms, with or without modification, are permitted
5 ! provided that the following conditions are met:
6 !
7 ! * Redistributions of source code must retain the above copyright notice, this list
8 ! of conditions and the following disclaimer.
9 ! * Redistributions in binary form must reproduce the above copyright notice, this list
10 ! of conditions and the following disclaimer in the documentation and/or other materials
11 ! provided with the distribution.
12 ! * Neither the name of the LMD/IPSL/CNRS/UPMC nor the names of its
13 ! contributors may be used to endorse or promote products derived from this software without
14 ! specific prior written permission.
15 !
16 ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
17 ! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
18 ! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
19 ! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20 ! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21 ! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
22 ! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
23 ! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24 
25 
26 !------------------------------------------------------------------------------------
27 ! Authors: Sandrine Bony and Helene Chepfer (LMD/IPSL, CNRS, UPMC, France).
28 !------------------------------------------------------------------------------------
29 MODULE mod_lmd_ipsl_stats
30  USE mod_llnl_stats
31  IMPLICIT NONE
32 
33 CONTAINS
34  SUBROUTINE diag_lidar(npoints,ncol,llm,max_bin,nrefl &
35  ,pnorm,pmol,refl,land,pplay,undef,ok_lidar_cfad &
36  ,cfad2,srbval &
37  ,ncat,lidarcld,cldlayer,parasolrefl)
38 !
39 ! -----------------------------------------------------------------------------------
40 ! Lidar outputs :
41 !
42 ! Diagnose cloud fraction (3D cloud fraction + low/middle/high/total cloud fraction
43 ! from the lidar signals (ATB and molecular ATB) computed from model outputs
44 ! +
45 ! Compute CFADs of lidar scattering ratio SR and of depolarization index
46 !
47 ! Authors: Sandrine Bony and Helene Chepfer (LMD/IPSL, CNRS, UPMC, France).
48 !
49 ! December 2008, S. Bony, H. Chepfer and J-L. Dufresne :
50 ! - change of the cloud detection threshold S_cld from 3 to 5, for better
51 ! with both day and night observations. The optical thinest clouds are missed.
52 ! - remove of the detection of the first fully attenuated layer encountered from above.
53 ! December 2008, A. Bodas-Salcedo:
54 ! - Dimensions of pmol reduced to (npoints,llm)
55 ! August 2009, A. Bodas-Salcedo:
56 ! - Warning message regarding PARASOL being valid only over ocean deleted.
57 ! February 2010, A. Bodas-Salcedo:
58 ! - Undef passed into cosp_cfad_sr
59 ! June 2010, T. Yokohata, T. Nishimura and K. Ogochi
60 ! Optimisation of COSP_CFAD_SR
61 !
62 ! Version 1.0 (June 2007)
63 ! Version 1.1 (May 2008)
64 ! Version 1.2 (June 2008)
65 ! Version 2.0 (October 2008)
66 ! Version 2.1 (December 2008)
67 ! c------------------------------------------------------------------------------------
68 
69 ! c inputs :
70  integer npoints
71  integer ncol
72  integer llm
73  integer max_bin ! nb of bins for SR CFADs
74  integer ncat ! nb of cloud layer types (low,mid,high,total)
75  integer nrefl ! nb of solar zenith angles for parasol reflectances
76 
77  real undef ! undefined value
78  real pnorm(npoints,ncol,llm) ! lidar ATB
79  real pmol(npoints,llm) ! molecular ATB
80  real land(npoints) ! Landmask [0 - Ocean, 1 - Land]
81  real pplay(npoints,llm) ! pressure on model levels (Pa)
82  logical ok_lidar_cfad ! true if lidar CFAD diagnostics need to be computed
83  real refl(npoints,ncol,nrefl) ! subgrid parasol reflectance ! parasol
84 
85 ! c outputs :
86  real lidarcld(npoints,llm) ! 3D "lidar" cloud fraction
87  real cldlayer(npoints,ncat) ! "lidar" cloud fraction (low, mid, high, total)
88  real cfad2(npoints,max_bin,llm) ! CFADs of SR
89  real srbval(max_bin) ! SR bins in CFADs
90  real parasolrefl(npoints,nrefl)! grid-averaged parasol reflectance
91 
92 ! c threshold for cloud detection :
93  real s_clr
94  parameter(s_clr = 1.2)
95  real s_cld
96 ! parameter (S_cld = 3.0) ! Previous thresold for cloud detection
97  parameter(s_cld = 5.0) ! New (dec 2008) thresold for cloud detection
98  real s_att
99  parameter(s_att = 0.01)
100 
101 ! c local variables :
102  integer ic,k
103  real x3d(npoints,ncol,llm)
104  real x3d_c(npoints,llm),pnorm_c(npoints,llm)
105  real xmax
106 !
107 ! c -------------------------------------------------------
108 ! c 0- Initializations
109 ! c -------------------------------------------------------
110 !
111 
112 ! Should be modified in future version
113  xmax=undef-1.0
114 
115 ! c -------------------------------------------------------
116 ! c 1- Lidar scattering ratio :
117 ! c -------------------------------------------------------
118 !
119 ! where ((pnorm.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
120 ! x3d = pnorm/pmol
121 ! elsewhere
122 ! x3d = undef
123 ! end where
124 ! A.B-S: pmol reduced to 2D (npoints,llm) (Dec 08)
125  do ic = 1, ncol
126  pnorm_c = pnorm(:,ic,:)
127  where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
128  x3d_c = pnorm_c/pmol
129  elsewhere
130  x3d_c = undef
131  end where
132  x3d(:,ic,:) = x3d_c
133  enddo
134 
135 ! c -------------------------------------------------------
136 ! c 2- Diagnose cloud fractions (3D, low, middle, high, total)
137 ! c from subgrid-scale lidar scattering ratios :
138 ! c -------------------------------------------------------
139 
140  CALL cosp_cldfrac(npoints,ncol,llm,ncat, &
141  x3d,pplay, s_att,s_cld,undef,lidarcld, &
142  cldlayer)
143 
144 ! c -------------------------------------------------------
145 ! c 3- CFADs
146 ! c -------------------------------------------------------
147  if (ok_lidar_cfad) then
148 !
149 ! c CFADs of subgrid-scale lidar scattering ratios :
150 ! c -------------------------------------------------------
151  CALL cosp_cfad_sr(npoints,ncol,llm,max_bin,undef, &
152  x3d, &
153  s_att,s_clr,xmax,cfad2,srbval)
154 
155  endif ! ok_lidar_cfad
156 ! c -------------------------------------------------------
157 
158 ! c -------------------------------------------------------
159 ! c 4- Compute grid-box averaged Parasol reflectances
160 ! c -------------------------------------------------------
161 
162  parasolrefl(:,:) = 0.0
163 
164  do k = 1, nrefl
165  do ic = 1, ncol
166  parasolrefl(:,k) = parasolrefl(:,k) + refl(:,ic,k)
167  enddo
168  enddo
169 
170  do k = 1, nrefl
171  parasolrefl(:,k) = parasolrefl(:,k) / float(ncol)
172 ! if land=1 -> parasolrefl=undef
173 ! if land=0 -> parasolrefl=parasolrefl
174  parasolrefl(:,k) = parasolrefl(:,k) * max(1.0-land(:),0.0) &
175  + (1.0 - max(1.0-land(:),0.0))*undef
176  enddo
177 
178  RETURN
179  END SUBROUTINE diag_lidar
180 
181 
182 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
183 !-------------------- FUNCTION COSP_CFAD_SR ------------------------
184 ! Author: Sandrine Bony (LMD/IPSL, CNRS, Paris)
185 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
186  SUBROUTINE cosp_cfad_sr(Npoints,Ncolumns,Nlevels,Nbins,undef, &
187  x,s_att,s_clr,xmax,cfad,srbval)
188  IMPLICIT NONE
189 
190 !--- Input arguments
191 ! Npoints: Number of horizontal points
192 ! Ncolumns: Number of subcolumns
193 ! Nlevels: Number of levels
194 ! Nbins: Number of x axis bins
195 ! xmax: maximum value allowed for x
196 ! S_att: Threshold for full attenuation
197 ! S_clr: Threshold for clear-sky layer
198 !
199 !--- Input-Outout arguments
200 ! x: variable to process (Npoints,Ncolumns,Nlevels), mofified where saturation occurs
201 !
202 ! -- Output arguments
203 ! srbval : values of the histogram bins
204 ! cfad: 2D histogram on each horizontal point
205 
206 ! Input arguments
207  integer npoints,ncolumns,nlevels,nbins
208  real xmax,s_att,s_clr,undef
209 ! Input-output arguments
210  real x(npoints,ncolumns,nlevels)
211 ! Output :
212  real cfad(npoints,nbins,nlevels)
213  real srbval(nbins)
214 ! Local variables
215  integer i, j, k, ib
216  real srbval_ext(0:nbins)
217 
218 ! c -------------------------------------------------------
219 ! c 0- Initializations
220 ! c -------------------------------------------------------
221  if ( nbins .lt. 6) return
222 
223  srbval(1) = s_att
224  srbval(2) = s_clr
225  srbval(3) = 3.0
226  srbval(4) = 5.0
227  srbval(5) = 7.0
228  srbval(6) = 10.0
229  do i = 7, min(10,nbins)
230  srbval(i) = srbval(i-1) + 5.0
231  enddo
232  DO i = 11, min(13,nbins)
233  srbval(i) = srbval(i-1) + 10.0
234  enddo
235  srbval(min(14,nbins)) = 80.0
236  srbval(nbins) = xmax
237  cfad(:,:,:) = 0.0
238 
239  srbval_ext(1:nbins) = srbval
240  srbval_ext(0) = -1.0
241 ! c -------------------------------------------------------
242 ! c c- Compute CFAD
243 ! c -------------------------------------------------------
244 
245  do j = 1, nlevels
246  do ib = 1, nbins
247  do k = 1, ncolumns
248  do i = 1, npoints
249  if (x(i,k,j) /= undef) then
250  if ((x(i,k,j).gt.srbval_ext(ib-1)).and.(x(i,k,j).le.srbval_ext(ib))) &
251  cfad(i,ib,j) = cfad(i,ib,j) + 1.0
252  else
253  cfad(i,ib,j) = undef
254  endif
255  enddo
256  enddo
257  enddo
258  enddo
259 
260  where (cfad .ne. undef) cfad = cfad / float(ncolumns)
261 
262 ! c -------------------------------------------------------
263  RETURN
264  END SUBROUTINE cosp_cfad_sr
265 
266 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
267 !-------------------- SUBROUTINE COSP_CLDFRAC -------------------
268 ! c Purpose: Cloud fraction diagnosed from lidar measurements
269 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
270  SUBROUTINE cosp_cldfrac(Npoints,Ncolumns,Nlevels,Ncat, &
271  x,pplay,s_att,s_cld,undef,lidarcld, &
272  cldlayer)
273  IMPLICIT NONE
274 ! Input arguments
275  integer npoints,ncolumns,nlevels,ncat
276  real x(npoints,ncolumns,nlevels)
277  real pplay(npoints,nlevels)
278  real s_att,s_cld
279  real undef
280 ! Output :
281  real lidarcld(npoints,nlevels) ! 3D cloud fraction
282  real cldlayer(npoints,ncat) ! low, middle, high, total cloud fractions
283 ! Local variables
284  integer ip, k, iz, ic
285  real p1
286  real cldy(npoints,ncolumns,nlevels)
287  real srok(npoints,ncolumns,nlevels)
288  real cldlay(npoints,ncolumns,ncat)
289  real nsublay(npoints,ncolumns,ncat), nsublayer(npoints,ncat)
290  real nsub(npoints,nlevels)
291 
292  real cldlay1(npoints,ncolumns)
293  real cldlay2(npoints,ncolumns)
294  real cldlay3(npoints,ncolumns)
295  real nsublay1(npoints,ncolumns)
296  real nsublay2(npoints,ncolumns)
297  real nsublay3(npoints,ncolumns)
298 
299 
300 ! ---------------------------------------------------------------
301 ! 1- initialization
302 ! ---------------------------------------------------------------
303 
304  if ( ncat .ne. 4 ) then
305  print *,'Error in lmd_ipsl_stats.cosp_cldfrac, Ncat must be 4, not',ncat
306  stop
307  endif
308 
309  lidarcld = 0.0
310  nsub = 0.0
311  cldlay = 0.0
312  nsublay = 0.0
313 
314 ! ---------------------------------------------------------------
315 ! 2- Cloud detection
316 ! ---------------------------------------------------------------
317 
318  do k = 1, nlevels
319 
320 ! cloud detection at subgrid-scale:
321  where ( (x(:,:,k).gt.s_cld) .and. (x(:,:,k).ne. undef) )
322  cldy(:,:,k)=1.0
323  elsewhere
324  cldy(:,:,k)=0.0
325  endwhere
326 
327 ! number of usefull sub-columns:
328  where ( (x(:,:,k).gt.s_att) .and. (x(:,:,k).ne. undef) )
329  srok(:,:,k)=1.0
330  elsewhere
331  srok(:,:,k)=0.0
332  endwhere
333 
334  enddo ! k
335 
336 ! ---------------------------------------------------------------
337 ! 3- grid-box 3D cloud fraction and layered cloud fractions (ISCCP pressure
338 ! categories) :
339 ! ---------------------------------------------------------------
340  lidarcld = 0.0
341  nsub = 0.0
342 
343 !! XXX: Use cldlay[1-3] and nsublay[1-3] to avoid bank-conflicts.
344  cldlay1 = 0.0
345  cldlay2 = 0.0
346  cldlay3 = 0.0
347  cldlay(:,:,4) = 0.0 !! XXX: Ncat == 4
348  nsublay1 = 0.0
349  nsublay2 = 0.0
350  nsublay3 = 0.0
351  nsublay(:,:,4) = 0.0
352  do k = nlevels, 1, -1
353  do ic = 1, ncolumns
354  do ip = 1, npoints
355  p1 = pplay(ip,k)
356 
357  if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
358  cldlay3(ip,ic) = max(cldlay3(ip,ic), cldy(ip,ic,k))
359  nsublay3(ip,ic) = max(nsublay3(ip,ic), srok(ip,ic,k))
360  else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid clouds
361  cldlay2(ip,ic) = max(cldlay2(ip,ic), cldy(ip,ic,k))
362  nsublay2(ip,ic) = max(nsublay2(ip,ic), srok(ip,ic,k))
363  else
364  cldlay1(ip,ic) = max(cldlay1(ip,ic), cldy(ip,ic,k))
365  nsublay1(ip,ic) = max(nsublay1(ip,ic), srok(ip,ic,k))
366  endif
367 
368  cldlay(ip,ic,4) = max(cldlay(ip,ic,4), cldy(ip,ic,k))
369  lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k)
370  nsublay(ip,ic,4) = max(nsublay(ip,ic,4),srok(ip,ic,k))
371  nsub(ip,k)=nsub(ip,k) + srok(ip,ic,k)
372  enddo
373  enddo
374  enddo
375  cldlay(:,:,1) = cldlay1
376  cldlay(:,:,2) = cldlay2
377  cldlay(:,:,3) = cldlay3
378  nsublay(:,:,1) = nsublay1
379  nsublay(:,:,2) = nsublay2
380  nsublay(:,:,3) = nsublay3
381 
382 ! -- grid-box 3D cloud fraction
383 
384  where ( nsub(:,:).gt.0.0 )
385  lidarcld(:,:) = lidarcld(:,:)/nsub(:,:)
386  elsewhere
387  lidarcld(:,:) = undef
388  endwhere
389 
390 ! -- layered cloud fractions
391 
392  cldlayer = 0.0
393  nsublayer = 0.0
394 
395  do iz = 1, ncat
396  do ic = 1, ncolumns
397 
398  cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz)
399  nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz)
400 
401  enddo
402  enddo
403  where ( nsublayer(:,:).gt.0.0 )
404  cldlayer(:,:) = cldlayer(:,:)/nsublayer(:,:)
405  elsewhere
406  cldlayer(:,:) = undef
407  endwhere
408 
409  RETURN
410  END SUBROUTINE cosp_cldfrac
411 ! ---------------------------------------------------------------
412 
413 END MODULE mod_lmd_ipsl_stats