source: LMDZ5/branches/IPSLCM6.0.8/libf/phylmd/cosp/lmd_ipsl_stats.F90 @ 5455

Last change on this file since 5455 was 2435, checked in by Laurent Fairhead, 9 years ago

Merged trunk changes r2396:2434 into testing branch
  • 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!------------------------------------------------------------------------------------
31MODULE MOD_LMD_IPSL_STATS
32  USE MOD_LLNL_STATS
33  IMPLICIT NONE
34
35CONTAINS
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^5*alpha50 + ATB^4*beta50 + ATB^3*gamma50 + ATB^2*delta50 + 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! ---------------------------------------------------------------
575do ncol=1,Ncolumns
576do 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^5*alpha50 + ATB^4*beta50 + ATB^3*gamma50 + ATB^2*delta50
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
80399 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!
811if(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
829endif
830     
831      toplvlsat=0
832
833enddo
834enddo
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
846lidarcldphasetmp(:,:)=lidarcldphase(:,:,1)+lidarcldphase(:,:,2);
847WHERE (lidarcldphasetmp(:,:).gt. 0.)
848   lidarcldphase(:,:,6)=lidarcldphase(:,:,1)/lidarcldphasetmp(:,:)
849ELSEWHERE
850   lidarcldphase(:,:,6) = undef
851ENDWHERE
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)
888cldlayerphasetmp(:,:)=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
928do nlev=1,Nlevels
929do ncol=1,Ncolumns     
930do i=1,Npoints
931do itemp=1,Ntemp
932if(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
936elseif(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
940elseif(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
944endif
945enddo
946enddo
947enddo
948enddo
949
950! Check temperature cloud fraction
951do i=1,Npoints
952   do itemp=1,Ntemp
953checktemp=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
961enddo
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
967WHERE(sumlidarcldtemp(:,:)>0.)
968  lidarcldtemp(:,:,5)=lidarcldtemp(:,:,2)/sumlidarcldtemp(:,:)
969ELSEWHERE
970  lidarcldtemp(:,:,5)=undef
971ENDWHERE
972
973do i=1,4
974  WHERE(lidarcldtempind(:,:).gt.0.)
975     lidarcldtemp(:,:,i) = lidarcldtemp(:,:,i)/lidarcldtempind(:,:)
976  ELSEWHERE
977     lidarcldtemp(:,:,i) = undef
978  ENDWHERE
979enddo
980
981       RETURN
982      END SUBROUTINE COSP_CLDFRAC
983! ---------------------------------------------------------------
984
985
986END MODULE MOD_LMD_IPSL_STATS
Note: See TracBrowser for help on using the repository browser.