| 1 |
c subroutine igrf_sub(xlat,xlong,year,height, |
| 2 |
c & xl,icode,dip,dec) |
| 3 |
c---------------------------------------------------------------- |
| 4 |
c INPUT: |
| 5 |
c xlat geodatic latitude in degrees |
| 6 |
c xlong geodatic longitude in degrees |
| 7 |
c year decimal year (year+month/12.0-0.5 or year+day-of-year/365 |
| 8 |
c or 366 if leap year) |
| 9 |
c height height in km |
| 10 |
c OUTPUT: |
| 11 |
c xl L value |
| 12 |
c icode =1 L is correct; =2 L is not correct; |
| 13 |
c =3 an approximation is used |
| 14 |
c dip geomagnetic inclination in degrees |
| 15 |
c dec geomagnetic declination in degress |
| 16 |
c---------------------------------------------------------------- |
| 17 |
c |
| 18 |
c REAL LATI,LONGI |
| 19 |
c COMMON/GENER/ UMR,ERA,AQUAD,BQUAD |
| 20 |
c SAVE /GENER/ |
| 21 |
C |
| 22 |
c CALL INITIZE |
| 23 |
c ibbb=0 |
| 24 |
c ALOG2=ALOG(2.) |
| 25 |
c ISTART=1 |
| 26 |
c lati=xlat |
| 27 |
c longi=xlong |
| 28 |
c |
| 29 |
C----------------CALCULATE PROFILES----------------------------------- |
| 30 |
c |
| 31 |
c CALL FELDCOF(YEAR,DIMO) |
| 32 |
c CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS) |
| 33 |
c CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1) |
| 34 |
c DIP=ASIN(BDOWN/BABS)/UMR |
| 35 |
c DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR |
| 36 |
c RETURN |
| 37 |
c END |
| 38 |
c |
| 39 |
c |
| 40 |
C SHELLIG.FOR, Version 2.0, January 1992 |
| 41 |
C |
| 42 |
C 11/01/91-DKB- SHELLG: lowest starting point for B0 search is 2 |
| 43 |
C 1/27/92-DKB- Adopted to IGRF-91 coeffcients model |
| 44 |
C 2/05/92-DKB- Reduce variable-names: INTER(P)SHC,EXTRA(P)SHC,INITI(ALI)ZE |
| 45 |
C 8/08/95-DKB- Updated to IGRF-45-95; new coeff. DGRF90, IGRF95, IGRF95S |
| 46 |
C 5/31/00-DKB- Updated to IGRF-45-00; new coeff.: IGRF00, IGRF00s |
| 47 |
C 3/24/05-DKB- Updated to IGRF-45-10; new coeff.: IGRF05, IGRF05s |
| 48 |
C 4/25/05-DKB- CALL FELDI and DO 1111 I=1,7 [Alexey Petrov] |
| 49 |
C |
| 50 |
C********************************************************************* |
| 51 |
C SUBROUTINES FINDB0, SHELLG, STOER, FELDG, FELDCOF, GETSHC, * |
| 52 |
C INTERSHC, EXTRASHC, INITIZE * |
| 53 |
C********************************************************************* |
| 54 |
C********************************************************************* |
| 55 |
C |
| 56 |
C |
| 57 |
SUBROUTINE FINDB0(STPS,BDEL,VALUE,BEQU,RR0) |
| 58 |
C-------------------------------------------------------------------- |
| 59 |
C FINDS SMALLEST MAGNETIC FIELD STRENGTH ON FIELD LINE |
| 60 |
C |
| 61 |
C INPUT: STPS STEP SIZE FOR FIELD LINE TRACING |
| 62 |
C COMMON/FIDB0/ |
| 63 |
C SP DIPOLE ORIENTED COORDINATES FORM SHELLG; P(1,*), |
| 64 |
C P(2,*), P(3,*) CLOSEST TO MAGNETIC EQUATOR |
| 65 |
C BDEL REQUIRED ACCURACY = [ B(LAST) - BEQU ] / BEQU |
| 66 |
C B(LAST) IS FIELD STRENGTH BEFORE BEQU |
| 67 |
C |
| 68 |
C OUTPUT: VALUE =.FALSE., IF BEQU IS NOT MINIMAL VALUE ON FIELD LINE |
| 69 |
C BEQU MAGNETIC FIELD STRENGTH AT MAGNETIC EQUATOR |
| 70 |
C RR0 EQUATORIAL RADIUS NORMALIZED TO EARTH RADIUS |
| 71 |
C BDEL FINAL ACHIEVED ACCURACY |
| 72 |
C-------------------------------------------------------------------- |
| 73 |
DIMENSION P(8,4),SP(3) |
| 74 |
LOGICAL VALUE |
| 75 |
COMMON/FIDB0/ SP |
| 76 |
SAVE /FIDB0/ |
| 77 |
C |
| 78 |
STEP=STPS |
| 79 |
IRUN=0 |
| 80 |
7777 IRUN=IRUN+1 |
| 81 |
IF(IRUN.GT.5) THEN |
| 82 |
VALUE=.FALSE. |
| 83 |
GOTO 8888 |
| 84 |
ENDIF |
| 85 |
C*********************FIRST THREE POINTS |
| 86 |
P(1,2)=SP(1) |
| 87 |
P(2,2)=SP(2) |
| 88 |
P(3,2)=SP(3) |
| 89 |
STEP=-SIGN(STEP,P(3,2)) |
| 90 |
CALL STOER(P(1,2),BQ2,R2) |
| 91 |
P(1,3)=P(1,2)+0.5*STEP*P(4,2) |
| 92 |
P(2,3)=P(2,2)+0.5*STEP*P(5,2) |
| 93 |
P(3,3)=P(3,2)+0.5*STEP |
| 94 |
CALL STOER(P(1,3),BQ3,R3) |
| 95 |
P(1,1)=P(1,2)-STEP*(2.*P(4,2)-P(4,3)) |
| 96 |
P(2,1)=P(2,2)-STEP*(2.*P(5,2)-P(5,3)) |
| 97 |
P(3,1)=P(3,2)-STEP |
| 98 |
CALL STOER(P(1,1),BQ1,R1) |
| 99 |
P(1,3)=P(1,2)+STEP*(20.*P(4,3)-3.*P(4,2)+P(4,1))/18. |
| 100 |
P(2,3)=P(2,2)+STEP*(20.*P(5,3)-3.*P(5,2)+P(5,1))/18. |
| 101 |
P(3,3)=P(3,2)+STEP |
| 102 |
CALL STOER(P(1,3),BQ3,R3) |
| 103 |
C******************INVERT SENSE IF REQUIRED |
| 104 |
IF(BQ3.LE.BQ1) GOTO 2 |
| 105 |
STEP=-STEP |
| 106 |
R3=R1 |
| 107 |
BQ3=BQ1 |
| 108 |
DO 1 I=1,5 |
| 109 |
ZZ=P(I,1) |
| 110 |
P(I,1)=P(I,3) |
| 111 |
1 P(I,3)=ZZ |
| 112 |
C******************INITIALIZATION |
| 113 |
2 STEP12=STEP/12. |
| 114 |
VALUE=.TRUE. |
| 115 |
BMIN=1.E4 |
| 116 |
BOLD=1.E4 |
| 117 |
C******************CORRECTOR (FIELD LINE TRACING) |
| 118 |
N=0 |
| 119 |
5555 P(1,3)=P(1,2)+STEP12*(5.*P(4,3)+8.*P(4,2)-P(4,1)) |
| 120 |
N=N+1 |
| 121 |
P(2,3)=P(2,2)+STEP12*(5.*P(5,3)+8.*P(5,2)-P(5,1)) |
| 122 |
C******************PREDICTOR (FIELD LINE TRACING) |
| 123 |
P(1,4)=P(1,3)+STEP12*(23.*P(4,3)-16.*P(4,2)+5.*P(4,1)) |
| 124 |
P(2,4)=P(2,3)+STEP12*(23.*P(5,3)-16.*P(5,2)+5.*P(5,1)) |
| 125 |
P(3,4)=P(3,3)+STEP |
| 126 |
CALL STOER(P(1,4),BQ3,R3) |
| 127 |
DO 1111 J=1,3 |
| 128 |
C DO 1111 I=1,8 |
| 129 |
DO 1111 I=1,7 |
| 130 |
1111 P(I,J)=P(I,J+1) |
| 131 |
B=SQRT(BQ3) |
| 132 |
IF(B.LT.BMIN) BMIN=B |
| 133 |
IF(B.LE.BOLD) THEN |
| 134 |
BOLD=B |
| 135 |
ROLD=1./R3 |
| 136 |
SP(1)=P(1,4) |
| 137 |
SP(2)=P(2,4) |
| 138 |
SP(3)=P(3,4) |
| 139 |
GOTO 5555 |
| 140 |
ENDIF |
| 141 |
IF(BOLD.NE.BMIN) THEN |
| 142 |
VALUE=.FALSE. |
| 143 |
ENDIF |
| 144 |
BDELTA=(B-BOLD)/BOLD |
| 145 |
IF(BDELTA.GT.BDEL) THEN |
| 146 |
STEP=STEP/10. |
| 147 |
GOTO 7777 |
| 148 |
ENDIF |
| 149 |
8888 RR0=ROLD |
| 150 |
BEQU=BOLD |
| 151 |
BDEL=BDELTA |
| 152 |
RETURN |
| 153 |
END |
| 154 |
C |
| 155 |
C |
| 156 |
SUBROUTINE SHELLG(GLAT,GLON,ALT,DIMO,FL,ICODE,B0) |
| 157 |
C-------------------------------------------------------------------- |
| 158 |
C CALCULATES L-VALUE FOR SPECIFIED GEODAETIC COORDINATES, ALTITUDE |
| 159 |
C AND GEMAGNETIC FIELD MODEL. |
| 160 |
C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTER, INTERNAL NOTE |
| 161 |
C NO. 67, 1970. |
| 162 |
C G. KLUGE, COMPUTER PHYSICS COMMUNICATIONS 3, 31-35, 1972 |
| 163 |
C-------------------------------------------------------------------- |
| 164 |
C CHANGES (D. BILITZA, NOV 87): |
| 165 |
C - USING CORRECT DIPOL MOMENT I.E.,DIFFERENT COMMON/MODEL/ |
| 166 |
C - USING IGRF EARTH MAGNETIC FIELD MODELS FROM 1945 TO 1990 |
| 167 |
C-------------------------------------------------------------------- |
| 168 |
C INPUT: ENTRY POINT SHELLG |
| 169 |
C GLAT GEODETIC LATITUDE IN DEGREES (NORTH) |
| 170 |
C GLON GEODETIC LONGITUDE IN DEGREES (EAST) |
| 171 |
C ALT ALTITUDE IN KM ABOVE SEA LEVEL |
| 172 |
C |
| 173 |
C ENTRY POINT SHELLC |
| 174 |
C V(3) CARTESIAN COORDINATES IN EARTH RADII (6371.2 KM) |
| 175 |
C X-AXIS POINTING TO EQUATOR AT 0 LONGITUDE |
| 176 |
C Y-AXIS POINTING TO EQUATOR AT 90 LONG. |
| 177 |
C Z-AXIS POINTING TO NORTH POLE |
| 178 |
C |
| 179 |
C DIMO DIPOL MOMENT IN GAUSS (NORMALIZED TO EARTH RADIUS) |
| 180 |
C |
| 181 |
C COMMON |
| 182 |
C X(3) NOT USED |
| 183 |
C H(144) FIELD MODEL COEFFICIENTS ADJUSTED FOR SHELLG |
| 184 |
C----------------------------------------------------------------------- |
| 185 |
C OUTPUT: FL L-VALUE |
| 186 |
C ICODE =1 NORMAL COMPLETION |
| 187 |
C =2 UNPHYSICAL CONJUGATE POINT (FL MEANINGLESS) |
| 188 |
C =3 SHELL PARAMETER GREATER THAN LIMIT UP TO |
| 189 |
C WHICH ACCURATE CALCULATION IS REQUIRED; |
| 190 |
C APPROXIMATION IS USED. |
| 191 |
C B0 MAGNETIC FIELD STRENGTH IN GAUSS |
| 192 |
C----------------------------------------------------------------------- |
| 193 |
DIMENSION V(3),U(3,3),P(8,100),SP(3) |
| 194 |
COMMON X(3),H(196) |
| 195 |
COMMON/FIDB0/ SP |
| 196 |
SAVE /FIDB0/ |
| 197 |
COMMON/GENER/ UMR,ERA,AQUAD,BQUAD |
| 198 |
SAVE /GENER/ |
| 199 |
REAL FLS |
| 200 |
C |
| 201 |
C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3 |
| 202 |
C-- STEP IS STEP SIZE FOR FIELD LINE TRACING |
| 203 |
C-- STEQ IS STEP SIZE FOR INTEGRATION |
| 204 |
C |
| 205 |
DATA RMIN,RMAX /0.05,1.01/ |
| 206 |
DATA STEP,STEQ /0.20,0.03/ |
| 207 |
BEQU=1.E10 |
| 208 |
FLS = FL |
| 209 |
C*****ENTRY POINT SHELLG TO BE USED WITH GEODETIC CO-ORDINATES |
| 210 |
RLAT=GLAT*UMR |
| 211 |
CT=SIN(RLAT) |
| 212 |
ST=COS(RLAT) |
| 213 |
D=SQRT(AQUAD-(AQUAD-BQUAD)*CT*CT) |
| 214 |
X(1)=(ALT+AQUAD/D)*ST/ERA |
| 215 |
X(3)=(ALT+BQUAD/D)*CT/ERA |
| 216 |
RLON=GLON*UMR |
| 217 |
X(2)=X(1)*SIN(RLON) |
| 218 |
X(1)=X(1)*COS(RLON) |
| 219 |
GOTO9 |
| 220 |
ENTRY SHELLC(V,FL,B0) |
| 221 |
C*****ENTRY POINT SHELLC TO BE USED WITH CARTESIAN CO-ORDINATES |
| 222 |
X(1)=V(1) |
| 223 |
X(2)=V(2) |
| 224 |
X(3)=V(3) |
| 225 |
C*****CONVERT TO DIPOL-ORIENTED CO-ORDINATES |
| 226 |
DATA U/ +0.3511737,-0.9148385,-0.1993679, |
| 227 |
A +0.9335804,+0.3583680,+0.0000000, |
| 228 |
B +0.0714471,-0.1861260,+0.9799247/ |
| 229 |
9 RQ=1./(X(1)*X(1)+X(2)*X(2)+X(3)*X(3)) |
| 230 |
R3H=SQRT(RQ*SQRT(RQ)) |
| 231 |
P(1,2)=(X(1)*U(1,1)+X(2)*U(2,1)+X(3)*U(3,1))*R3H |
| 232 |
P(2,2)=(X(1)*U(1,2)+X(2)*U(2,2) )*R3H |
| 233 |
P(3,2)=(X(1)*U(1,3)+X(2)*U(2,3)+X(3)*U(3,3))*RQ |
| 234 |
C*****FIRST THREE POINTS OF FIELD LINE |
| 235 |
STEP=-SIGN(STEP,P(3,2)) |
| 236 |
CALL STOER(P(1,2),BQ2,R2) |
| 237 |
B0=SQRT(BQ2) |
| 238 |
P(1,3)=P(1,2)+0.5*STEP*P(4,2) |
| 239 |
P(2,3)=P(2,2)+0.5*STEP*P(5,2) |
| 240 |
P(3,3)=P(3,2)+0.5*STEP |
| 241 |
CALL STOER(P(1,3),BQ3,R3) |
| 242 |
P(1,1)=P(1,2)-STEP*(2.*P(4,2)-P(4,3)) |
| 243 |
P(2,1)=P(2,2)-STEP*(2.*P(5,2)-P(5,3)) |
| 244 |
P(3,1)=P(3,2)-STEP |
| 245 |
CALL STOER(P(1,1),BQ1,R1) |
| 246 |
P(1,3)=P(1,2)+STEP*(20.*P(4,3)-3.*P(4,2)+P(4,1))/18. |
| 247 |
P(2,3)=P(2,2)+STEP*(20.*P(5,3)-3.*P(5,2)+P(5,1))/18. |
| 248 |
P(3,3)=P(3,2)+STEP |
| 249 |
CALL STOER(P(1,3),BQ3,R3) |
| 250 |
C*****INVERT SENSE IF REQUIRED |
| 251 |
IF(BQ3.LE.BQ1)GOTO2 |
| 252 |
STEP=-STEP |
| 253 |
R3=R1 |
| 254 |
BQ3=BQ1 |
| 255 |
DO 1 I=1,7 |
| 256 |
ZZ=P(I,1) |
| 257 |
P(I,1)=P(I,3) |
| 258 |
1 P(I,3)=ZZ |
| 259 |
C*****SEARCH FOR LOWEST MAGNETIC FIELD STRENGTH |
| 260 |
2 IF(BQ1.LT.BEQU) THEN |
| 261 |
BEQU=BQ1 |
| 262 |
IEQU=1 |
| 263 |
ENDIF |
| 264 |
IF(BQ2.LT.BEQU) THEN |
| 265 |
BEQU=BQ2 |
| 266 |
IEQU=2 |
| 267 |
ENDIF |
| 268 |
IF(BQ3.LT.BEQU) THEN |
| 269 |
BEQU=BQ3 |
| 270 |
IEQU=3 |
| 271 |
ENDIF |
| 272 |
C*****INITIALIZATION OF INTEGRATION LOOPS |
| 273 |
STEP12=STEP/12. |
| 274 |
STEP2=STEP+STEP |
| 275 |
STEQ=SIGN(STEQ,STEP) |
| 276 |
FI=0. |
| 277 |
ICODE=1 |
| 278 |
ORADIK=0. |
| 279 |
OTERM=0. |
| 280 |
STP=R2*STEQ |
| 281 |
Z=P(3,2)+STP |
| 282 |
STP=STP/0.75 |
| 283 |
P(8,1)=STEP2*(P(1,1)*P(4,1)+P(2,1)*P(5,1)) |
| 284 |
P(8,2)=STEP2*(P(1,2)*P(4,2)+P(2,2)*P(5,2)) |
| 285 |
C*****MAIN LOOP (FIELD LINE TRACING) |
| 286 |
DO 3 N=3,3333 |
| 287 |
C*****CORRECTOR (FIELD LINE TRACING) |
| 288 |
P(1,N)=P(1,N-1)+STEP12*(5.*P(4,N)+8.*P(4,N-1)-P(4,N-2)) |
| 289 |
P(2,N)=P(2,N-1)+STEP12*(5.*P(5,N)+8.*P(5,N-1)-P(5,N-2)) |
| 290 |
C*****PREPARE EXPANSION COEFFICIENTS FOR INTERPOLATION |
| 291 |
C*****OF SLOWLY VARYING QUANTITIES |
| 292 |
P(8,N)=STEP2*(P(1,N)*P(4,N)+P(2,N)*P(5,N)) |
| 293 |
C0=P(1,N-1)**2+P(2,N-1)**2 |
| 294 |
C1=P(8,N-1) |
| 295 |
C2=(P(8,N)-P(8,N-2))*0.25 |
| 296 |
C3=(P(8,N)+P(8,N-2)-C1-C1)/6.0 |
| 297 |
D0=P(6,N-1) |
| 298 |
D1=(P(6,N)-P(6,N-2))*0.5 |
| 299 |
D2=(P(6,N)+P(6,N-2)-D0-D0)*0.5 |
| 300 |
E0=P(7,N-1) |
| 301 |
E1=(P(7,N)-P(7,N-2))*0.5 |
| 302 |
E2=(P(7,N)+P(7,N-2)-E0-E0)*0.5 |
| 303 |
C*****INNER LOOP (FOR QUADRATURE) |
| 304 |
4 T=(Z-P(3,N-1))/STEP |
| 305 |
IF(T.GT.1.)GOTO5 |
| 306 |
HLI=0.5*(((C3*T+C2)*T+C1)*T+C0) |
| 307 |
ZQ=Z*Z |
| 308 |
R=HLI+SQRT(HLI*HLI+ZQ) |
| 309 |
IF(R.NE.R)THEN |
| 310 |
FL = FLS |
| 311 |
RETURN |
| 312 |
ENDIF |
| 313 |
IF(R.LE.RMIN)GOTO30 |
| 314 |
RQ=R*R |
| 315 |
FF=SQRT(1.+3.*ZQ/RQ) |
| 316 |
RADIK=B0-((D2*T+D1)*T+D0)*R*RQ*FF |
| 317 |
IF((R-RMAX).le.0.) goto 44 |
| 318 |
IF((R-RMAX).gt.0.) goto 45 |
| 319 |
45 ICODE=2 |
| 320 |
RADIK=RADIK-12.*(R-RMAX)**2 |
| 321 |
44 IF(RADIK+RADIK.LE.ORADIK) GOTO 10 |
| 322 |
TERM=SQRT(RADIK)*FF*((E2*T+E1)*T+E0)/(RQ+ZQ) |
| 323 |
FI=FI+STP*(OTERM+TERM) |
| 324 |
ORADIK=RADIK |
| 325 |
OTERM=TERM |
| 326 |
STP=R*STEQ |
| 327 |
Z=Z+STP |
| 328 |
GOTO4 |
| 329 |
C*****PREDICTOR (FIELD LINE TRACING) |
| 330 |
5 P(1,N+1)=P(1,N)+STEP12*(23.*P(4,N)-16.*P(4,N-1)+5.*P(4,N-2)) |
| 331 |
P(2,N+1)=P(2,N)+STEP12*(23.*P(5,N)-16.*P(5,N-1)+5.*P(5,N-2)) |
| 332 |
P(3,N+1)=P(3,N)+STEP |
| 333 |
CALL STOER(P(1,N+1),BQ3,R3) |
| 334 |
C*****SEARCH FOR LOWEST MAGNETIC FIELD STRENGTH |
| 335 |
IF(BQ3.LT.BEQU) THEN |
| 336 |
IEQU=N+1 |
| 337 |
BEQU=BQ3 |
| 338 |
ENDIF |
| 339 |
3 CONTINUE |
| 340 |
10 IF(IEQU.lt.2) IEQU=2 |
| 341 |
SP(1)=P(1,IEQU-1) |
| 342 |
SP(2)=P(2,IEQU-1) |
| 343 |
SP(3)=P(3,IEQU-1) |
| 344 |
IF(ORADIK.LT.1E-15)GOTO11 |
| 345 |
FI=FI+STP/0.75*OTERM*ORADIK/(ORADIK-RADIK) |
| 346 |
C |
| 347 |
C-- The minimal allowable value of FI was changed from 1E-15 to 1E-12, |
| 348 |
C-- because 1E-38 is the minimal allowable arg. for ALOG in our envir. |
| 349 |
C-- D. Bilitza, Nov 87. |
| 350 |
C |
| 351 |
11 FI=0.5*ABS(FI)/SQRT(B0)+1E-12 |
| 352 |
C |
| 353 |
C*****COMPUTE L FROM B AND I. SAME AS CARMEL IN INVAR. |
| 354 |
C |
| 355 |
C-- Correct dipole moment is used here. D. Bilitza, Nov 87. |
| 356 |
C |
| 357 |
DIMOB0=DIMO/B0 |
| 358 |
arg1=alog(FI) |
| 359 |
arg2=alog(DIMOB0) |
| 360 |
c arg = FI*FI*FI/DIMOB0 |
| 361 |
c if(abs(arg).gt.88.0) arg=88.0 |
| 362 |
XX=3*arg1-arg2 |
| 363 |
IF(XX.GT.23.0) GOTO 776 |
| 364 |
IF(XX.GT.11.7) GOTO 775 |
| 365 |
IF(XX.GT.+3.0) GOTO 774 |
| 366 |
IF(XX.GT.-3.0) GOTO 773 |
| 367 |
IF(XX.GT.-22.) GOTO 772 |
| 368 |
c 771 GG=3.33338E-1*XX+3.0062102E-1 |
| 369 |
GG=3.33338E-1*XX+3.0062102E-1 |
| 370 |
GOTO777 |
| 371 |
772 GG=((((((((-8.1537735E-14*XX+8.3232531E-13)*XX+1.0066362E-9)*XX+ |
| 372 |
18.1048663E-8)*XX+3.2916354E-6)*XX+8.2711096E-5)*XX+1.3714667E-3)* |
| 373 |
2XX+1.5017245E-2)*XX+4.3432642E-1)*XX+6.2337691E-1 |
| 374 |
GOTO777 |
| 375 |
773 GG=((((((((2.6047023E-10*XX+2.3028767E-9)*XX-2.1997983E-8)*XX- |
| 376 |
15.3977642E-7)*XX-3.3408822E-6)*XX+3.8379917E-5)*XX+1.1784234E-3)* |
| 377 |
2XX+1.4492441E-2)*XX+4.3352788E-1)*XX+6.228644E-1 |
| 378 |
GOTO777 |
| 379 |
774 GG=((((((((6.3271665E-10*XX-3.958306E-8)*XX+9.9766148E-07)*XX- |
| 380 |
11.2531932E-5)*XX+7.9451313E-5)*XX-3.2077032E-4)*XX+2.1680398E-3)* |
| 381 |
2XX+1.2817956E-2)*XX+4.3510529E-1)*XX+6.222355E-1 |
| 382 |
GOTO777 |
| 383 |
775 GG=(((((2.8212095E-8*XX-3.8049276E-6)*XX+2.170224E-4)*XX-6.7310339 |
| 384 |
1E-3)*XX+1.2038224E-1)*XX-1.8461796E-1)*XX+2.0007187E0 |
| 385 |
GOTO777 |
| 386 |
776 GG=XX-3.0460681E0 |
| 387 |
777 FL=EXP(ALOG((1.+EXP(GG))*DIMOB0)/3.0) |
| 388 |
RETURN |
| 389 |
C*****APPROXIMATION FOR HIGH VALUES OF L. |
| 390 |
30 ICODE=3 |
| 391 |
T=-P(3,N-1)/STEP |
| 392 |
FL=1./(ABS(((C3*T+C2)*T+C1)*T+C0)+1E-15) |
| 393 |
RETURN |
| 394 |
END |
| 395 |
C |
| 396 |
C |
| 397 |
SUBROUTINE STOER(P,BQ,R) |
| 398 |
C******************************************************************* |
| 399 |
C* SUBROUTINE USED FOR FIELD LINE TRACING IN SHELLG * |
| 400 |
C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG * |
| 401 |
C******************************************************************* |
| 402 |
DIMENSION P(7),U(3,3) |
| 403 |
COMMON XI(3),H(196) |
| 404 |
C*****XM,YM,ZM ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES |
| 405 |
ZM=P(3) |
| 406 |
FLI=P(1)*P(1)+P(2)*P(2)+1E-15 |
| 407 |
R=0.5*(FLI+SQRT(FLI*FLI+(ZM+ZM)**2)) |
| 408 |
RQ=R*R |
| 409 |
WR=SQRT(R) |
| 410 |
XM=P(1)*WR |
| 411 |
YM=P(2)*WR |
| 412 |
C*****TRANSFORM TO GEOGRAPHIC CO-ORDINATE SYSTEM |
| 413 |
DATA U/ +0.3511737,-0.9148385,-0.1993679, |
| 414 |
A +0.9335804,+0.3583680,+0.0000000, |
| 415 |
B +0.0714471,-0.1861260,+0.9799247/ |
| 416 |
XI(1)=XM*U(1,1)+YM*U(1,2)+ZM*U(1,3) |
| 417 |
XI(2)=XM*U(2,1)+YM*U(2,2)+ZM*U(2,3) |
| 418 |
XI(3)=XM*U(3,1) +ZM*U(3,3) |
| 419 |
C*****COMPUTE DERIVATIVES |
| 420 |
C CALL FELDI(XI,H) |
| 421 |
CALL FELDI |
| 422 |
Q=H(1)/RQ |
| 423 |
DX=H(3)+H(3)+Q*XI(1) |
| 424 |
DY=H(4)+H(4)+Q*XI(2) |
| 425 |
DZ=H(2)+H(2)+Q*XI(3) |
| 426 |
C*****TRANSFORM BACK TO GEOMAGNETIC CO-ORDINATE SYSTEM |
| 427 |
DXM=U(1,1)*DX+U(2,1)*DY+U(3,1)*DZ |
| 428 |
DYM=U(1,2)*DX+U(2,2)*DY |
| 429 |
DZM=U(1,3)*DX+U(2,3)*DY+U(3,3)*DZ |
| 430 |
DR=(XM*DXM+YM*DYM+ZM*DZM)/R |
| 431 |
C*****FORM SLOWLY VARYING EXPRESSIONS |
| 432 |
P(4)=(WR*DXM-0.5*P(1)*DR)/(R*DZM) |
| 433 |
P(5)=(WR*DYM-0.5*P(2)*DR)/(R*DZM) |
| 434 |
DSQ=RQ*(DXM*DXM+DYM*DYM+DZM*DZM) |
| 435 |
BQ=DSQ*RQ*RQ |
| 436 |
P(6)=SQRT(DSQ/(RQ+3.*ZM*ZM)) |
| 437 |
P(7)=P(6)*(RQ+ZM*ZM)/(RQ*DZM) |
| 438 |
RETURN |
| 439 |
END |
| 440 |
C |
| 441 |
C |
| 442 |
SUBROUTINE FELDG(GLAT,GLON,ALT,BNORTH,BEAST,BDOWN,BABS) |
| 443 |
C------------------------------------------------------------------- |
| 444 |
C CALCULATES EARTH MAGNETIC FIELD FROM SPHERICAL HARMONICS MODEL |
| 445 |
C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTRE, INTERNAL NOTE 61, |
| 446 |
C 1970. |
| 447 |
C-------------------------------------------------------------------- |
| 448 |
C CHANGES (D. BILITZA, NOV 87): |
| 449 |
C - FIELD COEFFICIENTS IN BINARY DATA FILES INSTEAD OF BLOCK DATA |
| 450 |
C - CALCULATES DIPOL MOMENT |
| 451 |
C-------------------------------------------------------------------- |
| 452 |
C INPUT: ENTRY POINT FELDG |
| 453 |
C GLAT GEODETIC LATITUDE IN DEGREES (NORTH) |
| 454 |
C GLON GEODETIC LONGITUDE IN DEGREES (EAST) |
| 455 |
C ALT ALTITUDE IN KM ABOVE SEA LEVEL |
| 456 |
C |
| 457 |
C ENTRY POINT FELDC |
| 458 |
C V(3) CARTESIAN COORDINATES IN EARTH RADII (6371.2 KM) |
| 459 |
C X-AXIS POINTING TO EQUATOR AT 0 LONGITUDE |
| 460 |
C Y-AXIS POINTING TO EQUATOR AT 90 LONG. |
| 461 |
C Z-AXIS POINTING TO NORTH POLE |
| 462 |
C |
| 463 |
C COMMON BLANK AND ENTRY POINT FELDI ARE NEEDED WHEN USED |
| 464 |
C IN CONNECTION WITH L-CALCULATION PROGRAM SHELLG. |
| 465 |
C |
| 466 |
C COMMON /MODEL/ AND /GENER/ |
| 467 |
C UMR = ATAN(1.0)*4./180. <DEGREE>*UMR=<RADIANT> |
| 468 |
C ERA EARTH RADIUS FOR NORMALIZATION OF CARTESIAN |
| 469 |
C COORDINATES (6371.2 KM) |
| 470 |
C AQUAD, BQUAD SQUARE OF MAJOR AND MINOR HALF AXIS FOR |
| 471 |
C EARTH ELLIPSOID AS RECOMMENDED BY INTERNATIONAL |
| 472 |
C ASTRONOMICAL UNION (6378.160, 6356.775 KM). |
| 473 |
C NMAX MAXIMUM ORDER OF SPHERICAL HARMONICS |
| 474 |
C TIME YEAR (DECIMAL: 1973.5) FOR WHICH MAGNETIC |
| 475 |
C FIELD IS TO BE CALCULATED |
| 476 |
C G(M) NORMALIZED FIELD COEFFICIENTS (SEE FELDCOF) |
| 477 |
C M=NMAX*(NMAX+2) |
| 478 |
C------------------------------------------------------------------------ |
| 479 |
C OUTPUT: BABS MAGNETIC FIELD STRENGTH IN GAUSS |
| 480 |
C BNORTH, BEAST, BDOWN COMPONENTS OF THE FIELD WITH RESPECT |
| 481 |
C TO THE LOCAL GEODETIC COORDINATE SYSTEM, WITH AXIS |
| 482 |
C POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST |
| 483 |
C AND DOWNWARD. |
| 484 |
C----------------------------------------------------------------------- |
| 485 |
DIMENSION V(3),B(3),G(196) |
| 486 |
CHARACTER*258 NAME |
| 487 |
INTEGER NMAX |
| 488 |
REAL TIME |
| 489 |
COMMON XI(3),H(196) |
| 490 |
COMMON/MODEL/ G,NMAX,TIME,NAME |
| 491 |
SAVE/MODEL/ |
| 492 |
COMMON/GENER/ UMR,ERA,AQUAD,BQUAD |
| 493 |
SAVE/GENER/ |
| 494 |
C |
| 495 |
C-- IS RECORDS ENTRY POINT |
| 496 |
C |
| 497 |
C*****ENTRY POINT FELDG TO BE USED WITH GEODETIC CO-ORDINATES |
| 498 |
IS=1 |
| 499 |
RLAT=GLAT*UMR |
| 500 |
CT=SIN(RLAT) |
| 501 |
ST=COS(RLAT) |
| 502 |
D=SQRT(AQUAD-(AQUAD-BQUAD)*CT*CT) |
| 503 |
RLON=GLON*UMR |
| 504 |
CP=COS(RLON) |
| 505 |
SP=SIN(RLON) |
| 506 |
ZZZ=(ALT+BQUAD/D)*CT/ERA |
| 507 |
RHO=(ALT+AQUAD/D)*ST/ERA |
| 508 |
XXX=RHO*CP |
| 509 |
YYY=RHO*SP |
| 510 |
GOTO10 |
| 511 |
ENTRY FELDC(V,B) |
| 512 |
C*****ENTRY POINT FELDC TO BE USED WITH CARTESIAN CO-ORDINATES |
| 513 |
IS=2 |
| 514 |
XXX=V(1) |
| 515 |
YYY=V(2) |
| 516 |
ZZZ=V(3) |
| 517 |
10 RQ=1./(XXX*XXX+YYY*YYY+ZZZ*ZZZ) |
| 518 |
XI(1)=XXX*RQ |
| 519 |
XI(2)=YYY*RQ |
| 520 |
XI(3)=ZZZ*RQ |
| 521 |
GOTO20 |
| 522 |
ENTRY FELDI |
| 523 |
C*****ENTRY POINT FELDI USED FOR L COMPUTATION |
| 524 |
IS=3 |
| 525 |
20 IHMAX=NMAX*NMAX+1 |
| 526 |
LAST=IHMAX+NMAX+NMAX |
| 527 |
IMAX=NMAX+NMAX-1 |
| 528 |
DO 8 I=IHMAX,LAST |
| 529 |
8 H(I)=G(I) |
| 530 |
DO 6 K=1,3,2 |
| 531 |
I=IMAX |
| 532 |
IH=IHMAX |
| 533 |
1 IL=IH-I |
| 534 |
F=2./FLOAT(I-K+2) |
| 535 |
X=XI(1)*F |
| 536 |
Y=XI(2)*F |
| 537 |
Z=XI(3)*(F+F) |
| 538 |
I=I-2 |
| 539 |
c print *,' I ',I |
| 540 |
IF((I-1).lt.0) goto 5 |
| 541 |
IF((I-1).eq.0) goto 4 |
| 542 |
IF((I-1).gt.0) goto 2 |
| 543 |
2 DO 3 M=3,I,2 |
| 544 |
H(IL+M+1)=G(IL+M+1)+Z*H(IH+M+1)+X*(H(IH+M+3)-H(IH+M-1)) |
| 545 |
A -Y*(H(IH+M+2)+H(IH+M-2)) |
| 546 |
3 H(IL+M)=G(IL+M)+Z*H(IH+M)+X*(H(IH+M+2)-H(IH+M-2)) |
| 547 |
A +Y*(H(IH+M+3)+H(IH+M-1)) |
| 548 |
4 H(IL+2)=G(IL+2)+Z*H(IH+2)+X*H(IH+4)-Y*(H(IH+3)+H(IH)) |
| 549 |
H(IL+1)=G(IL+1)+Z*H(IH+1)+Y*H(IH+4)+X*(H(IH+3)-H(IH)) |
| 550 |
5 H(IL)=G(IL)+Z*H(IH)+2.*(X*H(IH+1)+Y*H(IH+2)) |
| 551 |
IH=IL |
| 552 |
IF(I.GE.K)GOTO1 |
| 553 |
6 CONTINUE |
| 554 |
IF(IS.EQ.3)RETURN |
| 555 |
S=.5*H(1)+2.*(H(2)*XI(3)+H(3)*XI(1)+H(4)*XI(2)) |
| 556 |
T=(RQ+RQ)*SQRT(RQ) |
| 557 |
BXXX=T*(H(3)-S*XXX) |
| 558 |
BYYY=T*(H(4)-S*YYY) |
| 559 |
BZZZ=T*(H(2)-S*ZZZ) |
| 560 |
IF(IS.EQ.2)GOTO7 |
| 561 |
BABS=SQRT(BXXX*BXXX+BYYY*BYYY+BZZZ*BZZZ) |
| 562 |
BEAST=BYYY*CP-BXXX*SP |
| 563 |
BRHO=BYYY*SP+BXXX*CP |
| 564 |
BNORTH=BZZZ*ST-BRHO*CT |
| 565 |
BDOWN=-BZZZ*CT-BRHO*ST |
| 566 |
RETURN |
| 567 |
7 B(1)=BXXX |
| 568 |
B(2)=BYYY |
| 569 |
B(3)=BZZZ |
| 570 |
RETURN |
| 571 |
END |
| 572 |
C |
| 573 |
C |
| 574 |
SUBROUTINE FELDCOF(YEAR,DIMO) |
| 575 |
C------------------------------------------------------------------------ |
| 576 |
C DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS |
| 577 |
C |
| 578 |
C INPUT: YEAR DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO |
| 579 |
C BE CALCULATED |
| 580 |
C OUTPUT: DIMO GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED |
| 581 |
C TO EARTH'S RADIUS) AT THE TIME (YEAR) |
| 582 |
C D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771, |
| 583 |
C (301)286-9536 NOV 1987. |
| 584 |
C ### updated to IGRF-2000 version -dkb- 5/31/2000 |
| 585 |
C ### updated to IGRF-2005 version -dkb- 3/24/2000 |
| 586 |
C----------------------------------------------------------------------- |
| 587 |
CHARACTER*258 FIL1, FIL2 |
| 588 |
CHARACTER*258 FILMOD |
| 589 |
C ### FILMOD, DTEMOD arrays +1 |
| 590 |
c DIMENSION GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14) |
| 591 |
DIMENSION GH1(196),GH2(196),GHA(196),FILMOD(3),DTEMOD(3) |
| 592 |
DOUBLE PRECISION X,F0,F |
| 593 |
INTEGER L1,L2,L3 |
| 594 |
INTEGER NMAX |
| 595 |
REAL TIME |
| 596 |
CHARACTER *258 P1,P2,P3 |
| 597 |
COMMON/PPATH/ L1,L2,L3,P1, P2, P3 |
| 598 |
SAVE/PPATH/ |
| 599 |
COMMON/MODEL/ GH1,NMAX,TIME,FIL1 |
| 600 |
SAVE/MODEL/ |
| 601 |
COMMON/GENER/ UMR,ERAD,AQUAD,BQUAD |
| 602 |
SAVE/GENER/ |
| 603 |
C ### updated to 2005 |
| 604 |
C CHARACTER COEFPATH*80, COEF1*80, COEF2*80, COEF3*80 |
| 605 |
|
| 606 |
c COEFPATH = 'OrbitalInfo/src/' |
| 607 |
c COEF1 = 'dgrf00.dat' |
| 608 |
c COEF2 = 'igrf05.dat' |
| 609 |
c COEF3 = 'igrf05s.dat' |
| 610 |
c COEF1 = COEFPATH(1:16)//COEF1 |
| 611 |
c COEF2 = COEFPATH(1:16)//COEF2 |
| 612 |
c COEF3 = COEFPATH(1:16)//COEF3 |
| 613 |
c FILMOD(1) = COEF1 |
| 614 |
c FILMOD(2) = COEF2 |
| 615 |
c FILMOD(3) = COEF3 |
| 616 |
c print *, "qui" |
| 617 |
FILMOD(1) = P1(1:L1) |
| 618 |
FILMOD(2) = P2(1:L2) |
| 619 |
FILMOD(3) = P3(1:L3) |
| 620 |
c print *, "qua" |
| 621 |
c FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat' |
| 622 |
c FILMOD(2) = 'OrbitalInfo/src/igrf05.dat' |
| 623 |
c FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat' |
| 624 |
c WRITE(*,*) FILMOD(1) |
| 625 |
c WRITE(*,*) FILMOD(2) |
| 626 |
c WRITE(*,*) FILMOD(3) |
| 627 |
c DATA FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/ |
| 628 |
DATA DTEMOD / 2005., 2010., 2015./ |
| 629 |
c |
| 630 |
c DATA FILMOD /'dgrf45.dat', 'dgrf50.dat', |
| 631 |
c 1 'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat', |
| 632 |
c 2 'dgrf70.dat', 'dgrf75.dat', 'dgrf80.dat', |
| 633 |
c 3 'dgrf85.dat', 'dgrf90.dat', 'dgrf95.dat', |
| 634 |
c 4 'dgrf00.dat','igrf05.dat','igrf05s.dat'/ |
| 635 |
c DATA DTEMOD / 1945., 1950., 1955., 1960., 1965., 1970., |
| 636 |
c 1 1975., 1980., 1985., 1990., 1995., 2000.,2005.,2010./ |
| 637 |
C |
| 638 |
C ### numye = numye + 1 ; is number of years represented by IGRF |
| 639 |
C |
| 640 |
c NUMYE=13 |
| 641 |
NUMYE=2 |
| 642 |
c print *, "quo" |
| 643 |
|
| 644 |
C |
| 645 |
C IS=0 FOR SCHMIDT NORMALIZATION IS=1 GAUSS NORMALIZATION |
| 646 |
C IU IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS |
| 647 |
C |
| 648 |
IU = 10 |
| 649 |
IS = 0 |
| 650 |
C-- DETERMINE IGRF-YEARS FOR INPUT-YEAR |
| 651 |
TIME = YEAR |
| 652 |
IYEA = INT(YEAR/5.)*5 |
| 653 |
c L = (IYEA - 1945)/5 + 1 |
| 654 |
L = (IYEA - 2000)/5 + 1 |
| 655 |
IF(L.LT.1) L=1 |
| 656 |
IF(L.GT.NUMYE) L=NUMYE |
| 657 |
DTE1 = DTEMOD(L) |
| 658 |
FIL1 = FILMOD(L) |
| 659 |
DTE2 = DTEMOD(L+1) |
| 660 |
FIL2 = FILMOD(L+1) |
| 661 |
c WRITE(*,*) FIL1 |
| 662 |
c WRITE(*,*) FIL2 |
| 663 |
c print *, "que" |
| 664 |
C-- GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS |
| 665 |
CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER) |
| 666 |
IF (IER .NE. 0) STOP |
| 667 |
c print *, "quessss" |
| 668 |
CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER) |
| 669 |
IF (IER .NE. 0) STOP |
| 670 |
c print *, "quj" |
| 671 |
C-- DETERMINE IGRF COEFFICIENTS FOR YEAR |
| 672 |
IF (L .LE. NUMYE-1) THEN |
| 673 |
CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2, |
| 674 |
1 NMAX2, GH2, NMAX, GHA) |
| 675 |
ELSE |
| 676 |
CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2, |
| 677 |
1 GH2, NMAX, GHA) |
| 678 |
ENDIF |
| 679 |
c print *, "quw" |
| 680 |
C-- DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G |
| 681 |
F0=0.D0 |
| 682 |
DO 1234 J=1,3 |
| 683 |
F = GHA(J) * 1.D-5 |
| 684 |
F0 = F0 + F * F |
| 685 |
1234 CONTINUE |
| 686 |
DIMO = REAL(DSQRT(F0)) |
| 687 |
|
| 688 |
GH1(1) = 0.0 |
| 689 |
I=2 |
| 690 |
F0=1.D-5 |
| 691 |
IF(IS.EQ.0) F0=-F0 |
| 692 |
SQRT2=SQRT(2.) |
| 693 |
|
| 694 |
c print *, "quq" |
| 695 |
|
| 696 |
DO 9 N=1,NMAX |
| 697 |
X = N |
| 698 |
F0 = F0 * X * X / (4.D0 * X - 2.D0) |
| 699 |
IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X |
| 700 |
F = F0 * 0.5D0 |
| 701 |
IF(IS.EQ.0) F = F * SQRT2 |
| 702 |
GH1(I) = GHA(I-1) * REAL(F0) |
| 703 |
I = I+1 |
| 704 |
DO 9 M=1,N |
| 705 |
F = F * (X + M) / (X - M + 1.D0) |
| 706 |
IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M)) |
| 707 |
GH1(I) = GHA(I-1) * REAL(F) |
| 708 |
GH1(I+1) = GHA(I) * REAL(F) |
| 709 |
I=I+2 |
| 710 |
9 CONTINUE |
| 711 |
RETURN |
| 712 |
END |
| 713 |
C |
| 714 |
C |
| 715 |
SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER) |
| 716 |
|
| 717 |
C =============================================================== |
| 718 |
C |
| 719 |
C Version 1.01 |
| 720 |
C |
| 721 |
C Reads spherical harmonic coefficients from the specified |
| 722 |
C file into an array. |
| 723 |
C |
| 724 |
C Input: |
| 725 |
C IU - Logical unit number |
| 726 |
C FSPEC - File specification |
| 727 |
C |
| 728 |
C Output: |
| 729 |
C NMAX - Maximum degree and order of model |
| 730 |
C ERAD - Earth's radius associated with the spherical |
| 731 |
C harmonic coefficients, in the same units as |
| 732 |
C elevation |
| 733 |
C GH - Schmidt quasi-normal internal spherical |
| 734 |
C harmonic coefficients |
| 735 |
C IER - Error number: = 0, no error |
| 736 |
C = -2, records out of order |
| 737 |
C = FORTRAN run-time error number |
| 738 |
C |
| 739 |
C A. Zunde |
| 740 |
C USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 |
| 741 |
C |
| 742 |
C =============================================================== |
| 743 |
|
| 744 |
CHARACTER FSPEC*(*), FOUT*258 |
| 745 |
DIMENSION GH(*) |
| 746 |
C --------------------------------------------------------------- |
| 747 |
C Open coefficient file. Read past first header record. |
| 748 |
C Read degree and order of model and Earth's radius. |
| 749 |
C --------------------------------------------------------------- |
| 750 |
WRITE(FOUT,667) FSPEC |
| 751 |
c 667 FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12) |
| 752 |
667 FORMAT(A258) |
| 753 |
c print *," gui" |
| 754 |
OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999) |
| 755 |
c print *," gua" |
| 756 |
READ (IU, *, IOSTAT=IER, ERR=999) |
| 757 |
c print *," gue" |
| 758 |
READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD |
| 759 |
c print *," guo" |
| 760 |
C --------------------------------------------------------------- |
| 761 |
C Read the coefficient file, arranged as follows: |
| 762 |
C |
| 763 |
C N M G H |
| 764 |
C ---------------------- |
| 765 |
C / 1 0 GH(1) - |
| 766 |
C / 1 1 GH(2) GH(3) |
| 767 |
C / 2 0 GH(4) - |
| 768 |
C / 2 1 GH(5) GH(6) |
| 769 |
C NMAX*(NMAX+3)/2 / 2 2 GH(7) GH(8) |
| 770 |
C records \ 3 0 GH(9) - |
| 771 |
C \ . . . . |
| 772 |
C \ . . . . |
| 773 |
C NMAX*(NMAX+2) \ . . . . |
| 774 |
C elements in GH \ NMAX NMAX . . |
| 775 |
C |
| 776 |
C N and M are, respectively, the degree and order of the |
| 777 |
C coefficient. |
| 778 |
C --------------------------------------------------------------- |
| 779 |
|
| 780 |
I = 0 |
| 781 |
DO 2211 NN = 1, NMAX |
| 782 |
DO 2233 MM = 0, NN |
| 783 |
READ (IU, *, IOSTAT=IER, ERR=999) N, M, G, H |
| 784 |
IF (NN .NE. N .OR. MM .NE. M) THEN |
| 785 |
IER = -2 |
| 786 |
GOTO 999 |
| 787 |
ENDIF |
| 788 |
I = I + 1 |
| 789 |
GH(I) = G |
| 790 |
IF (M .NE. 0) THEN |
| 791 |
I = I + 1 |
| 792 |
GH(I) = H |
| 793 |
ENDIF |
| 794 |
2233 CONTINUE |
| 795 |
2211 CONTINUE |
| 796 |
c print *," guj" |
| 797 |
|
| 798 |
999 CLOSE (IU) |
| 799 |
c print *," guw IER",IER |
| 800 |
if ( IER .eq. -1 ) IER = 0 ! gfortran 4.1.2 bug workaround... hoping not to create problems with other versions |
| 801 |
|
| 802 |
RETURN |
| 803 |
END |
| 804 |
C |
| 805 |
C |
| 806 |
SUBROUTINE INTERSHC (DATE, DTE1, NMAX1, GH1, DTE2, |
| 807 |
1 NMAX2, GH2, NMAX, GH) |
| 808 |
|
| 809 |
C =============================================================== |
| 810 |
C |
| 811 |
C Version 1.01 |
| 812 |
C |
| 813 |
C Interpolates linearly, in time, between two spherical |
| 814 |
C harmonic models. |
| 815 |
C |
| 816 |
C Input: |
| 817 |
C DATE - Date of resulting model (in decimal year) |
| 818 |
C DTE1 - Date of earlier model |
| 819 |
C NMAX1 - Maximum degree and order of earlier model |
| 820 |
C GH1 - Schmidt quasi-normal internal spherical |
| 821 |
C harmonic coefficients of earlier model |
| 822 |
C DTE2 - Date of later model |
| 823 |
C NMAX2 - Maximum degree and order of later model |
| 824 |
C GH2 - Schmidt quasi-normal internal spherical |
| 825 |
C harmonic coefficients of later model |
| 826 |
C |
| 827 |
C Output: |
| 828 |
C GH - Coefficients of resulting model |
| 829 |
C NMAX - Maximum degree and order of resulting model |
| 830 |
C |
| 831 |
C A. Zunde |
| 832 |
C USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 |
| 833 |
C |
| 834 |
C =============================================================== |
| 835 |
|
| 836 |
DIMENSION GH1(*), GH2(*), GH(*) |
| 837 |
|
| 838 |
C --------------------------------------------------------------- |
| 839 |
C The coefficients (GH) of the resulting model, at date |
| 840 |
C DATE, are computed by linearly interpolating between the |
| 841 |
C coefficients of the earlier model (GH1), at date DTE1, |
| 842 |
C and those of the later model (GH2), at date DTE2. If one |
| 843 |
C model is smaller than the other, the interpolation is |
| 844 |
C performed with the missing coefficients assumed to be 0. |
| 845 |
C --------------------------------------------------------------- |
| 846 |
|
| 847 |
FACTOR = (DATE - DTE1) / (DTE2 - DTE1) |
| 848 |
|
| 849 |
IF (NMAX1 .EQ. NMAX2) THEN |
| 850 |
K = NMAX1 * (NMAX1 + 2) |
| 851 |
NMAX = NMAX1 |
| 852 |
ELSE IF (NMAX1 .GT. NMAX2) THEN |
| 853 |
K = NMAX2 * (NMAX2 + 2) |
| 854 |
L = NMAX1 * (NMAX1 + 2) |
| 855 |
DO 1122 I = K + 1, L |
| 856 |
1122 GH(I) = GH1(I) + FACTOR * (-GH1(I)) |
| 857 |
NMAX = NMAX1 |
| 858 |
ELSE |
| 859 |
K = NMAX1 * (NMAX1 + 2) |
| 860 |
L = NMAX2 * (NMAX2 + 2) |
| 861 |
DO 1133 I = K + 1, L |
| 862 |
1133 GH(I) = FACTOR * GH2(I) |
| 863 |
NMAX = NMAX2 |
| 864 |
ENDIF |
| 865 |
|
| 866 |
DO 1144 I = 1, K |
| 867 |
1144 GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I)) |
| 868 |
|
| 869 |
RETURN |
| 870 |
END |
| 871 |
C |
| 872 |
C |
| 873 |
SUBROUTINE EXTRASHC (DATE, DTE1, NMAX1, GH1, NMAX2, |
| 874 |
1 GH2, NMAX, GH) |
| 875 |
|
| 876 |
C =============================================================== |
| 877 |
C |
| 878 |
C Version 1.01 |
| 879 |
C |
| 880 |
C Extrapolates linearly a spherical harmonic model with a |
| 881 |
C rate-of-change model. |
| 882 |
C |
| 883 |
C Input: |
| 884 |
C DATE - Date of resulting model (in decimal year) |
| 885 |
C DTE1 - Date of base model |
| 886 |
C NMAX1 - Maximum degree and order of base model |
| 887 |
C GH1 - Schmidt quasi-normal internal spherical |
| 888 |
C harmonic coefficients of base model |
| 889 |
C NMAX2 - Maximum degree and order of rate-of-change |
| 890 |
C model |
| 891 |
C GH2 - Schmidt quasi-normal internal spherical |
| 892 |
C harmonic coefficients of rate-of-change model |
| 893 |
C |
| 894 |
C Output: |
| 895 |
C GH - Coefficients of resulting model |
| 896 |
C NMAX - Maximum degree and order of resulting model |
| 897 |
C |
| 898 |
C A. Zunde |
| 899 |
C USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 |
| 900 |
C |
| 901 |
C =============================================================== |
| 902 |
|
| 903 |
DIMENSION GH1(*), GH2(*), GH(*) |
| 904 |
|
| 905 |
C --------------------------------------------------------------- |
| 906 |
C The coefficients (GH) of the resulting model, at date |
| 907 |
C DATE, are computed by linearly extrapolating the coef- |
| 908 |
C ficients of the base model (GH1), at date DTE1, using |
| 909 |
C those of the rate-of-change model (GH2), at date DTE2. If |
| 910 |
C one model is smaller than the other, the extrapolation is |
| 911 |
C performed with the missing coefficients assumed to be 0. |
| 912 |
C --------------------------------------------------------------- |
| 913 |
|
| 914 |
FACTOR = (DATE - DTE1) |
| 915 |
|
| 916 |
IF (NMAX1 .EQ. NMAX2) THEN |
| 917 |
K = NMAX1 * (NMAX1 + 2) |
| 918 |
NMAX = NMAX1 |
| 919 |
ELSE IF (NMAX1 .GT. NMAX2) THEN |
| 920 |
K = NMAX2 * (NMAX2 + 2) |
| 921 |
L = NMAX1 * (NMAX1 + 2) |
| 922 |
DO 1155 I = K + 1, L |
| 923 |
1155 GH(I) = GH1(I) |
| 924 |
NMAX = NMAX1 |
| 925 |
ELSE |
| 926 |
K = NMAX1 * (NMAX1 + 2) |
| 927 |
L = NMAX2 * (NMAX2 + 2) |
| 928 |
DO 1166 I = K + 1, L |
| 929 |
1166 GH(I) = FACTOR * GH2(I) |
| 930 |
NMAX = NMAX2 |
| 931 |
ENDIF |
| 932 |
|
| 933 |
DO 1177 I = 1, K |
| 934 |
1177 GH(I) = GH1(I) + FACTOR * GH2(I) |
| 935 |
|
| 936 |
RETURN |
| 937 |
END |
| 938 |
C |
| 939 |
C |
| 940 |
SUBROUTINE INITIZE(TP1,TL1,TP2,TL2,TP3,TL3) |
| 941 |
C---------------------------------------------------------------- |
| 942 |
C Initializes the parameters in COMMON/GENER/ |
| 943 |
C |
| 944 |
C UMR = ATAN(1.0)*4./180. <DEGREE>*UMR=<RADIANT> |
| 945 |
C ERA EARTH RADIUS FOR NORMALIZATION OF CARTESIAN |
| 946 |
C COORDINATES (6371.2 KM) |
| 947 |
C EREQU MAJOR HALF AXIS FOR EARTH ELLIPSOID (6378.160 KM) |
| 948 |
C ERPOL MINOR HALF AXIS FOR EARTH ELLIPSOID (6356.775 KM) |
| 949 |
C AQUAD SQUARE OF MAJOR HALF AXIS FOR EARTH ELLIPSOID |
| 950 |
C BQUAD SQUARE OF MINOR HALF AXIS FOR EARTH ELLIPSOID |
| 951 |
C |
| 952 |
C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL |
| 953 |
C ASTRONOMICAL UNION . |
| 954 |
C----------------------------------------------------------------- |
| 955 |
INTEGER TL1,TL2,TL3 |
| 956 |
CHARACTER (len=*) :: TP1,TP2,TP3 |
| 957 |
INTEGER L1,L2,L3 |
| 958 |
CHARACTER *258 P1,P2,P3 |
| 959 |
COMMON/PPATH/ L1,L2,L3,P1, P2, P3 |
| 960 |
SAVE/PPATH/ |
| 961 |
|
| 962 |
COMMON/GENER/UMR,ERA,AQUAD,BQUAD |
| 963 |
SAVE/GENER/ |
| 964 |
|
| 965 |
L1 = TL1 |
| 966 |
L2 = TL2 |
| 967 |
L3 = TL3 |
| 968 |
P1 = TP1(1:L1) |
| 969 |
P2 = TP2(1:L2) |
| 970 |
P3 = TP3(1:L3) |
| 971 |
|
| 972 |
ERA=6371.2 |
| 973 |
EREQU=6378.16 |
| 974 |
ERPOL=6356.775 |
| 975 |
AQUAD=EREQU*EREQU |
| 976 |
BQUAD=ERPOL*ERPOL |
| 977 |
UMR=ATAN(1.0)*4./180. |
| 978 |
RETURN |
| 979 |
END |
Parent Directory
| 
Revision Log