1 |
kusanagi |
1.1 |
C BILCAL, VERSION 3.0, AUGUST 1995 |
2 |
|
|
C |
3 |
|
|
Cmm/dd/yy |
4 |
|
|
C 1/25/92-DKB-Modified for use with the IGRF-91 coefficients, which |
5 |
|
|
C were provided by R. Langel, GSFC. |
6 |
|
|
C 2/ 5/92-DKB-Reduce variable-name: INITI(ALI)ZE |
7 |
|
|
C 3/25/96-DKB-Modified for use with the IGRF-95 coefficients, which |
8 |
|
|
C were provided by R. Langel, GSFC. |
9 |
|
|
C 6/ 6/00-DKB-Modified for use with IGRF-2000 coefficients. |
10 |
|
|
C11/14/01-DKB-Add IMIN=0 above 4927 READ(...) [Rui Pereira] |
11 |
|
|
C04/25/05-DKB-IBBB instead of IBB in data statem. [Alexey Petrov] |
12 |
|
|
C |
13 |
|
|
C***************************************************************** |
14 |
|
|
C**************** IGRF MAGNETIC FIELD MODEL ********************* |
15 |
|
|
C**************** SHELLG L-VALUE CALCULATION ********************* |
16 |
|
|
C***************************************************************** |
17 |
|
|
C***************************************************************** |
18 |
|
|
C*** THIS PROGRAM PRODUCES PROFILES OF: *** |
19 |
|
|
C*** GEOMAGNETIC FIELD STRENGTH (GAUSS) *** |
20 |
|
|
C*** L-VALUE *** |
21 |
|
|
C***************************************************************** |
22 |
|
|
C*** FOR SPECIFIED: *** |
23 |
|
|
C*** YEAR (DECIMAL YEAR, E.G., 1995.5 FOR MID 1995) *** |
24 |
|
|
C*** GEODATIC LATITUDE AND LONGITUDE (DEGREE) *** |
25 |
|
|
C*** ALTITUDE (KM) *** |
26 |
|
|
C***************************************************************** |
27 |
|
|
C***************************************************************** |
28 |
|
|
C* --------------------ADDRESS-------------------------- * |
29 |
|
|
C* I DR. DIETER BILITZA (301)513-1664 I * |
30 |
|
|
C* I GSFC, NSSDC, CODE 933, GREENBELT, MD 20771, USA I * |
31 |
|
|
C* I SPAN: NSSDCA::BILITZA, NSSDC::BILITZA I * |
32 |
|
|
C* I BITNET: BILITZA%NSSDCA.SPAN@DFTNIC.BITNET I * |
33 |
|
|
C* ----------------------------------------------------- * |
34 |
|
|
C***************************************************************** |
35 |
|
|
C***************************************************************** |
36 |
|
|
C***************************************************************** |
37 |
|
|
INTEGER EGNR,AGNR,OGNR |
38 |
|
|
REAL LATI,LONGI,MLAT,MLONG |
39 |
|
|
CHARACTER*4 ITEXT(4),LTEX |
40 |
|
|
CHARACTER*7 ITB |
41 |
|
|
CHARACTER*11 NAME |
42 |
|
|
LOGICAL NOTBEG,VAL |
43 |
|
|
DIMENSION DEN(8),TEMP(2),XVAR(4),VARE(4),VARB(4) |
44 |
|
|
COMMON/GENER/ UMR,ERA,AQUAD,BQUAD |
45 |
|
|
C |
46 |
|
|
DATA ITEXT /'LATI','LONG','H/km','YEAR'/ |
47 |
|
|
DATA LATI,LONGI,HEIGHT,YEAR,IVAR,BVAR,EVAR,SVAR,IBBB,JAGNR |
48 |
|
|
& /45.1,293.1,100,1985.5,3,100,1000,100,0,2/ |
49 |
|
|
c### year limit modified |
50 |
|
|
DATA VARB /-90.0,-360.0,0.00000,1940.0/ |
51 |
|
|
DATA VARE /+90.0,+360.0,30000.0,2010.0/ |
52 |
|
|
C |
53 |
|
|
CALL INITIZE |
54 |
|
|
ALOG2=ALOG(2.) |
55 |
|
|
ISTART=1 |
56 |
|
|
C |
57 |
|
|
C FIRST SPECIFY YOUR COMPUTERS CHANNEL NUMBERS .................... |
58 |
|
|
C EGNR=INPUT, MONITO=MONITOR, KONSOL=MESSAGES...................... |
59 |
|
|
C AGNR=DISPLAY, OGNR=FILE OUTPUT................................... |
60 |
|
|
C |
61 |
|
|
EGNR=5 |
62 |
|
|
MONITO=6 |
63 |
|
|
AGNR=6 |
64 |
|
|
OGNR=16 |
65 |
|
|
WRITE(MONITO,5100) |
66 |
|
|
5100 FORMAT(1X/////4X,54('*')/4X, |
67 |
|
|
& '****** IGRF GEOMAGNETIC FIELD MODEL 1945 - 2005 ******'/4X, |
68 |
|
|
& '*********** SHELLG L-VALUE CALCULATION *************'/1X, |
69 |
|
|
& 60('*')/' This program allows you to produce B and L ', |
70 |
|
|
& 'profiles in '/' latitude, longitude, year or altitude.'/ |
71 |
|
|
& ' In each of the following windows you will be ', |
72 |
|
|
& 'asked to enter'/' one or more values, defining the conditions', |
73 |
|
|
& ' for your tables.'/' In each window the current value(s) is', |
74 |
|
|
& ' (are) shown in the right'/' upper corner (#...#). You can ', |
75 |
|
|
& 'choose the current values by'/' entering / at the prompt.'/ |
76 |
|
|
& ' If you enter a wrong character or a value outside the ', |
77 |
|
|
& 'allowed'/' parameter range, the program will ask you for a', |
78 |
|
|
& ' new entry.'/' After your tables are displayed, you can ', |
79 |
|
|
& 'change any parameter'/' you wish to change and create a ', |
80 |
|
|
& 'new profile.'/' You can leave the program at any point ', |
81 |
|
|
& 'by entering Ctrl Z.'/1X,25('*'),' GOOD LUCK ',25('*')) |
82 |
|
|
NOTBEG=.FALSE. |
83 |
|
|
GOTO 5508 |
84 |
|
|
C---------------START ENTERING PARAMETERS---------------------------- |
85 |
|
|
3293 CONTINUE |
86 |
|
|
ISTART=ISTART+1 |
87 |
|
|
C---------------WINDOW 1: WHICH PARAMETER CHANGE ?------------------- |
88 |
|
|
WRITE(MONITO,5602) LATI,LONGI,ITEXT(IVAR),HEIGHT,YEAR,BVAR, |
89 |
|
|
& EVAR,SVAR |
90 |
|
|
5602 FORMAT(1X//' **** WHICH PARAMETER DO YOU WANT TO CHANGE?'/ |
91 |
|
|
& 1X,60('-')/' 0 NO FURTHER CHANGES, CALCULATE PROFILE'/ |
92 |
|
|
& ' 1 LATITUDE #',F6.1,'#',7X,'5 DISPLAY OR STORE'/ |
93 |
|
|
& ' 2 LONGITUDE #',F6.1,'#',7X,'6 SELECTION OF VARIABLE #', |
94 |
|
|
& A4,'#'/' 3 ALTITUDE #',F8.1,'#',5X,'7 VARIABLE RANGE'/ |
95 |
|
|
& ' 4 YEAR #',F6.1,'#',11X,'#',F8.1,',',F8.1,',', |
96 |
|
|
& F8.1,'#'/29X,'8 B OR B/B0'/1X,60('-')/' ENTER NUMBER') |
97 |
|
|
IMIN=0 |
98 |
|
|
IMAX=8 |
99 |
|
|
4925 READ(EGNR,*,ERR=8600,END=6666) ISWIT |
100 |
|
|
IF((ISWIT.GE.IMIN).AND.(ISWIT.LE.IMAX)) GOTO 8601 |
101 |
|
|
8600 WRITE(MONITO,4924) IMIN,IMAX |
102 |
|
|
4924 FORMAT(' Your input is outside the value range:', |
103 |
|
|
& I2,' to',I2/' try again') |
104 |
|
|
GOTO 4925 |
105 |
|
|
8601 GOTO (5505,3329,3339,5502,6780,5508,5503,5504,9138) ISWIT+1 |
106 |
|
|
C--------------WINDOW 2: DISPLAY OPTIONS-------------------------- |
107 |
|
|
5508 WRITE(MONITO,5123) JAGNR |
108 |
|
|
5123 FORMAT(/' DO YOU WANT YOUR PROFILES',32X,'#',I1,'#'/5X, |
109 |
|
|
& 'DISPLAYED ON YOUR MONITOR: ENTER 0 AT PROMPT'/5X, |
110 |
|
|
& 'STORED IN FILE OUTPUT.IGR: ENTER 1 AT PROMPT'/5X, |
111 |
|
|
& 'DISPLAYED AND STORED: ENTER 2 AT PROMPT') |
112 |
|
|
WRITE(MONITO,8630) |
113 |
|
|
IMAX=2 |
114 |
|
|
IMIN=0 |
115 |
|
|
4927 READ(EGNR,*,ERR=8603,END=6666) JAGNR |
116 |
|
|
IF((JAGNR.GE.IMIN).AND.(JAGNR.LE.IMAX)) GOTO 8602 |
117 |
|
|
8603 WRITE(MONITO,4924) IMIN,IMAX |
118 |
|
|
GOTO 4927 |
119 |
|
|
8602 IVARNR=0 |
120 |
|
|
IF(JAGNR.GT.0) |
121 |
|
|
& OPEN(UNIT=OGNR,FILE='OUTPUT.IGR',STATUS='NEW',FORM='FORMATTED') |
122 |
|
|
IF(JAGNR.EQ.1) AGNR=OGNR |
123 |
|
|
IF(NOTBEG) GOTO 3293 |
124 |
|
|
C---------------WINDOW 3: SELECT VARIABLE------------------------ |
125 |
|
|
5503 WRITE(MONITO,5040) IVAR |
126 |
|
|
5040 FORMAT(1X//' SELECT YOUR VARIABLE:',31X,'#LAST:',I1,'#'// |
127 |
|
|
& ' 1 LATITUDE 3 ALTITUDE'/ |
128 |
|
|
& ' 2 LONGITUDE 4 YEAR') |
129 |
|
|
WRITE(MONITO,8630) |
130 |
|
|
8630 FORMAT(1X,60('-')/' Enter / to use previous value(s) ', |
131 |
|
|
& '(see # .. #); Ctrl Z to exit') |
132 |
|
|
IMIN=1 |
133 |
|
|
IMAX=4 |
134 |
|
|
4929 READ(EGNR,*,ERR=8605,END=6666) IVAR |
135 |
|
|
IF((IVAR.GE.IMIN).AND.(IVAR.LE.IMAX)) GOTO 5504 |
136 |
|
|
8605 WRITE(MONITO,4924) IMIN,IMAX |
137 |
|
|
GOTO 4929 |
138 |
|
|
C--------------WINDOW 4: SELECT VARIABLE RANGE--------------------- |
139 |
|
|
5504 WRITE(MONITO,5044) BVAR,EVAR,SVAR |
140 |
|
|
5044 FORMAT(1X//' CHOOSE YOUR VARIABLE RANGE:',5X,' BEGIN, END, ', |
141 |
|
|
& 'STEPWIDTH ?'/32X,'#',F8.1,',',F8.1,',',F8.1,'#') |
142 |
|
|
WRITE(MONITO,8630) |
143 |
|
|
VAMIN=VARB(IVAR) |
144 |
|
|
VAMAX=VARE(IVAR) |
145 |
|
|
4931 READ(EGNR,*,ERR=8606,END=6666) BVAR,EVAR,SVAR |
146 |
|
|
IF((BVAR.GE.VAMIN).AND.(EVAR.LE.VAMAX)) GOTO 8607 |
147 |
|
|
8606 WRITE(MONITO,4930) VAMIN,VAMAX |
148 |
|
|
4930 FORMAT(' Your input is outside the value range:', |
149 |
|
|
& F8.1,' to',F8.1/' try again') |
150 |
|
|
GOTO 4931 |
151 |
|
|
8607 LANZ=INT((EVAR-BVAR)/SVAR)+1 |
152 |
|
|
IF(NOTBEG) GOTO 3293 |
153 |
|
|
IVARNR=IVARNR+1 |
154 |
|
|
IF(IVARNR.EQ.IVAR) GOTO 7339 |
155 |
|
|
C--------------WINDOW 5: LATITUDE----------------------------------- |
156 |
|
|
3329 WRITE(MONITO,5000) LATI |
157 |
|
|
5000 FORMAT(1X//1X,'GEOD LATITUDE ? !NORTH! [DEGREE,DECIMAL]', |
158 |
|
|
& 8X,'#',F5.1,'#') |
159 |
|
|
WRITE(MONITO,8630) |
160 |
|
|
XMAX=VARE(1) |
161 |
|
|
XMIN=VARB(1) |
162 |
|
|
4933 READ(EGNR,*,ERR=8608,END=6666) LATI |
163 |
|
|
IF((LATI.GE.XMIN).AND.(LATI.LE.XMAX)) GOTO 8609 |
164 |
|
|
8608 WRITE(MONITO,4930) XMIN,XMAX |
165 |
|
|
GOTO 4933 |
166 |
|
|
8609 IF(NOTBEG) GOTO 3293 |
167 |
|
|
7339 IVARNR=IVARNR+1 |
168 |
|
|
IF(IVARNR.EQ.IVAR) GOTO 7500 |
169 |
|
|
C---------------WINDOW 6: LONGITUDE--------------------------------- |
170 |
|
|
3339 WRITE(MONITO,6001) LONGI |
171 |
|
|
6001 FORMAT(1X//1X,'GEOD LONGITUDE ? !EAST! [DEGREE,DECIMAL]', |
172 |
|
|
& 7X,'#',F6.1,'#') |
173 |
|
|
WRITE(MONITO,8630) |
174 |
|
|
XMAX=VARE(2) |
175 |
|
|
XMIN=VARB(2) |
176 |
|
|
4934 READ(EGNR,*,ERR=8610,END=6666) LONGI |
177 |
|
|
IF((LONGI.GE.XMIN).AND.(LONGI.LE.XMAX)) GOTO 8611 |
178 |
|
|
8610 WRITE(MONITO,4930) XMIN,XMAX |
179 |
|
|
GOTO 4934 |
180 |
|
|
8611 IF(NOTBEG) GOTO 3293 |
181 |
|
|
7500 IVARNR=IVARNR+1 |
182 |
|
|
IF(IVARNR.EQ.IVAR) GOTO 5551 |
183 |
|
|
C---------------WINDOW 7: ALTITUDE--------------------------------- |
184 |
|
|
5502 WRITE(MONITO,6002) HEIGHT |
185 |
|
|
6002 FORMAT(1X//1X,'ALTITUDE ? [KM]',33X,'#',F7.1,'#') |
186 |
|
|
WRITE(MONITO,8630) |
187 |
|
|
XMAX=VARE(3) |
188 |
|
|
XMIN=VARB(3) |
189 |
|
|
4936 READ(EGNR,*,ERR=8615,END=6666) HEIGHT |
190 |
|
|
IF((HEIGHT.GE.XMIN).AND.(HEIGHT.LE.XMAX)) GOTO 8616 |
191 |
|
|
8615 WRITE(MONITO,4930) XMIN,XMAX |
192 |
|
|
GOTO 4936 |
193 |
|
|
8616 IF(NOTBEG) GOTO 3293 |
194 |
|
|
5551 IVARNR=IVARNR+1 |
195 |
|
|
IF(IVARNR.EQ.IVAR) GOTO 9138 |
196 |
|
|
C----------------WINDOW 8: YEAR------------------------------------ |
197 |
|
|
6780 WRITE(MONITO,6004) YEAR |
198 |
|
|
6004 FORMAT(1X//' YEAR(EPOCH) ?',9X,'*decimal*',9X,'#',F6.1,'#') |
199 |
|
|
WRITE(MONITO,8630) |
200 |
|
|
XMAX=VARE(4) |
201 |
|
|
XMIN=VARB(4) |
202 |
|
|
4938 READ(EGNR,*,ERR=8617,END=6666) YEAR |
203 |
|
|
IF((YEAR.GE.XMIN).AND.(YEAR.LE.XMAX)) GOTO 8618 |
204 |
|
|
8617 WRITE(MONITO,4930) XMIN,XMAX |
205 |
|
|
GOTO 4938 |
206 |
|
|
8618 IF(NOTBEG) GOTO 3293 |
207 |
|
|
C----------------WINDOW 9: ABSOLUTE OR NORMALIZED B-------------- |
208 |
|
|
9138 WRITE(MONITO,6204) IBBB |
209 |
|
|
6204 FORMAT(1X//' OUTPUT OPTION: B OR B/B0 ?',19X,'#',I1,'#'// |
210 |
|
|
& 4X,'if you enter 0, the absolute magnetic field strength'/ |
211 |
|
|
& 4X,'will be listed, otherwise the field strength normalized'/ |
212 |
|
|
& 4X,'to the field strength at the magnetic equator is listed') |
213 |
|
|
WRITE(MONITO,8630) |
214 |
|
|
4738 READ(EGNR,*,ERR=8717,END=6666) IBBB |
215 |
|
|
IF(IBBB.NE.0) THEN |
216 |
|
|
ITB=' B/B0 ' |
217 |
|
|
ELSE |
218 |
|
|
ITB='B/Gauss' |
219 |
|
|
ENDIF |
220 |
|
|
GOTO 8718 |
221 |
|
|
8717 WRITE(MONITO,4630) |
222 |
|
|
4630 FORMAT(' Your input should be a integer value'/' try again') |
223 |
|
|
GOTO 4738 |
224 |
|
|
8718 IF(NOTBEG) GOTO 3293 |
225 |
|
|
C----------------CALCULATE PROFILES----------------------------------- |
226 |
|
|
5505 WRITE(AGNR,3910) ITEXT(IVAR),ITB |
227 |
|
|
IF(JAGNR.EQ.2) WRITE(OGNR,3910) ITEXT(IVAR),ITB |
228 |
|
|
3910 FORMAT(1X//////////// |
229 |
|
|
& 5X,A4,' DIMO ',A7,' B-NORTH B-EAST B-DOWN ', |
230 |
|
|
& ' DIP DEC L-VALUE C') |
231 |
|
|
XVAR(1)=LATI |
232 |
|
|
XVAR(2)=LONGI |
233 |
|
|
XVAR(3)=HEIGHT |
234 |
|
|
XVAR(4)=YEAR |
235 |
|
|
LFD=0 |
236 |
|
|
XVAR(IVAR)=BVAR-SVAR |
237 |
|
|
2123 XVAR(IVAR)=XVAR(IVAR)+SVAR |
238 |
|
|
LFD=LFD+1 |
239 |
|
|
LATI=XVAR(1) |
240 |
|
|
LONGI=XVAR(2) |
241 |
|
|
HEIGHT=XVAR(3) |
242 |
|
|
YEAR=XVAR(4) |
243 |
|
|
IF((IVAR.LT.4).AND.(LFD.GT.1)) GOTO 2910 |
244 |
|
|
CALL FELDCOF(YEAR,DIMO) |
245 |
|
|
2910 CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS) |
246 |
|
|
CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1) |
247 |
|
|
IF(IABS(ICODE).GT.9) THEN |
248 |
|
|
WRITE(MONITO,7227) ICODE |
249 |
|
|
7227 FORMAT(' ICODE=',I10,' is set to 2') |
250 |
|
|
ICODE=2 |
251 |
|
|
ENDIF |
252 |
|
|
IF(IBBB.EQ.0) GOTO 2299 |
253 |
|
|
BEQU=DIMO/(XL*XL*XL) |
254 |
|
|
IF(ICODE.EQ.1) THEN |
255 |
|
|
BDEL=1.E-3 |
256 |
|
|
CALL FINDB0(0.05,BDEL,VAL,BEQ,RR0) |
257 |
|
|
IF(VAL) BEQU=BEQ |
258 |
|
|
ENDIF |
259 |
|
|
2299 DIP=ASIN(BDOWN/BABS)/UMR |
260 |
|
|
DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR |
261 |
|
|
XCOR=XVAR(IVAR) |
262 |
|
|
IF(IBBB.EQ.0) THEN |
263 |
|
|
WRITE(AGNR,7117) XCOR,DIMO,BABS,BNORTH,BEAST,BDOWN, |
264 |
|
|
& DIP,DEC,XL,ICODE |
265 |
|
|
IF(JAGNR.EQ.2) WRITE(OGNR,7117) XCOR,DIMO,BABS,BNORTH, |
266 |
|
|
& BEAST,BDOWN,DIP,DEC,XL,ICODE |
267 |
|
|
7117 FORMAT(1X,F8.2,F8.4,4(1X,F7.5),2F7.1,F8.3,I3) |
268 |
|
|
ELSE |
269 |
|
|
BBX=BABS/BEQU |
270 |
|
|
IF(BBX.GT.9999.999) BBX=9999.999 |
271 |
|
|
WRITE(AGNR,7177) XCOR,DIMO,BBX,BNORTH,BEAST,BDOWN, |
272 |
|
|
& DIP,DEC,XL,ICODE |
273 |
|
|
IF(JAGNR.EQ.2) WRITE(OGNR,7177) XCOR,DIMO,BBX,BNORTH, |
274 |
|
|
& BEAST,BDOWN,DIP,DEC,XL,ICODE |
275 |
|
|
7177 FORMAT(1X,F8.2,F8.4,F8.3,3(1X,F7.5),2F7.1,F8.3,I3) |
276 |
|
|
ENDIF |
277 |
|
|
IF(XCOR.LT.EVAR) GOTO 2123 |
278 |
|
|
WRITE(AGNR,2193) LATI,LONGI,HEIGHT,YEAR |
279 |
|
|
IF(JAGNR.EQ.2) WRITE(OGNR,2193) LATI,LONGI,HEIGHT,YEAR |
280 |
|
|
C ### edition date corrected |
281 |
|
|
2193 FORMAT(1X,'------- International Geomagnetic Reference Field', |
282 |
|
|
& ' --- Edition 2000 -------'/' LATI=',F7.1,' LONGI=',F6.1, |
283 |
|
|
& ' I DIMO is Dipol I C=1 L and B0 correct'/ |
284 |
|
|
& ' ALT=',F7.1,' YEAR=',F6.1,' I Moment in Gauss', |
285 |
|
|
& ' I =2 wrong, =3 approx.'/1X,74('-')) |
286 |
|
|
IF(HEIGHT.GT.5000.0) THEN |
287 |
|
|
WRITE(AGNR,5611) |
288 |
|
|
IF(JAGNR.EQ.2) WRITE(OGNR,5611) |
289 |
|
|
ENDIF |
290 |
|
|
5611 FORMAT(' !! REMINDER: this field model does not', |
291 |
|
|
& ' include external sources !!') |
292 |
|
|
C ### year limits corrected |
293 |
|
|
IF((YEAR.LT.1945.0).OR.(YEAR.GT.2005.0)) THEN |
294 |
|
|
WRITE(AGNR,5612) |
295 |
|
|
IF(JAGNR.EQ.2) WRITE(OGNR,5612) |
296 |
|
|
ENDIF |
297 |
|
|
C ### timeperiod corrected |
298 |
|
|
5612 FORMAT(' !! REMINDER: Recommended time period is 1945', |
299 |
|
|
& ' to 2005 !!') |
300 |
|
|
C-----------------LAST WINDOW: CONTINUE ?----------------------- |
301 |
|
|
918 WRITE(MONITO,5600) |
302 |
|
|
5600 FORMAT(1X/' **** DO YOU WANT TO CONTINUE?'/1X,60('-')/ |
303 |
|
|
& ' "0" QUIT AND EXIT "1" NEW PARAMETERS'/ |
304 |
|
|
& 1X,60('-')) |
305 |
|
|
IMIN=0 |
306 |
|
|
IMAX=1 |
307 |
|
|
8651 READ(EGNR,*,ERR=8652,END=6666) IALL |
308 |
|
|
IF((IALL.GE.IMIN).AND.(IALL.LE.IMAX)) GOTO 8653 |
309 |
|
|
8652 WRITE(MONITO,4924) IMIN,IMAX |
310 |
|
|
GOTO 8651 |
311 |
|
|
8653 NOTBEG=.TRUE. |
312 |
|
|
IF(IALL.EQ.1) GOTO 3293 |
313 |
|
|
6666 CONTINUE |
314 |
|
|
STOP |
315 |
|
|
END |