source: LMDZ5/trunk/libf/phylmd/cosp/lidar_simulator.F90 @ 5314

Last change on this file since 5314 was 2428, checked in by idelkadi, 9 years ago

Mise a jour du simulateur COSP (passage de la version v3.2 a la version v1.4) :

  • mise a jour des sources pour ISCCP, CALIPSO et PARASOL
  • prise en compte des changements de phases pour les nuages (Calipso)
  • rajout de plusieurs diagnostiques (fraction nuageuse en fonction de la temperature, ...)

http://lmdz.lmd.jussieu.fr/Members/aidelkadi/cosp
  • 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.3 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      IMPLICIT NONE
124      REAL :: SRsat
125      PARAMETER (SRsat = 0.01) ! threshold full attenuation
126
127      LOGICAL ok_parasol
128      PARAMETER (ok_parasol=.true.)  ! set to .true. if you want to activate parasol reflectances
129
130      INTEGER i, k
131     
132      INTEGER INDX_LSLIQ,INDX_LSICE,INDX_CVLIQ,INDX_CVICE
133      PARAMETER (INDX_LSLIQ=1,INDX_LSICE=2,INDX_CVLIQ=3,INDX_CVICE=4)
134! inputs:
135      INTEGER npoints,nlev,npart,ice_type
136      INTEGER nrefl
137      real undef                 ! undefined value
138      REAL pres(npoints,nlev)    ! pressure full levels
139      REAL presf(npoints,nlev+1) ! pressure half levels
140      REAL q_lsliq(npoints,nlev), q_lsice(npoints,nlev)
141      REAL q_cvliq(npoints,nlev), q_cvice(npoints,nlev)
142      REAL ls_radliq(npoints,nlev), ls_radice(npoints,nlev)
143      REAL cv_radliq(npoints,nlev), cv_radice(npoints,nlev)
144
145! outputs (for each subcolumn):
146
147      REAL pmol(npoints,nlev)  ! molecular backscatter signal power (m^-1.sr^-1)
148      REAL pnorm(npoints,nlev) ! total lidar backscatter signal power (m^-1.sr^-1)
149      REAL tautot(npoints,nlev)! optical thickess integrated from top
150      REAL refl(npoints,nrefl)! parasol reflectance ! parasol
151
152! actsim variables:
153
154      REAL km, rdiffm, Qscat, Cmol
155      PARAMETER (Cmol = 6.2446e-32) ! depends on wavelength
156      PARAMETER (km = 1.38e-23)     ! Boltzmann constant (J/K)
157
158      PARAMETER (rdiffm = 0.7)      ! multiple scattering correction parameter
159      PARAMETER (Qscat = 2.0)       ! particle scattering efficiency at 532 nm
160
161      REAL rholiq, rhoice
162      PARAMETER (rholiq=1.0e+03)     ! liquid water (kg/m3)
163      PARAMETER (rhoice=0.5e+03)     ! ice (kg/m3)
164
165      REAL pi, rhopart(npart)
166      REAL polpart(npart,5)  ! polynomial coefficients derived for spherical and non spherical
167                             ! particules
168
169!   grid-box variables:
170      REAL rad_part(npoints,nlev,npart)
171      REAL rhoair(npoints,nlev), zheight(npoints,nlev+1)
172      REAL beta_mol(npoints,nlev), alpha_mol(npoints,nlev)
173      REAL kp_part(npoints,nlev,npart)
174
175!   sub-column variables:
176      REAL qpart(npoints,nlev,npart) ! mixing ratio particles in each subcolumn
177      REAL alpha_part(npoints,nlev,npart)
178      REAL tau_mol_lay(npoints)  ! temporary variable, moL. opt. thickness of layer k
179      REAL tau_mol(npoints,nlev) ! optical thickness between TOA and bottom of layer k
180      REAL tau_part(npoints,nlev,npart)
181      REAL betatot(npoints,nlev)
182      REAL tautot_lay(npoints)   ! temporary variable, total opt. thickness of layer k
183!     Optical thickness from TOA to surface for Parasol
184     REAL tautot_S_liq(npoints),tautot_S_ice(npoints)     ! for liq and ice clouds
185
186
187! Local variables
188      REAL Alpha, Beta, Gamma  ! Polynomial coefficient for ATBperp computation
189      REAL temp(npoints,nlev)                   ! temperature of layer k
190      REAL betatot_ice(npoints,nlev)    ! backscatter coefficient for ice particles
191      REAL beta_perp_ice(npoints,nlev)  ! perpendicular backscatter coefficient for ice
192      REAL betatot_liq(npoints,nlev)    ! backscatter coefficient for liquid particles
193      REAL beta_perp_liq(npoints,nlev)  ! perpendicular backscatter coefficient for liq
194      REAL tautot_ice(npoints,nlev)     ! total optical thickness of ice
195      REAL tautot_liq(npoints,nlev)     ! total optical thickness of liq
196      REAL tautot_lay_ice(npoints)    ! total optical thickness of ice in the layer k
197      REAL tautot_lay_liq(npoints)    ! total optical thickness of liq in the layer k
198      REAL pnorm_liq(npoints,nlev)    ! lidar backscattered signal power for liquid
199      REAL pnorm_ice(npoints,nlev)    ! lidar backscattered signal power for ice
200      REAL pnorm_perp_ice(npoints,nlev) ! perpendicular lidar backscattered signal power for ice
201      REAL pnorm_perp_liq(npoints,nlev) ! perpendicular lidar backscattered signal power for liq
202
203! Output variable
204      REAL pnorm_perp_tot (npoints,nlev) ! perpendicular lidar backscattered signal power
205
206!------------------------------------------------------------
207!---- 0. Initialisation :
208!------------------------------------------------------------
209betatot_ice(:,:)=0
210betatot_liq(:,:)=0
211beta_perp_ice(:,:)=0
212beta_perp_liq(:,:)=0
213tautot_ice(:,:)=0
214tautot_liq(:,:)=0
215tautot_lay_ice(:)=0;
216tautot_lay_liq(:)=0;
217pnorm_liq(:,:)=0
218pnorm_ice(:,:)=0
219pnorm_perp_ice(:,:)=0
220pnorm_perp_liq(:,:)=0
221pnorm_perp_tot(:,:)=0
222
223
224! Polynomial coefficients (Alpha, Beta, Gamma) which allow to compute the ATBperpendicular
225! as a function of the ATB for ice or liquid cloud particles derived from CALIPSO-GOCCP
226! observations at 120m vertical grid (Cesana and Chepfer, JGR, 2013).
227!
228! Relationship between ATBice and ATBperp,ice for ice particles
229!  ATBperp,ice = Alpha*ATBice
230         Alpha = 0.2904
231
232! Relationship between ATBice and ATBperp,ice for liquid particles
233!  ATBperp,ice = Beta*ATBice^2 + Gamma*ATBice
234         Beta = 0.4099
235         Gamma = 0.009
236
237!------------------------------------------------------------
238!---- 1. Preliminary definitions and calculations :
239!------------------------------------------------------------
240
241      if ( npart .ne. 4 ) then
242        print *,'Error in lidar_simulator, npart should be 4, not',npart
243        stop
244      endif
245
246      pi = dacos(-1.D0)
247
248! Polynomial coefficients for spherical liq/ice particles derived from Mie theory.
249! Polynomial coefficients for non spherical particles derived from a composite of
250! Ray-tracing theory for large particles (e.g. Noel et al., Appl. Opt., 2001)
251! and FDTD theory for very small particles (Yang et al., JQSRT, 2003).
252
253! We repeat the same coefficients for LS and CONV cloud to make code more readable
254!*     LS Liquid water coefficients:
255         polpart(INDX_LSLIQ,1) =  2.6980e-8
256         polpart(INDX_LSLIQ,2) = -3.7701e-6
257         polpart(INDX_LSLIQ,3) =  1.6594e-4
258         polpart(INDX_LSLIQ,4) = -0.0024
259         polpart(INDX_LSLIQ,5) =  0.0626
260!*     LS Ice coefficients:
261      if (ice_type.eq.0) then     
262         polpart(INDX_LSICE,1) = -1.0176e-8
263         polpart(INDX_LSICE,2) =  1.7615e-6
264         polpart(INDX_LSICE,3) = -1.0480e-4
265         polpart(INDX_LSICE,4) =  0.0019
266         polpart(INDX_LSICE,5) =  0.0460
267      endif
268!*     LS Ice NS coefficients:
269      if (ice_type.eq.1) then
270         polpart(INDX_LSICE,1) = 1.3615e-8
271         polpart(INDX_LSICE,2) = -2.04206e-6
272         polpart(INDX_LSICE,3) = 7.51799e-5
273         polpart(INDX_LSICE,4) = 0.00078213
274         polpart(INDX_LSICE,5) = 0.0182131
275      endif
276!*     CONV Liquid water coefficients:
277         polpart(INDX_CVLIQ,1) =  2.6980e-8
278         polpart(INDX_CVLIQ,2) = -3.7701e-6
279         polpart(INDX_CVLIQ,3) =  1.6594e-4
280         polpart(INDX_CVLIQ,4) = -0.0024
281         polpart(INDX_CVLIQ,5) =  0.0626
282!*     CONV Ice coefficients:
283      if (ice_type.eq.0) then
284         polpart(INDX_CVICE,1) = -1.0176e-8
285         polpart(INDX_CVICE,2) =  1.7615e-6
286         polpart(INDX_CVICE,3) = -1.0480e-4
287         polpart(INDX_CVICE,4) =  0.0019
288         polpart(INDX_CVICE,5) =  0.0460
289      endif
290      if (ice_type.eq.1) then
291         polpart(INDX_CVICE,1) = 1.3615e-8
292         polpart(INDX_CVICE,2) = -2.04206e-6
293         polpart(INDX_CVICE,3) = 7.51799e-5
294         polpart(INDX_CVICE,4) = 0.00078213
295         polpart(INDX_CVICE,5) = 0.0182131
296      endif
297
298! density:
299!*    clear-sky air:
300      rhoair = pres/(287.04*temp)
301
302!*    liquid/ice particules:
303      rhopart(INDX_LSLIQ) = rholiq
304      rhopart(INDX_LSICE) = rhoice
305      rhopart(INDX_CVLIQ) = rholiq
306      rhopart(INDX_CVICE) = rhoice
307
308! effective radius particles:
309      rad_part(:,:,INDX_LSLIQ) = ls_radliq(:,:)
310      rad_part(:,:,INDX_LSICE) = ls_radice(:,:)
311      rad_part(:,:,INDX_CVLIQ) = cv_radliq(:,:)
312      rad_part(:,:,INDX_CVICE) = cv_radice(:,:)
313      rad_part(:,:,:)=MAX(rad_part(:,:,:),0.)
314      rad_part(:,:,:)=MIN(rad_part(:,:,:),70.0e-6)
315     
316! altitude at half pressure levels:
317      zheight(:,1) = 0.0
318      do k = 2, nlev+1
319        zheight(:,k) = zheight(:,k-1) &
320                  -(presf(:,k)-presf(:,k-1))/(rhoair(:,k-1)*9.81)
321      enddo
322
323!------------------------------------------------------------
324!---- 2. Molecular alpha and beta:
325!------------------------------------------------------------
326
327      beta_mol = pres/km/temp * Cmol
328      alpha_mol = 8.0*pi/3.0 * beta_mol
329
330!------------------------------------------------------------
331!---- 3. Particles alpha and beta:
332!------------------------------------------------------------
333
334! polynomes kp_lidar derived from Mie theory:
335      do i = 1, npart
336       where ( rad_part(:,:,i).gt.0.0)
337         kp_part(:,:,i) = &
338            polpart(i,1)*(rad_part(:,:,i)*1e6)**4 &
339          + polpart(i,2)*(rad_part(:,:,i)*1e6)**3 &
340          + polpart(i,3)*(rad_part(:,:,i)*1e6)**2 &
341          + polpart(i,4)*(rad_part(:,:,i)*1e6) &
342          + polpart(i,5)
343        elsewhere
344         kp_part(:,:,i) = 0.
345        endwhere
346      enddo
347     
348! mixing ratio particules in each subcolumn:
349          qpart(:,:,INDX_LSLIQ) = q_lsliq(:,:) ! oct08
350          qpart(:,:,INDX_LSICE) = q_lsice(:,:) ! oct08
351          qpart(:,:,INDX_CVLIQ) = q_cvliq(:,:) ! oct08
352          qpart(:,:,INDX_CVICE) = q_cvice(:,:) ! oct08
353
354! alpha of particles in each subcolumn:
355      do i = 1, npart
356        where ( rad_part(:,:,i).gt.0.0)
357          alpha_part(:,:,i) = 3.0/4.0 * Qscat &
358                 * rhoair(:,:) * qpart(:,:,i) &
359                 / (rhopart(i) * rad_part(:,:,i) )
360        elsewhere
361          alpha_part(:,:,i) = 0.
362        endwhere
363      enddo
364
365!------------------------------------------------------------
366!---- 4.1 Total Backscatter signal:
367!------------------------------------------------------------
368
369! optical thickness (molecular):
370!     opt. thick of each layer
371      tau_mol(:,1:nlev) = alpha_mol(:,1:nlev) &
372         & *(zheight(:,2:nlev+1)-zheight(:,1:nlev))
373!     opt. thick from TOA
374      DO k = nlev-1, 1, -1
375        tau_mol(:,k) = tau_mol(:,k) + tau_mol(:,k+1)
376      ENDDO
377
378! optical thickness (particles):
379
380      tau_part = rdiffm * alpha_part
381      DO i = 1, npart
382!       opt. thick of each layer
383        tau_part(:,:,i) = tau_part(:,:,i) &
384           & * (zheight(:,2:nlev+1)-zheight(:,1:nlev) )
385!       opt. thick from TOA
386        DO k = nlev-1, 1, -1
387          tau_part(:,k,i) = tau_part(:,k,i) + tau_part(:,k+1,i)
388        ENDDO
389      ENDDO
390
391! molecular signal:
392!      Upper layer
393       pmol(:,nlev) = beta_mol(:,nlev) / (2.*tau_mol(:,nlev)) &
394            & * (1.-exp(-2.0*tau_mol(:,nlev)))
395!      Other layers
396       DO k= nlev-1, 1, -1
397        tau_mol_lay(:) = tau_mol(:,k)-tau_mol(:,k+1) ! opt. thick. of layer k
398        WHERE (tau_mol_lay(:).GT.0.)
399          pmol(:,k) = beta_mol(:,k) * EXP(-2.0*tau_mol(:,k+1)) / (2.*tau_mol_lay(:)) &
400            & * (1.-exp(-2.0*tau_mol_lay(:)))
401        ELSEWHERE
402!         This must never happend, but just in case, to avoid div. by 0
403          pmol(:,k) = beta_mol(:,k) * EXP(-2.0*tau_mol(:,k+1))
404        END WHERE
405      END DO
406
407! Total signal (molecular + particules):
408!
409!
410! For performance reason on vector computers, the 2 following lines should not be used
411! and should be replace by the later one.
412!      betatot(:,:) = beta_mol(:,:) + sum(kp_part*alpha_part,dim=3)
413!      tautot(:,:)  = tau_mol(:,:)  + sum(tau_part,dim=3)
414      betatot(:,:) = beta_mol(:,:)
415      tautot(:,:)  = tau_mol(:,:)
416      DO i = 1, npart
417           betatot(:,:) = betatot(:,:) + kp_part(:,:,i)*alpha_part(:,:,i)
418           tautot(:,:) = tautot(:,:)  + tau_part(:,:,i)
419      ENDDO ! i
420!
421!     Upper layer
422      pnorm(:,nlev) = betatot(:,nlev) / (2.*tautot(:,nlev)) &
423            & * (1.-exp(-2.0*tautot(:,nlev)))
424
425!     Other layers
426      DO k= nlev-1, 1, -1
427          tautot_lay(:) = tautot(:,k)-tautot(:,k+1) ! optical thickness of layer k
428        WHERE (tautot_lay(:).GT.0.)
429          pnorm(:,k) = betatot(:,k) * EXP(-2.0*tautot(:,k+1)) / (2.*tautot_lay(:)) &
430               & * (1.-EXP(-2.0*tautot_lay(:)))
431        ELSEWHERE
432!         This must never happend, but just in case, to avoid div. by 0
433          pnorm(:,k) = betatot(:,k) * EXP(-2.0*tautot(:,k+1))
434        END WHERE
435      END DO
436
437!------------------------------------------------------------
438!---- 4.2 Ice/Liq Backscatter signal:
439!------------------------------------------------------------
440
441! Contribution of the molecular to beta
442      betatot_ice(:,:) = beta_mol(:,:)
443      betatot_liq(:,:) = beta_mol(:,:)
444
445      tautot_ice(:,:) = tau_mol(:,:)
446      tautot_liq(:,:) = tau_mol(:,:)
447
448      DO i = 2, npart,2
449           betatot_ice(:,:) = betatot_ice(:,:)+ kp_part(:,:,i)*alpha_part(:,:,i)
450           tautot_ice(:,:) = tautot_ice(:,:)  + tau_part(:,:,i)
451      ENDDO ! i
452      DO i = 1, npart,2
453           betatot_liq(:,:) = betatot_liq(:,:)+ kp_part(:,:,i)*alpha_part(:,:,i)
454           tautot_liq(:,:) = tautot_liq(:,:)  + tau_part(:,:,i)
455      ENDDO ! i
456
457
458! Computation of the ice and liquid lidar backscattered signal (ATBice and ATBliq)
459!     Ice only
460!     Upper layer
461      pnorm_ice(:,nlev) = betatot_ice(:,nlev) / (2.*tautot_ice(:,nlev)) &
462            & * (1.-exp(-2.0*tautot_ice(:,nlev)))
463
464      DO k= nlev-1, 1, -1
465          tautot_lay_ice(:) = tautot_ice(:,k)-tautot_ice(:,k+1)
466        WHERE (tautot_lay_ice(:).GT.0.)
467         pnorm_ice(:,k)=betatot_ice(:,k)*EXP(-2.0*tautot_ice(:,k+1))/(2.*tautot_lay_ice(:)) &
468               & * (1.-EXP(-2.0*tautot_lay_ice(:)))
469        ELSEWHERE
470         pnorm_ice(:,k)=betatot_ice(:,k)*EXP(-2.0*tautot_ice(:,k+1))
471        END WHERE
472      ENDDO
473
474!     Liquid only
475!     Upper layer
476      pnorm_liq(:,nlev) = betatot_liq(:,nlev) / (2.*tautot_liq(:,nlev)) &
477            & * (1.-exp(-2.0*tautot_liq(:,nlev)))
478
479      DO k= nlev-1, 1, -1
480          tautot_lay_liq(:) = tautot_liq(:,k)-tautot_liq(:,k+1)
481        WHERE (tautot_lay_liq(:).GT.0.)
482          pnorm_liq(:,k)=betatot_liq(:,k)*EXP(-2.0*tautot_liq(:,k+1))/(2.*tautot_lay_liq(:)) &
483               & * (1.-EXP(-2.0*tautot_lay_liq(:)))
484        ELSEWHERE
485          pnorm_liq(:,k)=betatot_liq(:,k)*EXP(-2.0*tautot_liq(:,k+1))
486        END WHERE
487      ENDDO
488
489
490! Computation of ATBperp,ice/liq from ATBice/liq including the multiple scattering
491! contribution (Cesana and Chepfer 2013, JGR)
492!  ATBperp,ice = Alpha*ATBice
493!  ATBperp,liq = Beta*ATBliq^2 + Gamma*ATBliq
494
495      DO k= nlev, 1, -1
496              pnorm_perp_ice(:,k) = Alpha * pnorm_ice(:,k) ! Ice particles
497              pnorm_perp_liq(:,k) = 1000*Beta * pnorm_liq(:,k)**2 + Gamma * pnorm_liq(:,k) ! Liquid particles
498      ENDDO
499
500! Computation of beta_perp_ice/liq using the lidar equation
501!     Ice only
502!     Upper layer
503      beta_perp_ice(:,nlev) = pnorm_perp_ice(:,nlev) * (2.*tautot_ice(:,nlev)) &
504            & / (1.-exp(-2.0*tautot_ice(:,nlev)))
505
506      DO k= nlev-1, 1, -1
507        tautot_lay_ice(:) = tautot_ice(:,k)-tautot_ice(:,k+1)
508        WHERE (tautot_lay_ice(:).GT.0.)
509         beta_perp_ice(:,k) = pnorm_perp_ice(:,k)/ EXP(-2.0*tautot_ice(:,k+1)) * (2.*tautot_lay_ice(:)) &
510            & / (1.-exp(-2.0*tautot_lay_ice(:)))
511
512        ELSEWHERE
513         beta_perp_ice(:,k)=pnorm_perp_ice(:,k)/EXP(-2.0*tautot_ice(:,k+1))
514        END WHERE
515      ENDDO
516
517!     Liquid only
518!     Upper layer
519      beta_perp_liq(:,nlev) = pnorm_perp_liq(:,nlev) * (2.*tautot_liq(:,nlev)) &
520            & / (1.-exp(-2.0*tautot_liq(:,nlev)))
521
522      DO k= nlev-1, 1, -1
523          tautot_lay_liq(:) = tautot_liq(:,k)-tautot_liq(:,k+1)
524        WHERE (tautot_lay_liq(:).GT.0.)
525         beta_perp_liq(:,k) = pnorm_perp_liq(:,k)/ EXP(-2.0*tautot_liq(:,k+1)) * (2.*tautot_lay_liq(:)) &
526            & / (1.-exp(-2.0*tautot_lay_liq(:)))
527
528        ELSEWHERE
529         beta_perp_liq(:,k)=pnorm_perp_liq(:,k)/EXP(-2.0*tautot_liq(:,k+1))
530        END WHERE
531      ENDDO
532
533
534
535!------------------------------------------------------------
536!---- 4.3 Perpendicular Backscatter signal:
537!------------------------------------------------------------
538
539! Computation of the total perpendicular lidar signal (ATBperp for liq+ice)
540!     Upper layer
541    WHERE(tautot(:,nlev).GT.0)
542          pnorm_perp_tot(:,nlev) = &
543              (beta_perp_ice(:,nlev)+beta_perp_liq(:,nlev)-(beta_mol(:,nlev)/(1+1/0.0284))) / (2.*tautot(:,nlev)) &
544              & * (1.-exp(-2.0*tautot(:,nlev)))
545    ELSEWHERE
546    pnorm_perp_tot(:,nlev) = 0.
547    ENDWHERE
548
549!     Other layers
550      DO k= nlev-1, 1, -1
551          tautot_lay(:) = tautot(:,k)-tautot(:,k+1) ! optical thickness of layer k
552
553          ! The perpendicular component of the molecular backscattered signal (Betaperp) has been
554          ! taken into account two times (once for liquid and once for ice).
555          ! We remove one contribution using
556          ! Betaperp=beta_mol(:,k)/(1+1/0.0284)) [bodhaine et al. 1999] in the following equations:
557            WHERE (pnorm(:,k).eq.0)
558                  pnorm_perp_tot(:,k)=0.
559                  ELSEWHERE
560                    WHERE (tautot_lay(:).GT.0.)
561                      pnorm_perp_tot(:,k) = &
562                          (beta_perp_ice(:,k)+beta_perp_liq(:,k)-(beta_mol(:,k)/(1+1/0.0284))) * &
563                          EXP(-2.0*tautot(:,k+1)) / (2.*tautot_lay(:)) &
564                          & * (1.-EXP(-2.0*tautot_lay(:)))
565                    ELSEWHERE
566          !         This must never happen, but just in case, to avoid div. by 0
567                      pnorm_perp_tot(:,k) = &
568                           (beta_perp_ice(:,k)+beta_perp_liq(:,k)-(beta_mol(:,k)/(1+1/0.0284))) * &
569                          EXP(-2.0*tautot(:,k+1))
570                    END WHERE
571            ENDWHERE
572
573      END DO
574
575!-------- End computation Lidar --------------------------
576
577!---------------------------------------------------------
578!  Parasol/Polder module
579!
580!  Purpose : Compute reflectance for one particular viewing direction
581!  and 5 solar zenith angles (calculation valid only over ocean)
582! ---------------------------------------------------------
583
584! initialization:
585    refl(:,:) = 0.0
586
587! activate parasol calculations:
588    if (ok_parasol) then
589
590!     Optical thickness from TOA to surface
591      tautot_S_liq = 0.
592      tautot_S_ice = 0.
593      tautot_S_liq(:) = tautot_S_liq(:) &
594         + tau_part(:,1,1) + tau_part(:,1,3)
595      tautot_S_ice(:) = tautot_S_ice(:) &
596         + tau_part(:,1,2) + tau_part(:,1,4)
597
598      call parasol(npoints,nrefl,undef  &
599                 ,tautot_S_liq,tautot_S_ice &
600                 ,refl)
601
602    endif ! ok_parasol
603
604  END SUBROUTINE lidar_simulator
605!
606!---------------------------------------------------------------------------------
607!
608  SUBROUTINE parasol(npoints,nrefl,undef  &
609                       ,tautot_S_liq,tautot_S_ice  &
610                       ,refl)
611!---------------------------------------------------------------------------------
612! Purpose: To compute Parasol reflectance signal from model-simulated profiles
613!          of cloud water and cloud fraction in each sub-column of each model
614!          gridbox.
615!
616!
617! December 2008, S. Bony,  H. Chepfer and J-L. Dufresne :
618! - optimization for vectorization
619!
620! Version 2.0 (October 2008)
621! Version 2.1 (December 2008)
622!---------------------------------------------------------------------------------
623
624    IMPLICIT NONE
625
626! inputs
627    INTEGER npoints              ! Number of horizontal gridpoints
628    INTEGER nrefl                ! Number of angles for which the reflectance
629                                 ! is computed. Can not be greater then ntetas
630    REAL undef                   ! Undefined value. Currently not used
631    REAL tautot_S_liq(npoints)   ! liquid water cloud optical thickness,
632                                   ! integrated from TOA to surface
633    REAL tautot_S_ice(npoints)   ! same for ice water clouds only
634! outputs
635    REAL refl(npoints,nrefl)     ! Parasol reflectances
636!
637! Local variables
638    REAL tautot_S(npoints)       ! cloud optical thickness, from TOA to surface
639    REAL frac_taucol_liq(npoints), frac_taucol_ice(npoints)
640
641    REAL pi
642!   look up table variables:
643    INTEGER ny, it
644    INTEGER ntetas, nbtau        ! number of angle and of optical thickness
645                                   ! of the look-up table
646    PARAMETER (ntetas=5, nbtau=7)
647    REAL aa(ntetas,nbtau-1), ab(ntetas,nbtau-1)
648    REAL ba(ntetas,nbtau-1), bb(ntetas,nbtau-1) 
649    REAL tetas(ntetas),tau(nbtau)                       
650    REAL r_norm(ntetas)
651    REAL rlumA(ntetas,nbtau), rlumB(ntetas,nbtau)       
652    REAL rlumA_mod(npoints,5), rlumB_mod(npoints,5)
653
654    DATA tau   /0., 1., 5., 10., 20., 50., 100./
655    DATA tetas /0., 20., 40., 60., 80./
656   
657! Look-up table for spherical liquid particles:
658    DATA (rlumA(1,ny),ny=1,nbtau) /0.03, 0.090886, 0.283965, &
659     0.480587, 0.695235, 0.908229, 1.0 /
660    DATA (rlumA(2,ny),ny=1,nbtau) /0.03, 0.072185, 0.252596, &
661      0.436401,  0.631352, 0.823924, 0.909013 /
662    DATA (rlumA(3,ny),ny=1,nbtau) /0.03, 0.058410, 0.224707, &
663      0.367451,  0.509180, 0.648152, 0.709554 /
664    DATA (rlumA(4,ny),ny=1,nbtau) /0.03, 0.052498, 0.175844, &
665      0.252916,  0.326551, 0.398581, 0.430405 /
666    DATA (rlumA(5,ny),ny=1,nbtau) /0.03, 0.034730, 0.064488, &
667      0.081667,  0.098215, 0.114411, 0.121567 /
668
669! Look-up table for ice particles:
670    DATA (rlumB(1,ny),ny=1,nbtau) /0.03, 0.092170, 0.311941, &
671       0.511298, 0.712079 , 0.898243 , 0.976646 /
672    DATA (rlumB(2,ny),ny=1,nbtau) /0.03, 0.087082, 0.304293, &
673       0.490879,  0.673565, 0.842026, 0.912966 /
674    DATA (rlumB(3,ny),ny=1,nbtau) /0.03, 0.083325, 0.285193, &
675      0.430266,  0.563747, 0.685773,  0.737154 /
676    DATA (rlumB(4,ny),ny=1,nbtau) /0.03, 0.084935, 0.233450, &
677      0.312280, 0.382376, 0.446371, 0.473317 /
678    DATA (rlumB(5,ny),ny=1,nbtau) /0.03, 0.054157, 0.089911, &
679      0.107854, 0.124127, 0.139004, 0.145269 /
680
681!--------------------------------------------------------------------------------
682! Lum_norm=f(tetaS,tau_cloud) derived from adding-doubling calculations
683!        valid ONLY ABOVE OCEAN (albedo_sfce=5%)
684!        valid only in one viewing direction (theta_v=30�, phi_s-phi_v=320�)
685!        based on adding-doubling radiative transfer computation
686!        for tau values (0 to 100) and for tetas values (0 to 80)
687!        for 2 scattering phase functions: liquid spherical, ice non spherical
688
689    IF ( nrefl.GT. ntetas ) THEN
690        PRINT *,'Error in lidar_simulator, nrefl should be less then ',ntetas,' not',nrefl
691        STOP
692    ENDIF
693
694    rlumA_mod=0
695    rlumB_mod=0
696!
697    pi = ACOS(-1.0)
698    r_norm(:)=1./ COS(pi/180.*tetas(:))
699!
700    tautot_S_liq(:)=MAX(tautot_S_liq(:),tau(1))
701    tautot_S_ice(:)=MAX(tautot_S_ice(:),tau(1))
702    tautot_S(:) = tautot_S_ice(:) + tautot_S_liq(:)
703!
704! relative fraction of the opt. thick due to liquid or ice clouds
705    WHERE (tautot_S(:) .GT. 0.)
706        frac_taucol_liq(:) = tautot_S_liq(:) / tautot_S(:)
707        frac_taucol_ice(:) = tautot_S_ice(:) / tautot_S(:)
708    ELSEWHERE
709        frac_taucol_liq(:) = 1.
710        frac_taucol_ice(:) = 0.
711    END WHERE
712    tautot_S(:)=MIN(tautot_S(:),tau(nbtau))
713!
714! Linear interpolation :
715
716    DO ny=1,nbtau-1
717! microphysics A (liquid clouds)
718      aA(:,ny) = (rlumA(:,ny+1)-rlumA(:,ny))/(tau(ny+1)-tau(ny))
719      bA(:,ny) = rlumA(:,ny) - aA(:,ny)*tau(ny)
720! microphysics B (ice clouds)
721      aB(:,ny) = (rlumB(:,ny+1)-rlumB(:,ny))/(tau(ny+1)-tau(ny))
722      bB(:,ny) = rlumB(:,ny) - aB(:,ny)*tau(ny)
723    ENDDO
724!
725    DO it=1,ntetas
726      DO ny=1,nbtau-1
727        WHERE (tautot_S(:).GE.tau(ny).AND.tautot_S(:).LE.tau(ny+1))
728            rlumA_mod(:,it) = aA(it,ny)*tautot_S(:) + bA(it,ny)
729            rlumB_mod(:,it) = aB(it,ny)*tautot_S(:) + bB(it,ny)
730        END WHERE
731      END DO
732    END DO
733!
734    DO it=1,ntetas
735      refl(:,it) = frac_taucol_liq(:) * rlumA_mod(:,it) &
736         + frac_taucol_ice(:) * rlumB_mod(:,it)
737! normalized radiance -> reflectance:
738      refl(:,it) = refl(:,it) * r_norm(it)
739    ENDDO
740
741    RETURN
742  END SUBROUTINE parasol
Note: See TracBrowser for help on using the repository browser.