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

Last change on this file since 2176 was 2033, checked in by emillour, 6 years ago

Generic GCM:

  • Cleanup around aeroptproperties.90, remove arrays that were not initialized but still used, probably for nothing since removing them leads to no change in results. Possibly to revisit and further clean up later.

EM

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