| blob | f37d167a67ebe44d33e5c99dc226ecc69dda18d4 |
1 !#define SW_DEBUG
3 module module_phot_tuv
5 use params_mod, only : dp, m2km, ppm2vmr, o2vmr, km2cm, m2s
7 IMPLICIT NONE
9 private
10 public :: tuv_driver, tuv_init, tuv_timestep_init
12 real, parameter :: conv1 = km2cm*.5
13 real, parameter :: conv2 = conv1*ppm2vmr
14 real, parameter :: bext340 = 5.E-6
15 real, parameter :: bexth2o = 5.E-6
17 integer :: nj
18 integer :: nlambda_start = 1
19 integer :: nlambda_af_start = 1
20 integer :: nlambda_af_end = 1
21 integer :: nconc, ntemp, nwave
22 integer :: n_temp_data, n_o3_data, n_air_dens_data
23 integer :: j_o2_ndx = -1
24 integer :: last, next
25 integer :: curjulday = 0
26 integer, allocatable :: rxn_ndx(:)
27 real :: dels
28 real :: esfact
29 logical :: has_exo_coldens
30 logical :: tuv_is_initialized = .false.
32 real, allocatable :: temp_tab(:)
33 real, allocatable :: conc_tab(:)
34 real, allocatable :: del_temp_tab(:)
35 real, allocatable :: del_conc_tab(:)
36 real, allocatable :: wl(:)
37 real, allocatable :: wc(:)
38 real, allocatable :: dw(:)
39 real, allocatable :: w_fac(:)
40 real, allocatable :: etfl(:)
41 real, allocatable :: albedo(:)
42 real, allocatable :: o2_xs(:)
43 real, allocatable :: so2_xs(:)
44 real, allocatable :: par_wght(:)
45 real, allocatable :: ery_wght(:)
46 real, allocatable :: z_temp_data(:), temp_data(:)
47 real, allocatable :: z_o3_data(:), o3_data(:)
48 real, allocatable :: z_air_dens_data(:), air_dens_data(:)
49 real, allocatable :: o3_xs_tab(:,:)
50 real, allocatable :: no2_xs_tab(:,:)
51 real, allocatable :: xsqy_tab(:,:,:,:)
52 character(len=32), allocatable :: tuv_jname(:)
53 logical :: has_so2, has_no2
54 logical :: is_full_tuv = .true.
55 ! logical :: is_full_tuv = .false.
56 logical :: rxn_initialized
57 logical, allocatable :: xsqy_is_zdep(:)
59 type column_density
60 integer :: ncoldens_levs
61 integer :: ndays_of_year
62 real(8), pointer :: col_levs(:)
63 real(8), pointer :: day_of_year(:)
64 real(8), pointer :: o3_col_dens(:,:,:,:)
65 real(8), pointer :: o2_col_dens(:,:,:,:)
66 logical :: is_allocated
67 end type column_density
69 type(column_density), allocatable :: col_dens(:)
71 CONTAINS
73 subroutine tuv_driver( &
74 dm, curr_secs, ktau, config_flags, haveaer, &
75 gmt, julday, t_phy, moist, aerwrf, &
76 p8w, t8w, p_phy, chem, rho_phy, &
77 dz8w, xlat, xlong, z, z_at_w, gd_cloud, gd_cloud2, &
78 bl_cld,bl_cld2, &
79 tauaer1,tauaer2,tauaer3,tauaer4, &
80 gaer1,gaer2,gaer3,gaer4, &
81 waer1,waer2,waer3,waer4, &
82 ph_macr,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, &
83 ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
84 ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
85 ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob, &
86 ph_n2o5,ph_o2,ph_n2o,ph_pooh,ph_pan,ph_mvk,ph_hyac, &
87 ph_glyald,ph_mek,ph_gly, &
88 pm2_5_dry,pm2_5_water,uvrad, &
89 dt_cld,af_dir,af_dn,af_up,par,erythema, &
90 o3_gfs_du, &
91 ids,ide, jds,jde, kds,kde, &
92 ims,ime, jms,jme, kms,kme, &
93 its,ite, jts,jte, kts,kte )
95 USE module_configure
96 USE module_state_description
97 USE module_model_constants
98 USE module_data_radm2
99 USE tuv_subs, only : tuv_radfld, sundis, calc_zenith
100 USE module_params, only : kz
101 USE srb, only : sjo2
102 USE module_rxn, only : xsqy_table => xsqy_tab, the_subs, set_initialization
103 USE module_rxn, only : get_initialization
104 USE module_xsections, only : o3xs, no2xs_jpl06a
106 !---------------------------------------------------------------------
107 ! ... dummy arguments
108 !---------------------------------------------------------------------
109 INTEGER, INTENT(IN ) :: dm,julday, &
110 ids,ide, jds,jde, kds,kde, &
111 ims,ime, jms,jme, kms,kme, &
112 its,ite, jts,jte, kts,kte
113 INTEGER, INTENT(IN ) :: ktau
114 REAL(KIND=8), INTENT(IN ) :: curr_secs
115 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), &
116 INTENT(IN ) :: moist
117 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
118 INTENT(INOUT ) :: &
119 pm2_5_dry,pm2_5_water
120 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
121 INTENT(INOUT ) :: &
122 gd_cloud,gd_cloud2,bl_cld,bl_cld2
123 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
124 INTENT(IN ) :: tauaer1, tauaer2, tauaer3, tauaer4, &
125 waer1, waer2, waer3, waer4, &
126 gaer1, gaer2, gaer3, gaer4
127 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), &
128 INTENT(INOUT ) :: &
129 ph_macr,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2,&
130 ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, &
131 ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, &
132 ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,&
133 ph_n2o5,ph_o2,ph_n2o,ph_pooh,ph_pan,ph_mvk,ph_hyac, &
134 ph_glyald,ph_mek,ph_gly
135 REAL, INTENT(INOUT) :: dt_cld(ims:ime,kms:kme,jms:jme), &
136 af_dir(ims:ime,kms:kme,jms:jme), &
137 af_dn(ims:ime,kms:kme,jms:jme), &
138 af_up(ims:ime,kms:kme,jms:jme), &
139 par(ims:ime,kms:kme,jms:jme), &
140 erythema(ims:ime,kms:kme,jms:jme)
141 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), &
142 INTENT(IN ) :: chem
143 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , &
144 INTENT(IN ) :: &
145 z, &
146 t_phy, &
147 p_phy, &
148 dz8w, &
149 t8w,p8w,z_at_w , &
150 aerwrf , &
151 rho_phy
152 REAL, DIMENSION( ims:ime , jms:jme ) , &
153 INTENT(IN ) :: &
154 xlat, &
155 xlong
156 REAL, DIMENSION( ims:ime , jms:jme ) , &
157 INTENT(INOUT ) :: uvrad
158 REAL, INTENT(IN ) :: gmt
160 REAL, DIMENSION( ims:ime, jme:jme ), &
161 OPTIONAL,INTENT(IN ) :: o3_gfs_du
163 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
165 LOGICAL, INTENT(IN) :: haveaer
167 !---------------------------------------------------------------------
168 ! ... local variables
169 !---------------------------------------------------------------------
170 real, parameter :: boltz = 1.38065e-23 ! J/K/molecule
171 real, parameter :: xair_kg = 1.e3/28.97*6.023e23*1.e-6
172 real(dp), parameter :: Pa2hPa = .01_dp
174 integer :: astat, ierr
175 integer :: k, ki, kp1, i, j, n, wn
176 integer :: ktsm1, ktem1
177 integer :: nlev, nlyr
178 integer :: jndx
179 integer :: n_tuv_z, n_exo_z, last_n_tuv_z
180 integer :: min_ndx(2), max_ndx(2)
181 integer :: myproc
182 integer(kind=8) :: ixhour
183 integer :: minndx(2)
185 real :: xmin, gmtp, uvb_dd1, uvb_du1, uvb_dir1
186 real :: zenith, dobsi
187 real :: delz_top
188 real(kind=8) :: xtime, xhour
189 real :: max_tauaer
190 real :: Dobson(2)
191 real :: xsect(nwave)
192 real :: wrk(nwave), wrk_lam(nwave)
193 ! real :: tlev(kts-1:kte)
194 ! real :: tlyr(kts-1:kte)
195 ! real :: o33(kts-1:kte)
196 real :: rhoa(kts-1:kte)
197 ! real :: dens_air(kts-1:kte)
198 real :: aerext(kts-1:kte)
199 ! real :: qll(kts-1:kte)
200 ! real :: aircol(kts-1:kte)
201 ! real :: o3col(kts-1:kte)
202 ! real :: o2col(kts-1:kte)
203 ! real :: so2col(kts-1:kte)
204 ! real :: no2col(kts-1:kte)
205 real :: dtcld_col(kts:kte)
206 ! real :: cldfrac(kts-1:kte)
207 real :: par_col(kts:kte)
208 ! real :: tauaer300(kts-1:kte), tauaer400(kts-1:kte), &
209 ! tauaer600(kts-1:kte), tauaer999(kts-1:kte)
210 ! real :: waer300(kts-1:kte), waer400(kts-1:kte), &
211 ! waer600(kts-1:kte), waer999(kts-1:kte)
212 ! real :: gaer300(kts-1:kte), gaer400(kts-1:kte), &
213 ! gaer600(kts-1:kte), gaer999(kts-1:kte)
214 real :: zen_angle(ims:ime,jms:jme)
215 real :: z_top(its:ite,jts:jte)
216 real :: z_exo(1), dens_exo(1), temp_exo(1)
217 real :: dummy(kz)
218 real(dp) :: p_top(its:ite,jts:jte)
219 real(dp) :: o2_exo_col(its:ite,jts:jte)
220 real(dp) :: o3_exo_col(its:ite,jts:jte)
221 real(dp) :: o3_exo_col_at_grnd(its:ite,jts:jte)
223 logical :: do_alloc
224 logical :: has_aer_ra_feedback
226 real, allocatable :: tuv_z(:)
227 real, allocatable :: dtuv_z(:)
228 real, allocatable :: cldfrac(:)
229 real, allocatable :: qll(:)
230 real, allocatable :: tlev(:)
231 real, allocatable :: tlyr(:)
232 real, allocatable :: dens_air(:)
233 real, allocatable :: o33(:)
234 real, allocatable :: aircol(:)
235 real, allocatable :: o3col(:)
236 real, allocatable :: o2col(:)
237 real, allocatable :: so2col(:)
238 real, allocatable :: no2col(:)
239 real, allocatable :: tauaer300(:), tauaer400(:), tauaer600(:), tauaer999(:)
240 real, allocatable :: waer300(:), waer400(:), waer600(:), waer999(:)
241 real, allocatable :: gaer300(:), gaer400(:), gaer600(:), gaer999(:)
243 real, allocatable :: dtcld(:,:), dtaer(:,:)
244 real, allocatable :: omcld(:,:), omaer(:,:)
245 real, allocatable :: gcld(:,:), gaer(:,:)
247 real, allocatable :: rad_fld(:,:)
248 real, allocatable :: e_fld(:,:)
249 real, allocatable :: tuv_prate(:,:)
250 real, allocatable :: xsqy(:,:)
251 real, allocatable :: srb_o2_xs(:,:)
252 real, allocatable :: o3_xs(:,:), o3_xs_tpose(:,:)
253 real, allocatable :: no2_xs(:,:), no2_xs_tpose(:,:)
254 real, allocatable :: rad_fld_tpose(:,:)
255 real, allocatable :: dir_fld(:,:), dwn_fld(:,:), up_fld(:,:)
256 real, allocatable :: e_dir(:,:), e_dn(:,:), e_up(:,:)
258 REAL :: zexo_grd(2*kte)
259 CHARACTER(len=256) :: msg
261 #ifdef SW_DEBUG
262 logical :: tuv_diags
263 #endif
265 call wrf_get_myproc( myproc )
266 has_aer_ra_feedback = config_flags%aer_ra_feedback == 1
268 xtime = curr_secs/60._8 ! min since simulation start
269 ixhour = int( gmt+.01,8 ) + int( xtime/60._8,8 )
270 xhour = real( ixhour,8 )
271 xmin = 60.*gmt + real( xtime-xhour*60._8,4 )
272 gmtp = real( mod(xhour,24._8),4 )
273 gmtp = gmtp + xmin/60.
275 call calc_zenith( xlat, -xlong, &
276 julday, real(gmtp,8), zen_angle, &
277 its, ite, jts, jte, &
278 ims, ime, jms, jme )
280 #ifdef SW_DEBUG
281 if( myproc == 6 ) then
282 write(*,'(''tuv_drvr: its,jts = '',2i5)') its,jts
283 minndx(:) = minloc( zen_angle(its:ite,jts:jte) )
284 write(*,'(''tuv_drvr: min zen_angle ndx = '',2i5)') minndx(:)
285 write(*,'(''tuv_drvr: min zen_angle = '',1p,g15.8)') minval( zen_angle(its:ite,jts:jte) )
286 minndx(1) = minndx(1) + its - 1
287 minndx(2) = minndx(2) + jts - 1
288 write(*,'(''tuv_drvr: minndx = '',2i5)') minndx(:)
289 write(*,'(''tuv_drvr: min zen_angle = '',1p,g15.8)') zen_angle(minndx(1),minndx(2))
290 endif
291 #endif
294 do j = jts,jte
295 where( zen_angle(its:ite,j) == 90. )
296 zen_angle(its:ite,j) = 89.9
297 endwhere
298 end do
300 ktsm1 = kts - 1 ; ktem1 = kte - 1
302 allocate( tuv_prate(kts:kte,nj),stat=ierr )
303 if( ierr /= 0 ) then
304 call wrf_error_fatal( 'tuv_driver: failed to allocate tuv_prate' )
305 endif
307 any_daylight: &
308 if( any( zen_angle(its:ite,jts:jte) < 90. ) ) then
309 if( config_flags%has_o3_exo_coldens ) then
310 ! p_top(its:ite,jts:jte) = Pa2hPa * real( p_phy(its:ite,kte,jts:jte),dp )
311 z_exo(1) = 50.
312 call z_interp( z_air_dens_data, air_dens_data, n_air_dens_data, &
313 z_exo, dens_exo )
314 call z_interp( z_temp_data, temp_data, n_temp_data, z_exo, temp_exo )
315 p_top(its:ite,jts:jte) = temp_exo(1)*boltz*dens_exo(1)*1.e6*Pa2hPa
316 call tuv_timestep_init( dm, julday )
317 call p_interp( o2_exo_col, o3_exo_col, o3_exo_col_at_grnd, p_top, &
318 dm, its, ite, jts, jte )
319 endif
320 allocate( rad_fld(nwave,kts:kte),e_fld(kts:kte,nwave),stat=astat )
321 ierr = ierr + astat
322 allocate( dir_fld(kts:kte,nwave), dwn_fld(kts:kte,nwave), up_fld(kts:kte,nwave),stat=astat )
323 ierr = ierr + astat
324 allocate( e_dir(kts:kte,nwave), e_dn(kts:kte,nwave), e_up(kts:kte,nwave),stat=astat )
325 ierr = ierr + astat
326 if( .not. is_full_tuv ) then
327 if( any( .not. xsqy_is_zdep(:) ) ) then
328 allocate( rad_fld_tpose(kts:kte,nwave),stat=astat )
329 endif
330 elseif( any( xsqy_table(2:nj)%tpflag == 0 ) ) then
331 allocate( rad_fld_tpose(kts:kte,nwave),stat=astat )
332 endif
333 ierr = ierr + astat
334 allocate( srb_o2_xs(nwave,kts:kte) )
335 ierr = ierr + astat
336 allocate( xsqy(nwave,kts:kte) )
337 ierr = ierr + astat
338 if( ierr /= 0 ) then
339 call wrf_error_fatal( 'tuv_driver: failed to allocate rad_fld ... xsqy' )
340 endif
341 !-----------------------------------------------------------------------------
342 ! set solar distance factor
343 !-----------------------------------------------------------------------------
344 if( curjulday /= julday ) then
345 curjulday = julday
346 esfact = sundis( julday )
347 endif
348 if( .not. config_flags%scale_o3_to_grnd_exo_coldens ) then
349 if( config_flags%scale_o3_to_du_at_grnd ) then
350 dobsi = max( 0.,config_flags%du_at_grnd )
351 else
352 dobsi = 0.
353 endif
354 endif
355 endif any_daylight
357 if( ierr /= 0 ) then
358 call wrf_error_fatal( 'tuv_driver: failed to allocate module variables' )
359 endif
361 z_top(:,:) = 0.
362 do j = jts,jte
363 do i = its,ite
364 if( zen_angle(i,j) < 90. ) then
365 wrk(1) = m2km*(z(i,kte,j) - z_at_w(i,kts,j))
366 z_top(i,j) = real( ceiling(wrk(1)) )
367 delz_top = z_top(i,j) - wrk(1)
368 if( delz_top < .3 ) then
369 z_top(i,j) = z_top(i,j) + 1.
370 endif
371 endif
372 end do
373 end do
375 if( is_full_tuv ) then
376 rxn_initialized = .not. get_initialization()
377 if( .not. rxn_initialized ) then
378 do n = 1,nj
379 jndx = rxn_ndx(n)
380 if( jndx /= -1 ) then
381 call the_subs(jndx)%xsqy_sub(nwave+1,wl,wc,kz,dummy,dummy,jndx)
382 endif
383 enddo
384 call set_initialization( status=.false. )
385 endif
386 endif
388 last_n_tuv_z = -1
389 lat_loop : &
390 do j = jts,jte
391 long_loop : &
392 do i = its,ite
393 #ifdef SW_DEBUG
394 tuv_diags = i == minndx(1) .and. j == minndx(2)
395 #endif
396 do n = 1,nj
397 tuv_prate(:,n) = 0.
398 end do
399 zenith = zen_angle(i,j)
400 !---------------------------------------------------------------------
401 ! if night, skip radiative field calculation
402 !---------------------------------------------------------------------
403 has_daylight : &
404 if( zenith < 90. ) then
405 do k = kts,kte
406 rad_fld(:,k) = 0.
407 end do
408 do wn = 1,nwave
409 e_fld(:,wn) = 0.
410 end do
411 wrk(1) = m2km*(z(i,kte,j) - z_at_w(i,kts,j))
412 call setzgrid( wrk(1), n_exo_z, zexo_grd )
413 n_tuv_z = kte + n_exo_z
414 nlev = n_tuv_z - ktsm1 + 1
415 nlyr = nlev - 1
416 do_alloc = n_tuv_z /= last_n_tuv_z
417 last_n_tuv_z = n_tuv_z
418 !---------------------------------------------------------------------
419 ! allocate column variables
420 !---------------------------------------------------------------------
421 if( do_alloc ) then
422 call tuv_deallocate
423 call tuv_allocate
424 endif
426 !---------------------------------------------------------------------
427 ! column vertical grid including exo-model top
428 !---------------------------------------------------------------------
429 tuv_z(ktsm1) = 0.
430 tuv_z(kts:kte) = (z(i,kts:kte,j) - z_at_w(i,kts,j))*m2km
431 tuv_z(kte+1:kte+n_exo_z) = zexo_grd(1:n_exo_z)
432 ! tuv_z(kte+1:kte+n_exo_z) = real( (/ (n,n=k,k+n_exo_z-1) /) )
433 ! tuv_z(kte+n_exo_z+1:n_tuv_z) = real( (/ (25+5*n,n=1,5) /) )
434 !---------------------------------------------------------------------
435 ! cloud fraction
436 !---------------------------------------------------------------------
437 cldfrac(:) = 0.
438 if( config_flags%cld_od_opt > 1 ) then
439 if( config_flags%pht_cldfrc_opt == 1 ) then
440 call cldfrac_binary( cldfrac(kts:kte), &
441 moist(i,kts:kte,j,p_qc), moist(i,kts:kte,j,p_qi), &
442 moist(i,kts:kte,j,p_qs), kts, kte )
443 elseif( config_flags%pht_cldfrc_opt == 2 ) then
444 if( config_flags%cu_diag == 1 ) then
445 if ( config_flags%phot_blcld .and. config_flags%icloud_bl == 1 ) then
446 call cldfrac_fractional( cldfrac(kts:kte), &
447 moist(i,kts:kte,j,p_qv), &
448 moist(i,kts:kte,j,p_qc) + gd_cloud(i,kts:kte,j) + &
449 bl_cld(i,kts:kte,j), &
450 moist(i,kts:kte,j,p_qi) + gd_cloud2(i,kts:kte,j) + &
451 bl_cld2(i,kts:kte,j), &
452 moist(i,kts:kte,j,p_qs), &
453 p_phy(i,kts:kte,j), t_phy(i,kts:kte,j), kts, kte )
454 else
455 call cldfrac_fractional( cldfrac(kts:kte), &
456 moist(i,kts:kte,j,p_qv), &
457 moist(i,kts:kte,j,p_qc) + gd_cloud(i,kts:kte,j), &
458 moist(i,kts:kte,j,p_qi) + gd_cloud2(i,kts:kte,j), &
459 moist(i,kts:kte,j,p_qs), &
460 p_phy(i,kts:kte,j), t_phy(i,kts:kte,j), kts, kte )
461 endif
462 else
463 call cldfrac_fractional( cldfrac(kts:kte), &
464 moist(i,kts:kte,j,p_qv), moist(i,kts:kte,j,p_qc), &
465 moist(i,kts:kte,j,p_qi), moist(i,kts:kte,j,p_qs), &
466 p_phy(i,kts:kte,j), t_phy(i,kts:kte,j), kts, kte )
467 endif
468 endif
469 endif
470 !---------------------------------------------------------------------
471 ! aerosols
472 !---------------------------------------------------------------------
473 tauaer300(:) = 0.
474 tauaer400(:) = 0.
475 tauaer600(:) = 0.
476 tauaer999(:) = 0.
477 waer300(:) = 0.
478 waer400(:) = 0.
479 waer600(:) = 0.
480 waer999(:) = 0.
481 gaer300(:) = 0.
482 gaer400(:) = 0.
483 gaer600(:) = 0.
484 gaer999(:) = 0.
485 if( has_aer_ra_feedback ) then
486 tauaer300(kts:kte) = tauaer1(i,kts:kte,j)
487 tauaer400(kts:kte) = tauaer2(i,kts:kte,j)
488 tauaer600(kts:kte) = tauaer3(i,kts:kte,j)
489 tauaer999(kts:kte) = tauaer4(i,kts:kte,j)
490 waer300(kts:kte) = waer1(i,kts:kte,j)
491 waer400(kts:kte) = waer2(i,kts:kte,j)
492 waer600(kts:kte) = waer3(i,kts:kte,j)
493 waer999(kts:kte) = waer4(i,kts:kte,j)
494 gaer300(kts:kte) = gaer1(i,kts:kte,j)
495 gaer400(kts:kte) = gaer2(i,kts:kte,j)
496 gaer600(kts:kte) = gaer3(i,kts:kte,j)
497 gaer999(kts:kte) = gaer4(i,kts:kte,j)
498 tauaer300(ktsm1) = tauaer300(kts)
499 tauaer400(ktsm1) = tauaer400(kts)
500 tauaer600(ktsm1) = tauaer600(kts)
501 tauaer999(ktsm1) = tauaer999(kts)
502 waer300(ktsm1) = waer300(kts)
503 waer400(ktsm1) = waer400(kts)
504 waer600(ktsm1) = waer600(kts)
505 waer999(ktsm1) = waer999(kts)
506 gaer300(ktsm1) = gaer300(kts)
507 gaer400(ktsm1) = gaer400(kts)
508 gaer600(ktsm1) = gaer600(kts)
509 gaer999(ktsm1) = gaer999(kts)
510 endif
512 tlev(ktsm1) = t8w(i,kts,j)
513 tlev(kts:kte) = t_phy(i,kts:kte,j)
514 call z_interp( z_temp_data, temp_data, n_temp_data, &
515 tuv_z(kte+1:n_tuv_z), tlev(kte+1:n_tuv_z) )
516 tlyr(ktsm1:n_tuv_z-1) = .5*(tlev(ktsm1:n_tuv_z-1) + tlev(kts:n_tuv_z))
518 rhoa(ktsm1) = p8w(i,kts,j)/(t8w(i,kts,j)*r_d)
519 rhoa(kts:kte) = rho_phy(i,kts:kte,j)
520 dens_air(ktsm1:kte) = xair_kg*rhoa(ktsm1:kte) ! air num.(molecules/cm^3)
521 call z_interp( z_air_dens_data, air_dens_data, n_air_dens_data, &
522 tuv_z(kte+1:n_tuv_z), dens_air(kte+1:n_tuv_z) )
524 o33(kts:kte) = ppm2vmr*chem(i,kts:kte,j,p_o3)*dens_air(kts:kte)
525 o33(ktsm1) = o33(kts)
526 call z_interp( z_o3_data, o3_data, n_o3_data, &
527 tuv_z(kte+1:n_tuv_z), o33(kte+1:n_tuv_z) )
529 qll(:) = 0.
530 qll(kts:kte) = moist(i,kts:kte,j,p_qc) + moist(i,kts:kte,j,p_qi)
531 if( config_flags%cu_diag == 1 ) then
532 if ( config_flags%phot_blcld .and. config_flags%icloud_bl == 1 ) then
533 qll(kts:kte) = qll(kts:kte) + gd_cloud(i,kts:kte,j) + &
534 gd_cloud2(i,kts:kte,j) + &
535 bl_cld(i,kts:kte,j) + &
536 bl_cld2(i,kts:kte,j)
537 else
538 qll(kts:kte) = qll(kts:kte) + gd_cloud(i,kts:kte,j) + gd_cloud2(i,kts:kte,j)
539 endif
540 endif
541 qll(kts:kte) = 1.e3*rhoa(kts:kte)*qll(kts:kte)
542 where( qll(kts:kte) < 1.e-5 )
543 qll(kts:kte) = 0.
544 endwhere
546 if( haveaer .and. ktau > 1 )then
547 aerext(ktsm1) = aerext(kts)
548 else
549 aerext(ktsm1) = aerwrf(i,kts,j)
550 endif
552 if( .not. is_full_tuv ) then
553 call xs_int( o3_xs, tlyr, o3_xs_tab )
554 call xs_int( no2_xs, tlyr, no2_xs_tab )
555 else
556 call o3xs( nlyr,tlyr,nwave+1,wl,o3_xs_tpose )
557 call no2xs_jpl06a( nlyr,tlyr,nwave+1,wl,no2_xs_tpose )
558 o3_xs = transpose( o3_xs_tpose )
559 no2_xs = transpose( no2_xs_tpose )
560 endif
562 dtuv_z(ktsm1:n_tuv_z-1) = tuv_z(kts:n_tuv_z) - tuv_z(ktsm1:n_tuv_z-1)
563 aircol(ktsm1:n_tuv_z-1) = &
564 km2cm*dtuv_z(ktsm1:n_tuv_z-1) &
565 *real(sqrt(real(dens_air(kts:n_tuv_z),kind=8)*real(dens_air(ktsm1:n_tuv_z-1),kind=8)),kind=4)
566 ! km2cm*dtuv_z(ktsm1:n_tuv_z-1)*(dens_air(kts:n_tuv_z) - dens_air(ktsm1:n_tuv_z-1)) &
567 ! / log( dens_air(kts:n_tuv_z)/dens_air(ktsm1:n_tuv_z-1) )
568 o3col(ktsm1:n_tuv_z-1) = conv1*dtuv_z(ktsm1:n_tuv_z-1)*(o33(kts:n_tuv_z) + o33(ktsm1:n_tuv_z-1))
569 o2col(ktsm1:n_tuv_z-1) = o2vmr*aircol(ktsm1:n_tuv_z-1)
570 if( config_flags%has_o3_exo_coldens ) then
571 o3col(n_tuv_z-1) = o3col(n_tuv_z-1) + real( o3_exo_col(i,j),4 )
572 o2col(n_tuv_z-1) = o2col(n_tuv_z-1) + real( o2_exo_col(i,j),4 )
573 aircol(n_tuv_z-1) = aircol(n_tuv_z-1) + o2col(n_tuv_z-1)/o2vmr
574 endif
575 if( config_flags%scale_o3_to_grnd_exo_coldens ) then
576 dobsi = real( o3_exo_col_at_grnd(i,j),4 )/2.687e16
577 endif
578 if (config_flags%scale_o3_to_gfs_tot ) then
579 dobsi = o3_gfs_du(i,j)
580 endif
582 so2col(:) = 0.
583 if( has_so2 ) then
584 so2col(ktsm1) = conv2*dtuv_z(ktsm1) &
585 *chem(i,kts,j,p_so2)*(dens_air(ktsm1) + dens_air(kts))
586 so2col(kts:ktem1) = conv2*dtuv_z(kts:ktem1) &
587 *(chem(i,kts:ktem1,j,p_so2)*dens_air(kts:ktem1) &
588 + chem(i,kts+1:kte,j,p_so2)*dens_air(kts+1:kte))
589 endif
590 no2col(:) = 0.
591 if( has_no2 ) then
592 no2col(ktsm1) = conv2*dtuv_z(ktsm1) &
593 *chem(i,kts,j,p_no2)*(dens_air(ktsm1) + dens_air(kts))
594 no2col(kts:ktem1) = conv2*dtuv_z(kts:ktem1) &
595 *(chem(i,kts:ktem1,j,p_no2)*dens_air(kts:ktem1) &
596 + chem(i,kts+1:kte,j,p_no2)*dens_air(kts+1:kte))
597 endif
599 call tuv_radfld( nlambda_start, config_flags%cld_od_opt, cldfrac, &
600 nlyr, nwave, zenith, tuv_z, albedo, &
601 aircol, o2col, o3col, so2col, no2col, &
602 tauaer300, tauaer400, tauaer600, tauaer999, &
603 waer300, waer400, waer600, waer999, &
604 gaer300, gaer400, gaer600, gaer999, &
605 dtaer, omaer, gaer, dtcld, omcld, gcld, &
606 has_aer_ra_feedback, &
607 qll, dobsi, o3_xs, no2_xs, o2_xs, &
608 so2_xs, wl(1), wc, tlev, srb_o2_xs, rad_fld, e_fld, &
609 e_dir, e_dn, e_up, &
610 dir_fld, dwn_fld, up_fld, dtcld_col )
612 do k = kts,kte
613 af_dir(i,k,j) = dot_product( dir_fld(k,nlambda_af_start:nlambda_af_end),dw(nlambda_af_start:nlambda_af_end) )
614 af_dn(i,k,j) = dot_product( dwn_fld(k,nlambda_af_start:nlambda_af_end),dw(nlambda_af_start:nlambda_af_end) )
615 af_up(i,k,j) = dot_product( up_fld(k,nlambda_af_start:nlambda_af_end),dw(nlambda_af_start:nlambda_af_end) )
616 end do
618 dt_cld(i,kts:kte,j) = dtcld_col(kts:kte)
619 wrk(nlambda_start:nwave) = dw(nlambda_start:nwave)*etfl(nlambda_start:nwave)
620 wrk_lam(nlambda_start:nwave) = wrk(nlambda_start:nwave)*par_wght(nlambda_start:nwave)
621 par(i,kts:kte,j) = matmul( e_fld(kts:kte,nlambda_start:nwave), wrk_lam(nlambda_start:nwave) )
622 wrk_lam(nlambda_start:nwave) = wrk(nlambda_start:nwave)*ery_wght(nlambda_start:nwave)
623 erythema(i,kts:kte,j) = matmul( e_fld(kts:kte,nlambda_start:nwave), wrk_lam(nlambda_start:nwave) )
625 if( .not. is_full_tuv ) then
626 if( any( .not. xsqy_is_zdep(:) ) ) then
627 rad_fld_tpose = transpose( rad_fld )
628 endif
629 elseif( any( xsqy_table(1:nj)%tpflag == 0 ) ) then
630 rad_fld_tpose = transpose( rad_fld )
631 endif
633 #ifdef SW_DEBUG
634 if( tuv_diags ) then
635 Dobson(1) = sum( o3col(ktsm1:ktem1) )/2.687e16
636 Dobson(2) = o3col(kte)/2.687e16
637 write(*,'(''tuv_drvr: o3col = '',1p,2g15.8)') Dobson(:)
638 open(unit=33,file='wrfchm_inp')
639 write(33,*) nlev
640 write(33,*) kte
641 write(33,*) tuv_z(ktsm1:n_tuv_z)
642 write(33,*) tlev(ktsm1:n_tuv_z)
643 write(33,*) o3col(ktsm1:n_tuv_z-1)
644 write(33,*) so2col(ktsm1:n_tuv_z-1)
645 write(33,*) no2col(ktsm1:n_tuv_z-1)
646 write(33,*) dens_air(ktsm1:n_tuv_z)
647 write(33,*) aircol(ktsm1:n_tuv_z-1)
648 do wn = 1,nwave
649 write(33,*) dtaer(ktsm1:n_tuv_z-1,wn)
650 end do
651 do wn = 1,nwave
652 write(33,*) omaer(ktsm1:n_tuv_z-1,wn)
653 end do
654 do wn = 1,nwave
655 write(33,*) gaer(ktsm1:n_tuv_z-1,wn)
656 end do
657 do wn = 1,nwave
658 write(33,*) dtcld(ktsm1:n_tuv_z-1,wn)
659 end do
660 do wn = 1,nwave
661 write(33,*) omcld(ktsm1:n_tuv_z-1,wn)
662 end do
663 do wn = 1,nwave
664 write(33,*) gcld(ktsm1:n_tuv_z-1,wn)
665 end do
666 write(33,*) zen_angle(i,j)
667 write(33,*) esfact
668 write(33,*) dobsi
669 close( 33 )
671 open(unit=33,file='WRF-TUV.flx.out' )
672 do n = nlambda_start,nwave
673 write(33,*) dir_fld(kts:kte,n)
674 end do
675 write(33,*) ' '
676 do n = nlambda_start,nwave
677 write(33,*) dwn_fld(kts:kte,n)
678 end do
679 write(33,*) ' '
680 do n = nlambda_start,nwave
681 write(33,*) up_fld(kts:kte,n)
682 end do
683 write(33,*) ' '
684 do n = nlambda_start,nwave
685 write(33,*) rad_fld_tpose(kts:kte,n)
686 end do
687 close( 33 )
688 endif
689 #endif
691 rate_loop: &
692 do n = 1,nj
693 !---------------------------------------------------------------------
694 ! set cross-section x quantum yields
695 !---------------------------------------------------------------------
696 if( n /= j_o2_ndx ) then
697 if( .not. is_full_tuv ) then
698 if( xsqy_is_zdep(n) ) then
699 call xsqy_int( n, xsqy, tlev(kts:kte), dens_air(kts:kte) )
700 endif
701 else
702 jndx = rxn_ndx(n)
703 if( jndx /= -1 ) then
704 if( xsqy_table(jndx)%tpflag /= 0 ) then
705 call the_subs(jndx)%xsqy_sub(nwave+1,wl,wc,nlev,tlev,dens_air,jndx)
706 endif
707 endif
708 endif
709 elseif( .not. is_full_tuv ) then
710 call sjo2( kte, nwave, srb_o2_xs, xsqy )
711 endif
712 !---------------------------------------------------------------------
713 ! compute tuv photorates
714 !---------------------------------------------------------------------
715 if( .not. is_full_tuv ) then
716 if( xsqy_is_zdep(n) ) then
717 do k = kts,kte
718 xsect(nlambda_start:nwave) = xsqy(nlambda_start:nwave,k)*w_fac(nlambda_start:nwave)*esfact
719 tuv_prate(k,n) = m2s*dot_product( rad_fld(nlambda_start:nwave,k),xsect(nlambda_start:nwave) )
720 end do
721 else
722 xsect(nlambda_start:nwave) = xsqy_tab(nlambda_start:nwave,1,1,n)*w_fac(nlambda_start:nwave)*esfact
723 tuv_prate(:,n) = m2s*matmul( rad_fld_tpose(:,nlambda_start:nwave),xsect(nlambda_start:nwave) )
724 endif
725 else
726 if( n /= j_o2_ndx ) then
727 if( xsqy_table(jndx)%tpflag > 0 ) then
728 do k = kts,kte
729 xsect(nlambda_start:nwave) = xsqy_table(jndx)%sq(k+1,nlambda_start:nwave)*w_fac(nlambda_start:nwave)*esfact
730 tuv_prate(k,n) = m2s*dot_product( rad_fld(nlambda_start:nwave,k),xsect(nlambda_start:nwave) )
731 end do
732 else
733 xsect(nlambda_start:nwave) = xsqy_table(jndx)%sq(nlambda_start:nwave,1)*w_fac(nlambda_start:nwave)*esfact
734 tuv_prate(:,n) = m2s*matmul( rad_fld_tpose(:,nlambda_start:nwave),xsect(nlambda_start:nwave) )
735 endif
736 else
737 do k = kts,kte
738 xsect(nlambda_start:nwave) = srb_o2_xs(nlambda_start:nwave,k)*w_fac(nlambda_start:nwave)*esfact
739 tuv_prate(k,n) = m2s*dot_product( rad_fld(nlambda_start:nwave,k),xsect(nlambda_start:nwave) )
740 end do
741 endif
742 endif
743 end do rate_loop
745 #ifdef SW_DEBUG
746 if( myproc == 6 ) then
747 if( tuv_diags ) then
748 open(unit=33,file='WRF-TUV.j.out' )
749 do n = 1,nj
750 select case( n )
751 case( 2:4,7:10,12:13 )
752 do k = kts,kte,5
753 write(33,'(1p,5g15.7)') tuv_prate(k:min(k+4,kte),n)/m2s
754 end do
755 if( n < 13 ) then
756 write(33,*) ' '
757 endif
758 end select
759 end do
760 close(33)
761 call wrf_abort
762 endif
763 endif
764 #endif
765 endif has_daylight
766 #if ( WRF_KPP == 1 )
767 #include "tuv2wrf_jvals.inc"
768 #endif
770 end do long_loop
771 end do lat_loop
773 #ifdef SW_DEBUG
774 if( any( zen_angle(its:ite,jts:jte) < 90. ) ) then
775 min_ndx(:) = minloc( z_top,mask=z_top>0. )
776 min_ndx(:) = (/ min_ndx(1) + its - 1,min_ndx(2) + jts - 1 /)
777 max_ndx(:) = maxloc( z_top,mask=z_top>0. )
778 max_ndx(:) = (/ max_ndx(1) + its - 1,max_ndx(2) + jts - 1 /)
779 write(msg,'(''tuv_driver: z_top min indices = '',2i6)') min_ndx(:)
780 call wrf_debug( 0,trim(msg) )
781 write(msg,'(''tuv_driver: z_top max indices = '',2i6)') max_ndx(:)
782 call wrf_debug( 0,trim(msg) )
783 write(msg,'(''tuv_driver: z_top min,max = '',1p2g15.7)') z_top(min_ndx(1),min_ndx(2)),z_top(max_ndx(1),max_ndx(2))
784 call wrf_debug( 0,trim(msg) )
785 endif
786 #endif
788 call tuv_deallocate
790 ierr = 0
791 if( allocated( rad_fld ) ) then
792 deallocate( rad_fld,stat=astat )
793 ierr = ierr + astat
794 deallocate( dir_fld, dwn_fld, up_fld,stat=astat )
795 ierr = ierr + astat
796 deallocate( e_dir, e_dn, e_up,stat=astat )
797 ierr = ierr + astat
798 endif
799 if( allocated( rad_fld_tpose ) ) then
800 deallocate( rad_fld_tpose,stat=astat )
801 ierr = ierr + astat
802 endif
803 if( allocated( e_fld ) ) then
804 deallocate( e_fld,stat=astat )
805 ierr = ierr + astat
806 endif
807 if( allocated( xsqy ) ) then
808 deallocate( xsqy,stat=astat )
809 ierr = ierr + astat
810 endif
811 if( allocated( srb_o2_xs ) ) then
812 deallocate( srb_o2_xs,stat=astat )
813 ierr = ierr + astat
814 endif
816 if( allocated( tuv_prate ) ) then
817 deallocate( tuv_prate,stat=astat )
818 ierr = ierr + astat
819 endif
821 if( ierr /= 0 ) then
822 call wrf_error_fatal( 'tuv_driver: failed to deallocate local variables' )
823 endif
825 CONTAINS
827 subroutine setzgrid( ztop, nexo, zexo_grd )
828 !---------------------------------------------------------------------
829 ! set the vertical grid above model top up to 50 km for photorate
830 ! calculation
831 !---------------------------------------------------------------------
833 integer, intent(out) :: nexo
834 real, intent(in) :: ztop
835 real, intent(out) :: zexo_grd(:)
837 real, parameter :: ztrop = 25.
838 real, parameter :: dztrop = 1.
839 real, parameter :: zlid = 50.
840 real, parameter :: dzlid = 5.
842 real :: z
843 real :: zBase
845 nexo = 0
846 if( ztop >= zlid ) then
847 return
848 elseif( ztop <= ztrop ) then
849 zBase = ztrop
850 else
851 zBase = dzlid*real( int(ztop/dzlid) )
852 endif
854 z = int( ztop )
855 do while( z < ztrop )
856 z = z + dztrop
857 nexo = nexo + 1
858 zexo_grd(nexo) = z
859 enddo
861 z = zBase
862 do while( z < zlid )
863 z = z + dzlid
864 nexo = nexo + 1
865 zexo_grd(nexo) = z
866 end do
868 end subroutine setzgrid
870 subroutine tuv_allocate
871 !---------------------------------------------------------------------
872 ! allocate column variables
873 !---------------------------------------------------------------------
875 allocate( tuv_z(ktsm1:n_tuv_z),dtuv_z(ktsm1:n_tuv_z-1),stat=ierr )
876 allocate( cldfrac(ktsm1:n_tuv_z),stat=astat )
877 ierr = ierr + astat
878 allocate( qll(ktsm1:n_tuv_z),stat=astat )
879 ierr = ierr + astat
880 allocate( tlev(ktsm1:n_tuv_z),tlyr(ktsm1:n_tuv_z-1),stat=astat )
881 ierr = ierr + astat
882 allocate( dens_air(ktsm1:n_tuv_z),stat=astat )
883 ierr = ierr + astat
884 allocate( o33(ktsm1:n_tuv_z),stat=astat )
885 ierr = ierr + astat
886 allocate( o3_xs(nwave,ktsm1:n_tuv_z-1),stat=astat )
887 ierr = ierr + astat
888 if( is_full_tuv ) then
889 allocate( o3_xs_tpose(ktsm1:n_tuv_z-1,nwave),stat=astat )
890 ierr = ierr + astat
891 allocate( no2_xs_tpose(ktsm1:n_tuv_z-1,nwave),stat=astat )
892 ierr = ierr + astat
893 endif
894 allocate( no2_xs(nwave,ktsm1:n_tuv_z-1),stat=astat )
895 ierr = ierr + astat
896 allocate( aircol(ktsm1:n_tuv_z-1),stat=astat )
897 ierr = ierr + astat
898 allocate( o2col(ktsm1:n_tuv_z-1),stat=astat )
899 ierr = ierr + astat
900 allocate( o3col(ktsm1:n_tuv_z-1),stat=astat )
901 ierr = ierr + astat
902 allocate( no2col(ktsm1:n_tuv_z-1),stat=astat )
903 ierr = ierr + astat
904 allocate( so2col(ktsm1:n_tuv_z-1),stat=astat )
905 ierr = ierr + astat
906 allocate( tauaer300(ktsm1:n_tuv_z-1),tauaer400(ktsm1:n_tuv_z-1), &
907 tauaer600(ktsm1:n_tuv_z-1),tauaer999(ktsm1:n_tuv_z-1),stat=astat )
908 ierr = ierr + astat
909 allocate( waer300(ktsm1:n_tuv_z-1),waer400(ktsm1:n_tuv_z-1), &
910 waer600(ktsm1:n_tuv_z-1),waer999(ktsm1:n_tuv_z-1),stat=astat )
911 ierr = ierr + astat
912 allocate( gaer300(ktsm1:n_tuv_z-1),gaer400(ktsm1:n_tuv_z-1), &
913 gaer600(ktsm1:n_tuv_z-1),gaer999(ktsm1:n_tuv_z-1),stat=astat )
914 ierr = ierr + astat
915 allocate( dtaer(ktsm1:n_tuv_z-1,nwave),omaer(ktsm1:n_tuv_z-1,nwave), &
916 gaer(ktsm1:n_tuv_z-1,nwave),stat=astat )
917 ierr = ierr + astat
918 allocate( dtcld(ktsm1:n_tuv_z-1,nwave),omcld(ktsm1:n_tuv_z-1,nwave), &
919 gcld(ktsm1:n_tuv_z-1,nwave),stat=astat )
920 ierr = ierr + astat
922 if( ierr /= 0 ) then
923 call wrf_error_fatal( 'tuv_driver: failed to allocate column variables' )
924 endif
926 end subroutine tuv_allocate
928 subroutine tuv_deallocate
929 !---------------------------------------------------------------------
930 ! deallocate column variables
931 !---------------------------------------------------------------------
933 integer :: astat, ierr
935 ierr = 0
936 if( allocated( tuv_z ) ) then
937 deallocate( tuv_z,stat=astat )
938 ierr = ierr + astat
939 endif
940 if( allocated( dtuv_z ) ) then
941 deallocate( dtuv_z,stat=astat )
942 ierr = ierr + astat
943 endif
944 if( allocated( cldfrac ) ) then
945 deallocate( cldfrac,stat=astat )
946 ierr = ierr + astat
947 endif
948 if( allocated( qll ) ) then
949 deallocate( qll,stat=astat )
950 ierr = ierr + astat
951 endif
952 if( allocated( tlev ) ) then
953 deallocate( tlev,stat=astat )
954 ierr = ierr + astat
955 endif
956 if( allocated( tlyr ) ) then
957 deallocate( tlyr,stat=astat )
958 ierr = ierr + astat
959 endif
960 if( allocated( dens_air ) ) then
961 deallocate( dens_air,stat=astat )
962 ierr = ierr + astat
963 endif
964 if( allocated( o33 ) ) then
965 deallocate( o33,stat=astat )
966 ierr = ierr + astat
967 endif
968 if( allocated( o3_xs ) ) then
969 deallocate( o3_xs,stat=astat )
970 ierr = ierr + astat
971 endif
972 if( allocated( o3_xs_tpose ) ) then
973 deallocate( o3_xs_tpose,stat=astat )
974 ierr = ierr + astat
975 endif
976 if( allocated( no2_xs ) ) then
977 deallocate( no2_xs,stat=astat )
978 ierr = ierr + astat
979 endif
980 if( allocated( no2_xs_tpose ) ) then
981 deallocate( no2_xs_tpose,stat=astat )
982 ierr = ierr + astat
983 endif
984 if( allocated( aircol ) ) then
985 deallocate( aircol,stat=astat )
986 ierr = ierr + astat
987 endif
988 if( allocated( o2col ) ) then
989 deallocate( o2col,stat=astat )
990 ierr = ierr + astat
991 endif
992 if( allocated( o3col ) ) then
993 deallocate( o3col,stat=astat )
994 ierr = ierr + astat
995 endif
996 if( allocated( no2col ) ) then
997 deallocate( no2col,stat=astat )
998 ierr = ierr + astat
999 endif
1000 if( allocated( so2col ) ) then
1001 deallocate( so2col,stat=astat )
1002 ierr = ierr + astat
1003 endif
1004 if( allocated( tauaer300 ) ) then
1005 deallocate( tauaer300, waer300, gaer300,stat=astat )
1006 ierr = ierr + astat
1007 endif
1008 if( allocated( tauaer400 ) ) then
1009 deallocate( tauaer400, waer400, gaer400,stat=astat )
1010 ierr = ierr + astat
1011 endif
1012 if( allocated( tauaer600 ) ) then
1013 deallocate( tauaer600, waer600, gaer600,stat=astat )
1014 ierr = ierr + astat
1015 endif
1016 if( allocated( tauaer999 ) ) then
1017 deallocate( tauaer999, waer999, gaer999,stat=astat )
1018 ierr = ierr + astat
1019 endif
1020 if( allocated( dtcld ) ) then
1021 deallocate( dtcld, omcld, gcld,stat=astat )
1022 ierr = ierr + astat
1023 endif
1024 if( allocated( dtaer ) ) then
1025 deallocate( dtaer, omaer, gaer,stat=astat )
1026 ierr = ierr + astat
1027 endif
1029 if( ierr /= 0 ) then
1030 call wrf_error_fatal( 'tuv_deallocate: failed to deallocate all variables' )
1031 endif
1033 end subroutine tuv_deallocate
1035 end subroutine tuv_driver
1037 subroutine tuv_init( &
1038 domain, config_flags, z_at_w, aerwrf, g, &
1039 af_lambda_start, af_lambda_end, start_lambda, &
1040 ids,ide, jds,jde, kds,kde, &
1041 ims,ime, jms,jme, kms,kme, &
1042 its,ite, jts,jte, kts,kte )
1044 USE module_configure
1045 USE params_mod, only : lambda_cutoff
1046 USE srb, only : init_srb
1047 USE module_state_description, only : p_so2, p_no2, param_first_scalar
1048 USE module_rxn, only : rxn_init, xsqy_table => xsqy_tab, npht_tab
1049 USE module_rxn, only : get_initialization
1051 ! .. Scalar Arguments ..
1052 INTEGER, INTENT (IN) :: domain
1053 INTEGER, INTENT (IN) :: ide, ids, ime, ims, ite, its, jde, jds, &
1054 jme, jms, jte, jts, kde, kds, kme, kms, kte, kts
1055 REAL, INTENT (IN) :: g
1056 REAL, INTENT (IN) :: af_lambda_start, af_lambda_end, start_lambda
1057 ! .. Array Arguments ..
1058 REAL, INTENT (INOUT) :: aerwrf(ims:ime,kms:kme,jms:jme)
1059 REAL, INTENT (IN) :: z_at_w(ims:ime,kms:kme,jms:jme)
1061 TYPE(grid_config_rec_type), INTENT(IN) :: config_flags
1063 ! .. Local Scalars ..
1064 INTEGER :: astat
1065 INTEGER :: i, j, k, n, wn
1066 ! .. Local Arrays ..
1067 CHARACTER(len=256) :: msg
1069 #ifndef NETCDF
1070 call wrf_error_fatal( 'tuv_init: requires netcdf' )
1071 #endif
1073 DO j = jts, min(jte,jde-1)
1074 DO k = kts, min(kte,kde-1)
1075 DO i = its, min(ite,ide-1)
1076 aerwrf(i,k,j) = 0.
1077 END DO
1078 END DO
1079 END DO
1082 is_initialized: &
1083 if( .not. tuv_is_initialized ) then
1084 has_exo_coldens = config_flags%has_o3_exo_coldens .or. config_flags%scale_o3_to_grnd_exo_coldens
1085 is_full_tuv = config_flags%is_full_tuv
1086 #if ( WRF_KPP == 1 )
1087 #include "tuvdef_jvals.inc"
1088 #endif
1089 call get_xsqy_tab
1090 if( .not. is_full_tuv ) then
1091 allocate( xsqy_is_zdep(nj),stat=astat )
1092 if( astat /= 0 ) then
1093 call wrf_error_fatal( 'tuv_init: failed to allocate xsqy_is_zdep' )
1094 endif
1095 xsqy_is_zdep(:) = .false.
1096 if( j_o2_ndx > 0 ) then
1097 xsqy_is_zdep(j_o2_ndx) = .true.
1098 endif
1099 do n = 1,nj
1100 if( n /= j_o2_ndx ) then
1101 t_loop : do j = 1,nconc
1102 do i = 1,ntemp
1103 if( any( xsqy_tab(:,i,j,n) /= xsqy_tab(:,1,1,n) ) ) then
1104 xsqy_is_zdep(n) = .true.
1105 exit t_loop
1106 endif
1107 end do
1108 end do t_loop
1109 endif
1110 end do
1111 endif
1112 has_so2 = p_so2 >= param_first_scalar
1113 has_no2 = p_no2 >= param_first_scalar
1114 call init_srb
1115 !---------------------------------------------------------------------
1116 ! ... locate starting wave bin
1117 !---------------------------------------------------------------------
1118 lambda_cutoff = start_lambda
1119 do nlambda_start = 1,nwave
1120 if( wc(nlambda_start) >= lambda_cutoff ) then
1121 exit
1122 endif
1123 end do
1124 if( nlambda_start > nwave ) then
1125 write(msg,'(''tuv_init: '',1pg15.7,'' is not in photo wavelength interval ('',g15.7,'','',g15.7,'')'')') &
1126 lambda_cutoff,wl(1),wl(nwave+1)
1127 call wrf_error_fatal( trim(msg) )
1128 endif
1129 !---------------------------------------------------------------------
1130 ! ... locate af output wave bins
1131 !---------------------------------------------------------------------
1132 do nlambda_af_start = nlambda_start,nwave
1133 if( wl(nlambda_af_start) <= af_lambda_start .and. &
1134 wl(nlambda_af_start+1) >= af_lambda_start ) then
1135 exit
1136 endif
1137 end do
1138 do nlambda_af_end = nlambda_start,nwave
1139 if( wl(nlambda_af_end) <= af_lambda_end .and. &
1140 wl(nlambda_af_end+1) >= af_lambda_end ) then
1141 exit
1142 endif
1143 end do
1144 allocate( par_wght(nwave), ery_wght(nwave), stat=astat )
1145 if( astat /= 0 ) then
1146 call wrf_error_fatal( 'tuv_init: failed to allocate par_wght,ery_wght' )
1147 endif
1148 !---------------------------------------------------------------------
1149 ! ... set PAR weight
1150 !---------------------------------------------------------------------
1151 where (wc(:) > 400. .AND. wc(:) < 700.)
1152 par_wght(:) = 8.36e-3 * wc(:)
1153 elsewhere
1154 par_wght(:) = 0.
1155 end where
1156 !---------------------------------------------------------------------
1157 ! ... set erythema weight
1158 !---------------------------------------------------------------------
1159 call fery( nwave, wc, ery_wght )
1160 !---------------------------------------------------------------------
1161 ! ... set surface albedo
1162 !---------------------------------------------------------------------
1163 DO wn = 1,nwave
1164 IF (wl(wn)<400.) then
1165 albedo(wn) = 0.05
1166 elseIF (wl(wn)>=400. .AND. wl(wn)<450.) then
1167 albedo(wn) = 0.06
1168 elseIF (wl(wn)>=450. .AND. wl(wn)<500.) then
1169 albedo(wn) = 0.08
1170 elseIF (wl(wn)>=500. .AND. wl(wn)<550.) then
1171 albedo(wn) = 0.10
1172 elseIF (wl(wn)>=550. .AND. wl(wn)<600.) then
1173 albedo(wn) = 0.11
1174 elseIF (wl(wn)>=600. .AND. wl(wn)<640.) then
1175 albedo(wn) = 0.12
1176 elseIF (wl(wn)>=640. .AND. wl(wn)<660.) then
1177 albedo(wn) = 0.135
1178 elseIF (wl(wn)>=660.) then
1179 albedo(wn) = 0.15
1180 endIF
1181 END DO
1182 !---------------------------------------------------------------------
1183 ! ... if full TUV then initialize
1184 !---------------------------------------------------------------------
1185 if( is_full_tuv ) then
1186 call rxn_init( nwave+1,wl )
1187 allocate( rxn_ndx(nj),stat=astat )
1188 if( astat /= 0 ) then
1189 call wrf_error_fatal( 'tuv_init: failed to allocate rxn_ndx' )
1190 endif
1191 rxn_ndx(1:nj) = -1
1192 do j = 1,nj
1193 if( j /= j_o2_ndx ) then
1194 do n = 2,npht_tab
1195 if( trim(xsqy_table(n)%wrf_label) == trim(tuv_jname(j)) ) then
1196 rxn_ndx(j) = n
1197 exit
1198 endif
1199 enddo
1200 endif
1201 enddo
1202 rxn_initialized = .not. get_initialization()
1203 endif
1204 tuv_is_initialized = .true.
1205 endif is_initialized
1207 !---------------------------------------------------------------------
1208 ! ... get exo model column densities
1209 !---------------------------------------------------------------------
1210 if( has_exo_coldens ) then
1211 call get_exo_coldens( domain, config_flags%exo_coldens_inname, &
1212 ids,ide, jds,jde, kds,kde, &
1213 ims,ime, jms,jme, kms,kme, &
1214 its,ite, jts,jte, kts,kte )
1215 endif
1217 end subroutine tuv_init
1219 subroutine get_xsqy_tab
1220 !---------------------------------------------------------------------
1221 ! ... read in the cross section,quantum yield tables
1222 !---------------------------------------------------------------------
1224 use params_mod, only : hc
1225 use srb, only : ila, isrb
1226 use srb, only : nchebev_term, nchebev_wave
1227 use srb, only : chebev_ac, chebev_bc
1229 !---------------------------------------------------------------------
1230 ! ... local variables
1231 !---------------------------------------------------------------------
1232 integer :: astat, ierr
1233 integer :: m
1234 integer :: ncid, dimid, varid
1235 character(len=132) :: filename
1236 character(len=132) :: err_msg
1237 character(len=64) :: varname
1239 #ifdef NETCDF
1240 include 'netcdf.inc'
1242 !---------------------------------------------------------------------
1243 ! ... function declarations
1244 !---------------------------------------------------------------------
1245 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
1247 master_proc_a: &
1248 if( wrf_dm_on_monitor() ) then
1249 filename = 'wrf_tuv_xsqy.nc'
1250 err_msg = 'get_xsqy_tab: failed to open file ' // trim(filename)
1251 call handle_ncerr( nf_open( trim(filename), nf_noclobber, ncid ), trim(err_msg) )
1252 !---------------------------------------------------------------------
1253 ! ... get dimensions
1254 !---------------------------------------------------------------------
1255 err_msg = 'get_xsqy_tab: failed to get nwave id'
1256 call handle_ncerr( nf_inq_dimid( ncid, 'nwave', dimid ), trim(err_msg) )
1257 err_msg = 'get_xsqy_tab: failed to get nwave'
1258 call handle_ncerr( nf_inq_dimlen( ncid, dimid, nwave ), trim(err_msg) )
1259 err_msg = 'get_xsqy_tab: failed to get ntemp id'
1260 call handle_ncerr( nf_inq_dimid( ncid, 'ntemp', dimid ), trim(err_msg) )
1261 err_msg = 'get_xsqy_tab: failed to get ntemp'
1262 call handle_ncerr( nf_inq_dimlen( ncid, dimid, ntemp ), trim(err_msg) )
1263 err_msg = 'get_xsqy_tab: failed to get nconc id'
1264 call handle_ncerr( nf_inq_dimid( ncid, 'nconc', dimid ), trim(err_msg) )
1265 err_msg = 'get_xsqy_tab: failed to get nconc'
1266 call handle_ncerr( nf_inq_dimlen( ncid, dimid, nconc ), trim(err_msg) )
1267 err_msg = 'get_xsqy_tab: failed to get nchebev_term id'
1268 call handle_ncerr( nf_inq_dimid( ncid, 'nchebev_term', dimid ), trim(err_msg) )
1269 err_msg = 'get_xsqy_tab: failed to get nchebev'
1270 call handle_ncerr( nf_inq_dimlen( ncid, dimid, nchebev_term ), trim(err_msg) )
1271 err_msg = 'get_xsqy_tab: failed to get nchebev_wave id'
1272 call handle_ncerr( nf_inq_dimid( ncid, 'nchebev_wave', dimid ), trim(err_msg) )
1273 err_msg = 'get_xsqy_tab: failed to get nchebev'
1274 call handle_ncerr( nf_inq_dimlen( ncid, dimid, nchebev_wave ), trim(err_msg) )
1275 err_msg = 'get_xsqy_tab: failed to get n_temp_data id'
1276 call handle_ncerr( nf_inq_dimid( ncid, 'n_temp_data', dimid ), trim(err_msg) )
1277 err_msg = 'get_xsqy_tab: failed to get n_temp_data'
1278 call handle_ncerr( nf_inq_dimlen( ncid, dimid, n_temp_data ), trim(err_msg) )
1279 err_msg = 'get_xsqy_tab: failed to get n_o3_data id'
1280 call handle_ncerr( nf_inq_dimid( ncid, 'n_o3_data', dimid ), trim(err_msg) )
1281 err_msg = 'get_xsqy_tab: failed to get n_o3_data'
1282 call handle_ncerr( nf_inq_dimlen( ncid, dimid, n_o3_data ), trim(err_msg) )
1283 err_msg = 'get_xsqy_tab: failed to get n_air_dens_data id'
1284 call handle_ncerr( nf_inq_dimid( ncid, 'n_air_dens_data', dimid ), trim(err_msg) )
1285 err_msg = 'get_xsqy_tab: failed to get n_air_dens_data'
1286 call handle_ncerr( nf_inq_dimlen( ncid, dimid, n_air_dens_data ), trim(err_msg) )
1287 endif master_proc_a
1288 #ifdef DM_PARALLEL
1289 !---------------------------------------------------------------------
1290 ! ... bcast the dimensions
1291 !---------------------------------------------------------------------
1292 CALL wrf_dm_bcast_bytes ( nwave, IWORDSIZE )
1293 CALL wrf_dm_bcast_bytes ( ntemp, IWORDSIZE )
1294 CALL wrf_dm_bcast_bytes ( nconc, IWORDSIZE )
1295 CALL wrf_dm_bcast_bytes ( n_temp_data, IWORDSIZE )
1296 CALL wrf_dm_bcast_bytes ( n_o3_data, IWORDSIZE )
1297 CALL wrf_dm_bcast_bytes ( n_air_dens_data, IWORDSIZE )
1298 CALL wrf_dm_bcast_bytes ( nchebev_term, IWORDSIZE )
1299 CALL wrf_dm_bcast_bytes ( nchebev_wave, IWORDSIZE )
1300 #endif
1301 !---------------------------------------------------------------------
1302 ! ... allocate module arrays
1303 !---------------------------------------------------------------------
1304 ierr = 0
1305 allocate( z_temp_data(n_temp_data), z_o3_data(n_o3_data), &
1306 z_air_dens_data(n_air_dens_data),stat=astat )
1307 ierr = astat + ierr
1308 allocate( temp_data(n_temp_data), o3_data(n_o3_data), &
1309 air_dens_data(n_air_dens_data),stat=astat )
1310 ierr = astat + ierr
1311 allocate( wl(nwave+1), wc(nwave), dw(nwave), w_fac(nwave), &
1312 etfl(nwave), albedo(nwave), stat=astat )
1313 ierr = astat + ierr
1314 if( .not. is_full_tuv ) then
1315 allocate( temp_tab(ntemp), conc_tab(nconc), stat=astat )
1316 ierr = astat + ierr
1317 allocate( del_temp_tab(ntemp-1), del_conc_tab(nconc-1), stat=astat )
1318 ierr = astat + ierr
1319 endif
1320 allocate( chebev_ac(nchebev_term,nchebev_wave), stat=astat )
1321 ierr = astat + ierr
1322 allocate( chebev_bc(nchebev_term,nchebev_wave), stat=astat )
1323 ierr = astat + ierr
1324 allocate( o2_xs(nwave), so2_xs(nwave), stat=astat )
1325 ierr = astat + ierr
1326 allocate( o3_xs_tab(nwave,ntemp), no2_xs_tab(nwave,ntemp), stat=astat )
1327 ierr = astat + ierr
1328 if( .not. is_full_tuv ) then
1329 allocate( xsqy_tab(nwave,ntemp,nconc,nj), stat=astat )
1330 ierr = astat + ierr
1331 endif
1332 if( ierr /= 0 ) then
1333 call wrf_error_fatal( 'tuv_init: failed to allocate z_temp_data ... xsqy_tab' )
1334 endif
1335 !---------------------------------------------------------------------
1336 ! ... read arrays
1337 !---------------------------------------------------------------------
1338 master_proc_b: &
1339 if( wrf_dm_on_monitor() ) then
1340 err_msg = 'get_xsqy_tab: failed to get z_temp_data variable id'
1341 call handle_ncerr( nf_inq_varid( ncid, 'z_temp_data', varid ), trim(err_msg) )
1342 err_msg = 'get_xsqy_tab: failed to read z_temp_data variable'
1343 call handle_ncerr( nf_get_var_real( ncid, varid, z_temp_data ), trim(err_msg) )
1345 err_msg = 'get_xsqy_tab: failed to get z_o3_data variable id'
1346 call handle_ncerr( nf_inq_varid( ncid, 'z_o3_data', varid ), trim(err_msg) )
1347 err_msg = 'get_xsqy_tab: failed to read z_o3_data variable'
1348 call handle_ncerr( nf_get_var_real( ncid, varid, z_o3_data ), trim(err_msg) )
1350 err_msg = 'get_xsqy_tab: failed to get z_air_dens_data variable id'
1351 call handle_ncerr( nf_inq_varid( ncid, 'z_air_dens_data', varid ), trim(err_msg) )
1352 err_msg = 'get_xsqy_tab: failed to read z_air_dens_data variable'
1353 call handle_ncerr( nf_get_var_real( ncid, varid, z_air_dens_data ), trim(err_msg) )
1355 err_msg = 'get_xsqy_tab: failed to get temp_data variable id'
1356 call handle_ncerr( nf_inq_varid( ncid, 'temp_data', varid ), trim(err_msg) )
1357 err_msg = 'get_xsqy_tab: failed to read temp_data variable'
1358 call handle_ncerr( nf_get_var_real( ncid, varid, temp_data ), trim(err_msg) )
1360 err_msg = 'get_xsqy_tab: failed to get o3_data variable id'
1361 call handle_ncerr( nf_inq_varid( ncid, 'o3_data', varid ), trim(err_msg) )
1362 err_msg = 'get_xsqy_tab: failed to read o3_data variable'
1363 call handle_ncerr( nf_get_var_real( ncid, varid, o3_data ), trim(err_msg) )
1365 err_msg = 'get_xsqy_tab: failed to get air_dens_data variable id'
1366 call handle_ncerr( nf_inq_varid( ncid, 'air_dens_data', varid ), trim(err_msg) )
1367 err_msg = 'get_xsqy_tab: failed to read air_dens_data variable'
1368 call handle_ncerr( nf_get_var_real( ncid, varid, air_dens_data ), trim(err_msg) )
1370 err_msg = 'get_xsqy_tab: failed to get wl variable id'
1371 call handle_ncerr( nf_inq_varid( ncid, 'wl', varid ), trim(err_msg) )
1372 err_msg = 'get_xsqy_tab: failed to read wl variable'
1373 call handle_ncerr( nf_get_var_real( ncid, varid, wl ), trim(err_msg) )
1374 err_msg = 'get_xsqy_tab: failed to get wc variable id'
1375 call handle_ncerr( nf_inq_varid( ncid, 'wc', varid ), trim(err_msg) )
1376 err_msg = 'get_xsqy_tab: failed to read wc variable'
1377 call handle_ncerr( nf_get_var_real( ncid, varid, wc ), trim(err_msg) )
1378 err_msg = 'get_xsqy_tab: failed to get etfl variable id'
1379 call handle_ncerr( nf_inq_varid( ncid, 'etf', varid ), trim(err_msg) )
1380 err_msg = 'get_xsqy_tab: failed to read etfl variable'
1381 call handle_ncerr( nf_get_var_real( ncid, varid, etfl ), trim(err_msg) )
1383 err_msg = 'get_xsqy_tab: failed to get chebev_ac variable id'
1384 call handle_ncerr( nf_inq_varid( ncid, 'chebev_ac', varid ), trim(err_msg) )
1385 err_msg = 'get_xsqy_tab: failed to read chebev_ac variable'
1386 call handle_ncerr( nf_get_var_double( ncid, varid, chebev_ac ), trim(err_msg) )
1387 err_msg = 'get_xsqy_tab: failed to get chebev_bc variable id'
1388 call handle_ncerr( nf_inq_varid( ncid, 'chebev_bc', varid ), trim(err_msg) )
1389 err_msg = 'get_xsqy_tab: failed to read chebev_bc variable'
1390 call handle_ncerr( nf_get_var_double( ncid, varid, chebev_bc ), trim(err_msg) )
1392 err_msg = 'get_xsqy_tab: failed to get ila variable id'
1393 call handle_ncerr( nf_inq_varid( ncid, 'ila', varid ), trim(err_msg) )
1394 err_msg = 'get_xsqy_tab: failed to read ila variable'
1395 call handle_ncerr( nf_get_var_int( ncid, varid, ila ), trim(err_msg) )
1397 err_msg = 'get_xsqy_tab: failed to get isrb variable id'
1398 call handle_ncerr( nf_inq_varid( ncid, 'isrb', varid ), trim(err_msg) )
1399 err_msg = 'get_xsqy_tab: failed to read isrb variable'
1400 call handle_ncerr( nf_get_var_int( ncid, varid, isrb ), trim(err_msg) )
1402 if( .not. is_full_tuv ) then
1403 err_msg = 'get_xsqy_tab: failed to temp_tab variable id'
1404 call handle_ncerr( nf_inq_varid( ncid, 'temps', varid ), trim(err_msg) )
1405 err_msg = 'get_xsqy_tab: failed to read temp_tab variable'
1406 call handle_ncerr( nf_get_var_real( ncid, varid, temp_tab ), trim(err_msg) )
1407 err_msg = 'get_xsqy_tab: failed to conc_tab variable id'
1408 call handle_ncerr( nf_inq_varid( ncid, 'concs', varid ), trim(err_msg) )
1409 err_msg = 'get_xsqy_tab: failed to read conc_tab variable'
1410 call handle_ncerr( nf_get_var_real( ncid, varid, conc_tab ), trim(err_msg) )
1411 endif
1412 err_msg = 'get_xsqy_tab: failed to o2_xs variable id'
1413 call handle_ncerr( nf_inq_varid( ncid, 'o2_xs', varid ), trim(err_msg) )
1414 err_msg = 'get_xsqy_tab: failed to read o2_xs variable'
1415 call handle_ncerr( nf_get_var_real( ncid, varid, o2_xs ), trim(err_msg) )
1416 err_msg = 'get_xsqy_tab: failed to so2_xs variable id'
1417 call handle_ncerr( nf_inq_varid( ncid, 'so2_xs', varid ), trim(err_msg) )
1418 err_msg = 'get_xsqy_tab: failed to read so2_xs variable'
1419 call handle_ncerr( nf_get_var_real( ncid, varid, so2_xs ), trim(err_msg) )
1420 err_msg = 'get_xsqy_tab: failed to o3_xs_tab variable id'
1421 call handle_ncerr( nf_inq_varid( ncid, 'o3_xs', varid ), trim(err_msg) )
1422 err_msg = 'get_xsqy_tab: failed to read o3_xs_tab variable'
1423 call handle_ncerr( nf_get_var_real( ncid, varid, o3_xs_tab ), trim(err_msg) )
1424 err_msg = 'get_xsqy_tab: failed to no2_xs_tab variable id'
1425 call handle_ncerr( nf_inq_varid( ncid, 'no2_xs', varid ), trim(err_msg) )
1426 err_msg = 'get_xsqy_tab: failed to read no2_xs_tab variable'
1427 call handle_ncerr( nf_get_var_real( ncid, varid, no2_xs_tab ), trim(err_msg) )
1428 if( .not. is_full_tuv ) then
1429 do m = 1,nj
1430 varname = trim(tuv_jname(m)) // '_xsqy'
1431 err_msg = 'get_xsqy_tab: failed to ' // trim(varname) //' variable id'
1432 call handle_ncerr( nf_inq_varid( ncid, trim(varname), varid ), trim(err_msg) )
1433 err_msg = 'get_xsqy_tab: failed to read ' // trim(varname) // ' variable'
1434 call handle_ncerr( nf_get_var_real( ncid, varid, xsqy_tab(:,:,:,m) ), trim(err_msg) )
1435 end do
1436 endif
1437 endif master_proc_b
1439 #ifdef DM_PARALLEL
1440 !---------------------------------------------------------------------
1441 ! ... bcast the arrays
1442 !---------------------------------------------------------------------
1443 CALL wrf_dm_bcast_bytes ( ila, IWORDSIZE )
1444 CALL wrf_dm_bcast_bytes ( isrb, IWORDSIZE )
1445 CALL wrf_dm_bcast_bytes( wl, (nwave+1)*RWORDSIZE )
1446 CALL wrf_dm_bcast_bytes( wc, nwave*RWORDSIZE )
1447 CALL wrf_dm_bcast_bytes( etfl, nwave*RWORDSIZE )
1448 if( .not. is_full_tuv ) then
1449 CALL wrf_dm_bcast_bytes( temp_tab, ntemp*RWORDSIZE )
1450 CALL wrf_dm_bcast_bytes( conc_tab, nconc*RWORDSIZE )
1451 endif
1452 CALL wrf_dm_bcast_bytes( o2_xs, nwave*RWORDSIZE )
1453 CALL wrf_dm_bcast_bytes( so2_xs, nwave*RWORDSIZE )
1454 CALL wrf_dm_bcast_bytes( z_temp_data, n_temp_data*RWORDSIZE )
1455 CALL wrf_dm_bcast_bytes( z_o3_data, n_o3_data*RWORDSIZE )
1456 CALL wrf_dm_bcast_bytes( z_air_dens_data, n_air_dens_data*RWORDSIZE )
1457 CALL wrf_dm_bcast_bytes( temp_data, n_temp_data*RWORDSIZE )
1458 CALL wrf_dm_bcast_bytes( o3_data, n_o3_data*RWORDSIZE )
1459 CALL wrf_dm_bcast_bytes( air_dens_data, n_air_dens_data*RWORDSIZE )
1460 #if RWORDSIZE == 4
1461 CALL wrf_dm_bcast_bytes( chebev_ac, nchebev_term*nchebev_wave*2*RWORDSIZE )
1462 CALL wrf_dm_bcast_bytes( chebev_bc, nchebev_term*nchebev_wave*2*RWORDSIZE )
1463 #endif
1464 #if RWORDSIZE == 8
1465 CALL wrf_dm_bcast_bytes( chebev_ac, nchebev_term*nchebev_wave*RWORDSIZE )
1466 CALL wrf_dm_bcast_bytes( chebev_bc, nchebev_term*nchebev_wave*RWORDSIZE )
1467 #endif
1468 CALL wrf_dm_bcast_bytes( o3_xs_tab, nwave*ntemp*RWORDSIZE )
1469 CALL wrf_dm_bcast_bytes( no2_xs_tab, nwave*ntemp*RWORDSIZE )
1470 if( .not. is_full_tuv ) then
1471 CALL wrf_dm_bcast_bytes( xsqy_tab, nwave*ntemp*nconc*nj*RWORDSIZE )
1472 endif
1473 #endif
1475 if( .not. is_full_tuv ) then
1476 del_temp_tab(:ntemp-1) = 1./(temp_tab(2:ntemp) - temp_tab(1:ntemp-1))
1477 del_conc_tab(:nconc-1) = 1./(conc_tab(2:nconc) - conc_tab(1:nconc-1))
1478 endif
1479 dw(:nwave) = wl(2:nwave+1) - wl(1:nwave)
1480 w_fac(:nwave) = dw(:nwave)*etfl(:nwave)*1.e-13*wc(:nwave)/hc
1482 if( wrf_dm_on_monitor() ) then
1483 err_msg = 'get_xsqy_tab: failed to close file ' // trim(filename)
1484 call handle_ncerr( nf_close( ncid ),trim(err_msg) )
1485 endif
1486 #endif
1488 end subroutine get_xsqy_tab
1490 subroutine get_exo_coldens( dm, exo_coldens_filename, &
1491 ids,ide, jds,jde, kds,kde, &
1492 ims,ime, jms,jme, kms,kme, &
1493 its,ite, jts,jte, kts,kte )
1494 !---------------------------------------------------------------------
1495 ! ... read in the exo column o2,o3 densities
1496 !---------------------------------------------------------------------
1498 !---------------------------------------------------------------------
1499 ! ... dummy arguments
1500 !---------------------------------------------------------------------
1501 integer, intent(in) :: dm
1502 integer, intent(in) :: ide, ids, ime, ims, ite, its, jde, jds, &
1503 jme, jms, jte, jts, kde, kds, kme, kms, kte, kts
1504 character(len=*), intent(in) :: exo_coldens_filename
1506 !---------------------------------------------------------------------
1507 ! ... local variables
1508 !---------------------------------------------------------------------
1509 INTEGER :: i, j, k
1510 integer :: astat
1511 integer :: ncid
1512 integer :: dimid
1513 integer :: varid
1514 integer :: max_dom
1515 integer :: cpos
1516 integer :: iend, jend
1517 integer :: lon_e, lat_e
1518 integer :: ncoldens_levs
1519 integer :: ndays_of_year
1520 ! real, allocatable :: coldens(:,:,:,:)
1521 character(len=128) :: err_msg
1522 character(len=64) :: filename
1523 character(len=2) :: id_num
1525 !---------------------------------------------------------------------
1526 ! ... function declarations
1527 !---------------------------------------------------------------------
1528 LOGICAL, EXTERNAL :: wrf_dm_on_monitor
1530 #ifdef NETCDF
1531 include 'netcdf.inc'
1532 #define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes ( A , size ( A ) * DWORDSIZE )
1534 !---------------------------------------------------------------------
1535 ! ... allocate column_density type
1536 !---------------------------------------------------------------------
1537 if( .not. allocated(col_dens) ) then
1538 CALL nl_get_max_dom( 1,max_dom )
1539 allocate( col_dens(max_dom),stat=astat )
1540 if( astat /= 0 ) then
1541 call wrf_error_fatal( 'get_exo_coldens: failed to allocate col_dens' )
1542 endif
1543 write(err_msg,'(''get_exo_coldens: intializing '',i2,'' domains'')') max_dom
1544 call wrf_message( trim(err_msg) )
1545 col_dens(:)%is_allocated = .false.
1546 endif
1547 !---------------------------------------------------------------------
1548 ! ... open column density netcdf file
1549 !---------------------------------------------------------------------
1550 col_dens_allocated : &
1551 if( .not. col_dens(dm)%is_allocated ) then
1552 ! if( wrf_dm_on_monitor() ) then
1553 cpos = index( exo_coldens_filename, '<domain>' )
1554 if( cpos > 0 ) then
1555 write(id_num,'(i2.2)') dm
1556 filename = exo_coldens_filename(:cpos-1) // 'd' // id_num
1557 else
1558 filename = trim( exo_coldens_filename )
1559 endif
1560 err_msg = 'get_exo_coldens: intializing domain ' // id_num
1561 call wrf_message( trim(err_msg) )
1562 err_msg = 'get_exo_coldens: failed to open file ' // trim(filename)
1563 call handle_ncerr( nf_open( trim(filename), nf_noclobber, ncid ), trim(err_msg) )
1564 !---------------------------------------------------------------------
1565 ! ... get dimensions
1566 !---------------------------------------------------------------------
1567 err_msg = 'get_exo_coldens: failed to get col_dens levels id'
1568 call handle_ncerr( nf_inq_dimid( ncid, 'coldens_levs', dimid ), trim(err_msg) )
1569 err_msg = 'get_exo_coldens: failed to get col_dens levels'
1570 call handle_ncerr( nf_inq_dimlen( ncid, dimid, ncoldens_levs ), trim(err_msg) )
1571 err_msg = 'get_exo_coldens: failed to get number of days in year id'
1572 call handle_ncerr( nf_inq_dimid( ncid, 'ndays_of_year', dimid ), trim(err_msg) )
1573 err_msg = 'get_exo_coldens: failed to get number of days in year'
1574 call handle_ncerr( nf_inq_dimlen( ncid, dimid, ndays_of_year ), trim(err_msg) )
1575 err_msg = 'get_exo_coldens: failed to get west_east id'
1576 call handle_ncerr( nf_inq_dimid( ncid, 'west_east', dimid ), trim(err_msg) )
1577 err_msg = 'get_exo_coldens: failed to get west_east'
1578 call handle_ncerr( nf_inq_dimlen( ncid, dimid, lon_e ), trim(err_msg) )
1579 err_msg = 'get_exo_coldens: failed to get south_north id'
1580 call handle_ncerr( nf_inq_dimid( ncid, 'south_north', dimid ), trim(err_msg) )
1581 err_msg = 'get_exo_coldens: failed to get south_north'
1582 call handle_ncerr( nf_inq_dimlen( ncid, dimid, lat_e ), trim(err_msg) )
1583 ! end IF
1584 #ifdef DM_PARALLEL
1585 !---------------------------------------------------------------------
1586 ! ... bcast the dimensions
1587 !---------------------------------------------------------------------
1588 ! CALL wrf_dm_bcast_bytes ( ncoldens_levs , IWORDSIZE )
1589 ! CALL wrf_dm_bcast_bytes ( ndays_of_year , IWORDSIZE )
1590 ! CALL wrf_dm_bcast_bytes ( lon_e , IWORDSIZE )
1591 ! CALL wrf_dm_bcast_bytes ( lat_e , IWORDSIZE )
1592 #endif
1593 !---------------------------------------------------------------------
1594 ! ... allocate local arrays
1595 !---------------------------------------------------------------------
1596 iend = min( ite,ide-1 )
1597 jend = min( jte,jde-1 )
1598 ! allocate( coldens(lon_e,lat_e,ncoldens_levs,ndays_of_year), stat=astat )
1599 ! if( astat /= 0 ) then
1600 ! call wrf_error_fatal( 'get_exo_coldens: failed to allocate coldens' )
1601 ! end if
1602 !---------------------------------------------------------------------
1603 ! ... allocate column_density type component arrays
1604 !---------------------------------------------------------------------
1605 col_dens(dm)%ncoldens_levs = ncoldens_levs
1606 col_dens(dm)%ndays_of_year = ndays_of_year
1607 allocate( col_dens(dm)%col_levs(ncoldens_levs), &
1608 col_dens(dm)%day_of_year(ndays_of_year), stat=astat )
1609 if( astat /= 0 ) then
1610 call wrf_error_fatal( 'get_exo_coldens: failed to allocate col_levs,day_of_year' )
1611 end if
1612 allocate( col_dens(dm)%o3_col_dens(its:iend,jts:jend,ncoldens_levs,ndays_of_year), &
1613 col_dens(dm)%o2_col_dens(its:iend,jts:jend,ncoldens_levs,ndays_of_year), stat=astat )
1614 if( astat /= 0 ) then
1615 call wrf_error_fatal( 'get_exo_coldens: failed to allocate o3_col_dens,o2_col_dens' )
1616 end if
1617 col_dens(dm)%is_allocated = .true.
1618 !---------------------------------------------------------------------
1619 ! ... read arrays
1620 !---------------------------------------------------------------------
1621 ! IF ( wrf_dm_on_monitor() ) THEN
1622 err_msg = 'get_exo_coldens: failed to get col_levs variable id'
1623 call handle_ncerr( nf_inq_varid( ncid, 'coldens_levs', varid ), trim(err_msg) )
1624 err_msg = 'get_exo_coldens: failed to read col_levs variable'
1625 call handle_ncerr( nf_get_var_double( ncid, varid, col_dens(dm)%col_levs ), trim(err_msg) )
1626 err_msg = 'get_exo_coldens: failed to get days_of_year variable id'
1627 call handle_ncerr( nf_inq_varid( ncid, 'days_of_year', varid ), trim(err_msg) )
1628 err_msg = 'get_exo_coldens: failed to read days_of_year variable'
1629 call handle_ncerr( nf_get_var_double( ncid, varid, col_dens(dm)%day_of_year ), trim(err_msg) )
1630 err_msg = 'get_exo_coldens: failed to get o3 col_dens variable id'
1631 call handle_ncerr( nf_inq_varid( ncid, 'o3_column_density', varid ), trim(err_msg) )
1632 err_msg = 'get_exo_coldens: failed to read o3 col_dens variable'
1633 ! call handle_ncerr( nf_get_var_real( ncid, varid, coldens ), trim(err_msg) )
1634 call handle_ncerr( nf_get_vara_double( ncid, varid, (/its,jts,1,1/), &
1635 (/iend-its+1,jend-jts+1,ncoldens_levs,ndays_of_year/), &
1636 col_dens(dm)%o3_col_dens(its:iend,jts:jend,1:ncoldens_levs,1:ndays_of_year) ), trim(err_msg) )
1637 ! ENDIF
1639 #ifdef DM_PARALLEL
1640 ! call wrf_debug( 0,'get_exo_coldens: bcast col_levs' )
1641 ! DM_BCAST_MACRO(col_dens(dm)%col_levs)
1642 ! call wrf_debug( 0,'get_exo_coldens: bcast day_of_year' )
1643 ! DM_BCAST_MACRO(col_dens(dm)%day_of_year)
1644 ! call wrf_debug( 0,'get_exo_coldens: bcast o3_col_dens' )
1645 ! CALL wrf_dm_bcast_bytes ( coldens, size(coldens)*RWORDSIZE )
1646 #endif
1648 ! col_dens(dm)%o3_col_dens(its:iend,jts:jend,:ncoldens_levs,:ndays_of_year) = &
1649 ! coldens(its:iend,jts:jend,:ncoldens_levs,:ndays_of_year)
1651 ! IF ( wrf_dm_on_monitor() ) THEN
1652 err_msg = 'get_exo_coldens: failed to get o2 col_dens variable id'
1653 call handle_ncerr( nf_inq_varid( ncid, 'o2_column_density', varid ), trim(err_msg) )
1654 err_msg = 'get_exo_coldens: failed to read o2 col_dens variable'
1655 ! call handle_ncerr( nf_get_var_real( ncid, varid, coldens ), trim(err_msg) )
1656 call handle_ncerr( nf_get_vara_double( ncid, varid, (/its,jts,1,1/), &
1657 (/iend-its+1,jend-jts+1,ncoldens_levs,ndays_of_year/), &
1658 col_dens(dm)%o2_col_dens(its:iend,jts:jend,1:ncoldens_levs,1:ndays_of_year) ), trim(err_msg) )
1659 !---------------------------------------------------------------------
1660 ! ... close column density netcdf file
1661 !---------------------------------------------------------------------
1662 err_msg = 'get_exo_coldens: failed to close file ' // trim(filename)
1663 call handle_ncerr( nf_close( ncid ), trim(err_msg) )
1664 ! end if
1666 #ifdef DM_PARALLEL
1667 ! call wrf_debug( 0,'get_exo_coldens: bcast o2_col_dens' )
1668 ! CALL wrf_dm_bcast_bytes ( coldens, size(coldens)*RWORDSIZE )
1669 #endif
1671 ! col_dens(dm)%o2_col_dens(its:iend,jts:jend,:ncoldens_levs,:ndays_of_year) = &
1672 ! coldens(its:iend,jts:jend,:ncoldens_levs,:ndays_of_year)
1674 ! deallocate( coldens )
1676 !---------------------------------------------------------------------
1677 ! ... some diagnostics
1678 !---------------------------------------------------------------------
1679 call wrf_debug( 100,' ' )
1680 write(err_msg,'(''get_exo_coldens: coldens variables for domain '',i2)') dm
1681 call wrf_debug( 100,trim(err_msg) )
1682 call wrf_debug( 100,'get_exo_coldens: days_of_year' )
1683 do k = 1,ndays_of_year,5
1684 write(err_msg,'(1p,5g15.7)') col_dens(dm)%day_of_year(k:min(k+4,ndays_of_year))
1685 call wrf_debug( 100,trim(err_msg) )
1686 end do
1687 call wrf_debug( 100,'get_exo_coldens: coldens levels' )
1688 do k = 1,ncoldens_levs,5
1689 write(err_msg,'(1p,5g15.7)') col_dens(dm)%col_levs(k:min(k+4,ncoldens_levs))
1690 call wrf_debug( 100,trim(err_msg) )
1691 end do
1692 call wrf_debug( 100,' ' )
1693 endif col_dens_allocated
1694 #endif
1696 end subroutine get_exo_coldens
1698 subroutine tuv_timestep_init( dm, julday )
1699 !-----------------------------------------------------------------------------
1700 ! ... setup the exo column time interpolation
1701 !-----------------------------------------------------------------------------
1703 implicit none
1705 !-----------------------------------------------------------------------------
1706 ! ... dummy arguments
1707 !-----------------------------------------------------------------------------
1708 integer, intent(in) :: dm ! domain number
1709 integer, intent(in) :: julday ! day of year at present time step
1711 !-----------------------------------------------------------------------------
1712 ! ... local variables
1713 !-----------------------------------------------------------------------------
1714 integer :: m
1715 real(dp) :: calday
1717 calday = real( julday,kind=dp)
1718 if( calday < col_dens(dm)%day_of_year(1) ) then
1719 next = 1
1720 last = 12
1721 dels = (365._dp + calday - col_dens(dm)%day_of_year(12)) &
1722 / (365._dp + col_dens(dm)%day_of_year(1) - col_dens(dm)%day_of_year(12))
1723 else if( calday >= col_dens(dm)%day_of_year(12) ) then
1724 next = 1
1725 last = 12
1726 dels = (calday - col_dens(dm)%day_of_year(12)) &
1727 / (365. + col_dens(dm)%day_of_year(1) - col_dens(dm)%day_of_year(12))
1728 else
1729 do m = 11,1,-1
1730 if( calday >= col_dens(dm)%day_of_year(m) ) then
1731 exit
1732 end if
1733 end do
1734 last = m
1735 next = m + 1
1736 dels = (calday - col_dens(dm)%day_of_year(m)) / (col_dens(dm)%day_of_year(m+1) - col_dens(dm)%day_of_year(m))
1737 end if
1739 end subroutine tuv_timestep_init
1741 subroutine z_interp( z_tab, tab, ntab, z_out, out )
1743 integer, intent(in) :: ntab
1744 real, intent(in) :: z_tab(ntab)
1745 real, intent(in) :: tab(ntab)
1746 real, intent(in) :: z_out(:)
1747 real, intent(out) :: out(:)
1748 !---------------------------------------------------------------
1749 ! ... altitude interpolation
1750 !---------------------------------------------------------------
1752 integer :: k, kt, ktm1, n
1753 real :: delz
1755 n = size(out)
1756 do k = 1,n
1757 if( z_out(k) <= z_tab(1) ) then
1758 out(k) = z_tab(1)
1759 elseif( z_out(k) >= z_tab(ntab) ) then
1760 out(k) = z_tab(ntab)
1761 else
1762 do kt = 1,ntab
1763 if( z_tab(kt) >= z_out(k) ) then
1764 ktm1 = kt - 1
1765 delz = (z_out(k) - z_tab(ktm1))/(z_tab(kt) - z_tab(ktm1))
1766 out(k) = tab(ktm1) + delz*(tab(kt) - tab(ktm1))
1767 exit
1768 endif
1769 end do
1770 endif
1771 end do
1773 end subroutine z_interp
1775 subroutine p_interp( o2_exo_col, o3_exo_col, o3_exo_col_at_grnd, ptop, &
1776 dm, its, ite, jts, jte )
1777 !---------------------------------------------------------------
1778 ! ... pressure interpolation for exo col density
1779 !---------------------------------------------------------------
1781 implicit none
1783 !---------------------------------------------------------------
1784 ! ... dummy arguments
1785 !---------------------------------------------------------------
1786 integer, intent(in) :: dm
1787 integer, intent(in) :: its, ite
1788 integer, intent(in) :: jts, jte
1789 real(dp), intent(in) :: ptop(its:ite,jts:jte) ! pressure at photolysis top (hPa)
1790 real(dp), intent(out) :: o2_exo_col(its:ite,jts:jte) ! exo model o2 column density (molecules/cm^2)
1791 real(dp), intent(out) :: o3_exo_col(its:ite,jts:jte) ! exo model o3 column density (molecules/cm^2)
1792 real(dp), intent(out) :: o3_exo_col_at_grnd(its:ite,jts:jte) ! exo model o3 column density at grnd (molecules/cm^2)
1794 !---------------------------------------------------------------
1795 ! ... local variables
1796 !---------------------------------------------------------------
1797 integer :: i, j, k, ku, kl
1798 integer :: Kgrnd
1799 real(dp) :: pinterp
1800 real(dp) :: delp
1801 real(dp) :: tint_vals(2)
1803 Kgrnd = col_dens(dm)%ncoldens_levs
1805 lat_loop : &
1806 do j = jts,jte
1807 long_loop : &
1808 do i = its,ite
1809 pinterp = ptop(i,j)
1810 if( pinterp < col_dens(dm)%col_levs(1) ) then
1811 ku = 1
1812 kl = 1
1813 delp = 0._dp
1814 else
1815 do ku = 2,col_dens(dm)%ncoldens_levs
1816 if( pinterp <= col_dens(dm)%col_levs(ku) ) then
1817 kl = ku - 1
1818 delp = log( pinterp/col_dens(dm)%col_levs(kl) )/log( col_dens(dm)%col_levs(ku)/col_dens(dm)%col_levs(kl) )
1819 exit
1820 end if
1821 end do
1822 end if
1823 tint_vals(1) = col_dens(dm)%o2_col_dens(i,j,kl,last) &
1824 + delp * (col_dens(dm)%o2_col_dens(i,j,ku,last) &
1825 - col_dens(dm)%o2_col_dens(i,j,kl,last))
1826 tint_vals(2) = col_dens(dm)%o2_col_dens(i,j,kl,next) &
1827 + delp * (col_dens(dm)%o2_col_dens(i,j,ku,next) &
1828 - col_dens(dm)%o2_col_dens(i,j,kl,next))
1829 o2_exo_col(i,j) = tint_vals(1) + dels * (tint_vals(2) - tint_vals(1))
1830 tint_vals(1) = col_dens(dm)%o3_col_dens(i,j,kl,last) &
1831 + delp * (col_dens(dm)%o3_col_dens(i,j,ku,last) &
1832 - col_dens(dm)%o3_col_dens(i,j,kl,last))
1833 tint_vals(2) = col_dens(dm)%o3_col_dens(i,j,kl,next) &
1834 + delp * (col_dens(dm)%o3_col_dens(i,j,ku,next) &
1835 - col_dens(dm)%o3_col_dens(i,j,kl,next))
1836 o3_exo_col(i,j) = tint_vals(1) + dels * (tint_vals(2) - tint_vals(1))
1837 tint_vals(1) = col_dens(dm)%o3_col_dens(i,j,Kgrnd,last)
1838 tint_vals(2) = col_dens(dm)%o3_col_dens(i,j,Kgrnd,next)
1839 o3_exo_col_at_grnd(i,j) = tint_vals(1) + dels * (tint_vals(2) - tint_vals(1))
1840 end do long_loop
1841 end do lat_loop
1843 end subroutine p_interp
1845 subroutine xsqy_int( n, xsqy, tlev, dens_air )
1846 !---------------------------------------------------------------------
1847 ! ... interpolate m,t tables for xs * qy
1848 !---------------------------------------------------------------------
1850 integer, intent(in) :: n
1851 real, intent(in) :: tlev(:)
1852 real, intent(in) :: dens_air(:)
1853 real, intent(out) :: xsqy(:,:)
1855 real, parameter :: m0 = 2.45e19
1856 integer :: tndx, mndx, tndxp1, mndxp1
1857 integer :: k, ku
1858 real :: temp, dens
1859 real :: w(4)
1860 real :: del_t, del_d
1862 ku = size( tlev )
1863 do k = 1,ku
1864 temp = tlev(k)
1865 do tndx = 1,ntemp
1866 if( temp_tab(tndx) > temp ) then
1867 exit
1868 endif
1869 end do
1870 tndx = max( min( tndx,ntemp ) - 1,1 )
1871 tndxp1 = tndx + 1
1872 del_t = max( 0.,min( 1.,(temp - temp_tab(tndx))*del_temp_tab(tndx) ) )
1874 ! dens = dens_air(k)
1875 dens = dens_air(k)/m0
1876 do mndx = 1,nconc
1877 if( conc_tab(mndx) > dens ) then
1878 exit
1879 endif
1880 end do
1881 mndx = max( min( mndx,nconc ) - 1,1 )
1882 mndxp1 = mndx + 1
1883 del_d = max( 0.,min( 1.,(dens - conc_tab(mndx))*del_conc_tab(mndx) ) )
1885 w(1) = (1. - del_t)*(1. - del_d)
1886 w(2) = del_t*(1. - del_d)
1887 w(3) = (1. - del_t)*del_d
1888 w(4) = del_t*del_d
1890 xsqy(1:nwave,k) = w(1)*xsqy_tab(1:nwave,tndx,mndx,n) &
1891 + w(2)*xsqy_tab(1:nwave,tndxp1,mndx,n) &
1892 + w(3)*xsqy_tab(1:nwave,tndx,mndxp1,n) &
1893 + w(4)*xsqy_tab(1:nwave,tndxp1,mndxp1,n)
1894 end do
1896 end subroutine xsqy_int
1898 subroutine xs_int( xs, tlev, xs_tab )
1899 !---------------------------------------------------------------------
1900 ! ... interpolate tables for xs
1901 !---------------------------------------------------------------------
1903 real, intent(in) :: tlev(:)
1904 real, intent(in) :: xs_tab(:,:)
1905 real, intent(out) :: xs(:,:)
1907 integer :: tndx, tndxp1
1908 integer :: k, ku
1909 real :: temp
1910 real :: w(2)
1911 real :: del_t
1913 ku = size( tlev )
1914 do k = 1,ku
1915 temp = tlev(k)
1916 do tndx = 1,ntemp
1917 if( temp_tab(tndx) > temp ) then
1918 exit
1919 endif
1920 end do
1921 tndx = max( min( tndx,ntemp ) - 1,1 )
1922 tndxp1 = tndx + 1
1923 del_t = max( 0.,min( 1.,(temp - temp_tab(tndx))*del_temp_tab(tndx) ) )
1925 w(1) = (1. - del_t)
1926 w(2) = del_t
1928 xs(1:nwave,k) = w(1)*xs_tab(1:nwave,tndx) &
1929 + w(2)*xs_tab(1:nwave,tndxp1)
1930 end do
1932 end subroutine xs_int
1934 FUNCTION chap(zeta)
1935 ! chapman function is used when the solar zenith angle exceeds
1936 ! 75 deg.
1937 ! interpolates between values given in, e.g., mccartney (1976).
1938 ! .. Scalar Arguments ..
1939 REAL, intent(in) :: zeta
1940 ! .. Local Scalars ..
1941 REAL :: rm
1942 INTEGER :: i
1943 LOGICAL :: fnd
1944 ! .. Local Arrays ..
1945 REAL :: y(75:96)
1946 ! .. Function Return Value ..
1947 REAL :: chap
1948 ! .. Data Statements ..
1949 DATA (y(i),i=75,96) &
1950 /3.800, 4.055, 4.348, 4.687, 5.083, 5.551, 6.113, &
1951 6.799, 7.650, 8.732, 10.144, 12.051, 14.730, 18.686, 24.905, 35.466, &
1952 55.211, 96.753, 197., 485., 1476., 9999./
1954 fnd = .false.
1955 DO i = 75, 96
1956 rm = real(i)
1957 IF (zeta < rm) then
1958 chap = y(i) + (y(i+1) - y(i))*(zeta - (rm - 1.))
1959 fnd = .true.
1960 exit
1961 ENDIF
1962 END DO
1964 IF( .not. fnd ) then
1965 chap = y(96)
1966 ENDIF
1968 END FUNCTION chap
1970 SUBROUTINE fery( nwave, w, wght )
1971 !-----------------------------------------------------------------------------
1972 != PURPOSE:
1973 != Calculate the action spectrum value for erythema at a given wavelength
1974 != according to: McKinlay, A.F and B.L.Diffey, A reference action spectrum
1975 != for ultraviolet induced erythema in human skin, CIE Journal, vol 6,
1976 != pp 17-22, 1987.
1977 != Value at 300 nm = 0.6486
1978 !-----------------------------------------------------------------------------
1980 INTEGER, intent(in) :: nwave
1981 REAL, intent(in) :: w(:)
1982 REAL, intent(out) :: wght(:)
1984 WHERE( w(1:nwave) < 298. )
1985 wght(1:nwave) = 1.
1986 ELSEWHERE( w(1:nwave) >= 298. .AND. w(1:nwave) < 328. )
1987 wght(1:nwave) = 10.**(0.094*(298. - w(1:nwave)))
1988 ELSEWHERE( w(1:nwave) >= 328. .AND. w(1:nwave) < 400. )
1989 wght(1:nwave) = 10.**(0.015*(139. - w(1:nwave)))
1990 ELSEWHERE( w(1:nwave) >= 400. )
1991 wght(1:nwave) = 1.e-36
1992 ENDWHERE
1994 END SUBROUTINE fery
1996 SUBROUTINE cldfrac_binary( CLDFRA,QC,QI, QS, kts, kte )
1997 !---------------------------------------------------------------------
1998 ! !DESCRIPTION:
1999 ! Compute cloud fraction from input ice, snow, and cloud water fields
2000 ! if provided.
2001 !
2002 ! Whether QI or QC is active or not is determined from the indices of
2003 ! the fields into the 4D scalar arrays in WRF. These indices are
2004 ! P_QI and P_QC, respectively, and they are passed in to the routine
2005 ! to enable testing to see if QI and QC represent active fields in
2006 ! the moisture 4D scalar array carried by WRF.
2007 !
2008 ! If a field is active its index will have a value greater than or
2009 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
2010 ! this routine.
2011 !EOP
2012 !---------------------------------------------------------------------
2014 IMPLICIT NONE
2016 INTEGER, INTENT(IN) :: kts, kte
2018 REAL, INTENT(OUT ) :: CLDFRA(kts:kte)
2020 REAL, INTENT(IN) :: QI(kts:kte), &
2021 QC(kts:kte), &
2022 QS(kts:kte)
2023 !----------------------------------------------------------------------
2024 ! local variables
2025 !----------------------------------------------------------------------
2026 REAL, parameter :: thresh = 1.e-9
2028 INTEGER :: j
2030 where( (qc(kts:kte) + qi(kts:kte) + qs(kts:kte)) > thresh )
2031 cldfra(kts:kte) = 1.
2032 elsewhere
2033 cldfra(kts:kte) = 0.
2034 endwhere
2036 END SUBROUTINE cldfrac_binary
2038 SUBROUTINE cldfrac_fractional( CLDFRA, QV, QC, QI, QS, &
2039 p_phy, t_phy, &
2040 kts, kte )
2042 !----------------------------------------------------------------------
2043 ! dummy arguments
2044 !----------------------------------------------------------------------
2045 INTEGER, INTENT(IN) :: kts,kte
2047 REAL, INTENT(OUT) :: CLDFRA(kts:kte)
2049 REAL, INTENT(IN) :: QV(kts:kte), &
2050 QI(kts:kte), &
2051 QC(kts:kte), &
2052 QS(kts:kte), &
2053 t_phy(kts:kte), &
2054 p_phy(kts:kte)
2057 !----------------------------------------------------------------------
2058 ! local variables
2059 !----------------------------------------------------------------------
2060 REAL , PARAMETER :: ALPHA0 = 100.
2061 REAL , PARAMETER :: GAMMA = 0.49
2062 REAL , PARAMETER :: QCLDMIN = 1.E-12
2063 REAL , PARAMETER :: PEXP = 0.25
2064 REAL , PARAMETER :: RHGRID =1.0
2065 REAL , PARAMETER :: SVP1 = 0.61078
2066 REAL , PARAMETER :: SVP2 = 17.2693882
2067 REAL , PARAMETER :: SVPI2 = 21.8745584
2068 REAL , PARAMETER :: SVP3 = 35.86
2069 REAL , PARAMETER :: SVPI3 = 7.66
2070 REAL , PARAMETER :: SVPT0 = 273.15
2071 REAL , PARAMETER :: r_d = 287.
2072 REAL , PARAMETER :: r_v = 461.6
2073 REAL , PARAMETER :: ep_2 = r_d/r_v
2075 INTEGER :: i,j,k
2076 INTEGER :: imax, jmax, kmax
2077 REAL :: RHUM, tc, esw, esi, weight, qvsw, qvsi, qvs_weight
2078 REAL :: QCLD, DENOM, ARG, SUBSAT, wrk
2079 REAL :: relhum_max, wrk_max
2082 ! !DESCRIPTION:
2083 !----------------------------------------------------------------------
2084 ! Compute cloud fraction from input ice and cloud water fields
2085 ! if provided.
2086 !
2087 ! Whether QI or QC is active or not is determined from the indices of
2088 ! the fields into the 4D scalar arrays in WRF. These indices are
2089 ! P_QI and P_QC, respectively, and they are passed in to the routine
2090 ! to enable testing to see if QI and QC represent active fields in
2091 ! the moisture 4D scalar array carried by WRF.
2092 !
2093 ! If a field is active its index will have a value greater than or
2094 ! equal to PARAM_FIRST_SCALAR, which is also an input argument to
2095 ! this routine.
2096 !----------------------------------------------------------------------
2097 !EOP
2100 !-----------------------------------------------------------------------
2101 ! COMPUTE GRID-SCALE CLOUD COVER FOR RADIATION
2102 ! (modified by Ferrier, Feb '02)
2103 !
2104 ! Cloud fraction parameterization follows Randall, 1994
2105 ! (see Hong et al., 1998)
2106 !-----------------------------------------------------------------------
2107 ! Note: ep_2=287./461.6 Rd/Rv
2108 ! Note: R_D=287.
2110 ! Alternative calculation for critical RH for grid saturation
2111 ! RHGRID=0.90+.08*((100.-DX)/95.)**.5
2113 ! Calculate saturation mixing ratio weighted according to the fractions of
2114 ! water and ice.
2115 ! Following:
2116 ! Murray, F.W. 1966. ``On the computation of Saturation Vapor Pressure'' J. Appl. Meteor. 6 p.204
2117 ! es (in mb) = 6.1078 . exp[ a . (T-273.16)/ (T-b) ]
2118 !
2119 ! over ice over water
2120 ! a = 21.8745584 17.2693882
2121 ! b = 7.66 35.86
2122 !---------------------------------------------------------------------
2124 relhum_max = -100.
2125 wrk_max = -10000.
2126 imax = 0; kmax = 0; jmax = 0
2128 vert_loop: &
2129 DO k = kts,kte
2130 !---------------------------------------------------------------------
2131 ! Determine cloud fraction (modified from original algorithm)
2132 !---------------------------------------------------------------------
2133 QCLD = QI(k) + QC(k) + QS(k)
2134 has_cloud : &
2135 IF( QCLD >= QCLDMIN ) THEN
2136 tc = t_phy(k) - SVPT0
2137 esw = 1000.0 * SVP1 * EXP( SVP2 * tc / ( t_phy(k) - SVP3 ) )
2138 esi = 1000.0 * SVP1 * EXP( SVPI2 * tc / ( t_phy(k) - SVPI3 ) )
2139 QVSW = EP_2 * esw / ( p_phy(k) - esw )
2140 QVSI = EP_2 * esi / ( p_phy(k) - esi )
2142 weight = (QI(k) + QS(k)) / QCLD
2143 QVS_WEIGHT = (1. - weight)*QVSW + weight*QVSI
2144 RHUM = QV(k)/QVS_WEIGHT !--- Relative humidity
2145 !---------------------------------------------------------------------
2146 ! Assume zero cloud fraction if there is no cloud mixing ratio
2147 !---------------------------------------------------------------------
2148 IF( RHUM >= RHGRID )THEN
2149 !---------------------------------------------------------------------
2150 ! Assume cloud fraction of unity if near saturation and the cloud
2151 ! mixing ratio is at or above the minimum threshold
2152 !---------------------------------------------------------------------
2153 CLDFRA(k) = 1.
2154 ELSE
2155 !---------------------------------------------------------------------
2156 ! Adaptation of original algorithm (Randall, 1994; Zhao, 1995)
2157 ! modified based on assumed grid-scale saturation at RH=RHgrid.
2158 !---------------------------------------------------------------------
2159 SUBSAT = MAX( 1.E-10,RHGRID*QVS_WEIGHT - QV(k) )
2160 DENOM = SUBSAT**GAMMA
2161 ARG = MAX( -6.9,-ALPHA0*QCLD/DENOM ) ! <-- EXP(-6.9)=.001
2162 !---------------------------------------------------------------------
2163 ! prevent negative values (new)
2164 !---------------------------------------------------------------------
2165 RHUM = MAX( 1.E-10, RHUM )
2166 wrk = (RHUM/RHGRID)**PEXP*(1. - EXP( ARG ))
2167 IF( wrk >= .01 ) then
2168 CLDFRA(k) = wrk
2169 ENDIF
2170 ENDIF
2171 ENDIF has_cloud
2172 END DO vert_loop
2174 END SUBROUTINE cldfrac_fractional
2176 #ifdef NETCDF
2177 subroutine handle_ncerr( ret, mes )
2178 !---------------------------------------------------------------------
2179 ! ... netcdf error handling routine
2180 !---------------------------------------------------------------------
2182 implicit none
2184 !---------------------------------------------------------------------
2185 ! ... dummy arguments
2186 !---------------------------------------------------------------------
2187 integer, intent(in) :: ret
2188 character(len=*), intent(in) :: mes
2190 include 'netcdf.inc'
2192 if( ret /= nf_noerr ) then
2193 call wrf_message( trim(mes) )
2194 call wrf_message( trim(nf_strerror(ret)) )
2195 call wrf_abort
2196 end if
2198 end subroutine handle_ncerr
2199 #endif
2201 end module module_phot_tuv