C********************************************************** real function PLAS1(ALT) C T.L. Gulyaeva....................................................May 2014 CC CC==== IZMIRAN PLASMASPHERE ELECTRON DENSITY PROFILE ABOVE h05top CC Ref. Gulyaeva T.L., X.Huang, and B.W. Reinisch. Plasmaspheric extension of CC topside electron density profiles. Adv. Space Res., 29, No.6, 825-831, 2002. CC CC T.Gulyaeva............................., July 2004 - September 2013. CC Formulation is changed fitting plasmasphere model to topside ionosphere model at CC the Topside Transition Ionosphere Height CC HEITOP at, Electron Density ELTOP = 0.5*NmF2. CC Output : electron density, PLAS, [m^-3], at altitude X [km] above HEITOP CC C C Parameters used by function PLAS1: C C W = R12 = Rzs 12-months smoothed sunspot number C GLON - geographic longitude C HOUR - local time, UT - universal time, C NEW= GMLT Geomagnetic Local Time C LDAYNR - day-of-year C ALON, BLON - corrected geomagnetic coordinates C HEITOP = h05top - altitute in the topside at Ne05 = 0.5* NmF2 C X - altitude in the range: [h05top, 36,000 km] C INTEGER ND real N6370,L,LT real eltop COMMON /bplas/w,glat,glon,put,ldaynr,heitop,eltop PARAMETER ( PI=3.14159265, RE=6370., RAD=PI/180.) C C Specify corrected geomagnetic coordinates alon, blon C SMI SUBR. ABLON TO PREPARE FOR PLASMASPHERE MODEL alon=365. blon=365. if (((glon.ge.130.).and.(glon.le.178.)).or.((glon.ge.330.). &and.(glon.le.340.))) goto 209 CALL ABLON(GLAT,GLON,ALON,BLON) 209 continue HHUP=HEITOP ND=LDAYNR LT=PUT+GLON/15.0 if (LT.ge.24.) LT=LT-24.0 if (LT.lt.0.) LT=LT+24.0 L= 1.0 + ALT/RE c c- A1 = 1.-0.1*COS(PI*(LT-4.)/12.) C A1-cor - corrected A1(L,GMLT) xn = 0.9 xd = 1.-0.1*cos(PI/L**5) c- A1 = xn+(xd-xn)/(1.+exp(7.-GMLT))+(xn-xd)/(1.+exp(17.-GMLT)) A1 = xn+(xd-xn)/(1.+exp(11.-LT))+(xn-xd)/(1.+exp(21.-LT)) ! LT function A2 = 1.+0.2*SQRT(W) A3 = 1.+0.001*W AL = 1.-EXP(-0.04*L**4) AL2 = 1.-EXP(-0.04*2**4) C !!!! cold c1=3.0E+09 C !!New!! CORRECTION: C1=3E+09 AS UPPER LIMIT OF COEFFICIENT c- C1=.0002*ELTOP*(RE-HHUP) C New expression for C1(GMLT) c1=(1.0+cos(-pi*(LT-12.)/24.)/2.)*1.0E+09 C!!!! if (alon.le.360.) goto 10 A0 = (1.0+0.7*COS(RAD*(GLON+70.)))*(1.+COS(2*PI*(ND+16)/365.25)) BL = 0.5*A0*AL-0.7*L*EXP(0.1*A0*AL) BL2 = 0.5*A0*AL2-1.4*EXP(0.1*A0*AL2) c- N6370 = 3.0E9*A1*A2*EXP(A3*BL2) N6370 = c1*A1*A2*EXP(A3*BL2) c+!!! old formula: H = 1342.5/ALOG(CN1000/N6370) IS REPLACED BY : H=(RE-HHUP)*.25/ALOG(ELTOP/N6370) if (ALT.le.6370.) then C!!! old formula: DN = CN1000*EXP((1000.-ALT)/H/L/L) IS REPLACED BY : DN=ELTOP*EXP((HHUP-ALT)/H/L/L) else C!!! old formula: DN = 3.0E9*A1*A2*EXP(A3*BL) IS REPLACED BY : DN=C1*A1*A2*EXP(A3*BL) endif Plas1=DN RETURN 10 a0a= (1.0+0.7*COS(RAD*(ALON+70.)))*(1.+COS(2*PI*(ND+16)/365.25)) b0b= (1.0+0.7*COS(RAD*(BLON+70.)))*(1.+COS(2*PI*(ND+16)/365.25)) bla= 0.5*a0a*AL-0.7*L*EXP(0.1*a0a*AL) blb= 0.5*b0b*AL-0.7*L*EXP(0.1*b0b*AL) bl2a= 0.5*a0a*AL2-0.7*2*EXP(0.1*a0a*AL2) bl2b= 0.5*b0b*AL2-0.7*2*EXP(0.1*b0b*AL2) drea = 3.0E9*A1*A2*EXP(A3*bl2a) dreb = 3.0E9*A1*A2*EXP(A3*bl2b) drea=alog10(drea) dreb=alog10(dreb) delon=blon-alon if (delon.le.0.) delon=delon+360. deglon=glon-alon if (deglon.le.0.) deglon=deglon+360. dne=drea+(dreb-drea)*deglon/delon N6370=10.**dne C!!! H = 1342.5/ALOG(CN1000/N6370) H=(RE-HHUP)*.25/ALOG(ELTOP/N6370) if (ALT.le.6370.) then C!!! DN = CN1000*EXP((1000.-ALT)/H/L/L) DN=ELTOP*EXP((HHUP-ALT)/H/L/L) else C!!! dna = 3.0E9*A1*A2*EXP(A3*bla) dna = C1*A1*A2*EXP(A3*bla) dna=alog10(dna) C!!! dnb = 3.0E9*A1*A2*EXP(A3*blb) dnb = C1*A1*A2*EXP(A3*blb) dnb=alog10(dnb) dnx=dna+(dnb-dna)*deglon/delon dn=10.**dnx endif Plas1=DN RETURN END C C------------------------------------------------------------------ SUBROUTINE ablon(xlati,xlongi,alonn,blonn) C==== GRID POINTS FOR THE FUNCTION PLAS(alt,...) DIMENSION xlat(37),xlon(37),zlon(6) DATA xlat /10.62,10.43,10.06,9.60,9.10,8.76,8.60,8.74,9.29,9.82, &9.83,9.27,8.38,7.48,7.07,7.05,6.25,4.07,1.76,0.36,-0.42,-1.37, &-2.84,-4.01,-4.15,-4.37,-6.34,-9.73,-12.75,-14.46,-14.12,-10.41, &-3.90,2.66,7.44,9.95,10.62/ DATA xlon /0.,10.,20.,30.,40.,50.,60.,70.,80.,90.,100.,110.,120., &130.,140.,150.,160.,170.,180.,190.,200.,210.,220.,230.,240.,250., &260.,270.,280.,290.,300.,310.,320.,330.,340.,350.,360./ do 1 i=1,36 k=i if ((xlongi.ge.xlon(i)).and.(xlongi.lt.xlon(i+1))) goto 2 1 CONTINUE 2 fimeq=xlat(k)+(xlat(k+1)-xlat(k))*(xlongi-xlon(k))*0.1 if ((xlongi.lt.178.).or.(xlongi.gt.330.)) goto 3 if (xlati.le.(fimeq+1.)) goto 11 if (xlati.lt.2.6) goto 4 alonn=178. blonn=330. goto 11 4 nl=0 do 5 i=1,36 c k=i if (((xlati.ge.xlat(i)).and.(xlati.lt.xlat(i+1))).or. &((xlati.le.xlat(i)).and.(xlati.gt.xlat(i+1)))) goto 15 goto 5 15 nl=nl+1 zlon(nl)=xlon(i)+10.*(xlati-xlat(i))/(xlat(i+1)-xlat(i)) 5 continue alonn=zlon(1) blonn=zlon(2) goto 11 3 if (xlati.ge.(fimeq-1.)) goto 11 if (xlati.gt.7.6) goto 6 alonn=340. blonn=130. goto 11 6 nl=0 do 7 i=1,36 c k=i if (((xlati.ge.xlat(i)).and.(xlati.lt.xlat(i+1))).or. &((xlati.le.xlat(i)).and.(xlati.gt.xlat(i+1)))) goto 17 goto 7 17 nl=nl+1 zlon(nl)=xlon(i)+10.*(xlati-xlat(i))/(xlat(i+1)-xlat(i)) 7 continue if ((xlongi.ge.zlon(1)).and.(xlongi.le.130.)) goto 8 alonn=zlon(nl) blonn=zlon(1) goto 11 8 continue do 9 i=1,nl k=i if ((xlongi.ge.zlon(i)).and.(xlongi.le.zlon(i+1))) goto 10 9 continue 10 alonn=zlon(k) blonn=zlon(k+1) 11 RETURN END C C C********************************************************** C New: T.L. Gulyaeva.................................Sep 30, 2013 subroutine prepPlas(hm,aep,bb) C Determine height h05top(0.5*NmF2) for NeQuick Plasmasphere starting height: real NeNeQ dimension hm(3),aep(3),bb(6) COMMON /bplas/r12i,glat,ylon,put,ldaynr,heitop,eltop C C Start iterations to obtain h05top: c eltop=aep(1)*1E+11/8. ! 0.5*NmF2 htop=hm(1) ! start from hmF2 eltop=NeNeQ(htop,hm,aep,bb)/2. ! 0.5NmF2 dh=100. iflag=0 1 htop=htop+dh e1=NeNeQ(htop,hm,aep,bb) if (iflag.gt.3) goto 2 if (e1.lt.eltop) then htop=htop-dh dh=dh/10. iflag=iflag+1 goto 1 else goto 1 endif 2 heitop=htop RETURN end C*********************************************************************** real function NeNeQ(h,hm,aep,bb) real NeMdGR dimension hm(3),aep(3),bb(6) COMMON /bplas/r12i,glat,ylon,put,ldaynr,heitop,eltop if ((h.gt.heitop).and.(heitop.gt.hm(1))) then C choose option Plasmasphere: NeNeq=plas1(h) return endif if (h.gt.hm(1)) then C choose option topside Ionophere aNmax=NeMdGR(aep,hm,bb,hm(1)) NeNeQ=topq(h,aNmax,hm(1),bb(6)) return endif C choose option bottomside Ionophere NeNeQ=NeMdGR(aep,hm,bb,h) return end C------------------------------------------------------ real function NeMdGR(aep,hm,bb,h) dimension aep(3),hm(3),bb(6),B(3) save parameter (f1=10.0, f2=1.0) parameter (h0=90.0, Hd=5.0) B(1)=bb(5) B(2)=bb(3) B(3)=bb(1) if (h.gt.hm(3)) B(3)=bb(2) if (h.gt.hm(2)) B(2)=bb(4) sum=0.0 do jj=1,3 arg0=(h-hm(jj)) arg=arg0/B(jj) if (jj.gt.1) then d=abs(h-hm(1)) arg=arg*exp(f1/(1.0+f2*d)) endif if (h.lt.h0) arg=arg*(Hd+h0-h)/Hd if (abs(arg).gt.25.0) then s0=0.0 else ee=exp(arg) s0=aep(jj)*ee/(1.0+ee)**2 endif sum=sum+s0 enddo NeMdGR=sum*1.0E11 return end C------------------------------------------------------------- real function topq(h,aNo,hmax,Ho) parameter (g=0.125,rfac=100.0) dh=h-hmax gh=g*dh z=dh/(Ho*(1.0+rfac*gh/(rfac*Ho+gh))) ee=fexp(z) if (ee.gt.1.0E11) then ep=4.0/ee else ep=4.0*ee/(1.0+ee)**2 endif topq=aNo*ep return end C------------------------------------------------------------- REAL FUNCTION DINTEC(h1,h2,hm,aep,bb) C==== TECion+TECpl INTEGRATION PROCEDURE TEC/m^3 REAL NeNeQ dimension hm(3),aep(3),bb(6) common /bplas/r12i,glat,ylon,gmlt,ldaynr,heitop,eltop dnemin=NeNeQ(h1,hm,aep,bb) c k=ifix((ah-hmin)/20.+0.5) C if ((h2-h1).gt.1000.) then st=20. ! step=20 km above heitop else st=10.0 endif k=int((h2-h1)/st+0.5) x1=h1 y1=dnemin SUM=0.0 do 1 l=1,k x2=x1+st y2=NeNeQ(x2,hm,aep,bb) SUM=SUM+(y2+y1)*st x1=x2 y1=y2 1 continue dintec=sum*500. return end C C------------------------------------------------------------- subroutine subinteg(h2,tec1,hm,aep,bb) C integration along vertical profile C h2 = upper boundary for vTEC [80.0,hpl] dimension hm(3),aep(3),bb(6) COMMON /bplas/r12i,glat,ylon,put,ldaynr,heitop,eltop h0=80.0 h1=hm(1) crem h2=20200.0 if (h2.le.heitop) then tec1=dintec(h0,h2,hm,aep,bb) goto 25 else h1a=hm(1) h1b=heitop tec1=dintec(h0, h1a,hm,aep,bb) tec2=dintec(h1a,h1b,hm,aep,bb) tec3=dintec(h1b,h2, hm,aep,bb) tec4=tec1+tec2+tec3 endif 25 tec1=tec4/1.0E15 !TECU/10. equals to 10^15 m^-2 RETURN END C----------------------------------------------------------------------------