source: trunk/LMDZ.MARS/libf/phymars/aeroptproperties.F @ 259

Last change on this file since 259 was 89, checked in by aslmd, 14 years ago

mars: nouvelle version amelioree de aeroptproperties.F dans phymars. permet de limiter les incertitudes et laisse la possibilite d'augmenter ngau. la derniere version de cette routine sur les sources reference de JB [dans /u/jbmlmd/LMDZ.VERSO/11W02_SP40_LMDZ.BETA.RECHERCHE] n'avait pas ete incluse dans la version de reference SVN [et LMDZ.MARS.BETA].
File size: 52.9 KB
Line 
1      SUBROUTINE aeroptproperties(ngrid,nlayer,reffrad,nueffrad,
2     &                            QVISsQREF3d,omegaVIS3d,gVIS3d,
3     &                            QIRsQREF3d,omegaIR3d,gIR3d,
4     &                            QREFvis3d,QREFir3d,
5     &                            omegaREFvis3d,omegaREFir3d)
6      IMPLICIT NONE
7c     =============================================================
8c     Aerosol Optical Properties
9c
10c     Description:
11c       Compute the scattering parameters in each grid
12c       box, depending on aerosol grain sizes. Log-normal size
13c       distribution and Gauss-Legendre integration are used.
14
15c     Parameters:
16c       Don't forget to set the value of varyingnueff below; If
17c       the effective variance of the distribution for the given
18c       aerosol is considered homogeneous in the atmosphere, please
19c       set varyingnueff(iaer) to .false. Resulting computational
20c       time will be much better.
21
22c     Authors: J.-B. Madeleine, F. Forget, F. Montmessin
23c     =============================================================
24
25#include "dimensions.h"
26#include "dimphys.h"
27#include "callkeys.h"
28#include "dimradmars.h"
29#include "yomaer.h"
30
31c     Local variables
32c     ---------------
33
34c     =============================================================
35      LOGICAL, PARAMETER :: varyingnueff(naerkind) = .false.
36c     =============================================================
37
38c     Min. and max radius of the interpolation grid (in METERS)
39      REAL, PARAMETER :: refftabmin = 1e-8 !5e-8
40      REAL, PARAMETER :: refftabmax = 35e-6
41c     Log of the min and max variance of the interpolation grid
42      REAL, PARAMETER :: nuefftabmin = -4.6
43      REAL, PARAMETER :: nuefftabmax = 0.
44c     Number of effective radius of the interpolation grid
45      INTEGER, PARAMETER :: refftabsize = 100
46c     Number of effective variances of the interpolation grid
47      INTEGER, PARAMETER :: nuefftabsize = 100
48c     Interpolation grid indices (reff,nueff)
49      INTEGER :: grid_i,grid_j
50c     Intermediate variable
51      REAL :: var_tmp,var3d_tmp(ngridmx,nlayermx)
52c     Bilinear interpolation factors
53      REAL :: kx,ky,k1,k2,k3,k4
54c     Size distribution parameters
55      REAL :: sizedistk1,sizedistk2
56c     Pi!
57      REAL,SAVE :: pi
58c     Variables used by the Gauss-Legendre integration:
59      INTEGER radius_id,gausind
60      REAL kint
61      REAL drad
62      INTEGER, PARAMETER :: ngau = 10
63      REAL weightgaus(ngau),radgaus(ngau)
64      SAVE weightgaus,radgaus
65c     DATA weightgaus/.2955242247,.2692667193,.2190863625,
66c    &                .1494513491,.0666713443/
67c     DATA radgaus/.1488743389,.4333953941,.6794095682,
68c    &             .8650633666,.9739065285/
69      DATA    radgaus/0.07652652113350,0.22778585114165,
70     &                0.37370608871528,0.51086700195146,
71     &                0.63605368072468,0.74633190646476,
72     &                0.83911697181213,0.91223442826796,
73     &                0.96397192726078,0.99312859919241/
74
75      DATA weightgaus/0.15275338723120,0.14917298659407,
76     &                0.14209610937519,0.13168863843930,
77     &                0.11819453196154,0.10193011980823,
78     &                0.08327674160932,0.06267204829828,
79     &                0.04060142982019,0.01761400714091/
80
81c     Indices
82      INTEGER :: i,j,k,l,m,iaer,idomain
83      INTEGER :: ig,lg,chg
84
85c     Local saved variables
86c     ---------------------
87
88c     Radius axis of the interpolation grid
89      REAL,SAVE :: refftab(refftabsize)
90c     Variance axis of the interpolation grid
91      REAL,SAVE :: nuefftab(nuefftabsize)
92c     Volume ratio of the grid
93      REAL,SAVE :: logvratgrid,vratgrid
94c     Grid used to remember which calculation is done
95      LOGICAL,SAVE :: checkgrid(refftabsize,nuefftabsize,naerkind,2)
96     &                = .false.
97c     Optical properties of the grid (VISIBLE)
98      REAL,SAVE :: qsqrefVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
99      REAL,SAVE :: qextVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
100      REAL,SAVE :: qscatVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
101      REAL,SAVE :: omegVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
102      REAL,SAVE :: gVISgrid(refftabsize,nuefftabsize,nsun,naerkind)
103c     Optical properties of the grid (INFRARED)
104      REAL,SAVE :: qsqrefIRgrid(refftabsize,nuefftabsize,nir,naerkind)
105      REAL,SAVE :: qextIRgrid(refftabsize,nuefftabsize,nir,naerkind)
106      REAL,SAVE :: qscatIRgrid(refftabsize,nuefftabsize,nir,naerkind)
107      REAL,SAVE :: omegIRgrid(refftabsize,nuefftabsize,nir,naerkind)
108      REAL,SAVE :: gIRgrid(refftabsize,nuefftabsize,nir,naerkind)
109c     Optical properties of the grid (REFERENCE WAVELENGTHS)
110      REAL,SAVE :: qrefVISgrid(refftabsize,nuefftabsize,naerkind)
111      REAL,SAVE :: qscatrefVISgrid(refftabsize,nuefftabsize,naerkind)
112      REAL,SAVE :: qrefIRgrid(refftabsize,nuefftabsize,naerkind)
113      REAL,SAVE :: qscatrefIRgrid(refftabsize,nuefftabsize,naerkind)
114      REAL,SAVE :: omegrefVISgrid(refftabsize,nuefftabsize,naerkind)
115      REAL,SAVE :: omegrefIRgrid(refftabsize,nuefftabsize,naerkind)
116c     Firstcall
117      LOGICAL,SAVE :: firstcall = .true.
118c     Variables used by the Gauss-Legendre integration:
119      REAL,SAVE :: normd(refftabsize,nuefftabsize,naerkind,2)
120      REAL,SAVE :: dista(refftabsize,nuefftabsize,naerkind,2,ngau)
121      REAL,SAVE :: distb(refftabsize,nuefftabsize,naerkind,2,ngau)
122
123      REAL,SAVE :: radGAUSa(ngau,naerkind,2)
124      REAL,SAVE :: radGAUSb(ngau,naerkind,2)
125
126      REAL,SAVE :: qsqrefVISa(nsun,ngau,naerkind)
127      REAL,SAVE :: qrefVISa(ngau,naerkind)
128      REAL,SAVE :: qsqrefVISb(nsun,ngau,naerkind)
129      REAL,SAVE :: qrefVISb(ngau,naerkind)
130      REAL,SAVE :: omegVISa(nsun,ngau,naerkind)
131      REAL,SAVE :: omegrefVISa(ngau,naerkind)
132      REAL,SAVE :: omegVISb(nsun,ngau,naerkind)
133      REAL,SAVE :: omegrefVISb(ngau,naerkind)
134      REAL,SAVE :: gVISa(nsun,ngau,naerkind)
135      REAL,SAVE :: gVISb(nsun,ngau,naerkind)
136
137      REAL,SAVE :: qsqrefIRa(nir,ngau,naerkind)
138      REAL,SAVE :: qrefIRa(ngau,naerkind)
139      REAL,SAVE :: qsqrefIRb(nir,ngau,naerkind)
140      REAL,SAVE :: qrefIRb(ngau,naerkind)
141      REAL,SAVE :: omegIRa(nir,ngau,naerkind)
142      REAL,SAVE :: omegrefIRa(ngau,naerkind)
143      REAL,SAVE :: omegIRb(nir,ngau,naerkind)
144      REAL,SAVE :: omegrefIRb(ngau,naerkind)
145      REAL,SAVE :: gIRa(nir,ngau,naerkind)
146      REAL,SAVE :: gIRb(nir,ngau,naerkind)
147
148      REAL :: radiusm
149      REAL :: radiusr
150
151c     Inputs
152c     ------
153
154      INTEGER :: ngrid,nlayer
155c     Aerosol effective radius used for radiative transfer (meter)
156      REAL :: reffrad(ngridmx,nlayermx,naerkind)
157c     Aerosol effective variance used for radiative transfer (n.u.)
158      REAL :: nueffrad(ngridmx,nlayermx,naerkind)
159
160c     Outputs
161c     -------
162
163      REAL :: QVISsQREF3d(ngridmx,nlayermx,nsun,naerkind)
164      REAL :: omegaVIS3d(ngridmx,nlayermx,nsun,naerkind)
165      REAL :: gVIS3d(ngridmx,nlayermx,nsun,naerkind)
166
167      REAL :: QIRsQREF3d(ngridmx,nlayermx,nir,naerkind)
168      REAL :: omegaIR3d(ngridmx,nlayermx,nir,naerkind)
169      REAL :: gIR3d(ngridmx,nlayermx,nir,naerkind)
170
171      REAL :: QREFvis3d(ngridmx,nlayermx,naerkind)
172      REAL :: QREFir3d(ngridmx,nlayermx,naerkind)
173
174      REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind)
175      REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind)
176
177c     Tests
178c     -----
179
180      LOGICAL,SAVE :: out_qwg = .false.
181      INTEGER, PARAMETER :: out_iaer = 1
182      REAL :: out_qext(ngridmx,nlayermx)
183      REAL :: out_omeg(ngridmx,nlayermx)
184      REAL :: out_g(ngridmx,nlayermx)
185      INTEGER :: out_nchannel
186      CHARACTER*1 :: out_str
187
188      DO iaer = 1, naerkind ! Loop on aerosol kind
189        IF ( (nsize(iaer,1).EQ.1).AND.(nsize(iaer,2).EQ.1) ) THEN
190c==================================================================
191c       If there is one single particle size, optical
192c         properties of the considered aerosol are homogeneous
193          DO lg = 1, nlayer
194            DO ig = 1, ngrid
195              DO chg = 1, nsun
196                QVISsQREF3d(ig,lg,chg,iaer)=QVISsQREF(chg,iaer,1)
197                omegaVIS3d(ig,lg,chg,iaer)=omegaVIS(chg,iaer,1)
198                gVIS3d(ig,lg,chg,iaer)=gVIS(chg,iaer,1)
199              ENDDO
200              DO chg = 1, nir
201                QIRsQREF3d(ig,lg,chg,iaer)=QIRsQREF(chg,iaer,1)
202                omegaIR3d(ig,lg,chg,iaer)=omegaIR(chg,iaer,1)
203                gIR3d(ig,lg,chg,iaer)=gIR(chg,iaer,1)
204              ENDDO
205              QREFvis3d(ig,lg,iaer)=QREFvis(iaer,1)
206              QREFir3d(ig,lg,iaer)=QREFir(iaer,1)
207              omegaREFvis3d(ig,lg,iaer)=omegaREFvis(iaer,1)
208              omegaREFir3d(ig,lg,iaer)=omegaREFir(iaer,1)
209            ENDDO
210          ENDDO
211        ELSE ! Varying effective radius and variance
212      DO idomain = 1, 2 ! Loop on visible or infrared channel
213c==================================================================
214c     1. Creating the effective radius and variance grid
215c     --------------------------------------------------
216      IF (firstcall) THEN
217
218c       1.1 Pi!
219        pi = 2. * asin(1.e0)
220
221c       1.2 Effective radius
222        refftab(1)    = refftabmin
223        refftab(refftabsize) = refftabmax
224
225        logvratgrid = log(refftabmax/refftabmin) /
226     &                float(refftabsize-1)*3.
227        vratgrid = exp(logvratgrid)
228
229        do i = 2, refftabsize-1
230          refftab(i) = refftab(i-1)*vratgrid**(1./3.)
231        enddo
232
233c       1.3 Effective variance
234        do i = 0, nuefftabsize-1
235          nuefftab(i+1) = exp( nuefftabmin +
236     &             i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) )
237        enddo
238        firstcall = .false.
239      ENDIF
240
241c       1.4 Radius middle point and range for Gauss integration
242        radiusm=
243     &         0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) +
244     &              radiustab(iaer,idomain,1))
245        radiusr=
246     &         0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) -
247     &              radiustab(iaer,idomain,1))
248
249c       1.5 Interpolating data at the Gauss quadrature points:
250        DO gausind=1,ngau
251          drad=radiusr*radgaus(gausind)
252          radGAUSa(gausind,iaer,idomain)=radiusm-drad
253
254          radius_id=minloc(abs(radiustab(iaer,idomain,:) -
255     &                         (radiusm-drad)),DIM=1)
256          IF ((radiustab(iaer,idomain,radius_id) -
257     &         (radiusm-drad)).GT.0) THEN
258            radius_id=radius_id-1
259          ENDIF
260          IF (radius_id.GE.nsize(iaer,idomain)) THEN
261            radius_id=nsize(iaer,idomain)-1
262            kint = 1.
263          ELSEIF (radius_id.LT.1) THEN
264            radius_id=1
265            kint = 0.
266          ELSE
267          kint = ( (radiusm-drad) -
268     &             radiustab(iaer,idomain,radius_id) ) /
269     &           ( radiustab(iaer,idomain,radius_id+1) -
270     &             radiustab(iaer,idomain,radius_id) )
271          ENDIF
272          IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
273            DO m=1,nsun
274            qsqrefVISa(m,gausind,iaer)=
275     &              (1-kint)*QVISsQREF(m,iaer,radius_id) +
276     &              kint*QVISsQREF(m,iaer,radius_id+1)
277            omegVISa(m,gausind,iaer)=
278     &              (1-kint)*omegaVIS(m,iaer,radius_id) +
279     &              kint*omegaVIS(m,iaer,radius_id+1)
280            gVISa(m,gausind,iaer)=
281     &              (1-kint)*gVIS(m,iaer,radius_id) +
282     &              kint*gVIS(m,iaer,radius_id+1)
283            ENDDO
284            qrefVISa(gausind,iaer)=
285     &              (1-kint)*QREFvis(iaer,radius_id) +
286     &              kint*QREFvis(iaer,radius_id+1)
287            omegrefVISa(gausind,iaer)=
288     &              (1-kint)*omegaREFvis(iaer,radius_id) +
289     &              kint*omegaREFvis(iaer,radius_id+1)
290          ELSE ! INFRARED DOMAIN ----------------------------------
291            DO m=1,nir
292            qsqrefIRa(m,gausind,iaer)=
293     &              (1-kint)*QIRsQREF(m,iaer,radius_id) +
294     &              kint*QIRsQREF(m,iaer,radius_id+1)
295            omegIRa(m,gausind,iaer)=
296     &              (1-kint)*omegaIR(m,iaer,radius_id) +
297     &              kint*omegaIR(m,iaer,radius_id+1)
298            gIRa(m,gausind,iaer)=
299     &              (1-kint)*gIR(m,iaer,radius_id) +
300     &              kint*gIR(m,iaer,radius_id+1)
301            ENDDO
302            qrefIRa(gausind,iaer)=
303     &              (1-kint)*QREFir(iaer,radius_id) +
304     &              kint*QREFir(iaer,radius_id+1)
305            omegrefIRa(gausind,iaer)=
306     &              (1-kint)*omegaREFir(iaer,radius_id) +
307     &              kint*omegaREFir(iaer,radius_id+1)
308          ENDIF
309        ENDDO
310
311        DO gausind=1,ngau
312          drad=radiusr*radgaus(gausind)
313          radGAUSb(gausind,iaer,idomain)=radiusm+drad
314
315          radius_id=minloc(abs(radiustab(iaer,idomain,:) -
316     &                         (radiusm+drad)),DIM=1)
317          IF ((radiustab(iaer,idomain,radius_id) -
318     &         (radiusm+drad)).GT.0) THEN
319            radius_id=radius_id-1
320          ENDIF
321          IF (radius_id.GE.nsize(iaer,idomain)) THEN
322            radius_id=nsize(iaer,idomain)-1
323            kint = 1.
324          ELSEIF (radius_id.LT.1) THEN
325            radius_id=1
326            kint = 0.
327          ELSE
328            kint = ( (radiusm+drad) -
329     &               radiustab(iaer,idomain,radius_id) ) /
330     &             ( radiustab(iaer,idomain,radius_id+1) -
331     &               radiustab(iaer,idomain,radius_id) )
332          ENDIF
333          IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
334            DO m=1,nsun
335            qsqrefVISb(m,gausind,iaer)=
336     &              (1-kint)*QVISsQREF(m,iaer,radius_id) +
337     &              kint*QVISsQREF(m,iaer,radius_id+1)
338            omegVISb(m,gausind,iaer)=
339     &              (1-kint)*omegaVIS(m,iaer,radius_id) +
340     &              kint*omegaVIS(m,iaer,radius_id+1)
341            gVISb(m,gausind,iaer)=
342     &              (1-kint)*gVIS(m,iaer,radius_id) +
343     &              kint*gVIS(m,iaer,radius_id+1)
344            ENDDO
345            qrefVISb(gausind,iaer)=
346     &              (1-kint)*QREFvis(iaer,radius_id) +
347     &              kint*QREFvis(iaer,radius_id+1)
348            omegrefVISb(gausind,iaer)=
349     &              (1-kint)*omegaREFvis(iaer,radius_id) +
350     &              kint*omegaREFvis(iaer,radius_id+1)
351          ELSE ! INFRARED DOMAIN ----------------------------------
352            DO m=1,nir
353            qsqrefIRb(m,gausind,iaer)=
354     &              (1-kint)*QIRsQREF(m,iaer,radius_id) +
355     &              kint*QIRsQREF(m,iaer,radius_id+1)
356            omegIRb(m,gausind,iaer)=
357     &              (1-kint)*omegaIR(m,iaer,radius_id) +
358     &              kint*omegaIR(m,iaer,radius_id+1)
359            gIRb(m,gausind,iaer)=
360     &              (1-kint)*gIR(m,iaer,radius_id) +
361     &              kint*gIR(m,iaer,radius_id+1)
362            ENDDO
363            qrefIRb(gausind,iaer)=
364     &              (1-kint)*QREFir(iaer,radius_id) +
365     &              kint*QREFir(iaer,radius_id+1)
366            omegrefIRb(gausind,iaer)=
367     &              (1-kint)*omegaREFir(iaer,radius_id) +
368     &              kint*omegaREFir(iaer,radius_id+1)
369          ENDIF
370        ENDDO
371c==================================================================
372      IF ( .NOT.varyingnueff(iaer) ) THEN          ! CONSTANT NUEFF
373c==================================================================
374c     2. Compute the scattering parameters using linear
375c       interpolation over grain sizes and constant nueff
376c     ---------------------------------------------------
377
378      DO lg = 1,nlayer
379        DO ig = 1,ngrid
380c         2.1 Effective radius index and kx calculation
381          var_tmp=reffrad(ig,lg,iaer)/refftabmin
382          var_tmp=log(var_tmp)*3.
383          var_tmp=var_tmp/logvratgrid+1.
384          grid_i=floor(var_tmp)
385          IF (grid_i.GE.refftabsize) THEN
386c           WRITE(*,*) 'Warning: particle size in grid box #'
387c           WRITE(*,*) ig,' is too large to be used by the '
388c           WRITE(*,*) 'radiative transfer; please extend the '
389c           WRITE(*,*) 'interpolation grid to larger grain sizes.'
390            grid_i=refftabsize-1
391            kx = 1.
392          ELSEIF (grid_i.LT.1) THEN
393c           WRITE(*,*) 'Warning: particle size in grid box #'
394c           WRITE(*,*) ig,' is too small to be used by the '
395c           WRITE(*,*) 'radiative transfer; please extend the '
396c           WRITE(*,*) 'interpolation grid to smaller grain sizes.'
397            grid_i=1
398            kx = 0.
399          ELSE
400            kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) /
401     &           ( refftab(grid_i+1)-refftab(grid_i) )
402          ENDIF
403c         2.3 Integration
404          DO j=grid_i,grid_i+1
405c             2.3.1 Check if the calculation has been done
406              IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN
407c               2.3.2 Log-normal dist., r_g and sigma_g are defined
408c                 in [hansen_1974], "Light scattering in planetary
409c                 atmospheres", Space Science Reviews 16 527-610.
410c                 Here, sizedistk1=r_g and sizedistk2=sigma_g^2
411                sizedistk2 = log(1.+nueffrad(1,1,iaer))
412                sizedistk1 = exp(2.5*sizedistk2)
413                sizedistk1 = refftab(j) / sizedistk1
414
415                normd(j,1,iaer,idomain) = 1e-30
416                DO gausind=1,ngau
417                  drad=radiusr*radgaus(gausind)
418                  dista(j,1,iaer,idomain,gausind) =
419     &              LOG((radiusm-drad)/sizedistk1)
420                  dista(j,1,iaer,idomain,gausind) =
421     &              EXP(-dista(j,1,iaer,idomain,gausind) *
422     &              dista(j,1,iaer,idomain,gausind) *
423     &              0.5e0/sizedistk2)/(radiusm-drad)
424                  dista(j,1,iaer,idomain,gausind) =
425     &              dista(j,1,iaer,idomain,gausind) /
426     &              (sqrt(2e0*pi*sizedistk2))
427
428                  distb(j,1,iaer,idomain,gausind) =
429     &              LOG((radiusm+drad)/sizedistk1)
430                  distb(j,1,iaer,idomain,gausind) =
431     &              EXP(-distb(j,1,iaer,idomain,gausind) *
432     &              distb(j,1,iaer,idomain,gausind) *
433     &              0.5e0/sizedistk2)/(radiusm+drad)
434                  distb(j,1,iaer,idomain,gausind) =
435     &              distb(j,1,iaer,idomain,gausind) /
436     &              (sqrt(2e0*pi*sizedistk2))
437
438                  normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) +
439     &              weightgaus(gausind) *
440     &              (
441     &              distb(j,1,iaer,idomain,gausind) * pi *
442     &              radGAUSb(gausind,iaer,idomain) *
443     &              radGAUSb(gausind,iaer,idomain) +
444     &              dista(j,1,iaer,idomain,gausind) * pi *
445     &              radGAUSa(gausind,iaer,idomain) *
446     &              radGAUSa(gausind,iaer,idomain)
447     &              )
448                ENDDO
449                IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
450c                 2.3.3.vis Initialization
451                  qsqrefVISgrid(j,1,:,iaer)=0.
452                  qextVISgrid(j,1,:,iaer)=0.
453                  qscatVISgrid(j,1,:,iaer)=0.
454                  omegVISgrid(j,1,:,iaer)=0.
455                  gVISgrid(j,1,:,iaer)=0.
456                  qrefVISgrid(j,1,iaer)=0.
457                  qscatrefVISgrid(j,1,iaer)=0.
458                  omegrefVISgrid(j,1,iaer)=0.
459
460                  DO gausind=1,ngau
461                    DO m=1,nsun
462c                     Convolution:
463                      qextVISgrid(j,1,m,iaer) =
464     &                  qextVISgrid(j,1,m,iaer) +
465     &                  weightgaus(gausind) *
466     &                  (
467     &                  qsqrefVISb(m,gausind,iaer) *
468     &                  qrefVISb(gausind,iaer) *
469     &                  pi*radGAUSb(gausind,iaer,idomain) *
470     &                  radGAUSb(gausind,iaer,idomain) *
471     &                  distb(j,1,iaer,idomain,gausind) +
472     &                  qsqrefVISa(m,gausind,iaer) *
473     &                  qrefVISa(gausind,iaer) *
474     &                  pi*radGAUSa(gausind,iaer,idomain) *
475     &                  radGAUSa(gausind,iaer,idomain) *
476     &                  dista(j,1,iaer,idomain,gausind)
477     &                  )
478                      qscatVISgrid(j,1,m,iaer) =
479     &                  qscatVISgrid(j,1,m,iaer) +
480     &                  weightgaus(gausind) *
481     &                  (
482     &                  omegVISb(m,gausind,iaer) *
483     &                  qsqrefVISb(m,gausind,iaer) *
484     &                  qrefVISb(gausind,iaer) *
485     &                  pi*radGAUSb(gausind,iaer,idomain) *
486     &                  radGAUSb(gausind,iaer,idomain) *
487     &                  distb(j,1,iaer,idomain,gausind) +
488     &                  omegVISa(m,gausind,iaer) *
489     &                  qsqrefVISa(m,gausind,iaer) *
490     &                  qrefVISa(gausind,iaer) *
491     &                  pi*radGAUSa(gausind,iaer,idomain) *
492     &                  radGAUSa(gausind,iaer,idomain) *
493     &                  dista(j,1,iaer,idomain,gausind)
494     &                  )
495                      gVISgrid(j,1,m,iaer) =
496     &                  gVISgrid(j,1,m,iaer) +
497     &                  weightgaus(gausind) *
498     &                  (
499     &                  omegVISb(m,gausind,iaer) *
500     &                  qsqrefVISb(m,gausind,iaer) *
501     &                  qrefVISb(gausind,iaer) *
502     &                  gVISb(m,gausind,iaer) *
503     &                  pi*radGAUSb(gausind,iaer,idomain) *
504     &                  radGAUSb(gausind,iaer,idomain) *
505     &                  distb(j,1,iaer,idomain,gausind) +
506     &                  omegVISa(m,gausind,iaer) *
507     &                  qsqrefVISa(m,gausind,iaer) *
508     &                  qrefVISa(gausind,iaer) *
509     &                  gVISa(m,gausind,iaer) *
510     &                  pi*radGAUSa(gausind,iaer,idomain) *
511     &                  radGAUSa(gausind,iaer,idomain) *
512     &                  dista(j,1,iaer,idomain,gausind)
513     &                  )
514                    ENDDO
515                    qrefVISgrid(j,1,iaer) =
516     &                qrefVISgrid(j,1,iaer) +
517     &                weightgaus(gausind) *
518     &                (
519     &                qrefVISb(gausind,iaer) *
520     &                pi*radGAUSb(gausind,iaer,idomain) *
521     &                radGAUSb(gausind,iaer,idomain) *
522     &                distb(j,1,iaer,idomain,gausind) +
523     &                qrefVISa(gausind,iaer) *
524     &                pi*radGAUSa(gausind,iaer,idomain) *
525     &                radGAUSa(gausind,iaer,idomain) *
526     &                dista(j,1,iaer,idomain,gausind)
527     &                )
528                    qscatrefVISgrid(j,1,iaer) =
529     &                qscatrefVISgrid(j,1,iaer) +
530     &                weightgaus(gausind) *
531     &                (
532     &                omegrefVISb(gausind,iaer) *
533     &                qrefVISb(gausind,iaer) *
534     &                pi*radGAUSb(gausind,iaer,idomain) *
535     &                radGAUSb(gausind,iaer,idomain) *
536     &                distb(j,1,iaer,idomain,gausind) +
537     &                omegrefVISa(gausind,iaer) *
538     &                qrefVISa(gausind,iaer) *
539     &                pi*radGAUSa(gausind,iaer,idomain) *
540     &                radGAUSa(gausind,iaer,idomain) *
541     &                dista(j,1,iaer,idomain,gausind)
542     &                )
543                  ENDDO
544
545                  qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) /
546     &                          normd(j,1,iaer,idomain)
547                  qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /
548     &                          normd(j,1,iaer,idomain)
549                  omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /
550     &                         qrefVISgrid(j,1,iaer)
551                  DO m=1,nsun
552                    qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /
553     &                          normd(j,1,iaer,idomain)
554                    qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /
555     &                          normd(j,1,iaer,idomain)
556                    gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) /
557     &                          qscatVISgrid(j,1,m,iaer) /
558     &                          normd(j,1,iaer,idomain)
559
560                    qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /
561     &                          qrefVISgrid(j,1,iaer)
562                    omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /
563     &                          qextVISgrid(j,1,m,iaer)
564                  ENDDO
565                ELSE                   ! INFRARED DOMAIN ----------
566c                 2.3.3.ir Initialization
567                  qsqrefIRgrid(j,1,:,iaer)=0.
568                  qextIRgrid(j,1,:,iaer)=0.
569                  qscatIRgrid(j,1,:,iaer)=0.
570                  omegIRgrid(j,1,:,iaer)=0.
571                  gIRgrid(j,1,:,iaer)=0.
572                  qrefIRgrid(j,1,iaer)=0.
573                  qscatrefIRgrid(j,1,iaer)=0.
574                  omegrefIRgrid(j,1,iaer)=0.
575
576                  DO gausind=1,ngau
577                    DO m=1,nir
578c                     Convolution:
579                      qextIRgrid(j,1,m,iaer) =
580     &                  qextIRgrid(j,1,m,iaer) +
581     &                  weightgaus(gausind) *
582     &                  (
583     &                  qsqrefIRb(m,gausind,iaer) *
584     &                  qrefVISb(gausind,iaer) *
585     &                  pi*radGAUSb(gausind,iaer,idomain) *
586     &                  radGAUSb(gausind,iaer,idomain) *
587     &                  distb(j,1,iaer,idomain,gausind) +
588     &                  qsqrefIRa(m,gausind,iaer) *
589     &                  qrefVISa(gausind,iaer) *
590     &                  pi*radGAUSa(gausind,iaer,idomain) *
591     &                  radGAUSa(gausind,iaer,idomain) *
592     &                  dista(j,1,iaer,idomain,gausind)
593     &                  )
594                      qscatIRgrid(j,1,m,iaer) =
595     &                  qscatIRgrid(j,1,m,iaer) +
596     &                  weightgaus(gausind) *
597     &                  (
598     &                  omegIRb(m,gausind,iaer) *
599     &                  qsqrefIRb(m,gausind,iaer) *
600     &                  qrefVISb(gausind,iaer) *
601     &                  pi*radGAUSb(gausind,iaer,idomain) *
602     &                  radGAUSb(gausind,iaer,idomain) *
603     &                  distb(j,1,iaer,idomain,gausind) +
604     &                  omegIRa(m,gausind,iaer) *
605     &                  qsqrefIRa(m,gausind,iaer) *
606     &                  qrefVISa(gausind,iaer) *
607     &                  pi*radGAUSa(gausind,iaer,idomain) *
608     &                  radGAUSa(gausind,iaer,idomain) *
609     &                  dista(j,1,iaer,idomain,gausind)
610     &                  )
611                      gIRgrid(j,1,m,iaer) =
612     &                  gIRgrid(j,1,m,iaer) +
613     &                  weightgaus(gausind) *
614     &                  (
615     &                  omegIRb(m,gausind,iaer) *
616     &                  qsqrefIRb(m,gausind,iaer) *
617     &                  qrefVISb(gausind,iaer) *
618     &                  gIRb(m,gausind,iaer) *
619     &                  pi*radGAUSb(gausind,iaer,idomain) *
620     &                  radGAUSb(gausind,iaer,idomain) *
621     &                  distb(j,1,iaer,idomain,gausind) +
622     &                  omegIRa(m,gausind,iaer) *
623     &                  qsqrefIRa(m,gausind,iaer) *
624     &                  qrefVISa(gausind,iaer) *
625     &                  gIRa(m,gausind,iaer) *
626     &                  pi*radGAUSa(gausind,iaer,idomain) *
627     &                  radGAUSa(gausind,iaer,idomain) *
628     &                  dista(j,1,iaer,idomain,gausind)
629     &                  )
630                    ENDDO
631                    qrefIRgrid(j,1,iaer) =
632     &                qrefIRgrid(j,1,iaer) +
633     &                weightgaus(gausind) *
634     &                (
635     &                qrefIRb(gausind,iaer) *
636     &                pi*radGAUSb(gausind,iaer,idomain) *
637     &                radGAUSb(gausind,iaer,idomain) *
638     &                distb(j,1,iaer,idomain,gausind) +
639     &                qrefIRa(gausind,iaer) *
640     &                pi*radGAUSa(gausind,iaer,idomain) *
641     &                radGAUSa(gausind,iaer,idomain) *
642     &                dista(j,1,iaer,idomain,gausind)
643     &                )
644                    qscatrefIRgrid(j,1,iaer) =
645     &                qscatrefIRgrid(j,1,iaer) +
646     &                weightgaus(gausind) *
647     &                (
648     &                omegrefIRb(gausind,iaer) *
649     &                qrefIRb(gausind,iaer) *
650     &                pi*radGAUSb(gausind,iaer,idomain) *
651     &                radGAUSb(gausind,iaer,idomain) *
652     &                distb(j,1,iaer,idomain,gausind) +
653     &                omegrefIRa(gausind,iaer) *
654     &                qrefIRa(gausind,iaer) *
655     &                pi*radGAUSa(gausind,iaer,idomain) *
656     &                radGAUSa(gausind,iaer,idomain) *
657     &                dista(j,1,iaer,idomain,gausind)
658     &                )
659                  ENDDO
660
661                  qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) /
662     &                          normd(j,1,iaer,idomain)
663                  qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /
664     &                          normd(j,1,iaer,idomain)
665                  omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /
666     &                         qrefIRgrid(j,1,iaer)
667                  DO m=1,nir
668                    qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /
669     &                          normd(j,1,iaer,idomain)
670                    qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /
671     &                          normd(j,1,iaer,idomain)
672                    gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) /
673     &                          qscatIRgrid(j,1,m,iaer) /
674     &                          normd(j,1,iaer,idomain)
675
676                    qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /
677     &                          qrefVISgrid(j,1,iaer)
678                    omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /
679     &                          qextIRgrid(j,1,m,iaer)
680                  ENDDO
681                ENDIF                  ! --------------------------
682                checkgrid(j,1,iaer,idomain) = .true.
683              ENDIF !checkgrid
684          ENDDO !grid_i
685c         2.4 Linear interpolation
686          k1 = (1-kx)
687          k2 = kx
688          IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
689          DO m=1,nsun
690            QVISsQREF3d(ig,lg,m,iaer) =
691     &                  k1*qsqrefVISgrid(grid_i,1,m,iaer) +
692     &                  k2*qsqrefVISgrid(grid_i+1,1,m,iaer)
693            omegaVIS3d(ig,lg,m,iaer) =
694     &                  k1*omegVISgrid(grid_i,1,m,iaer) +
695     &                  k2*omegVISgrid(grid_i+1,1,m,iaer)
696            gVIS3d(ig,lg,m,iaer) =
697     &                  k1*gVISgrid(grid_i,1,m,iaer) +
698     &                  k2*gVISgrid(grid_i+1,1,m,iaer)
699          ENDDO !nsun
700          QREFvis3d(ig,lg,iaer) =
701     &                  k1*qrefVISgrid(grid_i,1,iaer) +
702     &                  k2*qrefVISgrid(grid_i+1,1,iaer)
703          omegaREFvis3d(ig,lg,iaer) =
704     &                  k1*omegrefVISgrid(grid_i,1,iaer) +
705     &                  k2*omegrefVISgrid(grid_i+1,1,iaer)
706          ELSE                   ! INFRARED -----------------------
707          DO m=1,nir
708            QIRsQREF3d(ig,lg,m,iaer) =
709     &                  k1*qsqrefIRgrid(grid_i,1,m,iaer) +
710     &                  k2*qsqrefIRgrid(grid_i+1,1,m,iaer)
711            omegaIR3d(ig,lg,m,iaer) =
712     &                  k1*omegIRgrid(grid_i,1,m,iaer) +
713     &                  k2*omegIRgrid(grid_i+1,1,m,iaer)
714            gIR3d(ig,lg,m,iaer) =
715     &                  k1*gIRgrid(grid_i,1,m,iaer) +
716     &                  k2*gIRgrid(grid_i+1,1,m,iaer)
717          ENDDO !nir
718          QREFir3d(ig,lg,iaer) =
719     &                  k1*qrefIRgrid(grid_i,1,iaer) +
720     &                  k2*qrefIRgrid(grid_i+1,1,iaer)
721          omegaREFir3d(ig,lg,iaer) =
722     &                  k1*omegrefIRgrid(grid_i,1,iaer) +
723     &                  k2*omegrefIRgrid(grid_i+1,1,iaer)
724          ENDIF                  ! --------------------------------
725        ENDDO !nlayermx
726      ENDDO !ngridmx
727c==================================================================
728      ELSE                                          ! VARYING NUEFF
729c==================================================================
730c     3. Compute the scattering parameters in each grid box
731c       using bilinear interpolation over the grain sizes
732c       and the effective variances;
733c     -----------------------------------------------------
734
735      DO lg = 1,nlayer
736        DO ig = 1,ngrid
737c         3.1 Effective variance index and ky calculation
738          var_tmp=log(nueffrad(ig,lg,iaer))
739          grid_j=floor( (nuefftabsize-1)/(nuefftabmax-nuefftabmin)*
740     &       (var_tmp-nuefftabmin)+1. )
741          IF (grid_j.GE.nuefftabsize) THEN
742c           WRITE(*,*) 'Warning: effective variance '
743c           WRITE(*,*) 'is too large to be used by the '
744c           WRITE(*,*) 'radiative transfer; please extend the '
745c           WRITE(*,*) 'interpolation grid to larger values.'
746            grid_j=nuefftabsize-1
747            ky = 1.
748          ELSEIF (grid_j.LT.1) THEN
749c           WRITE(*,*) 'Warning: effective variance '
750c           WRITE(*,*) 'is too small to be used by the '
751c           WRITE(*,*) 'radiative transfer; please extend the '
752c           WRITE(*,*) 'interpolation grid to smaller values.'
753            grid_j=1
754            ky = 0.
755          ELSE
756            ky = ( nueffrad(ig,lg,iaer)-nuefftab(grid_j) ) /
757     &           ( nuefftab(grid_j+1)-nuefftab(grid_j) )
758          ENDIF
759c         3.2 Effective radius index and kx calculation
760          var_tmp=reffrad(ig,lg,iaer)/refftabmin
761          var_tmp=log(var_tmp)*3.
762          var_tmp=var_tmp/logvratgrid+1.
763          grid_i=floor(var_tmp)
764          IF (grid_i.GE.refftabsize) THEN
765c           WRITE(*,*) 'Warning: particle size in grid box #'
766c           WRITE(*,*) ig,' is too large to be used by the '
767c           WRITE(*,*) 'radiative transfer; please extend the '
768c           WRITE(*,*) 'interpolation grid to larger grain sizes.'
769            grid_i=refftabsize-1
770            kx = 1.
771          ELSEIF (grid_i.LT.1) THEN
772c           WRITE(*,*) 'Warning: particle size in grid box #'
773c           WRITE(*,*) ig,' is too small to be used by the '
774c           WRITE(*,*) 'radiative transfer; please extend the '
775c           WRITE(*,*) 'interpolation grid to smaller grain sizes.'
776            grid_i=1
777            kx = 0.
778          ELSE
779            kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) /
780     &           ( refftab(grid_i+1)-refftab(grid_i) )
781          ENDIF
782c         3.3 Integration
783          DO j=grid_i,grid_i+1
784            DO k=grid_j,grid_j+1
785c             3.3.1 Check if the calculation has been done
786              IF (.NOT.checkgrid(j,k,iaer,idomain)) THEN
787
788c               3.3.2 Log-normal dist., r_g and sigma_g are defined
789c                 in [hansen_1974], "Light scattering in planetary
790c                 atmospheres", Space Science Reviews 16 527-610.
791c                 Here, sizedistk1=r_g and sizedistk2=sigma_g^2
792                sizedistk2 = log(1.+nuefftab(k))
793                sizedistk1 = exp(2.5*sizedistk2)
794                sizedistk1 = refftab(j) / sizedistk1
795
796                normd(j,k,iaer,idomain) = 1e-30
797                DO gausind=1,ngau
798                  drad=radiusr*radgaus(gausind)
799
800                  dista(j,k,iaer,idomain,gausind) =
801     &              LOG((radiusm-drad)/sizedistk1)
802                  dista(j,k,iaer,idomain,gausind) =
803     &              EXP(-dista(j,k,iaer,idomain,gausind) *
804     &              dista(j,k,iaer,idomain,gausind) *
805     &              0.5e0/sizedistk2)/(radiusm-drad)
806                  dista(j,k,iaer,idomain,gausind) =
807     &              dista(j,k,iaer,idomain,gausind) /
808     &              (sqrt(2e0*pi*sizedistk2))
809
810                  distb(j,k,iaer,idomain,gausind) =
811     &              LOG((radiusm+drad)/sizedistk1)
812                  distb(j,k,iaer,idomain,gausind) =
813     &              EXP(-distb(j,k,iaer,idomain,gausind) *
814     &              distb(j,k,iaer,idomain,gausind) *
815     &              0.5e0/sizedistk2)/(radiusm+drad)
816                  distb(j,k,iaer,idomain,gausind) =
817     &              distb(j,k,iaer,idomain,gausind) /
818     &              (sqrt(2e0*pi*sizedistk2))
819
820                  normd(j,k,iaer,idomain)=normd(j,k,iaer,idomain) +
821     &              weightgaus(gausind) *
822     &              (
823     &              distb(j,k,iaer,idomain,gausind) * pi *
824     &              radGAUSb(gausind,iaer,idomain) *
825     &              radGAUSb(gausind,iaer,idomain) +
826     &              dista(j,k,iaer,idomain,gausind) * pi *
827     &              radGAUSa(gausind,iaer,idomain) *
828     &              radGAUSa(gausind,iaer,idomain)
829     &              )
830                ENDDO
831                IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
832c                 2.3.3.vis Initialization
833                  qsqrefVISgrid(j,k,:,iaer)=0.
834                  qextVISgrid(j,k,:,iaer)=0.
835                  qscatVISgrid(j,k,:,iaer)=0.
836                  omegVISgrid(j,k,:,iaer)=0.
837                  gVISgrid(j,k,:,iaer)=0.
838                  qrefVISgrid(j,k,iaer)=0.
839                  qscatrefVISgrid(j,k,iaer)=0.
840                  omegrefVISgrid(j,k,iaer)=0.
841
842                  DO gausind=1,ngau
843                    DO m=1,nsun
844c                     Convolution:
845                      qextVISgrid(j,k,m,iaer) =
846     &                  qextVISgrid(j,k,m,iaer) +
847     &                  weightgaus(gausind) *
848     &                  (
849     &                  qsqrefVISb(m,gausind,iaer) *
850     &                  qrefVISb(gausind,iaer) *
851     &                  pi*radGAUSb(gausind,iaer,idomain) *
852     &                  radGAUSb(gausind,iaer,idomain) *
853     &                  distb(j,k,iaer,idomain,gausind) +
854     &                  qsqrefVISa(m,gausind,iaer) *
855     &                  qrefVISa(gausind,iaer) *
856     &                  pi*radGAUSa(gausind,iaer,idomain) *
857     &                  radGAUSa(gausind,iaer,idomain) *
858     &                  dista(j,k,iaer,idomain,gausind)
859     &                  )
860                      qscatVISgrid(j,k,m,iaer) =
861     &                  qscatVISgrid(j,k,m,iaer) +
862     &                  weightgaus(gausind) *
863     &                  (
864     &                  omegVISb(m,gausind,iaer) *
865     &                  qsqrefVISb(m,gausind,iaer) *
866     &                  qrefVISb(gausind,iaer) *
867     &                  pi*radGAUSb(gausind,iaer,idomain) *
868     &                  radGAUSb(gausind,iaer,idomain) *
869     &                  distb(j,k,iaer,idomain,gausind) +
870     &                  omegVISa(m,gausind,iaer) *
871     &                  qsqrefVISa(m,gausind,iaer) *
872     &                  qrefVISa(gausind,iaer) *
873     &                  pi*radGAUSa(gausind,iaer,idomain) *
874     &                  radGAUSa(gausind,iaer,idomain) *
875     &                  dista(j,k,iaer,idomain,gausind)
876     &                  )
877                      gVISgrid(j,k,m,iaer) =
878     &                  gVISgrid(j,k,m,iaer) +
879     &                  weightgaus(gausind) *
880     &                  (
881     &                  omegVISb(m,gausind,iaer) *
882     &                  qsqrefVISb(m,gausind,iaer) *
883     &                  qrefVISb(gausind,iaer) *
884     &                  gVISb(m,gausind,iaer) *
885     &                  pi*radGAUSb(gausind,iaer,idomain) *
886     &                  radGAUSb(gausind,iaer,idomain) *
887     &                  distb(j,k,iaer,idomain,gausind) +
888     &                  omegVISa(m,gausind,iaer) *
889     &                  qsqrefVISa(m,gausind,iaer) *
890     &                  qrefVISa(gausind,iaer) *
891     &                  gVISa(m,gausind,iaer) *
892     &                  pi*radGAUSa(gausind,iaer,idomain) *
893     &                  radGAUSa(gausind,iaer,idomain) *
894     &                  dista(j,k,iaer,idomain,gausind)
895     &                  )
896                    ENDDO
897                    qrefVISgrid(j,k,iaer) =
898     &                qrefVISgrid(j,k,iaer) +
899     &                weightgaus(gausind) *
900     &                (
901     &                qrefVISb(gausind,iaer) *
902     &                pi*radGAUSb(gausind,iaer,idomain) *
903     &                radGAUSb(gausind,iaer,idomain) *
904     &                distb(j,k,iaer,idomain,gausind) +
905     &                qrefVISa(gausind,iaer) *
906     &                pi*radGAUSa(gausind,iaer,idomain) *
907     &                radGAUSa(gausind,iaer,idomain) *
908     &                dista(j,k,iaer,idomain,gausind)
909     &                )
910                    qscatrefVISgrid(j,k,iaer) =
911     &                qscatrefVISgrid(j,k,iaer) +
912     &                weightgaus(gausind) *
913     &                (
914     &                omegrefVISb(gausind,iaer) *
915     &                qrefVISb(gausind,iaer) *
916     &                pi*radGAUSb(gausind,iaer,idomain) *
917     &                radGAUSb(gausind,iaer,idomain) *
918     &                distb(j,k,iaer,idomain,gausind) +
919     &                omegrefVISa(gausind,iaer) *
920     &                qrefVISa(gausind,iaer) *
921     &                pi*radGAUSa(gausind,iaer,idomain) *
922     &                radGAUSa(gausind,iaer,idomain) *
923     &                dista(j,k,iaer,idomain,gausind)
924     &                )
925                  ENDDO
926                  qrefVISgrid(j,k,iaer)=qrefVISgrid(j,k,iaer) /
927     &                          normd(j,k,iaer,idomain)
928                  qscatrefVISgrid(j,k,iaer)=qscatrefVISgrid(j,k,iaer) /
929     &                          normd(j,k,iaer,idomain)
930                  omegrefVISgrid(j,k,iaer)=qscatrefVISgrid(j,k,iaer) /
931     &                         qrefVISgrid(j,k,iaer)
932                  DO m=1,nsun
933                    qextVISgrid(j,k,m,iaer)=qextVISgrid(j,k,m,iaer) /
934     &                          normd(j,k,iaer,idomain)
935                    qscatVISgrid(j,k,m,iaer)=qscatVISgrid(j,k,m,iaer) /
936     &                          normd(j,k,iaer,idomain)
937                    gVISgrid(j,k,m,iaer)=gVISgrid(j,k,m,iaer) /
938     &                          qscatVISgrid(j,k,m,iaer) /
939     &                          normd(j,k,iaer,idomain)
940
941                    qsqrefVISgrid(j,k,m,iaer)=qextVISgrid(j,k,m,iaer) /
942     &                          qrefVISgrid(j,k,iaer)
943                    omegVISgrid(j,k,m,iaer)=qscatVISgrid(j,k,m,iaer) /
944     &                          qextVISgrid(j,k,m,iaer)
945                  ENDDO
946                ELSE                   ! INFRARED DOMAIN ----------
947c                 2.3.3.ir Initialization
948                  qsqrefIRgrid(j,k,:,iaer)=0.
949                  qextIRgrid(j,k,:,iaer)=0.
950                  qscatIRgrid(j,k,:,iaer)=0.
951                  omegIRgrid(j,k,:,iaer)=0.
952                  gIRgrid(j,k,:,iaer)=0.
953                  qrefIRgrid(j,k,iaer)=0.
954                  qscatrefIRgrid(j,k,iaer)=0.
955                  omegrefIRgrid(j,k,iaer)=0.
956
957                  DO gausind=1,ngau
958                    DO m=1,nir
959c                     Convolution:
960                      qextIRgrid(j,k,m,iaer) =
961     &                  qextIRgrid(j,k,m,iaer) +
962     &                  weightgaus(gausind) *
963     &                  (
964     &                  qsqrefIRb(m,gausind,iaer) *
965     &                  qrefVISb(gausind,iaer) *
966     &                  pi*radGAUSb(gausind,iaer,idomain) *
967     &                  radGAUSb(gausind,iaer,idomain) *
968     &                  distb(j,k,iaer,idomain,gausind) +
969     &                  qsqrefIRa(m,gausind,iaer) *
970     &                  qrefVISa(gausind,iaer) *
971     &                  pi*radGAUSa(gausind,iaer,idomain) *
972     &                  radGAUSa(gausind,iaer,idomain) *
973     &                  dista(j,k,iaer,idomain,gausind)
974     &                  )
975                      qscatIRgrid(j,k,m,iaer) =
976     &                  qscatIRgrid(j,k,m,iaer) +
977     &                  weightgaus(gausind) *
978     &                  (
979     &                  omegIRb(m,gausind,iaer) *
980     &                  qsqrefIRb(m,gausind,iaer) *
981     &                  qrefVISb(gausind,iaer) *
982     &                  pi*radGAUSb(gausind,iaer,idomain) *
983     &                  radGAUSb(gausind,iaer,idomain) *
984     &                  distb(j,k,iaer,idomain,gausind) +
985     &                  omegIRa(m,gausind,iaer) *
986     &                  qsqrefIRa(m,gausind,iaer) *
987     &                  qrefVISa(gausind,iaer) *
988     &                  pi*radGAUSa(gausind,iaer,idomain) *
989     &                  radGAUSa(gausind,iaer,idomain) *
990     &                  dista(j,k,iaer,idomain,gausind)
991     &                  )
992                      gIRgrid(j,k,m,iaer) =
993     &                  gIRgrid(j,k,m,iaer) +
994     &                  weightgaus(gausind) *
995     &                  (
996     &                  omegIRb(m,gausind,iaer) *
997     &                  qsqrefIRb(m,gausind,iaer) *
998     &                  qrefVISb(gausind,iaer) *
999     &                  gIRb(m,gausind,iaer) *
1000     &                  pi*radGAUSb(gausind,iaer,idomain) *
1001     &                  radGAUSb(gausind,iaer,idomain) *
1002     &                  distb(j,k,iaer,idomain,gausind) +
1003     &                  omegIRa(m,gausind,iaer) *
1004     &                  qsqrefIRa(m,gausind,iaer) *
1005     &                  qrefVISa(gausind,iaer) *
1006     &                  gIRa(m,gausind,iaer) *
1007     &                  pi*radGAUSa(gausind,iaer,idomain) *
1008     &                  radGAUSa(gausind,iaer,idomain) *
1009     &                  dista(j,k,iaer,idomain,gausind)
1010     &                  )
1011                    ENDDO
1012                    qrefIRgrid(j,k,iaer) =
1013     &                qrefIRgrid(j,k,iaer) +
1014     &                weightgaus(gausind) *
1015     &                (
1016     &                qrefIRb(gausind,iaer) *
1017     &                pi*radGAUSb(gausind,iaer,idomain) *
1018     &                radGAUSb(gausind,iaer,idomain) *
1019     &                distb(j,k,iaer,idomain,gausind) +
1020     &                qrefIRa(gausind,iaer) *
1021     &                pi*radGAUSa(gausind,iaer,idomain) *
1022     &                radGAUSa(gausind,iaer,idomain) *
1023     &                dista(j,k,iaer,idomain,gausind)
1024     &                )
1025                    qscatrefIRgrid(j,k,iaer) =
1026     &                qscatrefIRgrid(j,k,iaer) +
1027     &                weightgaus(gausind) *
1028     &                (
1029     &                omegrefIRb(gausind,iaer) *
1030     &                qrefIRb(gausind,iaer) *
1031     &                pi*radGAUSb(gausind,iaer,idomain) *
1032     &                radGAUSb(gausind,iaer,idomain) *
1033     &                distb(j,k,iaer,idomain,gausind) +
1034     &                omegrefIRa(gausind,iaer) *
1035     &                qrefIRa(gausind,iaer) *
1036     &                pi*radGAUSa(gausind,iaer,idomain) *
1037     &                radGAUSa(gausind,iaer,idomain) *
1038     &                dista(j,k,iaer,idomain,gausind)
1039     &                )
1040                  ENDDO
1041                  qrefIRgrid(j,k,iaer)=qrefIRgrid(j,k,iaer) /
1042     &                          normd(j,k,iaer,idomain)
1043                  qscatrefIRgrid(j,k,iaer)=qscatrefIRgrid(j,k,iaer) /
1044     &                          normd(j,k,iaer,idomain)
1045                  omegrefIRgrid(j,k,iaer)=qscatrefIRgrid(j,k,iaer) /
1046     &                         qrefIRgrid(j,k,iaer)
1047                  DO m=1,nir
1048                    qextIRgrid(j,k,m,iaer)=qextIRgrid(j,k,m,iaer) /
1049     &                          normd(j,k,iaer,idomain)
1050                    qscatIRgrid(j,k,m,iaer)=qscatIRgrid(j,k,m,iaer) /
1051     &                          normd(j,k,iaer,idomain)
1052                    gIRgrid(j,k,m,iaer)=gIRgrid(j,k,m,iaer) /
1053     &                          qscatIRgrid(j,k,m,iaer) /
1054     &                          normd(j,k,iaer,idomain)
1055
1056                    qsqrefIRgrid(j,k,m,iaer)=qextIRgrid(j,k,m,iaer) /
1057     &                          qrefVISgrid(j,k,iaer)
1058                    omegIRgrid(j,k,m,iaer)=qscatIRgrid(j,k,m,iaer) /
1059     &                          qextIRgrid(j,k,m,iaer)
1060                  ENDDO
1061                ENDIF                  ! --------------------------
1062                checkgrid(j,k,iaer,idomain) = .true.
1063              ENDIF !checkgrid
1064            ENDDO !grid_j
1065          ENDDO !grid_i
1066c         3.4 Bilinear interpolation
1067          k1 = (1-kx)*(1-ky)
1068          k2 = kx*(1-ky)
1069          k3 = kx*ky
1070          k4 = (1-kx)*ky
1071          IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
1072          DO m=1,nsun
1073            QVISsQREF3d(ig,lg,m,iaer) =
1074     &                  k1*qsqrefVISgrid(grid_i,grid_j,m,iaer) +
1075     &                  k2*qsqrefVISgrid(grid_i+1,grid_j,m,iaer) +
1076     &                  k3*qsqrefVISgrid(grid_i+1,grid_j+1,m,iaer) +
1077     &                  k4*qsqrefVISgrid(grid_i,grid_j+1,m,iaer)
1078            omegaVIS3d(ig,lg,m,iaer) =
1079     &                  k1*omegVISgrid(grid_i,grid_j,m,iaer) +
1080     &                  k2*omegVISgrid(grid_i+1,grid_j,m,iaer) +
1081     &                  k3*omegVISgrid(grid_i+1,grid_j+1,m,iaer) +
1082     &                  k4*omegVISgrid(grid_i,grid_j+1,m,iaer)
1083            gVIS3d(ig,lg,m,iaer) =
1084     &                  k1*gVISgrid(grid_i,grid_j,m,iaer) +
1085     &                  k2*gVISgrid(grid_i+1,grid_j,m,iaer) +
1086     &                  k3*gVISgrid(grid_i+1,grid_j+1,m,iaer) +
1087     &                  k4*gVISgrid(grid_i,grid_j+1,m,iaer)
1088          ENDDO !nsun
1089          QREFvis3d(ig,lg,iaer) =
1090     &                  k1*qrefVISgrid(grid_i,grid_j,iaer) +
1091     &                  k2*qrefVISgrid(grid_i+1,grid_j,iaer) +
1092     &                  k3*qrefVISgrid(grid_i+1,grid_j+1,iaer) +
1093     &                  k4*qrefVISgrid(grid_i,grid_j+1,iaer)
1094          omegaREFvis3d(ig,lg,iaer) =
1095     &                  k1*omegrefVISgrid(grid_i,grid_j,iaer) +
1096     &                  k2*omegrefVISgrid(grid_i+1,grid_j,iaer) +
1097     &                  k3*omegrefVISgrid(grid_i+1,grid_j+1,iaer) +
1098     &                  k4*omegrefVISgrid(grid_i,grid_j+1,iaer)
1099          ELSE                   ! INFRARED -----------------------
1100          DO m=1,nir
1101            QIRsQREF3d(ig,lg,m,iaer) =
1102     &                  k1*qsqrefIRgrid(grid_i,grid_j,m,iaer) +
1103     &                  k2*qsqrefIRgrid(grid_i+1,grid_j,m,iaer) +
1104     &                  k3*qsqrefIRgrid(grid_i+1,grid_j+1,m,iaer) +
1105     &                  k4*qsqrefIRgrid(grid_i,grid_j+1,m,iaer)
1106            omegaIR3d(ig,lg,m,iaer) =
1107     &                  k1*omegIRgrid(grid_i,grid_j,m,iaer) +
1108     &                  k2*omegIRgrid(grid_i+1,grid_j,m,iaer) +
1109     &                  k3*omegIRgrid(grid_i+1,grid_j+1,m,iaer) +
1110     &                  k4*omegIRgrid(grid_i,grid_j+1,m,iaer)
1111            gIR3d(ig,lg,m,iaer) =
1112     &                  k1*gIRgrid(grid_i,grid_j,m,iaer) +
1113     &                  k2*gIRgrid(grid_i+1,grid_j,m,iaer) +
1114     &                  k3*gIRgrid(grid_i+1,grid_j+1,m,iaer) +
1115     &                  k4*gIRgrid(grid_i,grid_j+1,m,iaer)
1116          ENDDO !nir
1117          QREFir3d(ig,lg,iaer) =
1118     &                  k1*qrefIRgrid(grid_i,grid_j,iaer) +
1119     &                  k2*qrefIRgrid(grid_i+1,grid_j,iaer) +
1120     &                  k3*qrefIRgrid(grid_i+1,grid_j+1,iaer) +
1121     &                  k4*qrefIRgrid(grid_i,grid_j+1,iaer)
1122          omegaREFir3d(ig,lg,iaer) =
1123     &                  k1*omegrefIRgrid(grid_i,grid_j,iaer) +
1124     &                  k2*omegrefIRgrid(grid_i+1,grid_j,iaer) +
1125     &                  k3*omegrefIRgrid(grid_i+1,grid_j+1,iaer) +
1126     &                  k4*omegrefIRgrid(grid_i,grid_j+1,iaer)
1127          ENDIF                  ! --------------------------------
1128        ENDDO !nlayermx
1129      ENDDO !ngridmx
1130
1131      ENDIF ! varyingnueff
1132c==================================================================
1133      ENDDO ! idomain
1134
1135      ENDIF ! nsize = 1
1136
1137      ENDDO ! iaer (loop on aerosol kind)
1138
1139c=====Radiative properties - TESTS=================================
1140      IF (out_qwg) THEN
1141      IF (ngrid.NE.1) THEN
1142        DO out_nchannel = 1, 2
1143c     -------------------------------------------------------------
1144          DO lg = 1, nlayer
1145            DO ig = 1, ngrid
1146              out_qext(ig,lg) =
1147     &           QVISsQREF3d(ig,lg,out_nchannel,out_iaer)*
1148     &           QREFvis3d(ig,lg,out_iaer)
1149              out_omeg(ig,lg) =
1150     &           omegaVIS3d(ig,lg,out_nchannel,out_iaer)
1151              out_g(ig,lg)    = gVIS3d(ig,lg,out_nchannel,out_iaer)
1152            ENDDO ! ig
1153          ENDDO ! lg
1154          write(out_str(1:1),'(i1.1)') out_nchannel
1155          call WRITEDIAGFI(ngrid,'qextvis'//out_str,"Ext.efficiency","",
1156     &                     3,out_qext)
1157          call WRITEDIAGFI(ngrid,'omegvis'//out_str,"Sing.Scat.Alb.","",
1158     &                     3,out_omeg)
1159          call WRITEDIAGFI(ngrid,'gvis'//out_str,"Asym.Factor","",
1160     &                     3,out_g)
1161c     -------------------------------------------------------------
1162        ENDDO ! out_nchannel
1163        DO out_nchannel = 2, 4
1164c     -------------------------------------------------------------
1165          DO lg = 1, nlayer
1166            DO ig = 1, ngrid
1167              out_qext(ig,lg) =
1168     &          QIRsQREF3d(ig,lg,out_nchannel,out_iaer)*
1169     &          QREFir3d(ig,lg,out_iaer)
1170              out_omeg(ig,lg) =
1171     &          omegaIR3d(ig,lg,out_nchannel,out_iaer)
1172              out_g(ig,lg)    = gIR3d(ig,lg,out_nchannel,out_iaer)
1173            ENDDO ! ig
1174          ENDDO ! lg
1175          write(out_str(1:1),'(i1.1)') out_nchannel
1176          call WRITEDIAGFI(ngrid,'qextir'//out_str,"Ext.efficiency","",
1177     &                     3,out_qext)
1178          call WRITEDIAGFI(ngrid,'omegir'//out_str,"Sing.Scat.Alb.","",
1179     &                     3,out_omeg)
1180          call WRITEDIAGFI(ngrid,'gir'//out_str,"Asym.Factor","",
1181     &                     3,out_g)
1182c     -------------------------------------------------------------
1183        ENDDO ! out_nchannel
1184        call WRITEDIAGFI(ngrid,"omegvisref","Sing.Scat.Alb.","",
1185     &                   3,omegaREFvis3d(1,1,out_iaer))
1186        call WRITEDIAGFI(ngrid,"omegirref","Sing.Scat.Alb.","",
1187     &                   3,omegaREFir3d(1,1,out_iaer))
1188      ENDIF ! ngrid.EQ.1
1189      ENDIF ! out_qwg
1190c==================================================================
1191
1192      RETURN
1193      END
Note: See TracBrowser for help on using the repository browser.