source: LMDZ6/branches/contrails/libf/phylmd/cospv2/lidar_simulator.f90 @ 5440

Last change on this file since 5440 was 5268, checked in by abarral, 2 months ago

.f90 <-> .F90 depending on cpp key use
File size: 68.8 KB
Line 
1! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2! Copyright (c) 2009, Centre National de la Recherche Scientifique
3! All rights reserved.
4!
5! Redistribution and use in source and binary forms, with or without modification, are
6! permitted provided that the following conditions are met:
7!
8! 1. Redistributions of source code must retain the above copyright notice, this list of
9!    conditions and the following disclaimer.
10!
11! 2. 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
13!    materials provided with the distribution.
14!
15! 3. Neither the name of the copyright holder nor the names of its contributors may be
16!    used to endorse or promote products derived from this software without specific prior
17!    written permission.
18!
19! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY
20! EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF
21! MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL
22! THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
23! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT
24! OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
25! INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
27! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
28!
29! History
30! May 2007: ActSim code of M. Chiriaco and H. Chepfer rewritten by S. Bony
31!
32! May 2008, H. Chepfer:
33! - Units of pressure inputs: Pa
34! - Non Spherical particles : LS Ice NS coefficients, CONV Ice NS coefficients
35! - New input: ice_type (0=ice-spheres ; 1=ice-non-spherical)
36!
37! June 2008, A. Bodas-Salcedo:
38! - Ported to Fortran 90 and optimisation changes
39!
40! August 2008, J-L Dufresne:
41! - Optimisation changes (sum instructions suppressed)
42!
43! October 2008, S. Bony,  H. Chepfer and J-L. Dufresne : 
44! - Interface with COSP v2.0:
45!      cloud fraction removed from inputs
46!      in-cloud condensed water now in input (instead of grid-averaged value)
47!      depolarisation diagnostic removed
48!      parasol (polder) reflectances (for 5 different solar zenith angles) added
49!
50! December 2008, S. Bony,  H. Chepfer and J-L. Dufresne :
51! - Modification of the integration of the lidar equation.
52! - change the cloud detection threshold
53!
54! April 2008, A. Bodas-Salcedo:
55! - Bug fix in computation of pmol and pnorm of upper layer
56!
57! April 2008, J-L. Dufresne
58! - Bug fix in computation of pmol and pnorm, thanks to Masaki Satoh: a factor 2
59! was missing. This affects the ATB values but not the cloud fraction.
60!
61! January 2013, G. Cesana and H. Chepfer:
62! - Add the perpendicular component of the backscattered signal (pnorm_perp_tot) in the arguments
63! - Add the temperature for each levels (temp) in the arguments
64! - Add the computation of the perpendicular component of the backscattered lidar signal
65! Reference: Cesana G. and H. Chepfer (2013): Evaluation of the cloud water phase
66! in a climate model using CALIPSO-GOCCP, J. Geophys. Res., doi: 10.1002/jgrd.50376
67!
68! May 2015 - D. Swales - Modified for COSPv2.0
69!
70! Mar 2018 - R. Guzman - Added OPAQ subroutines
71! References OPAQ:
72!
73!       Guzman et al. (2017): Direct atmosphere opacity observations from CALIPSO provide
74! new constraints on cloud-radiation interactions. JGR-Atmospheres, DOI: 10.1002/2016JD025946
75!       Vaillant de Guelis et al. (2017a): The link between outgoing longwave radiation and
76! the altitude at which a spaceborne lidar beam is fully attenuated. AMT, 10, 4659-4685,
77! https://doi.org/10.5194/amt-10-4659-2017
78!
79! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
80module mod_lidar_simulator
81  USE COSP_KINDS,         ONLY: wp
82  USE MOD_COSP_CONFIG,    ONLY: SR_BINS,S_CLD,S_ATT,S_CLD_ATT,R_UNDEF,calipso_histBsct,  &
83                                use_vgrid,vgrid_zl,vgrid_zu,vgrid_z,atlid_histBsct,      &
84                                grLidar532_histBsct,S_CLD_ATLID,S_ATT_ATLID,S_CLD_ATT_ATLID
85  USE MOD_COSP_STATS,     ONLY: COSP_CHANGE_VERTICAL_GRID,hist1d
86  implicit none
87 
88  ! Polynomial coefficients (Alpha, Beta, Gamma) which allow to compute the
89  ! ATBperpendicular as a function of the ATB for ice or liquid cloud particles
90  ! derived from CALIPSO-GOCCP observations at 120m vertical grid
91  ! (Cesana and Chepfer, JGR, 2013).
92  !
93  ! Relationship between ATBice and ATBperp,ice for ice particles:
94  !                ATBperp,ice = Alpha*ATBice
95  ! Relationship between ATBice and ATBperp,ice for liquid particles:
96  !          ATBperp,ice = Beta*ATBice^2 + Gamma*ATBice
97  real(wp) :: &
98       alpha,beta,gamma   
99
100contains
101  ! ######################################################################################
102  ! SUBROUTINE lidar_subcolumn
103  ! Inputs with a vertical dimensions (nlev) should ordered in along the vertical
104  ! dimension from TOA-2-SFC, for example: varIN(nlev) is varIN @ SFC.
105  ! ######################################################################################
106  subroutine lidar_subcolumn(npoints, ncolumns, nlev, lground, beta_mol, tau_mol,        &
107       betatot, tautot, pmol, pnorm, betatot_ice, tautot_ice, betatot_liq, tautot_liq,   &
108       pnorm_perp_tot)
109
110    ! INPUTS
111    INTEGER,intent(in) :: &
112         npoints,      & ! Number of gridpoints
113         ncolumns,     & ! Number of subcolumns
114         nlev            ! Number of levels
115    logical,intent(in) :: &
116         lground         ! True for ground-based lidar simulator
117    REAL(WP),intent(in),dimension(npoints,nlev) :: &
118         beta_mol,     & ! Molecular backscatter coefficient
119         tau_mol         ! Molecular optical depth
120    REAL(WP),intent(in),dimension(npoints,ncolumns,nlev)       :: &
121         betatot,      & !
122         tautot          ! Optical thickess integrated from top
123    ! Optional Inputs
124    REAL(WP),intent(in),dimension(npoints,ncolumns,nlev),optional       :: &         
125         betatot_ice,  & ! Backscatter coefficient for ice particles
126         betatot_liq,  & ! Backscatter coefficient for liquid particles
127         tautot_ice,   & ! Total optical thickness of ice
128         tautot_liq      ! Total optical thickness of liq
129
130    ! OUTPUTS
131    REAL(WP),intent(out),dimension(npoints,nlev) :: &
132         pmol            ! Molecular attenuated backscatter lidar signal power(m^-1.sr^-1)
133    REAL(WP),intent(out),dimension(npoints,ncolumns,nlev) :: &
134         pnorm           ! Molecular backscatter signal power (m^-1.sr^-1)
135    ! Optional outputs
136    REAL(WP),intent(out),dimension(npoints,ncolumns,nlev),optional :: &
137         pnorm_perp_tot  ! Perpendicular lidar backscattered signal power
138
139    ! LOCAL VARIABLES
140    INTEGER :: k,icol,zi,zf,zinc
141    logical :: lphaseoptics
142    REAL(WP),dimension(npoints) :: &
143         tautot_lay        !
144    REAL(WP),dimension(npoints,ncolumns,nlev) :: &
145         pnorm_liq,      & ! Lidar backscattered signal power for liquid
146         pnorm_ice,      & ! Lidar backscattered signal power for ice
147         pnorm_perp_ice, & ! Perpendicular lidar backscattered signal power for ice
148         pnorm_perp_liq, & ! Perpendicular lidar backscattered signal power for liq
149         beta_perp_ice,  & ! Perpendicular backscatter coefficient for ice
150         beta_perp_liq     ! Perpendicular backscatter coefficient for liquid
151
152    ! Phase optics?
153    lphaseoptics=.false.
154    if (present(betatot_ice) .and. present(betatot_liq) .and. present(tautot_liq) .and. &
155         present(tautot_ice)) lphaseoptics=.true.
156   
157    ! Is this lidar spaceborne (default) or ground-based?
158    if (lground) then
159       zi   = nlev
160       zf   = 1
161       zinc = -1
162    else
163       zi   = 1
164       zf   = nlev
165       zinc = 1
166    endif   
167
168    ! ####################################################################################
169    ! *) Molecular signal
170    ! ####################################################################################
171    call cmp_backsignal(nlev,npoints,beta_mol(1:npoints,zi:zf:zinc),&
172                        tau_mol(1:npoints,zi:zf:zinc),pmol(1:npoints,zi:zf:zinc))
173                       
174    ! ####################################################################################
175    ! PLANE PARRALLEL FIELDS
176    ! ####################################################################################
177    do icol=1,ncolumns
178       ! #################################################################################
179       ! *) Total Backscatter signal
180       ! #################################################################################
181       call cmp_backsignal(nlev,npoints,betatot(1:npoints,icol,zi:zf:zinc),&
182            tautot(1:npoints,icol,zi:zf:zinc),pnorm(1:npoints,icol,zi:zf:zinc))
183       
184       ! #################################################################################
185       ! *) Ice/Liq Backscatter signal
186       ! #################################################################################
187       if (lphaseoptics) then
188          ! Computation of the ice and liquid lidar backscattered signal (ATBice and ATBliq)
189          ! Ice only
190          call cmp_backsignal(nlev,npoints,betatot_ice(1:npoints,icol,zi:zf:zinc),&
191               tautot_ice(1:npoints,icol,zi:zf:zinc), pnorm_ice(1:npoints,icol,zi:zf:zinc))
192          ! Liquid only
193          call cmp_backsignal(nlev,npoints,betatot_liq(1:npoints,icol,zi:zf:zinc),&
194               tautot_liq(1:npoints,icol,zi:zf:zinc), pnorm_liq(1:npoints,icol,zi:zf:zinc))
195       endif
196    enddo
197
198    ! ####################################################################################
199    ! PERDENDICULAR FIELDS (Only needed if distinguishing by phase (ice/liquid))
200    ! ####################################################################################
201    if (lphaseoptics) then
202       do icol=1,ncolumns
203          ! #################################################################################
204          ! *) Ice/Liq Perpendicular Backscatter signal
205          ! #################################################################################
206          ! Computation of ATBperp,ice/liq from ATBice/liq including the multiple scattering
207          ! contribution (Cesana and Chepfer 2013, JGR)
208          do k=1,nlev
209             ! Ice particles
210             pnorm_perp_ice(1:npoints,icol,k) = Alpha * pnorm_ice(1:npoints,icol,k)
211             
212             ! Liquid particles
213             pnorm_perp_liq(1:npoints,icol,k) = 1000._wp*Beta*pnorm_liq(1:npoints,icol,k)**2+&
214                  Gamma*pnorm_liq(1:npoints,icol,k)
215          enddo
216         
217          ! #################################################################################
218          ! *) Computation of beta_perp_ice/liq using the lidar equation
219          ! #################################################################################
220          ! Ice only
221          call cmp_beta(nlev,npoints,pnorm_perp_ice(1:npoints,icol,zi:zf:zinc),&
222               tautot_ice(1:npoints,icol,zi:zf:zinc),beta_perp_ice(1:npoints,icol,zi:zf:zinc))       
223         
224          ! Liquid only
225          call cmp_beta(nlev,npoints,pnorm_perp_liq(1:npoints,icol,zi:zf:zinc),&
226               tautot_liq(1:npoints,icol,zi:zf:zinc),beta_perp_liq(1:npoints,icol,zi:zf:zinc))
227         
228          ! #################################################################################
229          ! *) Perpendicular Backscatter signal
230          ! #################################################################################
231          ! Computation of the total perpendicular lidar signal (ATBperp for liq+ice)
232          ! Upper layer
233          WHERE(tautot(1:npoints,icol,1) .gt. 0)
234             pnorm_perp_tot(1:npoints,icol,1) = (beta_perp_ice(1:npoints,icol,1)+           &
235                  beta_perp_liq(1:npoints,icol,1)-                                          &
236                  (beta_mol(1:npoints,1)/(1._wp+1._wp/0.0284_wp))) /                        &
237                  (2._wp*tautot(1:npoints,icol,1))*                                         &
238                  (1._wp-exp(-2._wp*tautot(1:npoints,icol,1)))
239          ELSEWHERE
240             pnorm_perp_tot(1:npoints,icol,1) = 0._wp
241          ENDWHERE
242         
243          ! Other layers
244          do k=2,nlev
245             ! Optical thickness of layer k
246             tautot_lay(1:npoints) = tautot(1:npoints,icol,k)-tautot(1:npoints,icol,k-1)
247             
248             ! The perpendicular component of the molecular backscattered signal (Betaperp)
249             ! has been taken into account two times (once for liquid and once for ice).
250             ! We remove one contribution using
251             ! Betaperp=beta_mol(:,k)/(1+1/0.0284)) [bodhaine et al. 1999] in the following
252             ! equations:
253             WHERE (pnorm(1:npoints,icol,k) .eq. 0)
254                pnorm_perp_tot(1:npoints,icol,k)=0._wp
255             ELSEWHERE
256                where(tautot_lay(1:npoints) .gt. 0.)
257                   pnorm_perp_tot(1:npoints,icol,k) = (beta_perp_ice(1:npoints,icol,k)+     &
258                        beta_perp_liq(1:npoints,icol,k)-(beta_mol(1:npoints,k)/(1._wp+1._wp/  &
259                        0.0284_wp)))*EXP(-2._wp*tautot(1:npoints,icol,k-1))/                  &
260                        (2._wp*tautot_lay(1:npoints))* (1._wp-EXP(-2._wp*tautot_lay(1:npoints)))
261                elsewhere
262                   pnorm_perp_tot(1:npoints,icol,k) = (beta_perp_ice(1:npoints,icol,k)+     &
263                        beta_perp_liq(1:npoints,icol,k)-(beta_mol(1:npoints,k)/(1._wp+1._wp/  &
264                        0.0284_wp)))*EXP(-2._wp*tautot(1:npoints,icol,k-1))
265                endwhere
266             ENDWHERE
267          END DO
268       enddo
269    end if
270  end subroutine lidar_subcolumn
271
272  ! ######################################################################################
273  ! SUBROUTINE lidar_column
274  ! ######################################################################################
275  subroutine lidar_column(npoints, ncol, nlevels, llm, max_bin, ntype, platform, pnorm, pmol,             &
276       pplay, zlev, zlev_half, vgrid_z, ok_lidar_cfad, ncat, cfad2, lidarcld, cldlayer,  &
277       ! Optional stuff below
278       tmp, pnorm_perp, surfelev,  lidarcldphase, lidarcldtype, cldtype, cldtypetemp, &
279       cldtypemeanz, cldtypemeanzse, cldthinemis, cldlayerphase, lidarcldtmp)
280
281    integer,parameter :: &
282         nphase = 6 ! Number of cloud layer phase types
283
284    ! Inputs
285    integer,intent(in) :: &
286         npoints, & ! Number of horizontal grid points
287         ncol,    & ! Number of subcolumns
288         nlevels, & ! Number of vertical layers (OLD grid)
289         llm,     & ! Number of vertical layers (NEW grid)
290         max_bin, & ! Number of bins for SR CFADs
291         ncat,    & ! Number of cloud layer types (low,mid,high,total)
292         ntype      ! Number of OPAQ products (opaque/thin cloud + z_opaque)
293    character(len=*),intent(in) :: &
294         platform   ! Name of platform (e.g. calipso,atlid,grLidar532)
295    real(wp),intent(in),dimension(npoints,ncol,Nlevels) :: &
296         pnorm      ! Lidar ATB
297    real(wp),intent(in),dimension(npoints,Nlevels) :: &
298         pmol,    & ! Molecular ATB
299         pplay      ! Pressure on model levels (Pa)
300    logical,intent(in) :: &
301         ok_lidar_cfad ! True if lidar CFAD diagnostics need to be computed
302    real(wp),intent(in),dimension(npoints,nlevels) :: &
303         zlev        ! Model full levels
304    real(wp),intent(in),dimension(npoints,nlevels+1) :: &
305         zlev_half   ! Model half levels
306    real(wp),intent(in),dimension(llm) :: &
307         vgrid_z     ! mid-level altitude of the output vertical grid
308    ! Optional Inputs
309    real(wp),intent(in),dimension(npoints,ncol,Nlevels),optional :: &
310         pnorm_perp ! Lidar perpendicular ATB
311    real(wp),intent(in),dimension(npoints),optional :: &
312         surfelev   ! Surface Elevation (m)
313    real(wp),intent(in),dimension(npoints,Nlevels),optional :: &
314         tmp        ! Temperature at each levels
315   
316    ! Outputs
317    real(wp),intent(inout),dimension(npoints,llm) :: &
318         lidarcld      ! 3D "lidar" cloud fraction
319    real(wp),intent(inout),dimension(npoints,ncat) :: &
320         cldlayer      ! "lidar" cloud layer fraction (low, mid, high, total)
321    real(wp),intent(inout),dimension(npoints,max_bin,llm) :: &
322         cfad2         ! CFADs of SR
323    ! Optional Outputs
324    real(wp),intent(out),dimension(npoints,ntype),optional :: &
325         cldtype,    & ! "lidar" OPAQ type covers (opaque/thin cloud + z_opaque)
326         cldtypetemp   ! Opaque and thin clouds + z_opaque temperature
327    real(wp),intent(out),dimension(npoints,2),optional :: & 
328         cldtypemeanz  ! Opaque and thin clouds altitude
329    real(wp),intent(out),dimension(npoints,3),optional :: & 
330         cldtypemeanzse ! Opaque, thin clouds and z_opaque altitude with respect to SE
331    real(wp),intent(out),dimension(npoints),optional :: &   
332         cldthinemis   ! Thin clouds emissivity computed from SR
333    real(wp),intent(out),dimension(npoints,llm,nphase),optional :: &
334         lidarcldphase ! 3D "lidar" phase cloud fraction
335    real(wp),intent(out),dimension(npoints,llm,ntype+1),optional :: &
336         lidarcldtype ! 3D "lidar" OPAQ type fraction
337    real(wp),intent(out),dimension(npoints,40,5),optional :: &
338         lidarcldtmp   ! 3D "lidar" phase cloud fraction as a function of temp
339    real(wp),intent(out),dimension(npoints,ncat,nphase),optional :: &
340         cldlayerphase ! "lidar" phase low mid high cloud fraction
341
342    ! Local Variables
343    integer :: ic,i,j
344    logical :: lcalipso,latlid,lgrlidar532
345    real(wp),dimension(npoints,ncol,llm) :: &
346         x3d
347    real(wp),dimension(npoints,llm) :: &
348         x3d_c,pnorm_c
349    real(wp)  :: &
350         xmax
351    real(wp),dimension(npoints,1,Nlevels) :: t_in,ph_in,betamol_in
352    real(wp),dimension(npoints,ncol,llm)  :: pnormFlip,pnorm_perpFlip
353    real(wp),dimension(npoints,1,llm)     :: tmpFlip,pplayFlip,betamolFlip
354    real(wp),dimension(SR_BINS+1)         :: histBsct
355   
356    ! Which lidar platform?
357    lcalipso = .false.
358    latlid = .false.
359    lgrlidar532 = .false.
360    if (platform .eq. 'calipso') lcalipso=.true.
361    if (platform .eq. 'atlid') latlid=.true.
362    if (platform .eq. 'grlidar532') lgrlidar532=.true.
363       
364    ! Vertically regrid input data
365    if (use_vgrid) then
366       ph_in(:,1,:) = pplay(:,nlevels:1:-1)
367       call cosp_change_vertical_grid(Npoints,1,Nlevels,zlev(:,nlevels:1:-1),zlev_half(:,nlevels:1:-1),&
368            ph_in,llm,vgrid_zl(llm:1:-1),vgrid_zu(llm:1:-1),pplayFlip(:,1,llm:1:-1))
369       betamol_in(:,1,:) = pmol(:,nlevels:1:-1)
370       call cosp_change_vertical_grid(Npoints,1,Nlevels,zlev(:,nlevels:1:-1),zlev_half(:,nlevels:1:-1),&
371            betamol_in,llm,vgrid_zl(llm:1:-1),vgrid_zu(llm:1:-1),betamolFlip(:,1,llm:1:-1))
372       call cosp_change_vertical_grid(Npoints,Ncol,Nlevels,zlev(:,nlevels:1:-1),zlev_half(:,nlevels:1:-1),&
373            pnorm(:,:,nlevels:1:-1),llm,vgrid_zl(llm:1:-1),vgrid_zu(llm:1:-1),pnormFlip(:,:,llm:1:-1))
374       if (lcalipso) then
375          t_in(:,1,:)=tmp(:,nlevels:1:-1)
376          call cosp_change_vertical_grid(Npoints,1,Nlevels,zlev(:,nlevels:1:-1),zlev_half(:,nlevels:1:-1),&
377               t_in,llm,vgrid_zl(llm:1:-1),vgrid_zu(llm:1:-1),tmpFlip(:,1,llm:1:-1))
378          call cosp_change_vertical_grid(Npoints,Ncol,Nlevels,zlev(:,nlevels:1:-1),zlev_half(:,nlevels:1:-1),&
379               pnorm_perp(:,:,nlevels:1:-1),llm,vgrid_zl(llm:1:-1),vgrid_zu(llm:1:-1),pnorm_perpFlip(:,:,llm:1:-1))
380       endif
381    endif
382
383    ! Initialization (The histogram bins, are set up during initialization and the
384    ! maximum value is used as the upper bounds.)
385    if (lcalipso)    then
386       xmax = maxval(calipso_histBsct)
387       histBsct = calipso_histBsct
388    endif
389    if (latlid)      then
390       xmax = maxval(atlid_histBsct)
391       histBsct = atlid_histBsct
392    endif
393    if (lgrlidar532) then
394       xmax = maxval(grLidar532_histBsct)
395       histBsct = grLidar532_histBsct
396    endif
397       
398    ! Compute LIDAR scattering ratio
399    if (use_vgrid) then
400       do ic = 1, ncol
401          pnorm_c = pnormFlip(:,ic,:)
402          where ((pnorm_c .lt. xmax) .and. (betamolFlip(:,1,:) .lt. xmax) .and.          &
403                (betamolFlip(:,1,:) .gt. 0.0 ))
404             x3d_c = pnorm_c/betamolFlip(:,1,:)
405          elsewhere
406             x3d_c = R_UNDEF
407          end where
408          x3d(:,ic,:) = x3d_c
409       enddo
410       if (lcalipso) then
411          ! Diagnose cloud fractions for subcolumn lidar scattering ratios
412          CALL COSP_CLDFRAC(npoints,ncol,llm,ncat,nphase,tmpFlip,x3d,pnormFlip,pnorm_perpFlip,&
413               pplayFlip,S_att,S_cld,S_cld_att,R_UNDEF,lidarcld,cldlayer,lidarcldphase,&
414               cldlayerphase,lidarcldtmp)                         
415
416          ! Calipso opaque cloud diagnostics
417          CALL COSP_OPAQ(npoints,ncol,llm,ntype,tmpFlip,x3d,S_att,S_cld,R_UNDEF,lidarcldtype, &
418               cldtype,cldtypetemp,cldtypemeanz,cldtypemeanzse,cldthinemis,vgrid_z,surfelev)
419       endif
420       if (latlid) then
421          CALL COSP_CLDFRAC_NOPHASE(npoints,ncol,llm,ncat,x3d,pnormFlip,pplayFlip,  &
422               S_att_atlid,S_cld_atlid,S_cld_att_atlid,R_UNDEF,lidarcld,cldlayer)   
423       endif
424       if (lgrLidar532) then
425          CALL COSP_CLDFRAC_NOPHASE(npoints,ncol,llm,ncat,x3d,pnormFlip,pplayFlip,  &
426               S_att,S_cld,S_cld_att,R_UNDEF,lidarcld,cldlayer)
427       endif
428    else
429       do ic = 1, ncol
430          pnorm_c = pnorm(:,ic,:)
431          where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
432             x3d_c = pnorm_c/pmol
433          elsewhere
434             x3d_c = R_UNDEF
435          end where
436          x3d(:,ic,:) = x3d_c
437       enddo
438       if (lcalipso) then
439          ! Diagnose cloud fractions for subcolumn lidar scattering ratios
440          CALL COSP_CLDFRAC(npoints,ncol,nlevels,ncat,nphase,tmp,x3d,pnorm,pnorm_perp,pplay,&
441               S_att,S_cld,S_cld_att,R_UNDEF,lidarcld,cldlayer,lidarcldphase,  &
442               cldlayerphase,lidarcldtmp)
443          ! Calipso opaque cloud diagnostics
444          CALL COSP_OPAQ(npoints,ncol,nlevels,ntype,tmp,x3d,S_att,S_cld,R_UNDEF,lidarcldtype, &
445               cldtype,cldtypetemp,cldtypemeanz,cldtypemeanzse,cldthinemis,vgrid_z,surfelev)
446       endif
447       if (latlid) then
448          CALL COSP_CLDFRAC_NOPHASE(npoints,ncol,nlevels,ncat,x3d,pnorm,pplay,  &
449               S_att_atlid,S_cld_atlid,S_cld_att_atlid, R_UNDEF,lidarcld,cldlayer)
450       endif
451       if (lgrlidar532) then
452          CALL COSP_CLDFRAC_NOPHASE(npoints,ncol,nlevels,ncat,x3d,pnorm,pplay,      &
453               S_att,S_cld,S_cld_att,R_UNDEF,lidarcld,cldlayer)
454       endif
455    endif
456
457    ! CFADs
458    if (ok_lidar_cfad) then
459       ! CFADs of subgrid-scale lidar scattering ratios
460       do i=1,Npoints
461          do j=1,llm
462             cfad2(i,:,j) = hist1D(ncol,x3d(i,:,j),SR_BINS,histBsct)
463          enddo
464       enddo
465       where(cfad2 .ne. R_UNDEF) cfad2=cfad2/ncol
466    endif
467   
468    ! Unit conversions
469    where(lidarcld /= R_UNDEF)      lidarcld      = lidarcld*100._wp
470    where(cldlayer /= R_UNDEF)      cldlayer      = cldlayer*100._wp
471    if (lcalipso) then
472       where(cldtype(:,1) /= R_UNDEF)  cldtype(:,1)  = cldtype(:,1)*100._wp
473       where(cldtype(:,2) /= R_UNDEF)  cldtype(:,2)  = cldtype(:,2)*100._wp
474       where(cldlayerphase /= R_UNDEF) cldlayerphase = cldlayerphase*100._wp
475       where(lidarcldphase /= R_UNDEF) lidarcldphase = lidarcldphase*100._wp
476       where(lidarcldtype /= R_UNDEF)  lidarcldtype  = lidarcldtype*100._wp
477       where(lidarcldtmp /= R_UNDEF)   lidarcldtmp   = lidarcldtmp*100._wp
478    endif
479  end subroutine lidar_column
480
481  ! ######################################################################################
482  ! The subroutines below compute the attenuated backscatter signal and the lidar
483  ! backscatter coefficients using eq (1) from doi:0094-8276/08/2008GL034207
484  ! ######################################################################################
485  subroutine cmp_backsignal(nlev,npoints,beta,tau,pnorm)
486    ! INPUTS
487    integer, intent(in) :: nlev,npoints
488    real(wp),intent(in),dimension(npoints,nlev) :: beta,tau
489
490    ! OUTPUTS
491    real(wp),intent(out),dimension(npoints,nlev) :: pnorm
492
493    ! Internal Variables
494    real(wp), dimension(npoints) :: tautot_lay
495    integer :: k
496
497    ! Uppermost layer
498    pnorm(:,1) = beta(:,1) / (2._wp*tau(:,1)) * (1._wp-exp(-2._wp*tau(:,1)))
499
500    ! Other layers
501    do k=2,nlev
502       tautot_lay(:) = tau(:,k)-tau(:,k-1)
503       WHERE (tautot_lay(:) .gt. 0.)
504          pnorm(:,k) = beta(:,k)*EXP(-2._wp*tau(:,k-1)) /&
505               (2._wp*tautot_lay(:))*(1._wp-EXP(-2._wp*tautot_lay(:)))
506       ELSEWHERE
507          ! This must never happen, but just in case, to avoid div. by 0
508          pnorm(:,k) = beta(:,k) * EXP(-2._wp*tau(:,k-1))
509       END WHERE
510
511    END DO
512  end subroutine cmp_backsignal
513
514  subroutine cmp_beta(nlev,npoints,pnorm,tau,beta)
515    ! INPUTS
516    integer, intent(in) :: nlev,npoints
517    real(wp),intent(in),dimension(npoints,nlev) :: pnorm,tau
518
519    ! OUTPUTS
520    real(wp),intent(out),dimension(npoints,nlev) :: beta
521
522    ! Internal Variables
523    real(wp), dimension(npoints) :: tautot_lay
524    integer :: k
525    real(wp) :: epsrealwp
526
527    epsrealwp = epsilon(1._wp)
528    beta(:,1) = pnorm(:,1) * (2._wp*tau(:,1))/(1._wp-exp(-2._wp*tau(:,1)))
529    do k=2,nlev
530       tautot_lay(:) = tau(:,k)-tau(:,k-1)       
531       WHERE ( EXP(-2._wp*tau(:,k-1)) .gt. epsrealwp )
532          WHERE (tautot_lay(:) .gt. 0.)
533             beta(:,k) = pnorm(:,k)/ EXP(-2._wp*tau(:,k-1))* &
534                  (2._wp*tautot_lay(:))/(1._wp-exp(-2._wp*tautot_lay(:)))
535          ELSEWHERE
536             beta(:,k)=pnorm(:,k)/EXP(-2._wp*tau(:,k-1))
537          END WHERE
538       ELSEWHERE
539          beta(:,k)=pnorm(:,k)/epsrealwp
540       END WHERE
541    ENDDO
542
543  end subroutine cmp_beta
544    ! ####################################################################################
545    ! SUBROUTINE cosp_cldfrac
546    ! Conventions: Ncat must be equal to 4
547    ! ####################################################################################
548    SUBROUTINE COSP_CLDFRAC(Npoints,Ncolumns,Nlevels,Ncat,Nphase,tmp,x,ATB,ATBperp,      &
549                               pplay,S_att,S_cld,S_cld_att,undef,lidarcld,cldlayer,      &
550                               lidarcldphase,cldlayerphase,lidarcldtemp)
551    ! Parameters
552    integer,parameter :: Ntemp=40 ! indice of the temperature vector
553    real(wp),parameter,dimension(Ntemp+1) :: &
554       tempmod = [0.0,   183.15,186.15,189.15,192.15,195.15,198.15,201.15,204.15,207.15, &
555                  210.15,213.15,216.15,219.15,222.15,225.15,228.15,231.15,234.15,237.15, &
556                  240.15,243.15,246.15,249.15,252.15,255.15,258.15,261.15,264.15,267.15, &
557                  270.15,273.15,276.15,279.15,282.15,285.15,288.15,291.15,294.15,297.15, &
558                  473.15]
559         
560    ! Polynomial coefficient of the phase discrimination line used to separate liquid from ice
561    ! (Cesana and Chepfer, JGR, 2013)
562    ! ATBperp = ATB^5*alpha50 + ATB^4*beta50 + ATB^3*gamma50 + ATB^2*delta50 + ATB*epsilon50 + zeta50
563    real(wp),parameter :: &
564       alpha50   = 9.0322e+15_wp,  & !
565       beta50    = -2.1358e+12_wp, & !
566       gamma50   = 173.3963e06_wp, & !
567       delta50   = -3.9514e03_wp,  & !
568       epsilon50 = 0.2559_wp,      & !
569       zeta50    = -9.4776e-07_wp    !
570       
571        ! Inputs
572    integer,intent(in) :: &
573       Npoints,  & ! Number of gridpoints
574       Ncolumns, & ! Number of subcolumns
575       Nlevels,  & ! Number of vertical levels
576       Ncat,     & ! Number of cloud layer types
577       Nphase      ! Number of cloud layer phase types
578                       ! [ice,liquid,undefined,false ice,false liquid,Percent of ice]
579    real(wp),intent(in) :: &
580       S_att,    & !
581       S_cld,    & !
582       S_cld_att,& ! New threshold for undefine cloud phase detection
583       undef       ! Undefined value
584    real(wp),intent(in),dimension(Npoints,Ncolumns,Nlevels) :: &
585       x,        & !
586       ATB,      & ! 3D attenuated backscatter
587       ATBperp     ! 3D attenuated backscatter (perpendicular)
588    real(wp),intent(in),dimension(Npoints,Nlevels) :: &
589       tmp,      & ! Temperature   
590       pplay       ! Pressure
591
592        ! Outputs
593    real(wp),intent(out),dimension(Npoints,Ntemp,5) :: &
594       lidarcldtemp  ! 3D Temperature 1=tot,2=ice,3=liq,4=undef,5=ice/ice+liq
595    real(wp),intent(out),dimension(Npoints,Nlevels,Nphase) :: &
596       lidarcldphase ! 3D cloud phase fraction
597    real(wp),intent(out),dimension(Npoints,Nlevels) :: &
598       lidarcld      ! 3D cloud fraction
599    real(wp),intent(out),dimension(Npoints,Ncat) :: &
600       cldlayer      ! Low, middle, high, total cloud fractions
601    real(wp),intent(out),dimension(Npoints,Ncat,Nphase) :: &
602       cldlayerphase ! Low, middle, high, total cloud fractions for ice liquid and undefine phase   
603   
604    ! Local variables
605    integer  :: &
606       ip, k, iz, ic, ncol, nlev, i, itemp, toplvlsat
607    real(wp) :: &
608       p1,checktemp, ATBperp_tmp,checkcldlayerphase, checkcldlayerphase2
609    real(wp),dimension(Npoints,Nlevels) :: &
610       nsub,lidarcldphasetmp   
611    real(wp),dimension(Npoints,Ntemp) :: &
612       sumlidarcldtemp,lidarcldtempind
613    real(wp),dimension(Npoints,Ncolumns,Ncat) :: &
614       cldlay,nsublay   
615    real(wp),dimension(Npoints,Ncat) :: &
616       nsublayer,cldlayerphasetmp,cldlayerphasesum
617    real(wp),dimension(Npoints,Ncolumns,Nlevels) :: &   
618       tmpi, & ! Temperature of ice cld
619       tmpl, & ! Temperature of liquid cld
620       tmpu, & ! Temperature of undef cld
621       cldy, & !
622       srok    !
623    real(wp),dimension(Npoints,Ncolumns,Ncat,Nphase) :: &
624       cldlayphase ! subgrided low mid high phase cloud fraction
625             
626    ! ####################################################################################
627        ! 1) Initialize   
628    ! ####################################################################################
629    lidarcld              = 0._wp
630    nsub                  = 0._wp
631    cldlay                = 0._wp
632    nsublay               = 0._wp
633    ATBperp_tmp           = 0._wp
634    lidarcldphase(:,:,:)  = 0._wp
635    cldlayphase(:,:,:,:)  = 0._wp
636    cldlayerphase(:,:,:)  = 0._wp
637    tmpi(:,:,:)           = 0._wp
638    tmpl(:,:,:)           = 0._wp
639    tmpu(:,:,:)           = 0._wp
640    cldlayerphasesum(:,:) = 0._wp
641    lidarcldtemp(:,:,:)   = 0._wp
642    lidarcldtempind(:,:)  = 0._wp
643    sumlidarcldtemp(:,:)  = 0._wp
644    lidarcldphasetmp(:,:) = 0._wp
645    toplvlsat             = 0
646
647    ! ####################################################################################
648    ! 2) Cloud detection
649    ! ####################################################################################
650    do k=1,Nlevels
651       ! Cloud detection at subgrid-scale:
652       where ((x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
653          cldy(:,:,k)=1._wp
654       elsewhere
655          cldy(:,:,k)=0._wp
656       endwhere
657       
658       ! Number of usefull sub-columns:
659       where ((x(:,:,k) .gt. S_att) .and. (x(:,:,k) .ne. undef) )
660          srok(:,:,k)=1._wp
661       elsewhere
662          srok(:,:,k)=0._wp
663       endwhere
664    enddo   
665   
666    ! ####################################################################################
667    ! 3) Grid-box 3D cloud fraction and layered cloud fractions(ISCCP pressure categories)
668    ! ####################################################################################
669    lidarcld = 0._wp
670    nsub     = 0._wp
671    cldlay   = 0._wp
672    nsublay  = 0._wp
673    do k=1,Nlevels
674       do ic = 1, Ncolumns
675          do ip = 1, Npoints
676         
677             ! Computation of the cloud fraction as a function of the temperature instead
678             ! of height, for ice,liquid and all clouds
679             if(srok(ip,ic,k).gt.0.)then
680                do itemp=1,Ntemp
681                   if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
682                      lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1._wp
683                   endif
684                enddo
685             endif
686             
687             if(cldy(ip,ic,k).eq.1.)then
688                do itemp=1,Ntemp
689                   if( (tmp(ip,k) .ge. tempmod(itemp)).and.(tmp(ip,k) .lt. tempmod(itemp+1)) )then
690                      lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1._wp
691                   endif
692                enddo
693             endif
694
695             iz=1
696             p1 = pplay(ip,k)
697             if ( p1.gt.0. .and. p1.lt.(440._wp*100._wp)) then ! high clouds
698                iz=3
699             else if(p1.ge.(440._wp*100._wp) .and. p1.lt.(680._wp*100._wp)) then ! mid clouds
700                iz=2
701             endif
702             
703             cldlay(ip,ic,iz) = MAX(cldlay(ip,ic,iz),cldy(ip,ic,k))
704             cldlay(ip,ic,4)  = MAX(cldlay(ip,ic,4),cldy(ip,ic,k))
705             lidarcld(ip,k)   = lidarcld(ip,k) + cldy(ip,ic,k)
706             
707             nsublay(ip,ic,iz) = MAX(nsublay(ip,ic,iz),srok(ip,ic,k))
708             nsublay(ip,ic,4)  = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
709             nsub(ip,k)        = nsub(ip,k) + srok(ip,ic,k)
710             
711          enddo
712       enddo
713    enddo   
714   
715    ! Grid-box 3D cloud fraction
716    where ( nsub(:,:).gt.0.0 )
717       lidarcld(:,:) = lidarcld(:,:)/nsub(:,:)
718    elsewhere
719       lidarcld(:,:) = undef
720    endwhere
721   
722    ! Layered cloud fractions
723    cldlayer  = 0._wp
724    nsublayer = 0._wp
725    do iz = 1, Ncat
726       do ic = 1, Ncolumns
727          cldlayer(:,iz)  = cldlayer(:,iz)  + cldlay(:,ic,iz)
728          nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz)
729       enddo
730    enddo
731    where (nsublayer(:,:) .gt. 0.0)
732       cldlayer(:,:) = cldlayer(:,:)/nsublayer(:,:)
733    elsewhere
734       cldlayer(:,:) = undef
735    endwhere
736             
737    ! ####################################################################################
738    ! 4) Grid-box 3D cloud Phase
739    ! ####################################################################################
740   
741    ! ####################################################################################
742    ! 4.1) For Cloudy pixels with 8.16km < z < 19.2km
743    ! ####################################################################################
744    do ncol=1,Ncolumns
745       do i=1,Npoints         
746          do nlev=1,23 ! from 19.2km until 8.16km
747               p1 = pplay(1,nlev)
748
749             ! Avoid zero values
750             if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
751                ! Computation of the ATBperp along the phase discrimination line
752                ATBperp_tmp = (ATB(i,ncol,nlev)**5)*alpha50 + (ATB(i,ncol,nlev)**4)*beta50 + &
753                     (ATB(i,ncol,nlev)**3)*gamma50 + (ATB(i,ncol,nlev)**2)*delta50 + &
754                     ATB(i,ncol,nlev)*epsilon50 + zeta50     
755                ! ########################################################################
756                ! 4.1.a) Ice: ATBperp above the phase discrimination line
757                ! ########################################################################
758                if((ATBperp(i,ncol,nlev)-ATBperp_tmp) .ge. 0.)then ! Ice clouds
759
760                   ! ICE with temperature above 273,15°K = Liquid (false ice)
761                   if(tmp(i,nlev) .gt. 273.15) then ! Temperature above 273,15 K
762                     ! Liquid: False ice corrected by the temperature to Liquid
763                      lidarcldphase(i,nlev,2) = lidarcldphase(i,nlev,2)+1._wp ! False ice detection ==> added to Liquid
764                                   
765                      tmpl(i,ncol,nlev)       = tmp(i,nlev)
766                      lidarcldphase(i,nlev,5) = lidarcldphase(i,nlev,5)+1._wp ! Keep the information "temperature criterium used"                     
767                                                                              ! to classify the phase cloud
768                      cldlayphase(i,ncol,4,2) = 1. ! tot cloud
769                      if (p1 .gt. 0. .and. p1.lt.(440._wp*100._wp)) then ! high cloud
770                         cldlayphase(i,ncol,3,2) = 1._wp
771                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then ! mid cloud
772                         cldlayphase(i,ncol,2,2) = 1._wp
773                      else ! low cloud
774                         cldlayphase(i,ncol,1,2) = 1._wp
775                      endif
776                      cldlayphase(i,ncol,4,5) = 1._wp ! tot cloud
777                      ! High cloud
778                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
779                         cldlayphase(i,ncol,3,5) = 1._wp
780                      ! Middle cloud
781                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
782                         cldlayphase(i,ncol,2,5) = 1._wp
783                      ! Low cloud
784                      else
785                         cldlayphase(i,ncol,1,5) = 1._wp
786                      endif
787                   else
788                      ! ICE with temperature below 273,15°K
789                      lidarcldphase(i,nlev,1) = lidarcldphase(i,nlev,1)+1._wp
790                      tmpi(i,ncol,nlev)       = tmp(i,nlev)
791                      cldlayphase(i,ncol,4,1) = 1._wp ! tot cloud
792                      ! High cloud
793                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
794                         cldlayphase(i,ncol,3,1) = 1._wp
795                      ! Middle cloud   
796                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
797                         cldlayphase(i,ncol,2,1) = 1._wp
798                      ! Low cloud
799                      else
800                         cldlayphase(i,ncol,1,1) = 1._wp
801                      endif
802                   endif
803                ! ########################################################################
804                ! 4.1.b) Liquid: ATBperp below the phase discrimination line
805                ! ########################################################################
806                else
807                   ! Liquid with temperature above 231,15°K
808                   if(tmp(i,nlev) .gt. 231.15_wp) then
809                      lidarcldphase(i,nlev,2) = lidarcldphase(i,nlev,2)+1._wp
810                      tmpl(i,ncol,nlev)       = tmp(i,nlev)
811                      cldlayphase(i,ncol,4,2) = 1._wp ! tot cloud
812                      ! High cloud
813                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
814                         cldlayphase(i,ncol,3,2) = 1._wp
815                      ! Middle cloud   
816                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
817                         cldlayphase(i,ncol,2,2) = 1._wp
818                      ! Low cloud   
819                      else
820                         cldlayphase(i,ncol,1,2) = 1._wp
821                      endif
822                   else
823                      ! Liquid with temperature below 231,15°K = Ice (false liquid)
824                      tmpi(i,ncol,nlev)       = tmp(i,nlev)
825                      lidarcldphase(i,nlev,1) = lidarcldphase(i,nlev,1)+1._wp ! false liquid detection ==> added to ice
826                      lidarcldphase(i,nlev,4) = lidarcldphase(i,nlev,4)+1._wp
827                      cldlayphase(i,ncol,4,4) = 1._wp ! tot cloud
828                      ! High cloud
829                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
830                         cldlayphase(i,ncol,3,4) = 1._wp
831                      ! Middle cloud   
832                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
833                         cldlayphase(i,ncol,2,4) = 1._wp
834                      ! Low cloud
835                      else
836                         cldlayphase(i,ncol,1,4) = 1._wp
837                      endif
838                      cldlayphase(i,ncol,4,1) = 1._wp ! tot cloud
839                      ! High cloud
840                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
841                         cldlayphase(i,ncol,3,1) = 1._wp
842                      ! Middle cloud   
843                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
844                         cldlayphase(i,ncol,2,1) = 1._wp
845                      ! Low cloud   
846                      else
847                         cldlayphase(i,ncol,1,1) = 1._wp
848                      endif
849                   endif
850                endif ! end of discrimination condition
851             endif ! end of cloud condition
852          enddo ! end of altitude loop
853
854          ! ##############################################################################
855          ! 4.2) For Cloudy pixels with 0km < z < 8.16km
856          ! ##############################################################################
857          toplvlsat = 0
858          do nlev=24,Nlevels! from 8.16km until 0km
859             p1 = pplay(i,nlev)
860
861             if((cldy(i,ncol,nlev) .eq. 1.) .and. (ATBperp(i,ncol,nlev) .gt. 0.) )then
862                ! Computation of the ATBperp of the phase discrimination line
863                ATBperp_tmp = (ATB(i,ncol,nlev)**5)*alpha50 + (ATB(i,ncol,nlev)**4)*beta50 + &
864                     (ATB(i,ncol,nlev)**3)*gamma50 + (ATB(i,ncol,nlev)**2)*delta50 + &
865                     ATB(i,ncol,nlev)*epsilon50 + zeta50
866                ! ########################################################################
867                ! 4.2.a) Ice: ATBperp above the phase discrimination line
868                ! ########################################################################
869                ! ICE with temperature above 273,15°K = Liquid (false ice)
870                if((ATBperp(i,ncol,nlev)-ATBperp_tmp) .ge. 0.)then ! Ice clouds
871                   if(tmp(i,nlev) .gt. 273.15)then
872                      lidarcldphase(i,nlev,2) = lidarcldphase(i,nlev,2)+1._wp ! false ice ==> liq
873                      tmpl(i,ncol,nlev)       = tmp(i,nlev)
874                      lidarcldphase(i,nlev,5) = lidarcldphase(i,nlev,5)+1._wp
875                      cldlayphase(i,ncol,4,2) = 1._wp ! tot cloud
876                      ! High cloud
877                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
878                         cldlayphase(i,ncol,3,2) = 1._wp
879                      ! Middle cloud   
880                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
881                         cldlayphase(i,ncol,2,2) = 1._wp
882                      ! Low cloud
883                      else
884                         cldlayphase(i,ncol,1,2) = 1._wp
885                      endif
886                     
887                      cldlayphase(i,ncol,4,5) = 1. ! tot cloud
888                      ! High cloud
889                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
890                         cldlayphase(i,ncol,3,5) = 1._wp
891                      ! Middle cloud   
892                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
893                         cldlayphase(i,ncol,2,5) = 1._wp
894                      ! Low cloud   
895                      else
896                         cldlayphase(i,ncol,1,5) = 1._wp
897                      endif
898                   else
899                      ! ICE with temperature below 273,15°K
900                      lidarcldphase(i,nlev,1) = lidarcldphase(i,nlev,1)+1._wp
901                     tmpi(i,ncol,nlev)       = tmp(i,nlev)
902                      cldlayphase(i,ncol,4,1) = 1._wp ! tot cloud
903                      ! High cloud
904                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
905                         cldlayphase(i,ncol,3,1) = 1._wp
906                      ! Middle cloud   
907                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt.(680._wp*100._wp)) then
908                         cldlayphase(i,ncol,2,1) = 1._wp
909                      ! Low cloud   
910                      else
911                         cldlayphase(i,ncol,1,1) = 1._wp
912                      endif
913                   endif
914                   
915                ! ########################################################################
916                ! 4.2.b) Liquid: ATBperp below the phase discrimination line
917                ! ########################################################################
918                else
919                   ! Liquid with temperature above 231,15°K
920                   if(tmp(i,nlev) .gt. 231.15)then
921                      lidarcldphase(i,nlev,2) = lidarcldphase(i,nlev,2)+1._wp
922                      tmpl(i,ncol,nlev)       = tmp(i,nlev)
923                      cldlayphase(i,ncol,4,2) = 1._wp ! tot cloud
924                      ! High cloud
925                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
926                         cldlayphase(i,ncol,3,2) = 1._wp
927                      ! Middle cloud   
928                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
929                         cldlayphase(i,ncol,2,2) = 1._wp
930                      ! Low cloud   
931                      else
932                         cldlayphase(i,ncol,1,2) = 1._wp
933                      endif
934                   else
935                      ! Liquid with temperature below 231,15°K = Ice (false liquid)
936                      tmpi(i,ncol,nlev)       = tmp(i,nlev)
937                      lidarcldphase(i,nlev,1) = lidarcldphase(i,nlev,1)+1._wp ! false liq ==> ice
938                      lidarcldphase(i,nlev,4) = lidarcldphase(i,nlev,4)+1._wp ! false liq ==> ice
939                      cldlayphase(i,ncol,4,4) = 1._wp ! tot cloud
940                      ! High cloud
941                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
942                         cldlayphase(i,ncol,3,4) = 1._wp
943                      ! Middle   
944                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
945                         cldlayphase(i,ncol,2,4) = 1._wp
946                      ! Low cloud   
947                      else
948                         cldlayphase(i,ncol,1,4) = 1._wp
949                      endif
950                     
951                      cldlayphase(i,ncol,4,1) = 1._wp ! tot cloud
952                      ! High cloud
953                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
954                         cldlayphase(i,ncol,3,1) = 1._wp
955                      ! Middle cloud   
956                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
957                         cldlayphase(i,ncol,2,1) = 1._wp
958                      ! Low cloud   
959                      else
960                         cldlayphase(i,ncol,1,1) = 1._wp
961                      endif
962                   endif
963                endif ! end of discrimination condition
964               
965                toplvlsat=0
966               
967                ! Find the level of the highest cloud with SR>30
968                if(x(i,ncol,nlev) .gt. S_cld_att) then ! SR > 30.
969                    toplvlsat = nlev+1
970                    goto 99
971                endif
972             endif ! end of cloud condition
973          enddo ! end of altitude loop
97499        continue
975         
976          ! ##############################################################################
977          ! Undefined phase: For a cloud located below another cloud with SR>30
978          ! see Cesana and Chepfer 2013 Sect.III.2
979          ! ##############################################################################
980          if(toplvlsat.ne.0) then
981             do nlev = toplvlsat,Nlevels
982                p1 = pplay(i,nlev)
983                if(cldy(i,ncol,nlev).eq.1.)then
984                   lidarcldphase(i,nlev,3) = lidarcldphase(i,nlev,3)+1._wp
985                   tmpu(i,ncol,nlev)       = tmp(i,nlev)
986                   cldlayphase(i,ncol,4,3) = 1._wp ! tot cloud
987                   ! High cloud
988                   if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
989                      cldlayphase(i,ncol,3,3) = 1._wp
990                   ! Middle cloud   
991                   else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
992                      cldlayphase(i,ncol,2,3) = 1._wp
993                   ! Low cloud   
994                   else
995                      cldlayphase(i,ncol,1,3) = 1._wp
996                   endif
997                endif
998             enddo
999          endif
1000          toplvlsat=0
1001       enddo
1002    enddo
1003     
1004    ! ####################################################################################
1005    ! Computation of final cloud phase diagnosis
1006    ! ####################################################################################
1007
1008    ! Compute the Ice percentage in cloud = ice/(ice+liq) as a function of the occurrences
1009    lidarcldphasetmp(:,:) = lidarcldphase(:,:,1)+lidarcldphase(:,:,2);
1010    WHERE (lidarcldphasetmp(:,:) .gt. 0.)
1011       lidarcldphase(:,:,6)=lidarcldphase(:,:,1)/lidarcldphasetmp(:,:)
1012    ELSEWHERE
1013       lidarcldphase(:,:,6) = undef
1014    ENDWHERE
1015   
1016    ! Compute Phase 3D Cloud Fraction
1017    !WHERE (nsub(:,Nlevels:1:-1) .gt. 0.0 )
1018    WHERE (nsub(:,:) .gt. 0.0 ) 
1019       lidarcldphase(:,:,1)=lidarcldphase(:,:,1)/nsub(:,:)
1020       lidarcldphase(:,:,2)=lidarcldphase(:,:,2)/nsub(:,:)
1021       lidarcldphase(:,:,3)=lidarcldphase(:,:,3)/nsub(:,:)
1022       lidarcldphase(:,:,4)=lidarcldphase(:,:,4)/nsub(:,:)
1023       lidarcldphase(:,:,5)=lidarcldphase(:,:,5)/nsub(:,:)
1024    ELSEWHERE
1025       lidarcldphase(:,:,1) = undef
1026       lidarcldphase(:,:,2) = undef
1027       lidarcldphase(:,:,3) = undef
1028       lidarcldphase(:,:,4) = undef
1029       lidarcldphase(:,:,5) = undef
1030    ENDWHERE
1031
1032    ! Compute Phase low mid high cloud fractions
1033    do iz = 1, Ncat
1034       do i=1,Nphase-3
1035          do ic = 1, Ncolumns
1036             cldlayerphase(:,iz,i)  = cldlayerphase(:,iz,i)  + cldlayphase(:,ic,iz,i)
1037             cldlayerphasesum(:,iz) = cldlayerphasesum(:,iz) + cldlayphase(:,ic,iz,i)
1038          enddo
1039       enddo
1040    enddo
1041    do iz = 1, Ncat
1042       do i=4,5
1043          do ic = 1, Ncolumns
1044             cldlayerphase(:,iz,i) = cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i)
1045          enddo
1046       enddo
1047    enddo
1048   
1049    ! Compute the Ice percentage in cloud = ice/(ice+liq)
1050    cldlayerphasetmp(:,:)=cldlayerphase(:,:,1)+cldlayerphase(:,:,2)
1051    WHERE (cldlayerphasetmp(:,:).gt. 0.)
1052       cldlayerphase(:,:,6)=cldlayerphase(:,:,1)/cldlayerphasetmp(:,:)
1053    ELSEWHERE
1054       cldlayerphase(:,:,6) = undef
1055    ENDWHERE
1056   
1057    do i=1,Nphase-1
1058       WHERE ( cldlayerphasesum(:,:).gt.0.0 )
1059          cldlayerphase(:,:,i) = (cldlayerphase(:,:,i)/cldlayerphasesum(:,:)) * cldlayer(:,:)
1060       ENDWHERE
1061    enddo
1062   
1063    do i=1,Npoints
1064       do iz=1,Ncat
1065          checkcldlayerphase=0.
1066          checkcldlayerphase2=0.
1067          if (cldlayerphasesum(i,iz) .gt. 0.0 )then
1068             do ic=1,Nphase-3
1069                checkcldlayerphase = checkcldlayerphase+cldlayerphase(i,iz,ic)
1070             enddo
1071             checkcldlayerphase2 = cldlayer(i,iz)-checkcldlayerphase
1072             if((checkcldlayerphase2 .gt. 0.01) .or. (checkcldlayerphase2 .lt. -0.01) ) print *, checkcldlayerphase,cldlayer(i,iz)
1073          endif
1074       enddo
1075    enddo
1076   
1077    do i=1,Nphase-1
1078       WHERE (nsublayer(:,:) .eq. 0.0)
1079          cldlayerphase(:,:,i) = undef
1080       ENDWHERE
1081    enddo
1082 
1083    ! Compute Phase 3D as a function of temperature
1084    do nlev=1,Nlevels
1085       do ncol=1,Ncolumns
1086          do i=1,Npoints
1087             do itemp=1,Ntemp
1088                if(tmpi(i,ncol,nlev).gt.0.)then
1089                   if((tmpi(i,ncol,nlev) .ge. tempmod(itemp)) .and. (tmpi(i,ncol,nlev) .lt. tempmod(itemp+1)) )then
1090                      lidarcldtemp(i,itemp,2)=lidarcldtemp(i,itemp,2)+1._wp
1091                   endif
1092                elseif(tmpl(i,ncol,nlev) .gt. 0.)then
1093                   if((tmpl(i,ncol,nlev) .ge. tempmod(itemp)) .and. (tmpl(i,ncol,nlev) .lt. tempmod(itemp+1)) )then
1094                      lidarcldtemp(i,itemp,3)=lidarcldtemp(i,itemp,3)+1._wp
1095                   endif
1096                elseif(tmpu(i,ncol,nlev) .gt. 0.)then
1097                   if((tmpu(i,ncol,nlev) .ge. tempmod(itemp)) .and. (tmpu(i,ncol,nlev) .lt. tempmod(itemp+1)) )then
1098                      lidarcldtemp(i,itemp,4)=lidarcldtemp(i,itemp,4)+1._wp
1099                   endif
1100                endif
1101             enddo
1102          enddo
1103       enddo
1104    enddo
1105   
1106    ! Check temperature cloud fraction
1107    do i=1,Npoints
1108       do itemp=1,Ntemp
1109          checktemp=lidarcldtemp(i,itemp,2)+lidarcldtemp(i,itemp,3)+lidarcldtemp(i,itemp,4)
1110          !if(checktemp .NE. lidarcldtemp(i,itemp,1))then
1111          !   print *, i,itemp
1112          !   print *, lidarcldtemp(i,itemp,1:4)
1113          !endif
1114         
1115       enddo
1116    enddo
1117   
1118    ! Compute the Ice percentage in cloud = ice/(ice+liq)
1119    sumlidarcldtemp(:,:)=lidarcldtemp(:,:,2)+lidarcldtemp(:,:,3)   
1120    WHERE(sumlidarcldtemp(:,:) .gt. 0.)
1121       lidarcldtemp(:,:,5)=lidarcldtemp(:,:,2)/sumlidarcldtemp(:,:)
1122    ELSEWHERE
1123       lidarcldtemp(:,:,5)=undef
1124    ENDWHERE
1125   
1126    do i=1,4
1127       WHERE(lidarcldtempind(:,:) .gt. 0.)
1128          lidarcldtemp(:,:,i) = lidarcldtemp(:,:,i)/lidarcldtempind(:,:)
1129       ELSEWHERE
1130          lidarcldtemp(:,:,i) = undef
1131       ENDWHERE
1132    enddo
1133   
1134    RETURN
1135  END SUBROUTINE COSP_CLDFRAC
1136 
1137  ! ####################################################################################
1138  ! SUBROUTINE cosp_cldfrac_nophase
1139  ! Conventions: Ncat must be equal to 4
1140  ! ####################################################################################
1141  SUBROUTINE COSP_CLDFRAC_NOPHASE(Npoints,Ncolumns,Nlevels,Ncat,x,ATB,pplay,      &
1142       S_att,S_cld,S_cld_att,undef,lidarcld,cldlayer)
1143   
1144    ! Inputs
1145    integer,intent(in) :: &
1146         Npoints,  & ! Number of gridpoints
1147         Ncolumns, & ! Number of subcolumns
1148         Nlevels,  & ! Number of vertical levels
1149         Ncat        ! Number of cloud layer types
1150    real(wp),intent(in) :: &
1151         S_att,    & !
1152         S_cld,    & !
1153         S_cld_att,& ! New threshold for undefine cloud phase detection
1154         undef       ! Undefined value
1155    real(wp),intent(in),dimension(Npoints,Ncolumns,Nlevels) :: &
1156         x,        & !
1157         ATB         ! 3D attenuated backscatter
1158    real(wp),intent(in),dimension(Npoints,Nlevels) :: &
1159         pplay       ! Pressure
1160   
1161    ! Outputs
1162    real(wp),intent(out),dimension(Npoints,Nlevels) :: &
1163         lidarcld      ! 3D cloud fraction
1164    real(wp),intent(out),dimension(Npoints,Ncat) :: &
1165         cldlayer      ! Low, middle, high, total cloud fractions
1166   
1167    ! Local variables
1168    integer  :: &
1169         ip, k, iz, ic, ncol, nlev, i
1170    real(wp) :: &
1171         p1
1172    real(wp),dimension(Npoints,Nlevels) :: &
1173         nsub
1174    real(wp),dimension(Npoints,Ncolumns,Ncat) :: &
1175         cldlay,nsublay   
1176    real(wp),dimension(Npoints,Ncat) :: &
1177         nsublayer
1178    real(wp),dimension(Npoints,Ncolumns,Nlevels) :: &   
1179         cldy, & !
1180         srok    !
1181   
1182    ! ####################################################################################
1183    ! 1) Initialize   
1184    ! ####################################################################################
1185    lidarcld           = 0._wp
1186    nsub               = 0._wp
1187    cldlay             = 0._wp
1188    nsublay            = 0._wp
1189   
1190    ! ####################################################################################
1191    ! 2) Cloud detection
1192    ! ####################################################################################
1193    do k=1,Nlevels
1194       ! Cloud detection at subgrid-scale:
1195       where ((x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
1196          cldy(:,:,k)=1._wp
1197       elsewhere
1198          cldy(:,:,k)=0._wp
1199       endwhere
1200       
1201       ! Number of usefull sub-columns:
1202       where ((x(:,:,k) .gt. S_att) .and. (x(:,:,k) .ne. undef) )
1203          srok(:,:,k)=1._wp
1204       elsewhere
1205          srok(:,:,k)=0._wp
1206       endwhere
1207    enddo   
1208
1209    ! ####################################################################################
1210    ! 3) Grid-box 3D cloud fraction and layered cloud fractions(ISCCP pressure categories)
1211    ! ####################################################################################
1212    do k=1,Nlevels
1213       do ic = 1, Ncolumns
1214          do ip = 1, Npoints
1215
1216             iz=1
1217             p1 = pplay(ip,k)
1218             if ( p1.gt.0. .and. p1.lt.(440._wp*100._wp)) then ! high clouds
1219                iz=3
1220             else if(p1.ge.(440._wp*100._wp) .and. p1.lt.(680._wp*100._wp)) then ! mid clouds
1221                iz=2
1222             endif
1223             
1224             cldlay(ip,ic,iz) = MAX(cldlay(ip,ic,iz),cldy(ip,ic,k))
1225             cldlay(ip,ic,4)  = MAX(cldlay(ip,ic,4),cldy(ip,ic,k))
1226             lidarcld(ip,k)   = lidarcld(ip,k) + cldy(ip,ic,k)
1227             
1228             nsublay(ip,ic,iz) = MAX(nsublay(ip,ic,iz),srok(ip,ic,k))
1229             nsublay(ip,ic,4)  = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
1230             nsub(ip,k)        = nsub(ip,k) + srok(ip,ic,k)
1231             
1232          enddo
1233       enddo
1234    enddo   
1235   
1236    ! Grid-box 3D cloud fraction
1237    where ( nsub(:,:).gt.0.0 )
1238       lidarcld(:,:) = lidarcld(:,:)/nsub(:,:)
1239    elsewhere
1240       lidarcld(:,:) = undef
1241    endwhere
1242   
1243    ! Layered cloud fractions
1244    cldlayer  = 0._wp
1245    nsublayer = 0._wp
1246    do iz = 1, Ncat
1247       do ic = 1, Ncolumns
1248          cldlayer(:,iz)  = cldlayer(:,iz)  + cldlay(:,ic,iz)
1249          nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz)
1250       enddo
1251    enddo
1252    where (nsublayer(:,:) .gt. 0.0)
1253       cldlayer(:,:) = cldlayer(:,:)/nsublayer(:,:)
1254    elsewhere
1255       cldlayer(:,:) = undef
1256    endwhere
1257
1258    RETURN
1259  END SUBROUTINE COSP_CLDFRAC_NOPHASE
1260
1261    ! ####################################################################################
1262    ! SUBROUTINE cosp_opaq
1263    ! Conventions: Ntype must be equal to 3
1264    ! ####################################################################################
1265    SUBROUTINE COSP_OPAQ(Npoints,Ncolumns,Nlevels,Ntype,tmp,x,S_att,S_cld,undef,lidarcldtype,   &
1266                         cldtype,cldtypetemp,cldtypemeanz,cldtypemeanzse,cldthinemis,vgrid_z,   &
1267                         surfelev)
1268
1269    ! Local parameter
1270    real(wp),parameter  :: &
1271       S_att_opaq = 0.06_wp, & ! Fully Attenuated threshold (Guzman et al. 2017, JGR-Atmospheres)
1272       eta = 0.6_wp            ! Multiple-scattering factor (Vaillant de Guelis et al. 2017a, AMT)
1273
1274        ! Inputs
1275    integer,intent(in) :: &
1276       Npoints,  & ! Number of gridpoints
1277       Ncolumns, & ! Number of subcolumns
1278       Nlevels,  & ! Number of vertical levels
1279       Ntype       ! Number of OPAQ cloud types (opaque, thin clouds and z_opaque)
1280    real(wp),intent(in) :: &
1281       S_att,    & ! Fully Attenuated legacy threshold
1282       S_cld,    & ! Cloud detection threshold
1283       undef       ! Undefined value
1284    real(wp),intent(in),dimension(Nlevels) :: &
1285       vgrid_z     ! mid-level vertical profile altitude (subcolumns)
1286    real(wp),intent(in),dimension(Npoints,Ncolumns,Nlevels) :: &
1287       x           ! SR profiles (subcolumns)
1288    real(wp),intent(in),dimension(Npoints,Nlevels) :: &
1289       tmp         ! Temperature profiles
1290    real(wp),intent(in),dimension(Npoints) :: &
1291       surfelev    ! Surface Elevation (SE)
1292
1293        ! Outputs
1294    real(wp),intent(out),dimension(Npoints,Nlevels,Ntype+1) :: &
1295       lidarcldtype   ! 3D OPAQ product fraction (opaque clouds, thin clouds, z_opaque, opacity)
1296    real(wp),intent(out),dimension(Npoints,Ntype) :: &
1297       cldtype,     & ! Opaque/thin cloud covers + z_opaque altitude
1298       cldtypetemp    ! Opaque and thin clouds + z_opaque temperature
1299    real(wp),intent(out),dimension(Npoints,2) :: &
1300       cldtypemeanz   ! Opaque and thin clouds altitude
1301    real(wp),intent(out),dimension(Npoints,3) :: &
1302       cldtypemeanzse ! Opaque, thin clouds and z_opaque altitude with respect to SE
1303    real(wp),intent(out),dimension(Npoints) :: &
1304       cldthinemis    ! Thin clouds emissivity
1305   
1306    ! Local variables
1307    integer  :: &
1308       ip, k, zopac, ic, iz, z_top, z_base, topcloud
1309    real(wp)  :: &
1310       srmean, srcount, trans2, tau_app, tau_vis, tau_ir, cloudemis
1311    real(wp),dimension(Npoints) :: &
1312       count_emis
1313    real(wp),dimension(Npoints,Nlevels) :: &
1314       nsub, nsubopaq
1315    real(wp),dimension(Npoints,Ncolumns,Ntype+1) :: & ! Opaque, thin, z_opaque and all cloud cover
1316       cldlay, nsublay   
1317    real(wp),dimension(Npoints,Ntype) :: &
1318       nsublayer
1319    real(wp),dimension(Npoints,Ncolumns,Nlevels) :: &   
1320       cldy,     & !
1321       cldyopaq, & !
1322       srok,     & !
1323       srokopaq    !
1324
1325    ! ####################################################################################
1326        ! 1) Initialize   
1327    ! ####################################################################################
1328    cldtype(:,:)          = 0._wp
1329    cldtypetemp(:,:)      = 0._wp
1330    cldtypemeanz(:,:)     = 0._wp
1331    cldtypemeanzse(:,:)   = 0._wp
1332    cldthinemis(:)        = 0._wp
1333    count_emis(:)         = 0._wp
1334    lidarcldtype(:,:,:)   = 0._wp
1335    nsub                  = 0._wp
1336    nsubopaq              = 0._wp
1337    cldlay                = 0._wp
1338    nsublay               = 0._wp
1339    nsublayer             = 0._wp
1340
1341    ! ####################################################################################
1342    ! 2) Cloud detection and Fully attenuated layer detection
1343    ! ####################################################################################
1344    do k=1,Nlevels
1345       ! Cloud detection at subgrid-scale:
1346       where ( (x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
1347          cldy(:,:,k)=1._wp
1348       elsewhere
1349          cldy(:,:,k)=0._wp
1350       endwhere
1351       ! Fully attenuated layer detection at subgrid-scale:
1352       where ( (x(:,:,k) .lt. S_att_opaq) .and. (x(:,:,k) .ge. 0.) .and. (x(:,:,k) .ne. undef) ) !DEBUG
1353          cldyopaq(:,:,k)=1._wp
1354       elsewhere
1355          cldyopaq(:,:,k)=0._wp
1356       endwhere
1357
1358
1359       ! Number of usefull sub-column layers:
1360       where ( (x(:,:,k) .gt. S_att) .and. (x(:,:,k) .ne. undef) )
1361          srok(:,:,k)=1._wp
1362       elsewhere
1363          srok(:,:,k)=0._wp
1364       endwhere
1365       ! Number of usefull sub-columns layers for z_opaque 3D fraction:
1366       where ( (x(:,:,k) .ge. 0.) .and. (x(:,:,k) .ne. undef) ) !DEBUG
1367          srokopaq(:,:,k)=1._wp
1368       elsewhere
1369          srokopaq(:,:,k)=0._wp
1370       endwhere
1371    enddo
1372
1373    ! ####################################################################################
1374    ! 3) Grid-box 3D OPAQ product fraction and cloud type cover (opaque/thin) + mean z_opaque
1375    ! ####################################################################################
1376
1377    do k=1,Nlevels
1378       do ic = 1, Ncolumns
1379          do ip = 1, Npoints
1380
1381             cldlay(ip,ic,1)   = MAX(cldlay(ip,ic,1),cldyopaq(ip,ic,k)) ! Opaque cloud
1382             cldlay(ip,ic,4)   = MAX(cldlay(ip,ic,4),cldy(ip,ic,k))     ! All cloud
1383
1384             nsublay(ip,ic,1)  = MAX(nsublay(ip,ic,1),srok(ip,ic,k))
1385             nsublay(ip,ic,2)  = MAX(nsublay(ip,ic,2),srok(ip,ic,k))
1386!             nsublay(ip,ic,4)  = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
1387             nsub(ip,k)        = nsub(ip,k) + srok(ip,ic,k)
1388             nsubopaq(ip,k)    = nsubopaq(ip,k) + srokopaq(ip,ic,k)
1389
1390          enddo
1391       enddo
1392    enddo   
1393
1394! OPAQ variables
1395     do ic = 1, Ncolumns
1396        do ip = 1, Npoints
1397
1398     ! Declaring non-opaque cloudy profiles as thin cloud profiles
1399           if ( cldlay(ip,ic,4).gt. 0. .and. cldlay(ip,ic,1) .eq. 0. ) then
1400              cldlay(ip,ic,2)  =  1._wp
1401           endif
1402
1403     ! Filling in 3D and 2D variables
1404
1405     ! Opaque cloud profiles
1406           if ( cldlay(ip,ic,1) .eq. 1. ) then
1407              zopac = 0._wp
1408              z_top = 0._wp
1409              do k=1,Nlevels-1
1410     ! Declaring z_opaque altitude and opaque cloud fraction for 3D and 2D variables
1411     ! From SFC-2-TOA ( actually from vgrid_z(SFC+1) = vgrid_z(Nlevels-1) )
1412                 if ( cldy(ip,ic,Nlevels-k) .eq. 1. .and. zopac .eq. 0. ) then
1413                    lidarcldtype(ip,Nlevels-k + 1,3) = lidarcldtype(ip,Nlevels-k + 1,3) + 1._wp
1414                    cldlay(ip,ic,3)                  = vgrid_z(Nlevels-k+1)      ! z_opaque altitude
1415                    nsublay(ip,ic,3)                 = 1._wp
1416                    zopac = Nlevels-k+1                        ! z_opaque vertical index on vgrid_z
1417                 endif
1418                 if ( cldy(ip,ic,Nlevels-k) .eq. 1. ) then
1419                    lidarcldtype(ip,Nlevels-k ,1)    = lidarcldtype(ip,Nlevels-k ,1) + 1._wp
1420                    z_top = Nlevels-k    ! top cloud layer vertical index on vgrid_z
1421                 endif
1422              enddo
1423     ! Summing opaque cloud mean temperatures and altitudes
1424     ! as defined in Vaillant de Guelis et al. 2017a, AMT
1425              if (zopac .ne. 0) then
1426                 cldtypetemp(ip,1) = cldtypetemp(ip,1) + ( tmp(ip,zopac) + tmp(ip,z_top) )/2.
1427                 cldtypetemp(ip,3) = cldtypetemp(ip,3) + tmp(ip,zopac)                 ! z_opaque
1428                 cldtypemeanz(ip,1) = cldtypemeanz(ip,1) + ( vgrid_z(zopac) + vgrid_z(z_top) )/2.
1429                 cldtypemeanzse(ip,1) = cldtypemeanzse(ip,1) + (( vgrid_z(zopac) + vgrid_z(z_top) )/2.) - surfelev(ip)
1430                 cldtypemeanzse(ip,3) = cldtypemeanzse(ip,3) + ( vgrid_z(zopac) - surfelev(ip) )
1431              else
1432                 cldlay(ip,ic,1) = 0
1433              endif
1434           endif
1435
1436     ! Thin cloud profiles
1437           if ( cldlay(ip,ic,2) .eq. 1. ) then
1438              topcloud = 0._wp
1439              z_top = 0._wp
1440              z_base = 0._wp
1441              do k=1,Nlevels
1442     ! Declaring thin cloud fraction for 3D variable
1443     ! From TOA-2-SFC
1444                 if ( cldy(ip,ic,k) .eq. 1. .and. topcloud .eq. 1. ) then
1445                    lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1._wp
1446                    z_base = k ! bottom cloud layer
1447                 endif
1448                 if ( cldy(ip,ic,k) .eq. 1. .and. topcloud .eq. 0. ) then
1449                    lidarcldtype(ip,k,2) = lidarcldtype(ip,k,2) + 1._wp
1450                    z_top = k  ! top cloud layer
1451                    z_base = k ! bottom cloud layer
1452                    topcloud = 1._wp
1453                 endif
1454              enddo
1455     ! Computing mean emissivity using layers below the bottom cloud layer to the surface
1456              srmean = 0._wp
1457              srcount = 0._wp
1458              cloudemis = 0._wp
1459              do k=z_base+1,Nlevels
1460                 if (  (x(ip,ic,k) .gt. S_att_opaq) .and. (x(ip,ic,k) .lt. 1.0) .and. (x(ip,ic,k) .ne. undef)  ) then
1461                    srmean = srmean + x(ip,ic,k)
1462                    srcount = srcount + 1.
1463                 endif
1464              enddo
1465              ! If clear sky layers exist below bottom cloud layer
1466              if ( srcount .gt. 0. ) then
1467                 trans2 = srmean/srcount              ! thin cloud transmittance**2
1468                 tau_app = -(log(trans2))/2.          ! apparent cloud optical depth
1469                 tau_vis = tau_app/eta                ! cloud visible optical depth (multiple scat.)
1470                 tau_ir = tau_vis/2.                  ! approx. relation between visible and IR ODs
1471                 cloudemis = 1. - exp(-tau_ir)        ! no diffusion in IR considered : emis = 1-T
1472                 count_emis(ip) = count_emis(ip) + 1.
1473              endif
1474     ! Summing thin cloud mean temperatures and altitudes
1475     ! as defined in Vaillant de Guelis et al. 2017a, AMT
1476              cldtypetemp(ip,2) = cldtypetemp(ip,2) + ( tmp(ip,z_base) + tmp(ip,z_top) )/2.
1477              cldtypemeanz(ip,2) = cldtypemeanz(ip,2) + ( vgrid_z(z_base) + vgrid_z(z_top) )/2.
1478              cldtypemeanzse(ip,2) = cldtypemeanzse(ip,2) + (( vgrid_z(z_base) + vgrid_z(z_top) )/2.) - surfelev(ip)
1479              cldthinemis(ip) = cldthinemis(ip) + cloudemis
1480           endif
1481
1482       enddo
1483    enddo   
1484
1485    ! 3D cloud types fraction (opaque=1 and thin=2 clouds)
1486    where ( nsub(:,:) .gt. 0. )
1487       lidarcldtype(:,:,1) = lidarcldtype(:,:,1)/nsub(:,:)
1488       lidarcldtype(:,:,2) = lidarcldtype(:,:,2)/nsub(:,:)
1489    elsewhere
1490       lidarcldtype(:,:,1) = undef
1491       lidarcldtype(:,:,2) = undef
1492    endwhere
1493    ! 3D z_opaque fraction (=3)
1494    where ( nsubopaq(:,:) .gt. 0. )
1495       lidarcldtype(:,:,3) = lidarcldtype(:,:,3)/nsubopaq(:,:)
1496    elsewhere
1497       lidarcldtype(:,:,3) = undef
1498       lidarcldtype(:,:,4) = undef !declaring undef for opacity as well
1499    endwhere
1500    ! 3D opacity fraction (=4) !Summing z_opaque fraction from TOA(k=1) to SFC(k=Nlevels)
1501       lidarcldtype(:,1,4) = lidarcldtype(:,1,3) !top layer equal to 3D z_opaque fraction
1502    do ip = 1, Npoints
1503        do k = 2, Nlevels
1504            if ( (lidarcldtype(ip,k,3) .ne. undef) .and. (lidarcldtype(ip,k-1,4) .ne. undef) ) then
1505                lidarcldtype(ip,k,4) = lidarcldtype(ip,k,3) + lidarcldtype(ip,k-1,4)
1506            else
1507                lidarcldtype(ip,k,4) = undef
1508            endif
1509        enddo
1510    enddo
1511
1512    ! Layered cloud types (opaque, thin and z_opaque 2D variables)
1513
1514    do iz = 1, Ntype
1515       do ic = 1, Ncolumns
1516          cldtype(:,iz)  = cldtype(:,iz)  + cldlay(:,ic,iz)
1517          nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz)
1518       enddo
1519    enddo
1520
1521    ! Mean temperature and altitude
1522    where (cldtype(:,1) .gt. 0.)
1523       cldtypetemp(:,1) = cldtypetemp(:,1)/cldtype(:,1) ! opaque cloud temp
1524       cldtypetemp(:,3) = cldtypetemp(:,3)/cldtype(:,1) ! z_opaque
1525       cldtypemeanz(:,1) = cldtypemeanz(:,1)/cldtype(:,1) ! opaque cloud alt
1526       cldtypemeanzse(:,1) = cldtypemeanzse(:,1)/cldtype(:,1) ! opaque cloud alt - SE
1527       cldtypemeanzse(:,3) = cldtypemeanzse(:,3)/cldtype(:,1) ! z_opaque - SE
1528    elsewhere
1529       cldtypetemp(:,1) = undef
1530       cldtypetemp(:,3) = undef
1531       cldtypemeanz(:,1) = undef
1532       cldtypemeanzse(:,1) = undef
1533       cldtypemeanzse(:,3) = undef
1534    endwhere
1535
1536    where (cldtype(:,2) .gt. 0.) ! thin cloud
1537       cldtypetemp(:,2) = cldtypetemp(:,2)/cldtype(:,2)
1538       cldtypemeanz(:,2) = cldtypemeanz(:,2)/cldtype(:,2)
1539       cldtypemeanzse(:,2) = cldtypemeanzse(:,2)/cldtype(:,2)
1540    elsewhere
1541       cldtypetemp(:,2) = undef
1542       cldtypemeanz(:,2) = undef
1543       cldtypemeanzse(:,2) = undef
1544    endwhere
1545
1546    ! Mean thin cloud emissivity
1547    where (count_emis(:) .gt. 0.) ! thin cloud
1548       cldthinemis(:) = cldthinemis(:)/count_emis(:)
1549    elsewhere
1550       cldthinemis(:) = undef
1551    endwhere
1552
1553    where (nsublayer(:,:) .gt. 0.)
1554       cldtype(:,:) = cldtype(:,:)/nsublayer(:,:)
1555    elsewhere
1556       cldtype(:,:) = undef
1557    endwhere
1558
1559  END SUBROUTINE COSP_OPAQ
1560
1561end module mod_lidar_simulator
Note: See TracBrowser for help on using the repository browser.