Update version info for release v4.6.1 (#2122)
[WRF.git] / phys / module_ra_aerosol.F
blob13b769e7dc5fcff00ad97e24d6ec9bd1d0bee548
1 ! Aerosol optical properties from AOD at 550 nm, aerosol type and relative humidity
2 ! Author: Jose A. Ruiz-Arias
3 ! Init date: 11/2012
5 ! Changes:
6 !    12/2012, jararias: re-inception using aer_*_opt and aer_*_val options
7 ! 16/07/2013, jararias: Angstrom exponent limited to values greater than -0.3. Personal
8 !                       communication con Chris Gueymard
9 ! 16/07/2013, jararias: Added generation of mean values of angexp, ssa and asy between 0.4 um
10 !                       and 1 um to be promoted to the wrfout when parameterization is used (*_opt=3)
11 ! 12/2013, jararias: integration in the official WRFv3.6
13 module module_ra_aerosol
15 ! global variables
17 contains
20 !--------------------------------------------------------------
21 !       INDICES CONVENTION
22 !--------------------------------------------------------------       
23 !       kms:kme define the range for full-level indices
24 !       kts:kte define the range for half-level indices
25 !       
26 !       kms=1 is the first full level at surface
27 !       kts=1 is the first half level at surface
28 !       
29 !       kme is the last full level at toa
30 !       kte is the last half level at toa
31 !       
32 !       There is one more full level than half levels.
33 !       Therefore, kme=kte+1. I checked it in one of my
34 !       simulations:
35 !       
36 !         namelist.input:
37 !              s_vert=1  e_vert=28        
38 !         code:
39 !              kms= 1    kts= 1
40 !              kms=28    kte=27
41 !       
42 !       In the vertical dimension there is no tiling for
43 !       parallelization as in the horizontal dimensions.
44 !       For i-dim and j-dim, the t-indices define the
45 !       range of indices over which each tile runs.
46 !--------------------------------------------------------------
48 !  namelist options:
49 !    aer_aod550_opt = [1,2] : 
50 !        1 = input constant value for AOD at 550 nm from namelist. 
51 !            In this case, the value is read from aer_aod550_val; 
52 !        2 = input value from auxiliary input 15. It is a time-varying 2D grid in netcdf wrf-compatible 
53 !            format. The default operation is aer_aod550_opt=1 and aer_aod550_val=0.12
54 !    aer_angexp_opt = [1,2,3] : 
55 !        1 = input constant value for Angstrom exponent from namelist. In this case, the value is read 
56 !            from aer_angexp_val; 
57 !        2 = input value from auxiliary input 15, as in aer_aod550_opt; 
58 !        3 = Angstrom exponent value estimated from the aerosol type defined in aer_type, and modulated 
59 !            with the RH in WRF. Default operation is aer_angexp_opt = 1, and aer_angexp_val=1.3.
60 !    aer_ssa_opt and aer_asy_opt are similar to aer_angexp_opt.
62 !    aer_type = [1,2,3] : 1 for rural, 2 is urban and 3 is maritime.
63 !--------------------------------------------------------------
65 subroutine calc_aerosol_goddard_sw(ht,dz8w,p,t3d,qv3d,aer_type,                               &
66                                    aer_aod550_opt, aer_angexp_opt, aer_ssa_opt, aer_asy_opt,  &
67                                    aer_aod550_val, aer_angexp_val, aer_ssa_val, aer_asy_val,  &
68                                    aod5502d, angexp2d, aerssa2d, aerasy2d,                    &
69                                    ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte,           &
70                                    tauaer, ssaaer, asyaer                                     )
72     USE module_wrf_error , ONLY : wrf_err_message
73     
74     implicit none
76     ! constants
77     integer, parameter :: N_BANDS=11
79     ! mid-band wavelength / 550. nm
80     real :: lower_wvl(N_BANDS),upper_wvl(N_BANDS)
81     integer :: kk
82     data (lower_wvl(kk),kk=1,N_BANDS) /0.175,0.225,0.245,0.280,0.295,0.310,0.325,0.400,0.700,1.220,2.270/
83     data (upper_wvl(kk),kk=1,N_BANDS) /0.225,0.245,0.280,0.295,0.310,0.325,0.400,0.700,1.220,2.270,10.00/
85     ! I/O variables
86     integer, intent(in) :: ims,ime,jms,jme,kms,kme, &
87                            its,ite,jts,jte,kts,kte
88     real, dimension(ims:ime, kms:kme, jms:jme), intent(in) :: p,     & ! pressure (Pa)
89                                                               t3d,   & ! temperature (K)
90                                                               dz8w,  & ! dz between full levels (m)
91                                                               qv3d     ! water vapor mixing ratio (kg/kg)
92     integer, intent(in) :: aer_type
93     integer, intent(in) :: aer_aod550_opt, aer_angexp_opt, aer_ssa_opt, aer_asy_opt
94     real, intent(in)    :: aer_aod550_val, aer_angexp_val, aer_ssa_val, aer_asy_val
96     real, dimension(ims:ime, jms:jme),                     intent(in)    :: ht
97     real, dimension(ims:ime, jms:jme),           optional, intent(inout) :: aod5502d, angexp2d, aerssa2d, aerasy2d
98     real, dimension(ims:ime, kms:kme, jms:jme, 1:N_BANDS), intent(inout) :: tauaer, ssaaer, asyaer
100     ! local variables
101     integer :: i,j,k,nb
102     real :: aod_rate,angexp_val,x,xy,xx
103     real, dimension(ims:ime, jms:jme, 1:N_BANDS) :: aod550spc
104     real, dimension(ims:ime, kms:kme, jms:jme)   :: rh  ! relative humidity
106     call calc_relative_humidity(p,t3d,qv3d,                 &
107                                 ims,ime,jms,jme,kms,kme,    &
108                                 its,ite,jts,jte,kts,kte,    &
109                                 rh                          )
111     aer_aod550_opt_select: select case(aer_aod550_opt)
112        !case(0)
113        ! reserved for climatology
114        case(1)
115           if (aer_aod550_val .lt. 0 ) then
116              write(wrf_err_message,'("aer_aod550_val must be positive. Negative value ",F7.4," found")') aer_aod550_val
117              call wrf_error_fatal(wrf_err_message)
118           end if
119           write( wrf_err_message, '("aer_aod550_opt=",I1,": AOD@550 nm fixed to value ",F6.3)') aer_aod550_opt,aer_aod550_val
120           call wrf_debug(100, wrf_err_message )
121           do j=jts,jte
122              do i=its,ite
123                 aod5502d(i,j)=aer_aod550_val
124              end do
125           end do
127        case(2) 
128           if (.not.(present(aod5502d))) then
129              write(wrf_err_message,*) 'Expected gridded total AOD@550 nm, but it is not in the radiation driver'
130              call wrf_error_fatal(wrf_err_message)
131           end if
132           if (minval(aod5502d) .lt. 0) then
133              call wrf_error_fatal('AOD@550 must be positive. Negative value(s) found in auxinput')
134           end if
135           write( wrf_err_message, '("aer_aod550_opt=",I1,": AOD@550 nm read from auxinput (min=",F6.3," max=",F6.3,")")') &
136                                   aer_aod550_opt,minval(aod5502d),maxval(aod5502d)
137           call wrf_debug(100, wrf_err_message )
139        case default
140           write(wrf_err_message,*) 'Expected aer_aod550_opt=[1,2]. Got',aer_aod550_opt
141           call wrf_error_fatal(wrf_err_message)
142     end select aer_aod550_opt_select
145     ! here, the 3d aod550 is calculated according to the aer_angexp_opt case
146     aer_angexp_opt_select: select case(aer_angexp_opt)
147        !case(0)
148        ! reserved for climatology
149        case(1)
150           if (aer_angexp_val .lt. -0.3) then
151              write(wrf_err_message,'("WARNING: aer_angexp_val limited to -0.3. Illegal value ",F7.4," found")') aer_angexp_val
152              call wrf_debug(100,wrf_err_message)
153           end if
154           if (aer_angexp_val .gt. 2.5) then
155              write(wrf_err_message,'("WARNING: aer_angexp_val limited to 2.5. Illegal value ",F7.4," found")') aer_angexp_val
156              call wrf_debug(100,wrf_err_message)
157           end if
158           write( wrf_err_message , '("aer_angexp_opt=",I1,": Aerosol Angstrom exponent fixed to value ",F6.3)') &
159                                       aer_angexp_opt,aer_angexp_val
160           call wrf_debug(100, wrf_err_message )
161           angexp_val=min(2.5,max(-0.3,aer_angexp_val))
162           do nb=1,N_BANDS
163              if ((angexp_val .lt. 0.999) .or. (angexp_val .gt. 1.001)) then
164                 aod_rate=((0.55**angexp_val)*(upper_wvl(nb)**(1.-angexp_val)- &
165                           lower_wvl(nb)**(1.-angexp_val)))/((1.-angexp_val)*(upper_wvl(nb)-lower_wvl(nb)))
166              else
167                 aod_rate=(0.55/(upper_wvl(nb)-lower_wvl(nb)))*log(upper_wvl(nb)/lower_wvl(nb))
168              end if
169              do j=jts,jte
170                 do i=its,ite
171                    aod550spc(i,j,nb)=aod5502d(i,j)*aod_rate
172                 end do
173              end do
174           end do
175           do j=jts,jte
176              do i=its,ite
177                 angexp2d(i,j)=angexp_val
178              end do
179           end do
180        case(2)
181           if (.not.(present(angexp2d))) then
182              write(wrf_err_message,*) 'Expected gridded aerosol Angstrom exponent, but it is not in the radiation driver'
183              call wrf_error_fatal(wrf_err_message)
184           end if
185           write( wrf_err_message, '("aer_angexp_opt=",I1,": Angstrom exponent read from auxinput (min=",F6.3," max=",F6.3,")")') &
186                                   aer_angexp_opt,minval(angexp2d),maxval(angexp2d)
187           call wrf_debug(100, wrf_err_message )
188           do j=jts,jte
189              do i=its,ite
190                 angexp_val=min(2.5,max(-0.3,angexp2d(i,j)))
191                 do nb=1,N_BANDS
192                    if ((angexp_val .lt. 0.999) .or. (angexp_val .gt. 1.001)) then
193                       aod_rate=((0.55**angexp_val)*(upper_wvl(nb)**(1.-angexp_val)- &
194                                 lower_wvl(nb)**(1.-angexp_val)))/((1.-angexp_val)*(upper_wvl(nb)-lower_wvl(nb)))
195                    else
196                       aod_rate=(0.55/(upper_wvl(nb)-lower_wvl(nb)))*log(upper_wvl(nb)/lower_wvl(nb))
197                    end if
198                    aod550spc(i,j,nb)=aod5502d(i,j)*aod_rate
199                 end do
200              end do
201           end do
203        case(3)
204           ! spectral disaggregation based on prescribed aerosol type and relative humidity
205           write( wrf_err_message, '("aer_angexp_opt=",I1,": Angstrom exponent calculated from RH and aer_type ",I1)') &
206                                      aer_angexp_opt,aer_type
207           call wrf_debug(100, wrf_err_message )
208           call calc_spectral_aod_goddard_sw(ims,ime,jms,jme,kms,kme,  &
209                                             its,ite,jts,jte,kts,kte,  &
210                                             rh,aer_type,aod5502d,     &
211                                             aod550spc                 )
212           ! added July, 16th, 2013: angexp2d is in the wrfout when aer_angexp_opt=3. It is the average
213           ! value in the spectral bands between 0.4 and 1. um
214           do j=jts,jte
215              do i=its,ite
216                 xy=0
217                 xx=0
218                 do nb=8,9  ! bands between 0.4 and  1.0 um
219                    ! the slope of a linear regression with intercept=0 is m=E(xy)/E(x^2), where y=m*x
220                    x=log(0.5*(lower_wvl(nb)+upper_wvl(nb))/0.55)
221                    xy=xy+x*log(aod550spc(i,j,nb)/aod5502d(i,j))
222                    xx=xx+x*x
223                 end do
224                 angexp2d(i,j)=-xy/xx
225              end do
226           end do
228        case default
229           write(wrf_err_message,*) 'Expected aer_angexp_opt=[1,2,3]. Got',aer_angexp_opt
230           call wrf_error_fatal(wrf_err_message)
231     end select aer_angexp_opt_select
233     ! exponental -vertical- profile
234     call aod_profiler(ht,dz8w,aod550spc,n_bands,ims,ime,jms,jme,kms,kme, &
235                       its,ite,jts,jte,kts,kte,tauaer                     )
237     aer_ssa_opt_select: select case(aer_ssa_opt)
238        !case(0)
239        ! reserved for climatology
240        case(1)
241           if ((aer_ssa_val .lt. 0) .or. (aer_ssa_val .gt. 1)) then
242              write(wrf_err_message,'("aer_ssa_val must be within [0,1]. Illegal value ",F7.4," found")') aer_ssa_val
243              call wrf_error_fatal(wrf_err_message)
244           end if
245           write( wrf_err_message , '("aer_ssa_opt=",I1,": single-scattering albedo fixed to value ",F6.3)') &
246                                       aer_ssa_opt,aer_ssa_val
247           call wrf_debug(100, wrf_err_message )
248           do j=jts,jte
249              do i=its,ite
250                 do k=kts,kte
251                    do nb=1,N_BANDS
252                       ! no spectral disaggregation
253                       ssaaer(i,k,j,nb)=aer_ssa_val
254                    end do
255                 end do
256              end do
257           end do
259           do j = jts, jte
260              do i = its, ite
261                 aerssa2d(i, j) = aer_ssa_val
262              end do
263           end do
265        case(2)
266           if (.not.(present(aerssa2d))) then
267              write(wrf_err_message,*) &
268                   'Expected gridded aerosol single-scattering albedo, but it is not in the radiation driver'
269              call wrf_error_fatal(wrf_err_message)
270           end if
271           if ((minval(aerssa2d) .lt. 0) .or. (maxval(aerssa2d) .gt. 1)) then
272              call wrf_error_fatal('Aerosol single-scattering albedo must be within [0,1]. ' // &
273                                   'Out of bounds value(s) found in auxinput')
274           end if
275           write( wrf_err_message, &
276               '("aer_ssa_opt=",I1,": single-scattering albedo read from auxinput (min=",F6.3," max=",F6.3,")")') &
277                                   aer_ssa_opt,minval(aerssa2d),maxval(aerssa2d)
278           call wrf_debug(100, wrf_err_message )
279           do j=jts,jte
280              do i=its,ite
281                 do k=kts,kte
282                    do nb=1,N_BANDS
283                       ! no spectral disaggregation
284                       ssaaer(i,k,j,nb)=aerssa2d(i,j)
285                    end do
286                 end do
287              end do
288           end do
290        case(3)
291           ! spectral disaggregation based on prescribed aerosol type and relative humidity
292           write( wrf_err_message, '("aer_ssa_opt=",I1,": single-scattering albedo calculated from RH and aer_type ",I1)') &
293                                      aer_ssa_opt,aer_type
294           call wrf_debug(100, wrf_err_message )
295           call calc_spectral_ssa_goddard_sw(ims,ime,jms,jme,kms,kme,  &
296                                             its,ite,jts,jte,kts,kte,  &
297                                             rh,aer_type,ssaaer        )
298           ! added July, 16th, 2013: aerssa2d is in the wrfout when aer_ssa_opt=3. It is the average
299           ! value in the spectral bands between 0.4 and 1. um
300           do j=jts,jte
301              do i=its,ite
302                 aerssa2d(i,j)=0
303              end do
304           end do
305           do j=jts,jte
306              do i=its,ite
307                 do nb=8,9  ! bands between 0.4 and 1.0 um
308                    aerssa2d(i,j)=aerssa2d(i,j)+ssaaer(i,kts,j,nb)
309                 end do
310                 aerssa2d(i,j)=aerssa2d(i,j)/2.
311              end do
312           end do
314        case default
315           write(wrf_err_message,*) 'Expected aer_ssa_opt=[1,2,3]. Got',aer_ssa_opt
316           call wrf_error_fatal(wrf_err_message)
317     end select aer_ssa_opt_select
319     aer_asy_opt_select: select case(aer_asy_opt)
320        !case(0)
321        ! reserved for climatology
322        case(1)
323           if ((aer_asy_val .lt. -1) .or. (aer_asy_val .gt. 1)) then
324              write(wrf_err_message,'("aer_asy_val must be withing [-1,1]. Illegal value ",F7.4," found")') aer_asy_val
325              call wrf_error_fatal(wrf_err_message)
326           end if
327           write( wrf_err_message , '("aer_asy_opt=",I1,": asymmetry parameter fixed to value ",F6.3)') aer_asy_opt,aer_asy_val
328           call wrf_debug(100, wrf_err_message )
329           do j=jts,jte
330              do i=its,ite
331                 do k=kts,kte
332                    do nb=1,N_BANDS
333                       asyaer(i,k,j,nb)=aer_asy_val
334                    end do
335                 end do
336              end do
337           end do
339           do j = jts, jte
340              do i = its, ite
341                 aerasy2d(i, j) = aer_asy_val
342              end do
343           end do
345        case(2)
346           if (.not.(present(aerasy2d))) then
347              write(wrf_err_message,*) 'Expected gridded aerosol asymmetry parameter, but it is not in the radiation driver'
348              call wrf_error_fatal(wrf_err_message)
349           end if
350           if ((minval(aerasy2d) .lt. -1) .or. (maxval(aerasy2d) .gt. 1)) then
351              call wrf_error_fatal('Aerosol asymmetry parameter must be within [-1,1]. Out of bounds value(s) found in auxinput')
352           end if
353           write( wrf_err_message, '("aer_asy_opt=",I1,": asymmetry parameter read from auxinput (min=",F6.3," max=",F6.3,")")') &
354                                   aer_asy_opt,minval(aerasy2d),maxval(aerasy2d)
355           call wrf_debug(100, wrf_err_message )
356           do j=jts,jte
357              do i=its,ite
358                 do k=kts,kte
359                    do nb=1,N_BANDS
360                       asyaer(i,k,j,nb)=aerasy2d(i,j)
361                    end do
362                 end do
363              end do
364           end do
366        case(3)
367           ! spectral disaggregation based on a prescribed aerosol type and relative humidity
368           write( wrf_err_message, '("aer_asy_opt=",I1,": asymmetry parameter calculated from RH and aer_type ",I1)') &
369                                      aer_asy_opt,aer_type
370           call wrf_debug(100, wrf_err_message )
371           call calc_spectral_asy_goddard_sw(ims,ime,jms,jme,kms,kme,  &
372                                             its,ite,jts,jte,kts,kte,  &
373                                             rh,aer_type,asyaer        )
374           ! added July, 16th, 2013: aerasy2d is in the wrfout when aer_asy_opt=3. It is the average
375           ! value in the spectral bands between 0.4 and 1. um
376           do j=jts,jte
377              do i=its,ite
378                 aerasy2d(i,j)=0
379              end do
380           end do
381           do j=jts,jte
382              do i=its,ite
383                 do nb=8,9  ! bands between 0.4 and 1.0 um
384                    aerasy2d(i,j)=aerasy2d(i,j)+asyaer(i,kts,j,nb)
385                 end do
386                 aerasy2d(i,j)=aerasy2d(i,j)/2.
387              end do
388           end do
390        case default
391           write(wrf_err_message,*) 'Expected aer_asy_opt=[1,2,3]. Got',aer_asy_opt
392           call wrf_error_fatal(wrf_err_message)
393     end select aer_asy_opt_select
394 end subroutine calc_aerosol_goddard_sw
397 subroutine calc_aerosol_rrtmg_sw(ht,dz8w,p,t3d,qv3d,aer_type,                               &
398                                  aer_aod550_opt, aer_angexp_opt, aer_ssa_opt, aer_asy_opt,  &
399                                  aer_aod550_val, aer_angexp_val, aer_ssa_val, aer_asy_val,  &
400                                  aod5502d, angexp2d, aerssa2d, aerasy2d,                    &
401                                  ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte,           &
402                                  tauaer, ssaaer, asyaer, aod5503d                           )
404     USE module_wrf_error , ONLY : wrf_err_message
405     
406     implicit none
408     ! constants
409     integer, parameter :: N_BANDS=14
410     ! local index variables
411     integer :: i,j,k,nb
413     real :: lower_wvl(N_BANDS),upper_wvl(N_BANDS)
414     data (lower_wvl(i),i=1,N_BANDS) /3.077,2.500,2.150,1.942,1.626,1.299,1.242,0.7782,0.6250,0.4415,0.3448,0.2632,0.2000,3.846/
415     data (upper_wvl(i),i=1,N_BANDS) /3.846,3.077,2.500,2.150,1.942,1.626,1.299,1.2420,0.7782,0.6250,0.4415,0.3448,0.2632,12.195/
417     ! I/O variables
418     real, dimension(ims:ime, kms:kme, jms:jme), intent(in) :: p,     & ! pressure (Pa)
419                                                               t3d,   & ! temperature (K)
420                                                               dz8w,  & ! dz between full levels (m)
421                                                               qv3d     ! water vapor mixing ratio (kg/kg)
422     integer, intent(in) :: ims,ime,jms,jme,kms,kme, &
423                            its,ite,jts,jte,kts,kte
424     integer, intent(in) :: aer_type
425     integer, intent(in) :: aer_aod550_opt, aer_angexp_opt, aer_ssa_opt, aer_asy_opt
426     real, intent(in)    :: aer_aod550_val, aer_angexp_val, aer_ssa_val, aer_asy_val
428     real, dimension(ims:ime, jms:jme), intent(in)    :: ht
429     real, dimension(ims:ime, jms:jme), optional, intent(inout) :: aod5502d, angexp2d, aerssa2d, aerasy2d
430     real, dimension(ims:ime, kms:kme, jms:jme, 1:N_BANDS), intent(inout) :: tauaer, ssaaer, asyaer
432     real, dimension(ims:ime, kms:kme, jms:jme), optional, intent(in) :: aod5503d   ! trude
434     ! local variables
435     real :: angexp_val,aod_rate,x,xy,xx
436     real, dimension(ims:ime, jms:jme, 1:N_BANDS) :: aod550spc
437     real, dimension(ims:ime, kms:kme, jms:jme, 1:N_BANDS) :: aod550spc3d    !  trude
438     real, dimension(ims:ime, kms:kme, jms:jme)   :: rh  ! relative humidity
440     call calc_relative_humidity(p,t3d,qv3d,                 &
441                                 ims,ime,jms,jme,kms,kme,    &
442                                 its,ite,jts,jte,kts,kte,    &
443                                 rh                          )
445     aer_aod550_opt_select: select case(aer_aod550_opt)
446        !case(0)
447        ! reserved for climatology
448        case(1)
449           if (aer_aod550_val .lt. 0) then
450              write(wrf_err_message,'("aer_aod550_val must be positive. Negative value ",F7.4," found")') aer_aod550_val
451              call wrf_error_fatal(wrf_err_message)
452           end if
453           write( wrf_err_message, '("aer_aod550_opt=",I1,": AOD@550 nm fixed to value ",F6.3)') aer_aod550_opt,aer_aod550_val
454           call wrf_debug(100, wrf_err_message )
455           do j=jts,jte
456              do i=its,ite
457                 aod5502d(i,j)=aer_aod550_val
458              end do
459           end do
461        case(2) 
462           if (.not.(present(aod5502d))) then
463              write(wrf_err_message,*) 'Expected gridded total AOD@550 nm, but it is not in the radiation driver'
464              call wrf_error_fatal(wrf_err_message)
465           end if
466           if (minval(aod5502d) .lt. 0) then
467              call wrf_error_fatal('AOD@550 must be positive. Negative value(s) found in auxinput')
468           end if
469           write( wrf_err_message, '("aer_aod550_opt=",I1,": AOD@550 nm read from auxinput (min=",F6.3," max=",F6.3,")")') &
470                                   aer_aod550_opt,minval(aod5502d),maxval(aod5502d)
471           call wrf_debug(100, wrf_err_message )
473        case default
474           write(wrf_err_message,*) 'Expected aer_aod550_opt=[1,2]. Got',aer_aod550_opt
475           call wrf_error_fatal(wrf_err_message)
476     end select aer_aod550_opt_select
479     ! here, the 3d aod550 is calculated according to the aer_angexp_opt case
480     aer_angexp_opt_select: select case(aer_angexp_opt)
481        !case(0)
482        ! reserved for climatology
483        case(1)
484           if (aer_angexp_val .lt. -0.3) then
485              write(wrf_err_message,'("WARNING: aer_angexp_val limited to -0.3. Illegal value ",F7.4," found")') aer_angexp_val
486              call wrf_debug(100,wrf_err_message)
487           end if
488           if (aer_angexp_val .gt. 2.5) then
489              write(wrf_err_message,'("WARNING: aer_angexp_val limited to 2.5. Illegal value ",F7.4," found")') aer_angexp_val
490              call wrf_debug(100,wrf_err_message)
491           end if
492           write( wrf_err_message , '("aer_angexp_opt=",I1,": Aerosol Angstrom exponent fixed to value ",F6.3)') &
493                                       aer_angexp_opt,aer_angexp_val
494           call wrf_debug(100, wrf_err_message )
495           angexp_val=min(2.5,max(-0.3,aer_angexp_val))
496           do nb=1,N_BANDS
497              if ((angexp_val .lt. 0.999) .or. (angexp_val .gt. 1.001)) then
498                 aod_rate=((0.55**angexp_val)*(upper_wvl(nb)**(1.-angexp_val)- &
499                           lower_wvl(nb)**(1.-angexp_val)))/((1.-angexp_val)*(upper_wvl(nb)-lower_wvl(nb)))
500              else
501                 aod_rate=(0.55/(upper_wvl(nb)-lower_wvl(nb)))*log(upper_wvl(nb)/lower_wvl(nb))
502              end if
503              do j=jts,jte
504                 do i=its,ite
505                    aod550spc(i,j,nb)=aod5502d(i,j)*aod_rate
506                 end do
507              end do
508           end do
509           do j=jts,jte
510              do i=its,ite
511                 angexp2d(i,j)=angexp_val
512              end do
513           end do
514        case(2)
515           if (.not.(present(angexp2d))) then
516              write(wrf_err_message,*) 'Expected gridded aerosol Angstrom exponent, but it is not in the radiation driver'
517              call wrf_error_fatal(wrf_err_message)
518           end if
519           write( wrf_err_message, '("aer_angexp_opt=",I1,": Angstrom exponent read from auxinput (min=",F6.3," max=",F6.3,")")') &
520                                   aer_angexp_opt,minval(angexp2d),maxval(angexp2d)
521           call wrf_debug(100, wrf_err_message )
522           do j=jts,jte
523              do i=its,ite
524                 angexp_val=min(2.5,max(-0.3,angexp2d(i,j)))
525                 do nb=1,N_BANDS
526                    if ((angexp_val .lt. 0.999) .or. (angexp_val .gt. 1.001)) then
527                       aod_rate=((0.55**angexp_val)*(upper_wvl(nb)**(1.-angexp_val)- &
528                                 lower_wvl(nb)**(1.-angexp_val)))/((1.-angexp_val)*(upper_wvl(nb)-lower_wvl(nb)))
529                    else
530                       aod_rate=(0.55/(upper_wvl(nb)-lower_wvl(nb)))*log(upper_wvl(nb)/lower_wvl(nb))
531                    end if
532                    aod550spc(i,j,nb)=aod5502d(i,j)*aod_rate
533                 end do
534              end do
535           end do
537        case(3)
538           ! spectral disaggregation based on a prescribed aerosol type and relative humidity
539           write( wrf_err_message, '("aer_angexp_opt=",I1,": Angstrom exponent calculated from RH and aer_type ",I1)') &
540                                      aer_angexp_opt,aer_type
541           call wrf_debug(100, wrf_err_message )
542           call calc_spectral_aod_rrtmg_sw(ims,ime,jms,jme,kms,kme,      &
543                                           its,ite,jts,jte,kts,kte,      &
544                                           rh,aer_type,aod5502d,         &
545                                           aod550spc,                    &
546                                           aod5503d, aod550spc3d)         ! trude
548           do j=jts,jte
549             do i=its,ite
550               angexp2d(i,j) = 0.0
551             enddo
552           enddo
554           if (present(aod5503d)) then
555             do j=jts,jte
556              do k=kts,kte
557               do i=its,ite
558                 xy=0
559                 xx=0
560                 do nb=8,N_BANDS-3  ! bands between 0.4 and  1.0 um
561                    ! the slope of a linear regression with intercept=0 is m=E(xy)/E(x^2), where y=m*x
562                    x=log(0.5*(lower_wvl(nb)+upper_wvl(nb))/0.55)
563                    xy=xy+x*log(aod550spc3d(i,k,j,nb)/aod5503d(i,k,j))
564                    xx=xx+x*x
565                 end do
566                 angexp2d(i,j) = angexp2d(i,j) - (xy/(N_BANDS-3-8+1))/(xx/(N_BANDS-3-8+1))*(p(i,k,j)-p(i,k+1,j))/(p(i,kts,j)-p(i,kte,j))
567               enddo
568              enddo
569             enddo
570           else
572           ! added July, 16th, 2013: angexp2d is in the wrfout when aer_angexp_opt=3. It is the average
573           ! value in the spectral bands between 0.4 and 1. um
574           do j=jts,jte
575              do i=its,ite
576                 xy=0
577                 xx=0
578                 do nb=8,N_BANDS-3  ! bands between 0.4 and  1.0 um
579                    ! the slope of a linear regression with intercept=0 is m=E(xy)/E(x^2), where y=m*x
580                    x=log(0.5*(lower_wvl(nb)+upper_wvl(nb))/0.55)
581                    xy=xy+x*log(aod550spc(i,j,nb)/aod5502d(i,j))
582                    xx=xx+x*x
583                 end do
584                 angexp2d(i,j)=-(xy/(N_BANDS-3-8+1))/(xx/(N_BANDS-3-8+1))
585              end do
586           end do
587           endif
589        case default
590           write(wrf_err_message,*) 'Expected aer_angexp_opt=[1,2,3]. Got',aer_angexp_opt
591           call wrf_error_fatal(wrf_err_message)
592     end select aer_angexp_opt_select
594 !..If 3D AOD (at 550nm) was provided explicitly, then no need to assume a
595 !.. vertical distribution, just use what was provided.  (Trude)
597       if (present(aod5503d)) then
598          do nb=1,N_BANDS
599           do j=jts,jte
600            do k=kts,kte
601             do i=its,ite
602                tauaer(i,k,j,nb) = aod550spc3d(i,k,j,nb)
603             enddo
604            enddo
605           enddo
606          enddo
607       else
608          ! exponental -vertical- profile
609          call aod_profiler(ht,dz8w,aod550spc,n_bands,ims,ime,jms,jme,kms,kme, &
610                       its,ite,jts,jte,kts,kte,tauaer                     )
611       endif
613     aer_ssa_opt_select: select case(aer_ssa_opt)
614        !case(0)
615        ! reserved for climatology
616        case(1)
617           if ((aer_ssa_val .lt. 0) .or. (aer_ssa_val .gt. 1)) then
618              write(wrf_err_message,'("aer_ssa_val must be within [0,1]. Illegal value ",F7.4," found")') aer_ssa_val
619              call wrf_error_fatal(wrf_err_message)
620           end if
621           write( wrf_err_message, &
622                 '("aer_ssa_opt=",I1,": single-scattering albedo fixed to value ",F6.3)') aer_ssa_opt,aer_ssa_val
623           call wrf_debug(100, wrf_err_message )
624           do j=jts,jte
625              do i=its,ite
626                 do k=kts,kte
627                    do nb=1,N_BANDS
628                       ! no spectral disaggregation
629                       ssaaer(i,k,j,nb)=aer_ssa_val
630                    end do
631                 end do
632              end do
633           end do
634           do j=jts,jte
635              do i=its,ite
636                 aerssa2d(i,j)=aer_ssa_val
637              end do
638           end do
640        case(2)
641           if (.not.(present(aerssa2d))) then
642              write(wrf_err_message,*) 'Expected gridded aerosol single-scattering albedo, but it is not in the radiation driver'
643              call wrf_error_fatal(wrf_err_message)
644           end if
645           if ((minval(aerssa2d) .lt. 0) .or. (maxval(aerssa2d) .gt. 1)) then
646              call wrf_error_fatal('Aerosol single-scattering albedo must be within [0,1]. ' // &
647                                   'Out of bounds value(s) found in auxinput')
648           end if
649           write( wrf_err_message, '("aer_ssa_opt=",I1,": single-scattering albedo from auxinput (min=",F6.3," max=",F6.3,")")') &
650                                   aer_ssa_opt,minval(aerssa2d),maxval(aerssa2d)
651           call wrf_debug(100, wrf_err_message )
652           do j=jts,jte
653              do i=its,ite
654                 do k=kts,kte
655                    do nb=1,N_BANDS
656                       ! no spectral disaggregation
657                       ssaaer(i,k,j,nb)=aerssa2d(i,j)
658                    end do
659                 end do
660              end do
661           end do
663        case(3)
664           ! spectral disaggregation based on a prescribed aerosol type and relative humidity
665           write( wrf_err_message, '("aer_ssa_opt=",I1,": single-scattering albedo calculated from RH and aer_type ",I1)') &
666                                      aer_ssa_opt,aer_type
667           call wrf_debug(100, wrf_err_message )
668           call calc_spectral_ssa_rrtmg_sw(ims,ime,jms,jme,kms,kme,  &
669                                           its,ite,jts,jte,kts,kte,  &
670                                           rh,aer_type,ssaaer        )
671           ! added July, 16th, 2013: aerssa2d is in the wrfout when aer_ssa_opt=3. It is the average
672           ! value in the spectral bands between 0.4 and 1. um
673           do j=jts,jte
674              do i=its,ite
675                 aerssa2d(i,j)=0
676              end do
677           end do
678           do j=jts,jte
679              do i=its,ite
680                 do nb=8,N_BANDS-3  ! bands between 0.4 and 1.0 um
681                    aerssa2d(i,j)=aerssa2d(i,j)+ssaaer(i,kts,j,nb)
682                 end do
683                 aerssa2d(i,j)=aerssa2d(i,j)/(N_BANDS-3-8+1)
684              end do
685           end do
687        case default
688           write(wrf_err_message,*) 'Expected aer_ssa_opt=[1,2,3]. Got',aer_ssa_opt
689           call wrf_error_fatal(wrf_err_message)
690     end select aer_ssa_opt_select
692     aer_asy_opt_select: select case(aer_asy_opt)
693        !case(0)
694        ! reserved for climatology
695        case(1)
696           if ((aer_asy_val .lt. 0) .or. (aer_asy_val .gt. 1)) then
697              write(wrf_err_message,'("aer_asy_val must be withing [-1,1]. Illegal value ",F7.4," found")') aer_asy_val
698              call wrf_error_fatal(wrf_err_message)
699           end if
700           write( wrf_err_message , '("aer_asy_opt=",I1,": asymmetry parameter fixed to value ",F6.3)') aer_asy_opt,aer_asy_val
701           call wrf_debug(100, wrf_err_message )
702           do j=jts,jte
703              do i=its,ite
704                 do k=kts,kte
705                    do nb=1,N_BANDS
706                       asyaer(i,k,j,nb)=aer_asy_val
707                    end do
708                 end do
709              end do
710           end do
711           do j=jts,jte
712              do i=its,ite
713                 aerasy2d(i,j)=aer_asy_val
714              end do
715           end do
717        case(2)
718           if (.not.(present(aerasy2d))) then
719              write(wrf_err_message,*) 'Expected gridded aerosol asymmetry parameter, but it is not in the radiation driver'
720              call wrf_error_fatal(wrf_err_message)
721           end if
722           if ((minval(aerasy2d) .lt. -1) .or. (maxval(aerasy2d) .gt. 1)) then
723              call wrf_error_fatal('Aerosol asymmetry parameter must be within [-1,1]. Out of bounds value(s) found in auxinput')
724           end if
725           write( wrf_err_message, '("aer_asy_opt=",I1,": asymmetry parameter read from auxinput (min=",F6.3," max=",F6.3,")")') &
726                                   aer_asy_opt,minval(aerasy2d),maxval(aerasy2d)
727           call wrf_debug(100, wrf_err_message )
728           do j=jts,jte
729              do i=its,ite
730                 do k=kts,kte
731                    do nb=1,N_BANDS
732                       asyaer(i,k,j,nb)=aerasy2d(i,j)
733                    end do
734                 end do
735              end do
736           end do
738        case(3)
739           ! spectral disaggregation based on a prescribed aerosol type and relative humidity
740           write( wrf_err_message, '("aer_asy_opt=",I1,": asymmetry parameter calculated from RH and aer_type ",I1)') &
741                                      aer_asy_opt,aer_type
742           call wrf_debug(100, wrf_err_message )
743           call calc_spectral_asy_rrtmg_sw(ims,ime,jms,jme,kms,kme,  &
744                                           its,ite,jts,jte,kts,kte,  &
745                                           rh,aer_type,asyaer        )
746           ! added July, 16th, 2013: aerasy2d is in the wrfout when aer_asy_opt=3. It is the average
747           ! value in the spectral bands between 0.4 and 1. um
748           do j=jts,jte
749              do i=its,ite
750                 aerasy2d(i,j)=0
751              end do
752           end do
753           do j=jts,jte
754              do i=its,ite
755                 do nb=8,N_BANDS-3  ! bands between 0.4 and 1.0 um
756                    aerasy2d(i,j)=aerasy2d(i,j)+asyaer(i,kts,j,nb)
757                 end do
758                 aerasy2d(i,j)=aerasy2d(i,j)/(N_BANDS-3-8+1)
759              end do
760           end do
762        case default
763           write(wrf_err_message,*) 'Expected aer_asy_opt=[1,2,3]. Got',aer_asy_opt
764           call wrf_error_fatal(wrf_err_message)
765     end select aer_asy_opt_select
766 end subroutine calc_aerosol_rrtmg_sw
769 subroutine calc_spectral_aod_goddard_sw(ims,ime,jms,jme,kms,kme,  &
770                                         its,ite,jts,jte,kts,kte,  &
771                                         rh,aer_type,aod550,       &
772                                         tauaer                    )
773    implicit none
774    
775    ! constants
776    integer, parameter :: N_AER_TYPES=3
777    integer, parameter :: N_RH=8
778    integer, parameter :: N_BANDS=11
779    integer, parameter :: N_INT_POINTS=4
781    ! I/O variables
782    integer, intent(in) :: ims,ime,jms,jme,kms,kme, &
783                           its,ite,jts,jte,kts,kte
784    integer, intent(in) :: aer_type
785    real, dimension(ims:ime, kms:kme, jms:jme),   intent(in) :: rh        ! relative humidity
786    real, dimension(ims:ime, jms:jme),            intent(in) :: aod550    ! Total AOD at 550 nm at surface
787    real, dimension(ims:ime, jms:jme, 1:N_BANDS), intent(inout) :: tauaer ! Total spectral aerosol optical depth at surface
789    ! local variables
790    integer :: i,j,k,ib,imax,imin,ii,jj
791    real    :: rhs(N_RH),lj
792    real    :: raod_lut(N_AER_TYPES,N_BANDS,N_RH)
794    ! relative humidity steps
795    data (rhs(i),i=1,8) /0.,50.,70.,80.,90.,95.,98.,99./
797    ! aer_type = 1 : rural (SF79)
798    data (raod_lut(1,ib,1),ib=1,N_BANDS) /2.5849,2.2217,2.0768,1.8320,1.7528,1.6889,1.4802,1.0127,0.5051,0.2292,0.0924/
799    data (raod_lut(1,ib,2),ib=1,N_BANDS) /2.5645,2.2070,2.0642,1.8228,1.7446,1.6816,1.4753,1.0121,0.5061,0.2302,0.0930/
800    data (raod_lut(1,ib,3),ib=1,N_BANDS) /2.5285,2.1809,2.0419,1.8064,1.7301,1.6685,1.4667,1.0111,0.5080,0.2321,0.0943/
801    data (raod_lut(1,ib,4),ib=1,N_BANDS) /2.4861,2.1501,2.0155,1.7871,1.7129,1.6530,1.4565,1.0100,0.5104,0.2345,0.0959/
802    data (raod_lut(1,ib,5),ib=1,N_BANDS) /2.3824,2.0745,1.9505,1.7392,1.6703,1.6146,1.4310,1.0072,0.5172,0.2415,0.1005/
803    data (raod_lut(1,ib,6),ib=1,N_BANDS) /2.2463,1.9744,1.8642,1.6752,1.6132,1.5630,1.3965,1.0039,0.5293,0.2540,0.1091/
804    data (raod_lut(1,ib,7),ib=1,N_BANDS) /2.0619,1.8373,1.7452,1.5861,1.5336,1.4908,1.3479,1.0007,0.5559,0.2827,0.1298/
805    data (raod_lut(1,ib,8),ib=1,N_BANDS) /1.9550,1.7568,1.6751,1.5332,1.4861,1.4477,1.3185,1.0004,0.5825,0.3131,0.1532/
807    ! aer_type = 2 : urban (SF79)
808    data (raod_lut(2,ib,1),ib=1,N_BANDS) /2.3427,2.0454,1.9254,1.7206,1.6538,1.5997,1.4210,1.0163,0.5774,0.3072,0.1486/
809    data (raod_lut(2,ib,2),ib=1,N_BANDS) /2.3288,2.0352,1.9166,1.7141,1.6480,1.5944,1.4175,1.0139,0.5671,0.2953,0.1393/
810    data (raod_lut(2,ib,3),ib=1,N_BANDS) /2.3033,2.0164,1.9004,1.7021,1.6373,1.5848,1.4111,1.0115,0.5588,0.2859,0.1322/
811    data (raod_lut(2,ib,4),ib=1,N_BANDS) /2.2717,1.9931,1.8803,1.6872,1.6239,1.5727,1.4030,1.0095,0.5546,0.2814,0.1288/
812    data (raod_lut(2,ib,5),ib=1,N_BANDS) /2.1862,1.9299,1.8256,1.6464,1.5875,1.5398,1.3809,1.0055,0.5523,0.2788,0.1269/
813    data (raod_lut(2,ib,6),ib=1,N_BANDS) /2.0524,1.8301,1.7390,1.5814,1.5294,1.4870,1.3453,1.0001,0.5550,0.2818,0.1291/
814    data (raod_lut(2,ib,7),ib=1,N_BANDS) /1.8164,1.6516,1.5830,1.4630,1.4229,1.3901,1.2790,0.9909,0.5657,0.2937,0.1381/
815    data (raod_lut(2,ib,8),ib=1,N_BANDS) /1.6384,1.5144,1.4622,1.3699,1.3388,1.3132,1.2257,0.9839,0.5778,0.3077,0.1489/
817    ! aer_type = 3 : maritime (SF79)
818    data (raod_lut(3,ib,1),ib=1,N_BANDS) /1.6325,1.5098,1.4581,1.3667,1.3359,1.3106,1.2238,1.0070,0.7288,0.5088,0.3361/
819    data (raod_lut(3,ib,2),ib=1,N_BANDS) /1.5338,1.4327,1.3898,1.3135,1.2876,1.2662,1.1927,1.0058,0.7593,0.5556,0.3875/
820    data (raod_lut(3,ib,3),ib=1,N_BANDS) /1.4227,1.3449,1.3117,1.2520,1.2316,1.2148,1.1563,1.0049,0.8006,0.6225,0.4656/
821    data (raod_lut(3,ib,4),ib=1,N_BANDS) /1.3439,1.2821,1.2554,1.2073,1.1908,1.1772,1.1295,1.0046,0.8359,0.6827,0.5405/
822    data (raod_lut(3,ib,5),ib=1,N_BANDS) /1.2451,1.2023,1.1837,1.1499,1.1383,1.1286,1.0944,1.0050,0.8897,0.7800,0.6701/
823    data (raod_lut(3,ib,6),ib=1,N_BANDS) /1.1867,1.1548,1.1408,1.1153,1.1064,1.0991,1.0730,1.0056,0.9280,0.8533,0.7745/
824    data (raod_lut(3,ib,7),ib=1,N_BANDS) /1.1485,1.1234,1.1124,1.0923,1.0852,1.0794,1.0586,1.0062,0.9565,0.9099,0.8589/
825    data (raod_lut(3,ib,8),ib=1,N_BANDS) /1.1352,1.1124,1.1025,1.0842,1.0778,1.0725,1.0536,1.0065,0.9671,0.9315,0.8920/
827    ! this spectral disaggregation is only at surface (k=kts)
828    do j=jts,jte
829       do i=its,ite
830          ! common part of the Lagrange's interpolator
831          !   only depends on the relative humidity value
832          ii=1
833          do while ( (ii.le.N_RH) .and. (rh(i,kts,j).gt.rhs(ii)) )
834             ii=ii+1
835          end do
836          imin=max(1,ii-N_INT_POINTS/2-1)
837          imax=min(N_RH,ii+N_INT_POINTS/2)
839          do ib=1,N_BANDS
840             tauaer(i,j,ib)=0.
841             do jj=imin,imax
842                lj=1.
843                do k=imin,imax
844                   if (k.ne.jj) lj=lj*(rh(i,kts,j)-rhs(k))/(rhs(jj)-rhs(k))
845                end do
846                tauaer(i,j,ib)=tauaer(i,j,ib)+lj*raod_lut(aer_type,ib,jj)*aod550(i,j)
847             end do
848          end do
849       end do
850    end do
851 end subroutine calc_spectral_aod_goddard_sw
853 subroutine calc_spectral_ssa_goddard_sw(ims,ime,jms,jme,kms,kme,  &
854                                         its,ite,jts,jte,kts,kte,  &
855                                         rh,aer_type,              &
856                                         ssaaer                    )
857    implicit none
858    
859    ! constants
860    integer, parameter :: N_AER_TYPES=3
861    integer, parameter :: N_RH=8
862    integer, parameter :: N_BANDS=11
863    integer, parameter :: N_INT_POINTS=4
865    ! I/O variables
866    integer, intent(in) :: ims,ime,jms,jme,kms,kme, &
867                           its,ite,jts,jte,kts,kte
868    integer, intent(in) :: aer_type
869    real, dimension(ims:ime, kms:kme, jms:jme),            intent(in)    :: rh     ! surface relative humidity
870    real, dimension(ims:ime, kms:kme, jms:jme, 1:N_BANDS), intent(inout) :: ssaaer ! aerosol single-scattering albedo at surface
872    ! local variables
873    integer :: i,j,k,kk,ib,imax,imin,ii,jj
874    real    :: rhs(N_RH),lj
875    real    :: ssa_lut(N_AER_TYPES,N_BANDS,N_RH)
877    ! relative humidity steps
878    data (rhs(i),i=1,8) /0.,50.,70.,80.,90.,95.,98.,99./
880    ! aer_type = 1 : rural (SF79)
881    data (ssa_lut(1,ib,1),ib=1,N_BANDS) /0.6719,0.8192,0.8657,0.9223,0.9333,0.9400,0.9509,0.9428,0.8932,0.8165,0.7855/
882    data (ssa_lut(1,ib,2),ib=1,N_BANDS) /0.6813,0.8251,0.8703,0.9251,0.9357,0.9421,0.9524,0.9447,0.8963,0.8217,0.7665/
883    data (ssa_lut(1,ib,3),ib=1,N_BANDS) /0.6929,0.8322,0.8760,0.9292,0.9395,0.9457,0.9557,0.9482,0.9048,0.8338,0.7431/
884    data (ssa_lut(1,ib,4),ib=1,N_BANDS) /0.7437,0.8626,0.8997,0.9445,0.9530,0.9581,0.9661,0.9607,0.9263,0.8739,0.7215/
885    data (ssa_lut(1,ib,5),ib=1,N_BANDS) /0.7987,0.8951,0.9249,0.9602,0.9666,0.9704,0.9760,0.9730,0.9498,0.9118,0.7200/
886    data (ssa_lut(1,ib,6),ib=1,N_BANDS) /0.8232,0.9089,0.9353,0.9665,0.9721,0.9754,0.9802,0.9780,0.9595,0.9327,0.7229/
887    data (ssa_lut(1,ib,7),ib=1,N_BANDS) /0.8472,0.9205,0.9435,0.9712,0.9765,0.9797,0.9850,0.9838,0.9710,0.9482,0.7340/
888    data (ssa_lut(1,ib,8),ib=1,N_BANDS) /0.8645,0.9322,0.9530,0.9771,0.9814,0.9839,0.9874,0.9871,0.9780,0.9589,0.7363/
890    ! aer_type = 2: urban (SF79)
891    data (ssa_lut(2,ib,1),ib=1,N_BANDS) /0.5832,0.6144,0.6244,0.6368,0.6393,0.6409,0.6436,0.6375,0.5856,0.4748,0.3908/
892    data (ssa_lut(2,ib,2),ib=1,N_BANDS) /0.5933,0.6244,0.6343,0.6467,0.6492,0.6508,0.6536,0.6477,0.5961,0.4861,0.3957/
893    data (ssa_lut(2,ib,3),ib=1,N_BANDS) /0.6444,0.6823,0.6917,0.6999,0.7010,0.7019,0.7145,0.7552,0.6365,0.5559,0.4342/
894    data (ssa_lut(2,ib,4),ib=1,N_BANDS) /0.7195,0.7490,0.7587,0.7714,0.7743,0.7762,0.7805,0.7796,0.7419,0.6545,0.4987/
895    data (ssa_lut(2,ib,5),ib=1,N_BANDS) /0.7791,0.8071,0.8163,0.8287,0.8316,0.8336,0.8384,0.8413,0.8156,0.7477,0.5632/
896    data (ssa_lut(2,ib,6),ib=1,N_BANDS) /0.8230,0.8486,0.8573,0.8690,0.8718,0.8738,0.8789,0.8843,0.8691,0.8202,0.6178/
897    data (ssa_lut(2,ib,7),ib=1,N_BANDS) /0.8653,0.8877,0.8953,0.9059,0.9084,0.9103,0.9155,0.9232,0.9180,0.8895,0.6733/
898    data (ssa_lut(2,ib,8),ib=1,N_BANDS) /0.8861,0.9064,0.9134,0.9233,0.9257,0.9276,0.9327,0.9413,0.9410,0.9234,0.6987/
900    ! aer_type = 3 : maritime (SF79)
901    data (ssa_lut(3,ib,1),ib=1,N_BANDS) /0.7849,0.8895,0.9221,0.9610,0.9683,0.9727,0.9802,0.9829,0.9793,0.9742,0.9493/
902    data (ssa_lut(3,ib,2),ib=1,N_BANDS) /0.7965,0.8962,0.9271,0.9639,0.9707,0.9748,0.9816,0.9843,0.9814,0.9765,0.9043/
903    data (ssa_lut(3,ib,3),ib=1,N_BANDS) /0.8250,0.9116,0.9384,0.9701,0.9760,0.9795,0.9852,0.9876,0.9861,0.9819,0.8586/
904    data (ssa_lut(3,ib,4),ib=1,N_BANDS) /0.8907,0.9464,0.9635,0.9834,0.9869,0.9890,0.9922,0.9939,0.9935,0.9892,0.8278/
905    data (ssa_lut(3,ib,5),ib=1,N_BANDS) /0.9147,0.9595,0.9730,0.9884,0.9911,0.9925,0.9944,0.9956,0.9954,0.9905,0.8233/
906    data (ssa_lut(3,ib,6),ib=1,N_BANDS) /0.9336,0.9689,0.9795,0.9915,0.9935,0.9945,0.9959,0.9969,0.9967,0.9908,0.8178/
907    data (ssa_lut(3,ib,7),ib=1,N_BANDS) /0.9543,0.9788,0.9861,0.9944,0.9958,0.9965,0.9974,0.9980,0.9978,0.9904,0.8122/
908    data (ssa_lut(3,ib,8),ib=1,N_BANDS) /0.9669,0.9848,0.9902,0.9962,0.9971,0.9976,0.9982,0.9986,0.9985,0.9892,0.8062/
910    do j=jts,jte
911       do i=its,ite
912          do k=kts,kte
913             ! common part of the Lagrange's interpolator
914             !   only depends on the relative humidity value
915             ii=1
916             do while ( (ii.le.N_RH) .and. (rh(i,k,j).gt.rhs(ii)) )
917                ii=ii+1
918             end do
919             imin=max(1,ii-N_INT_POINTS/2-1)
920             imax=min(N_RH,ii+N_INT_POINTS/2)
922             do ib=1,N_BANDS
923                ssaaer(i,k,j,ib)=0.
924                do jj=imin,imax
925                   lj=1.
926                   do kk=imin,imax
927                      if (kk.ne.jj) lj=lj*(rh(i,k,j)-rhs(kk))/(rhs(jj)-rhs(kk))
928                   end do
929                   ssaaer(i,k,j,ib)=ssaaer(i,k,j,ib)+lj*ssa_lut(aer_type,ib,jj)
930                end do
931             end do
932          end do
933       end do
934    end do
935 end subroutine calc_spectral_ssa_goddard_sw
937 subroutine calc_spectral_asy_goddard_sw(ims,ime,jms,jme,kms,kme,  &
938                                         its,ite,jts,jte,kts,kte,  &
939                                         rh,aer_type,              &
940                                         asyaer                    )
941    implicit none
942    
943    ! constants
944    integer, parameter :: N_AER_TYPES=3
945    integer, parameter :: N_RH=8
946    integer, parameter :: N_BANDS=11
947    integer, parameter :: N_INT_POINTS=4
949    ! I/O variables
950    integer, intent(in) :: ims,ime,jms,jme,kms,kme, &
951                           its,ite,jts,jte,kts,kte
952    integer, intent(in) :: aer_type
953    real, dimension(ims:ime, kms:kme, jms:jme),            intent(in)    :: rh ! surface relative humidity
954    real, dimension(ims:ime, kms:kme, jms:jme, 1:N_BANDS), intent(inout) :: asyaer ! aerosol asymmetry parameter at surface
956    ! local variables
957    integer :: i,j,k,kk,ib,imax,imin,ii,jj
958    real    :: rhs(N_RH),lj
959    real    :: asy_lut(N_AER_TYPES,N_BANDS,N_RH)
961    ! relative humidity steps
962    data (rhs(i),i=1,8) /0.,50.,70.,80.,90.,95.,98.,99./
964    ! aer_type = 1 : rural (SF79)
965    data (asy_lut(1,ib,1),ib=1,N_BANDS) /0.7602,0.7152,0.7007,0.6820,0.6778,0.6750,0.6675,0.6485,0.6231,0.6481,0.7524/
966    data (asy_lut(1,ib,2),ib=1,N_BANDS) /0.7617,0.7191,0.7052,0.6870,0.6829,0.6800,0.6724,0.6538,0.6277,0.6509,0.7543/
967    data (asy_lut(1,ib,3),ib=1,N_BANDS) /0.7651,0.7250,0.7122,0.6958,0.6922,0.6898,0.6831,0.6645,0.6377,0.6575,0.7592/
968    data (asy_lut(1,ib,4),ib=1,N_BANDS) /0.7738,0.7460,0.7372,0.7260,0.7236,0.7219,0.7172,0.7003,0.6725,0.6806,0.7635/
969    data (asy_lut(1,ib,5),ib=1,N_BANDS) /0.7765,0.7627,0.7581,0.7518,0.7502,0.7490,0.7453,0.7312,0.7026,0.6949,0.7539/
970    data (asy_lut(1,ib,6),ib=1,N_BANDS) /0.7761,0.7683,0.7656,0.7615,0.7604,0.7595,0.7564,0.7438,0.7157,0.7019,0.7491/
971    data (asy_lut(1,ib,7),ib=1,N_BANDS) /0.7773,0.7754,0.7746,0.7730,0.7725,0.7720,0.7702,0.7610,0.7383,0.7279,0.7747/
972    data (asy_lut(1,ib,8),ib=1,N_BANDS) /0.7777,0.7792,0.7795,0.7794,0.7793,0.7790,0.7780,0.7714,0.7514,0.7405,0.7815/
974    ! aer_type = 2: urban (SF79)
975    data (asy_lut(2,ib,1),ib=1,N_BANDS) /0.7794,0.7504,0.7395,0.7223,0.7173,0.7132,0.6996,0.6640,0.6250,0.6406,0.7319/
976    data (asy_lut(2,ib,2),ib=1,N_BANDS) /0.7820,0.7543,0.7439,0.7275,0.7227,0.7188,0.7057,0.6708,0.6317,0.6458,0.7351/
977    data (asy_lut(2,ib,3),ib=1,N_BANDS) /0.7912,0.7711,0.7634,0.7507,0.7469,0.7437,0.7327,0.7016,0.6631,0.6691,0.7477/
978    data (asy_lut(2,ib,4),ib=1,N_BANDS) /0.7951,0.7851,0.7809,0.7733,0.7708,0.7686,0.7607,0.7354,0.6979,0.6898,0.7486/
979    data (asy_lut(2,ib,5),ib=1,N_BANDS) /0.7923,0.7905,0.7892,0.7858,0.7844,0.7831,0.7778,0.7581,0.7238,0.7056,0.7460/
980    data (asy_lut(2,ib,6),ib=1,N_BANDS) /0.7879,0.7924,0.7932,0.7927,0.7920,0.7913,0.7877,0.7733,0.7431,0.7201,0.7458/
981    data (asy_lut(2,ib,7),ib=1,N_BANDS) /0.7828,0.7925,0.7951,0.7973,0.7974,0.7972,0.7957,0.7872,0.7638,0.7396,0.7523/
982    data (asy_lut(2,ib,8),ib=1,N_BANDS) /0.7808,0.7921,0.7954,0.7989,0.7994,0.7995,0.7993,0.7944,0.7754,0.7677,0.7622/
984    ! aer_type = 3 : maritime (SF79)
985    data (asy_lut(3,ib,1),ib=1,N_BANDS) /0.7532,0.7204,0.7102,0.6980,0.6956,0.6940,0.6896,0.6780,0.6814,0.6953,0.6927/
986    data (asy_lut(3,ib,2),ib=1,N_BANDS) /0.7604,0.7321,0.7225,0.7092,0.7060,0.7036,0.6978,0.6918,0.6945,0.7119,0.7217/
987    data (asy_lut(3,ib,3),ib=1,N_BANDS) /0.7717,0.7491,0.7415,0.7309,0.7284,0.7266,0.7226,0.7210,0.7273,0.7466,0.7677/
988    data (asy_lut(3,ib,4),ib=1,N_BANDS) /0.7956,0.7875,0.7844,0.7794,0.7779,0.7768,0.7738,0.7717,0.7753,0.7773,0.8198/
989    data (asy_lut(3,ib,5),ib=1,N_BANDS) /0.8009,0.7963,0.7948,0.7927,0.7922,0.7918,0.7907,0.7868,0.7897,0.7975,0.8295/
990    data (asy_lut(3,ib,6),ib=1,N_BANDS) /0.8081,0.8051,0.8043,0.8036,0.8035,0.8035,0.8031,0.7990,0.7935,0.8014,0.8386/
991    data (asy_lut(3,ib,7),ib=1,N_BANDS) /0.8143,0.8179,0.8186,0.8185,0.8181,0.8176,0.8156,0.8107,0.8041,0.8071,0.8482/
992    data (asy_lut(3,ib,8),ib=1,N_BANDS) /0.8205,0.8252,0.8264,0.8271,0.8269,0.8267,0.8251,0.8202,0.8131,0.8119,0.8547/
994    do j=jts,jte
995       do i=its,ite
996          do k=kts,kte
997             ! common part of the Lagrange's interpolator
998             !   only depends on the relative humidity value
999             ii=1
1000             do while ( (ii.le.N_RH) .and. (rh(i,k,j).gt.rhs(ii)) )
1001                ii=ii+1
1002             end do
1003             imin=max(1,ii-N_INT_POINTS/2-1)
1004             imax=min(N_RH,ii+N_INT_POINTS/2)
1006             do ib=1,N_BANDS
1007                asyaer(i,k,j,ib)=0.
1008                do jj=imin,imax
1009                   lj=1.
1010                   do kk=imin,imax
1011                      if (kk.ne.jj) lj=lj*(rh(i,k,j)-rhs(kk))/(rhs(jj)-rhs(kk))
1012                   end do
1013                   asyaer(i,k,j,ib)=asyaer(i,k,j,ib)+lj*asy_lut(aer_type,ib,jj)
1014                end do
1015             end do
1016          end do
1017       end do
1018    end do
1019 end subroutine calc_spectral_asy_goddard_sw
1021 subroutine calc_spectral_aod_rrtmg_sw(ims,ime,jms,jme,kms,kme,          &
1022                                       its,ite,jts,jte,kts,kte,          &
1023                                       rh,aer_type,aod550,               &
1024                                       tauaer,                           &
1025                                       aod550_3d, tauaer3d)               ! trude
1027    implicit none
1028    
1029    ! constants
1030    integer, parameter :: N_AER_TYPES=3
1031    integer, parameter :: N_RH=8
1032    integer, parameter :: N_BANDS=14
1033    integer, parameter :: N_INT_POINTS=4
1035    ! I/O variables
1036    integer, intent(in) :: ims,ime,jms,jme,kms,kme, &
1037                           its,ite,jts,jte,kts,kte
1038    integer, intent(in) :: aer_type
1039    real, dimension(ims:ime, kms:kme, jms:jme),   intent(in) :: rh        ! relative humidity
1040    real, dimension(ims:ime, jms:jme),            intent(in) :: aod550    ! Total AOD at 550 nm at surface
1041    real, dimension(ims:ime, jms:jme, 1:N_BANDS), intent(inout) :: tauaer ! Total spectral aerosol optical depth at surface
1043    ! ++ Trude
1044    real, dimension(ims:ime, kms:kme, jms:jme), optional, intent(in) :: aod550_3d ! 3D  AOD at 550 nm 
1045    real, dimension(ims:ime, kms:kme, jms:jme, 1:N_BANDS), optional, intent(inout) :: tauaer3d ! 
1046    ! -- Trude
1048    ! local variables
1049    integer :: i,j,k,ib,imax,imin,ii,jj,kk
1050    real    :: rhs(N_RH),lj
1051    real    :: raod_lut(N_AER_TYPES,N_BANDS,N_RH)
1053    ! relative humidity steps
1054    data (rhs(i),i=1,8) /0.,50.,70.,80.,90.,95.,98.,99./
1056    ! aer_type = 1 : rural (SF79)
1057    data (raod_lut(1,ib,1),ib=1,N_BANDS) /0.0735,0.0997,0.1281,0.1529,0.1882,0.2512,0.3010,0.4550,0.7159,1.0357, &
1058                                          1.3582,1.6760,2.2523,0.0582/
1059    data (raod_lut(1,ib,2),ib=1,N_BANDS) /0.0741,0.1004,0.1289,0.1537,0.1891,0.2522,0.3021,0.4560,0.7166,1.0351, &
1060                                          1.3547,1.6687,2.2371,0.0587/
1061    data (raod_lut(1,ib,3),ib=1,N_BANDS) /0.0752,0.1017,0.1304,0.1554,0.1909,0.2542,0.3042,0.4580,0.7179,1.0342, &
1062                                          1.3485,1.6559,2.2102,0.0596/
1063    data (raod_lut(1,ib,4),ib=1,N_BANDS) /0.0766,0.1034,0.1323,0.1575,0.1932,0.2567,0.3068,0.4605,0.7196,1.0332, &
1064                                          1.3411,1.6407,2.1785,0.0608/
1065    data (raod_lut(1,ib,5),ib=1,N_BANDS) /0.0807,0.1083,0.1379,0.1635,0.1998,0.2639,0.3143,0.4677,0.7244,1.0305, &
1066                                          1.3227,1.6031,2.1006,0.0644/
1067    data (raod_lut(1,ib,6),ib=1,N_BANDS) /0.0884,0.1174,0.1482,0.1746,0.2118,0.2769,0.3277,0.4805,0.7328,1.0272, &
1068                                          1.2977,1.5525,1.9976,0.0712/
1069    data (raod_lut(1,ib,7),ib=1,N_BANDS) /0.1072,0.1391,0.1724,0.2006,0.2396,0.3066,0.3581,0.5087,0.7510,1.0231, &
1070                                          1.2622,1.4818,1.8565,0.0878/
1071    data (raod_lut(1,ib,8),ib=1,N_BANDS) /0.1286,0.1635,0.1991,0.2288,0.2693,0.3377,0.3895,0.5372,0.7686,1.0213, &
1072                                          1.2407,1.4394,1.7739,0.1072/
1074    ! aer_type = 2 : urban (SF79)
1075    data (raod_lut(2,ib,1),ib=1,N_BANDS) /0.1244,0.1587,0.1939,0.2233,0.2635,0.3317,0.3835,0.5318,0.7653,1.0344, &
1076                                          1.3155,1.5885,2.0706,0.1033/
1077    data (raod_lut(2,ib,2),ib=1,N_BANDS) /0.1159,0.1491,0.1834,0.2122,0.2518,0.3195,0.3712,0.5207,0.7585,1.0331, &
1078                                          1.3130,1.5833,2.0601,0.0956/
1079    data (raod_lut(2,ib,3),ib=1,N_BANDS) /0.1093,0.1416,0.1752,0.2035,0.2427,0.3099,0.3615,0.5118,0.7529,1.0316, &
1080                                          1.3083,1.5739,2.0408,0.0898/
1081    data (raod_lut(2,ib,4),ib=1,N_BANDS) /0.1062,0.1381,0.1712,0.1993,0.2382,0.3052,0.3567,0.5074,0.7501,1.0302, &
1082                                          1.3025,1.5620,2.0168,0.0870/
1083    data (raod_lut(2,ib,5),ib=1,N_BANDS) /0.1045,0.1361,0.1690,0.1970,0.2357,0.3025,0.3540,0.5049,0.7486,1.0271, &
1084                                          1.2864,1.5297,1.9518,0.0854/
1085    data (raod_lut(2,ib,6),ib=1,N_BANDS) /0.1065,0.1384,0.1716,0.1997,0.2386,0.3056,0.3571,0.5078,0.7504,1.0227, &
1086                                          1.2603,1.4780,1.8492,0.0872/
1087    data (raod_lut(2,ib,7),ib=1,N_BANDS) /0.1147,0.1478,0.1820,0.2107,0.2503,0.3179,0.3696,0.5192,0.7575,1.0146, &
1088                                          1.2116,1.3830,1.6658,0.0946/
1089    data (raod_lut(2,ib,8),ib=1,N_BANDS) /0.1247,0.1590,0.1943,0.2237,0.2639,0.3322,0.3840,0.5322,0.7656,1.0082, &
1090                                          1.1719,1.3075,1.5252,0.1036/
1092    ! aer_type = 3 : maritime (SF79)
1093    data (raod_lut(3,ib,1),ib=1,N_BANDS) /0.3053,0.3507,0.3932,0.4261,0.4681,0.5334,0.5797,0.6962,0.8583,1.0187, &
1094                                          1.1705,1.3049,1.5205,0.2748/
1095    data (raod_lut(3,ib,2),ib=1,N_BANDS) /0.3566,0.4023,0.4443,0.4765,0.5170,0.5792,0.6227,0.7298,0.8756,1.0162, &
1096                                          1.1472,1.2614,1.4415,0.3256/
1097    data (raod_lut(3,ib,3),ib=1,N_BANDS) /0.4359,0.4803,0.5203,0.5505,0.5879,0.6441,0.6828,0.7756,0.8985,1.0135, &
1098                                          1.1198,1.2109,1.3518,0.4051/
1099    data (raod_lut(3,ib,4),ib=1,N_BANDS) /0.5128,0.5544,0.5913,0.6187,0.6523,0.7020,0.7358,0.8149,0.9174,1.0115, &
1100                                          1.0995,1.1740,1.2875,0.4835/
1101    data (raod_lut(3,ib,5),ib=1,N_BANDS) /0.6479,0.6816,0.7108,0.7320,0.7575,0.7946,0.8193,0.8752,0.9455,1.0092, &
1102                                          1.0728,1.1263,1.2061,0.6236/
1103    data (raod_lut(3,ib,6),ib=1,N_BANDS) /0.7582,0.7831,0.8043,0.8196,0.8377,0.8636,0.8806,0.9184,0.9649,1.0080, &
1104                                          1.0564,1.0973,1.1576,0.7399/
1105    data (raod_lut(3,ib,7),ib=1,N_BANDS) /0.8482,0.8647,0.8785,0.8884,0.9000,0.9164,0.9272,0.9506,0.9789,1.0072, &
1106                                          1.0454,1.0780,1.1256,0.8360/
1107    data (raod_lut(3,ib,8),ib=1,N_BANDS) /0.8836,0.8965,0.9073,0.9149,0.9239,0.9365,0.9448,0.9626,0.9841,1.0069, &
1108                                          1.0415,1.0712,1.1145,0.8741/
1111 !   ++ Trude    ;  if 3D AOD, disaggreaget at all levels.
1112    if (present(aod550_3d)) then
1113       do j=jts,jte
1114          do i=its,ite
1115             !   common part of the Lagrange's interpolator
1116             !   only depends on the relative humidity value
1117             do kk = kts,kte
1118                ii=1
1119                do while ( (ii.le.N_RH) .and. (rh(i,kk,j).gt.rhs(ii)) )
1120                   ii=ii+1
1121                end do
1122                imin=max(1,ii-N_INT_POINTS/2-1)
1123                imax=min(N_RH,ii+N_INT_POINTS/2)
1125                do ib=1,N_BANDS
1126                   tauaer3d(i,kk,j,ib)=0.
1127                   do jj=imin,imax
1128                      lj=1.
1129                      do k=imin,imax
1130                         if (k.ne.jj) lj=lj*(rh(i,kk,j)-rhs(k))/(rhs(jj)-rhs(k))
1131                      end do
1132                      tauaer3d(i,kk,j,ib)=tauaer3d(i,kk,j,ib)+lj*raod_lut(aer_type,ib,jj)*aod550_3d(i,kk,j)
1133                   end do
1134                end do
1135             end do
1136          end do
1137       end do
1138     else
1139 ! -- Trude
1141    do j=jts,jte
1142       do i=its,ite
1143          ! common part of the Lagrange's interpolator
1144          !   only depends on the relative humidity value
1145          ii=1
1146          do while ( (ii.le.N_RH) .and. (rh(i,kts,j).gt.rhs(ii)) )
1147             ii=ii+1
1148          end do
1149          imin=max(1,ii-N_INT_POINTS/2-1)
1150          imax=min(N_RH,ii+N_INT_POINTS/2)
1152          do ib=1,N_BANDS
1153             tauaer(i,j,ib)=0.
1154             do jj=imin,imax
1155                lj=1.
1156                do k=imin,imax
1157                   if (k.ne.jj) lj=lj*(rh(i,kts,j)-rhs(k))/(rhs(jj)-rhs(k))
1158                end do
1159                tauaer(i,j,ib)=tauaer(i,j,ib)+lj*raod_lut(aer_type,ib,jj)*aod550(i,j)
1160             end do
1161          end do
1162       end do
1163    end do
1164    endif
1166 end subroutine calc_spectral_aod_rrtmg_sw
1168 subroutine calc_spectral_ssa_rrtmg_sw(ims,ime,jms,jme,kms,kme,  &
1169                                       its,ite,jts,jte,kts,kte,  &
1170                                       rh,aer_type,              &
1171                                       ssaaer                    )
1172    implicit none
1173    
1174    ! constants
1175    integer, parameter :: N_AER_TYPES=3
1176    integer, parameter :: N_RH=8
1177    integer, parameter :: N_BANDS=14
1178    integer, parameter :: N_INT_POINTS=4
1180    ! I/O variables
1181    integer, intent(in) :: ims,ime,jms,jme,kms,kme, &
1182                           its,ite,jts,jte,kts,kte
1183    integer, intent(in) :: aer_type
1184    real, dimension(ims:ime, kms:kme, jms:jme),            intent(in)    :: rh     ! surface relative humidity
1185    real, dimension(ims:ime, kms:kme, jms:jme, 1:N_BANDS), intent(inout) :: ssaaer ! aerosol single-scattering albedo at surface
1187    ! local variables
1188    integer :: i,j,k,kk,ib,imax,imin,ii,jj
1189    real    :: rhs(N_RH),lj
1190    real    :: ssa_lut(N_AER_TYPES,N_BANDS,N_RH)
1192    ! relative humidity steps
1193    data (rhs(i),i=1,8) /0.,50.,70.,80.,90.,95.,98.,99./
1195    ! aer_type = 1 : rural (SF79)
1196    data (ssa_lut(1,ib,1),ib=1,N_BANDS) /0.8730,0.6695,0.8530,0.8601,0.8365,0.7949,0.8113,0.8810,0.9305,0.9436, &
1197                                         0.9532,0.9395,0.8007,0.8634/
1198    data (ssa_lut(1,ib,2),ib=1,N_BANDS) /0.8428,0.6395,0.8571,0.8645,0.8408,0.8007,0.8167,0.8845,0.9326,0.9454, &
1199                                         0.9545,0.9416,0.8070,0.8589/
1200    data (ssa_lut(1,ib,3),ib=1,N_BANDS) /0.8000,0.6025,0.8668,0.8740,0.8503,0.8140,0.8309,0.8943,0.9370,0.9489, &
1201                                         0.9577,0.9451,0.8146,0.8548/
1202    data (ssa_lut(1,ib,4),ib=1,N_BANDS) /0.7298,0.5666,0.9030,0.9049,0.8863,0.8591,0.8701,0.9178,0.9524,0.9612, &
1203                                         0.9677,0.9576,0.8476,0.8578/
1204    data (ssa_lut(1,ib,5),ib=1,N_BANDS) /0.7010,0.5606,0.9312,0.9288,0.9183,0.9031,0.9112,0.9439,0.9677,0.9733, &
1205                                         0.9772,0.9699,0.8829,0.8590/
1206    data (ssa_lut(1,ib,6),ib=1,N_BANDS) /0.6933,0.5620,0.9465,0.9393,0.9346,0.9290,0.9332,0.9549,0.9738,0.9782, &
1207                                         0.9813,0.9750,0.8980,0.8594/
1208    data (ssa_lut(1,ib,7),ib=1,N_BANDS) /0.6842,0.5843,0.9597,0.9488,0.9462,0.9470,0.9518,0.9679,0.9808,0.9839, &
1209                                         0.9864,0.9794,0.9113,0.8648/
1210    data (ssa_lut(1,ib,8),ib=1,N_BANDS) /0.6786,0.5897,0.9658,0.9522,0.9530,0.9610,0.9651,0.9757,0.9852,0.9871, &
1211                                         0.9883,0.9835,0.9236,0.8618/
1213    ! aer_type = 2: urban (SF79)
1214    data (ssa_lut(2,ib,1),ib=1,N_BANDS) /0.4063,0.3663,0.4093,0.4205,0.4487,0.4912,0.5184,0.5743,0.6233,0.6392, &
1215                                         0.6442,0.6408,0.6105,0.4094/
1216    data (ssa_lut(2,ib,2),ib=1,N_BANDS) /0.4113,0.3654,0.4215,0.4330,0.4604,0.5022,0.5293,0.5848,0.6336,0.6493, &
1217                                         0.6542,0.6507,0.6205,0.4196/
1218    data (ssa_lut(2,ib,3),ib=1,N_BANDS) /0.4500,0.3781,0.4924,0.5050,0.5265,0.5713,0.6048,0.6274,0.6912,0.7714, &
1219                                         0.7308,0.7027,0.6772,0.4820/
1220    data (ssa_lut(2,ib,4),ib=1,N_BANDS) /0.5075,0.4139,0.5994,0.6127,0.6350,0.6669,0.6888,0.7333,0.7704,0.7809, &
1221                                         0.7821,0.7762,0.7454,0.5709/
1222    data (ssa_lut(2,ib,5),ib=1,N_BANDS) /0.5596,0.4570,0.7009,0.7118,0.7317,0.7583,0.7757,0.8093,0.8361,0.8422, &
1223                                         0.8406,0.8337,0.8036,0.6525/
1224    data (ssa_lut(2,ib,6),ib=1,N_BANDS) /0.6008,0.4971,0.7845,0.7906,0.8075,0.8290,0.8418,0.8649,0.8824,0.8849, &
1225                                         0.8815,0.8739,0.8455,0.7179/
1226    data (ssa_lut(2,ib,7),ib=1,N_BANDS) /0.6401,0.5407,0.8681,0.8664,0.8796,0.8968,0.9043,0.9159,0.9244,0.9234, &
1227                                         0.9182,0.9105,0.8849,0.7796/
1228    data (ssa_lut(2,ib,8),ib=1,N_BANDS) /0.6567,0.5618,0.9073,0.9077,0.9182,0.9279,0.9325,0.9398,0.9440,0.9413, &
1229                                         0.9355,0.9278,0.9039,0.8040/
1231    ! aer_type = 3 : maritime (SF79)
1232    data (ssa_lut(3,ib,1),ib=1,N_BANDS) /0.9697,0.9183,0.9749,0.9820,0.9780,0.9712,0.9708,0.9778,0.9831,0.9827, &
1233                                         0.9826,0.9723,0.8763,0.9716/
1234    data (ssa_lut(3,ib,2),ib=1,N_BANDS) /0.9070,0.8491,0.9730,0.9816,0.9804,0.9742,0.9738,0.9802,0.9847,0.9841, &
1235                                         0.9838,0.9744,0.8836,0.9546/
1236    data (ssa_lut(3,ib,3),ib=1,N_BANDS) /0.8378,0.7761,0.9797,0.9827,0.9829,0.9814,0.9812,0.9852,0.9882,0.9875, &
1237                                         0.9871,0.9791,0.9006,0.9348/
1238    data (ssa_lut(3,ib,4),ib=1,N_BANDS) /0.7866,0.7249,0.9890,0.9822,0.9856,0.9917,0.9924,0.9932,0.9943,0.9938, &
1239                                         0.9933,0.9887,0.9393,0.9204/
1240    data (ssa_lut(3,ib,5),ib=1,N_BANDS) /0.7761,0.7164,0.9959,0.9822,0.9834,0.9941,0.9955,0.9952,0.9960,0.9956, &
1241                                         0.9951,0.9922,0.9538,0.9152/
1242    data (ssa_lut(3,ib,6),ib=1,N_BANDS) /0.7671,0.7114,0.9902,0.9786,0.9838,0.9954,0.9970,0.9965,0.9971,0.9968, &
1243                                         0.9964,0.9943,0.9644,0.9158/
1244    data (ssa_lut(3,ib,7),ib=1,N_BANDS) /0.7551,0.7060,0.9890,0.9743,0.9807,0.9966,0.9989,0.9978,0.9982,0.9980, &
1245                                         0.9978,0.9964,0.9757,0.9122/
1246    data (ssa_lut(3,ib,8),ib=1,N_BANDS) /0.7439,0.7000,0.9870,0.9695,0.9769,0.9970,1.0000,0.9984,0.9988,0.9986, &
1247                                         0.9984,0.9975,0.9825,0.9076/
1249    do j=jts,jte
1250       do i=its,ite
1251          do k=kts,kte
1252             ! common part of the Lagrange's interpolator
1253             !   only depends on the relative humidity value
1254             ii=1
1255             do while ( (ii.le.N_RH) .and. (rh(i,k,j).gt.rhs(ii)) )
1256                ii=ii+1
1257             end do
1258             imin=max(1,ii-N_INT_POINTS/2-1)
1259             imax=min(N_RH,ii+N_INT_POINTS/2)
1261             do ib=1,N_BANDS
1262                ssaaer(i,k,j,ib)=0.
1263                do jj=imin,imax
1264                   lj=1.
1265                   do kk=imin,imax
1266                      if (kk.ne.jj) lj=lj*(rh(i,k,j)-rhs(kk))/(rhs(jj)-rhs(kk))
1267                   end do
1268                   ssaaer(i,k,j,ib)=ssaaer(i,k,j,ib)+lj*ssa_lut(aer_type,ib,jj)
1269                end do
1270             end do
1271          end do
1272       end do
1273    end do
1274 end subroutine calc_spectral_ssa_rrtmg_sw
1276 subroutine calc_spectral_asy_rrtmg_sw(ims,ime,jms,jme,kms,kme,  &
1277                                       its,ite,jts,jte,kts,kte,  &
1278                                       rh,aer_type,              &
1279                                       asyaer                    )
1280    implicit none
1281    
1282    ! constants
1283    integer, parameter :: N_AER_TYPES=3
1284    integer, parameter :: N_RH=8
1285    integer, parameter :: N_BANDS=14
1286    integer, parameter :: N_INT_POINTS=4
1288    ! I/O variables
1289    integer, intent(in) :: ims,ime,jms,jme,kms,kme, &
1290                           its,ite,jts,jte,kts,kte
1291    integer, intent(in) :: aer_type
1292    real, dimension(ims:ime, kms:kme, jms:jme),            intent(in)    :: rh ! surface relative humidity
1293    real, dimension(ims:ime, kms:kme, jms:jme, 1:N_BANDS), intent(inout) :: asyaer ! aerosol asymmetry parameter at surface
1295    ! local variables
1296    integer :: i,j,k,kk,ib,imax,imin,ii,jj
1297    real    :: rhs(N_RH),lj
1298    real    :: asy_lut(N_AER_TYPES,N_BANDS,N_RH)
1300    ! relative humidity steps
1301    data (rhs(i),i=1,8) /0.,50.,70.,80.,90.,95.,98.,99./
1303    ! aer_type = 1 : rural (SF79)
1304    data (asy_lut(1,ib,1),ib=1,N_BANDS) /0.7444,0.7711,0.7306,0.7103,0.6693,0.6267,0.6169,0.6207,0.6341,0.6497, &
1305                                         0.6630,0.6748,0.7208,0.7419/
1306    data (asy_lut(1,ib,2),ib=1,N_BANDS) /0.7444,0.7747,0.7314,0.7110,0.6711,0.6301,0.6210,0.6251,0.6392,0.6551, &
1307                                         0.6680,0.6799,0.7244,0.7436/
1308    data (asy_lut(1,ib,3),ib=1,N_BANDS) /0.7438,0.7845,0.7341,0.7137,0.6760,0.6381,0.6298,0.6350,0.6497,0.6657, &
1309                                         0.6790,0.6896,0.7300,0.7477/
1310    data (asy_lut(1,ib,4),ib=1,N_BANDS) /0.7336,0.7934,0.7425,0.7217,0.6925,0.6665,0.6616,0.6693,0.6857,0.7016, &
1311                                         0.7139,0.7218,0.7495,0.7574/
1312    data (asy_lut(1,ib,5),ib=1,N_BANDS) /0.7111,0.7865,0.7384,0.7198,0.6995,0.6864,0.6864,0.6987,0.7176,0.7326, &
1313                                         0.7427,0.7489,0.7644,0.7547/
1314    data (asy_lut(1,ib,6),ib=1,N_BANDS) /0.7009,0.7828,0.7366,0.7196,0.7034,0.6958,0.6979,0.7118,0.7310,0.7452, &
1315                                         0.7542,0.7593,0.7692,0.7522/
1316    data (asy_lut(1,ib,7),ib=1,N_BANDS) /0.7226,0.8127,0.7621,0.7434,0.7271,0.7231,0.7248,0.7351,0.7506,0.7622, &
1317                                         0.7688,0.7719,0.7756,0.7706/
1318    data (asy_lut(1,ib,8),ib=1,N_BANDS) /0.7296,0.8219,0.7651,0.7513,0.7404,0.7369,0.7386,0.7485,0.7626,0.7724, &
1319                                         0.7771,0.7789,0.7790,0.7760/
1321    ! aer_type = 2: urban (SF79)
1322    data (asy_lut(2,ib,1),ib=1,N_BANDS) /0.7399,0.7372,0.7110,0.6916,0.6582,0.6230,0.6147,0.6214,0.6412,0.6655, &
1323                                         0.6910,0.7124,0.7538,0.7395/
1324    data (asy_lut(2,ib,2),ib=1,N_BANDS) /0.7400,0.7419,0.7146,0.6952,0.6626,0.6287,0.6209,0.6280,0.6481,0.6723, &
1325                                         0.6974,0.7180,0.7575,0.7432/
1326    data (asy_lut(2,ib,3),ib=1,N_BANDS) /0.7363,0.7614,0.7303,0.7100,0.6815,0.6550,0.6498,0.6590,0.6802,0.7032, &
1327                                         0.7255,0.7430,0.7735,0.7580/
1328    data (asy_lut(2,ib,4),ib=1,N_BANDS) /0.7180,0.7701,0.7358,0.7163,0.6952,0.6807,0.6801,0.6935,0.7160,0.7370, &
1329                                         0.7553,0.7681,0.7862,0.7623/
1330    data (asy_lut(2,ib,5),ib=1,N_BANDS) /0.7013,0.7733,0.7374,0.7203,0.7057,0.7006,0.7035,0.7192,0.7415,0.7596, &
1331                                         0.7739,0.7827,0.7906,0.7596/
1332    data (asy_lut(2,ib,6),ib=1,N_BANDS) /0.6922,0.7773,0.7404,0.7264,0.7170,0.7179,0.7228,0.7389,0.7595,0.7746, &
1333                                         0.7851,0.7909,0.7918,0.7562/
1334    data (asy_lut(2,ib,7),ib=1,N_BANDS) /0.6928,0.7875,0.7491,0.7393,0.7345,0.7397,0.7455,0.7602,0.7773,0.7883, &
1335                                         0.7944,0.7970,0.7912,0.7555/
1336    data (asy_lut(2,ib,8),ib=1,N_BANDS) /0.7021,0.7989,0.7590,0.7512,0.7613,0.7746,0.7718,0.7727,0.7867,0.7953, &
1337                                         0.7988,0.7994,0.7906,0.7600/
1339    ! aer_type = 3 : maritime (SF79)
1340    data (asy_lut(3,ib,1),ib=1,N_BANDS) /0.6620,0.7011,0.7111,0.7068,0.6990,0.6918,0.6883,0.6827,0.6768,0.6773, &
1341                                         0.6863,0.6940,0.7245,0.6719/
1342    data (asy_lut(3,ib,2),ib=1,N_BANDS) /0.6880,0.7394,0.7297,0.7240,0.7162,0.7083,0.7038,0.6957,0.6908,0.6917, &
1343                                         0.6952,0.7035,0.7356,0.6977/
1344    data (asy_lut(3,ib,3),ib=1,N_BANDS) /0.7266,0.7970,0.7666,0.7593,0.7505,0.7427,0.7391,0.7293,0.7214,0.7210, &
1345                                         0.7212,0.7265,0.7519,0.7340/
1346    data (asy_lut(3,ib,4),ib=1,N_BANDS) /0.7683,0.8608,0.8120,0.8030,0.7826,0.7679,0.7713,0.7760,0.7723,0.7716, &
1347                                         0.7726,0.7767,0.7884,0.7768/
1348    data (asy_lut(3,ib,5),ib=1,N_BANDS) /0.7776,0.8727,0.8182,0.8083,0.7985,0.7939,0.7953,0.7913,0.7846,0.7870, &
1349                                         0.7899,0.7918,0.7969,0.7870/
1350    data (asy_lut(3,ib,6),ib=1,N_BANDS) /0.7878,0.8839,0.8231,0.8130,0.8050,0.7977,0.7945,0.7932,0.7955,0.7992, &
1351                                         0.8025,0.8035,0.8055,0.7956/
1352    data (asy_lut(3,ib,7),ib=1,N_BANDS) /0.8005,0.8957,0.8273,0.8179,0.8105,0.8035,0.8010,0.8030,0.8081,0.8108, &
1353                                         0.8143,0.8174,0.8174,0.8042/
1354    data (asy_lut(3,ib,8),ib=1,N_BANDS) /0.8104,0.9034,0.8294,0.8212,0.8144,0.8087,0.8077,0.8118,0.8175,0.8202, &
1355                                         0.8239,0.8265,0.8246,0.8095/
1357    do j=jts,jte
1358       do i=its,ite
1359          do k=kts,kte
1360             ! common part of the Lagrange's interpolator
1361             !   only depends on the relative humidity value
1362             ii=1
1363             do while ( (ii.le.N_RH) .and. (rh(i,k,j).gt.rhs(ii)) )
1364                ii=ii+1
1365             end do
1366             imin=max(1,ii-N_INT_POINTS/2-1)
1367             imax=min(N_RH,ii+N_INT_POINTS/2)
1369             do ib=1,N_BANDS
1370                asyaer(i,k,j,ib)=0.
1371                do jj=imin,imax
1372                   lj=1.
1373                   do kk=imin,imax
1374                      if (kk.ne.jj) lj=lj*(rh(i,k,j)-rhs(kk))/(rhs(jj)-rhs(kk))
1375                   end do
1376                   asyaer(i,k,j,ib)=asyaer(i,k,j,ib)+lj*asy_lut(aer_type,ib,jj)
1377                end do
1378             end do
1379          end do
1380       end do
1381    end do
1382 end subroutine calc_spectral_asy_rrtmg_sw
1384 subroutine aod_profiler(ht,dz8w,taod550,n_bands,   &
1385                         ims,ime,jms,jme,kms,kme,   &
1386                         its,ite,jts,jte,kts,kte,   &
1387                         aod550                     &
1388                                                    )
1389    use module_wrf_error , only : wrf_err_message
1390    
1391    implicit none
1392    
1393    ! constants
1394    real, parameter :: scale_height=2500. ! meters
1396    ! I/O variables
1397    integer, intent(in) :: n_bands
1398    integer, intent(in) :: ims,ime,jms,jme,kms,kme, &
1399                           its,ite,jts,jte,kts,kte
1400    real, dimension( ims:ime, jms:jme),            intent(in) :: ht
1401    real, dimension( ims:ime, kms:kme, jms:jme ),  intent(in) :: dz8w
1402    real, dimension( ims:ime, jms:jme, 1:n_bands), intent(in) :: taod550
1403    real, dimension( ims:ime, kms:kme, jms:jme, 1:n_bands ), intent(inout) :: aod550
1405    ! local variables
1406    real, dimension(its:ite,kts:kte) :: z2d,aod5502d
1407    real, dimension(its:ite)         :: htoa
1408    real :: aod_scale
1409    real :: aod_acum
1410    integer :: i,j,k,nb
1412    ! input variables from driver are defined such as kms is sfc and
1413    ! kme is toa. Equivalently, kts is sfc and kte is toa
1414    do j=jts,jte
1415       ! heigth profile
1416       !   kts=surface, kte=toa
1417       do i=its,ite
1418          z2d(i,kts)=ht(i,j)+0.5*dz8w(i,kts,j)
1419          do k=kts+1,kte
1420             z2d(i,k)=z2d(i,k-1)+0.5*(dz8w(i,k-1,j)+dz8w(i,k,j))
1421          end do
1422          htoa(i)=z2d(i,kte)+0.5*dz8w(i,kte,j)
1423       end do
1425       do nb=1,n_bands
1426          ! AOD exponential profile
1427          do i=its,ite
1428             aod_scale=taod550(i,j,nb)/(scale_height*(exp(-ht(i,j)/scale_height)-exp(-htoa(i)/scale_height)))
1429             do k=kts,kte
1430                aod550(i,k,j,nb)=aod_scale*dz8w(i,k,j)*exp(-z2d(i,k)/scale_height)
1431             end do
1432          end do 
1433       end do ! nb-loop
1434    end do ! j-loop
1435 end subroutine aod_profiler
1437 subroutine calc_relative_humidity(p,t3d,qv3d,                 &
1438                                   ims,ime,jms,jme,kms,kme,    &
1439                                   its,ite,jts,jte,kts,kte,    &
1440                                   rh                          )
1441    implicit none
1442    
1443    ! I/O variables
1444    integer, intent(in) :: ims,ime,jms,jme,kms,kme,its,ite,jts,jte,kts,kte
1445    ! Naming convention: 8~at => p8w reads as "p-at-w" (w=full levels)
1446    real, dimension(ims:ime, kms:kme, jms:jme), intent(in)    :: p,    & ! pressure (Pa)
1447                                                                 t3d,  & ! temperature (K)
1448                                                                 qv3d    ! water vapor mixing ratio (kg/kg)
1449    real, dimension(ims:ime, kms:kme, jms:jme), intent(inout) :: rh      ! relative humidity at surface
1451    ! local variables
1452    real    :: tc,rv,es,e
1453    integer :: i,j,k
1455    do j=jts,jte
1456       do i=its,ite
1457          do k=kts,kte                               ! only calculations at surface level
1458             tc=t3d(i,k,j)-273.15                    ! temperature (C)
1459             rv=max(0.,qv3d(i,k,j))                  ! water vapor mixing ration (kg kg-1)
1460             es=6.112*exp((17.6*tc)/(tc+243.5))      ! saturation vapor pressure, hPa, Bolton (1980)
1461             e =0.01*rv*p(i,k,j)/(rv+0.62197)        ! vapor pressure, hPa, (ECMWF handouts, page 6, Atmosph. Thermdyn.)
1462                                                     ! rv=eps * e/(p-e) -> e=p * rv/(rv+eps), eps=0.62197
1463             rh(i,k,j)=min(99.,max(0.,100.*e/es))    ! relative humidity (%)
1464          end do
1465       end do
1466    end do
1467 end subroutine calc_relative_humidity
1469 end module module_ra_aerosol