Line data Source code
1 : !=======================================================================================================
2 : ! extracted from TIEGCM model
3 : ! see section 6.3 in https://www.hao.ucar.edu/modeling/tgcm/doc/description/model_description.pdf
4 : !=======================================================================================================
5 : module tei_mod
6 : !
7 : ! Calculate electron and ion temperatures.
8 : !
9 : use shr_kind_mod, only : r8 => shr_kind_r8
10 : use shr_const_mod, only : pi => shr_const_pi ! Boltzmann constant and pi
11 : use mo_constants, only : gask => rgas_cgs
12 : use mo_constants, only : boltz=>boltz_cgs
13 : use physconst, only : avogad ! molecules/kmole
14 :
15 : implicit none
16 : private
17 : public :: settei
18 :
19 : ! g = 8.7 m/s at 400 km altitude
20 : real(r8), parameter :: grav = 870._r8 ! (cm/s^2)
21 : real(r8), parameter :: dipmin = 0.24_r8 ! minimum mag dip angle (2.5 deg horizontal res)
22 : real(r8), parameter :: rtd = 180._r8/pi
23 : real(r8), parameter :: evergs = 1.602e-12_r8 ! 1 eV = 1.602e-12 ergs
24 :
25 : real(r8), parameter :: rmassinv_o2 = 1._r8/32._r8
26 : real(r8), parameter :: rmassinv_o1 = 1._r8/16._r8
27 : real(r8), parameter :: rmassinv_n2 = 1._r8/28._r8
28 : real(r8), parameter :: avo = avogad*1.e-3_r8 ! molecules/mole
29 : contains
30 :
31 : !
32 0 : subroutine settei(tn,o2,o1,n2,ne,te,ti,op,o2p,nop,barm,qji_ti, &
33 0 : f107, chi, qtot, qteaur, rlatm, dipmag, pmid, pint, lev0, lev1,ncol, &
34 0 : te_out, ti_out, qtotal )
35 :
36 : use trsolv_mod, only : trsolv
37 : !
38 : ! Args:
39 : integer,intent(in) :: lev0,lev1,ncol
40 : real(r8),dimension(lev0:lev1,ncol),intent(in) :: &
41 : tn, & ! neutral temperature (deg K)
42 : o2, & ! molecular oxygen (mmr)
43 : o1, & ! atomic oxygen (mmr)
44 : n2, & ! molecular nitrogen (mmr)
45 : ne, & ! electron density (cm3)
46 : te, & ! electron temperature (from previous time step) (K)
47 : ti, & ! ion temperature (from previous time step) (K)
48 : op, & ! O+ number dens (/cm^3)
49 : o2p, & ! O2+ number dens (/cm^3)
50 : nop, & ! NO+ number dens (/cm^3)
51 : barm, & ! mean molecular weight (g/mole)
52 : qji_ti ! joule heating from qjoule_ti (used ui,vi) ergs/s/g
53 :
54 : real(r8), intent(in) :: f107
55 : real(r8), intent(in) :: chi(ncol) ! solar zenith angle (radians)
56 : real(r8), intent(in) :: qtot(lev0:lev1,ncol) ! total ionization rate (s-1 cm-3)
57 : real(r8), intent(in) :: qteaur(ncol) ! (units ?? )
58 : real(r8), intent(in) :: rlatm(ncol) ! geo mag latitude (radians)
59 : real(r8), intent(in) :: dipmag(ncol) ! geo mag dip angle (radians)
60 : real(r8), intent(in) :: pmid(lev0:lev1,ncol) ! mid-level press (dyne/cm^2)
61 : real(r8), intent(in) :: pint(lev0:lev1+1,ncol) ! interface press (dyne/cm^2)
62 :
63 : !
64 : ! Output args:
65 : real(r8),dimension(lev0:lev1,ncol),intent(out) :: &
66 : te_out, & ! output electron temperature (deg K)
67 : ti_out ! output ion temperature (deg K)
68 : real(r8), intent(out) :: qtotal(lev0:lev1,ncol)! (ergs/sec/gm)
69 :
70 : !
71 : ! Local:
72 : integer :: i,k,n,ier
73 : integer :: nk,nkm1
74 : real(r8),dimension(lev0:lev1,ncol) :: &
75 0 : xnmbar, & ! p0*e(-z)*barm/kT (minpoints or interfaces)
76 0 : p_coef, & ! coefficient for trisolv (s1)
77 0 : q_coef, & ! coefficient for trisolv (s2)
78 0 : r_coef, & ! coefficient for trisolv (s3)
79 0 : rhs ! right-hand-side for trisolv (s4)
80 : real(r8),dimension(lev0:lev1) :: &
81 0 : te_int, & ! electron temperature (interfaces)
82 0 : tn_int, & ! neutral temperature (interfaces)
83 0 : o2n, & ! O2 number density (interfaces)
84 0 : o1n, & ! O1 number density (interfaces)
85 0 : n2n, & ! N2 number density (interfaces)
86 0 : root_te, & ! sqrt(te)
87 0 : root_tn, & ! sqrt(tn)
88 0 : root_ne, & ! sqrt(ne)
89 0 : tek0, & ! ke/te**2.5 (s15)
90 0 : h_mid,h_int,&
91 0 : barm_int,&
92 : fki, & ! work array
93 0 : qe, & ! source term (s10)
94 0 : q_eni, & ! heating from electron/neutral and electron/ion collisions
95 0 : coll_en2v, & ! electron/N2vib collision (s9)
96 : !
97 : ! Cooling rates (heat loss):
98 0 : loss_en2v, & ! electron/N2vib loss term (s10)
99 0 : loss_en2, & ! electron/N2 loss
100 0 : loss_eo2, & ! electron/O2 loss
101 0 : loss_eo1d, & ! electron/O(1d) loss
102 0 : loss_eo1, & ! electron/O loss
103 0 : loss_xen, & ! L0*(E,N) (s8)
104 0 : loss_en ! electrons/neutrals loss (s11)
105 : real(r8),dimension(lev0:lev1,ncol) :: &
106 0 : loss_ei, & ! electron/ion loss (s10)
107 0 : loss_in ! ion/neutral loss (s9)
108 : real(r8),parameter :: &
109 : fpolar = -3.0e+9_r8, & ! polar te flux
110 : del = 1.e-6_r8 , & ! some small value
111 : root2 = sqrt(2._r8), &
112 : !
113 : ! Correction factors for neutral heating due to L(E,O1D)
114 : alam = 0.0069_r8 , &
115 : ad = 0.0091_r8 , &
116 : sd = 2.3e-11_r8
117 : real(r8) :: &
118 : f107te ! solar flux
119 : !
120 : ! a,fed,fen,fe,sindipmag have a z dimension only for diagnostic plotting:
121 : real(r8) :: &
122 : a,fed,fen,& ! day/night
123 : fe, & ! heat flux at upper boundary
124 : sindipmag ! sin(dipmag)
125 : !
126 0 : real(r8) :: dellnp(lev0:lev1) ! delta of log interface pressures
127 :
128 0 : qtotal(:,:) = 0._r8
129 :
130 0 : col_loop: do i=1,ncol
131 0 : do k=lev0,lev1
132 0 : dellnp(k) = abs(log(pint(k,i)) - log(pint(k+1,i)))
133 : end do
134 : !
135 0 : f107te = f107
136 0 : if (f107te > 235._r8) f107te = 235._r8
137 0 : nk = lev1-lev0+1
138 0 : nkm1 = nk-1
139 : !
140 0 : if (abs(rlatm(i))-pi/4.5_r8 >= 0._r8) then
141 : a = 1._r8
142 : else
143 0 : a = .5_r8*(1._r8+sin(pi*(abs(rlatm(i))-pi/9._r8)/(pi/4.5_r8)))
144 : endif
145 : !
146 : ! Increased heat flux for TE fom protonosphere.
147 0 : fed = ( -5.0e+7_r8*f107te*a-4.0e+7_r8*f107te)*1.2_r8
148 0 : fen = fed/2._r8
149 0 : fed = fed+qteaur(i) ! t4
150 0 : fen = fen+qteaur(i) ! t5
151 0 : if (chi(i)-.5_r8*pi >= 0._r8) then ! chi==t2
152 : fe = fen ! t1
153 : else
154 0 : fe = fed
155 : endif
156 0 : if ((chi(i)*rtd-80._r8)*(chi(i)*rtd-100._r8)>=0._r8) then
157 0 : fe = fe*evergs
158 : else
159 : fe = (.5_r8*(fed+fen)+.5_r8*(fed-fen)* &
160 0 : cos(pi*(chi(i)*rtd-80._r8)/20._r8))*evergs
161 : endif
162 : !
163 : ! Add fpolar if magnetic latitude >= 60 degrees:
164 0 : if (abs(rlatm(i)-pi/3._r8)>=0._r8) then
165 0 : fe = fe+fpolar*evergs
166 : end if
167 :
168 : !
169 : ! te,o2,o,n2,tn at interfaces:
170 0 : do k=lev0+1,lev1-1
171 0 : te_int(k) = .5_r8*(te(k,i)+te(k-1,i))
172 0 : o2n(k) = .5_r8*(o2(k,i)+o2(k-1,i))
173 0 : o1n(k) = .5_r8*(o1(k,i)+o1(k-1,i))
174 0 : n2n(k) = .5_r8*(n2(k,i)+n2(k-1,i))
175 0 : tn_int(k) = .5_r8*(tn(k,i)+tn(k-1,i))
176 0 : barm_int(k) = .5_r8*(barm(k,i)+barm(k-1,i))
177 : enddo ! k=lev0+1,lev1-2
178 : !
179 : ! Bottom:
180 0 : te_int(lev0) = 1.5_r8*te(lev0,i)-.5_r8*te(lev0+1,i)
181 0 : o2n(lev0) = 1.5_r8*o2(lev0,i)-.5_r8*o2(lev0+1,i)
182 0 : o1n(lev0) = 1.5_r8*o1(lev0,i)-.5_r8*o1(lev0+1,i)
183 0 : n2n(lev0) = 1.5_r8*n2(lev0,i)-.5_r8*n2(lev0+1,i)
184 0 : tn_int(lev0) = 1.5_r8*tn(lev0,i)-.5_r8*tn(lev0+1,i)
185 0 : barm_int(lev0) = 1.5_r8*barm(lev0,i)-.5_r8*barm(lev0+1,i)
186 : !
187 : ! Top:
188 0 : te_int(lev1) = 1.5_r8*te(lev1-1,i)-.5_r8*te(lev1-2,i)
189 0 : o2n(lev1) = 1.5_r8*o2(lev1-1,i)-.5_r8*o2(lev1-2,i)
190 0 : o1n(lev1) = 1.5_r8*o1(lev1-1,i)-.5_r8*o1(lev1-2,i)
191 0 : n2n(lev1) = 1.5_r8*n2(lev1-1,i)-.5_r8*n2(lev1-2,i)
192 0 : tn_int(lev1) = 1.5_r8*tn(lev1-1,i)-.5_r8*tn(lev1-2,i)
193 0 : barm_int(lev1) = 1.5_r8*barm(lev1-1,i)-.5_r8*barm(lev1-2,i)
194 : !
195 : ! N2:
196 0 : do k=lev0,lev1
197 0 : if (n2n(k) < 0._r8) n2n(k) = 0._r8
198 : enddo ! k=lev0,lev1
199 :
200 : !
201 : ! Convert o2,o,n2 to number density (interfaces):
202 :
203 0 : do k=lev0,lev1
204 : ! xnmbar at interfaces:
205 0 : xnmbar(k,i) = pint(k,i)*barm_int(k)/(boltz*tn_int(k)) ! s8
206 0 : o2n(k) = xnmbar(k,i)*o2n(k)*rmassinv_o2 ! s13
207 0 : o1n(k) = xnmbar(k,i)*o1n(k)*rmassinv_o1 ! s12
208 0 : n2n(k) = xnmbar(k,i)*n2n(k)*rmassinv_n2 ! s11
209 0 : root_te(k) = sqrt(te_int(k))
210 : !
211 : tek0(k) = 7.5e5_r8/ &
212 : (1._r8+3.22e4_r8*te_int(k)**2/ &
213 : ne(k,i)*(root_te(k)* &
214 : (2.82e-17_r8 - 3.41e-21_r8 * te_int (k))*n2n(k)+ &
215 : (2.20e-16_r8 + 7.92e-18_r8 * root_te(k))*o2n(k)+ &
216 0 : 1.10e-16_r8 * (1._r8+5.7e-4_r8 * te_int (k))*o1n(k)))*evergs
217 :
218 : enddo ! k=lev0,lev1
219 :
220 0 : do k=lev0,lev1-1
221 0 : h_mid(k) = gask*tn(k,i)/(barm(k,i)*grav) ! s7
222 : enddo ! k=lev0,lev1-1
223 0 : do k=lev0,lev1
224 0 : h_int(k) = gask*tn_int(k)/(barm_int(k)*grav) ! s6
225 : enddo ! k=lev0,lev1
226 :
227 0 : if (abs(dipmag(i)) >= dipmin) then
228 0 : sindipmag = (sin(dipmag(i)))**2 ! t2,s2
229 : else
230 : sindipmag = (sin(dipmin))**2
231 : endif
232 0 : if (sindipmag < .10_r8) sindipmag = .10_r8
233 : !
234 : ! Start coefficients and rhs for trsolv:
235 0 : do k=lev0,lev1-1
236 0 : p_coef(k,i) = 2._r8/7._r8*sindipmag/(h_mid(k)*dellnp(k)**2) ! s1
237 0 : r_coef(k,i) = p_coef(k,i)*tek0(k+1)/h_int(k+1) ! s3
238 0 : p_coef(k,i) = p_coef(k,i)*tek0(k )/h_int(k ) ! s1
239 0 : q_coef(k,i) = -(p_coef(k,i)+r_coef(k,i)) ! s2
240 0 : rhs(k,i) = 0._r8 ! s4
241 : enddo ! k=lev0,lev1-1
242 : !
243 : ! Bottom boundary:
244 0 : q_coef(lev0,i) = q_coef(lev0,i)-p_coef(lev0,i)
245 0 : rhs(lev0,i) = rhs(lev0,i)-2._r8*p_coef(lev0,i)*tn_int(lev0)**3.5_r8
246 0 : p_coef(lev0,i) = 0._r8
247 : !
248 : ! Upper boundary:
249 0 : q_coef(lev1-1,i) = q_coef(lev1-1,i)+r_coef(lev1-1,i)
250 : rhs(lev1-1,i) = rhs(lev1-1,i)+r_coef(lev1-1,i)*dellnp(lev1-1)*3.5_r8* &
251 0 : h_int(lev1)*fe/tek0(lev1)
252 0 : r_coef(lev1-1,i) = 0._r8
253 :
254 : !
255 : !
256 : ! Set Ne (midpoints "(K+1/2)"):
257 : !
258 0 : do k=lev0,lev1-1
259 0 : root_ne(k) = ne(k,i)*ne(k+1,i)
260 0 : if (root_ne(k) < 1.e4_r8) root_ne(k) = 1.e4_r8
261 0 : root_ne(k) = sqrt(root_ne(k))
262 : enddo ! k=lev0,lev1-1
263 :
264 : !
265 : ! Set up o2,o,n2 number densities at midpoints:
266 : !
267 :
268 0 : do k=lev0,lev1-1
269 0 : xnmbar(k,i) = pmid(k,i)*barm(k,i)/(boltz*tn(k,i))
270 0 : o2n(k) = xnmbar(k,i)*o2(k,i)*rmassinv_o2 ! s14
271 0 : o1n(k) = xnmbar(k,i)*o1(k,i)*rmassinv_o1 ! s13
272 0 : n2n(k) = n2(k,i)
273 0 : if (n2n(k) < 0._r8) n2n(k) = 0._r8
274 0 : n2n(k) = xnmbar(k,i)*n2n(k)*rmassinv_n2 ! s12
275 : !
276 : ! Calculate source term qe (s10)
277 : ! Comment from earlier version (maybe the *1.0 below was once *2.0):
278 : ! "Correction facor of 2 increase in TE heating rate"
279 : !
280 0 : qe(k) = log(root_ne(k)/(o2n(k)+n2n(k)+0.1_r8*o1n(k)))
281 : qe(k) = exp(-((((0.001996_r8*qe(k)+0.08034_r8)*qe(k)+1.166_r8)* &
282 0 : qe(k)+6.941_r8)*qe(k)+12.75_r8))*1.0_r8
283 : !
284 : ! Subtract qe from right-hand-side:
285 0 : rhs(k,i) = rhs(k,i)-qe(k)*qtot(k,i)*evergs
286 0 : root_te(k) = sqrt(te(k,i))
287 : !
288 : ! Electron/N2 collision A(E,N2,VIB) (s9):
289 : !
290 0 : if (te(k,i) >= 1000._r8) then
291 0 : coll_en2v(k) = 2.e-7_r8*exp(-4605.2_r8/te(k,i))
292 : else
293 0 : coll_en2v(k) = 5.71e-8_r8*exp(-3352.6_r8/te(k,i))
294 : endif
295 0 : if (te(k,i) > 2000._r8) then
296 0 : coll_en2v(k) = 2.53e-6_r8*root_te(k)*exp(-17620._r8/te(k,i))
297 : end if
298 : !
299 : ! Loss due to electron/n2 collision L0(E,N2,VIB)/(NE*N(N2)) (s10)
300 : !
301 0 : loss_en2v(k) = 3200._r8*(1._r8/te(k,i)-1._r8/tn(k,i))
302 0 : loss_en2v(k) = sign(abs(loss_en2v(k))+del,loss_en2v(k)) ! avoid divide by zero in the next line when te = tn
303 : ! by adding a small value to loss_en2v
304 : ! this assumes te >= tn
305 0 : loss_en2v(k) = -3200._r8/(te(k,i)*tn(k,i))*(1._r8-exp(loss_en2v(k)))/loss_en2v(k)
306 0 : loss_en2v(k) = 1.3e-4_r8*loss_en2v(k)*coll_en2v(k)
307 : enddo ! k=lev0,lev1-1
308 :
309 : !
310 : ! Calculate and sum cooling rates (heat loss) due to interactions between
311 : ! electrons/neutrals, electrons/ions, ions/neutrals
312 : !
313 0 : do k=lev0,lev1-1
314 : !
315 : ! Electron/N2 loss rate:
316 : ! loss_en2 = (L0(E,N2)+L0(E,N2,ROT)+L0(E,N2,VIB))/NE (s11)
317 : !
318 0 : loss_en2(k) = n2n(k)*(1.77E-19_r8*(1._r8-1.21E-4_r8*te(k,i))* te(k,i) + 2.9e-14_r8/root_te(k) + loss_en2v(k))
319 : !
320 : ! Start total of electron/neutral loss rate (s11):
321 : !
322 0 : loss_en(k) = loss_en2(k)
323 : !
324 : ! Electron/O2 loss rates: (L0(E,O2)+L0(E,O2,ROT)+L0(E,O2,VIB)/NE
325 : !
326 : loss_eo2(k) = o2n(k)*(1.21e-18_r8*(1._r8+3.6e-2_r8*root_te(k))* &
327 0 : root_te(k)+6.9e-14_r8/root_te(k)+3.125e-21_r8*te(k,i)**2)
328 0 : loss_en(k) = loss_en(k)+loss_eo2(k)
329 : !
330 : ! Electron/O(1d) loss rates: L0(E,O,1D)/(NE*N(O))
331 : !
332 0 : loss_eo1d(k) = 22713._r8*(1._r8/te(k,i)-1._r8/tn(k,i))
333 0 : loss_eo1d(k) = sign(abs(loss_eo1d(k))+del,loss_eo1d(k))
334 0 : loss_eo1d(k) = 22713._r8/(te(k,i)*tn(k,i))*(1._r8-exp(loss_eo1d(k)))/loss_eo1d(k)
335 : !
336 : ! loss_eo1d function often fails here with bad argument to exp()
337 : ! due to high te and/or high loss_eo1d from above.
338 : ! loss_eo1d(k) = 1.57e-12*exp((2.4e4+0.3*(te(k)-1500.)-
339 : ! 1.947e-5*(te(k)-1500.)*(te(k)-4000.))*(te(k)-3000.)/
340 : ! (3000.*te(k)))*loss_eo1d(k)
341 0 : loss_eo1d(k) = 0._r8
342 : !
343 : ! Electron/O1 loss rates: (L0(E,O)+L0(E,O,F))/NE
344 : !
345 : loss_eo1(k) = o1n(k)*(7.9e-19_r8*(1._r8+5.7e-4_r8*te(k,i))* &
346 : root_te(k)+3.4e-12_r8*(1._r8-7.e-5_r8*te(k,i))/tn(k,i)* &
347 0 : (150._r8/te(k,i)+0.4_r8))
348 :
349 0 : loss_en(k) = loss_en(k)+loss_eo1(k)
350 : !
351 : ! loss_xen = L0*(E,N) (s8)
352 : !
353 : loss_xen(k) = (loss_en(k)+o1n(k)*(1._r8-alam/(ad+sd* &
354 0 : n2n(k)))*loss_eo1d(k))*root_ne(k)*evergs
355 : !
356 : ! Complete total electron/neutral loss rate L0(E,N) (s11):
357 : !
358 : loss_en(k) = (loss_en(k)+o1n(k)*loss_eo1d(k))* &
359 0 : root_ne(k)*evergs
360 : !
361 : ! Calculate L0(E) = L(E)/(TE-TI), where L(E) is loss due to
362 : ! interactions between electrons and ions.
363 : !
364 : loss_ei(k,i) = 3.2e-8_r8*root_ne(k)/(root_te(k)*te(k,i))* &
365 0 : 15._r8*(op(k,i)+0.5_r8*o2p(k,i)+0.53_r8*nop(k,i))*evergs
366 :
367 0 : root_tn(k) = sqrt(tn(k,i))
368 : !
369 : ! loss_in = ion/neutral cooling = L0(I,N) =L(I,N)/(TI-TN)
370 : !
371 : loss_in(k,i) = ((6.6e-14_r8*n2n(k)+5.8e-14_r8*o2n(k)+0.21e-14_r8* &
372 : o1n(k)*root2*root_tn(k))*op(k,i)+(5.45e-14_r8*o2n(k)+ &
373 : 5.9e-14_r8*n2n(k)+4.5e-14_r8*o1n(k))*nop(k,i)+(5.8e-14_r8* &
374 : n2n(k)+4.4e-14_r8*o1n(k)+0.14e-14_r8*o2n(k)*root_tn(k))* &
375 0 : o2p(k,i))*evergs
376 : !
377 : ! Complete tridiagonal matrix coefficients and rhs:
378 : !
379 : ! q_coef = q_coef-(L0(E,N)+L0(E,I))/TE**2.5 = Q
380 : !
381 0 : q_coef(k,i) = q_coef(k,i)-(loss_en(k)+loss_ei(k,i))/te(k,i)**2.5_r8
382 : !
383 : ! rhs = rhs-L0(E,N)*TN-L0(E)*TI
384 : !
385 0 : rhs(k,i) = rhs(k,i)-loss_en(k)*tn(k,i)-loss_ei(k,i)*ti(k,i)
386 :
387 : enddo ! k=lev0,lev1-1
388 :
389 : ! Calculate heating due to electron/neutral and electron/ion collisions
390 : ! (ergs/sec/gm):
391 : !
392 0 : do k=lev0,lev1-1
393 0 : if (te(k,i)-ti(k,i) >= 0._r8) then
394 0 : q_eni(k)=loss_ei(k,i)*(te(k,i)-ti(k,i))
395 : else
396 0 : q_eni(k) = 0._r8
397 : endif
398 0 : q_eni(k) = (loss_xen(k)*(te(k,i)-tn(k,i))+q_eni(k)) * avo/xnmbar(k,i)
399 : enddo ! k=lev0,lev1-1
400 :
401 : !
402 : ! Add collisional heating to Q for use in thermodynamic equation.
403 0 : do k=lev0,lev1-2
404 0 : qtotal(k+1,i) = qtotal(k+1,i)+.5_r8*(q_eni(k)+q_eni(k+1))
405 : enddo ! k=lev0,lev1-2
406 : !
407 : ! Upper and lower boundaries:
408 0 : qtotal(lev0,i) = qtotal(lev0,i)+1.5_r8*q_eni(lev0)-0.5_r8*q_eni(lev0+1)
409 0 : qtotal(lev1,i) = qtotal(lev1,i)+1.5_r8*q_eni(lev1-1)-0.5_r8*q_eni(lev1-2)
410 : end do col_loop
411 :
412 : !
413 : ! Solve tridiagonal system:
414 : !
415 0 : call trsolv(p_coef,q_coef,r_coef,rhs,te_out, lev0,lev1, lev0,lev1-1, 1,ncol )
416 :
417 0 : col_loop2: do i=1,ncol
418 : !
419 : ! Te = Te**(2./7.):
420 0 : do k=lev0,lev1-1
421 0 : te_out(k,i) = te_out(k,i)**(2._r8/7._r8)
422 : enddo
423 : !
424 : ! 10/21/03 btf: make this check after te*(2/7), rather than before.
425 : !
426 : ! Te must be >= Tn:
427 0 : do k=lev0,lev1-1
428 0 : if (te_out(k,i) < tn(k,i)) te_out(k,i) = tn(k,i)
429 : enddo
430 : !
431 : ! 1/9/08 btf: put spval in top level of te:
432 : ! te_out(lev1) = spval
433 0 : te_out(lev1,i) = te_out(lev1-1,i)
434 : !
435 : ! Set ion temperature output. Use joule heating qji_ti from sub
436 : ! qjoule_ti (see qjoule.F). lev1 not calculated.
437 : !
438 0 : do k=lev0,lev1-1
439 0 : ti_out(k,i) = (qji_ti(k,i)*(xnmbar(k,i)/avo)+ &
440 : loss_ei(k,i)*te_out(k,i)+loss_in(k,i)*tn(k,i))/&
441 0 : (loss_ei(k,i)+loss_in(k,i))
442 : !
443 : ! ti must be at least as large as tn:
444 0 : if (ti_out(k,i) < tn(k,i)) ti_out(k,i) = tn(k,i)
445 : enddo
446 : !
447 : ! 1/9/08 btf: put spval in top level of ti:
448 : ! ti_out(lev1) = spval
449 0 : ti_out(lev1,i) = ti_out(lev1-1,i)
450 : end do col_loop2
451 :
452 0 : end subroutine settei
453 :
454 : end module tei_mod
455 :
|
