source: trunk/mars/libf/phymars/aeroptproperties.F @ 73

Last change on this file since 73 was 38, checked in by emillour, 15 years ago

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