MODULE iono_h IMPLICIT NONE CONTAINS subroutine phdisrate(ig,nlayer,chemthermod,zenit,i) ! Calculates photoionization and photodissociation rates from the ! photoabsorption rates calculated in jthermcalc_e107 and the ! ionization/dissociation branching ratios in param_read_e107 ! apr 2002 fgg first version !********************************************************************** use param_v4_h, only: ninter,nabs, & jfotsout,fluxtop, & jion,jdistot,jdistot_b, & efdisco2,efdiso2,efdish2o, & efdish2o2,efdish2,efdiso3, & efdiso,efdisn,efdish, & efdisno,efdisn2,efdisno2, & efdisco,efionco2,efiono2,efionn2, & efionco,efiono3p,efionn, & efionno,efionh ! & implicit none ! arguments integer i !altitude integer ig,chemthermod,nlayer real zenit ! local variables integer inter,iz,j real lambda real jdis(nabs,ninter,nlayer) character*1 dn !********************************************************************** ! photodissociation and photoionization rates initialization jion(:,i,:) = 0.d0 jdistot(:,i) = 0.d0 jdistot_b(:,i) = 0.d0 ! jion(1,i,1) = 0.d0 ! CO2 channel 1 ( --> CO2+ + e- ) ! jion(1,i,2) = 0.d0 ! CO2 channel 2 ( --> O+ + CO + e- ) ! jion(1,i,3) = 0.d0 ! CO2 channel 3 ( --> CO+ + O + e- ) ! jion(1,i,4) = 0.d0 ! CO2 channel 4 ( --> C+ + O2 + e- ) ! jion(2,i,1) = 0.d0 ! O2 (only one ionization channel) ! jion(3,i,1) = 0.d0 ! O3P (only one ionization channel) ! jion(4,i,1) = 0.d0 ! H2O (no ionization) ! jion(5,i,1) = 0.d0 ! H2 (no ionization) ! jion(6,i,1) = 0.d0 ! H2O2 (no ionization) ! jion(7,i,1) = 0.d0 ! O3 (no ionization) ! jion(8,i,1) = 0.d0 ! N2 channel 1 ( --> n2+ + e- ) ! jion(8,i,2) = 0.d0 ! N2 channel 2 ( --> n+ + n + e- ) ! jion(9,i,1) = 0.d0 ! N ( --> N+ + e- ) ! jion(10,i,1)= 0.d0 ! We do not mind its ionization ! jion(11,i,1) = 0.d0 ! CO channel 1 ( --> co+ + e- ) ! jion(11,i,2) = 0.d0 ! CO channel 2 ( --> c+ + o + e- ) ! jion(11,i,3) = 0.d0 ! CO channel 3 ( --> o+ + c + e- ) ! jion(12,i,1) = 0.d0 ! H ( --> H+ + e- ) ! jion(13,i,1) = 0.d0 ! do j=1,nabs ! jdistot(j,i) = 0. ! jdistot_b(j,i) = 0. ! end do if(zenit.gt.140.) then dn='n' else dn='d' end if if(dn.eq.'n') return ! photodissociation and photoionization rates for each species do inter=1,ninter ! CO2 jdis(1,inter,i) = jfotsout(inter,1,i) * fluxtop(inter) & * efdisco2(inter) if(inter.gt.29.and.inter.le.32) then jdistot(1,i) = jdistot(1,i) + jdis(1,inter,i) else if(inter.le.29) then jdistot_b(1,i) = jdistot_b(1,i) + jdis(1,inter,i) end if jion(1,i,1)=jion(1,i,1) + & jfotsout(inter,1,i)*fluxtop(inter)*efionco2(inter,1) jion(1,i,2)=jion(1,i,2) + & jfotsout(inter,1,i)*fluxtop(inter)*efionco2(inter,2) jion(1,i,3)=jion(1,i,3) + & jfotsout(inter,1,i)*fluxtop(inter)*efionco2(inter,3) jion(1,i,4)=jion(1,i,4) + & jfotsout(inter,1,i)*fluxtop(inter)*efionco2(inter,4) ! O2 jdis(2,inter,i) = jfotsout(inter,2,i) * fluxtop(inter) & * efdiso2(inter) if(inter.ge.31) then jdistot(2,i) = jdistot(2,i) + jdis(2,inter,i) else if(inter.eq.30) then jdistot(2,i)=jdistot(2,i)+0.02*jdis(2,inter,i) jdistot_b(2,i)=jdistot_b(2,i)+0.98*jdis(2,inter,i) else if(inter.lt.31) then jdistot_b(2,i) = jdistot_b(2,i) + jdis(2,inter,i) end if jion(2,i,1)=jion(2,i,1) + & jfotsout(inter,2,i) * fluxtop(inter) * efiono2(inter,1) jion(2,1,2)=jion(2,1,2) + & jfotsout(inter,2,i) * fluxtop(inter) * efiono2(inter,2) !(1.-efdiso2(inter)) ! O3P jion(3,i,1)=jion(3,i,1) + & jfotsout(inter,3,i) * fluxtop(inter) * efiono3p(inter) ! H2O jdis(4,inter,i) = jfotsout(inter,4,i) * fluxtop(inter) & * efdish2o(inter) jdistot(4,i) = jdistot(4,i) + jdis(4,inter,i) ! H2 jdis(5,inter,i) = jfotsout(inter,5,i) * fluxtop(inter) & * efdish2(inter) jdistot(5,i) = jdistot(5,i) + jdis(5,inter,i) ! H2O2 jdis(6,inter,i) = jfotsout(inter,6,i) * fluxtop(inter) & * efdish2o2(inter) jdistot(6,i) = jdistot(6,i) + jdis(6,inter,i) !Only if O3, N or ion chemistry requested if(chemthermod.ge.1) then ! O3 jdis(7,inter,i) = jfotsout(inter,7,i) * fluxtop(inter) & * efdiso3(inter) if(inter.eq.34) then jdistot(7,i) = jdistot(7,i) + jdis(7,inter,i) else if (inter.eq.35) then jdistot(7,i) = jdistot(7,i) + 0.997 * jdis(7,inter,i) jdistot_b(7,i) = jdistot_b(7,i) + 0.003 * jdis(7,inter,i) else if (inter.eq.36) then jdistot_b(7,i) = jdistot_b(7,i) + jdis(7,inter,i) endif endif !Of chemthermod.ge.1 !Only if N or ion chemistry requested if(chemthermod.ge.2) then ! N2 jdis(8,inter,i) = jfotsout(inter,8,i) * fluxtop(inter) & * efdisn2(inter) jdistot(8,i) = jdistot(8,i) + jdis(8,inter,i) jion(8,i,1) = jion(8,i,1) + jfotsout(inter,8,i) * & fluxtop(inter) * efionn2(inter,1) jion(8,i,2) = jion(8,i,2) + jfotsout(inter,8,i) * & fluxtop(inter) * efionn2(inter,2) ! N jion(9,i,1) = jion(9,i,1) + jfotsout(inter,9,i) * & fluxtop(inter) * efionn(inter) ! NO jdis(10,inter,i) = jfotsout(inter,10,i) * fluxtop(inter) & * efdisno(inter) jdistot(10,i) = jdistot(10,i) + jdis(10,inter,i) jion(10,i,1) = jion(10,i,1) + jfotsout(inter,10,i) * & fluxtop(inter) * efionno(inter) ! NO2 jdis(13,inter,i) = jfotsout(inter,13,i) * fluxtop(inter) & * efdisno2(inter) jdistot(13,i) = jdistot(13,i) + jdis(13,inter,i) endif !Of chemthermod.ge.2 ! CO jdis(11,inter,i) = jfotsout(inter,11,i) * fluxtop(inter) & * efdisco(inter) jdistot(11,i) = jdistot(11,i) + jdis(11,inter,i) jion(11,i,1) = jion(11,i,1) + jfotsout(inter,11,i) * & fluxtop(inter) * efionco(inter,1) jion(11,i,2) = jion(11,i,2) + jfotsout(inter,11,i) * & fluxtop(inter) * efionco(inter,2) jion(11,i,3) = jion(11,i,3) + jfotsout(inter,11,i) * & fluxtop(inter) * efionco(inter,3) ! H jion(12,i,1) = jion(12,i,1) + jfotsout(inter,12,i) * & fluxtop(inter) * efionh(inter) end do return end !*********************************************************************** function P_indice_factor(sza_local,indice) ! Computes the P_indice_factor for the electronic temperature of ! Theis et al. 1984 model !*********************************************************************** ! Arguments integer :: indice ! the indiceth factor real :: sza_local ! solar zenith angle in degree ! local variables: logical, save :: firstcall = .true. real,parameter :: pi = 4.*atan(1.0) real,parameter :: deg2rad = pi/180. real, save :: a(4,4), b(4,4) ! Fourrier coefficient integer :: jj real :: P_indice_factor_local ! output real :: P_indice_factor if (firstcall) then a(1,1:4) = (/ 0.3567296E+01, 0.1508654E-01, & 0.4030668E-02, 0.7808922E-02 /) a(2,1:4) = (/ -0.2775023E+04, 0.8828382E+01, & 0.1584634E+02,-0.1397511E+03 /) a(3,1:4) = (/ -0.8403277E+02,-0.1444215E+01, & -0.2192531E+01, 0.6054179E+00 /) a(4,1:4) = (/ 0.1056939E-03, 0.1387796E-04, & -0.1125676E-04,-0.1435086E-04 /) b(1,1:4) = (/ 0.0 , 0.1383608E-01, & 0.6072980E-02,-0.1321784E-01 /) b(2,1:4) = (/ 0.0 ,-0.1237310E+03, & -0.9694563E+02, 0.1389812E+03 /) b(3,1:4) = (/ 0.0 , 0.2649297E+00, & -0.6347786E+00,-0.5396207E+00 /) b(4,1:4) = (/ 0.0 ,-0.8090800E-05, & -0.1232726E-04, 0.1383023E-04 /) firstcall = .false. endif P_indice_factor_local = 0. do jj=1,4 P_indice_factor_local = P_indice_factor_local + & ( a(indice,jj)*cos((jj-1)*sza_local*deg2rad)+ & b(indice,jj)*sin((jj-1)*sza_local*deg2rad) ) enddo P_indice_factor = P_indice_factor_local return end function P_indice_factor !*********************************************************************** function temp_elect(zkm,tt,sza_local,origin) ! Computes the electronic temperature, either from Theis et al 1980 ! model (origin=1) or Theis et al. 1984 (origin=2) ! or NEUTRAL (origin .ge. 2) <=== NOT USED !*********************************************************************** ! Arguments real :: tt ! Temperature real :: zkm ! Altitude in km integer :: origin ! Theis et al 1980 model (origin=1) or Theis et al 1984 model (origin=2) or NEUTRAL (origin>2) real :: sza_local ! solar zenith angle in degree ! local variables: logical, save :: firstcall = .true. real,parameter :: pi = 4.*atan(1.0) real,parameter :: deg2rad = pi/180. real :: temp_elect ! electronic temperatures real :: z_incre, sza_incre integer :: ii, iz1, iz2, jj, jsza1, jsza2 real :: dfx, dfy, dfxy, dT ! origin = 1 integer, parameter :: nsza_Theis = 13 integer, parameter :: nz_Theis = 10 real, save :: t_Theis(nsza_Theis,nz_Theis) real, save :: z_Theis(nz_Theis), sza_Theis(nsza_Theis) ! origin = 2 real :: sza_security ! degree real, parameter :: Altimini = 150. ! km real :: log10_temp, factor_e if (firstcall) then ! 2022/09/17: IPLS-Venus GCM ayant zmax < 300 km, on a mis les donnees pour z < 350 km z_Theis(1:10)= (/ 160.,170.,180.,190.,200.,225.,250., & 275.,300.,350. /) sza_Theis(1:13)= (/ 0,15,30,45,60,75,90,105,120,135,150,165,180 /) t_Theis(1,1:10) = (/ 1136.,1514.,1893.,2255.,2589.,3275.,3757., & 4083.,4304.,4604. /) t_Theis(2,1:10) = (/ 1157.,1530.,1900.,2253.,2577.,3240.,3708., & 4027.,4246.,4516. /) t_Theis(3,1:10) = (/ 1210.,1567.,1915.,2243.,2542.,3150.,3581., & 3883.,4099.,4391. /) t_Theis(4,1:10) = (/ 1262.,1598.,1921.,2221.,2491.,3040.,3431., & 3712.,3923.,4232. /) t_Theis(5,1:10) = (/ 1279.,1598.,1902.,2181.,2433.,2491.,3305., & 3507.,3773.,4083. /) t_Theis(6,1:10) = (/ 1249.,1560.,1856.,2129.,2375.,2871.,3226., & 3482.,3675.,3967. /) t_Theis(7,1:10) = (/ 1197.,1505.,1801.,2076.,2324.,2827.,3184., & 3436.,3621.,3881. /) t_Theis(8,1:10) = (/ 1163.,1467.,1760.,2034.,2283.,2790.,3149., & 3400.,3576.,3808. /) t_Theis(9,1:10) = (/ 1180.,1472.,1753.,2015.,2245.,2743.,3092., & 3337.,3509.,3726. /) t_Theis(10,1:10)= (/ 1259.,1528.,1783.,2020.,2235.,2679.,3004., & 3239.,3409.,3631. /) t_Theis(11,1:10)= (/ 1383.,1618.,1838.,2041.,2225.,2609.,2902., & 3124.,3295.,3535. /) t_Theis(12,1:10)= (/ 1502.,1703.,1890.,2062.,2220.,2555.,2820., & 3032.,3204.,3463. /) t_Theis(13,1:10)= (/ 1552.,1739.,1912.,2071.,2218.,2534.,2789., & 2997.,3169.,3436. /) if (origin.gt.2) then write(*,*)'Error in function temp_elect:' write(*,*)'Origin must be either 1 or 2' write(*,*)'Using neutral temperature instead' endif firstcall = .false. endif if(origin.eq.1) then if ( zkm .lt. 130. ) then temp_elect = tt else if(zkm .lt. 160.) then if (sza_local .lt. 180.) then do jj=nsza_Theis,2,-1 if ( sza_local .lt. sza_Theis(jj) ) then jsza1 = jj - 1 jsza2 = jj endif enddo dfx = (t_Theis(jsza2,1) - t_Theis(jsza1,1)) & / (cos(sza_Theis(jsza2)*deg2rad) & -cos(sza_Theis(jsza1)*deg2rad)) dT = t_Theis(jsza1,1) + & dfx*(cos(sza_local*deg2rad) & -cos(sza_Theis(jsza1)*deg2rad)) dfy = (tt-dT) / (z_Theis(1)-130.) temp_elect = dT + dfy*(z_Theis(1)-zkm) else dfy = (tt-t_Theis(nsza_Theis,1)) / (z_Theis(1)-130.) temp_elect = t_Theis(13,1) + dfy*(z_Theis(1)-zkm) endif else if(zkm .lt. 350. .and. sza_local .lt. 180.) then do ii=nz_Theis,2,-1 if ( zkm .lt. z_Theis(ii) ) then iz1 = ii - 1 iz2 = ii endif enddo do jj=nsza_Theis,2,-1 if ( sza_local .lt. sza_Theis(jj) ) then jsza1 = jj - 1 jsza2 = jj endif enddo dfx = t_Theis(jsza2,iz1) - t_Theis(jsza1,iz1) dfy = t_Theis(jsza1,iz2) - t_Theis(jsza1,iz1) dfxy= t_Theis(jsza1,iz1) + t_Theis(jsza2,iz2) & -(t_Theis(jsza2,iz1) + t_Theis(jsza1,iz2)) z_incre =(zkm-z_Theis(iz1))/(z_Theis(iz2)-z_Theis(iz1)) sza_incre =(cos(sza_local*deg2rad) & -cos(sza_Theis(jsza1)*deg2rad)) & /(cos(sza_Theis(jsza2)*deg2rad) & -cos(sza_Theis(jsza1)*deg2rad)) temp_elect = dfx*sza_incre + dfy*z_incre & - dfxy*sza_incre*z_incre + t_Theis(jsza1,iz1) else if (sza_local .ge. 180.) then if (zkm .lt. 350.) then do ii=nz_Theis,2,-1 if ( zkm .lt. z_Theis(ii) ) then iz1 = ii - 1 iz2 = ii endif enddo dfy = (t_Theis(nsza_Theis,iz2)-t_Theis(nsza_Theis,iz1)) & / (z_Theis(iz2)-z_Theis(iz1)) temp_elect = t_Theis(nsza_Theis,iz1)+dfy*(zkm-z_Theis(iz1)) else temp_elect = t_Theis(nsza_Theis,nz_Theis) endif endif if (zkm .ge. 350.) then do jj=nsza_Theis,2,-1 if ( sza_local .lt. sza_Theis(jj) ) then jsza1 = jj - 1 jsza2 = jj endif enddo dfx = (t_Theis(jsza2,nz_Theis) - t_Theis(jsza1,nz_Theis)) & / (cos(sza_Theis(jsza2)*deg2rad) & -cos(sza_Theis(jsza1)*deg2rad)) temp_elect = t_Theis(jsza1,nz_Theis) + & dfx*(cos(sza_local*deg2rad) & -cos(sza_Theis(jsza1)*deg2rad)) endif endif endif else if(origin.eq.2) then if (sza_local .gt. 180.) then sza_security = 180. else sza_security = sza_local endif if (zkm .lt. 130.) then temp_elect = tt else if (zkm .lt. Altimini) then log10_temp = P_indice_factor(sza_security,1) & + ( P_indice_factor(sza_security,2) & / ((Altimini + P_indice_factor(sza_security,3))**2.) ) & + Altimini * P_indice_factor(sza_security,4) factor_e = exp((zkm-Altimini)/10.) !temp_elect = 10.**(log10_temp) + (tt-10.**(log10_temp))* & ! (Altimini-zkm)/(Altimini-130.) temp_elect = 10.**(log10_temp)*factor_e + tt*(1.-factor_e) else log10_temp = P_indice_factor(sza_security,1) & + ( P_indice_factor(sza_security,2) & / ((zkm + P_indice_factor(sza_security,3))**2.) ) & + zkm * P_indice_factor(sza_security,4) temp_elect = 10.**(log10_temp) endif else !MAVEN measured electronic temperature (Ergun et al., GRL 2015) !Note that the Langmuir probe is not sensitive below ~500K, so !electronic temperatures in the lower thermosphere (<150 km) may !be overestimated by this formula !if(zkm.le.120) then ! temp_elect = tt !else ! temp_elect=((3140.+120.)/2.)+((3140.-120.)/2.)*tanh((zkm-241.)/60.) !endif !else !write(*,*)'Error in function temp_elect:' !write(*,*)'Origin must be either 1 or 2' !write(*,*)'Using neutral temperature instead' temp_elect = tt endif return end function temp_elect !*********************************************************************** function temp_ion(zkm,tt,sza_local,origin) ! Computes the ion temperature, from VIRA Model based on the model ! of Miller et al. 1980 & 1984. ! or NEUTRAL (origin=2) <=== NOT USED !*********************************************************************** ! Arguments real :: tt ! Temperature real :: zkm ! Altitude in km integer :: origin ! Miller et al 1984 model (origin=1) or NEUTRAL (origin=2) real :: sza_local ! solar zenith angle in degree ! local variables: logical, save :: firstcall = .true. real,parameter :: pi = 4.*atan(1.0) real,parameter :: deg2rad = pi/180. real :: temp_ion ! electronic temperatures real :: z_incre, sza_incre integer :: ii, iz1, iz2, jj, jsza1, jsza2 real :: dfx, dfy, dfxy, dT ! origin = 1 real :: sza_security ! degree real, parameter :: Altimini = 150. ! km - doit etre egal a z_Miller(1) real :: log10_temp, factor_e real :: zkm_security ! km integer, parameter :: nsza_Miller = 9 integer, parameter :: nz_Miller = 11 real, save :: t_Miller(nsza_Miller,nz_Miller) real, save :: z_Miller(nz_Miller), sza_Miller(nsza_Miller) if (firstcall) then ! 2022/09/17: IPLS-Venus GCM ayant zmax < 300 km, on a mis les donnees pour z < 350 km z_Miller(1:11) = (/ 150.,160.,170.,180.,190.,200.,225., & 250.,275.,300.,350. /) sza_Miller(1:9) = (/ 20.,40.,60.,75.,85.,95.,110.,135.,160. /) t_Miller(1,1:11)= (/ 300., 380., 470., 560., 580., 600., 660., & 930.,1400.,1900.,2700. /) t_Miller(2,1:11)= (/ 310., 390., 490., 580., 640., 620., 660., & 940.,1500.,1900.,2400. /) ! Il y a une erreur dans la temperature ion dans CHAPTER IV de VIRA ! du coup, j'ai pris l'article MILLER et al. 1984 pour 65 SZA t_Miller(3,1:11)= (/ 330., 410., 490., 540., 600., 610., 670., & 880.,1300.,1700.,2100. /) t_Miller(4,1:11)= (/ 310., 390., 510., 580., 600., 600., 670., & 920.,1400.,1700.,2000. /) t_Miller(5,1:11)= (/ 320., 430., 530., 600., 640., 660., 840., & 1100.,1300.,1500.,1800. /) t_Miller(6,1:11)= (/ 290., 440., 530., 580., 640., 750.,1100., & 1300.,1400.,1500.,1700. /) t_Miller(7,1:11)= (/ 230., 410., 600., 800.,1100.,1300.,1500., & 1700.,1700.,1800.,1800. /) t_Miller(8,1:11)= (/ 210., 310., 540., 910.,1500.,2000.,2400., & 2700.,2600.,2500.,2400. /) t_Miller(9,1:11)= (/ 230., 350., 500., 800.,1200.,1700.,2800., & 3500.,4100.,4700.,5500. /) if (origin.gt.1) then write(*,*)'Error in function temp_ion:' write(*,*)'Origin must be either 1' write(*,*)'Using neutral temperature instead' endif firstcall = .false. endif if(origin.eq.1) then ! Altitude security if (zkm .gt. z_Miller(nz_Miller)) then zkm_security = z_Miller(nz_Miller) else zkm_security = zkm endif if (zkm_security .lt. 130.) then temp_ion = tt else ! SZA security if (sza_local .gt. sza_Miller(nsza_Miller)) then sza_security = sza_Miller(nsza_Miller) jsza1 = nsza_Miller - 1 jsza2 = nsza_Miller else if (sza_local .lt. sza_Miller(1)) then sza_security = sza_Miller(1) jsza1 = 2 - 1 jsza2 = 2 else sza_security = sza_local do jj=nsza_Miller,2,-1 if ( sza_security .lt. sza_Miller(jj) ) then jsza1 = jj - 1 jsza2 = jj endif enddo endif if(zkm_security .lt. Altimini) then dfx = (t_Miller(jsza2,1) - t_Miller(jsza1,1)) & / (cos(sza_Miller(jsza2)*deg2rad) & -cos(sza_Miller(jsza1)*deg2rad)) dT = t_Miller(jsza1,1) + & dfx*(cos(sza_security*deg2rad) & -cos(sza_Miller(jsza1)*deg2rad)) dfy = (tt-dT) / (z_Miller(1)-130.) temp_ion = dT + dfy*(z_Miller(1)-zkm_security) else do ii=nz_Miller,2,-1 if ( zkm_security .lt. z_Miller(ii) ) then iz1 = ii - 1 iz2 = ii endif enddo dfx = t_Miller(jsza2,iz1) - t_Miller(jsza1,iz1) dfy = t_Miller(jsza1,iz2) - t_Miller(jsza1,iz1) dfxy= t_Miller(jsza1,iz1) + t_Miller(jsza2,iz2) & -(t_Miller(jsza2,iz1) + t_Miller(jsza1,iz2)) z_incre =(zkm_security-z_Miller(iz1)) & /(z_Miller(iz2)-z_Miller(iz1)) sza_incre =(cos(sza_security*deg2rad) & -cos(sza_Miller(jsza1)*deg2rad)) & /(cos(sza_Miller(jsza2)*deg2rad) & -cos(sza_Miller(jsza1)*deg2rad)) temp_ion = dfx*sza_incre + dfy*z_incre & - dfxy*sza_incre*z_incre + t_Miller(jsza1,iz1) endif endif else !MAVEN measured electronic temperature (Ergun et al., GRL 2015) !Note that the Langmuir probe is not sensitive below ~500K, so !electronic temperatures in the lower thermosphere (<150 km) may !be overestimated by this formula !if(zkm.le.120) then ! temp_elect = tt !else ! temp_elect=((3140.+120.)/2.)+((3140.-120.)/2.)*tanh((zkm-241.)/60.) !endif !else temp_ion = tt endif return end function temp_ion !*********************************************************************** function temp_mars_elect(zkm,tt,sza_local,origin) ! Computes the electronic temperature, either from Viking (origin=1) ! or MAVEN (origin=2) !*********************************************************************** ! Arguments real tt ! Temperature real zkm ! Altitude in km integer origin ! Viking (origin=1) or MAVEN (origin=2) real sza_local ! solar zenith angle ! local variables: real temp_mars_elect ! electronic temperatures real zhanson(9),tehanson(9) real incremento integer ii, i1, i2 zhanson(1:9) = (/ 120.,130.,150.,175.,200.,225.,250.,275.,300. /) tehanson(2:9) = (/ 200.,300.,500.,1250.,2000.,2200.,2400.,2500. /) tehanson(1) = tt if(origin.eq.1) then if ( zkm .le. 120. ) then temp_mars_elect = tt else if(zkm .ge.300.) then temp_mars_elect=tehanson(9) else do ii=9,2,-1 if ( zkm .lt. zhanson(ii) ) then i1 = ii - 1 i2 = ii endif enddo incremento=(tehanson(i2)-tehanson(i1)) / & (zhanson(i2)-zhanson(i1)) temp_mars_elect = tehanson(i1) + & (zkm-zhanson(i1)) * incremento endif else if(origin.eq.2) then !MAVEN measured electronic temperature (Ergun et al., GRL 2015) !Note that the Langmuir probe is not sensitive below ~500K, so !electronic temperatures in the lower thermosphere (<150 km) may !be overestimated by this formula if(zkm.le.120) then temp_mars_elect = tt else temp_mars_elect=((3140.+120.)/2.) + & ((3140.-120.)/2.)*tanh((zkm-241.)/60.) endif else write(*,*)'Error in function temp_elect:' write(*,*)'Origin must be either 1 or 2' write(*,*)'Using neutral temperature instead' temp_mars_elect = tt endif return end function temp_mars_elect END MODULE iono_h