--- DarthVader/OrbitalInfo/src/igrf_sub.for 2007/04/27 10:35:35 1.2 +++ DarthVader/OrbitalInfo/src/igrf_sub.for 2014/01/17 09:54:31 1.9 @@ -1,40 +1,40 @@ - subroutine igrf_sub(xlat,xlong,year,height, - & xl,icode,dip,dec) +c subroutine igrf_sub(xlat,xlong,year,height, +c & xl,icode,dip,dec) c---------------------------------------------------------------- c INPUT: -c xlat geodatic latitude in degrees -c xlong geodatic longitude in degrees -c year decimal year (year+month/12.0-0.5 or year+day-of-year/365 -c or 366 if leap year) -c height height in km +c xlat geodatic latitude in degrees +c xlong geodatic longitude in degrees +c year decimal year (year+month/12.0-0.5 or year+day-of-year/365 +c or 366 if leap year) +c height height in km c OUTPUT: -c xl L value -c icode =1 L is correct; =2 L is not correct; -c =3 an approximation is used -c dip geomagnetic inclination in degrees -c dec geomagnetic declination in degress +c xl L value +c icode =1 L is correct; =2 L is not correct; +c =3 an approximation is used +c dip geomagnetic inclination in degrees +c dec geomagnetic declination in degress c---------------------------------------------------------------- - - REAL LATI,LONGI - COMMON/GENER/ UMR,ERA,AQUAD,BQUAD - SAVE /GENER/ -C - CALL INITIZE - ibbb=0 - ALOG2=ALOG(2.) - ISTART=1 - lati=xlat - longi=xlong +c +c REAL LATI,LONGI +c COMMON/GENER/ UMR,ERA,AQUAD,BQUAD +c SAVE /GENER/ +C +c CALL INITIZE +c ibbb=0 +c ALOG2=ALOG(2.) +c ISTART=1 +c lati=xlat +c longi=xlong c C----------------CALCULATE PROFILES----------------------------------- c - CALL FELDCOF(YEAR,DIMO) - CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS) - CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1) - DIP=ASIN(BDOWN/BABS)/UMR - DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR - RETURN - END +c CALL FELDCOF(YEAR,DIMO) +c CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS) +c CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1) +c DIP=ASIN(BDOWN/BABS)/UMR +c DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR +c RETURN +c END c c C SHELLIG.FOR, Version 2.0, January 1992 @@ -191,7 +191,7 @@ C B0 MAGNETIC FIELD STRENGTH IN GAUSS C----------------------------------------------------------------------- DIMENSION V(3),U(3,3),P(8,100),SP(3) - COMMON X(3),H(144) + COMMON X(3),H(196) COMMON/FIDB0/ SP SAVE /FIDB0/ COMMON/GENER/ UMR,ERA,AQUAD,BQUAD @@ -308,7 +308,8 @@ RQ=R*R FF=SQRT(1.+3.*ZQ/RQ) RADIK=B0-((D2*T+D1)*T+D0)*R*RQ*FF - IF(R-RMAX)44,44,45 + IF((R-RMAX).le.0.) goto 44 + IF((R-RMAX).gt.0.) goto 45 45 ICODE=2 RADIK=RADIK-12.*(R-RMAX)**2 44 IF(RADIK+RADIK.LE.ORADIK) GOTO 10 @@ -341,24 +342,25 @@ C-- because 1E-38 is the minimal allowable arg. for ALOG in our envir. C-- D. Bilitza, Nov 87. C -11 FI=0.5*ABS(FI)/SQRT(B0)+1E-12 +11 FI=0.5*ABS(FI)/SQRT(B0)+1E-12 C C*****COMPUTE L FROM B AND I. SAME AS CARMEL IN INVAR. C C-- Correct dipole moment is used here. D. Bilitza, Nov 87. C DIMOB0=DIMO/B0 - arg1=alog(FI) - arg2=alog(DIMOB0) -c arg = FI*FI*FI/DIMOB0 -c if(abs(arg).gt.88.0) arg=88.0 + arg1=alog(FI) + arg2=alog(DIMOB0) +c arg = FI*FI*FI/DIMOB0 +c if(abs(arg).gt.88.0) arg=88.0 XX=3*arg1-arg2 IF(XX.GT.23.0) GOTO 776 IF(XX.GT.11.7) GOTO 775 IF(XX.GT.+3.0) GOTO 774 IF(XX.GT.-3.0) GOTO 773 IF(XX.GT.-22.) GOTO 772 - 771 GG=3.33338E-1*XX+3.0062102E-1 +c 771 GG=3.33338E-1*XX+3.0062102E-1 + GG=3.33338E-1*XX+3.0062102E-1 GOTO777 772 GG=((((((((-8.1537735E-14*XX+8.3232531E-13)*XX+1.0066362E-9)*XX+ 18.1048663E-8)*XX+3.2916354E-6)*XX+8.2711096E-5)*XX+1.3714667E-3)* @@ -392,7 +394,7 @@ C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG * C******************************************************************* DIMENSION P(7),U(3,3) - COMMON XI(3),H(144) + COMMON XI(3),H(196) C*****XM,YM,ZM ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES ZM=P(3) FLI=P(1)*P(1)+P(2)*P(2)+1E-15 @@ -474,11 +476,11 @@ C POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST C AND DOWNWARD. C----------------------------------------------------------------------- - DIMENSION V(3),B(3),G(144) + DIMENSION V(3),B(3),G(196) CHARACTER*258 NAME INTEGER NMAX REAL TIME - COMMON XI(3),H(144) + COMMON XI(3),H(196) COMMON/MODEL/ G,NMAX,TIME,NAME SAVE/MODEL/ COMMON/GENER/ UMR,ERA,AQUAD,BQUAD @@ -528,7 +530,10 @@ Y=XI(2)*F Z=XI(3)*(F+F) I=I-2 - IF(I-1)5,4,2 +c print *,' I ',I + IF((I-1).lt.0) goto 5 + IF((I-1).eq.0) goto 4 + IF((I-1).gt.0) goto 2 2 DO 3 M=3,I,2 H(IL+M+1)=G(IL+M+1)+Z*H(IH+M+1)+X*(H(IH+M+3)-H(IH+M-1)) A -Y*(H(IH+M+2)+H(IH+M-2)) @@ -560,177 +565,192 @@ END C C - SUBROUTINE FELDCOF(YEAR,DIMO) + SUBROUTINE FELDCOF(YEAR,DIMO) C------------------------------------------------------------------------ -C DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS -C -C INPUT: YEAR DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO -C BE CALCULATED -C OUTPUT: DIMO GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED -C TO EARTH'S RADIUS) AT THE TIME (YEAR) -C D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771, -C (301)286-9536 NOV 1987. -C ### updated to IGRF-2000 version -dkb- 5/31/2000 -C ### updated to IGRF-2005 version -dkb- 3/24/2000 +C DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS +C +C INPUT: YEAR DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO +C BE CALCULATED +C OUTPUT: DIMO GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED +C TO EARTH'S RADIUS) AT THE TIME (YEAR) +C D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771, +C (301)286-9536 NOV 1987. +C ### updated to IGRF-2000 version -dkb- 5/31/2000 +C ### updated to IGRF-2005 version -dkb- 3/24/2000 C----------------------------------------------------------------------- - CHARACTER*258 FIL1, FIL2 - CHARACTER*258 FILMOD -C ### FILMOD, DTEMOD arrays +1 -c DIMENSION GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14) - DIMENSION GH1(144),GH2(120),GHA(144),FILMOD(3),DTEMOD(3) - DOUBLE PRECISION X,F0,F - INTEGER L1,L2,L3 - INTEGER NMAX - REAL TIME - CHARACTER *258 P1,P2,P3 - COMMON/PPATH/ L1,L2,L3,P1, P2, P3 - SAVE/PPATH/ - COMMON/MODEL/ GH1,NMAX,TIME,FIL1 - SAVE/MODEL/ - COMMON/GENER/ UMR,ERAD,AQUAD,BQUAD - SAVE/GENER/ -C ### updated to 2005 -C CHARACTER COEFPATH*80, COEF1*80, COEF2*80, COEF3*80 - -c COEFPATH = 'OrbitalInfo/src/' -c COEF1 = 'dgrf00.dat' -c COEF2 = 'igrf05.dat' -c COEF3 = 'igrf05s.dat' -c COEF1 = COEFPATH(1:16)//COEF1 -c COEF2 = COEFPATH(1:16)//COEF2 -c COEF3 = COEFPATH(1:16)//COEF3 -c FILMOD(1) = COEF1 -c FILMOD(2) = COEF2 -c FILMOD(3) = COEF3 + CHARACTER*258 FIL1, FIL2 + CHARACTER*258 FILMOD +C ### FILMOD, DTEMOD arrays +1 +c DIMENSION GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14) + DIMENSION GH1(196),GH2(196),GHA(196),FILMOD(3),DTEMOD(3) + DOUBLE PRECISION X,F0,F + INTEGER L1,L2,L3 + INTEGER NMAX + REAL TIME + CHARACTER *258 P1,P2,P3 + COMMON/PPATH/ L1,L2,L3,P1, P2, P3 + SAVE/PPATH/ + COMMON/MODEL/ GH1,NMAX,TIME,FIL1 + SAVE/MODEL/ + COMMON/GENER/ UMR,ERAD,AQUAD,BQUAD + SAVE/GENER/ +C ### updated to 2005 +C CHARACTER COEFPATH*80, COEF1*80, COEF2*80, COEF3*80 + +c COEFPATH = 'OrbitalInfo/src/' +c COEF1 = 'dgrf00.dat' +c COEF2 = 'igrf05.dat' +c COEF3 = 'igrf05s.dat' +c COEF1 = COEFPATH(1:16)//COEF1 +c COEF2 = COEFPATH(1:16)//COEF2 +c COEF3 = COEFPATH(1:16)//COEF3 +c FILMOD(1) = COEF1 +c FILMOD(2) = COEF2 +c FILMOD(3) = COEF3 +c print *, "qui" FILMOD(1) = P1(1:L1) FILMOD(2) = P2(1:L2) FILMOD(3) = P3(1:L3) -c FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat' -c FILMOD(2) = 'OrbitalInfo/src/igrf05.dat' -c FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat' +c print *, "qua" +c FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat' +c FILMOD(2) = 'OrbitalInfo/src/igrf05.dat' +c FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat' c WRITE(*,*) FILMOD(1) c WRITE(*,*) FILMOD(2) c WRITE(*,*) FILMOD(3) c DATA FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/ - DATA DTEMOD / 2000., 2005., 2010./ -c -c DATA FILMOD /'dgrf45.dat', 'dgrf50.dat', + DATA DTEMOD / 2005., 2010., 2015./ +c +c DATA FILMOD /'dgrf45.dat', 'dgrf50.dat', c 1 'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat', c 2 'dgrf70.dat', 'dgrf75.dat', 'dgrf80.dat', c 3 'dgrf85.dat', 'dgrf90.dat', 'dgrf95.dat', c 4 'dgrf00.dat','igrf05.dat','igrf05s.dat'/ -c DATA DTEMOD / 1945., 1950., 1955., 1960., 1965., 1970., +c DATA DTEMOD / 1945., 1950., 1955., 1960., 1965., 1970., c 1 1975., 1980., 1985., 1990., 1995., 2000.,2005.,2010./ -C -C ### numye = numye + 1 ; is number of years represented by IGRF -C -c NUMYE=13 - NUMYE=2 - -C -C IS=0 FOR SCHMIDT NORMALIZATION IS=1 GAUSS NORMALIZATION -C IU IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS -C - IU = 10 - IS = 0 -C-- DETERMINE IGRF-YEARS FOR INPUT-YEAR - TIME = YEAR - IYEA = INT(YEAR/5.)*5 -c L = (IYEA - 1945)/5 + 1 - L = (IYEA - 2000)/5 + 1 - IF(L.LT.1) L=1 - IF(L.GT.NUMYE) L=NUMYE - DTE1 = DTEMOD(L) - FIL1 = FILMOD(L) - DTE2 = DTEMOD(L+1) - FIL2 = FILMOD(L+1) -C-- GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS - CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER) - IF (IER .NE. 0) STOP - CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER) - IF (IER .NE. 0) STOP -C-- DETERMINE IGRF COEFFICIENTS FOR YEAR - IF (L .LE. NUMYE-1) THEN - CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2, - 1 NMAX2, GH2, NMAX, GHA) - ELSE - CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2, - 1 GH2, NMAX, GHA) - ENDIF -C-- DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G - F0=0.D0 - DO 1234 J=1,3 - F = GHA(J) * 1.D-5 - F0 = F0 + F * F -1234 CONTINUE - DIMO = DSQRT(F0) - - GH1(1) = 0.0 - I=2 - F0=1.D-5 - IF(IS.EQ.0) F0=-F0 - SQRT2=SQRT(2.) +C +C ### numye = numye + 1 ; is number of years represented by IGRF +C +c NUMYE=13 + NUMYE=2 +c print *, "quo" + +C +C IS=0 FOR SCHMIDT NORMALIZATION IS=1 GAUSS NORMALIZATION +C IU IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS +C + IU = 10 + IS = 0 +C-- DETERMINE IGRF-YEARS FOR INPUT-YEAR + TIME = YEAR + IYEA = INT(YEAR/5.)*5 +c L = (IYEA - 1945)/5 + 1 + L = (IYEA - 2000)/5 + 1 + IF(L.LT.1) L=1 + IF(L.GT.NUMYE) L=NUMYE + DTE1 = DTEMOD(L) + FIL1 = FILMOD(L) + DTE2 = DTEMOD(L+1) + FIL2 = FILMOD(L+1) +c WRITE(*,*) FIL1 +c WRITE(*,*) FIL2 +c print *, "que" +C-- GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS + CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER) + IF (IER .NE. 0) STOP +c print *, "quessss" + CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER) + IF (IER .NE. 0) STOP +c print *, "quj" +C-- DETERMINE IGRF COEFFICIENTS FOR YEAR + IF (L .LE. NUMYE-1) THEN + CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2, + 1 NMAX2, GH2, NMAX, GHA) + ELSE + CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2, + 1 GH2, NMAX, GHA) + ENDIF +c print *, "quw" +C-- DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G + F0=0.D0 + DO 1234 J=1,3 + F = GHA(J) * 1.D-5 + F0 = F0 + F * F + 1234 CONTINUE + DIMO = REAL(DSQRT(F0)) + + GH1(1) = 0.0 + I=2 + F0=1.D-5 + IF(IS.EQ.0) F0=-F0 + SQRT2=SQRT(2.) +c print *, "quq" + DO 9 N=1,NMAX - X = N - F0 = F0 * X * X / (4.D0 * X - 2.D0) - IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X - F = F0 * 0.5D0 - IF(IS.EQ.0) F = F * SQRT2 - GH1(I) = GHA(I-1) * F0 - I = I+1 - DO 9 M=1,N - F = F * (X + M) / (X - M + 1.D0) - IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M)) - GH1(I) = GHA(I-1) * F - GH1(I+1) = GHA(I) * F - I=I+2 -9 CONTINUE - RETURN - END -C -C - SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER) + X = N + F0 = F0 * X * X / (4.D0 * X - 2.D0) + IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X + F = F0 * 0.5D0 + IF(IS.EQ.0) F = F * SQRT2 + GH1(I) = GHA(I-1) * REAL(F0) + I = I+1 + DO 9 M=1,N + F = F * (X + M) / (X - M + 1.D0) + IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M)) + GH1(I) = GHA(I-1) * REAL(F) + GH1(I+1) = GHA(I) * REAL(F) + I=I+2 + 9 CONTINUE + RETURN + END +C +C + SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER) C =============================================================== C -C Version 1.01 -C -C Reads spherical harmonic coefficients from the specified -C file into an array. -C -C Input: -C IU - Logical unit number -C FSPEC - File specification -C -C Output: -C NMAX - Maximum degree and order of model -C ERAD - Earth's radius associated with the spherical -C harmonic coefficients, in the same units as -C elevation -C GH - Schmidt quasi-normal internal spherical -C harmonic coefficients -C IER - Error number: = 0, no error -C = -2, records out of order -C = FORTRAN run-time error number -C -C A. Zunde -C USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 -C -C =============================================================== - - CHARACTER FSPEC*(*), FOUT*258 - DIMENSION GH(*) -C --------------------------------------------------------------- -C Open coefficient file. Read past first header record. -C Read degree and order of model and Earth's radius. -C --------------------------------------------------------------- +C Version 1.01 +C +C Reads spherical harmonic coefficients from the specified +C file into an array. +C +C Input: +C IU - Logical unit number +C FSPEC - File specification +C +C Output: +C NMAX - Maximum degree and order of model +C ERAD - Earth's radius associated with the spherical +C harmonic coefficients, in the same units as +C elevation +C GH - Schmidt quasi-normal internal spherical +C harmonic coefficients +C IER - Error number: = 0, no error +C = -2, records out of order +C = FORTRAN run-time error number +C +C A. Zunde +C USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 +C +C =============================================================== + + CHARACTER FSPEC*(*), FOUT*258 + DIMENSION GH(*) +C --------------------------------------------------------------- +C Open coefficient file. Read past first header record. +C Read degree and order of model and Earth's radius. +C --------------------------------------------------------------- WRITE(FOUT,667) FSPEC -c 667 FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12) -667 FORMAT(A258) +c 667 FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12) + 667 FORMAT(A258) +c print *," gui" OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999) +c print *," gua" READ (IU, *, IOSTAT=IER, ERR=999) +c print *," gue" READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD +c print *," guo" C --------------------------------------------------------------- C Read the coefficient file, arranged as follows: C @@ -752,8 +772,8 @@ C --------------------------------------------------------------- I = 0 - DO 2211 NN = 1, NMAX - DO 2233 MM = 0, NN + DO 2211 NN = 1, NMAX + DO 2233 MM = 0, NN READ (IU, *, IOSTAT=IER, ERR=999) N, M, G, H IF (NN .NE. N .OR. MM .NE. M) THEN IER = -2 @@ -767,9 +787,11 @@ ENDIF 2233 CONTINUE 2211 CONTINUE +c print *," guj" 999 CLOSE (IU) - +c print *," guw IER",IER + if ( IER .eq. -1 ) IER = 0 ! gfortran 4.1.2 bug workaround... hoping not to create problems with other versions RETURN END @@ -824,18 +846,18 @@ ELSE IF (NMAX1 .GT. NMAX2) THEN K = NMAX2 * (NMAX2 + 2) L = NMAX1 * (NMAX1 + 2) - DO 1122 I = K + 1, L + DO 1122 I = K + 1, L 1122 GH(I) = GH1(I) + FACTOR * (-GH1(I)) NMAX = NMAX1 ELSE K = NMAX1 * (NMAX1 + 2) L = NMAX2 * (NMAX2 + 2) - DO 1133 I = K + 1, L + DO 1133 I = K + 1, L 1133 GH(I) = FACTOR * GH2(I) NMAX = NMAX2 ENDIF - DO 1144 I = 1, K + DO 1144 I = 1, K 1144 GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I)) RETURN @@ -891,18 +913,18 @@ ELSE IF (NMAX1 .GT. NMAX2) THEN K = NMAX2 * (NMAX2 + 2) L = NMAX1 * (NMAX1 + 2) - DO 1155 I = K + 1, L + DO 1155 I = K + 1, L 1155 GH(I) = GH1(I) NMAX = NMAX1 ELSE K = NMAX1 * (NMAX1 + 2) L = NMAX2 * (NMAX2 + 2) - DO 1166 I = K + 1, L + DO 1166 I = K + 1, L 1166 GH(I) = FACTOR * GH2(I) NMAX = NMAX2 ENDIF - DO 1177 I = 1, K + DO 1177 I = 1, K 1177 GH(I) = GH1(I) + FACTOR * GH2(I) RETURN @@ -925,7 +947,7 @@ C ASTRONOMICAL UNION . C----------------------------------------------------------------- INTEGER TL1,TL2,TL3 - CHARACTER *258 TP1,TP2,TP3 + CHARACTER (len=258) TP1,TP2,TP3 INTEGER L1,L2,L3 CHARACTER *258 P1,P2,P3 COMMON/PPATH/ L1,L2,L3,P1, P2, P3