Context Navigation

iono_h.F90 @ 2889

Last change on this file since 2889 was 2836, checked in by abierjon, 3 years ago

VENUS GCM:

INCLUDING THE IONOSPHERE CODE IN VENUS GCM

ATTENTION: INCREASED MEMORY DEMAND

NEEDS AT LEAST 135 GB of allocated memory

==============================================================
===> LIST OF APPENDED FILES AND INTERNAL ADDITIONS <==========
==============================================================

NEW MODEL SPECIES

deftank/traceur-chemistry-IONOSPHERE.def

Neutrals: 4 Species: N, NO, NO2, N(2D)

Ions: 15 species:

CO2+, CO+, O+, O2+, H+,

N2+, H2O+, OH+, C+, HCO+,

H3O+, HCO2+, N+, NO+, elec

FROM 36 to 55 chemical species

NEW KEYWORD OF PHYSIQ.DEF

deftank/physiq-96x96x90-chemistry-IONOSPHERE.def

ok_ionchem: keyword supposed to activate ion chemistry.

(be careful that n, no, no2, n2d and the ion species are in the deftank/traceur-chemistry-IONOSPHERE.def)

ok_jonline: keyword supposed to activate the online photochemistry

(be careful that n, no, no2 and n2d are in the traceur-chemistry-IONOSPHERE.def)

==============================================================
===> LIST OF MODIFIED PROGRAMS <=========================
==============================================================

nonoro_gwd_ran_mod.F90

Change EPFLUXMAX value from 5.E-3 to 1.E-3

photochemistry_venus.F90 (krates; photolysis_online; indices;)

Import of keywords: ok_ionchem, tuneupperatm & ok_jonline
addition of ion species in order
Forcing electroneutrality
Update of the reactions a001 and a002 with taking into account the other species
Change of formula for a002 with the formulation of Baulch et al., 1976 (confirmed by Smith and Robertson, 2008)

photolysis_mod.F90

Modification of the rdsolarflux subroutine to include interpolation with Atlas1 and Atlas3 in connection with E10.7

concentration2.F90

Added chemical species n2d, no, no2 and n in the conductivity calculation.

The ions have been excluded because their sum is 10^{5 times less dense than the neutrals and

their thermal conductivity is unknown}

iono_h.F90

Addition of the phdisrate routine
replace the electronic temperature of Mars by that of

origin = 1: Theis et al. 1980 (Venus) with bilinear interpolation altitude/cos(SZA)
origin = 2: Theis et al. 1984 (Venus) with the formula of the electronic temperature at high solar activity

addition of an ion temperature model based on VIRA

cleshphys.h

added ok_ionchem & ok_jonline in COMMON/clesphys_l/

conf_phys.f90

add ok_ionchem & ok_jonline as parameters to read from physiq.def file set to .false. by default

chemparam_mod.F90

Add species in M_tr and corresponding i_X. Set all i_X to zero before reading traceur-chemistry-IONOSPHERE.def
Added Type_tr table to differentiate species: 1 == neutral, 2 == ION, 3 == liquid, 10 == others

euvheat.F90; hrtherm.F; jthermcalc_e107.F; param_read_e107.F

Normalization with Mars

A.M

File size: 26.8 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, save :: a(4,4), b(4,4) ! Fourrier coefficient
225	integer :: jj
226
227	real :: P_indice_factor_local
228
229	! output
230
231	real :: P_indice_factor
232
233	if (firstcall) then
234
235	a(1,1:4) = (/ 0.3567296E+01, 0.1508654E-01, &
236	0.4030668E-02, 0.7808922E-02 /)
237	a(2,1:4) = (/ -0.2775023E+04, 0.8828382E+01, &
238	0.1584634E+02,-0.1397511E+03 /)
239	a(3,1:4) = (/ -0.8403277E+02,-0.1444215E+01, &
240	-0.2192531E+01, 0.6054179E+00 /)
241	a(4,1:4) = (/ 0.1056939E-03, 0.1387796E-04, &
242	-0.1125676E-04,-0.1435086E-04 /)
243
244	b(1,1:4) = (/ 0.0 , 0.1383608E-01, &
245	0.6072980E-02,-0.1321784E-01 /)
246	b(2,1:4) = (/ 0.0 ,-0.1237310E+03, &
247	-0.9694563E+02, 0.1389812E+03 /)
248	b(3,1:4) = (/ 0.0 , 0.2649297E+00, &
249	-0.6347786E+00,-0.5396207E+00 /)
250	b(4,1:4) = (/ 0.0 ,-0.8090800E-05, &
251	-0.1232726E-04, 0.1383023E-04 /)
252
253	firstcall = .false.
254	endif
255
256	P_indice_factor_local = 0.
257	do jj=1,4
258	P_indice_factor_local = P_indice_factor_local + &
259	( a(indice,jj)cosd((jj-1)sza_local) + &
260	b(indice,jj)sind((jj-1)sza_local) )
261	enddo
262
263	P_indice_factor = P_indice_factor_local
264
265	return
266
267	end function P_indice_factor
268
269	!***********************************************************************
270	function temp_elect(zkm,tt,sza_local,origin)
271
272	! Computes the electronic temperature, either from Theis et al 1980
273	! model (origin=1) or Theis et al. 1984 (origin=2)
274	! or NEUTRAL (origin .ge. 2) <=== NOT USED
275
276	!***********************************************************************
277
278	! Arguments
279
280	real :: tt ! Temperature
281	real :: zkm ! Altitude in km
282	integer :: origin ! Theis et al 1980 model (origin=1) or Theis et al 1984 model (origin=2) or NEUTRAL (origin>2)
283	real :: sza_local ! solar zenith angle in degree
284
285	! local variables:
286
287	logical, save :: firstcall = .true.
288	real :: temp_elect ! electronic temperatures
289
290	real :: z_incre, sza_incre
291	integer :: ii, iz1, iz2, jj, jsza1, jsza2
292	real :: dfx, dfy, dfxy, dT
293
294	! origin = 1
295	integer, parameter :: nsza_Theis = 13
296	integer, parameter :: nz_Theis = 10
297	real, save :: t_Theis(nsza_Theis,nz_Theis)
298	real, save :: z_Theis(nz_Theis), sza_Theis(nsza_Theis)
299
300	! origin = 2
301	real :: sza_security ! degree
302	real, parameter :: Altimini = 150. ! km
303	real :: log10_temp, factor_e
304
305
306	if (firstcall) then
307
308	! 2022/09/17: IPLS-Venus GCM ayant zmax < 300 km, on a mis les donnees pour z < 350 km
309
310	z_Theis(1:10)= (/ 160.,170.,180.,190.,200.,225.,250., &
311	275.,300.,350. /)
312	sza_Theis(1:13)= (/ 0,15,30,45,60,75,90,105,120,135,150,165,180 /)
313
314	t_Theis(1,1:10) = (/ 1136.,1514.,1893.,2255.,2589.,3275.,3757., &
315	4083.,4304.,4604. /)
316	t_Theis(2,1:10) = (/ 1157.,1530.,1900.,2253.,2577.,3240.,3708., &
317	4027.,4246.,4516. /)
318	t_Theis(3,1:10) = (/ 1210.,1567.,1915.,2243.,2542.,3150.,3581., &
319	3883.,4099.,4391. /)
320	t_Theis(4,1:10) = (/ 1262.,1598.,1921.,2221.,2491.,3040.,3431., &
321	3712.,3923.,4232. /)
322	t_Theis(5,1:10) = (/ 1279.,1598.,1902.,2181.,2433.,2491.,3305., &
323	3507.,3773.,4083. /)
324	t_Theis(6,1:10) = (/ 1249.,1560.,1856.,2129.,2375.,2871.,3226., &
325	3482.,3675.,3967. /)
326	t_Theis(7,1:10) = (/ 1197.,1505.,1801.,2076.,2324.,2827.,3184., &
327	3436.,3621.,3881. /)
328	t_Theis(8,1:10) = (/ 1163.,1467.,1760.,2034.,2283.,2790.,3149., &
329	3400.,3576.,3808. /)
330	t_Theis(9,1:10) = (/ 1180.,1472.,1753.,2015.,2245.,2743.,3092., &
331	3337.,3509.,3726. /)
332	t_Theis(10,1:10)= (/ 1259.,1528.,1783.,2020.,2235.,2679.,3004., &
333	3239.,3409.,3631. /)
334	t_Theis(11,1:10)= (/ 1383.,1618.,1838.,2041.,2225.,2609.,2902., &
335	3124.,3295.,3535. /)
336	t_Theis(12,1:10)= (/ 1502.,1703.,1890.,2062.,2220.,2555.,2820., &
337	3032.,3204.,3463. /)
338	t_Theis(13,1:10)= (/ 1552.,1739.,1912.,2071.,2218.,2534.,2789., &
339	2997.,3169.,3436. /)
340
341	if (origin.gt.2) then
342	write(,)'Error in function temp_elect:'
343	write(,)'Origin must be either 1 or 2'
344	write(,)'Using neutral temperature instead'
345	endif
346
347	firstcall = .false.
348	endif
349
350	if(origin.eq.1) then
351	if ( zkm .lt. 130. ) then
352	temp_elect = tt
353	else if(zkm .lt. 160.) then
354	if (sza_local .lt. 180.) then
355	do jj=nsza_Theis,2,-1
356	if ( sza_local .lt. sza_Theis(jj) ) then
357	jsza1 = jj - 1
358	jsza2 = jj
359	endif
360	enddo
361	dfx = (t_Theis(jsza2,1) - t_Theis(jsza1,1)) &
362	/ (cosd(sza_Theis(jsza2))-cosd(sza_Theis(jsza1)))
363
364	dT = t_Theis(jsza1,1) + &
365	dfx*(cosd(sza_local)-cosd(sza_Theis(jsza1)))
366
367	dfy = (tt-dT) / (z_Theis(1)-130.)
368
369	temp_elect = dT + dfy*(z_Theis(1)-zkm)
370
371	else
372	dfy = (tt-t_Theis(nsza_Theis,1)) / (z_Theis(1)-130.)
373
374	temp_elect = t_Theis(13,1) + dfy*(z_Theis(1)-zkm)
375	endif
376	else
377	if(zkm .lt. 350. .and. sza_local .lt. 180.) then
378	do ii=nz_Theis,2,-1
379	if ( zkm .lt. z_Theis(ii) ) then
380	iz1 = ii - 1
381	iz2 = ii
382	endif
383	enddo
384	do jj=nsza_Theis,2,-1
385	if ( sza_local .lt. sza_Theis(jj) ) then
386	jsza1 = jj - 1
387	jsza2 = jj
388	endif
389	enddo
390	dfx = t_Theis(jsza2,iz1) - t_Theis(jsza1,iz1)
391	dfy = t_Theis(jsza1,iz2) - t_Theis(jsza1,iz1)
392	dfxy= t_Theis(jsza1,iz1) + t_Theis(jsza2,iz2) &
393	-(t_Theis(jsza2,iz1) + t_Theis(jsza1,iz2))
394
395	z_incre =(zkm-z_Theis(iz1))/(z_Theis(iz2)-z_Theis(iz1))
396	sza_incre =(cosd(sza_local)-cosd(sza_Theis(jsza1))) &
397	/(cosd(sza_Theis(jsza2))-cosd(sza_Theis(jsza1)))
398
399	temp_elect = dfxsza_incre + dfyz_incre &
400	- dfxysza_increz_incre + t_Theis(jsza1,iz1)
401	else
402	if (sza_local .ge. 180.) then
403	if (zkm .lt. 350.) then
404	do ii=nz_Theis,2,-1
405	if ( zkm .lt. z_Theis(ii) ) then
406	iz1 = ii - 1
407	iz2 = ii
408	endif
409	enddo
410	dfy = (t_Theis(nsza_Theis,iz2)-t_Theis(nsza_Theis,iz1)) &
411	/ (z_Theis(iz2)-z_Theis(iz1))
412
413	temp_elect = t_Theis(nsza_Theis,iz1)+dfy*(zkm-z_Theis(iz1))
414	else
415	temp_elect = t_Theis(nsza_Theis,nz_Theis)
416	endif
417	endif
418	if (zkm .ge. 350.) then
419	do jj=nsza_Theis,2,-1
420	if ( sza_local .lt. sza_Theis(jj) ) then
421	jsza1 = jj - 1
422	jsza2 = jj
423	endif
424	enddo
425	dfx = (t_Theis(jsza2,nz_Theis) - t_Theis(jsza1,nz_Theis)) &
426	/ (cosd(sza_Theis(jsza2))-cosd(sza_Theis(jsza1)))
427
428	temp_elect = t_Theis(jsza1,nz_Theis) + &
429	dfx*(cosd(sza_local)-cosd(sza_Theis(jsza1)))
430	endif
431	endif
432	endif
433	else if(origin.eq.2) then
434	if (sza_local .gt. 180.) then
435	sza_security = 180.
436	else
437	sza_security = sza_local
438	endif
439
440	if (zkm .lt. 130.) then
441	temp_elect = tt
442	else if (zkm .lt. Altimini) then
443	log10_temp = P_indice_factor(sza_security,1) &
444	+ ( P_indice_factor(sza_security,2) &
445	/ ((Altimini + P_indice_factor(sza_security,3))**2.) ) &
446	+ Altimini * P_indice_factor(sza_security,4)
447
448	factor_e = exp((zkm-Altimini)/10.)
449
450	!temp_elect = 10.(log10_temp) + (tt-10.(log10_temp))* &
451	! (Altimini-zkm)/(Altimini-130.)
452	temp_elect = 10.*(log10_temp)factor_e + tt*(1.-factor_e)
453
454	else
455	log10_temp = P_indice_factor(sza_security,1) &
456	+ ( P_indice_factor(sza_security,2) &
457	/ ((zkm + P_indice_factor(sza_security,3))**2.) ) &
458	+ zkm * P_indice_factor(sza_security,4)
459
460	temp_elect = 10.**(log10_temp)
461	endif
462	else
463	!MAVEN measured electronic temperature (Ergun et al., GRL 2015)
464	!Note that the Langmuir probe is not sensitive below ~500K, so
465	!electronic temperatures in the lower thermosphere (<150 km) may
466	!be overestimated by this formula
467	!if(zkm.le.120) then
468	! temp_elect = tt
469	!else
470	! temp_elect=((3140.+120.)/2.)+((3140.-120.)/2.)*tanh((zkm-241.)/60.)
471	!endif
472	!else
473	!write(,)'Error in function temp_elect:'
474	!write(,)'Origin must be either 1 or 2'
475	!write(,)'Using neutral temperature instead'
476	temp_elect = tt
477	endif
478
479	return
480
481	end function temp_elect
482
483	!***********************************************************************
484	function temp_ion(zkm,tt,sza_local,origin)
485
486	! Computes the ion temperature, from VIRA Model based on the model
487	! of Miller et al. 1980 & 1984.
488	! or NEUTRAL (origin=2) <=== NOT USED
489
490	!***********************************************************************
491
492	! Arguments
493
494	real :: tt ! Temperature
495	real :: zkm ! Altitude in km
496	integer :: origin ! Miller et al 1984 model (origin=1) or NEUTRAL (origin=2)
497	real :: sza_local ! solar zenith angle in degree
498
499	! local variables:
500
501	logical, save :: firstcall = .true.
502	real :: temp_ion ! electronic temperatures
503
504	real :: z_incre, sza_incre
505	integer :: ii, iz1, iz2, jj, jsza1, jsza2
506	real :: dfx, dfy, dfxy, dT
507
508	! origin = 1
509	real :: sza_security ! degree
510	real, parameter :: Altimini = 150. ! km - doit etre egal a z_Miller(1)
511	real :: log10_temp, factor_e
512	real :: zkm_security ! km
513	integer, parameter :: nsza_Miller = 9
514	integer, parameter :: nz_Miller = 11
515	real, save :: t_Miller(nsza_Miller,nz_Miller)
516	real, save :: z_Miller(nz_Miller), sza_Miller(nsza_Miller)
517
518
519	if (firstcall) then
520
521	! 2022/09/17: IPLS-Venus GCM ayant zmax < 300 km, on a mis les donnees pour z < 350 km
522
523	z_Miller(1:11) = (/ 150.,160.,170.,180.,190.,200.,225., &
524	250.,275.,300.,350. /)
525
526	sza_Miller(1:9) = (/ 20.,40.,60.,75.,85.,95.,110.,135.,160. /)
527
528	t_Miller(1,1:11)= (/ 300., 380., 470., 560., 580., 600., 660., &
529	930.,1400.,1900.,2700. /)
530	t_Miller(2,1:11)= (/ 310., 390., 490., 580., 640., 620., 660., &
531	940.,1500.,1900.,2400. /)
532	! Il y a une erreur dans la temperature ion dans CHAPTER IV de VIRA
533	! du coup, j'ai pris l'article MILLER et al. 1984 pour 65 SZA
534	t_Miller(3,1:11)= (/ 330., 410., 490., 540., 600., 610., 670., &
535	880.,1300.,1700.,2100. /)
536	t_Miller(4,1:11)= (/ 310., 390., 510., 580., 600., 600., 670., &
537	920.,1400.,1700.,2000. /)
538	t_Miller(5,1:11)= (/ 320., 430., 530., 600., 640., 660., 840., &
539	1100.,1300.,1500.,1800. /)
540	t_Miller(6,1:11)= (/ 290., 440., 530., 580., 640., 750.,1100., &
541	1300.,1400.,1500.,1700. /)
542	t_Miller(7,1:11)= (/ 230., 410., 600., 800.,1100.,1300.,1500., &
543	1700.,1700.,1800.,1800. /)
544	t_Miller(8,1:11)= (/ 210., 310., 540., 910.,1500.,2000.,2400., &
545	2700.,2600.,2500.,2400. /)
546	t_Miller(9,1:11)= (/ 230., 350., 500., 800.,1200.,1700.,2800., &
547	3500.,4100.,4700.,5500. /)
548
549	if (origin.gt.1) then
550	write(,)'Error in function temp_ion:'
551	write(,)'Origin must be either 1'
552	write(,)'Using neutral temperature instead'
553	endif
554
555	firstcall = .false.
556	endif
557
558	if(origin.eq.1) then
559	! Altitude security
560	if (zkm .gt. z_Miller(nz_Miller)) then
561	zkm_security = z_Miller(nz_Miller)
562	else
563	zkm_security = zkm
564	endif
565
566	if (zkm_security .lt. 130.) then
567	temp_ion = tt
568	else
569	! SZA security
570	if (sza_local .gt. sza_Miller(nsza_Miller)) then
571	sza_security = sza_Miller(nsza_Miller)
572	jsza1 = nsza_Miller - 1
573	jsza2 = nsza_Miller
574	else if (sza_local .lt. sza_Miller(1)) then
575	sza_security = sza_Miller(1)
576	jsza1 = 2 - 1
577	jsza2 = 2
578	else
579	sza_security = sza_local
580	do jj=nsza_Miller,2,-1
581	if ( sza_security .lt. sza_Miller(jj) ) then
582	jsza1 = jj - 1
583	jsza2 = jj
584	endif
585	enddo
586	endif
587
588
589	if(zkm_security .lt. Altimini) then
590	dfx = (t_Miller(jsza2,1) - t_Miller(jsza1,1)) &
591	/ (cosd(sza_Miller(jsza2))-cosd(sza_Miller(jsza1)))
592
593	dT = t_Miller(jsza1,1) + &
594	dfx*(cosd(sza_security)-cosd(sza_Miller(jsza1)))
595
596	dfy = (tt-dT) / (z_Miller(1)-130.)
597
598	temp_ion = dT + dfy*(z_Miller(1)-zkm_security)
599	else
600	do ii=nz_Miller,2,-1
601	if ( zkm_security .lt. z_Miller(ii) ) then
602	iz1 = ii - 1
603	iz2 = ii
604	endif
605	enddo
606	dfx = t_Miller(jsza2,iz1) - t_Miller(jsza1,iz1)
607	dfy = t_Miller(jsza1,iz2) - t_Miller(jsza1,iz1)
608	dfxy= t_Miller(jsza1,iz1) + t_Miller(jsza2,iz2) &
609	-(t_Miller(jsza2,iz1) + t_Miller(jsza1,iz2))
610
611	z_incre =(zkm_security-z_Miller(iz1)) &
612	/(z_Miller(iz2)-z_Miller(iz1))
613	sza_incre =(cosd(sza_security)-cosd(sza_Miller(jsza1))) &
614	/(cosd(sza_Miller(jsza2))-cosd(sza_Miller(jsza1)))
615
616	temp_ion = dfxsza_incre + dfyz_incre &
617	- dfxysza_increz_incre + t_Miller(jsza1,iz1)
618	endif
619	endif
620	else
621	!MAVEN measured electronic temperature (Ergun et al., GRL 2015)
622	!Note that the Langmuir probe is not sensitive below ~500K, so
623	!electronic temperatures in the lower thermosphere (<150 km) may
624	!be overestimated by this formula
625	!if(zkm.le.120) then
626	! temp_elect = tt
627	!else
628	! temp_elect=((3140.+120.)/2.)+((3140.-120.)/2.)*tanh((zkm-241.)/60.)
629	!endif
630	!else
631	temp_ion = tt
632	endif
633
634	return
635
636	end function temp_ion
637
638
639	!***********************************************************************
640	function temp_mars_elect(zkm,tt,sza_local,origin)
641
642	! Computes the electronic temperature, either from Viking (origin=1)
643	! or MAVEN (origin=2)
644
645	!***********************************************************************
646
647	! Arguments
648
649	real tt ! Temperature
650	real zkm ! Altitude in km
651	integer origin ! Viking (origin=1) or MAVEN (origin=2)
652	real sza_local ! solar zenith angle
653
654	! local variables:
655	real temp_mars_elect ! electronic temperatures
656	real zhanson(9),tehanson(9)
657	real incremento
658	integer ii, i1, i2
659
660	zhanson(1:9) = (/ 120.,130.,150.,175.,200.,225.,250.,275.,300. /)
661	tehanson(2:9) = (/ 200.,300.,500.,1250.,2000.,2200.,2400.,2500. /)
662	tehanson(1) = tt
663
664	if(origin.eq.1) then
665	if ( zkm .le. 120. ) then
666	temp_mars_elect = tt
667	else if(zkm .ge.300.) then
668	temp_mars_elect=tehanson(9)
669	else
670	do ii=9,2,-1
671	if ( zkm .lt. zhanson(ii) ) then
672	i1 = ii - 1
673	i2 = ii
674	endif
675	enddo
676	incremento=(tehanson(i2)-tehanson(i1)) / &
677	(zhanson(i2)-zhanson(i1))
678	temp_mars_elect = tehanson(i1) + &
679	(zkm-zhanson(i1)) * incremento
680	endif
681	else if(origin.eq.2) then
682	!MAVEN measured electronic temperature (Ergun et al., GRL 2015)
683	!Note that the Langmuir probe is not sensitive below ~500K, so
684	!electronic temperatures in the lower thermosphere (<150 km) may
685	!be overestimated by this formula
686	if(zkm.le.120) then
687	temp_mars_elect = tt
688	else
689	temp_mars_elect=((3140.+120.)/2.) + &
690	((3140.-120.)/2.)*tanh((zkm-241.)/60.)
691	endif
692	else
693	write(,)'Error in function temp_elect:'
694	write(,)'Origin must be either 1 or 2'
695	write(,)'Using neutral temperature instead'
696	temp_mars_elect = tt
697	endif
698
699	return
700
701	end function temp_mars_elect
702
703	END MODULE iono_h

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

Download in other formats:

Original Format