C IRIFO7-A, OCTOBER 1979
C INTERNATIONAL REFERENCE IONOSPHERE (IRI).
C THIS PROGRAM PRODUCES HEIGHT PROFILES OF TEMPERATURE
C AND DENSITY OF THE ELECTRONS AND IONS IN THE IONOSPHERE FOR
C SPECIFIED LOCATION,TIME AND SOLAR ACTIVITY.
C ADDRESS@ PROF.K.RAWER, HERRENSTR.43
C D 7801 MARCH , F.R.G.
C INPUT@
C 1.) LATITUDE(F6.1),LONGITUDE(F6.1),SOLAR-ACTIVITY(R)(F6.1),
C MONTH(F4.1),HOUR(F4.1,1X),XHI(F6.1),BEGIN,END,STEPWIDTH(OF THE
C REQUIRED HEIGHT RANGE)(3(F6.1,1X)),KOBE(I1,1X),JNE,JNEMAX,
C JTN,JTE,JTI,JTETI,JO,JHHE,JO2,JNO,JMAG(11I1)
C 2.) HMF2/KM AND NMF2/M-3 OR FOF2/MHZ (F5.1,E9.3)
C IF HMF2 LESS 10.0 THEN THE CCIR VALUES FOR FOF2 AND HMF2
C ARE TAKEN BY USE OF THE CCIR-TAPE FOR FOF2 AND M3000.
C THE USER SHOULD ADAPT PROCEDURE CCIRCA TO HIS OWN
C CCIR-TAPE.
C KOBE=0 SHOULD BE USED FOR PRINTER OUTPUT AND KOBE=1 FOR
C PUNCHER OR TAPE OUTPUT.THE TEN INPUT VARIABLES
C (JNE...JNO) ARE SWITCHED (1=YES,0=NO) TO CHOOSE YOUR
C PARAMETERS.JMAG=0 MEANS GEOGRAFIC JMAG=1 GEOMAGNETIC
C LATITUDE AND LONGITUDE.
C
C-WDC-A THIS PROGRAM HAS BEEN MODIFIED AT WDC-A TO BE MORE MACHINE
C-WDC-A INDEPENDENT.
C-WDC-A NOTE 1) INIALIZATION OF VARIBLES IS NOT NECESSARY.
C-WDC-A NOTE 2) THIS VERSION OF THE PROGRAM DOES NOT ALLOW YOU
C-WDC-A TO CHOSE WHICH PRAMETERS ARE OUTPUT. A FIXED OUTPUT
C-WDC-A WAS USED TO ASSURE EXECUTION ON ANY MACHINE.
C-WDC-A NOTE 3) YOU WILL HAVE TO OPEN YOUR CHANNEL NUMBERS, ON A
C-WDC-A CONTROL DATA MACHINE THIS WOULD BE A PROGRAM CARD,ON
C-WDC-A A DATA GENERAL THIS WOULD BE AN OPEN CARD, ETC. THE
C-WDC-A CHANNEL VARIBLES SPECIFIED BY THE VARIBLES EGNR,AGNR
C-WDC-A AND KONSOL CAN BE CHANGED TO ANY VALUE YOU LIKE,
C-WDC-A SO LONG AS YOU REPLACE THOSE CHANGES IN YOUR CHANNEL
C-WDC-A OPENING STATEMENT.
C-WDC-A IN THIS VERSION:
C-WDC-A EGNR = 1 --INPUT DATA
C-WDC-A AGNR = 4 --OUTPUT- PRINTOUT OR DISK FILE
C-WDC-A KONSOL = 2 --INTERACTIVE DISPLAY OF INFORMATIVE STATEMENTS
C-WDC-A OR OUTPUT TO LINE PRINTER
C-WDC-A TAPE5 = 5 --CCIR COEFFICENTS
C-WDC-A NOTE 4) THIS VERSION OF IRIFO7 IS INTENDED TO USE THE CCIR
C-WDC-A COEFFICENTS PROVIDED WITH THIS PROGRAM TAPE. THE READ
C-WDC-A STATEMENTS IN CCIRCA MUST BE CHANGED IF YOU ARE USING
C-WDC-A A DIFFERENT SET OF COEFFICENTS.
C
PROGRAM IRI(INPUT,OUTPUT,TAPE1=INPUT,TAPE2=OUTPUT,TAPE4,
XTAPE5)
INTEGER EGNR,AGNR,SMONTH,SHOUR,DAYNR,DDO(4)
1,IIF(4),DO2(2)
REAL LATI,LONGI,MO2(3),MO(5),LOGE,MONTH,MODIP,NMF2,
2NDEL0,NMF1,NME,NMD,K,MAGBR,
1NHABR,NDEL,NDX,NEI,MM,MLAT,MLONG,NOBO2
DIMENSION F(3),B0F(2,2,8),OUTF(50,12),JF(10),RIF(4),
1PF1O(12),PF3O(12),HO(4),
2PF2O(4),HO2(2),CTNN(3),PG1O(80),PG2O(32),PG3O(80)
LOGICAL WINTER,SUMMER,SCHALT,EXT,NIGHT,F1REG,VALLEY,
1IRDUPP,CCIREI,ANF,TEMP,TEVAL,TIMP
COMMON/BLOCK1/HMF2,NMF2,HMF1/BLOCK2/B0,B1,C1,HZ,T,G(144),HST,STR
4/BLOCK3/HDX,HME,
1NME,HMD,NMD,HEF,D1,K,FP30,FP3U,FP1,FP2/BLOCK4/HB,HC,P9,P10,
2P11,P12,IRDL,OW,AGNR,IRDUPP,P(8)/BLOCK5/ZX,TNX,ATN,CTN(3)/
3BLOCK8/HS,TNHS,XSM(2),MM(3)/BLO10/BETA,ETA,DELTA,ZETA/
XBLO11/MONTH,LATI/BLOCK6/NIGHT,E(4)
EXTERNAL SSRSU,XE1,XE2,XE3,XE4,XE5,XE6,TEDER
C
C-WDC-A THIS DATA STATEMENT IS NECESSARY IF ALL PARAMETERS ARE TO BE
C-WDC-A OUTPUT. THE OPTION OF VARIBLE TYPES OF OUTPUT IS NOT
C-WDC-A ACTIVE IN THIS VERSION, THUS ALL OUTPUT TYPES ARE SUPPLIES
C-WDC-A BY MAKING ALL ELEMENTS OF THE ARRAY JF EQUAL 1.
C
DATA JF /10*1/
LOGE=1/ALOG(10.0)
UMR=0.01745329
PHI=3.1415927
CCIREI=.FALSE.
C "COMMENT" FIRST SPECIFY YOUR COMPUTERS CHANNEL NUMBERS, E.G.
EGNR=1
AGNR=4
KONSOL=2
C
C-WDC-1 IF YOU WISH TO READ MORE THAN ONE HOUR AT A TIME ADD AN
C-WDC-1 END OF FILE TEST BELOW. IF EOF TRUE THEN GO TO STATEMENT 999
C
2000 READ(EGNR,2) LATI,LONGI,R,MONTH,HOUR,XHI,AH,EH,SH,KOBE,
1JMAG,HMF2,FOF2
2 FORMAT(3F6.1,2F4.1,F6.1,1X,3(F6.1,1X),I1,1X,I1,F5.1,E9.3)
IOND=JF(8)+JF(9)+JF(10)+JF(7)
C
C INPUT OF TRANSFORMATION COEFFICIENTS (PROCEDURE FIELDG) ..................
C
CALL KOEFFI(B0F)
CALL KOEFFB(G)
C
C CALCULATION OF XHI,SUNSET,SUNRISE,COV,XHIM ...............................
C
IF(JMAG.LE.0) GOTO 888
MLAT=LATI
MLONG=LONGI
888 CALL GGM(JMAG,LONGI,LATI,MLONG,MLAT)
DAYNR=IFIX(MONTH)*30-15
SUNDEC=-0.40915*COS(2.0*PHI/365.25*FLOAT(DAYNR+8))
Z1=SIN(SUNDEC)*SIN(LATI*UMR)
Z2=COS(SUNDEC)*COS(LATI*UMR)
IF(ABS(Z2).GT.0.0) GOTO 120
SAX=24.0
IF(Z1.GE.0.0) SAX=0.0
GOTO 140
120 IF(ABS(Z1/Z2).LE.1.0) GOTO 510
SAX=24.0
IF((Z1+Z2).GE.0.0) SAX=0.0
GOTO 140
510 SAX=12.0-ACOS(-Z1/Z2)/(UMR*15.0)
140 SUX=24.0-SAX
XLSTA=15.0*(HOUR-12.0)
COV=63.75+R*(0.728+R*0.00089)
IF(XHI.GT.(-10.0)) GOTO 520
COSXHI=SIN(LATI*UMR)*SIN(SUNDEC)+COS(LATI*UMR)*COS(SUNDEC)
1*COS(XLSTA*UMR)
XHI=ACOS(COSXHI)/UMR
GOTO 503
520 COSXHI=COS(XHI*UMR)
503 XHI0=0.87+0.0061*ABS(LATI)
XHI100=0.68+0.0089*ABS(LATI)
XHIM=(XHI0+(XHI100-XHI0)*R/100.0)/UMR
CALL FIELDG(LATI,LONGI,300.0,XMA,YMA,ZMA,BET,DIP,MODIP)
MAGBR=ATAN(0.5*TAN(DIP*UMR))/UMR
SUNDEC=SUNDEC/UMR
C
C CLASSIFICATION OF DIFFERENT TIMES AND REGIONS ............................
C
F1REG=.TRUE.
VALLEY=.FALSE.
SUMMER=.FALSE.
WINTER=.FALSE.
NIGHT=.FALSE.
TEMP=.TRUE.
TIMP=.TRUE.
TEVAL=.FALSE.
IF((MONTH.GT.10.0).OR.(MONTH.LT.3.0)) WINTER=.TRUE.
IF((MONTH.GT.4.0).AND.(MONTH.LT.9.0)) SUMMER=.TRUE.
EXT=SUMMER
IF(LATI.GT.0.0) GOTO 1111
EXT=WINTER
WINTER=SUMMER
1111 SUMMER=EXT
IF((HOUR.GT.SUX).OR.(HOUR.LT.SAX)) NIGHT=.TRUE.
IF(XHI.GT.XHIM) F1REG=.FALSE.
IF(NIGHT) F1REG=.FALSE.
IF(WINTER) F1REG=.FALSE.
C
C INPUT OF THE ION DENSITY PARAMETER ARRAYS PF1O,PF2O AND PF3O......
C
IF(IOND.LT.1) GOTO 141
IIF(1)=2
IF(ABS(LATI).LT.30.0) IIF(1)=1
IIF(2)=2
IF(COV.LT.100.0) IIF(2)=1
IIF(3)=3
IF(WINTER) IIF(3)=4
IF(SUMMER) IIF(3)=2
IIF(4)=1
IF(NIGHT) IIF(4)=2
CALL KOEFP1(PG1O)
CALL KOEFP2 (PG2O)
CALL KOEFP3 (PG3O)
DO 6000 I=1, 4
RIF(I)=FLOAT(IIF(I))
6000 CONTINUE
CALL SUFE (PG1O,RIF,12,PF1O)
CALL SUFE(PG2O,RIF,4,PF2O)
CALL SUFE(PG3O,RIF,12,PF3O)
C
C CALCULATION OF ELECTRON DENSITY PARAMETERS...............................
C
141 DELL=4.32
IF(ABS(MODIP).GE.18.0) DELL=1.0+EXP(-(ABS(MODIP)-30.0)/10.0)
IF(.NOT.F1REG) GOTO 150
C1=0.1244-(4.44E-4)*R+0.09/DELL
FOF1=FOF1ED(MAGBR,R,COSXHI)
NMF1=1.24E10*FOF1*FOF1
150 XHINON=XHIS(MONTH,12.0,LATI)
FOE=FOEEDI(R,XHI,XHINON,LATI,HOUR,SUX,SAX)
IF(HMF2.GE.10.0) GOTO 501
CALL F2OUT(MODIP,LATI,LONGI,MAGBR,R,MONTH,HOUR,FOE,HMF2,FOF2)
GOTO 502
501 IF(FOF2.GT.100.0) FOF2=SQRT(FOF2/1.24E+10)
502 NMF2=1.24E10*FOF2*FOF2
NME=1.24E10*FOE*FOE
NMD=XMDED(XHI,R,4.0E8)
7000 HME=HPOL(HOUR,110.0,105.0,SAX,SUX)
HMD=HPOL(HOUR,81.0,88.0,SAX,SUX)
SMONTH=IFIX(MONTH/3.0)
IF(SMONTH.LT.1) SMONTH=4
B0=XPOL(B0F,HOUR,SAX,SUX,DELL,SMONTH,R)
COS2=COS(MLAT*UMR)
COS2=COS2*COS2
FLU=(COV-40.0)/30.0
ETA=0.058798-0.0070305*COS2+FLU*(-0.014065+0.0069724*COS2)+
1(0.0024287+0.0042810*COS2-0.00015280*FOF2)*FOF2
ZETA=0.078922-0.0046702*COS2+FLU*(-0.019132+0.0076545*COS2)+
1(0.0032513+0.0060290*COS2-0.00020872*FOF2)*FOF2
BETA=-128.03+20.253*COS2+FLU*(-8.0755-0.65896*COS2)+(0.44041
1+0.71458*COS2-0.042966*FOF2)*FOF2
XXX=EP1ST(-94.45,BETA)
XXXX=DEP1ST(-94.45,BETA)
DELTA=(ETA*XXX-ZETA/2.0)/(ETA*XXXX+ZETA/400.0)
B1=3.0
F(1)=HPOL(HOUR,0.02+0.03/DELL,0.05,SAX,SUX)
F(2)=HPOL(HOUR,4.6,4.5,SAX,SUX)
F(3)=HPOL(HOUR,-11.5,-4.0,SAX,SUX)
NDEL0=5.0
IF(WINTER) NDEL0=10.0
DNDH0=0.016
IF(SUMMER) DNDH0=0.01
NHABR=HPOL(HOUR,10.5/DELL,28.0,SAX,SUX)
XXX=EPSTEP(45.0,67.0,-10.0,20.0,ABS(LATI))
HBR=HPOL(HOUR,17.8/DELL,XXX,SAX,SUX)
NDEL=HPOL(HOUR,NDEL0/DELL,81.0,SAX,SUX)
DNDHBR=HPOL(HOUR,DNDH0/DELL,0.06,SAX,SUX)
IF(NDEL.GT.1.0) VALLEY=.TRUE.
IF(.NOT.VALLEY) GOTO 600
CALL TAL(HME,NME,NHABR,NDEL,HBR,1.0,DNDHBR,EXT,E)
IF(.NOT.EXT) GOTO 600
WRITE(KONSOL,650)
650 FORMAT(1X,"NO ELECTRON DENSITY E-VALLEY,AS THE
1 MODEL FUNCTION HAS A SECOND EXTREMUM IN THE VALLEY-REGION")
VALLEY=.FALSE.
600 FP1=F(1)
FP2=-FP1*FP1/2.0
FP30=(-F(2)*FP2-FP1+1.0/F(2))/(F(2)*F(2))
FP3U=(-F(3)*FP2-FP1-1.0/F(3))/(F(3)*F(3))
HDX=HMD+F(2)
X=HDX-HMD
NDX=NMD*EXP(X*(FP1+X*(FP2+X*FP30)))
DNDX=NDX*(FP1+X*(2.0*FP2+X*3.0*FP30))
X=HME-HDX
K=-DNDX*X/(NDX*ALOG(NDX/NME))
D1=DNDX/(NDX*K*X**(K-1.0))
IF(.NOT.VALLEY) HBR=0.0
IF(F1REG) GOTO 700
HMF1=0.0
NMF1=0.0
C1=0.0
700 HEF=HME+HBR
C
C SEARCH FOR HMF1,HST,HEF,AND HA.......................................
C
7107 IF(.NOT.F1REG) NMF1=(NME+NMF2)/2.0
H=0.0
133 H=H+10.0
IF(H.GT.(HMF2-HME)) GOTO 135
CALL REGFA1(HMF2-10.0-H,HMF2,1.0,NMF1,XE2,SCHALT,HMF1)
IF(SCHALT) GOTO 133
GOTO 137
135 WRITE(KONSOL,11)
11 FORMAT(1H ,1X," HMF1 IS
1 NOT EVALUATED BY THE FUNCTION XE2")
HMF1=HMF2
C1=0.0
137 IF(HMF1.GT.(HEF+10.0)) GOTO 380
C1=0.0
HMF1=HMF2
WRITE(KONSOL,9)
9 FORMAT(1H ,1X," HEF GREATER AS HMF1")
380 H=0.0
125 H=H+3.0
IF(H.GT.(HMF1-HME)) GOTO 900
CALL REGFA1(HMF1,HME+H,1.0,NME,XE3,SCHALT,HST)
STR=HST
IF(SCHALT) GOTO 125
GOTO 360
900 B1=B1+0.5
IF(B1.GE.10.0) GOTO 901
H=0.0
GOTO 133
901 WRITE(KONSOL,100)
100 FORMAT(1H ,1X," HST IS NOT
1EVALUATED BY THE FUNCTION XE3")
HST=(HMF1+HME-HBR)/2.0
360 HZ=(HST+HMF1)/2.0
WRITE (KONSOL,902) B1
902 FORMAT (1X, "WE PUT B1= ", F3.1, "TO GET HST")
IF(HZ.GT.(HEF+10.0)) GOTO 950
HST=(HMF1+HEF)/2.0
HZ=(HST+HMF1)/2.0
950 D=HZ-HST
HEF=HME+HBR
T=D*D/(HZ-HEF-D)
C
C CALCULATION OF NEUTRAL TEMPERATURE PARAMETER...........................
C
HTA=90.0
TNA=183.0
ZX=125.0
ANF=.TRUE.
Z1=-180.0
882 TUN=TUNCAL(COV,LATI,SUNDEC,Z1)
TNX=371.6678+0.0518806*TUN-294.3505*EXP(-0.00216222*TUN)
ATN=0.63662*(TUN-TNX)
HDEL=ZX-HTA
TDEL=TNX-TNA
HD2=HDEL*HDEL
CTN(1)=1.9*TDEL/HDEL
CTN(3)=3.0*TDEL/(HD2*HD2)
CTN(2)=CTN(3)*1.333333*HDEL-CTN(1)/HD2
IF(.NOT.ANF) GOTO 881
ANF=.FALSE.
TX=TN(130.0,TNX,ATN,CTN)
TNXN=TNX
ATNN=ATN
CTNN(1)=CTN(1)
CTNN(2)=CTN(2)
CTNN(3)=CTN(3)
Z1=XLSTA
GOTO 882
C
C CALCULATION OF ELECTRON TEMPERATURE PARAMETER......................
C
881 HTA=120.0
H0T=200.0
H0N=400.0
IF(ABS(MLAT).GT.40.0) H0T=350.0
H0=H0T
IF(NIGHT) H0=H0N
HMAX=70.0*EXP((-1.4E-3)*MLAT*MLAT)+200.0
CALL TELAT(MLAT,1600.0,700.0,0.5,0.47,0.024,0.0,F1N,F2N)
TX=TN(H0N,TNXN,ATNN,CTNN)
TEN=TE(H0N,F1N,F2N,HMAX,0.0,H0N,0.0,0.0,0.0)
DTEN=(TEN-TX)/(H0N-HTA)
CALL TELAT(MLAT,2325.0,725.0,1.0,3.4,-0.014,2.56E4,F1T,F2T)
TX=TN(H0,TNXN,ATNN,CTNN)
XXX=TE(H0,F1T,F2T,HMAX,2.0,H0T,0.0,0.0,0.0)
XXXX=TE(H0,F1N,F2N,0.0,0.0,H0N,TX,HTA,DTEN)
TE0=HPOL(HOUR,XXX,XXXX,SAX,SUX)
TX=TN(H0,TNX,ATN,CTN)
QUO=(TE0-TX)/(H0-HTA)
C
C CALCULATION OF ION TEMPERATURE PARAMETERS
C
887 XSM(1)=430.0
MM(2)=HPOL(HOUR,3.0,0.0,SAX,SUX)
Z1=EXP(-0.09*MLAT)
YSM2=1240.0-1400.0*Z1/((1.0+Z1)*(1.0+Z1))
X2=ABS(MLAT)
X1=X2*(0.47+X2*0.024)*UMR
X3=COS(X1)
X4=YSM2
TIN=1200.0-300.0*SIGN(1.0,X3)*SQRT(ABS(X3))
YSM2=TIN
IF(X4.GE.TIN) YSM2=HPOL(HOUR,X4,TIN,SAX,SUX)
TX=TN(XSM(1),TNXN,ATNN,CTNN)
XXX=TE(XSM(1),F1T,F2T,HMAX,2.0,H0T,0.0,0.0,0.0)
XXXX=TE(XSM(1),F1N,F2N,0.0,0.0,H0N,TX,HTA,DTEN)
Z1=HPOL(HOUR,XXX,XXXX,SAX,SUX)
Z2=TN(XSM(1),TNX,ATN,CTN)
IF(YSM2.LE.Z2) YSM2=(Z1+Z2)/2.0
CALL REGFA1(130.0,500.0,0.1,YSM2,TEDER,SCHALT,HS)
IF(.NOT.SCHALT) GOTO 250
HS=200.0
250 TNHS=TN(HS,TNX,ATN,CTN)
MM(1)=DTNDH(HS,ATN,CTN)
MM(3)=MM(2)
XSM(2)=XSM(1)
C
C SMOOTHING OF THE ION TEMPERATURE TO KEEP IT LESS THAN
C THE ELECTRON TEMPERATURE AT HSM
C
HSM=1000.0
XXX=TE(HSM,F1T,F2T,HMAX,2.0,H0T,0.0,0.0,0.0)
XXXX=TE(HSM,F1N,F2N,HMAX,0.0,H0N,0.0,0.0,0.0)
TESM=HPOL(HOUR,XXX,XXXX,SAX,SUX)
IF(TESM.GE.TI(HSM)) GOTO 240
XXX=DTEDH(HSM,F2T,HMAX,2.0,H0T,0.0,0.0)
XXXX=DTEDH(HSM,F2N,HMAX,0.0,H0N,0.0,0.0)
MM(3)=HPOL(HOUR,XXX,XXXX,SAX,SUX)
XSM(2)=(TESM-10.0-YSM2+MM(2)*XSM(1)-MM(3)*HSM)/(MM(2)-MM(3))
C
C CALCULATION OF ION DENSITY PARAMETER...............................
C
240 IF(IOND.LT.1) GOTO 189
Z1=0.0
IF(XHI.LE.90.0) Z1=COSXHI
HFIXO=300.0
IF((IIF(2).EQ.2).AND.(IIF(3).EQ.2)) HFIXO=249.0
MO(1)=EPSTEP(PF1O(1),PF1O(2),PF1O(3),PF1O(4),Z1)
MO(2)=EPSTEP(PF1O(5),PF1O(6),PF1O(7),PF1O(8),Z1)
MO(3)=0.0
HO(1)=EPSTEP(PF1O(9),PF1O(10),PF1O(11),PF1O(12),Z1)
HO(4)=PF2O(1)
MO(4)=PF2O(2)
MO(5)=PF2O(3)
DDO(1)=9
DDO(2)=5
DDO(3)=5
DDO(4)=50
7100 HO(2)=290.0
IF((IIF(2).EQ.2).AND.(IIF(3).EQ.2)) HO(2)=237.0
HO(3)=(4.60517-MO(5)*(HO(4)-PF2O(4)))/MO(4)+HO(4)
IF(HO(2).LT.HO(3)) GOTO 7101
MO(4)=MO(4)-0.001
GOTO 7100
7101 Z2=5.0
X=HO(2)
Z1=0.0
7102 X=X+Z2
Y=RPID(X,HFIXO,98.0,4,MO,DDO,HO)
IF(Y.LE.Z1) GOTO 7103
Z1=Y
GOTO 7102
7103 IF(Z2.LE.1.0) GOTO 7104
X=X-Z2
Z2=1.0
GOTO 7102
7104 H0O=X-0.5
IF(XHI.GE.90.0) COSXHI=0.0
DO 7105 I=2,4,2
L=I/2
HO2(L)=PF3O(1+I)+PF3O(2+I)*COSXHI
7105 MO2(L+1)=PF3O(7+I)+PF3O(8+I)*COSXHI
DO2(1)=5
DO2(2)=5
MO2(1)=PF3O(7)+PF3O(8)*COSXHI
7106 Y=RPID(H0O,PF3O(1),PF3O(2),2,MO2,DO2,HO2)
IF(Y.LE.0.1) GOTO 189
MO2(3)=MO2(3)-0.02
GOTO 7106
C
C CALCULATION FOR THE REQUIRED HEIGHT RANGE.......................
C
189 HA=65.0
IF(NIGHT) HA=80.0
IF(AH.LT.HA) AH=HA
IF(EH.GT.1000.0) EH=1000.0
KOMB=JF(4)+JF(5)+JF(6)
IF(((EH-AH)/SH+1.0).LT.50.0) GOTO 230
EH=AH+49.0*SH
WRITE(AGNR,190) EH
190 FORMAT(1H ,1X," TOO MANY HEIGHT STEPS,EH IS REDUCED
1 TO ",F4.0)
230 I=0
X=EH
300 I=I+1
IN=1
OUTF(I,IN)=X
IF((JF(1)+JF(2)).LE.0) GOTO 330
NEI=XE(X)
IF(JF(1).LE.0) GOTO 340
IN=IN+1
OUTF(I,IN)=NEI
340 IF(JF(2).LE.0) GOTO 330
IN=IN+1
OUTF(I,IN)=NEI/NMF2
330 IF(KOMB.LE.0) GOTO 7108
X1=TN(X,TNX,ATN,CTN)
IF(X.GE.HTA) GOTO 7109
Z2=-1.0
TEH=-1.0
GOTO 7110
7109 TX=TN(X,TNXN,ATNN,CTNN)
Z1=TE(X,F1N,F2N,HMAX,0.0,H0N,TX,HTA,DTEN)
Z2=X1
IF(X.GE.HS) Z2=TI(X)
IF(X.LT.H0) GOTO 7119
XXX=TE(X,F1T,F2T,HMAX,2.0,H0T,0.0,0.0,0.0)
TEH=HPOL(HOUR,XXX,Z1,SAX,SUX)
GOTO 7110
7119 TEH=TE(X,0.0,0.0,0.0,0.0,H0,X1,HTA,QUO)
7110 IF(JF(3).LE.0) GOTO 7112
IN=IN+1
OUTF(I,IN)=X1
7112 IF(JF(4).LE.0) GOTO 7113
IN=IN+1
OUTF(I,IN)=TEH
7113 IF(JF(5).LE.0) GOTO 7120
IN=IN+1
OUTF(I,IN)=Z2
7120 IF(JF(6).LE.0) GOTO 7108
IN=IN+1
OUTF(I,IN)=SIGN(1.0,TEH)*TEH/Z2
7108 IF(IOND.LE.0) GOTO 7118
Z1=RPID(X,HFIXO,98.0,4,MO,DDO,HO)
Z2=RPID(X,PF3O(1),PF3O(2),2,MO2,DO2,HO2)
C FIRST COMPUTE NOBO2, RATIO OF NO+TOO2+
C AT HEIGHT HB. THEN THE SAME RATIO AT X
Z1B= RPID (H0O,HFIXO,98.0,4,MO,DDO,HO)
Z2B= RPID (H0O,PF3O(1),PF3O(2),2,MO2,DO2,HO2)
NOBO2= (100.0-Z1B-Z2B)/Z2B
IF(JF(7).LE.0) GOTO 7114
IN=IN+1
OUTF(I,IN)=Z1
7114 IF(JF(8).LE.0) GOTO 7116
IN=IN+1
CALL RDHHE(X,H0O,Z1,Z2,NOBO2,10.0,OUTF(I,IN),OUTF(I,IN+1))
7116 IF(JF(9).LE.0) GOTO 7117
IN=IN+2
OUTF(I,IN)=Z2
7117 IF(JF(10).LE.0) GOTO 7118
IN=IN+1
OUTF(I,IN)=RDNO(X,H0O,Z2,Z1,NOBO2)
7118 X=X-SH
IF(X.GE.AH) GOTO 300
IEI=I
C
C OUTPUT ON THE SPECIFIED DEVICE.................................
C
IF(KOBE.NE.1) GOTO 7020
WRITE(AGNR,7030) (JF(I),I=1,10),IEI,LATI,LONGI,R,MONTH,HOUR
7030 FORMAT(1X,10I1,2X,I3,2(2X,F6.1),2X,F5.1,2X,F4.1,2X,F5.2)
7020 IF(KOBE.NE.0) GOTO 7041
WRITE (AGNR,7051)
7051 FORMAT (1X,"INPUT@")
IF (JMAG.LT.1) GOTO 7053
WRITE (AGNR,7052) MLAT,MLONG
7052 FORMAT (1X,"MLAT=",F6.1,2X,"MLONG=",F6.1)
GOTO 7055
7053 WRITE (AGNR,7054) LATI,LONGI
7054 FORMAT (1X,"LATI=",F6.1,2X,"LONGI=",F6.1)
7055 WRITE(AGNR,7050) R,MONTH,HOUR
7050 FORMAT(1X,"R=",F5.0," MONTH=",F4.1," HOUR=",F5.2)
WRITE (AGNR,7061)
7061 FORMAT (1X,"CALCULATED VALUES@")
IF (JMAG.LT.1) GOTO 7063
WRITE (AGNR,7062) LATI,LONGI
7062 FORMAT (1X,"LATI=",F6.1,2X,"LONGI=",F6.1)
GOTO 7065
7063 WRITE (AGNR,7064)MLAT,MLONG
7064 FORMAT (1X,"MLAT=",F6.1,2X,"MLONG=",F6.1)
7065 WRITE(AGNR,7060) DIP,MODIP,MAGBR,XHI
7060 FORMAT(1X,"DIP=",F6.1," MODIP=",F6.1," MAGLA=",F6.1," XHI=",
XF6.1)
WRITE(AGNR,7066) SAX,SUX,SUNDEC
7066 FORMAT(1X,"SUNRISE",F4.1,"L.T. SUNSET@",F4.1,"L.T.",5X,
X"SUNDEC.=",F6.1)
Z1=0.0
IF(F1REG) Z1=HMF1
Z2=0.0
IF(F1REG) Z2=NMF1
WRITE(AGNR,7070) NMF2,Z2,NME,NMD
7070 FORMAT(1X,"NMF2=",E7.2," NMF1=",E7.2," NME=",E7.2," NMD=",E7.2)
WRITE(AGNR,7080) HMF2,Z1,HME,HMD
7080 FORMAT(1X,"HMF2=",F5.1," HMF1=",F5.1," HME=",F5.1," HMD=",F5.1)
WRITE(AGNR,7998)
7998 FORMAT(1H0,3X,1HH,9X,2HNE,7X,6HN/NMAX,5X,2HTN,7X,2HTE,7X,2HTI,
16X,5HTE/TI,5X,4HRDO+,4X,4HRDH+,3X,5HRDHE+,3X,5HRDO2+,3X,5HRDNO+)
7041 DO 7040 I=1,IEI
7040 WRITE(AGNR,7999) (OUTF(I,L),L=1,IN)
7999 FORMAT(1X,F6.1,3X,E10.4,3X,F6.4,3X,3(F6.1,3X),F7.4,5(3X,F5.1))
IF(KOBE.EQ.1) WRITE(AGNR,35)
35 FORMAT(1X,"111111")
C
C-WDC-A ADD A "GO TO 2000" IF YOU WISH TO READ THE DATA FOR MORE THAN
C-WDC-A ONE HOUR TO A RUN.
C
999 STOP
END
C D.BILITZA,IPW FREIBURG,15.1.77 *******************************************
C "COMMENT" D. BILITZA,H.THIEMANN,K.RAWER IPW FREIBURG,
C AUG.79 **********************
C IRI-PROCEDURES ***************************************
C FREF.@ K. RAWER, S. RAMAKRISHNAN, D. BILITZA,
C \INTERNATIONAL REFERENCE IONOSPHERE 1978, U.R.S.I. BRUSSELS''
C FOR CALCULATING ELECTRON DENSITY AND TEMPERATURE AND ION TEMPERATURE
C AND RELATIVE DENSITY IN THE HEIGHT RANGE 70 - 1000 KM.?
C IRI-PROCEDURES
C \REF.@ K. RAWER, S. RAMAKRISHNAN, D. BILITZA'
C TO CALCULATR ELECTRON DENSITY AND TEMPERATURE AND ION TEMPERATURE
C AND RELATIVE DENSITY IN THE HEIGHT RANGE 70 - 1000 KM.
C I. ELECTRON DENSITY ------------------------------------------------------
FUNCTION EP1ST(X,BET)
EP1ST=1.0/(1.0+EXP(-X/BET))
RETURN
END
FUNCTION DEP1ST(X,BET)
U=EXP(-X/BET)
DEP1ST=U/(BET*(1.0+U)*(1.0+U))
RETURN
END
FUNCTION XE1(H)
C REPRESENTS THE ELECTRON DENSITY FOR HEIGHTS NOT GREATER 1000 KM
C AND NOT LESS HMF2,BY USING A HARMONIZED BENT-MODEL.
C GLOBAL PARAMETERS@ ETA,ZETA,BETA,DELTA,NMF2,HMF2
C "COMMENT" K.RAWER, S. RAMAKRISHNAN 1978. REPRESENTING ELECTRON DENSITY
C PROFILE
C FOR HEIGHTS NOT GREATER 1000 KM AND NOT LESS HMF2
C BY HARMONIZED BENT-MODEL ADMITTING VARIABILITY OF
C GLOBAL PARAMETERS@ ETA, ZETA, BETA, DELTA WITH GEOMAGNETIC LATITUDE
C MLAT, SOLAR FLUXFLU AND CRITICAL FREQUENCY FOF2.
C ALO GLOBAL ARE PEAK-PARAMETERS NMF2,HMF2?
REAL NMF2
COMMON/BLOCK1/HMF2,NMF2,HMF1/BLO10/BETA,ETA,DELTA,ZETA
X=(H-HMF2)/(1000.0-HMF2)*700.0+300.0-DELTA
Y=(1000.0-HMF2)/700.0*(BETA*ETA*ALOG((1.0+EXP((X-394.50)/BETA))/
1(1.0+EXP((-94.5-DELTA)/BETA)))+ZETA*(100.0*ALOG((1.0+EXP
2((X-300.0)/100.0))/(1.0+EXP(-DELTA/100.0)))-X+300.0-DELTA))
XE1=NMF2*EXP(-Y)
RETURN
END
REAL FUNCTION XE2(H)
C ELECTRON DENSITY FOR HEIGHTS LESS HMF2 AND NOT LESS HMF1
REAL NMF2
COMMON/BLOCK1/HMF2,NMF2,HMF1/BLOCK2/B0,B1,C1,HZ,T,G(144),HST,STR
X=(HMF2-H)/B0
IF(ABS(X).LT.1.0E-10) GOTO 100
XE2=NMF2*EXP(-X**B1)/COSH(X)
GOTO 200
100 XE2=NMF2
200 RETURN
END
REAL FUNCTION XE3(H)
C ELECTRON DENSITY FOR HEIGHTS LESS HMF1 AND NOT LESS HZ
REAL NMF2
COMMON/BLOCK1/HMF2,NMF2,HMF1/BLOCK2/B0,B1,C1,HZ,T,G(144),HST,STR
XE3=XE2(H)+NMF2*C1*SQRT(ABS(HMF1-H)/B0)
RETURN
END
REAL FUNCTION XE4(H)
C ELECTRON DENSITY FOR HEIGHTS LESS HZ AND NOT LESS HEF
COMMON/BLOCK2/B0,B1,C1,HZ,T,G(144),HST,STR/BLOCK3/
1HDX,HME,NME,HMD,NMD,HEF,D1,K,FP30,FP3U,FP1,FP2
A=((STR/HST-1)*(H-HZ))/(HEF-HZ)+1
XE4=XE3 (A*(HZ+T/2.0-SIGN(1.0,T)*SQRT(T*(HZ-H+T/4.0))))
RETURN
END
REAL FUNCTION XE5(H)
C ELECTRON DENSITY FOR HEIGHTS LESS HEF AND NOT LESS HME(VALLEY-REGION)
REAL NME,NMD,K
LOGICAL NIGHT
COMMON/BLOCK3/HDX,HME,NME,HMD,NMD,HEF,D1,K,FP30,FP3U,FP1,
1FP2/BLOCK6/NIGHT,E(4)
T3=H-HME
T1=T3*T3*(E(1)+T3*(E(2)+T3*(E(3)+T3*E(4))))
IF(NIGHT) GOTO 100
XE5=NME*(T1+1.0)
GOTO 200
100 XE5=NME*EXP(T1)
200 RETURN
END
REAL FUNCTION XE6(H)
C ELECTRON DENSITY FOR HEIGHTS LESS HME AND NOT LESS HA
REAL NME,NMD,K
COMMON/BLOCK3/HDX,HME,NME,HMD,NMD,HEF,D1,K,FP30,FP3U,FP1,FP2
IF(H.GT.HDX) GOTO 100
Z=H-HMD
FP3=FP3U
IF(Z.GT.0.0) FP3=FP30
XE6=NMD*EXP(Z*(FP1+Z*(FP2+Z*FP3)))
GOTO 200
100 Z=HME-H
XE6=NME*EXP(-D1*Z**K)
200 RETURN
END
REAL FUNCTION XE(H)
C ELECTRON DENSITY BEETWEEN HA(KM) AND 1000 KM.FUNCTIONS NE1...6 ARE USED.
C "COMMENT" ELECTRON DENSITY BEETWEEN HA(KM) AND 1000 KM
C SUMMARIZING PROCEDURES NE1....6?
REAL NMF2,NME,NMD,K
COMMON/BLOCK1/HMF2,NMF2,HMF1/BLOCK2/B0,B1,C1,
1HZ,T,G(144),HST,STR/BLOCK3/HDX,H
1ME,NME,HMD,NMD,HEF,D1,K,FP30,FP3U,FP1,FP2
IF(H.LT.HMF2) GOTO 100
XE=XE1(H)
GOTO 200
100 IF(H.LT.HMF1) GOTO 300
XE=XE2(H)
GOTO 200
300 IF(H.LT.HZ) GOTO 400
XE=XE3(H)
GOTO 200
400 IF(H.LT.HEF) GOTO 500
XE=XE4(H)
GOTO 200
500 IF(H.LT.HME) GOTO 600
XE=XE5(H)
GOTO 200
600 XE=XE6(H)
200 RETURN
END
C II. ELECTRON TEMPERATURE -------------------------------------------------
FUNCTION DTEDH(H,F2,HMAX,D,H0,DTNDH,A)
C D.BILITZA,1.7.78,CALCULATES THE HEIGHT DERIVATE OF THE
C ELECTRON TEMPERATURE MODEL PROCEDURE TE
C "COMMENT" D.BILITZA 1978, CALCULATES THE HEIGHT DERIVAT IVE
C ELECTRON TEMPERATURE PROFILE AFTER PROCEDURE TE?
IF(H.LT.H0) GOTO 100
Y=EXP(-0.03*(H-HMAX))
F=1.0+Y
DTEDH=D-F2*0.03*Y*(1.0-Y)/(F*F*F)
RETURN
100 DTEDH=DTNDH+A
RETURN
END
SUBROUTINE TELAT(PHI,A,B,XN,P1,P2,C,F1,F2)
C D.BILITZA,1.7.78,CALCULATES THE TEMPERATURE MODEL PARTS F1 AND F2
C DEPENDING ON GEOMAGNETIC LATITUDE PHI
C "COMMENT" D. BILITZA 1978, CALCULATES THE PARAMETERS F1 AND F2
C NEEDED IN THE ELECTRON TEMPERATURE MODEL PROCEDURE TE,
C BOTH DEPENDING ON GEOMAGNETIC LATITUDE (MLAT,INTERNAL,PHI)?
XPHI=ABS(PHI)
F=(P1+P2*XPHI)*XPHI
F=F*0.01745329
Y=COS(F)
F1=A-B*SIGN(1.0,Y)*ABS(Y)**XN
Y=EXP(-0.1*XPHI)
F=1.0+Y
F2=C*Y/(F*F)
RETURN
END
FUNCTION TE(H,F1,F2,HMAX,D,HH0,TNH,HTA,A)
C D.BILITZA,1.7.78,MODEL FOR THE ELECTRON TEMPERATURE IN THE HEIGHT RANGE
C 120-1000 KM.F1,F2 ARE THE LATITUDE DEPENDENT PARAMETERS CALCULATED BY
C TELAT.HMAX IS THE HEIGHT OF THE RELATIVE MAXIMUM,D THE TEMPERATURE
C HEIGHT GRADIENT,HH0 THE INTERSECTION HEIGHT FOR THE EXTRAPOLATION TO THE
C CIRA 72 VALUES TNA AT HTA AND A,TNA ARE THE COEFFICIENTS FOR THE
C EXTRAPOLATION FUNCTION
C "COMMENT" D. BILITZA 1978. ELECTRON TEMPERATURE PROFILE
C IN THE HEIGHT RANGE 120-1000 KM. SIMPLIFIED MODEL EQUATIONS
C FITTED WITH PLUGGE AND SPENNERS AEROS MODEL POTENTIAL OF 197
C \SPACE RES. ,197'. F1, F2 ARE LATITUDE DEPENDENT
C PARAMETERS GIVEN BY TELAT. HMAX IS THE HEIGHT OF A
C MAXIMUM, D THE TEMPERATURE HEIGHT GRADIENT, H0 THE INTERSECTION HEIGHT
C BELOW WHICH THE PROFILE IS BUILT TO APPROACH THE CIRA 72
C TEMPERATURE TNA AT HTA . A TNH ARE COEFFICIENTS OF THE FITTING FUNCTION?
IF(H.LT.HH0) GOTO 10
Y=EXP(-0.03*(H-HMAX))
F=1.0+Y
TE=F1+F2*Y/(F*F)+D*(H-700.0)
RETURN
10 TE=TNH+A*(H-HTA)
RETURN
END
C III. ION TEMPERATURE ----------------------------------------------------
FUNCTION TUNCAL(COV,XLATI,SD,SLSTA)
C CALCULATES THE EXOSPHERIC TEMPERATURE FOR COVINGTON-INDEX
C COV,LATITUDE XLATI,SOLARDEKLINATION SD AND LOCAL SOLAR TIME
C ANGLE SLSTA USING CIRA72-MODEL
UMR=0.01745329
TC=379.0+3.24*COV
ETA=ABS(XLATI-SD)/2.0
THETA=ABS(XLATI+SD)/2.0
H1=COS(ETA*UMR)**2.2
H2=SIN(THETA*UMR)**2.2
TD=1.0+0.3*H1
TN=1.0+0.3*H2
A=(TD-TN)/TN
TAU=SLSTA-37.0+6.0*SIN((SLSTA+43.0)*UMR)
X=COS(TAU/2.0*UMR)
TUNCAL=TC*TN*(1.0+A*SIGN(1.0,X)*ABS(X)**3.0)
RETURN
END
REAL FUNCTION TN(X,TTNX,ATN,CCTN)
C NEUTRAL TEMPERATURE
C "COMMENT" D. BILITZA, 1978, NEUTRAL TEMPERATURE PROFILE
C APPROXIMATING ROUGHLY CIRA72-MODEL?
DIMENSION CCTN(3)
Z=X-125.0
IF(Z.LE.0.0) GOTO 100
Y=Z**2.5
Y=(1.0+(4.5E-6)*Y)*Z
YY=CCTN(1)/ATN
Y=YY*Y
Y=ATAN(Y)*ATN
TN=TTNX+Y
GOTO 200
100 TN=TTNX+Z*(CCTN(1)+Z*Z*(CCTN(2)+Z*CCTN(3)))
200 RETURN
END
FUNCTION DTNDH(H,ATN,CCTN)
C DERIVATIVE OF NEUTRAL TEMPERATURE
C "COMMENT" D.BILITZA, 1978. APPROXIMATE HEIGHT DERIVATIVE OF NEUTRAL
C CIRA-72 MODEL TEMPERATURE PROFILE?
DIMENSION CCTN(3)
Z=H-125.0
IF(Z.GT.0.0) GOTO 100
DTNDH=CCTN(1)+Z*Z*(3.0*CCTN(2)+Z*4.0*CCTN(3))
RETURN
100 H1=CCTN(1)/ATN
H2=Z**2.5
H3=H1*Z*(1.0+(4.5E-6)*H2)
DTNDH=ATN/(1.0+H3*H3)*H1*(1.0+(15.75E-6)*H2)
RETURN
END
REAL FUNCTION TI(H)
C ION TEMPERATURE FOR HEIGHTS NOT GREATER 1000 KM AND NOT LESS HS
REAL MM,G(2)
COMMON/BLOCK8/HS,TNHS,XSM(2),MM(3)
G(1)=20.0
G(2)=50.0
SUM=MM(1)*(H-HS)+TNHS
DO 100 I=1,2
A=H-XSM(I)
IF((H-XSM(I)).LT.(100.0*G(I))) A=G(I)*ALOG(1.0+EXP((H-XSM(I))
1/G(I)))
B=HS-XSM(I)
IF((HS-XSM(I)).LT.(100.0*G(I))) B=G(I)*ALOG(1.0+EXP((HS-XSM(I))
1/G(I)))
100 SUM=SUM+(MM(I+1)-MM(I))*(A-B)
TI=SUM
RETURN
END
REAL FUNCTION TEDER(H)
C THIS FUNCTION TOGETHER WITH THE SUBROUTINE REGFA1 IS USED TO FIND THE
C HEIGHT WHERE DIFFERENCE BETWEEN TN AND TI BEGINS
C "COMMENT" THIS FUNCTION ALONG WITH PROCEDURE REGFA1
C ALLOWS TO FIND THE HEIGHT ABOVE WHICH TN BEGINS TO BE DIFERENT FROM TI?
REAL MM
COMMON/BLOCK5/ZX,TNX,ATN,CTN(3)/BLOCK8/HS,TNHS,XSM(2),MM(3)
TNH=TN(H,TNX,ATN,CTN)
DTDX=DTNDH(H,ATN,CTN)
TEDER=DTDX*(XSM(1)-H)+TNH
RETURN
END
C IV. ION RELATIVE PRECENTAGE DENSITY
REAL FUNCTION RPID (H, H0, N0, M, ST, ID, XS)
C THIS ANALYTIC FUNCTION IS USED TO REPRESENT THE RELATIVE
C PRECENTAGE DENSITY OF ATOMAR AND MOLECULAR OXYGEN IONS
C "COMMENT" H. THIEMANN, 1979. ANALYTIC FUNCTION DESCRIBING RELATIVE
C PERCENTAGE DENSITY OF ATOMIC AND MOLECULAR OXYGEN IONS.
C COEFFICIENTS CHOSEN SO AS TO APPROACH BELOW 200 KM
C DANILOV-SEMENOV"S DATA AND DUMBS-SPENNER"S ABOVE?
REAL N0
DIMENSION ID(4), ST(5), XS(4)
SUM=(H-H0)*ST(1)
DO 100 I=1,M
A=H-XS(I)
IF (A.LT.FLOAT(100*ID(I))) A=FLOAT(ID(I))*ALOG(1.0+EXP((H-XS(I))/
1FLOAT(ID(I))))
B=H0-XS(I)
IF (B.LT.FLOAT(100*ID(I))) B=FLOAT(ID(I))*ALOG(1.0+EXP((H0-XS(I))/
1FLOAT(ID(I))))
100 SUM=SUM+(ST(I+1)-ST(I))*(A-B)
SUM=N0*EXP(SUM)
IF (SUM.LT.1.0E-10) SUM=0.0
RPID=SUM
RETURN
END
SUBROUTINE RDHHE (H,HB,RDOH,RDO2H,RNO,PEHE,RDH,RDHE)
C RAWER, OCT.79,H+ AND HE+ RELATIVE PERECENTAGE
C DENSITY FOR HEIGHTS NOT GREATER THAN 1000 KM.
C RNO IS THE RATIO OF NO+TO O2+DENSITY AT H=HB.
IF (H-HB)100,100,200
100 RDH=0.0
GOTO 300
200 RDH=(100.0-RDOH-(1.0+RNO)*RDO2H)*(1.0-PEHE/100.0)
300 RDHE=RDH*PEHE/(100.0-PEHE)
RETURN
END
REAL FUNCTION RDNO(H,HB,RDO2H,RDOH,RNO)
C D.BILTIZA, 1978. NO+RELATIVE PERCENTAGE
C DENSITY FOR HEIGHTS NOT LESS THAN 100 KM.
IF (H-HB) 100,100,200
100 Y=100.0-RDO2H-RDOH
GOTO 300
200 Y=RNO*RDO2H
300 IF(Y.LT.0.0005) Y=0.0
RDNO=Y
RETURN
END
C END IRI-FUNCTIONS******************************************************
REAL FUNCTION XMDED (XHI,R,XW)
C BILITZA, 24.3.78, CALCULATES THE ELECTRON DENSITY OF
C THE D REGION RELATIVE MAXIMUM. XHI IS THE SUN ZENITH ANGLE,
C R THE ZUERICH SUNSPOT NUMBER AND XW THE DESIRED NIGHT VALUE
C "COMMENT" D. BILITZA, 1978, CALCULATES ELECTRON DENSITY OF
C D RELATIVE MAXIMUM. XHI IS SOLAR ZENITH ANGLE, R ZURICH
C SUNSPOT NUMBER AND NW THE ESTIMATED NIGHT VALUE.
C \SOURCE@MECHTLY-BILITZA,1974'.?
Y=6.05+0.088*R
YW=XW/1.0E8
Z=(-0.1/(ALOG(YW/Y)))**0.3704
SUXHI=ACOS(Z)
IF (SUXHI.LT.1.0472) SUXHI=1.0472
XXHI=XHI/57.2957795
IF (XXHI.GT.SUXHI) GOTO 100
X=COS(XXHI)
XMDED=Y*1.0E8*EXP(-0.1/X**2.7)
RETURN
100 XMDED=YW*1.0E8
RETURN
END
REAL FUNCTION FOEEDI(R,XHI,XHIM,SLATI,ST,SU,SA)
C BILITZA,17.5.1977,CALCULATES FOE BY THE ENDINBURGH-METHODE.
C INPUT@ COVINGTON-INDEX (COV),SUN-ZENITH-ANGLE ACTUAL (XHI) AND FOR
C MIDDAY (XHIM),LATITUDE (SLATI),SUNSET,AND SUNRISE (SU,SA,HOUR)
C "COMMENT" D.BILITZA, 1977, CALCULATES FOE BY THE EDINBURGH-METHOD.
C INPUT@ZUERICH SUNSPOT NUMBER(R),SOLAR ZENITH+ANGLE@ ACTUAL AT HOUR T @XHI
C AT NOON@ XHIM , LATITUDE (LATI),SUNSET AND SUNRISE
C (SU, SA/HOUR). \REF.@ KOURIS-MUGGELETON, CCIR DOC. 6/3/07 SEPT.73'?
XLATI=ABS(SLATI)
COV=63.75+R*(0.728+0.00089*R)
A=1.0+0.0094*(COV-66.0)
SL=COS(XLATI*0.0174533)
IF(XLATI.LT.32.0) GOTO 100
SM=0.11-0.49*SL
C=92.0+35.0*SL
GOTO 400
100 SM=-1.93+1.92*SL
C=23.0+116.0*SL
400 IF(XHIM.GT.89.0) XHIM=89.0
B=COS(XHIM*0.0174533)**SM
SP=1.31
IF(XLATI.GT.12.0) SP=1.2
DX=0.0
IF(XHI.GT.73.0.AND.XHI.LE.90.0) DX=(6.27E-13)*(XHI-50.0)**8.0
IF(ST.GT.SU) GOTO 300
IF(ST.GT.SA) GOTO 600
D=0.077**SP*EXP(-1.68*(SA-ST))
GOTO 200
600 D=COS((XHI-DX)*0.0174533)**SP
GOTO 200
300 D=0.077**SP*EXP(-1.01*(ST-SU))
200 FOE=A*B*C*D
SP=1.0+0.0098*R
SMIN=0.017*SP*SP
IF(FOE.LT.SMIN) FOE=SMIN
FOEEDI=(FOE)**0.25
RETURN
END
REAL FUNCTION FOF1ED(YLATI,R,COSXHI)
C D. BILITZA,CALCULATES FOF1 FOR THE LATITUDE LATI THE ZUERICH
C SUNSPOT NUMBER R AND THE COSINUS OF THE SUN ZENITH ANGLE
C COSXHI.
C \REF. E. D. DURCHARME E. A.,RADIO SCIENCE,8,10,837-839'
C "COMMENT" R. EYFRIG, 1979. CALCULATES FOF1 FOR DIP- LATITUDE LATI ,
C ZURICH SUNSPOT NUMBER R AND COSINUS OF SOLAR ZENITH
C ANGLE COSXHI. \REF. E.D.DUCHARME E. A. , RADIO SCIENCE,
C 8, 10, 837-839, HOWEVER WITH MAGNETIC INSTEAD OF GEOMAGNETIC LATITUDE?
XLATI=ABS(YLATI)
F0=4.35+XLATI*(0.0058-(1.2E-4)*XLATI)
F100=5.348+XLATI*(0.011-(2.3E-4)*XLATI)
FS=F0+(F100-F0)*R/100.0
XMUE=0.093+XLATI*(0.0046-(5.4E-5)*XLATI)+(3.0E-4)*R
FOF1ED=FS*COSXHI**XMUE
RETURN
END
REAL FUNCTION HMF2ED(XMAGBR,R,X,XM3)
C D. BILITZA,CALCULATES HMF2 FOR THE MAGNETIC LATITUDE XMAGBR
C AND THE ZUERICH SUNSPOT NUMBER R BY USING XM3 AND THE RATIO
C FOF2/FOE.
C \REF. D. BILITZA ET. AL. TO BE PUBLISHED IN@ TELECOMM.J.'
F1=(2.32E-3)*R+0.222
F2=1.2-(1.16E-2)*EXP((2.39E-2)*R)
F3=0.096*(R-25.0)/150.0
DELM=F1*(1.0-R/150.0*EXP(-XMAGBR*XMAGBR/1600.0))/(X-F2)+F3
HMF2ED=1490.0/(XM3+DELM)-176.0
RETURN
END
SUBROUTINE REGFA1(X11,X22,EPS,FW,F,SCHALT,X)
C REGULA-FALSI-PROCEDURE TO FIND X WITH F(X)-FW=0.X1,X2 ARE THE STARTING
C VALUES.IF THE X-INTERVALL IS LESS EPS,THE PROCEDURE ENDS.
C IF SIGN(1.0,F(X1)-FW)=SIGN(1.0,F(X2)-FW) SCHALT=TRUE.
C "COMMENT" REGULA-FALSI-PROCEDURE TO FIND X WITH F(X)-FW=0.
C X1, X2 ARE THE STARTING VALUES. THE COMOUTATION ENDS WHEN THE X-INTERVAL
C HAS BECOME LESS THAN EPS . IF SIGN(F(X1)-FW)= SIGN(F(X2)-FW)
C THE PROCEDURE JUMPS TO LABEL IRRTUM?
LOGICAL L1,LINKS,K,SCHALT
SCHALT=.FALSE.
X1=X11
X2=X22
F1=F(X1)-FW
F2=F(X2)-FW
IF(ABS(SIGN(1.0,F1)-SIGN(1.0,F2))-10E-10) 100,110,110
110 K=.FALSE.
NG=2
200 X=(X1*F2-X2*F1)/(F2-F1)
GOTO 400
300 L1=LINKS
IF(LINKS) GOTO 500
X=X1+FLOAT(NG-1)*(X2-X1)/FLOAT(NG)
GOTO 400
500 X=X1+(X2-X1)/FLOAT(NG)
400 FX=F(X)-FW
LINKS=SIGN(1.0,F1)*SIGN(1.0,FX).GT.0
K=.NOT.K
IF(LINKS) GOTO 600
X2=X
F2=FX
GOTO 700
600 X1=X
F1=FX
700 IF(ABS(X2-X1).LE.EPS) GOTO 800
IF(K) GOTO 300
IF((LINKS.AND.(.NOT.L1)).OR.(.NOT.LINKS.AND.L1)) NG=2*NG
GOTO 200
100 X=0.0
SCHALT=.TRUE.
800 RETURN
END
SUBROUTINE TAL(SHMAX,STMAX,SHABR,SDELTA,SHBR,STE0,SDTDH0,AUS6,SPT)
C D. BILITZA,15.1.77,CALCULATION OF THE COEF. SPT(4) FOR A POLYNOM OF
C 5TH DEGREE (K1=1,K2=SPT(1),...,K5=SPT(4),WHICH FITS TO THE
C VALLEY WITH THE PARAMETERS@ SHABR=HMIN-SHMAX,SDELTA IS THE
C PERCENTAGE DEPTH,SHBR THE WIDTH,STE0 THE UPPER FIXPOINT
C NORMALISED TO STMAX AND SDTDH0 IS THE LOGARITHMIC DEVIATION AT
C THIS FIXPOINT.IF THERE IS A MINIMUM OR A SECOND MAXIMUM IN THE
C VALLEY REGION AUS6=.TRUE.
DIMENSION SPT(4)
LOGICAL AUS6
X=SHABR
Y=SHBR
IF(SDELTA.LE.60.0) GOTO 100
Z1=ALOG(1.0-SDELTA/100.0)/(X*X)
Z2=ALOG(STE0)/(Y*Y)
Z3=SDTDH0/(2.0*Y)
GOTO 200
100 Z1=-SDELTA/(X*X*100.0)
Z2=(STE0-1.0)/(Y*Y)
Z3=SDTDH0*STE0/(2.0*Y)
200 Z4=X-Y
SPT(4)=2.0*(Z1*(Y-2.0*X)*Y-Z2*(X-2.0*Y)*X+Z3*Z4*X)/(X*Y*Z4*Z4*Z4)
SPT(3)=(Z1*(2.0*Y-3.0*X)+Z2*X)/(X*Z4*Z4)-(2.0*X+Y)*SPT(4)
SPT(2)=-2.0*Z1/X-2.0*X*SPT(3)-3.0*X*X*SPT(4)
SPT(1)=Z1-X*(SPT(2)+X*(SPT(3)+X*SPT(4)))
AUS6=.FALSE.
Z1=(4.0*SPT(3)+5.0*SPT(4)*X)/(10.0*SPT(4))
Z2=Z1*Z1+2.0*SPT(1)/(5.0*SPT(4)*X)
IF(Z2.LT.0.0) GOTO 300
Z3=SQRT(Z2)
Z2=-1.0*Z1+Z3
IF(Z2.GT.0.0.AND.Z2.LT.Y) AUS6=.TRUE.
Z2=-1.0*Z1-Z3
IF(Z2.GT.0.0.AND.Z2.LT.Y) AUS6=.TRUE.
300 RETURN
END
REAL FUNCTION XHIS(XMONTH,XHOUR,XLATI)
C CALCULATES THE SOLAR-ZENITH-ANGLE
DAYNR=XMONTH*30.0-15.0
SUNDEC=-0.40915*COS(1.7202E-2*(DAYNR+8.0))
UMR=1.7453E-2
COSXHI=SIN(XLATI*UMR)*SIN(SUNDEC) +COS(XLATI*UMR)*COS(SUNDEC)*
1COS(2.6180E-1*(XHOUR-12.0))
XHIS=ACOS(COSXHI)/UMR
RETURN
END
SUBROUTINE GGM(ART,XLG,BG,XLM,BM)
C CALCULATES GEOMAGNETIC COORDINATES (XLM,BM) FROM GEOGRAFIC COORDINATES
C (XLG,BG) FOR ART=0 AND REVERSE FOR ART=1
INTEGER ART
ZPI=6.2831853
FAKTOR=0.0174533
CBG=11.4*FAKTOR
CI=COS(CBG)
SI=SIN(CBG)
IF(ART.EQ.0) GOTO 10
CBM=COS(BM*FAKTOR)
SBM=SIN(BM*FAKTOR)
CLM=COS(XLM*FAKTOR)
SLM=SIN(XLM*FAKTOR)
SBG=SBM*CI-CBM*CLM*SI
BG=ASIN(SBG)
CBG=COS(BG)
SLG=(CBM*SLM)/CBG
CLG=(SBM*SI+CBM*CLM*CI)/CBG
XLG=ACOS(CLG)
IF(SLG.LT.0.0) XLG=ZPI-ACOS(CLG)
BG=BG/FAKTOR
XLG=XLG/FAKTOR
XLG=XLG-69.8
IF(XLG.LT.0.0) XLG=XLG+360.0
GOTO 20
10 YLG=XLG+69.8
CBG=COS(BG*FAKTOR)
SBG=SIN(BG*FAKTOR)
CLG=COS(YLG*FAKTOR)
SLG=SIN(YLG*FAKTOR)
SBM=SBG*CI+CBG*CLG*SI
BM=ASIN(SBM)
CBM=COS(BM)
SLM=(CBG*SLG)/CBM
CLM=(-SBG*SI+CBG*CLG*CI)/CBM
XLM=ACOS(CLM)
IF(SLM.LT.0.0) XLM=ZPI-ACOS(CLM)
BM=BM/FAKTOR
XLM=XLM/FAKTOR
20 RETURN
END
REAL FUNCTION SSRSU(HO)
REAL LATI,MONTH
COMMON/BLO11/MONTH,LATI
SSRSU=XHIS(MONTH,HO,LATI)
RETURN
END
SUBROUTINE KOEFFB(FIELD)
C TRANSFORMATION COEFFICIENTS G(1@144) VALID FOR 1973 FOR CALCULATING
C THE MAGNETIC FIELD ACCORDING TO POGO 68/10
DIMENSION FIELD (144)
REAL FELD (144)
DATA (FELD(KK),KK=1,80)/0.0, 0.1506723,0.0101742, -0.0286519, 0.009
12606, -0.0130846, 0.0089594, -0.0136808,-0.0001508, -0.0093977,
2 0.0130650, 0.0020520, -0.0121956, -0.0023451, -0.0208555,
3 0.0068416,-0.0142659, -0.0093322, -0.0021364, -0.0078910,
4 0.0045586, 0.0128904, -0.0002951, -0.0237245,0.0289493,
5 0.0074605, -0.0105741, -0.0005116, -0.0105732, -0.0058542,
60.0033268, 0.0078164,0.0211234, 0.0099309, 0.0362792, -0.0201070,
7 -0.0046350, -0.0058722, 0.0011147, -0.0013949,
8 -0.0108838, 0.0322263, -0.0147390, 0.0031247, 0.0111986,
9 -0.0109394,0.0058112, 0.2739046, -0.0155682, -0.0253272,
1 0.0163782, 0.0205730, 0.0022081, 0.0112749,-0.0098427,
2 0.0072705, 0.0195189, -0.0081132, -0.0071889, -0.0579970,
3 -0.0856642, 0.1884260,-0.7391512, 0.1210288, -0.0241888,
4 -0.0052464, -0.0096312, -0.0044834, 0.0201764, 0.0258343,
50.0083033, 0.0077187,0.0586055,0.0102236, -0.0396107,
6 -0.0167860, -0.2019911, -0.5810815,0.0379916, 3.7508268/
DATA (FELD(KZ),KZ=81,144)/1.8133030,
7 -0.0564250, -0.0557352, 0.1335347, -0.0142641,
8 -0.1024618,0.0970994, -0.0751830,-0.1274948, 0.0402073,
9 0.0386290, 0.1883088, 0.1838960, -0.7848989,0.7591817,
1 -0.9302389,-0.8560960, 0.6633250, -4.6363869, -13.2599277,
2 0.1002136, 0.0855714,-0.0991981, -0.0765378,-0.0455264,
3 0.1169326, -0.2604067, 0.1800076, -0.2223685, -0.6347679,
40.5334222, -0.3459502,-0.1573697, 0.8589464, 1.7815990,
5-6.3347645, -3.1513653, -9.9927750,13.3327637, -35.4897308,
637.3466339, -0.5257398, 0.0571474, -0.5421217, 0.2404770,
7 -0.1747774,-0.3433644, 0.4829708,0.3935944, 0.4885033,
8 0.8488121, -0.7640999, -1.8884945, 3.2930784,-7.3497229,
9 0.1672821,-0.2306652, 10.5782146, 12.6031065, 8.6579742,
1 215.5209961, -27.1419220,22.3405762,1108.6394043/
K=0
DO 10 I=1,144
K=K+1
10 FIELD(K)=FELD(I)
RETURN
END
SUBROUTINE FIELDG(DLAT,DLONG,ALT,X,Y,Z,F,DIP,SMODIP)
C IPW,28.12.77,THIS IS A SPECIAL VERSION OF THE POGO 68/10 MODEL
C TO ECONOMIZE COMPUTER TIME AND MEMORY.X,Y,Z,F,DIP ARE THE COORDINATES
C MAGNETUDE (GAUSS) AND ANGLE (VERSUS THE HORICONTAL AXIS) OF THE
C EARTH MAGNETIC FIELD IN THE HEIGHT ALT/KM AT THE LATITUDE
C DLAT AND THE LONGITUDE DLONG
C "COMMENT" THIS IS A SPECIAL VERSION OF THE POGO68/107
C MAGNETIC FIELD LEGENDRE MODEL (IN GAUSS UNITS).
C F IS TOTAL FIELD, Z THE (DOWNWARD) VERTICAL COMPONENT WHILE
C X AND Y ARE COMPONENTS IN THE EQUATORIAL PLANE X TO ZERO LONGITUDE).
C DIP=INCLINATION ANGLE/DEGREE. MODIP=RAWER"S MODFIED DIP.
C REFERENCE@ METEOROLOGICAL AND ASTRONOMICAL INFLUENCES ON
C RADIO WAVE PROPAGATION. PERGAMON PRESS, 1963
C INPUT@ DLAT, DLONG=GEOGRAPHIC COORDINATES/DEGREE, ALT=ALTITUDE/KM.
C THE SET OF COEFFICIENTS IS GIVEN IN THE FOLLOWING PROCEDURE?
DIMENSION H(144),XI(3)
COMMON/BLOCK2/B0,B1,C1,HZ,T,G(144),HST,STR
RLAT=DLAT*0.0174533
CT=SIN(RLAT)
ST=COS(RLAT)
NMAX=11
D=SQRT(40680925.0-272336.0*CT*CT)
RLONG=DLONG*0.0174533
CP=COS(RLONG)
SP=SIN(RLONG)
ZZZ=(ALT+40408589.0/D)*CT/6371.2
RHO=(ALT+40680925.0/D)*ST/6371.2
XXX=RHO*CP
YYY=RHO*SP
RQ=1.0/(XXX*XXX+YYY*YYY+ZZZ*ZZZ)
XI(1)=XXX*RQ
XI(2)=YYY*RQ
XI(3)=ZZZ*RQ
IHMAX=NMAX*NMAX+1
LAST=IHMAX+NMAX+NMAX
IMAX=NMAX+NMAX-1
DO 100 I=IHMAX,LAST
100 H(I)=G(I)
DO 200 K=1,3,2
I=IMAX
IH=IHMAX
300 IL=IH-I
F1=2.0/FLOAT(I-K+2)
X1=XI(1)*F1
Y1=XI(2)*F1
Z1=XI(3)*(F1+F1)
I=I-2
IF((I-1).LT.0) GOTO 400
IF((I-1).EQ.0) GOTO 500
DO 600 M=3,I,2
H(IL+M+1)=G(IL+M+1)+Z1*H(IH+M+1)+X1*(H(IH+M+3)-H(IH+M-1))-Y1*
1(H(IH+M+2)+H(IH+M-2))
H(IL+M)=G(IL+M)+Z1*H(IH+M)+X1*(H(IH+M+2)-H(IH+M-2))+Y1*(H(IH+
1M+3)+H(IH+M-1))
600 CONTINUE
500 H(IL+2)=G(IL+2)+Z1*H(IH+2)+X1*H(IH+4)-Y1*(H(IH+3)+H(IH))
H(IL+1)=G(IL+1)+Z1*H(IH+1)+Y1*H(IH+4)+X1*(H(IH+3)-H(IH))
400 H(IL)=G(IL)+Z1*H(IH)+2.0*(X1*H(IH+1)+Y1*H(IH+2))
700 IH=IL
IF(I.GE.K) GOTO 300
200 CONTINUE
S=0.5*H(1)+2.0*(H(2)*XI(3)+H(3)*XI(1)+H(4)*XI(2))
XT=(RQ+RQ)*SQRT(RQ)
X=XT*(H(3)-S*XXX)
Y=XT*(H(4)-S*YYY)
Z=XT*(H(2)-S*ZZZ)
F=SQRT(X*X+Y*Y+Z*Z)
Z=-Z*CT-(Y*SP+X*CP)*ST
DIP=ASIN(Z/F)
SMODIP=ASIN(DIP/SQRT(DIP*DIP+ST))
DIP=DIP*57.2957795
9 SMODIP=SMODIP*57.2957795
RETURN
END
SUBROUTINE F2OUT(XMODIP,XLATI,XLONGI,XMAGBR,R,XMONTH,HOUR,FOE,
1XHMF2,FOF2)
C D.BILITZA,20.8.78,CALCULATES XHMF2 AND NMF2 BY USING THE
C CCIR-MAPS FOR FOF2 AND M3000.YOU HAVE TO ADJUST THE TAPE
C READING PROCEDURE CCIRCA TO YOUR COMPUTER SYSTEM AND CCIR-TAPE
INTEGER QF(9),QM(7)
DIMENSION FF0(988)
REAL MM0(441)
REWIND 5
LMONTH=IFIX(XMONTH)
CALL CCIRCA(R,LMONTH,FF0,MM0)
QF(1)=11
QF(2)=11
QF(3)=8
QF(4)=4
QF(5)=1
QF(6)=0
QF(7)=0
QF(8)=0
QF(9)=0
QM(1)=6
QM(2)=7
QM(3)=5
QM(4)=2
QM(5)=1
QM(6)=0
QM(7)=0
FOF2=GAMMA1(XMODIP,XLATI,XLONGI,.FALSE.,HOUR,6,QF,9,76,13,988,
1FF0)
XM3000=GAMMA1(XMODIP,XLATI,XLONGI,.FALSE.,HOUR,4,QM,7,49,9,441,
1MM0)
XHMF2=HMF2ED(XMAGBR,R,FOF2/FOE,XM3000)
RETURN
END
REAL FUNCTION GAMMA1(SMODIP,SLAT,SLONG,UT,XHOUR,IHARM,NQ,K1,M,
1MM,M3,SFE)
C SHEIKH,4.3.77,CALCULATES A QUANTITY (FOF2/M3000) USING CCIR
C NUMERICAL MAP COEFFICIENTS SFE(M,2*IHARM+1).
C NQ(K1) IS AN INTEGER ARRAY GIVING NQ1,NQ2,NQ3,...AS HIGHEST
C DEGREE OF LATITUDE FOR EACH LONGITUDE HARMONIC (IHARM).
C M=1+NQ1+2(NQ2+1)+2(NQ3+1)+... .IF UT IS TRUE,THEN HOUR IS
C TAKEN AS UNIVERSAL TIME,OTHERWISE IT IS ASSUMED AS LOCAL
C TIME.
DIMENSION NQ(K1),XSINX(13),COEF(100),C(12),S(12),SFE(M3)
LOGICAL UT
RD=57.2957795
IF(UT) GOTO 100
TIME=(15.0*XHOUR-180.0-SLONG)/RD
GOTO 150
100 TIME=(15.0*XHOUR-180.0)/RD
150 S(1)=SIN(TIME)
C(1)=COS(TIME)
DO 250 I=2,IHARM
C(I)=C(1)*C(I-1)-S(1)*S(I-1)
S(I)=C(1)*S(I-1)+S(1)*C(I-1)
250 CONTINUE
DO 300 I=1,M
COEF(I)=SFE((I-1)*MM+1)
DO 300 J=1,IHARM
COEF(I)=COEF(I)+SFE((I-1)*MM+2*J)*S(J)+SFE((I-1)*MM+2*J+1)*C(J)
300 CONTINUE
SUM=COEF(1)
SS=SIN(SMODIP/RD)
S3=SS
XSINX(1)=1.0
INDEX=NQ(1)
DO 350 J=1,INDEX
SUM=SUM+COEF(1+J)*SS
XSINX(J+1)=SS
SS=SS*S3
350 CONTINUE
XSINX(NQ(1)+2)=SS*S3
NP=NQ(1)+1
SS=COS(SLAT/RD)
S3=SS
DO 400 J=2,K1
S0=SLONG*FLOAT(J-1)/RD
S1=COS(S0)
S2=SIN(S0)
INDEX=NQ(J)+1
DO 450 L=1,INDEX
NP=NP+1
SUM=SUM+COEF(NP)*XSINX(L)*SS*S1
NP=NP+1
SUM=SUM+COEF(NP)*XSINX(L)*SS*S2
450 CONTINUE
SS=SS*S3
400 CONTINUE
GAMMA1=SUM
RETURN
END
SUBROUTINE CCIRCA(R,IMONTH,FOF2,SM3000)
C THIS SUBROUTINE CALCULATES THE COEFFICIENTS ARRAYS FOF2(76,13) AND
C SM3000(49,9) FOR SOLARACTIVITY (R) AND MONTH USING OUR SPECIAL CCIRFO-TAPE
DIMENSION FOF2(988),SM3000(441),F2(13,76,2),FM3(9,49,2)
DO 10 I=1,IMONTH
READ(5,1)
1 FORMAT( )
READ(5,2) ((F2(L,M,1),M=1,76),L=1,13)
READ(5,1)
READ(5,2) ((F2(L,M,2),M=1,76),L=1,13)
READ(5,1)
READ(5,2) ((FM3(L,M,1),M=1,49),L=1,9)
READ(5,1)
READ(5,2) ((FM3(L,M,2),M=1,49),L=1,9)
2 FORMAT(4(5X,E13.7))
10 CONTINUE
DO 20 I=1,76
DO 20 J=1,13
K=J+13*(I-1)
20 FOF2(K)=(F2(J,I,1)*(100.0-R)+F2(J,I,2)*R)/100.0
DO 30 I=1,49
DO 30 J=1,9
K=J+9*(I-1)
30 SM3000(K)=(FM3(J,I,1)*(100.0-R)+FM3(J,I,2)*R)/100.0
RETURN
END
SUBROUTINE KOEFFI(BOF)
C COEFFICIENTS FOR BOTTOMSIDE F2 REGION ELECTRON DENSITY
DIMENSION BOF(2,2,8)
REAL FELD(32)
DATA FELD/114.0,64.0,134.0,77.0,128.0,66.0,75.0,73.0,113.0,115.0,
1150.0,116.0,138.0,123.0,94.0,132.0,72.0,84.0,83.0,89.0,75.0,85.0,
257.0,76.0,102.0,100.0,120.0,110.0,107.0,103.0,76.0,86.0/
L=0
K=0
DO 10 I=1,2
DO 10 J=1,2
DO 10 K=1,8
L=L+1
10 BOF(I,J,K)=FELD(L)
RETURN
END
SUBROUTINE KOEFP1(PG1O)
C COEFFICIENTS PG1O FOR CALCULATING O+ PROFILES
DIMENSION PG1O(80)
REAL FELD (80)
DATA FELD/-11.0,-11.0,4.0,-11.0,0.08018,
60.13027,0.04216,0.25 ,-0.00686,0.00999,
15.113,0.1 ,170.0,180.0,0.1175,0.15,-11.0,
71.0 ,2.0,-11.0,0.069,0.161,0.254,0.18,0.0161,
20.0216,0.03014,0.1,152.0,167.0,0.04916,
80.17,-11.0,2.0,2.0,-11.0,0.072,0.092,0.014,0.21,
30.01389,0.03863,0.05762,0.12,165.0,168.0,0.008,
90.258,-11.0,1.0,3.0,-11.0,0.091,0.088,
40.008,0.34,0.0067,0.0195,0.04,0.1,158.0,172.0,
10.01,0.24,-11.0,2.0,3.0, -11.0,0.083,0.102,
50.045,0.03,0.00127,0.01,0.05,0.09,167.0,185.0,
20.015,0.18/
K=0
DO 10 I=1,80
K=K+1
10 PG1O(K)=FELD(I)
RETURN
END
SUBROUTINE KOEFP2(PG2O)
C COEFFICIENTS FOR CALCULATING O+ PROFILES
DIMENSION PG2O(32)
REAL FELD(32)
DATA FELD/1.0,-11.0,-11.0,1.0,695.0,-.000781,
3-.00264,2177.0,1.0,-11.0,-11.0,2.0,570.0,
1-.002,-.0052,1040.0,2.0,-11.0,-11.0,1.0,695.0,
4-.000786,-.00165,3367.0,2.0,-11.0,-11.0,2.0,
2575.0,-.00126,-.00524,1380.0/
K=0
DO 10 I=1,32
K=K+1
10 PG2O(K)=FELD(I)
RETURN
END
SUBROUTINE KOEFP3(PG3O)
C COEFFICIENTS FOR CALCULATING O2+ PROFILES
DIMENSION PG3O(80)
REAL FELD(80)
DATA FELD/-11.0,1.0,2.0,-11.0,160.0,31.0,130.0,
9-10.0,198.0,0.0,0.05922,-0.07983,
1-0.00397,0.00085,-0.00313,0.0,-11.0,2.0,2.0,-11.0,
8140.0,30.0,130.0,-10.0,
2190.0,0.0,0.05107,-0.07964,0.00097,-0.01118,-0.02614,
7-0.09537,
3-11.0,1.0,3.0,-11.0,140.0,37.0,125.0,0.0,182.0,
60.0,0.0307,-0.04968,-0.00248,
4-0.02451,-0.00313,0.0,-11.0,2.0,3.0,-11.0,
5140.0,37.0,125.0,0.0,170.0,0.0,
50.02806,-0.04716,0.00066,-0.02763,-0.02247,-0.01919,
6-11.0,-11.0,4.0,-11.0,140.0,45.0,136.0,-9.0,
7181.0,-26.0,0.02994,-0.04879,
7-0.01396,0.00089,-0.09929,0.05589/
K=0
DO 10 I=1,80
K=K+1
10 PG3O(K)=FELD(I)
RETURN
END
SUBROUTINE SUFE (FIELD,RFE,M,FE)
C SPECIAL SUBROUTINE FOR CHOOSING THE DESIRED
C ION DENSITY PARAMETER ARRAYS
DIMENSION RFE(4),FE(12),FIELD(80),EFE(4)
K=0
100 DO 101 I=1,4
K=K+1
101 EFE(I)=FIELD(K)
DO 111 I=1,M
K=K+1
111 FE(I)=FIELD(K)
DO 120 I=1,4
IF((EFE(I).GT.-10.0).AND.(RFE(I).NE.EFE(I))) GOTO 100
120 CONTINUE
RETURN
END
FUNCTION HPOL(XHOUR,TW,XNW,SA,SU)
C SPECIAL FUNCTION FOR TIME-INTERPOLATION.TW,XNW ARE THE DAY
C AND NIGHT VALUE OF THE PARAMETER.SA,SU ARE TIME (HOURS) OF
C SUNRISE AND SUNSET.
C "COMMENT" D.BILITZA,1978 SPECIAL PROCEDURE FOR TIME-INTERPOLATION
C USING EPSTEIN STEP FUNKTION OF ONE HOUR WIDTH (DAY/NIGHT FUNCTION).
C TW,NW ARE THE DAY AND NIGHT VALUE OF THE PARAMETER.
C SA, SU ARE TIME (HOURS) OF SUNRISE AND SUNSET/HOUR?
IF(ABS(SU-SA).GE.24.0) GOTO 100
HPOL=XNW+(TW-XNW)/(1.0+EXP(-(XHOUR-SA)/1.0))+(XNW-TW)
1/(1.0+EXP(-(XHOUR-SU)/1.0))
GOTO 200
100 HPOL=XNW
IF(SA.LT.1.0) HPOL=TW
200 RETURN
END
REAL FUNCTION XPOL(AA,XHOUR,SA,SU,DELL,IMONTH,R)
C SPECIAL FUNCTION TO INTERPOLATE ARRAY AA(LATI,R,MONTH,HOUR)=
C A(2,2,4,2) IN LATITUDE (LATI), SOLAR ACTIVITY (R) AND TIME (HOUR)
DIMENSION AA(2,2,8),SIPH(2),SIPL(2)
JMONTH=IMONTH*2
DO 10 ISR=1,2
DO 15 ISLAT=1,2
15 SIPH(ISLAT)=HPOL(XHOUR,AA(ISLAT,ISR,JMONTH-1),AA(ISLAT,ISR,
1JMONTH),SA,SU)
SIPL(ISR)=SIPH(1)+(SIPH(2)-SIPH(1))/DELL
10 CONTINUE
XPOL=SIPL(1)+(SIPL(2)-SIPL(1))/90.0*(R-10.0)
RETURN
END
FUNCTION EPSTEP(BE,AB,TH,FIX,COSI)
C EPSTEIN STEP FUNCTION = TRANSITION WITH THICKNESS TH BETWEEN
C BE AND AB CENTERED WHERE VARIABLE COSI EQUALS FIX
EPSTEP=AB+(BE-AB)/(1.0+EXP(-(COSI-FIX)/TH))
RETURN
END