source: LMDZ4/branches/LMDZ4-dev-20091210/libf/cosp/lidar_simulator.F90 @ 2570

Last change on this file since 2570 was 1262, checked in by idelkadi, 15 years ago

-Rajout des fichiers cosp_input_nl.txt et cosp_output_nl.txt pour controler les parametres entrees sorties COSP
-Rajout du repertoire cosp contenant les sources de COSP

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.0*pi/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.0*tau_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.0*tautot(:,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.