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

Last change on this file since 3436 was 3353, checked in by afalco, 6 months ago

Pluto PCM:
Added zrecast & old sedim ;
Choose haze file ;
AF

File size: 38.2 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).LT.minrad).OR.(MAXVAL(reffrad).GT.maxrad)) then
260           WRITE(*,*) 'Warning: particle size in grid box #'
261           WRITE(*,*) ig,' is larger than optical properties. '
262           WRITE(*,*) 'reffrad=',MINVAL(reffrad),'-',MAXVAL(reffrad)
263           WRITE(*,*) 'radiustab=',minrad,'-',maxrad
264
265           ! ensure reffrad is within bounds of radiustab
266          !  WHERE(reffrad.LT.minrad)
267          !     reffrad=minrad
268          !  ENDWHERE
269          !  WHERE(reffrad.GT.maxrad)
270          !     reffrad=maxrad
271          !  ENDWHERE
272          !  WRITE(*,*) 'Truncated reffrad within radiustab bounds:'
273          !  WRITE(*,*) 'New reffrad=',MINVAL(reffrad),'-',MAXVAL(reffrad)
274        ENDIF
275
276
277      ENDIF
278
279! refftab est un tableau de rayon donne par les parametres refftabmin et
280! max de cette routine. cest la grille sur laquelle on va interpoler
281! Idem pour nuefftab
282! radiustab est le tableau de rayon du .dat (1,1,nb rayon)
283
284!       1.4 Radius middle point and range for Gauss integration
285        radiusm=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) + radiustab(iaer,idomain,1)) ! moyenne entre 1er rayon et rayons
286        radiusr=0.5*(radiustab(iaer,idomain,nsize(iaer,idomain)) - radiustab(iaer,idomain,1)) ! ecart avec premier rayon / 2
287
288!       1.5 Interpolating data at the Gauss quadrature points:
289        ! sur [-1:0]
290        DO gausind=1,ngau
291          drad=radiusr*radgaus(gausind)
292          radGAUSa(gausind,iaer,idomain)=radiusm-drad
293
294          ! formule 2.33 dans these JB
295          ! il ne faut pas tomber dans cas change1 ou change2 : il faut
296          ! un tableau de donnee dune certaine taille... bug!
297          radius_id=minloc(abs(radiustab(iaer,idomain,:) - (radiusm-drad)),DIM=1)
298          IF ((radiustab(iaer,idomain,radius_id) - (radiusm-drad)).GT.0) THEN
299            radius_id=radius_id-1
300          ENDIF
301          IF (radius_id.GE.nsize(iaer,idomain)) THEN
302            radius_id=nsize(iaer,idomain)-1
303            kint = 1.
304          ELSEIF (radius_id.LT.1) THEN
305            radius_id=1
306            kint = 0.
307          ELSE
308            kint = ( (radiusm-drad) -                           &
309                   radiustab(iaer,idomain,radius_id) ) /        &
310                 ( radiustab(iaer,idomain,radius_id+1) -        &
311                   radiustab(iaer,idomain,radius_id) )
312          ENDIF
313          IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
314            DO m=1,L_NSPECTV
315               qsqrefVISa(m,gausind,iaer)=                      &
316                    (1-kint)*QVISsQREF(m,iaer,radius_id) +      &
317                    kint*QVISsQREF(m,iaer,radius_id+1)
318            omegVISa(m,gausind,iaer)=                           &
319                    (1-kint)*omegaVIS(m,iaer,radius_id) +       &
320                    kint*omegaVIS(m,iaer,radius_id+1)
321            gVISa(m,gausind,iaer)=                              &
322                    (1-kint)*gVIS(m,iaer,radius_id) +           &
323                    kint*gVIS(m,iaer,radius_id+1)
324            ENDDO
325            qrefVISa(gausind,iaer)=                             &
326                    (1-kint)*QREFvis(iaer,radius_id) +          &
327                    kint*QREFvis(iaer,radius_id+1)
328            omegrefVISa(gausind,iaer)=                          &
329                    (1-kint)*omegaREFvis(iaer,radius_id) +      &
330                    kint*omegaREFvis(iaer,radius_id+1)
331          ELSE ! INFRARED DOMAIN ----------------------------------
332            DO m=1,L_NSPECTI
333            qsqrefIRa(m,gausind,iaer)=                          &
334                    (1-kint)*QIRsQREF(m,iaer,radius_id) +       &
335                    kint*QIRsQREF(m,iaer,radius_id+1)
336            omegIRa(m,gausind,iaer)=                            &
337                    (1-kint)*omegaIR(m,iaer,radius_id) +        &
338                    kint*omegaIR(m,iaer,radius_id+1)
339            gIRa(m,gausind,iaer)=                               &
340                    (1-kint)*gIR(m,iaer,radius_id) +            &
341                    kint*gIR(m,iaer,radius_id+1)
342            ENDDO
343            qrefIRa(gausind,iaer)=                              &
344                    (1-kint)*QREFir(iaer,radius_id) +           &
345                    kint*QREFir(iaer,radius_id+1)
346            omegrefIRa(gausind,iaer)=                           &
347                    (1-kint)*omegaREFir(iaer,radius_id) +       &
348                    kint*omegaREFir(iaer,radius_id+1)
349          ENDIF
350        ENDDO
351
352        ! sur [0:1] radgaus
353        DO gausind=1,ngau
354          drad=radiusr*radgaus(gausind)
355          radGAUSb(gausind,iaer,idomain)=radiusm+drad
356
357          radius_id=minloc(abs(radiustab(iaer,idomain,:) -      &
358                               (radiusm+drad)),DIM=1)
359          IF ((radiustab(iaer,idomain,radius_id) -              &
360               (radiusm+drad)).GT.0) THEN
361            radius_id=radius_id-1
362          ENDIF
363          IF (radius_id.GE.nsize(iaer,idomain)) THEN
364            radius_id=nsize(iaer,idomain)-1
365            kint = 1.
366          ELSEIF (radius_id.LT.1) THEN
367            radius_id=1
368            kint = 0.
369          ELSE
370            kint = ( (radiusm+drad) -                           &
371                     radiustab(iaer,idomain,radius_id) ) /      &
372                   ( radiustab(iaer,idomain,radius_id+1) -      &
373                     radiustab(iaer,idomain,radius_id) )
374          ENDIF
375          IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------------
376            DO m=1,L_NSPECTV
377            qsqrefVISb(m,gausind,iaer)=                         &
378                    (1-kint)*QVISsQREF(m,iaer,radius_id) +      &
379                    kint*QVISsQREF(m,iaer,radius_id+1)
380            omegVISb(m,gausind,iaer)=                           &
381                    (1-kint)*omegaVIS(m,iaer,radius_id) +       &
382                    kint*omegaVIS(m,iaer,radius_id+1)
383            gVISb(m,gausind,iaer)=                              &
384                    (1-kint)*gVIS(m,iaer,radius_id) +           &
385                    kint*gVIS(m,iaer,radius_id+1)
386            ENDDO
387            qrefVISb(gausind,iaer)=                             &
388                    (1-kint)*QREFvis(iaer,radius_id) +          &
389                    kint*QREFvis(iaer,radius_id+1)
390            omegrefVISb(gausind,iaer)=                          &
391                    (1-kint)*omegaREFvis(iaer,radius_id) +      &
392                    kint*omegaREFvis(iaer,radius_id+1)
393          ELSE ! INFRARED DOMAIN ----------------------------------
394            DO m=1,L_NSPECTI
395            qsqrefIRb(m,gausind,iaer)=                          &
396                    (1-kint)*QIRsQREF(m,iaer,radius_id) +       &
397                    kint*QIRsQREF(m,iaer,radius_id+1)
398            omegIRb(m,gausind,iaer)=                            &
399                    (1-kint)*omegaIR(m,iaer,radius_id) +        &
400                    kint*omegaIR(m,iaer,radius_id+1)
401            gIRb(m,gausind,iaer)=                               &
402                    (1-kint)*gIR(m,iaer,radius_id) +            &
403                    kint*gIR(m,iaer,radius_id+1)
404            ENDDO
405            qrefIRb(gausind,iaer)=                              &
406                    (1-kint)*QREFir(iaer,radius_id) +           &
407                    kint*QREFir(iaer,radius_id+1)
408            omegrefIRb(gausind,iaer)=                           &
409                    (1-kint)*omegaREFir(iaer,radius_id) +       &
410                    kint*omegaREFir(iaer,radius_id+1)
411          ENDIF
412        ENDDO
413
414!==================================================================
415! CONSTANT NUEFF FROM HERE
416!==================================================================
417
418!     2. Compute the scattering parameters using linear
419!       interpolation over grain sizes and constant nueff
420!     ---------------------------------------------------
421      DO lg = 1,nlayer
422        DO ig = 1, ngrid
423!         2.1 Effective radius index and kx calculation
424          var_tmp=reffrad(ig,lg,iaer)/refftabmin
425          var_tmp=log(var_tmp)*3.
426          var_tmp=var_tmp/logvratgrid+1.
427          grid_i=floor(var_tmp)
428!         grid_i: get index of reffrad in refftab
429
430          IF (grid_i.GE.refftabsize) THEN
431           WRITE(*,*) 'Warning: particle size in grid box #'
432           WRITE(*,*) ig,' is too large to be used by the '
433           WRITE(*,*) 'radiative transfer; please extend the '
434           WRITE(*,*) 'interpolation grid to larger grain sizes.'
435           grid_i=refftabsize-1
436           kx = 1.
437           stop
438          ELSEIF (grid_i.LE.1) THEN
439           WRITE(*,*) 'Warning: particle size in grid box #'
440           WRITE(*,*) ig,' is too small to be used by the '
441           WRITE(*,*) 'radiative transfer; please extend the '
442           WRITE(*,*) 'interpolation grid to smaller grain sizes.'
443           grid_i=1
444           kx = 0.
445           stop
446          ELSE
447            kx = ( reffrad(ig,lg,iaer)-refftab(grid_i) ) /            &
448                 ( refftab(grid_i+1)-refftab(grid_i) )
449            ! kx is the coeff of interpolation between the 2 radius adjacent to reffrad in refftab
450          ENDIF
451!         2.3 Integration
452          DO j=grid_i,grid_i+1
453!             2.3.1 Check if the calculation has been done
454              IF (.NOT.checkgrid(j,1,iaer,idomain)) THEN
455!               2.3.2 Log-normal dist., r_g and sigma_g are defined
456!                 in [hansen_1974], "Light scattering in planetary
457!                 atmospheres", Space Science Reviews 16 527-610.
458!                 Here, sizedistk1=r_g and sizedistk2=sigma_g^2
459                sizedistk2 = log(1.+nueffrad(1,1,iaer))
460                sizedistk1 = exp(2.5*sizedistk2)
461                sizedistk1 = refftab(j) / sizedistk1
462
463                normd(j,1,iaer,idomain) = 1e-30
464                DO gausind=1,ngau
465                  drad=radiusr*radgaus(gausind)
466                  dista(j,1,iaer,idomain,gausind) = LOG((radiusm-drad)/sizedistk1)
467                  dista(j,1,iaer,idomain,gausind) =                   &
468                    EXP(-dista(j,1,iaer,idomain,gausind) *            &
469                    dista(j,1,iaer,idomain,gausind) *                 &
470                    0.5e0/sizedistk2)/(radiusm-drad)
471                  dista(j,1,iaer,idomain,gausind) =                   &
472                    dista(j,1,iaer,idomain,gausind) /                 &
473                    (sqrt(2e0*pi*sizedistk2))
474
475                  distb(j,1,iaer,idomain,gausind) = LOG((radiusm+drad)/sizedistk1)
476                  distb(j,1,iaer,idomain,gausind) =                   &
477                    EXP(-distb(j,1,iaer,idomain,gausind) *            &
478                    distb(j,1,iaer,idomain,gausind) *                 &
479                    0.5e0/sizedistk2)/(radiusm+drad)
480                  distb(j,1,iaer,idomain,gausind) =                   &
481                    distb(j,1,iaer,idomain,gausind) /                 &
482                    (sqrt(2e0*pi*sizedistk2))
483
484                  normd(j,1,iaer,idomain)=normd(j,1,iaer,idomain) +   &
485                    weightgaus(gausind) *                             &
486                    (                                                 &
487                    distb(j,1,iaer,idomain,gausind) * pi *            &
488                    radGAUSb(gausind,iaer,idomain) *                  &
489                    radGAUSb(gausind,iaer,idomain) +                  &
490                    dista(j,1,iaer,idomain,gausind) * pi *            &
491                    radGAUSa(gausind,iaer,idomain) *                  &
492                    radGAUSa(gausind,iaer,idomain)                    &
493                    )
494                ENDDO
495
496                IF (idomain.EQ.1) THEN ! VISIBLE DOMAIN -----------
497!                 2.3.3.vis Initialization
498                  qsqrefVISgrid(j,1,:,iaer)=0.
499                  qextVISgrid(j,1,:,iaer)=0.
500                  qscatVISgrid(j,1,:,iaer)=0.
501                  omegVISgrid(j,1,:,iaer)=0.
502                  gVISgrid(j,1,:,iaer)=0.
503                  qrefVISgrid(j,1,iaer)=0.
504                  qscatrefVISgrid(j,1,iaer)=0.
505                  omegrefVISgrid(j,1,iaer)=0.
506
507                  DO gausind=1,ngau
508                    DO m=1,L_NSPECTV
509!                     Convolution:
510                      qextVISgrid(j,1,m,iaer) =              &
511                        qextVISgrid(j,1,m,iaer) +            &
512                        weightgaus(gausind) *                &
513                        (                                    &
514                        qsqrefVISb(m,gausind,iaer) *         &
515                        qrefVISb(gausind,iaer) *             &
516                        pi*radGAUSb(gausind,iaer,idomain) *  &
517                        radGAUSb(gausind,iaer,idomain) *     &
518                        distb(j,1,iaer,idomain,gausind) +    &
519                        qsqrefVISa(m,gausind,iaer) *         &
520                        qrefVISa(gausind,iaer) *             &
521                        pi*radGAUSa(gausind,iaer,idomain) *  &
522                        radGAUSa(gausind,iaer,idomain) *     &
523                        dista(j,1,iaer,idomain,gausind)      &
524                        )
525                      qscatVISgrid(j,1,m,iaer) =             &
526                        qscatVISgrid(j,1,m,iaer) +           &
527                        weightgaus(gausind) *                &
528                        (                                    &
529                        omegVISb(m,gausind,iaer) *           &
530                        qsqrefVISb(m,gausind,iaer) *         &
531                        qrefVISb(gausind,iaer) *             &
532                        pi*radGAUSb(gausind,iaer,idomain) *  &
533                        radGAUSb(gausind,iaer,idomain) *     &
534                        distb(j,1,iaer,idomain,gausind) +    &
535                        omegVISa(m,gausind,iaer) *           &
536                        qsqrefVISa(m,gausind,iaer) *         &
537                        qrefVISa(gausind,iaer) *             &
538                        pi*radGAUSa(gausind,iaer,idomain) *  &
539                        radGAUSa(gausind,iaer,idomain) *     &
540                        dista(j,1,iaer,idomain,gausind)      &
541                        )
542                      gVISgrid(j,1,m,iaer) =                 &
543                        gVISgrid(j,1,m,iaer) +               &
544                        weightgaus(gausind) *                &
545                        (                                    &
546                        omegVISb(m,gausind,iaer) *           &
547                        qsqrefVISb(m,gausind,iaer) *         &
548                        qrefVISb(gausind,iaer) *             &
549                        gVISb(m,gausind,iaer) *              &
550                        pi*radGAUSb(gausind,iaer,idomain) *  &
551                        radGAUSb(gausind,iaer,idomain) *     &
552                        distb(j,1,iaer,idomain,gausind) +    &
553                        omegVISa(m,gausind,iaer) *           &
554                        qsqrefVISa(m,gausind,iaer) *         &
555                        qrefVISa(gausind,iaer) *             &
556                        gVISa(m,gausind,iaer) *              &
557                        pi*radGAUSa(gausind,iaer,idomain) *  &
558                        radGAUSa(gausind,iaer,idomain) *     &
559                        dista(j,1,iaer,idomain,gausind)      &
560                        )
561                    ENDDO
562                    qrefVISgrid(j,1,iaer) =                  &
563                      qrefVISgrid(j,1,iaer) +                &
564                      weightgaus(gausind) *                  &
565                      (                                      &
566                      qrefVISb(gausind,iaer) *               &
567                      pi*radGAUSb(gausind,iaer,idomain) *    &
568                      radGAUSb(gausind,iaer,idomain) *       &
569                      distb(j,1,iaer,idomain,gausind) +      &
570                      qrefVISa(gausind,iaer) *               &
571                      pi*radGAUSa(gausind,iaer,idomain) *    &
572                      radGAUSa(gausind,iaer,idomain) *       &
573                      dista(j,1,iaer,idomain,gausind)        &
574                      )
575                    qscatrefVISgrid(j,1,iaer) =              &
576                      qscatrefVISgrid(j,1,iaer) +            &
577                      weightgaus(gausind) *                  &
578                      (                                      &
579                      omegrefVISb(gausind,iaer) *            &
580                      qrefVISb(gausind,iaer) *               &
581                      pi*radGAUSb(gausind,iaer,idomain) *    &
582                      radGAUSb(gausind,iaer,idomain) *       &
583                      distb(j,1,iaer,idomain,gausind) +      &
584                      omegrefVISa(gausind,iaer) *            &
585                      qrefVISa(gausind,iaer) *               &
586                      pi*radGAUSa(gausind,iaer,idomain) *    &
587                      radGAUSa(gausind,iaer,idomain) *       &
588                      dista(j,1,iaer,idomain,gausind)        &
589                      )
590                  ENDDO
591
592                  qrefVISgrid(j,1,iaer)=qrefVISgrid(j,1,iaer) /          &
593                                normd(j,1,iaer,idomain)
594                  qscatrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /  &
595                                normd(j,1,iaer,idomain)
596                  omegrefVISgrid(j,1,iaer)=qscatrefVISgrid(j,1,iaer) /   &
597                               qrefVISgrid(j,1,iaer)
598                  DO m=1,L_NSPECTV
599                    qextVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /    &
600                                normd(j,1,iaer,idomain)
601                    qscatVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /  &
602                                normd(j,1,iaer,idomain)
603                    gVISgrid(j,1,m,iaer)=gVISgrid(j,1,m,iaer) /          &
604                                qscatVISgrid(j,1,m,iaer) /               &
605                                normd(j,1,iaer,idomain)
606
607                    qsqrefVISgrid(j,1,m,iaer)=qextVISgrid(j,1,m,iaer) /  &
608                                qrefVISgrid(j,1,iaer)
609                    omegVISgrid(j,1,m,iaer)=qscatVISgrid(j,1,m,iaer) /   &
610                                qextVISgrid(j,1,m,iaer)
611                  ENDDO
612
613                ELSE                   ! INFRARED DOMAIN ----------
614!                 2.3.3.ir Initialization
615                  qsqrefIRgrid(j,1,:,iaer)=0.
616                  qextIRgrid(j,1,:,iaer)=0.
617                  qscatIRgrid(j,1,:,iaer)=0.
618                  omegIRgrid(j,1,:,iaer)=0.
619                  gIRgrid(j,1,:,iaer)=0.
620                  qrefIRgrid(j,1,iaer)=0.
621                  qscatrefIRgrid(j,1,iaer)=0.
622                  omegrefIRgrid(j,1,iaer)=0.
623
624                  DO gausind=1,ngau
625                    DO m=1,L_NSPECTI
626!                     Convolution:
627                      qextIRgrid(j,1,m,iaer) =                  &
628                        qextIRgrid(j,1,m,iaer) +                &
629                        weightgaus(gausind) *                   &
630                        (                                       &
631                        qsqrefIRb(m,gausind,iaer) *             &
632                        qrefVISb(gausind,iaer) *                &
633                        pi*radGAUSb(gausind,iaer,idomain) *     &
634                        radGAUSb(gausind,iaer,idomain) *        &
635                        distb(j,1,iaer,idomain,gausind) +       &
636                        qsqrefIRa(m,gausind,iaer) *             &
637                        qrefVISa(gausind,iaer) *                &
638                        pi*radGAUSa(gausind,iaer,idomain) *     &
639                        radGAUSa(gausind,iaer,idomain) *        &
640                        dista(j,1,iaer,idomain,gausind)         &
641                        )
642                      qscatIRgrid(j,1,m,iaer) =                 &
643                        qscatIRgrid(j,1,m,iaer) +               &
644                        weightgaus(gausind) *                   &
645                        (                                       &
646                        omegIRb(m,gausind,iaer) *               &
647                        qsqrefIRb(m,gausind,iaer) *             &
648                        qrefVISb(gausind,iaer) *                &
649                        pi*radGAUSb(gausind,iaer,idomain) *     &
650                        radGAUSb(gausind,iaer,idomain) *        &
651                        distb(j,1,iaer,idomain,gausind) +       &
652                        omegIRa(m,gausind,iaer) *               &
653                        qsqrefIRa(m,gausind,iaer) *             &
654                        qrefVISa(gausind,iaer) *                &
655                        pi*radGAUSa(gausind,iaer,idomain) *     &
656                        radGAUSa(gausind,iaer,idomain) *        &
657                        dista(j,1,iaer,idomain,gausind)         &
658                        )
659                      gIRgrid(j,1,m,iaer) =                     &
660                        gIRgrid(j,1,m,iaer) +                   &
661                        weightgaus(gausind) *                   &
662                        (                                       &
663                        omegIRb(m,gausind,iaer) *               &
664                        qsqrefIRb(m,gausind,iaer) *             &
665                        qrefVISb(gausind,iaer) *                &
666                        gIRb(m,gausind,iaer) *                  &
667                        pi*radGAUSb(gausind,iaer,idomain) *     &
668                        radGAUSb(gausind,iaer,idomain) *        &
669                        distb(j,1,iaer,idomain,gausind) +       &
670                        omegIRa(m,gausind,iaer) *               &
671                        qsqrefIRa(m,gausind,iaer) *             &
672                        qrefVISa(gausind,iaer) *                &
673                        gIRa(m,gausind,iaer) *                  &
674                        pi*radGAUSa(gausind,iaer,idomain) *     &
675                        radGAUSa(gausind,iaer,idomain) *        &
676                        dista(j,1,iaer,idomain,gausind)         &
677                        )
678                    ENDDO
679                    qrefIRgrid(j,1,iaer) =                      &
680                      qrefIRgrid(j,1,iaer) +                    &
681                      weightgaus(gausind) *                     &
682                      (                                         &
683                      qrefIRb(gausind,iaer) *                   &
684                      pi*radGAUSb(gausind,iaer,idomain) *       &
685                      radGAUSb(gausind,iaer,idomain) *          &
686                      distb(j,1,iaer,idomain,gausind) +         &
687                      qrefIRa(gausind,iaer) *                   &
688                      pi*radGAUSa(gausind,iaer,idomain) *       &
689                      radGAUSa(gausind,iaer,idomain) *          &
690                      dista(j,1,iaer,idomain,gausind)           &
691                      )
692                    qscatrefIRgrid(j,1,iaer) =                  &
693                      qscatrefIRgrid(j,1,iaer) +                &
694                      weightgaus(gausind) *                     &
695                      (                                         &
696                      omegrefIRb(gausind,iaer) *                &
697                      qrefIRb(gausind,iaer) *                   &
698                      pi*radGAUSb(gausind,iaer,idomain) *       &
699                      radGAUSb(gausind,iaer,idomain) *          &
700                      distb(j,1,iaer,idomain,gausind) +         &
701                      omegrefIRa(gausind,iaer) *                &
702                      qrefIRa(gausind,iaer) *                   &
703                      pi*radGAUSa(gausind,iaer,idomain) *       &
704                      radGAUSa(gausind,iaer,idomain) *          &
705                      dista(j,1,iaer,idomain,gausind)           &
706                      )
707                  ENDDO
708
709                  qrefIRgrid(j,1,iaer)=qrefIRgrid(j,1,iaer) /          &
710                                normd(j,1,iaer,idomain)
711                  qscatrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /  &
712                                normd(j,1,iaer,idomain)
713                  omegrefIRgrid(j,1,iaer)=qscatrefIRgrid(j,1,iaer) /   &
714                               qrefIRgrid(j,1,iaer)
715                  DO m=1,L_NSPECTI
716                    qextIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /    &
717                                normd(j,1,iaer,idomain)
718                    qscatIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /  &
719                                normd(j,1,iaer,idomain)
720                    gIRgrid(j,1,m,iaer)=gIRgrid(j,1,m,iaer) /          &
721                                qscatIRgrid(j,1,m,iaer) /              &
722                                normd(j,1,iaer,idomain)
723
724                    qsqrefIRgrid(j,1,m,iaer)=qextIRgrid(j,1,m,iaer) /  &
725                                qrefVISgrid(j,1,iaer)
726                    omegIRgrid(j,1,m,iaer)=qscatIRgrid(j,1,m,iaer) /   &
727                                qextIRgrid(j,1,m,iaer)
728                  ENDDO
729                ENDIF                  ! --------------------------
730                checkgrid(j,1,iaer,idomain) = .true.
731              ENDIF !checkgrid
732          ENDDO !grid_i
733!         2.4 Linear interpolation
734          k1 = (1-kx)
735          k2 = kx
736          IF (idomain.EQ.1) THEN ! VISIBLE ------------------------
737          DO m=1,L_NSPECTV
738             QVISsQREF3d(ig,lg,m,iaer) =                           &
739                        k1*qsqrefVISgrid(grid_i,1,m,iaer) +        &
740                        k2*qsqrefVISgrid(grid_i+1,1,m,iaer)
741            omegaVIS3d(ig,lg,m,iaer) =                             &
742                        k1*omegVISgrid(grid_i,1,m,iaer) +          &
743                        k2*omegVISgrid(grid_i+1,1,m,iaer)
744            gVIS3d(ig,lg,m,iaer) =                                 &
745                        k1*gVISgrid(grid_i,1,m,iaer) +             &
746                        k2*gVISgrid(grid_i+1,1,m,iaer)
747          ENDDO !L_NSPECTV
748          QREFvis3d(ig,lg,iaer) =                                  &
749                        k1*qrefVISgrid(grid_i,1,iaer) +            &
750                        k2*qrefVISgrid(grid_i+1,1,iaer)
751          omegaREFvis3d(ig,lg,iaer) =                              &
752                        k1*omegrefVISgrid(grid_i,1,iaer) +         &
753                        k2*omegrefVISgrid(grid_i+1,1,iaer)
754          ELSE                   ! INFRARED -----------------------
755          DO m=1,L_NSPECTI
756            QIRsQREF3d(ig,lg,m,iaer) =                             &
757                        k1*qsqrefIRgrid(grid_i,1,m,iaer) +         &
758                        k2*qsqrefIRgrid(grid_i+1,1,m,iaer)
759            omegaIR3d(ig,lg,m,iaer) =                              &
760                        k1*omegIRgrid(grid_i,1,m,iaer) +           &
761                        k2*omegIRgrid(grid_i+1,1,m,iaer)
762            gIR3d(ig,lg,m,iaer) =                                  &
763                        k1*gIRgrid(grid_i,1,m,iaer) +              &
764                        k2*gIRgrid(grid_i+1,1,m,iaer)
765          ENDDO !L_NSPECTI
766          QREFir3d(ig,lg,iaer) =                                   &
767                        k1*qrefIRgrid(grid_i,1,iaer) +             &
768                        k2*qrefIRgrid(grid_i+1,1,iaer)
769          omegaREFir3d(ig,lg,iaer) =                               &
770                        k1*omegrefIRgrid(grid_i,1,iaer) +          &
771                        k2*omegrefIRgrid(grid_i+1,1,iaer)
772          ENDIF                  ! --------------------------------
773        ENDDO !nlayer
774      ENDDO !ngrid
775
776!==================================================================
777
778
779
780      ENDDO ! idomain
781
782      ENDIF ! nsize = 1
783
784      ENDDO ! iaer (loop on aerosol kind)
785
786      RETURN
787      END SUBROUTINE aeroptproperties
788
789
790
791
Note: See TracBrowser for help on using the repository browser.