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