Line data Source code
1 : #ifdef AIX
2 : #define USE_ESSL
3 : #endif
4 : #define USE_BDE
5 :
6 : module mo_jshort
7 :
8 : use shr_kind_mod, only : r8 => shr_kind_r8
9 : use physconst, only : pi
10 : use mo_constants, only : d2r
11 : use cam_abortutils, only : endrun
12 : use cam_logfile, only : iulog
13 : use spmd_utils, only : masterproc
14 : use ppgrid, only : pver
15 : use phys_control, only : waccmx_is
16 :
17 : implicit none
18 :
19 : interface jshort
20 : module procedure jshort_photo
21 : module procedure jshort_hrates
22 : end interface
23 :
24 : private
25 : public :: jshort_init
26 : public :: jshort_timestep_init
27 : public :: jshort
28 : public :: sphers
29 : public :: slant_col
30 : public :: nj
31 :
32 : save
33 :
34 : !------------------------------------------------------------------------------
35 : ! ... define altitude and wavelength parameters
36 : !------------------------------------------------------------------------------
37 : integer, parameter :: num_ms93tuv = 4 ! wavelength bins for MS, 93
38 : integer, parameter :: nw_ms93 = 4 ! wavelength bins for MS, 93
39 : integer, parameter :: nsrc_tot = 19 ! total bins for SRC
40 : integer, parameter :: nsrb_tot = 14 ! total bins <200nm for SRB
41 : integer, parameter :: nsrbtuv = 17 ! total SRB bins in TUV
42 : real(r8), parameter :: hc = 6.62608e-34_r8 * 2.9979e8_r8 / 1.e-9_r8
43 : real(r8), parameter :: wc_o2_a = 175.05_r8 ! (nm)
44 : real(r8), parameter :: wc_o2_b = 242.37_r8 ! (nm)
45 : real(r8), parameter :: wc_o3_a = 310.32_r8 ! (nm)
46 : real(r8), parameter :: wc_o3_b = 1179.87_r8 ! (nm)
47 : real(r8), parameter :: we_ms(nw_ms93+1) = (/ 181.6_r8, 183.1_r8, 184.6_r8, 190.2_r8, 192.5_r8 /)
48 :
49 : integer :: nw ! Number of wavelength bins <200nm
50 : integer :: nj ! Number of photorates
51 : real(r8) :: wtno50(6,2)
52 : real(r8) :: wtno90(6,2)
53 : real(r8) :: wtno100(6,2)
54 : real(r8) :: csno50(6,2)
55 : real(r8) :: csno90(6,2)
56 : real(r8) :: csno100(6,2)
57 : real(r8) :: ac(20,nsrbtuv)
58 : real(r8) :: bc(20,nsrbtuv) ! chebyshev polynomial coeffs
59 : real(r8) :: wave_num(nsrbtuv)
60 : real(r8), allocatable :: wc(:)
61 : real(r8), allocatable :: we(:)
62 : real(r8), allocatable :: wlintv(:)
63 : real(r8), allocatable :: bde_o2_a(:)
64 : real(r8), allocatable :: bde_o2_b(:)
65 : real(r8), allocatable :: bde_o3_a(:)
66 : real(r8), allocatable :: bde_o3_b(:)
67 : real(r8), allocatable :: etfphot(:)
68 : real(r8), allocatable :: etfphot_ms93(:)
69 : real(r8), allocatable :: xs_o2src(:)
70 : real(r8), allocatable :: xs_o3a(:)
71 : real(r8), allocatable :: xs_o3b(:)
72 : real(r8), allocatable :: xs_wl(:,:)
73 :
74 : real(r8), parameter :: lno2_llimit = 38._r8 ! ln(NO2) lower limit
75 : real(r8), parameter :: lno2_ulimit = 56._r8 ! ln(NO2) upper limit
76 :
77 : contains
78 :
79 0 : subroutine jshort_init( xs_coef_file, xs_short_file, sht_indexer )
80 : !------------------------------------------------------------------------------
81 : ! ... initialize photorates for 120nm <= lambda <= 200nm
82 : !------------------------------------------------------------------------------
83 :
84 : use mo_util, only : rebin
85 : use solar_irrad_data, only : data_nbins=>nbins, data_we => we, data_etf => sol_etf
86 :
87 : implicit none
88 :
89 : !------------------------------------------------------------------------------
90 : ! ... dummy arguments
91 : !------------------------------------------------------------------------------
92 : character(len=*), intent(in) :: xs_coef_file, xs_short_file
93 : integer, intent(inout) :: sht_indexer(:)
94 :
95 : !------------------------------------------------------------------------------
96 : ! ... set the wavelength grid, exoatmospheric flux,
97 : ! and cross sections (for <200nm) - contained in
98 : ! a NetCDF file
99 : !------------------------------------------------------------------------------
100 0 : call get_crs( xs_short_file, sht_indexer )
101 0 : if(masterproc) then
102 0 : write(iulog,*) ' '
103 0 : write(iulog,*) '============================================'
104 0 : write(iulog,*) 'jshort_init: finished get_crs'
105 0 : write(iulog,*) 'jshort_init: nj, nw = ',nj,nw
106 0 : write(iulog,*) 'jshort_init: wc'
107 0 : write(iulog,*) wc(:)
108 0 : write(iulog,*) '============================================'
109 0 : write(iulog,*) ' '
110 : end if
111 0 : we(:nw) = wc(:nw) - .5_r8*wlintv(:nw)
112 0 : we(nw+1) = wc(nw) + .5_r8*wlintv(nw)
113 0 : if(masterproc) then
114 0 : write(iulog,*) ' '
115 0 : write(iulog,*) '-------------------------------------------'
116 0 : write(iulog,*) 'jshort_init: diagnostics before rebin'
117 0 : write(iulog,*) 'jshort_init: data_nbins, nw = ',data_nbins,nw
118 0 : write(iulog,*) 'jshort_init: data_we range'
119 0 : write(iulog,'(1p,5g15.7)') minval(data_we(:)),maxval(data_we(:))
120 0 : write(iulog,*) 'jshort_init: we range'
121 0 : write(iulog,'(1p,5g15.7)') minval(we(:)),maxval(we(:))
122 : end if
123 0 : call rebin( data_nbins, nw, data_we, we, data_etf, etfphot )
124 0 : if(masterproc) then
125 0 : write(iulog,*) 'jshort_init: etfphot'
126 0 : write(iulog,'(1p,5g15.7)') etfphot(:)
127 0 : write(iulog,*) '-------------------------------------------'
128 0 : write(iulog,*) ' '
129 0 : write(iulog,*) 'jshort_init: diagnostics for ms93'
130 0 : call rebin( data_nbins, nw_ms93, data_we, we_ms, data_etf, etfphot_ms93 )
131 0 : write(iulog,'(1p,5g15.7)') etfphot_ms93(:)
132 0 : write(iulog,*) '-------------------------------------------'
133 0 : write(iulog,*) ' '
134 : end if
135 : !------------------------------------------------------------------------------
136 : ! ... loads Chebyshev polynomial Coeff
137 : !------------------------------------------------------------------------------
138 0 : call xs_init(xs_coef_file)
139 :
140 : !------------------------------------------------------------------------------
141 : ! ... initialize no photolysis
142 : !------------------------------------------------------------------------------
143 0 : call jno_init
144 :
145 0 : end subroutine jshort_init
146 :
147 0 : subroutine get_crs( xs_short_file, sht_indexer )
148 : !=============================================================================
149 : ! PURPOSE:
150 : ! Reads a NetCDF file that contains:
151 : ! Cross_sections*quantum yield data for species <200nm.
152 : ! Exoatmospheric flux on the wavelength grid of the crs
153 : !=============================================================================
154 : ! PARAMETERS:
155 : ! Input:
156 : ! xs_short_file.... NetCDF file that contains the crs*QY for wavenum species
157 : !
158 : ! Output:
159 : ! xs_species.. Cross Sections * QY data for each species
160 : ! etfphot..... Exoatmospheric flux in photons cm-2 s-1 nm-1
161 : ! etfphot_ms93.Minshwanner and Siskind JNO SRB etf (on MS93 grid)
162 : ! wc.......... wavelength center (nm)
163 : ! numwl ...... Number of wavelengths < 200nm in NetCDF input file
164 : ! wlintv ..... Wavelength inteval for grid, nm
165 : !=============================================================================
166 : ! EDIT HISTORY:
167 : ! Created by Doug Kinnison, 1/14/2002
168 : ! Modified by S. Walters, 4/2/2003
169 : !=============================================================================
170 :
171 0 : use chem_mods, only : phtcnt, pht_alias_lst, rxt_tag_lst
172 : use ioFileMod, only : getfil
173 : use error_messages, only : alloc_err
174 : use pio, only : file_desc_t, pio_get_var, pio_closefile, pio_noerr, &
175 : pio_seterrorhandling, pio_bcast_error, pio_internal_error, pio_inq_varid, &
176 : pio_inq_dimlen, pio_nowrite, pio_inq_dimid
177 : use cam_pio_utils, only : cam_pio_openfile
178 : implicit none
179 :
180 : !------------------------------------------------------------------------------
181 : ! ... dummy arguments
182 : !------------------------------------------------------------------------------
183 : integer, intent(inout) :: sht_indexer(:)
184 : character(len=*), intent(in) :: xs_short_file
185 :
186 : !------------------------------------------------------------------------------
187 : ! ... local variables
188 : !------------------------------------------------------------------------------
189 : integer :: wn
190 : type(file_desc_t) :: ncid
191 : integer :: ierr
192 : integer :: i, m, ndx
193 : integer :: varid, dimid
194 : integer :: wrk_ndx(phtcnt)
195 0 : real(r8), allocatable :: xs_species(:)
196 : character(len=256) :: locfn
197 :
198 : !------------------------------------------------------------------------------
199 : ! ... open NetCDF File
200 : !------------------------------------------------------------------------------
201 0 : call getfil(xs_short_file, locfn, 0)
202 0 : call cam_pio_openfile(ncid, trim(locfn), PIO_NOWRITE)
203 :
204 : !------------------------------------------------------------------------------
205 : ! ... get dimensions
206 : !------------------------------------------------------------------------------
207 0 : ierr = pio_inq_dimid( ncid, 'numwl', dimid )
208 0 : ierr = pio_inq_dimlen( ncid, dimid, nw )
209 :
210 : !------------------------------------------------------------------------------
211 : ! ... check for cross section in dataset
212 : !------------------------------------------------------------------------------
213 0 : call pio_seterrorhandling(ncid, pio_bcast_error)
214 : do m = 1,phtcnt
215 : if( pht_alias_lst(m,1) == ' ' ) then
216 : ierr = pio_inq_varid( ncid, rxt_tag_lst(m), varid )
217 : if( ierr == PIO_noerr ) then
218 : sht_indexer(m) = varid
219 : end if
220 : else if( pht_alias_lst(m,1) == 'userdefined' ) then
221 : sht_indexer(m) = -1
222 : else
223 : ierr = pio_inq_varid( ncid, pht_alias_lst(m,1), varid )
224 : if( ierr == PIO_noerr ) then
225 : sht_indexer(m) = varid
226 : else
227 : write(iulog,*) 'get_crs : ',rxt_tag_lst(m)(:len_trim(rxt_tag_lst(m))),' alias ', &
228 : pht_alias_lst(m,1)(:len_trim(pht_alias_lst(m,1))),' not in dataset'
229 : call endrun
230 : end if
231 : end if
232 : end do
233 0 : call pio_seterrorhandling(ncid, pio_internal_error)
234 :
235 0 : if (masterproc) then
236 0 : write(iulog,*) ' '
237 0 : write(iulog,*) '###############################################'
238 0 : write(iulog,*) 'get_crs: sht_indexer'
239 0 : write(iulog,'(10i6)') sht_indexer(:)
240 0 : write(iulog,*) '###############################################'
241 0 : write(iulog,*) ' '
242 : endif
243 :
244 0 : nj = 0
245 : do m = 1,phtcnt
246 : if( sht_indexer(m) > 0 ) then
247 : if( any( sht_indexer(:m-1) == sht_indexer(m) ) ) then
248 : cycle
249 : end if
250 : nj = nj + 1
251 : end if
252 : end do
253 :
254 : !------------------------------------------------------------------------------
255 : ! ... allocate arrays
256 : !------------------------------------------------------------------------------
257 0 : allocate( wc(nw),stat=ierr )
258 0 : if( ierr /= 0 ) then
259 0 : call alloc_err( ierr, 'get_crs', 'wc', nw )
260 : end if
261 0 : allocate( we(nw+1),stat=ierr )
262 0 : if( ierr /= 0 ) then
263 0 : call alloc_err( ierr, 'get_crs', 'we', nw+1 )
264 : end if
265 0 : allocate( wlintv(nw),stat=ierr )
266 0 : if( ierr /= 0 ) then
267 0 : call alloc_err( ierr, 'get_crs', 'wlintv', nw )
268 : end if
269 0 : allocate( etfphot(nw),stat=ierr )
270 0 : if( ierr /= 0 ) then
271 0 : call alloc_err( ierr, 'get_crs', 'etfphot', nw )
272 : end if
273 0 : allocate( bde_o2_a(nw),bde_o2_b(nw),bde_o3_a(nw),bde_o3_b(nw),stat=ierr )
274 0 : if( ierr /= 0 ) then
275 0 : call alloc_err( ierr, 'get_crs', 'bde_o2_a ... bde_o3_b', nw )
276 : end if
277 0 : allocate( etfphot_ms93(nw_ms93),stat=ierr )
278 0 : if( ierr /= 0 ) then
279 0 : call alloc_err( ierr, 'get_crs', 'etfphot_ms93', nw_ms93 )
280 : end if
281 0 : allocate( xs_o2src(nw),stat=ierr )
282 0 : if( ierr /= 0 ) then
283 0 : call alloc_err( ierr, 'get_crs', 'xs_o2src', nw )
284 : end if
285 0 : allocate( xs_o3a(nw),xs_o3b(nw),stat=ierr )
286 0 : if( ierr /= 0 ) then
287 0 : call alloc_err( ierr, 'get_crs', 'xs_o3a,xs_o3b', nw )
288 : end if
289 0 : allocate( xs_species(nw),xs_wl(nw,nj),stat=ierr )
290 0 : if( ierr /= 0 ) then
291 0 : call alloc_err( ierr, 'get_crs', 'xs_species,xs_wl', nw*nj )
292 : end if
293 :
294 : !------------------------------------------------------------------------------
295 : ! ... read arrays
296 : !------------------------------------------------------------------------------
297 0 : ierr = pio_inq_varid( ncid, 'wc', varid )
298 0 : ierr = pio_get_var( ncid, varid, wc )
299 0 : ierr = pio_inq_varid( ncid, 'wlintv', varid )
300 0 : ierr = pio_get_var( ncid, varid, wlintv )
301 0 : ierr = pio_inq_varid( ncid, 'xs_o2src', varid )
302 0 : ierr = pio_get_var( ncid, varid, xs_o2src )
303 : ndx = 0
304 : do m = 1,phtcnt
305 : if( sht_indexer(m) > 0 ) then
306 : if( any( sht_indexer(:m-1) == sht_indexer(m) ) ) then
307 : cycle
308 : end if
309 : ierr = pio_get_var( ncid, sht_indexer(m), xs_species )
310 : ndx = ndx + 1
311 : xs_wl(:,ndx) = xs_species(:)*wlintv(:)
312 : end if
313 : end do
314 0 : deallocate( xs_species )
315 0 : if( ndx /= nj ) then
316 0 : write(iulog,*) 'get_crs : ndx count /= cross section count'
317 0 : call endrun
318 : end if
319 : !------------------------------------------------------------------------------
320 : ! ... get jo3 cross sections
321 : !------------------------------------------------------------------------------
322 0 : ierr = pio_inq_varid( ncid, 'jo3_a', varid )
323 0 : ierr = pio_get_var( ncid, varid, xs_o3a )
324 0 : ierr = pio_inq_varid( ncid, 'jo3_b', varid )
325 0 : ierr = pio_get_var( ncid, varid, xs_o3b )
326 : !------------------------------------------------------------------------------
327 : ! ... setup final sht_indexer
328 : !------------------------------------------------------------------------------
329 0 : ndx = 0
330 0 : wrk_ndx(:) = sht_indexer(:)
331 : do m = 1,phtcnt
332 : if( wrk_ndx(m) > 0 ) then
333 : ndx = ndx + 1
334 : i = wrk_ndx(m)
335 : where( wrk_ndx(:) == i )
336 : sht_indexer(:) = ndx
337 : wrk_ndx(:) = -100000
338 : end where
339 : end if
340 : end do
341 :
342 0 : call pio_closefile( ncid )
343 :
344 : #ifdef USE_BDE
345 0 : if (masterproc) write(iulog,*) 'Jshort using bdes'
346 : #else
347 : if (masterproc) write(iulog,*) 'Jshort not using bdes'
348 : #endif
349 0 : do wn = 1,nw
350 : #ifdef USE_BDE
351 0 : bde_o2_a(wn) = max(0._r8, hc*(wc_o2_a - wc(wn))/(wc_o2_a*wc(wn)) )
352 0 : bde_o2_b(wn) = max(0._r8, hc*(wc_o2_b - wc(wn))/(wc_o2_b*wc(wn)) )
353 0 : bde_o3_a(wn) = max(0._r8, hc*(wc_o3_a - wc(wn))/(wc_o3_a*wc(wn)) )
354 0 : bde_o3_b(wn) = max(0._r8, hc*(wc_o3_b - wc(wn))/(wc_o3_b*wc(wn)) )
355 : #else
356 : bde_o2_a(wn) = hc/wc(wn)
357 : bde_o2_b(wn) = hc/wc(wn)
358 : bde_o3_a(wn) = hc/wc(wn)
359 : bde_o3_b(wn) = hc/wc(wn)
360 : #endif
361 : end do
362 :
363 0 : end subroutine get_crs
364 :
365 0 : subroutine xs_init(xs_coef_file)
366 : !-------------------------------------------------------------
367 : ! ... Loads XS_COEFFS containing the Chebyshev
368 : ! polynomial coeffs necessary to calculate O2 effective
369 : ! cross-sections
370 : !-------------------------------------------------------------
371 :
372 0 : use ioFileMod, only : getfil
373 : use units, only : getunit, freeunit
374 :
375 : !------------------------------------------------------------------------------
376 : ! ... Dummy arguments
377 : !------------------------------------------------------------------------------
378 : character(len=*), intent(in) :: xs_coef_file
379 :
380 : !-------------------------------------------------------------
381 : ! ... Local variables
382 : !-------------------------------------------------------------
383 : integer :: unit ! file unit number
384 : integer :: istat ! i/o status
385 : integer :: i, j
386 : character(len=256) :: locfn
387 :
388 : !----------------------------------------------------------------------
389 : ! ... Get first strato photo rate file
390 : !----------------------------------------------------------------------
391 0 : call getfil(xs_coef_file, locfn, 0)
392 : !----------------------------------------------------------------------
393 : ! ... open file
394 : !----------------------------------------------------------------------
395 0 : unit = getunit()
396 : open( unit = unit, &
397 : file = trim(locfn), &
398 : status = 'old', &
399 : form = 'formatted', &
400 0 : iostat = istat )
401 0 : if( istat /= 0 ) then
402 : !----------------------------------------------------------------------
403 : ! ... Open error exit
404 : !----------------------------------------------------------------------
405 0 : write(iulog,*) 'xs_init: error ',istat,' opening file ',trim(locfn)
406 0 : call endrun
407 : end if
408 : !----------------------------------------------------------------------
409 : ! ... read file
410 : !----------------------------------------------------------------------
411 0 : read(unit,901)
412 0 : do i = 1,20
413 0 : read(unit,903,iostat=istat) ac(i,:)
414 0 : if( istat /= 0 ) then
415 0 : write(iulog,*) 'xs_init: error ',istat,' reading ac'
416 0 : call endrun
417 : end if
418 : end do
419 :
420 0 : read(unit,901)
421 0 : do i = 1,20
422 0 : read(unit,903,iostat=istat) bc(i,:)
423 0 : if( istat /= 0 ) then
424 0 : write(iulog,*) 'xs_init: error ',istat,' reading bc'
425 0 : call endrun
426 : end if
427 : end do
428 0 : close( unit )
429 0 : call freeunit( unit )
430 :
431 0 : wave_num(17:1:-1) = (/ (48250._r8 + (500._r8*i),i=1,17) /)
432 :
433 : 901 format( / )
434 : 903 format( 17(e23.14,1x))
435 :
436 0 : end subroutine xs_init
437 :
438 0 : subroutine jno_init
439 : !-------------------------------------------------------------
440 : ! ... Initialization for no photolysis
441 : !-------------------------------------------------------------
442 :
443 : implicit none
444 :
445 : !-------------------------------------------------------------
446 : ! ... Local variables
447 : !-------------------------------------------------------------
448 : real(r8), dimension(24) :: a, b, c
449 :
450 : !------------------------------------------------------------------------------
451 : ! ... 6 sub-intervals for O2 5-0 at 265K,
452 : ! 2 sub-sub-intervals for NO 0-0 at 250K
453 : !------------------------------------------------------------------------------
454 : a(:) = (/ 0._r8, 0._r8, 0._r8, 0._r8, &
455 : 5.12e-02_r8, 5.68e-03_r8, 1.32e-18_r8, 4.41e-17_r8, &
456 : 1.36e-01_r8, 1.52e-02_r8, 6.35e-19_r8, 4.45e-17_r8, &
457 : 1.65e-01_r8, 1.83e-02_r8, 7.09e-19_r8, 4.50e-17_r8, &
458 : 1.41e-01_r8, 1.57e-02_r8, 2.18e-19_r8, 2.94e-17_r8, &
459 0 : 4.50e-02_r8, 5.00e-03_r8, 4.67e-19_r8, 4.35e-17_r8 /)
460 :
461 : !------------------------------------------------------------------------------
462 : ! ... sub-intervals for o2 9-0 band,
463 : ! 2 sub-sub-intervals for no 1-0 at 250 k
464 : !------------------------------------------------------------------------------
465 : b(:) = (/ 0._r8, 0._r8, 0._r8, 0._r8, &
466 : 0._r8, 0._r8, 0._r8, 0._r8, &
467 : 1.93e-03_r8, 2.14e-04_r8, 3.05e-21_r8, 3.20e-21_r8, &
468 : 9.73e-02_r8, 1.08e-02_r8, 5.76e-19_r8, 5.71e-17_r8, &
469 : 9.75e-02_r8, 1.08e-02_r8, 2.29e-18_r8, 9.09e-17_r8, &
470 0 : 3.48e-02_r8, 3.86e-03_r8, 2.21e-18_r8, 6.00e-17_r8 /)
471 :
472 : !------------------------------------------------------------------------------
473 : ! ... sub-intervals for o2 10-0 band,
474 : ! 2 sub-sub-intervals for no 1-0 at 250 k
475 : !------------------------------------------------------------------------------
476 : c(:) = (/ 4.50e-02_r8, 5.00e-03_r8, 1.80e-18_r8, 1.40e-16_r8, &
477 : 1.80e-01_r8, 2.00e-02_r8, 1.50e-18_r8, 1.52e-16_r8, &
478 : 2.25e-01_r8, 2.50e-02_r8, 5.01e-19_r8, 7.00e-17_r8, &
479 : 2.25e-01_r8, 2.50e-02_r8, 7.20e-20_r8, 2.83e-17_r8, &
480 : 1.80e-01_r8, 2.00e-02_r8, 6.72e-20_r8, 2.73e-17_r8, &
481 0 : 4.50e-02_r8, 5.00e-03_r8, 1.49e-21_r8, 6.57e-18_r8 /)
482 :
483 0 : wtno50 (1:6,1) = a(1:24:4)
484 0 : wtno50 (1:6,2) = a(2:24:4)
485 0 : csno50 (1:6,1) = a(3:24:4)
486 0 : csno50 (1:6,2) = a(4:24:4)
487 0 : wtno90 (1:6,1) = b(1:24:4)
488 0 : wtno90 (1:6,2) = b(2:24:4)
489 0 : csno90 (1:6,1) = b(3:24:4)
490 0 : csno90 (1:6,2) = b(4:24:4)
491 0 : wtno100(1:6,1) = c(1:24:4)
492 0 : wtno100(1:6,2) = c(2:24:4)
493 0 : csno100(1:6,1) = c(3:24:4)
494 0 : csno100(1:6,2) = c(4:24:4)
495 :
496 0 : end subroutine jno_init
497 :
498 0 : subroutine jshort_timestep_init
499 : !---------------------------------------------------------------
500 : ! ... set etfphot if required
501 : !---------------------------------------------------------------
502 :
503 : use mo_util, only : rebin
504 : use solar_irrad_data, only : data_nbins=>nbins, data_we => we, data_etf => sol_etf
505 :
506 : implicit none
507 :
508 0 : call rebin( data_nbins, nw, data_we, we, data_etf, etfphot )
509 0 : call rebin( data_nbins, nw_ms93, data_we, we_ms, data_etf, etfphot_ms93 )
510 :
511 0 : end subroutine jshort_timestep_init
512 :
513 0 : subroutine jshort_hrates( nlev, zen, o2_vmr, o3_vmr, o2cc, &
514 0 : o3cc, tlev, zkm, mw, qrs, cparg, &
515 0 : lchnk, long, co2cc, scco2, do_diag )
516 : !==============================================================================!
517 : ! Subroutine Jshort !
518 : !==============================================================================!
519 : ! Purpose: !
520 : ! To calculate thermal heating rates for lamba < 200 nm
521 : !==============================================================================!
522 : ! This routine uses JO2 parameterizations based on: !
523 : ! Lyman alpha... Chabrillat and Kockarts, GRL, 25, 2659, 1998 !
524 : ! SRC .......... Brasseur and Solomon, 1986 (from TUV) !
525 : ! SRB .......... Koppers and Murtagh, Ann. Geophys., 14, 68-79, 1996 !
526 : ! (supplied by Dan Marsh, NCAR ACD !
527 : ! and JNO: !
528 : ! SRB .......... Minschwanner and Siskind, JGR< 98, 20401, 1993. !
529 : !==============================================================================!
530 : ! Input: !
531 : ! o2cc....... O2 concentration, molecule cm-3 !
532 : ! o3cc....... O3 concentration, molecule cm-3 !
533 : ! zen........ zenith angle, units = degrees !
534 : ! tlev....... Temperature Profile (K) !
535 : ! zkm ....... Altitude, km !
536 : ! !
537 : ! Output: !
538 : ! qrs ... short wavelength thermal heating rates
539 : !==============================================================================!
540 : ! !
541 : ! Approach: !
542 : ! !
543 : ! 1) Call sphers (taken from TUV) !
544 : ! -> derives dsdh and nid used in slant column routines !
545 : ! -> zenith angle dependent !
546 : ! !
547 : ! 2) Call slant_col (taken from TUV) !
548 : ! -> derives the slant column for each species !
549 : ! !
550 : ! 3) Calls get_crs !
551 : ! -> read a NetCDF file !
552 : ! -> returns cross sections*quantum yields for all species that !
553 : ! have absorption below 200nm. !
554 : ! !
555 : ! 4) Derives transmission and photolysis rates for selective species !
556 : ! !
557 : !==============================================================================!
558 : ! EDIT HISTORY: !
559 : ! Created by Doug Kinnison, 3/14/2002 !
560 : !==============================================================================!
561 :
562 0 : use physconst, only : avogad
563 : use error_messages, only : alloc_err
564 :
565 : implicit none
566 :
567 : integer, parameter :: branch = 2 ! two photolysis pathways for JO2
568 : real(r8), parameter :: km2cm = 1.e5_r8
569 :
570 : !------------------------------------------------------------------------------
571 : ! ... dummy arguments
572 : !------------------------------------------------------------------------------
573 : integer, intent(in) :: nlev ! model vertical levels
574 : integer, intent(in) :: lchnk ! chunk index
575 : integer, intent(in) :: long ! chunk index
576 : real(r8), intent(in) :: zen ! Zenith angle (degrees)
577 : real(r8), intent(in) :: o2_vmr(nlev) ! o2 conc (mol/mol)
578 : real(r8), intent(in) :: o3_vmr(nlev) ! o3 conc (mol/mol)
579 : real(r8), intent(in) :: o2cc(nlev) ! o2 conc (mol/cm^3)
580 : real(r8), intent(in) :: co2cc(nlev) ! co2 conc (mol/cm^3)
581 : real(r8), intent(in) :: o3cc(nlev) ! o3 conc (mol/cm^3)
582 : real(r8), intent(in) :: tlev(nlev) ! Temperature profile
583 : real(r8), intent(in) :: zkm(nlev) ! Altitude, km
584 : real(r8), intent(in) :: mw(nlev) ! atms molecular weight
585 : real(r8), intent(in) :: cparg(nlev) ! column specific heat capacity
586 : real(r8), intent(inout) :: qrs(:,:) ! sw heating rates
587 : real(r8), intent(out) :: scco2(nlev) ! co2 column concentration (molec/cm^2)
588 : logical, intent(in) :: do_diag
589 :
590 : !------------------------------------------------------------------------------
591 : ! ... local variables
592 : !------------------------------------------------------------------------------
593 : integer :: k, k1 ! Altitude indicies
594 : integer :: wn ! Wavelength index
595 : integer :: astat
596 0 : integer :: nid(0:nlev) ! Number of layers crossed by the direct
597 : ! beam when travelling from the top of the
598 : ! atmosphere to layer i; NID(i), i = 0..NZ-1
599 : real(r8) :: hfactor
600 0 : real(r8) :: dsdh(0:nlev,nlev) ! Slant path of direct beam through each
601 : ! layer crossed when travelling from the top of
602 : ! the atmosphere to layer i; DSDH(i,j), i = 0.
603 : ! NZ-1, j = 1..NZ-1 (see sphers.f)
604 0 : real(r8), allocatable :: fnorm(:,:) ! Normalized ETF
605 0 : real(r8), allocatable :: trans_o2(:,:) ! Transmission o2 (total)
606 0 : real(r8), allocatable :: trans_o3(:,:) ! Transmission, ozone
607 0 : real(r8), allocatable :: wrk(:) ! wrk array
608 0 : real(r8) :: jo2_lya(nlev) ! Total photolytic rate constant for Ly alpha
609 : real(r8) :: jo2_srb(nlev) ! Total JO2 for SRB
610 : real(r8) :: jo2_src(nlev) ! Total JO2 for SRC
611 0 : real(r8) :: delz(nlev) ! layer thickness (cm)
612 0 : real(r8) :: o2scol(nlev) ! O2 Slant Column
613 0 : real(r8) :: o3scol(nlev) ! O3 Slant Column
614 0 : real(r8) :: rmla(nlev) ! Transmission, Lyman Alpha (other species)
615 0 : real(r8) :: ro2la(nlev) ! Transmission, Lyman Alpha (for JO2)
616 0 : real(r8) :: tlevmin(nlev)
617 0 : real(r8) :: abs_col(nlev)
618 0 : real(r8) :: tsrb(nlev,nsrbtuv) ! Transmission in the SRB
619 0 : real(r8) :: xs_o2srb(nlev,nsrbtuv) ! Cross section * QY for O2 in SRB
620 :
621 0 : allocate( fnorm(nlev,nw),stat=astat )
622 0 : if( astat /= 0 ) then
623 0 : call alloc_err( astat, 'jshort_hrates', 'fnorm', nw*nlev )
624 : end if
625 0 : allocate( trans_o2(nlev,nw),stat=astat )
626 0 : if( astat /= 0 ) then
627 0 : call alloc_err( astat, 'jshort_hrates', 'trans_o2', nw*nlev )
628 : end if
629 0 : allocate( trans_o3(nlev,nw),stat=astat )
630 0 : if( astat /= 0 ) then
631 0 : call alloc_err( astat, 'jshort_hrates', 'trans_o3', nw*nlev )
632 : end if
633 0 : allocate( wrk(nw),stat=astat )
634 0 : if( astat /= 0 ) then
635 0 : call alloc_err( astat, 'jshort_hrates', 'wrk', nw )
636 : end if
637 :
638 : !------------------------------------------------------------------------------
639 : ! ... derive Slant Path for Spherical Atmosphere
640 : !------------------------------------------------------------------------------
641 0 : call sphers( nlev, zkm, zen, dsdh, nid )
642 :
643 : !------------------------------------------------------------------------------
644 : ! ... derive O2, O3, and CO2 slant columns
645 : !------------------------------------------------------------------------------
646 0 : delz(1:nlev-1) = km2cm*(zkm(1:nlev-1) - zkm(2:nlev))
647 0 : call slant_col( nlev, delz, dsdh, nid, o2cc, o2scol )
648 0 : call slant_col( nlev, delz, dsdh, nid, o3cc, o3scol )
649 0 : call slant_col( nlev, delz, dsdh, nid, co2cc, scco2 )
650 :
651 : !------------------------------------------------------------------------------
652 : ! ... transmission due to ozone
653 : !------------------------------------------------------------------------------
654 0 : do wn = 1,nw
655 0 : abs_col(:) = min( (xs_o3a(wn) + xs_o3b(wn))*o3scol(:),100._r8 )
656 0 : trans_o3(:,wn) = exp( -abs_col(:) )
657 : end do
658 :
659 : !------------------------------------------------------------------------------
660 : ! ... derive the cross section and transmission for
661 : ! molecular oxygen Lya, SRC, and SRB's
662 : !------------------------------------------------------------------------------
663 : ! ... transmission due to molecular oxygen in the SRC
664 : !------------------------------------------------------------------------------
665 0 : do wn = 1,nsrc_tot
666 0 : abs_col(:) = min( xs_o2src(wn)*o2scol(:),100._r8 )
667 0 : trans_o2(:,wn) = exp( -abs_col(:) )
668 : end do
669 :
670 : !------------------------------------------------------------------------------
671 : ! ... transmission and cross sections due to O2 at lyman alpha
672 : !------------------------------------------------------------------------------
673 0 : call lymana( nlev, o2scol, rmla, ro2la )
674 :
675 : !------------------------------------------------------------------------------
676 : ! ... place lya reduction faction in transmission array
677 : ! this must follow the SRC placement (above)
678 : !------------------------------------------------------------------------------
679 0 : trans_o2(:,1) = rmla(:)
680 :
681 : !------------------------------------------------------------------------------
682 : ! ... Molecular Oxygen, SRB
683 : !------------------------------------------------------------------------------
684 : ! ... Koppers Grid (see Koppers and Murtagh, Ann. Geophys., 14, 68-79, 1996)
685 : ! # wl(i) wl(i+1)
686 : ! 1 174.4 177.0
687 : ! 2 177.0 178.6
688 : ! 3 178.6 180.2
689 : ! 4 180.2 181.8
690 : ! 5 181.8 183.5
691 : ! 6 183.5 185.2
692 : ! 7 185.2 186.9
693 : ! 8 186.9 188.7
694 : ! 9 188.7 190.5
695 : ! 10 190.5 192.3
696 : ! 11 192.3 194.2
697 : ! 12 194.2 196.1
698 : ! 13 196.1 198.0
699 : ! 14 198.0 200.0 <<last wl bin for <200nm
700 : ! ----------------------
701 : ! 15 200.0 202.0
702 : ! 16 202.0 204.1
703 : ! 17 204.1 205.8
704 : !------------------------------------------------------------------------------
705 :
706 : !------------------------------------------------------------------------------
707 : ! ... Kopper SRB parameterization is used here
708 : !------------------------------------------------------------------------------
709 0 : tlevmin(nlev:1:-1) = min( max( tlev(:),150._r8 ),350._r8 )
710 0 : call calc_o2srb( nlev, nid, o2scol, tlevmin, tsrb, xs_o2srb )
711 :
712 : !------------------------------------------------------------------------------
713 : ! ... Place Koppers SRB transmission (1-14) on user
714 : ! grid #20 (wc=175.7nm)
715 : !------------------------------------------------------------------------------
716 0 : do wn = 1,nsrb_tot
717 0 : trans_o2(:,wn+nsrc_tot) = tsrb(:,wn)
718 : end do
719 :
720 : !------------------------------------------------------------------------------
721 : ! ... Derive the normalize flux at each altitude,
722 : ! corrected for O2 and O3 absorption
723 : !------------------------------------------------------------------------------
724 0 : do wn = 1,nw ! nw = 33 (nsrb_tot+nsrc_tot)
725 0 : fnorm(:,wn) = etfphot(wn)*trans_o2(:,wn)*trans_o3(:,wn)
726 : end do
727 :
728 : !------------------------------------------------------------------------------
729 : ! ... Derive the O2 rate constant and apply branching ratio (QY)
730 : !------------------------------------------------------------------------------
731 : ! ... SRC and SRB QY
732 : ! Longward of 174.65 the product is O2 + hv => O(3P) + O(3P)
733 : ! Shortward of 174.65 the product is O2 + hv => O(3P) + O(1D)
734 : ! The QY is assumed to be unity in both wavelength ranges.
735 : !
736 : ! ... Lyman Alpha QY
737 : ! O2 + hv -> O(3P) + O(3P) at Lyman Alpha has a QY = 0.47
738 : ! O2 + hv -> O(3P) + O(1D) at Lyman Alpha has a QY = 0.53
739 : ! Lacoursiere et al., J. Chem. Phys. 110., 1949-1958, 1999.
740 : !------------------------------------------------------------------------------
741 : ! ... lyman alpha
742 : !------------------------------------------------------------------------------
743 0 : jo2_lya(:) = etfphot(1)*ro2la(:)*wlintv(1)
744 :
745 0 : wrk(1:nsrc_tot) = xs_o2src(1:nsrc_tot)*wlintv(1:nsrc_tot) &
746 0 : *bde_o2_a(1:nsrc_tot)
747 0 : wrk(1) = 0._r8
748 : !------------------------------------------------------------------------------
749 : ! ... o2 src heating
750 : !------------------------------------------------------------------------------
751 0 : if( do_diag ) then
752 0 : write(iulog,*) '-------------------------------------------------'
753 0 : write(iulog,*) 'jshort_hrates: fnorm,wrk at long,lchnk = ',long,lchnk
754 0 : write(iulog,'(1p,5g12.5)') fnorm(nlev,1:nsrc_tot)
755 0 : write(iulog,*) ' '
756 0 : write(iulog,'(1p,5g12.5)') wrk(1:nsrc_tot)
757 0 : write(iulog,*) '-------------------------------------------------'
758 : end if
759 : #ifdef USE_ESSL
760 : call dgemm( 'N', 'N', nlev, 1, nsrc_tot, &
761 : 1._r8, fnorm, nlev, wrk, nw, &
762 : 0._r8, qrs, nlev )
763 : #else
764 0 : qrs(:,1) = matmul( fnorm(:,1:nsrc_tot),wrk(1:nsrc_tot) )
765 : #endif
766 : !------------------------------------------------------------------------------
767 : ! ... o2 srb heating
768 : !------------------------------------------------------------------------------
769 0 : do k = 1,nlev
770 0 : wrk(1:nsrb_tot) = xs_o2srb(k,1:nsrb_tot)*wlintv(nsrc_tot+1:nsrc_tot+nsrb_tot) &
771 0 : *bde_o2_b(nsrc_tot+1:nsrc_tot+nsrb_tot)
772 0 : qrs(k,2) = dot_product( fnorm(k,nsrc_tot+1:nsrc_tot+nsrb_tot),wrk(1:nsrb_tot) )
773 : end do
774 :
775 0 : if( do_diag ) then
776 0 : write(iulog,*) '-------------------------------------------------'
777 0 : write(iulog,*) 'jshort_hrates: lya,bde_o2_a,qrs(nlev) at long,lchnk = ',long,lchnk
778 0 : write(iulog,'(1p,5g12.5)') jo2_lya(nlev),bde_o2_a(2),qrs(nlev,1)
779 0 : write(iulog,*) '-------------------------------------------------'
780 : end if
781 : !------------------------------------------------------------------------------
782 : ! ... total o2 heating
783 : !------------------------------------------------------------------------------
784 : !------------------------------------------------------------------------------
785 : ! ... Branch 1, O2 + hv => O(3P) + O(3P); wavelengths >175nm
786 : !------------------------------------------------------------------------------
787 0 : qrs(:,2) = qrs(:,2) + jo2_lya(:)*.47_r8*bde_o2_b(2)
788 : !------------------------------------------------------------------------------
789 : ! ... Branch 2, O2 + hv => O(3P) + O(1D); wavelengths <175nm
790 : !------------------------------------------------------------------------------
791 0 : qrs(:,1) = qrs(:,1) + jo2_lya(:)*.53_r8*bde_o2_a(2)
792 0 : if( do_diag ) then
793 0 : write(iulog,*) '-------------------------------------------------'
794 0 : write(iulog,*) 'jshort_hrates: o2(1),qrs(nlev) at long,lchnk = ',long,lchnk
795 0 : write(iulog,'(1p,5g12.5)') o2_vmr(1),qrs(nlev,1)
796 0 : write(iulog,*) '-------------------------------------------------'
797 : end if
798 :
799 : !------------------------------------------------------------------------------
800 : ! ... the o3 heating rates
801 : !------------------------------------------------------------------------------
802 0 : wrk(:) = xs_o3a(:)*wlintv(:)*bde_o3_a(:)
803 : #ifdef USE_ESSL
804 : call dgemm( 'N', 'N', nlev, 1, nw, &
805 : 1._r8, fnorm, nlev, wrk, nw, &
806 : 0._r8, qrs(1,3), nlev )
807 : #else
808 0 : qrs(:,3) = matmul( fnorm,wrk )
809 : #endif
810 0 : wrk(:) = xs_o3b(:)*wlintv(:)*bde_o3_b(:)
811 : #ifdef USE_ESSL
812 : call dgemm( 'N', 'N', nlev, 1, nw, &
813 : 1._r8, fnorm, nlev, wrk, nw, &
814 : 0._r8, qrs(1,4), nlev )
815 : #else
816 0 : qrs(:,4) = matmul( fnorm,wrk )
817 : #endif
818 :
819 : !------------------------------------------------------------------------------
820 : ! ... form actual heating rates (k/s)
821 : !------------------------------------------------------------------------------
822 0 : do k = 1,nlev
823 0 : k1 = nlev - k + 1
824 0 : hfactor = avogad/(cparg(k1)*mw(k1))
825 0 : qrs(k,1) = qrs(k,1)*hfactor*o2_vmr(k1)
826 0 : qrs(k,2) = qrs(k,2)*hfactor*o2_vmr(k1)
827 0 : qrs(k,3) = qrs(k,3)*hfactor*o3_vmr(k1)
828 0 : qrs(k,4) = qrs(k,4)*hfactor*o3_vmr(k1)
829 : end do
830 :
831 0 : deallocate( fnorm, trans_o2, trans_o3, wrk )
832 :
833 0 : end subroutine jshort_hrates
834 :
835 0 : subroutine jshort_photo( nlev, zen, n2cc, o2cc, o3cc, &
836 0 : nocc, tlev, zkm, jo2_sht, jno_sht, jsht )
837 : !==============================================================================!
838 : ! Subroutine Jshort !
839 : ! !
840 : !==============================================================================!
841 : ! Purpose: !
842 : ! To calculate the total J for JO2, JNO, and selective species below 200nm.!
843 : ! !
844 : !==============================================================================!
845 : ! This routine uses JO2 parameterizations based on: !
846 : ! Lyman alpha... Chabrillat and Kockarts, GRL, 25, 2659, 1998 !
847 : ! SRC .......... Brasseur and Solomon, 1986 (from TUV) !
848 : ! SRB .......... Koppers and Murtagh, Ann. Geophys., 14, 68-79, 1996 !
849 : ! (supplied by Dan Marsh, NCAR ACD !
850 : ! and JNO: !
851 : ! SRB .......... Minschwanner and Siskind, JGR< 98, 20401, 1993. !
852 : ! !
853 : !==============================================================================!
854 : ! Input: !
855 : ! n2cc....... N2 concentration, molecule cm-3 !
856 : ! o2cc....... O2 concentration, molecule cm-3 !
857 : ! o3cc....... O3 concentration, molecule cm-3 !
858 : ! nocc....... NO concentration, molecule cm-3 !
859 : ! n2cc....... N2 concentration, molecule cm-3 !
860 : ! zen........ zenith angle, units = degrees !
861 : ! tlev....... Temperature Profile (K) !
862 : ! zkm ....... Altitude, km !
863 : ! !
864 : ! Output: !
865 : ! jo2_sht ... O2 photolytic rate constant, sec-1, <200nm !
866 : ! jno_sht ... NO photolytic rate constant, sec-1, SRB !
867 : ! jsht. Photolytic rate constant for other species below 200nm !
868 : ! !
869 : !==============================================================================!
870 : ! !
871 : ! Approach: !
872 : ! !
873 : ! 1) Call sphers (taken from TUV) !
874 : ! -> derives dsdh and nid used in slant column routines !
875 : ! -> zenith angle dependent !
876 : ! !
877 : ! 2) Call slant_col (taken from TUV) !
878 : ! -> derives the slant column for each species !
879 : ! !
880 : ! 3) Calls get_crs !
881 : ! -> read a NetCDF file !
882 : ! -> returns cross sections*quantum yields for all species that !
883 : ! have absorption below 200nm. !
884 : ! !
885 : ! 4) Derives transmission and photolysis rates for selective species !
886 : ! !
887 : !==============================================================================!
888 : ! EDIT HISTORY: !
889 : ! Created by Doug Kinnison, 3/14/2002 !
890 : !==============================================================================!
891 :
892 : use error_messages, only : alloc_err
893 :
894 : implicit none
895 :
896 : integer, parameter :: branch = 2 ! two photolysis pathways for JO2
897 : real(r8), parameter :: km2cm = 1.e5_r8
898 :
899 : !------------------------------------------------------------------------------
900 : ! ... dummy arguments
901 : !------------------------------------------------------------------------------
902 : integer, intent(in) :: nlev ! model vertical levels
903 : real(r8), intent(in) :: zen ! Zenith angle (degrees)
904 : real(r8), intent(in) :: n2cc(nlev) ! Molecular Nitrogen conc (mol/cm^3)
905 : real(r8), intent(in) :: o2cc(nlev) ! Molecular Oxygen conc (mol/cm^3)
906 : real(r8), intent(in) :: o3cc(nlev) ! Ozone concentration (mol/cm^3)
907 : real(r8), intent(in) :: nocc(nlev) ! Nitric Oxide conc (mol/cm^3)
908 : real(r8), intent(in) :: tlev(nlev) ! Temperature profile
909 : real(r8), intent(in) :: zkm(nlev) ! Altitude, km
910 : real(r8), intent(out) :: jo2_sht(nlev,branch) ! JO2, sec-1, <200nm
911 : real(r8), intent(out) :: jno_sht(nlev) ! JNO, sec-1, SRB
912 : real(r8), intent(out) :: jsht(:,:) ! Additional J's
913 :
914 : !------------------------------------------------------------------------------
915 : ! ... local variables
916 : !------------------------------------------------------------------------------
917 : integer :: k ! Altitude index
918 : integer :: wn ! Wavelength index
919 : integer :: astat
920 0 : integer :: nid(0:nlev) ! Number of layers crossed by the direct
921 : ! beam when travelling from the top of the
922 : ! atmosphere to layer i; NID(i), i = 0..NZ-1
923 : real(r8) :: hfactor
924 0 : real(r8) :: dsdh(0:nlev,nlev) ! Slant path of direct beam through each
925 : ! layer crossed when travelling from the top of
926 : ! the atmosphere to layer i; DSDH(i,j), i = 0.
927 : ! NZ-1, j = 1..NZ-1 (see sphers.f)
928 0 : real(r8), allocatable :: fnorm(:,:) ! Normalized ETF
929 0 : real(r8), allocatable :: trans_o2(:,:) ! Transmission o2 (total)
930 0 : real(r8), allocatable :: trans_o3(:,:) ! Transmission, ozone
931 0 : real(r8), allocatable :: wrk(:) ! wrk array
932 0 : real(r8) :: jo2_lya(nlev) ! Total photolytic rate constant for Ly alpha
933 0 : real(r8) :: jo2_srb(nlev) ! Total JO2 for SRB
934 0 : real(r8) :: jo2_src(nlev) ! Total JO2 for SRC
935 0 : real(r8) :: delz(nlev) ! layer thickness (cm)
936 0 : real(r8) :: o2scol(nlev) ! O2 Slant Column
937 0 : real(r8) :: o3scol(nlev) ! O3 Slant Column
938 0 : real(r8) :: noscol(nlev) ! NO Slant Column
939 0 : real(r8) :: rmla(nlev) ! Transmission, Lyman Alpha (other species)
940 0 : real(r8) :: ro2la(nlev) ! Transmission, Lyman Alpha (for JO2)
941 0 : real(r8) :: tlevmin(nlev)
942 0 : real(r8) :: abs_col(nlev)
943 0 : real(r8) :: tsrb(nlev,nsrbtuv) ! Transmission in the SRB
944 0 : real(r8) :: xs_o2srb(nlev,nsrbtuv) ! Cross section * QY for O2 in SRB
945 :
946 0 : allocate( fnorm(nlev,nw),stat=astat )
947 0 : if( astat /= 0 ) then
948 0 : call alloc_err( astat, 'jshort_photo', 'fnorm', nw*nlev )
949 : end if
950 0 : allocate( trans_o2(nlev,nw),stat=astat )
951 0 : if( astat /= 0 ) then
952 0 : call alloc_err( astat, 'jshort_photo', 'trans_o2', nw*nlev )
953 : end if
954 0 : allocate( trans_o3(nlev,nw),stat=astat )
955 0 : if( astat /= 0 ) then
956 0 : call alloc_err( astat, 'jshort_photo', 'trans_o3', nw*nlev )
957 : end if
958 0 : allocate( wrk(nw),stat=astat )
959 0 : if( astat /= 0 ) then
960 0 : call alloc_err( astat, 'jshort_photo', 'wrk', nw )
961 : end if
962 :
963 : !------------------------------------------------------------------------------
964 : ! ... Derive Slant Path for Spherical Atmosphere
965 : !------------------------------------------------------------------------------
966 0 : call sphers( nlev, zkm, zen, dsdh, nid )
967 :
968 : !------------------------------------------------------------------------------
969 : ! ... Derive O2, O3, and NO Slant Column
970 : !------------------------------------------------------------------------------
971 0 : delz(1:nlev-1) = km2cm*(zkm(1:nlev-1) - zkm(2:nlev))
972 0 : call slant_col( nlev, delz, dsdh, nid, o2cc, o2scol )
973 0 : call slant_col( nlev, delz, dsdh, nid, o3cc, o3scol )
974 0 : call slant_col( nlev, delz, dsdh, nid, nocc, noscol )
975 :
976 : !------------------------------------------------------------------------------
977 : ! ... Transmission due to ozone
978 : !------------------------------------------------------------------------------
979 0 : do wn = 1,nw
980 0 : abs_col(:) = min( (xs_o3a(wn) + xs_o3b(wn))*o3scol(:),100._r8 )
981 0 : trans_o3(:,wn) = exp( -abs_col(:) )
982 : end do
983 :
984 : !------------------------------------------------------------------------------
985 : ! ... Derive the cross section and transmission for
986 : ! molecular oxygen Lya, SRC, and SRB's
987 : !------------------------------------------------------------------------------
988 : !------------------------------------------------------------------------------
989 : ! ... Transmission due to molecular oxygen in the SRC
990 : !------------------------------------------------------------------------------
991 0 : do wn = 1,nsrc_tot
992 0 : abs_col(:) = min( xs_o2src(wn)*o2scol(:),100._r8 )
993 0 : trans_o2(:,wn) = exp( -abs_col(:) )
994 : end do
995 :
996 : !------------------------------------------------------------------------------
997 : ! ... Transmission and cross sections due to O2 at lyman alpha
998 : !------------------------------------------------------------------------------
999 0 : call lymana( nlev, o2scol, rmla, ro2la )
1000 :
1001 : !------------------------------------------------------------------------------
1002 : ! ... Place lya reduction faction in transmission array
1003 : ! ... This must follow the SRC placement (above)
1004 : !------------------------------------------------------------------------------
1005 0 : trans_o2(:,1) = rmla(:)
1006 :
1007 : !------------------------------------------------------------------------------
1008 : ! ... Molecular Oxygen, SRB
1009 : !------------------------------------------------------------------------------
1010 : ! ... Koppers Grid (see Koppers and Murtagh, Ann. Geophys., 14, 68-79, 1996)
1011 : ! # wl(i) wl(i+1)
1012 : ! 1 174.4 177.0
1013 : ! 2 177.0 178.6
1014 : ! 3 178.6 180.2
1015 : ! 4 180.2 181.8
1016 : ! 5 181.8 183.5
1017 : ! 6 183.5 185.2
1018 : ! 7 185.2 186.9
1019 : ! 8 186.9 188.7
1020 : ! 9 188.7 190.5
1021 : ! 10 190.5 192.3
1022 : ! 11 192.3 194.2
1023 : ! 12 194.2 196.1
1024 : ! 13 196.1 198.0
1025 : ! 14 198.0 200.0 <<last wl bin for <200nm
1026 : ! ----------------------
1027 : ! 15 200.0 202.0
1028 : ! 16 202.0 204.1
1029 : ! 17 204.1 205.8
1030 : !------------------------------------------------------------------------------
1031 :
1032 : !------------------------------------------------------------------------------
1033 : ! ... Kopper SRB parameterization is used here
1034 : !------------------------------------------------------------------------------
1035 0 : tlevmin(nlev:1:-1) = min( max( tlev(:),150._r8 ),350._r8 )
1036 0 : call calc_o2srb( nlev, nid, o2scol, tlevmin, tsrb, xs_o2srb )
1037 :
1038 : !------------------------------------------------------------------------------
1039 : ! ... Place Koppers SRB transmission (1-14) on user
1040 : ! grid #20 (wc=175.7nm)
1041 : !------------------------------------------------------------------------------
1042 0 : do wn = 1,nsrb_tot
1043 0 : trans_o2(:,wn+nsrc_tot) = tsrb(:,wn)
1044 : end do
1045 :
1046 : !------------------------------------------------------------------------------
1047 : ! ... Derive the normalize flux at each altitude,
1048 : ! corrected for O2 and O3 absorption
1049 : !------------------------------------------------------------------------------
1050 0 : do wn = 1,nw ! nw = 33 (nsrb_tot+nsrc_tot)
1051 0 : fnorm(:,wn) = etfphot(wn)*trans_o2(:,wn)*trans_o3(:,wn)
1052 : end do
1053 :
1054 : !------------------------------------------------------------------------------
1055 : ! ... Derive the O2 rate constant and apply branching ratio (QY)
1056 : !------------------------------------------------------------------------------
1057 : ! ... SRC and SRB QY
1058 : ! Longward of 174.65 the product is O2 + hv => O(3P) + O(3P)
1059 : ! Shortward of 174.65 the product is O2 + hv => O(3P) + O(1D)
1060 : ! The QY is assumed to be unity in both wavelength ranges.
1061 : !
1062 : ! ... Lyman Alpha QY
1063 : ! O2 + hv -> O(3P) + O(3P) at Lyman Alpha has a QY = 0.47
1064 : ! O2 + hv -> O(3P) + O(1D) at Lyman Alpha has a QY = 0.53
1065 : ! Lacoursiere et al., J. Chem. Phys. 110., 1949-1958, 1999.
1066 : !------------------------------------------------------------------------------
1067 : ! ... Lyman Alpha
1068 : !------------------------------------------------------------------------------
1069 0 : jo2_lya(:) = etfphot(1)*ro2la(:)*wlintv(1)
1070 :
1071 0 : wrk(1:nsrc_tot) = xs_o2src(1:nsrc_tot)*wlintv(1:nsrc_tot)
1072 0 : wrk(1) = 0._r8
1073 : !------------------------------------------------------------------------------
1074 : ! ... o2 src photolysis
1075 : !------------------------------------------------------------------------------
1076 : #ifdef USE_ESSL
1077 : call dgemm( 'N', 'N', nlev, 1, nsrc_tot, &
1078 : 1._r8, fnorm, nlev, wrk, nw, &
1079 : 0._r8, jo2_src, nlev )
1080 : #else
1081 0 : jo2_src(:) = matmul( fnorm(:,1:nsrc_tot),wrk(1:nsrc_tot) )
1082 : #endif
1083 : !------------------------------------------------------------------------------
1084 : ! ... o2 srb photolysis
1085 : !------------------------------------------------------------------------------
1086 0 : do k = 1,nlev
1087 0 : wrk(1:nsrb_tot) = xs_o2srb(k,1:nsrb_tot)*wlintv(nsrc_tot+1:nsrc_tot+nsrb_tot)
1088 0 : jo2_srb(k) = dot_product( fnorm(k,nsrc_tot+1:nsrc_tot+nsrb_tot),wrk(1:nsrb_tot) )
1089 : end do
1090 :
1091 : !------------------------------------------------------------------------------
1092 : ! ... total o2 photolysis
1093 : !------------------------------------------------------------------------------
1094 : !------------------------------------------------------------------------------
1095 : ! ... Branch 1, O2 + hv => O(3P) + O(3P); wavelengths >175nm
1096 : !------------------------------------------------------------------------------
1097 0 : jo2_sht(:,1) = jo2_lya(:)*.47_r8 + jo2_srb(:)
1098 : !------------------------------------------------------------------------------
1099 : ! ... Branch 2, O2 + hv => O(3P) + O(1D); wavelengths <175nm
1100 : !------------------------------------------------------------------------------
1101 0 : jo2_sht(:,2) = jo2_lya(:)*.53_r8 + jo2_src(:)
1102 :
1103 : !------------------------------------------------------------------------------
1104 : ! ... Derive the NO rate constant Minsch. and Siskind, JGR, 98, 20401, 1993
1105 : !------------------------------------------------------------------------------
1106 : call calc_jno( nlev, etfphot_ms93, n2cc, o2scol, o3scol, &
1107 0 : noscol, jno_sht )
1108 :
1109 : !------------------------------------------------------------------------------
1110 : ! ... Derive addtional rate constants for species with wl < 200 nm.!
1111 : ! Temperature dependence of the cross sections are not included in this
1112 : ! version.
1113 : !------------------------------------------------------------------------------
1114 : #if defined USE_ESSL
1115 : call dgemm( 'N', 'N', nlev, nj, nw, &
1116 : 1._r8, fnorm, nlev, xs_wl, nw, &
1117 : 0._r8, jsht, nlev )
1118 : #else
1119 0 : jsht(:,:) = matmul( fnorm,xs_wl )
1120 : #endif
1121 :
1122 0 : deallocate( fnorm, trans_o2, trans_o3, wrk )
1123 :
1124 0 : end subroutine jshort_photo
1125 :
1126 0 : subroutine sphers( nlev, z, zenith_angle, dsdh, nid )
1127 : !=============================================================================!
1128 : ! Subroutine sphers !
1129 : !=============================================================================!
1130 : ! PURPOSE: !
1131 : ! Calculate slant path over vertical depth ds/dh in spherical geometry. !
1132 : ! Calculation is based on: A.Dahlback, and K.Stamnes, A new spheric model !
1133 : ! for computing the radiation field available for photolysis and heating !
1134 : ! at twilight, Planet.Space Sci., v39, n5, pp. 671-683, 1991 (Appendix B) !
1135 : !=============================================================================!
1136 : ! PARAMETERS: !
1137 : ! NZ - INTEGER, number of specified altitude levels in the working (I) !
1138 : ! grid !
1139 : ! Z - REAL, specified altitude working grid (km) (I) !
1140 : ! ZEN - REAL, solar zenith angle (degrees) (I) !
1141 : ! DSDH - REAL, slant path of direct beam through each layer crossed (O) !
1142 : ! when travelling from the top of the atmosphere to layer i; !
1143 : ! DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 !
1144 : ! NID - INTEGER, number of layers crossed by the direct beam when (O) !
1145 : ! travelling from the top of the atmosphere to layer i; !
1146 : ! NID(i), i = 0..NZ-1 !
1147 : !=============================================================================!
1148 : ! EDIT HISTORY: !
1149 : ! Original: Taken By Doug Kinnison from Sasha Madronich, TUV Code, V4.1a, !
1150 : ! on 1/1/02 !
1151 : !=============================================================================!
1152 :
1153 : use physconst, only : rearth
1154 :
1155 : implicit none
1156 :
1157 : !------------------------------------------------------------------------------
1158 : ! ... Dummy arguments
1159 : !------------------------------------------------------------------------------
1160 : integer, intent(in) :: nlev ! number model vertical levels
1161 : integer, intent(out) :: nid(0:nlev) ! see above
1162 : real(r8), intent (in) :: zenith_angle ! zenith_angle
1163 : real(r8), intent (in) :: z(nlev) ! geometric altitude (km)
1164 : real(r8), intent (out) :: dsdh(0:nlev,nlev) ! see above
1165 :
1166 :
1167 : !------------------------------------------------------------------------------
1168 : ! ... Local variables
1169 : !------------------------------------------------------------------------------
1170 : real(r8) :: radius
1171 : real(r8) :: re
1172 : real(r8) :: zenrad
1173 : real(r8) :: rpsinz
1174 : real(r8) :: const0
1175 : real(r8) :: rj
1176 : real(r8) :: rjp1
1177 : real(r8) :: dsj
1178 : real(r8) :: dhj
1179 : real(r8) :: ga
1180 : real(r8) :: gb
1181 : real(r8) :: sm
1182 0 : real(r8) :: zd(0:nlev-1)
1183 :
1184 : integer :: i
1185 : integer :: j
1186 : integer :: k
1187 : integer :: id
1188 : integer :: nlayer
1189 :
1190 :
1191 0 : radius = rearth*1.e-3_r8 ! radius earth (km)
1192 :
1193 : !------------------------------------------------------------------------------
1194 : ! ... set zenith angle in radians
1195 : !------------------------------------------------------------------------------
1196 0 : zenrad = zenith_angle*d2r
1197 0 : const0 = sin( zenrad )
1198 :
1199 : !------------------------------------------------------------------------------
1200 : ! ... set number of layers:
1201 : !------------------------------------------------------------------------------
1202 0 : nlayer = nlev - 1
1203 :
1204 : !------------------------------------------------------------------------------
1205 : ! ... include the elevation above sea level to the radius of the earth:
1206 : !------------------------------------------------------------------------------
1207 0 : re = radius + z(nlev)
1208 :
1209 : !------------------------------------------------------------------------------
1210 : ! ... inverse coordinate of z
1211 : !------------------------------------------------------------------------------
1212 0 : do k = 0,nlayer
1213 0 : zd(k) = z(k+1) - z(nlev)
1214 : end do
1215 :
1216 : !------------------------------------------------------------------------------
1217 : ! ... initialize dsdh(i,j), nid(i)
1218 : !------------------------------------------------------------------------------
1219 0 : nid(:) = 0
1220 0 : do j = 1,nlev
1221 0 : dsdh(:,j) = 0._r8
1222 : end do
1223 :
1224 : !------------------------------------------------------------------------------
1225 : ! ... calculate ds/dh of every layer
1226 : !------------------------------------------------------------------------------
1227 0 : do i = 0,nlayer
1228 0 : rpsinz = (re + zd(i)) * const0
1229 0 : if( zenith_angle <= 90._r8 .or. rpsinz >= re ) then
1230 : !------------------------------------------------------------------------------
1231 : ! Find index of layer in which the screening height lies
1232 : !------------------------------------------------------------------------------
1233 0 : id = i
1234 0 : if( zenith_angle > 90._r8 ) then
1235 0 : do j = 1,nlayer
1236 0 : if( rpsinz < (zd(j-1) + re) .and. rpsinz >= (zd(j) + re) ) then
1237 : id = j
1238 : exit
1239 : end if
1240 : end do
1241 : end if
1242 :
1243 0 : do j = 1,id
1244 0 : sm = 1._r8
1245 0 : if( j == id .and. id == i .and. zenith_angle > 90._r8 ) then
1246 0 : sm = -1._r8
1247 : end if
1248 0 : rj = re + zd(j-1)
1249 0 : rjp1 = re + zd(j)
1250 0 : dhj = zd(j-1) - zd(j)
1251 0 : ga = max( rj*rj - rpsinz*rpsinz,0._r8 )
1252 0 : gb = max( rjp1*rjp1 - rpsinz*rpsinz,0._r8 )
1253 0 : if( id > i .and. j == id ) then
1254 0 : dsj = sqrt( ga )
1255 : else
1256 0 : dsj = sqrt( ga ) - sm*sqrt( gb )
1257 : end if
1258 0 : dsdh(i,j) = dsj / dhj
1259 : end do
1260 0 : nid(i) = id
1261 : else
1262 0 : nid(i) = -1
1263 : end if
1264 : end do
1265 :
1266 0 : end subroutine sphers
1267 :
1268 0 : subroutine slant_col( nlev, delz, dsdh, nid, absden, scol )
1269 : !=============================================================================!
1270 : ! PURPOSE: !
1271 : ! Derive Column
1272 : !=============================================================================!
1273 : ! PARAMETERS: !
1274 : ! NLEV - INTEGER, number of specified altitude levels in the working (I) !
1275 : ! grid !
1276 : ! DELZ - REAL, specified altitude working grid (km) (I) !
1277 : ! DSDH - REAL, slant path of direct beam through each layer crossed (O) !
1278 : ! when travelling from the top of the atmosphere to layer i; !
1279 : ! DSDH(i,j), i = 0..NZ-1, j = 1..NZ-1 !
1280 : ! NID - INTEGER, number of layers crossed by the direct beam when (O) !
1281 : ! travelling from the top of the atmosphere to layer i; !
1282 : ! NID(i), i = 0..NZ-1 !
1283 : ! specified altitude at each specified wavelength !
1284 : ! absden - REAL, absorber concentration, molecules cm-3 !
1285 : ! SCOL - REAL, absorber Slant Column, molecules cm-2 !
1286 : !=============================================================================!
1287 : ! EDIT HISTORY: !
1288 : ! 09/01 Read in profile from an input file, DEK !
1289 : ! 01/02 Taken from Sasha Madronich's TUV code !
1290 : !=============================================================================!
1291 :
1292 : implicit none
1293 :
1294 : !------------------------------------------------------------------------------
1295 : ! ... Dummy arguments
1296 : !------------------------------------------------------------------------------
1297 : integer, intent(in) :: nlev
1298 : integer, intent(in) :: nid(0:nlev) ! see above
1299 : real(r8), intent(in) :: delz(nlev) ! layer thickness (cm)
1300 : real(r8), intent(in) :: dsdh(0:nlev,nlev) ! see above
1301 : real(r8), intent(in) :: absden(nlev) ! absorber concentration (molec. cm-3)
1302 : real(r8), intent(out) :: scol(nlev) ! absorber Slant Column (molec. cm-2)
1303 :
1304 : !------------------------------------------------------------------------------
1305 : ! ... Local variables
1306 : !------------------------------------------------------------------------------
1307 : real(r8), parameter :: largest = 1.e+36_r8
1308 :
1309 : real(r8) :: sum
1310 : real(r8) :: hscale
1311 : real(r8) :: numer, denom
1312 0 : real(r8) :: cz(nlev)
1313 :
1314 : integer :: id
1315 : integer :: j
1316 : integer :: k
1317 :
1318 : !------------------------------------------------------------------------------
1319 : ! ... compute column increments (logarithmic integrals)
1320 : !------------------------------------------------------------------------------
1321 0 : do k = 1,nlev-1
1322 0 : if( absden(k) /= 0._r8 .and. absden(k+1) /= 0._r8 ) then
1323 0 : cz(nlev-k) = (absden(k) - absden(k+1))/log( absden(k)/absden(k+1) ) * delz(k)
1324 : else
1325 0 : cz(nlev-k) = .5_r8*(absden(k) + absden(k+1)) * delz(k)
1326 : end if
1327 : end do
1328 :
1329 : !------------------------------------------------------------------------------
1330 : ! ... Include exponential tail integral from infinity to model top
1331 : ! specify scale height near top of data.For WACCM-X model, scale
1332 : ! height needs to be increased for higher model top
1333 : !------------------------------------------------------------------------------
1334 0 : if (nlev==pver) then
1335 0 : if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then
1336 : hscale = 20.e5_r8
1337 : else
1338 0 : hscale = 10.e5_r8
1339 : endif
1340 0 : cz(nlev-1) = cz(nlev-1) + hscale * absden(1)
1341 : endif
1342 :
1343 : !------------------------------------------------------------------------------
1344 : ! ... Calculate vertical and slant column from each level:
1345 : ! work downward
1346 : !------------------------------------------------------------------------------
1347 0 : do id = 0,nlev-1
1348 0 : sum = 0._r8
1349 0 : if( nid(id) >= 0 ) then
1350 : !------------------------------------------------------------------------------
1351 : ! ... Single pass layers:
1352 : !------------------------------------------------------------------------------
1353 0 : do j = 1, min(nid(id), id)
1354 0 : sum = sum + cz(nlev-j)*dsdh(id,j)
1355 : end do
1356 : !------------------------------------------------------------------------------
1357 : ! ... Double pass layers:
1358 : !------------------------------------------------------------------------------
1359 0 : do j = min(nid(id),id)+1, nid(id)
1360 0 : sum = sum + 2._r8*cz(nlev-j)*dsdh(id,j)
1361 : end do
1362 : else
1363 : sum = largest
1364 : end if
1365 0 : scol(nlev-id) = sum
1366 : end do
1367 0 : scol(nlev) = .95_r8*scol(nlev-1)
1368 :
1369 0 : end subroutine slant_col
1370 :
1371 0 : subroutine lymana( nlev, o2scol, rm, ro2 )
1372 : !-----------------------------------------------------------------------------!
1373 : ! PURPOSE: !
1374 : ! Calculate the effective absorption cross section of O2 in the Lyman-Alpha !
1375 : ! bands and an effective O2 optical depth at all altitudes. Parameterized !
1376 : ! after: Chabrillat, S., and G. Kockarts, Simple parameterization of the !
1377 : ! absorption of the solar Lyman-Alpha line, Geophysical Research Letters, !
1378 : ! Vol.24, No.21, pp 2659-2662, 1997. doi:10.1029/97GL52690 (note there is a !
1379 : ! correction to this paper - the table was missing minuses in the exponents)!
1380 : !-----------------------------------------------------------------------------!
1381 : ! PARAMETERS: !
1382 : ! nz - INTEGER, number of specified altitude levels in the working (I) !
1383 : ! grid !
1384 : ! o2scol - REAL, slant overhead O2 column (molec/cc) at each specified (I) !
1385 : ! altitude !
1386 : ! dto2la - REAL, optical depth due to O2 absorption at each specified (O) !
1387 : ! vertical layer !
1388 : ! xso2la - REAL, molecular absorption cross section in LA bands (O) !
1389 : !-----------------------------------------------------------------------------!
1390 : ! EDIT HISTORY: !
1391 : ! 01/15/2002 Taken from Sasha Madronich's TUV Version 4.1a, Doug Kinnison !
1392 : ! 01/15/2002 Upgraded to F90, DK !
1393 : !-----------------------------------------------------------------------------!
1394 :
1395 : implicit none
1396 :
1397 : !------------------------------------------------------------------------------
1398 : ! ... Dummy arguments
1399 : !------------------------------------------------------------------------------
1400 : integer, intent(in) :: nlev
1401 : real(r8), intent(in) :: o2scol(nlev)
1402 : real(r8), intent(out) :: ro2(nlev)
1403 : real(r8), intent(out) :: rm(nlev)
1404 :
1405 : !------------------------------------------------------------------------------
1406 : ! ... Local variables
1407 : !------------------------------------------------------------------------------
1408 : real(r8),save :: b(3)
1409 : real(r8),save :: c(3)
1410 : real(r8),save :: d(3)
1411 : real(r8),save :: e(3)
1412 :
1413 : data b / 6.8431e-01_r8, 2.29841e-01_r8, 8.65412e-02_r8 /, &
1414 : c / 8.22114e-21_r8, 1.77556e-20_r8, 8.22112e-21_r8 /, &
1415 : d / 6.0073e-21_r8, 4.28569e-21_r8, 1.28059e-20_r8 /, &
1416 : e / 8.21666e-21_r8, 1.63296e-20_r8, 4.85121e-17_r8 /
1417 :
1418 : integer :: i, k
1419 : real(r8) :: wrk, term
1420 :
1421 : !------------------------------------------------------------------------------
1422 : ! ... Calculate reduction factors at every altitude
1423 : !------------------------------------------------------------------------------
1424 0 : do k = 1,nlev
1425 : wrk = 0._r8
1426 0 : do i = 1,2 ! pc Dan Marsh
1427 0 : term = e(i)*o2scol(k)
1428 0 : if( term < 100._r8 ) then
1429 0 : wrk = wrk + d(i) * exp( -term )
1430 : end if
1431 : end do
1432 0 : ro2(k) = wrk
1433 0 : wrk = 0._r8
1434 0 : do i = 1,3
1435 0 : term = c(i)*o2scol(k)
1436 0 : if( term < 100._r8 ) then
1437 0 : wrk = wrk + b(i) * exp( -term )
1438 : end if
1439 : end do
1440 0 : rm(k) = wrk
1441 : end do
1442 :
1443 0 : end subroutine lymana
1444 :
1445 0 : subroutine calc_o2srb( nlev, nid, o2col, tlev, tsrb, xscho2 )
1446 : !-----------------------------------------------------------------------------!
1447 : ! PURPOSE: !
1448 : ! Calculate the equivalent absorption cross section of O2 in the SR bands. !
1449 : ! The algorithm is based on parameterization of G.A. Koppers, and !
1450 : ! D.P. Murtagh [ref. Ann.Geophys., 14 68-79, 1996] !
1451 : ! Final values do include effects from the Herzberg continuum. !
1452 : !-----------------------------------------------------------------------------!
1453 : ! PARAMETERS: !
1454 : ! NZ - INTEGER, number of specified altitude levels in the working (I) !
1455 : ! grid !
1456 : ! O2COL - REAL, slant overhead O2 column (molec/cc) at each specified (I) !
1457 : ! altitude !
1458 : ! TLEV - tmeperature at each level (I) !
1459 : ! TSRB - REAL, transmission for the SRB !
1460 : ! XSCHO2 - REAL, molecular absorption cross section in SR bands at (O) !
1461 : ! each specified wavelength. Includes Herzberg continuum !
1462 : !-----------------------------------------------------------------------------!
1463 : ! EDIT HISTORY: Taken from TUV, 1/17/2002 !
1464 : ! This code was supplied to TUV by Dan Marsh. !
1465 : !-----------------------------------------------------------------------------!
1466 :
1467 : implicit none
1468 :
1469 : !------------------------------------------------------------------------------
1470 : ! ... Dummy arguments
1471 : !------------------------------------------------------------------------------
1472 : integer, intent(in) :: nlev
1473 : integer, intent(in) :: nid(0:nlev)
1474 : real(r8), intent (in) :: o2col(nlev)
1475 : real(r8), intent (in) :: tlev(nlev)
1476 : real(r8), intent (out) :: tsrb(nlev,nsrbtuv)
1477 : real(r8), intent (out) :: xscho2(nlev,nsrbtuv)
1478 :
1479 : !------------------------------------------------------------------------------
1480 : ! ... Local variables
1481 : !------------------------------------------------------------------------------
1482 : integer :: i, k, ktop, ktop1, kbot
1483 : real(r8) :: x, dto2
1484 : real(r8) :: den, num
1485 : real(r8) :: term1, term2
1486 0 : real(r8) :: dtsrb(nlev)
1487 : real(r8) :: tsrb_rev(nlev,nsrbtuv)
1488 : real(r8) :: xs(nsrbtuv)
1489 :
1490 : !------------------------------------------------------------------------------
1491 : ! ... Calculate cross sections
1492 : !------------------------------------------------------------------------------
1493 0 : ktop = nlev
1494 0 : kbot = 0
1495 :
1496 0 : do k = 1,nlev
1497 0 : x = log( o2col(k) )
1498 0 : if( x >= lno2_llimit .and. x <= lno2_ulimit ) then
1499 : call effxs( x, tlev(k), xs )
1500 0 : xscho2(k,:) = xs(:)
1501 0 : else if( x < lno2_llimit ) then
1502 0 : ktop1 = k-1
1503 0 : ktop = min( ktop1,ktop )
1504 0 : else if( x > lno2_ulimit ) then
1505 0 : kbot = k
1506 : end if
1507 : end do
1508 :
1509 0 : if ( kbot == nlev ) then
1510 0 : tsrb(:,:) = 0._r8
1511 0 : xscho2(:,:) = 0._r8
1512 : return
1513 : endif
1514 : !------------------------------------------------------
1515 : ! ... Fill in cross section where X is out of range
1516 : ! by repeating edge table values
1517 : !-------------------------------------------------------
1518 0 : if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then
1519 :
1520 : ! Need to be careful with nlev values for kbot and ktop.
1521 : ! This was handled by Hanli Liu fix.
1522 0 : if ( kbot < nlev ) then
1523 0 : do k = 1,kbot
1524 0 : xscho2(k,:) = xscho2(kbot+1,:)
1525 : end do
1526 0 : if (ktop < nlev) then
1527 0 : do k = ktop+1,nlev
1528 0 : xscho2(k,:) = xscho2(ktop,:)
1529 : end do
1530 : else
1531 0 : xscho2(nlev,:) = 2.0e-19_r8
1532 : endif
1533 : else
1534 0 : do k = 1,kbot
1535 0 : xscho2(k,:) = 2.0e-19_r8
1536 : enddo
1537 : endif
1538 :
1539 : else
1540 0 : do k = 1,kbot
1541 0 : xscho2(k,:) = xscho2(kbot+1,:)
1542 : end do
1543 0 : do k = ktop+1,nlev
1544 0 : xscho2(k,:) = xscho2(ktop,:)
1545 : end do
1546 : endif
1547 :
1548 : !-------------------------------------------------------
1549 : ! ... Calculate incremental optical depths
1550 : !-------------------------------------------------------
1551 0 : do i = 1,nsrbtuv
1552 0 : do k = 1,nlev-1
1553 0 : if( nid(nlev-k) /= -1 ) then
1554 : !-------------------------------------------------------
1555 : ! ... Calculate an optical depth weighted by density
1556 : !-------------------------------------------------------
1557 0 : num = xscho2(k+1,i)*o2col(k+1) - xscho2(k,i)*o2col(k)
1558 0 : if( num == 0._r8 ) then
1559 0 : write(iulog,*) 'calc_o2srb : o2col(k:k+1),xscho2(k:k+1,i) = ',o2col(k:k+1),xscho2(k:k+1,i),' @ i,k = ',i,k
1560 : end if
1561 0 : term1 = log( xscho2(k+1,i)/xscho2(k,i) )
1562 0 : term2 = log( o2col(k+1)/o2col(k) )
1563 0 : if( term2 == 0._r8 ) then
1564 0 : write(iulog,*) 'calc_o2srb : o2col(k:k+1),xscho2(k:k+1,i) = ',o2col(k:k+1),xscho2(k:k+1,i),' @ i,k = ',i,k
1565 0 : call endrun
1566 : end if
1567 0 : den = 1._r8 + log( xscho2(k+1,i)/xscho2(k,i) )/log( o2col(k+1)/o2col(k) )
1568 0 : dto2 = abs(num/den)
1569 0 : if( dto2 < 100._r8 ) then
1570 0 : dtsrb(k) = exp( -dto2 )
1571 : else
1572 0 : dtsrb(k) = 0._r8
1573 : end if
1574 : else
1575 0 : dtsrb(k) = 0._r8
1576 : end if
1577 : end do
1578 : !-----------------------------------------------
1579 : ! ... Calculate Transmission for SRB
1580 : !-----------------------------------------------
1581 0 : if (nlev==pver) then ! waccm
1582 0 : tsrb(nlev,i) = 1._r8
1583 0 : do k = nlev-1,1,-1
1584 0 : tsrb(k,i) = tsrb(k+1,i)*dtsrb(k)
1585 : end do
1586 : else ! cam-chem
1587 0 : tsrb(nlev,i) = exp(-xscho2(nlev,i)*o2col(nlev))
1588 0 : do k = nlev-1,1,-1
1589 0 : tsrb(k,i) = tsrb(k+1,i)*dtsrb(k)
1590 : end do
1591 : endif
1592 :
1593 : end do
1594 :
1595 : end subroutine calc_o2srb
1596 :
1597 0 : subroutine effxs( x, t, xs )
1598 : !-------------------------------------------------------------
1599 : ! Subroutine for evaluating the effective cross section
1600 : ! of O2 in the Schumann-Runge bands using parameterization
1601 : ! of G.A. Koppers, and D.P. Murtagh [ref. Ann.Geophys., 14
1602 : ! 68-79, 1996]
1603 : !
1604 : ! method:
1605 : ! ln(xs) = A(X)[T-220]+B(X)
1606 : ! X = log of slant column of O2
1607 : ! A,B are calculated from Chebyshev polynomial coeffs
1608 : ! AC and BC using Clenshaw summation algorithm within
1609 : ! the interval 38<ln(NO2)<56.
1610 : !
1611 : ! Revision History:
1612 : !
1613 : ! drm 2/97 initial coding
1614 : !-------------------------------------------------------------
1615 :
1616 : implicit none
1617 :
1618 : !-------------------------------------------------------------
1619 : ! ... Dummy arguments
1620 : !-------------------------------------------------------------
1621 : real(r8), intent(in) :: x
1622 : real(r8), intent(in) :: t
1623 : real(r8), intent(out) :: xs(nsrbtuv)
1624 :
1625 : !-------------------------------------------------------------
1626 : ! ... Local variables
1627 : !-------------------------------------------------------------
1628 : real(r8) :: a(nsrbtuv), b(nsrbtuv)
1629 :
1630 0 : call calc_params( x, a, b )
1631 :
1632 0 : xs(:) = exp( a(:)*(t - 220._r8) + b(:) )
1633 :
1634 0 : end subroutine effxs
1635 :
1636 :
1637 0 : subroutine calc_params( x, a, b )
1638 : !-------------------------------------------------------------
1639 : ! calculates coefficients (A,B), used in calculating the
1640 : ! effective cross section, for 17 wavelength intervals
1641 : ! as a function of log O2 column density (X)
1642 : ! Wavelength intervals are defined in WMO1985
1643 : !-------------------------------------------------------------
1644 :
1645 : !-------------------------------------------------------------
1646 : ! ... Dummy arguments
1647 : !-------------------------------------------------------------
1648 : real(r8), intent(in) :: x
1649 : real(r8), intent(out) :: a(nsrbtuv), b(nsrbtuv)
1650 :
1651 : !-------------------------------------------------------------
1652 : ! ... Local variables
1653 : !-------------------------------------------------------------
1654 : integer :: i
1655 :
1656 0 : if (x<lno2_llimit .or. x>lno2_ulimit) then
1657 0 : call endrun('mo_jshort::calc_params of O2 abs xs: x is not in the valid range. ')
1658 : end if
1659 :
1660 : !-------------------------------------------------------------
1661 : ! ... evaluate at each wavelength
1662 : ! for a set of 20 Chebyshev coeficients
1663 : !-------------------------------------------------------------
1664 0 : do i = 1,nsrbtuv
1665 0 : a(i) = evalchebpoly( ac(:,i), x )
1666 0 : b(i) = evalchebpoly( bc(:,i), x )
1667 : end do
1668 :
1669 : contains
1670 :
1671 : ! Use Clenshaw summation algorithm to evaluate Chebyshev polynomial at point
1672 : ! [pnt - (lno2_ulimit + lno2_llimit)/2]/[(lno2_ulimit - lno2_llimit)/2]
1673 : ! given coefficients coefs within limits lim1 and lim2
1674 0 : pure function evalchebpoly( coefs, pnt ) result(cval)
1675 : real(r8), intent(in) :: coefs(:)
1676 : real(r8), intent(in) :: pnt
1677 :
1678 : real(r8) :: cval
1679 : real(r8) :: fac(2)
1680 : real(r8) :: csum(2) ! Clenshaw summation
1681 : integer :: ndx
1682 : integer :: ncoef
1683 :
1684 0 : ncoef = size(coefs)
1685 :
1686 0 : fac(1) = (2._r8*pnt-lno2_llimit-lno2_ulimit)/(lno2_ulimit-lno2_llimit)
1687 0 : fac(2) = 2._r8*fac(1)
1688 :
1689 : ! Clenshaw recurrence summation
1690 0 : csum(:) = 0.0_r8
1691 0 : do ndx = ncoef, 2, -1
1692 0 : cval = csum(1)
1693 0 : csum(1) = fac(2)*csum(1) - csum(2) + coefs(ndx)
1694 0 : csum(2) = cval
1695 : end do
1696 :
1697 0 : cval = fac(1)*csum(1) - csum(2) + 0.5_r8*coefs(1)
1698 :
1699 0 : end function evalchebpoly
1700 :
1701 : end subroutine calc_params
1702 :
1703 0 : subroutine calc_jno( nlev, etfphot_ms93, n2cc, o2scol, o3scol, &
1704 0 : noscol, jno )
1705 : !-----------------------------------------------------------------------------!
1706 : ! PURPOSE: !
1707 : ! Compute the total photolytic rate constant for NO in the SR bands !
1708 : ! - following the approach of Minshwanner and Siskind, JGR, !
1709 : ! 98, D11, 20401-20412, 1993. !
1710 : ! !
1711 : !-----------------------------------------------------------------------------!
1712 : ! PARAMETERS: !
1713 : ! NZ - INTEGER, number of specified altitude levels !
1714 : ! !
1715 : ! etfphot_ms93 - Extraterrestrial Flux, within the MS 1993 Grid !
1716 : ! units of photons cm-2 sec-1 nm-1 !
1717 : ! n2cc - N2 conc (molecules cm-3) !
1718 : ! o3scol - Ozone Slant Column (molecules cm-2) !
1719 : ! o2scol - Oxygen Slant Column (molecules cm-2) !
1720 : ! noscol - Nitric Oxide Slant Column(molecules cm-2) !
1721 : ! !
1722 : ! LOCAL VARIABLES: !
1723 : ! tauo3 - Transmission factor in the Hartley Band of O3 !
1724 : ! etfphot_ms93 - Solar Irr. on Minschwaner and Siskind 1993 (MS93) Grid !
1725 : ! xs_o3ms93 - O3 cross section on the MS93 Grid !
1726 : ! !
1727 : ! OUTPUT VARIABLES: !
1728 : ! jno - photolytic rate constant !
1729 : ! each specified altitude !
1730 : ! !
1731 : !-----------------------------------------------------------------------------!
1732 : ! EDIT HISTORY: !
1733 : ! 08/01 Created, Doug Kinnison, NCAR, ACD !
1734 : !-----------------------------------------------------------------------------!
1735 :
1736 : implicit none
1737 :
1738 : !------------------------------------------------------------------------------
1739 : ! ... Dummy arguments
1740 : !------------------------------------------------------------------------------
1741 : integer, intent(in) :: nlev
1742 : real(r8), intent(in) :: etfphot_ms93(num_ms93tuv)
1743 : real(r8), intent(in) :: n2cc(nlev)
1744 : real(r8), intent(in) :: o3scol(nlev)
1745 : real(r8), intent(in) :: o2scol(nlev)
1746 : real(r8), intent(in) :: noscol(nlev)
1747 : real(r8), intent(out) :: jno(nlev)
1748 :
1749 : !------------------------------------------------------------------------------
1750 : ! ... Local variables
1751 : !------------------------------------------------------------------------------
1752 : integer :: i, iw, lev
1753 : real(r8) :: jno50
1754 : real(r8) :: jno90
1755 : real(r8) :: jno100
1756 0 : real(r8) :: tauo3(nlev,num_ms93tuv)
1757 :
1758 : !------------------------------------------------------------------------------
1759 : ! ... O3 SRB Cross Sections from WMO 1985, interpolated onto MS, 1993 grid
1760 : !------------------------------------------------------------------------------
1761 : real(r8), save :: xso3_ms93(num_ms93tuv) = (/ 7.3307600e-19_r8, 6.9660105E-19_r8, 5.9257699E-19_r8, 4.8372219E-19_r8 /)
1762 :
1763 : !------------------------------------------------------------------------------
1764 : ! ... delta wavelength of the MS, 1993 grid
1765 : !------------------------------------------------------------------------------
1766 : real(r8), save :: wlintv_ms93(num_ms93tuv) = (/ 1.50_r8, 1.50_r8, 5.6_r8, 2.3_r8 /)
1767 :
1768 : !------------------------------------------------------------------------------
1769 : ! ... O2 SRB Cross Sections for the six ODF regions, MS, 1993
1770 : !------------------------------------------------------------------------------
1771 : real(r8), save :: cs250(6) = (/ 1.117e-23_r8, 2.447e-23_r8, 7.188e-23_r8, 3.042e-22_r8, 1.748e-21_r8, 1.112e-20_r8 /)
1772 : real(r8), save :: cs290(6) = (/ 1.350e-22_r8, 2.991e-22_r8, 7.334e-22_r8, 3.074e-21_r8, 1.689e-20_r8, 1.658e-19_r8 /)
1773 : real(r8), save :: cs2100(6) = (/ 2.968e-22_r8, 5.831e-22_r8, 2.053e-21_r8, 8.192e-21_r8, 4.802e-20_r8, 2.655e-19_r8 /)
1774 :
1775 : !------------------------------------------------------------------------------
1776 : ! ... derive tauo3 for the three o2 srb
1777 : ! ... iw = 1,2, and 4 are used below for jno
1778 : !------------------------------------------------------------------------------
1779 0 : do iw = 1,num_ms93tuv
1780 0 : tauo3(:,iw) = exp( -xso3_ms93(iw)*o3scol(:) )
1781 : end do
1782 :
1783 : !------------------------------------------------------------------------------
1784 : ! ... Call PJNO Function to derive SR Band JNO contributions
1785 : ! Called in order of wavelength interval (shortest firs)
1786 : !------------------------------------------------------------------------------
1787 0 : do lev = 1,nlev
1788 0 : jno100 = pjno( 1, cs2100, wtno100, csno100 )
1789 0 : jno90 = pjno( 2, cs290, wtno90, csno90 )
1790 0 : jno50 = pjno( 4, cs250, wtno50, csno50 )
1791 0 : jno(lev) = jno50 + jno90 + jno100
1792 : end do
1793 :
1794 : contains
1795 :
1796 0 : function pjno( w, cso2, wtno, csno )
1797 : !------------------------------------------------------------------------------
1798 : ! ... uses xsec at center of g subinterval for o2
1799 : ! uses mean values for no
1800 : !------------------------------------------------------------------------------
1801 : implicit none
1802 :
1803 : !------------------------------------------------------------------------------
1804 : ! ... parameters
1805 : !------------------------------------------------------------------------------
1806 : integer, parameter :: ngint = 6
1807 : integer, parameter :: nno = 2
1808 :
1809 : !----------------------------------------------------------------
1810 : ! ... Dummy arguments
1811 : !----------------------------------------------------------------
1812 : integer, intent(in) :: w
1813 : real(r8), intent(in) :: cso2(ngint)
1814 : real(r8), intent(in) :: csno(ngint,nno)
1815 : real(r8), intent(in) :: wtno(ngint,nno)
1816 :
1817 : !----------------------------------------------------------------
1818 : ! ... Function declarations
1819 : !----------------------------------------------------------------
1820 : real(r8) :: pjno
1821 :
1822 : !----------------------------------------------------------------
1823 : ! ... Local variables
1824 : !----------------------------------------------------------------
1825 : integer :: jj, i, k
1826 : real(r8) :: tauno
1827 : real(r8) :: transno
1828 : real(r8) :: transo2
1829 : real(r8) :: tauo2
1830 : real(r8) :: jno
1831 : real(r8) :: jno1
1832 :
1833 : !----------------------------------------------------------------
1834 : ! ... derive the photolysis frequency for no within a given
1835 : ! srb (i.e., 5-0, 9-0, 10-0)
1836 : !----------------------------------------------------------------
1837 0 : jno = 0._r8
1838 0 : do k = 1,ngint
1839 0 : tauo2 = o2scol(lev) * cso2(k)
1840 0 : if( tauo2 < 50._r8 ) then
1841 0 : transo2 = exp( -tauo2 )
1842 : else
1843 : transo2 = 0._r8
1844 : end if
1845 0 : jno1 = 0._r8
1846 0 : do jj = 1,nno
1847 0 : tauno = noscol(lev)*csno(k,jj)
1848 0 : if( tauno < 50._r8 ) then
1849 0 : transno = exp( -tauno )
1850 : else
1851 : transno = 0._r8
1852 : end if
1853 0 : jno1 = jno1 + csno(k,jj) * wtno(k,jj) * transno
1854 : end do
1855 0 : jno = jno + jno1*transo2
1856 : end do
1857 :
1858 0 : pjno = wlintv_ms93(w)*etfphot_ms93(w)*tauo3(lev,w)*jno
1859 :
1860 : !----------------------------------------------------------------
1861 : ! ... correct for the predissociation of the deltq 1-0
1862 : ! transition in the srb (5-0)
1863 : !----------------------------------------------------------------
1864 0 : if( w == 4 ) then
1865 0 : pjno = 1.65e9_r8/(5.1e7_r8 + 1.65e9_r8 + (1.5e-9_r8*n2cc(nlev-lev+1)))*pjno
1866 : end if
1867 :
1868 0 : end function pjno
1869 :
1870 : end subroutine calc_jno
1871 :
1872 : end module mo_jshort
|
