/[PAMELA software]/yodaUtility/igrf/bilcal.for
ViewVC logotype

Contents of /yodaUtility/igrf/bilcal.for

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (show annotations) (download)
Sat Jun 10 06:33:31 2006 UTC (18 years, 6 months ago) by kusanagi
Branch: MAIN
CVS Tags: yodaUtility2_0/00, yodaUtility2_2/00, yodaUtility2_1/00, HEAD
Software for geomagnetic field calculation. The validity time interval is 2000-2010.

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

  ViewVC Help
Powered by ViewVC 1.1.23