C BILCAL, VERSION 3.0, AUGUST 1995 C Cmm/dd/yy C 1/25/92-DKB-Modified for use with the IGRF-91 coefficients, which C were provided by R. Langel, GSFC. C 2/ 5/92-DKB-Reduce variable-name: INITI(ALI)ZE C 3/25/96-DKB-Modified for use with the IGRF-95 coefficients, which C were provided by R. Langel, GSFC. C 6/ 6/00-DKB-Modified for use with IGRF-2000 coefficients. C11/14/01-DKB-Add IMIN=0 above 4927 READ(...) [Rui Pereira] C04/25/05-DKB-IBBB instead of IBB in data statem. [Alexey Petrov] C C***************************************************************** C**************** IGRF MAGNETIC FIELD MODEL ********************* C**************** SHELLG L-VALUE CALCULATION ********************* C***************************************************************** C***************************************************************** C*** THIS PROGRAM PRODUCES PROFILES OF: *** C*** GEOMAGNETIC FIELD STRENGTH (GAUSS) *** C*** L-VALUE *** C***************************************************************** C*** FOR SPECIFIED: *** C*** YEAR (DECIMAL YEAR, E.G., 1995.5 FOR MID 1995) *** C*** GEODATIC LATITUDE AND LONGITUDE (DEGREE) *** C*** ALTITUDE (KM) *** C***************************************************************** C***************************************************************** C* --------------------ADDRESS-------------------------- * C* I DR. DIETER BILITZA (301)513-1664 I * C* I GSFC, NSSDC, CODE 933, GREENBELT, MD 20771, USA I * C* I SPAN: NSSDCA::BILITZA, NSSDC::BILITZA I * C* I BITNET: BILITZA%NSSDCA.SPAN@DFTNIC.BITNET I * C* ----------------------------------------------------- * C***************************************************************** C***************************************************************** C***************************************************************** INTEGER EGNR,AGNR,OGNR REAL LATI,LONGI,MLAT,MLONG CHARACTER*4 ITEXT(4),LTEX CHARACTER*7 ITB CHARACTER*11 NAME LOGICAL NOTBEG,VAL DIMENSION DEN(8),TEMP(2),XVAR(4),VARE(4),VARB(4) COMMON/GENER/ UMR,ERA,AQUAD,BQUAD C DATA ITEXT /'LATI','LONG','H/km','YEAR'/ DATA LATI,LONGI,HEIGHT,YEAR,IVAR,BVAR,EVAR,SVAR,IBBB,JAGNR & /45.1,293.1,100,1985.5,3,100,1000,100,0,2/ c### year limit modified DATA VARB /-90.0,-360.0,0.00000,1940.0/ DATA VARE /+90.0,+360.0,30000.0,2010.0/ C CALL INITIZE ALOG2=ALOG(2.) ISTART=1 C C FIRST SPECIFY YOUR COMPUTERS CHANNEL NUMBERS .................... C EGNR=INPUT, MONITO=MONITOR, KONSOL=MESSAGES...................... C AGNR=DISPLAY, OGNR=FILE OUTPUT................................... C EGNR=5 MONITO=6 AGNR=6 OGNR=16 WRITE(MONITO,5100) 5100 FORMAT(1X/////4X,54('*')/4X, & '****** IGRF GEOMAGNETIC FIELD MODEL 1945 - 2005 ******'/4X, & '*********** SHELLG L-VALUE CALCULATION *************'/1X, & 60('*')/' This program allows you to produce B and L ', & 'profiles in '/' latitude, longitude, year or altitude.'/ & ' In each of the following windows you will be ', & 'asked to enter'/' one or more values, defining the conditions', & ' for your tables.'/' In each window the current value(s) is', & ' (are) shown in the right'/' upper corner (#...#). You can ', & 'choose the current values by'/' entering / at the prompt.'/ & ' If you enter a wrong character or a value outside the ', & 'allowed'/' parameter range, the program will ask you for a', & ' new entry.'/' After your tables are displayed, you can ', & 'change any parameter'/' you wish to change and create a ', & 'new profile.'/' You can leave the program at any point ', & 'by entering Ctrl Z.'/1X,25('*'),' GOOD LUCK ',25('*')) NOTBEG=.FALSE. GOTO 5508 C---------------START ENTERING PARAMETERS---------------------------- 3293 CONTINUE ISTART=ISTART+1 C---------------WINDOW 1: WHICH PARAMETER CHANGE ?------------------- WRITE(MONITO,5602) LATI,LONGI,ITEXT(IVAR),HEIGHT,YEAR,BVAR, & EVAR,SVAR 5602 FORMAT(1X//' **** WHICH PARAMETER DO YOU WANT TO CHANGE?'/ & 1X,60('-')/' 0 NO FURTHER CHANGES, CALCULATE PROFILE'/ & ' 1 LATITUDE #',F6.1,'#',7X,'5 DISPLAY OR STORE'/ & ' 2 LONGITUDE #',F6.1,'#',7X,'6 SELECTION OF VARIABLE #', & A4,'#'/' 3 ALTITUDE #',F8.1,'#',5X,'7 VARIABLE RANGE'/ & ' 4 YEAR #',F6.1,'#',11X,'#',F8.1,',',F8.1,',', & F8.1,'#'/29X,'8 B OR B/B0'/1X,60('-')/' ENTER NUMBER') IMIN=0 IMAX=8 4925 READ(EGNR,*,ERR=8600,END=6666) ISWIT IF((ISWIT.GE.IMIN).AND.(ISWIT.LE.IMAX)) GOTO 8601 8600 WRITE(MONITO,4924) IMIN,IMAX 4924 FORMAT(' Your input is outside the value range:', & I2,' to',I2/' try again') GOTO 4925 8601 GOTO (5505,3329,3339,5502,6780,5508,5503,5504,9138) ISWIT+1 C--------------WINDOW 2: DISPLAY OPTIONS-------------------------- 5508 WRITE(MONITO,5123) JAGNR 5123 FORMAT(/' DO YOU WANT YOUR PROFILES',32X,'#',I1,'#'/5X, & 'DISPLAYED ON YOUR MONITOR: ENTER 0 AT PROMPT'/5X, & 'STORED IN FILE OUTPUT.IGR: ENTER 1 AT PROMPT'/5X, & 'DISPLAYED AND STORED: ENTER 2 AT PROMPT') WRITE(MONITO,8630) IMAX=2 IMIN=0 4927 READ(EGNR,*,ERR=8603,END=6666) JAGNR IF((JAGNR.GE.IMIN).AND.(JAGNR.LE.IMAX)) GOTO 8602 8603 WRITE(MONITO,4924) IMIN,IMAX GOTO 4927 8602 IVARNR=0 IF(JAGNR.GT.0) & OPEN(UNIT=OGNR,FILE='OUTPUT.IGR',STATUS='NEW',FORM='FORMATTED') IF(JAGNR.EQ.1) AGNR=OGNR IF(NOTBEG) GOTO 3293 C---------------WINDOW 3: SELECT VARIABLE------------------------ 5503 WRITE(MONITO,5040) IVAR 5040 FORMAT(1X//' SELECT YOUR VARIABLE:',31X,'#LAST:',I1,'#'// & ' 1 LATITUDE 3 ALTITUDE'/ & ' 2 LONGITUDE 4 YEAR') WRITE(MONITO,8630) 8630 FORMAT(1X,60('-')/' Enter / to use previous value(s) ', & '(see # .. #); Ctrl Z to exit') IMIN=1 IMAX=4 4929 READ(EGNR,*,ERR=8605,END=6666) IVAR IF((IVAR.GE.IMIN).AND.(IVAR.LE.IMAX)) GOTO 5504 8605 WRITE(MONITO,4924) IMIN,IMAX GOTO 4929 C--------------WINDOW 4: SELECT VARIABLE RANGE--------------------- 5504 WRITE(MONITO,5044) BVAR,EVAR,SVAR 5044 FORMAT(1X//' CHOOSE YOUR VARIABLE RANGE:',5X,' BEGIN, END, ', & 'STEPWIDTH ?'/32X,'#',F8.1,',',F8.1,',',F8.1,'#') WRITE(MONITO,8630) VAMIN=VARB(IVAR) VAMAX=VARE(IVAR) 4931 READ(EGNR,*,ERR=8606,END=6666) BVAR,EVAR,SVAR IF((BVAR.GE.VAMIN).AND.(EVAR.LE.VAMAX)) GOTO 8607 8606 WRITE(MONITO,4930) VAMIN,VAMAX 4930 FORMAT(' Your input is outside the value range:', & F8.1,' to',F8.1/' try again') GOTO 4931 8607 LANZ=INT((EVAR-BVAR)/SVAR)+1 IF(NOTBEG) GOTO 3293 IVARNR=IVARNR+1 IF(IVARNR.EQ.IVAR) GOTO 7339 C--------------WINDOW 5: LATITUDE----------------------------------- 3329 WRITE(MONITO,5000) LATI 5000 FORMAT(1X//1X,'GEOD LATITUDE ? !NORTH! [DEGREE,DECIMAL]', & 8X,'#',F5.1,'#') WRITE(MONITO,8630) XMAX=VARE(1) XMIN=VARB(1) 4933 READ(EGNR,*,ERR=8608,END=6666) LATI IF((LATI.GE.XMIN).AND.(LATI.LE.XMAX)) GOTO 8609 8608 WRITE(MONITO,4930) XMIN,XMAX GOTO 4933 8609 IF(NOTBEG) GOTO 3293 7339 IVARNR=IVARNR+1 IF(IVARNR.EQ.IVAR) GOTO 7500 C---------------WINDOW 6: LONGITUDE--------------------------------- 3339 WRITE(MONITO,6001) LONGI 6001 FORMAT(1X//1X,'GEOD LONGITUDE ? !EAST! [DEGREE,DECIMAL]', & 7X,'#',F6.1,'#') WRITE(MONITO,8630) XMAX=VARE(2) XMIN=VARB(2) 4934 READ(EGNR,*,ERR=8610,END=6666) LONGI IF((LONGI.GE.XMIN).AND.(LONGI.LE.XMAX)) GOTO 8611 8610 WRITE(MONITO,4930) XMIN,XMAX GOTO 4934 8611 IF(NOTBEG) GOTO 3293 7500 IVARNR=IVARNR+1 IF(IVARNR.EQ.IVAR) GOTO 5551 C---------------WINDOW 7: ALTITUDE--------------------------------- 5502 WRITE(MONITO,6002) HEIGHT 6002 FORMAT(1X//1X,'ALTITUDE ? [KM]',33X,'#',F7.1,'#') WRITE(MONITO,8630) XMAX=VARE(3) XMIN=VARB(3) 4936 READ(EGNR,*,ERR=8615,END=6666) HEIGHT IF((HEIGHT.GE.XMIN).AND.(HEIGHT.LE.XMAX)) GOTO 8616 8615 WRITE(MONITO,4930) XMIN,XMAX GOTO 4936 8616 IF(NOTBEG) GOTO 3293 5551 IVARNR=IVARNR+1 IF(IVARNR.EQ.IVAR) GOTO 9138 C----------------WINDOW 8: YEAR------------------------------------ 6780 WRITE(MONITO,6004) YEAR 6004 FORMAT(1X//' YEAR(EPOCH) ?',9X,'*decimal*',9X,'#',F6.1,'#') WRITE(MONITO,8630) XMAX=VARE(4) XMIN=VARB(4) 4938 READ(EGNR,*,ERR=8617,END=6666) YEAR IF((YEAR.GE.XMIN).AND.(YEAR.LE.XMAX)) GOTO 8618 8617 WRITE(MONITO,4930) XMIN,XMAX GOTO 4938 8618 IF(NOTBEG) GOTO 3293 C----------------WINDOW 9: ABSOLUTE OR NORMALIZED B-------------- 9138 WRITE(MONITO,6204) IBBB 6204 FORMAT(1X//' OUTPUT OPTION: B OR B/B0 ?',19X,'#',I1,'#'// & 4X,'if you enter 0, the absolute magnetic field strength'/ & 4X,'will be listed, otherwise the field strength normalized'/ & 4X,'to the field strength at the magnetic equator is listed') WRITE(MONITO,8630) 4738 READ(EGNR,*,ERR=8717,END=6666) IBBB IF(IBBB.NE.0) THEN ITB=' B/B0 ' ELSE ITB='B/Gauss' ENDIF GOTO 8718 8717 WRITE(MONITO,4630) 4630 FORMAT(' Your input should be a integer value'/' try again') GOTO 4738 8718 IF(NOTBEG) GOTO 3293 C----------------CALCULATE PROFILES----------------------------------- 5505 WRITE(AGNR,3910) ITEXT(IVAR),ITB IF(JAGNR.EQ.2) WRITE(OGNR,3910) ITEXT(IVAR),ITB 3910 FORMAT(1X//////////// & 5X,A4,' DIMO ',A7,' B-NORTH B-EAST B-DOWN ', & ' DIP DEC L-VALUE C') XVAR(1)=LATI XVAR(2)=LONGI XVAR(3)=HEIGHT XVAR(4)=YEAR LFD=0 XVAR(IVAR)=BVAR-SVAR 2123 XVAR(IVAR)=XVAR(IVAR)+SVAR LFD=LFD+1 LATI=XVAR(1) LONGI=XVAR(2) HEIGHT=XVAR(3) YEAR=XVAR(4) IF((IVAR.LT.4).AND.(LFD.GT.1)) GOTO 2910 CALL FELDCOF(YEAR,DIMO) 2910 CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS) CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1) IF(IABS(ICODE).GT.9) THEN WRITE(MONITO,7227) ICODE 7227 FORMAT(' ICODE=',I10,' is set to 2') ICODE=2 ENDIF IF(IBBB.EQ.0) GOTO 2299 BEQU=DIMO/(XL*XL*XL) IF(ICODE.EQ.1) THEN BDEL=1.E-3 CALL FINDB0(0.05,BDEL,VAL,BEQ,RR0) IF(VAL) BEQU=BEQ ENDIF 2299 DIP=ASIN(BDOWN/BABS)/UMR DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR XCOR=XVAR(IVAR) IF(IBBB.EQ.0) THEN WRITE(AGNR,7117) XCOR,DIMO,BABS,BNORTH,BEAST,BDOWN, & DIP,DEC,XL,ICODE IF(JAGNR.EQ.2) WRITE(OGNR,7117) XCOR,DIMO,BABS,BNORTH, & BEAST,BDOWN,DIP,DEC,XL,ICODE 7117 FORMAT(1X,F8.2,F8.4,4(1X,F7.5),2F7.1,F8.3,I3) ELSE BBX=BABS/BEQU IF(BBX.GT.9999.999) BBX=9999.999 WRITE(AGNR,7177) XCOR,DIMO,BBX,BNORTH,BEAST,BDOWN, & DIP,DEC,XL,ICODE IF(JAGNR.EQ.2) WRITE(OGNR,7177) XCOR,DIMO,BBX,BNORTH, & BEAST,BDOWN,DIP,DEC,XL,ICODE 7177 FORMAT(1X,F8.2,F8.4,F8.3,3(1X,F7.5),2F7.1,F8.3,I3) ENDIF IF(XCOR.LT.EVAR) GOTO 2123 WRITE(AGNR,2193) LATI,LONGI,HEIGHT,YEAR IF(JAGNR.EQ.2) WRITE(OGNR,2193) LATI,LONGI,HEIGHT,YEAR C ### edition date corrected 2193 FORMAT(1X,'------- International Geomagnetic Reference Field', & ' --- Edition 2000 -------'/' LATI=',F7.1,' LONGI=',F6.1, & ' I DIMO is Dipol I C=1 L and B0 correct'/ & ' ALT=',F7.1,' YEAR=',F6.1,' I Moment in Gauss', & ' I =2 wrong, =3 approx.'/1X,74('-')) IF(HEIGHT.GT.5000.0) THEN WRITE(AGNR,5611) IF(JAGNR.EQ.2) WRITE(OGNR,5611) ENDIF 5611 FORMAT(' !! REMINDER: this field model does not', & ' include external sources !!') C ### year limits corrected IF((YEAR.LT.1945.0).OR.(YEAR.GT.2005.0)) THEN WRITE(AGNR,5612) IF(JAGNR.EQ.2) WRITE(OGNR,5612) ENDIF C ### timeperiod corrected 5612 FORMAT(' !! REMINDER: Recommended time period is 1945', & ' to 2005 !!') C-----------------LAST WINDOW: CONTINUE ?----------------------- 918 WRITE(MONITO,5600) 5600 FORMAT(1X/' **** DO YOU WANT TO CONTINUE?'/1X,60('-')/ & ' "0" QUIT AND EXIT "1" NEW PARAMETERS'/ & 1X,60('-')) IMIN=0 IMAX=1 8651 READ(EGNR,*,ERR=8652,END=6666) IALL IF((IALL.GE.IMIN).AND.(IALL.LE.IMAX)) GOTO 8653 8652 WRITE(MONITO,4924) IMIN,IMAX GOTO 8651 8653 NOTBEG=.TRUE. IF(IALL.EQ.1) GOTO 3293 6666 CONTINUE STOP END