| 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 |
Parent Directory
| 
Revision Log