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

Last change on this file since 2156 was 1775, checked in by aslmd, 7 years ago

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