source: LMDZ6/trunk/libf/phylmd/cosp/lidar_simulator.f90 @ 5422

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

.f90 <-> .F90 depending on cpp key use
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 29.4 KB
Line 
1! Copyright (c) 2009, Centre National de la Recherche Scientifique
2! All rights reserved.
3! $Revision: 88 $, $Date: 2013-11-13 15:08:38 +0100 (mer. 13 nov. 2013) $
4! $URL: http://cfmip-obs-sim.googlecode.com/svn/stable/v1.4.0/actsim/lidar_simulator.F90 $
5!
6! Redistribution and use in source and binary forms, with or without modification, are permitted
7! provided that the following conditions are met:
8!
9!     * Redistributions of source code must retain the above copyright notice, this list
10!       of conditions and the following disclaimer.
11!     * Redistributions in binary form must reproduce the above copyright notice, this list
12!       of conditions and the following disclaimer in the documentation and/or other materials
13!       provided with the distribution.
14!     * Neither the name of the LMD/IPSL/CNRS/UPMC nor the names of its
15!       contributors may be used to endorse or promote products derived from this software without
16!       specific prior written permission.
17!
18! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
19! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
20! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
21! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
22! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
23! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
24! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
25! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26     
27      SUBROUTINE lidar_simulator(npoints,nlev,npart,nrefl &
28                , undef &
29                , pres, presf, temp &
30                , q_lsliq, q_lsice, q_cvliq, q_cvice &
31                , ls_radliq, ls_radice, cv_radliq, cv_radice &
32                , ice_type, pmol, pnorm, pnorm_perp_tot,tautot, refl )
33!
34!---------------------------------------------------------------------------------
35! Purpose: To compute lidar signal from model-simulated profiles of cloud water
36!          and cloud fraction in each sub-column of each model gridbox.
37!
38! References:
39! Chepfer H., S. Bony, D. Winker, M. Chiriaco, J.-L. Dufresne, G. Seze (2008),
40! Use of CALIPSO lidar observations to evaluate the cloudiness simulated by a
41! climate model, Geophys. Res. Lett., 35, L15704, doi:10.1029/2008GL034207.
42!
43! Previous references:
44! Chiriaco et al, MWR, 2006; Chepfer et al., MWR, 2007
45!
46! Contacts: Helene Chepfer (chepfer@lmd.polytechnique.fr), Sandrine Bony (bony@lmd.jussieu.fr)
47!
48! May 2007: ActSim code of M. Chiriaco and H. Chepfer rewritten by S. Bony
49!
50! May 2008, H. Chepfer:
51! - Units of pressure inputs: Pa
52! - Non Spherical particles : LS Ice NS coefficients, CONV Ice NS coefficients
53! - New input: ice_type (0=ice-spheres ; 1=ice-non-spherical)
54!
55! June 2008, A. Bodas-Salcedo:
56! - Ported to Fortran 90 and optimisation changes
57!
58! August 2008, J-L Dufresne:
59! - Optimisation changes (sum instructions suppressed)
60!
61! October 2008, S. Bony,  H. Chepfer and J-L. Dufresne : 
62! - Interface with COSP v2.0:
63!      cloud fraction removed from inputs
64!      in-cloud condensed water now in input (instead of grid-averaged value)
65!      depolarisation diagnostic removed
66!      parasol (polder) reflectances (for 5 different solar zenith angles) added
67!
68! December 2008, S. Bony,  H. Chepfer and J-L. Dufresne :
69! - Modification of the integration of the lidar equation.
70! - change the cloud detection threshold
71!
72! April 2008, A. Bodas-Salcedo:
73! - Bug fix in computation of pmol and pnorm of upper layer
74!
75! April 2008, J-L. Dufresne
76! - Bug fix in computation of pmol and pnorm, thanks to Masaki Satoh: a factor 2
77! was missing. This affects the ATB values but not the cloud fraction.
78!
79! January 2013, G. Cesana and H. Chepfer:
80! - Add the perpendicular component of the backscattered signal (pnorm_perp_tot) in the arguments
81! - Add the temperature for each levels (temp) in the arguments
82! - Add the computation of the perpendicular component of the backscattered lidar signal
83! Reference: Cesana G. and H. Chepfer (2013): Evaluation of the cloud water phase
84! in a climate model using CALIPSO-GOCCP, J. Geophys. Res., doi: 10.1002/jgrd.50376
85!
86!---------------------------------------------------------------------------------
87!
88! Inputs:
89!  npoints  : number of horizontal points
90!  nlev : number of vertical levels
91!  npart: numberb of cloud meteors (stratiform_liq, stratiform_ice, conv_liq, conv_ice).
92!        Currently npart must be 4
93!  nrefl: number of solar zenith angles for parasol reflectances
94!  pres : pressure in the middle of atmospheric layers (full levels): Pa
95!  presf: pressure in the interface of atmospheric layers (half levels): Pa
96!     presf(..,1) : surface pressure ; presf(..,nlev+1)= TOA pressure
97!  temp : temperature of atmospheric layers: K
98!  q_lsliq: LS sub-column liquid water mixing ratio (kg/kg)
99!  q_lsice: LS sub-column ice water mixing ratio (kg/kg)
100!  q_cvliq: CONV sub-column liquid water mixing ratio (kg/kg)
101!  q_cvice: CONV sub-column ice water mixing ratio (kg/kg)
102!  ls_radliq: effective radius of LS liquid particles (meters)
103!  ls_radice: effective radius of LS ice particles (meters)
104!  cv_radliq: effective radius of CONV liquid particles (meters)
105!  cv_radice: effective radius of CONV ice particles (meters)
106!  ice_type : ice particle shape hypothesis (ice_type=0 for spheres, ice_type=1
107!             for non spherical particles)
108!
109! Outputs:
110!  pmol : molecular attenuated backscatter lidar signal power (m^-1.sr^-1)
111!  pnorm: total attenuated backscatter lidar signal power (m^-1.sr^-1)
112!  pnorm_perp_tot: perpendicular attenuated backscatter lidar signal power (m^-1.sr^-1)
113!  tautot: optical thickess integrated from top to level z
114!  refl : parasol(polder) reflectance
115!
116! Version 1.0 (June 2007)
117! Version 1.1 (May 2008)
118! Version 1.2 (June 2008)
119! Version 2.0 (October 2008)
120! Version 2.1 (December 2008)
121!---------------------------------------------------------------------------------
122
123      USE MOD_COSP_CONSTANTS, only : ok_debug_cosp
124      IMPLICIT NONE
125      REAL :: SRsat
126      PARAMETER (SRsat = 0.01) ! threshold full attenuation
127
128      LOGICAL ok_parasol
129      PARAMETER (ok_parasol=.true.)  ! set to .true. if you want to activate parasol reflectances
130
131      INTEGER i, k
132     
133      INTEGER INDX_LSLIQ,INDX_LSICE,INDX_CVLIQ,INDX_CVICE
134      PARAMETER (INDX_LSLIQ=1,INDX_LSICE=2,INDX_CVLIQ=3,INDX_CVICE=4)
135! inputs:
136      INTEGER npoints,nlev,npart,ice_type
137      INTEGER nrefl
138      real undef                 ! undefined value
139      REAL pres(npoints,nlev)    ! pressure full levels
140      REAL presf(npoints,nlev+1) ! pressure half levels
141      REAL q_lsliq(npoints,nlev), q_lsice(npoints,nlev)
142      REAL q_cvliq(npoints,nlev), q_cvice(npoints,nlev)
143      REAL ls_radliq(npoints,nlev), ls_radice(npoints,nlev)
144      REAL cv_radliq(npoints,nlev), cv_radice(npoints,nlev)
145
146! outputs (for each subcolumn):
147
148      REAL pmol(npoints,nlev)  ! molecular backscatter signal power (m^-1.sr^-1)
149      REAL pnorm(npoints,nlev) ! total lidar backscatter signal power (m^-1.sr^-1)
150      REAL tautot(npoints,nlev)! optical thickess integrated from top
151      REAL refl(npoints,nrefl)! parasol reflectance ! parasol
152
153! actsim variables:
154
155      REAL km, rdiffm, Qscat, Cmol
156      PARAMETER (Cmol = 6.2446e-32) ! depends on wavelength
157      PARAMETER (km = 1.38e-23)     ! Boltzmann constant (J/K)
158
159      PARAMETER (rdiffm = 0.7)      ! multiple scattering correction parameter
160      PARAMETER (Qscat = 2.0)       ! particle scattering efficiency at 532 nm
161
162      REAL rholiq, rhoice
163      PARAMETER (rholiq=1.0e+03)     ! liquid water (kg/m3)
164      PARAMETER (rhoice=0.5e+03)     ! ice (kg/m3)
165
166      REAL pi, rhopart(npart)
167      REAL polpart(npart,5)  ! polynomial coefficients derived for spherical and non spherical
168                             ! particules
169
170!   grid-box variables:
171      REAL rad_part(npoints,nlev,npart)
172      REAL rhoair(npoints,nlev), zheight(npoints,nlev+1)
173      REAL beta_mol(npoints,nlev), alpha_mol(npoints,nlev)
174      REAL kp_part(npoints,nlev,npart)
175
176!   sub-column variables:
177      REAL qpart(npoints,nlev,npart) ! mixing ratio particles in each subcolumn
178      REAL alpha_part(npoints,nlev,npart)
179      REAL tau_mol_lay(npoints)  ! temporary variable, moL. opt. thickness of layer k
180      REAL tau_mol(npoints,nlev) ! optical thickness between TOA and bottom of layer k
181      REAL tau_part(npoints,nlev,npart)
182      REAL betatot(npoints,nlev)
183      REAL tautot_lay(npoints)   ! temporary variable, total opt. thickness of layer k
184!     Optical thickness from TOA to surface for Parasol
185     REAL tautot_S_liq(npoints),tautot_S_ice(npoints)     ! for liq and ice clouds
186
187
188! Local variables
189      REAL Alpha, Beta, Gamma  ! Polynomial coefficient for ATBperp computation
190      REAL temp(npoints,nlev)                   ! temperature of layer k
191      REAL betatot_ice(npoints,nlev)    ! backscatter coefficient for ice particles
192      REAL beta_perp_ice(npoints,nlev)  ! perpendicular backscatter coefficient for ice
193      REAL betatot_liq(npoints,nlev)    ! backscatter coefficient for liquid particles
194      REAL beta_perp_liq(npoints,nlev)  ! perpendicular backscatter coefficient for liq
195      REAL tautot_ice(npoints,nlev)     ! total optical thickness of ice
196      REAL tautot_liq(npoints,nlev)     ! total optical thickness of liq
197      REAL tautot_lay_ice(npoints)    ! total optical thickness of ice in the layer k
198      REAL tautot_lay_liq(npoints)    ! total optical thickness of liq in the layer k
199      REAL pnorm_liq(npoints,nlev)    ! lidar backscattered signal power for liquid
200      REAL pnorm_ice(npoints,nlev)    ! lidar backscattered signal power for ice
201      REAL pnorm_perp_ice(npoints,nlev) ! perpendicular lidar backscattered signal power for ice
202      REAL pnorm_perp_liq(npoints,nlev) ! perpendicular lidar backscattered signal power for liq
203
204      REAL :: seuil
205
206! Output variable
207      REAL pnorm_perp_tot (npoints,nlev) ! perpendicular lidar backscattered signal power
208
209!------------------------------------------------------------
210!---- 0. Initialisation :
211!------------------------------------------------------------
212betatot_ice(:,:)=0
213betatot_liq(:,:)=0
214beta_perp_ice(:,:)=0
215beta_perp_liq(:,:)=0
216tautot_ice(:,:)=0
217tautot_liq(:,:)=0
218tautot_lay_ice(:)=0;
219tautot_lay_liq(:)=0;
220pnorm_liq(:,:)=0
221pnorm_ice(:,:)=0
222pnorm_perp_ice(:,:)=0
223pnorm_perp_liq(:,:)=0
224pnorm_perp_tot(:,:)=0
225
226
227! Polynomial coefficients (Alpha, Beta, Gamma) which allow to compute the ATBperpendicular
228! as a function of the ATB for ice or liquid cloud particles derived from CALIPSO-GOCCP
229! observations at 120m vertical grid (Cesana and Chepfer, JGR, 2013).
230!
231! Relationship between ATBice and ATBperp,ice for ice particles
232!  ATBperp,ice = Alpha*ATBice
233         Alpha = 0.2904
234
235! Relationship between ATBice and ATBperp,ice for liquid particles
236!  ATBperp,ice = Beta*ATBice^2 + Gamma*ATBice
237         Beta = 0.4099
238         Gamma = 0.009
239
240  if (ok_debug_cosp) then
241     seuil=1.e-18
242  else
243     seuil=0.0
244  endif
245!------------------------------------------------------------
246!---- 1. Preliminary definitions and calculations :
247!------------------------------------------------------------
248
249      if ( npart .ne. 4 ) then
250        print *,'Error in lidar_simulator, npart should be 4, not',npart
251        stop
252      endif
253
254      pi = dacos(-1.D0)
255
256! Polynomial coefficients for spherical liq/ice particles derived from Mie theory.
257! Polynomial coefficients for non spherical particles derived from a composite of
258! Ray-tracing theory for large particles (e.g. Noel et al., Appl. Opt., 2001)
259! and FDTD theory for very small particles (Yang et al., JQSRT, 2003).
260
261! We repeat the same coefficients for LS and CONV cloud to make code more readable
262!*     LS Liquid water coefficients:
263         polpart(INDX_LSLIQ,1) =  2.6980e-8
264         polpart(INDX_LSLIQ,2) = -3.7701e-6
265         polpart(INDX_LSLIQ,3) =  1.6594e-4
266         polpart(INDX_LSLIQ,4) = -0.0024
267         polpart(INDX_LSLIQ,5) =  0.0626
268!*     LS Ice coefficients:
269      if (ice_type.eq.0) then     
270         polpart(INDX_LSICE,1) = -1.0176e-8
271         polpart(INDX_LSICE,2) =  1.7615e-6
272         polpart(INDX_LSICE,3) = -1.0480e-4
273         polpart(INDX_LSICE,4) =  0.0019
274         polpart(INDX_LSICE,5) =  0.0460
275      endif
276!*     LS Ice NS coefficients:
277      if (ice_type.eq.1) then
278         polpart(INDX_LSICE,1) = 1.3615e-8
279         polpart(INDX_LSICE,2) = -2.04206e-6
280         polpart(INDX_LSICE,3) = 7.51799e-5
281         polpart(INDX_LSICE,4) = 0.00078213
282         polpart(INDX_LSICE,5) = 0.0182131
283      endif
284!*     CONV Liquid water coefficients:
285         polpart(INDX_CVLIQ,1) =  2.6980e-8
286         polpart(INDX_CVLIQ,2) = -3.7701e-6
287         polpart(INDX_CVLIQ,3) =  1.6594e-4
288         polpart(INDX_CVLIQ,4) = -0.0024
289         polpart(INDX_CVLIQ,5) =  0.0626
290!*     CONV Ice coefficients:
291      if (ice_type.eq.0) then
292         polpart(INDX_CVICE,1) = -1.0176e-8
293         polpart(INDX_CVICE,2) =  1.7615e-6
294         polpart(INDX_CVICE,3) = -1.0480e-4
295         polpart(INDX_CVICE,4) =  0.0019
296         polpart(INDX_CVICE,5) =  0.0460
297      endif
298      if (ice_type.eq.1) then
299         polpart(INDX_CVICE,1) = 1.3615e-8
300         polpart(INDX_CVICE,2) = -2.04206e-6
301         polpart(INDX_CVICE,3) = 7.51799e-5
302         polpart(INDX_CVICE,4) = 0.00078213
303         polpart(INDX_CVICE,5) = 0.0182131
304      endif
305
306! density:
307!*    clear-sky air:
308      rhoair = pres/(287.04*temp)
309
310!*    liquid/ice particules:
311      rhopart(INDX_LSLIQ) = rholiq
312      rhopart(INDX_LSICE) = rhoice
313      rhopart(INDX_CVLIQ) = rholiq
314      rhopart(INDX_CVICE) = rhoice
315
316! effective radius particles:
317      rad_part(:,:,INDX_LSLIQ) = ls_radliq(:,:)
318      rad_part(:,:,INDX_LSICE) = ls_radice(:,:)
319      rad_part(:,:,INDX_CVLIQ) = cv_radliq(:,:)
320      rad_part(:,:,INDX_CVICE) = cv_radice(:,:)
321      rad_part(:,:,:)=MAX(rad_part(:,:,:),0.)
322      rad_part(:,:,:)=MIN(rad_part(:,:,:),70.0e-6)
323     
324! altitude at half pressure levels:
325      zheight(:,1) = 0.0
326      do k = 2, nlev+1
327        zheight(:,k) = zheight(:,k-1) &
328                  -(presf(:,k)-presf(:,k-1))/(rhoair(:,k-1)*9.81)
329      enddo
330
331!------------------------------------------------------------
332!---- 2. Molecular alpha and beta:
333!------------------------------------------------------------
334
335      beta_mol = pres/km/temp * Cmol
336      alpha_mol = 8.0*pi/3.0 * beta_mol
337
338!------------------------------------------------------------
339!---- 3. Particles alpha and beta:
340!------------------------------------------------------------
341
342! polynomes kp_lidar derived from Mie theory:
343      do i = 1, npart
344       where ( rad_part(:,:,i).gt.0.0)
345         kp_part(:,:,i) = &
346            polpart(i,1)*(rad_part(:,:,i)*1e6)**4 &
347          + polpart(i,2)*(rad_part(:,:,i)*1e6)**3 &
348          + polpart(i,3)*(rad_part(:,:,i)*1e6)**2 &
349          + polpart(i,4)*(rad_part(:,:,i)*1e6) &
350          + polpart(i,5)
351        elsewhere
352         kp_part(:,:,i) = 0.
353        endwhere
354      enddo
355     
356! mixing ratio particules in each subcolumn:
357          qpart(:,:,INDX_LSLIQ) = q_lsliq(:,:) ! oct08
358          qpart(:,:,INDX_LSICE) = q_lsice(:,:) ! oct08
359          qpart(:,:,INDX_CVLIQ) = q_cvliq(:,:) ! oct08
360          qpart(:,:,INDX_CVICE) = q_cvice(:,:) ! oct08
361
362! alpha of particles in each subcolumn:
363      do i = 1, npart
364        where ( rad_part(:,:,i).gt.0.0)
365          alpha_part(:,:,i) = 3.0/4.0 * Qscat &
366                 * rhoair(:,:) * qpart(:,:,i) &
367                 / (rhopart(i) * rad_part(:,:,i) )
368        elsewhere
369          alpha_part(:,:,i) = 0.
370        endwhere
371      enddo
372
373!------------------------------------------------------------
374!---- 4.1 Total Backscatter signal:
375!------------------------------------------------------------
376
377! optical thickness (molecular):
378!     opt. thick of each layer
379      tau_mol(:,1:nlev) = alpha_mol(:,1:nlev) &
380         & *(zheight(:,2:nlev+1)-zheight(:,1:nlev))
381!     opt. thick from TOA
382      DO k = nlev-1, 1, -1
383        tau_mol(:,k) = tau_mol(:,k) + tau_mol(:,k+1)
384      ENDDO
385
386! optical thickness (particles):
387
388      tau_part = rdiffm * alpha_part
389      DO i = 1, npart
390!       opt. thick of each layer
391        tau_part(:,:,i) = tau_part(:,:,i) &
392           & * (zheight(:,2:nlev+1)-zheight(:,1:nlev) )
393!       opt. thick from TOA
394        DO k = nlev-1, 1, -1
395          tau_part(:,k,i) = tau_part(:,k,i) + tau_part(:,k+1,i)
396        ENDDO
397      ENDDO
398
399! molecular signal:
400!      Upper layer
401       pmol(:,nlev) = beta_mol(:,nlev) / (2.*tau_mol(:,nlev)) &
402            & * (1.-exp(-2.0*tau_mol(:,nlev)))
403!      Other layers
404       DO k= nlev-1, 1, -1
405        tau_mol_lay(:) = tau_mol(:,k)-tau_mol(:,k+1) ! opt. thick. of layer k
406        WHERE (tau_mol_lay(:).GT.0.)
407          pmol(:,k) = beta_mol(:,k) * EXP(-2.0*tau_mol(:,k+1)) / (2.*tau_mol_lay(:)) &
408            & * (1.-exp(-2.0*tau_mol_lay(:)))
409        ELSEWHERE
410!         This must never happend, but just in case, to avoid div. by 0
411          pmol(:,k) = beta_mol(:,k) * EXP(-2.0*tau_mol(:,k+1))
412        END WHERE
413      END DO
414
415! Total signal (molecular + particules):
416!
417!
418! For performance reason on vector computers, the 2 following lines should not be used
419! and should be replace by the later one.
420!      betatot(:,:) = beta_mol(:,:) + sum(kp_part*alpha_part,dim=3)
421!      tautot(:,:)  = tau_mol(:,:)  + sum(tau_part,dim=3)
422      betatot(:,:) = beta_mol(:,:)
423      tautot(:,:)  = tau_mol(:,:)
424      DO i = 1, npart
425           betatot(:,:) = betatot(:,:) + kp_part(:,:,i)*alpha_part(:,:,i)
426           tautot(:,:) = tautot(:,:)  + tau_part(:,:,i)
427      ENDDO ! i
428!
429!     Upper layer
430      pnorm(:,nlev) = betatot(:,nlev) / (2.*tautot(:,nlev)) &
431            & * (1.-exp(-2.0*tautot(:,nlev)))
432
433!     Other layers
434      DO k= nlev-1, 1, -1
435          tautot_lay(:) = tautot(:,k)-tautot(:,k+1) ! optical thickness of layer k
436        WHERE (tautot_lay(:).GT.0.)
437          pnorm(:,k) = betatot(:,k) * EXP(-2.0*tautot(:,k+1)) / (2.*tautot_lay(:)) &
438               & * (1.-EXP(-2.0*tautot_lay(:)))
439        ELSEWHERE
440!         This must never happend, but just in case, to avoid div. by 0
441          pnorm(:,k) = betatot(:,k) * EXP(-2.0*tautot(:,k+1))
442        END WHERE
443      END DO
444
445!------------------------------------------------------------
446!---- 4.2 Ice/Liq Backscatter signal:
447!------------------------------------------------------------
448
449! Contribution of the molecular to beta
450      betatot_ice(:,:) = beta_mol(:,:)
451      betatot_liq(:,:) = beta_mol(:,:)
452
453      tautot_ice(:,:) = tau_mol(:,:)
454      tautot_liq(:,:) = tau_mol(:,:)
455
456      DO i = 2, npart,2
457           betatot_ice(:,:) = betatot_ice(:,:)+ kp_part(:,:,i)*alpha_part(:,:,i)
458           tautot_ice(:,:) = tautot_ice(:,:)  + tau_part(:,:,i)
459      ENDDO ! i
460      DO i = 1, npart,2
461           betatot_liq(:,:) = betatot_liq(:,:)+ kp_part(:,:,i)*alpha_part(:,:,i)
462           tautot_liq(:,:) = tautot_liq(:,:)  + tau_part(:,:,i)
463      ENDDO ! i
464
465
466! Computation of the ice and liquid lidar backscattered signal (ATBice and ATBliq)
467!     Ice only
468!     Upper layer
469      pnorm_ice(:,nlev) = betatot_ice(:,nlev) / (2.*tautot_ice(:,nlev)) &
470            & * (1.-exp(-2.0*tautot_ice(:,nlev)))
471
472      DO k= nlev-1, 1, -1
473          tautot_lay_ice(:) = tautot_ice(:,k)-tautot_ice(:,k+1)
474        WHERE (tautot_lay_ice(:).GT.0.)
475         pnorm_ice(:,k)=betatot_ice(:,k)*EXP(-2.0*tautot_ice(:,k+1))/(2.*tautot_lay_ice(:)) &
476               & * (1.-EXP(-2.0*tautot_lay_ice(:)))
477        ELSEWHERE
478         pnorm_ice(:,k)=betatot_ice(:,k)*EXP(-2.0*tautot_ice(:,k+1))
479        END WHERE
480      ENDDO
481
482!     Liquid only
483!     Upper layer
484      pnorm_liq(:,nlev) = betatot_liq(:,nlev) / (2.*tautot_liq(:,nlev)) &
485            & * (1.-exp(-2.0*tautot_liq(:,nlev)))
486
487      DO k= nlev-1, 1, -1
488          tautot_lay_liq(:) = tautot_liq(:,k)-tautot_liq(:,k+1)
489        WHERE (tautot_lay_liq(:).GT.0.)
490          pnorm_liq(:,k)=betatot_liq(:,k)*EXP(-2.0*tautot_liq(:,k+1))/(2.*tautot_lay_liq(:)) &
491               & * (1.-EXP(-2.0*tautot_lay_liq(:)))
492        ELSEWHERE
493          pnorm_liq(:,k)=betatot_liq(:,k)*EXP(-2.0*tautot_liq(:,k+1))
494        END WHERE
495      ENDDO
496
497
498! Computation of ATBperp,ice/liq from ATBice/liq including the multiple scattering
499! contribution (Cesana and Chepfer 2013, JGR)
500!  ATBperp,ice = Alpha*ATBice
501!  ATBperp,liq = Beta*ATBliq^2 + Gamma*ATBliq
502
503      DO k= nlev, 1, -1
504         pnorm_perp_ice(:,k) = Alpha * pnorm_ice(:,k) ! Ice particles
505         pnorm_perp_liq(:,k) = 1000*Beta * pnorm_liq(:,k)**2 + Gamma * pnorm_liq(:,k) ! Liquid particles
506      ENDDO
507
508! Computation of beta_perp_ice/liq using the lidar equation
509!     Ice only
510!     Upper layer
511      beta_perp_ice(:,nlev) = pnorm_perp_ice(:,nlev) * (2.*tautot_ice(:,nlev)) &
512            & / (1.-exp(-2.0*tautot_ice(:,nlev)))
513
514      DO k= nlev-1, 1, -1
515        tautot_lay_ice(:) = tautot_ice(:,k)-tautot_ice(:,k+1)
516        WHERE (tautot_lay_ice(:).GT.0.)
517         beta_perp_ice(:,k) = pnorm_perp_ice(:,k)/ EXP(-2.0*tautot_ice(:,k+1)) * (2.*tautot_lay_ice(:)) &
518            & / (1.-exp(-2.0*tautot_lay_ice(:)))
519
520        ELSEWHERE
521         beta_perp_ice(:,k)=pnorm_perp_ice(:,k)/EXP(-2.0*tautot_ice(:,k+1))
522        END WHERE
523      ENDDO
524
525!     Liquid only
526!     Upper layer
527      beta_perp_liq(:,nlev) = pnorm_perp_liq(:,nlev) * (2.*tautot_liq(:,nlev)) &
528            & / (1.-exp(-2.0*tautot_liq(:,nlev)))
529
530      DO k= nlev-1, 1, -1
531          tautot_lay_liq(:) = tautot_liq(:,k)-tautot_liq(:,k+1)
532        WHERE (tautot_lay_liq(:).GT.0.)
533         beta_perp_liq(:,k) = pnorm_perp_liq(:,k)/ max(seuil,EXP(-2.0*tautot_liq(:,k+1))) &
534            & * (2.*tautot_lay_liq(:)) / (1.-exp(-2.0*tautot_lay_liq(:)))
535
536        ELSEWHERE
537         beta_perp_liq(:,k)=pnorm_perp_liq(:,k)/EXP(-2.0*tautot_liq(:,k+1))
538        END WHERE
539      ENDDO
540
541
542
543!------------------------------------------------------------
544!---- 4.3 Perpendicular Backscatter signal:
545!------------------------------------------------------------
546
547! Computation of the total perpendicular lidar signal (ATBperp for liq+ice)
548!     Upper layer
549    WHERE(tautot(:,nlev).GT.0)
550          pnorm_perp_tot(:,nlev) = &
551              (beta_perp_ice(:,nlev)+beta_perp_liq(:,nlev)-(beta_mol(:,nlev)/(1+1/0.0284))) / (2.*tautot(:,nlev)) &
552              & * (1.-exp(-2.0*tautot(:,nlev)))
553    ELSEWHERE
554    pnorm_perp_tot(:,nlev) = 0.
555    ENDWHERE
556
557!     Other layers
558      DO k= nlev-1, 1, -1
559          tautot_lay(:) = tautot(:,k)-tautot(:,k+1) ! optical thickness of layer k
560
561          ! The perpendicular component of the molecular backscattered signal (Betaperp) has been
562          ! taken into account two times (once for liquid and once for ice).
563          ! We remove one contribution using
564          ! Betaperp=beta_mol(:,k)/(1+1/0.0284)) [bodhaine et al. 1999] in the following equations:
565            WHERE (pnorm(:,k).eq.0)
566                  pnorm_perp_tot(:,k)=0.
567                  ELSEWHERE
568                    WHERE (tautot_lay(:).GT.0.)
569                      pnorm_perp_tot(:,k) = &
570                          (beta_perp_ice(:,k)+beta_perp_liq(:,k)-(beta_mol(:,k)/(1+1/0.0284))) * &
571                          EXP(-2.0*tautot(:,k+1)) / (2.*tautot_lay(:)) &
572                          & * (1.-EXP(-2.0*tautot_lay(:)))
573                    ELSEWHERE
574          !         This must never happen, but just in case, to avoid div. by 0
575                      pnorm_perp_tot(:,k) = &
576                           (beta_perp_ice(:,k)+beta_perp_liq(:,k)-(beta_mol(:,k)/(1+1/0.0284))) * &
577                          EXP(-2.0*tautot(:,k+1))
578                    END WHERE
579            ENDWHERE
580
581      END DO
582
583!-------- End computation Lidar --------------------------
584
585!---------------------------------------------------------
586!  Parasol/Polder module
587!
588!  Purpose : Compute reflectance for one particular viewing direction
589!  and 5 solar zenith angles (calculation valid only over ocean)
590! ---------------------------------------------------------
591
592! initialization:
593    refl(:,:) = 0.0
594
595! activate parasol calculations:
596    if (ok_parasol) then
597
598!     Optical thickness from TOA to surface
599      tautot_S_liq = 0.
600      tautot_S_ice = 0.
601      tautot_S_liq(:) = tautot_S_liq(:) &
602         + tau_part(:,1,1) + tau_part(:,1,3)
603      tautot_S_ice(:) = tautot_S_ice(:) &
604         + tau_part(:,1,2) + tau_part(:,1,4)
605
606      call parasol(npoints,nrefl,undef  &
607                 ,tautot_S_liq,tautot_S_ice &
608                 ,refl)
609
610    endif ! ok_parasol
611
612  END SUBROUTINE lidar_simulator
613!
614!---------------------------------------------------------------------------------
615!
616  SUBROUTINE parasol(npoints,nrefl,undef  &
617                       ,tautot_S_liq,tautot_S_ice  &
618                       ,refl)
619!---------------------------------------------------------------------------------
620! Purpose: To compute Parasol reflectance signal from model-simulated profiles
621!          of cloud water and cloud fraction in each sub-column of each model
622!          gridbox.
623!
624!
625! December 2008, S. Bony,  H. Chepfer and J-L. Dufresne :
626! - optimization for vectorization
627!
628! Version 2.0 (October 2008)
629! Version 2.1 (December 2008)
630!---------------------------------------------------------------------------------
631
632    IMPLICIT NONE
633
634! inputs
635    INTEGER npoints              ! Number of horizontal gridpoints
636    INTEGER nrefl                ! Number of angles for which the reflectance
637                                 ! is computed. Can not be greater then ntetas
638    REAL undef                   ! Undefined value. Currently not used
639    REAL tautot_S_liq(npoints)   ! liquid water cloud optical thickness,
640                                   ! integrated from TOA to surface
641    REAL tautot_S_ice(npoints)   ! same for ice water clouds only
642! outputs
643    REAL refl(npoints,nrefl)     ! Parasol reflectances
644!
645! Local variables
646    REAL tautot_S(npoints)       ! cloud optical thickness, from TOA to surface
647    REAL frac_taucol_liq(npoints), frac_taucol_ice(npoints)
648
649    REAL pi
650!   look up table variables:
651    INTEGER ny, it
652    INTEGER ntetas, nbtau        ! number of angle and of optical thickness
653                                   ! of the look-up table
654    PARAMETER (ntetas=5, nbtau=7)
655    REAL aa(ntetas,nbtau-1), ab(ntetas,nbtau-1)
656    REAL ba(ntetas,nbtau-1), bb(ntetas,nbtau-1) 
657    REAL tetas(ntetas),tau(nbtau)                       
658    REAL r_norm(ntetas)
659    REAL rlumA(ntetas,nbtau), rlumB(ntetas,nbtau)       
660    REAL rlumA_mod(npoints,5), rlumB_mod(npoints,5)
661
662    DATA tau   /0., 1., 5., 10., 20., 50., 100./
663    DATA tetas /0., 20., 40., 60., 80./
664   
665! Look-up table for spherical liquid particles:
666    DATA (rlumA(1,ny),ny=1,nbtau) /0.03, 0.090886, 0.283965, &
667     0.480587, 0.695235, 0.908229, 1.0 /
668    DATA (rlumA(2,ny),ny=1,nbtau) /0.03, 0.072185, 0.252596, &
669      0.436401,  0.631352, 0.823924, 0.909013 /
670    DATA (rlumA(3,ny),ny=1,nbtau) /0.03, 0.058410, 0.224707, &
671      0.367451,  0.509180, 0.648152, 0.709554 /
672    DATA (rlumA(4,ny),ny=1,nbtau) /0.03, 0.052498, 0.175844, &
673      0.252916,  0.326551, 0.398581, 0.430405 /
674    DATA (rlumA(5,ny),ny=1,nbtau) /0.03, 0.034730, 0.064488, &
675      0.081667,  0.098215, 0.114411, 0.121567 /
676
677! Look-up table for ice particles:
678    DATA (rlumB(1,ny),ny=1,nbtau) /0.03, 0.092170, 0.311941, &
679       0.511298, 0.712079 , 0.898243 , 0.976646 /
680    DATA (rlumB(2,ny),ny=1,nbtau) /0.03, 0.087082, 0.304293, &
681       0.490879,  0.673565, 0.842026, 0.912966 /
682    DATA (rlumB(3,ny),ny=1,nbtau) /0.03, 0.083325, 0.285193, &
683      0.430266,  0.563747, 0.685773,  0.737154 /
684    DATA (rlumB(4,ny),ny=1,nbtau) /0.03, 0.084935, 0.233450, &
685      0.312280, 0.382376, 0.446371, 0.473317 /
686    DATA (rlumB(5,ny),ny=1,nbtau) /0.03, 0.054157, 0.089911, &
687      0.107854, 0.124127, 0.139004, 0.145269 /
688
689!--------------------------------------------------------------------------------
690! Lum_norm=f(tetaS,tau_cloud) derived from adding-doubling calculations
691!        valid ONLY ABOVE OCEAN (albedo_sfce=5%)
692!        valid only in one viewing direction (theta_v=30�, phi_s-phi_v=320�)
693!        based on adding-doubling radiative transfer computation
694!        for tau values (0 to 100) and for tetas values (0 to 80)
695!        for 2 scattering phase functions: liquid spherical, ice non spherical
696
697    IF ( nrefl.GT. ntetas ) THEN
698        PRINT *,'Error in lidar_simulator, nrefl should be less then ',ntetas,' not',nrefl
699        STOP
700    ENDIF
701
702    rlumA_mod=0
703    rlumB_mod=0
704!
705    pi = ACOS(-1.0)
706    r_norm(:)=1./ COS(pi/180.*tetas(:))
707!
708    tautot_S_liq(:)=MAX(tautot_S_liq(:),tau(1))
709    tautot_S_ice(:)=MAX(tautot_S_ice(:),tau(1))
710    tautot_S(:) = tautot_S_ice(:) + tautot_S_liq(:)
711!
712! relative fraction of the opt. thick due to liquid or ice clouds
713    WHERE (tautot_S(:) .GT. 0.)
714        frac_taucol_liq(:) = tautot_S_liq(:) / tautot_S(:)
715        frac_taucol_ice(:) = tautot_S_ice(:) / tautot_S(:)
716    ELSEWHERE
717        frac_taucol_liq(:) = 1.
718        frac_taucol_ice(:) = 0.
719    END WHERE
720    tautot_S(:)=MIN(tautot_S(:),tau(nbtau))
721!
722! Linear interpolation :
723
724    DO ny=1,nbtau-1
725! microphysics A (liquid clouds)
726      aA(:,ny) = (rlumA(:,ny+1)-rlumA(:,ny))/(tau(ny+1)-tau(ny))
727      bA(:,ny) = rlumA(:,ny) - aA(:,ny)*tau(ny)
728! microphysics B (ice clouds)
729      aB(:,ny) = (rlumB(:,ny+1)-rlumB(:,ny))/(tau(ny+1)-tau(ny))
730      bB(:,ny) = rlumB(:,ny) - aB(:,ny)*tau(ny)
731    ENDDO
732!
733    DO it=1,ntetas
734      DO ny=1,nbtau-1
735        WHERE (tautot_S(:).GE.tau(ny).AND.tautot_S(:).LE.tau(ny+1))
736            rlumA_mod(:,it) = aA(it,ny)*tautot_S(:) + bA(it,ny)
737            rlumB_mod(:,it) = aB(it,ny)*tautot_S(:) + bB(it,ny)
738        END WHERE
739      END DO
740    END DO
741!
742    DO it=1,ntetas
743      refl(:,it) = frac_taucol_liq(:) * rlumA_mod(:,it) &
744         + frac_taucol_ice(:) * rlumB_mod(:,it)
745! normalized radiance -> reflectance:
746      refl(:,it) = refl(:,it) * r_norm(it)
747    ENDDO
748
749    RETURN
750  END SUBROUTINE parasol
Note: See TracBrowser for help on using the repository browser.