Context Navigation

iono_h.F90 @ 3556

Last change on this file since 3556 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

Rev	Line
[2836]	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.
[2969]	224	real,parameter :: pi = 4.*atan(1.0)
	225	real,parameter :: deg2rad = pi/180.
[2836]	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 + &
[2969]	261	( a(indice,jj)cos((jj-1)sza_local*deg2rad)+ &
	262	b(indice,jj)sin((jj-1)sza_local*deg2rad) )
[2836]	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
[2969]	289	logical, save :: firstcall = .true.
	290	real,parameter :: pi = 4.*atan(1.0)
	291	real,parameter :: deg2rad = pi/180.
[2836]	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)) &
[2969]	366	/ (cos(sza_Theis(jsza2)*deg2rad) &
	367	-cos(sza_Theis(jsza1)*deg2rad))
[2836]	368
	369	dT = t_Theis(jsza1,1) + &
[2969]	370	dfx(cos(sza_localdeg2rad) &
	371	-cos(sza_Theis(jsza1)*deg2rad))
[2836]	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))
[2969]	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))
[2836]	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)) &
[2969]	434	/ (cos(sza_Theis(jsza2)*deg2rad) &
	435	-cos(sza_Theis(jsza1)*deg2rad))
[2836]	436
	437	temp_elect = t_Theis(jsza1,nz_Theis) + &
[2969]	438	dfx(cos(sza_localdeg2rad) &
	439	-cos(sza_Theis(jsza1)*deg2rad))
[2836]	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.
[2969]	512	real,parameter :: pi = 4.*atan(1.0)
	513	real,parameter :: deg2rad = pi/180.
[2836]	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)) &
[2969]	603	/ (cos(sza_Miller(jsza2)*deg2rad) &
	604	-cos(sza_Miller(jsza1)*deg2rad))
[2836]	605
	606	dT = t_Miller(jsza1,1) + &
[2969]	607	dfx(cos(sza_securitydeg2rad) &
	608	-cos(sza_Miller(jsza1)*deg2rad))
[2836]	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))
[2969]	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))
[2836]	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