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

Last change on this file since 1242 was 1212, checked in by aslmd, 11 years ago

LMDZ.MARS + MESOSCALE

A quite major commit, at least for MESOSCALE.
In a word: any ngrid deserves to be free.

  • no need to recompile when changing number of horizontal grid points or number of processors
  • latest version of LMDZ.MARS physics can be used
  • WARNING! Nesting is still yet to be fixed (since r1027)

Also some small bug fixes to LMDZ.MARS.

Changes in LMDZ.MARS

--> fixed a potential bug in thermal plume model because zlmax was computed both in thermcell_main_mars and calltherm_interface... so made it an OUT argument of calltherm_interface. also: changed the name to limz. and added precompiling flags to avoid the use of planetwide in MESOSCALE. in MESOSCALE we just go high enough (nlayer-5) and do not care about computational cost (although we certainly gain from not using MAXVAL).
--> moved allocations upward in inifis. does not change anything for GCM, but make MESOSCALE modifications simpler, and overall make inifis better organized: first allocations, then reading callphys.def file.
--> added precompiling flags around lines that are both useless for MESOSCALE (notably I/O) and recently adapted to parallel computations in the GCM
--> tidied up what is MESOSCALE vs. GCM in surfini

Changes in MESOSCALE

--> changed makemeso to allow dynamically set nx ny nprocs
--> changed makemeso to remove links to Fortran code adapted to parallel GCM and useless for mesoscale
--> changed ngridmx to ngrid in inifis includes

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