source: codes/icosagcm/branches/SATURN_DYNAMICO/LMDZ.COMMON/libf/phystd/aeroptproperties.F90 @ 222

Last change on this file since 222 was 222, checked in by ymipsl, 10 years ago

Creating temporary dynamico/lmdz/saturn branche

YM

File size: 35.7 KB
Line 
1      SUBROUTINE aeroptproperties(ngrid,nlayer,reffrad,nueffrad,   &
2                                  QVISsQREF3d,omegaVIS3d,gVIS3d,   &
3                                  QIRsQREF3d,omegaIR3d,gIR3d,      &
4                                  QREFvis3d,QREFir3d)!,            &
5!                                  omegaREFvis3d,omegaREFir3d)
6
7      use radinc_h,    only: L_NSPECTI,L_NSPECTV,nsizemax,naerkind
8      use radcommon_h, only: QVISsQREF,omegavis,gvis,QIRsQREF,omegair,gir
9      use radcommon_h, only: qrefvis,qrefir,omegarefvis,omegarefir
10      use radcommon_h, only: radiustab,nsize
11
12      implicit none
13
14!     =============================================================
15!     Aerosol Optical Properties
16!
17!     Description:
18!       Compute the scattering parameters in each grid
19!       box, depending on aerosol grain sizes. Log-normal size
20!       distribution and Gauss-Legendre integration are used.
21
22!     Parameters:
23!       Don't forget to set the value of varyingnueff below; If
24!       the effective variance of the distribution for the given
25!       aerosol is considered homogeneous in the atmosphere, please
26!       set varyingnueff(iaer) to .false. Resulting computational
27!       time will be much better.
28
29!     Authors: J.-B. Madeleine, F. Forget, F. Montmessin
30!     Slightly modified and converted to F90 by R. Wordsworth (2009)
31!     Varying nueff section removed by R. Wordsworth for simplicity
32!     ==============================================================
33
34#include "dimensions.h"
35#include "dimphys.h"
36#include "callkeys.h"
37
38!     Local variables
39!     ---------------
40
41
42
43!     =============================================================
44      LOGICAL, PARAMETER :: varyingnueff(naerkind) = .false.
45!     =============================================================
46
47!     Min. and max radius of the interpolation grid (in METERS)
48      REAL, PARAMETER :: refftabmin = 2e-8 !2e-8
49!      REAL, PARAMETER :: refftabmax = 35e-6
50      REAL, PARAMETER :: refftabmax = 1e-3
51!     Log of the min and max variance of the interpolation grid
52      REAL, PARAMETER :: nuefftabmin = -4.6
53      REAL, PARAMETER :: nuefftabmax = 0.
54!     Number of effective radius of the interpolation grid
55      INTEGER, PARAMETER :: refftabsize = 200
56!     Number of effective variances of the interpolation grid
57!      INTEGER, PARAMETER :: nuefftabsize = 100
58      INTEGER, PARAMETER :: nuefftabsize = 1
59!     Interpolation grid indices (reff,nueff)
60      INTEGER :: grid_i,grid_j
61!     Intermediate variable
62      REAL :: var_tmp,var3d_tmp(ngrid,nlayermx)
63!     Bilinear interpolation factors
64      REAL :: kx,ky,k1,k2,k3,k4
65!     Size distribution parameters
66      REAL :: sizedistk1,sizedistk2
67!     Pi!
68      REAL,SAVE :: pi
69!     Variables used by the Gauss-Legendre integration:
70      INTEGER radius_id,gausind
71      REAL kint
72      REAL drad
73      INTEGER, PARAMETER :: ngau = 10
74      REAL weightgaus(ngau),radgaus(ngau)
75      SAVE weightgaus,radgaus
76!      DATA weightgaus/.2955242247,.2692667193,.2190863625,.1494513491,.0666713443/
77!      DATA radgaus/.1488743389,.4333953941,.6794095682,.8650633666,.9739065285/
78      DATA    radgaus/0.07652652113350,0.22778585114165, &
79                      0.37370608871528,0.51086700195146, &
80                      0.63605368072468,0.74633190646476, &
81                      0.83911697181213,0.91223442826796, &
82                      0.96397192726078,0.99312859919241/
83
84      DATA weightgaus/0.15275338723120,0.14917298659407, &
85                      0.14209610937519,0.13168863843930, &
86                      0.11819453196154,0.10193011980823, &
87                      0.08327674160932,0.06267204829828, &
88                      0.04060142982019,0.01761400714091/
89!     Indices
90      INTEGER :: i,j,k,l,m,iaer,idomain
91      INTEGER :: ig,lg,chg
92
93!     Local saved variables
94!     ---------------------
95
96!     Radius axis of the interpolation grid
97      REAL,SAVE :: refftab(refftabsize)
98!     Variance axis of the interpolation grid
99      REAL,SAVE :: nuefftab(nuefftabsize)
100!     Volume ratio of the grid
101      REAL,SAVE :: logvratgrid,vratgrid
102!     Grid used to remember which calculation is done
103      LOGICAL,SAVE :: checkgrid(refftabsize,nuefftabsize,naerkind,2) = .false.
104!     Optical properties of the grid (VISIBLE)
105      REAL,SAVE :: qsqrefVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
106      REAL,SAVE :: qextVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
107      REAL,SAVE :: qscatVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
108      REAL,SAVE :: omegVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
109      REAL,SAVE :: gVISgrid(refftabsize,nuefftabsize,L_NSPECTV,naerkind)
110!     Optical properties of the grid (INFRARED)
111      REAL,SAVE :: qsqrefIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
112      REAL,SAVE :: qextIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
113      REAL,SAVE :: qscatIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
114      REAL,SAVE :: omegIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
115      REAL,SAVE :: gIRgrid(refftabsize,nuefftabsize,L_NSPECTI,naerkind)
116!     Optical properties of the grid (REFERENCE WAVELENGTHS)
117      REAL,SAVE :: qrefVISgrid(refftabsize,nuefftabsize,naerkind)
118      REAL,SAVE :: qscatrefVISgrid(refftabsize,nuefftabsize,naerkind)
119      REAL,SAVE :: qrefIRgrid(refftabsize,nuefftabsize,naerkind)
120      REAL,SAVE :: qscatrefIRgrid(refftabsize,nuefftabsize,naerkind)
121      REAL,SAVE :: omegrefVISgrid(refftabsize,nuefftabsize,naerkind)
122      REAL,SAVE :: omegrefIRgrid(refftabsize,nuefftabsize,naerkind)
123!     Firstcall
124      LOGICAL,SAVE :: firstcall = .true.
125!     Variables used by the Gauss-Legendre integration:
126      REAL,SAVE :: normd(refftabsize,nuefftabsize,naerkind,2)
127      REAL,SAVE :: dista(refftabsize,nuefftabsize,naerkind,2,ngau)
128      REAL,SAVE :: distb(refftabsize,nuefftabsize,naerkind,2,ngau)
129
130      REAL,SAVE :: radGAUSa(ngau,naerkind,2)
131      REAL,SAVE :: radGAUSb(ngau,naerkind,2)
132
133      REAL,SAVE :: qsqrefVISa(L_NSPECTV,ngau,naerkind)
134      REAL,SAVE :: qrefVISa(ngau,naerkind)
135      REAL,SAVE :: qsqrefVISb(L_NSPECTV,ngau,naerkind)
136      REAL,SAVE :: qrefVISb(ngau,naerkind)
137      REAL,SAVE :: omegVISa(L_NSPECTV,ngau,naerkind)
138      REAL,SAVE :: omegrefVISa(ngau,naerkind)
139      REAL,SAVE :: omegVISb(L_NSPECTV,ngau,naerkind)
140      REAL,SAVE :: omegrefVISb(ngau,naerkind)
141      REAL,SAVE :: gVISa(L_NSPECTV,ngau,naerkind)
142      REAL,SAVE :: gVISb(L_NSPECTV,ngau,naerkind)
143
144      REAL,SAVE :: qsqrefIRa(L_NSPECTI,ngau,naerkind)
145      REAL,SAVE :: qrefIRa(ngau,naerkind)
146      REAL,SAVE :: qsqrefIRb(L_NSPECTI,ngau,naerkind)
147      REAL,SAVE :: qrefIRb(ngau,naerkind)
148      REAL,SAVE :: omegIRa(L_NSPECTI,ngau,naerkind)
149      REAL,SAVE :: omegrefIRa(ngau,naerkind)
150      REAL,SAVE :: omegIRb(L_NSPECTI,ngau,naerkind)
151      REAL,SAVE :: omegrefIRb(ngau,naerkind)
152      REAL,SAVE :: gIRa(L_NSPECTI,ngau,naerkind)
153      REAL,SAVE :: gIRb(L_NSPECTI,ngau,naerkind)
154
155      REAL :: radiusm
156      REAL :: radiusr
157
158!     Inputs
159!     ------
160
161      INTEGER :: ngrid,nlayer
162!     Aerosol effective radius used for radiative transfer (meter)
163      REAL :: reffrad(ngrid,nlayermx,naerkind)
164!     Aerosol effective variance used for radiative transfer (n.u.)
165      REAL :: nueffrad(ngrid,nlayermx,naerkind)
166
167!     Outputs
168!     -------
169
170      REAL :: QVISsQREF3d(ngrid,nlayermx,L_NSPECTV,naerkind)
171      REAL :: omegaVIS3d(ngrid,nlayermx,L_NSPECTV,naerkind)
172      REAL :: gVIS3d(ngrid,nlayermx,L_NSPECTV,naerkind)
173
174      REAL :: QIRsQREF3d(ngrid,nlayermx,L_NSPECTI,naerkind)
175      REAL :: omegaIR3d(ngrid,nlayermx,L_NSPECTI,naerkind)
176      REAL :: gIR3d(ngrid,nlayermx,L_NSPECTI,naerkind)
177
178      REAL :: QREFvis3d(ngrid,nlayermx,naerkind)
179      REAL :: QREFir3d(ngrid,nlayermx,naerkind)
180
181      REAL :: omegaREFvis3d(ngrid,nlayermx,naerkind)
182      REAL :: omegaREFir3d(ngrid,nlayermx,naerkind)
183
184      DO iaer = 1, naerkind ! Loop on aerosol kind
185        IF ( (nsize(iaer,1).EQ.1).AND.(nsize(iaer,2).EQ.1) ) THEN
186!==================================================================
187!       If there is one single particle size, optical
188!         properties of the considered aerosol are homogeneous
189          DO lg = 1, nlayer
190            DO ig = 1, ngrid
191              DO chg = 1, L_NSPECTV
192                QVISsQREF3d(ig,lg,chg,iaer)=QVISsQREF(chg,iaer,1)
193                omegaVIS3d(ig,lg,chg,iaer)=omegaVIS(chg,iaer,1)
194                gVIS3d(ig,lg,chg,iaer)=gVIS(chg,iaer,1)
195              ENDDO
196              DO chg = 1, L_NSPECTI
197                QIRsQREF3d(ig,lg,chg,iaer)=QIRsQREF(chg,iaer,1)
198                omegaIR3d(ig,lg,chg,iaer)=omegaIR(chg,iaer,1)
199                gIR3d(ig,lg,chg,iaer)=gIR(chg,iaer,1)
200              ENDDO
201              QREFvis3d(ig,lg,iaer)=QREFvis(iaer,1)
202              QREFir3d(ig,lg,iaer)=QREFir(iaer,1)
203              omegaREFvis3d(ig,lg,iaer)=omegaREFvis(iaer,1)
204              omegaREFir3d(ig,lg,iaer)=omegaREFir(iaer,1)
205            ENDDO
206          ENDDO
207
208
209          if (firstcall) then
210             print*,'Optical prop. of the aerosol are homogenous for:'
211             print*,'iaer = ',iaer
212          endif
213
214        ELSE ! Varying effective radius and variance
215      DO idomain = 1, 2 ! Loop on visible or infrared channel
216!==================================================================
217!     1. Creating the effective radius and variance grid
218!     --------------------------------------------------
219      IF (firstcall) THEN
220
221!       1.1 Pi!
222        pi = 2. * asin(1.e0)
223
224!       1.2 Effective radius
225        refftab(1)    = refftabmin
226        refftab(refftabsize) = refftabmax
227
228        logvratgrid = log(refftabmax/refftabmin) / float(refftabsize-1)*3.
229        vratgrid = exp(logvratgrid)
230
231        do i = 2, refftabsize-1
232          refftab(i) = refftab(i-1)*vratgrid**(1./3.)
233        enddo
234
235!       1.3 Effective variance
236        if(nuefftabsize.eq.1)then ! addded by RDW
237           print*,'Warning: no variance range in aeroptproperties'
238           nuefftab(1)=0.2
239        else
240           do i = 0, nuefftabsize-1
241              nuefftab(i+1) = exp( nuefftabmin + i*(nuefftabmax-nuefftabmin)/(nuefftabsize-1) )
242           enddo
243        endif
244
245        firstcall = .false.
246      ENDIF
247
248!       1.4 Radius middle point and range for Gauss integration
249        radiusm=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) + radiustab(iaer,idomain,1))
250        radiusr=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) - radiustab(iaer,idomain,1))
251
252!       1.5 Interpolating data at the Gauss quadrature points:
253        DO gausind=1,ngau
254          drad=radiusr*radgaus(gausind)
255          radGAUSa(gausind,iaer,idomain)=radiusm-drad
256
257          radius_id=minloc(abs(radiustab(iaer,idomain,:) - (radiusm-drad)),DIM=1)
258          IF ((radiustab(iaer,idomain,radius_id) - (radiusm-drad)).GT.0) THEN
259            radius_id=radius_id-1
260          ENDIF
261          IF (radius_id.GE.nsize(iaer,idomain)) THEN
262            radius_id=nsize(iaer,idomain)-1
263            kint = 1.
264          ELSEIF (radius_id.LT.1) THEN
265            radius_id=1
266            kint = 0.
267          ELSE
268          kint = ( (radiusm-drad) -                             &
269                   radiustab(iaer,idomain,radius_id) ) /        &
270                 ( radiustab(iaer,idomain,radius_id+1) -        &
271                   radiustab(iaer,idomain,radius_id) )
272          ENDIF
273          IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
274            DO m=1,L_NSPECTV
275               qsqrefVISa(m,gausind,iaer)=                      &
276                    (1-kint)*QVISsQREF(m,iaer,radius_id) +      &
277                    kint*QVISsQREF(m,iaer,radius_id+1)
278            omegVISa(m,gausind,iaer)=                           &
279                    (1-kint)*omegaVIS(m,iaer,radius_id) +       &
280                    kint*omegaVIS(m,iaer,radius_id+1)
281            gVISa(m,gausind,iaer)=                              &
282                    (1-kint)*gVIS(m,iaer,radius_id) +           &
283                    kint*gVIS(m,iaer,radius_id+1)
284            ENDDO
285            qrefVISa(gausind,iaer)=                             &
286                    (1-kint)*QREFvis(iaer,radius_id) +          &
287                    kint*QREFvis(iaer,radius_id+1)
288            omegrefVISa(gausind,iaer)=                          &
289                    (1-kint)*omegaREFvis(iaer,radius_id) +      &
290                    kint*omegaREFvis(iaer,radius_id+1)
291          ELSE ! INFRARED DOMAIN ----------------------------------
292            DO m=1,L_NSPECTI
293            qsqrefIRa(m,gausind,iaer)=                          &
294                    (1-kint)*QIRsQREF(m,iaer,radius_id) +       &
295                    kint*QIRsQREF(m,iaer,radius_id+1)
296            omegIRa(m,gausind,iaer)=                            &
297                    (1-kint)*omegaIR(m,iaer,radius_id) +        &
298                    kint*omegaIR(m,iaer,radius_id+1)
299            gIRa(m,gausind,iaer)=                               &
300                    (1-kint)*gIR(m,iaer,radius_id) +            &
301                    kint*gIR(m,iaer,radius_id+1)
302            ENDDO
303            qrefIRa(gausind,iaer)=                              &
304                    (1-kint)*QREFir(iaer,radius_id) +           &
305                    kint*QREFir(iaer,radius_id+1)
306            omegrefIRa(gausind,iaer)=                           &
307                    (1-kint)*omegaREFir(iaer,radius_id) +       &
308                    kint*omegaREFir(iaer,radius_id+1)
309          ENDIF
310        ENDDO
311
312        DO gausind=1,ngau
313          drad=radiusr*radgaus(gausind)
314          radGAUSb(gausind,iaer,idomain)=radiusm+drad
315
316          radius_id=minloc(abs(radiustab(iaer,idomain,:) -      &
317                               (radiusm+drad)),DIM=1)
318          IF ((radiustab(iaer,idomain,radius_id) -              &
319               (radiusm+drad)).GT.0) THEN
320            radius_id=radius_id-1
321          ENDIF
322          IF (radius_id.GE.nsize(iaer,idomain)) THEN
323            radius_id=nsize(iaer,idomain)-1
324            kint = 1.
325          ELSEIF (radius_id.LT.1) THEN
326            radius_id=1
327            kint = 0.
328          ELSE
329            kint = ( (radiusm+drad) -                           &
330                     radiustab(iaer,idomain,radius_id) ) /      &
331                   ( radiustab(iaer,idomain,radius_id+1) -      &
332                     radiustab(iaer,idomain,radius_id) )
333          ENDIF
334          IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
335            DO m=1,L_NSPECTV
336            qsqrefVISb(m,gausind,iaer)=                         &
337                    (1-kint)*QVISsQREF(m,iaer,radius_id) +      &
338                    kint*QVISsQREF(m,iaer,radius_id+1)   
339            omegVISb(m,gausind,iaer)=                           &
340                    (1-kint)*omegaVIS(m,iaer,radius_id) +       &
341                    kint*omegaVIS(m,iaer,radius_id+1)
342            gVISb(m,gausind,iaer)=                              &
343                    (1-kint)*gVIS(m,iaer,radius_id) +           &
344                    kint*gVIS(m,iaer,radius_id+1)
345            ENDDO
346            qrefVISb(gausind,iaer)=                             &
347                    (1-kint)*QREFvis(iaer,radius_id) +          &
348                    kint*QREFvis(iaer,radius_id+1)
349            omegrefVISb(gausind,iaer)=                          &
350                    (1-kint)*omegaREFvis(iaer,radius_id) +      &
351                    kint*omegaREFvis(iaer,radius_id+1)
352          ELSE ! INFRARED DOMAIN ----------------------------------
353            DO m=1,L_NSPECTI
354            qsqrefIRb(m,gausind,iaer)=                          &
355                    (1-kint)*QIRsQREF(m,iaer,radius_id) +       &
356                    kint*QIRsQREF(m,iaer,radius_id+1)
357            omegIRb(m,gausind,iaer)=                            &
358                    (1-kint)*omegaIR(m,iaer,radius_id) +        &
359                    kint*omegaIR(m,iaer,radius_id+1)
360            gIRb(m,gausind,iaer)=                               &
361                    (1-kint)*gIR(m,iaer,radius_id) +            &
362                    kint*gIR(m,iaer,radius_id+1)
363            ENDDO
364            qrefIRb(gausind,iaer)=                              &
365                    (1-kint)*QREFir(iaer,radius_id) +           &
366                    kint*QREFir(iaer,radius_id+1)
367            omegrefIRb(gausind,iaer)=                           &
368                    (1-kint)*omegaREFir(iaer,radius_id) +       &
369                    kint*omegaREFir(iaer,radius_id+1)
370          ENDIF
371        ENDDO
372
373!==================================================================
374! CONSTANT NUEFF FROM HERE
375!==================================================================
376
377!     2. Compute the scattering parameters using linear
378!       interpolation over grain sizes and constant nueff
379!     ---------------------------------------------------
380
381      DO lg = 1,nlayer
382        DO ig = 1, ngrid
383!         2.1 Effective radius index and kx calculation
384          var_tmp=reffrad(ig,lg,iaer)/refftabmin
385          var_tmp=log(var_tmp)*3.
386          var_tmp=var_tmp/logvratgrid+1.
387          grid_i=floor(var_tmp)
388          IF (grid_i.GE.refftabsize) THEN
389!           WRITE(*,*) 'Warning: particle size in grid box #'
390!           WRITE(*,*) ig,' is too large to be used by the '
391!           WRITE(*,*) 'radiative transfer; please extend the '
392!           WRITE(*,*) 'interpolation grid to larger grain sizes.'
393            grid_i=refftabsize-1
394            kx = 1.
395          ELSEIF (grid_i.LT.1) THEN
396!           WRITE(*,*) 'Warning: particle size in grid box #'
397!           WRITE(*,*) ig,' is too small to be used by the '
398!           WRITE(*,*) 'radiative transfer; please extend the '
399!           WRITE(*,*) 'interpolation grid to smaller grain sizes.'
400            grid_i=1
401            kx = 0.
402          ELSE
403            kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) /            &
404                 ( refftab(grid_i+1)-refftab(grid_i) )
405          ENDIF
406!         2.3 Integration
407          DO j=grid_i,grid_i+1
408!             2.3.1 Check if the calculation has been done
409              IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN
410!               2.3.2 Log-normal dist., r_g and sigma_g are defined
411!                 in [hansen_1974], "Light scattering in planetary
412!                 atmospheres", Space Science Reviews 16 527-610.
413!                 Here, sizedistk1=r_g and sizedistk2=sigma_g^2
414                sizedistk2 = log(1.+nueffrad(1,1,iaer))
415                sizedistk1 = exp(2.5*sizedistk2)
416                sizedistk1 = refftab(j) / sizedistk1
417
418                normd(j,1,iaer,idomain) = 1e-30
419                DO gausind=1,ngau
420                  drad=radiusr*radgaus(gausind)
421                  dista(j,1,iaer,idomain,gausind) = LOG((radiusm-drad)/sizedistk1)
422                  dista(j,1,iaer,idomain,gausind) =                   &
423                    EXP(-dista(j,1,iaer,idomain,gausind) *            &
424                    dista(j,1,iaer,idomain,gausind) *                 &
425                    0.5e0/sizedistk2)/(radiusm-drad)                 
426                  dista(j,1,iaer,idomain,gausind) =                   &
427                    dista(j,1,iaer,idomain,gausind) /                 &
428                    (sqrt(2e0*pi*sizedistk2))
429
430                  distb(j,1,iaer,idomain,gausind) = LOG((radiusm+drad)/sizedistk1)
431                  distb(j,1,iaer,idomain,gausind) =                   &
432                    EXP(-distb(j,1,iaer,idomain,gausind) *            &
433                    distb(j,1,iaer,idomain,gausind) *                 &
434                    0.5e0/sizedistk2)/(radiusm+drad)
435                  distb(j,1,iaer,idomain,gausind) =                   &
436                    distb(j,1,iaer,idomain,gausind) /                 &
437                    (sqrt(2e0*pi*sizedistk2))
438
439                  normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) +   &
440                    weightgaus(gausind) *                             &
441                    (                                                 &
442                    distb(j,1,iaer,idomain,gausind) * pi *            &
443                    radGAUSb(gausind,iaer,idomain) *                  &
444                    radGAUSb(gausind,iaer,idomain) +                  &
445                    dista(j,1,iaer,idomain,gausind) * pi *            &
446                    radGAUSa(gausind,iaer,idomain) *                  &
447                    radGAUSa(gausind,iaer,idomain)                    &
448                    )
449                ENDDO
450                IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
451!                 2.3.3.vis Initialization
452                  qsqrefVISgrid(j,1,:,iaer)=0.
453                  qextVISgrid(j,1,:,iaer)=0.
454                  qscatVISgrid(j,1,:,iaer)=0.
455                  omegVISgrid(j,1,:,iaer)=0.
456                  gVISgrid(j,1,:,iaer)=0.
457                  qrefVISgrid(j,1,iaer)=0.
458                  qscatrefVISgrid(j,1,iaer)=0.
459                  omegrefVISgrid(j,1,iaer)=0.
460
461                  DO gausind=1,ngau
462                    DO m=1,L_NSPECTV
463!                     Convolution:
464                      qextVISgrid(j,1,m,iaer) =              &
465                        qextVISgrid(j,1,m,iaer) +            & 
466                        weightgaus(gausind) *                &
467                        (                                    &
468                        qsqrefVISb(m,gausind,iaer) *         &
469                        qrefVISb(gausind,iaer) *             &
470                        pi*radGAUSb(gausind,iaer,idomain) *  &
471                        radGAUSb(gausind,iaer,idomain) *     &
472                        distb(j,1,iaer,idomain,gausind) +    &
473                        qsqrefVISa(m,gausind,iaer) *         &
474                        qrefVISa(gausind,iaer) *             &
475                        pi*radGAUSa(gausind,iaer,idomain) *  &
476                        radGAUSa(gausind,iaer,idomain) *     &
477                        dista(j,1,iaer,idomain,gausind)      &
478                        )
479                      qscatVISgrid(j,1,m,iaer) =             &
480                        qscatVISgrid(j,1,m,iaer) +           &
481                        weightgaus(gausind) *                &
482                        (                                    &
483                        omegVISb(m,gausind,iaer) *           &
484                        qsqrefVISb(m,gausind,iaer) *         &
485                        qrefVISb(gausind,iaer) *             &
486                        pi*radGAUSb(gausind,iaer,idomain) *  &
487                        radGAUSb(gausind,iaer,idomain) *     &
488                        distb(j,1,iaer,idomain,gausind) +    &
489                        omegVISa(m,gausind,iaer) *           &
490                        qsqrefVISa(m,gausind,iaer) *         &
491                        qrefVISa(gausind,iaer) *             &
492                        pi*radGAUSa(gausind,iaer,idomain) *  &
493                        radGAUSa(gausind,iaer,idomain) *     &
494                        dista(j,1,iaer,idomain,gausind)      &
495                        )
496                      gVISgrid(j,1,m,iaer) =                 &
497                        gVISgrid(j,1,m,iaer) +               &
498                        weightgaus(gausind) *                &
499                        (                                    &
500                        omegVISb(m,gausind,iaer) *           &
501                        qsqrefVISb(m,gausind,iaer) *         &
502                        qrefVISb(gausind,iaer) *             &
503                        gVISb(m,gausind,iaer) *              &
504                        pi*radGAUSb(gausind,iaer,idomain) *  &
505                        radGAUSb(gausind,iaer,idomain) *     &
506                        distb(j,1,iaer,idomain,gausind) +    &
507                        omegVISa(m,gausind,iaer) *           &
508                        qsqrefVISa(m,gausind,iaer) *         &
509                        qrefVISa(gausind,iaer) *             &
510                        gVISa(m,gausind,iaer) *              &
511                        pi*radGAUSa(gausind,iaer,idomain) *  &
512                        radGAUSa(gausind,iaer,idomain) *     &
513                        dista(j,1,iaer,idomain,gausind)      &
514                        )
515                    ENDDO
516                    qrefVISgrid(j,1,iaer) =                  &
517                      qrefVISgrid(j,1,iaer) +                &
518                      weightgaus(gausind) *                  &
519                      (                                      &
520                      qrefVISb(gausind,iaer) *               &
521                      pi*radGAUSb(gausind,iaer,idomain) *    &
522                      radGAUSb(gausind,iaer,idomain) *       &
523                      distb(j,1,iaer,idomain,gausind) +      &
524                      qrefVISa(gausind,iaer) *               &
525                      pi*radGAUSa(gausind,iaer,idomain) *    &
526                      radGAUSa(gausind,iaer,idomain) *       &
527                      dista(j,1,iaer,idomain,gausind)        &
528                      )
529                    qscatrefVISgrid(j,1,iaer) =              &
530                      qscatrefVISgrid(j,1,iaer) +            &
531                      weightgaus(gausind) *                  &
532                      (                                      &
533                      omegrefVISb(gausind,iaer) *            &
534                      qrefVISb(gausind,iaer) *               & 
535                      pi*radGAUSb(gausind,iaer,idomain) *    &
536                      radGAUSb(gausind,iaer,idomain) *       &
537                      distb(j,1,iaer,idomain,gausind) +      &
538                      omegrefVISa(gausind,iaer) *            &
539                      qrefVISa(gausind,iaer) *               &
540                      pi*radGAUSa(gausind,iaer,idomain) *    &
541                      radGAUSa(gausind,iaer,idomain) *       &
542                      dista(j,1,iaer,idomain,gausind)        &
543                      )
544                  ENDDO
545
546                  qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) /          &
547                                normd(j,1,iaer,idomain)       
548                  qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /  &
549                                normd(j,1,iaer,idomain)
550                  omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /   &
551                               qrefVISgrid(j,1,iaer)
552                  DO m=1,L_NSPECTV
553                    qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /    &
554                                normd(j,1,iaer,idomain)
555                    qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /  &
556                                normd(j,1,iaer,idomain)
557                    gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) /          &
558                                qscatVISgrid(j,1,m,iaer) /               &
559                                normd(j,1,iaer,idomain)
560
561                    qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /  &
562                                qrefVISgrid(j,1,iaer)
563                    omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /   &
564                                qextVISgrid(j,1,m,iaer)
565                  ENDDO
566                ELSE                   ! INFRARED DOMAIN ----------
567!                 2.3.3.ir Initialization
568                  qsqrefIRgrid(j,1,:,iaer)=0.
569                  qextIRgrid(j,1,:,iaer)=0.
570                  qscatIRgrid(j,1,:,iaer)=0.
571                  omegIRgrid(j,1,:,iaer)=0.
572                  gIRgrid(j,1,:,iaer)=0.
573                  qrefIRgrid(j,1,iaer)=0.
574                  qscatrefIRgrid(j,1,iaer)=0.
575                  omegrefIRgrid(j,1,iaer)=0.
576
577                  DO gausind=1,ngau
578                    DO m=1,L_NSPECTI
579!                     Convolution:
580                      qextIRgrid(j,1,m,iaer) =                  &
581                        qextIRgrid(j,1,m,iaer) +                &
582                        weightgaus(gausind) *                   &
583                        (                                       &
584                        qsqrefIRb(m,gausind,iaer) *             &
585                        qrefVISb(gausind,iaer) *                &
586                        pi*radGAUSb(gausind,iaer,idomain) *     &
587                        radGAUSb(gausind,iaer,idomain) *        &
588                        distb(j,1,iaer,idomain,gausind) +       &
589                        qsqrefIRa(m,gausind,iaer) *             &
590                        qrefVISa(gausind,iaer) *                &
591                        pi*radGAUSa(gausind,iaer,idomain) *     &
592                        radGAUSa(gausind,iaer,idomain) *        &
593                        dista(j,1,iaer,idomain,gausind)         &
594                        )
595                      qscatIRgrid(j,1,m,iaer) =                 &
596                        qscatIRgrid(j,1,m,iaer) +               &
597                        weightgaus(gausind) *                   &
598                        (                                       &
599                        omegIRb(m,gausind,iaer) *               &
600                        qsqrefIRb(m,gausind,iaer) *             &
601                        qrefVISb(gausind,iaer) *                &
602                        pi*radGAUSb(gausind,iaer,idomain) *     &
603                        radGAUSb(gausind,iaer,idomain) *        &
604                        distb(j,1,iaer,idomain,gausind) +       &
605                        omegIRa(m,gausind,iaer) *               &
606                        qsqrefIRa(m,gausind,iaer) *             &
607                        qrefVISa(gausind,iaer) *                &
608                        pi*radGAUSa(gausind,iaer,idomain) *     &
609                        radGAUSa(gausind,iaer,idomain) *        &
610                        dista(j,1,iaer,idomain,gausind)         &
611                        )
612                      gIRgrid(j,1,m,iaer) =                     &
613                        gIRgrid(j,1,m,iaer) +                   &
614                        weightgaus(gausind) *                   &
615                        (                                       &
616                        omegIRb(m,gausind,iaer) *               &
617                        qsqrefIRb(m,gausind,iaer) *             &
618                        qrefVISb(gausind,iaer) *                &
619                        gIRb(m,gausind,iaer) *                  &
620                        pi*radGAUSb(gausind,iaer,idomain) *     &
621                        radGAUSb(gausind,iaer,idomain) *        &
622                        distb(j,1,iaer,idomain,gausind) +       &
623                        omegIRa(m,gausind,iaer) *               &
624                        qsqrefIRa(m,gausind,iaer) *             &
625                        qrefVISa(gausind,iaer) *                &
626                        gIRa(m,gausind,iaer) *                  &
627                        pi*radGAUSa(gausind,iaer,idomain) *     &
628                        radGAUSa(gausind,iaer,idomain) *        &
629                        dista(j,1,iaer,idomain,gausind)         &
630                        )
631                    ENDDO
632                    qrefIRgrid(j,1,iaer) =                      &
633                      qrefIRgrid(j,1,iaer) +                    &
634                      weightgaus(gausind) *                     &
635                      (                                         &
636                      qrefIRb(gausind,iaer) *                   &
637                      pi*radGAUSb(gausind,iaer,idomain) *       &
638                      radGAUSb(gausind,iaer,idomain) *          &
639                      distb(j,1,iaer,idomain,gausind) +         &
640                      qrefIRa(gausind,iaer) *                   &
641                      pi*radGAUSa(gausind,iaer,idomain) *       &
642                      radGAUSa(gausind,iaer,idomain) *          &
643                      dista(j,1,iaer,idomain,gausind)           &
644                      )
645                    qscatrefIRgrid(j,1,iaer) =                  &
646                      qscatrefIRgrid(j,1,iaer) +                &
647                      weightgaus(gausind) *                     &
648                      (                                         &
649                      omegrefIRb(gausind,iaer) *                &
650                      qrefIRb(gausind,iaer) *                   &
651                      pi*radGAUSb(gausind,iaer,idomain) *       &
652                      radGAUSb(gausind,iaer,idomain) *          &
653                      distb(j,1,iaer,idomain,gausind) +         &
654                      omegrefIRa(gausind,iaer) *                &
655                      qrefIRa(gausind,iaer) *                   &
656                      pi*radGAUSa(gausind,iaer,idomain) *       &
657                      radGAUSa(gausind,iaer,idomain) *          &
658                      dista(j,1,iaer,idomain,gausind)           &
659                      )
660                  ENDDO
661 
662                  qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) /          &
663                                normd(j,1,iaer,idomain)
664                  qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /  &
665                                normd(j,1,iaer,idomain)
666                  omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /   &
667                               qrefIRgrid(j,1,iaer)
668                  DO m=1,L_NSPECTI
669                    qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /    &
670                                normd(j,1,iaer,idomain)
671                    qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /  &
672                                normd(j,1,iaer,idomain)
673                    gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) /          &
674                                qscatIRgrid(j,1,m,iaer) /              &
675                                normd(j,1,iaer,idomain)
676
677                    qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /  &
678                                qrefVISgrid(j,1,iaer)
679                    omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /   &
680                                qextIRgrid(j,1,m,iaer)
681                  ENDDO
682                ENDIF                  ! --------------------------
683                checkgrid(j,1,iaer,idomain) = .true.
684              ENDIF !checkgrid
685          ENDDO !grid_i
686!         2.4 Linear interpolation
687          k1 = (1-kx)
688          k2 = kx
689          IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
690          DO m=1,L_NSPECTV
691             QVISsQREF3d(ig,lg,m,iaer) =                           &
692                        k1*qsqrefVISgrid(grid_i,1,m,iaer) +        &
693                        k2*qsqrefVISgrid(grid_i+1,1,m,iaer)
694            omegaVIS3d(ig,lg,m,iaer) =                             &
695                        k1*omegVISgrid(grid_i,1,m,iaer) +          &
696                        k2*omegVISgrid(grid_i+1,1,m,iaer)
697            gVIS3d(ig,lg,m,iaer) =                                 &
698                        k1*gVISgrid(grid_i,1,m,iaer) +             &
699                        k2*gVISgrid(grid_i+1,1,m,iaer)
700          ENDDO !L_NSPECTV
701          QREFvis3d(ig,lg,iaer) =                                  &
702                        k1*qrefVISgrid(grid_i,1,iaer) +            &
703                        k2*qrefVISgrid(grid_i+1,1,iaer)
704          omegaREFvis3d(ig,lg,iaer) =                              &
705                        k1*omegrefVISgrid(grid_i,1,iaer) +         &
706                        k2*omegrefVISgrid(grid_i+1,1,iaer)
707          ELSE                   ! INFRARED -----------------------
708          DO m=1,L_NSPECTI
709            QIRsQREF3d(ig,lg,m,iaer) =                             &
710                        k1*qsqrefIRgrid(grid_i,1,m,iaer) +         &
711                        k2*qsqrefIRgrid(grid_i+1,1,m,iaer)
712            omegaIR3d(ig,lg,m,iaer) =                              &
713                        k1*omegIRgrid(grid_i,1,m,iaer) +           &
714                        k2*omegIRgrid(grid_i+1,1,m,iaer) 
715            gIR3d(ig,lg,m,iaer) =                                  & 
716                        k1*gIRgrid(grid_i,1,m,iaer) +              &
717                        k2*gIRgrid(grid_i+1,1,m,iaer)
718          ENDDO !L_NSPECTI
719          QREFir3d(ig,lg,iaer) =                                   &
720                        k1*qrefIRgrid(grid_i,1,iaer) +             &
721                        k2*qrefIRgrid(grid_i+1,1,iaer)
722          omegaREFir3d(ig,lg,iaer) =                               &
723                        k1*omegrefIRgrid(grid_i,1,iaer) +          &
724                        k2*omegrefIRgrid(grid_i+1,1,iaer)
725          ENDIF                  ! --------------------------------
726        ENDDO !nlayermx
727      ENDDO !ngrid
728
729!==================================================================
730
731
732
733      ENDDO ! idomain
734
735      ENDIF ! nsize = 1
736
737      ENDDO ! iaer (loop on aerosol kind)
738
739      RETURN
740    END SUBROUTINE aeroptproperties
741
742
743
744     
Note: See TracBrowser for help on using the repository browser.