source: LMDZ4/trunk/libf/cosp/lmd_ipsl_stats.F90 @ 3782

Last change on this file since 3782 was 1414, checked in by idelkadi, 14 years ago

Passage a la version cosp.v1.3 pour le Lidar et ISCCP
Corrections de bugs pour ISCCP et optimisation pour le Lidar

File size: 14.3 KB
Line 
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!------------------------------------------------------------------------------------
29MODULE MOD_LMD_IPSL_STATS
30  USE MOD_LLNL_STATS
31  IMPLICIT NONE
32
33CONTAINS
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
413END MODULE MOD_LMD_IPSL_STATS
Note: See TracBrowser for help on using the repository browser.