Line data Source code
1 : !-------------------------------------------------------------------------------
2 : ! Energetic Particle Precipitation (EPP) forcings module
3 : ! Manages ionization of the atmosphere due to energetic particles, which consists of
4 : ! solar protons events (SPE), galactic cosmic rays(GCR), medium energy electrons (MEE)
5 : !-------------------------------------------------------------------------------
6 : module epp_ionization
7 : use shr_kind_mod, only : r8 => shr_kind_r8, cs => shr_kind_cs, cl=> shr_kind_cl
8 : use spmd_utils, only : masterproc
9 : use cam_abortutils, only : endrun
10 : use cam_logfile, only : iulog
11 : use ppgrid, only : pcols, pver, begchunk, endchunk
12 : use phys_grid, only : get_ncols_p
13 : use pio, only : var_desc_t, file_desc_t
14 : use pio, only : pio_get_var, pio_inq_varid, pio_get_att
15 : use pio, only : pio_inq_varndims, pio_inq_vardimid, pio_inq_dimname, pio_inq_dimlen
16 : use pio, only : PIO_NOWRITE
17 : use cam_pio_utils, only : cam_pio_openfile
18 : use ioFileMod, only : getfil
19 : use input_data_utils, only : time_coordinate
20 :
21 : implicit none
22 : private
23 :
24 : public :: epp_ionization_readnl ! read namelist variables
25 : public :: epp_ionization_init ! initialization
26 : public :: epp_ionization_adv ! read and time/space interpolate the data
27 : public :: epp_ionization_ionpairs! ion pairs production rates
28 : public :: epp_ionization_setmag ! update geomagnetic coordinates mapping
29 : public :: epp_ionization_active
30 :
31 : character(len=cl) :: epp_all_filepath = 'NONE'
32 : character(len=cs) :: epp_all_varname = 'epp_ion_rates'
33 : character(len=cl) :: epp_mee_filepath = 'NONE'
34 : character(len=cs) :: epp_mee_varname = 'iprm'
35 : character(len=cl) :: epp_spe_filepath = 'NONE'
36 : character(len=cs) :: epp_spe_varname = 'iprp'
37 : character(len=cl) :: epp_gcr_filepath = 'NONE'
38 : character(len=cs) :: epp_gcr_varname = 'iprg'
39 :
40 : logical, protected :: epp_ionization_active = .false.
41 :
42 : type input_obj_t
43 : type(file_desc_t) :: fid
44 : type(var_desc_t) :: vid
45 : character(len=32) :: units
46 : integer :: nlevs = 0
47 : integer :: nglats = 0
48 : real(r8), allocatable :: press(:)
49 : real(r8), allocatable :: glats(:)
50 : real(r8), allocatable :: gwght(:,:) ! (pcol, begchunk:endchunk)
51 : integer, allocatable :: glatn(:,:) ! (pcol, begchunk:endchunk)
52 : real(r8), allocatable :: indata(:,:,:,:) ! (pcol,nlevs,begchunk:endchunk,2) inputs at indexm and indexp
53 : type(time_coordinate) :: time_coord
54 : endtype input_obj_t
55 :
56 : type(input_obj_t), pointer :: epp_in => null()
57 : type(input_obj_t), pointer :: spe_in => null()
58 : type(input_obj_t), pointer :: mee_in => null()
59 : type(input_obj_t), pointer :: gcr_in => null()
60 :
61 : contains
62 :
63 : !-----------------------------------------------------------------------------
64 : !-----------------------------------------------------------------------------
65 1536 : subroutine epp_ionization_readnl(nlfile)
66 :
67 : use namelist_utils, only: find_group_name
68 : use units, only: getunit, freeunit
69 : use spmd_utils, only: mpicom, mpi_character, masterprocid
70 :
71 : character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input
72 :
73 : ! Local variables
74 : integer :: unitn, ierr
75 : character(len=*), parameter :: subname = 'epp_ionization_readnl'
76 :
77 : namelist /epp_ionization_nl/ epp_all_filepath, epp_all_varname, &
78 : epp_mee_filepath, epp_mee_varname, epp_spe_filepath, epp_spe_varname, epp_gcr_filepath, epp_gcr_varname
79 :
80 : ! Read namelist
81 1536 : if (masterproc) then
82 2 : unitn = getunit()
83 2 : open( unitn, file=trim(nlfile), status='old' )
84 2 : call find_group_name(unitn, 'epp_ionization_nl', status=ierr)
85 2 : if (ierr == 0) then
86 0 : read(unitn, epp_ionization_nl, iostat=ierr)
87 0 : if (ierr /= 0) then
88 0 : call endrun(subname // ':: ERROR reading namelist')
89 : end if
90 : end if
91 2 : close(unitn)
92 2 : call freeunit(unitn)
93 : end if
94 :
95 : ! Broadcast namelist variables
96 1536 : call mpi_bcast(epp_all_filepath, len(epp_all_filepath), mpi_character, masterprocid, mpicom, ierr)
97 1536 : call mpi_bcast(epp_mee_filepath, len(epp_mee_filepath), mpi_character, masterprocid, mpicom, ierr)
98 1536 : call mpi_bcast(epp_spe_filepath, len(epp_spe_filepath), mpi_character, masterprocid, mpicom, ierr)
99 1536 : call mpi_bcast(epp_gcr_filepath, len(epp_gcr_filepath), mpi_character, masterprocid, mpicom, ierr)
100 :
101 1536 : call mpi_bcast(epp_all_varname, len(epp_all_varname), mpi_character, masterprocid, mpicom, ierr)
102 1536 : call mpi_bcast(epp_mee_varname, len(epp_mee_varname), mpi_character, masterprocid, mpicom, ierr)
103 1536 : call mpi_bcast(epp_spe_varname, len(epp_spe_varname), mpi_character, masterprocid, mpicom, ierr)
104 1536 : call mpi_bcast(epp_gcr_varname, len(epp_gcr_varname), mpi_character, masterprocid, mpicom, ierr)
105 :
106 1536 : epp_ionization_active = epp_all_filepath /= 'NONE'
107 1536 : epp_ionization_active = epp_mee_filepath /= 'NONE' .or. epp_ionization_active
108 1536 : epp_ionization_active = epp_spe_filepath /= 'NONE' .or. epp_ionization_active
109 1536 : epp_ionization_active = epp_gcr_filepath /= 'NONE' .or. epp_ionization_active
110 :
111 1536 : if ( epp_ionization_active .and. masterproc ) then
112 0 : write(iulog,*) subname//':: epp_all_filepath = '//trim(epp_all_filepath)
113 0 : write(iulog,*) subname//':: epp_mee_filepath = '//trim(epp_mee_filepath)
114 0 : write(iulog,*) subname//':: epp_spe_filepath = '//trim(epp_spe_filepath)
115 0 : write(iulog,*) subname//':: epp_gcr_filepath = '//trim(epp_gcr_filepath)
116 : endif
117 :
118 1536 : end subroutine epp_ionization_readnl
119 :
120 : !-----------------------------------------------------------------------------
121 : !-----------------------------------------------------------------------------
122 0 : subroutine epp_ionization_init()
123 : use cam_history, only : addfld
124 :
125 : character(len=32) :: fldunits
126 0 : fldunits = ''
127 :
128 0 : if (epp_all_filepath /= 'NONE') then
129 0 : epp_in => create_input_obj(epp_all_filepath,epp_all_varname)
130 0 : fldunits = trim(epp_in%units)
131 : else
132 0 : if (epp_mee_filepath /= 'NONE') then
133 0 : mee_in => create_input_obj(epp_mee_filepath,epp_mee_varname)
134 0 : fldunits = trim(mee_in%units)
135 : endif
136 0 : if (epp_spe_filepath /= 'NONE') then
137 0 : spe_in => create_input_obj(epp_spe_filepath,epp_spe_varname)
138 0 : fldunits = trim(spe_in%units)
139 : endif
140 0 : if (epp_gcr_filepath /= 'NONE') then
141 0 : gcr_in => create_input_obj(epp_gcr_filepath,epp_gcr_varname)
142 0 : fldunits = trim(gcr_in%units)
143 : endif
144 : endif
145 0 : call addfld( 'EPPions', (/ 'lev' /), 'A', fldunits, 'EPP ionization data' )
146 :
147 0 : end subroutine epp_ionization_init
148 :
149 : !-----------------------------------------------------------------------------
150 : !-----------------------------------------------------------------------------
151 0 : subroutine epp_ionization_setmag( maglat )
152 : real(r8), intent(in) :: maglat(pcols,begchunk:endchunk)
153 :
154 0 : if (.not.epp_ionization_active) return
155 :
156 0 : if ( associated(epp_in) ) then
157 0 : call set_wghts(maglat,epp_in)
158 : else
159 0 : if ( associated(mee_in) ) then
160 0 : call set_wghts(maglat,mee_in)
161 : endif
162 0 : if ( associated(spe_in) ) then
163 0 : call set_wghts(maglat,spe_in)
164 : endif
165 0 : if ( associated(gcr_in) ) then
166 0 : call set_wghts(maglat,gcr_in)
167 : endif
168 : endif
169 :
170 0 : end subroutine epp_ionization_setmag
171 :
172 : !-----------------------------------------------------------------------------
173 : !-----------------------------------------------------------------------------
174 370944 : subroutine epp_ionization_adv
175 :
176 370944 : if (.not.epp_ionization_active) return
177 :
178 0 : if ( associated(epp_in) ) then
179 0 : call update_input(epp_in)
180 : else
181 0 : if ( associated(spe_in) ) then
182 0 : call update_input(spe_in)
183 : endif
184 0 : if ( associated(gcr_in) ) then
185 0 : call update_input(gcr_in)
186 : endif
187 0 : if ( associated(mee_in) ) then
188 0 : call update_input(mee_in)
189 : endif
190 : endif
191 :
192 : end subroutine epp_ionization_adv
193 :
194 : !-----------------------------------------------------------------------------
195 : !-----------------------------------------------------------------------------
196 1489176 : subroutine epp_ionization_ionpairs( ncol, lchnk, pmid, temp, ionpairs )
197 :
198 : integer, intent(in) :: ncol, lchnk
199 : real(r8), intent(in) :: pmid(:,:), temp(:,:)
200 : real(r8), intent(out) :: ionpairs(:,:) ! ion pair production rate
201 :
202 2314006344 : ionpairs = 0._r8
203 1489176 : if (.not.epp_ionization_active) return
204 :
205 0 : if ( associated(epp_in) ) then
206 0 : ionpairs(:ncol,:) = ionpairs(:ncol,:) + interp_ionpairs( ncol, lchnk, pmid, temp, epp_in )
207 : else
208 0 : if ( associated(spe_in) ) then
209 0 : ionpairs(:ncol,:) = ionpairs(:ncol,:) + interp_ionpairs( ncol, lchnk, pmid, temp, spe_in )
210 : endif
211 0 : if ( associated(gcr_in) ) then
212 0 : ionpairs(:ncol,:) = ionpairs(:ncol,:) + interp_ionpairs( ncol, lchnk, pmid, temp, gcr_in )
213 : endif
214 0 : if ( associated(mee_in) ) then
215 0 : ionpairs(:ncol,:) = ionpairs(:ncol,:) + interp_ionpairs( ncol, lchnk, pmid, temp, mee_in )
216 : endif
217 : endif
218 :
219 : end subroutine epp_ionization_ionpairs
220 :
221 : ! private methods
222 : !-----------------------------------------------------------------------------
223 : !-----------------------------------------------------------------------------
224 0 : subroutine update_input( input )
225 : type(input_obj_t), pointer :: input
226 :
227 0 : if ( input%time_coord%read_more() ) then
228 0 : call input%time_coord%advance()
229 0 : call read_next_data( input )
230 : else
231 0 : call input%time_coord%advance()
232 : endif
233 :
234 0 : end subroutine update_input
235 :
236 : !-----------------------------------------------------------------------------
237 : !-----------------------------------------------------------------------------
238 0 : subroutine read_next_data( input )
239 : type(input_obj_t), pointer :: input
240 :
241 : ! read data corresponding surrounding time indices
242 0 : if ( input%nglats > 0 ) then
243 0 : call read_2d_profile( input )
244 : else
245 0 : call read_1d_profile( input )
246 : endif
247 :
248 0 : end subroutine read_next_data
249 :
250 : !-----------------------------------------------------------------------------
251 : !-----------------------------------------------------------------------------
252 0 : function interp_ionpairs( ncol, lchnk, pmid, temp, input ) result( ionpairs )
253 : use interpolate_data, only : lininterp, lininterp_init, lininterp_finish, extrap_method_zero, interp_type
254 : use air_composition, only : rairv
255 : use cam_history, only : outfld
256 :
257 : integer, intent(in) :: ncol, lchnk
258 : real(r8), intent(in) :: pmid(:,:) ! Pa
259 : real(r8), intent(in) :: temp(:,:) ! K
260 : type(input_obj_t), pointer :: input
261 : real(r8) :: ionpairs(ncol,pver)
262 :
263 : real(r8) :: fctr1, fctr2
264 0 : real(r8) :: wrk(ncol,input%nlevs)
265 0 : real(r8) :: ions_diags(ncol,pver) ! for diagnostics
266 : integer :: i
267 : type (interp_type) :: interp_wgts
268 :
269 0 : if (input%time_coord%time_interp) then
270 : ! time interpolate
271 0 : fctr1 = input%time_coord%wghts(1)
272 0 : fctr2 = input%time_coord%wghts(2)
273 0 : wrk(:ncol,:) = fctr1*input%indata(:ncol,:,lchnk,1) + fctr2*input%indata(:ncol,:,lchnk,2)
274 : else
275 0 : wrk(:ncol,:) = input%indata(:ncol,:,lchnk,1)
276 : endif
277 :
278 : ! vertical interpolate ...
279 : ! interpolate to model levels
280 0 : do i = 1,ncol
281 :
282 : ! interpolate over log pressure
283 0 : call lininterp_init( log(input%press(:input%nlevs)*1.e2_r8), input%nlevs, &
284 0 : log(pmid(i,:pver)), pver, extrap_method_zero, interp_wgts )
285 0 : call lininterp( wrk(i,:input%nlevs), input%nlevs, &
286 0 : ionpairs(i,:pver), pver, interp_wgts )
287 0 : call lininterp_finish(interp_wgts)
288 :
289 0 : ions_diags(i,:pver) = ionpairs(i,:pver)
290 :
291 0 : if ( index(trim(input%units), 'g^-1') > 0 ) then
292 : ! convert to ionpairs/cm3/sec
293 0 : ionpairs(i,:pver) = ionpairs(i,:pver) *(1.e-3_r8*pmid(i,:pver)/(rairv(i,:pver,lchnk)*temp(i,:pver)))
294 : endif
295 : enddo
296 :
297 0 : call outfld( 'EPPions', ions_diags(:ncol,:), ncol, lchnk )
298 :
299 0 : end function interp_ionpairs
300 :
301 : !-----------------------------------------------------------------------------
302 : ! read 2D profile (geomag-lat vs press) and transfer to geographic grid
303 : !-----------------------------------------------------------------------------
304 0 : subroutine read_2d_profile( input )
305 :
306 : type(input_obj_t), pointer :: input
307 :
308 : ! local vars
309 0 : real(r8) :: wrk2d( input%nglats, input%nlevs, 2 )
310 : integer :: t, c, i, ntimes, ncols, ierr
311 : real(r8) :: wght1, wght2
312 : integer :: gndx1, gndx2
313 : integer :: cnt(3), strt(3)
314 :
315 0 : if (input%time_coord%time_interp) then
316 : ntimes = 2
317 : else
318 0 : ntimes = 1
319 : endif
320 :
321 0 : cnt(1) = input%nglats
322 0 : cnt(2) = input%nlevs
323 0 : cnt(3) = ntimes
324 :
325 0 : strt(:) = 1
326 0 : strt(3) = input%time_coord%indxs(1)
327 :
328 0 : ierr = pio_get_var( input%fid, input%vid, strt, cnt, wrk2d )
329 :
330 0 : do t = 1,ntimes
331 0 : do c=begchunk,endchunk
332 0 : ncols = get_ncols_p(c)
333 0 : do i = 1,ncols
334 0 : gndx1 = input%glatn(i,c)
335 0 : if (gndx1>0) then
336 0 : wght1 = input%gwght(i,c)
337 0 : gndx2 = gndx1+1
338 0 : if (gndx2.le.input%nglats) then
339 0 : wght2 = 1._r8-wght1
340 0 : input%indata(i,:,c,t) = wght1*wrk2d(gndx1,:,t) &
341 0 : + wght2*wrk2d(gndx2,:,t)
342 : else
343 0 : input%indata(i,:,c,t) = wght1*wrk2d(gndx1,:,t)
344 : endif
345 : else
346 0 : input%indata(i,:,c,t) = 0._r8
347 : endif
348 : end do
349 : end do
350 : end do
351 :
352 0 : end subroutine read_2d_profile
353 :
354 : !-----------------------------------------------------------------------------
355 : ! read 1D vertical profile and transfer to geographic grid poleward of 60 degrees geomag-lat
356 : !-----------------------------------------------------------------------------
357 0 : subroutine read_1d_profile( input )
358 :
359 : type(input_obj_t), pointer :: input
360 :
361 : ! local vars
362 0 : real(r8) :: wrk( input%nlevs, 2 )
363 : integer :: t, c, i, ntimes, ncols, ierr
364 : integer :: cnt(2), strt(2)
365 :
366 0 : if (input%time_coord%time_interp) then
367 : ntimes = 2
368 : else
369 0 : ntimes = 1
370 : endif
371 :
372 0 : cnt(1) = input%nlevs
373 0 : cnt(2) = ntimes
374 :
375 0 : strt(:) = 1
376 0 : strt(2) = input%time_coord%indxs(1)
377 :
378 0 : ierr = pio_get_var( input%fid, input%vid, strt, cnt, wrk )
379 :
380 0 : do t = 1,ntimes
381 0 : do c=begchunk,endchunk
382 0 : ncols = get_ncols_p(c)
383 0 : do i = 1,ncols
384 0 : input%indata(i,:,c,t) = input%gwght(i,c)*wrk(:,t)
385 : end do
386 : end do
387 : end do
388 :
389 0 : end subroutine read_1d_profile
390 :
391 : !-----------------------------------------------------------------------------
392 : !-----------------------------------------------------------------------------
393 0 : function create_input_obj( path, varname ) result(in_obj)
394 : use infnan, only : nan, assignment(=)
395 :
396 : character(*), intent(in) :: path
397 : character(*), intent(in) :: varname
398 : type(input_obj_t), pointer :: in_obj
399 :
400 : character(len=cl) :: filen
401 : character(len=cl) :: data_units
402 : character(len=cs) :: dimname
403 : integer :: i, ierr
404 0 : integer, allocatable :: dimids(:)
405 : integer :: pres_did, pres_vid, glat_did, glat_vid, ndims
406 :
407 0 : if (path .eq. 'NONE') return
408 :
409 0 : allocate(in_obj)
410 :
411 0 : call in_obj%time_coord%initialize( path )
412 :
413 0 : call getfil( path, filen, 0 )
414 0 : call cam_pio_openfile( in_obj%fid, filen, PIO_NOWRITE )
415 :
416 0 : ierr = pio_inq_varid( in_obj%fid, varname, in_obj%vid )
417 :
418 0 : ierr = pio_get_att( in_obj%fid, in_obj%vid, 'units', data_units)
419 0 : in_obj%units = trim(data_units(1:32))
420 :
421 0 : ierr = pio_inq_varndims( in_obj%fid, in_obj%vid, ndims )
422 0 : allocate( dimids(ndims) )
423 :
424 0 : ierr = pio_inq_vardimid( in_obj%fid, in_obj%vid, dimids)
425 0 : pres_did = -1
426 0 : glat_did = -1
427 0 : do i = 1,ndims
428 0 : ierr = pio_inq_dimname( in_obj%fid, dimids(i), dimname )
429 0 : select case( trim(dimname(1:4)) )
430 : case ( 'pres', 'lev', 'plev' )
431 0 : pres_did = dimids(i)
432 0 : ierr = pio_inq_varid( in_obj%fid, dimname, pres_vid)
433 : case ( 'glat' )
434 0 : glat_did = dimids(i)
435 0 : ierr = pio_inq_varid( in_obj%fid, dimname, glat_vid)
436 : case default
437 : end select
438 : end do
439 :
440 0 : deallocate( dimids )
441 :
442 0 : if (pres_did>0) then
443 0 : ierr = pio_inq_dimlen( in_obj%fid, pres_did, in_obj%nlevs )
444 0 : allocate( in_obj%press(in_obj%nlevs) )
445 0 : ierr = pio_get_var( in_obj%fid, pres_vid, in_obj%press )
446 : endif
447 0 : if (glat_did>0) then
448 0 : ierr = pio_inq_dimlen( in_obj%fid, glat_did, in_obj%nglats )
449 0 : allocate( in_obj%glats(in_obj%nglats) )
450 0 : ierr = pio_get_var( in_obj%fid, glat_vid, in_obj%glats )
451 0 : allocate( in_obj%glatn(pcols,begchunk:endchunk) )
452 : endif
453 :
454 0 : allocate( in_obj%gwght(pcols,begchunk:endchunk) )
455 :
456 0 : if (in_obj%time_coord%time_interp) then
457 0 : allocate( in_obj%indata(pcols,in_obj%nlevs,begchunk:endchunk,2) )
458 : else
459 0 : allocate( in_obj%indata(pcols,in_obj%nlevs,begchunk:endchunk,1) )
460 : endif
461 0 : in_obj%indata = nan
462 :
463 0 : end function create_input_obj
464 :
465 : !-----------------------------------------------------------------------
466 : !-----------------------------------------------------------------------
467 0 : subroutine set_wghts( maglat, input )
468 :
469 : real(r8), intent(in) :: maglat(pcols,begchunk:endchunk)
470 : type(input_obj_t), pointer :: input
471 :
472 : integer :: i, c, ncols, imag
473 :
474 0 : if (input%nglats>1) then ! read in general EPP 2D ionpairs production rates
475 0 : do c = begchunk,endchunk
476 0 : ncols = get_ncols_p(c)
477 0 : col_loop: do i = 1,ncols
478 0 : if ( maglat(i,c) .lt. input%glats(1) ) then
479 0 : input%glatn(i,c) = 1
480 0 : input%gwght(i,c) = 1._r8
481 0 : elseif ( maglat(i,c) .gt. input%glats(input%nglats) ) then
482 0 : input%glatn(i,c) = input%nglats
483 0 : input%gwght(i,c) = 1._r8
484 : else
485 0 : mag_loop: do imag = 1,input%nglats-1
486 0 : if ( maglat(i,c) .ge. input%glats(imag) .and. &
487 0 : maglat(i,c) .lt. input%glats(imag+1) ) then
488 0 : input%gwght(i,c) = (input%glats(imag+1)-maglat(i,c) ) &
489 0 : / (input%glats(imag+1)-input%glats(imag))
490 0 : input%glatn(i,c) = imag
491 0 : exit mag_loop
492 : endif
493 : enddo mag_loop
494 : endif
495 : enddo col_loop
496 : enddo
497 : else ! read in 1D SPE ionpairs profile ...
498 0 : do c = begchunk,endchunk
499 0 : ncols = get_ncols_p(c)
500 0 : do i = 1,ncols
501 0 : if ( abs(maglat(i,c)) .ge. 60._r8 ) then ! poleward of 60 degrees
502 0 : input%gwght(i,c) = 1._r8
503 : else
504 0 : input%gwght(i,c) = 0._r8
505 : endif
506 : enddo
507 : enddo
508 : endif
509 :
510 0 : call read_next_data( input ) ! update the inputs when wghts are updated
511 :
512 0 : end subroutine set_wghts
513 :
514 0 : end module epp_ionization
|