source: trunk/LMDZ.GENERIC/libf/phystd/aeroptpropertiesJBnew.F90 @ 135

Last change on this file since 135 was 135, checked in by aslmd, 14 years ago

CHANGEMENT ARBORESCENCE ETAPE 2 -- NON COMPLET

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