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

Annotation of /yodaUtility/igrf/bilcal.for

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (hide 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 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

  ViewVC Help
Powered by ViewVC 1.1.23