source: trunk/LMDZ.PLUTO.old/libf/phypluto/aeroptproperties.F90 @ 3237

Last change on this file since 3237 was 3175, checked in by emillour, 23 months ago

Pluto PCM:
Add the old Pluto LMDZ for reference (required prior step to making
an LMDZ.PLUTO using the same framework as the other physics packages).
TB+EM

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