C INDEX SUBROUTINES:___INDX2018.FOR___T.L.Gulyaeva___May. 2018 C C INDEX SUBROUTINES:___INDX2017.FOR___T.L.Gulyaeva___May. 2017 C C INDEX SUBROUTINES:___INDX2015c.FOR___T.L.Gulyaeva___Feb. 2016 C C INDEX SUBROUTINES:___INDX2015.FOR___T.L.Gulyaeva___Dec. 2014 C C == AMENDMENTS ...................................Dec. 2014 C STORM RELATED SUBROUTINE FROM IRI2001: = = APF == is replaced C by: == SUBAPF == which enters geomagnetic Kp and Ap 3-h indices C and solar daily Rz and F107 indices from standard = kpYEAR = file C for given 'YEAR' and preceding 'YEAR-1' for calculations of 81-days C indices R81 and F81 instead of former APRZ.DAT file (cancelled input) C C INDEX SUBROUTINES:___INDX2014.FOR___T.L.Gulyaeva___Apr. 2014 C C INDEX SUBROUTINES:___INDX2013.FOR___T.L.Gulyaeva___Dec. 2013 C C INDEX SUBROUTINES:___INDX2011.FOR___T.L.Gulyaeva___Dec. 2011 C C INDEX SUBROUTINES:___INDX2007.FOR___T.L.Gulyaeva___June 2009 C C ******************* ISO_IRI_SMI PROJECT ******************* C C ________________________ISOMAIN3.FOR_______________________ C C UPDATED FOR ISO_IRI CORRECTIONS____ C C Changes from ISOMAIN2.FOR_:________________________Sep. 2005 C C New SUBROUTINE TOPH05 is composed for modeling h05top height C at Ne05=0.5*NmF2 deduced from ISIS and IK19 topside electron C density profiles. Topside IRI-Bent model is corrected with C QF-factor obtained from h05top, NmF2 and hmF2. C Dependence of h05top(SSN) introduced C C == AMENDMENTS : T.L.Gulyaeva Sep. 2001 C C == STORM RELATED SUBROUTINES FROM IRI2001: APF, STORM == C Ref.[IRIFUN.FOR, version 2001.1, May 7, 2001 D. Bilitza] C C 2001.1 05/07/01 Include storm subroutine STORM and Ap access s/w C The subroutine APF uses the input unit 13 for reading 3-hourly C magnetic Ap indices. The subroutine TCON uses input unit 12 for C reading the 12-month-running means of ionospheric index IG and C of sunspot number Rz. C C == GKPM-subroutine using APF-read ap_indices produces: C == instantaneous kp-index from ap-index C == time-weighted accumulation Kpm index from kp-indices C == for plasmasphere model C == F81=COV81-subroutine produces daily F10.7 (Covington) solar C == radio flux index averaged for 81 days preceding given day C == R81=Rz81-subroutine produces daily sunspot number Rz C == averaged for 81 days preceding given day C************************************************************** C************************************************************** C C********** INTERNATIONAL REFERENCE IONOSPHERE **************** C**********RUSSIAN STANDARD MODEL OF IONOSPHERE**************** C**************** FUNCTIONS,SUBROUTINES ********************* C************************************************************** C** initialize: INITIALIZE (needs to be called before using C** subroutines or functions) C** NE: XE1,ZERO,DXE1N,XE2,XE3_1,XE4_1,XE5,XE6,XE_1 C** TE/TI: ELTEIK,SPHARM_IK,TEBA,SPHARM,ELTE,TEDE,TI,TEDER, C** TN,DTNDH C** NI: RPID,RDHHE,RDNO,KOEFP1,KOEFP2,KOEFP3,SUFE C** IONCO2, APROK,IONCOM_NEW,IONCO1 C** PEAKS: FOUT,XMOUT,HMF2ED,FOF1ED,f1_c1,f1_prob,FOEEDI,XMDED, C** GAMMA1 C** PROFILE PAR:B0_98,TAL,VALGUL C** MAG. FIELD: GGM,FIELDG,CONVER(corrected latitude) C** FUNCTIONS: REGFA1 C** TIME: SOCO,HPOL,MODA,UT_LT C** EPSTEIN: RLAY,D1LAY,D2LAY,EPTR,EPST,EPSTEP,EPLA C** LAY: XE2TO5,XEN,VALGUL,ROGUL,LNGLSN,LSKNM,INILAY C** INDICES: TCON,TCONGEC, SUBAPF,GKPM C** Storm: LSTID,STORM, SUBAPF,GKPM C** ion drift: vdrift C************************************************************** C C************************************************************** C*** -------------------ADDRESSES------------------------ *** C*** I PROF. K. RAWER DR. D. BILITZA I *** C*** I HERRENSTR. 43 GSFC CODE 933 I *** C*** I 7801 MARCH 1 GREENBELT MD 20771 I *** C*** I F.R.G. USA I *** C*** ---------------------------------------------------- *** C*** I DR. T.L. GULYAEVA I *** C*** I IZMIRAN 142190 I *** C*** I TROITSK, MOSCOW I *** C*** I RUSSIA I *** C*** I gulyaeva@izmiran.ru I *** C*** ---------------------------------------------------- *** C************************************************************** C C SUBROUTINE GKPM(IAS,DKP,RKP) C C T.L. Gulyaeva, Dec. 2013: C New: put real kp instead of Kpm calculation: C C T.L. Gulyaeva, Sep. 2001: C =================================================================== C calculates geomagnetic 3h kp-index, dkp(13), and Kpm index, rkp, C used by the plasmasphere model C for given year, month, day and UT hour using ap-indices C collected by APF-subroutine in the array ias(13). C The ap-indices are transformed to kp-indices with conventional C thresholds for ap - kp . C Kpm presents forecast of kp-index 3h in advance based on accumulation C of instantaneous kp-indices for 12 preceding hours ranked by C decreasing order. C C Ref. T.L. Gulyaeva, G. De Franceschi, and L. Perrone. 'Electron C temperature variations at the F2 layer peak height durung the space C weather month of September 1999'. Adv. Space Res.31(4), 965-970,2003. C DIMENSION ias(13),lap(0:27),tkp(0:27),dkp(13),X1(4),X2(4) &,icurkp(0:7),iprekp(0:7) DATA lap/0,2,3,4,5,6,7,9,12,15,18,22,27,32,39,48,56,67,80,94,111 &,132,154,179,207,236,300,400/ DATA tkp/0.0,.3,.7,1.0,1.3,1.7,2.0,2.3,2.7,3.0,3.3,3.7,4.0,4.3 &,4.7,5.0,5.3,5.7,6.0,6.3,6.7,7.0,7.3,7.7,8.0,8.3,8.7,9.0/ COMMON /FIKP/icurkp,iprekp,iyear_cur,imn_cur,idy_cur C C calculate 3h kp-indices : kp(13) presents kp-index for current 3h UT bin C DO J=1,13 IF (IAS(J).EQ.0) THEN DKP(J)=0.0 GOTO 2 ELSE DO I=1,27 IF ((IAS(j).GT.LAP(i-1)).AND.(IAS(j).LE.LAP(i))) THEN DKP(j)=TKP(i) GOTO 2 ENDIF ENDDO ENDIF 2 ENDDO C New: put real kp instead of Kpm calculation: RKP=dkp(13) goto 870 C calculate Kpm-index RKP=0.0 STAU=1.0 C K1=13 K1=8 DO N=1,4 NK=K1+N X1(N)=DKP(NK) X2(N)=0.0 ENDDO MM=0 C rank 4 preceding 3h kp-indices by decreasing order DO M=1,4 C find naximum of 4 preceding 3h kp-indices (XM) XM=0.0 DO N=1,4 IF (XM.LT.X1(N)) THEN XM=X1(N) ENDIF ENDDO IF (XM.EQ.0.0) GOTO 860 C ranking 4 preceding kp by decreasing order DO N=1,4 IF (X1(N).LT.XM) GOTO 830 MM=MM+1 X2(MM)=XM X1(N)=0.0 830 CONTINUE ENDDO ENDDO C 860 DO K=1,4 RKP=RKP+STAU*X2(K) STAU=STAU*.7788 ENDDO RKP=RKP*0.35 870 RETURN END C************************************************************** C C subroutine initialize dimension indap(13),akp(13) COMMON /CONST/UMR,PI/ARGEXP/ARGMAX & /const1/humr,dumr,xfap,indap,akp,cov,rzs ARGMAX=88.0 pi=atan(1.0)*4. UMR=pi/180. humr=pi/12. dumr = pi / 182.5 return end C C************************************************************* CC SUBROUTINE SUBAPF(IYYYY,IMN,IDY,HOUR,IAP,F81,R81) CC------------------------------------------------------------------- CC T.L. Gulyaeva........................................May. 2017 C CC Source: former IRI-Plas subroutine: CC SUBROUTINE APF(IYYYY,IMN,ID,HOUR,IAP,F81,R81) CC Modified to produce sunspot number R81, and solar radio flux F81 CC from two successive kpYEAR annual files CC Input: CC kpYEAR (preceding year: IYYYY-1) CC kpYEAR (given year: IYYYY) c-------------------------------------------------------------------- c D. Bilitza .......................................... May, 2001. c c finds Ap indices for IRI-storm for given year IYYYY (yyyy), month c (IMN) and day (ID) and UT hour (HOUR decimal hours). The indices are c stored in IAP(13). IAP(13) is ap index for UT=hour. Intervals are UT: c (0-3),(3-6),(6-9),(9-12),(12-15),(15-18),(18-21),(21-24). c c T.L. Gulyaeva ...................................... Sep. 2001. c Modified for producing solar radio flux F10.7, F81, averaged for c 81 days preceeding given day C C------------------------------------------------------------------------ DIMENSION iap(13),iap0(8),iap1(8),iap2(8),iapt(24),lm(12) &,icv(82),irz(82),icurkp(0:7),iprekp(0:7) CHARACTER*4 YEAR,YEAR_pre CHARACTER*80 infilekp COMMON /FIKP/icurkp,iprekp,iyear_cur,imn_cur,idy_cur DATA LM/31,28,31,30,31,30,31,31,30,31,30,31/ C C START: ipcnt=0 ! keep cnt -82,...,-1 do i=1,8 iap0(i)=0 iap1(i)=0 iap2(i)=0 enddo do k=1,82 irz(k)=0 icv(k)=0 enddo C Input file kpYEAR starts from 1948: if((iyyyy.lt.1948).or.(iyyyy.gt.iyear_cur)) goto 21 !* IF ((iyyyy.eq.iyear_cur).and.(imn.gt.imn_cur)) goto 21 !* IF ((iyyyy.eq.iyear_cur).and.(imn.eq.imn_cur).and. !* * (idy.gt.idy_cur)) goto 21 !*TEST C iy=iyyyy-1900 770 z1=iy/4.0 lm(2)=28 if(iyyyy/4*4.eq.iyyyy) lm(2)=29 if (iy.ge.100) then iyr=iy-100 else iyr=iy endif call blet4(iyyyy,YEAR) iYYYY_pre=iYYYY-1 call blet4(iYYYY_pre,YEAR_pre) C idend=idy C C Day-of-year: ldaend=ndoy(iyyyy,imn,idend) ! current day-of-year C ipre=0 if (ldaend.lt.81) then ipre=1 ! Include data for preceding year endif C infilekp='kpYEAR' C 10 FORMAT(3I2,6X,8(I2),3X,8(I3),7X,I3,I3,1X,I1) C C------------------------------------------------------------------------------- C Go to read current year file kpYEAR (avoid preceding year input) if current day IDY>81 day-of-year: C if (ipre.eq.0) GOTO 203 !avoid input of preceding year file kpYEAR>>>>> C C------------------------------------------------------------------------------- C Include Input file 'kpYEAR' for preceding year C infilekp(3:6)=YEAR_pre !!! Tamara OPEN(13,FILE=infilekp,STATUS='OLD',ERR=21) 201 continue do jj=1,8 iap2(jj)=iap1(jj) iap1(jj)=iap0(jj) enddo READ(13,10,END=202) JYR,JMN,JDY,(icurkp(k),k=0,7) &,(iap0(i),i=1,8),ir1,if1,if2 C C Replace missed F10.7 by regression model: if (if1.eq.0) then cov=63.8255+0.8868*ir1 ! daily F10.7(Rz) model icv(82)=nint(cov*10.) else icv(82)=if1*10+if2 !F107*10 endif irz(82)=ir1 C Move data 1 line up: do k=1,81 irz(k)=irz(k+1) icv(k)=icv(k+1) enddo GOTO 201 202 close(unit=13) ipcnt=81-ldaend ! C C For current year: 203 infilekp(3:6)=YEAR ! C OPEN(13,FILE=infilekp,ERR=21) C C Read file kpyear for current year: C 204 continue do jj=1,8 iprekp(jj-1)=icurkp(jj-1) iap2(jj)=iap1(jj) iap1(jj)=iap0(jj) enddo READ(13,10,err=29,END=205) JYR,JMN,JDY,(icurkp(k),k=0,7) &,(iap0(i),i=1,8),ir1,if1,if2 irz(82)=ir1 if (if1.eq.0) then cov=63.8255+0.8868*ir1 ! daily F10.7(Rz) model icv(82)=nint(cov*10.) else icv(82)=if1*10+if2 !F107*10 endif C Move data 1 line up: 140 do k=1,81 irz(k)=irz(k+1) icv(k)=icv(k+1) enddo C ipcnt=ipcnt+1 if (ipcnt.lt.81) goto 204 IF ((JYR.eq.IYR).and.(JMN.eq.IMN).and.(JDY.eq.idend)) then goto 207 ELSE ipcnt=ipcnt-1 goto 204 ENDIF C data for 81 preceding days are collected 207 F81=0.0 R81=0.0 icnv=0 do k=1,81 R81=R81+irz(k) if (icv(k).gt.0) then F81=F81+icv(k) icnv=icnv+1 endif enddo F81=F81/(icnv*10.) R81=R81/81.0 C C === collection of ap-indices for 13 3h UT bins: 205 do jj=1,8 iapt(jj)=iap2(jj) iapt(8+jj)=iap1(jj) iapt(16+jj)=iap0(jj) enddo if(hour.ge.24.0) hour=0.0 ih=int(hour/3.)+1 ihb=16+ih ilauf=ihb do jj=13,1,-1 iap(jj)=iapt(ilauf) ilauf=ilauf-1 enddo C goto 30 21 write(*,100) 100 format(1X,'Date is outside range of Kp-Ap indices file.', & ' STORM model is turned off.') IAP(1)=-iap(1) goto 30 29 write(*,*) 'Error in input k-index file',infilekp C pause ' ' stop 30 CLOSE(unit=13) RETURN END C __________________________________________________________________________ C C integer function ndoy(nyear,nmn,ndy) C To define day-of-year C DIMENSION IM(0:12) DATA IM/31,31,28,31,30,31,30,31,31,30,31,30,31/ IF(nyear/4*4.EQ.nyear) THEN IM(2)=29 ELSE IM(2)=28 ENDIF C Day-of-year LDA1 mosum=0 C NEW++++++++++++++++++++ if (nmn.eq.0) then nmn=12 nyear=nyear-1 endif if(nmn.gt.1) then do 1234 i=1,nmn-1 1234 mosum=mosum+im(i) endif ndoy=mosum+ndy END FUNCTION C ----------------------------------------------------------- subroutine blet2(jlet,alet2) c CHARACTER*1 IN(0:9) C integer*2 ilet DATA IN/'0','1','2','3','4','5','6','7','8','9'/ CHARACTER*2 alet2 ilet=jlet alet2='00' j10=ilet/10 j1=ilet-j10*10 5 do i=0,9 if (j10.eq.i) then alet2(1:1)=IN(i) endif enddo 7 do i=0,9 if (j1.eq.i) then alet2(2:2)=IN(i) endif enddo return end C-------------------------------------------------------------------------------- C subroutine blet3(ilet,alet3) c nn=1,2,3,4 CHARACTER*1 IN(0:9) DATA IN/'0','1','2','3','4','5','6','7','8','9'/ CHARACTER*3 alet3 C j100=ilet/100 j10=(ilet-j100*100)/10 j1=ilet-j100*100-j10*10 3 do i=0,9 if (j100.eq.i) then alet3(1:1)=IN(i) endif enddo 5 do i=0,9 if (j10.eq.i) then alet3(2:2)=IN(i) endif enddo 7 do i=0,9 if (j1.eq.i) then alet3(3:3)=IN(i) endif enddo return end C C--------------------------------------------------------------- subroutine blet4(ilet,alet4) c nn=1,2,3,4 INTEGER ilet CHARACTER*4 alet4 CHARACTER*1 IN(0:9) DATA IN/'0','1','2','3','4','5','6','7','8','9'/ alet4='0000' C j1000=ilet/1000 j100=(ilet-j1000*1000)/100 j10=(ilet-j1000*1000-j100*100)/10 j1=ilet-j1000*1000-j100*100-j10*10 3 do i=0,9 if (j1000.eq.i) then alet4(1:1)=IN(i) endif enddo 4 do i=0,9 if (j100.eq.i) then alet4(2:2)=IN(i) endif enddo 5 do i=0,9 if (j10.eq.i) then alet4(3:3)=IN(i) endif enddo 7 do i=0,9 if (j1.eq.i) then alet4(4:4)=IN(i) endif enddo return end C =============================================================== c SUBROUTINE STORM(ap,rga,rgo,coor,rgma,ut,doy,cf,rap) C E.A. Araujo-Pradere and T.J. Fuller-Rowell .............May, 2001 C C Refs. C (i) T.J.Fuller-Rowell, M.V.Codrescu, and E.A.Araujo-Pradere, C Capturing the storm-time ionospheric response in an empirical C model. AGU Geophys. Monograph, 125, 393-401, 2001. C (ii)E.A. Araujo-Pradere, T.J. Fuller-Rowell, and M.V. Codrescu, C STORM: An empirical storm-time ionospheric correction model, C Radio Sci.,37(5),1070, DOI:10.1029/2001RS002467, 2002. C C Fortran code to obtain the foF2 storm-time correction factor at a C given location and time, using the current and the 12 previous C ap values as input. C C ap ---> (13 elements integer array). Array with the preceeding C 13 value of the 3-hourly ap index. The 13th value C in the array will contain the ap at the UT of interest, C the 12th value will contain the 1st three hourly interval C preceeding the time of interest, and so on to the 1st C ap value at the earliest time. C coor --> (integer). If coor = 2, rga should contain the geomagnetic C latitude. C If coor = 1, rga should contain the geographic C latitude. C rga ---> (real, -90 to 90) geographic or geomagnetic latitude. C rgo ---> (real, 0 to 360, positive east from Greenwich.) C geographic longitude, only used if coor=1. C rgma --> (corrected geomagnetic latitude) C ut ---> (integer, hours 00 to 23) Universal Time of interest. C doy ---> (integer, 1 to 366)Day of the year. C cf ---> (real) The output; the storm-time correction factor used C to scale foF2, foF2 * cf. C rap ---> integrated ap-index used for producing storm factor cf C DIMENSIONS AND COEFFICIENTS VALUES cc- COMMON cc- & /CONST1/HUMR,DUMR,XFAP,INDAP,AKP,COV,RZS DIMENSION c4(20) DATA c4/0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00, +0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00, +0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00/ DIMENSION c3(20) DATA c3/0.00E+00,0.00E+00,0.00E+00,0.00E+00,0.00E+00,-9.44E-12, +0.00E+00,3.04E-12,0.00E+00,9.32E-12,-1.07E-11,0.00E+00,0.00E+00, +0.00E+00,1.09E-11,0.00E+00,0.00E+00,0.00E+00,0.00E+00,-1.01E-11/ DIMENSION c2(20) DATA c2/1.16E-08,0.00E+00,0.00E+00,-1.46E-08,0.00E+00,9.86E-08, +2.25E-08,-1.67E-08,-1.62E-08,-9.42E-08,1.17E-07,4.32E-08,3.97E-08, +3.13E-08,-8.04E-08,3.91E-08,2.58E-08,3.45E-08,4.76E-08,1.13E-07/ DIMENSION c1(20) DATA c1/-9.17E-05,-1.37E-05,0.00E+00,7.14E-05,0.00E+00,-3.21E-04, +-1.66E-04,-4.10E-05,1.36E-04,2.29E-04,-3.89E-04,-3.08E-04, +-2.81E-04,-1.90E-04,4.76E-05,-2.80E-04,-2.07E-04,-2.91E-04, +-3.30E-04,-4.04E-04/ DIMENSION c0(20) DATA c0/1.0136E+00,1.0478E+00,1.00E+00,1.0258E+00,1.00E+00, +1.077E+00,1.0543E+00,1.0103E+00,9.9927E-01,9.6876E-01,1.0971E+00, +1.0971E+00,1.0777E+00,1.1134E+00,1.0237E+00,1.0703E+00,1.0248E+00, +1.0945E+00,1.1622E+00,1.1393E+00/ DIMENSION fap(36) DATA fap/0.,0.,0.037037037,0.074074074,0.111111111,0.148148148, 10.185185185,0.222222222,0.259259259,0.296296296,0.333333333, 20.37037037,0.407407407,0.444444444,0.481481481,0.518518519, 30.555555556,0.592592593,0.62962963,0.666666667,0.703703704, 40.740740741,0.777777778,0.814814815,0.851851852,0.888888889, 50.925925926,0.962962963,1.,0.66666667,0.33333334,0.,0.333333, 60.666666,1.,0.7/ integer code(8,6) data code/3,4,5,4,3,2,1,2,3,2,1,2,3,4,5,4,8,7,6,7,8,9,10,9, *13,12,11,12,13,14,15,14,18,17,16,17,18,19,20,19,18,17,16,17, *18,19,20,19/ INTEGER ape(39) INTEGER ap(13) INTEGER ut,doy,dayno,coor,s1,s2,l1,l2 REAL rgma, rap, rga, rgo, rs, rl C C ----------------------------------------------------------------------- rgo1=rgo C CALLING THE PROGRAM TO CONVERT TO GEOMAGNETIC COORDINATES IF (coor .EQ. 1) THEN CALL CONVER (rga,rgo,rgma) ELSE IF (coor .EQ. 2) THEN rgma = rga ELSE WRITE (6,*)' ' WRITE (6,*)' ' WRITE (6,*)' Wrong Coordinates Selection -------- >>', coor WRITE (6,*)' ' GOTO 100 ENDIF IF (COOR.EQ.2) THEN RGMA=RGA C FROM 3-HOURLY TO HOURLY ap i = 1 DO 10 k = 1,13 DO j = 1,3 ape(i) = ap(k) i = i + 1 END DO 10 CONTINUE C TO OBTAIN THE INTEGRAL OF ap. C INTEGRAL OF ap if(ut.eq.24) ut=0 IF (ut .EQ. 0 .OR. ut .EQ. 3 .OR. ut .EQ. 6 .OR. ut .EQ. 9 .OR. 1ut .EQ. 12 .OR. ut .EQ. 15 .OR. ut .EQ. 18 .OR. ut .EQ. 21) THEN k = 1 ELSE IF (ut .EQ. 1 .OR. ut .EQ. 4 .OR. ut .EQ. 7 .OR. ut .EQ. 10 1.OR.ut .EQ. 13 .OR. ut .EQ. 16 .OR. ut .EQ. 19 .OR. ut .EQ. 22) 2THEN k = 2 ELSE IF (ut .EQ. 2 .OR. ut .EQ. 5 .OR. ut .EQ. 8 .OR. ut .EQ. 11 1.OR. ut .EQ. 14 .OR. ut .EQ. 17 .OR. ut .EQ. 20 .OR. ut .EQ. 23) 2THEN k = 3 ELSE WRITE (6,*)' ' WRITE (6,*)' ' WRITE (6,*)' Wrong Universal Time value -------- >>', ut WRITE (6,*)' ' GOTO 100 END IF rap = 0 DO j = 1,36 rap = rap + fap(j) * ape(k+j) END DO if(rap.le.200.)then cf=1.0 goto 100 end if C Test TLG: 111 continue if(doy.gt.366.or.doy.lt.1)then WRITE (6,*)' ' WRITE (6,*)' ' WRITE (6,*)' ' WRITE (6,*)' Wrong Day of Year value --- >>', doy WRITE (6,*)' ' GOTO 100 end if if(rgma.gt.90.0.or.rgma.lt.-90.0)then WRITE (6,*)' ' WRITE (6,*)' ' WRITE (6,*)' ' WRITE (6,*)' Wrong GEOMAGNETIC LATITUDE value --- >>', rgma WRITE (6,*)' ' GOTO 100 end if c write(6,*)rgma dayno=doy if(rgma.lt.0.0)then dayno=doy+172 if(dayno.gt.365)dayno=dayno-365 end if if (dayno.ge.82) rs=(dayno-82.)/45.6+1. if (dayno.lt.82) rs=(dayno+283.)/45.6+1. s1=rs facs=rs-s1 s2=s1+1 if(s2.eq.9) s2=1 c write(6,*)s1,s2,rs rgma = abs(rgma) rl=(rgma+10.)/20.+1 if(rl.eq.6.0)rl=5.9 l1=rl facl=rl-l1 l2=l1+1 c write(6,*)l1,l2,rl C FACTORS CALCULATIONS if(rap.lt.300.)then rapf=300. n1=code(s1,l1) cf1=c4(n1)*(rapf**4)+c3(n1) * (rapf**3) + c2(n1) * (rapf**2) + 1c1(n1) * rapf + c0(n1) n2=code(s1,l2) cf2=c4(n2)*(rapf**4)+c3(n2) * (rapf**3) + c2(n2) * (rapf**2) + 1c1(n2) * rapf + c0(n2) n3=code(s2,l1) cf3=c4(n3)*(rapf**4)+c3(n3) * (rapf**3) + c2(n3) * (rapf**2) + 1c1(n3) * rapf + c0(n3) n4=code(s2,l2) cf4=c4(n4)*(rapf**4)+c3(n4) * (rapf**3) + c2(n4) * (rapf**2) + 1c1(n4) * rapf + c0(n4) C INTERPOLATION cf300=cf1*(1 - facs) * (1 - facl) + cf2 * (1 - facs) * (facl) + *cf3 * (facs) * (1 - facl) + cf4 * (facs) * (facl) cf = (cf300-1.0)*rap/100.-2.*cf300+3. goto 100 end if n1=code(s1,l1) c write(6,*)n1 cf1 = c4(n1) * (rap**4) + c3(n1) * (rap**3) + c2(n1) * (rap**2) + 1c1(n1) * rap + c0(n1) n2=code(s1,l2) cf2 = c4(n2) * (rap**4) + c3(n2) * (rap**3) + c2(n2) * (rap**2) + 1c1(n2) * rap + c0(n2) n3=code(s2,l1) cf3 = c4(n3) * (rap**4) + c3(n3) * (rap**3) + c2(n3) * (rap**2) + 1c1(n3) * rap + c0(n3) n4=code(s2,l2) cf4 = c4(n4) * (rap**4) + c3(n4) * (rap**3) + c2(n4) * (rap**2) + 1c1(n4) * rap + c0(n4) c INTERPOLATION cf = cf1 * (1 - facs) * (1 - facl) + cf2 * (1 - facs) * (facl) + *cf3 * (facs) * (1 - facl) + cf4 * (facs) * (facl) 100 CONTINUE RETURN END C T.L. Gulyaeva ............................................ ........ 2007 C C New subroutine TOP05 to correct Topside IRI-Bent model using C QF-factor obtained from h05top height at Ne=0.5*NmF2, NmF2 and hmF2 . C CREM SUBROUTINE TOPH05(covi,AMLAT,TIME,HMAX,HT05,sg) SUBROUTINE TOPH05(rz,AMLAT,TIME,HMAX,HT05,sg) C C C// Gulyaeva T.L. (2003) Variations in the half-width of the topside ionosphere C// according to the observations by space ionosondes ISIS 1,ISIS 2, and IK19. C// International J. of Geomagnetism and Aeronomy, 4(3), 201-207. C// Gulyaeva T.L., Titheridge J.E. (2006) Advanced specification of electron density C// and temperature in the IRI ionosphere-plasmasphere model. C// Adv. Space Res. 38(11), 2587-2595, doi:10.1016/j.asr.2005.08.045. C C Implementation of empirical RAT=(h05top-hmF2)/hmF2 derived from ISIS and IK19 C topside electron density profiles to obtain half peak density topside height C h05top from the Chebishev polinomial coefficients given for C (1) 4 levels of solar activity: Rz= 0, 50, 100, 150 C solar radio flux covi=60, 106, 152, 198 (F10.7 option) C (2) 10 selected grids of geomagnetic latitude (N=S):0,10,20,30,40,50,60,70,80,90 C (3) 5 selected grids of local time: 0, 6, 12, 18, 24. C (4) 4 seasonal grids: 1 equinox(sg=90deg), 2 summer (sg=180), C 3 equinox (sg=270), 4 winter(sg=360) C DIMENSION CVLEV(4) COMMON /BLOCK1/HMF2,XNMF2,XHMF1,F1REG & /QTOP/Y05,H05TOP,QF,XNETOP,XM3000,HHALF,TAU,HTOP,HCM,DHTOP,TCM crem DATA CVLEV/60.,106.,152.,198./ DATA CVLEV/0.,50.,100.,150./ ! Rz levels LOGICAL F1REG ABMLAT=ABS(AMLAT) crem IR=IFIX((covi-60.)/46.)+1 IR=IFIX(rz/50.)+1 M1=IFIX(ABMLAT/10.)+1 L1=IFIX(TIME/6.)+1 M2=M1+1 IF(M1.EQ.10) M2=10 L2=L1+1 IF(L1.EQ.5) L2=5 C C INTERPOLATE RAT FOR GIVEN RZI C C Call Chebishev approximation to interpolate for given ABMLAT, HRLT call CHEBISH(CVLEV(IR),TIME,ABMLAT,XX,SG) IF (IR.EQ.4) THEN RAT05=XX GOTO 10 ENDIF call CHEBISH(CVLEV(IR+1),TIME,ABMLAT,YY,sg) crem RAT05=XX+(YY-XX)*(COVI-CVLEV(IR))/46. RAT05=XX+(YY-XX)*(RZ-CVLEV(IR))/50. 10 HT05=HMAX*(1.+RAT05) C Add include update by DHTOP HT05=HT05+DHTOP RETURN END C******************************************************************** SUBROUTINE CHEBISH(COVS,HOURLT,ABMLAT,RATCH,SG) C CHEBISHEV POLINOMIALS FOR ABMLAT(10),HOURLT(5) C CR((C0...C5),(LT=0,6,...24),(SG=season grids=90,180,270,360) C (COV=60,106,152,198) C Ref. T.L.Gulyaeva, J.E.Titheridge. Advanced specification of electron density and C temperature in the IRI ionosphere-plasmasphere model. Adv. Space Res., 2005. c REAL UK(0:10),CR(0:5,5,3,4),YI(5),YY(5,3) REAL BR(6,5,3,4),YI(5),YY(5,3) REAL PL1(5),PL2(5),PL3(5),CL(0:3) C DATA rad/0.01745329/ DATA PL1/-2.,-1.,0.,1.,2./ DATA PL2/2.,-1.,-2.,-1.,2./ DATA PL3/-1.,2.,0.,-2.,1./ DATA BR/ C Polinomial Coefficients B1,B2,B3,B4,B5,B6 for COV=60 (mlat/10=0,1,...,9) C Equinox B0MLAT: * -1.5859, 3.5789, -3.7884, 2.7094,-1.2962,.6759 *, -5.3449, 12.8379,-12.0165, 5.9746,-1.9084,.7669 *,-12.8000, 35.3084,-38.0043,19.6004,-4.4974,.6975 *, 5.8282,-13.3538, 9.1674,-0.9593,-0.8909,.6062 *, -1.5859, 3.5789, -3.7884, 2.7094,-1.2962,.6759 C Summer B0MLAT *, -7.1103, 21.0389,-24.5539,14.1607,-3.8537,.7266 *, 5.5333,-10.6242, 4.8751, 1.0587,-1.0821,.7527 *,-15.4487, 42.9269,-45.0314,21.4718,-4.2116,.6026 *, -6.6436, 16.4533,-15.5142, 6.8287,-1.2871,.4976 *, -7.1103, 21.0389,-24.5539,14.1607,-3.8537,.7266 C Winter B0MLAT *, 14.9103,-35.2337, 27.3078,-6.5362,-0.6265,.7509 *, 2.3846, -5.8840, 3.7023, 0.8525,-1.2663,.7086 *, -9.8846, 26.6649,-27.0173,12.6959,-2.6536,.6295 *, 1.7692, -2.3578, -0.7945, 2.2477,-0.9691,.5719 *, 14.9103,-35.2337, 27.3078,-6.5362,-0.6265,.7509 C Polinomial Coefficients B1,B2,B3,B4,B5,B6 for COV=106 (mlat=0,10,...,90) C Equinox B1MLAT *, -4.1218, 10.6136,-11.4922, 6.0470,-1.3620,.5563 *, 0.9077, 2.9562, -8.9880, 6.8680,-1.9621,.7737 *,-16.2744, 42.8047,-43.7009,20.7965,-4.0697,.6619 *,-17.3038, 44.3336,-40.9249,15.9042,-2.1554,.4796 *, -4.1218, 10.6136,-11.4922, 6.0470,-1.3620,.5563 C Summer B1MLAT *, -4.9692, 16.5753,-21.3543,12.7061,-3.1758,.6446 *, 1.9000, -2.8167, -0.9962, 3.0687,-1.3454,.6859 *, 7.6769,-14.8343, 6.7030, 1.5578,-1.0626,.4291 *, 5.4833,-10.6322, 4.7571, 1.2178,-0.8223,.4615 *, -4.9692, 16.5753,-21.3543,12.7061,-3.1758,.6446 C Winter B1MLAT *, -4.7282, 13.4491,-15.6931, 8.8388,-1.9732,.5874 *, 5.6756,-14.8458, 11.8927,-2.2632,-0.6122,.6948 *,-14.2872, 40.0829,-41.2716,18.1696,-2.7203,.4916 *,-13.6128, 33.4657,-29.7231,11.0972,-1.2884,.5034 *, -4.7282, 13.4491,-15.6931, 8.8388,-1.9732,.5874 C Polinomial Coefficients B1,B2,B3,B4,B5,B6 for COV=152 (mlat=0,10,...,90) C Equinox B2MLAT *, -3.3282, 10.4296,-12.4722, 6.7623,-1.5172,.4931 *, -8.9744, 20.1311,-17.4552, 7.6518,-1.7371,.6702 *, 12.0462,-27.8932, 20.6241,-4.5781, 0.0814,.3501 *,-17.0551, 42.3258,-37.1874,13.3608,-1.4804,.4216 *, -3.3282, 10.4296,-12.4722, 6.7623,-1.5172,.4931 C Summer B2MLAT *, 7.3077,-17.1579, 11.6872,-0.7405,-1.0298,.5754 *, 19.2641,-45.1886, 34.3297,-8.1879,-0.1875,.6562 *, 6.0987,-11.0903, 4.3569, 1.4001,-0.7309,.3885 *, 5.9295,-13.9205, 10.2347,-2.2818, 0.0853,.3915 *, 7.3077,-17.1579, 11.6872,-0.7405,-1.0298,.5754 C Winter B2MLAT *, -1.6821, 8.6010,-13.6570, 8.6307,-1.9846,.5635 *, 5.4679,-12.3750, 7.5620, 0.5394,-1.4415,.6659 *, -8.0821, 21.9288,-21.8597, 9.3455,-1.4644,.3599 *, -8.3000, 19.3076,-16.3295, 6.1619,-0.9144,.3846 *, -1.6821, 8.6010,-13.6570, 8.6307,-1.9846,.5635 C Polinomial Coefficients B1,B2,B3,B4,B5,B6 for COV=198 (mlat=0,10,...,90) C Equinox B3MLAT *,-16.4051, 28.2880,-16.0982, 4.6328,-1.0405,.5486 *, 13.0846,-34.8291, 30.0074,-8.6402, 0.1529,.6165 *, 19.7474,-42.7116, 28.9430,-6.0487, 0.1492,.3748 *, 16.2795,-36.6982, 26.5094,-6.3492, 0.2926,.3946 *,-16.4051, 28.2880,-16.0982, 4.6328,-1.0405,.5486 C Summer B3MLAT *, 4.6410,-13.7931, 11.6548,-1.9248,-0.7246,.5264 *, -2.4090, 3.1805, -2.8423, 2.8861,-0.9937,.5234 *, 6.3410,-13.9643, 8.2461,-0.0186,-0.7009,.3582 *, 9.0987,-20.8618, 14.7262,-2.8798,-0.0512,.3662 *, 4.6410,-13.7931, 11.6548,-1.9248,-0.7246,.5264 C Winter B3MLAT *, -4.6526, 12.1878,-14.4047, 8.5226,-2.0493,.5903 *, 3.9821, -6.9477, 0.8382, 3.4898,-1.5694,.6283 *, -7.0474, 17.3974,-17.3465, 8.3671,-1.5708,.3759 *, 4.2782, -9.9880, 5.9834, 0.0975,-0.4900,.3842 *, -4.6526, 12.1878,-14.4047, 8.5226,-2.0493,.5903/ C DATA UL/-2.,-1.,0.,1.,2./ do k=0,3 cl(k)=0. enddo C c,, IR=IFIX((covs-60.)/46.)+1 IR=IFIX(covs/50.)+1 C Given geomagnetic latitude parameter: xi=abmlat/100. DO LS=1,3 DO LL=1,5 B1=BR(6,LL,LS,IR) B2=BR(5,LL,LS,IR) B3=BR(4,LL,LS,IR) B4=BR(3,LL,LS,IR) B5=BR(2,LL,LS,IR) B6=BR(1,LL,LS,IR) HLT=(LL-1)*6.0 YY(LL,LS)=B1+xi*(B2+xi*(B3+xi*(B4+xi*(B5+xi*B6)))) ENDDO ENDDO ! end of season/day cycle C Apply seasonal interpolation do i=1,5 p0=(2.*YY(i,1)+YY(i,2)+YY(i,3))/4. p1=(YY(i,3)-YY(i,2))/2. p2=(YY(i,2)+YY(i,3)-2.*YY(i,1))/4. YI(i)=p0+p1*cos(sg*rad)+p2*cos(2.*sg*rad) enddo DO K=1,5 CL(0)=CL(0)+YI(K) CL(1)=CL(1)+YI(K)*PL1(K) CL(2)=CL(2)+YI(K)*PL2(K) CL(3)=CL(3)+YI(K)*PL3(K) ENDDO CL(0)=CL(0)/5. CL(1)=CL(1)/10. CL(2)=CL(2)/14. CL(3)=CL(3)/12. ULL=(HOURLT-12.)/6. ZA=CL(0)-2.*CL(2) RATCH=ZA+ULL*(CL(1)-3.4*CL(3)+ULL*(CL(2)+ULL*CL(3))) RETURN END C real function fmlt(ut,xlat,xlong,day) C Calculats geomagnetic local time in hours for given C universal time UT(sec) C XLAT/XLONG=geodetic Latitude/Longitude (deg.) C DAY= day of year C D. Bilitza, COSPAR-1980, paper 7.3.1 C rad=57.29578 flat=xlat/rad flong=xlong/rad delta=0.409207*sin((day-80.0)*3.14159/184.0) beta=(180.-ut/240.0)/rad px=cos(flat)*cos(flong) py=cos(flat)*sin(flong) pz=sin(flat) sx=cos(delta)*cos(beta) sy=cos(delta)*sin(beta) sz=sin(delta) p1=0.35117*px-0.91483*py-0.19937*pz p2=0.93358*px+0.35837*py s1=0.35117*sx-0.91483*sy-0.19937*sz s2=0.93358*sx+0.35837*sy thetp=atan2(p2,p1) thets=atan2(s2,s1) fmlt=(thetp-thets+3.141593)/.2617994 afmlt=anint(fmlt*10.) fmlt=afmlt/10. if(fmlt.ge.24.) fmlt=fmlt-24. if(fmlt.lt.0.) fmlt=fmlt+24. return end C C --------------------------------------------------------------------- C SUBROUTINE CONVER(rga,rgo,rgma) C This subroutine converts a geographic latitude and longitude C location to a corrected geomagnetic latitude. C C INPUT: C geographic latitude -90. to +90. C geographic longitude 0. to 360. positive east from Greenwich. C C OUTPUT: C corrected geomagnetic latitude -90. to +90. DIMENSION CORMAG(20,91) DATA ((CORMAG(i,j),i=1,20),j=1,31)/ +163.68,163.68,163.68,163.68,163.68,163.68, +163.68,163.68,163.68,163.68,163.68,163.68,163.68,163.68, +163.68,163.68,163.68,163.68,163.68,163.68,162.60,163.12, +163.64,164.18,164.54,164.90,165.16,165.66,166.00,165.86, +165.20,164.38,163.66,162.94,162.42,162.00,161.70,161.70, +161.80,162.14,161.20,162.18,163.26,164.44,165.62,166.60, +167.42,167.80,167.38,166.82,166.00,164.66,163.26,162.16, +161.18,160.40,159.94,159.80,159.98,160.44,159.80,161.14, +162.70,164.50,166.26,167.90,169.18,169.72,169.36,168.24, +166.70,164.80,162.90,161.18,159.74,158.60,157.94,157.80, +157.98,158.72,158.40,160.10,162.02,164.28,166.64,169.00, +170.80,171.72,171.06,169.46,167.10,164.64,162.18,160.02, +158.20,156.80,156.04,155.80,156.16,157.02,157.00,158.96, +161.24,163.86,166.72,169.80,172.42,173.72,172.82,170.34, +167.30,164.22,161.34,158.74,156.60,155.00,154.08,153.90, +154.36,155.36,155.50,157.72,160.36,163.32,166.60,170.20, +173.70,175.64,174.18,170.80,167.10,163.56,160.24,157.36, +154.96,153.10,152.08,151.92,152.46,153.76,154.10,156.52, +159.36,162.52,166.24,170.30,174.62,177.48,175.04,170.82, +166.60,162.70,159.02,155.88,153.22,151.20,150.08,149.92, +150.64,152.20,152.80,155.32,158.28,161.70,165.58,170.00, +174.84,178.46,175.18,170.38,165.80,161.64,157.80,154.38, +151.52,149.30,148.18,148.02,148.92,150.60,151.40,154.08, +157.18,160.68,164.78,169.40,174.34,177.44,174.28,169.44, +164.70,160.34,156.30,152.78,149.72,147.40,146.18,146.04, +147.12,149.04,150.10,152.88,156.00,159.58,163.78,168.50, +173.28,175.60,172.86,168.14,163.40,158.98,154.88,151.10, +147.98,145.50,144.18,144.14,145.40,147.48,148.80,151.68, +154.88,158.48,162.68,167.40,171.76,173.60,171.12,166.68, +162.00,157.48,153.28,149.50,146.18,143.50,142.18,142.24, +143.68,145.98,147.50,150.54,153.68,157.28,161.42,166.10, +170.10,171.48,169.22,164.98,160.40,155.88,151.68,147.80, +144.34,141.60,140.18,140.26,141.98,144.62,146.30,149.34, +152.48,155.98,160.08,164.60,168.34,169.38,167.20,163.18, +158.60,154.18,149.98,146.02,142.54,139.70,138.18,138.46, +140.26,143.16,145.10,148.14,151.18,154.60,158.68,163.10, +166.48,167.28,165.18,161.32,156.90,152.48,148.28,144.32, +140.74,137.80,136.22,136.48,138.64,141.76,143.90,146.98, +149.98,153.30,157.24,161.40,164.52,165.16,162.86,159.42, +155.00,150.68,146.48,142.52,138.94,135.90,134.22,134.68, +137.02,140.40,142.70,145.84,148.76,151.92,155.74,159.70, +162.52,162.96,160.98,157.42,153.10,148.84,144.68,140.82, +137.20,134.00,132.32,132.80,135.42,139.10,141.60,144.74, +147.46,150.52,154.20,158.00,160.46,160.76,158.86,155.36, +151.20,146.94,142.88,139.02,135.40,132.10,130.32,131.00, +133.80,137.74,140.50,143.58,146.24,149.12,152.60,156.20, +158.40,158.66,156.76,153.36,149.30,145.04,141.08,137.30, +133.60,130.30,128.42,129.12,132.28,136.44,139.30,142.48, +144.94,147.64,150.48,154.30,156.34,156.36,154.56,151.26, +147.30,143.14,139.20,135.50,131.90,128.40,126.52,127.32, +130.76,135.18,138.20,141.28,143.72,146.24,149.26,152.40, +154.24,154.16,152.36,149.16,145.30,141.24,137.30,133.70, +130.10,126.60,124.62,125.54,129.16,133.92,137.10,140.18, +142.42,144.66,147.62,150.50,152.18,151.96,150.16,147.10, +143.30,139.24,135.50,131.90,128.36,124.80,122.72,123.74, +127.64,132.62,135.90,139.02,141.12,143.18,145.92,148.60, +149.98,149.76,148.04,145.00,141.20,137.30,133.60,130.10, +126.60,123.00,120.86,121.96,126.12,131.36,134.80,137.88, +139.80,141.68,144.08,146.60,147.88,147.56,145.84,142.90, +139.20,135.30,131.70,128.28,124.86,121.30,118.96,120.18, +124.70,130.16,133.60,136.72,138.48,140.10,142.38,144.60, +145.72,145.34,143.64,140.80,137.10,133.30,129.72,126.48, +123.10,119.50,117.16,118.48,123.18,128.86,132.40,135.42, +137.08,138.50,140.54,142.60,143.52,143.06,141.44,138.70, +135.10,131.30,127.82,124.58,121.40,117.70,115.26,116.70, +121.66,127.60,131.20,134.22,135.66,136.82,138.70,140.60, +141.36,140.86,139.24,136.50,133.00,129.30,125.92,122.78, +119.60,116.00,113.40,114.92,120.16,126.30,130.00,132.92, +134.24,135.14,136.80,138.60,139.16,138.64,137.12,134.40, +130.90,127.20,123.92,120.96,117.90,114.20,111.56,113.12, +118.64,124.90,128.70,131.56,132.74,133.44,134.90,136.50, +137.00,136.36,134.82,132.30,128.70,125.16,121.94,119.06, +116.10,112.50,109.70,111.42,117.14,123.60,127.30,130.16, +131.22,131.66,133.00,134.50,134.80,134.14,132.62,130.14, +126.60,123.06,119.94,117.16,114.30,110.70,107.80,109.64, +115.62,122.24,125.90,128.76,129.62,129.96,131.06,132.40, +132.60,131.86,130.42,128.00,124.50,120.96,117.96,115.26, +112.54,108.90,105.94,107.86,114.02,120.84/ DATA ((CORMAG(i,j),i=1,20),j=32,61)/ +124.05,126.79, +127.55,127.83,128.90,130.21,130.41,129.71,128.33,125.96, +122.49,118.96,115.97,113.26,110.52,106.89,104.01,106.00, +112.21,119.06,122.19,124.82,125.48,125.69,126.73,128.03, +128.22,127.55,126.23,123.92,120.47,116.97,113.97,111.26, +108.50,104.89,102.08,104.14,110.41,117.29,120.34,122.85, +123.40,123.56,124.57,125.84,126.03,125.40,124.14,121.88, +118.46,114.97,111.98,109.26,106.48,102.88,100.15,102.28, +108.60,115.51,118.49,120.88,121.33,121.42,122.40,123.65, +123.84,123.24,122.04,119.83,116.45,112.97,109.98,107.26, +104.46,100.87,098.22,100.42,106.79,113.74,116.63,118.91, +119.26,119.29,120.24,121.47,121.65,121.09,119.95,117.79, +114.43,110.98,107.99,105.26,102.44,098.87,096.29,098.56, +104.98,111.96,114.78,116.94,117.19,117.15,118.07,119.28, +119.46,118.93,117.86,115.75,112.42,108.98,106.00,103.26, +100.42,096.86,094.36,096.70,103.18,110.19,112.93,114.97, +115.12,115.02,115.91,117.09,117.27,116.78,115.76,113.71, +110.41,106.98,104.00,101.26,098.40,094.85,092.43,094.84, +101.37,108.41,111.07,113.00,113.04,112.88,113.74,114.91, +115.08,114.62,113.67,111.67,108.39,104.99,102.01,099.26, +096.38,092.85,090.51,092.97,099.56,106.64,109.22,111.03, +110.97,110.75,111.58,112.72,112.89,112.47,111.57,109.63, +106.38,102.99,100.01,097.26,094.36,090.84,088.58,091.11, +097.75,104.86,107.37,109.06,108.90,108.61,109.41,110.53, +110.70,110.31,109.48,107.59,104.37,100.99,098.02,095.26, +092.34,088.83,086.65,089.25,095.95,103.09,105.51,107.09, +106.83,106.48,107.25,108.35,108.51,108.16,107.39,105.55, +102.35,099.00,096.03,093.26,090.32,086.83,084.72,087.39, +094.14,101.31,103.66,105.12,104.76,104.34,105.08,106.16, +106.32,106.00,105.29,103.50,100.34,097.00,094.03,091.26, +088.30,084.82,082.79,085.53,092.33,099.54,101.81,103.15, +102.68,102.21,102.92,103.97,104.13,103.85,103.20,101.46, +098.33,095.00,092.04,089.26,086.28,082.81,080.86,083.67, +090.52,097.76,099.95,101.18,100.61,100.07,100.75,101.79, +101.94,101.69,101.10,099.42,096.31,093.01,090.04,087.26, +084.26,080.81,078.93,081.81,088.72,095.99,098.10,099.21, +098.54,097.94,098.59,099.60,099.75,099.54,099.01,097.38, +094.30,091.01,088.05,085.26,082.24,078.80,077.00,079.95, +086.91,094.21,096.25,097.24,096.47,095.81,096.43,097.41, +097.56,097.39,096.92,095.34,092.29,089.01,086.06,083.26, +080.22,076.79,075.07,078.09,085.10,092.43,094.39,095.27, +094.40,093.67,094.26,095.23,095.37,095.23,094.82,093.30, +090.27,087.02,084.06,081.26,078.20,074.79,073.14,076.23, +083.30,090.66,092.54,093.30,092.32,091.54,092.10,093.04, +093.18,093.08,092.73,091.26,088.26,085.02,082.07,079.26, +076.18,072.78,071.21,074.37,081.49,088.88,090.69,091.33, +090.25,089.40,089.93,090.85,090.99,090.92,090.63,089.21, +086.25,083.02,080.07,077.26,074.16,070.77,069.28,072.51, +079.68,087.11,088.83,089.36,088.18,087.27,087.77,088.67, +088.80,088.77,088.54,087.17,084.23,081.03,078.08,075.26, +072.14,068.77,067.35,070.65,077.87,085.33,086.98,087.39, +086.11,085.13,085.60,086.48,086.61,086.61,086.45,085.13, +082.22,079.03,076.09,073.26,070.12,066.76,065.42,068.79, +076.07,083.56,085.13,085.42,084.04,083.00,083.44,084.29, +084.42,084.46,084.35,083.09,080.21,077.03,074.09,071.26, +068.10,064.75,063.49,066.93,074.26,081.78,083.27,083.45, +081.96,080.86,081.27,082.11,082.23,082.30,082.26,081.05, +078.19,075.04,072.10,069.26,066.08,062.75,061.57,065.06, +072.45,080.01,081.42,081.48,079.89,078.73,079.11,079.92, +080.04,080.15,080.16,079.01,076.18,073.04,070.10,067.26, +064.06,060.74,059.64,063.20,070.64,078.23,079.57,079.51, +077.82,076.59,076.94,077.73,077.85,077.99,078.07,076.97, +074.17,071.04,068.11,065.26,062.04,058.73,057.71,061.34, +068.84,076.46,077.71,077.54,075.75,074.46,074.78,075.55, +075.66,075.84,075.98,074.93,072.15,069.05,066.12,063.26, +060.02,056.73,055.78,059.48,067.03,074.68,075.86,075.57, +073.68,072.32,072.61,073.36,073.47,073.68,073.88,072.88, +070.14,067.05,064.12,061.26,058.00,054.72,053.85,057.62, +065.22,072.91,074.01,073.60,071.60,070.19,070.45,071.17, +071.28,071.53,071.79,070.84,068.13,065.05,062.13,059.26, +055.98,052.71,051.92,055.76,063.41,071.13,072.15,071.63, +069.53,068.05,068.28,068.99,069.09,069.37,069.69,068.80, +066.11,063.06,060.13,057.26,053.96,050.71,049.99,053.90, +061.61,069.36,070.30,069.66,067.46,065.92,066.12,066.80, +066.90,067.22,067.60,066.76,064.10,061.06,058.14,055.26, +051.94,048.70,048.06,052.04,059.80,067.58/ DATA ((CORMAG(i,j),i=1,20),j=62,91)/ +067.70,067.06, +065.08,063.72,063.98,064.60,064.80,065.12,065.60,064.86, +062.40,059.26,056.24,053.18,049.84,046.60,046.12,050.12, +057.52,064.80,064.90,064.42,062.70,061.62,061.78,062.40, +062.60,063.04,063.58,063.00,060.60,057.46,054.42,051.18, +047.70,044.60,044.22,048.02,055.06,061.92,062.10,061.72, +060.32,059.50,059.68,060.20,060.46,060.94,061.58,061.00, +058.70,055.66,052.52,049.18,045.60,042.50,042.22,046.00, +052.60,058.98,059.20,059.18,058.12,057.32,057.48,058.00, +058.30,058.84,059.48,059.04,056.90,053.86,050.62,047.10, +043.50,040.50,040.28,043.98,050.22,056.18,056.40,056.64, +055.84,055.20,055.38,055.80,056.16,056.84,057.48,057.04, +055.10,052.06,048.70,045.10,041.40,038.40,038.28,041.88, +047.94,053.44,053.70,054.14,053.56,053.10,053.24,053.70, +054.06,054.74,055.38,055.14,053.20,050.26,046.80,043.10, +039.34,036.40,036.38,039.96,045.56,050.84,051.10,051.70, +051.36,051.00,051.14,051.50,051.96,052.64,053.38,053.08, +051.30,048.36,044.90,041.02,037.24,034.40,034.38,037.86, +043.28,048.20,048.50,049.26,049.18,048.90,049.04,049.40, +049.86,050.64,051.28,051.08,049.40,046.46,042.98,039.02, +035.14,032.40,032.48,035.72,041.00,045.70,046.00,046.96, +046.98,046.80,046.94,047.30,047.76,048.54,049.28,049.08, +047.40,044.56,041.08,037.02,033.14,030.40,030.58,033.84, +038.72,043.20,043.50,044.62,044.80,044.80,044.94,045.20, +045.76,046.54,047.18,046.98,045.50,042.66,039.08,035.02, +031.14,028.40,028.58,031.82,036.52,040.80,041.20,042.32, +042.54,042.70,042.84,043.20,043.66,044.44,045.08,044.98, +043.50,040.76,037.08,033.04,029.04,026.40,026.68,029.82, +034.34,038.40,038.80,040.12,040.60,040.70,040.84,041.10, +041.62,042.34,042.98,042.88,041.50,038.76,035.18,031.04, +027.14,024.50,024.78,027.70,032.14,036.06,036.50,037.88, +038.50,038.68,038.84,039.10,039.56,040.34,040.88,040.82, +039.40,036.76,033.18,029.12,025.14,022.50,022.88,025.90, +029.96,033.86,034.30,035.68,036.42,036.68,036.84,037.10, +037.56,038.24,038.88,038.72,037.40,034.76,031.18,027.12, +023.14,020.60,020.98,023.90,027.88,031.66,032.10,033.58, +034.32,034.68,034.84,035.10,035.56,036.24,036.78,036.62, +035.30,032.72,029.18,025.14,021.24,018.70,019.08,021.90, +025.88,029.42,029.90,031.48,032.32,032.68,032.84,033.10, +033.56,034.22,034.68,034.42,033.20,030.72,027.28,023.22, +019.34,016.80,017.24,020.00,023.78,027.32,027.70,029.38, +030.24,030.68,030.94,031.20,031.66,032.22,032.58,032.32, +031.10,028.62,025.28,021.32,017.48,015.00,015.38,018.18, +021.80,025.22,025.70,027.28,028.24,028.78,029.04,029.30, +029.66,030.22,030.50,030.22,029.00,026.62,023.30,019.42, +015.64,013.10,013.54,016.28,019.80,023.12,023.60,025.24, +026.24,026.78,027.14,027.40,027.76,028.22,028.40,028.12, +026.80,024.52,021.30,017.52,013.78,011.30,011.74,014.48, +017.90,021.12,021.60,023.24,024.34,024.88,025.24,025.50, +025.86,026.22,026.40,025.98,024.70,022.48,019.40,015.72, +012.04,009.50,009.94,012.58,016.02,019.12,019.60,021.24, +022.34,022.98,023.34,023.70,024.00,024.30,024.40,023.88, +022.60,020.48,017.52,014.00,010.34,007.80,008.18,010.88, +014.22,017.18,017.60,019.34,020.44,021.16,021.54,021.90, +022.16,022.40,022.32,021.78,020.60,018.48,015.62,012.20, +008.68,006.00,006.44,009.18,012.42,015.28,015.80,017.44, +018.54,019.26,019.74,020.10,020.30,020.50,020.32,019.72, +018.50,016.54,013.84,010.68,007.14,004.40,004.74,007.58, +010.74,013.48,014.00,015.54,016.74,017.46,017.94,018.30, +018.50,018.58,018.32,017.72,016.50,014.64,012.24,009.18, +005.84,002.90,003.30,006.16,009.14,011.84,012.30,013.78, +014.94,015.66,016.24,016.50,016.70,016.70,016.42,015.78, +014.60,012.90,010.66,007.86,004.88,001.60,001.72,004.96, +007.84,010.24,010.70,012.14,013.24,013.96,014.44,014.80, +014.90,014.88,014.52,013.92,012.80,011.30,009.28,006.94, +004.32,001.80,001.94,004.34,006.78,008.94,009.40,010.58, +011.64,012.36,012.74,013.10,013.20,013.08,012.72,012.12, +011.10,009.86,008.30,006.50,004.60,003.10,003.16,004.50, +006.20,007.90,008.40,009.42,010.14,010.76,011.14,011.40, +011.40,011.38,011.02,010.46,009.70,008.72,007.64,006.46, +005.42,004.60,004.70,005.34,006.24,007.36,007.90,008.46, +008.92,009.28,009.54,009.70,009.70,009.68,009.42,009.06, +008.60,008.08,007.56,007.02,006.56,006.30,006.30,006.52, +006.96,007.38,008.15,008.15,008.15,008.15,008.15,008.15, +008.15,008.15,008.15,008.15,008.15,008.15,008.15,008.15, +008.15,008.15,008.15,008.15,008.15,008.15/ C Data Input rlan = rga rlo = rgo C From "normal" geographic latitude C to angle from South Pole. rla = rlan + 90 IF (rlo .EQ. 360) THEN rlo = 0 END IF C PROXIMITY C coefficients of the latitudinal points LA1 = (INT(rla/2)+1) LA2 = LA1 + 1 if(la2.gt.91) la2=91 C coefficients of the longitudinal points LO1 = (INT(rlo/18)+1) corr LO2 = LO1 + 1 LO2 = MOD(LO1,20) + 1 C Four points of Geomagnetic Coordinates gm1 = CORMAG(LO1,LA1) gm2 = CORMAG(LO1,LA2) gm3 = CORMAG(LO2,LA1) gm4 = CORMAG(LO2,LA2) C latitudinal points X1 = ABS(rla - (INT(rla))) X2 = 2. - X1 C longitudinal points Y1 = ABS(rlo - (INT(rlo))) Y2 = 18. - Y1 C X AND Y VALUES x = X1 / (X1 + X2) y = Y1 / (Y1 + Y2) C INTERPOLATION gmla = gm1 * (1 - x) * (1 - y) + gm2 * (1 - y) * (x) + gm3 * (y) 1 * (1 - x) + gm4 * (x) * (y) C OUTPUT OF THE PROGRAM C From corrected geomagnetic latitude from North Pole C to "normal" geomagnetic latitude. rgma = 90. - gmla END c C-------------------------------------------------- real function peakh(foE,foF2,M3000) real MF,M3000 sqM=M3000*M3000 MF=M3000*sqrt((0.0196*sqM+1.)/(1.2967*sqM-1.0)) If(foE.ge.1.0E-30) then ratio=foF2/foE ratio=djoin(ratio,1.75,20.0,ratio-1.75) dM=0.253/(ratio-1.215)-0.012 else dM=-0.012 endif peakh=1490.0*MF/(M3000+dM)-176.0 return end C----------------------------------------------------------------- real function djoin(f1,f2,alpha,x) real f1,f2,alpha,x,ee,fexp ee=fexp(alpha*x) djoin=(f1*ee+f2)/(ee+1.0) return end C----------------------------------------------------------------- real function fexp(a) real a if(a.gt.80.0) then fexp=5.5406E34 return endif if(a.lt.-80.0D0) then fexp=1.8049E-35 return endif fexp=exp(a) return end C-----------------------------------------------------------------