Context Navigation

iono_h.F90 @ 3567

Last change on this file since 3567 was 2969, checked in by emillour, 19 months ago

Venus PCM:
Some adaptations for gfortran:

the cosd()/sind() functions (arguments in degrees) are depreciated use cos()/sin() instead with a conversion to radians of input arguments
comparisons between logicals is done with .eqv. not .eq.

File size: 27.7 KB

Line
1	MODULE iono_h
2
3	IMPLICIT NONE
4
5	CONTAINS
6
7	subroutine phdisrate(ig,nlayer,chemthermod,zenit,i)
8
9	! Calculates photoionization and photodissociation rates from the
10	! photoabsorption rates calculated in jthermcalc_e107 and the
11	! ionization/dissociation branching ratios in param_read_e107
12
13	! apr 2002 fgg first version
14	!**********************************************************************
15
16	use param_v4_h, only: ninter,nabs, &
17	jfotsout,fluxtop, &
18	jion,jdistot,jdistot_b, &
19	efdisco2,efdiso2,efdish2o, &
20	efdish2o2,efdish2,efdiso3, &
21	efdiso,efdisn,efdish, &
22	efdisno,efdisn2,efdisno2, &
23	efdisco,efionco2,efiono2,efionn2, &
24	efionco,efiono3p,efionn, &
25	efionno,efionh
26	! &
27	implicit none
28
29	! arguments
30
31	integer i !altitude
32	integer ig,chemthermod,nlayer
33	real zenit
34
35	! local variables
36
37	integer inter,iz,j
38	real lambda
39	real jdis(nabs,ninter,nlayer)
40	character*1 dn
41
42
43	!**********************************************************************
44
45	! photodissociation and photoionization rates initialization
46
47	jion(:,i,:) = 0.d0
48	jdistot(:,i) = 0.d0
49	jdistot_b(:,i) = 0.d0
50
51	! jion(1,i,1) = 0.d0 ! CO2 channel 1 ( --> CO2+ + e- )
52	! jion(1,i,2) = 0.d0 ! CO2 channel 2 ( --> O+ + CO + e- )
53	! jion(1,i,3) = 0.d0 ! CO2 channel 3 ( --> CO+ + O + e- )
54	! jion(1,i,4) = 0.d0 ! CO2 channel 4 ( --> C+ + O2 + e- )
55	! jion(2,i,1) = 0.d0 ! O2 (only one ionization channel)
56	! jion(3,i,1) = 0.d0 ! O3P (only one ionization channel)
57	! jion(4,i,1) = 0.d0 ! H2O (no ionization)
58	! jion(5,i,1) = 0.d0 ! H2 (no ionization)
59	! jion(6,i,1) = 0.d0 ! H2O2 (no ionization)
60	! jion(7,i,1) = 0.d0 ! O3 (no ionization)
61	! jion(8,i,1) = 0.d0 ! N2 channel 1 ( --> n2+ + e- )
62	! jion(8,i,2) = 0.d0 ! N2 channel 2 ( --> n+ + n + e- )
63	! jion(9,i,1) = 0.d0 ! N ( --> N+ + e- )
64	! jion(10,i,1)= 0.d0 ! We do not mind its ionization
65	! jion(11,i,1) = 0.d0 ! CO channel 1 ( --> co+ + e- )
66	! jion(11,i,2) = 0.d0 ! CO channel 2 ( --> c+ + o + e- )
67	! jion(11,i,3) = 0.d0 ! CO channel 3 ( --> o+ + c + e- )
68	! jion(12,i,1) = 0.d0 ! H ( --> H+ + e- )
69	! jion(13,i,1) = 0.d0
70	! do j=1,nabs
71	! jdistot(j,i) = 0.
72	! jdistot_b(j,i) = 0.
73	! end do
74
75	if(zenit.gt.140.) then
76	dn='n'
77	else
78	dn='d'
79	end if
80
81	if(dn.eq.'n') return
82
83	! photodissociation and photoionization rates for each species
84
85	do inter=1,ninter
86	! CO2
87	jdis(1,inter,i) = jfotsout(inter,1,i) * fluxtop(inter) &
88	* efdisco2(inter)
89	if(inter.gt.29.and.inter.le.32) then
90	jdistot(1,i) = jdistot(1,i) + jdis(1,inter,i)
91	else if(inter.le.29) then
92	jdistot_b(1,i) = jdistot_b(1,i) + jdis(1,inter,i)
93	end if
94	jion(1,i,1)=jion(1,i,1) + &
95	jfotsout(inter,1,i)fluxtop(inter)efionco2(inter,1)
96	jion(1,i,2)=jion(1,i,2) + &
97	jfotsout(inter,1,i)fluxtop(inter)efionco2(inter,2)
98	jion(1,i,3)=jion(1,i,3) + &
99	jfotsout(inter,1,i)fluxtop(inter)efionco2(inter,3)
100	jion(1,i,4)=jion(1,i,4) + &
101	jfotsout(inter,1,i)fluxtop(inter)efionco2(inter,4)
102
103
104	! O2
105	jdis(2,inter,i) = jfotsout(inter,2,i) * fluxtop(inter) &
106	* efdiso2(inter)
107	if(inter.ge.31) then
108	jdistot(2,i) = jdistot(2,i) + jdis(2,inter,i)
109	else if(inter.eq.30) then
110	jdistot(2,i)=jdistot(2,i)+0.02*jdis(2,inter,i)
111	jdistot_b(2,i)=jdistot_b(2,i)+0.98*jdis(2,inter,i)
112	else if(inter.lt.31) then
113	jdistot_b(2,i) = jdistot_b(2,i) + jdis(2,inter,i)
114	end if
115	jion(2,i,1)=jion(2,i,1) + &
116	jfotsout(inter,2,i) * fluxtop(inter) * efiono2(inter,1)
117	jion(2,1,2)=jion(2,1,2) + &
118	jfotsout(inter,2,i) * fluxtop(inter) * efiono2(inter,2)
119	!(1.-efdiso2(inter))
120
121	! O3P
122	jion(3,i,1)=jion(3,i,1) + &
123	jfotsout(inter,3,i) * fluxtop(inter) * efiono3p(inter)
124
125	! H2O
126	jdis(4,inter,i) = jfotsout(inter,4,i) * fluxtop(inter) &
127	* efdish2o(inter)
128	jdistot(4,i) = jdistot(4,i) + jdis(4,inter,i)
129
130	! H2
131	jdis(5,inter,i) = jfotsout(inter,5,i) * fluxtop(inter) &
132	* efdish2(inter)
133	jdistot(5,i) = jdistot(5,i) + jdis(5,inter,i)
134
135	! H2O2
136	jdis(6,inter,i) = jfotsout(inter,6,i) * fluxtop(inter) &
137	* efdish2o2(inter)
138	jdistot(6,i) = jdistot(6,i) + jdis(6,inter,i)
139
140	!Only if O3, N or ion chemistry requested
141	if(chemthermod.ge.1) then
142	! O3
143	jdis(7,inter,i) = jfotsout(inter,7,i) * fluxtop(inter) &
144	* efdiso3(inter)
145	if(inter.eq.34) then
146	jdistot(7,i) = jdistot(7,i) + jdis(7,inter,i)
147	else if (inter.eq.35) then
148	jdistot(7,i) = jdistot(7,i) + 0.997 * jdis(7,inter,i)
149	jdistot_b(7,i) = jdistot_b(7,i) + 0.003 * jdis(7,inter,i)
150	else if (inter.eq.36) then
151	jdistot_b(7,i) = jdistot_b(7,i) + jdis(7,inter,i)
152	endif
153	endif !Of chemthermod.ge.1
154
155	!Only if N or ion chemistry requested
156	if(chemthermod.ge.2) then
157	! N2
158	jdis(8,inter,i) = jfotsout(inter,8,i) * fluxtop(inter) &
159	* efdisn2(inter)
160	jdistot(8,i) = jdistot(8,i) + jdis(8,inter,i)
161	jion(8,i,1) = jion(8,i,1) + jfotsout(inter,8,i) * &
162	fluxtop(inter) * efionn2(inter,1)
163	jion(8,i,2) = jion(8,i,2) + jfotsout(inter,8,i) * &
164	fluxtop(inter) * efionn2(inter,2)
165
166	! N
167	jion(9,i,1) = jion(9,i,1) + jfotsout(inter,9,i) * &
168	fluxtop(inter) * efionn(inter)
169
170	! NO
171	jdis(10,inter,i) = jfotsout(inter,10,i) * fluxtop(inter) &
172	* efdisno(inter)
173	jdistot(10,i) = jdistot(10,i) + jdis(10,inter,i)
174	jion(10,i,1) = jion(10,i,1) + jfotsout(inter,10,i) * &
175	fluxtop(inter) * efionno(inter)
176
177	! NO2
178	jdis(13,inter,i) = jfotsout(inter,13,i) * fluxtop(inter) &
179	* efdisno2(inter)
180	jdistot(13,i) = jdistot(13,i) + jdis(13,inter,i)
181
182	endif !Of chemthermod.ge.2
183
184	! CO
185	jdis(11,inter,i) = jfotsout(inter,11,i) * fluxtop(inter) &
186	* efdisco(inter)
187	jdistot(11,i) = jdistot(11,i) + jdis(11,inter,i)
188	jion(11,i,1) = jion(11,i,1) + jfotsout(inter,11,i) * &
189	fluxtop(inter) * efionco(inter,1)
190	jion(11,i,2) = jion(11,i,2) + jfotsout(inter,11,i) * &
191	fluxtop(inter) * efionco(inter,2)
192	jion(11,i,3) = jion(11,i,3) + jfotsout(inter,11,i) * &
193	fluxtop(inter) * efionco(inter,3)
194
195	! H
196	jion(12,i,1) = jion(12,i,1) + jfotsout(inter,12,i) * &
197	fluxtop(inter) * efionh(inter)
198
199
200	end do
201
202
203	return
204
205
206	end
207
208	!***********************************************************************
209	function P_indice_factor(sza_local,indice)
210
211	! Computes the P_indice_factor for the electronic temperature of
212	! Theis et al. 1984 model
213
214	!***********************************************************************
215
216	! Arguments
217
218	integer :: indice ! the indiceth factor
219	real :: sza_local ! solar zenith angle in degree
220
221	! local variables:
222
223	logical, save :: firstcall = .true.
224	real,parameter :: pi = 4.*atan(1.0)
225	real,parameter :: deg2rad = pi/180.
226	real, save :: a(4,4), b(4,4) ! Fourrier coefficient
227	integer :: jj
228
229	real :: P_indice_factor_local
230
231	! output
232
233	real :: P_indice_factor
234
235	if (firstcall) then
236
237	a(1,1:4) = (/ 0.3567296E+01, 0.1508654E-01, &
238	0.4030668E-02, 0.7808922E-02 /)
239	a(2,1:4) = (/ -0.2775023E+04, 0.8828382E+01, &
240	0.1584634E+02,-0.1397511E+03 /)
241	a(3,1:4) = (/ -0.8403277E+02,-0.1444215E+01, &
242	-0.2192531E+01, 0.6054179E+00 /)
243	a(4,1:4) = (/ 0.1056939E-03, 0.1387796E-04, &
244	-0.1125676E-04,-0.1435086E-04 /)
245
246	b(1,1:4) = (/ 0.0 , 0.1383608E-01, &
247	0.6072980E-02,-0.1321784E-01 /)
248	b(2,1:4) = (/ 0.0 ,-0.1237310E+03, &
249	-0.9694563E+02, 0.1389812E+03 /)
250	b(3,1:4) = (/ 0.0 , 0.2649297E+00, &
251	-0.6347786E+00,-0.5396207E+00 /)
252	b(4,1:4) = (/ 0.0 ,-0.8090800E-05, &
253	-0.1232726E-04, 0.1383023E-04 /)
254
255	firstcall = .false.
256	endif
257
258	P_indice_factor_local = 0.
259	do jj=1,4
260	P_indice_factor_local = P_indice_factor_local + &
261	( a(indice,jj)cos((jj-1)sza_local*deg2rad)+ &
262	b(indice,jj)sin((jj-1)sza_local*deg2rad) )
263	enddo
264
265	P_indice_factor = P_indice_factor_local
266
267	return
268
269	end function P_indice_factor
270
271	!***********************************************************************
272	function temp_elect(zkm,tt,sza_local,origin)
273
274	! Computes the electronic temperature, either from Theis et al 1980
275	! model (origin=1) or Theis et al. 1984 (origin=2)
276	! or NEUTRAL (origin .ge. 2) <=== NOT USED
277
278	!***********************************************************************
279
280	! Arguments
281
282	real :: tt ! Temperature
283	real :: zkm ! Altitude in km
284	integer :: origin ! Theis et al 1980 model (origin=1) or Theis et al 1984 model (origin=2) or NEUTRAL (origin>2)
285	real :: sza_local ! solar zenith angle in degree
286
287	! local variables:
288
289	logical, save :: firstcall = .true.
290	real,parameter :: pi = 4.*atan(1.0)
291	real,parameter :: deg2rad = pi/180.
292	real :: temp_elect ! electronic temperatures
293
294	real :: z_incre, sza_incre
295	integer :: ii, iz1, iz2, jj, jsza1, jsza2
296	real :: dfx, dfy, dfxy, dT
297
298	! origin = 1
299	integer, parameter :: nsza_Theis = 13
300	integer, parameter :: nz_Theis = 10
301	real, save :: t_Theis(nsza_Theis,nz_Theis)
302	real, save :: z_Theis(nz_Theis), sza_Theis(nsza_Theis)
303
304	! origin = 2
305	real :: sza_security ! degree
306	real, parameter :: Altimini = 150. ! km
307	real :: log10_temp, factor_e
308
309
310	if (firstcall) then
311
312	! 2022/09/17: IPLS-Venus GCM ayant zmax < 300 km, on a mis les donnees pour z < 350 km
313
314	z_Theis(1:10)= (/ 160.,170.,180.,190.,200.,225.,250., &
315	275.,300.,350. /)
316	sza_Theis(1:13)= (/ 0,15,30,45,60,75,90,105,120,135,150,165,180 /)
317
318	t_Theis(1,1:10) = (/ 1136.,1514.,1893.,2255.,2589.,3275.,3757., &
319	4083.,4304.,4604. /)
320	t_Theis(2,1:10) = (/ 1157.,1530.,1900.,2253.,2577.,3240.,3708., &
321	4027.,4246.,4516. /)
322	t_Theis(3,1:10) = (/ 1210.,1567.,1915.,2243.,2542.,3150.,3581., &
323	3883.,4099.,4391. /)
324	t_Theis(4,1:10) = (/ 1262.,1598.,1921.,2221.,2491.,3040.,3431., &
325	3712.,3923.,4232. /)
326	t_Theis(5,1:10) = (/ 1279.,1598.,1902.,2181.,2433.,2491.,3305., &
327	3507.,3773.,4083. /)
328	t_Theis(6,1:10) = (/ 1249.,1560.,1856.,2129.,2375.,2871.,3226., &
329	3482.,3675.,3967. /)
330	t_Theis(7,1:10) = (/ 1197.,1505.,1801.,2076.,2324.,2827.,3184., &
331	3436.,3621.,3881. /)
332	t_Theis(8,1:10) = (/ 1163.,1467.,1760.,2034.,2283.,2790.,3149., &
333	3400.,3576.,3808. /)
334	t_Theis(9,1:10) = (/ 1180.,1472.,1753.,2015.,2245.,2743.,3092., &
335	3337.,3509.,3726. /)
336	t_Theis(10,1:10)= (/ 1259.,1528.,1783.,2020.,2235.,2679.,3004., &
337	3239.,3409.,3631. /)
338	t_Theis(11,1:10)= (/ 1383.,1618.,1838.,2041.,2225.,2609.,2902., &
339	3124.,3295.,3535. /)
340	t_Theis(12,1:10)= (/ 1502.,1703.,1890.,2062.,2220.,2555.,2820., &
341	3032.,3204.,3463. /)
342	t_Theis(13,1:10)= (/ 1552.,1739.,1912.,2071.,2218.,2534.,2789., &
343	2997.,3169.,3436. /)
344
345	if (origin.gt.2) then
346	write(,)'Error in function temp_elect:'
347	write(,)'Origin must be either 1 or 2'
348	write(,)'Using neutral temperature instead'
349	endif
350
351	firstcall = .false.
352	endif
353
354	if(origin.eq.1) then
355	if ( zkm .lt. 130. ) then
356	temp_elect = tt
357	else if(zkm .lt. 160.) then
358	if (sza_local .lt. 180.) then
359	do jj=nsza_Theis,2,-1
360	if ( sza_local .lt. sza_Theis(jj) ) then
361	jsza1 = jj - 1
362	jsza2 = jj
363	endif
364	enddo
365	dfx = (t_Theis(jsza2,1) - t_Theis(jsza1,1)) &
366	/ (cos(sza_Theis(jsza2)*deg2rad) &
367	-cos(sza_Theis(jsza1)*deg2rad))
368
369	dT = t_Theis(jsza1,1) + &
370	dfx(cos(sza_localdeg2rad) &
371	-cos(sza_Theis(jsza1)*deg2rad))
372
373	dfy = (tt-dT) / (z_Theis(1)-130.)
374
375	temp_elect = dT + dfy*(z_Theis(1)-zkm)
376
377	else
378	dfy = (tt-t_Theis(nsza_Theis,1)) / (z_Theis(1)-130.)
379
380	temp_elect = t_Theis(13,1) + dfy*(z_Theis(1)-zkm)
381	endif
382	else
383	if(zkm .lt. 350. .and. sza_local .lt. 180.) then
384	do ii=nz_Theis,2,-1
385	if ( zkm .lt. z_Theis(ii) ) then
386	iz1 = ii - 1
387	iz2 = ii
388	endif
389	enddo
390	do jj=nsza_Theis,2,-1
391	if ( sza_local .lt. sza_Theis(jj) ) then
392	jsza1 = jj - 1
393	jsza2 = jj
394	endif
395	enddo
396	dfx = t_Theis(jsza2,iz1) - t_Theis(jsza1,iz1)
397	dfy = t_Theis(jsza1,iz2) - t_Theis(jsza1,iz1)
398	dfxy= t_Theis(jsza1,iz1) + t_Theis(jsza2,iz2) &
399	-(t_Theis(jsza2,iz1) + t_Theis(jsza1,iz2))
400
401	z_incre =(zkm-z_Theis(iz1))/(z_Theis(iz2)-z_Theis(iz1))
402	sza_incre =(cos(sza_local*deg2rad) &
403	-cos(sza_Theis(jsza1)*deg2rad)) &
404	/(cos(sza_Theis(jsza2)*deg2rad) &
405	-cos(sza_Theis(jsza1)*deg2rad))
406
407	temp_elect = dfxsza_incre + dfyz_incre &
408	- dfxysza_increz_incre + t_Theis(jsza1,iz1)
409	else
410	if (sza_local .ge. 180.) then
411	if (zkm .lt. 350.) then
412	do ii=nz_Theis,2,-1
413	if ( zkm .lt. z_Theis(ii) ) then
414	iz1 = ii - 1
415	iz2 = ii
416	endif
417	enddo
418	dfy = (t_Theis(nsza_Theis,iz2)-t_Theis(nsza_Theis,iz1)) &
419	/ (z_Theis(iz2)-z_Theis(iz1))
420
421	temp_elect = t_Theis(nsza_Theis,iz1)+dfy*(zkm-z_Theis(iz1))
422	else
423	temp_elect = t_Theis(nsza_Theis,nz_Theis)
424	endif
425	endif
426	if (zkm .ge. 350.) then
427	do jj=nsza_Theis,2,-1
428	if ( sza_local .lt. sza_Theis(jj) ) then
429	jsza1 = jj - 1
430	jsza2 = jj
431	endif
432	enddo
433	dfx = (t_Theis(jsza2,nz_Theis) - t_Theis(jsza1,nz_Theis)) &
434	/ (cos(sza_Theis(jsza2)*deg2rad) &
435	-cos(sza_Theis(jsza1)*deg2rad))
436
437	temp_elect = t_Theis(jsza1,nz_Theis) + &
438	dfx(cos(sza_localdeg2rad) &
439	-cos(sza_Theis(jsza1)*deg2rad))
440	endif
441	endif
442	endif
443	else if(origin.eq.2) then
444	if (sza_local .gt. 180.) then
445	sza_security = 180.
446	else
447	sza_security = sza_local
448	endif
449
450	if (zkm .lt. 130.) then
451	temp_elect = tt
452	else if (zkm .lt. Altimini) then
453	log10_temp = P_indice_factor(sza_security,1) &
454	+ ( P_indice_factor(sza_security,2) &
455	/ ((Altimini + P_indice_factor(sza_security,3))**2.) ) &
456	+ Altimini * P_indice_factor(sza_security,4)
457
458	factor_e = exp((zkm-Altimini)/10.)
459
460	!temp_elect = 10.(log10_temp) + (tt-10.(log10_temp))* &
461	! (Altimini-zkm)/(Altimini-130.)
462	temp_elect = 10.*(log10_temp)factor_e + tt*(1.-factor_e)
463
464	else
465	log10_temp = P_indice_factor(sza_security,1) &
466	+ ( P_indice_factor(sza_security,2) &
467	/ ((zkm + P_indice_factor(sza_security,3))**2.) ) &
468	+ zkm * P_indice_factor(sza_security,4)
469
470	temp_elect = 10.**(log10_temp)
471	endif
472	else
473	!MAVEN measured electronic temperature (Ergun et al., GRL 2015)
474	!Note that the Langmuir probe is not sensitive below ~500K, so
475	!electronic temperatures in the lower thermosphere (<150 km) may
476	!be overestimated by this formula
477	!if(zkm.le.120) then
478	! temp_elect = tt
479	!else
480	! temp_elect=((3140.+120.)/2.)+((3140.-120.)/2.)*tanh((zkm-241.)/60.)
481	!endif
482	!else
483	!write(,)'Error in function temp_elect:'
484	!write(,)'Origin must be either 1 or 2'
485	!write(,)'Using neutral temperature instead'
486	temp_elect = tt
487	endif
488
489	return
490
491	end function temp_elect
492
493	!***********************************************************************
494	function temp_ion(zkm,tt,sza_local,origin)
495
496	! Computes the ion temperature, from VIRA Model based on the model
497	! of Miller et al. 1980 & 1984.
498	! or NEUTRAL (origin=2) <=== NOT USED
499
500	!***********************************************************************
501
502	! Arguments
503
504	real :: tt ! Temperature
505	real :: zkm ! Altitude in km
506	integer :: origin ! Miller et al 1984 model (origin=1) or NEUTRAL (origin=2)
507	real :: sza_local ! solar zenith angle in degree
508
509	! local variables:
510
511	logical, save :: firstcall = .true.
512	real,parameter :: pi = 4.*atan(1.0)
513	real,parameter :: deg2rad = pi/180.
514	real :: temp_ion ! electronic temperatures
515
516	real :: z_incre, sza_incre
517	integer :: ii, iz1, iz2, jj, jsza1, jsza2
518	real :: dfx, dfy, dfxy, dT
519
520	! origin = 1
521	real :: sza_security ! degree
522	real, parameter :: Altimini = 150. ! km - doit etre egal a z_Miller(1)
523	real :: log10_temp, factor_e
524	real :: zkm_security ! km
525	integer, parameter :: nsza_Miller = 9
526	integer, parameter :: nz_Miller = 11
527	real, save :: t_Miller(nsza_Miller,nz_Miller)
528	real, save :: z_Miller(nz_Miller), sza_Miller(nsza_Miller)
529
530
531	if (firstcall) then
532
533	! 2022/09/17: IPLS-Venus GCM ayant zmax < 300 km, on a mis les donnees pour z < 350 km
534
535	z_Miller(1:11) = (/ 150.,160.,170.,180.,190.,200.,225., &
536	250.,275.,300.,350. /)
537
538	sza_Miller(1:9) = (/ 20.,40.,60.,75.,85.,95.,110.,135.,160. /)
539
540	t_Miller(1,1:11)= (/ 300., 380., 470., 560., 580., 600., 660., &
541	930.,1400.,1900.,2700. /)
542	t_Miller(2,1:11)= (/ 310., 390., 490., 580., 640., 620., 660., &
543	940.,1500.,1900.,2400. /)
544	! Il y a une erreur dans la temperature ion dans CHAPTER IV de VIRA
545	! du coup, j'ai pris l'article MILLER et al. 1984 pour 65 SZA
546	t_Miller(3,1:11)= (/ 330., 410., 490., 540., 600., 610., 670., &
547	880.,1300.,1700.,2100. /)
548	t_Miller(4,1:11)= (/ 310., 390., 510., 580., 600., 600., 670., &
549	920.,1400.,1700.,2000. /)
550	t_Miller(5,1:11)= (/ 320., 430., 530., 600., 640., 660., 840., &
551	1100.,1300.,1500.,1800. /)
552	t_Miller(6,1:11)= (/ 290., 440., 530., 580., 640., 750.,1100., &
553	1300.,1400.,1500.,1700. /)
554	t_Miller(7,1:11)= (/ 230., 410., 600., 800.,1100.,1300.,1500., &
555	1700.,1700.,1800.,1800. /)
556	t_Miller(8,1:11)= (/ 210., 310., 540., 910.,1500.,2000.,2400., &
557	2700.,2600.,2500.,2400. /)
558	t_Miller(9,1:11)= (/ 230., 350., 500., 800.,1200.,1700.,2800., &
559	3500.,4100.,4700.,5500. /)
560
561	if (origin.gt.1) then
562	write(,)'Error in function temp_ion:'
563	write(,)'Origin must be either 1'
564	write(,)'Using neutral temperature instead'
565	endif
566
567	firstcall = .false.
568	endif
569
570	if(origin.eq.1) then
571	! Altitude security
572	if (zkm .gt. z_Miller(nz_Miller)) then
573	zkm_security = z_Miller(nz_Miller)
574	else
575	zkm_security = zkm
576	endif
577
578	if (zkm_security .lt. 130.) then
579	temp_ion = tt
580	else
581	! SZA security
582	if (sza_local .gt. sza_Miller(nsza_Miller)) then
583	sza_security = sza_Miller(nsza_Miller)
584	jsza1 = nsza_Miller - 1
585	jsza2 = nsza_Miller
586	else if (sza_local .lt. sza_Miller(1)) then
587	sza_security = sza_Miller(1)
588	jsza1 = 2 - 1
589	jsza2 = 2
590	else
591	sza_security = sza_local
592	do jj=nsza_Miller,2,-1
593	if ( sza_security .lt. sza_Miller(jj) ) then
594	jsza1 = jj - 1
595	jsza2 = jj
596	endif
597	enddo
598	endif
599
600
601	if(zkm_security .lt. Altimini) then
602	dfx = (t_Miller(jsza2,1) - t_Miller(jsza1,1)) &
603	/ (cos(sza_Miller(jsza2)*deg2rad) &
604	-cos(sza_Miller(jsza1)*deg2rad))
605
606	dT = t_Miller(jsza1,1) + &
607	dfx(cos(sza_securitydeg2rad) &
608	-cos(sza_Miller(jsza1)*deg2rad))
609
610	dfy = (tt-dT) / (z_Miller(1)-130.)
611
612	temp_ion = dT + dfy*(z_Miller(1)-zkm_security)
613	else
614	do ii=nz_Miller,2,-1
615	if ( zkm_security .lt. z_Miller(ii) ) then
616	iz1 = ii - 1
617	iz2 = ii
618	endif
619	enddo
620	dfx = t_Miller(jsza2,iz1) - t_Miller(jsza1,iz1)
621	dfy = t_Miller(jsza1,iz2) - t_Miller(jsza1,iz1)
622	dfxy= t_Miller(jsza1,iz1) + t_Miller(jsza2,iz2) &
623	-(t_Miller(jsza2,iz1) + t_Miller(jsza1,iz2))
624
625	z_incre =(zkm_security-z_Miller(iz1)) &
626	/(z_Miller(iz2)-z_Miller(iz1))
627	sza_incre =(cos(sza_security*deg2rad) &
628	-cos(sza_Miller(jsza1)*deg2rad)) &
629	/(cos(sza_Miller(jsza2)*deg2rad) &
630	-cos(sza_Miller(jsza1)*deg2rad))
631
632	temp_ion = dfxsza_incre + dfyz_incre &
633	- dfxysza_increz_incre + t_Miller(jsza1,iz1)
634	endif
635	endif
636	else
637	!MAVEN measured electronic temperature (Ergun et al., GRL 2015)
638	!Note that the Langmuir probe is not sensitive below ~500K, so
639	!electronic temperatures in the lower thermosphere (<150 km) may
640	!be overestimated by this formula
641	!if(zkm.le.120) then
642	! temp_elect = tt
643	!else
644	! temp_elect=((3140.+120.)/2.)+((3140.-120.)/2.)*tanh((zkm-241.)/60.)
645	!endif
646	!else
647	temp_ion = tt
648	endif
649
650	return
651
652	end function temp_ion
653
654
655	!***********************************************************************
656	function temp_mars_elect(zkm,tt,sza_local,origin)
657
658	! Computes the electronic temperature, either from Viking (origin=1)
659	! or MAVEN (origin=2)
660
661	!***********************************************************************
662
663	! Arguments
664
665	real tt ! Temperature
666	real zkm ! Altitude in km
667	integer origin ! Viking (origin=1) or MAVEN (origin=2)
668	real sza_local ! solar zenith angle
669
670	! local variables:
671	real temp_mars_elect ! electronic temperatures
672	real zhanson(9),tehanson(9)
673	real incremento
674	integer ii, i1, i2
675
676	zhanson(1:9) = (/ 120.,130.,150.,175.,200.,225.,250.,275.,300. /)
677	tehanson(2:9) = (/ 200.,300.,500.,1250.,2000.,2200.,2400.,2500. /)
678	tehanson(1) = tt
679
680	if(origin.eq.1) then
681	if ( zkm .le. 120. ) then
682	temp_mars_elect = tt
683	else if(zkm .ge.300.) then
684	temp_mars_elect=tehanson(9)
685	else
686	do ii=9,2,-1
687	if ( zkm .lt. zhanson(ii) ) then
688	i1 = ii - 1
689	i2 = ii
690	endif
691	enddo
692	incremento=(tehanson(i2)-tehanson(i1)) / &
693	(zhanson(i2)-zhanson(i1))
694	temp_mars_elect = tehanson(i1) + &
695	(zkm-zhanson(i1)) * incremento
696	endif
697	else if(origin.eq.2) then
698	!MAVEN measured electronic temperature (Ergun et al., GRL 2015)
699	!Note that the Langmuir probe is not sensitive below ~500K, so
700	!electronic temperatures in the lower thermosphere (<150 km) may
701	!be overestimated by this formula
702	if(zkm.le.120) then
703	temp_mars_elect = tt
704	else
705	temp_mars_elect=((3140.+120.)/2.) + &
706	((3140.-120.)/2.)*tanh((zkm-241.)/60.)
707	endif
708	else
709	write(,)'Error in function temp_elect:'
710	write(,)'Origin must be either 1 or 2'
711	write(,)'Using neutral temperature instead'
712	temp_mars_elect = tt
713	endif
714
715	return
716
717	end function temp_mars_elect
718
719	END MODULE iono_h

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format