Context Navigation

lidar_simulator.F90 @ 5157

Last change on this file since 5157 was 5099, checked in by abarral, 4 months ago
Replace most uses of CPP_DUST by the corresponding logical defined in lmdz_cppkeys_wrapper.F90 Convert several files from .F to .f90 to allow Dust to compile w/o rrtm/ecrad Create lmdz_yoerad.f90 (lint) Remove "!" on otherwise empty line
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	!------------------------------------------------------------
212	betatot_ice(:,:)=0
213	betatot_liq(:,:)=0
214	beta_perp_ice(:,:)=0
215	beta_perp_liq(:,:)=0
216	tautot_ice(:,:)=0
217	tautot_liq(:,:)=0
218	tautot_lay_ice(:)=0;
219	tautot_lay_liq(:)=0;
220	pnorm_liq(:,:)=0
221	pnorm_ice(:,:)=0
222	pnorm_perp_ice(:,:)=0
223	pnorm_perp_liq(:,:)=0
224	pnorm_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 = BetaATBice^2 + GammaATBice
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.0pi/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.0tau_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.0tautot(:,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.0tautot_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.0tautot_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.0tautot_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.0tautot_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 = BetaATBliq^2 + GammaATBliq
502
503	DO k= nlev, 1, -1
504	pnorm_perp_ice(:,k) = Alpha * pnorm_ice(:,k) ! Ice particles
505	pnorm_perp_liq(:,k) = 1000Beta 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.0tautot_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.0tautot_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.0tautot(:,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.

Download in other formats:

Original Format