source: LMDZ5/trunk/libf/cosp/lidar_simulator.F90 @ 2114

Last change on this file since 2114 was 1907, checked in by lguez, 11 years ago

Added a copyright property to every file of the distribution, except
for the fcm files (which have their own copyright). Use svn propget on
a file to see the copyright. For instance:

$ svn propget copyright libf/phylmd/physiq.F90
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

Also added the files defining the CeCILL version 2 license, in French
and English, at the top of the LMDZ tree.

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