c
c The small main program below is an example of how to compute field
c components with T89D_SP.
c See GEOPACK-2008.DOC for a sample field line tracing program.
c
dimension parmod(10)
1 print *, ' enter x,y,z,ps,iopt'
read*, x,y,z,ps,iopt
call T89D_SP (iopt,parmod,ps,x,y,z,bx,by,bz)
print *, bx,by,bz
goto 1
end
C=======================================================================
C
SUBROUTINE T89D_SP (IOPT,PARMOD,PS,X,Y,Z,BX,BY,BZ)
C (single precision version)
C
C COMPUTES GSM COMPONENTS OF THE MAGNETIC FIELD PRODUCED BY EXTRA-
C TERRESTRIAL CURRENT SYSTEMS IN THE GEOMAGNETOSPHERE. THE MODEL IS
C VALID UP TO GEOCENTRIC DISTANCES OF 70 RE AND IS BASED ON THE MER-
C GED IMP-A,C,D,E,F,G,H,I,J (1966-1974), HEOS-1 AND -2 (1969-1974),
C AND ISEE-1 AND -2 SPACECRAFT DATA SET.
C
c NOTE OF NOV 12, 2014:
C
C THIS IS THIRD UPGRADE OF THE OLD T89 MODEL, IN WHICH ALL COEFFICIENTS
C (HENCE, THE MODEL OUTPUT) REMAINED IDENTICAL TO THOSE OF T89C.
C ALL MODIFICATIONS WERE MADE WITH RESPECT TO THE CODE STRUCTURE, WHICH
C MADE IT POSSIBLE TO SIGNIFICANTLY BOOST THE SPEED.
C--------------------------------------------------------------------------
C ----------- OLD COMMENTS OF 1992: ------------------------------------
C
C THIS IS A MODIFIED VERSION (T89c), WHICH REPLACED THE ORIGINAL ONE
C IN 1992 AND DIFFERS FROM IT IN THE FOLLOWING:
C
C (1) ISEE-1,2 DATA WERE ADDED TO THE ORIGINAL IMP-HEOS DATASET
C (2) TWO TERMS WERE ADDED TO THE ORIGINAL TAIL FIELD MODES, ALLOWING
C A MODULATION OF THE CURRENT BY THE GEODIPOLE TILT ANGLE
C
C REFERENCE FOR THE ORIGINAL MODEL: N.A. TSYGANENKO, A MAGNETOSPHERIC MAGNETIC
C FIELD MODEL WITH A WARPED TAIL CURRENT SHEET: PLANET.SPACE SCI., V.37,
C PP.5-20, 1989.
C
C----INPUT PARAMETERS: IOPT - SPECIFIES THE GROUND DISTURBANCE LEVEL:
C
C IOPT= 1 2 3 4 5 6 7
C CORRESPOND TO:
C KP= 0,0+ 1-,1,1+ 2-,2,2+ 3-,3,3+ 4-,4,4+ 5-,5,5+ > =6-
C
C PS - GEODIPOLE TILT ANGLE IN RADIANS
C X, Y, Z - GSM COORDINATES OF THE POINT IN EARTH RADII
C
C----OUTPUT PARAMETERS: BX,BY,BZ - GSM COMPONENTS OF THE MODEL MAGNETIC
C FIELD IN NANOTESLAS
c
c THE PARAMETER PARMOD(10) IS A DUMMY ARRAY. IT IS NOT USED IN THIS
C SUBROUTINE AND IS PROVIDED JUST FOR MAKING IT COMPATIBLE WITH THE
C NEW VERSION (4/16/96) OF THE GEOPACK SOFTWARE.
C
C THIS RELEASE OF T89C/D IS DATED FEB 12, 1996; UPDATED APR 29, 2013
C--------------------------------------------------------------------------
C
C
C AUTHOR: NIKOLAI A. TSYGANENKO
C HSTX CORP./NASA GSFC (1992-2007)
C SPB STATE UNIVERSITY (2007-present)
C
DIMENSION PARAM(30,7),A(30),PARMOD(10)
DATA A02,XLW2,YN,RPI,RT/25.,170.,30.,0.31830989D0,30./
DATA XD,XLD2/0.,40./
C
C The last 2 quantities define variation of tail sheet thickness along X
C
DATA SXC,XLWC2/4.,50./
C
C The two quantities belong to the function WC which confines tail closure
c current in X- and Y- direction
C
DATA DXL/20./
DATA PARAM/-116.53,-10719.,42.375,59.753,-11363.,
* 1.7844,30.268,-.35372E-01,-0.66832E-01,0.16456E-01,-1.3024,
* 0.16529E-02,0.20293E-02,20.289,-0.25203E-01,224.91,-9234.8,
* 22.788,7.8813,1.8362,-0.27228,8.8184,2.8714,14.468,
* 32.177,0.01,0.0,7.0459,4.0,20.0,
* -55.553,-13198.,
* 60.647,61.072,-16064.,2.2534,34.407,-0.38887E-01,
* -0.94571E-01,0.27154E-01,-1.3901,0.13460E-02,0.13238E-02,
* 23.005,-0.30565E-01,55.047,-3875.7,20.178,7.9693,
* 1.4575,0.89471,9.4039,3.5215,14.474,36.555,0.01,
* 0.0,7.0787,4.0,20.,
*-101.34,-13480.,111.35,12.386,
* -24699.,2.6459,38.948,-0.34080E-01,-0.12404,0.29702E-01,
* -1.4052,0.12103E-02,0.16381E-02,24.49,-0.37705E-01,-298.32,
* 4400.9,18.692,7.9064,1.3047,2.4541,9.7012,7.1624,
* 14.288,33.822,0.01,0.0,6.7442,4.0,20.0,
*-181.69,
* -12320.,173.79,-96.664,-39051.,3.2633,44.968,
* -0.46377E-01,-0.16686,0.048298,-1.5473,0.10277E-02,
* 0.31632E-02,27.341,-0.50655E-01,-514.10,12482.,16.257,
* 8.5834,1.0194,3.6148,8.6042,5.5057,13.778,32.373,
* 0.01,0.0,7.3195,4.0,20.0,
*-436.54,-9001.0,323.66,
* -410.08,-50340.,3.9932,58.524,-0.38519D-01,-0.26822,
* 0.74528E-01,-1.4268,-0.10985E-02,0.96613E-02,27.557,
* -0.56522E-01,-867.03,20652.,14.101,8.3501,0.72996,
* 3.8149,9.2908,6.4674,13.729,28.353,0.01,0.0,
* 7.4237,4.0,20.0,-707.77,-4471.9,432.81,-435.51,
* -60400.,4.6229,68.178,-0.88245E-01,-0.21002,0.11846,
* -2.6711,0.22305E-02,.10910E-01,27.547,-0.54080E-01,-424.23,
* 1100.2,13.954,7.5337,0.89714,3.7813,8.2945,5.174,
* 14.213,25.237,0.01,0.0,7.0037,4.0,20.0,-1190.4,
* 2749.9,742.56,-1110.3,-77193.,7.6727,102.05,
*-0.96015E-01,-0.74507,0.11214,-1.3614,0.15157E-02,0.22283E-01,
*23.164,-0.74146E-01,-2219.1,48253.,12.714,7.6777,.57138,
* 2.9633,9.3909,9.7263,11.123,21.558,0.01,0.0,
* 4.4518,4.0,20.0/
DATA IOP/10/
SAVE
C
IF (IOP.NE.IOPT) THEN
C
IOP=IOPT
DO 1 I=1,30
1 A(I)=PARAM(I,IOPT)
C
DYC=A(30)
DYC2=DYC**2
DX=A(18)
HA02=0.5*A02
RDX2M=-1./DX**2
RDX2=-RDX2M
RDYC2=1./DYC2
HLWC2M=-0.5*XLWC2
DRDYC2=-2.*RDYC2
DRDYC3=2.*RDYC2*SQRT(RDYC2)
HXLW2M=-0.5*XLW2
ADR=A(19)
D0=A(20)
DD=A(21)
RC=A(22)
G=A(23)
AT=A(24)
DT=D0
DEL=A(26)
P=A(25)
Q=A(27)
SX=A(28)
GAM=A(29)
HXLD2M=-0.5*XLD2
ADSL=0.
XGHS=0.
H=0.
HS=0.
GAMH=0.
W1=-0.5/DX
DBLDEL=2.*DEL
W2=W1*2.
W4=-1./3.
W3=W4/DX
W5=-0.5
W6=-3.
AK1=A(1)
AK2=A(2)
AK3=A(3)
AK4=A(4)
AK5=A(5)
AK6=A(6)
AK7=A(7)
AK8=A(8)
AK9=A(9)
AK10=A(10)
AK11=A(11)
AK12=A(12)
AK13=A(13)
AK14=A(14)
AK15=A(15)
AK16=A(16)
AK17=A(17)
SXA=0.
SYA=0.
SZA=0.
AK610=AK6*W1+AK10*W5
AK711=AK7*W2-AK11
AK812=AK8*W2+AK12*W6
AK913=AK9*W3+AK13*W4
RDXL=1./DXL
HRDXL=0.5*RDXL
A6H=AK6*0.5
A9T=AK9/3.
YNP=RPI/YN*0.5
YND=2.*YN
C
ENDIF
C
SPS = SIN(PS)
CPS = COS(PS)
C
X2=X*X
Y2=Y*Y
Z2=Z*Z
TPS=SPS/CPS
HTP=TPS*0.5
GSP=G*SPS
XSM=X*CPS-Z*SPS
ZSM=X*SPS+Z*CPS
C
C CALCULATE THE FUNCTION ZS DEFINING THE SHAPE OF THE TAIL CURRENT SHEET
C AND ITS SPATIAL DERIVATIVES:
C
XRC=XSM+RC
XRC16=XRC**2+16.
SXRC=SQRT(XRC16)
Y4=Y2*Y2
Y410=Y4+1.E4
SY4=SPS/Y410
GSY4=G*SY4
ZS1=HTP*(XRC-SXRC)
DZSX=-ZS1/SXRC
ZS=ZS1-GSY4*Y4
D2ZSGY=-SY4/Y410*4.E4*Y2*Y
DZSY=G*D2ZSGY
C
C CALCULATE THE COMPONENTS OF THE RING CURRENT CONTRIBUTION:
C
XSM2=XSM**2
DSQT=SQRT(XSM2+A02)
FA0=0.5*(1.+XSM/DSQT)
DDR=D0+DD*FA0
DFA0=HA02/DSQT**3
ZR=ZSM-ZS
TR=SQRT(ZR**2+DDR**2)
RTR=1./TR
RO2=XSM2+Y2
ADRT=ADR+TR
ADRT2=ADRT**2
FK=1./(ADRT2+RO2)
DSFC=SQRT(FK)
FC=FK**2*DSFC
FACXY=3.0*ADRT*FC*RTR
XZR=XSM*ZR
YZR=Y*ZR
DBXDP=FACXY*XZR
DER25=FACXY*YZR
XZYZ=XSM*DZSX+Y*DZSY
FAQ=ZR*XZYZ-DDR*DD*DFA0*XSM
DBZDP=FC*(2.*ADRT2-RO2)+FACXY*FAQ
DER15=DBXDP*CPS+DBZDP*SPS
DER35=DBZDP*CPS-DBXDP*SPS
C
C CALCULATE THE TAIL CURRENT SHEET CONTRIBUTION:
C
DELY2=DEL*Y2
D=DT+DELY2
IF (ABS(GAM).LT.1.E-6) GOTO 8
XXD=XSM-XD
RQD=1./(XXD**2+XLD2)
RQDS=SQRT(RQD)
H=0.5*(1.+XXD*RQDS)
HS=-HXLD2M*RQD*RQDS
GAMH=GAM*H
D=D+GAMH
XGHS=XSM*GAM*HS
ADSL=-D*XGHS
8 D2=D**2
T=SQRT(ZR**2+D2)
XSMX=XSM-SX
RDSQ2=1./(XSMX**2+XLW2)
RDSQ=SQRT(RDSQ2)
V=0.5*(1.-XSMX*RDSQ)
DVX=HXLW2M*RDSQ*RDSQ2
OM=SQRT(SQRT(XSM2+16.)-XSM)
OMS=-OM/(OM*OM+XSM)*0.5
RDY=1./(P+Q*OM)
OMSV=OMS*V
RDY2=RDY**2
FY=1./(1.+Y2*RDY2)
W=V*FY
YFY1=2.*FY*Y2*RDY2
FYPR=YFY1*RDY
FYDY=FYPR*FY
DWX=DVX*FY+FYDY*Q*OMSV
YDWY=-V*YFY1*FY
DDY=DBLDEL*Y
ATT=AT+T
S1=SQRT(ATT**2+RO2)
F5=1./S1
F7=1./(S1+ATT)
F1=F5*F7
F3=F5**3
F9=ATT*F3
FS=ZR*XZYZ-D*Y*DDY+ADSL
XDWX=XSM*DWX+YDWY
RTT=1./T
WT=W*RTT
BRRZ1=WT*F1
BRRZ2=WT*F3
DBXC1=BRRZ1*XZR
DBXC2=BRRZ2*XZR
TLT2=PS**2
DER21=BRRZ1*YZR
DER22=BRRZ2*YZR
DER216=DER21*TLT2
DER217=DER22*TLT2
WTFS=WT*FS
DBZC1=W*F5+XDWX*F7+WTFS*F1
DBZC2=W*F9+XDWX*F1+WTFS*F3
DER11=DBXC1*CPS+DBZC1*SPS
DER12=DBXC2*CPS+DBZC2*SPS
DER31=DBZC1*CPS-DBXC1*SPS
DER32=DBZC2*CPS-DBXC2*SPS
DER116=DER11*TLT2
DER117=DER12*TLT2
DER316=DER31*TLT2
DER317=DER32*TLT2
C
C CALCULATE CONTRIBUTION FROM THE CLOSURE CURRENTS
C
ZPL=Z+RT
ZMN=Z-RT
ROGSM2=X2+Y2
SPL=SQRT(ZPL**2+ROGSM2)
SMN=SQRT(ZMN**2+ROGSM2)
XSXC=X-SXC
RQC2=1./(XSXC**2+XLWC2)
RQC=SQRT(RQC2)
FYC=1./(1.+Y2*RDYC2)
WC=0.5*(1.-XSXC*RQC)*FYC
DWCX=HLWC2M*RQC2*RQC*FYC
DWCY=DRDYC2*WC*FYC*Y
SZRP=1./(SPL+ZPL)
SZRM=1./(SMN-ZMN)
XYWC=X*DWCX+Y*DWCY
WCSP=WC/SPL
WCSM=WC/SMN
FXYP=WCSP*SZRP
FXYM=WCSM*SZRM
FXPL=X*FXYP
FXMN=-X*FXYM
FYPL=Y*FXYP
FYMN=-Y*FXYM
FZPL=WCSP+XYWC*SZRP
FZMN=WCSM+XYWC*SZRM
DER13=FXPL+FXMN
DER14=(FXPL-FXMN)*SPS
DER23=FYPL+FYMN
DER24=(FYPL-FYMN)*SPS
DER33=FZPL+FZMN
DER34=(FZPL-FZMN)*SPS
C
C NOW CALCULATE CONTRIBUTION FROM CHAPMAN-FERRARO SOURCES + ALL OTHER
C
EX=EXP(X/DX)
EC=EX*CPS
ES=EX*SPS
ECZ=EC*Z
ESZ=ES*Z
ESZY2=ESZ*Y2
ESZZ2=ESZ*Z2
ECZ2=ECZ*Z
ESY=ES*Y
C
C FINALLY, CALCULATE NET EXTERNAL MAGNETIC FIELD COMPONENTS,
C BUT FIRST OF ALL THOSE FOR C.-F. FIELD:
C
SX1=AK6*ECZ+AK7*ES+AK8*ESY*Y+AK9*ESZ*Z
SY1=AK10*ECZ*Y+AK11*ESY+AK12*ESY*Y2+AK13*ESY*Z2
SZ1=AK14*EC+AK15*EC*Y2+AK610*ECZ2+AK711*ESZ+AK812
* *ESZY2+AK913*ESZZ2
BXCL=AK3*DER13+AK4*DER14
BYCL=AK3*DER23+AK4*DER24
BZCL=AK3*DER33+AK4*DER34
BXT=AK1*DER11+AK2*DER12+BXCL +AK16*DER116+AK17*DER117
BYT=AK1*DER21+AK2*DER22+BYCL +AK16*DER216+AK17*DER217
BZT=AK1*DER31+AK2*DER32+BZCL +AK16*DER316+AK17*DER317
BX=BXT+AK5*DER15+SX1+SXA
BY=BYT+AK5*DER25+SY1+SYA
BZ=BZT+AK5*DER35+SZ1+SZA
C
RETURN
END
c
c=======================================================================