C IRIFO7, DECEMBER 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 IF XHI LESS -10.0 THEN XHI WILL BE CALCULATED.
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.
INTEGER EGNR,AGNR,SMONTH,DAYNR,DDO,IIF,DO2
REAL LATI,LONGI,MO2,MO,LOGE,MONTH,MODIP,NMF2,M3000,MAGBR,
&NDEL0,NMF1,MUE,NME,NMD,K,NHABR,NDEL,NDX,NE0,NEI,MM,MLAT,
&MLONG,NOBO2
DIMENSION F(3),B0F(2,2,8),OUTF(50,12),UF(38),RIF(4),XSM(2),
&FF(38),FA(3,11),FZ(3,11),JF(10),DDO(4),IIF(4),DO2(2),MM(3),
&Q(3),QT(3),PF1O(12),PF3O(12),HO(4),MO2(3),MO(5),G(144),E(4),
&PF2O(4),HO2(2),CTNN(3),PG1O(80),PG2O(32),PG3O(80),CTN(3)
LOGICAL WINTER,SUMMER,SCHALT,EXT,NIGHT,F1REG,VALLEY,
&IRDLOW,IRDUPP,CCIREI,ANF
COMMON/BLOCK1/HMF2,NMF2,HMF1/BLOCK2/B0,B1,C1,HZ,T,HST,STR
&/BLOCK3/HDX,HME,NME,HMD,NMD,HEF,D1,K,FP30,FP3U,FP1,FP2
COMMON/BLOCK1/HMF2,NMF2,HMF1/BLOCK2/B0,B1,C1,HZ,T,HST,STR
&/BLOCK3/HDX,HME,NME,HMD,NMD,HEF,D1,K,FP30,FP3U,FP1,FP2
&/BLOCK4/G/BLOCK5/ZX,TNX,ATN,CTN/BLOCK6/NIGHT,E
&/BLOCK8/HS,TNHS,XSM,MM/BLO10/BETA,ETA,DELTA,ZETA
EXTERNAL XE1,XE2,XE3,XE4,XE5,XE6,TEDER
DATA FF,UF/4H(1X,,4HF6.1,1H,,35*1H ,4H(1H0,4H,3X,,4H1HH,,
&3H2X,,34*1H /,FA/4H3X,E,4H10.4,1H,,4H3X,F,4H6.4,,1H ,
&4H3X,F,4H6.1,,1H ,4H3X,F,4H6.1,,1H ,4H3X,F,4H6.1,,1H ,
&4H3X,F,4H7.4,,1H ,4H3X,F,4H5.1,,1H ,4H3X,F,4H5.1,,1H ,
&4H3X,F,4H5.1,,1H ,4H3X,F,4H5.1,,1H ,4H3X,F,4H5.1,,1H /,FZ/
&4H7X,2,4HHNE,,3H3X,,4H4X,6,4HHN/N,4HMAX,,4H5X,2,4HHTN,,3H2X,,
&4H5X,2,4HHTE,,3H2X,,4H5X,2,4HHTI,,3H2X,,4H4X,5,4HHTE/,3HTI,,
&4H5X,4,4HHRDO,2H+,,4H4X,4,4HHRDH,2H+,,4H3X,5,4HHRDH,3HE+,,
&4H3X,5,4HHRDO,3H2+,,4H3X,5,4HHRDN,3H0+,/,ENFF,ENUF/2*1H)/
LOGE=1/ALOG(10.0)
UMR=0.01745329
PHI=3.1415927
CCIREI=.FALSE.
C FIRST SPECIFY YOUR COMPUTERS CHANNEL NUMBERS,
C EGNR=INPUT, AGNR=OUTPUT, KONSOL=MONITOR OUTPUT
EGNR=5
AGNR=6
KONSOL=6
READ(EGNR,2) LATI,LONGI,R,MONTH,HOUR,XHI,AH,EH,SH,KOBE,(JF(I),
&I=1,10),JMAG,HMF2,FOF2
WRITE(KONSOL,2122) LATI,LONGI,R,MONTH,HOUR,XHI,AH,EH,SH,
&KOBE,(JF(I),I=1,10),JMAG,HMF2,FOF2
2 FORMAT(3F6.1,2F4.1,F6.1,1X,3(F6.1,1X),I1,1X,11I1/F5.1,E9.3)
2122 FORMAT(1H*,3F6.1,2F4.1,F6.1,1X,3(F6.1,1X),I1,1X,11I1/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-ARCOS(-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=Z1+Z2*COS(XLSTA*UMR)
XHI=ARCOS(COSXHI)/UMR
GOTO 503
520 COSXHI=COS(XHI*UMR)
503 XHI0=49.84733+0.349504*ABS(LATI)
XHI100=38.96113+0.509932*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.
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).OR.NIGHT.OR.WINTER) 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)+
&(0.0024287+0.0042810*COS2-0.00015280*FOF2)*FOF2
ZETA=0.078922-0.0046702*COS2+FLU*(-0.019132+0.0076545*COS2)+
&(0.0032513+0.0060290*COS2-0.00020872*FOF2)*FOF2
BETA=-128.03+20.253*COS2+FLU*(-8.0755-0.65896*COS2)+(0.44041
&+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
& 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
924 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
& 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 IF (HMF2-HEF.GE.B0) GOTO 922
B0=(HMF2-HEF)/1.1
WRITE(KONSOL,923) B0
923 FORMAT(1X,"CORRECTION BY DECREASING OF B0 TO",F7.1)
GOTO 924
922 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
&EVALUATED 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
& 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,THE RATIO OF NO+ TO O2+ AT HEIGHT OF REL.
C O+ MAX.(H0O).SAME RATIO IS CONSERVED AT ALL HEIGHTS ABOVE
IF(JF(7).LE.0) GOTO 7114
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=",F6.1)
WRITE(AGNR,7066) SAX,SUX,SUNDEC
7066 FORMAT(1X,"SUNRISE",F4.1,"L.T. SUNSET@",F4.1,"L.T.",
&5X,"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)
KARL=4
DO 7091 J=1,7
IF(JF(J).LT.1) GOTO 7091
DO 7092 I=1,3
FF(KARL+I)=FA(I,J)
7092 UF(KARL+I)=FZ(I,J)
KARL=KARL+3
7091 CONTINUE
J=8
IF(JF(J).LT.1) GOTO 7093
DO 7094 I=1,3
FF(KARL+I)=FA(I,J)
UF(KARL+I)=FZ(I,J)
II=I+3
JJ=J+1
FF(KARL+II)=FA(I,JJ)
7094 UF(KARL+II)=FZ(I,JJ)
KARL=KARL+6
7093 CONTINUE
DO 7095 J=9, 10
IF(JF(J).LT.1) GOTO 7095
JJ=J+1
DO 7096 I=1,3
FF(KARL+I)=FA(I,JJ)
7096 UF(KARL+I)=FZ(I,JJ)
KARL=KARL+3
7095 CONTINUE
FF(KARL+1)=ENFF
UF(KARL+1)=ENUF
WRITE(AGNR,UF)
7041 DO 7040 I=1,IEI
7040 WRITE(AGNR,FF) (OUTF(I,L),L=1,IN)
IF(KOBE.EQ.1) WRITE(AGNR,35)
35 FORMAT(1X,"111111")
STOP
END
C D.BILITZA,H.THIEMANN,K.RAWER,IPW FREIBURG,
C AUG.79 **********************
C IRI-PROCEDURES ***************************************
C \REF.@ K. RAWER, S. RAMAKRISHNAN, D. BILITZA,
C INTERNATIONAL REFERENCE IONOSPHERE 1978, U.R.S.I. BRUSSELS'
C FOR CALCULATION OF ELECTRON DENSITY, TEMPERATURE, ION TEMPERATURE
C AND ION RELATIVE DENSITY IN THE HEIGHT RANGE 70 - 1000 KM.?
C I. ELECTRON DENSITY ----------------------------------
FUNCTION EP1ST(X,BET)
C EPSTEIN-STEP-FUNCTION WITH TRANSITION FROM Y=0 TO Y=1 AT X=0
C AND TRANSITION THICKNESS 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 REPRESENTING ELECTRON DENSITY
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 GEOM. LATITUDE
C (MLAT), SOLAR FLUX(FLU) AND CRITICAL FREQUENCY(FOF2).
C ALO GLOBAL ARE PEAK-PARAMETERS NMF2,HMF2?
C \REF.@K.RAWER,S.RAMAKRISHNAN,1978'
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.0+EXP((-94.5-DELTA)/BETA)))+ZETA*(100.0*ALOG((1.0+
&EXP((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,HST,STR
X=(HMF2-H)/B0
IF(ABS(X).LT.1.0E-10) GOTO 100
XE2=NMF2*EXP(-X**B1)/COSH(X)
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,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,HST,STR/BLOCK3/
&HDX,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)
REAL NME,NMD,K
LOGICAL NIGHT
COMMON/BLOCK3/HDX,HME,NME,HMD,NMD,HEF,D1,K,FP30,FP3U,FP1,
&FP2/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
C SUMMARIZING PROCEDURES NE1....6?
REAL NMF2,NME,NMD,K
COMMON/BLOCK1/HMF2,NMF2,HMF1/BLOCK2/B0,B1,C1,HZ,T,HST,STR
&/BLOCK3/HDX,HME,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 1978, CALCULATES THE HEIGHT DERIVATIVE
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 1978, CALCULATES THE PARAMETERS F1 AND F2
C NEEDED IN THE ELECTRON TEMPERATURE MODEL PROCEDURE TE,
C BOTH DEPENDING ON GEOMAGNETIC LATITUDE (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 1978. ELECTRON TEMPERATURE PROFILE
C IN THE HEIGHT RANGE 120-1000 KM. F1,F2 ARE LATITUDE DEPENDENT AND
C ARE GIVEN BY TELAT. HMAX IS THE HEIGHT OF A POSSIBLE MAXIMUM
C D THE TEMPERATURE HEIGHT GRADIENT, H0 THE INTERSECTION HEIGHT
C BELOW WHICH THE PROFILE IS BUILT TO APPROACH THE CIRA 72 TEMP.
C 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 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 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))/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))/G(I)))
100 SUM=SUM+(MM(I+1)-MM(I))*(A-B)
TI=SUM
RETURN
END
REAL FUNCTION TEDER(H)
C THIS FUNCTION ALONG WITH PROCEDURE REGFA1 ALLOWS TO FIND
C THE HEIGHT ABOVE WHICH TN BEGINS TO BE DIFFERENT 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 D.BILITZA,1977,THIS ANALYTIC FUNCTION IS USED TO REPRESENT THE
C RELATIVE PRECENTAGE DENSITY OF ATOMAR AND MOLECULAR OXYGEN IONS.
C THE M+1 HEIGHT GRADIENTS ST(M+1) ARE CONNECTED WITH EPSTEIN-
C STEP-FUNCTIONS AT THE STEP HEIGHTS XS(M) WITH TRANSITION
C THICKNESSES ID(M). RPID(H0,H0,N0,....)=N0.
C INSTEAD OF 88.0 YOU CAN USE THE HIGHEST ALLOWED ARGUMENT
C FOR EXP AT YOUR COMPUTER.
REAL N0
DIMENSION ID(4), ST(5), XS(4)
SUM=(H-H0)*ST(1)
DO 100 I=1,M
A=H-XS(I)
XI=FLOAT(ID(I))
IF (A.LT.88.0*XI) A=XI*ALOG(1.0+EXP((H-XS(I))/XI))
B=H0-XS(I)
IF (B.LT.88.0*XI) B=XI*ALOG(1.0+EXP((H0-XS(I))/XI))
100 SUM=SUM+(ST(I+1)-ST(I))*(A-B)
SUM=N0*EXP(SUM)
IF (SUM.LT.0.00005) 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 PERECENTAGEDENSITY FOR HEIGHTS
C NOT GREATER THAN 1000 KM.THE O+ AND O2+ REL. PER. DENSITIES SHOULD
C BE GIVEN (RDOH,RDO2H). HB IS THE HEIGHT WERE N(O+)=100_.
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.BILITZA, 1978. NO+ RELATIVE PERCENTAGE DENSITY FOR HEIGHTS
C NOT LESS THAN 100 KM.FOR MORE INFORMATION SEE SUBROUTINE RDHHE.
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 D. BILITZA, 1978, CALCULATES ELECTRON DENSITY OF D MAXIMUM.
C XHI IS SOLAR ZENITH ANGLE, R ZURICH SUNSPOT NUMBER AND XW THE
C REQUIRED NIGHT VALUE \SOURCE@MECHTLY-BILITZA,1974'.
Y=6.05+0.088*R
YW=XW/1.0E8
Z=(-0.1/(ALOG(YW/Y)))**0.3704
SUXHI=ARCOS(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 D.BILITZA, 1977, CALCULATES FOE BY THE EDINBURGH-METHOD.
C INPUT@ZUERICH SUNSPOT NUMBER(R),SOLAR ZENITH ANGLE AT LOCAL
C HOUR ST (XHI) AND AT NOON (XHIM), LATITUDE (SLATI),SUNSET
C AND SUNRISE (SU, SA/HOUR).
C \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,1977, CALCULATES FOF1 FOR DIP- LATITUDE YLATI ,
C ZURICH SUNSPOT NUMBER R AND COSINUS OF SOLAR ZENITH
C ANGLE COSXHI. \REF@E.D.DUCHARME ET.AL.,RADIO SCI.,8,10,837-839,
C HOWEVER WITH MAG. INSTEAD OF GEOM. LATITUDE,R.EYFRIG,1979'
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 THE PEAK HEIGHT HMF2 FOR THE MAGNETIC
C LATITUDE XMAGBR AND THE ZUERICH SUNSPOT NUMBER R
C USING M3000(XM3) AND THE RATIO FOF2/FOE.
C \REF. D. BILITZA ET. AL. TELECOMM.J.,46,549-553,1979'
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
C STARTING VALUES. THE COMUTATION ENDS WHEN THE X-INTERVAL
C HAS BECOME LESS THAN EPS . IF SIGN(F(X1)-FW)= SIGN(F(X2)-FW)
C THEN SCHALT=.TRUE.
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,CAL. OF THE COEF. SPT(4) FOR A POLYNOM OF
C 5TH DEGREE (SPT(0)=1) TO FIT THE VALLEY REPR. BY THE PARAMETER@
C SHMAX IS THE HEIGHT OF THE LOWER MAXIMUM WHERE THE VALLEY STARTS,
C STMAX IS THE DENSITY AT SHMAX, SHABR IS THE HEIGHT DIFFERENCE
C BETWEEN THE HEIGHT OF THE DEEPEST VALLEY DENSITY AND SHMAX, SDELTA
C IS THE PERCENTAGE DEPTH,SHBR THE WIDTH,STE0 THE UPPER FIXPOINT
C NORMALISED TO STMAX AND SDTDH0 IS THE LOGARITHMIC DERIV. 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)*
&COS(2.6180E-1*(XHOUR-12.0))
XHIS=ARCOS(COSXHI)/UMR
RETURN
END
SUBROUTINE GGM(ART,XLG,BG,XLM,BM)
C CALCULATES GEOM. COORDINATES (XLM,BM) FROM GEOG. 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=ARSIN(SBG)
CBG=COS(BG)
SLG=(CBM*SLM)/CBG
CLG=(SBM*SI+CBM*CLM*CI)/CBG
XLG=ARCOS(CLG)
IF(SLG.LT.0.0) XLG=ZPI-ARCOS(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=ARSIN(SBM)
CBM=COS(BM)
SLM=(CBG*SLG)/CBM
CLM=(-SBG*SI+CBG*CLG*CI)/CBM
XLM=ARCOS(CLM)
IF(SLM.LT.0.0) XLM=ZPI-ARCOS(CLM)
BM=BM/FAKTOR
XLM=XLM/FAKTOR
20 RETURN
END
SUBROUTINE KOEFFB(FIELD)
C SHEIK,1977,TRANSFORMATION COEFF. G(144) VALID FOR 1973 FOR
C CALCULATING THE MAGNETIC FIELD ACCORDING TO POGO 68/10
DIMENSION FIELD (144)
REAL FEL1 (72),FEL2(72)
DATA FEL1/0.0, 0.1506723,0.0101742, -0.0286519, 0.0092606,
& -0.0130846, 0.0089594, -0.0136808,-0.0001508, -0.0093977,
& 0.0130650, 0.0020520, -0.0121956, -0.0023451, -0.0208555,
& 0.0068416,-0.0142659, -0.0093322, -0.0021364, -0.0078910,
& 0.0045586, 0.0128904, -0.0002951, -0.0237245,0.0289493,
& 0.0074605, -0.0105741, -0.0005116, -0.0105732, -0.0058542,
&0.0033268, 0.0078164,0.0211234, 0.0099309, 0.0362792,
&-0.0201070,-0.0046350,-0.0058722,0.0011147,-0.0013949,
& -0.0108838, 0.0322263, -0.0147390, 0.0031247, 0.0111986,
& -0.0109394,0.0058112, 0.2739046, -0.0155682, -0.0253272,
& 0.0163782, 0.0205730, 0.0022081, 0.0112749,-0.0098427,
& 0.0072705, 0.0195189, -0.0081132, -0.0071889, -0.0579970,
& -0.0856642, 0.1884260,-0.7391512, 0.1210288, -0.0241888,
& -0.0052464, -0.0096312, -0.0044834, 0.0201764, 0.0258343,
&0.0083033, 0.0077187/
DATA FEL2/0.0586055,0.0102236,-0.0396107,
& -0.0167860, -0.2019911, -0.5810815,0.0379916, 3.7508268,
& 1.8133030, -0.0564250, -0.0557352, 0.1335347, -0.0142641,
& -0.1024618,0.0970994, -0.0751830,-0.1274948, 0.0402073,
& 0.0386290, 0.1883088, 0.1838960, -0.7848989,0.7591817,
& -0.9302389,-0.8560960, 0.6633250, -4.6363869, -13.2599277,
& 0.1002136, 0.0855714,-0.0991981, -0.0765378,-0.0455264,
& 0.1169326, -0.2604067, 0.1800076, -0.2223685, -0.6347679,
&0.5334222, -0.3459502,-0.1573697, 0.8589464, 1.7815990,
&-6.3347645, -3.1513653, -9.9927750,13.3327637, -35.4897308,
&37.3466339, -0.5257398, 0.0571474, -0.5421217, 0.2404770,
& -0.1747774,-0.3433644, 0.4829708,0.3935944, 0.4885033,
& 0.8488121, -0.7640999, -1.8884945, 3.2930784,-7.3497229,
& 0.1672821,-0.2306652, 10.5782146, 12.6031065, 8.6579742,
& 215.5209961, -27.1419220,22.3405762,1108.6394043/
K=0
DO 10 I=1,72
K=K+1
FIELD(K)=FEL1(I)
10 FIELD(72+K)=FEL2(I)
RETURN
END
SUBROUTINE FIELDG(DLAT,DLONG,ALT,X,Y,Z,F,DIP,SMODIP)
C THIS IS A SPECIAL VERSION OF THE POGO68/107
C MAGNETIC FIELD LEGENDRE MODEL (IN UNITS OF GAUSS).
C F IS TOTAL FIELD, Z THE (DOWNWARD) VERTICAL COMPONENT WHILE
C X AND Y ARE COMPONENTS IN THE EQUAT. PLANE X TO ZERO LONGITUDE).
C DIP=INCLINATION ANGLE/DEGREE. MODIP=RAWER"S MODFIED DIP.
C REFERENCE@ METEOROLOGICAL AND ASTRONOMICAL INFLUENCES ON RADIO
C WAVE PROPAGATION. PERGAMON PRESS,1963,ALSO RADIO SCI. 10,669,75
C INPUT@ DLAT, DLONG=GEOGRAPHIC COORDINATES/DEG., ALT=ALTITUDE/KM.
C THE SET OF COEFFICIENTS IS GIVEN BY SUBROUTINE KOEFFB.
DIMENSION H(144),XI(3)
COMMON/BLOCK4/G(144)
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*(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+M+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=ARSIN(Z/F)
SMODIP=ARSIN(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,
&XHMF2,FOF2)
C D.BILITZA,20.8.78,CALCULATES XHMF2 AND NMF2 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)
DATA QF/11,11,8,4,1,0,0,0,0/,QM/6,7,5,2,1,0,0/
REWIND 5
LMONTH=IFIX(XMONTH)
CALL CCIRCA(R,LMONTH,FF0,MM0)
FOF2=GAMMA1(XMODIP,XLATI,XLONGI,.FALSE.,HOUR,6,QF,9,76,13,988,
&FF0)
XM3000=GAMMA1(XMODIP,XLATI,XLONGI,.FALSE.,HOUR,4,QM,7,49,9,441,
&MM0)
XHMF2=HMF2ED(XMAGBR,R,FOF2/FOE,XM3000)
RETURN
END
REAL FUNCTION GAMMA1(SMODIP,SLAT,SLONG,UT,XHOUR,IHARM,NQ,K1,M,
&MM,M3,SFE)
C SHEIKH,4.3.77,CALCULATES THE RATIO (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
C CCIR-FORTRAN-TAPE.
DIMENSION FOF2(988),SM3000(441),F2(13,76,2),FM3(9,49,2)
DO 10 I=1,IMONTH
READ(10,1) F2
10 READ(10,1) FM3
1 FORMAT(8E15.8)
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 DATA FOR THICKNESS OF 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,
&150.0,116.0,138.0,123.0,94.0,132.0,72.0,84.0,83.0,89.0,75.0,85.0,
&57.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 THIEMANN,1979,COEFFICIENTS PG1O FOR CALCULATING O+ PROFILES
C BELOW THE F2-MAXIMUM. CHOSEN TO APPROACH DANILOV-
C SEMENOV"S COMPILATION.
DIMENSION PG1O(80)
REAL FELD (80)
DATA FELD/-11.0,-11.0,4.0,-11.0,0.08018,
&0.13027,0.04216,0.25 ,-0.00686,0.00999,
&5.113,0.1 ,170.0,180.0,0.1175,0.15,-11.0,
&1.0 ,2.0,-11.0,0.069,0.161,0.254,0.18,0.0161,
&0.0216,0.03014,0.1,152.0,167.0,0.04916,
&0.17,-11.0,2.0,2.0,-11.0,0.072,0.092,0.014,0.21,
&0.01389,0.03863,0.05762,0.12,165.0,168.0,0.008,
&0.258,-11.0,1.0,3.0,-11.0,0.091,0.088,
&0.008,0.34,0.0067,0.0195,0.04,0.1,158.0,172.0,
&0.01,0.24,-11.0,2.0,3.0, -11.0,0.083,0.102,
&0.045,0.03,0.00127,0.01,0.05,0.09,167.0,185.0,
&0.015,0.18/
K=0
DO 10 I=1,80
K=K+1
10 PG1O(K)=FELD(I)
RETURN
END
SUBROUTINE KOEFP2(PG2O)
C THIEMANN,1979,COEFFICIENTS FOR CALCULATION OF O+ PROFILES
C ABOVE THE F2-MAXIMUM (DUMBS,SPENNER@AEROS-COMPILATION)
DIMENSION PG2O(32)
REAL FELD(32)
DATA FELD/1.0,-11.0,-11.0,1.0,695.0,-.000781,
&-.00264,2177.0,1.0,-11.0,-11.0,2.0,570.0,
&-.002,-.0052,1040.0,2.0,-11.0,-11.0,1.0,695.0,
&-.000786,-.00165,3367.0,2.0,-11.0,-11.0,2.0,
&575.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 THIEMANN,1979,COEFFICIENTS FOR CALCULATING O2+ PROFILES.
C CHOSEN AS TO APPROACH DANILOV-SEMENOV"S COMPILATION.
DIMENSION PG3O(80)
REAL FELD(80)
DATA FELD/-11.0,1.0,2.0,-11.0,160.0,31.0,130.0,
&-10.0,198.0,0.0,0.05922,-0.07983,
&-0.00397,0.00085,-0.00313,0.0,-11.0,2.0,2.0,-11.0,
&140.0,30.0,130.0,-10.0,
&190.0,0.0,0.05107,-0.07964,0.00097,-0.01118,-0.02614,
&-0.09537,
&-11.0,1.0,3.0,-11.0,140.0,37.0,125.0,0.0,182.0,
&0.0,0.0307,-0.04968,-0.00248,
&-0.02451,-0.00313,0.0,-11.0,2.0,3.0,-11.0,
&140.0,37.0,125.0,0.0,170.0,0.0,
&0.02806,-0.04716,0.00066,-0.02763,-0.02247,-0.01919,
&-11.0,-11.0,4.0,-11.0,140.0,45.0,136.0,-9.0,
&181.0,-26.0,0.02994,-0.04879,
&-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 SUBR. FOR CHOOSING THE REQUIRED ION DENSITY PARAMETER
C ARRAY. THE INPUT FIELD INCLUDES DIFF. SETS OF DIMENSION M EACH
C CARACTERISED BY 4 HEADER NUMBERS. RFE(4) SHOULD CONTAIN THE
C CHOSEN HEADER NUMBERS.FE(M) IS THE CORRESPONDING SET.
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 D.BILITZA,1978 PROCEDURE FOR TIME-INTERPOLATION
C USING EPSTEIN STEP FUNCTION OF ONE HOUR WIDTH (DAY/NIGHT
C FUNCTION). 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.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,JMONTH),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
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