1 ! Aerosol optical properties from AOD at 550 nm, aerosol type and relative humidity
2 ! Author: Jose A. Ruiz-Arias
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
20 !--------------------------------------------------------------
22 !--------------------------------------------------------------
23 ! kms:kme define the range for full-level indices
24 ! kts:kte define the range for half-level indices
26 ! kms=1 is the first full level at surface
27 ! kts=1 is the first half level at surface
29 ! kme is the last full level at toa
30 ! kte is the last half level at toa
32 ! There is one more full level than half levels.
33 ! Therefore, kme=kte+1. I checked it in one of my
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 !--------------------------------------------------------------
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
77 integer, parameter :: N_BANDS=11
79 ! mid-band wavelength / 550. nm
80 real :: lower_wvl(N_BANDS),upper_wvl(N_BANDS)
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/
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
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, &
111 aer_aod550_opt_select: select case(aer_aod550_opt)
113 ! reserved for climatology
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)
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 )
123 aod5502d(i,j)=aer_aod550_val
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)
132 if (minval(aod5502d) .lt. 0) then
133 call wrf_error_fatal('AOD@550 must be positive. Negative value(s) found in auxinput')
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 )
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)
148 ! reserved for climatology
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)
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)
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))
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)))
167 aod_rate=(0.55/(upper_wvl(nb)-lower_wvl(nb)))*log(upper_wvl(nb)/lower_wvl(nb))
171 aod550spc(i,j,nb)=aod5502d(i,j)*aod_rate
177 angexp2d(i,j)=angexp_val
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)
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 )
190 angexp_val=min(2.5,max(-0.3,angexp2d(i,j)))
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)))
196 aod_rate=(0.55/(upper_wvl(nb)-lower_wvl(nb)))*log(upper_wvl(nb)/lower_wvl(nb))
198 aod550spc(i,j,nb)=aod5502d(i,j)*aod_rate
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, &
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
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))
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)
239 ! reserved for climatology
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)
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 )
252 ! no spectral disaggregation
253 ssaaer(i,k,j,nb)=aer_ssa_val
261 aerssa2d(i, j) = aer_ssa_val
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)
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')
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 )
283 ! no spectral disaggregation
284 ssaaer(i,k,j,nb)=aerssa2d(i,j)
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)') &
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, &
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
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)
310 aerssa2d(i,j)=aerssa2d(i,j)/2.
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)
321 ! reserved for climatology
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)
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 )
333 asyaer(i,k,j,nb)=aer_asy_val
341 aerasy2d(i, j) = aer_asy_val
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)
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')
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 )
360 asyaer(i,k,j,nb)=aerasy2d(i,j)
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)') &
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, &
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
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)
386 aerasy2d(i,j)=aerasy2d(i,j)/2.
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
409 integer, parameter :: N_BANDS=14
410 ! local index variables
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/
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
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, &
445 aer_aod550_opt_select: select case(aer_aod550_opt)
447 ! reserved for climatology
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)
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 )
457 aod5502d(i,j)=aer_aod550_val
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)
466 if (minval(aod5502d) .lt. 0) then
467 call wrf_error_fatal('AOD@550 must be positive. Negative value(s) found in auxinput')
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 )
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)
482 ! reserved for climatology
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)
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)
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))
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)))
501 aod_rate=(0.55/(upper_wvl(nb)-lower_wvl(nb)))*log(upper_wvl(nb)/lower_wvl(nb))
505 aod550spc(i,j,nb)=aod5502d(i,j)*aod_rate
511 angexp2d(i,j)=angexp_val
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)
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 )
524 angexp_val=min(2.5,max(-0.3,angexp2d(i,j)))
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)))
530 aod_rate=(0.55/(upper_wvl(nb)-lower_wvl(nb)))*log(upper_wvl(nb)/lower_wvl(nb))
532 aod550spc(i,j,nb)=aod5502d(i,j)*aod_rate
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, &
546 aod5503d, aod550spc3d) ! trude
554 if (present(aod5503d)) then
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))
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))
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
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))
584 angexp2d(i,j)=-(xy/(N_BANDS-3-8+1))/(xx/(N_BANDS-3-8+1))
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
602 tauaer(i,k,j,nb) = aod550spc3d(i,k,j,nb)
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 )
613 aer_ssa_opt_select: select case(aer_ssa_opt)
615 ! reserved for climatology
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)
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 )
628 ! no spectral disaggregation
629 ssaaer(i,k,j,nb)=aer_ssa_val
636 aerssa2d(i,j)=aer_ssa_val
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)
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')
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 )
656 ! no spectral disaggregation
657 ssaaer(i,k,j,nb)=aerssa2d(i,j)
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)') &
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, &
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
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)
683 aerssa2d(i,j)=aerssa2d(i,j)/(N_BANDS-3-8+1)
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)
694 ! reserved for climatology
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)
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 )
706 asyaer(i,k,j,nb)=aer_asy_val
713 aerasy2d(i,j)=aer_asy_val
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)
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')
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 )
732 asyaer(i,k,j,nb)=aerasy2d(i,j)
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)') &
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, &
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
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)
758 aerasy2d(i,j)=aerasy2d(i,j)/(N_BANDS-3-8+1)
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, &
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
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
790 integer :: i,j,k,ib,imax,imin,ii,jj
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)
830 ! common part of the Lagrange's interpolator
831 ! only depends on the relative humidity value
833 do while ( (ii.le.N_RH) .and. (rh(i,kts,j).gt.rhs(ii)) )
836 imin=max(1,ii-N_INT_POINTS/2-1)
837 imax=min(N_RH,ii+N_INT_POINTS/2)
844 if (k.ne.jj) lj=lj*(rh(i,kts,j)-rhs(k))/(rhs(jj)-rhs(k))
846 tauaer(i,j,ib)=tauaer(i,j,ib)+lj*raod_lut(aer_type,ib,jj)*aod550(i,j)
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, &
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
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
873 integer :: i,j,k,kk,ib,imax,imin,ii,jj
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/
913 ! common part of the Lagrange's interpolator
914 ! only depends on the relative humidity value
916 do while ( (ii.le.N_RH) .and. (rh(i,k,j).gt.rhs(ii)) )
919 imin=max(1,ii-N_INT_POINTS/2-1)
920 imax=min(N_RH,ii+N_INT_POINTS/2)
927 if (kk.ne.jj) lj=lj*(rh(i,k,j)-rhs(kk))/(rhs(jj)-rhs(kk))
929 ssaaer(i,k,j,ib)=ssaaer(i,k,j,ib)+lj*ssa_lut(aer_type,ib,jj)
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, &
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
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
957 integer :: i,j,k,kk,ib,imax,imin,ii,jj
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/
997 ! common part of the Lagrange's interpolator
998 ! only depends on the relative humidity value
1000 do while ( (ii.le.N_RH) .and. (rh(i,k,j).gt.rhs(ii)) )
1003 imin=max(1,ii-N_INT_POINTS/2-1)
1004 imax=min(N_RH,ii+N_INT_POINTS/2)
1011 if (kk.ne.jj) lj=lj*(rh(i,k,j)-rhs(kk))/(rhs(jj)-rhs(kk))
1013 asyaer(i,k,j,ib)=asyaer(i,k,j,ib)+lj*asy_lut(aer_type,ib,jj)
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, &
1025 aod550_3d, tauaer3d) ! trude
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
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
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 !
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
1115 ! common part of the Lagrange's interpolator
1116 ! only depends on the relative humidity value
1119 do while ( (ii.le.N_RH) .and. (rh(i,kk,j).gt.rhs(ii)) )
1122 imin=max(1,ii-N_INT_POINTS/2-1)
1123 imax=min(N_RH,ii+N_INT_POINTS/2)
1126 tauaer3d(i,kk,j,ib)=0.
1130 if (k.ne.jj) lj=lj*(rh(i,kk,j)-rhs(k))/(rhs(jj)-rhs(k))
1132 tauaer3d(i,kk,j,ib)=tauaer3d(i,kk,j,ib)+lj*raod_lut(aer_type,ib,jj)*aod550_3d(i,kk,j)
1143 ! common part of the Lagrange's interpolator
1144 ! only depends on the relative humidity value
1146 do while ( (ii.le.N_RH) .and. (rh(i,kts,j).gt.rhs(ii)) )
1149 imin=max(1,ii-N_INT_POINTS/2-1)
1150 imax=min(N_RH,ii+N_INT_POINTS/2)
1157 if (k.ne.jj) lj=lj*(rh(i,kts,j)-rhs(k))/(rhs(jj)-rhs(k))
1159 tauaer(i,j,ib)=tauaer(i,j,ib)+lj*raod_lut(aer_type,ib,jj)*aod550(i,j)
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, &
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
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
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/
1252 ! common part of the Lagrange's interpolator
1253 ! only depends on the relative humidity value
1255 do while ( (ii.le.N_RH) .and. (rh(i,k,j).gt.rhs(ii)) )
1258 imin=max(1,ii-N_INT_POINTS/2-1)
1259 imax=min(N_RH,ii+N_INT_POINTS/2)
1266 if (kk.ne.jj) lj=lj*(rh(i,k,j)-rhs(kk))/(rhs(jj)-rhs(kk))
1268 ssaaer(i,k,j,ib)=ssaaer(i,k,j,ib)+lj*ssa_lut(aer_type,ib,jj)
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, &
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
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
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/
1360 ! common part of the Lagrange's interpolator
1361 ! only depends on the relative humidity value
1363 do while ( (ii.le.N_RH) .and. (rh(i,k,j).gt.rhs(ii)) )
1366 imin=max(1,ii-N_INT_POINTS/2-1)
1367 imax=min(N_RH,ii+N_INT_POINTS/2)
1374 if (kk.ne.jj) lj=lj*(rh(i,k,j)-rhs(kk))/(rhs(jj)-rhs(kk))
1376 asyaer(i,k,j,ib)=asyaer(i,k,j,ib)+lj*asy_lut(aer_type,ib,jj)
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, &
1389 use module_wrf_error , only : wrf_err_message
1394 real, parameter :: scale_height=2500. ! meters
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
1406 real, dimension(its:ite,kts:kte) :: z2d,aod5502d
1407 real, dimension(its:ite) :: htoa
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
1416 ! kts=surface, kte=toa
1418 z2d(i,kts)=ht(i,j)+0.5*dz8w(i,kts,j)
1420 z2d(i,k)=z2d(i,k-1)+0.5*(dz8w(i,k-1,j)+dz8w(i,k,j))
1422 htoa(i)=z2d(i,kte)+0.5*dz8w(i,kte,j)
1426 ! AOD exponential profile
1428 aod_scale=taod550(i,j,nb)/(scale_height*(exp(-ht(i,j)/scale_height)-exp(-htoa(i)/scale_height)))
1430 aod550(i,k,j,nb)=aod_scale*dz8w(i,k,j)*exp(-z2d(i,k)/scale_height)
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, &
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
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 (%)
1467 end subroutine calc_relative_humidity
1469 end module module_ra_aerosol