--- DarthVader/OrbitalInfo/src/igrf_sub.for 2007/04/27 10:35:35 1.2 +++ DarthVader/OrbitalInfo/src/igrf_sub.for 2009/08/04 14:01:27 1.3 @@ -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 @@ -358,7 +359,8 @@ 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)* @@ -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,189 @@ 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(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 +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 WRITE(*,*) FILMOD(1) -c WRITE(*,*) FILMOD(2) -c WRITE(*,*) FILMOD(3) +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', +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 print *, "que" +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 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 = 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) * 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) 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 +769,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 +784,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 +843,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 +910,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