Line data Source code
1 : module mo_setext
2 :
3 : use shr_kind_mod, only : r8 => shr_kind_r8
4 : use shr_const_mod,only : pi => shr_const_pi
5 : use cam_logfile, only : iulog
6 :
7 : private
8 : public :: setext_inti, setext, has_ions
9 :
10 : save
11 :
12 : integer :: co_ndx, no_ndx, xno_ndx, o_ndx
13 : integer :: op_ndx, o2p_ndx, np_ndx, n2p_ndx, n2d_ndx, n_ndx, e_ndx, oh_ndx
14 : logical :: has_ions = .false.
15 : logical :: has_dregion_ions = .false.
16 : integer :: aoa_nh_ndx
17 :
18 : contains
19 :
20 1490712 : subroutine setext_inti
21 : !--------------------------------------------------------
22 : ! ... Initialize the external forcing module
23 : !--------------------------------------------------------
24 :
25 : use mo_chem_utls, only : get_extfrc_ndx, get_spc_ndx
26 : use cam_history, only : addfld
27 : use spmd_utils, only : masterproc
28 : use mee_ionization,only : mee_ion_init
29 :
30 : implicit none
31 :
32 1536 : co_ndx = get_extfrc_ndx( 'CO' )
33 1536 : no_ndx = get_extfrc_ndx( 'NO' )
34 1536 : xno_ndx = get_extfrc_ndx( 'XNO' )
35 1536 : aoa_nh_ndx = get_extfrc_ndx('AOA_NH')
36 :
37 1536 : op_ndx = get_extfrc_ndx( 'Op' )
38 1536 : o2p_ndx = get_extfrc_ndx( 'O2p' )
39 1536 : np_ndx = get_extfrc_ndx( 'Np' )
40 1536 : n2p_ndx = get_extfrc_ndx( 'N2p' )
41 1536 : n2d_ndx = get_extfrc_ndx( 'N2D' )
42 1536 : n_ndx = get_extfrc_ndx( 'N' )
43 1536 : e_ndx = get_extfrc_ndx( 'e' )
44 1536 : oh_ndx = get_extfrc_ndx( 'OH' )
45 1536 : o_ndx = get_extfrc_ndx( 'O' )
46 :
47 1536 : has_ions = op_ndx > 0 .and. o2p_ndx > 0 .and. np_ndx > 0 .and. n2p_ndx > 0 .and. e_ndx > 0
48 1536 : has_dregion_ions = has_ions .and. (get_spc_ndx( 'O2m' )>0)
49 :
50 1536 : if (masterproc) then
51 2 : write(iulog,*) ' '
52 2 : write(iulog,*) 'setext_inti: diagnostics: co_ndx, no_ndx, xno_ndx'
53 2 : write(iulog,'(10i5)') co_ndx, no_ndx, xno_ndx
54 : endif
55 :
56 3072 : call addfld( 'NO_Lightning', (/ 'lev' /), 'A','molec/cm3/s', 'lightning NO source' )
57 3072 : call addfld( 'NO_Aircraft', (/ 'lev' /), 'A', 'molec/cm3/s', 'aircraft NO source' )
58 3072 : call addfld( 'CO_Aircraft', (/ 'lev' /), 'A', 'molec/cm3/s', 'aircraft CO source' )
59 :
60 3072 : call addfld( 'EPP_ionpairs', (/ 'lev' /), 'A', 'pairs/cm3/s', 'EPP ionization forcing' )
61 3072 : call addfld( 'GCR_ionpairs', (/ 'lev' /), 'A', 'pairs/cm3/s', 'GCR ionization forcing' )
62 :
63 1536 : if (.not.has_dregion_ions) then
64 1536 : if ( n2d_ndx > 0 .and. n_ndx>0 ) then
65 0 : call addfld( 'N4S_EPP', (/ 'lev' /), 'I', 'molec/cm3/s', 'solar proton event N(4S) source' )
66 0 : call addfld( 'N2D_EPP', (/ 'lev' /), 'I', 'molec/cm3/s', 'solar proton event N(2D) source' )
67 1536 : elseif ( no_ndx > 0 .and. n_ndx>0 ) then
68 0 : call addfld( 'N4S_EPP', (/ 'lev' /), 'I', 'molec/cm3/s', 'solar proton event N(4S) source' )
69 0 : call addfld( 'NO_EPP', (/ 'lev' /), 'I', 'molec/cm3/s', 'solar proton event NO source' )
70 : endif
71 1536 : if ( oh_ndx > 0 ) then
72 0 : call addfld( 'OH_EPP', (/ 'lev' /), 'I', 'molec/cm3/s', 'solar proton event OH source' )
73 : endif
74 : endif
75 1536 : if ( has_ions ) then
76 0 : call addfld( 'P_Op', (/ 'lev' /), 'I', '/s', 'production o+' )
77 0 : call addfld( 'P_O2p', (/ 'lev' /), 'I', '/s', 'production o2+' )
78 0 : call addfld( 'P_N2p', (/ 'lev' /), 'I', '/s', 'production n2+' )
79 0 : call addfld( 'P_Np', (/ 'lev' /), 'I', '/s', 'production n+' )
80 0 : call addfld( 'P_IONS', (/ 'lev' /), 'I', '/s', 'total ion production' )
81 : endif
82 :
83 1536 : if ( aoa_nh_ndx > 0 ) then
84 0 : call addfld('AOA_NH_XFRC', (/ 'lev' /), 'A', 'molec/cm3/s', 'external forcing for AOA_NH' )
85 : endif
86 :
87 1536 : call mee_ion_init()
88 :
89 1536 : end subroutine setext_inti
90 :
91 1489176 : subroutine setext( extfrc, zint_abs, zint_rel, cldtop, &
92 1489176 : zmid, lchnk, tfld, o2mmr, ommr, &
93 1489176 : pmid, mbar, rlats, calday, ncol, rlons, pbuf )
94 : !--------------------------------------------------------
95 : ! ... for this latitude slice:
96 : ! - form the production from datasets
97 : ! - form the nox (xnox) production from lighing
98 : ! - form the nox (xnox) production from airplanes
99 : ! - form the co production from airplanes
100 : !--------------------------------------------------------
101 :
102 1536 : use cam_history, only : outfld
103 : use ppgrid, only : pver, pcols
104 : use mo_airplane, only : airpl_set
105 : use chem_mods, only : extcnt
106 : use mo_lightning, only : prod_no
107 :
108 : use mo_extfrc, only : extfrc_set
109 : use hco_cc_emissions, only : hco_set_extfrc
110 : use chem_mods, only : extcnt
111 : use tracer_srcs, only : num_tracer_srcs, tracer_src_flds, get_srcs_data
112 : use mo_chem_utls, only : get_extfrc_ndx
113 :
114 : use mo_aurora, only : aurora
115 : use gcr_ionization, only : gcr_ionization_ionpairs
116 : use epp_ionization, only : epp_ionization_ionpairs
117 : use mee_ionization, only : mee_ionpairs
118 : use spehox, only : hox_prod_factor
119 :
120 : use physics_buffer, only : physics_buffer_desc
121 :
122 : use phys_control, only : use_hemco ! Use Harmonized Emissions Component (HEMCO)
123 :
124 : implicit none
125 :
126 : !--------------------------------------------------------
127 : ! ... dummy arguments
128 : !--------------------------------------------------------
129 : !--------------------------------------------------------
130 : integer, intent(in) :: lchnk ! chunk id
131 : integer, intent(in) :: ncol ! columns in chunk
132 : real(r8), intent(in) :: zint_abs(ncol,pver+1) ! interface geopot height ( km )
133 : real(r8), intent(in) :: zint_rel(ncol,pver+1) ! interface geopot height ( km )
134 : real(r8), intent(in) :: cldtop(ncol) ! cloud top index
135 : real(r8), intent(out) :: extfrc(ncol,pver,extcnt) ! the "extraneous" forcing
136 :
137 : real(r8), intent(in) :: calday ! calendar day of year
138 : real(r8), intent(in) :: rlats(ncol) ! column latitudes (radians)
139 : real(r8), intent(in) :: rlons(ncol) ! column longitudes (radians)
140 : real(r8), intent(in) :: zmid(ncol,pver) ! midpoint geopot height ( km )
141 : real(r8), intent(in) :: pmid(pcols,pver) ! midpoint pressure (Pa)
142 : real(r8), intent(in) :: tfld(pcols,pver) ! midpoint temperature (K)
143 : real(r8), intent(in) :: mbar(ncol,pver) ! mean molecular mass (g/mole)
144 : real(r8), intent(in) :: o2mmr(ncol,pver) ! o2 concentration (kg/kg)
145 : real(r8), intent(in) :: ommr(ncol,pver) ! o concentration (kg/kg)
146 :
147 : type(physics_buffer_desc),pointer :: pbuf(:)
148 :
149 : !--------------------------------------------------------
150 : ! ... local variables
151 : !--------------------------------------------------------
152 : integer :: i, k, cldind
153 2978352 : real(r8) :: srcs_offline( ncol, pver )
154 : integer :: ndx
155 :
156 2978352 : real(r8), dimension(ncol,pver) :: no_lgt
157 :
158 2978352 : real(r8) :: gcr_ipr(ncol,pver) ! ion pairs production rate associated with galactic comsic rays
159 2978352 : real(r8) :: epp_ipr(ncol,pver) ! ion pairs production rate associated with energetic particles
160 2978352 : real(r8) :: epp_hox(ncol,pver) ! HOx production rate associated with energetic particles
161 :
162 : real(r8), parameter :: rad2deg = 180._r8/pi ! radians to degrees conversion factor
163 : real(r8) :: xlat
164 :
165 2978352 : real(r8) :: mee_ap_ipr(ncol,pver) ! ion pairs production rate from Ap formulation
166 :
167 2314006344 : call mee_ionpairs(ncol, lchnk, pmid, zmid*1.e3_r8, tfld, mee_ap_ipr)
168 :
169 25455558960 : extfrc(:,:,:) = 0._r8
170 :
171 2314006344 : no_lgt(:,:) = 0._r8
172 :
173 1489176 : if(use_hemco) then
174 : !--------------------------------------------------------
175 : ! ... set frcing from datasets (HEMCO)
176 : !--------------------------------------------------------
177 0 : call hco_set_extfrc( lchnk, zint_rel, extfrc, ncol, pbuf )
178 : else
179 : !--------------------------------------------------------
180 : ! ... set frcing from datasets
181 : !--------------------------------------------------------
182 1489176 : call extfrc_set( lchnk, zint_rel, extfrc, ncol )
183 : endif
184 :
185 : !--------------------------------------------------------
186 : ! ... set nox production from lighting
187 : ! note: from ground to cloud top production is c shaped
188 : !--------------------------------------------------------
189 1489176 : if ( no_ndx > 0 ) then
190 0 : do i = 1,ncol
191 0 : cldind = nint( cldtop(i) )
192 0 : if( cldind < pver .and. cldind > 0 ) then
193 0 : extfrc(i,cldind:pver,no_ndx) = extfrc(i,cldind:pver,no_ndx) &
194 0 : + prod_no(i,cldind:pver,lchnk)
195 0 : no_lgt(i,cldind:pver) = prod_no(i,cldind:pver,lchnk)
196 : end if
197 : end do
198 : endif
199 1489176 : if ( xno_ndx > 0 ) then
200 0 : do i = 1,ncol
201 0 : cldind = nint( cldtop(i) )
202 0 : if( cldind < pver .and. cldind > 0 ) then
203 0 : extfrc(i,cldind:pver,xno_ndx) = extfrc(i,cldind:pver,xno_ndx) &
204 0 : + prod_no(i,cldind:pver,lchnk)
205 : end if
206 : end do
207 : endif
208 :
209 1489176 : call outfld( 'NO_Lightning', no_lgt(:ncol,:), ncol, lchnk )
210 :
211 1489176 : call airpl_set( lchnk, ncol, no_ndx, co_ndx, xno_ndx, cldtop, zint_abs, extfrc)
212 :
213 1489176 : do i = 1,num_tracer_srcs
214 :
215 0 : ndx = get_extfrc_ndx( tracer_src_flds(i) )
216 0 : call get_srcs_data( tracer_src_flds(i), srcs_offline, ncol, lchnk, pbuf )
217 1489176 : do k = 1,pver
218 0 : extfrc(:ncol,k,ndx) = extfrc(:ncol,k,ndx) + srcs_offline(:ncol,k)
219 : enddo
220 :
221 : enddo
222 :
223 : !
224 : ! CCMI : external forcing for AOA_NH
225 : ! set to 100 ppb per year
226 : !
227 1489176 : if ( aoa_nh_ndx > 0 ) then
228 0 : extfrc(:ncol,:,aoa_nh_ndx) = 1.e-7_r8/(86400._r8*365._r8)
229 0 : do i=1,ncol
230 0 : xlat = rlats(i)*rad2deg ! convert to degrees
231 0 : if ( xlat >= 30._r8 .and. xlat <= 50._r8 ) then
232 0 : extfrc(i,pver,aoa_nh_ndx) = 0._r8
233 : end if
234 : end do
235 0 : call outfld('AOA_NH_XFRC',extfrc(:ncol,:,aoa_nh_ndx), ncol, lchnk )
236 : end if
237 :
238 1489176 : if ( has_ions ) then
239 : !---------------------------------------------------------------------
240 : ! ... set ion auroral production
241 : !---------------------------------------------------------------------
242 :
243 : call aurora( tfld, o2mmr, ommr, mbar, rlats, &
244 : extfrc(:,:,o2p_ndx), extfrc(:,:,op_ndx), extfrc(:,:,n2p_ndx), extfrc(:,:,np_ndx), pmid, &
245 0 : lchnk, calday, ncol, rlons, pbuf )
246 : !---------------------------------------------------------------------
247 : ! ... set n(2d) and n(4s) auroral production
248 : ! Stan Solomon HAO
249 : ! include production of N by secondary auroral hot electrons (e_s*):
250 : ! this is not a "real" reaction; instead, the production is parameterized in terms
251 : ! of the production rate of N2+ by primary electrons, QN2P (which is in the model),
252 : ! as follows:
253 : !---------------------------------------------------------------------
254 0 : do k = 1,pver
255 0 : extfrc(:,k,n2d_ndx) = 1.57_r8*.6_r8*extfrc(:,k,n2p_ndx)
256 0 : extfrc(:,k,n_ndx) = 1.57_r8*.4_r8*extfrc(:,k,n2p_ndx)
257 : end do
258 :
259 : endif
260 :
261 : !---------------------------------------------------------------------
262 : ! ... set SPE NOx and HOx production
263 : ! Jackman et al., JGR, 2005
264 : ! production of 1.25 Nitrogen atoms/ion pair with branching ratios
265 : ! of 0.55 N(4S) and 0.7 N(2D).
266 : !---------------------------------------------------------------------
267 : !---------------------------------------------------------------------
268 : ! ion pairs production due to Galactic Cosmic Rays
269 : !---------------------------------------------------------------------
270 1489176 : call gcr_ionization_ionpairs( ncol, lchnk, gcr_ipr )
271 1489176 : call outfld( 'GCR_ionpairs', gcr_ipr, ncol, lchnk )
272 :
273 : !---------------------------------------------------------------------
274 : ! ion pairs production due to Energetic Particle Precipitation
275 : !---------------------------------------------------------------------
276 1489176 : call epp_ionization_ionpairs( ncol, lchnk, pmid, tfld, epp_ipr )
277 1489176 : call outfld( 'EPP_ionpairs', epp_ipr, ncol, lchnk )
278 :
279 2314006344 : epp_ipr(:ncol,:pver) = epp_ipr(:ncol,:) + gcr_ipr(:ncol,:) + mee_ap_ipr(:ncol,:)
280 :
281 1489176 : if (has_dregion_ions) then
282 : ! D-region ion chemistry is active ...
283 : ! N2p production
284 0 : extfrc(:ncol,:pver,n2p_ndx) = extfrc(:ncol,:pver,n2p_ndx) + 0.585_r8 * epp_ipr(:ncol,:pver)
285 : ! O2p production
286 0 : extfrc(:ncol,:pver,o2p_ndx) = extfrc(:ncol,:pver,o2p_ndx) + 0.15_r8 * epp_ipr(:ncol,:pver)
287 : ! Np
288 0 : extfrc(:ncol,:pver,np_ndx) = extfrc(:ncol,:pver,np_ndx) + 0.185_r8 * epp_ipr(:ncol,:pver)
289 : ! Op
290 0 : extfrc(:ncol,:pver,op_ndx) = extfrc(:ncol,:pver,op_ndx) + 0.076_r8 * epp_ipr(:ncol,:pver)
291 : ! N2D/N4S branching
292 : ! new initial rates
293 0 : extfrc(:ncol,:pver,n2d_ndx) = extfrc(:ncol,:pver,n2d_ndx) + 0.583_r8 * epp_ipr(:ncol,:pver)
294 0 : extfrc(:ncol,:pver,n_ndx) = extfrc(:ncol,:pver,n_ndx) + 0.502_r8 * epp_ipr(:ncol,:pver)
295 : ! O
296 0 : extfrc(:ncol,:pver,o_ndx) = extfrc(:ncol,:pver,o_ndx) + 1.074_r8 * epp_ipr(:ncol,:pver)
297 :
298 : else
299 : ! D-region ion chemistry is NOT active
300 1489176 : if ( n2d_ndx>0 .and. n_ndx>0 ) then
301 0 : extfrc(:ncol,:pver,n2d_ndx) = extfrc(:ncol,:pver,n2d_ndx) + 0.7_r8*epp_ipr(:ncol,:pver)
302 0 : extfrc(:ncol,:pver, n_ndx) = extfrc(:ncol,:pver, n_ndx) + 0.55_r8*epp_ipr(:ncol,:pver)
303 0 : call outfld( 'N2D_EPP', 0.7_r8*epp_ipr(:ncol,:), ncol, lchnk ) ! N(2D) produciton (molec/cm3/s)
304 0 : call outfld( 'N4S_EPP',0.55_r8*epp_ipr(:ncol,:), ncol, lchnk ) ! N(4S) produciton (molec/cm3/s)
305 1489176 : elseif ( no_ndx>0 .and. n_ndx>0 ) then
306 : ! for mechanisms that do not include N2D -- the EPP produce NO
307 0 : extfrc(:ncol,:pver, no_ndx) = extfrc(:ncol,:pver, no_ndx) + 0.7_r8*epp_ipr(:ncol,:pver)
308 0 : extfrc(:ncol,:pver, n_ndx) = extfrc(:ncol,:pver, n_ndx) + 0.55_r8*epp_ipr(:ncol,:pver)
309 0 : call outfld( 'NO_EPP', 0.7_r8*epp_ipr(:ncol,:), ncol, lchnk ) ! NO produciton (molec/cm3/s)
310 0 : call outfld( 'N4S_EPP',0.55_r8*epp_ipr(:ncol,:), ncol, lchnk ) ! N(4S) produciton (molec/cm3/s)
311 : endif
312 1489176 : if ( oh_ndx > 0 ) then
313 0 : do i = 1,ncol
314 0 : epp_hox(i,:pver) = epp_ipr(i,:pver)*hox_prod_factor( epp_ipr(i,:pver), zmid(i,:pver) )
315 : end do
316 0 : extfrc(:ncol,:pver, oh_ndx) = extfrc(:ncol,:pver, oh_ndx) + epp_hox(:ncol,:pver)
317 0 : call outfld( 'OH_EPP' , epp_hox(:ncol,:), ncol, lchnk ) ! HOX produciton (molec/cm3/s)
318 : endif
319 : endif
320 :
321 1489176 : if ( has_ions ) then
322 : !---------------------------------------------------------------------
323 : ! ... set total electron production
324 : !---------------------------------------------------------------------
325 0 : do k = 1,pver
326 0 : extfrc(:,k,e_ndx) = extfrc(:,k,op_ndx) + extfrc(:,k,o2p_ndx) &
327 0 : + extfrc(:,k,np_ndx) + extfrc(:,k,n2p_ndx)
328 : end do
329 0 : call outfld( 'P_Op', extfrc(:,:,op_ndx), ncol, lchnk )
330 0 : call outfld( 'P_O2p', extfrc(:,:,o2p_ndx), ncol, lchnk )
331 0 : call outfld( 'P_Np', extfrc(:,:,np_ndx), ncol, lchnk )
332 0 : call outfld( 'P_N2p', extfrc(:,:,n2p_ndx), ncol, lchnk )
333 0 : call outfld( 'P_IONS',extfrc(:,:,e_ndx), ncol, lchnk )
334 : end if
335 :
336 1489176 : end subroutine setext
337 :
338 : end module mo_setext
|
