source: LMDZ4/branches/LMDZ4-dev-20091210/libf/cosp/lmd_ipsl_stats.F90 @ 1280

Last change on this file since 1280 was 1262, checked in by idelkadi, 15 years ago

-Rajout des fichiers cosp_input_nl.txt et cosp_output_nl.txt pour controler les parametres entrees sorties COSP
-Rajout du repertoire cosp contenant les sources de COSP

File size: 13.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!
56! Version 1.0 (June 2007)
57! Version 1.1 (May 2008)
58! Version 1.2 (June 2008)
59! Version 2.0 (October 2008)
60! Version 2.1 (December 2008)
61! c------------------------------------------------------------------------------------
62
63! c inputs :
64      integer npoints
65      integer ncol
66      integer llm
67      integer max_bin               ! nb of bins for SR CFADs
68      integer ncat                  ! nb of cloud layer types (low,mid,high,total)
69      integer nrefl                 ! nb of solar zenith angles for parasol reflectances
70
71      real undef                    ! undefined value
72      real pnorm(npoints,ncol,llm)  ! lidar ATB
73      real pmol(npoints,llm)        ! molecular ATB
74      real land(npoints)            ! Landmask [0 - Ocean, 1 - Land]   
75      real pplay(npoints,llm)       ! pressure on model levels (Pa)
76      logical ok_lidar_cfad         ! true if lidar CFAD diagnostics need to be computed
77      real refl(npoints,ncol,nrefl) ! subgrid parasol reflectance ! parasol
78
79! c outputs :
80      real lidarcld(npoints,llm)     ! 3D "lidar" cloud fraction
81      real cldlayer(npoints,ncat)    ! "lidar" cloud fraction (low, mid, high, total)
82      real cfad2(npoints,max_bin,llm) ! CFADs of SR 
83      real srbval(max_bin)           ! SR bins in CFADs 
84      real parasolrefl(npoints,nrefl)! grid-averaged parasol reflectance
85
86! c threshold for cloud detection :
87      real S_clr
88      parameter (S_clr = 1.2)
89      real S_cld
90!      parameter (S_cld = 3.0)  ! Previous thresold for cloud detection
91      parameter (S_cld = 5.0)  ! New (dec 2008) thresold for cloud detection
92      real S_att
93      parameter (S_att = 0.01)
94
95! c local variables :
96      integer ic,k
97      real x3d(npoints,ncol,llm)
98      real x3d_c(npoints,llm),pnorm_c(npoints,llm)
99      real xmax
100!
101! c -------------------------------------------------------
102! c 0- Initializations
103! c -------------------------------------------------------
104! Parasol reflectance algorithm is not valid over land. Write
105! a warning if there is no land. Landmask [0 - Ocean, 1 - Land]
106      IF ( MAXVAL(land(:)) .EQ. 0.0) THEN
107          WRITE (*,*) 'WARNING. PARASOL reflectance is not valid over land' &
108             & ,' and there is only land'
109      END IF
110
111!  Should be modified in future version
112      xmax=undef-1.0
113
114! c -------------------------------------------------------
115! c 1- Lidar scattering ratio :
116! c -------------------------------------------------------
117!
118!       where ((pnorm.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
119!          x3d = pnorm/pmol
120!       elsewhere
121!           x3d = undef
122!       end where
123! A.B-S: pmol reduced to 2D (npoints,llm) (Dec 08)
124      do ic = 1, ncol
125        pnorm_c = pnorm(:,ic,:)
126        where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
127            x3d_c = pnorm_c/pmol
128        elsewhere
129            x3d_c = undef
130        end where
131        x3d(:,ic,:) = x3d_c
132      enddo
133
134! c -------------------------------------------------------
135! c 2- Diagnose cloud fractions (3D, low, middle, high, total)
136! c from subgrid-scale lidar scattering ratios :
137! c -------------------------------------------------------
138
139      CALL COSP_CLDFRAC(npoints,ncol,llm,ncat,  &
140              x3d,pplay, S_att,S_cld,undef,lidarcld, &
141              cldlayer)
142
143! c -------------------------------------------------------
144! c 3- CFADs
145! c -------------------------------------------------------
146      if (ok_lidar_cfad) then
147!
148! c CFADs of subgrid-scale lidar scattering ratios :
149! c -------------------------------------------------------
150      CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin, &
151                 x3d, &
152                 S_att,S_clr,xmax,cfad2,srbval)
153
154      endif   ! ok_lidar_cfad
155! c -------------------------------------------------------
156
157! c -------------------------------------------------------
158! c 4- Compute grid-box averaged Parasol reflectances
159! c -------------------------------------------------------
160
161      parasolrefl(:,:) = 0.0
162
163      do k = 1, nrefl
164       do ic = 1, ncol
165         parasolrefl(:,k) = parasolrefl(:,k) + refl(:,ic,k)
166       enddo
167      enddo
168
169      do k = 1, nrefl
170        parasolrefl(:,k) = parasolrefl(:,k) / float(ncol)
171! if land=1 -> parasolrefl=undef
172! if land=0 -> parasolrefl=parasolrefl
173        parasolrefl(:,k) = parasolrefl(:,k) * MAX(1.0-land(:),0.0) &
174                           + (1.0 - MAX(1.0-land(:),0.0))*undef
175      enddo
176
177      RETURN
178      END SUBROUTINE diag_lidar
179         
180         
181!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
182!-------------------- FUNCTION COSP_CFAD_SR ------------------------
183! Author: Sandrine Bony (LMD/IPSL, CNRS, Paris)
184!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
185      SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins, &
186                      x,S_att,S_clr,xmax,cfad,srbval)
187      IMPLICIT NONE
188
189!--- Input arguments
190! Npoints: Number of horizontal points
191! Ncolumns: Number of subcolumns
192! Nlevels: Number of levels
193! Nbins: Number of x axis bins
194! xmax: maximum value allowed for x
195! S_att: Threshold for full attenuation
196! S_clr: Threshold for clear-sky layer
197!
198!--- Input-Outout arguments
199! x: variable to process (Npoints,Ncolumns,Nlevels), mofified where saturation occurs
200!
201! -- Output arguments
202! srbval : values of the histogram bins
203! cfad: 2D histogram on each horizontal point
204
205! Input arguments
206      integer Npoints,Ncolumns,Nlevels,Nbins
207      real xmax,S_att,S_clr
208! Input-output arguments
209      real x(Npoints,Ncolumns,Nlevels)
210! Output :
211      real cfad(Npoints,Nbins,Nlevels)
212      real srbval(Nbins)
213! Local variables
214      integer i, j, k, ib
215
216! c -------------------------------------------------------
217! c 0- Initializations
218! c -------------------------------------------------------
219      if ( Nbins .lt. 6) return
220
221      srbval(1) =  S_att
222      srbval(2) =  S_clr
223      srbval(3) =  3.0
224      srbval(4) =  5.0
225      srbval(5) =  7.0
226      srbval(6) = 10.0
227      do i = 7, MIN(10,Nbins)
228       srbval(i) = srbval(i-1) + 5.0
229      enddo
230      DO i = 11, MIN(13,Nbins)
231       srbval(i) = srbval(i-1) + 10.0
232      enddo
233      srbval(MIN(14,Nbins)) = 80.0
234      srbval(Nbins) = xmax
235      cfad(:,:,:) = 0.0
236
237
238! c -------------------------------------------------------
239! c c- Compute CFAD
240! c -------------------------------------------------------
241
242        do j = Nlevels, 1, -1
243          do k = 1, Ncolumns
244              where ( x(:,k,j).le.srbval(1) ) &
245                        cfad(:,1,j) = cfad(:,1,j) + 1.0
246          enddo  !k
247        enddo  !j
248
249      do ib = 2, Nbins
250        do j = Nlevels, 1, -1
251          do k = 1, Ncolumns
252              where ( ( x(:,k,j).gt.srbval(ib-1) ) .and. ( x(:,k,j).le.srbval(ib) ) ) &
253                        cfad(:,ib,j) = cfad(:,ib,j) + 1.0
254          enddo  !k
255        enddo  !j
256      enddo  !ib
257
258      cfad(:,:,:) = cfad(:,:,:) / float(Ncolumns)
259
260! c -------------------------------------------------------
261      RETURN
262      END SUBROUTINE COSP_CFAD_SR
263
264!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
265!-------------------- SUBROUTINE COSP_CLDFRAC -------------------
266! c Purpose: Cloud fraction diagnosed from lidar measurements
267!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
268      SUBROUTINE COSP_CLDFRAC(Npoints,Ncolumns,Nlevels,Ncat, &
269                  x,pplay,S_att,S_cld,undef,lidarcld, &
270                  cldlayer)
271      IMPLICIT NONE
272! Input arguments
273      integer Npoints,Ncolumns,Nlevels,Ncat
274      real x(Npoints,Ncolumns,Nlevels)
275      real pplay(Npoints,Nlevels)
276      real S_att,S_cld
277      real undef
278! Output :
279      real lidarcld(Npoints,Nlevels) ! 3D cloud fraction
280      real cldlayer(Npoints,Ncat)    ! low, middle, high, total cloud fractions
281! Local variables
282      integer ip, k, iz, ic
283      real p1
284      real cldy(Npoints,Ncolumns,Nlevels)
285      real srok(Npoints,Ncolumns,Nlevels)
286      real cldlay(Npoints,Ncolumns,Ncat)
287      real nsublay(Npoints,Ncolumns,Ncat), nsublayer(Npoints,Ncat)
288      real nsub(Npoints,Nlevels)
289
290
291! ---------------------------------------------------------------
292! 1- initialization
293! ---------------------------------------------------------------
294
295      if ( Ncat .ne. 4 ) then
296         print *,'Error in lmd_ipsl_stats.cosp_cldfrac, Ncat must be 4, not',Ncat
297         stop
298      endif
299
300      lidarcld = 0.0
301      nsub = 0.0
302      cldlay = 0.0
303      nsublay = 0.0
304
305! ---------------------------------------------------------------
306! 2- Cloud detection
307! ---------------------------------------------------------------
308
309      do k = 1, Nlevels
310
311! cloud detection at subgrid-scale:
312         where ( (x(:,:,k).gt.S_cld) .and. (x(:,:,k).ne. undef) )
313           cldy(:,:,k)=1.0
314         elsewhere
315           cldy(:,:,k)=0.0
316         endwhere
317
318! number of usefull sub-columns:
319         where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef)  )
320           srok(:,:,k)=1.0
321         elsewhere
322           srok(:,:,k)=0.0
323         endwhere
324
325      enddo ! k
326
327! ---------------------------------------------------------------
328! 3- grid-box 3D cloud fraction and layered cloud fractions (ISCCP pressure
329! categories) :
330! ---------------------------------------------------------------
331! Test abderr
332      do k = Nlevels, 1, -1
333       do ic = 1, Ncolumns
334        do ip = 1, Npoints
335          iz=1
336          p1 = pplay(ip,k)
337          if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
338            iz=3
339          else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid clouds
340            iz=2
341         endif
342
343         cldlay(ip,ic,iz) = MAX(cldlay(ip,ic,iz),cldy(ip,ic,k))
344         cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4),cldy(ip,ic,k))
345         lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k)
346
347         nsublay(ip,ic,iz) = MAX(nsublay(ip,ic,iz),srok(ip,ic,k))
348         nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
349         nsub(ip,k)=nsub(ip,k) + srok(ip,ic,k)
350
351        enddo
352       enddo
353      enddo
354
355! -- grid-box 3D cloud fraction
356
357      where ( nsub(:,:).gt.0.0 )
358         lidarcld(:,:) = lidarcld(:,:)/nsub(:,:)
359      elsewhere
360         lidarcld(:,:) = undef
361      endwhere
362
363! -- layered cloud fractions
364
365      cldlayer = 0.0
366      nsublayer = 0.0
367
368      do iz = 1, Ncat
369       do ic = 1, Ncolumns
370
371          cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz)   
372          nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz)
373
374       enddo
375      enddo
376      where ( nsublayer(:,:).gt.0.0 )
377         cldlayer(:,:) = cldlayer(:,:)/nsublayer(:,:)
378      elsewhere
379         cldlayer(:,:) = undef
380      endwhere
381
382      RETURN
383      END SUBROUTINE COSP_CLDFRAC
384! ---------------------------------------------------------------
385         
386END MODULE MOD_LMD_IPSL_STATS
Note: See TracBrowser for help on using the repository browser.