Context Navigation

lidar_simulator.F90 @ 6

Last change on this file since 6 was 1, checked in by emillour, 14 years ago
Import initial LMDZ5
File size: 21.7 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
182	!------------------------------------------------------------
183	!---- 1. Preliminary definitions and calculations :
184	!------------------------------------------------------------
185
186	if ( npart .ne. 4 ) then
187	print *,'Error in lidar_simulator, npart should be 4, not',npart
188	stop
189	endif
190
191	pi = dacos(-1.D0)
192
193	! Polynomial coefficients for spherical liq/ice particles derived from Mie theory.
194	! Polynomial coefficients for non spherical particles derived from a composite of
195	! Ray-tracing theory for large particles (e.g. Noel et al., Appl. Opt., 2001)
196	! and FDTD theory for very small particles (Yang et al., JQSRT, 2003).
197
198	! We repeat the same coefficients for LS and CONV cloud to make code more readable
199	!* LS Liquid water coefficients:
200	polpart(INDX_LSLIQ,1) = 2.6980e-8
201	polpart(INDX_LSLIQ,2) = -3.7701e-6
202	polpart(INDX_LSLIQ,3) = 1.6594e-4
203	polpart(INDX_LSLIQ,4) = -0.0024
204	polpart(INDX_LSLIQ,5) = 0.0626
205	!* LS Ice coefficients:
206	if (ice_type.eq.0) then
207	polpart(INDX_LSICE,1) = -1.0176e-8
208	polpart(INDX_LSICE,2) = 1.7615e-6
209	polpart(INDX_LSICE,3) = -1.0480e-4
210	polpart(INDX_LSICE,4) = 0.0019
211	polpart(INDX_LSICE,5) = 0.0460
212	endif
213	!* LS Ice NS coefficients:
214	if (ice_type.eq.1) then
215	polpart(INDX_LSICE,1) = 1.3615e-8
216	polpart(INDX_LSICE,2) = -2.04206e-6
217	polpart(INDX_LSICE,3) = 7.51799e-5
218	polpart(INDX_LSICE,4) = 0.00078213
219	polpart(INDX_LSICE,5) = 0.0182131
220	endif
221	!* CONV Liquid water coefficients:
222	polpart(INDX_CVLIQ,1) = 2.6980e-8
223	polpart(INDX_CVLIQ,2) = -3.7701e-6
224	polpart(INDX_CVLIQ,3) = 1.6594e-4
225	polpart(INDX_CVLIQ,4) = -0.0024
226	polpart(INDX_CVLIQ,5) = 0.0626
227	!* CONV Ice coefficients:
228	if (ice_type.eq.0) then
229	polpart(INDX_CVICE,1) = -1.0176e-8
230	polpart(INDX_CVICE,2) = 1.7615e-6
231	polpart(INDX_CVICE,3) = -1.0480e-4
232	polpart(INDX_CVICE,4) = 0.0019
233	polpart(INDX_CVICE,5) = 0.0460
234	endif
235	if (ice_type.eq.1) then
236	polpart(INDX_CVICE,1) = 1.3615e-8
237	polpart(INDX_CVICE,2) = -2.04206e-6
238	polpart(INDX_CVICE,3) = 7.51799e-5
239	polpart(INDX_CVICE,4) = 0.00078213
240	polpart(INDX_CVICE,5) = 0.0182131
241	endif
242
243	! density:
244	!* clear-sky air:
245	rhoair = pres/(287.04*temp)
246
247	!* liquid/ice particules:
248	rhopart(INDX_LSLIQ) = rholiq
249	rhopart(INDX_LSICE) = rhoice
250	rhopart(INDX_CVLIQ) = rholiq
251	rhopart(INDX_CVICE) = rhoice
252
253	! effective radius particles:
254	rad_part(:,:,INDX_LSLIQ) = ls_radliq(:,:)
255	rad_part(:,:,INDX_LSICE) = ls_radice(:,:)
256	rad_part(:,:,INDX_CVLIQ) = cv_radliq(:,:)
257	rad_part(:,:,INDX_CVICE) = cv_radice(:,:)
258	rad_part(:,:,:)=MAX(rad_part(:,:,:),0.)
259	rad_part(:,:,:)=MIN(rad_part(:,:,:),70.0e-6)
260
261	! altitude at half pressure levels:
262	zheight(:,1) = 0.0
263	do k = 2, nlev+1
264	zheight(:,k) = zheight(:,k-1) &
265	-(presf(:,k)-presf(:,k-1))/(rhoair(:,k-1)*9.81)
266	enddo
267
268	! cloud fraction (0 or 1) in each sub-column:
269	! (if frac_out=1or2 -> frac_sub=1; if frac_out=0 -> frac_sub=0)
270	frac_sub = MIN( frac_out, 1.0 )
271
272	!------------------------------------------------------------
273	!---- 2. Molecular alpha and beta:
274	!------------------------------------------------------------
275
276	beta_mol = pres/km/temp * Cmol
277	alpha_mol = 8.0pi/3.0 beta_mol
278
279	!------------------------------------------------------------
280	!---- 3. Particles alpha and beta:
281	!------------------------------------------------------------
282
283	! polynomes kp_lidar derived from Mie theory:
284	do i = 1, npart
285	where ( rad_part(:,:,i).gt.0.0)
286	kp_part(:,:,i) = &
287	polpart(i,1)(rad_part(:,:,i)1e6)**4 &
288	+ polpart(i,2)(rad_part(:,:,i)1e6)**3 &
289	+ polpart(i,3)(rad_part(:,:,i)1e6)**2 &
290	+ polpart(i,4)(rad_part(:,:,i)1e6) &
291	+ polpart(i,5)
292	elsewhere
293	kp_part(:,:,i) = 0.
294	endwhere
295	enddo
296
297	! mixing ratio particules in each subcolumn:
298	qpart(:,:,INDX_LSLIQ) = q_lsliq(:,:) ! oct08
299	qpart(:,:,INDX_LSICE) = q_lsice(:,:) ! oct08
300	qpart(:,:,INDX_CVLIQ) = q_cvliq(:,:) ! oct08
301	qpart(:,:,INDX_CVICE) = q_cvice(:,:) ! oct08
302
303	! alpha of particles in each subcolumn:
304	do i = 1, npart
305	where ( rad_part(:,:,i).gt.0.0)
306	alpha_part(:,:,i) = 3.0/4.0 * Qscat &
307	* rhoair(:,:) * qpart(:,:,i) &
308	/ (rhopart(i) * rad_part(:,:,i) )
309	elsewhere
310	alpha_part(:,:,i) = 0.
311	endwhere
312	enddo
313
314	!------------------------------------------------------------
315	!---- 4. Backscatter signal:
316	!------------------------------------------------------------
317
318	! optical thickness (molecular):
319	! opt. thick of each layer
320	tau_mol(:,1:nlev) = alpha_mol(:,1:nlev) &
321	& *(zheight(:,2:nlev+1)-zheight(:,1:nlev))
322	! opt. thick from TOA
323	DO k = nlev-1, 1, -1
324	tau_mol(:,k) = tau_mol(:,k) + tau_mol(:,k+1)
325	ENDDO
326
327	! optical thickness (particles):
328
329	tau_part = rdiffm * alpha_part
330	DO i = 1, npart
331	! opt. thick of each layer
332	tau_part(:,:,i) = tau_part(:,:,i) &
333	& * (zheight(:,2:nlev+1)-zheight(:,1:nlev) )
334	! opt. thick from TOA
335	DO k = nlev-1, 1, -1
336	tau_part(:,k,i) = tau_part(:,k,i) + tau_part(:,k+1,i)
337	ENDDO
338	ENDDO
339
340	! molecular signal:
341	! Upper layer
342	pmol(:,nlev) = beta_mol(:,nlev) / (2.*tau_mol(:,nlev)) &
343	& * (1.-exp(-2.0*tau_mol(:,nlev)))
344	! Other layers
345	DO k= nlev-1, 1, -1
346	tau_mol_lay(:) = tau_mol(:,k)-tau_mol(:,k+1) ! opt. thick. of layer k
347	WHERE (tau_mol_lay(:).GT.0.)
348	pmol(:,k) = beta_mol(:,k) * EXP(-2.0tau_mol(:,k+1)) / (2.tau_mol_lay(:)) &
349	& * (1.-exp(-2.0*tau_mol_lay(:)))
350	ELSEWHERE
351	! This must never happend, but just in case, to avoid div. by 0
352	pmol(:,k) = beta_mol(:,k) * EXP(-2.0*tau_mol(:,k+1))
353	END WHERE
354	END DO
355	!
356	! Total signal (molecular + particules):
357	!
358	! For performance reason on vector computers, the 2 following lines should not be used
359	! and should be replace by the later one.
360	! betatot(:,:) = beta_mol(:,:) + sum(kp_part*alpha_part,dim=3)
361	! tautot(:,:) = tau_mol(:,:) + sum(tau_part,dim=3)
362	betatot(:,:) = beta_mol(:,:)
363	tautot(:,:) = tau_mol(:,:)
364	DO i = 1, npart
365	betatot(:,:) = betatot(:,:) + kp_part(:,:,i)*alpha_part(:,:,i)
366	tautot(:,:) = tautot(:,:) + tau_part(:,:,i)
367	ENDDO ! i
368	!
369	! Upper layer
370	pnorm(:,nlev) = betatot(:,nlev) / (2.*tautot(:,nlev)) &
371	& * (1.-exp(-2.0*tautot(:,nlev)))
372	! Other layers
373	DO k= nlev-1, 1, -1
374	tautot_lay(:) = tautot(:,k)-tautot(:,k+1) ! optical thickness of layer k
375	WHERE (tautot_lay(:).GT.0.)
376	pnorm(:,k) = betatot(:,k) * EXP(-2.0tautot(:,k+1)) / (2.tautot_lay(:)) &
377	!correc pnorm(:,k) = betatot(:,k) * EXP(-2.0*tautot(:,k+1)) & ! correc Satoh
378	!correc & / (2.0*tautot_lay(:)) & ! correc Satoh
379	& * (1.-EXP(-2.0*tautot_lay(:)))
380	ELSEWHERE
381	! This must never happend, but just in case, to avoid div. by 0
382	pnorm(:,k) = betatot(:,k) * EXP(-2.0*tautot(:,k+1))
383	END WHERE
384	END DO
385
386	!-------- End computation Lidar --------------------------
387
388	!---------------------------------------------------------
389	! Parasol/Polder module
390	!
391	! Purpose : Compute reflectance for one particular viewing direction
392	! and 5 solar zenith angles (calculation valid only over ocean)
393	! ---------------------------------------------------------
394
395	! initialization:
396	refl(:,:) = 0.0
397
398	! activate parasol calculations:
399	if (ok_parasol) then
400
401	! Optical thickness from TOA to surface
402	tautot_S_liq = 0.
403	tautot_S_ice = 0.
404	tautot_S_liq(:) = tautot_S_liq(:) &
405	+ tau_part(:,1,1) + tau_part(:,1,3)
406	tautot_S_ice(:) = tautot_S_ice(:) &
407	+ tau_part(:,1,2) + tau_part(:,1,4)
408
409	call parasol(npoints,nrefl,undef &
410	,tautot_S_liq,tautot_S_ice &
411	,refl)
412
413	endif ! ok_parasol
414
415	END SUBROUTINE lidar_simulator
416	!
417	!---------------------------------------------------------------------------------
418	!
419	SUBROUTINE parasol(npoints,nrefl,undef &
420	,tautot_S_liq,tautot_S_ice &
421	,refl)
422	!---------------------------------------------------------------------------------
423	! Purpose: To compute Parasol reflectance signal from model-simulated profiles
424	! of cloud water and cloud fraction in each sub-column of each model
425	! gridbox.
426	!
427	!
428	! December 2008, S. Bony, H. Chepfer and J-L. Dufresne :
429	! - optimization for vectorization
430	!
431	! Version 2.0 (October 2008)
432	! Version 2.1 (December 2008)
433	!---------------------------------------------------------------------------------
434
435	IMPLICIT NONE
436
437	! inputs
438	INTEGER npoints ! Number of horizontal gridpoints
439	INTEGER nrefl ! Number of angles for which the reflectance
440	! is computed. Can not be greater then ntetas
441	REAL undef ! Undefined value. Currently not used
442	REAL tautot_S_liq(npoints) ! liquid water cloud optical thickness,
443	! integrated from TOA to surface
444	REAL tautot_S_ice(npoints) ! same for ice water clouds only
445	! outputs
446	REAL refl(npoints,nrefl) ! Parasol reflectances
447	!
448	! Local variables
449	REAL tautot_S(npoints) ! cloud optical thickness, from TOA to surface
450	REAL frac_taucol_liq(npoints), frac_taucol_ice(npoints)
451
452	REAL pi
453	! look up table variables:
454	INTEGER ny, it
455	INTEGER ntetas, nbtau ! number of angle and of optical thickness
456	! of the look-up table
457	PARAMETER (ntetas=5, nbtau=7)
458	REAL aa(ntetas,nbtau-1), ab(ntetas,nbtau-1)
459	REAL ba(ntetas,nbtau-1), bb(ntetas,nbtau-1)
460	REAL tetas(ntetas),tau(nbtau)
461	REAL r_norm(ntetas)
462	REAL rlumA(ntetas,nbtau), rlumB(ntetas,nbtau)
463	REAL rlumA_mod(npoints,5), rlumB_mod(npoints,5)
464
465	DATA tau /0., 1., 5., 10., 20., 50., 100./
466	DATA tetas /0., 20., 40., 60., 80./
467
468	! Look-up table for spherical liquid particles:
469	DATA (rlumA(1,ny),ny=1,nbtau) /0.03, 0.090886, 0.283965, &
470	0.480587, 0.695235, 0.908229, 1.0 /
471	DATA (rlumA(2,ny),ny=1,nbtau) /0.03, 0.072185, 0.252596, &
472	0.436401, 0.631352, 0.823924, 0.909013 /
473	DATA (rlumA(3,ny),ny=1,nbtau) /0.03, 0.058410, 0.224707, &
474	0.367451, 0.509180, 0.648152, 0.709554 /
475	DATA (rlumA(4,ny),ny=1,nbtau) /0.03, 0.052498, 0.175844, &
476	0.252916, 0.326551, 0.398581, 0.430405 /
477	DATA (rlumA(5,ny),ny=1,nbtau) /0.03, 0.034730, 0.064488, &
478	0.081667, 0.098215, 0.114411, 0.121567 /
479
480	! Look-up table for ice particles:
481	DATA (rlumB(1,ny),ny=1,nbtau) /0.03, 0.092170, 0.311941, &
482	0.511298, 0.712079 , 0.898243 , 0.976646 /
483	DATA (rlumB(2,ny),ny=1,nbtau) /0.03, 0.087082, 0.304293, &
484	0.490879, 0.673565, 0.842026, 0.912966 /
485	DATA (rlumB(3,ny),ny=1,nbtau) /0.03, 0.083325, 0.285193, &
486	0.430266, 0.563747, 0.685773, 0.737154 /
487	DATA (rlumB(4,ny),ny=1,nbtau) /0.03, 0.084935, 0.233450, &
488	0.312280, 0.382376, 0.446371, 0.473317 /
489	DATA (rlumB(5,ny),ny=1,nbtau) /0.03, 0.054157, 0.089911, &
490	0.107854, 0.124127, 0.139004, 0.145269 /
491
492	!--------------------------------------------------------------------------------
493	! Lum_norm=f(tetaS,tau_cloud) derived from adding-doubling calculations
494	! valid ONLY ABOVE OCEAN (albedo_sfce=5%)
495	! valid only in one viewing direction (theta_v=30�, phi_s-phi_v=320�)
496	! based on adding-doubling radiative transfer computation
497	! for tau values (0 to 100) and for tetas values (0 to 80)
498	! for 2 scattering phase functions: liquid spherical, ice non spherical
499
500	IF ( nrefl.GT. ntetas ) THEN
501	PRINT *,'Error in lidar_simulator, nrefl should be less then ',ntetas,' not',nrefl
502	STOP
503	ENDIF
504
505	rlumA_mod=0
506	rlumB_mod=0
507	!
508	pi = ACOS(-1.0)
509	r_norm(:)=1./ COS(pi/180.*tetas(:))
510	!
511	tautot_S_liq(:)=MAX(tautot_S_liq(:),tau(1))
512	tautot_S_ice(:)=MAX(tautot_S_ice(:),tau(1))
513	tautot_S(:) = tautot_S_ice(:) + tautot_S_liq(:)
514	!
515	! relative fraction of the opt. thick due to liquid or ice clouds
516	WHERE (tautot_S(:) .GT. 0.)
517	frac_taucol_liq(:) = tautot_S_liq(:) / tautot_S(:)
518	frac_taucol_ice(:) = tautot_S_ice(:) / tautot_S(:)
519	ELSEWHERE
520	frac_taucol_liq(:) = 1.
521	frac_taucol_ice(:) = 0.
522	END WHERE
523	tautot_S(:)=MIN(tautot_S(:),tau(nbtau))
524	!
525	! Linear interpolation :
526
527	DO ny=1,nbtau-1
528	! microphysics A (liquid clouds)
529	aA(:,ny) = (rlumA(:,ny+1)-rlumA(:,ny))/(tau(ny+1)-tau(ny))
530	bA(:,ny) = rlumA(:,ny) - aA(:,ny)*tau(ny)
531	! microphysics B (ice clouds)
532	aB(:,ny) = (rlumB(:,ny+1)-rlumB(:,ny))/(tau(ny+1)-tau(ny))
533	bB(:,ny) = rlumB(:,ny) - aB(:,ny)*tau(ny)
534	ENDDO
535	!
536	DO it=1,ntetas
537	DO ny=1,nbtau-1
538	WHERE (tautot_S(:).GE.tau(ny).AND.tautot_S(:).LE.tau(ny+1))
539	rlumA_mod(:,it) = aA(it,ny)*tautot_S(:) + bA(it,ny)
540	rlumB_mod(:,it) = aB(it,ny)*tautot_S(:) + bB(it,ny)
541	END WHERE
542	END DO
543	END DO
544	!
545	DO it=1,ntetas
546	refl(:,it) = frac_taucol_liq(:) * rlumA_mod(:,it) &
547	+ frac_taucol_ice(:) * rlumB_mod(:,it)
548	! normalized radiance -> reflectance:
549	refl(:,it) = refl(:,it) * r_norm(it)
550	ENDDO
551
552	RETURN
553	END SUBROUTINE parasol

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

Download in other formats:

Original Format