igrf10-02.f
PROGRAM IGRF10
C
C This is a program for synthesising geomagnetic field values from the
C International Geomagnetic Reference Field series of models as agreed
c in Decmember 2004 by IAGA Working Group V-MOD.
C It is the 10th generation IGRF, ie the 9th revision.
C The main-field models for 1900.0, 1905.0,..1940.0 and 2005.0 are
C non-definitive, those for 1945.0, 1950.0,...2000.0 are definitive and
C the secular-variation model for 2005.0 to 2010.0 is non-definitive.
C
C Main-field models are to degree and order 10 (ie 120 coefficients)
C for 1900.0-1995.0 and to 13 (ie 195 coefficients) for 2000.0 onwards.
C The predictive secular-variation model is to degree and order 8 (ie 80
C coefficients).
C
C Options include values at different locations at different
C times (spot), values at same location at one year intervals
C (time series), grid of values at one time (grid); geodetic or
C geocentric coordinates, latitude & longitude entered as decimal
C degrees or degrees & minutes (not in grid), choice of main field
C or secular variation or both (grid only).
C
c Adapted from 8th generation version to include new maximum degree for
c main-field models for 2000.0 and onwards and use WGS84 spheroid instead
c of International Astronomical Union 1966 spheroid as recommended by IAGA
c in July 2003. Reference radius remains as 6371.2 km - it is NOT the mean
c radius (= 6371.0 km) but 6371.2 km is what is used in determining the
c coefficients. Adaptation by Susan Macmillan, August 2003 (for
c 9th generation) and December 2004.
c 1995.0 coefficients as published in igrf9coeffs.xls and igrf10coeffs.xls
c used - (Kimmo Korhonen spotted 1 nT difference in 11 coefficients)
c Susan Macmillan July 2005
C
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
CHARACTER*11 TYPE
dimension wwlat(30000),wwlon(30000)
dimension yr0(30000),gr0(30000),gto0(30000),ga0(30000)
dimension iyr1(30000),jdy1(30000),min1(30000)
dimension speed(30000),head(30000),kdep(30000),kflg(30000)
dimension tmag(30000),nn(30000)
dtmn=1900.0
dtmx=2015.0
C
WRITE(6,*)
WRITE(6,*)'******************************************************'
WRITE(6,*)'* IGRF SYNTHESIS PROGRAM *'
WRITE(6,*)'* *'
WRITE(6,*)'* A program for the computation of geomagnetic *'
WRITE(6,*)'* field elements from the International Geomagnetic *'
WRITE(6,*)'* Reference Field (9th generation) as revised in *'
WRITE(6,*)'* December 2004 by the IAGA Working Group V-MOD. *'
WRITE(6,*)'* *'
WRITE(6,*)'* It is valid for dates from 1900.0 to 2010.0, *'
WRITE(6,*)'* values up to 2015.0 will be computed but with *'
WRITE(6,*)'* reduced accuracy. Values for dates before 1945.0 *'
WRITE(6,*)'* and after 2000.0 are non-definitive, otherwise the *'
WRITE(6,*)'* values are definitive. *'
WRITE(6,*)'* *'
WRITE(6,*)'* Susan Macmillan British Geological Survey *'
WRITE(6,*)'* Chair IAGA Working Group V-MOD *'
WRITE(6,*)'******************************************************'
WRITE(6,*)
c
open(20,file='H:\Myhome\gh90a.bs1',
+ status='old')
open(30,file='H:\Myhome\gh90a.bs2')
c if the year is URU, ik=1;others ik=0
ik=0
c
FACT = 180.0/3.141592654
NCOUNT = 0
C
itype=1
C
iopt=1
idm=2
c
i=0
1 continue
read(20,*,end=99) no,iyr,jdy,min,alt,aln,
+ spd0,hdg0,idp,mgt,mga,
+ gr,gto,ga,iflg
c write(6,*) no
i=i+1
nn(i)=no
iyr1(i)=iyr
jdy1(i)=jdy
min1(i)=min
wwlat(i)=dble(alt)
wwlon(i)=dble(aln)
speed(i)=spd0
head(i)=hdg0
kdep(i)=idp
tmag(i)=float(mgt)
gr0(i)=gr
gto0(i)=gto
ga0(i)=ga
kflg(i)=iflg
iyr=iyr+1900
yr0(i)=float(iyr)+float(jdy)/float(365+ik)
go to 1
99 close(20)
jnum=i
write(6,*) 'read end-30',jnum
k=0
C
50 ncount=1
ALT=0
C
do 235 is=1,jnum
xlt=wwlat(is)
xln=wwlon(is)
date=yr0(is)
IF (XLT.LT.-90.0.OR.XLT.GT.90.0) then
GO TO 202
end if
IF (XLN.LT.-360.0.OR.XLN.GT.360.0) then
GO TO 203
ENDIF
C
CLT = 90.0 - XLT
IF (CLT.LT.0.0.OR.CLT.GT.180.0) GO TO 204
IF (XLN.LE.-360.0.OR.XLN.GE.360.0) GO TO 205
C
CALL IGRF10SYN (0,DATE,ITYPE,ALT,CLT,XLN,X,Y,Z,F)
D = FACT*ATAN2(Y,X)
H = SQRT(X*X + Y*Y)
S = FACT*ATAN2(Z,H)
CALL DDECDM (D,IDEC,IDECM)
CALL DDECDM (S,INC,INCM)
C
CALL IGRF10SYN (1,DATE,ITYPE,ALT,CLT,XLN,DX,DY,DZ,F1)
DD = (60.0*FACT*(X*DY - Y*DX))/(H*H)
DH = (X*DX + Y*DY)/H
DS = (60.0*FACT*(H*DZ - Z*DH))/(F*F)
DF = (H*DH + Z*DZ)/F
C
amag=tmag(is)-f
c
write(6,*) is,iyr1(is),jdy1(is),min1(is),tmag(is),amag
c
write(30,207) nn(is),iyr1(is),jdy1(is),min1(is),xlt,xln,
+ speed(is),head(is),kdep(is),tmag(is),amag,
+ gr0(is),gto0(is),ga0(is),kflg(is)
207 format(i7,i3,i4,i5,f10.5,f11.5,f7.2,f6.1,i6,f9.1,f8.1,f7.1,
+ f10.1,f8.1,i2)
c
C
235 continue
C
60 CONTINUE
C
C
159 continue
close(30)
STOP
C
209 WRITE(6,972) DATE
972 FORMAT (' ***** Error *****'/' DATE =',F9.3,
1 ' - out of range')
STOP
C
210 WRITE(6,973) ALT,ITYPE
973 FORMAT (' ***** Error *****'/' A value of ALT =',F10.3,
1 ' is not allowed when ITYPE =',I2)
STOP
C
202 WRITE(6,966) XLT
966 FORMAT (' ***** Error *****'/' XLT =',F9.3,
1 ' - out of range')
STOP
C
203 WRITE(6,967) XLN
967 FORMAT (' ***** Error *****'/' XLN =',F10.3,
1 ' - out of range')
STOP
C
204 WRITE(6,968) LTD,LTM
968 FORMAT (' ***** Error *****'/' Latitude out of range',
1 ' - LTD =',I6,5X,'LTM =',I4)
STOP
C
205 WRITE(6,969) LND,LNM
969 FORMAT (' ***** Error *****'/' Longitude out of range',
1 ' - LND =',I8,5X,'LNM =',I4)
STOP
C
END
C
SUBROUTINE DMDDEC (I,M,X)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
DE = I
EM = M
IF (I.LT.0) EM = -EM
X = DE + EM/60.0
RETURN
END
C
SUBROUTINE DDECDM (X,I,M)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)
if(x.ge.0) then
sig=1.1
else
sig=-1.1
end if
c SIG = SIGN(1.1,X)
DR = ABS(X)
I = INT(DR)
T = I
M = NINT(60.*(DR - T))
IF (M.EQ.60) THEN
M = 0
I = I + 1
ENDIF
ISIG = INT(SIG)
IF (I.NE.0) THEN
I = I * ISIG
ELSE
IF (M.NE.0) M = M * ISIG
ENDIF
RETURN
END
c
subroutine igrf10syn (isv,date,itype,alt,colat,elong,x,y,z,f)
c
c This is a synthesis routine for the 10th generation IGRF as agreed
c in December 2004 by IAGA Working Group V-MOD. It is valid 1900.0 to
c 2010.0 inclusive. Values for dates from 1945.0 to 2000.0 inclusive are
c definitve, otherwise they are non-definitive.
c INPUT
c isv = 0 if main-field values are required
c isv = 1 if secular variation values are required
c date = year A.D. Must be greater than or equal to 1900.0 and
c less than or equal to 2015.0. Warning message is given
c for dates greater than 2010.0. Must be double precision.
c itype = 1 if geodetic (spheroid)
c itype = 2 if geocentric (sphere)
c alt = height in km above sea level if itype = 1
c = distance from centre of Earth in km if itype = 2 (>3485 km)
c colat = colatitude (0-180)
c elong = east-longitude (0-360)
c alt, colat and elong must be double precision.
c OUTPUT
c x = north component (nT) if isv = 0, nT/year if isv = 1
c y = east component (nT) if isv = 0, nT/year if isv = 1
c z = vertical component (nT) if isv = 0, nT/year if isv = 1
c f = total intensity (nT) if isv = 0, rubbish if isv = 1
c
c To get the other geomagnetic elements (D, I, H and secular
c variations dD, dH, dI and dF) use routines ptoc and ptocsv.
c
c Adapted from 8th generation version to include new maximum degree for
c main-field models for 2000.0 and onwards and use WGS84 spheroid instead
c of International Astronomical Union 1966 spheroid as recommended by IAGA
c in July 2003. Reference radius remains as 6371.2 km - it is NOT the mean
c radius (= 6371.0 km) but 6371.2 km is what is used in determining the
c coefficients. Adaptation by Susan Macmillan, August 2003 (for
c 9th generation) and December 2004.
c 1995.0 coefficients as published in igrf9coeffs.xls and igrf10coeffs.xls
c used - (Kimmo Korhonen spotted 1 nT difference in 11 coefficients)
c Susan Macmillan July 2005
c
implicit double precision (a-h,o-z)
dimension gh(3060),g0(120),g1(120),g2(120),g3(120),g4(120),
1 g5(120),g6(120),g7(120),g8(120),g9(120),ga(120),
2 gb(120),gc(120),gd(120),ge0(120),gf(120),gg(120),
3 gi(120),gj(120),gk(195),gl(195),gm(195),gp(195),
4 p(105),q(105),cl(13),sl(13)
equivalence (g0,gh( 1)),(g1,gh(121)),(g2,gh(241)),(g3,gh(361)),
1 (g4,gh(481)),(g5,gh(601)),(g6,gh(721)),(g7,gh(841)),
2 (g8,gh(961)),(g9,gh(1081)),(ga,gh(1201)),
3 (gb,gh(1321)),(gc,gh(1441)),(gd,gh(1561)),
4 (ge0,gh(1681)),(gf,gh(1801)),(gg,gh(1921)),
5 (gi,gh(2041)),(gj,gh(2161)),(gk,gh(2281)),
6 (gl,gh(2476)),(gm,gh(2671)),(gp,gh(2866))
c
open(10,file='H:\Myhome\programs\igrf10coeffs.data',status='old')
do 120 kkk=1,120
read(10,*) g0(kkk)
120 continue
do 121 kkk=1,120
read(10,*) g1(kkk)
121 continue
do 122 kkk=1,120
read(10,*) g2(kkk)
122 continue
do 123 kkk=1,120
read(10,*) g3(kkk)
123 continue
do 124 kkk=1,120
read(10,*) g4(kkk)
124 continue
do 125 kkk=1,120
read(10,*) g5(kkk)
125 continue
do 126 kkk=1,120
read(10,*) g6(kkk)
126 continue
do 127 kkk=1,120
read(10,*) g7(kkk)
127 continue
do 128 kkk=1,120
read(10,*) g8(kkk)
128 continue
do 129 kkk=1,120
read(10,*) g9(kkk)
129 continue
do 130 kkk=1,120
read(10,*) ga(kkk)
130 continue
do 131 kkk=1,120
read(10,*) gb(kkk)
131 continue
do 132 kkk=1,120
read(10,*) gc(kkk)
132 continue
do 133 kkk=1,120
read(10,*) gd(kkk)
133 continue
do 134 kkk=1,120
read(10,*) ge0(kkk)
134 continue
do 135 kkk=1,120
read(10,*) gf(kkk)
135 continue
do 136 kkk=1,120
read(10,*) gg(kkk)
136 continue
do 137 kkk=1,120
read(10,*) gi(kkk)
137 continue
do 138 kkk=1,120
read(10,*) gj(kkk)
138 continue
do 139 kkk=1,120
read(10,*) gk(kkk)
139 continue
do 140 kkk=1,195
read(10,*) gl(kkk)
140 continue
do 141 kkk=1,195
read(10,*) gm(kkk)
141 continue
do 142 kkk=1,80
read(10,*) gp(kkk)
142 continue
do 143 lll=121, 195
gk(lll)=0
143 continue
do 144 lll=81, 195
gp(lll)=0.0
144 continue
c
close(10)
c write(6,*) gp(1),g0(1)
c
c
c set initial values
c
x = 0.0
y = 0.0
z = 0.0
if (date.lt.1900.0.or.date.gt.2015.0) go to 11
if (date.gt.2010.0) write (6,960) date
960 format (/' This version of the IGRF is intended for use up',
1 ' to 2010.0.'/' values for',f9.3,' will be computed',
2 ' but may be of reduced accuracy'/)
if (date.ge.2005.0) go to 1
t = 0.2*(date - 1900.0)
ll = t
one = ll
t = t - one
if (date.lt.1995.0) then
nmx = 10
nc = nmx*(nmx+2)
ll = nc*ll
kmx = (nmx+1)*(nmx+2)/2
else
nmx = 13
nc = nmx*(nmx+2)
ll = 0.2*(date - 1995.0)
ll = 120*19 + nc*ll
kmx = (nmx+1)*(nmx+2)/2
endif
tc = 1.0 - t
if (isv.eq.1) then
tc = -0.2
t = 0.2
end if
go to 2
c
1 t = date - 2005.0
tc = 1.0
if (isv.eq.1) then
t = 1.0
tc = 0.0
end if
ll = 2670
nmx = 13
nc = nmx*(nmx+2)
kmx = (nmx+1)*(nmx+2)/2
2 r = alt
one = colat*0.017453292
ct = cos(one)
st = sin(one)
one = elong*0.017453292
cl(1) = cos(one)
sl(1) = sin(one)
cd = 1.0
sd = 0.0
l = 1
m = 1
n = 0
if (itype.eq.2) go to 3
c
c conversion from geodetic to geocentric coordinates
c (using the WGS84 spheroid)
c
a2 = 40680631.6
b2 = 40408296.0
one = a2*st*st
two = b2*ct*ct
three = one + two
rho = sqrt(three)
r = sqrt(alt*(alt + 2.0*rho) + (a2*one + b2*two)/three)
cd = (alt + rho)/r
sd = (a2 - b2)/rho*ct*st/r
one = ct
ct = ct*cd - st*sd
st = st*cd + one*sd
c
3 ratio = 6371.2/r
rr = ratio*ratio
c
c computation of Schmidt quasi-normal coefficients p and x(=q)
c
p(1) = 1.0
p(3) = st
q(1) = 0.0
q(3) = ct
do 10 k=2,kmx
if (n.ge.m) go to 4
m = 0
n = n + 1
rr = rr*ratio
fn = n
gn = n - 1
4 fm = m
if (m.ne.n) go to 5
if (k.eq.3) go to 6
one = sqrt(1.0 - 0.5/fm)
j = k - n - 1
p(k) = one*st*p(j)
q(k) = one*(st*q(j) + ct*p(j))
cl(m) = cl(m-1)*cl(1) - sl(m-1)*sl(1)
sl(m) = sl(m-1)*cl(1) + cl(m-1)*sl(1)
go to 6
5 gmm = m*m
one = sqrt(fn*fn - gmm)
two = sqrt(gn*gn - gmm)/one
three = (fn + gn)/one
i = k - n
j = i - n + 1
p(k) = three*ct*p(i) - two*p(j)
q(k) = three*(ct*q(i) - st*p(i)) - two*q(j)
c
c synthesis of x, y and z in geocentric coordinates
c
6 lm = ll + l
one = (tc*gh(lm) + t*gh(lm+nc))*rr
if (m.eq.0) go to 9
two = (tc*gh(lm+1) + t*gh(lm+nc+1))*rr
three = one*cl(m) + two*sl(m)
x = x + three*q(k)
z = z - (fn + 1.0)*three*p(k)
if (st.eq.0.0) go to 7
y = y + (one*sl(m) - two*cl(m))*fm*p(k)/st
go to 8
7 y = y + (one*sl(m) - two*cl(m))*q(k)*ct
8 l = l + 2
go to 10
9 x = x + one*q(k)
z = z - (fn + 1.0)*one*p(k)
l = l + 1
10 m = m + 1
c
c conversion to coordinate system specified by itype
c
one = x
x = x*cd + z*sd
z = z*cd - one*sd
f = sqrt(x*x + y*y + z*z)
c
return
c
c error return if date out of bounds
c
11 f = 1.0d8
write (6,961) date
961 format (/' This subroutine will not work with a date of',
1 f9.3,'. Date must be in the range 1900.0.ge.date',
2 '.le.2015.0. On return f = 1.0d8., x = y = z = 0.')
return
end