source: trunk/LMDZ.GENERIC/libf/phystd/aeroptproperties.F90 @ 1477

Last change on this file since 1477 was 1397, checked in by milmd, 10 years ago

In LMDZ.GENERIC replacement of all phystd .h files by module files.

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