Line data Source code
1 : module mo_jeuv
2 :
3 : use shr_kind_mod, only : r8 => shr_kind_r8
4 : use cam_abortutils, only : endrun
5 : use cam_logfile, only : iulog
6 : use euvac, only : euvac_etf
7 : use solar_euv_data, only : solar_euv_data_etf, solar_euv_data_active
8 : use spmd_utils, only : masterproc, masterprocid, mpicom, mpi_real8
9 :
10 : implicit none
11 :
12 : private
13 : public :: jeuv_init, jeuv, heuv, neuv
14 : public :: nIonRates
15 :
16 : save
17 :
18 : !------------------------------------------------------------------------------
19 : ! define EUV photolysis cross sections, branching ratios,
20 : ! wavelength parameters,etc
21 : !------------------------------------------------------------------------------
22 : integer, parameter :: neuv = 26 ! number of photolysis/ionization reactions
23 : integer, parameter :: nIonRates = 11 ! number of photo-ionizations rates needed for waccmx
24 : integer, parameter :: nmaj = 3 ! number of major neutral species (O,O2,N2)
25 : integer, parameter :: nspecies = 5 ! number of neutral species(O,N,O2,N2,CO2)
26 : integer, parameter :: nstat = 6 ! maximum number of ionization/dissociation
27 : ! /excitation states for each speies
28 : integer, parameter :: lmax = 23 ! number of wavelength bins in EUV
29 :
30 : real(r8), parameter :: heat_eff_fac = .08_r8 ! heating efficiency factor
31 : ! Solar EUV direct heating was increased from 5% to 8% to bring it closer to the
32 : ! TIE-GCM value (5% applied twice for a total of 10%) -- Hanli Liu & Stan Solomon
33 :
34 : real(r8), parameter :: hc = 6.62608e-34_r8 * 2.9979e8_r8 / 1.e-9_r8
35 :
36 : real(r8) :: sigabs(lmax,nspecies) ! absorption cross sections of major species
37 : real(r8) :: branch_p(lmax,nstat,nmaj) = 0._r8 ! branching ratios for photoionization/dissociation
38 : real(r8) :: branch_e(lmax,nstat,nmaj) = 0._r8 ! branching ratios for photoelectron ionization/dissociation/excitation
39 : real(r8) :: energy(lmax) ! solar energy
40 :
41 : real(r8), pointer :: solar_etf(:)
42 : logical :: do_heating(13)
43 :
44 : contains
45 :
46 1536 : subroutine jeuv_init (photon_file, electron_file, indexer)
47 : !==============================================================================
48 : ! Purpose:
49 : ! read tabulated data:
50 : ! (1) thermosphere neutral species' absorption cross sections,
51 : ! photoionization/dissociation branching ratios
52 : ! (2) read photoelectron enhancement factor, photoelectron ionization/
53 : ! dissociation/excitation branching ratios
54 : ! (3) read solar flux
55 : !==============================================================================
56 :
57 : use units, only : getunit, freeunit
58 : use ioFileMod, only : getfil
59 : use mo_chem_utls, only : get_rxt_ndx
60 :
61 : implicit none
62 :
63 : character(len=*), intent(in) :: photon_file
64 : character(len=*), intent(in) :: electron_file
65 : integer, optional,intent(inout) :: indexer(:)
66 :
67 : !------------------------------------------------------------------------------
68 : ! local variables
69 : !------------------------------------------------------------------------------
70 : integer :: i, j, m ! loop indicies
71 : integer :: unit ! fortran i/o unit number
72 : integer :: istat ! file op status
73 : real(r8) :: waves(lmax) ! wavelength array for short bound of bins (A)
74 : real(r8) :: wavel(lmax) ! wavelength array for long bound of bins (A)
75 : real(r8) :: wc(lmax) ! wavelength bin center (nm)
76 : character(len=200) :: str,fmt ! string for comments in data file
77 : character(len=256) :: locfn
78 :
79 : integer :: jeuv_1_ndx, ierr, jndx
80 : character(len=2) :: mstring
81 : character(len=7) :: jstring
82 : logical :: do_jeuv
83 :
84 1536 : do_jeuv=.false.
85 1536 : do_heating=.false.
86 :
87 41472 : do m = 1,neuv
88 39936 : write ( mstring, '(I2)' ) m
89 39936 : jstring = 'jeuv_'//trim(adjustl(mstring))
90 39936 : jndx = get_rxt_ndx( jstring )
91 39936 : if (jndx>0) then
92 0 : do_jeuv=.true.
93 0 : indexer(m) = jndx
94 : endif
95 41472 : if (jndx>0 .and. m<14) then
96 0 : do_heating(m) = .true.
97 : endif
98 : enddo
99 :
100 1536 : if (.not.do_jeuv) return
101 :
102 0 : if (solar_euv_data_active) then
103 0 : solar_etf => solar_euv_data_etf
104 : else
105 0 : solar_etf => euvac_etf
106 : endif
107 :
108 0 : if ( size(solar_etf) /= lmax ) then
109 0 : write(iulog,*) 'jeuv_init: the size of solar_etf is incorrect '
110 0 : write(iulog,*) ' ... size(solar_etf) = ',size(solar_etf)
111 0 : write(iulog,*) ' .............. lmax = ',lmax
112 0 : call endrun('jeuv_init: the size of solar_etf is incorrect ')
113 : endif
114 :
115 : !------------------------------------------------------------------------------
116 : ! read from data file the absorption cross sections for neutral species,
117 : ! braching ratios for photoionization/dissociation, and braching ratios
118 : ! for photoelectron ionization/dissociation/excitation
119 : !------------------------------------------------------------------------------
120 0 : master: if (masterproc) then ! read in ascii data only on master task then b-cast
121 : !------------------------------------------------------------------------------
122 : ! read neutral species' absorption cross section and
123 : ! photoionization/dissociation branching ratio
124 : !------------------------------------------------------------------------------
125 0 : unit = getunit()
126 0 : call getfil( photon_file, locfn, 0 )
127 0 : open( unit, file = trim(locfn), status='UNKNOWN', iostat=istat )
128 0 : if( istat /= 0 ) then
129 0 : write(iulog,*) 'jeuv_init: failed to open ',trim(locfn),'; error = ',istat
130 0 : call endrun
131 : end if
132 : !------------------------------------------------------------------------------
133 : ! read O
134 : !------------------------------------------------------------------------------
135 0 : read(unit,*,iostat=istat) str
136 0 : if( istat /= 0 ) then
137 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
138 0 : call endrun
139 : end if
140 0 : read(unit,*,iostat=istat) str
141 0 : if( istat /= 0 ) then
142 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
143 0 : call endrun
144 : end if
145 0 : fmt = "(f7.2,2x,f7.2,2x,f9.3,4(2x,f7.3))"
146 0 : do i = 1,lmax
147 0 : read(unit,fmt,iostat=istat) waves(i), wavel(i), sigabs(i,1), (branch_p(i,j,1),j=1,4)
148 0 : if( istat /= 0 ) then
149 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
150 0 : call endrun
151 : end if
152 : end do
153 : !------------------------------------------------------------------------------
154 : ! read O2
155 : !------------------------------------------------------------------------------
156 0 : read(unit,*,iostat=istat) str
157 0 : if( istat /= 0 ) then
158 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
159 0 : call endrun
160 : end if
161 0 : read(unit,*,iostat=istat) str
162 0 : if( istat /= 0 ) then
163 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
164 0 : call endrun
165 : end if
166 0 : fmt = "(f7.2,2x,f7.2,2x,f9.3,5(2x,f7.3))"
167 0 : do i = 1,lmax
168 0 : read(unit,fmt,iostat=istat) waves(i), wavel(i), sigabs(i,2), (branch_p(i,j,2),j=1,5)
169 0 : if( istat /= 0 ) then
170 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
171 0 : call endrun
172 : end if
173 : end do
174 : !------------------------------------------------------------------------------
175 : ! read N2
176 : !------------------------------------------------------------------------------
177 0 : read(unit,*,iostat=istat) str
178 0 : if( istat /= 0 ) then
179 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
180 0 : call endrun
181 : end if
182 0 : read(unit,*,iostat=istat) str
183 0 : if( istat /= 0 ) then
184 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
185 0 : call endrun
186 : end if
187 0 : do i = 1,lmax
188 0 : read(unit,fmt,iostat=istat) waves(i), wavel(i), sigabs(i,3), (branch_p(i,j,3),j=1,5)
189 0 : if( istat /= 0 ) then
190 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
191 0 : call endrun
192 : end if
193 : end do
194 : !------------------------------------------------------------------------------
195 : ! read N
196 : !------------------------------------------------------------------------------
197 0 : read(unit,*,iostat=istat) str
198 0 : if( istat /= 0 ) then
199 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
200 0 : call endrun
201 : end if
202 0 : read(unit,*,iostat=istat) str
203 0 : if( istat /= 0 ) then
204 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
205 0 : call endrun
206 : end if
207 0 : fmt = "(f7.2,2x,f7.2,2x,f9.3)"
208 0 : do i = 1,lmax
209 0 : read(unit,fmt,iostat=istat) waves(i), wavel(i), sigabs(i,4)
210 0 : if( istat /= 0 ) then
211 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
212 0 : call endrun
213 : end if
214 : end do
215 :
216 : !------------------------------------------------------------------------------
217 : ! read CO2
218 : !------------------------------------------------------------------------------
219 0 : read(unit,*,iostat=istat) str
220 0 : if( istat /= 0 ) then
221 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
222 0 : call endrun
223 : end if
224 0 : read(unit,*,iostat=istat) str
225 0 : if( istat /= 0 ) then
226 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
227 0 : call endrun
228 : end if
229 0 : fmt = "(f7.2,2x,f7.2,2x,f9.3)"
230 0 : do i = 1,lmax
231 0 : read(unit,fmt,iostat=istat) waves(i), wavel(i), sigabs(i,5)
232 0 : if( istat /= 0 ) then
233 0 : write(iulog,*) 'jeuv_init: failed to read CO2 data ',trim(locfn),'; error = ',istat
234 0 : call endrun
235 : end if
236 : end do
237 :
238 0 : close( unit )
239 :
240 : !------------------------------------------------------------------------------
241 : ! unit transfer for absorption cross sections
242 : ! from Megabarns to cm^2
243 : !------------------------------------------------------------------------------
244 0 : sigabs(:,:) = sigabs(:,:)*1.e-18_r8
245 :
246 : !------------------------------------------------------------------------------
247 : ! read photoelectron ionization/dissociation/excitation branching ratio
248 : !------------------------------------------------------------------------------
249 0 : call getfil( electron_file, locfn, 0 )
250 0 : open( unit, file = trim(locfn), status='UNKNOWN', iostat=istat )
251 0 : if( istat /= 0 ) then
252 0 : write(iulog,*) 'jeuv_init: failed to open ',trim(locfn),'; error = ',istat
253 0 : call endrun
254 : end if
255 : !------------------------------------------------------------------------------
256 : ! read O
257 : !------------------------------------------------------------------------------
258 0 : read(unit,*,iostat=istat) str
259 0 : if( istat /= 0 ) then
260 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
261 0 : call endrun
262 : end if
263 0 : read(unit,*,iostat=istat) str
264 0 : if( istat /= 0 ) then
265 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
266 0 : call endrun
267 : end if
268 0 : fmt="(f7.2,2x,f7.2,5(2x,f8.3))"
269 0 : do i = 1,lmax
270 0 : read(unit,fmt,iostat=istat) waves(i), wavel(i), (branch_e(i,j,1),j=1,5)
271 0 : if( istat /= 0 ) then
272 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
273 0 : call endrun
274 : end if
275 : end do
276 : !------------------------------------------------------------------------------
277 : ! read O2
278 : !------------------------------------------------------------------------------
279 0 : read(unit,*,iostat=istat) str
280 0 : if( istat /= 0 ) then
281 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
282 0 : call endrun
283 : end if
284 0 : read(unit,*,iostat=istat) str
285 0 : if( istat /= 0 ) then
286 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
287 0 : call endrun
288 : end if
289 0 : fmt = "(f7.2,2x,f7.2,6(2x,f8.3))"
290 0 : do i = 1,lmax
291 0 : read(unit,fmt,iostat=istat) waves(i), wavel(i), (branch_e(i,j,2),j=1,6)
292 0 : if( istat /= 0 ) then
293 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
294 0 : call endrun
295 : end if
296 : end do
297 : !------------------------------------------------------------------------------
298 : ! read N2
299 : !------------------------------------------------------------------------------
300 0 : read(unit,*,iostat=istat) str
301 0 : if( istat /= 0 ) then
302 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
303 0 : call endrun
304 : end if
305 0 : read(unit,*,iostat=istat) str
306 0 : if( istat /= 0 ) then
307 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
308 0 : call endrun
309 : end if
310 0 : do i = 1,lmax
311 0 : read(unit,fmt,iostat=istat) waves(i), wavel(i), (branch_e(i,j,3),j=1,6)
312 0 : if( istat /= 0 ) then
313 0 : write(iulog,*) 'jeuv_init: failed to read ',trim(locfn),'; error = ',istat
314 0 : call endrun
315 : end if
316 : end do
317 :
318 0 : close( unit )
319 :
320 0 : call freeunit( unit )
321 : end if master
322 :
323 0 : call mpi_bcast(sigabs, lmax*nspecies, mpi_real8, masterprocid, mpicom, ierr)
324 0 : call mpi_bcast(branch_p, lmax*nstat*nmaj, mpi_real8, masterprocid, mpicom, ierr)
325 0 : call mpi_bcast(branch_e, lmax*nstat*nmaj, mpi_real8, masterprocid, mpicom, ierr)
326 0 : call mpi_bcast(waves, lmax, mpi_real8, masterprocid, mpicom, ierr)
327 0 : call mpi_bcast(wavel, lmax, mpi_real8, masterprocid, mpicom, ierr)
328 :
329 0 : wc(:) = .1_r8*(waves(:) + wavel(:))*0.5_r8 ! A to nm
330 0 : energy(:) = heat_eff_fac*hc/wc(:)
331 :
332 : end subroutine jeuv_init
333 :
334 0 : subroutine jeuv( nlev, zen, occ, o2cc, n2cc, &
335 0 : zkm, euv_prates)
336 : !==============================================================================
337 : ! Purpose:
338 : ! Calculate euv photolysis/ionization rates (in s-1)
339 : !==============================================================================
340 : ! Arguments:
341 : ! nlev: number of model vertical levels
342 : ! zen: solar zenith angle in degree
343 : ! occ: atmoic oxygen number density (#/cm3)
344 : ! o2cc: molecular oxygen number density (#/cm3)
345 : ! n2cc: molecular nitrogen number density (#/cm3)
346 : ! zkm: altitude of model levels in KM
347 : ! euv_prates: array for EUV photolysis/ionization rates
348 : !==============================================================================
349 : ! Approach:
350 : ! call sphers
351 : ! input: zenith angle
352 : ! output: dsdh and nid used in slant column routine
353 : ! call slant_col
354 : ! input: dsdh and nid, etc
355 : ! output: slant column density
356 : ! calculate photon production rates and photoelectron production rates
357 : !==============================================================================
358 : ! Photolysis/ionization considered in EUV calculation
359 : !
360 : ! O + hv --> O+ (4S) + e* ! J1
361 : ! O + hv --> O+ (2D) + e* ! J2
362 : ! O + hv --> O+ (2P) + e* ! J3
363 : ! N (4S) + hv --> N+ + e* ! J4
364 : ! O2 + hv --> O2+ + e* ! J5
365 : ! N2 + hv --> N2+ + e* ! J6
366 : ! O2 + hv --> O + O+(4S) + e* ! J7
367 : ! O2 + hv --> O + O+(2D) + e* ! J8
368 : ! O2 + hv --> O + O+(2P) + e* ! J9
369 : ! N2 + hv --> N (4S) + N+ + e* ! J10
370 : ! N2 + hv --> N (2D) + N+ + e* ! J11
371 : ! O2 + hv --> O (3P) + O (3P) ! J12
372 : ! N2 + hv --> N (4S) + N (2D) ! J13
373 : !
374 : ! O + e* --> O+ (4S) + e* + e ! J14
375 : ! O + e* --> O+ (2D) + e*+ e ! J15
376 : ! O + e* --> O+ (2P) + e*+ e ! J16
377 : ! O2 + e* --> O2+ + e*+ e ! J17
378 : ! N2 + e*--> N2+ + e*+ e ! J18
379 : ! O2 + e* --> O + O+(4S) + e*+ e ! J19
380 : ! O2 + e*--> O + O+(2D) + e*+ e ! J20
381 : ! O2 + e* --> O + O+(2P) + e*+ e ! J21
382 : ! N2 + e* --> N (4S) + N+ + e*+ e ! J22
383 : ! N2 + e* --> N (2D) + N+ + e*+ e ! J23
384 : ! O2 + e* --> O (3P) + O (3P) + e* ! J24
385 : ! N2 + e* --> N (4S) + N (2D) + e* ! J25
386 : !==============================================================================
387 :
388 : use mo_jshort, only : sphers, slant_col
389 :
390 : implicit none
391 :
392 : !------------------------------------------------------------------------------
393 : ! dummy arguments
394 : !------------------------------------------------------------------------------
395 : integer, intent(in) :: nlev ! model vertical levels
396 : real(r8), intent(in) :: zen ! Zenith angle in degree
397 : real(r8), intent(in) :: occ(nlev) ! atmic oxygen number density (#/cm3)
398 : real(r8), intent(in) :: o2cc(nlev) ! Molecular oxygen number density (#/cm3)
399 : real(r8), intent(in) :: n2cc(nlev) ! molecular nitrogen number density(#/cm3)
400 : real(r8), intent(in) :: zkm(nlev) ! Altitude, km,from top to bottom
401 : real(r8), intent(out) :: euv_prates(:,:) ! EUV photolysis/ionization rates (s-1)
402 :
403 : !------------------------------------------------------------------------------
404 : ! local variables
405 : !------------------------------------------------------------------------------
406 : real(r8), parameter :: km2cm = 1.e5_r8
407 : integer :: l, k, m, n ! loop indecies
408 : real(r8) :: tau(lmax) ! wavelength dependant optical depth
409 0 : real(r8) :: delz(nlev) ! layer thickness (cm)
410 0 : real(r8) :: scol(nlev,nmaj) ! major species' (O,O2,N2) Slant Column Density
411 0 : real(r8) :: nflux(nlev,lmax)
412 : real(r8) :: wrk(nmaj) ! temporary array for photoabsorption rate
413 0 : real(r8) :: absorp(nlev,lmax) ! temporary array for photoabsorption rate
414 0 : real(r8) :: ioniz(nlev,lmax) ! temporary array for photoionization rate
415 0 : real(r8) :: dsdh(0:nlev,nlev) ! Slant path of direct beam through each layer
416 : ! crossed when travelling from the top of the atmosphere
417 : ! to layer i; dsdh(i,j), i = 0..NZ-1, j = 1..NZ-1
418 0 : integer :: nid(0:nlev) ! Number of layers crossed by the direct
419 : ! beam when travelling from the top of the
420 : ! atmosphere to layer i; NID(i), i = 0..NZ-1
421 0 : real(r8) :: p_photon(nlev,nstat,nspecies) ! photoionization/dissociation rates(s-1) (O,O2,N2,N)
422 0 : real(r8) :: p_electron(nlev,nstat,nmaj) ! photoelectron ionization/dissociation rates(s-1) (O,O2,N2)
423 :
424 0 : real(r8) :: prates(nlev,neuv)
425 : !------------------------------------------------------------------------------
426 : ! zero arrays
427 : !------------------------------------------------------------------------------
428 0 : p_photon(:,:,:) = 0._r8
429 0 : p_electron(:,:,:) = 0._r8
430 :
431 0 : call sphers( nlev, zkm, zen, dsdh, nid )
432 0 : delz(1:nlev-1) = km2cm*(zkm(1:nlev-1) - zkm(2:nlev))
433 0 : call slant_col( nlev, delz, dsdh, nid, occ, scol )
434 0 : call slant_col( nlev, delz, dsdh, nid, o2cc, scol(1,2) )
435 0 : call slant_col( nlev, delz, dsdh, nid, n2cc, scol(1,3) )
436 :
437 : !------------------------------------------------------------------------------
438 : ! The calculation is in the order from model bottom to model top
439 : ! because scol array is in this order.
440 : !------------------------------------------------------------------------------
441 0 : do k = 1,nlev
442 0 : wrk(:) = scol(k,:)
443 0 : tau(:) = matmul( sigabs(:,:nmaj),wrk )
444 0 : where( tau(:) < 20._r8 )
445 0 : nflux(k,:) = solar_etf(:) * exp( -tau(:) )
446 : elsewhere
447 : nflux(k,:) = 0._r8
448 : endwhere
449 : end do
450 :
451 : !------------------------------------------------------------------------------
452 : ! remember occ,o2cc and n2cc is from top to bottom
453 : !------------------------------------------------------------------------------
454 0 : do m = 1,nspecies
455 0 : do l = 1,lmax
456 0 : absorp(:,l) = sigabs(l,m) * nflux(:,l)
457 : end do
458 0 : if( m <= nmaj ) then
459 0 : do l = 1,lmax
460 0 : ioniz(:,l) = absorp(:,l) * branch_p(l,1,m)
461 : end do
462 0 : do n = 1,nstat
463 0 : p_photon(:,n,m) = matmul( absorp,branch_p(:,n,m) )
464 0 : p_electron(:,n,m) = matmul( ioniz,branch_e(:,n,m) )
465 : end do
466 : else
467 0 : p_photon(:,1,m) = matmul( nflux,sigabs(:,m) )
468 : end if
469 : end do
470 :
471 : !------------------------------------------------------------------------------
472 : ! set photolysis/ionization rate for each reaction
473 : !------------------------------------------------------------------------------
474 0 : prates(:,1) = p_photon(:,2,1)
475 0 : prates(:,2) = p_photon(:,3,1)
476 0 : prates(:,3) = p_photon(:,4,1)
477 0 : prates(:,4) = p_photon(:,1,4)
478 0 : prates(:,5) = p_photon(:,2,2) + p_photon(:,3,2)
479 0 : prates(:,6) = p_photon(:,2,3) + p_photon(:,3,3)
480 0 : prates(:,7) = .54_r8*p_photon(:,4,2)
481 0 : prates(:,8) = .24_r8*p_photon(:,4,2)
482 0 : prates(:,9) = .22_r8*p_photon(:,4,2)
483 0 : prates(:,10) = .2_r8*p_photon(:,4,3)
484 0 : prates(:,11) = .8_r8*p_photon(:,4,3)
485 0 : prates(:,12) = p_photon(:,5,2)
486 0 : prates(:,13) = p_photon(:,5,3)
487 0 : prates(:,14) = p_electron(:,2,1)
488 0 : prates(:,15) = p_electron(:,3,1)
489 0 : prates(:,16) = p_electron(:,4,1)
490 0 : prates(:,17) = p_electron(:,2,2) + p_electron(:,3,2)
491 0 : prates(:,18) = p_electron(:,2,3) + p_electron(:,3,3)
492 0 : prates(:,19) = .54_r8*p_electron(:,4,2)
493 0 : prates(:,20) = .24_r8*p_electron(:,4,2)
494 0 : prates(:,21) = .22_r8*p_electron(:,4,2)
495 0 : prates(:,22) = .2_r8*p_electron(:,4,3)
496 0 : prates(:,23) = .8_r8*p_electron(:,4,3)
497 0 : prates(:,24) = p_electron(:,5,2)
498 0 : prates(:,25) = p_electron(:,5,3)
499 0 : prates(:,26) = p_photon(:,1,5)
500 :
501 0 : do m = 1,neuv
502 0 : euv_prates(:,m) = prates(nlev:1:-1,m)
503 : enddo
504 :
505 0 : end subroutine jeuv
506 :
507 0 : subroutine heuv( nlev, zen, occ, o2cc, n2cc, &
508 0 : o_vmr, o2_vmr, n2_vmr, cparg, mw, &
509 0 : zkm, euv_hrates, kbot )
510 : !==============================================================================
511 : ! Purpose:
512 : ! Calculate euv photolysis heating rates
513 : !==============================================================================
514 : ! Arguments:
515 : ! nlev: number of model vertical levels
516 : ! zen: solar zenith angle in degree
517 : ! occ: atmoic oxygen number density (#/cm3)
518 : ! o2cc: molecular oxygen number density (#/cm3)
519 : ! n2cc: molecular nitrogen number density (#/cm3)
520 : ! zkm: altitude of model levels in KM
521 : ! euv_prates: array for EUV photolysis/ionization rates
522 : !==============================================================================
523 : ! Approach:
524 : ! call sphers
525 : ! input: zenith angle
526 : ! output: dsdh and nid used in slant column routine
527 : ! call slant_col
528 : ! input: dsdh and nid, etc
529 : ! output: slant column density
530 : ! calculate photon production rates and photoelectron production rates
531 : !==============================================================================
532 : ! Photolysis/ionization considered in EUV heating rate calculation
533 : ! O + hv --> O+ (4S) + e* ! J1
534 : ! O + hv --> O+ (2D) + e* ! J2
535 : ! O + hv --> O+ (2P) + e* ! J3
536 : ! N (4S) + hv --> N+ + e* ! J4
537 : ! O2 + hv --> O2+ + e* ! J5
538 : ! N2 + hv --> N2+ + e* ! J6
539 : ! O2 + hv --> O + O+(4S) + e* ! J7
540 : ! O2 + hv --> O + O+(2D) + e* ! J8
541 : ! O2 + hv --> O + O+(2P) + e* ! J9
542 : ! N2 + hv --> N (4S) + N+ + e* ! J10
543 : ! N2 + hv --> N (2D) + N+ + e* ! J11
544 : ! O2 + hv --> O (3P) + O (3P) ! J12
545 : ! N2 + hv --> N (4S) + N (2D) ! J13
546 : !==============================================================================
547 :
548 : use mo_jshort, only : sphers, slant_col
549 : use physconst, only : avogad
550 :
551 : implicit none
552 :
553 : !------------------------------------------------------------------------------
554 : ! dummy arguments
555 : !------------------------------------------------------------------------------
556 : integer, intent(in) :: nlev ! model vertical levels
557 : integer, intent(in) :: kbot ! heating vertical levels
558 : real(r8), intent(in) :: zen ! zenith angle (degrees)
559 : real(r8), intent(in) :: occ(nlev) ! atomic oxygen number density (molecules/cm^3)
560 : real(r8), intent(in) :: o2cc(nlev) ! molecular oxygen concentration (molecules/cm^3)
561 : real(r8), intent(in) :: n2cc(nlev) ! molecular nitrogen concentration (molecules/cm^3)
562 : real(r8), intent(in) :: o_vmr(nlev) ! atomic oxygen concentration (mol/mol)
563 : real(r8), intent(in) :: o2_vmr(nlev) ! molecular oxygen concentration (mol/mol)
564 : real(r8), intent(in) :: n2_vmr(nlev) ! molecular nitrogen concentration (mol/mol)
565 : real(r8), intent(in) :: zkm(nlev) ! midpoint geopotential (km)
566 : real(r8), intent(in) :: cparg(nlev) ! specific heat capacity
567 : real(r8), intent(in) :: mw(nlev) ! atm mean mass
568 : real(r8), intent(out) :: euv_hrates(:) ! euv heating rates
569 :
570 : !------------------------------------------------------------------------------
571 : ! local variables
572 : !------------------------------------------------------------------------------
573 : real(r8), parameter :: km2cm = 1.e5_r8
574 : integer :: l, k, k1, m, n ! indicies
575 : real(r8) :: tau(lmax) ! wavelength dependant optical depth
576 0 : real(r8) :: delz(nlev) ! layer thickness (cm)
577 0 : real(r8) :: hfactor(kbot)
578 0 : real(r8) :: scol(nlev,nmaj) ! major species' (O,O2,N2) Slant Column Density
579 0 : real(r8) :: nflux(kbot,lmax)
580 0 : real(r8) :: prates(kbot,13) ! working photorates array
581 : real(r8) :: wrk(nmaj) ! temporary array for photoabsorption rate
582 0 : real(r8) :: absorp(kbot,lmax) ! temporary array for photoabsorption rate
583 0 : real(r8) :: dsdh(0:nlev,nlev) ! Slant path of direct beam through each layer
584 : ! crossed when travelling from the top of the atmosphere
585 : ! to layer i; dsdh(i,j), i = 0..NZ-1, j = 1..NZ-1
586 0 : integer :: nid(0:nlev) ! Number of layers crossed by the direct
587 : ! beam when travelling from the top of the
588 : ! atmosphere to layer i; NID(i), i = 0..NZ-1
589 0 : real(r8) :: p_photon(kbot,nstat,nspecies) ! photoionization/dissociation rates(s-1) (O,O2,N2,N)
590 :
591 : !------------------------------------------------------------------------------
592 : ! zero arrays
593 : !------------------------------------------------------------------------------
594 0 : p_photon(:,:,:) = 0._r8
595 :
596 0 : call sphers( nlev, zkm, zen, dsdh, nid )
597 0 : delz(1:nlev-1) = km2cm*(zkm(1:nlev-1) - zkm(2:nlev))
598 0 : call slant_col( nlev, delz, dsdh, nid, occ, scol )
599 0 : call slant_col( nlev, delz, dsdh, nid, o2cc, scol(1,2) )
600 0 : call slant_col( nlev, delz, dsdh, nid, n2cc, scol(1,3) )
601 :
602 : !------------------------------------------------------------------------------
603 : ! The calculation is in the order from model bottom to model top
604 : ! because scol array is in this order.
605 : !------------------------------------------------------------------------------
606 0 : do k = 1,kbot
607 0 : k1 = nlev - kbot + k
608 0 : wrk(:) = scol(k1,:)
609 0 : tau(:) = matmul( sigabs(:,:nmaj),wrk )
610 0 : where( tau(:) < 20._r8 )
611 0 : nflux(k,:) = energy(:) * solar_etf(:) * exp( -tau(:) )
612 : elsewhere
613 0 : nflux(k,:) = 0._r8
614 : endwhere
615 : end do
616 :
617 : !------------------------------------------------------------------------------
618 : ! remember occ,o2cc and n2cc is from top to bottom
619 : !------------------------------------------------------------------------------
620 0 : do m = 1,nspecies
621 0 : do l = 1,lmax
622 0 : absorp(:,l) = sigabs(l,m) * nflux(:,l)
623 : end do
624 0 : if( m <= nmaj ) then
625 0 : do n = 1,nstat
626 0 : p_photon(:,n,m) = matmul( absorp,branch_p(:,n,m) )
627 : end do
628 : else
629 0 : p_photon(:,1,m) = matmul( nflux,sigabs(:,m) )
630 : end if
631 : end do
632 :
633 : !------------------------------------------------------------------------------
634 : ! use photolysis rates to compute heating rates
635 : ! only EUV photolysis reactions in the mechanism contribute to heating
636 : !------------------------------------------------------------------------------
637 0 : prates = 0._r8
638 :
639 0 : if (do_heating(1)) prates(:,1) = p_photon(:,2,1)
640 0 : if (do_heating(2)) prates(:,2) = p_photon(:,3,1)
641 0 : if (do_heating(3)) prates(:,3) = p_photon(:,4,1)
642 :
643 0 : if (do_heating(5)) prates(:,5) = p_photon(:,2,2) + p_photon(:,3,2)
644 0 : if (do_heating(6)) prates(:,6) = p_photon(:,2,3) + p_photon(:,3,3)
645 0 : if (do_heating(7)) prates(:,7) = .54_r8*p_photon(:,4,2)
646 0 : if (do_heating(8)) prates(:,8) = .24_r8*p_photon(:,4,2)
647 0 : if (do_heating(9)) prates(:,9) = .22_r8*p_photon(:,4,2)
648 0 : if (do_heating(10)) prates(:,10) = .2_r8*p_photon(:,4,3)
649 0 : if (do_heating(11)) prates(:,11) = .8_r8*p_photon(:,4,3)
650 0 : if (do_heating(12)) prates(:,12) = p_photon(:,5,2)
651 0 : if (do_heating(13)) prates(:,13) = p_photon(:,5,3)
652 0 : hfactor(:) = avogad/(cparg(:kbot)*mw(:kbot))
653 0 : euv_hrates(kbot+1:nlev) = 0._r8
654 0 : euv_hrates(:kbot) = ((prates(kbot:1:-1,1) + prates(kbot:1:-1,2) + prates(kbot:1:-1,3))*o_vmr(:kbot) &
655 : + (prates(kbot:1:-1,5) + prates(kbot:1:-1,7) + prates(kbot:1:-1,8) &
656 : + prates(kbot:1:-1,9) + prates(kbot:1:-1,12))*o2_vmr(:kbot) &
657 : + (prates(kbot:1:-1,6) + prates(kbot:1:-1,10) &
658 0 : + prates(kbot:1:-1,11) + prates(kbot:1:-1,13))*n2_vmr(:kbot))*hfactor(:)
659 :
660 0 : end subroutine heuv
661 :
662 : end module mo_jeuv
|
