source: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/cosp/mod_lmd_ipsl_stats.F90 @ 3373

Last change on this file since 3373 was 3242, checked in by idelkadi, 7 years ago

Nettoyage du code :

  • changement de nomes des routines pour avoir les memes nome des modules
  • corrections
  • dos2unix appliquee aux fichiers
  • Rajouter "implicit none" pour les routines ne le comportant pas
  • Suppression des commentaires unitils

Corrections pour tourner en mode debug :

controle par la cles ok_debug_cosp dans mod_cosp_constants.F90
la valeur par defaut etant .false.

  • 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: 44.6 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,ntype,lidarcld,lidarcldtype,lidarcldphase,cldlayer & !OPAQ
39                  ,cldtype,cldlayerphase,lidarcldtmp,parasolrefl,vgrid_z,profSR)          !OPAQ !TIBO
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 ntype                 ! nb of OPAQ products (opaque and thin clouds, z_opaque) !OPAQ
85      integer nrefl                 ! nb of solar zenith angles for parasol reflectances
86
87      real undef                    ! undefined value
88      real pnorm(npoints,ncol,llm)  ! lidar ATB
89      real pmol(npoints,llm)        ! molecular ATB
90      real land(npoints)            ! Landmask [0 - Ocean, 1 - Land]
91      real pplay(npoints,llm)       ! pressure on model levels (Pa)
92      logical ok_lidar_cfad         ! true if lidar CFAD diagnostics need to be computed
93      real refl(npoints,ncol,nrefl) ! subgrid parasol reflectance ! parasol
94      real tmp(npoints,llm)         ! temp at each levels
95      real pnorm_perp(npoints,ncol,llm)  ! lidar perpendicular ATB
96      real vgrid_z(llm)             ! mid-level altitude of the output vertical grid         !OPAQ
97
98! c outputs :
99      real lidarcld(npoints,llm)     ! 3D "lidar" cloud fraction
100      real lidarcldtype(npoints,llm,ntype+1)   ! 3D "lidar" OPAQ type fraction + opacity     !OPAQ
101      real sub(npoints,llm)     ! 3D "lidar" indice
102      real cldlayer(npoints,ncat)    ! "lidar" cloud layer fraction (low, mid, high, total)
103      real cldtype(npoints,ntype)  ! "lidar" OPAQ type covers (opaque/thin cloud + z_opaque) !OPAQ
104
105      real cfad2(npoints,max_bin,llm) ! CFADs of SR
106      real srbval(max_bin)           ! SR bins in CFADs
107      real parasolrefl(npoints,nrefl)! grid-averaged parasol reflectance
108!     real profSR(npoints,ncol,llm)  ! tableau avec les subcolumns SR !TIBO
109      real profSR(npoints,llm,ncol)  ! tableau avec les subcolumns SR !TIBO2
110
111! c threshold for cloud detection :
112      real S_clr
113      parameter (S_clr = 1.2)
114      real S_cld
115      parameter (S_cld = 5.)  ! Thresold for cloud detection
116      real S_att
117      parameter (S_att = 0.01)
118!      parameter (S_att = 0.06)  !OPAQ ! Threshold for "surface detection" equivalent
119
120! c local variables :
121      integer ic,k,i,j
122      real x3d(npoints,ncol,llm)
123      real x3d_c(npoints,llm),pnorm_c(npoints,llm)
124      real xmax
125
126! Output variables
127      integer,parameter :: nphase = 6 ! nb of cloud layer phase types (ice,liquid,undefined,false ice,false liquid,Percent of ice)
128      real lidarcldphase(npoints,llm,nphase)   ! 3D "lidar" phase cloud fraction
129      real lidarcldtmp(npoints,40,5)          ! 3D "lidar" phase cloud fraction as a function of temp
130      real cldlayerphase(npoints,ncat,nphase)  ! "lidar" phase low mid high cloud fraction
131
132! SR detection threshold
133      real, parameter  ::  S_cld_att = 30. ! New threshold for undefine cloud phase detection   
134
135
136!
137! c -------------------------------------------------------
138! c 0- Initializations
139! c -------------------------------------------------------
140!
141!  Should be modified in future version
142      xmax=undef-1.0
143
144! c -------------------------------------------------------
145! c 1- Lidar scattering ratio :
146! c -------------------------------------------------------
147
148      do ic = 1, ncol
149        pnorm_c = pnorm(:,ic,:)
150        where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
151            x3d_c = pnorm_c/pmol
152        elsewhere
153            x3d_c = undef
154        end where
155         x3d(:,ic,:) = x3d_c
156!       profSR(:,ic,:) = x3d(:,ic,:) !TIBO
157        profSR(:,:,ic) = x3d(:,ic,:) !TIBO2
158      enddo
159
160! c -------------------------------------------------------
161! c 2- Diagnose cloud fractions (3D, low, middle, high, total)
162! c from subgrid-scale lidar scattering ratios :
163! c -------------------------------------------------------
164
165    CALL COSP_CLDFRAC(npoints,ncol,llm,ncat,nphase,  &
166              tmp,x3d,pnorm,pnorm_perp,pplay, S_att,S_cld,S_cld_att,undef,lidarcld, &
167              cldlayer,lidarcldphase,sub,cldlayerphase,lidarcldtmp)
168
169    CALL COSP_OPAQ(npoints,ncol,llm,ntype,x3d,S_cld,undef,lidarcldtype,            & !OPAQ
170                   cldtype,vgrid_z)                                                  !OPAQ
171
172! c -------------------------------------------------------
173! c 3- CFADs
174! c -------------------------------------------------------
175      if (ok_lidar_cfad) then
176!
177! c CFADs of subgrid-scale lidar scattering ratios :
178! c -------------------------------------------------------
179      CALL COSP_CFAD_SR(npoints,ncol,llm,max_bin,undef, &
180                 x3d, &
181                 S_att,S_clr,xmax,cfad2,srbval)
182
183      endif   ! ok_lidar_cfad
184! c -------------------------------------------------------
185
186! c -------------------------------------------------------
187! c 4- Compute grid-box averaged Parasol reflectances
188! c -------------------------------------------------------
189
190      parasolrefl(:,:) = 0.0
191
192      do k = 1, nrefl
193       do ic = 1, ncol
194         parasolrefl(:,k) = parasolrefl(:,k) + refl(:,ic,k)
195       enddo
196      enddo
197
198      do k = 1, nrefl
199        parasolrefl(:,k) = parasolrefl(:,k) / float(ncol)
200! if land=1 -> parasolrefl=undef
201! if land=0 -> parasolrefl=parasolrefl
202        parasolrefl(:,k) = parasolrefl(:,k) * MAX(1.0-land(:),0.0) &
203                           + (1.0 - MAX(1.0-land(:),0.0))*undef
204      enddo
205
206      RETURN
207      END SUBROUTINE diag_lidar
208
209
210!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
211!-------------------- FUNCTION COSP_CFAD_SR ------------------------
212! Author: Sandrine Bony (LMD/IPSL, CNRS, Paris)
213!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
214      SUBROUTINE COSP_CFAD_SR(Npoints,Ncolumns,Nlevels,Nbins,undef, &
215                      x,S_att,S_clr,xmax,cfad,srbval)
216      IMPLICIT NONE
217
218!--- Input arguments
219! Npoints: Number of horizontal points
220! Ncolumns: Number of subcolumns
221! Nlevels: Number of levels
222! Nbins: Number of x axis bins
223! xmax: maximum value allowed for x
224! S_att: Threshold for full attenuation
225! S_clr: Threshold for clear-sky layer
226!
227!--- Input-Outout arguments
228! x: variable to process (Npoints,Ncolumns,Nlevels), mofified where saturation occurs
229!
230! -- Output arguments
231! srbval : values of the histogram bins
232! cfad: 2D histogram on each horizontal point
233
234! Input arguments
235      integer Npoints,Ncolumns,Nlevels,Nbins
236      real xmax,S_att,S_clr,undef
237! Input-output arguments
238      real x(Npoints,Ncolumns,Nlevels)
239! Output :
240      real cfad(Npoints,Nbins,Nlevels)
241      real srbval(Nbins)
242! Local variables
243      integer i, j, k, ib
244      real srbval_ext(0:Nbins)
245
246! c -------------------------------------------------------
247! c 0- Initializations
248! c -------------------------------------------------------
249      if ( Nbins .lt. 6) return
250
251      srbval(1) =  S_att
252      srbval(2) =  S_clr
253      srbval(3) =  3.0
254      srbval(4) =  5.0
255      srbval(5) =  7.0
256      srbval(6) = 10.0
257      do i = 7, MIN(10,Nbins)
258       srbval(i) = srbval(i-1) + 5.0
259      enddo
260      DO i = 11, MIN(13,Nbins)
261       srbval(i) = srbval(i-1) + 10.0
262      enddo
263      srbval(MIN(14,Nbins)) = 80.0
264      srbval(Nbins) = xmax
265      cfad(:,:,:) = 0.0
266
267      srbval_ext(1:Nbins) = srbval
268      srbval_ext(0) = -1.0
269! c -------------------------------------------------------
270! c c- Compute CFAD
271! c -------------------------------------------------------
272      do j = 1, Nlevels
273         do ib = 1, Nbins
274            do k = 1, Ncolumns
275               do i = 1, Npoints
276                  if (x(i,k,j) /= undef) then
277                     if ((x(i,k,j).gt.srbval_ext(ib-1)).and.(x(i,k,j).le.srbval_ext(ib))) &
278                          cfad(i,ib,j) = cfad(i,ib,j) + 1.0
279                  else
280                     cfad(i,ib,j) = undef
281                  endif
282               enddo
283            enddo
284         enddo
285      enddo
286
287      where (cfad .ne. undef)  cfad = cfad / float(Ncolumns)
288
289! c -------------------------------------------------------
290      RETURN
291      END SUBROUTINE COSP_CFAD_SR
292
293
294!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
295!-------------------- SUBROUTINE COSP_CLDFRAC -------------------
296! c Purpose: Cloud fraction diagnosed from lidar measurements
297!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
298      SUBROUTINE COSP_CLDFRAC(Npoints,Ncolumns,Nlevels,Ncat,Nphase, &
299                  tmp,x,ATB,ATBperp,pplay,S_att,S_cld,S_cld_att,undef,lidarcld, &
300                  cldlayer,lidarcldphase,nsub,cldlayerphase,lidarcldtemp)
301
302
303      IMPLICIT NONE
304! Input arguments
305      integer Npoints,Ncolumns,Nlevels,Ncat
306      real x(Npoints,Ncolumns,Nlevels)
307
308
309! Local parameters
310      integer nphase ! nb of cloud layer phase types
311                                      ! (ice,liquid,undefined,false ice,false liquid,Percent of ice)
312      integer,parameter  ::  Ntemp=40 ! indice of the temperature vector
313      integer ip, k, iz, ic, ncol, nlev, i, itemp  ! loop indice
314      real  S_cld_att ! New threshold for undefine cloud phase detection (SR=30)       
315      integer toplvlsat  ! level of the first cloud with SR>30
316      real alpha50, beta50, gamma50, delta50, epsilon50, zeta50 ! Polynomial Coef of the phase
317                                                                ! discrimination line   
318
319! Input variables
320      real tmp(Npoints,Nlevels)                 ! temperature
321      real ATB(Npoints,Ncolumns,Nlevels) ! 3D Attenuated backscatter
322      real ATBperp(Npoints,Ncolumns,Nlevels) ! 3D perpendicular attenuated backscatter
323      real pplay(Npoints,Nlevels)
324      real S_att,S_cld
325      real undef
326
327! Output variables
328      real lidarcldtemp(Npoints,Ntemp,5) ! 3D Temperature 1=tot,2=ice,3=liq,4=undef,5=ice/ice+liq
329      real tempmod(Ntemp+1)     ! temperature bins
330      real lidarcldphase(Npoints,Nlevels,Nphase)    ! 3D cloud phase fraction
331      real cldlayerphase(Npoints,Ncat,Nphase) ! low, middle, high, total cloud fractions for ice liquid and undefine phase
332      real lidarcld(Npoints,Nlevels) ! 3D cloud fraction
333      real cldlayer(Npoints,Ncat)    ! low, middle, high, total cloud fractions
334
335! Local variables
336      real tmpi(Npoints,Ncolumns,Nlevels)       ! temperature of ice cld
337      real tmpl(Npoints,Ncolumns,Nlevels)       ! temperature of liquid cld
338      real tmpu(Npoints,Ncolumns,Nlevels)       ! temperature of undef cld
339
340      real checktemp, ATBperp_tmp ! temporary variable
341      real checkcldlayerphase, checkcldlayerphase2 ! temporary variable
342      real sumlidarcldtemp(Npoints,Ntemp) ! temporary variable
343
344      real cldlayphase(Npoints,Ncolumns,Ncat,Nphase) ! subgrided low mid high phase cloud fraction
345      real cldlayerphasetmp(Npoints,Ncat) ! temporary variable
346      real cldlayerphasesum(Npoints,Ncat) ! temporary variable
347      real lidarcldtempind(Npoints,Ntemp) ! 3D Temperature indice
348      real lidarcldphasetmp(Npoints,Nlevels)  ! 3D sum of ice and liquid cloud occurences
349
350
351! Local variables
352      real p1
353      real cldy(Npoints,Ncolumns,Nlevels)
354      real srok(Npoints,Ncolumns,Nlevels)
355      real cldlay(Npoints,Ncolumns,Ncat)
356      real nsublay(Npoints,Ncolumns,Ncat), nsublayer(Npoints,Ncat)
357      real nsub(Npoints,Nlevels)
358
359#ifdef SYS_SX
360      real cldlay1(Npoints,Ncolumns)
361      real cldlay2(Npoints,Ncolumns)
362      real cldlay3(Npoints,Ncolumns)
363      real nsublay1(Npoints,Ncolumns)
364      real nsublay2(Npoints,Ncolumns)
365      real nsublay3(Npoints,Ncolumns)
366#endif
367
368
369
370
371! ---------------------------------------------------------------
372! 1- initialization
373! ---------------------------------------------------------------
374
375      if ( Ncat .ne. 4 ) then
376         print *,'Error in lmd_ipsl_stats.cosp_cldfrac, Ncat must be 4, not',Ncat
377         stop
378      endif
379
380      lidarcld = 0.0
381      nsub = 0.0
382      cldlay = 0.0
383      nsublay = 0.0
384
385      ATBperp_tmp = 0.
386      lidarcldphase(:,:,:) = 0.
387      cldlayphase(:,:,:,:) = 0.
388      cldlayerphase(:,:,:) = 0.
389      tmpi(:,:,:) = 0.
390      tmpl(:,:,:) = 0.
391      tmpu(:,:,:) = 0.
392      cldlayerphasesum(:,:) = 0.
393      lidarcldtemp(:,:,:) = 0.
394      lidarcldtempind(:,:) = 0.
395      sumlidarcldtemp(:,:) = 0.
396      toplvlsat=0
397      lidarcldphasetmp(:,:) = 0.
398
399! temperature bins
400      tempmod=(/-273.15,-90.,-87.,-84.,-81.,-78.,-75.,-72.,-69.,-66.,-63.,-60.,-57., &
401                -54.,-51.,-48.,-45.,-42.,-39.,-36.,-33.,-30.,-27.,-24.,-21.,-18.,  &
402                -15.,-12.,-9.,-6.,-3.,0.,3.,6.,9.,12.,15.,18.,21.,24.,200. /)
403       
404! convert C to K
405      tempmod=tempmod+273.15
406
407! Polynomial coefficient of the phase discrimination line used to separate liquid from ice
408! (Cesana and Chepfer, JGR, 2013)
409! ATBperp = ATB^5*alpha50 + ATB^4*beta50 + ATB^3*gamma50 + ATB^2*delta50 + ATB*epsilon50 + zeta50
410      alpha50   = 9.0322e+15
411      beta50    = -2.1358e+12
412      gamma50   = 173.3963e06
413      delta50   = -3.9514e03
414      epsilon50 = 0.2559
415      zeta50    = -9.4776e-07
416
417
418! ---------------------------------------------------------------
419! 2- Cloud detection
420! ---------------------------------------------------------------
421
422      do k = 1, Nlevels
423
424! cloud detection at subgrid-scale:
425         where ( (x(:,:,k).gt.S_cld) .and. (x(:,:,k).ne. undef) )
426           cldy(:,:,k)=1.0
427         elsewhere
428           cldy(:,:,k)=0.0
429         endwhere
430
431! number of usefull sub-columns:
432         where ( (x(:,:,k).gt.S_att) .and. (x(:,:,k).ne. undef)  )
433           srok(:,:,k)=1.0
434         elsewhere
435           srok(:,:,k)=0.0
436         endwhere
437
438      enddo ! k
439
440
441! ---------------------------------------------------------------
442! 3- grid-box 3D cloud fraction and layered cloud fractions (ISCCP pressure
443! categories) :
444! ---------------------------------------------------------------
445      lidarcld = 0.0
446      nsub = 0.0
447#ifdef SYS_SX
448!! XXX: Use cldlay[1-3] and nsublay[1-3] to avoid bank-conflicts.
449      cldlay1 = 0.0
450      cldlay2 = 0.0
451      cldlay3 = 0.0
452      cldlay(:,:,4) = 0.0 !! XXX: Ncat == 4
453      nsublay1 = 0.0
454      nsublay2 = 0.0
455      nsublay3 = 0.0
456      nsublay(:,:,4) = 0.0
457
458      do k = Nlevels, 1, -1
459       do ic = 1, Ncolumns
460        do ip = 1, Npoints
461
462         if(srok(ip,ic,k).gt.0.)then
463           ! Computation of the cloud fraction as a function of the temperature
464           ! instead of height, for ice,liquid and all clouds
465           do itemp=1,Ntemp
466             if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
467               lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1.
468             endif
469           enddo
470         endif
471
472         if (cldy(ip,ic,k).eq.1.) then
473           do itemp=1,Ntemp
474             if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
475               lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1.
476             endif
477           enddo
478         endif
479
480         p1 = pplay(ip,k)
481
482         if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
483            cldlay3(ip,ic) = MAX(cldlay3(ip,ic), cldy(ip,ic,k))
484            nsublay3(ip,ic) = MAX(nsublay3(ip,ic), srok(ip,ic,k))
485         else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid clouds
486            cldlay2(ip,ic) = MAX(cldlay2(ip,ic), cldy(ip,ic,k))
487            nsublay2(ip,ic) = MAX(nsublay2(ip,ic), srok(ip,ic,k))
488         else
489            cldlay1(ip,ic) = MAX(cldlay1(ip,ic), cldy(ip,ic,k))
490            nsublay1(ip,ic) = MAX(nsublay1(ip,ic), srok(ip,ic,k))
491         endif
492
493         cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4), cldy(ip,ic,k))
494         lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k)
495         nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
496         nsub(ip,k)=nsub(ip,k) + srok(ip,ic,k)
497        enddo
498       enddo
499      enddo
500      cldlay(:,:,1) = cldlay1
501      cldlay(:,:,2) = cldlay2
502      cldlay(:,:,3) = cldlay3
503      nsublay(:,:,1) = nsublay1
504      nsublay(:,:,2) = nsublay2
505      nsublay(:,:,3) = nsublay3
506#else
507      cldlay = 0.0
508      nsublay = 0.0
509      do k = Nlevels, 1, -1
510       do ic = 1, Ncolumns
511        do ip = 1, Npoints
512
513          ! Computation of the cloud fraction as a function of the temperature
514          ! instead of height, for ice,liquid and all clouds
515          if(srok(ip,ic,k).gt.0.)then
516          do itemp=1,Ntemp
517            if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
518              lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1.
519            endif
520          enddo
521          endif
522
523          if(cldy(ip,ic,k).eq.1.)then
524          do itemp=1,Ntemp
525            if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
526              lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1.
527            endif
528          enddo
529          endif
530!
531
532          iz=1
533          p1 = pplay(ip,k)
534          if ( p1.gt.0. .and. p1.lt.(440.*100.)) then ! high clouds
535            iz=3
536          else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid clouds
537            iz=2
538         endif
539
540         cldlay(ip,ic,iz) = MAX(cldlay(ip,ic,iz),cldy(ip,ic,k))
541         cldlay(ip,ic,4) = MAX(cldlay(ip,ic,4),cldy(ip,ic,k))
542         lidarcld(ip,k)=lidarcld(ip,k) + cldy(ip,ic,k)
543
544         nsublay(ip,ic,iz) = MAX(nsublay(ip,ic,iz),srok(ip,ic,k))
545         nsublay(ip,ic,4) = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
546         nsub(ip,k)=nsub(ip,k) + srok(ip,ic,k)
547
548        enddo
549       enddo
550      enddo
551#endif
552
553
554! -- grid-box 3D cloud fraction
555
556      where ( nsub(:,:).gt.0.0 )
557         lidarcld(:,:) = lidarcld(:,:)/nsub(:,:)
558      elsewhere
559         lidarcld(:,:) = undef
560      endwhere
561
562! -- layered cloud fractions
563
564      cldlayer = 0.0
565      nsublayer = 0.0
566
567      do iz = 1, Ncat
568       do ic = 1, Ncolumns
569
570          cldlayer(:,iz)=cldlayer(:,iz) + cldlay(:,ic,iz)
571          nsublayer(:,iz)=nsublayer(:,iz) + nsublay(:,ic,iz)
572
573       enddo
574      enddo
575      where ( nsublayer(:,:).gt.0.0 )
576         cldlayer(:,:) = cldlayer(:,:)/nsublayer(:,:)
577      elsewhere
578         cldlayer(:,:) = undef
579      endwhere
580
581! ---------------------------------------------------------------
582! 4- grid-box 3D cloud Phase :
583! ---------------------------------------------------------------
584! ---------------------------------------------------------------
585! 4.1 - For Cloudy pixels with 8.16km < z < 19.2km
586! ---------------------------------------------------------------
587do ncol=1,Ncolumns
588do i=1,Npoints
589
590      do nlev=Nlevels,18,-1  ! from 19.2km until 8.16km
591         p1 = pplay(i,nlev)
592
593
594! Avoid zero values
595        if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
596! Computation of the ATBperp along the phase discrimination line
597           ATBperp_tmp = (ATB(i,ncol,nlev)**5)*alpha50 + (ATB(i,ncol,nlev)**4)*beta50 + &
598                         (ATB(i,ncol,nlev)**3)*gamma50 + (ATB(i,ncol,nlev)**2)*delta50 + &
599                          ATB(i,ncol,nlev)*epsilon50 + zeta50
600
601!____________________________________________________________________________________________________
602!
603!4.1.a Ice: ATBperp above the phase discrimination line
604!____________________________________________________________________________________________________
605!
606           if( (ATBperp(i,ncol,nlev)-ATBperp_tmp).ge.0. )then   ! Ice clouds
607             ! ICE with temperature above 273,15°K = Liquid (false ice)
608            if(tmp(i,nlev).gt.273.15)then                ! Temperature above 273,15 K
609              ! Liquid: False ice corrected by the temperature to Liquid
610               lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.   ! false ice detection ==> added to Liquid
611               tmpl(i,ncol,nlev)=tmp(i,nlev)
612               lidarcldphase(i,nlev,5)=lidarcldphase(i,nlev,5)+1.   ! keep the information "temperature criterium used"
613                                                    ! to classify the phase cloud
614                   cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
615                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
616                   cldlayphase(i,ncol,3,2) = 1.
617                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
618                   cldlayphase(i,ncol,2,2) = 1.
619                else                                                    ! low cloud
620                   cldlayphase(i,ncol,1,2) = 1.
621                endif
622                   cldlayphase(i,ncol,4,5) = 1.                         ! tot cloud
623                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
624                   cldlayphase(i,ncol,3,5) = 1.
625                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
626                   cldlayphase(i,ncol,2,5) = 1.
627                else                                                    ! low cloud
628                   cldlayphase(i,ncol,1,5) = 1.
629                endif
630
631             else
632             ! ICE with temperature below 273,15°K
633              lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1.
634              tmpi(i,ncol,nlev)=tmp(i,nlev)
635                   cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
636                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
637                   cldlayphase(i,ncol,3,1) = 1.
638                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
639                   cldlayphase(i,ncol,2,1) = 1.
640                else                                                    ! low cloud
641                   cldlayphase(i,ncol,1,1) = 1.
642                endif
643
644              endif
645
646!____________________________________________________________________________________________________
647!
648! 4.1.b Liquid: ATBperp below the phase discrimination line
649!____________________________________________________________________________________________________
650!
651             else                                        ! Liquid clouds
652              ! Liquid with temperature above 231,15°K
653            if(tmp(i,nlev).gt.231.15)then
654               lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.
655               tmpl(i,ncol,nlev)=tmp(i,nlev)
656                   cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
657                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
658                   cldlayphase(i,ncol,3,2) = 1. 
659                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
660                   cldlayphase(i,ncol,2,2) = 1.
661                else                                                    ! low cloud
662                   cldlayphase(i,ncol,1,2) = 1.
663                endif
664
665             else
666             ! Liquid with temperature below 231,15°K = Ice (false liquid)
667               tmpi(i,ncol,nlev)=tmp(i,nlev)
668               lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1.   ! false liquid detection ==> added to ice
669               lidarcldphase(i,nlev,4)=lidarcldphase(i,nlev,4)+1.   ! keep the information "temperature criterium used"
670                                                    ! to classify the phase cloud
671                   cldlayphase(i,ncol,4,4) = 1.                         ! tot cloud
672                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
673                   cldlayphase(i,ncol,3,4) = 1. 
674                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
675                   cldlayphase(i,ncol,2,4) = 1.
676                else                                                    ! low cloud
677                   cldlayphase(i,ncol,1,4) = 1.
678                endif
679                   cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
680                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
681                   cldlayphase(i,ncol,3,1) = 1. 
682                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
683                   cldlayphase(i,ncol,2,1) = 1.
684                else                                                    ! low cloud
685                   cldlayphase(i,ncol,1,1) = 1.
686                endif
687
688             endif
689
690            endif  ! end of discrimination condition
691         endif  ! end of cloud condition
692      enddo ! end of altitude loop
693
694
695
696! ---------------------------------------------------------------
697! 4.2 - For Cloudy pixels with 0km < z < 8.16km
698! ---------------------------------------------------------------
699
700      toplvlsat=0
701      do nlev=17,1,-1  ! from 8.16km until 0km
702         p1 = pplay(i,nlev)
703
704        if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
705! Phase discrimination line : ATBperp = ATB^5*alpha50 + ATB^4*beta50 + ATB^3*gamma50 + ATB^2*delta50
706!                                  + ATB*epsilon50 + zeta50
707! Computation of the ATBperp of the phase discrimination line
708           ATBperp_tmp = (ATB(i,ncol,nlev)**5)*alpha50 + (ATB(i,ncol,nlev)**4)*beta50 + &
709                         (ATB(i,ncol,nlev)**3)*gamma50 + (ATB(i,ncol,nlev)**2)*delta50 + &
710                          ATB(i,ncol,nlev)*epsilon50 + zeta50
711!____________________________________________________________________________________________________
712!
713! 4.2.a Ice: ATBperp above the phase discrimination line
714!____________________________________________________________________________________________________
715!
716            ! ICE with temperature above 273,15°K = Liquid (false ice)
717          if( (ATBperp(i,ncol,nlev)-ATBperp_tmp).ge.0. )then   ! Ice clouds
718            if(tmp(i,nlev).gt.273.15)then
719               lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.  ! false ice ==> liq
720               tmpl(i,ncol,nlev)=tmp(i,nlev)
721               lidarcldphase(i,nlev,5)=lidarcldphase(i,nlev,5)+1.
722
723                   cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
724               if ( p1.gt.0. .and. p1.lt.(440.*100.)) then              ! high cloud
725                   cldlayphase(i,ncol,3,2) = 1.
726                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
727                   cldlayphase(i,ncol,2,2) = 1.
728                else                                                    ! low cloud
729                   cldlayphase(i,ncol,1,2) = 1.
730                endif
731
732                   cldlayphase(i,ncol,4,5) = 1.                         ! tot cloud
733                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
734                   cldlayphase(i,ncol,3,5) = 1.
735                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
736                   cldlayphase(i,ncol,2,5) = 1.
737                else                                                    ! low cloud
738                   cldlayphase(i,ncol,1,5) = 1.
739                endif
740
741             else
742              ! ICE with temperature below 273,15°K
743              lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1.
744              tmpi(i,ncol,nlev)=tmp(i,nlev)
745
746                   cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
747                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
748                   cldlayphase(i,ncol,3,1) = 1.
749                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
750                   cldlayphase(i,ncol,2,1) = 1.
751                else                                                    ! low cloud
752                   cldlayphase(i,ncol,1,1) = 1.
753                endif
754
755              endif
756
757!____________________________________________________________________________________________________
758!
759! 4.2.b Liquid: ATBperp below the phase discrimination line
760!____________________________________________________________________________________________________
761!
762          else 
763             ! Liquid with temperature above 231,15°K
764            if(tmp(i,nlev).gt.231.15)then
765               lidarcldphase(i,nlev,2)=lidarcldphase(i,nlev,2)+1.
766               tmpl(i,ncol,nlev)=tmp(i,nlev)
767
768                   cldlayphase(i,ncol,4,2) = 1.                         ! tot cloud
769                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
770                   cldlayphase(i,ncol,3,2) = 1. 
771                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
772                   cldlayphase(i,ncol,2,2) = 1.
773                else                                                    ! low cloud
774                   cldlayphase(i,ncol,1,2) = 1.
775                endif
776
777             else
778             ! Liquid with temperature below 231,15°K = Ice (false liquid)
779               tmpi(i,ncol,nlev)=tmp(i,nlev)
780               lidarcldphase(i,nlev,1)=lidarcldphase(i,nlev,1)+1.  ! false liq ==> ice
781               lidarcldphase(i,nlev,4)=lidarcldphase(i,nlev,4)+1.  ! false liq ==> ice
782
783                   cldlayphase(i,ncol,4,4) = 1.                         ! tot cloud
784                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
785                   cldlayphase(i,ncol,3,4) = 1. 
786                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
787                   cldlayphase(i,ncol,2,4) = 1.
788                else                                                    ! low cloud
789                   cldlayphase(i,ncol,1,4) = 1.
790                endif
791
792                   cldlayphase(i,ncol,4,1) = 1.                         ! tot cloud
793                if ( p1.gt.0. .and. p1.lt.(440.*100.)) then             ! high cloud
794                   cldlayphase(i,ncol,3,1) = 1. 
795                else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then ! mid cloud
796                   cldlayphase(i,ncol,2,1) = 1.
797                else                                                    ! low cloud
798                   cldlayphase(i,ncol,1,1) = 1.
799                endif
800
801             endif
802           endif  ! end of discrimination condition
803
804            toplvlsat=0
805
806           ! Find the level of the highest cloud with SR>30
807            if(x(i,ncol,nlev).gt.S_cld_att)then  ! SR > 30.
808                toplvlsat=nlev-1
809                goto 99
810            endif
811
812        endif  ! end of cloud condition
813       enddo  ! end of altitude loop
814
81599 continue
816
817!____________________________________________________________________________________________________
818!
819! Undefined phase: For a cloud located below another cloud with SR>30
820! see Cesana and Chepfer 2013 Sect.III.2
821!____________________________________________________________________________________________________
822!
823if(toplvlsat.ne.0)then         
824      do nlev=toplvlsat,1,-1
825         p1 = pplay(i,nlev)
826        if(cldy(i,ncol,nlev).eq.1.)then
827           lidarcldphase(i,nlev,3)=lidarcldphase(i,nlev,3)+1.
828           tmpu(i,ncol,nlev)=tmp(i,nlev)
829
830                   cldlayphase(i,ncol,4,3) = 1.                         ! tot cloud
831          if ( p1.gt.0. .and. p1.lt.(440.*100.)) then              ! high cloud
832             cldlayphase(i,ncol,3,3) = 1.
833          else if(p1.ge.(440.*100.) .and. p1.lt.(680.*100.)) then  ! mid cloud
834             cldlayphase(i,ncol,2,3) = 1.
835          else                                                     ! low cloud
836             cldlayphase(i,ncol,1,3) = 1.
837          endif
838
839        endif   
840      enddo
841endif
842     
843      toplvlsat=0
844
845enddo
846enddo
847
848
849
850!____________________________________________________________________________________________________
851!
852! Computation of final cloud phase diagnosis
853!____________________________________________________________________________________________________
854!
855
856! Compute the Ice percentage in cloud = ice/(ice+liq) as a function
857! of the occurrences
858lidarcldphasetmp(:,:)=lidarcldphase(:,:,1)+lidarcldphase(:,:,2);
859WHERE (lidarcldphasetmp(:,:).gt. 0.)
860   lidarcldphase(:,:,6)=lidarcldphase(:,:,1)/lidarcldphasetmp(:,:)
861ELSEWHERE
862   lidarcldphase(:,:,6) = undef
863ENDWHERE
864
865! Compute Phase 3D Cloud Fraction
866     WHERE ( nsub(:,:).gt.0.0 )
867       lidarcldphase(:,:,1)=lidarcldphase(:,:,1)/nsub(:,:)
868       lidarcldphase(:,:,2)=lidarcldphase(:,:,2)/nsub(:,:)
869       lidarcldphase(:,:,3)=lidarcldphase(:,:,3)/nsub(:,:)
870       lidarcldphase(:,:,4)=lidarcldphase(:,:,4)/nsub(:,:)
871       lidarcldphase(:,:,5)=lidarcldphase(:,:,5)/nsub(:,:)
872     ELSEWHERE
873       lidarcldphase(:,:,1) = undef
874       lidarcldphase(:,:,2) = undef
875       lidarcldphase(:,:,3) = undef
876       lidarcldphase(:,:,4) = undef
877       lidarcldphase(:,:,5) = undef
878     ENDWHERE
879
880
881! Compute Phase low mid high cloud fractions
882    do iz = 1, Ncat
883       do i=1,Nphase-3
884       do ic = 1, Ncolumns
885          cldlayerphase(:,iz,i)=cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i)
886          cldlayerphasesum(:,iz)=cldlayerphasesum(:,iz)+cldlayphase(:,ic,iz,i)
887       enddo
888      enddo
889    enddo
890
891    do iz = 1, Ncat
892       do i=4,5
893       do ic = 1, Ncolumns
894          cldlayerphase(:,iz,i)=cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i)         
895       enddo
896       enddo
897    enddo
898   
899! Compute the Ice percentage in cloud = ice/(ice+liq)
900cldlayerphasetmp(:,:)=cldlayerphase(:,:,1)+cldlayerphase(:,:,2)
901    WHERE (cldlayerphasetmp(:,:).gt. 0.)
902       cldlayerphase(:,:,6)=cldlayerphase(:,:,1)/cldlayerphasetmp(:,:)
903    ELSEWHERE
904       cldlayerphase(:,:,6) = undef
905    ENDWHERE
906
907    do i=1,Nphase-1
908      WHERE ( cldlayerphasesum(:,:).gt.0.0 )
909         cldlayerphase(:,:,i) = (cldlayerphase(:,:,i)/cldlayerphasesum(:,:)) * cldlayer(:,:)
910      ENDWHERE
911    enddo
912
913
914    do i=1,Npoints
915       do iz=1,Ncat
916          checkcldlayerphase=0.
917          checkcldlayerphase2=0.
918
919          if (cldlayerphasesum(i,iz).gt.0.0 )then
920             do ic=1,Nphase-3
921                checkcldlayerphase=checkcldlayerphase+cldlayerphase(i,iz,ic) 
922             enddo
923             checkcldlayerphase2=cldlayer(i,iz)-checkcldlayerphase
924             if( (checkcldlayerphase2.gt.0.01).or.(checkcldlayerphase2.lt.-0.01) ) print *, checkcldlayerphase,cldlayer(i,iz)
925
926          endif
927
928       enddo
929    enddo
930
931    do i=1,Nphase-1
932      WHERE ( nsublayer(:,:).eq.0.0 )
933         cldlayerphase(:,:,i) = undef
934      ENDWHERE
935   enddo
936
937
938
939! Compute Phase 3D as a function of temperature
940do nlev=1,Nlevels
941do ncol=1,Ncolumns     
942do i=1,Npoints
943do itemp=1,Ntemp
944if(tmpi(i,ncol,nlev).gt.0.)then
945      if( (tmpi(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpi(i,ncol,nlev).lt.tempmod(itemp+1)) )then
946        lidarcldtemp(i,itemp,2)=lidarcldtemp(i,itemp,2)+1.
947      endif
948elseif(tmpl(i,ncol,nlev).gt.0.)then
949      if( (tmpl(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpl(i,ncol,nlev).lt.tempmod(itemp+1)) )then
950        lidarcldtemp(i,itemp,3)=lidarcldtemp(i,itemp,3)+1.
951      endif
952elseif(tmpu(i,ncol,nlev).gt.0.)then
953      if( (tmpu(i,ncol,nlev).ge.tempmod(itemp)).and.(tmpu(i,ncol,nlev).lt.tempmod(itemp+1)) )then
954        lidarcldtemp(i,itemp,4)=lidarcldtemp(i,itemp,4)+1.
955      endif
956endif
957enddo
958enddo
959enddo
960enddo
961
962! Check temperature cloud fraction
963do i=1,Npoints
964   do itemp=1,Ntemp
965checktemp=lidarcldtemp(i,itemp,2)+lidarcldtemp(i,itemp,3)+lidarcldtemp(i,itemp,4)
966
967        if(checktemp.NE.lidarcldtemp(i,itemp,1))then
968          print *, i,itemp
969          print *, lidarcldtemp(i,itemp,1:4)
970        endif
971
972   enddo
973enddo
974
975! Compute the Ice percentage in cloud = ice/(ice+liq)
976!   sumlidarcldtemp=sum(lidarcldtemp(:,:,2:3),3)
977   sumlidarcldtemp(:,:)=lidarcldtemp(:,:,2)+lidarcldtemp(:,:,3)
978
979WHERE(sumlidarcldtemp(:,:)>0.)
980  lidarcldtemp(:,:,5)=lidarcldtemp(:,:,2)/sumlidarcldtemp(:,:)
981ELSEWHERE
982  lidarcldtemp(:,:,5)=undef
983ENDWHERE
984
985do i=1,4
986  WHERE(lidarcldtempind(:,:).gt.0.)
987     lidarcldtemp(:,:,i) = lidarcldtemp(:,:,i)/lidarcldtempind(:,:)
988  ELSEWHERE
989     lidarcldtemp(:,:,i) = undef
990  ENDWHERE
991enddo
992
993       RETURN
994      END SUBROUTINE COSP_CLDFRAC
995! ---------------------------------------------------------------
996
997! BEGINNING OF OPAQ CHANGES
998    ! ####################################################################################
999    ! SUBROUTINE cosp_opaq
1000    ! Conventions: Ntype must be equal to 3 (opaque cloud, thin cloud, z_opaque)
1001    ! ####################################################################################
1002    SUBROUTINE COSP_OPAQ(Npoints,Ncolumns,Nlevels,Ntype,x,S_cld,undef,lidarcldtype,   &
1003                         cldtype,vgrid_z)
1004
1005      IMPLICIT NONE
1006! Input arguments
1007      integer Npoints,Ncolumns,Nlevels,Ntype
1008      real x(Npoints,Ncolumns,Nlevels)
1009      real S_cld
1010      real undef
1011      real vgrid_z(Nlevels)
1012! Output :
1013      real lidarcldtype(Npoints,Nlevels,Ntype+1) ! 3D "lidar" OPAQ type + opacity fraction
1014      real cldtype(Npoints,Ntype)              ! opaque and thin cloud covers, z_opaque
1015! Local variables
1016      integer ip, k, iz, ic, zopac
1017      real p1
1018      real cldy(Npoints,Ncolumns,Nlevels)
1019      real cldyopaq(Npoints,Ncolumns,Nlevels)
1020      real srok(Npoints,Ncolumns,Nlevels)
1021      real srokopaq(Npoints,Ncolumns,Nlevels)
1022      real cldlay(Npoints,Ncolumns,Ntype+1)  ! opaque, thin, z_opaque and all cloud cover
1023      real nsublay(Npoints,Ncolumns,Ntype+1) ! opaque, thin, z_opaque and all cloud cover
1024      real nsublayer(Npoints,Ntype)
1025      real nsub(Npoints,Nlevels)
1026      real nsubopaq(Npoints,Nlevels)
1027      real S_att_opaq
1028      real S_att
1029 
1030    ! ####################################################################################
1031        ! 1) Initialize   
1032    ! ####################################################################################
1033    cldtype               = 0.0
1034    lidarcldtype          = 0.0
1035    nsub                  = 0.0
1036    nsubopaq              = 0.0
1037    cldlay                = 0.0
1038    nsublay               = 0.0
1039    nsublayer             = 0.0
1040    S_att_opaq            = 0.06 ! Fully Attenuated threshold, from Guzman et al. 2017, JGR-A
1041    S_att                 = 0.01
1042
1043    ! ####################################################################################
1044    ! 2) Cloud detection and Fully attenuated layer detection
1045    ! ####################################################################################
1046    do k=1,Nlevels
1047       ! Cloud detection at subgrid-scale:
1048       where ( (x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
1049          cldy(:,:,k)=1.0
1050       elsewhere
1051          cldy(:,:,k)=0.0
1052       endwhere
1053       ! Fully attenuated layer detection at subgrid-scale:
1054       where ( (x(:,:,k) .gt. 0.0) .and. (x(:,:,k) .lt. S_att_opaq) .and. (x(:,:,k) .ne. undef) )
1055          cldyopaq(:,:,k)=1.0
1056       elsewhere
1057          cldyopaq(:,:,k)=0.0
1058       endwhere
1059
1060       ! Number of useful sub-column layers:
1061       where ( (x(:,:,k) .gt. S_att) .and. (x(:,:,k) .ne. undef) )
1062          srok(:,:,k)=1.0
1063       elsewhere
1064          srok(:,:,k)=0.0
1065       endwhere
1066       ! Number of useful sub-columns layers for z_opaque 3D fraction:
1067       where ( (x(:,:,k) .gt. 0.0) .and. (x(:,:,k) .ne. undef) )
1068          srokopaq(:,:,k)=1.0
1069       elsewhere
1070          srokopaq(:,:,k)=0.0
1071       endwhere
1072    enddo
1073
1074    ! ####################################################################################
1075    ! 3) Grid-box 3D OPAQ product fraction and cloud type cover (opaque/thin) + mean z_opaque
1076    ! ####################################################################################
1077
1078    do k= Nlevels,1,-1
1079       do ic = 1, Ncolumns
1080          do ip = 1, Npoints
1081
1082             cldlay(ip,ic,1)   = MAX(cldlay(ip,ic,1),cldyopaq(ip,ic,k)) ! Opaque clouds
1083             cldlay(ip,ic,4)   = MAX(cldlay(ip,ic,4),cldy(ip,ic,k))     ! All clouds
1084
1085             nsublay(ip,ic,1)  = MAX(nsublay(ip,ic,1),srok(ip,ic,k))
1086             nsublay(ip,ic,2)  = MAX(nsublay(ip,ic,2),srok(ip,ic,k))
1087!             nsublay(ip,ic,4)  = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
1088             nsub(ip,k)        = nsub(ip,k) + srok(ip,ic,k)
1089             nsubopaq(ip,k)    = nsubopaq(ip,k) + srokopaq(ip,ic,k)
1090
1091          enddo
1092       enddo
1093    enddo   
1094
1095! OPAQ variables
1096     do ic = 1, Ncolumns
1097        do ip = 1, Npoints
1098
1099     ! Declaring non-opaque cloudy profiles as thin cloud profiles
1100           if ( (cldlay(ip,ic,4) .eq. 1.0) .and. (cldlay(ip,ic,1) .eq. 0.0) ) then
1101              cldlay(ip,ic,2)  =  1.0
1102           endif
1103
1104     ! Filling in 3D and 2D variables
1105
1106     ! Opaque cloud profiles
1107           if ( cldlay(ip,ic,1) .eq. 1.0 ) then
1108              zopac = 0.0
1109              do k=2,Nlevels
1110     ! Declaring opaque cloud fraction and z_opaque altitude for 3D and 2D variables
1111                 if ( (cldy(ip,ic,k) .eq. 1.0) .and. (zopac .eq. 0.0) ) then
1112                    lidarcldtype(ip,k-1,3) = lidarcldtype(ip,k-1,3) + 1.0
1113                    cldlay(ip,ic,3)        = vgrid_z(k-1) !z_opaque altitude
1114                    nsublay(ip,ic,3)       = 1.0
1115                    zopac = 1.0
1116                 endif
1117                 if ( cldy(ip,ic,k) .eq. 1.0 ) then
1118                    lidarcldtype(ip,k,1)   = lidarcldtype(ip,k,1) + 1.0
1119                 endif
1120              enddo
1121           endif
1122
1123     ! Thin cloud profiles
1124           if ( cldlay(ip,ic,2) .eq. 1.0 ) then
1125              do k=1,Nlevels
1126     ! Declaring thin cloud fraction for 3D variable
1127                 if ( cldy(ip,ic,k) .eq. 1.0 ) then
1128                    lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1.0
1129                 endif
1130              enddo
1131           endif
1132
1133       enddo
1134    enddo   
1135
1136    ! 3D cloud types fraction (opaque=1 and thin=2)
1137    where ( nsub(:,:) .gt. 0.0 )
1138       lidarcldtype(:,:,1) = lidarcldtype(:,:,1)/nsub(:,:)
1139       lidarcldtype(:,:,2) = lidarcldtype(:,:,2)/nsub(:,:)
1140    elsewhere
1141       lidarcldtype(:,:,1) = undef
1142       lidarcldtype(:,:,2) = undef
1143    endwhere
1144    ! 3D z_opaque fraction (=3)
1145    where ( nsubopaq(:,:) .gt. 0.0 )
1146       lidarcldtype(:,:,3) = lidarcldtype(:,:,3)/nsubopaq(:,:)
1147    elsewhere
1148       lidarcldtype(:,:,3) = undef
1149    endwhere
1150    ! 3D opacity fraction (=4) !Summing z_opaque fraction from TOA(k=Nlevels) to SFC(k=1)
1151       lidarcldtype(:,Nlevels,4) = lidarcldtype(:,Nlevels,3)
1152    do ip = 1, Npoints
1153        do k = Nlevels-1, 1, -1
1154           if ( lidarcldtype(ip,k,3) .ne. undef ) then
1155              lidarcldtype(ip,k,4) = lidarcldtype(ip,k+1,4) + lidarcldtype(ip,k,3)
1156           endif
1157        enddo
1158    enddo
1159    where ( nsubopaq(:,:) .eq. 0.0 )
1160       lidarcldtype(:,:,4) = undef
1161    endwhere
1162
1163    ! Layered cloud types (opaque, thin and z_opaque 2D variables)
1164
1165    do iz = 1, Ntype
1166       do ic = 1, Ncolumns
1167          cldtype(:,iz)   = cldtype(:,iz)   + cldlay(:,ic,iz)
1168          nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz)
1169       enddo
1170    enddo
1171    where (nsublayer(:,:) .gt. 0.0)
1172       cldtype(:,:) = cldtype(:,:)/nsublayer(:,:)
1173    elsewhere
1174       cldtype(:,:) = undef
1175    endwhere
1176
1177  END SUBROUTINE COSP_OPAQ
1178! END OF OPAQ CHANGES
1179
1180
1181END MODULE MOD_LMD_IPSL_STATS
Note: See TracBrowser for help on using the repository browser.