Context Navigation

lidar_simulator.F90 @ 1865

Last change on this file since 1865 was 1526, checked in by idelkadi, 13 years ago

Les modificaions concerne la partie Lidar du simulateur COSP :

Correction dans phys_cosp.F90 du bug sur le calcul des altitudes aux milieux et aux interfaces des couches
Rajout du teste "Lidar parfait" (iflag_testlidar=.true. dans lidar_simulator.F90)
Correction d'une erreur sur la dimension du tableau frac_out dans cosp_lidar.F90
Sories sur les niveaux presnivs du modele dans le cas ou l'on desactive l'interpolation sur les 40 niveaux COSP

File size: 22.2 KB

Line
1	! Copyright (c) 2009, Centre National de la Recherche Scientifique
2	! All rights reserved.
3	!
4	! Redistribution and use in source and binary forms, with or without modification, are permitted
5	! provided that the following conditions are met:
6	!
7	! * Redistributions of source code must retain the above copyright notice, this list
8	! of conditions and the following disclaimer.
9	! * Redistributions in binary form must reproduce the above copyright notice, this list
10	! of conditions and the following disclaimer in the documentation and/or other materials
11	! provided with the distribution.
12	! * Neither the name of the LMD/IPSL/CNRS/UPMC nor the names of its
13	! contributors may be used to endorse or promote products derived from this software without
14	! specific prior written permission.
15	!
16	! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR
17	! IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND
18	! FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR
19	! CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
20	! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
21	! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER
22	! IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
23	! OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
24
25	SUBROUTINE lidar_simulator(npoints,nlev,npart,nrefl &
26	, undef &
27	, pres, presf, temp &
28	, q_lsliq, q_lsice, q_cvliq, q_cvice &
29	, ls_radliq, ls_radice, cv_radliq, cv_radice &
30	, frac_out, ice_type &
31	, pmol, pnorm, tautot, refl )
32	!
33	!---------------------------------------------------------------------------------
34	! Purpose: To compute lidar signal from model-simulated profiles of cloud water
35	! and cloud fraction in each sub-column of each model gridbox.
36	!
37	! References:
38	! Chepfer H., S. Bony, D. Winker, M. Chiriaco, J.-L. Dufresne, G. Seze (2008),
39	! Use of CALIPSO lidar observations to evaluate the cloudiness simulated by a
40	! climate model, Geophys. Res. Lett., 35, L15704, doi:10.1029/2008GL034207.
41	!
42	! Previous references:
43	! Chiriaco et al, MWR, 2006; Chepfer et al., MWR, 2007
44	!
45	! Contacts: Helene Chepfer (chepfer@lmd.polytechnique.fr), Sandrine Bony (bony@lmd.jussieu.fr)
46	!
47	! May 2007: ActSim code of M. Chiriaco and H. Chepfer rewritten by S. Bony
48	!
49	! May 2008, H. Chepfer:
50	! - Units of pressure inputs: Pa
51	! - Non Spherical particles : LS Ice NS coefficients, CONV Ice NS coefficients
52	! - New input: ice_type (0=ice-spheres ; 1=ice-non-spherical)
53	!
54	! June 2008, A. Bodas-Salcedo:
55	! - Ported to Fortran 90 and optimisation changes
56	!
57	! August 2008, J-L Dufresne:
58	! - Optimisation changes (sum instructions suppressed)
59	!
60	! October 2008, S. Bony, H. Chepfer and J-L. Dufresne :
61	! - Interface with COSP v2.0:
62	! cloud fraction removed from inputs
63	! in-cloud condensed water now in input (instead of grid-averaged value)
64	! depolarisation diagnostic removed
65	! parasol (polder) reflectances (for 5 different solar zenith angles) added
66	!
67	! December 2008, S. Bony, H. Chepfer and J-L. Dufresne :
68	! - Modification of the integration of the lidar equation.
69	! - change the cloud detection threshold
70	!
71	! April 2008, A. Bodas-Salcedo:
72	! - Bug fix in computation of pmol and pnorm of upper layer
73	!
74	! April 2008, J-L. Dufresne
75	! - Bug fix in computation of pmol and pnorm, thanks to Masaki Satoh: a factor 2
76	! was missing. This affects the ATB values but not the cloud fraction.
77	!
78	!---------------------------------------------------------------------------------
79	!
80	! Inputs:
81	! npoints : number of horizontal points
82	! nlev : number of vertical levels
83	! npart: numberb of cloud meteors (stratiform_liq, stratiform_ice, conv_liq, conv_ice).
84	! Currently npart must be 4
85	! nrefl: number of solar zenith angles for parasol reflectances
86	! pres : pressure in the middle of atmospheric layers (full levels): Pa
87	! presf: pressure in the interface of atmospheric layers (half levels): Pa
88	! presf(..,1) : surface pressure ; presf(..,nlev+1)= TOA pressure
89	! temp : temperature of atmospheric layers: K
90	! q_lsliq: LS sub-column liquid water mixing ratio (kg/kg)
91	! q_lsice: LS sub-column ice water mixing ratio (kg/kg)
92	! q_cvliq: CONV sub-column liquid water mixing ratio (kg/kg)
93	! q_cvice: CONV sub-column ice water mixing ratio (kg/kg)
94	! ls_radliq: effective radius of LS liquid particles (meters)
95	! ls_radice: effective radius of LS ice particles (meters)
96	! cv_radliq: effective radius of CONV liquid particles (meters)
97	! cv_radice: effective radius of CONV ice particles (meters)
98	! frac_out : cloud cover in each sub-column of the gridbox (output from scops)
99	! ice_type : ice particle shape hypothesis (ice_type=0 for spheres, ice_type=1
100	! for non spherical particles)
101	!
102	! Outputs:
103	! pmol : molecular attenuated backscatter lidar signal power (m^-1.sr^-1)
104	! pnorm: total attenuated backscatter lidar signal power (m^-1.sr^-1)
105	! tautot: optical thickess integrated from top to level z
106	! refl : parasol(polder) reflectance
107	!
108	! Version 1.0 (June 2007)
109	! Version 1.1 (May 2008)
110	! Version 1.2 (June 2008)
111	! Version 2.0 (October 2008)
112	! Version 2.1 (December 2008)
113	!---------------------------------------------------------------------------------
114
115	IMPLICIT NONE
116	REAL :: SRsat
117	PARAMETER (SRsat = 0.01) ! threshold full attenuation
118
119	LOGICAL ok_parasol
120	PARAMETER (ok_parasol=.true.) ! set to .true. if you want to activate parasol reflectances
121
122	INTEGER i, k
123
124	INTEGER INDX_LSLIQ,INDX_LSICE,INDX_CVLIQ,INDX_CVICE
125	PARAMETER (INDX_LSLIQ=1,INDX_LSICE=2,INDX_CVLIQ=3,INDX_CVICE=4)
126	! inputs:
127	INTEGER npoints,nlev,npart,ice_type
128	INTEGER nrefl
129	real undef ! undefined value
130	REAL pres(npoints,nlev) ! pressure full levels
131	REAL presf(npoints,nlev+1) ! pressure half levels
132	REAL temp(npoints,nlev)
133	REAL q_lsliq(npoints,nlev), q_lsice(npoints,nlev)
134	REAL q_cvliq(npoints,nlev), q_cvice(npoints,nlev)
135	REAL ls_radliq(npoints,nlev), ls_radice(npoints,nlev)
136	REAL cv_radliq(npoints,nlev), cv_radice(npoints,nlev)
137	REAL frac_out(npoints,nlev)
138
139	! outputs (for each subcolumn):
140
141	REAL pmol(npoints,nlev) ! molecular backscatter signal power (m^-1.sr^-1)
142	REAL pnorm(npoints,nlev) ! total lidar backscatter signal power (m^-1.sr^-1)
143	REAL tautot(npoints,nlev)! optical thickess integrated from top
144	REAL refl(npoints,nrefl)! parasol reflectance ! parasol
145
146	! actsim variables:
147
148	REAL km, rdiffm, Qscat, Cmol
149	PARAMETER (Cmol = 6.2446e-32) ! depends on wavelength
150	PARAMETER (km = 1.38e-23) ! Boltzmann constant (J/K)
151
152	PARAMETER (rdiffm = 0.7) ! multiple scattering correction parameter
153	PARAMETER (Qscat = 2.0) ! particle scattering efficiency at 532 nm
154
155	REAL rholiq, rhoice
156	PARAMETER (rholiq=1.0e+03) ! liquid water (kg/m3)
157	PARAMETER (rhoice=0.5e+03) ! ice (kg/m3)
158
159	REAL pi, rhopart(npart)
160	REAL polpart(npart,5) ! polynomial coefficients derived for spherical and non spherical
161	! particules
162
163	! grid-box variables:
164	REAL rad_part(npoints,nlev,npart)
165	REAL rhoair(npoints,nlev), zheight(npoints,nlev+1)
166	REAL beta_mol(npoints,nlev), alpha_mol(npoints,nlev)
167	REAL kp_part(npoints,nlev,npart)
168
169	! sub-column variables:
170	REAL frac_sub(npoints,nlev)
171	REAL qpart(npoints,nlev,npart) ! mixing ratio particles in each subcolumn
172	REAL alpha_part(npoints,nlev,npart)
173	REAL tau_mol_lay(npoints) ! temporary variable, moL. opt. thickness of layer k
174	REAL tau_mol(npoints,nlev) ! optical thickness between TOA and bottom of layer k
175	REAL tau_part(npoints,nlev,npart)
176	REAL betatot(npoints,nlev)
177	REAL tautot_lay(npoints) ! temporary variable, total opt. thickness of layer k
178	! Optical thickness from TOA to surface for Parasol
179	REAL tautot_S_liq(npoints),tautot_S_ice(npoints) ! for liq and ice clouds
180
181	! Abderrahmane 8-2-2011
182	Logical iflag_testlidar
183	PARAMETER (iflag_testlidar=.false.)
184
185	!------------------------------------------------------------
186	!---- 1. Preliminary definitions and calculations :
187	!------------------------------------------------------------
188
189	if ( npart .ne. 4 ) then
190	print *,'Error in lidar_simulator, npart should be 4, not',npart
191	stop
192	endif
193
194	pi = dacos(-1.D0)
195
196	! Polynomial coefficients for spherical liq/ice particles derived from Mie theory.
197	! Polynomial coefficients for non spherical particles derived from a composite of
198	! Ray-tracing theory for large particles (e.g. Noel et al., Appl. Opt., 2001)
199	! and FDTD theory for very small particles (Yang et al., JQSRT, 2003).
200
201	! We repeat the same coefficients for LS and CONV cloud to make code more readable
202	!* LS Liquid water coefficients:
203	polpart(INDX_LSLIQ,1) = 2.6980e-8
204	polpart(INDX_LSLIQ,2) = -3.7701e-6
205	polpart(INDX_LSLIQ,3) = 1.6594e-4
206	polpart(INDX_LSLIQ,4) = -0.0024
207	polpart(INDX_LSLIQ,5) = 0.0626
208	!* LS Ice coefficients:
209	if (ice_type.eq.0) then
210	polpart(INDX_LSICE,1) = -1.0176e-8
211	polpart(INDX_LSICE,2) = 1.7615e-6
212	polpart(INDX_LSICE,3) = -1.0480e-4
213	polpart(INDX_LSICE,4) = 0.0019
214	polpart(INDX_LSICE,5) = 0.0460
215	endif
216	!* LS Ice NS coefficients:
217	if (ice_type.eq.1) then
218	polpart(INDX_LSICE,1) = 1.3615e-8
219	polpart(INDX_LSICE,2) = -2.04206e-6
220	polpart(INDX_LSICE,3) = 7.51799e-5
221	polpart(INDX_LSICE,4) = 0.00078213
222	polpart(INDX_LSICE,5) = 0.0182131
223	endif
224	!* CONV Liquid water coefficients:
225	polpart(INDX_CVLIQ,1) = 2.6980e-8
226	polpart(INDX_CVLIQ,2) = -3.7701e-6
227	polpart(INDX_CVLIQ,3) = 1.6594e-4
228	polpart(INDX_CVLIQ,4) = -0.0024
229	polpart(INDX_CVLIQ,5) = 0.0626
230	!* CONV Ice coefficients:
231	if (ice_type.eq.0) then
232	polpart(INDX_CVICE,1) = -1.0176e-8
233	polpart(INDX_CVICE,2) = 1.7615e-6
234	polpart(INDX_CVICE,3) = -1.0480e-4
235	polpart(INDX_CVICE,4) = 0.0019
236	polpart(INDX_CVICE,5) = 0.0460
237	endif
238	if (ice_type.eq.1) then
239	polpart(INDX_CVICE,1) = 1.3615e-8
240	polpart(INDX_CVICE,2) = -2.04206e-6
241	polpart(INDX_CVICE,3) = 7.51799e-5
242	polpart(INDX_CVICE,4) = 0.00078213
243	polpart(INDX_CVICE,5) = 0.0182131
244	endif
245
246	! density:
247	!* clear-sky air:
248	rhoair = pres/(287.04*temp)
249
250	!* liquid/ice particules:
251	rhopart(INDX_LSLIQ) = rholiq
252	rhopart(INDX_LSICE) = rhoice
253	rhopart(INDX_CVLIQ) = rholiq
254	rhopart(INDX_CVICE) = rhoice
255
256	! effective radius particles:
257	rad_part(:,:,INDX_LSLIQ) = ls_radliq(:,:)
258	rad_part(:,:,INDX_LSICE) = ls_radice(:,:)
259	rad_part(:,:,INDX_CVLIQ) = cv_radliq(:,:)
260	rad_part(:,:,INDX_CVICE) = cv_radice(:,:)
261	rad_part(:,:,:)=MAX(rad_part(:,:,:),0.)
262	rad_part(:,:,:)=MIN(rad_part(:,:,:),70.0e-6)
263
264	! altitude at half pressure levels:
265	zheight(:,1) = 0.0
266	do k = 2, nlev+1
267	zheight(:,k) = zheight(:,k-1) &
268	-(presf(:,k)-presf(:,k-1))/(rhoair(:,k-1)*9.81)
269	enddo
270
271	! cloud fraction (0 or 1) in each sub-column:
272	! (if frac_out=1or2 -> frac_sub=1; if frac_out=0 -> frac_sub=0)
273	frac_sub = MIN( frac_out, 1.0 )
274
275	!------------------------------------------------------------
276	!---- 2. Molecular alpha and beta:
277	!------------------------------------------------------------
278
279	beta_mol = pres/km/temp * Cmol
280	alpha_mol = 8.0pi/3.0 beta_mol
281
282	!------------------------------------------------------------
283	!---- 3. Particles alpha and beta:
284	!------------------------------------------------------------
285
286	! polynomes kp_lidar derived from Mie theory:
287	do i = 1, npart
288	where ( rad_part(:,:,i).gt.0.0)
289	kp_part(:,:,i) = &
290	polpart(i,1)(rad_part(:,:,i)1e6)**4 &
291	+ polpart(i,2)(rad_part(:,:,i)1e6)**3 &
292	+ polpart(i,3)(rad_part(:,:,i)1e6)**2 &
293	+ polpart(i,4)(rad_part(:,:,i)1e6) &
294	+ polpart(i,5)
295	elsewhere
296	kp_part(:,:,i) = 0.
297	endwhere
298	enddo
299
300	! mixing ratio particules in each subcolumn:
301	qpart(:,:,INDX_LSLIQ) = q_lsliq(:,:) ! oct08
302	qpart(:,:,INDX_LSICE) = q_lsice(:,:) ! oct08
303	qpart(:,:,INDX_CVLIQ) = q_cvliq(:,:) ! oct08
304	qpart(:,:,INDX_CVICE) = q_cvice(:,:) ! oct08
305
306	! alpha of particles in each subcolumn:
307	do i = 1, npart
308	where ( rad_part(:,:,i).gt.0.0)
309	alpha_part(:,:,i) = 3.0/4.0 * Qscat &
310	* rhoair(:,:) * qpart(:,:,i) &
311	/ (rhopart(i) * rad_part(:,:,i) )
312	elsewhere
313	alpha_part(:,:,i) = 0.
314	endwhere
315	enddo
316
317	!------------------------------------------------------------
318	!---- 4. Backscatter signal:
319	!------------------------------------------------------------
320
321	! optical thickness (molecular):
322	! opt. thick of each layer
323	tau_mol(:,1:nlev) = alpha_mol(:,1:nlev) &
324	& *(zheight(:,2:nlev+1)-zheight(:,1:nlev))
325	! opt. thick from TOA
326	DO k = nlev-1, 1, -1
327	tau_mol(:,k) = tau_mol(:,k) + tau_mol(:,k+1)
328	ENDDO
329
330	! optical thickness (particles):
331
332	tau_part = rdiffm * alpha_part
333	DO i = 1, npart
334	! opt. thick of each layer
335	tau_part(:,:,i) = tau_part(:,:,i) &
336	& * (zheight(:,2:nlev+1)-zheight(:,1:nlev) )
337	! opt. thick from TOA
338	DO k = nlev-1, 1, -1
339	tau_part(:,k,i) = tau_part(:,k,i) + tau_part(:,k+1,i)
340	ENDDO
341	ENDDO
342
343	! molecular signal:
344	! Upper layer
345	pmol(:,nlev) = beta_mol(:,nlev) / (2.*tau_mol(:,nlev)) &
346	& * (1.-exp(-2.0*tau_mol(:,nlev)))
347	! Other layers
348	DO k= nlev-1, 1, -1
349	tau_mol_lay(:) = tau_mol(:,k)-tau_mol(:,k+1) ! opt. thick. of layer k
350	WHERE (tau_mol_lay(:).GT.0.)
351	pmol(:,k) = beta_mol(:,k) * EXP(-2.0tau_mol(:,k+1)) / (2.tau_mol_lay(:)) &
352	& * (1.-exp(-2.0*tau_mol_lay(:)))
353	ELSEWHERE
354	! This must never happend, but just in case, to avoid div. by 0
355	pmol(:,k) = beta_mol(:,k) * EXP(-2.0*tau_mol(:,k+1))
356	END WHERE
357	END DO
358	!
359	! Total signal (molecular + particules):
360	!
361	! For performance reason on vector computers, the 2 following lines should not be used
362	! and should be replace by the later one.
363	! betatot(:,:) = beta_mol(:,:) + sum(kp_part*alpha_part,dim=3)
364	! tautot(:,:) = tau_mol(:,:) + sum(tau_part,dim=3)
365	betatot(:,:) = beta_mol(:,:)
366	tautot(:,:) = tau_mol(:,:)
367	DO i = 1, npart
368	betatot(:,:) = betatot(:,:) + kp_part(:,:,i)*alpha_part(:,:,i)
369	tautot(:,:) = tautot(:,:) + tau_part(:,:,i)
370	ENDDO ! i
371	!
372	! Upper layer
373	pnorm(:,nlev) = betatot(:,nlev) / (2.*tautot(:,nlev)) &
374	& * (1.-exp(-2.0*tautot(:,nlev)))
375	! Other layers
376	DO k= nlev-1, 1, -1
377	tautot_lay(:) = tautot(:,k)-tautot(:,k+1) ! optical thickness of layer k
378	WHERE (tautot_lay(:).GT.0.)
379	pnorm(:,k) = betatot(:,k) * EXP(-2.0tautot(:,k+1)) / (2.tautot_lay(:)) &
380	!correc pnorm(:,k) = betatot(:,k) * EXP(-2.0*tautot(:,k+1)) & ! correc Satoh
381	!correc & / (2.0*tautot_lay(:)) & ! correc Satoh
382	& * (1.-EXP(-2.0*tautot_lay(:)))
383	ELSEWHERE
384	! This must never happend, but just in case, to avoid div. by 0
385	pnorm(:,k) = betatot(:,k) * EXP(-2.0*tautot(:,k+1))
386	END WHERE
387	END DO
388
389	if (iflag_testlidar) then
390	!+JLD test
391	! do k=1,nlev
392	! print*,'Min val de frac_out=',k,minval(frac_out(:,k))
393	! print*,'Max val de frac_out=',k,maxval(frac_out(:,k))
394	! enddo
395	where ( frac_out(:,:).ge.0.5)
396	! Correction AI 9 5 11 pnorm(:,:) = pmol(:,:)*10.
397	pnorm(:,:) = pmol(:,:)*50.
398	elsewhere
399	pnorm(:,:) = pmol(:,:)
400	endwhere
401	!-JLD test
402	endif
403
404	!-------- End computation Lidar --------------------------
405
406	!---------------------------------------------------------
407	! Parasol/Polder module
408	!
409	! Purpose : Compute reflectance for one particular viewing direction
410	! and 5 solar zenith angles (calculation valid only over ocean)
411	! ---------------------------------------------------------
412
413	! initialization:
414	refl(:,:) = 0.0
415
416	! activate parasol calculations:
417	if (ok_parasol) then
418
419	! Optical thickness from TOA to surface
420	tautot_S_liq = 0.
421	tautot_S_ice = 0.
422	tautot_S_liq(:) = tautot_S_liq(:) &
423	+ tau_part(:,1,1) + tau_part(:,1,3)
424	tautot_S_ice(:) = tautot_S_ice(:) &
425	+ tau_part(:,1,2) + tau_part(:,1,4)
426
427	call parasol(npoints,nrefl,undef &
428	,tautot_S_liq,tautot_S_ice &
429	,refl)
430
431	endif ! ok_parasol
432
433	END SUBROUTINE lidar_simulator
434	!
435	!---------------------------------------------------------------------------------
436	!
437	SUBROUTINE parasol(npoints,nrefl,undef &
438	,tautot_S_liq,tautot_S_ice &
439	,refl)
440	!---------------------------------------------------------------------------------
441	! Purpose: To compute Parasol reflectance signal from model-simulated profiles
442	! of cloud water and cloud fraction in each sub-column of each model
443	! gridbox.
444	!
445	!
446	! December 2008, S. Bony, H. Chepfer and J-L. Dufresne :
447	! - optimization for vectorization
448	!
449	! Version 2.0 (October 2008)
450	! Version 2.1 (December 2008)
451	!---------------------------------------------------------------------------------
452
453	IMPLICIT NONE
454
455	! inputs
456	INTEGER npoints ! Number of horizontal gridpoints
457	INTEGER nrefl ! Number of angles for which the reflectance
458	! is computed. Can not be greater then ntetas
459	REAL undef ! Undefined value. Currently not used
460	REAL tautot_S_liq(npoints) ! liquid water cloud optical thickness,
461	! integrated from TOA to surface
462	REAL tautot_S_ice(npoints) ! same for ice water clouds only
463	! outputs
464	REAL refl(npoints,nrefl) ! Parasol reflectances
465	!
466	! Local variables
467	REAL tautot_S(npoints) ! cloud optical thickness, from TOA to surface
468	REAL frac_taucol_liq(npoints), frac_taucol_ice(npoints)
469
470	REAL pi
471	! look up table variables:
472	INTEGER ny, it
473	INTEGER ntetas, nbtau ! number of angle and of optical thickness
474	! of the look-up table
475	PARAMETER (ntetas=5, nbtau=7)
476	REAL aa(ntetas,nbtau-1), ab(ntetas,nbtau-1)
477	REAL ba(ntetas,nbtau-1), bb(ntetas,nbtau-1)
478	REAL tetas(ntetas),tau(nbtau)
479	REAL r_norm(ntetas)
480	REAL rlumA(ntetas,nbtau), rlumB(ntetas,nbtau)
481	REAL rlumA_mod(npoints,5), rlumB_mod(npoints,5)
482
483	DATA tau /0., 1., 5., 10., 20., 50., 100./
484	DATA tetas /0., 20., 40., 60., 80./
485
486	! Look-up table for spherical liquid particles:
487	DATA (rlumA(1,ny),ny=1,nbtau) /0.03, 0.090886, 0.283965, &
488	0.480587, 0.695235, 0.908229, 1.0 /
489	DATA (rlumA(2,ny),ny=1,nbtau) /0.03, 0.072185, 0.252596, &
490	0.436401, 0.631352, 0.823924, 0.909013 /
491	DATA (rlumA(3,ny),ny=1,nbtau) /0.03, 0.058410, 0.224707, &
492	0.367451, 0.509180, 0.648152, 0.709554 /
493	DATA (rlumA(4,ny),ny=1,nbtau) /0.03, 0.052498, 0.175844, &
494	0.252916, 0.326551, 0.398581, 0.430405 /
495	DATA (rlumA(5,ny),ny=1,nbtau) /0.03, 0.034730, 0.064488, &
496	0.081667, 0.098215, 0.114411, 0.121567 /
497
498	! Look-up table for ice particles:
499	DATA (rlumB(1,ny),ny=1,nbtau) /0.03, 0.092170, 0.311941, &
500	0.511298, 0.712079 , 0.898243 , 0.976646 /
501	DATA (rlumB(2,ny),ny=1,nbtau) /0.03, 0.087082, 0.304293, &
502	0.490879, 0.673565, 0.842026, 0.912966 /
503	DATA (rlumB(3,ny),ny=1,nbtau) /0.03, 0.083325, 0.285193, &
504	0.430266, 0.563747, 0.685773, 0.737154 /
505	DATA (rlumB(4,ny),ny=1,nbtau) /0.03, 0.084935, 0.233450, &
506	0.312280, 0.382376, 0.446371, 0.473317 /
507	DATA (rlumB(5,ny),ny=1,nbtau) /0.03, 0.054157, 0.089911, &
508	0.107854, 0.124127, 0.139004, 0.145269 /
509
510	!--------------------------------------------------------------------------------
511	! Lum_norm=f(tetaS,tau_cloud) derived from adding-doubling calculations
512	! valid ONLY ABOVE OCEAN (albedo_sfce=5%)
513	! valid only in one viewing direction (theta_v=30�, phi_s-phi_v=320�)
514	! based on adding-doubling radiative transfer computation
515	! for tau values (0 to 100) and for tetas values (0 to 80)
516	! for 2 scattering phase functions: liquid spherical, ice non spherical
517
518	IF ( nrefl.GT. ntetas ) THEN
519	PRINT *,'Error in lidar_simulator, nrefl should be less then ',ntetas,' not',nrefl
520	STOP
521	ENDIF
522
523	rlumA_mod=0
524	rlumB_mod=0
525	!
526	pi = ACOS(-1.0)
527	r_norm(:)=1./ COS(pi/180.*tetas(:))
528	!
529	tautot_S_liq(:)=MAX(tautot_S_liq(:),tau(1))
530	tautot_S_ice(:)=MAX(tautot_S_ice(:),tau(1))
531	tautot_S(:) = tautot_S_ice(:) + tautot_S_liq(:)
532	!
533	! relative fraction of the opt. thick due to liquid or ice clouds
534	WHERE (tautot_S(:) .GT. 0.)
535	frac_taucol_liq(:) = tautot_S_liq(:) / tautot_S(:)
536	frac_taucol_ice(:) = tautot_S_ice(:) / tautot_S(:)
537	ELSEWHERE
538	frac_taucol_liq(:) = 1.
539	frac_taucol_ice(:) = 0.
540	END WHERE
541	tautot_S(:)=MIN(tautot_S(:),tau(nbtau))
542	!
543	! Linear interpolation :
544
545	DO ny=1,nbtau-1
546	! microphysics A (liquid clouds)
547	aA(:,ny) = (rlumA(:,ny+1)-rlumA(:,ny))/(tau(ny+1)-tau(ny))
548	bA(:,ny) = rlumA(:,ny) - aA(:,ny)*tau(ny)
549	! microphysics B (ice clouds)
550	aB(:,ny) = (rlumB(:,ny+1)-rlumB(:,ny))/(tau(ny+1)-tau(ny))
551	bB(:,ny) = rlumB(:,ny) - aB(:,ny)*tau(ny)
552	ENDDO
553	!
554	DO it=1,ntetas
555	DO ny=1,nbtau-1
556	WHERE (tautot_S(:).GE.tau(ny).AND.tautot_S(:).LE.tau(ny+1))
557	rlumA_mod(:,it) = aA(it,ny)*tautot_S(:) + bA(it,ny)
558	rlumB_mod(:,it) = aB(it,ny)*tautot_S(:) + bB(it,ny)
559	END WHERE
560	END DO
561	END DO
562	!
563	DO it=1,ntetas
564	refl(:,it) = frac_taucol_liq(:) * rlumA_mod(:,it) &
565	+ frac_taucol_ice(:) * rlumB_mod(:,it)
566	! normalized radiance -> reflectance:
567	refl(:,it) = refl(:,it) * r_norm(it)
568	ENDDO
569
570	RETURN
571	END SUBROUTINE parasol

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format