source: TOOLS/CMIP6_FORCING/AER_OPTICS/mie_6bands_sulfate.f @ 3385

Last change on this file since 3385 was 3385, checked in by oboucher, 4 years ago

Uploading package to prepare tropospheric aerosol optics LUT in LMDz fo rCMIP6

File size: 17.2 KB
Line 
1        PROGRAM mie
2        IMPLICIT NONE
3C
4C-------Mie computations for a size distribution 
5C       of homogeneous spheres.
6c
7C==========================================================
8C--Ref : Toon and Ackerman, Applied Optics, 1981
9C        Stephens, CSIRO, 1979
10C Attention : surdimensionement des tableaux
11C to be compiled with double precision option (-r8 on Sun)
12C AUTHOR: Olivier Boucher
13C-------SIZE distribution properties----------------
14C--sigma_g : geometric standard deviation 
15C--r_0     : geometric number mean radius (um)/modal radius
16C--Ntot    : total concentration in m-3 
17c
18        REAL rmin, rmax    !----integral bounds in  m
19c
20c--Nmode= 4 modes for SS in INCA
21c-- AS, CS, CI, SS
22c--NDis=1 distribution per mode
23        INTEGER Nmode, Ndis, mode, dis
24        PARAMETER (Nmode=2, Ndis=1)
25        REAL sigma_g(Ndis,Nmode), r_0(Ndis,Nmode), Ntot(Ndis,Nmode)
26c--sulfate Coarse Soluble (CS) & Accumulation Soluble (AS)
27c--becareful to list in the correct order if more than 1 dis per mode
28        DATA r_0    /0.433E-6,0.1E-6/   !--meters
29        DATA sigma_g/2.0,1.8/
30        DATA Ntot   /1.0,1.0/
31        CHARACTER*33 chmode(Nmode)
32        DATA chmode/'Sulfate Coarse Soluble (CS)', 
33     .              'Sulfate Accumulation Soluble (AS)'/
34c
35        REAL masse,volume,surface,rho
36        PARAMETER (rho=1.769E3)  !--dry density kg/m3 Tang and Munkelwitz
37c
38c---------- RH growth parameters----------------
39c
40        INTEGER Nrh,irh
41        PARAMETER(Nrh=12)
42        REAL RH_tab(Nrh),RH
43        DATA RH_tab/0.,10.,20,30.,40.,50.,60.,70.,80.,85.,90.,95./
44        REAL rh_growth(Nrh)
45c
46c--growth factor normalised to 0% relative humidity for sulfate
47        DATA rh_growth/1.000, 1.000, 1.000, 1.000, 1.169, 1.220, 
48     .                 1.282, 1.363, 1.485, 1.580, 1.735, 2.085/
49c
50c-------------------------------------
51c
52        COMPLEX m          !----refractive index m=n_r-i*n_i
53        INTEGER Nmax,Nstart !--number of iterations
54        COMPLEX k2, k3, z1, z2
55        COMPLEX u1,u5,u6,u8
56        COMPLEX a(1:21000), b(1:21000)
57        COMPLEX I 
58        INTEGER n  !--loop index
59        REAL pi, nnn
60        COMPLEX nn
61        REAL Q_ext, Q_abs, Q_sca, g, omega   !--parameters for radius r
62        REAL x  !--size parameter
63        REAL r  !--radius
64        REAL sigma_sca, sigma_ext, sigma_abs
65        REAL omegatot,  gtot !--averaged parameters
66        COMPLEX ksiz2(-1:21000), psiz2(1:21000)
67        COMPLEX nu1z1(1:21010), nu1z2(1:21010)
68        COMPLEX nu3z2(0:21000)
69        REAL number, deltar
70        INTEGER bin, Nbin, k
71        PARAMETER (Nbin=700)
72c
73C---wavelengths STREAMER
74        INTEGER Nwv, NwvmaxSW
75        PARAMETER (NwvmaxSW=25)
76        REAL lambda(1:NwvmaxSW)
77        DATA lambda/0.28E-6, 0.30E-6, 0.33E-6, 0.36E-6, 0.40E-6,
78     .              0.44E-6, 0.48E-6, 0.52E-6, 0.57E-6, 0.64E-6,
79     .              0.69E-6, 0.75E-6, 0.78E-6, 0.87E-6, 1.00E-6,
80     .              1.10E-6, 1.19E-6, 1.28E-6, 1.53E-6, 1.64E-6,
81     .              2.13E-6, 2.38E-6, 2.91E-6, 3.42E-6, 4.00E-6/
82c
83        INTEGER nb, nb_lambda
84        PARAMETER (nb_lambda=5)
85        REAL lambda_ref(nb_lambda)
86        DATA lambda_ref /0.443E-6,0.550E-6,0.670E-6,0.765E-6,0.865E-6/ 
87c
88c--read refractive index from old file with different wavelengths 
89c
90        REAL n_r_tab(1:Nrh,1:NwvmaxSW-1+nb_lambda)
91        REAL n_i_tab(1:Nrh,1:NwvmaxSW-1+nb_lambda)
92c
93C---TOA fluxes - Streamer Cs
94        REAL weight(1:NwvmaxSW-1)
95c        DATA weight/0.839920E1, 0.231208E2, 0.322393E2, 0.465058E2,
96c     .              0.678199E2, 0.798964E2, 0.771359E2, 0.888472E2,
97c     .              0.115281E3, 0.727565E2, 0.816992E2, 0.336172E2,
98c     .              0.914603E2, 0.112706E3, 0.658840E2, 0.524470E2,
99c     .              0.391067E2, 0.883864E2, 0.276672E2, 0.681812E2,
100c     .              0.190966E2, 0.250766E2, 0.128704E2, 0.698720E1/
101C---TOA fluxes - Tad
102c        DATA weight/ 4.20, 11.56, 16.12, 23.25, 33.91, 39.95, 38.57,
103c     .              44.42, 57.64, 29.36, 47.87, 16.81, 45.74, 56.35,
104c     .              32.94, 26.22, 19.55, 44.19, 13.83, 34.09,  9.55,
105c      .              12.54,  6.44,  3.49/
106C---BOA fluxes - Tad
107        DATA weight/ 0.01,  4.05, 9.51,  15.99, 26.07, 33.10, 33.07,
108     .              39.91, 52.67, 27.89, 43.60, 13.67, 42.22, 40.12,
109     .              32.70, 14.44, 19.48, 14.23, 13.43, 16.42,  8.33,
110     .               0.95,  0.65, 2.76/
111c
112        REAL lambda_int(1:NwvmaxSW-1+nb_lambda)
113c
114        REAL final_a(1:NwvmaxSW-1+nb_lambda)
115        REAL final_g(1:NwvmaxSW-1+nb_lambda)
116        REAL final_w(1:NwvmaxSW-1+nb_lambda)
117        REAL final_qext(1:NwvmaxSW-1+nb_lambda)
118c
119        INTEGER band, NbandSW, NbandLW
120        PARAMETER (NbandSW=6, NbandLW=5)
121c
122        REAL gcm_a(NbandSW+NbandLW)
123        REAL gcm_g(NbandSW+NbandLW)
124        REAL gcm_w(NbandSW+NbandLW)
125        REAL gcm_qext(NbandSW+NbandLW)
126        REAL gcm_weight_a(NbandSW+NbandLW)
127        REAL gcm_weight_g(NbandSW+NbandLW)
128        REAL gcm_weight_w(NbandSW+NbandLW)
129        REAL gcm_weight_qext(NbandSW+NbandLW)
130c
131        REAL ss_a(NbandSW+NbandLW+nb_lambda,Nrh)
132        REAL ss_w(NbandSW+NbandLW+nb_lambda,Nrh)
133        REAL ss_g(NbandSW+NbandLW+nb_lambda,Nrh)
134        REAL ss_qext(NbandSW+NbandLW+nb_lambda,Nrh)
135c
136        INTEGER NwvmaxLW
137        PARAMETER (NwvmaxLW=100)
138        REAL Tb, Planck
139        PARAMETER (Tb=260.0)
140c
141        INTEGER wv, nb_wv
142        PARAMETER (nb_wv=100)
143        REAL wv_SUL(1:nb_wv)
144        REAL index_r_SUL(1:Nrh,1:nb_wv)
145        REAL index_i_SUL(1:Nrh,1:nb_wv)
146        REAL rh_dummy
147        REAL count_n_r, count_n_i
148c
149c------opening output files
150c
151        OPEN (unit=14, file='SEXT_sulfate_6bands.txt')
152        OPEN (unit=15, file='G_sulfate_6bands.txt')
153        OPEN (unit=16, file='SSA_sulfate_6bands.txt')
154        OPEN (unit=17, file='QEXT_sulfate_6bands.txt')
155
156        OPEN (unit=24, file='SEXT_sulfate_5wave.txt')
157        OPEN (unit=25, file='QEXT_sulfate_5wave.txt')
158        OPEN (unit=26, file='SABS_sulfate_5wave.txt')
159c
160c--initializing wavelengths for calculations
161c
162        DO Nwv=1, NwvmaxSW-1
163        lambda_int(Nwv)=( lambda(Nwv)+lambda(Nwv+1) ) /2.
164        ENDDO
165c
166        DO nb=1, nb_lambda
167        lambda_int(NwvmaxSW-1+nb)=lambda_ref(nb)
168        ENDDO
169c
170c--lecture des indices from a high-resolution spectral file
171c-------RH dependent RI from file
172c
173       OPEN(unit=20,file='ri_sul_v2')
174c
175       DO irh=1,Nrh
176        DO wv=1, nb_wv
177          READ (20,*) rh_dummy, wv_SUL(wv),
178     .                index_r_SUL(irh,wv), index_i_SUL(irh,wv)
179        ENDDO
180       ENDDO
181c
182       CLOSE(20)
183c
184c---interpolate refractive index to tab values
185c--take average or closest neighbour wavelength
186c
187        DO irh=1, Nrh
188c
189        DO Nwv=1, NwvmaxSW-1+nb_lambda
190           n_r_tab(irh,Nwv)=0.0
191           n_i_tab(irh,Nwv)=0.0
192        ENDDO
193c
194c--interpolating on our wavelengths
195c
196        DO Nwv=1, NwvmaxSW-1
197c
198c--first real part
199          count_n_r=0.0
200          DO wv=1, nb_wv
201            IF (wv_SUL(wv).GT.lambda(Nwv).AND.
202     .          wv_SUL(wv).LT.lambda(Nwv+1)) THEN
203              n_r_tab(irh,Nwv)=n_r_tab(irh,Nwv)+index_r_SUL(irh,wv)
204              count_n_r=count_n_r+1.0
205            ENDIF 
206          ENDDO
207c
208          IF (count_n_r.GT.0.5) THEN
209c--averaging
210            n_r_tab(irh,Nwv)=n_r_tab(irh,Nwv)/count_n_r
211          ELSE
212c--otherwise closest neighbour
213          DO wv=1, nb_wv
214            IF (wv_SUL(wv).LT.lambda_int(Nwv)) THEN
215              n_r_tab(irh,Nwv)=index_r_SUL(irh,wv)
216            ENDIF 
217          ENDDO
218          ENDIF
219c
220c--then imaginary part
221          count_n_i=0.0
222          DO wv=1, nb_wv
223            IF (wv_SUL(wv).GT.lambda(Nwv).AND.
224     .          wv_SUL(wv).LT.lambda(Nwv+1)) THEN
225              n_i_tab(irh,Nwv)=n_i_tab(irh,Nwv)+index_i_SUL(irh,wv)
226              count_n_i=count_n_i+1.0
227            ENDIF 
228          ENDDO
229c
230          IF (count_n_i.GT.0.5) THEN
231c--averaging
232            n_i_tab(irh,Nwv)=n_i_tab(irh,Nwv)/count_n_i
233          ELSE
234c--otherwise closest neighbour
235          DO wv=1, nb_wv
236            IF (wv_SUL(wv).LT.lambda_int(Nwv)) THEN
237              n_i_tab(irh,Nwv)=index_i_SUL(irh,wv)
238            ENDIF 
239          ENDDO
240          ENDIF
241c
242        ENDDO !--wv
243c
244        ENDDO !--irh
245c
246c-----------------------------------------------------------
247c
248c--now defining nr and ni for my set of reference wavelengths for SUL
249        DO nb=1, nb_lambda
250c
251          DO irh=1,Nrh
252          DO wv=1, nb_wv
253            IF (wv_SUL(wv).LT.lambda_ref(nb)) THEN
254              n_r_tab(irh,NwvmaxSW-1+nb)=index_r_SUL(irh,wv)
255              n_i_tab(irh,NwvmaxSW-1+nb)=index_i_SUL(irh,wv)
256            ENDIF 
257          ENDDO
258          ENDDO
259c
260        ENDDO !--nb
261
262c        OPEN(unit=33,file='output_refr_index_sulfate_for_check.dat') 
263c        DO Nwv=1, NwvmaxSW-1+nb_lambda
264c          WRITE(33,*) lambda_int(Nwv), n_r_tab(1,Nwv), n_i_tab(1,Nwv),
265c     .                                 n_r_tab(12,Nwv), n_i_tab(12,Nwv)
266c        ENDDO
267c        CLOSE(33)
268         
269c
270c---now start calculations 
271c
272        DO mode=1, Nmode
273c
274        DO irh=1,Nrh       !---------LOOP OVER RH
275c
276c-map extra wavelengths to those from input file 
277c
278        rmin=0.002E-6*rh_growth(irh)
279        rmax=30.E-6*rh_growth(irh)
280c
281        DO Nwv=1, NwvmaxSW-1+nb_lambda
282c
283        m=CMPLX(n_r_tab(irh,Nwv),-n_i_tab(irh,Nwv))
284c
285        pi=4.*ATAN(1.)
286c
287        I=CMPLX(0.,1.)
288c
289          sigma_sca=0.0
290          sigma_ext=0.0
291          sigma_abs=0.0
292          gtot=0.0
293          omegatot=0.0
294          masse = 0.0
295          volume=0.0
296          surface=0.0
297c
298        DO bin=0, Nbin !---loop on size bins
299 
300        r=exp(log(rmin)+FLOAT(bin)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
301        x=2.*pi*r/lambda_int(Nwv)
302        deltar=1./FLOAT(Nbin)*(log(rmax)-log(rmin))
303c
304        number=0
305        DO dis=1, Ndis
306          number=number+
307     .       Ntot(dis,mode)/SQRT(2.*pi)/log(sigma_g(dis,mode))*
308     .       exp(-0.5*(log(r/(r_0(dis,mode)*rh_growth(irh)))/
309     .               log(sigma_g(dis,mode)))**2)
310        ENDDO
311c--80% RH mass
312c        masse=masse  +4./3.*pi*
313c     .          ((r/rh_growth(irh)*rh_growth(9))**3)*
314c     .           number*deltar*rho*1.E3          !--g/m3
315c--dry aerosol mass
316        masse=masse  +4./3.*pi*
317     .          ((r/rh_growth(irh))**3)*
318     .           number*deltar*rho*1.E3          !--g/m3
319c
320        volume=volume+4./3.*pi*(r**3)*number*deltar
321        surface=surface+4.*pi*r**2*number*deltar
322c
323        k2=m
324        k3=CMPLX(1.0,0.0)
325
326        z2=CMPLX(x,0.0)
327        z1=m*z2
328 
329        IF (0.0.LE.x.AND.x.LE.8.) THEN
330        Nmax=INT(x+4*x**(1./3.)+1.)+2
331        ELSEIF (8..LT.x.AND.x.LT.4200.) THEN
332        Nmax=INT(x+4.05*x**(1./3.)+2.)+1
333        ELSEIF (4200..LE.x.AND.x.LE.20000.) THEN
334        Nmax=INT(x+4*x**(1./3.)+2.)+1
335        ELSE
336        WRITE(10,*) 'x out of bound, x=', x
337        STOP
338        ENDIF
339
340        Nstart=Nmax+10
341
342C-----------loop for nu1z1, nu1z2
343
344       nu1z1(Nstart)=CMPLX(0.0,0.0)
345       nu1z2(Nstart)=CMPLX(0.0,0.0)
346       DO n=Nstart-1, 1 , -1
347       nn=CMPLX(FLOAT(n),0.0)
348       nu1z1(n)=(nn+1.)/z1 - 1./( (nn+1.)/z1 + nu1z1(n+1) ) 
349       nu1z2(n)=(nn+1.)/z2 - 1./( (nn+1.)/z2 + nu1z2(n+1) )
350       ENDDO
351
352C------------loop for nu3z2
353
354       nu3z2(0)=-I
355       DO n=1, Nmax
356       nn=CMPLX(FLOAT(n),0.0)
357       nu3z2(n)=-nn/z2 + 1./ (nn/z2 - nu3z2(n-1) ) 
358       ENDDO
359
360C-----------loop for psiz2 and ksiz2 (z2)
361       ksiz2(-1)=COS(REAL(z2))-I*SIN(REAL(z2))
362       ksiz2(0)=SIN(REAL(z2))+I*COS(REAL(z2))
363       DO n=1,Nmax
364       nn=CMPLX(FLOAT(n),0.0)
365       ksiz2(n)=(2.*nn-1.)/z2 * ksiz2(n-1) - ksiz2(n-2)
366       psiz2(n)=CMPLX(REAL(ksiz2(n)),0.0)
367       ENDDO
368
369C-----------loop for a(n) and b(n)
370   
371       DO n=1, Nmax
372        u1=k3*nu1z1(n) - k2*nu1z2(n)
373        u5=k3*nu1z1(n) - k2*nu3z2(n)
374        u6=k2*nu1z1(n) - k3*nu1z2(n)
375        u8=k2*nu1z1(n) - k3*nu3z2(n)
376        a(n)=psiz2(n)/ksiz2(n) * u1/u5
377        b(n)=psiz2(n)/ksiz2(n) * u6/u8
378       ENDDO
379
380C-----------------final loop--------------
381        Q_ext=0.0
382        Q_sca=0.0
383        g=0.0
384        DO n=Nmax-1,1,-1
385          nnn=FLOAT(n)
386          Q_ext=Q_ext+ (2.*nnn+1.) * REAL( a(n)+b(n) ) 
387          Q_sca=Q_sca+ (2.*nnn+1.) * 
388     .               REAL( a(n)*CONJG(a(n)) + b(n)*CONJG(b(n)) )
389          g=g + nnn*(nnn+2.)/(nnn+1.) * 
390     .       REAL( a(n)*CONJG(a(n+1))+b(n)*CONJG(b(n+1)) )  + 
391     .          (2.*nnn+1.)/nnn/(nnn+1.) * REAL(a(n)*CONJG(b(n)))
392        ENDDO
393        Q_ext=2./x**2 * Q_ext
394        Q_sca=2./x**2 * Q_sca
395        Q_abs=Q_ext-Q_sca
396        IF (AIMAG(m).EQ.0.0) Q_abs=0.0
397        omega=Q_sca/Q_ext
398        g=g*4./x**2/Q_sca
399c
400        sigma_sca=sigma_sca+r**2*Q_sca*number*deltar
401        sigma_abs=sigma_abs+r**2*Q_abs*number*deltar
402        sigma_ext=sigma_ext+r**2*Q_ext*number*deltar
403        omegatot=omegatot+r**2*Q_ext*omega*number*deltar
404        gtot    =gtot+r**2*Q_sca*g*number*deltar 
405c
406        ENDDO   !---bin
407C------------------------------------------------------------------
408
409        sigma_sca=pi*sigma_sca
410        sigma_abs=pi*sigma_abs
411        sigma_ext=pi*sigma_ext
412        gtot=pi*gtot/sigma_sca
413        omegatot=pi*omegatot/sigma_ext
414c
415        final_g(Nwv)=gtot
416        final_w(Nwv)=omegatot
417        final_a(Nwv)=sigma_ext/masse
418c
419c--averaged Qext for Yves
420        final_qext(Nwv)=sigma_ext/(surface/4.)
421c
422        ENDDO  !--loop on wavelength
423c
424c---averaging over LMDZ wavebands
425c
426        DO band=1, NbandSW
427          gcm_a(band)=0.0
428          gcm_g(band)=0.0
429          gcm_w(band)=0.0
430          gcm_qext(band)=0.0
431          gcm_weight_a(band)=0.0
432          gcm_weight_g(band)=0.0
433          gcm_weight_w(band)=0.0
434          gcm_weight_qext(band)=0.0
435        ENDDO
436c
437c---band 1 is now in the UV, so we take the first wavelength as being representative 
438c---it doesn't matter anyway because all radiation is absorbed in the stratosphere
439       DO Nwv=1,1
440          band=1
441          gcm_a(band)=gcm_a(band)+final_a(Nwv)*weight(Nwv)
442          gcm_weight_a(band)=gcm_weight_a(band)+weight(Nwv)
443c
444          gcm_w(band)=gcm_w(band)+
445     .                final_w(Nwv)*final_a(Nwv)*weight(Nwv)
446          gcm_weight_w(band)=gcm_weight_w(band)+
447     .                       final_a(Nwv)*weight(Nwv)
448c
449          gcm_g(band)=gcm_g(band)+
450     .                final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
451          gcm_weight_g(band)=gcm_weight_g(band)+
452     .                       final_a(Nwv)*final_w(Nwv)*weight(Nwv)
453c
454          gcm_qext(band)=gcm_qext(band)+final_qext(Nwv)*weight(Nwv)
455          gcm_weight_qext(band)=gcm_weight_qext(band)+weight(Nwv)
456        ENDDO
457c
458       DO Nwv=1,NwvmaxSW-1
459c
460        IF (Nwv.LE.5) THEN          !--RRTM spectral interval 2
461          band=2
462        ELSEIF (Nwv.LE.10) THEN     !--RRTM spectral interval 3
463          band=3
464        ELSEIF (Nwv.LE.16) THEN     !--RRTM spectral interval 4
465          band=4
466        ELSEIF (Nwv.LE.21) THEN     !--RRTM spectral interval 5
467          band=5
468        ELSE                        !--RRTM spectral interval 6
469          band=6
470        ENDIF
471c
472        gcm_a(band)=gcm_a(band)+final_a(Nwv)*weight(Nwv)
473        gcm_weight_a(band)=gcm_weight_a(band)+weight(Nwv)
474c
475        gcm_w(band)=gcm_w(band)+
476     .              final_w(Nwv)*final_a(Nwv)*weight(Nwv)
477        gcm_weight_w(band)=gcm_weight_w(band)+
478     .                     final_a(Nwv)*weight(Nwv)
479c
480        gcm_g(band)=gcm_g(band)+
481     .              final_g(Nwv)*final_a(Nwv)*final_w(Nwv)*weight(Nwv)
482        gcm_weight_g(band)=gcm_weight_g(band)+
483     .                     final_a(Nwv)*final_w(Nwv)*weight(Nwv)
484c
485        gcm_qext(band)=gcm_qext(band)+final_qext(Nwv)*weight(Nwv)
486        gcm_weight_qext(band)=gcm_weight_qext(band)+weight(Nwv)
487
488        ENDDO
489c
490        DO band=1, NbandSW
491          gcm_a(band)=gcm_a(band)/gcm_weight_a(band)
492          gcm_w(band)=gcm_w(band)/gcm_weight_w(band)
493          gcm_g(band)=gcm_g(band)/gcm_weight_g(band)
494          gcm_qext(band)=gcm_qext(band)/gcm_weight_qext(band)
495          ss_a(band,irh)=gcm_a(band)
496          ss_w(band,irh)=gcm_w(band)
497          ss_g(band,irh)=gcm_g(band)
498          ss_qext(band,irh)=gcm_qext(band)
499        ENDDO
500c
501        DO nb=NbandSW+1, NbandSW+nb_lambda
502          ss_a(nb,irh)=final_a(NwvmaxSW-1+nb-NbandSW)
503          ss_w(nb,irh)=final_w(NwvmaxSW-1+nb-NbandSW)
504          ss_g(nb,irh)=final_g(NwvmaxSW-1+nb-NbandSW)
505          ss_qext(nb,irh)=final_qext(NwvmaxSW-1+nb-NbandSW)
506        ENDDO
507c
508        ENDDO   !--fin boucle sur RH
509c
510c--Outputs
511C
512        IF (mode.EQ.1) THEN !--only CS because AS treated by internal mixing routine
513
514        WRITE(14,*) ' ! '//chmode(mode)
515        DO k=1, NbandSW
516        WRITE(14,951) (ss_a(k,irh),irh=1,Nrh)
517        ENDDO
518        WRITE(15,*) ' ! '//chmode(mode)
519        DO k=1, NbandSW
520        WRITE(15,951) (ss_g(k,irh),irh=1,Nrh)
521        ENDDO
522        WRITE(16,*) ' ! '//chmode(mode)
523        DO k=1, NbandSW
524        WRITE(16,951) (ss_w(k,irh),irh=1,Nrh)
525        ENDDO
526        DO k=1, NbandSW
527        WRITE(17,951) (ss_qext(k,irh),irh=1,Nrh)
528        ENDDO
529c
530        WRITE(24,*) ' ! extinction '//chmode(mode)
531        DO k=NbandSW+1,NbandSW+nb_lambda
532        WRITE(24,951) (ss_a(k,irh),irh=1,Nrh) 
533        ENDDO
534        WRITE(26,*) ' ! absorption '//chmode(mode)
535        DO k=NbandSW+1,NbandSW+nb_lambda
536        WRITE(26,951) ((1.0-ss_w(k,irh))*ss_a(k,irh),irh=1,Nrh) 
537        ENDDO
538c
539        WRITE(25,*) ' ! '//chmode(mode)
540        DO k=NbandSW+1,NbandSW+nb_lambda
541        WRITE(25,951) (ss_qext(k,irh),irh=1,Nrh) 
542        ENDDO
543c
544        ENDIF
545c
546        ENDDO !--boucle sur les modes
547
548951     FORMAT(1X,12(F6.3,','),' &')
549c
550        CLOSE(14)
551        CLOSE(15)
552        CLOSE(16)
553        CLOSE(17)
554c
555        CLOSE(24)
556        CLOSE(25)
557        CLOSE(26)
558c
559        END
Note: See TracBrowser for help on using the repository browser.