source: LMDZ6/trunk/libf/phylmd/cosp2/lidar_simulator.F90 @ 4021

Last change on this file since 4021 was 3358, checked in by idelkadi, 6 years ago

Implementation de la nouvelle version COSPv2 dans LMDZ.
Pour compiler avec makelmdz_fcma utiliser l'option "-cosp2 true"
File size: 47.2 KB
RevLine 
[3358]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! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
70module mod_lidar_simulator
71  USE COSP_KINDS,         ONLY: wp
72  USE MOD_COSP_CONFIG,    ONLY: SR_BINS,S_CLD,S_ATT,S_CLD_ATT,R_UNDEF,calipso_histBsct,  &
73                                use_vgrid,vgrid_zl,vgrid_zu
74  USE MOD_COSP_STATS,     ONLY: COSP_CHANGE_VERTICAL_GRID,hist1d
75  implicit none
76 
77  ! Polynomial coefficients (Alpha, Beta, Gamma) which allow to compute the
78  ! ATBperpendicular as a function of the ATB for ice or liquid cloud particles
79  ! derived from CALIPSO-GOCCP observations at 120m vertical grid
80  ! (Cesana and Chepfer, JGR, 2013).
81  !
82  ! Relationship between ATBice and ATBperp,ice for ice particles:
83  !                ATBperp,ice = Alpha*ATBice
84  ! Relationship between ATBice and ATBperp,ice for liquid particles:
85  !          ATBperp,ice = Beta*ATBice^2 + Gamma*ATBice
86  real(wp) :: &
87       alpha,beta,gamma   
88
89contains
90  ! ######################################################################################
91  ! SUBROUTINE lidar_subcolumn
92  ! Inputs with a vertical dimensions (nlev) should ordered in along the vertical
93  ! dimension from TOA-2-SFC, for example: varIN(nlev) is varIN @ SFC.
94  ! ######################################################################################
95  subroutine lidar_subcolumn(npoints,ncolumns,nlev,beta_mol,tau_mol,betatot,tautot,    &
96                         betatot_ice,tautot_ice,betatot_liq,tautot_liq,                  &
97                         pmol,pnorm,pnorm_perp_tot)
98
99    ! INPUTS
100    INTEGER,intent(in) :: &
101         npoints,      & ! Number of gridpoints
102         ncolumns,     & ! Number of subcolumns
103         nlev            ! Number of levels
104    REAL(WP),intent(in),dimension(npoints,nlev) :: &
105         beta_mol,     & ! Molecular backscatter coefficient
106         tau_mol         ! Molecular optical depth
107
108    REAL(WP),intent(in),dimension(npoints,ncolumns,nlev)       :: &
109         betatot,      & !
110         tautot,       & ! Optical thickess integrated from top
111         betatot_ice,  & ! Backscatter coefficient for ice particles
112         betatot_liq,  & ! Backscatter coefficient for liquid particles
113         tautot_ice,   & ! Total optical thickness of ice
114         tautot_liq      ! Total optical thickness of liq
115
116    ! OUTPUTS
117    REAL(WP),intent(out),dimension(npoints,nlev) :: &
118         pmol            ! Molecular attenuated backscatter lidar signal power(m^-1.sr^-1)
119    REAL(WP),intent(out),dimension(npoints,ncolumns,nlev) :: &
120         pnorm,        & ! Molecular backscatter signal power (m^-1.sr^-1)
121         pnorm_perp_tot  ! Perpendicular lidar backscattered signal power
122
123    ! LOCAL VARIABLES
124    INTEGER :: k,icol
125    REAL(WP),dimension(npoints) :: &
126         tautot_lay        !
127    REAL(WP),dimension(npoints,ncolumns,nlev) :: &
128         pnorm_liq,      & ! Lidar backscattered signal power for liquid
129         pnorm_ice,      & ! Lidar backscattered signal power for ice
130         pnorm_perp_ice, & ! Perpendicular lidar backscattered signal power for ice
131         pnorm_perp_liq, & ! Perpendicular lidar backscattered signal power for liq
132         beta_perp_ice,  & ! Perpendicular backscatter coefficient for ice
133         beta_perp_liq     ! Perpendicular backscatter coefficient for liquid   
134
135    ! ####################################################################################
136    ! *) Molecular signal
137    ! ####################################################################################
138    call cmp_backsignal(nlev,npoints,beta_mol(1:npoints,1:nlev),&
139                        tau_mol(1:npoints,1:nlev),pmol(1:npoints,1:nlev))
140                       
141    ! ####################################################################################
142    ! PLANE PARRALLEL FIELDS
143    ! ####################################################################################
144    do icol=1,ncolumns
145       ! #################################################################################
146       ! *) Total Backscatter signal
147       ! #################################################################################
148       call cmp_backsignal(nlev,npoints,betatot(1:npoints,icol,1:nlev),&
149            tautot(1:npoints,icol,1:nlev),pnorm(1:npoints,icol,1:nlev))
150       ! #################################################################################
151       ! *) Ice/Liq Backscatter signal
152       ! #################################################################################
153       ! Computation of the ice and liquid lidar backscattered signal (ATBice and ATBliq)
154       ! Ice only
155       call cmp_backsignal(nlev,npoints,betatot_ice(1:npoints,icol,1:nlev),&
156                      tautot_ice(1:npoints,icol,1:nlev),&
157                      pnorm_ice(1:npoints,icol,1:nlev))
158       ! Liquid only
159       call cmp_backsignal(nlev,npoints,betatot_liq(1:npoints,icol,1:nlev),&
160                      tautot_liq(1:npoints,icol,1:nlev),&
161                      pnorm_liq(1:npoints,icol,1:nlev))
162    enddo
163
164    ! ####################################################################################
165    ! PERDENDICULAR FIELDS
166    ! ####################################################################################
167    do icol=1,ncolumns
168
169       ! #################################################################################
170       ! *) Ice/Liq Perpendicular Backscatter signal
171       ! #################################################################################
172       ! Computation of ATBperp,ice/liq from ATBice/liq including the multiple scattering
173       ! contribution (Cesana and Chepfer 2013, JGR)
174       do k=1,nlev
175          ! Ice particles
176          pnorm_perp_ice(1:npoints,icol,k) = Alpha * pnorm_ice(1:npoints,icol,k)
177
178          ! Liquid particles
179          pnorm_perp_liq(1:npoints,icol,k) = 1000._wp*Beta*pnorm_liq(1:npoints,icol,k)**2+&
180               Gamma*pnorm_liq(1:npoints,icol,k)
181       enddo
182 
183       ! #################################################################################
184       ! *) Computation of beta_perp_ice/liq using the lidar equation
185       ! #################################################################################
186       ! Ice only
187       call cmp_beta(nlev,npoints,pnorm_perp_ice(1:npoints,icol,1:nlev),&
188              tautot_ice(1:npoints,icol,1:nlev),beta_perp_ice(1:npoints,icol,1:nlev))       
189 
190       ! Liquid only
191       call cmp_beta(nlev,npoints,pnorm_perp_liq(1:npoints,icol,1:nlev),&
192            tautot_liq(1:npoints,icol,1:nlev),beta_perp_liq(1:npoints,icol,1:nlev))
193         
194       ! #################################################################################
195       ! *) Perpendicular Backscatter signal
196       ! #################################################################################
197       ! Computation of the total perpendicular lidar signal (ATBperp for liq+ice)
198       ! Upper layer
199       WHERE(tautot(1:npoints,icol,1) .gt. 0)
200          pnorm_perp_tot(1:npoints,icol,1) = (beta_perp_ice(1:npoints,icol,1)+           &
201               beta_perp_liq(1:npoints,icol,1)-                                          &
202               (beta_mol(1:npoints,1)/(1._wp+1._wp/0.0284_wp))) /                        &
203               (2._wp*tautot(1:npoints,icol,1))*                                         &
204               (1._wp-exp(-2._wp*tautot(1:npoints,icol,1)))
205       ELSEWHERE
206          pnorm_perp_tot(1:npoints,icol,1) = 0._wp
207       ENDWHERE                                                 
208             
209       ! Other layers
210       do k=2,nlev
211          ! Optical thickness of layer k
212          tautot_lay(1:npoints) = tautot(1:npoints,icol,k)-tautot(1:npoints,icol,k-1)
213
214          ! The perpendicular component of the molecular backscattered signal (Betaperp)
215          ! has been taken into account two times (once for liquid and once for ice).
216          ! We remove one contribution using
217          ! Betaperp=beta_mol(:,k)/(1+1/0.0284)) [bodhaine et al. 1999] in the following
218          ! equations:
219          WHERE (pnorm(1:npoints,icol,k) .eq. 0)
220             pnorm_perp_tot(1:npoints,icol,k)=0._wp
221          ELSEWHERE
222             where(tautot_lay(1:npoints) .gt. 0.)
223                pnorm_perp_tot(1:npoints,icol,k) = (beta_perp_ice(1:npoints,icol,k)+     &
224                   beta_perp_liq(1:npoints,icol,k)-(beta_mol(1:npoints,k)/(1._wp+1._wp/  &
225                   0.0284_wp)))*EXP(-2._wp*tautot(1:npoints,icol,k-1))/                  &
226                   (2._wp*tautot_lay(1:npoints))* (1._wp-EXP(-2._wp*tautot_lay(1:npoints)))
227             elsewhere
228                pnorm_perp_tot(1:npoints,icol,k) = (beta_perp_ice(1:npoints,icol,k)+     &
229                   beta_perp_liq(1:npoints,icol,k)-(beta_mol(1:npoints,k)/(1._wp+1._wp/  &
230                   0.0284_wp)))*EXP(-2._wp*tautot(1:npoints,icol,k-1))
231             endwhere
232          ENDWHERE
233       END DO
234    enddo
235  end subroutine lidar_subcolumn
236
237  ! ######################################################################################
238  ! SUBROUTINE lidar_column
239  ! ######################################################################################
240  subroutine lidar_column(npoints,ncol,nlevels,llm,max_bin,tmp, pnorm,                   &
241                           pnorm_perp, pmol, pplay, ok_lidar_cfad, ncat, cfad2,    &
242                           lidarcld, lidarcldphase, cldlayer, zlev, zlev_half,           &
243                           cldlayerphase, lidarcldtmp)
244    integer,parameter :: &
245         nphase = 6 ! Number of cloud layer phase types
246
247    ! Inputs
248    integer,intent(in) :: &
249         npoints, & ! Number of horizontal grid points
250         ncol,    & ! Number of subcolumns
251         nlevels, & ! Number of vertical layers (OLD grid)
252         llm,     & ! Number of vertical layers (NEW grid)
253         max_bin, & ! Number of bins for SR CFADs
254         ncat       ! Number of cloud layer types (low,mid,high,total)
255    real(wp),intent(in),dimension(npoints,ncol,Nlevels) :: &
256         pnorm,   & ! Lidar ATB
257         pnorm_perp ! Lidar perpendicular ATB
258    real(wp),intent(in),dimension(npoints,Nlevels) :: &
259         pmol,    & ! Molecular ATB
260         pplay,   & ! Pressure on model levels (Pa)
261         tmp        ! Temperature at each levels
262    logical,intent(in) :: &
263         ok_lidar_cfad ! True if lidar CFAD diagnostics need to be computed
264    real(wp),intent(in),dimension(npoints,nlevels) :: &
265         zlev        ! Model full levels
266    real(wp),intent(in),dimension(npoints,nlevels+1) :: &
267         zlev_half   ! Model half levels
268         
269    ! Outputs
270    real(wp),intent(inout),dimension(npoints,llm) :: &
271         lidarcld      ! 3D "lidar" cloud fraction
272    real(wp),intent(inout),dimension(npoints,ncat) :: &
273         cldlayer      ! "lidar" cloud layer fraction (low, mid, high, total)
274    real(wp),intent(inout),dimension(npoints,llm,nphase) :: &
275         lidarcldphase ! 3D "lidar" phase cloud fraction
276    real(wp),intent(inout),dimension(npoints,40,5) :: &
277         lidarcldtmp   ! 3D "lidar" phase cloud fraction as a function of temp
278    real(wp),intent(inout),dimension(npoints,ncat,nphase) :: &
279         cldlayerphase ! "lidar" phase low mid high cloud fraction
280    real(wp),intent(inout),dimension(npoints,max_bin,llm) :: &
281         cfad2         ! CFADs of SR
282
283    ! Local Variables
284    integer :: ic,i,j
285    real(wp),dimension(npoints,ncol,llm) :: &
286         x3d
287    real(wp),dimension(npoints,llm) :: &
288         x3d_c,pnorm_c
289    real(wp)  :: &
290         xmax
291    real(wp),dimension(npoints,1,Nlevels) :: t_in,ph_in,betamol_in
292    real(wp),dimension(npoints,ncol,llm)  :: pnormFlip,pnorm_perpFlip
293    real(wp),dimension(npoints,1,llm)     :: tmpFlip,pplayFlip,betamolFlip
294
295    ! Vertically regrid input data
296    if (use_vgrid) then
297       t_in(:,1,:)=tmp(:,nlevels:1:-1)
298       call cosp_change_vertical_grid(Npoints,1,Nlevels,zlev(:,nlevels:1:-1),zlev_half(:,nlevels:1:-1),&
299            t_in,llm,vgrid_zl(llm:1:-1),vgrid_zu(llm:1:-1),tmpFlip(:,1,llm:1:-1))
300       ph_in(:,1,:) = pplay(:,nlevels:1:-1)
301       call cosp_change_vertical_grid(Npoints,1,Nlevels,zlev(:,nlevels:1:-1),zlev_half(:,nlevels:1:-1),&
302            ph_in,llm,vgrid_zl(llm:1:-1),vgrid_zu(llm:1:-1),pplayFlip(:,1,llm:1:-1))
303       betamol_in(:,1,:) = pmol(:,nlevels:1:-1)
304       call cosp_change_vertical_grid(Npoints,1,Nlevels,zlev(:,nlevels:1:-1),zlev_half(:,nlevels:1:-1),&
305            betamol_in,llm,vgrid_zl(llm:1:-1),vgrid_zu(llm:1:-1),betamolFlip(:,1,llm:1:-1))
306       call cosp_change_vertical_grid(Npoints,Ncol,Nlevels,zlev(:,nlevels:1:-1),zlev_half(:,nlevels:1:-1),&
307            pnorm(:,:,nlevels:1:-1),llm,vgrid_zl(llm:1:-1),vgrid_zu(llm:1:-1),pnormFlip(:,:,llm:1:-1))
308       call cosp_change_vertical_grid(Npoints,Ncol,Nlevels,zlev(:,nlevels:1:-1),zlev_half(:,nlevels:1:-1),&
309            pnorm_perp(:,:,nlevels:1:-1),llm,vgrid_zl(llm:1:-1),vgrid_zu(llm:1:-1),pnorm_perpFlip(:,:,llm:1:-1))
310    endif
311
312    ! Initialization (The histogram bins, are set up during initialization and the
313    ! maximum value is used as the upper bounds.)
314    xmax = maxval(calipso_histBsct)
315
316    ! Compute LIDAR scattering ratio
317    if (use_vgrid) then
318       do ic = 1, ncol
319          pnorm_c = pnormFlip(:,ic,:)
320          where ((pnorm_c .lt. xmax) .and. (betamolFlip(:,1,:) .lt. xmax) .and.          &
321                (betamolFlip(:,1,:) .gt. 0.0 ))
322             x3d_c = pnorm_c/betamolFlip(:,1,:)
323          elsewhere
324             x3d_c = R_UNDEF
325          end where
326          x3d(:,ic,:) = x3d_c
327       enddo
328       ! Diagnose cloud fractions for subcolumn lidar scattering ratios
329       CALL COSP_CLDFRAC(npoints,ncol,llm,ncat,nphase,tmpFlip,x3d,pnormFlip,             &
330                         pnorm_perpFlip,pplayFlip,S_att,S_cld,S_cld_att,R_UNDEF,         &
331                         lidarcld,cldlayer,lidarcldphase,cldlayerphase,lidarcldtmp)                           
332    else
333       do ic = 1, ncol
334          pnorm_c = pnorm(:,ic,:)
335          where ((pnorm_c.lt.xmax) .and. (pmol.lt.xmax) .and. (pmol.gt. 0.0 ))
336             x3d_c = pnorm_c/pmol
337          elsewhere
338             x3d_c = R_UNDEF
339          end where
340          x3d(:,ic,:) = x3d_c
341       enddo
342       ! Diagnose cloud fractions for subcolumn lidar scattering ratios
343       CALL COSP_CLDFRAC(npoints,ncol,nlevels,ncat,nphase,tmp,x3d,pnorm,pnorm_perp,pplay,&
344                         S_att,S_cld,S_cld_att,R_UNDEF,lidarcld,cldlayer,lidarcldphase,  &
345                         cldlayerphase,lidarcldtmp)
346    endif
347
348    ! CFADs
349    if (ok_lidar_cfad) then
350       ! CFADs of subgrid-scale lidar scattering ratios
351       do i=1,Npoints
352          do j=1,llm
353             cfad2(i,:,j) = hist1D(ncol,x3d(i,:,j),SR_BINS,calipso_histBsct)
354          enddo
355       enddo
356       where(cfad2 .ne. R_UNDEF) cfad2=cfad2/ncol
357
358    endif
359   
360    ! Unit conversions
361    where(lidarcld /= R_UNDEF)      lidarcld      = lidarcld*100._wp
362    where(cldlayer /= R_UNDEF)      cldlayer      = cldlayer*100._wp
363    where(cldlayerphase /= R_UNDEF) cldlayerphase = cldlayerphase*100._wp
364    where(lidarcldphase /= R_UNDEF) lidarcldphase = lidarcldphase*100._wp
365    where(lidarcldtmp /= R_UNDEF)   lidarcldtmp   = lidarcldtmp*100._wp
366   
367  end subroutine lidar_column
368
369  ! ######################################################################################
370  ! The subroutines below compute the attenuated backscatter signal and the lidar
371  ! backscatter coefficients using eq (1) from doi:0094-8276/08/2008GL034207
372  ! ######################################################################################
373  subroutine cmp_backsignal(nlev,npoints,beta,tau,pnorm)
374    ! INPUTS
375    integer, intent(in) :: nlev,npoints
376    real(wp),intent(in),dimension(npoints,nlev) :: beta,tau
377
378    ! OUTPUTS
379    real(wp),intent(out),dimension(npoints,nlev) :: pnorm
380
381    ! Internal Variables
382    real(wp), dimension(npoints) :: tautot_lay
383    integer :: k
384
385    ! Uppermost layer
386    pnorm(:,1) = beta(:,1) / (2._wp*tau(:,1)) * (1._wp-exp(-2._wp*tau(:,1)))
387
388    ! Other layers
389    do k=2,nlev
390       tautot_lay(:) = tau(:,k)-tau(:,k-1)
391       WHERE ( EXP(-2._wp*tau(:,k-1)) .gt. 0. )
392          WHERE (tautot_lay(:) .gt. 0.)
393             pnorm(:,k) = beta(:,k)*EXP(-2._wp*tau(:,k-1)) /&
394                  (2._wp*tautot_lay(:))*(1._wp-EXP(-2._wp*tautot_lay(:)))
395          ELSEWHERE
396             ! This must never happen, but just in case, to avoid div. by 0
397             pnorm(:,k) = beta(:,k) * EXP(-2._wp*tau(:,k-1))
398          END WHERE
399       ELSEWHERE
400          pnorm(:,k) = 0._wp!beta(:,k)
401       END WHERE
402    END DO
403  end subroutine cmp_backsignal
404
405  subroutine cmp_beta(nlev,npoints,pnorm,tau,beta)
406    ! INPUTS
407    integer, intent(in) :: nlev,npoints
408    real(wp),intent(in),dimension(npoints,nlev) :: pnorm,tau
409
410    ! OUTPUTS
411    real(wp),intent(out),dimension(npoints,nlev) :: beta
412
413    ! Internal Variables
414    real(wp), dimension(npoints) :: tautot_lay
415    integer :: k
416
417    beta(:,1) = pnorm(:,1) * (2._wp*tau(:,1))/(1._wp-exp(-2._wp*tau(:,1)))
418    do k=2,nlev
419       tautot_lay(:) = tau(:,k)-tau(:,k-1)       
420       WHERE ( EXP(-2._wp*tau(:,k-1)) .gt. 0. )
421          WHERE (tautot_lay(:) .gt. 0.)
422             beta(:,k) = pnorm(:,k)/ EXP(-2._wp*tau(:,k-1))* &
423                  (2._wp*tautot_lay(:))/(1._wp-exp(-2._wp*tautot_lay(:)))
424          ELSEWHERE
425             beta(:,k)=pnorm(:,k)/EXP(-2._wp*tau(:,k-1))
426          END WHERE
427       ELSEWHERE
428          beta(:,k)=pnorm(:,k)
429       END WHERE
430    ENDDO
431
432  end subroutine cmp_beta
433    ! ####################################################################################
434    ! SUBROUTINE cosp_cldfrac
435    ! Conventions: Ncat must be equal to 4
436    ! ####################################################################################
437    SUBROUTINE COSP_CLDFRAC(Npoints,Ncolumns,Nlevels,Ncat,Nphase,tmp,x,ATB,ATBperp,      &
438                               pplay,S_att,S_cld,S_cld_att,undef,lidarcld,cldlayer,      &
439                               lidarcldphase,cldlayerphase,lidarcldtemp)
440    ! Parameters
441    integer,parameter :: Ntemp=40 ! indice of the temperature vector
442    real(wp),parameter,dimension(Ntemp+1) :: &
443       tempmod = [0.0,   183.15,186.15,189.15,192.15,195.15,198.15,201.15,204.15,207.15, &
444                  210.15,213.15,216.15,219.15,222.15,225.15,228.15,231.15,234.15,237.15, &
445                  240.15,243.15,246.15,249.15,252.15,255.15,258.15,261.15,264.15,267.15, &
446                  270.15,273.15,276.15,279.15,282.15,285.15,288.15,291.15,294.15,297.15, &
447                  473.15]
448         
449    ! Polynomial coefficient of the phase discrimination line used to separate liquid from ice
450    ! (Cesana and Chepfer, JGR, 2013)
451    ! ATBperp = ATB^5*alpha50 + ATB^4*beta50 + ATB^3*gamma50 + ATB^2*delta50 + ATB*epsilon50 + zeta50
452    real(wp),parameter :: &
453       alpha50   = 9.0322e+15_wp,  & !
454       beta50    = -2.1358e+12_wp, & !
455       gamma50   = 173.3963e06_wp, & !
456       delta50   = -3.9514e03_wp,  & !
457       epsilon50 = 0.2559_wp,      & !
458       zeta50    = -9.4776e-07_wp    !
459       
460        ! Inputs
461    integer,intent(in) :: &
462       Npoints,  & ! Number of gridpoints
463       Ncolumns, & ! Number of subcolumns
464       Nlevels,  & ! Number of vertical levels
465       Ncat,     & ! Number of cloud layer types
466       Nphase      ! Number of cloud layer phase types
467                       ! [ice,liquid,undefined,false ice,false liquid,Percent of ice]
468    real(wp),intent(in) :: &
469       S_att,    & !
470       S_cld,    & !
471       S_cld_att,& ! New threshold for undefine cloud phase detection
472       undef       ! Undefined value
473    real(wp),intent(in),dimension(Npoints,Ncolumns,Nlevels) :: &
474       x,        & !
475       ATB,      & ! 3D attenuated backscatter
476       ATBperp     ! 3D attenuated backscatter (perpendicular)
477    real(wp),intent(in),dimension(Npoints,Nlevels) :: &
478       tmp,      & ! Temperature   
479       pplay       ! Pressure
480
481        ! Outputs
482    real(wp),intent(out),dimension(Npoints,Ntemp,5) :: &
483       lidarcldtemp  ! 3D Temperature 1=tot,2=ice,3=liq,4=undef,5=ice/ice+liq
484    real(wp),intent(out),dimension(Npoints,Nlevels,Nphase) :: &
485       lidarcldphase ! 3D cloud phase fraction
486    real(wp),intent(out),dimension(Npoints,Nlevels) :: &
487       lidarcld      ! 3D cloud fraction
488    real(wp),intent(out),dimension(Npoints,Ncat) :: &
489       cldlayer      ! Low, middle, high, total cloud fractions
490    real(wp),intent(out),dimension(Npoints,Ncat,Nphase) :: &
491       cldlayerphase ! Low, middle, high, total cloud fractions for ice liquid and undefine phase   
492   
493    ! Local variables
494    integer  :: &
495       ip, k, iz, ic, ncol, nlev, i, itemp, toplvlsat
496    real(wp) :: &
497       p1,checktemp, ATBperp_tmp,checkcldlayerphase, checkcldlayerphase2
498    real(wp),dimension(Npoints,Nlevels) :: &
499       nsub,lidarcldphasetmp   
500    real(wp),dimension(Npoints,Ntemp) :: &
501       sumlidarcldtemp,lidarcldtempind
502    real(wp),dimension(Npoints,Ncolumns,Ncat) :: &
503       cldlay,nsublay   
504    real(wp),dimension(Npoints,Ncat) :: &
505       nsublayer,cldlayerphasetmp,cldlayerphasesum
506    real(wp),dimension(Npoints,Ncolumns,Nlevels) :: &   
507       tmpi, & ! Temperature of ice cld
508       tmpl, & ! Temperature of liquid cld
509       tmpu, & ! Temperature of undef cld
510       cldy, & !
511       srok    !
512    real(wp),dimension(Npoints,Ncolumns,Ncat,Nphase) :: &
513       cldlayphase ! subgrided low mid high phase cloud fraction
514             
515    ! ####################################################################################
516        ! 1) Initialize   
517    ! ####################################################################################
518    lidarcld              = 0._wp
519    nsub                  = 0._wp
520    cldlay                = 0._wp
521    nsublay               = 0._wp
522    ATBperp_tmp           = 0._wp
523    lidarcldphase(:,:,:)  = 0._wp
524    cldlayphase(:,:,:,:)  = 0._wp
525    cldlayerphase(:,:,:)  = 0._wp
526    tmpi(:,:,:)           = 0._wp
527    tmpl(:,:,:)           = 0._wp
528    tmpu(:,:,:)           = 0._wp
529    cldlayerphasesum(:,:) = 0._wp
530    lidarcldtemp(:,:,:)   = 0._wp
531    lidarcldtempind(:,:)  = 0._wp
532    sumlidarcldtemp(:,:)  = 0._wp
533    lidarcldphasetmp(:,:) = 0._wp
534    toplvlsat             = 0
535
536    ! ####################################################################################
537    ! 2) Cloud detection
538    ! ####################################################################################
539    do k=1,Nlevels
540       ! Cloud detection at subgrid-scale:
541       where ((x(:,:,k) .gt. S_cld) .and. (x(:,:,k) .ne. undef) )
542          cldy(:,:,k)=1._wp
543       elsewhere
544          cldy(:,:,k)=0._wp
545       endwhere
546       
547       ! Number of usefull sub-columns:
548       where ((x(:,:,k) .gt. S_att) .and. (x(:,:,k) .ne. undef) )
549          srok(:,:,k)=1._wp
550       elsewhere
551          srok(:,:,k)=0._wp
552       endwhere
553    enddo   
554   
555    ! ####################################################################################
556    ! 3) Grid-box 3D cloud fraction and layered cloud fractions(ISCCP pressure categories)
557    ! ####################################################################################
558    lidarcld = 0._wp
559    nsub     = 0._wp
560    cldlay   = 0._wp
561    nsublay  = 0._wp
562    do k=1,Nlevels
563       do ic = 1, Ncolumns
564          do ip = 1, Npoints
565         
566             ! Computation of the cloud fraction as a function of the temperature instead
567             ! of height, for ice,liquid and all clouds
568             if(srok(ip,ic,k).gt.0.)then
569                do itemp=1,Ntemp
570                   if( (tmp(ip,k).ge.tempmod(itemp)).and.(tmp(ip,k).lt.tempmod(itemp+1)) )then
571                      lidarcldtempind(ip,itemp)=lidarcldtempind(ip,itemp)+1._wp
572                   endif
573                enddo
574             endif
575             
576             if(cldy(ip,ic,k).eq.1.)then
577                do itemp=1,Ntemp
578                   if( (tmp(ip,k) .ge. tempmod(itemp)).and.(tmp(ip,k) .lt. tempmod(itemp+1)) )then
579                      lidarcldtemp(ip,itemp,1)=lidarcldtemp(ip,itemp,1)+1._wp
580                   endif
581                enddo
582             endif
583
584             iz=1
585             p1 = pplay(ip,k)
586             if ( p1.gt.0. .and. p1.lt.(440._wp*100._wp)) then ! high clouds
587                iz=3
588             else if(p1.ge.(440._wp*100._wp) .and. p1.lt.(680._wp*100._wp)) then ! mid clouds
589                iz=2
590             endif
591             
592             cldlay(ip,ic,iz) = MAX(cldlay(ip,ic,iz),cldy(ip,ic,k))
593             cldlay(ip,ic,4)  = MAX(cldlay(ip,ic,4),cldy(ip,ic,k))
594             lidarcld(ip,k)   = lidarcld(ip,k) + cldy(ip,ic,k)
595             
596             nsublay(ip,ic,iz) = MAX(nsublay(ip,ic,iz),srok(ip,ic,k))
597             nsublay(ip,ic,4)  = MAX(nsublay(ip,ic,4),srok(ip,ic,k))
598             nsub(ip,k)        = nsub(ip,k) + srok(ip,ic,k)
599             
600          enddo
601       enddo
602    enddo   
603   
604    ! Grid-box 3D cloud fraction
605    where ( nsub(:,:).gt.0.0 )
606       lidarcld(:,:) = lidarcld(:,:)/nsub(:,:)
607    elsewhere
608       lidarcld(:,:) = undef
609    endwhere
610   
611    ! Layered cloud fractions
612    cldlayer  = 0._wp
613    nsublayer = 0._wp
614    do iz = 1, Ncat
615       do ic = 1, Ncolumns
616          cldlayer(:,iz)  = cldlayer(:,iz)  + cldlay(:,ic,iz)
617          nsublayer(:,iz) = nsublayer(:,iz) + nsublay(:,ic,iz)
618       enddo
619    enddo
620    where (nsublayer(:,:) .gt. 0.0)
621       cldlayer(:,:) = cldlayer(:,:)/nsublayer(:,:)
622    elsewhere
623       cldlayer(:,:) = undef
624    endwhere
625             
626    ! ####################################################################################
627    ! 4) Grid-box 3D cloud Phase
628    ! ####################################################################################
629   
630    ! ####################################################################################
631    ! 4.1) For Cloudy pixels with 8.16km < z < 19.2km
632    ! ####################################################################################
633    do ncol=1,Ncolumns
634       do i=1,Npoints         
635          do nlev=1,23 ! from 19.2km until 8.16km
636               p1 = pplay(1,nlev)
637
638             ! Avoid zero values
639             if( (cldy(i,ncol,nlev).eq.1.) .and. (ATBperp(i,ncol,nlev).gt.0.) )then
640                ! Computation of the ATBperp along the phase discrimination line
641                ATBperp_tmp = (ATB(i,ncol,nlev)**5)*alpha50 + (ATB(i,ncol,nlev)**4)*beta50 + &
642                     (ATB(i,ncol,nlev)**3)*gamma50 + (ATB(i,ncol,nlev)**2)*delta50 + &
643                     ATB(i,ncol,nlev)*epsilon50 + zeta50     
644                ! ########################################################################
645                ! 4.1.a) Ice: ATBperp above the phase discrimination line
646                ! ########################################################################
647                if((ATBperp(i,ncol,nlev)-ATBperp_tmp) .ge. 0.)then ! Ice clouds
648
649                   ! ICE with temperature above 273,15°K = Liquid (false ice)
650                   if(tmp(i,nlev) .gt. 273.15) then ! Temperature above 273,15 K
651                     ! Liquid: False ice corrected by the temperature to Liquid
652                      lidarcldphase(i,nlev,2) = lidarcldphase(i,nlev,2)+1._wp ! False ice detection ==> added to Liquid
653                                   
654                      tmpl(i,ncol,nlev)       = tmp(i,nlev)
655                      lidarcldphase(i,nlev,5) = lidarcldphase(i,nlev,5)+1._wp ! Keep the information "temperature criterium used"                     
656                                                                              ! to classify the phase cloud
657                      cldlayphase(i,ncol,4,2) = 1. ! tot cloud
658                      if (p1 .gt. 0. .and. p1.lt.(440._wp*100._wp)) then ! high cloud
659                         cldlayphase(i,ncol,3,2) = 1._wp
660                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then ! mid cloud
661                         cldlayphase(i,ncol,2,2) = 1._wp
662                      else ! low cloud
663                         cldlayphase(i,ncol,1,2) = 1._wp
664                      endif
665                      cldlayphase(i,ncol,4,5) = 1._wp ! tot cloud
666                      ! High cloud
667                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
668                         cldlayphase(i,ncol,3,5) = 1._wp
669                      ! Middle cloud
670                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
671                         cldlayphase(i,ncol,2,5) = 1._wp
672                      ! Low cloud
673                      else
674                         cldlayphase(i,ncol,1,5) = 1._wp
675                      endif
676                   else
677                      ! ICE with temperature below 273,15°K
678                      lidarcldphase(i,nlev,1) = lidarcldphase(i,nlev,1)+1._wp
679                      tmpi(i,ncol,nlev)       = tmp(i,nlev)
680                      cldlayphase(i,ncol,4,1) = 1._wp ! tot cloud
681                      ! High cloud
682                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
683                         cldlayphase(i,ncol,3,1) = 1._wp
684                      ! Middle cloud   
685                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
686                         cldlayphase(i,ncol,2,1) = 1._wp
687                      ! Low cloud
688                      else
689                         cldlayphase(i,ncol,1,1) = 1._wp
690                      endif
691                   endif
692                ! ########################################################################
693                ! 4.1.b) Liquid: ATBperp below the phase discrimination line
694                ! ########################################################################
695                else
696                   ! Liquid with temperature above 231,15°K
697                   if(tmp(i,nlev) .gt. 231.15_wp) then
698                      lidarcldphase(i,nlev,2) = lidarcldphase(i,nlev,2)+1._wp
699                      tmpl(i,ncol,nlev)       = tmp(i,nlev)
700                      cldlayphase(i,ncol,4,2) = 1._wp ! tot cloud
701                      ! High cloud
702                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
703                         cldlayphase(i,ncol,3,2) = 1._wp
704                      ! Middle cloud   
705                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
706                         cldlayphase(i,ncol,2,2) = 1._wp
707                      ! Low cloud   
708                      else
709                         cldlayphase(i,ncol,1,2) = 1._wp
710                      endif
711                   else
712                      ! Liquid with temperature below 231,15°K = Ice (false liquid)
713                      tmpi(i,ncol,nlev)       = tmp(i,nlev)
714                      lidarcldphase(i,nlev,1) = lidarcldphase(i,nlev,1)+1._wp ! false liquid detection ==> added to ice
715                      lidarcldphase(i,nlev,4) = lidarcldphase(i,nlev,4)+1._wp
716                      cldlayphase(i,ncol,4,4) = 1._wp ! tot cloud
717                      ! High cloud
718                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
719                         cldlayphase(i,ncol,3,4) = 1._wp
720                      ! Middle cloud   
721                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
722                         cldlayphase(i,ncol,2,4) = 1._wp
723                      ! Low cloud
724                      else
725                         cldlayphase(i,ncol,1,4) = 1._wp
726                      endif
727                      cldlayphase(i,ncol,4,1) = 1._wp ! tot cloud
728                      ! High cloud
729                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
730                         cldlayphase(i,ncol,3,1) = 1._wp
731                      ! Middle cloud   
732                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
733                         cldlayphase(i,ncol,2,1) = 1._wp
734                      ! Low cloud   
735                      else
736                         cldlayphase(i,ncol,1,1) = 1._wp
737                      endif
738                   endif
739                endif ! end of discrimination condition
740             endif ! end of cloud condition
741          enddo ! end of altitude loop
742
743          ! ##############################################################################
744          ! 4.2) For Cloudy pixels with 0km < z < 8.16km
745          ! ##############################################################################
746          toplvlsat = 0
747          do nlev=24,Nlevels! from 8.16km until 0km
748             p1 = pplay(i,nlev)
749
750             if((cldy(i,ncol,nlev) .eq. 1.) .and. (ATBperp(i,ncol,nlev) .gt. 0.) )then
751                ! Computation of the ATBperp of 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.2.a) Ice: ATBperp above the phase discrimination line
757                ! ########################################################################
758                ! ICE with temperature above 273,15°K = Liquid (false ice)
759                if((ATBperp(i,ncol,nlev)-ATBperp_tmp) .ge. 0.)then ! Ice clouds
760                   if(tmp(i,nlev) .gt. 273.15)then
761                      lidarcldphase(i,nlev,2) = lidarcldphase(i,nlev,2)+1._wp ! false ice ==> liq
762                      tmpl(i,ncol,nlev)       = tmp(i,nlev)
763                      lidarcldphase(i,nlev,5) = lidarcldphase(i,nlev,5)+1._wp
764                      cldlayphase(i,ncol,4,2) = 1._wp ! tot cloud
765                      ! High cloud
766                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
767                         cldlayphase(i,ncol,3,2) = 1._wp
768                      ! Middle cloud   
769                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
770                         cldlayphase(i,ncol,2,2) = 1._wp
771                      ! Low cloud
772                      else
773                         cldlayphase(i,ncol,1,2) = 1._wp
774                      endif
775                     
776                      cldlayphase(i,ncol,4,5) = 1. ! 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                ! ########################################################################
805                ! 4.2.b) Liquid: ATBperp below the phase discrimination line
806                ! ########################################################################
807                else
808                   ! Liquid with temperature above 231,15°K
809                   if(tmp(i,nlev) .gt. 231.15)then
810                      lidarcldphase(i,nlev,2) = lidarcldphase(i,nlev,2)+1._wp
811                      tmpl(i,ncol,nlev)       = tmp(i,nlev)
812                      cldlayphase(i,ncol,4,2) = 1._wp ! tot cloud
813                      ! High cloud
814                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
815                         cldlayphase(i,ncol,3,2) = 1._wp
816                      ! Middle cloud   
817                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
818                         cldlayphase(i,ncol,2,2) = 1._wp
819                      ! Low cloud   
820                      else
821                         cldlayphase(i,ncol,1,2) = 1._wp
822                      endif
823                   else
824                      ! Liquid with temperature below 231,15°K = Ice (false liquid)
825                      tmpi(i,ncol,nlev)       = tmp(i,nlev)
826                      lidarcldphase(i,nlev,1) = lidarcldphase(i,nlev,1)+1._wp ! false liq ==> ice
827                      lidarcldphase(i,nlev,4) = lidarcldphase(i,nlev,4)+1._wp ! false liq ==> ice
828                      cldlayphase(i,ncol,4,4) = 1._wp ! tot cloud
829                      ! High cloud
830                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
831                         cldlayphase(i,ncol,3,4) = 1._wp
832                      ! Middle   
833                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
834                         cldlayphase(i,ncol,2,4) = 1._wp
835                      ! Low cloud   
836                      else
837                         cldlayphase(i,ncol,1,4) = 1._wp
838                      endif
839                     
840                      cldlayphase(i,ncol,4,1) = 1._wp ! tot cloud
841                      ! High cloud
842                      if (p1 .gt. 0. .and. p1 .lt. (440._wp*100._wp)) then
843                         cldlayphase(i,ncol,3,1) = 1._wp
844                      ! Middle cloud   
845                      else if(p1 .ge. (440._wp*100._wp) .and. p1 .lt. (680._wp*100._wp)) then
846                         cldlayphase(i,ncol,2,1) = 1._wp
847                      ! Low cloud   
848                      else
849                         cldlayphase(i,ncol,1,1) = 1._wp
850                      endif
851                   endif
852                endif ! end of discrimination condition
853               
854                toplvlsat=0
855               
856                ! Find the level of the highest cloud with SR>30
857                if(x(i,ncol,nlev) .gt. S_cld_att) then ! SR > 30.
858                    toplvlsat = nlev+1
859                    goto 99
860                endif
861             endif ! end of cloud condition
862          enddo ! end of altitude loop
86399        continue
864         
865          ! ##############################################################################
866          ! Undefined phase: For a cloud located below another cloud with SR>30
867          ! see Cesana and Chepfer 2013 Sect.III.2
868          ! ##############################################################################
869          if(toplvlsat.ne.0) then
870             do nlev = toplvlsat,Nlevels
871                p1 = pplay(i,nlev)
872                if(cldy(i,ncol,nlev).eq.1.)then
873                   lidarcldphase(i,nlev,3) = lidarcldphase(i,nlev,3)+1._wp
874                   tmpu(i,ncol,nlev)       = tmp(i,nlev)
875                   cldlayphase(i,ncol,4,3) = 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,3) = 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,3) = 1._wp
882                   ! Low cloud   
883                   else
884                      cldlayphase(i,ncol,1,3) = 1._wp
885                   endif
886                endif
887             enddo
888          endif
889          toplvlsat=0
890       enddo
891    enddo
892     
893    ! ####################################################################################
894    ! Computation of final cloud phase diagnosis
895    ! ####################################################################################
896
897    ! Compute the Ice percentage in cloud = ice/(ice+liq) as a function of the occurrences
898    lidarcldphasetmp(:,:) = lidarcldphase(:,:,1)+lidarcldphase(:,:,2);
899    WHERE (lidarcldphasetmp(:,:) .gt. 0.)
900       lidarcldphase(:,:,6)=lidarcldphase(:,:,1)/lidarcldphasetmp(:,:)
901    ELSEWHERE
902       lidarcldphase(:,:,6) = undef
903    ENDWHERE
904   
905    ! Compute Phase 3D Cloud Fraction
906    !WHERE (nsub(:,Nlevels:1:-1) .gt. 0.0 )
907    WHERE (nsub(:,:) .gt. 0.0 ) 
908       lidarcldphase(:,:,1)=lidarcldphase(:,:,1)/nsub(:,:)
909       lidarcldphase(:,:,2)=lidarcldphase(:,:,2)/nsub(:,:)
910       lidarcldphase(:,:,3)=lidarcldphase(:,:,3)/nsub(:,:)
911       lidarcldphase(:,:,4)=lidarcldphase(:,:,4)/nsub(:,:)
912       lidarcldphase(:,:,5)=lidarcldphase(:,:,5)/nsub(:,:)
913    ELSEWHERE
914       lidarcldphase(:,:,1) = undef
915       lidarcldphase(:,:,2) = undef
916       lidarcldphase(:,:,3) = undef
917       lidarcldphase(:,:,4) = undef
918       lidarcldphase(:,:,5) = undef
919    ENDWHERE
920
921    ! Compute Phase low mid high cloud fractions
922    do iz = 1, Ncat
923       do i=1,Nphase-3
924          do ic = 1, Ncolumns
925             cldlayerphase(:,iz,i)  = cldlayerphase(:,iz,i)  + cldlayphase(:,ic,iz,i)
926             cldlayerphasesum(:,iz) = cldlayerphasesum(:,iz) + cldlayphase(:,ic,iz,i)
927          enddo
928       enddo
929    enddo
930    do iz = 1, Ncat
931       do i=4,5
932          do ic = 1, Ncolumns
933             cldlayerphase(:,iz,i) = cldlayerphase(:,iz,i) + cldlayphase(:,ic,iz,i)
934          enddo
935       enddo
936    enddo
937   
938    ! Compute the Ice percentage in cloud = ice/(ice+liq)
939    cldlayerphasetmp(:,:)=cldlayerphase(:,:,1)+cldlayerphase(:,:,2)
940    WHERE (cldlayerphasetmp(:,:).gt. 0.)
941       cldlayerphase(:,:,6)=cldlayerphase(:,:,1)/cldlayerphasetmp(:,:)
942    ELSEWHERE
943       cldlayerphase(:,:,6) = undef
944    ENDWHERE
945   
946    do i=1,Nphase-1
947       WHERE ( cldlayerphasesum(:,:).gt.0.0 )
948          cldlayerphase(:,:,i) = (cldlayerphase(:,:,i)/cldlayerphasesum(:,:)) * cldlayer(:,:)
949       ENDWHERE
950    enddo
951   
952    do i=1,Npoints
953       do iz=1,Ncat
954          checkcldlayerphase=0.
955          checkcldlayerphase2=0.
956          if (cldlayerphasesum(i,iz) .gt. 0.0 )then
957             do ic=1,Nphase-3
958                checkcldlayerphase = checkcldlayerphase+cldlayerphase(i,iz,ic)
959             enddo
960             checkcldlayerphase2 = cldlayer(i,iz)-checkcldlayerphase
961             if((checkcldlayerphase2 .gt. 0.01) .or. (checkcldlayerphase2 .lt. -0.01) ) print *, checkcldlayerphase,cldlayer(i,iz)
962          endif
963       enddo
964    enddo
965   
966    do i=1,Nphase-1
967       WHERE (nsublayer(:,:) .eq. 0.0)
968          cldlayerphase(:,:,i) = undef
969       ENDWHERE
970    enddo
971 
972    ! Compute Phase 3D as a function of temperature
973    do nlev=1,Nlevels
974       do ncol=1,Ncolumns
975          do i=1,Npoints
976             do itemp=1,Ntemp
977                if(tmpi(i,ncol,nlev).gt.0.)then
978                   if((tmpi(i,ncol,nlev) .ge. tempmod(itemp)) .and. (tmpi(i,ncol,nlev) .lt. tempmod(itemp+1)) )then
979                      lidarcldtemp(i,itemp,2)=lidarcldtemp(i,itemp,2)+1._wp
980                   endif
981                elseif(tmpl(i,ncol,nlev) .gt. 0.)then
982                   if((tmpl(i,ncol,nlev) .ge. tempmod(itemp)) .and. (tmpl(i,ncol,nlev) .lt. tempmod(itemp+1)) )then
983                      lidarcldtemp(i,itemp,3)=lidarcldtemp(i,itemp,3)+1._wp
984                   endif
985                elseif(tmpu(i,ncol,nlev) .gt. 0.)then
986                   if((tmpu(i,ncol,nlev) .ge. tempmod(itemp)) .and. (tmpu(i,ncol,nlev) .lt. tempmod(itemp+1)) )then
987                      lidarcldtemp(i,itemp,4)=lidarcldtemp(i,itemp,4)+1._wp
988                   endif
989                endif
990             enddo
991          enddo
992       enddo
993    enddo
994   
995    ! Check temperature cloud fraction
996    do i=1,Npoints
997       do itemp=1,Ntemp
998          checktemp=lidarcldtemp(i,itemp,2)+lidarcldtemp(i,itemp,3)+lidarcldtemp(i,itemp,4)
999          !if(checktemp .NE. lidarcldtemp(i,itemp,1))then
1000          !   print *, i,itemp
1001          !   print *, lidarcldtemp(i,itemp,1:4)
1002          !endif
1003         
1004       enddo
1005    enddo
1006   
1007    ! Compute the Ice percentage in cloud = ice/(ice+liq)
1008    sumlidarcldtemp(:,:)=lidarcldtemp(:,:,2)+lidarcldtemp(:,:,3)   
1009    WHERE(sumlidarcldtemp(:,:) .gt. 0.)
1010       lidarcldtemp(:,:,5)=lidarcldtemp(:,:,2)/sumlidarcldtemp(:,:)
1011    ELSEWHERE
1012       lidarcldtemp(:,:,5)=undef
1013    ENDWHERE
1014   
1015    do i=1,4
1016       WHERE(lidarcldtempind(:,:) .gt. 0.)
1017          lidarcldtemp(:,:,i) = lidarcldtemp(:,:,i)/lidarcldtempind(:,:)
1018       ELSEWHERE
1019          lidarcldtemp(:,:,i) = undef
1020       ENDWHERE
1021    enddo
1022   
1023    RETURN
1024  END SUBROUTINE COSP_CLDFRAC
1025
1026end module mod_lidar_simulator
Note: See TracBrowser for help on using the repository browser.