Context Navigation

lidar_simulator.F90 @ 5451

Last change on this file since 5451 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	!------------------------------------------------------------
209	betatot_ice(:,:)=0
210	betatot_liq(:,:)=0
211	beta_perp_ice(:,:)=0
212	beta_perp_liq(:,:)=0
213	tautot_ice(:,:)=0
214	tautot_liq(:,:)=0
215	tautot_lay_ice(:)=0;
216	tautot_lay_liq(:)=0;
217	pnorm_liq(:,:)=0
218	pnorm_ice(:,:)=0
219	pnorm_perp_ice(:,:)=0
220	pnorm_perp_liq(:,:)=0
221	pnorm_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 = BetaATBice^2 + GammaATBice
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.0pi/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.0tau_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.0tautot(:,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.0tautot_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.0tautot_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.0tautot_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.0tautot_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 = BetaATBliq^2 + GammaATBliq
494
495	DO k= nlev, 1, -1
496	pnorm_perp_ice(:,k) = Alpha * pnorm_ice(:,k) ! Ice particles
497	pnorm_perp_liq(:,k) = 1000Beta 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.0tautot_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.0tautot_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.0tautot(:,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.

Download in other formats:

Original Format