Line data Source code
1 : module mo_photo
2 : !----------------------------------------------------------------------
3 : ! ... photolysis interp table and related arrays
4 : !----------------------------------------------------------------------
5 :
6 : use shr_kind_mod, only : r8 => shr_kind_r8
7 : use ppgrid, only : pcols, pver, begchunk, endchunk
8 : use cam_abortutils, only : endrun
9 : use mo_constants, only : r2d,d2r
10 : use ref_pres, only : num_pr_lev, ptop_ref
11 : use pio
12 : use cam_pio_utils, only : cam_pio_openfile
13 : use spmd_utils, only : masterproc
14 : use cam_logfile, only : iulog
15 : use solar_parms_data, only : f107=>solar_parms_f107, f107a=>solar_parms_f107a
16 :
17 : implicit none
18 :
19 : private
20 :
21 : public :: photo_inti, table_photo
22 : public :: set_ub_col
23 : public :: setcol
24 : public :: photo_timestep_init
25 : public :: photo_register
26 :
27 : save
28 :
29 : real(r8), parameter :: kg2g = 1.e3_r8
30 : integer, parameter :: pverm = pver - 1
31 :
32 : integer :: jno_ndx
33 : integer :: jonitr_ndx
34 : integer :: jho2no2_ndx
35 : integer :: jo2_a_ndx, jo2_b_ndx
36 : integer :: ox_ndx, o3_ndx, o3_inv_ndx, o3rad_ndx
37 : integer :: n2_ndx, no_ndx, o2_ndx, o_ndx
38 : integer, allocatable :: lng_indexer(:)
39 : integer, allocatable :: sht_indexer(:)
40 : integer, allocatable :: euv_indexer(:)
41 :
42 : integer :: ki
43 : integer :: last
44 : integer :: next
45 : integer :: n_exo_levs
46 : real(r8) :: delp
47 : real(r8) :: dels
48 : real(r8), allocatable :: days(:)
49 : real(r8), allocatable :: levs(:)
50 : real(r8), allocatable :: o2_exo_coldens(:,:,:,:)
51 : real(r8), allocatable :: o3_exo_coldens(:,:,:,:)
52 : logical :: o_is_inv
53 : logical :: o2_is_inv
54 : logical :: n2_is_inv
55 : logical :: o3_is_inv
56 : logical :: no_is_inv
57 : logical :: has_o2_col
58 : logical :: has_o3_col
59 : logical :: has_fixed_press
60 : real(r8) :: max_zen_angle ! degrees
61 :
62 : integer :: jo1d_ndx, jo3p_ndx, jno2_ndx, jn2o5_ndx
63 : integer :: jhno3_ndx, jno3_ndx, jpan_ndx, jmpan_ndx
64 :
65 : integer :: jo1da_ndx, jo3pa_ndx, jno2a_ndx, jn2o5a_ndx, jn2o5b_ndx
66 : integer :: jhno3a_ndx, jno3a_ndx, jpana_ndx, jmpana_ndx, jho2no2a_ndx
67 : integer :: jonitra_ndx
68 :
69 : integer :: jppi_ndx, jepn1_ndx, jepn2_ndx, jepn3_ndx, jepn4_ndx, jepn6_ndx
70 : integer :: jepn7_ndx, jpni1_ndx, jpni2_ndx, jpni3_ndx, jpni4_ndx, jpni5_ndx
71 : logical :: do_jeuv = .false.
72 : logical :: do_jshort = .false.
73 : integer :: ion_rates_idx = -1
74 :
75 : contains
76 :
77 :
78 : !----------------------------------------------------------------------
79 : !----------------------------------------------------------------------
80 0 : subroutine photo_register
81 : use mo_jeuv, only : nIonRates
82 : use physics_buffer,only : pbuf_add_field, dtype_r8
83 :
84 : ! add photo-ionization rates to phys buffer for waccmx ionosphere module
85 :
86 0 : call pbuf_add_field('IonRates' , 'physpkg', dtype_r8, (/pcols,pver,nIonRates/), ion_rates_idx) ! Ionization rates for O+,O2+,N+,N2+,NO+
87 :
88 0 : endsubroutine photo_register
89 :
90 : !----------------------------------------------------------------------
91 : !----------------------------------------------------------------------
92 0 : subroutine photo_inti( xs_coef_file, xs_short_file, xs_long_file, rsf_file, &
93 : photon_file, electron_file, &
94 : exo_coldens_file, maxzen )
95 : !----------------------------------------------------------------------
96 : ! ... initialize photolysis module
97 : !----------------------------------------------------------------------
98 :
99 0 : use interpolate_data, only: lininterp_init, lininterp, lininterp_finish, interp_type
100 : use chem_mods, only : phtcnt
101 : use chem_mods, only : ncol_abs => nabscol
102 : use chem_mods, only : rxt_tag_lst, pht_alias_lst, pht_alias_mult
103 : use time_manager, only : get_calday
104 : use ioFileMod, only : getfil
105 : use mo_chem_utls, only : get_spc_ndx, get_rxt_ndx, get_inv_ndx
106 : use mo_jlong, only : jlong_init
107 : use mo_jshort, only : jshort_init
108 : use mo_jeuv, only : jeuv_init, neuv
109 : use phys_grid, only : get_ncols_p, get_rlat_all_p
110 : use solar_irrad_data,only : has_spectrum
111 : use photo_bkgrnd, only : photo_bkgrnd_init
112 : use cam_history, only : addfld
113 :
114 : implicit none
115 :
116 : !----------------------------------------------------------------------
117 : ! ... dummy arguments
118 : !----------------------------------------------------------------------
119 : character(len=*), intent(in) :: xs_long_file, rsf_file
120 : character(len=*), intent(in) :: exo_coldens_file
121 : real(r8), intent(in) :: maxzen
122 : ! waccm
123 : character(len=*), intent(in) :: xs_coef_file
124 : character(len=*), intent(in) :: xs_short_file
125 : character(len=*), intent(in) :: photon_file
126 : character(len=*), intent(in) :: electron_file
127 :
128 : !----------------------------------------------------------------------
129 : ! ... local variables
130 : !----------------------------------------------------------------------
131 : real(r8), parameter :: hPa2Pa = 100._r8
132 : integer :: k, n
133 : type(file_desc_t) :: ncid
134 : type(var_desc_t) :: vid
135 : type(interp_type) :: lat_wgts
136 : integer :: dimid
137 : integer :: nlat
138 : integer :: ntimes
139 : integer :: astat
140 : integer :: ndx
141 : integer :: spc_ndx
142 : integer :: ierr
143 : integer :: c, ncols
144 0 : integer, allocatable :: dates(:)
145 : real(r8) :: pinterp
146 0 : real(r8), allocatable :: lats(:)
147 0 : real(r8), allocatable :: coldens(:,:,:)
148 : character(len=256) :: locfn
149 : character(len=256) :: filespec
150 : real(r8), parameter :: trop_thrshld = 1._r8 ! Pa
151 : real(r8) :: to_lats(pcols)
152 :
153 :
154 : if( phtcnt < 1 ) then
155 : return
156 : end if
157 :
158 : ! maximum solar zenith angle for which photo-chemical rates are computed
159 : max_zen_angle = maxzen
160 : if (.not. max_zen_angle>0._r8) then
161 : write(iulog,*) 'photo_inti: max_zen_angle = ',max_zen_angle
162 : call endrun ('photo_inti: max_zen_angle must be greater then zero')
163 : end if
164 :
165 :
166 : ! jeuv_1,,, jeuv_25 --> need euv calculations --> must be waccm
167 : ! how to determine if shrt calc is needed ?? -- use top level pressure => waccm = true ? false
168 :
169 : if ( .not. has_spectrum ) then
170 : write(iulog,*) 'photo_inti: solar_irrad_data file needs to contain irradiance spectrum'
171 : call endrun('photo_inti: ERROR -- solar irradiance spectrum is missing')
172 : endif
173 :
174 : !----------------------------------------------------------------------
175 : ! ... allocate indexers
176 : !----------------------------------------------------------------------
177 : allocate( lng_indexer(phtcnt),stat=astat )
178 : if( astat /= 0 ) then
179 : write(iulog,*) 'photo_inti: lng_indexer allocation error = ',astat
180 : call endrun
181 : end if
182 : lng_indexer(:) = 0
183 : allocate( sht_indexer(phtcnt),stat=astat )
184 : if( astat /= 0 ) then
185 : write(iulog,*) 'photo_inti: Failed to allocate sht_indexer; error = ',astat
186 : call endrun
187 : end if
188 : sht_indexer(:) = 0
189 : allocate( euv_indexer(neuv),stat=astat )
190 : if( astat /= 0 ) then
191 : write(iulog,*) 'photo_inti: Failed to allocate euv_indexer; error = ',astat
192 : call endrun
193 : end if
194 : euv_indexer(:) = 0
195 :
196 : jno_ndx = get_rxt_ndx( 'jno' )
197 : jo2_a_ndx = get_rxt_ndx( 'jo2_a' )
198 : jo2_b_ndx = get_rxt_ndx( 'jo2_b' )
199 :
200 : jo1da_ndx = get_rxt_ndx( 'jo1da' )
201 : jo3pa_ndx = get_rxt_ndx( 'jo3pa' )
202 : jno2a_ndx = get_rxt_ndx( 'jno2a' )
203 : jn2o5a_ndx = get_rxt_ndx( 'jn2o5a' )
204 : jn2o5b_ndx = get_rxt_ndx( 'jn2o5b' )
205 : jhno3a_ndx = get_rxt_ndx( 'jhno3a' )
206 : jno3a_ndx = get_rxt_ndx( 'jno3a' )
207 : jpana_ndx = get_rxt_ndx( 'jpana' )
208 : jmpana_ndx = get_rxt_ndx( 'jmpana' )
209 : jho2no2a_ndx = get_rxt_ndx( 'jho2no2a' )
210 : jonitra_ndx = get_rxt_ndx( 'jonitra' )
211 :
212 : jo1d_ndx = get_rxt_ndx( 'jo1d' )
213 : jo3p_ndx = get_rxt_ndx( 'jo3p' )
214 : jno2_ndx = get_rxt_ndx( 'jno2' )
215 : jn2o5_ndx = get_rxt_ndx( 'jn2o5' )
216 : jn2o5_ndx = get_rxt_ndx( 'jn2o5' )
217 : jhno3_ndx = get_rxt_ndx( 'jhno3' )
218 : jno3_ndx = get_rxt_ndx( 'jno3' )
219 : jpan_ndx = get_rxt_ndx( 'jpan' )
220 : jmpan_ndx = get_rxt_ndx( 'jmpan' )
221 : jho2no2_ndx = get_rxt_ndx( 'jho2no2' )
222 : jonitr_ndx = get_rxt_ndx( 'jonitr' )
223 :
224 : jppi_ndx = get_rxt_ndx( 'jppi' )
225 : jepn1_ndx = get_rxt_ndx( 'jepn1' )
226 : jepn2_ndx = get_rxt_ndx( 'jepn2' )
227 : jepn3_ndx = get_rxt_ndx( 'jepn3' )
228 : jepn4_ndx = get_rxt_ndx( 'jepn4' )
229 : jepn6_ndx = get_rxt_ndx( 'jepn6' )
230 : jepn7_ndx = get_rxt_ndx( 'jepn7' )
231 : jpni1_ndx = get_rxt_ndx( 'jpni1' )
232 : jpni2_ndx = get_rxt_ndx( 'jpni2' )
233 : jpni3_ndx = get_rxt_ndx( 'jpni3' )
234 : jpni4_ndx = get_rxt_ndx( 'jpni4' )
235 : ! added to v02
236 : jpni5_ndx = get_rxt_ndx( 'jpni5' )
237 : ox_ndx = get_spc_ndx( 'OX' )
238 : if( ox_ndx < 1 ) then
239 : ox_ndx = get_spc_ndx( 'O3' )
240 : end if
241 : o3_ndx = get_spc_ndx( 'O3' )
242 : o3rad_ndx = get_spc_ndx( 'O3RAD' )
243 : o3_inv_ndx = get_inv_ndx( 'O3' )
244 :
245 : n2_ndx = get_inv_ndx( 'N2' )
246 : n2_is_inv = n2_ndx > 0
247 : if( .not. n2_is_inv ) then
248 : n2_ndx = get_spc_ndx( 'N2' )
249 : end if
250 : o2_ndx = get_inv_ndx( 'O2' )
251 : o2_is_inv = o2_ndx > 0
252 : if( .not. o2_is_inv ) then
253 : o2_ndx = get_spc_ndx( 'O2' )
254 : end if
255 : no_ndx = get_spc_ndx( 'NO' )
256 : no_is_inv = no_ndx < 1
257 : if( no_is_inv ) then
258 : no_ndx = get_inv_ndx( 'NO' )
259 : end if
260 : o3_is_inv = o3_ndx < 1
261 :
262 : o_ndx = get_spc_ndx( 'O' )
263 : o_is_inv = o_ndx < 1
264 : if( o_is_inv ) then
265 : o_ndx = get_inv_ndx( 'O' )
266 : end if
267 :
268 : do_jshort = o_ndx>0 .and. o2_ndx>0 .and. (o3_ndx>0.or.o3_inv_ndx>0) .and. n2_ndx>0 .and. no_ndx>0
269 :
270 : call jeuv_init( photon_file, electron_file, euv_indexer )
271 : do_jeuv = any(euv_indexer(:)>0)
272 :
273 : !----------------------------------------------------------------------
274 : ! ... call module initializers
275 : !----------------------------------------------------------------------
276 : call jlong_init( xs_long_file, rsf_file, lng_indexer )
277 : if (do_jeuv) then
278 : call photo_bkgrnd_init()
279 : call addfld('Qbkgndtot', (/ 'lev' /), 'A','cm-3 sec-1', 'background ionization rate ' )
280 : call addfld('Qbkgnd_o1', (/ 'lev' /), 'A','cm-3 sec-1', 'background ionization rate ' )
281 : call addfld('Qbkgnd_o2', (/ 'lev' /), 'A','cm-3 sec-1', 'background ionization rate ' )
282 : call addfld('Qbkgnd_n2', (/ 'lev' /), 'A','cm-3 sec-1', 'background ionization rate ' )
283 : call addfld('Qbkgnd_n1', (/ 'lev' /), 'A','cm-3 sec-1', 'background ionization rate ' )
284 : call addfld('Qbkgnd_no', (/ 'lev' /), 'A','cm-3 sec-1', 'background ionization rate ' )
285 : endif
286 : if (do_jshort) then
287 : call jshort_init( xs_coef_file, xs_short_file, sht_indexer )
288 : endif
289 : jho2no2_ndx = get_rxt_ndx( 'jho2no2_b' )
290 :
291 : !----------------------------------------------------------------------
292 : ! ... check that each photorate is in short or long datasets
293 : !----------------------------------------------------------------------
294 : if( any( ( abs(sht_indexer(:)) + abs(lng_indexer(:)) ) == 0 ) ) then
295 : write(iulog,*) ' '
296 : write(iulog,*) 'photo_inti: the following photorate(s) are not in'
297 : write(iulog,*) ' either the short or long datasets'
298 : write(iulog,*) ' '
299 : do ndx = 1,phtcnt
300 : if( abs(sht_indexer(ndx)) + abs(lng_indexer(ndx)) == 0 ) then
301 : write(iulog,*) ' ',trim( rxt_tag_lst(ndx) )
302 : end if
303 : end do
304 : call endrun
305 : end if
306 :
307 : !----------------------------------------------------------------------
308 : ! ... output any aliased photorates
309 : !----------------------------------------------------------------------
310 : if( masterproc ) then
311 : if( any( pht_alias_lst(:,1) /= ' ' ) ) then
312 : write(iulog,*) ' '
313 : write(iulog,*) 'photo_inti: the following short photorate(s) are aliased'
314 : write(iulog,*) ' '
315 : do ndx = 1,phtcnt
316 : if( pht_alias_lst(ndx,1) /= ' ' ) then
317 : if( pht_alias_mult(ndx,1) == 1._r8 ) then
318 : write(iulog,*) ' ',trim(rxt_tag_lst(ndx)),' -> ',trim(pht_alias_lst(ndx,1))
319 : else
320 : write(iulog,*) ' ',trim(rxt_tag_lst(ndx)),' -> ',pht_alias_mult(ndx,1),'*',trim(pht_alias_lst(ndx,1))
321 : end if
322 : end if
323 : end do
324 : end if
325 : if( any( pht_alias_lst(:,2) /= ' ' ) ) then
326 : write(iulog,*) ' '
327 : write(iulog,*) 'photo_inti: the following long photorate(s) are aliased'
328 : write(iulog,*) ' '
329 : do ndx = 1,phtcnt
330 : if( pht_alias_lst(ndx,2) /= ' ' ) then
331 : if( pht_alias_mult(ndx,2) == 1._r8 ) then
332 : write(iulog,*) ' ',trim(rxt_tag_lst(ndx)),' -> ',trim(pht_alias_lst(ndx,2))
333 : else
334 : write(iulog,*) ' ',trim(rxt_tag_lst(ndx)),' -> ',pht_alias_mult(ndx,2),'*',trim(pht_alias_lst(ndx,2))
335 : end if
336 : end if
337 : end do
338 : end if
339 :
340 : write(iulog,*) ' '
341 : write(iulog,*) '*********************************************'
342 : write(iulog,*) 'photo_inti: euv_indexer'
343 : write(iulog,'(10i6)') euv_indexer(:)
344 : write(iulog,*) 'photo_inti: sht_indexer'
345 : write(iulog,'(10i6)') sht_indexer(:)
346 : write(iulog,*) 'photo_inti: lng_indexer'
347 : write(iulog,'(10i6)') lng_indexer(:)
348 : write(iulog,*) '*********************************************'
349 : write(iulog,*) ' '
350 : endif
351 :
352 : !----------------------------------------------------------------------
353 : ! ... check for o2, o3 absorber columns
354 : !----------------------------------------------------------------------
355 : if( ncol_abs > 0 ) then
356 : spc_ndx = ox_ndx
357 : if( spc_ndx < 1 ) then
358 : spc_ndx = o3_ndx
359 : end if
360 : if( spc_ndx > 0 ) then
361 : has_o3_col = .true.
362 : else
363 : has_o3_col = .false.
364 : end if
365 : if( ncol_abs > 1 ) then
366 : if( o2_ndx > 1 ) then
367 : has_o2_col = .true.
368 : else
369 : has_o2_col = .false.
370 : end if
371 : else
372 : has_o2_col = .false.
373 : end if
374 : else
375 : has_o2_col = .false.
376 : has_o3_col = .false.
377 : end if
378 :
379 : if ( len_trim(exo_coldens_file) == 0 ) then
380 : has_o2_col = .false.
381 : has_o3_col = .false.
382 : endif
383 :
384 : has_abs_columns : if( has_o2_col .or. has_o3_col ) then
385 : !-----------------------------------------------------------------------
386 : ! ... open exo coldens file
387 : !-----------------------------------------------------------------------
388 : filespec = trim( exo_coldens_file )
389 : call getfil( filespec, locfn, 0 )
390 : call cam_pio_openfile( ncid, trim(locfn), PIO_NOWRITE )
391 :
392 : !-----------------------------------------------------------------------
393 : ! ... get grid dimensions from file
394 : !-----------------------------------------------------------------------
395 : ! ... timing
396 : !-----------------------------------------------------------------------
397 : ierr = pio_inq_dimid( ncid, 'month', dimid )
398 : ierr = pio_inq_dimlen( ncid, dimid, ntimes )
399 :
400 : if( ntimes /= 12 ) then
401 : call endrun('photo_inti: exo coldens is not annual period')
402 : end if
403 : allocate( dates(ntimes),days(ntimes),stat=astat )
404 : if( astat /= 0 ) then
405 : write(iulog,*) 'photo_inti: dates,days allocation error = ',astat
406 : call endrun
407 : end if
408 : dates(:) = (/ 116, 214, 316, 415, 516, 615, &
409 : 716, 816, 915, 1016, 1115, 1216 /)
410 : !-----------------------------------------------------------------------
411 : ! ... initialize the monthly day of year times
412 : !-----------------------------------------------------------------------
413 : do n = 1,ntimes
414 : days(n) = get_calday( dates(n), 0 )
415 : end do
416 : deallocate( dates )
417 : !-----------------------------------------------------------------------
418 : ! ... latitudes
419 : !-----------------------------------------------------------------------
420 : ierr = pio_inq_dimid( ncid, 'lat', dimid )
421 : ierr = pio_inq_dimlen( ncid, dimid, nlat )
422 : allocate( lats(nlat), stat=astat )
423 : if( astat /= 0 ) then
424 : write(iulog,*) 'photo_inti: lats allocation error = ',astat
425 : call endrun
426 : end if
427 : ierr = pio_inq_varid( ncid, 'lat', vid )
428 : ierr = pio_get_var( ncid, vid, lats )
429 : lats(:nlat) = lats(:nlat) * d2r
430 : !-----------------------------------------------------------------------
431 : ! ... levels
432 : !-----------------------------------------------------------------------
433 : ierr = pio_inq_dimid( ncid, 'lev', dimid )
434 : ierr = pio_inq_dimlen( ncid, dimid, n_exo_levs )
435 : allocate( levs(n_exo_levs), stat=astat )
436 : if( astat /= 0 ) then
437 : write(iulog,*) 'photo_inti: levs allocation error = ',astat
438 : call endrun
439 : end if
440 : ierr = pio_inq_varid( ncid, 'lev', vid )
441 : ierr = pio_get_var( ncid, vid, levs )
442 : levs(:n_exo_levs) = levs(:n_exo_levs) * hPa2Pa
443 : !-----------------------------------------------------------------------
444 : ! ... set up regridding
445 : !-----------------------------------------------------------------------
446 :
447 : allocate( coldens(nlat,n_exo_levs,ntimes),stat=astat )
448 : if( astat /= 0 ) then
449 : write(iulog,*) 'photo_inti: coldens allocation error = ',astat
450 : call endrun
451 : end if
452 : if( has_o2_col ) then
453 : allocate( o2_exo_coldens(n_exo_levs,pcols,begchunk:endchunk,ntimes),stat=astat )
454 : if( astat /= 0 ) then
455 : write(iulog,*) 'photo_inti: o2_exo_coldens allocation error = ',astat
456 : call endrun
457 : end if
458 : ierr = pio_inq_varid( ncid, 'O2_column_density', vid )
459 : ierr = pio_get_var( ncid, vid,coldens )
460 :
461 : do c=begchunk,endchunk
462 : ncols = get_ncols_p(c)
463 : call get_rlat_all_p(c, pcols, to_lats)
464 : call lininterp_init(lats, nlat, to_lats, ncols, 1, lat_wgts)
465 : do n=1,ntimes
466 : do k = 1,n_exo_levs
467 : call lininterp(coldens(:,k,n), nlat, o2_exo_coldens(k,:,c,n), ncols, lat_wgts)
468 : end do
469 : end do
470 : call lininterp_finish(lat_wgts)
471 : enddo
472 :
473 :
474 : end if
475 : if( has_o3_col ) then
476 : allocate( o3_exo_coldens(n_exo_levs,pcols,begchunk:endchunk,ntimes),stat=astat )
477 : if( astat /= 0 ) then
478 : write(iulog,*) 'photo_inti: o3_exo_coldens allocation error = ',astat
479 : call endrun
480 : end if
481 : ierr = pio_inq_varid( ncid, 'O3_column_density', vid )
482 : ierr = pio_get_var( ncid, vid,coldens )
483 :
484 : do c=begchunk,endchunk
485 : ncols = get_ncols_p(c)
486 : call get_rlat_all_p(c, pcols, to_lats)
487 : call lininterp_init(lats, nlat, to_lats, ncols, 1, lat_wgts)
488 : do n=1,ntimes
489 : do k = 1,n_exo_levs
490 : call lininterp(coldens(:,k,n), nlat, o3_exo_coldens(k,:,c,n), ncols, lat_wgts)
491 : end do
492 : end do
493 : call lininterp_finish(lat_wgts)
494 : enddo
495 : end if
496 : call pio_closefile (ncid)
497 : deallocate( coldens,stat=astat )
498 : if( astat /= 0 ) then
499 : write(iulog,*) 'photo_inti: failed to deallocate coldens; error = ',astat
500 : call endrun
501 : end if
502 : has_fixed_press = (num_pr_lev .ne. 0)
503 : !-----------------------------------------------------------------------
504 : ! ... setup the pressure interpolation
505 : !-----------------------------------------------------------------------
506 : if( has_fixed_press ) then
507 : pinterp = ptop_ref
508 : if( pinterp <= levs(1) ) then
509 : ki = 1
510 : delp = 0._r8
511 : else
512 : do ki = 2,n_exo_levs
513 : if( pinterp <= levs(ki) ) then
514 : delp = log( pinterp/levs(ki-1) )/log( levs(ki)/levs(ki-1) )
515 : exit
516 : end if
517 : end do
518 : end if
519 : #ifdef DEBUG
520 : if (masterproc) then
521 : write(iulog,*) '-----------------------------------'
522 : write(iulog,*) 'photo_inti: diagnostics'
523 : write(iulog,*) 'ki, delp = ',ki,delp
524 : if (ki>1) then
525 : write(iulog,*) 'pinterp,levs(ki-1:ki) = ',pinterp,levs(ki-1:ki)
526 : else
527 : write(iulog,*) 'pinterp,levs(ki) = ',pinterp,levs(ki)
528 : end if
529 : write(iulog,*) '-----------------------------------'
530 : endif
531 : #endif
532 : end if
533 : end if has_abs_columns
534 :
535 0 : end subroutine photo_inti
536 :
537 0 : subroutine table_photo( photos, pmid, pdel, temper, zmid, zint, &
538 0 : col_dens, zen_angle, srf_alb, lwc, clouds, &
539 0 : esfact, vmr, invariants, ncol, lchnk, pbuf )
540 : !-----------------------------------------------------------------
541 : ! ... table photorates for wavelengths > 200nm
542 : !-----------------------------------------------------------------
543 :
544 : use chem_mods, only : ncol_abs => nabscol, phtcnt, gas_pcnst, nfs
545 0 : use chem_mods, only : pht_alias_mult, indexm
546 : use mo_jshort, only : nsht => nj, jshort
547 : use mo_jlong, only : nlng => numj, jlong
548 : use mo_jeuv, only : neuv, jeuv, nIonRates
549 : use physics_buffer, only : physics_buffer_desc, pbuf_get_field
550 : use photo_bkgrnd, only : photo_bkgrnd_calc
551 : use cam_history, only : outfld
552 : use infnan, only : nan, assignment(=)
553 :
554 : implicit none
555 :
556 : !-----------------------------------------------------------------
557 : ! ... dummy arguments
558 : !-----------------------------------------------------------------
559 : integer, intent(in) :: lchnk
560 : integer, intent(in) :: ncol
561 : real(r8), intent(in) :: esfact ! earth sun distance factor
562 : real(r8), intent(in) :: vmr(ncol,pver,max(1,gas_pcnst)) ! vmr
563 : real(r8), intent(in) :: col_dens(ncol,pver,ncol_abs) ! column densities (molecules/cm^2)
564 : real(r8), intent(in) :: zen_angle(ncol) ! solar zenith angle (radians)
565 : real(r8), intent(in) :: srf_alb(pcols) ! surface albedo
566 : real(r8), intent(in) :: lwc(ncol,pver) ! liquid water content (kg/kg)
567 : real(r8), intent(in) :: clouds(ncol,pver) ! cloud fraction
568 : real(r8), intent(in) :: pmid(pcols,pver) ! midpoint pressure (Pa)
569 : real(r8), intent(in) :: pdel(pcols,pver) ! pressure delta about midpoint (Pa)
570 : real(r8), intent(in) :: temper(pcols,pver) ! midpoint temperature (K)
571 : real(r8), intent(in) :: zmid(ncol,pver) ! midpoint height (km)
572 : real(r8), intent(in) :: zint(ncol,pver) ! interface height (km)
573 : real(r8), intent(in) :: invariants(ncol,pver,max(1,nfs)) ! invariant densities (molecules/cm^3)
574 : real(r8), intent(inout) :: photos(ncol,pver,phtcnt) ! photodissociation rates (1/s)
575 : type(physics_buffer_desc),pointer :: pbuf(:)
576 :
577 : !-----------------------------------------------------------------
578 : ! ... local variables
579 : !-----------------------------------------------------------------
580 : real(r8), parameter :: Pa2mb = 1.e-2_r8 ! pascals to mb
581 :
582 : integer :: i, k, m ! indicies
583 : integer :: astat
584 : real(r8) :: sza
585 : real(r8) :: alias_factor
586 : real(r8) :: fac1(pver) ! work space for j(no) calc
587 : real(r8) :: fac2(pver) ! work space for j(no) calc
588 : real(r8) :: colo3(pver) ! vertical o3 column density
589 : real(r8) :: parg(pver) ! vertical pressure array (hPa)
590 :
591 : real(r8) :: cld_line(pver) ! vertical cloud array
592 : real(r8) :: lwc_line(pver) ! vertical lwc array
593 : real(r8) :: eff_alb(pver) ! effective albedo from cloud modifications
594 : real(r8) :: cld_mult(pver) ! clould multiplier
595 0 : real(r8), allocatable :: lng_prates(:,:) ! photorates matrix (1/s)
596 0 : real(r8), allocatable :: sht_prates(:,:) ! photorates matrix (1/s)
597 0 : real(r8), allocatable :: euv_prates(:,:) ! photorates matrix (1/s)
598 :
599 :
600 0 : real(r8), allocatable :: zarg(:)
601 0 : real(r8), allocatable :: tline(:) ! vertical temperature array
602 0 : real(r8), allocatable :: o_den(:) ! o density (molecules/cm^3)
603 0 : real(r8), allocatable :: o2_den(:) ! o2 density (molecules/cm^3)
604 0 : real(r8), allocatable :: o3_den(:) ! o3 density (molecules/cm^3)
605 0 : real(r8), allocatable :: no_den(:) ! no density (molecules/cm^3)
606 0 : real(r8), allocatable :: n2_den(:) ! n2 density (molecules/cm^3)
607 0 : real(r8), allocatable :: jno_sht(:) ! no short photorate
608 0 : real(r8), allocatable :: jo2_sht(:,:) ! o2 short photorate
609 :
610 0 : real(r8), pointer :: ionRates(:,:,:) ! Pointer to ionization rates for O+,O2+,N+,N2+,NO+ in pbuf (s-1 from modules mo_jeuv and mo_jshort)
611 :
612 : integer :: n_jshrt_levs, p1, p2
613 : real(r8) :: ideltaZkm
614 :
615 0 : real(r8) :: qbktot(ncol,pver)
616 0 : real(r8) :: qbko1(ncol,pver)
617 0 : real(r8) :: qbko2(ncol,pver)
618 0 : real(r8) :: qbkn2(ncol,pver)
619 0 : real(r8) :: qbkn1(ncol,pver)
620 0 : real(r8) :: qbkno(ncol,pver)
621 :
622 0 : qbktot(:,:) = nan
623 0 : qbko1(:,:) = nan
624 0 : qbko2(:,:) = nan
625 0 : qbkn2(:,:) = nan
626 0 : qbkn1(:,:) = nan
627 0 : qbkno(:,:) = nan
628 :
629 : if( phtcnt < 1 ) then
630 : return
631 : end if
632 :
633 : if ((.not.do_jshort) .or. (ptop_ref < 10._r8)) then
634 : n_jshrt_levs = pver
635 : p1 = 1
636 : p2 = pver
637 : else
638 : n_jshrt_levs = pver+1
639 : p1 = 2
640 : p2 = pver+1
641 : endif
642 :
643 : allocate( zarg(n_jshrt_levs) )
644 : allocate( tline(n_jshrt_levs) )
645 : if (do_jshort) then
646 : allocate( o_den(n_jshrt_levs) )
647 : allocate( o2_den(n_jshrt_levs) )
648 : allocate( o3_den(n_jshrt_levs) )
649 : allocate( no_den(n_jshrt_levs) )
650 : allocate( n2_den(n_jshrt_levs) )
651 : allocate( jno_sht(n_jshrt_levs) )
652 : allocate( jo2_sht(n_jshrt_levs,2) )
653 : endif
654 :
655 : !-----------------------------------------------------------------
656 : ! ... allocate short, long temp arrays
657 : !-----------------------------------------------------------------
658 : if ( do_jeuv ) then
659 : if (neuv>0) then
660 : allocate( euv_prates(pver,neuv),stat=astat )
661 : if( astat /= 0 ) then
662 : write(iulog,*) 'photo: Failed to allocate euv_prates; error = ',astat
663 : call endrun
664 : end if
665 : endif
666 : endif
667 :
668 : if (nsht>0) then
669 : allocate( sht_prates(n_jshrt_levs,nsht),stat=astat )
670 : if( astat /= 0 ) then
671 : write(iulog,*) 'photo: Failed to allocate sht_prates; error = ',astat
672 : call endrun
673 : end if
674 : endif
675 :
676 : if (nlng>0) then
677 : allocate( lng_prates(nlng,pver),stat=astat )
678 : if( astat /= 0 ) then
679 : write(iulog,*) 'photo: Failed to allocate lng_prates; error = ',astat
680 : call endrun
681 : end if
682 : endif
683 :
684 : !-----------------------------------------------------------------
685 : ! ... zero all photorates
686 : !-----------------------------------------------------------------
687 : do m = 1,max(1,phtcnt)
688 : do k = 1,pver
689 : photos(:,k,m) = 0._r8
690 : end do
691 : end do
692 :
693 : !------------------------------------------------------------------------------------------------------------
694 : ! Point to production rates array in physics buffer where rates will be stored for ionosphere module
695 : ! access. Also, initialize rates to zero before column loop since only daylight values are filled
696 : !------------------------------------------------------------------------------------------------------------
697 : if (ion_rates_idx>0) then
698 : call pbuf_get_field(pbuf, ion_rates_idx, ionRates)
699 : ionRates(:,:,:) = 0._r8
700 : endif
701 :
702 : col_loop : do i = 1,ncol
703 : if (do_jshort) then
704 :
705 : if ( o_is_inv ) then
706 : o_den(p1:p2) = invariants(i,:pver,o_ndx)
707 : else
708 : o_den(p1:p2) = vmr(i,:pver,o_ndx) * invariants(i,:pver,indexm)
709 : endif
710 : if ( o2_is_inv ) then
711 : o2_den(p1:p2) = invariants(i,:pver,o2_ndx)
712 : else
713 : o2_den(p1:p2) = vmr(i,:pver,o2_ndx) * invariants(i,:pver,indexm)
714 : endif
715 : if ( o3_is_inv ) then
716 : o3_den(p1:p2) = invariants(i,:pver,o3_inv_ndx)
717 : else
718 : o3_den(p1:p2) = vmr(i,:,o3_ndx) * invariants(i,:pver,indexm)
719 : endif
720 : if ( n2_is_inv ) then
721 : n2_den(p1:p2) = invariants(i,:,n2_ndx)
722 : else
723 : n2_den(p1:p2) = vmr(i,:pver,n2_ndx) * invariants(i,:pver,indexm)
724 : endif
725 : if ( no_is_inv ) then
726 : no_den(p1:p2) = invariants(i,:pver,no_ndx)
727 : else
728 : no_den(p1:p2) = vmr(i,:pver,no_ndx) * invariants(i,:pver,indexm)
729 : endif
730 :
731 : endif
732 : sza = zen_angle(i)*r2d
733 : daylight : if( sza >= 0._r8 .and. sza < max_zen_angle ) then
734 : parg(:) = Pa2mb*pmid(i,:)
735 : colo3(:) = col_dens(i,:,1)
736 : fac1(:) = pdel(i,:)
737 : lwc_line(:) = lwc(i,:)
738 : cld_line(:) = clouds(i,:)
739 :
740 :
741 : tline(p1:p2) = temper(i,:pver)
742 :
743 : zarg(p1:p2) = zmid(i,:pver)
744 :
745 : if ( ptop_ref > 10._r8 ) then
746 : if (jppi_ndx > 0 ) photos(i,:,jppi_ndx) = photos(i,:,jppi_ndx) + esfact * 0.42_r8
747 : if (jepn1_ndx > 0 ) photos(i,:,jepn1_ndx) = photos(i,:,jepn1_ndx) + esfact * 1.4_r8
748 : if (jepn2_ndx > 0 ) photos(i,:,jepn2_ndx) = photos(i,:,jepn2_ndx) + esfact * 3.8e-1_r8
749 : if (jepn3_ndx > 0 ) photos(i,:,jepn3_ndx) = photos(i,:,jepn3_ndx) + esfact * 4.7e-2_r8
750 : if (jepn4_ndx > 0 ) photos(i,:,jepn4_ndx) = photos(i,:,jepn4_ndx) + esfact * 1.1_r8
751 : if (jepn6_ndx > 0 ) photos(i,:,jepn6_ndx) = photos(i,:,jepn6_ndx) + esfact * 8.0e-4_r8
752 : if (jepn7_ndx > 0 ) photos(i,:,jepn7_ndx) = photos(i,:,jepn7_ndx) + esfact * 5.2e-2_r8
753 : if (jpni1_ndx > 0 ) photos(i,:,jpni1_ndx) = photos(i,:,jpni1_ndx) + esfact * 0.47_r8
754 : if (jpni2_ndx > 0 ) photos(i,:,jpni2_ndx) = photos(i,:,jpni2_ndx) + esfact * 0.24_r8
755 : if (jpni3_ndx > 0 ) photos(i,:,jpni3_ndx) = photos(i,:,jpni3_ndx) + esfact * 0.15_r8
756 : if (jpni4_ndx > 0 ) photos(i,:,jpni4_ndx) = photos(i,:,jpni4_ndx) + esfact * 6.2e-3_r8
757 : ! added to v02
758 : if (jpni5_ndx > 0 ) photos(i,:,jpni5_ndx) = photos(i,:,jpni5_ndx) + esfact * 1.0_r8
759 : endif
760 : if (do_jshort) then
761 : if ( ptop_ref > 10._r8 ) then
762 : !-----------------------------------------------------------------
763 : ! Only for lower lid versions of CAM (i.e., not for WACCM)
764 : ! Column O3 and O2 above the top of the model
765 : ! DEK 20110224
766 : !-----------------------------------------------------------------
767 : ideltaZkm = 1._r8/(zint(i,1) - zint(i,2))
768 :
769 : !-----------------------------------------------------------------
770 : !... set density (units: molecules cm-3)
771 : !... used for jshort
772 : !....... Assuming a scale height of 7km for ozone
773 : !....... Assuming a scale height of 7km for O2
774 : !-----------------------------------------------------------------
775 : o3_den(1) = o3_den(2)*7.0_r8 * ideltaZkm
776 :
777 : o2_den(1) = o2_den(2)*7.0_r8 * ideltaZkm
778 :
779 : no_den(1) = no_den(2)*0.9_r8
780 :
781 : n2_den(1) = n2_den(2)*0.9_r8
782 :
783 : tline(1) = tline(2) + 5.0_r8
784 :
785 : zarg(1) = zarg(2) + (zint(i,1) - zint(i,2))
786 :
787 : endif
788 :
789 : !-----------------------------------------------------------------
790 : ! ... short wave length component
791 : !-----------------------------------------------------------------
792 : call jshort( n_jshrt_levs, sza, n2_den, o2_den, o3_den, &
793 : no_den, tline, zarg, jo2_sht, jno_sht, sht_prates )
794 :
795 : do m = 1,phtcnt
796 : if( sht_indexer(m) > 0 ) then
797 : alias_factor = pht_alias_mult(m,1)
798 : if( alias_factor == 1._r8 ) then
799 : photos(i,pver:1:-1,m) = sht_prates(1:pver,sht_indexer(m))
800 : else
801 : photos(i,pver:1:-1,m) = alias_factor * sht_prates(1:pver,sht_indexer(m))
802 : end if
803 : end if
804 : end do
805 :
806 : if( jno_ndx > 0 ) photos(i,pver:1:-1,jno_ndx) = jno_sht(1:pver)
807 : if( jo2_a_ndx > 0 ) photos(i,pver:1:-1,jo2_a_ndx) = jo2_sht(1:pver,2)
808 : if( jo2_b_ndx > 0 ) photos(i,pver:1:-1,jo2_b_ndx) = jo2_sht(1:pver,1)
809 : endif
810 :
811 : if ( do_jeuv ) then
812 : !-----------------------------------------------------------------
813 : ! ... euv photorates do not include cloud effects ??
814 : !-----------------------------------------------------------------
815 : call jeuv( pver, sza, o_den, o2_den, n2_den, zarg, euv_prates )
816 : do m = 1,neuv
817 : if( euv_indexer(m) > 0 ) then
818 : photos(i,:,euv_indexer(m)) = esfact * euv_prates(:,m)
819 : endif
820 : enddo
821 : endif
822 :
823 : !-----------------------------------------------------------------
824 : ! ... compute eff_alb and cld_mult -- needs to be before jlong
825 : !-----------------------------------------------------------------
826 : call cloud_mod( zen_angle(i), cld_line, lwc_line, fac1, srf_alb(i), &
827 : eff_alb, cld_mult )
828 : cld_mult(:) = esfact * cld_mult(:)
829 :
830 : !-----------------------------------------------------------------
831 : ! ... long wave length component
832 : !-----------------------------------------------------------------
833 : call jlong( pver, sza, eff_alb, parg, tline, colo3, lng_prates )
834 : do m = 1,phtcnt
835 : if( lng_indexer(m) > 0 ) then
836 : alias_factor = pht_alias_mult(m,2)
837 : if( alias_factor == 1._r8 ) then
838 : photos(i,:,m) = (photos(i,:,m) + lng_prates(lng_indexer(m),:))*cld_mult(:)
839 : else
840 : photos(i,:,m) = (photos(i,:,m) + alias_factor * lng_prates(lng_indexer(m),:))*cld_mult(:)
841 : end if
842 : end if
843 : end do
844 :
845 : !-----------------------------------------------------------------
846 : ! ... calculate j(no) from formula
847 : !-----------------------------------------------------------------
848 : if( (jno_ndx > 0) .and. (.not.do_jshort)) then
849 : if( has_o2_col .and. has_o3_col ) then
850 : fac1(:) = 1.e-8_r8 * (abs(col_dens(i,:,2)/cos(zen_angle(i))))**.38_r8
851 : fac2(:) = 5.e-19_r8 * abs(col_dens(i,:,1)/cos(zen_angle(i)))
852 : photos(i,:,jno_ndx) = photos(i,:,jno_ndx) + 4.5e-6_r8 * exp( -(fac1(:) + fac2(:)) )
853 : end if
854 : end if
855 :
856 : !-----------------------------------------------------------------
857 : ! ... add near IR correction to ho2no2
858 : !-----------------------------------------------------------------
859 : if( jho2no2_ndx > 0 ) then
860 : photos(i,:,jho2no2_ndx) = photos(i,:,jho2no2_ndx) + 1.e-5_r8*cld_mult(:)
861 : endif
862 :
863 : ! Save photo-ionization rates to physics buffer accessed in ionosphere module for WACCMX
864 : if (ion_rates_idx>0) then
865 :
866 : ionRates(i,1:pver,1:nIonRates) = esfact * euv_prates(1:pver,1:nIonRates)
867 :
868 : endif
869 :
870 : end if daylight
871 :
872 : if (do_jeuv) then
873 : !-----------------------------------------------------------------
874 : ! include background ionization ...
875 : ! outside daylight block so this is applied in all columns
876 : !-----------------------------------------------------------------
877 : call photo_bkgrnd_calc( f107, o_den, o2_den, n2_den, no_den, zint(i,:),&
878 : photos(i,:,:), qbko1_out=qbko1(i,:), qbko2_out=qbko2(i,:), &
879 : qbkn2_out=qbkn2(i,:), qbkn1_out=qbkn1(i,:), qbkno_out=qbkno(i,:) )
880 : endif
881 :
882 : end do col_loop
883 :
884 : if ( do_jeuv ) then
885 : qbktot(:ncol,:) = qbko1(:ncol,:)+qbko2(:ncol,:)+qbkn2(:ncol,:)+qbkn1(:ncol,:)+qbkno(:ncol,:)
886 : call outfld( 'Qbkgndtot', qbktot(:ncol,:),ncol, lchnk )
887 : call outfld( 'Qbkgnd_o1', qbko1(:ncol,:), ncol, lchnk )
888 : call outfld( 'Qbkgnd_o2', qbko2(:ncol,:), ncol, lchnk )
889 : call outfld( 'Qbkgnd_n2', qbkn2(:ncol,:), ncol, lchnk )
890 : call outfld( 'Qbkgnd_n1', qbkn1(:ncol,:), ncol, lchnk )
891 : call outfld( 'Qbkgnd_no', qbkno(:ncol,:), ncol, lchnk )
892 : endif
893 :
894 : if ( allocated(lng_prates) ) deallocate( lng_prates )
895 : if ( allocated(sht_prates) ) deallocate( sht_prates )
896 : if ( allocated(euv_prates) ) deallocate( euv_prates )
897 :
898 : if ( allocated(zarg) ) deallocate( zarg )
899 : if ( allocated(tline) ) deallocate( tline )
900 : if ( allocated(o_den) ) deallocate( o_den )
901 : if ( allocated(o2_den) ) deallocate( o2_den )
902 : if ( allocated(o3_den) ) deallocate( o3_den )
903 : if ( allocated(no_den) ) deallocate( no_den )
904 : if ( allocated(n2_den) ) deallocate( n2_den )
905 : if ( allocated(jno_sht) ) deallocate( jno_sht )
906 : if ( allocated(jo2_sht) ) deallocate( jo2_sht )
907 :
908 : call set_xnox_photo( photos, ncol )
909 :
910 0 : end subroutine table_photo
911 :
912 : subroutine cloud_mod( zen_angle, clouds, lwc, delp, srf_alb, &
913 : eff_alb, cld_mult )
914 : !-----------------------------------------------------------------------
915 : ! ... cloud alteration factors for photorates and albedo
916 : !-----------------------------------------------------------------------
917 :
918 : implicit none
919 :
920 : !-----------------------------------------------------------------------
921 : ! ... dummy arguments
922 : !-----------------------------------------------------------------------
923 : real(r8), intent(in) :: zen_angle ! zenith angle
924 : real(r8), intent(in) :: srf_alb ! surface albedo
925 : real(r8), intent(in) :: clouds(pver) ! cloud fraction
926 : real(r8), intent(in) :: lwc(pver) ! liquid water content (mass mr)
927 : real(r8), intent(in) :: delp(pver) ! del press about midpoint in pascals
928 : real(r8), intent(out) :: eff_alb(pver) ! effective albedo
929 : real(r8), intent(out) :: cld_mult(pver) ! photolysis mult factor
930 :
931 : !-----------------------------------------------------------------------
932 : ! ... local variables
933 : !-----------------------------------------------------------------------
934 : integer :: k
935 : real(r8) :: coschi
936 : real(r8) :: del_lwp(pver)
937 : real(r8) :: del_tau(pver)
938 : real(r8) :: above_tau(pver)
939 : real(r8) :: below_tau(pver)
940 : real(r8) :: above_cld(pver)
941 : real(r8) :: below_cld(pver)
942 : real(r8) :: above_tra(pver)
943 : real(r8) :: below_tra(pver)
944 : real(r8) :: fac1(pver)
945 : real(r8) :: fac2(pver)
946 :
947 : real(r8), parameter :: rgrav = 1._r8/9.80616_r8
948 :
949 : !---------------------------------------------------------
950 : ! ... modify lwc for cloud fraction and form
951 : ! liquid water path for each layer
952 : !---------------------------------------------------------
953 : where( clouds(:) /= 0._r8 )
954 : del_lwp(:) = rgrav * lwc(:) * delp(:) * 1.e3_r8 / clouds(:)
955 : elsewhere
956 : del_lwp(:) = 0._r8
957 : endwhere
958 : !---------------------------------------------------------
959 : ! ... form tau for each model layer
960 : !---------------------------------------------------------
961 : where( clouds(:) /= 0._r8 )
962 : del_tau(:) = del_lwp(:) *.155_r8 * clouds(:)**1.5_r8
963 : elsewhere
964 : del_tau(:) = 0._r8
965 : end where
966 : !---------------------------------------------------------
967 : ! ... form integrated tau from top down
968 : !---------------------------------------------------------
969 : above_tau(1) = 0._r8
970 : do k = 1,pverm
971 : above_tau(k+1) = del_tau(k) + above_tau(k)
972 : end do
973 : !---------------------------------------------------------
974 : ! ... form integrated tau from bottom up
975 : !---------------------------------------------------------
976 : below_tau(pver) = 0._r8
977 : do k = pverm,1,-1
978 : below_tau(k) = del_tau(k+1) + below_tau(k+1)
979 : end do
980 : !---------------------------------------------------------
981 : ! ... form vertically averaged cloud cover above and below
982 : !---------------------------------------------------------
983 : above_cld(1) = 0._r8
984 : do k = 1,pverm
985 : above_cld(k+1) = clouds(k) * del_tau(k) + above_cld(k)
986 : end do
987 : do k = 2,pver
988 : if( above_tau(k) /= 0._r8 ) then
989 : above_cld(k) = above_cld(k) / above_tau(k)
990 : else
991 : above_cld(k) = above_cld(k-1)
992 : end if
993 : end do
994 : below_cld(pver) = 0._r8
995 : do k = pverm,1,-1
996 : below_cld(k) = clouds(k+1) * del_tau(k+1) + below_cld(k+1)
997 : end do
998 : do k = pverm,1,-1
999 : if( below_tau(k) /= 0._r8 ) then
1000 : below_cld(k) = below_cld(k) / below_tau(k)
1001 : else
1002 : below_cld(k) = below_cld(k+1)
1003 : end if
1004 : end do
1005 : !---------------------------------------------------------
1006 : ! ... modify above_tau and below_tau via jfm
1007 : !---------------------------------------------------------
1008 : where( above_cld(2:pver) /= 0._r8 )
1009 : above_tau(2:pver) = above_tau(2:pver) / above_cld(2:pver)
1010 : end where
1011 : where( below_cld(:pverm) /= 0._r8 )
1012 : below_tau(:pverm) = below_tau(:pverm) / below_cld(:pverm)
1013 : end where
1014 : where( above_tau(2:pver) < 5._r8 )
1015 : above_cld(2:pver) = 0._r8
1016 : end where
1017 : where( below_tau(:pverm) < 5._r8 )
1018 : below_cld(:pverm) = 0._r8
1019 : end where
1020 : !---------------------------------------------------------
1021 : ! ... form transmission factors
1022 : !---------------------------------------------------------
1023 : above_tra(:) = 11.905_r8 / (9.524_r8 + above_tau(:))
1024 : below_tra(:) = 11.905_r8 / (9.524_r8 + below_tau(:))
1025 : !---------------------------------------------------------
1026 : ! ... form effective albedo
1027 : !---------------------------------------------------------
1028 : where( below_cld(:) /= 0._r8 )
1029 : eff_alb(:) = srf_alb + below_cld(:) * (1._r8 - below_tra(:)) &
1030 : * (1._r8 - srf_alb)
1031 : elsewhere
1032 : eff_alb(:) = srf_alb
1033 : end where
1034 : coschi = max( cos( zen_angle ),.5_r8 )
1035 : where( del_lwp(:)*.155_r8 < 5._r8 )
1036 : fac1(:) = 0._r8
1037 : elsewhere
1038 : fac1(:) = 1.4_r8 * coschi - 1._r8
1039 : end where
1040 : fac2(:) = min( 0._r8,1.6_r8*coschi*above_tra(:) - 1._r8 )
1041 : cld_mult(:) = 1._r8 + fac1(:) * clouds(:) + fac2(:) * above_cld(:)
1042 : cld_mult(:) = max( .05_r8,cld_mult(:) )
1043 :
1044 0 : end subroutine cloud_mod
1045 :
1046 0 : subroutine set_ub_col( col_delta, vmr, invariants, ptop, pdel, ncol, lchnk )
1047 : !---------------------------------------------------------------
1048 : ! ... set the column densities at the upper boundary
1049 : !---------------------------------------------------------------
1050 :
1051 : use chem_mods, only : nfs, ncol_abs=>nabscol, indexm
1052 : use chem_mods, only : nabscol, gas_pcnst, indexm
1053 : use chem_mods, only : gas_pcnst
1054 :
1055 : implicit none
1056 :
1057 : !---------------------------------------------------------------
1058 : ! ... dummy args
1059 : !---------------------------------------------------------------
1060 : real(r8), intent(in) :: ptop(pcols) ! top pressure (Pa)
1061 : integer, intent(in) :: ncol ! number of columns in current chunk
1062 : integer, intent(in) :: lchnk ! latitude indicies in chunk
1063 : real(r8), intent(in) :: vmr(ncol,pver,gas_pcnst) ! xported species vmr
1064 : real(r8), intent(in) :: pdel(pcols,pver) ! pressure delta about midpoints (Pa)
1065 : real(r8), intent(in) :: invariants(ncol,pver,nfs)
1066 : real(r8), intent(out) :: col_delta(ncol,0:pver,max(1,nabscol)) ! /cm**2 o2,o3 col dens above model
1067 :
1068 : !---------------------------------------------------------------
1069 : ! ... local variables
1070 : !---------------------------------------------------------------
1071 : !---------------------------------------------------------------
1072 : ! note: xfactor = 10.*r/(k*g) in cgs units.
1073 : ! the factor 10. is to convert pdel
1074 : ! from pascals to dyne/cm**2.
1075 : !---------------------------------------------------------------
1076 : real(r8), parameter :: xfactor = 2.8704e21_r8/(9.80616_r8*1.38044_r8)
1077 : integer :: k, kl, spc_ndx
1078 : real(r8) :: tint_vals(2)
1079 0 : real(r8) :: o2_exo_col(ncol)
1080 0 : real(r8) :: o3_exo_col(ncol)
1081 : integer :: i
1082 :
1083 : !---------------------------------------------------------------
1084 : ! ... assign column density at the upper boundary
1085 : ! the first column is o3 and the second is o2.
1086 : ! add 10 du o3 column above top of model.
1087 : !---------------------------------------------------------------
1088 : !---------------------------------------------------------------
1089 : ! ... set exo absorber columns
1090 : !---------------------------------------------------------------
1091 0 : has_abs_cols : if( has_o2_col .and. has_o3_col ) then
1092 0 : if( has_fixed_press ) then
1093 0 : kl = ki - 1
1094 0 : if( has_o2_col ) then
1095 0 : do i = 1,ncol
1096 0 : if ( kl > 0 ) then
1097 0 : tint_vals(1) = o2_exo_coldens(kl,i,lchnk,last) &
1098 0 : + delp * (o2_exo_coldens(ki,i,lchnk,last) &
1099 0 : - o2_exo_coldens(kl,i,lchnk,last))
1100 0 : tint_vals(2) = o2_exo_coldens(kl,i,lchnk,next) &
1101 : + delp * (o2_exo_coldens(ki,i,lchnk,next) &
1102 0 : - o2_exo_coldens(kl,i,lchnk,next))
1103 : else
1104 0 : tint_vals(1) = o2_exo_coldens( 1,i,lchnk,last)
1105 0 : tint_vals(2) = o2_exo_coldens( 1,i,lchnk,next)
1106 : endif
1107 0 : o2_exo_col(i) = tint_vals(1) + dels * (tint_vals(2) - tint_vals(1))
1108 : end do
1109 : else
1110 0 : o2_exo_col(:) = 0._r8
1111 : end if
1112 0 : if( has_o3_col ) then
1113 0 : do i = 1,ncol
1114 0 : if ( kl > 0 ) then
1115 0 : tint_vals(1) = o3_exo_coldens(kl,i,lchnk,last) &
1116 0 : + delp * (o3_exo_coldens(ki,i,lchnk,last) &
1117 0 : - o3_exo_coldens(kl,i,lchnk,last))
1118 0 : tint_vals(2) = o3_exo_coldens(kl,i,lchnk,next) &
1119 : + delp * (o3_exo_coldens(ki,i,lchnk,next) &
1120 0 : - o3_exo_coldens(kl,i,lchnk,next))
1121 : else
1122 0 : tint_vals(1) = o3_exo_coldens( 1,i,lchnk,last)
1123 0 : tint_vals(2) = o3_exo_coldens( 1,i,lchnk,next)
1124 : endif
1125 0 : o3_exo_col(i) = tint_vals(1) + dels * (tint_vals(2) - tint_vals(1))
1126 : end do
1127 : else
1128 0 : o3_exo_col(:) = 0._r8
1129 : end if
1130 : #ifdef DEBUG
1131 0 : write(iulog,*) '-----------------------------------'
1132 0 : write(iulog,*) 'o2_exo_col'
1133 0 : write(iulog,'(1p,5g15.7)') o2_exo_col(:)
1134 0 : write(iulog,*) 'o3_exo_col'
1135 0 : write(iulog,'(1p,5g15.7)') o3_exo_col(:)
1136 0 : write(iulog,*) '-----------------------------------'
1137 : #endif
1138 : else
1139 : !---------------------------------------------------------------
1140 : ! ... do pressure interpolation
1141 : !---------------------------------------------------------------
1142 0 : call p_interp( lchnk, ncol, ptop, o2_exo_col, o3_exo_col )
1143 : end if
1144 : else
1145 0 : o2_exo_col(:) = 0._r8
1146 0 : o3_exo_col(:) = 0._r8
1147 : end if has_abs_cols
1148 :
1149 0 : if( o3rad_ndx > 0 ) then
1150 : spc_ndx = o3rad_ndx
1151 : else
1152 0 : spc_ndx = ox_ndx
1153 : end if
1154 0 : if( spc_ndx < 1 ) then
1155 0 : spc_ndx = o3_ndx
1156 : end if
1157 0 : if( spc_ndx > 0 ) then
1158 0 : col_delta(:,0,1) = o3_exo_col(:)
1159 : do k = 1,pver
1160 0 : col_delta(:ncol,k,1) = xfactor * pdel(:ncol,k) * vmr(:ncol,k,spc_ndx)
1161 : end do
1162 0 : else if( o3_inv_ndx > 0 ) then
1163 0 : col_delta(:,0,1) = o3_exo_col(:)
1164 0 : do k = 1,pver
1165 0 : col_delta(:ncol,k,1) = xfactor * pdel(:ncol,k) * invariants(:ncol,k,o3_inv_ndx)/invariants(:ncol,k,indexm)
1166 : end do
1167 : else
1168 0 : col_delta(:,:,1) = 0._r8
1169 : end if
1170 : if( ncol_abs > 1 ) then
1171 0 : if( o2_ndx > 1 ) then
1172 0 : col_delta(:,0,2) = o2_exo_col(:)
1173 0 : if( o2_is_inv ) then
1174 0 : do k = 1,pver
1175 0 : col_delta(:ncol,k,2) = xfactor * pdel(:ncol,k) * invariants(:ncol,k,o2_ndx)/invariants(:ncol,k,indexm)
1176 : end do
1177 : else
1178 : do k = 1,pver
1179 0 : col_delta(:ncol,k,2) = xfactor * pdel(:ncol,k) * vmr(:ncol,k,o2_ndx)
1180 : end do
1181 : endif
1182 : else
1183 0 : col_delta(:,:,2) = 0._r8
1184 : end if
1185 : end if
1186 :
1187 0 : end subroutine set_ub_col
1188 :
1189 0 : subroutine p_interp( lchnk, ncol, ptop, o2_exo_col, o3_exo_col )
1190 : !---------------------------------------------------------------
1191 : ! ... pressure interpolation for exo col density
1192 : !---------------------------------------------------------------
1193 :
1194 : implicit none
1195 :
1196 : !---------------------------------------------------------------
1197 : ! ... dummy arguments
1198 : !---------------------------------------------------------------
1199 : integer, intent(in) :: ncol ! no. of columns
1200 : real(r8), intent(out) :: o2_exo_col(ncol) ! exo model o2 column density (molecules/cm^2)
1201 : real(r8), intent(out) :: o3_exo_col(ncol) ! exo model o3 column density (molecules/cm^2)
1202 : integer, intent(in) :: lchnk ! latitude indicies in chunk
1203 : real(r8) :: ptop(pcols) ! top pressure (Pa)
1204 :
1205 : !---------------------------------------------------------------
1206 : ! ... local variables
1207 : !---------------------------------------------------------------
1208 : integer :: i, ki, kl
1209 : real(r8) :: pinterp
1210 : real(r8) :: delp
1211 : real(r8) :: tint_vals(2)
1212 :
1213 0 : do i = 1,ncol
1214 0 : pinterp = ptop(i)
1215 0 : if( pinterp < levs(1) ) then
1216 : ki = 0
1217 : delp = 0._r8
1218 : else
1219 0 : do ki = 2,n_exo_levs
1220 0 : if( pinterp <= levs(ki) ) then
1221 0 : delp = log( pinterp/levs(ki-1) )/log( levs(ki)/levs(ki-1) )
1222 0 : exit
1223 : end if
1224 : end do
1225 : end if
1226 0 : kl = ki - 1
1227 0 : if( has_o2_col ) then
1228 0 : tint_vals(1) = o2_exo_coldens(kl,i,lchnk,last) &
1229 0 : + delp * (o2_exo_coldens(ki,i,lchnk,last) &
1230 0 : - o2_exo_coldens(kl,i,lchnk,last))
1231 0 : tint_vals(2) = o2_exo_coldens(kl,i,lchnk,next) &
1232 : + delp * (o2_exo_coldens(ki,i,lchnk,next) &
1233 0 : - o2_exo_coldens(kl,i,lchnk,next))
1234 0 : o2_exo_col(i) = tint_vals(1) + dels * (tint_vals(2) - tint_vals(1))
1235 : else
1236 0 : o2_exo_col(i) = 0._r8
1237 : end if
1238 0 : if( has_o3_col ) then
1239 0 : tint_vals(1) = o3_exo_coldens(kl,i,lchnk,last) &
1240 0 : + delp * (o3_exo_coldens(ki,i,lchnk,last) &
1241 0 : - o3_exo_coldens(kl,i,lchnk,last))
1242 0 : tint_vals(2) = o3_exo_coldens(kl,i,lchnk,next) &
1243 : + delp * (o3_exo_coldens(ki,i,lchnk,next) &
1244 0 : - o3_exo_coldens(kl,i,lchnk,next))
1245 0 : o3_exo_col(i) = tint_vals(1) + dels * (tint_vals(2) - tint_vals(1))
1246 : else
1247 0 : o3_exo_col(i) = 0._r8
1248 : end if
1249 : end do
1250 :
1251 0 : end subroutine p_interp
1252 :
1253 0 : subroutine setcol( col_delta, col_dens )
1254 : !---------------------------------------------------------------
1255 : ! ... set the column densities
1256 : !---------------------------------------------------------------
1257 :
1258 : use chem_mods, only : ncol_abs=>nabscol
1259 :
1260 : implicit none
1261 :
1262 : !---------------------------------------------------------------
1263 : ! ... dummy arguments
1264 : !---------------------------------------------------------------
1265 : real(r8), intent(in) :: col_delta(:,0:,:) ! layer column densities (molecules/cm^2)
1266 : real(r8), intent(out) :: col_dens(:,:,:) ! column densities ( /cm**2 )
1267 :
1268 : !---------------------------------------------------------------
1269 : ! the local variables
1270 : !---------------------------------------------------------------
1271 : integer :: k, km1, m ! long, alt indicies
1272 :
1273 : !---------------------------------------------------------------
1274 : ! note: xfactor = 10.*r/(k*g) in cgs units.
1275 : ! the factor 10. is to convert pdel
1276 : ! from pascals to dyne/cm**2.
1277 : !---------------------------------------------------------------
1278 : real(r8), parameter :: xfactor = 2.8704e21_r8/(9.80616_r8*1.38044_r8)
1279 :
1280 : !---------------------------------------------------------------
1281 : ! ... compute column densities down to the
1282 : ! current eta index in the calling routine.
1283 : ! the first column is o3 and the second is o2.
1284 : !---------------------------------------------------------------
1285 0 : do m = 1,ncol_abs
1286 0 : col_dens(:,1,m) = col_delta(:,0,m) + .5_r8 * col_delta(:,1,m)
1287 0 : do k = 2,pver
1288 0 : km1 = k - 1
1289 0 : col_dens(:,k,m) = col_dens(:,km1,m) + .5_r8 * (col_delta(:,km1,m) + col_delta(:,k,m))
1290 : end do
1291 : enddo
1292 :
1293 0 : end subroutine setcol
1294 :
1295 0 : subroutine photo_timestep_init( calday )
1296 : use euvac, only : euvac_set_etf
1297 : use mo_jshort, only : jshort_timestep_init
1298 : use mo_jlong, only : jlong_timestep_init
1299 : use solar_euv_data, only : solar_euv_data_active
1300 :
1301 : !-----------------------------------------------------------------------------
1302 : ! ... setup the time interpolation
1303 : !-----------------------------------------------------------------------------
1304 :
1305 : implicit none
1306 :
1307 : !-----------------------------------------------------------------------------
1308 : ! ... dummy arguments
1309 : !-----------------------------------------------------------------------------
1310 : real(r8), intent(in) :: calday ! day of year at end of present time step
1311 :
1312 : !-----------------------------------------------------------------------------
1313 : ! ... local variables
1314 : !-----------------------------------------------------------------------------
1315 : integer :: m
1316 :
1317 0 : if ( do_jeuv ) then
1318 0 : if (.not.solar_euv_data_active) then
1319 0 : call euvac_set_etf( f107, f107a )
1320 : end if
1321 : endif
1322 :
1323 0 : if( has_o2_col .or. has_o3_col ) then
1324 0 : if( calday < days(1) ) then
1325 0 : next = 1
1326 0 : last = 12
1327 0 : dels = (365._r8 + calday - days(12)) / (365._r8 + days(1) - days(12))
1328 0 : else if( calday >= days(12) ) then
1329 0 : next = 1
1330 0 : last = 12
1331 0 : dels = (calday - days(12)) / (365._r8 + days(1) - days(12))
1332 : else
1333 0 : do m = 11,1,-1
1334 0 : if( calday >= days(m) ) then
1335 : exit
1336 : end if
1337 : end do
1338 0 : last = m
1339 0 : next = m + 1
1340 0 : dels = (calday - days(m)) / (days(m+1) - days(m))
1341 : end if
1342 : #ifdef DEBUG
1343 0 : if (masterproc) then
1344 0 : write(iulog,*) '-----------------------------------'
1345 0 : write(iulog,*) 'photo_timestep_init: diagnostics'
1346 0 : write(iulog,*) 'calday, last, next, dels = ',calday,last,next,dels
1347 0 : write(iulog,*) '-----------------------------------'
1348 : endif
1349 : #endif
1350 : end if
1351 :
1352 : !-----------------------------------------------------------------------
1353 : ! Set jlong etf
1354 : !-----------------------------------------------------------------------
1355 0 : call jlong_timestep_init
1356 :
1357 0 : if ( do_jshort ) then
1358 : !-----------------------------------------------------------------------
1359 : ! Set jshort etf
1360 : !-----------------------------------------------------------------------
1361 0 : call jshort_timestep_init
1362 : endif
1363 :
1364 0 : end subroutine photo_timestep_init
1365 :
1366 : !--------------------------------------------------------------------------
1367 : !--------------------------------------------------------------------------
1368 : subroutine set_xnox_photo( photos, ncol )
1369 0 : use chem_mods, only : phtcnt
1370 : implicit none
1371 : integer, intent(in) :: ncol
1372 : real(r8), intent(inout) :: photos(ncol,pver,phtcnt) ! photodissociation rates (1/s)
1373 :
1374 : if( jno2a_ndx > 0 .and. jno2_ndx > 0 ) then
1375 : photos(:,:,jno2a_ndx) = photos(:,:,jno2_ndx)
1376 : end if
1377 : if( jn2o5a_ndx > 0 .and. jn2o5_ndx > 0 ) then
1378 : photos(:,:,jn2o5a_ndx) = photos(:,:,jn2o5_ndx)
1379 : end if
1380 : if( jn2o5b_ndx > 0 .and. jn2o5_ndx > 0 ) then
1381 : photos(:,:,jn2o5b_ndx) = photos(:,:,jn2o5_ndx)
1382 : end if
1383 : if( jhno3a_ndx > 0 .and. jhno3_ndx > 0 ) then
1384 : photos(:,:,jhno3a_ndx) = photos(:,:,jhno3_ndx)
1385 : end if
1386 :
1387 : if( jno3a_ndx > 0 .and. jno3_ndx > 0 ) then
1388 : photos(:,:,jno3a_ndx) = photos(:,:,jno3_ndx)
1389 : end if
1390 : if( jho2no2a_ndx > 0 .and. jho2no2_ndx > 0 ) then
1391 : photos(:,:,jho2no2a_ndx) = photos(:,:,jho2no2_ndx)
1392 : end if
1393 : if( jmpana_ndx > 0 .and. jmpan_ndx > 0 ) then
1394 : photos(:,:,jmpana_ndx) = photos(:,:,jmpan_ndx)
1395 : end if
1396 : if( jpana_ndx > 0 .and. jpan_ndx > 0 ) then
1397 : photos(:,:,jpana_ndx) = photos(:,:,jpan_ndx)
1398 : end if
1399 : if( jonitra_ndx > 0 .and. jonitr_ndx > 0 ) then
1400 : photos(:,:,jonitra_ndx) = photos(:,:,jonitr_ndx)
1401 : end if
1402 : if( jo1da_ndx > 0 .and. jo1d_ndx > 0 ) then
1403 : photos(:,:,jo1da_ndx) = photos(:,:,jo1d_ndx)
1404 : end if
1405 : if( jo3pa_ndx > 0 .and. jo3p_ndx > 0 ) then
1406 : photos(:,:,jo3pa_ndx) = photos(:,:,jo3p_ndx)
1407 : end if
1408 :
1409 : endsubroutine set_xnox_photo
1410 :
1411 : end module mo_photo
|
