!***************************************************************************************************
! Copyright 2000 D. Boscher, 2004 S. Bourdarie
!
! This file is part of IRBEM-LIB.
!
! IRBEM-LIB is free software: you can redistribute it and/or modify
! it under the terms of the GNU Lesser General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! IRBEM-LIB is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public License
! along with IRBEM-LIB. If not, see .
!
SUBROUTINE CHAMP(xGEO,BxGEO,Bl,Ifail)
C
IMPLICIT NONE
C
INTEGER*4 I,IOPT,Ifail,kint
INTEGER*4 Activ,k_ext,k_l
REAL*8 xGEO(3)
REAL*8 xSM(3)
REAL*8 BxSM(3),BIMF(3)
REAL*8 Bxext(3)
REAL*8 Bxint(3),den,vel,dst
REAL*8 BxGEO(3),Bl,density,speed,dst_nt,fkp,fkp_osta
REAL*8 Pdyn_nPa,ByIMF_nt,BzIMF_nt,G1_tsy01,G2_tsy01,G3_tsy01
REAL*8 W1_tsy04,W2_tsy04,W3_tsy04,W4_tsy04,W5_tsy04,W6_tsy04
REAL*8 Al,Al_ind,BxIMF_nt
REAL*8 PARMOD(10),tilt,sn,Pdyn,BzIMF
C
REAL*8 Bo,xc,yc,zc,ct,st,cp,sp
C
COMMON /dipigrf/Bo,xc,yc,zc,ct,st,cp,sp
COMMON /magmod/k_ext,k_l,kint
COMMON /drivers/density,speed,dst_nt,Pdyn_nPa,BxIMF_nt,ByIMF_nt,
& BzIMF_nt,G1_tsy01,G2_tsy01,fkp,G3_tsy01,W1_tsy04,W2_tsy04,
& W3_tsy04,W4_tsy04,W5_tsy04,W6_tsy04,Al
COMMON /index/activ
COMMON /dip_ang/tilt
C
Ifail=0
IF(kint.EQ.0)THEN
CALL IGRF(xGEO(1),xGEO(2),xGEO(3),Bxint(1),Bxint(2),Bxint(3))
ENDIF
IF(kint.EQ.1) THEN
CALL DTD(xGEO(1),xGEO(2),xGEO(3),Bxint(1),Bxint(2),Bxint(3))
ENDIF
IF(kint.GE.2 .and. kint .le. 3) THEN
CALL get_intfield(xGEO(1),xGEO(2),xGEO(3),
& Bxint(1),Bxint(2),Bxint(3))
ENDIF
IF(kint .eq. 4) THEN
CALL myOwnMagField(xGEO,Bxint)
ENDIF
IF(kint .eq. 5) THEN
CALL centered_dipole(xGEO(1),xGEO(2),xGEO(3),Bxint(1),Bxint(2)
& ,Bxint(3))
ENDIF
C
Bxext(1) = 0.D0
Bxext(2) = 0.D0
Bxext(3) = 0.D0
c
c Case for Mead dynamic mag model (use of SM coordinates)
c inputs is magnetic activity index deduced from Kp
c This field model is only valid for rSM<15 Re
if (k_ext .eq. 1) then
IOPT=Activ
IF (xGEO(1)*xGEO(1)+xGEO(2)*xGEO(2)+xGEO(3)*xGEO(3)
& .GT.289.D0) THEN
Ifail=-1
RETURN
ENDIF
CALL GEO_SM(xGEO,xSM)
CALL MEAD(xSM(1),xSM(2),xSM(3),IOPT,BxSM(1),BxSM(2),BxSM(3))
CALL SM_GEO(BxSM,Bxext)
endif
c
c Case for Tsyganenko 87 short dynamic mag model (use of GSM coordinates)
c inputs is magnetic activity index deduced from Kp
c This field model is only valid for rGEO<30 Re
if (k_ext .eq. 2) then
IOPT=Activ
IF (xGEO(1)*xGEO(1)+xGEO(2)*xGEO(2)+xGEO(3)*xGEO(3)
& .GT.900.D0) THEN
Ifail=-1
RETURN
ENDIF
CALL GEO_GSM(xGEO,xSM)
CALL TSY87S(IOPT,xSM(1),xSM(2),xSM(3),BxSM(1),BxSM(2),BxSM(3))
CALL GSM_GEO(BxSM,Bxext)
endif
c
c Case for Tsyganenko 87 long dynamic mag model (use of GSM coordinates)
c inputs is magnetic activity index deduced from Kp
c This field model is only valid for rGEO<70 Re
if (k_ext .eq. 3) then
IOPT=Activ
IF (xGEO(1)*xGEO(1)+xGEO(2)*xGEO(2)+xGEO(3)*xGEO(3)
& .GT.4900.D0) THEN
Ifail=-1
RETURN
ENDIF
CALL GEO_GSM(xGEO,xSM)
CALL TSY87l(IOPT,xSM(1),xSM(2),xSM(3),BxSM(1),BxSM(2),BxSM(3))
CALL GSM_GEO(BxSM,Bxext)
endif
c
c Case for Tsyganenko 89 dynamic mag model (use of GSM coordinates)
c inputs is magnetic activity index deduced from Kp
c This field model is only valid for rGEO<70 Re
if (k_ext .eq. 4) then
IOPT=Activ
IF (xGEO(1)*xGEO(1)+xGEO(2)*xGEO(2)+xGEO(3)*xGEO(3)
& .GT.4900.D0) THEN
Ifail=-1
RETURN
ENDIF
CALL GEO_GSM(xGEO,xSM)
CALL T89C(IOPT,xSM(1),xSM(2),xSM(3),BxSM(1),BxSM(2),BxSM(3))
CALL GSM_GEO(BxSM,Bxext)
endif
c Case for Olson-Pfitzer quiet mag model (use of SM coordinates)
c no parameters as input
c This field model is only valid for rGEO<15 Re
if (k_ext .eq. 5) then
IF (xGEO(1)*xGEO(1)+xGEO(2)*xGEO(2)+xGEO(3)*xGEO(3)
& .GT.225.D0) THEN
Ifail=-1
RETURN
ENDIF
CALL GEO_SM(xGEO,xSM)
CALL BXYZMU(xSM(1),xSM(2),xSM(3),BxSM(1),BxSM(2),BxSM(3))
CALL SM_GEO(BxSM,Bxext)
endif
c Case for Olson-Pfitzer dynamic mag model (use of SM coordinates)
c inputs are SW density, speed, and Dst index
c This field model is only valid for rGEO<60 Re
if (k_ext .eq. 6) then
IF (xGEO(1)*xGEO(1)+xGEO(2)*xGEO(2)+xGEO(3)*xGEO(3)
& .GT.3600.D0) THEN
Ifail=-1
RETURN
ENDIF
CALL GEO_SM(xGEO,xSM)
den=density
vel=speed
dst=dst_nt
CALL BDYN(den,vel,dst,xSM(1),xSM(2),xSM(3),BxSM(1),
& BxSM(2),BxSM(3))
CALL SM_GEO(BxSM,Bxext)
endif
c
c Case for rescent Tsyganenko dynamic mag model 96 (use of GSM coordinates)
c inputs is SW pressure, (nPa), DST (nT), ByIMF and BzIMF (nT)
c
if (k_ext .eq. 7) then
IF (xGEO(1)*xGEO(1)+xGEO(2)*xGEO(2)+xGEO(3)*xGEO(3)
& .GT.1600.D0) THEN
Ifail=-1
RETURN
ENDIF
PARMOD(1)=Pdyn_nPa
PARMOD(2)=dst_nt
PARMOD(3)=ByIMF_nt
PARMOD(4)=BzIMF_nt
CALL GEO_GSM(xGEO,xSM)
CALL T96_01(PARMOD,xSM(1),xSM(2),xSM(3),BxSM(1),BxSM(2),
& BxSM(3))
CALL GSM_GEO(BxSM,Bxext)
endif
c
c Case for Ostapenko dynamic mag model 97 (use of SM coordinates)
c inputs is SW pressure, (nPa), DST (nT), f(Kp) and BzIMF (nT)
if (k_ext .eq. 8) then
dst=dst_nt
Pdyn=Pdyn_nPa
BzIMF=BzIMF_nt
fkp_osta=fkp
sn =SIN(tilt*4.D0*ATAN(1.D0)/180.d0)
CALL GEO_SM(xGEO,xSM)
CALL set_a(dst,Pdyn,fkp_osta,BzIMF,sn)
CALL BOM97(xSM,BxSM)
CALL SM_GEO(BxSM,Bxext)
endif
c
c Case for Tsyganenko dynamic mag model 01 (use of GSM coordinates)
c inputs is SW pressure, (nPa), DST (nT), ByIMF and BzIMF (nT), G1,G2
c
if (k_ext .eq. 9) then
PARMOD(1)=Pdyn_nPa
PARMOD(2)=dst_nt
PARMOD(3)=ByIMF_nt
PARMOD(4)=BzIMF_nt
PARMOD(5)=G1_tsy01
PARMOD(6)=G2_tsy01
CALL GEO_GSM(xGEO,xSM)
IF (xSM(1).LT.-15.D0) THEN
Ifail=-1
RETURN
ENDIF
CALL T01_01(PARMOD,xSM(1),xSM(2),xSM(3),BxSM(1),BxSM(2),
& BxSM(3))
CALL GSM_GEO(BxSM,Bxext)
endif
c
c Case for Tsyganenko dynamic mag model 01 storm (use of GSM coordinates)
c inputs is SW pressure, (nPa), DST (nT), ByIMF and BzIMF (nT), G1,G2
c
if (k_ext .eq. 10) then
PARMOD(1)=Pdyn_nPa
PARMOD(2)=dst_nt
PARMOD(3)=ByIMF_nt
PARMOD(4)=BzIMF_nt
PARMOD(5)=G2_tsy01
PARMOD(6)=G3_tsy01
CALL GEO_GSM(xGEO,xSM)
IF (xSM(1).LT.-15.D0) THEN
Ifail=-1
RETURN
ENDIF
CALL T01_01(PARMOD,xSM(1),xSM(2),xSM(3),BxSM(1),BxSM(2),
& BxSM(3))
CALL GSM_GEO(BxSM,Bxext)
endif
c
c Case for Tsyganenko dynamic mag model 04 storm (use of GSM coordinates)
c inputs is SW pressure, (nPa), DST (nT), ByIMF and BzIMF (nT), W1-6
c
if (k_ext .eq. 11) then
PARMOD(1)=Pdyn_nPa
PARMOD(2)=dst_nt
PARMOD(3)=ByIMF_nt
PARMOD(4)=BzIMF_nt
PARMOD(5)=W1_tsy04
PARMOD(6)=W2_tsy04
PARMOD(7)=W3_tsy04
PARMOD(8)=W4_tsy04
PARMOD(9)=W5_tsy04
PARMOD(10)=W6_tsy04
CALL GEO_GSM(xGEO,xSM)
IF (xSM(1).LT.-15.D0) THEN
Ifail=-1
RETURN
ENDIF
CALL T04_s(PARMOD,xSM(1),xSM(2),xSM(3),BxSM(1),BxSM(2),
& BxSM(3))
CALL GSM_GEO(BxSM,Bxext)
endif
c
c Case for Alexeev dynamic mag model 2000 (use of GSM coordinates)
c inputs is SW pressure, (nPa), DST (nT), ByIMF and BzIMF (nT), G1,G2
c
if (k_ext .eq. 12) then
IF (xGEO(1)*xGEO(1)+xGEO(2)*xGEO(2)+xGEO(3)*xGEO(3)
& .GT.50.D0) THEN
Ifail=-1
RETURN
ENDIF
den=density
vel=speed
dst=dst_nt
BIMF(1)=BxIMF_nt
BIMF(2)=ByIMF_nt
BIMF(3)=BzIMF_nt
Al_ind=Al
CALL GEO_GSM(xGEO,xSM)
CALL a2000(den,vel,BIMF,
& dst,al_ind,xSM,BxSM,ifail)
if (Ifail .eq. 0) CALL GSM_GEO(BxSM,Bxext)
endif
c
C
c Case for Tsyganenko and Sitnov 07d mag model (use of GSM coordinates)
c input is Pdyn
c
c 2015 version, new bess calls in external file TS07j_2015.f (previously TS07j.f)
IF (k_ext .eq. 13) THEN
c M. Sitnov and G. Stephens recommend limiting the valid output to 20 Re and
c less based on the current (at the time) expansion:
IF (xGEO(1)*xGEO(1)+xGEO(2)*xGEO(2)+xGEO(3)*xGEO(3)
& .GT.400.D0) THEN
Ifail=-1
RETURN
ENDIF
CALL GEO_GSM(xGEO,xSM)
CALL TS07D_2015 (xSM(1),xSM(2),xSM(3),BxSM(1),BxSM(2),
& BxSM(3)) ! note that no solar wind input is necessary as it is
! included through the TS07D common block
CALL GSM_GEO(BxSM,Bxext)
endif
c 2017 version, new bess calls implemented directly in Tsy_and_Sit_Jul2017.f code
c Code is double precision by design
IF (k_ext .eq. 14) THEN
c M. Sitnov and G. Stephens recommend limiting the valid output to 20 Re and
c less based on the current (at the time) expansion:
IF (xGEO(1)*xGEO(1)+xGEO(2)*xGEO(2)+xGEO(3)*xGEO(3)
& .GT.400.D0) THEN
Ifail=-1
RETURN
ENDIF
CALL GEO_GSM(xGEO,xSM)
CALL TS07D_JULY_2017 (xSM(1),xSM(2),xSM(3),BxSM(1),BxSM(2),
& BxSM(3)) ! note that no solar wind input is necessary as it is
! included through the TS07D common block
CALL GSM_GEO(BxSM,Bxext)
endif
Bl = 0.D0
DO I = 1,3
!WRITE(6,*)Bxext(I),Bxint(I)
BxGEO(I) = Bxext(I) + Bxint(I)
Bl = Bl + BxGEO(I)*BxGEO(I)
ENDDO
C
Bl = SQRT(Bl)
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE INIT_DTD(year)
C
IMPLICIT NONE
C
INTEGER*4 ierr
REAL*8 year
REAL*8 g(66),h(66)
REAL*8 xc,yc,zc !Re
REAL*8 thet,phit !radian
REAL*8 ct,st,cp,sp
REAL*8 Bo
C
COMMON /dipigrf/Bo,xc,yc,zc,ct,st,cp,sp
COMMON /dgrf/g,h
C
call get_igrf_coeffs(year,g,h,ierr)
c
call get_terms(g,h,thet,phit,xc,yc,zc,Bo)
c write(6,*)Bo
c
CALL calcul_GH1
C
ct = COS(thet)
st = SIN(thet)
cp = COS(phit)
sp = SIN(phit)
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE INIT_GSM(iyr,iday,secs,psi)
C
IMPLICIT NONE
C
INTEGER*4 iyr,iday
REAL*8 secs
REAL*8 gst,Slong,Srasn,Sdec
REAL*8 xS,yS,zS !GEI
REAL*8 xDip,yDip,zDip !GEO
REAL*8 xD,yD,zD !GEI
REAL*8 xSD,ySD,zSD,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
REAL*8 xc,yc,zc
REAL*8 ct,st,cp,sp
REAL*8 cgst,sgst
REAL*8 Bo
REAL*8 norm
REAL*8 psi
C
COMMON /Soleil/xS,yS,zS,cgst,sgst
COMMON /dipigrf/Bo,xc,yc,zc,ct,st,cp,sp
COMMON /sundip/xD,yD,zD,xSD,ySD,zSD
& ,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
CALL SUN(iyr,iday,secs,gst,Slong,Srasn,Sdec)
C
xS = COS(Srasn)*COS(Sdec)
yS = SIN(Srasn)*COS(Sdec)
zS = SIN(Sdec)
C
cgst = COS(gst)
sgst = SIN(gst)
C
xDip = st*cp
yDip = st*sp
zDip = ct
C
xD = cgst*xDip - sgst*yDip
yD = sgst*xDip + cgst*yDip
zD = zDip
C
xSD = yD*zS - yS*zD
ySD = zD*xS - zS*xD
zSD = xD*yS - xS*yD
C
norm = SQRT(xSD*xSD + ySD*ySD + zSD*zSD)
C
xSD = xSD/norm
ySD = ySD/norm
zSD = zSD/norm
C
xSSD = yS*zSD - ySD*zS
ySSD = zS*xSD - zSD*xS
zSSD = xS*ySD - xSD*yS
C
norm = SQRT(xSSD*xSSD + ySSD*ySSD + zSSD*zSSD)
C
xSSD = xSSD/norm
ySSD = ySSD/norm
zSSD = zSSD/norm
C
xSDD = ySD*zD - yD*zSD
ySDD = zSD*xD - zD*xSD
zSDD = xSD*yD - xD*ySD
C
norm = SQRT(xSDD*xSDD + ySDD*ySDD + zSDD*zSDD)
C
xSDD = xSDD/norm
ySDD = ySDD/norm
zSDD = zSDD/norm
C
psi = ASIN(xD*xS + yD*yS + zD*zS)
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE
& GEO_GSM(xGEO,xGSM)
C
IMPLICIT NONE
C
REAL*8 xGEO(3)
REAL*8 xGSM(3)
REAL*8 xGEI,yGEI,zGEI
REAL*8 xS,yS,zS,cgst,sgst
REAL*8 xD,yD,zD !GEI
REAL*8 xSD,ySD,zSD,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
COMMON /Soleil/xS,yS,zS,cgst,sgst
COMMON /sundip/xD,yD,zD,xSD,ySD,zSD
& ,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
xGEI = cgst*xGEO(1) - sgst*xGEO(2)
yGEI = sgst*xGEO(1) + cgst*xGEO(2)
zGEI = xGEO(3)
C
xGSM(1) = xS*xGEI + yS*yGEI + zS*zGEI
xGSM(2) = xSD*xGEI + ySD*yGEI + zSD*zGEI
xGSM(3) = xSSD*xGEI + ySSD*yGEI + zSSD*zGEI
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE
& GSM_GEO(xGSM,xGEO)
C
IMPLICIT NONE
C
REAL*8 xGSM(3)
REAL*8 xGEO(3)
REAL*8 xGEI,yGEI,zGEI
REAL*8 xS,yS,zS,cgst,sgst
REAL*8 xD,yD,zD !GEI
REAL*8 xSD,ySD,zSD,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
COMMON /Soleil/xS,yS,zS,cgst,sgst
COMMON /sundip/xD,yD,zD,xSD,ySD,zSD
& ,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
xGEI = xS*xGSM(1) + xSD*xGSM(2) + xSSD*xGSM(3)
yGEI = yS*xGSM(1) + ySD*xGSM(2) + ySSD*xGSM(3)
zGEI = zS*xGSM(1) + zSD*xGSM(2) + zSSD*xGSM(3)
C
xGEO(1) = cgst*xGEI + sgst*yGEI
xGEO(2) = -sgst*xGEI + cgst*yGEI
xGEO(3) = zGEI
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE
& GEO_GEI(xGEO,xGEI)
C
IMPLICIT NONE
C
REAL*8 xGEO(3)
REAL*8 xGEI(3)
REAL*8 xS,yS,zS,cgst,sgst
C
COMMON /Soleil/xS,yS,zS,cgst,sgst
C
xGEI(1) = cgst*xGEO(1) - sgst*xGEO(2)
xGEI(2) = sgst*xGEO(1) + cgst*xGEO(2)
xGEI(3) = xGEO(3)
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE
& GEI_GEO(xGEI,xGEO)
C
IMPLICIT NONE
C
REAL*8 xGEI(3)
REAL*8 xGEO(3)
REAL*8 xS,yS,zS,cgst,sgst
C
COMMON /Soleil/xS,yS,zS,cgst,sgst
C
xGEO(1) = cgst*xGEI(1) + sgst*xGEI(2)
xGEO(2) = -sgst*xGEI(1) + cgst*xGEI(2)
xGEO(3) = xGEI(3)
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE
& GEO_SM(xGEO,xSM)
C
IMPLICIT NONE
C
REAL*8 xGEO(3)
REAL*8 xSM(3)
REAL*8 xGEI,yGEI,zGEI
REAL*8 xS,yS,zS,cgst,sgst
REAL*8 xD,yD,zD !GEI
REAL*8 xSD,ySD,zSD,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
COMMON /Soleil/xS,yS,zS,cgst,sgst
COMMON /sundip/xD,yD,zD,xSD,ySD,zSD
& ,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
xGEI = cgst*xGEO(1) - sgst*xGEO(2)
yGEI = sgst*xGEO(1) + cgst*xGEO(2)
zGEI = xGEO(3)
C
xSM(1) = xSDD*xGEI + ySDD*yGEI + zSDD*zGEI
xSM(2) = xSD*xGEI + ySD*yGEI + zSD*zGEI
xSM(3) = xD*xGEI + yD*yGEI + zD*zGEI
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE
& SM_GEO(xSM,xGEO)
C
IMPLICIT NONE
C
REAL*8 xSM(3)
REAL*8 xGEO(3)
REAL*8 xGEI,yGEI,zGEI
REAL*8 xS,yS,zS,cgst,sgst
REAL*8 xD,yD,zD !GEI
REAL*8 xSD,ySD,zSD,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
COMMON /Soleil/xS,yS,zS,cgst,sgst
COMMON /sundip/xD,yD,zD,xSD,ySD,zSD
& ,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
xGEI = xSDD*xSM(1) + xSD*xSM(2) + xD*xSM(3)
yGEI = ySDD*xSM(1) + ySD*xSM(2) + yD*xSM(3)
zGEI = zSDD*xSM(1) + zSD*xSM(2) + zD*xSM(3)
C
xGEO(1) = cgst*xGEI + sgst*yGEI
xGEO(2) = -sgst*xGEI + cgst*yGEI
xGEO(3) = zGEI
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE
& SM_GSM(xSM,xGSM)
C
IMPLICIT NONE
C
REAL*8 xSM(3)
REAL*8 xGSM(3)
REAL*8 xGEI,yGEI,zGEI
REAL*8 xS,yS,zS,cgst,sgst
REAL*8 xD,yD,zD !GEI
REAL*8 xSD,ySD,zSD,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
COMMON /Soleil/xS,yS,zS,cgst,sgst
COMMON /sundip/xD,yD,zD,xSD,ySD,zSD
& ,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
xGEI = xSDD*xSM(1) + xSD*xSM(2) + xD*xSM(3)
yGEI = ySDD*xSM(1) + ySD*xSM(2) + yD*xSM(3)
zGEI = zSDD*xSM(1) + zSD*xSM(2) + zD*xSM(3)
C
xGSM(1) = xS*xGEI + yS*yGEI + zS*zGEI
xGSM(2) = xSD*xGEI + ySD*yGEI + zSD*zGEI
xGSM(3) = xSSD*xGEI + ySSD*yGEI + zSSD*zGEI
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE
& GSM_SM(xGSM,xSM)
C
IMPLICIT NONE
C
REAL*8 xGSM(3)
REAL*8 xSM(3)
REAL*8 xGEI,yGEI,zGEI
REAL*8 xS,yS,zS,cgst,sgst
REAL*8 xD,yD,zD !GEI
REAL*8 xSD,ySD,zSD,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
COMMON /Soleil/xS,yS,zS,cgst,sgst
COMMON /sundip/xD,yD,zD,xSD,ySD,zSD
& ,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
xGEI = xS*xGSM(1) + xSD*xGSM(2) + xSSD*xGSM(3)
yGEI = yS*xGSM(1) + ySD*xGSM(2) + ySSD*xGSM(3)
zGEI = zS*xGSM(1) + zSD*xGSM(2) + zSSD*xGSM(3)
C
xSM(1) = xSDD*xGEI + ySDD*yGEI + zSDD*zGEI
xSM(2) = xSD*xGEI + ySD*yGEI + zSD*zGEI
xSM(3) = xD*xGEI + yD*yGEI + zD*zGEI
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE GEO_MAG(xGEO,xMAG)
C
IMPLICIT NONE
C
REAL*8 xGEO(3)
REAL*8 xMAG(3)
REAL*8 xc,yc,zc !Re
REAL*8 ct,st,cp,sp
REAL*8 Bo
C
COMMON /dipigrf/Bo,xc,yc,zc,ct,st,cp,sp
C
xMAG(1) = xGEO(1)*ct*cp + xGEO(2)*ct*sp - xGEO(3)*st
xMAG(2) = -xGEO(1)*sp + xGEO(2)*cp
xMAG(3) = xGEO(1)*st*cp + xGEO(2)*st*sp + xGEO(3)*ct
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE MAG_GEO(xMAG,xGEO)
C
IMPLICIT NONE
C
REAL*8 xGEO(3)
REAL*8 xMAG(3)
REAL*8 xc,yc,zc !Re
REAL*8 ct,st,cp,sp
REAL*8 Bo
C
COMMON /dipigrf/Bo,xc,yc,zc,ct,st,cp,sp
C
xGEO(1) = xMAG(1)*ct*cp - xMAG(2)*sp + xMAG(3)*st*cp
xGEO(2) = xMAG(1)*ct*sp + xMAG(2)*cp + xMAG(3)*st*sp
xGEO(3) = -xMAG(1)*st + xMAG(3)*ct
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE GEO_GSE(xGEO,xGSE)
C
IMPLICIT NONE
C
REAL*8 xGSE(3)
REAL*8 xGEO(3)
REAL*8 xGEI,yGEI,zGEI
REAL*8 xS,yS,zS,cgst,sgst
REAL*8 xD,yD,zD !GEI
REAL*8 xSD,ySD,zSD,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
REAL*8 aa,bb,y1,y2,y3
C
COMMON /Soleil/xS,yS,zS,cgst,sgst
COMMON /sundip/xD,yD,zD,xSD,ySD,zSD
& ,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
aa = -0.3978D0
bb = 0.9175D0
C
y1 = aa*zS-bb*yS
y2 = bb*xS
y3 = -aa*xS
C
xGEI = cgst*xGEO(1) - sgst*xGEO(2)
yGEI = sgst*xGEO(1) + cgst*xGEO(2)
zGEI = xGEO(3)
C
xGSE(1) = xS*xGEI + yS*yGEI + zS*zGEI
xGSE(2) = y1*xGEI + y2*yGEI + y3*zGEI
xGSE(3) = aa*yGEI + bb*zGEI
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE GSE_GEO(xGSE,xGEO)
C
IMPLICIT NONE
C
REAL*8 xGSE(3)
REAL*8 xGEO(3)
REAL*8 xGEI,yGEI,zGEI
REAL*8 xS,yS,zS,cgst,sgst
REAL*8 xD,yD,zD !GEI
REAL*8 xSD,ySD,zSD,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
REAL*8 aa,bb,y1,y2,y3,det
C
COMMON /Soleil/xS,yS,zS,cgst,sgst
COMMON /sundip/xD,yD,zD,xSD,ySD,zSD
& ,xSSD,ySSD,zSSD,xSDD,ySDD,zSDD
C
aa = -0.3978D0
bb = 0.9175D0
C
y1 = aa*zS-bb*yS
y2 = bb*xS
y3 = -aa*xS
C
det = xS*y2*bb+zS*y1*aa-xS*y3*aa-yS*y1*bb
C
xGEI = (y2*bb-y3*aa)*xGSE(1) - (ys*bb-zS*aa)*xGSE(2)
& + (yS*y3-zS*y2)*xGSE(3)
yGEI = - y1*bb*xGSE(1) + xS*bb*xGSE(2) - (xS*y3-zS*y1)*xGSE(3)
zGEI = y1*aa*xGSE(1) - xS*aa*xGSE(2) + (xS*y2-yS*y1)*xGSE(3)
C
xGEI = xGEI/det
yGEI = yGEI/det
zGEI = zGEI/det
C
xGEO(1) = cgst*xGEI + sgst*yGEI
xGEO(2) = -sgst*xGEI + cgst*yGEI
xGEO(3) = zGEI
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE SPH_CAR(r,lati,longi,x)
C
IMPLICIT NONE
C
REAL*8 r,lati,longi
REAL*8 colati
REAL*8 x(3)
REAL*8 pi,rad
common /rconst/rad,pi
C
CALL INITIZE
colati = pi/2.d0 - lati*rad
C
x(1) = r*SIN(colati)*COS(longi*rad)
x(2) = r*SIN(colati)*SIN(longi*rad)
x(3) = r*COS(colati)
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE CAR_SPH(x,r,lati,longi)
C
IMPLICIT NONE
REAL*8 r,lati,longi
REAL*8 SQ
REAL*8 x(3)
REAL*8 pi,rad
common /rconst/rad,pi
C
CALL INITIZE
SQ = x(1)*x(1) + x(2)*x(2)
r = SQRT(SQ + x(3)*x(3))
IF (SQ.NE.0.d0) GOTO 2
longi = 0.d0
IF (x(3).LT.0.d0) THEN
lati = -90.d0
ELSE
lati = 90.d0
ENDIF
RETURN
2 SQ = SQRT(SQ)
longi = ATAN2(x(2),x(1))/rad
lati = 90.d0 - ATAN2(SQ,x(3))/rad
IF (longi.LT.0.d0) longi = longi + 360.d0
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE SPH_CAR_VECT(r,lati,longi,sph,x)
C
IMPLICIT NONE
C
REAL*8 r,lati,longi
REAL*8 colati
REAL*8 st,ct,sp,cp
REAL*8 sph(3),x(3)
REAL*8 pi,rad
common /rconst/rad,pi
C
CALL INITIZE
colati = pi/2.d0 - lati*rad
C
st = sin(colati)
ct = cos(colati)
sp = sin(longi*rad)
cp = cos(longi*rad)
x(1) = sph(1) * st*cp + sph(2) * ct*cp - sph(3) * sp
x(2) = sph(1) * st*sp + sph(2) * ct*sp + sph(3) * cp
x(3) = sph(1) * ct - sph(2) * st
C
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE CAR_SPH_VECT(r,lati,longi,x,sph)
C
IMPLICIT NONE
C
REAL*8 r,lati,longi
REAL*8 colati
REAL*8 st,ct,sp,cp
REAL*8 sph(3),x(3)
REAL*8 pi,rad
common /rconst/rad,pi
C
CALL INITIZE
colati = pi/2.d0 - lati*rad
C
st = sin(colati)
ct = cos(colati)
sp = sin(longi*rad)
cp = cos(longi*rad)
sph(1) = x(3) * ct + x(1) * st*cp + x(2) * st*sp
sph(2) = -st * x(3) + x(1) * ct*cp + x(2) * ct*sp
sph(3) = - sp * x(1) + x(2) * cp
C
C
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
SUBROUTINE SUN(iyr,iday,secs,gst,Slong,Srasn,Sdec)
C
IMPLICIT NONE
C
C program to calculate sideral, time and position of the Sun
C good for years 1901 through 2099. accuracy 0.006 degree
C input is iyr,iday (INTEGER*4s), and secs, defining universal time
C output is Greenwich mean sidereal time (gst) in degrees,
C longitude along ecliptic (Slong), and apparent right ascension
C and declination (Srasn,Sdec) of the Sun, all in degrees.
C
INTEGER*4 iyr,iday
INTEGER*4 iaux
REAL*8 secs,gst,Slong,Srasn,Sdec
C
REAL*8 aux
REAL*8 dj,fday
REAL*8 t,vl,g,obliq,slp,sind,cosd
REAL*8 pi,rad
common /rconst/rad,pi
C
IF (iyr.LT.1901 .OR. iyr.GT.2099) RETURN
C
fday = secs/86400.D0
dj = 365.D0*(iyr-1900.D0) + INT((iyr-1901)/4)
& + iday + fday - 0.5D0
t = dj/36525.D0
aux = 279.696678D0 + 0.9856473354D0*dj
iaux = INT(aux/360.D0)
vl = aux-360.D0*iaux
aux = 279.690983D0 + 0.9856473354D0*dj + 360.D0*fday + 180.D0
iaux = INT(aux/360.D0)
gst = aux-360.D0*iaux
aux = 358.475845D0 + 0.985600267D0*dj
iaux = INT(aux/360.D0)
g = (aux-360.D0*iaux)*rad
C
Slong = vl + (1.91946D0-0.004789D0*t)*SIN(g)
& + 0.020094D0*SIN(2.D0*g)
obliq = (23.45229D0 - 0.0130125D0*t) * rad
slp = (Slong - 0.005686D0) * rad
sind = SIN(obliq) * SIN(slp)
cosd = SQRT(1.D0 - sind*sind)
Sdec = ATAN(sind/cosd) / rad
Srasn = 180.D0
& - ATAN2(sind/(cosd*TAN(obliq)), -COS(slp)/cosd) / rad
C
gst = gst * rad
Slong = Slong * rad
Sdec = Sdec * rad
Srasn = Srasn * rad
C
RETURN
END
C
!C----------------------------------------------------------------------
!C
! subroutine get_igrf_coeffs(year,g,h,ierr)
!c
!C SET UP TO ACCEPT DATES BETWEEN 1965 AND 2005; COEFFICIENTS
!C THROUGH 1985 ARE FROM DGRF MODELS COEFFICIENTS FOR 1990
!C FROM IGRF
!C [EOS TRANS. AGU APRIL 21, 1992, P. 182]. INTERPOLATION IS USED
!C FOR YEARS BETWEEN DGRF MODELS, AND EXTRAPOLATION FOR YEARS
!C BETWEEN 1990 AND 2000.
!C Coefficients for 1995 were added by D. Boscher
!C Modified by Paul O'Brien: ierr set, but uses nearest date limit
!c Date limits are 1965 to 2015
!c
!c Mauricio Peredo
!c Hughes STX at NASA/GSFC
!c June 29, 1995
!c bidouille Daniel Boscher Fev 20,1997
!c bidouille Daniel Boscher May 02,2001
!c
!C Additional IGRF References:
!C (J. GEOMAG. GEOELECTR.(1982), V.34, P.313-315,
!C GEOMAGN. AND AERONOMY (1986), V.26, P.523-525).
!C
!C------INPUT PARAMETERS:
!C iyr - YEAR NUMBER (FROM 1965 UP TO 2000)
!C----- OUTPUT PARAMETERS:
!C g,h - coefficients for the igrf model interpolated (or extrapolated)
!c to the epoch of iyr
!C
!C
! IMPLICIT NONE
!C
! REAL*8 F1,F2,DT,G(66),H(66),G65(66),H65(66),G70(66),H70(66),
! *G75(66),H75(66),G80(66),H80(66),G85(66),H85(66),G90(66),
! *H90(66),G95(66),H95(66),G00(66),H00(66),G05(66),H05(66),
! *G2010(66),H2010(66),
! *DG2010(45),DH2010(45)
!
! INTEGER*4 N,ierr
! REAL*8 year
! INTEGER*4 year_error_reported ! only report year error once
! save year_error_reported
! DATA year_error_reported/0/ ! initially not reported
!c
! DATA G65/0.D0,-30334.D0,-2119.D0,-1662.D0,2997.D0,1594.D0,
! *1297.D0,-2038.D0,1292.D0,856.D0,957.D0,804.D0,479.D0,-390.D0,
! *252.D0,-219.D0,358.D0,254.D0,-31.D0,-157.D0,-62.D0,
! *45.D0,61.D0,8.D0,-228.D0,4.D0,1.D0,-111.D0,75.D0,-57.D0,4.D0,
! *13.D0,-26.D0,-6.D0,13.D0,1.D0,13.D0,5.D0,-4.D0,-14.D0,0.D0,
! *8.D0,-1.D0,11.D0,4.D0,8.D0,10.D0,2.D0,-13.D0,10.D0,-1.D0,-1.D0,
! *5.D0,1.D0,-2.D0,
! *-2.D0,-3.D0,2.D0,-5.D0,-2.D0,4.D0,4.D0,0.D0,2.D0,2.D0,0.D0/
! DATA H65/0.D0,0.D0,5776.D0,0.D0,-2016.D0,114.D0,0.D0,-404.D0,
! *240.D0,-165.D0,0.D0,148.D0,-269.D0,13.D0,-269.D0,0.D0,19.D0,
! *128.D0,-126.D0,-97.D0,81.D0,0.D0,-11.D0,100.D0,68.D0,-32.D0,
! *-8.D0,-7.D0,0.D0,-61.D0,-27.D0,-2.D0,6.D0,26.D0,-23.D0,-12.D0,
! *0.D0,7.D0,-12.D0,9.D0,-16.D0,4.D0,24.D0,-3.D0,-17.D0,0.D0,-22.D0,
! *15.D0,7.D0,-4.D0,-5.D0,10.D0,10.D0,-4.D0,1.D0,0.D0,2.D0,1.D0,
! *2.D0,6.D0,-4.D0,0.D0,-2.D0,3.D0,0.D0,-6.D0/
!c
! DATA G70/0.D0,-30220.D0,-2068.D0,-1781.D0,3000.D0,1611.D0,
! *1287.D0,-2091.D0,1278.D0,838.D0,952.D0,800.D0,461.D0,-395.D0,
! *234.D0,-216.D0,359.D0,262.D0,-42.D0,-160.D0,-56.D0,
! *43.D0,64.D0,15.D0,-212.D0,2.D0,3.D0,-112.D0,72.D0,-57.D0,1.D0,
! *14.D0,-22.D0,-2.D0,13.D0,-2.D0,14.D0,6.D0,-2.D0,-13.D0,-3.D0,
! *5.D0,0.D0,11.D0,3.D0,8.D0,10.D0,2.D0,-12.D0,10.D0,-1.D0,0.D0,
! *3.D0,1.D0,-1.D0,-3.D0,-3.D0,2.D0,-5.D0,-1.D0,6.D0,4.D0,1.D0,
! *0.D0,3.D0,-1.D0/
! DATA H70/0.D0,0.D0,5737.D0,0.D0,-2047.D0,25.D0,0.D0,-366.D0,
! *251.D0,-196.D0,0.D0,167.D0,-266.D0,26.D0,-279.D0,0.D0,26.D0,
! *139.D0,-139.D0,-91.D0,83.D0,0.D0,-12.D0,100.D0,72.D0,-37.D0,
! *-6.D0,1.D0,0.D0,-70.D0,-27.D0,-4.D0,8.D0,23.D0,-23.D0,-11.D0,
! *0.D0,7.D0,-15.D0,6.D0,-17.D0,6.D0,21.D0,-6.D0,-16.D0,0.D0,
! *-21.D0,16.D0,6.D0,-4.D0,-5.D0,10.D0,11.D0,-2.D0,1.D0,0.D0,
! *1.D0,1.D0,3.D0,4.D0,-4.D0,0.D0,-1.D0,3.D0,1.D0,-4.D0/
!c
! DATA G75/0.D0,-30100.D0,-2013.D0,-1902.D0,3010.D0,1632.D0,
! *1276.D0,-2144.D0,1260.D0,830.D0,946.D0,791.D0,438.D0,-405.D0,
! *216.D0,-218.D0,356.D0,264.D0,-59.D0,-159.D0,-49.D0,
! *45.D0,66.D0,28.D0,-198.D0,1.D0,6.D0,-111.D0,71.D0,-56.D0,1.D0,
! *16.D0,-14.D0,0.D0,12.D0,-5.D0,14.D0,6.D0,-1.D0,-12.D0,-8.D0,
! *4.D0,0.D0,10.D0,1.D0,7.D0,10.D0,2.D0,-12.D0,10.D0,-1.D0,-1.D0,
! *4.D0,1.D0,-2.D0,-3.D0,-3.D0,2.D0,-5.D0,-2.D0,5.D0,4.D0,1.D0,
! *0.D0,3.D0,-1.D0/
! DATA H75/0.D0,0.D0,5675.D0,0.D0,-2067.D0,-68.D0,0.D0,-333.D0,
! *262.D0,-223.D0,0.D0,191.D0,-265.D0,39.D0,-288.D0,0.D0,31.D0,
! *148.D0,-152.D0,-83.D0,88.D0,0.D0,-13.D0,99.D0,75.D0,-41.D0,
! *-4.D0,11.D0,0.D0,-77.D0,-26.D0,-5.D0,10.D0,22.D0,-23.D0,-12.D0,
! *0.D0,6.D0,-16.D0,4.D0,-19.D0,6.D0,18.D0,-10.D0,-17.D0,0.D0,
! *-21.D0,16.D0,7.D0,-4.D0,-5.D0,10.D0,11.D0,-3.D0,1.D0,0.D0,1.D0,
! *1.D0,3.D0,4.D0,-4.D0,-1.D0,-1.D0,3.D0,1.D0,-5.D0/
!c
! DATA G80/0.D0,-29992.D0,-1956.D0,-1997.D0,3027.D0,1663.D0,
! *1281.D0,-2180.D0,1251.D0,833.D0,938.D0,782.D0,398.D0,-419.D0,
! *199.D0,-218.D0,357.D0,261.D0,-74.D0,-162.D0,-48.D0,48.D0,66.D0,
! *42.D0,-192.D0,4.D0,14.D0,-108.D0,72.D0,-59.D0,2.D0,21.D0,-12.D0,
! *1.D0,11.D0,-2.D0,18.D0,6.D0,0.D0,-11.D0,-7.D0,4.D0,3.D0,6.D0,
! *-1.D0,5.D0,10.D0,1.D0,-12.D0,9.D0,-3.D0,-1.D0,7.D0,2.D0,
! *-5.D0,-4.D0,-4.D0,2.D0,-5.D0,-2.D0,5.D0,3.D0,1.D0,2.D0,3.D0,0.D0/
! DATA H80/0.D0,0.D0,5604.D0,0.D0,-2129.D0,-200.D0,0.D0,-336.D0,
! *271.D0,-252.D0,0.D0,212.D0,-257.D0,53.D0,-297.D0,0.D0,46.D0,
! *150.D0,-151.D0,-78.D0,92.D0,0.D0,-15.D0,93.D0,71.D0,-43.D0,
! *-2.D0,17.D0,0.D0,-82.D0,-27.D0,-5.D0,16.D0,18.D0,-23.D0,-10.D0,
! *0.D0,7.D0,-18.D0,4.D0,-22.D0,9.D0,16.D0,-13.D0,-15.D0,0.D0,
! *-21.D0,16.D0,9.D0,-5.D0,-6.D0,9.D0,10.D0,-6.D0,2.D0,0.D0,1.D0,
! *0.D0,3.D0,6.D0,-4.D0,0.D0,-1.D0,4.D0,0.D0,-6.D0/
!c
! DATA G85/0.D0,-29873.D0,-1905.D0,-2072.D0,3044.D0,1687.D0,
! *1296.D0,-2208.D0,1247.D0,829.D0,936.D0,780.D0,361.D0,-424.D0,
! *170.D0,-214.D0,355.D0,253.D0,-93.D0,-164.D0,-46.D0,
! *53.D0,65.D0,51.D0,-185.D0,4.D0,16.D0,-102.D0,74.D0,-62.D0,3.D0,
! *24.D0,-6.D0,4.D0,10.D0,0.D0,21.D0,6.D0,0.D0,-11.D0,-9.D0,4.D0,
! *4.D0,4.D0,-4.D0,5.D0,10.D0,1.D0,-12.D0,9.D0,-3.D0,-1.D0,7.D0,
! *1.D0,-5.D0,-4.D0,-4.D0,3.D0,-5.D0,-2.D0,5.D0,3.D0,1.D0,2.D0,
! *3.D0,0.D0/
! DATA H85/0.D0,0.D0,5500.D0,0.D0,-2197.D0,-306.D0,0.D0,-310.D0,
! *284.D0,-297.D0,0.D0,232.D0,-249.D0,69.D0,-297.D0,0.D0,47.D0,
! *150.D0,-154.D0,-75.D0,95.D0,0.D0,-16.D0,88.D0,69.D0,-48.D0,
! *-1.D0,21.D0,0.D0,-83.D0,-27.D0,-2.D0,20.D0,17.D0,-23.D0,-7.D0,
! *0.D0,8.D0,-19.D0,5.D0,-23.D0,11.D0,14.D0,-15.D0,-11.D0,0.D0,
! *-21.D0,15.D0,9.D0,-6.D0,-6.D0,9.D0,9.D0,-7.D0,2.D0,0.D0,1.D0,
! *0.D0,3.D0,6.D0,-4.D0,0.D0,-1.D0,4.D0,0.D0,-6.D0/
!c
! DATA G90/0.D0,-29775.D0, -1848.D0, -2131.D0, 3059.D0, 1686.D0,
! * 1314.D0, -2239.D0, 1248.D0, 802.D0, 939.D0, 780.D0,
! * 325.D0, -423.D0, 141.D0, -214.D0, 353.D0, 245.D0,
! * -109.D0, -165.D0, -36.D0, 61.D0, 65.D0, 59.D0,
! * -178.D0, 3.D0, 18.D0, -96.D0, 77.D0, -64.D0,
! * 2.D0, 26.D0, -1.D0, 5.D0, 9.D0, 0.D0,
! * 23.D0, 5.D0, -1.D0, -10.D0, -12.D0, 3.D0,
! * 4.D0, 2.D0, -6.D0, 4.D0, 9.D0, 1.D0,
! * -12.D0, 9.D0, -4.D0, -2.D0, 7.D0, 1.D0,
! * -6.D0, -3.D0, -4.D0, 2.D0, -5.D0, -2.D0,
! * 4.D0, 3.D0, 1.D0, 3.D0, 3.D0, 0.D0/
! DATA H90/0.D0, 0.D0, 5406.D0, 0.D0, -2279.D0, -373.D0,
! * 0.D0, -284.D0, 293.D0, -352.D0, 0.D0, 247.D0,
! * -240.D0, 84.D0, -299.D0, 0.D0, 46.D0, 154.D0,
! * -153.D0, -69.D0, 97.D0, 0.D0, -16.D0, 82.D0,
! * 69.D0, -52.D0, 1.D0, 24.D0, 0.D0, -80.D0,
! * -26.D0, 0.D0, 21.D0, 17.D0, -23.D0, -4.D0,
! * 0.D0, 10.D0, -19.D0, 6.D0, -22.D0, 12.D0,
! * 12.D0, -16.D0, -10.D0, 0.D0, -20.D0, 15.D0,
! * 11.D0, -7.D0, -7.D0, 9.D0, 8.D0, -7.D0,
! * 2.D0, 0.D0, 2.D0, 1.D0, 3.D0, 6.D0,
! * -4.D0, 0.D0, -2.D0, 3.D0, -1.D0, -6.D0/
!C
! DATA G95/0.D0,-29692.D0, -1784.D0, -2200.D0, 3070.D0, 1681.D0,
! * 1335.D0, -2267.D0, 1249.D0, 759.D0, 940.D0, 780.D0,
! * 290.D0, -418.D0, 122.D0, -214.D0, 352.D0, 235.D0,
! * -118.D0, -166.D0, -17.D0, 68.D0, 67.D0, 68.D0,
! * -170.D0, -1.D0, 19.D0, -93.D0, 77.D0, -72.D0,
! * 1.D0, 28.D0, 5.D0, 4.D0, 8.D0, -2.D0,
! * 25.D0, 6.D0, -6.D0, -9.D0, -14.D0, 9.D0,
! * 6.D0, -5.D0, -7.D0, 4.D0, 9.D0, 3.D0,
! * -10.D0, 8.D0, -8.D0, -1.D0, 10.D0, -2.D0,
! * -8.D0, -3.D0, -6.D0, 2.D0, -4.D0, -1.D0,
! * 4.D0, 2.D0, 2.D0, 5.D0, 1.D0, 0.D0/
! DATA H95/0.D0, 0.D0, 5306.D0, 0.D0, -2366.D0, -413.D0,
! * 0.D0, -262.D0, 302.D0, -427.D0, 0.D0, 262.D0,
! * -236.D0, 97.D0, -306.D0, 0.D0, 46.D0, 165.D0,
! * -143.D0, -55.D0, 107.D0, 0.D0, -17.D0, 72.D0,
! * 67.D0, -58.D0, 1.D0, 36.D0, 0.D0, -69.D0,
! * -25.D0, 4.D0, 24.D0, 17.D0, -24.D0, -6.D0,
! * 0.D0, 11.D0, -21.D0, 8.D0, -23.D0, 15.D0,
! * 11.D0, -16.D0, -4.D0, 0.D0, -20.D0, 15.D0,
! * 12.D0, -6.D0, -8.D0, 8.D0, 5.D0, -8.D0,
! * 3.D0, 0.D0, 1.D0, 0.D0, 4.D0, 5.D0,
! * -5.D0, -1.D0, -2.D0, 1.D0, -2.D0, -7.D0/
!C
! DATA G00/0.0D0,-29619.4D0,-1728.2D0,-2267.7D0,3068.4D0,1670.9D0,
! * 1339.6D0, -2288.0D0, 1252.1D0, 714.5D0, 932.3D0, 786.8D0,
! * 250.0D0, -403.0D0, 111.3D0, -218.8D0, 351.4D0, 222.3D0,
! * -130.4D0, -168.6D0, -12.9D0, 72.3D0, 68.2D0, 74.2D0,
! * -160.9D0, -5.9D0, 16.9D0, -90.4D0, 79.0D0, -74.0D0,
! * 0.0D0, 33.3D0, 9.1D0, 6.9D0, 7.3D0, -1.2D0,
! * 24.4D0, 6.6D0, -9.2D0, -7.9D0, -16.6D0, 9.1D0,
! * 7.0D0, -7.9D0, -7.0D0, 5.0D0, 9.4D0, 3.0D0,
! * -8.4D0, 6.3D0, -8.9D0, -1.5D0, 9.3D0, -4.3D0,
! * -8.2D0, -2.6D0, -6.0D0, 1.7D0, -3.1D0, -0.5D0,
! * 3.7D0, 1.0D0, 2.0D0, 4.2D0, 0.3D0, -1.1D0/
! DATA H00/0.0D0, 0.0D0, 5186.1D0, 0.0D0,-2481.6D0,-458.0D0,
! * 0.0D0, -227.6D0, 293.4D0, -491.1D0, 0.0D0, 272.6D0,
! * -231.9D0, 119.8D0, -303.8D0, 0.0D0, 43.8D0, 171.9D0,
! * -133.1D0, -39.3D0, 106.3D0, 0.0D0, -17.4D0, 63.7D0,
! * 65.1D0, -61.2D0, 0.7D0, 43.8D0, 0.0D0, -64.6D0,
! * -24.2D0, 6.2D0, 24.0D0, 14.8D0, -25.4D0, -5.8D0,
! * 0.0D0, 11.9D0, -21.5D0, 8.5D0, -21.5D0, 15.5D0,
! * 8.9D0, -14.9D0, -2.1D0, 0.0D0, -19.7D0, 13.4D0,
! * 12.5D0, -6.2D0, -8.4D0, 8.4D0, 3.8D0, -8.2D0,
! * 4.8D0, 0.0D0, 1.7D0, 0.0D0, 4.0D0, 4.9D0,
! * -5.9D0, -1.2D0, -2.9D0, 0.2D0, -2.2D0, -7.4D0/
!C
! DATA G05/0.0D0,-29554.63D0,-1669.05D0,-2337.24D0,3047.69D0,
! * 1657.76D0, 1336.3D0,-2305.83D0, 1246.39D0, 672.51D0,
! * 920.55D0, 797.96D0, 210.65D0, -379.86D0, 100.0D0,
! * -227.D0, 354.41D0, 208.95D0, -136.54D0,-168.05D0,
! * -13.55D0, 73.60D0, 69.56D0, 76.74D0,-151.34D0,
! * -14.58D0, 14.58D0, -86.36D0, 79.88D0, -74.46D0,
! * -1.65D0, 38.73D0, 12.3D0, 9.37D0, 5.42D0,
! * 1.94D0, 24.8D0, 7.62D0, -11.73D0, -6.88D0,
! * -18.11D0, 10.17D0, 9.36D0, -11.25D0, -4.87D0,
! * 5.58D0, 9.76D0, 3.58D0, -6.94D0, 5.01D0,
! * -10.76D0, -1.25D0, 8.76D0, -6.66D0, -9.22D0,
! * -2.17D0, -6.12D0, 1.42D0, -2.35D0, -0.15D0,
! * 3.06D0, 0.29D0, 2.06D0, 3.77D0, -0.21D0,
! * -2.09D0/
! DATA H05/0.0D0, 0.0D0, 5077.99D0, 0.0D0, -2594.5D0,
! * -515.43D0, 0.0D0, -198.86D0, 269.72D0, -524.72D0,
! * 0.0D0, 282.07D0,-225.23D0, 145.15D0, -305.36D0,
! * 0.0D0, 42.72D0, 180.25D0,-123.45D0, -19.57D0,
! * 103.85D0, 0.0D0, -20.33D0, 54.75D0, 63.63D0,
! * -63.53D0, 0.24D0, 50.94D0, 0.0D0, -61.14D0,
! * -22.57D0, 6.82D0, 25.35D0, 10.93D0, -26.32D0,
! * -4.64D0, 0.0D0, 11.20D0, -20.88D0, 9.83D0,
! * -19.71D0, 16.22D0, 7.61D0, -12.76D0, -0.06D0,
! * 0.0D0, -20.11D0, 12.67D0, 12.7D0, -6.72D0,
! * -8.16D0, 8.10D0, 2.92D0, -7.73D0, 6.01D0,
! * 0.0D0, 2.19D0, 0.10D0, 4.46D0, 4.76D0,
! * -6.58D0, -1.01D0, -3.47D0, -0.86D0, -2.31D0,
! * -7.93D0/
!C
! DATA G2010/
! * 0.00d0, !n=00 m=00-00
! * -29496.50d0, -1585.90d0, !n=01 m=00-01
! * -2396.60d0, 3026.00d0, 1668.60d0, !n=02 m=00-02
! * 1339.70d0, -2326.30d0, 1231.70d0, 634.20d0, !n=03 m=00-03
! * 912.60d0, 809.00d0, 166.60d0, -357.10d0, 89.70d0, !n=04 m=00-04
! * -231.10d0, 357.20d0, 200.30d0, -141.20d0, -163.10d0,
! * -7.70d0, !n=05 m=00-05
! * 72.80d0, 68.60d0, 76.00d0, -141.40d0, -22.90d0,
! * 13.10d0, -77.90d0, !n=06 m=00-06
! * 80.40d0, -75.00d0, -4.70d0, 45.30d0, 14.00d0,
! * 10.40d0, 1.60d0, 4.90d0, !n=07 m=00-07
! * 24.30d0, 8.20d0, -14.50d0, -5.70d0, -19.30d0,
! * 11.60d0, 10.90d0, -14.10d0, -3.70d0, !n=08 m=00-08
! * 5.40d0, 9.40d0, 3.40d0, -5.30d0, 3.10d0,
! * -12.40d0, -0.80d0, 8.40d0, -8.40d0, -10.10d0, !n=09 m=00-09
! * -2.00d0, -6.30d0, 0.90d0, -1.10d0, -0.20d0,
! * 2.50d0, -0.30d0, 2.20d0, 3.10d0, -1.00d0,
! * -2.80d0 !n=10 m=00-10
! * /
!
! DATA H2010/
! * 0.00d0, !n=00 m=00-00
! * 0.00d0, 4945.10d0, !n=01 m=00-01
! * 0.00d0, -2707.70d0, -575.40d0, !n=02 m=00-02
! * 0.00d0, -160.50d0, 251.70d0, -536.80d0, !n=03 m=00-03
! * 0.00d0, 286.40d0, -211.20d0, 164.40d0, -309.20d0, !n=04 m=00-04
! * 0.00d0, 44.70d0, 188.90d0, -118.10d0, 0.10d0,
! * 100.90d0, !n=05 m=00-05
! * 0.00d0, -20.80d0, 44.20d0, 61.50d0, -66.30d0,
! * 3.10d0, 54.90d0, !n=06 m=00-06
! * 0.00d0, -57.80d0, -21.20d0, 6.60d0, 24.90d0,
! * 7.00d0, -27.70d0, -3.40d0, !n=07 m=00-07
! * 0.00d0, 10.90d0, -20.00d0, 11.90d0, -17.40d0,
! * 16.70d0, 7.10d0, -10.80d0, 1.70d0, !n=08 m=00-08
! * 0.00d0, -20.50d0, 11.60d0, 12.80d0, -7.20d0,
! * -7.40d0, 8.00d0, 2.20d0, -6.10d0, 7.00d0, !n=09 m=00-09
! * 0.00d0, 2.80d0, -0.10d0, 4.70d0, 4.40d0,
! * -7.20d0, -1.00d0, -4.00d0, -2.00d0, -2.00d0,
! * -8.30d0 !n=10 m=00-10
! * /
!
! DATA DG2010/
! * 0.00d0, !n=00 m=00-00
! * 11.40d0, 16.70d0, !n=01 m=00-01
! * -11.30d0, -3.90d0, 2.70d0, !n=02 m=00-02
! * 1.30d0, -3.90d0, -2.90d0, -8.10d0, !n=03 m=00-03
! * -1.40d0, 2.00d0, -8.90d0, 4.40d0, -2.30d0, !n=04 m=00-04
! * -0.50d0, 0.50d0, -1.50d0, -0.70d0, 1.30d0,
! * 1.40d0, !n=05 m=00-05
! * -0.30d0, -0.30d0, -0.30d0, 1.90d0, -1.60d0,
! * -0.20d0, 1.80d0, !n=06 m=00-06
! * 0.20d0, -0.10d0, -0.60d0, 1.40d0, 0.30d0,
! * 0.10d0, -0.80d0, 0.40d0, !n=07 m=00-07
! * -0.10d0, 0.10d0, -0.50d0, 0.30d0, -0.30d0,
! * 0.30d0, 0.20d0, -0.50d0, 0.20d0 !n=08 m=00-08
! * /
!
! DATA DH2010/
! * 0.00d0, !n=00 m=00-00
! * 0.00d0, -28.80d0, !n=01 m=00-01
! * 0.00d0, -23.00d0, -12.90d0, !n=02 m=00-02
! * 0.00d0, 8.60d0, -2.90d0, -2.10d0, !n=03 m=00-03
! * 0.00d0, 0.40d0, 3.20d0, 3.60d0, -0.80d0, !n=04 m=00-04
! * 0.00d0, 0.50d0, 1.50d0, 0.90d0, 3.70d0,
! * -0.60d0, !n=05 m=00-05
! * 0.00d0, -0.10d0, -2.10d0, -0.40d0, -0.50d0,
! * 0.80d0, 0.50d0, !n=06 m=00-06
! * 0.00d0, 0.60d0, 0.30d0, -0.20d0, -0.10d0,
! * -0.80d0, -0.30d0, 0.20d0, !n=07 m=00-07
! * 0.00d0, 0.00d0, 0.20d0, 0.50d0, 0.40d0,
! * 0.10d0, -0.10d0, 0.40d0, 0.40d0 !n=08 m=00-08
! * /
!c
!C
! if ( (year.lt.1965.D0).or.(year.gt.2015.D0) ) then
! ierr=1
! if (year_error_reported .eq. 0) then
! write(*,999) year
! year_error_reported = 1
! endif
!c goto 300 ! continue, just throw an error
! endif
!c
! IF (year.LT.1970.D0) GOTO 50 !INTERPOLATE BETWEEN 1965 - 1970
! IF (year.LT.1975.D0) GOTO 60 !INTERPOLATE BETWEEN 1970 - 1975
! IF (year.LT.1980.D0) GOTO 70 !INTERPOLATE BETWEEN 1975 - 1980
! IF (year.LT.1985.D0) GOTO 80 !INTERPOLATE BETWEEN 1980 - 1985
! IF (year.LT.1990.D0) GOTO 90 !INTERPOLATE BETWEEN 1985 - 1990
! IF (year.LT.1995.D0) GOTO 100 !INTERPOLATE BETWEEN 1990 - 1995
! IF (year.LT.2000.D0) GOTO 110 !INTERPOLATE BETWEEN 1995 - 2000
! IF (year.LT.2005.D0) GOTO 120 !INTERPOLATE BETWEEN 2000 - 2005
! IF (year.LT.2010.D0) GOTO 130 !INTERPOLATE BETWEEN 2005 - 2010
!C
!C EXTRAPOLATE BETWEEN 2010 - 2020
!C
! DT=year-2010.D0
! if (DT.gt.10.D0) then
! DT = 10.D0 ! effectively 2020
! endif
! DO 40 N=1,66
! G(N)=G2010(N)
! H(N)=H2010(N)
! IF (N.GT.45) GOTO 40
! G(N)=G(N)+DG2010(N)*DT
! H(N)=H(N)+DH2010(N)*DT
!40 CONTINUE
! GOTO 300
!C
!C
!C INTERPOLATE BETWEEEN 1965 - 1970
!C
!50 if (year.gt.1965.D0) then
! F2=(year-1965.D0)/5.D0
! else
! F2=0.D0 ! effectively 1965
! endif
! F1=1.D0-F2
! DO 55 N=1,66
! G(N)=G65(N)*F1+G70(N)*F2
!55 H(N)=H65(N)*F1+H70(N)*F2
! GOTO 300
!C
!C
!C INTERPOLATE BETWEEN 1970 - 1975
!C
!60 F2=(year-1970.D0)/5.D0
! F1=1.D0-F2
! DO 65 N=1,66
! G(N)=G70(N)*F1+G75(N)*F2
!65 H(N)=H70(N)*F1+H75(N)*F2
! GOTO 300
!C
!C
!C INTERPOLATE BETWEEN 1975 - 1980
!C
!70 F2=(year-1975.D0)/5.D0
! F1=1.D0-F2
! DO 75 N=1,66
! G(N)=G75(N)*F1+G80(N)*F2
!75 H(N)=H75(N)*F1+H80(N)*F2
! GOTO 300
!C
!C
!C INTERPOLATE BETWEEN 1980 - 1985
!C
!80 F2=(year-1980.D0)/5.D0
! F1=1.D0-F2
! DO 85 N=1,66
! G(N)=G80(N)*F1+G85(N)*F2
!85 H(N)=H80(N)*F1+H85(N)*F2
! GOTO 300
!C
!C
!C INTERPOLATE BETWEEN 1985 - 1990
!C
!90 F2=(year-1985.D0)/5.D0
! F1=1.D0-F2
! DO 95 N=1,66
! G(N)=G85(N)*F1+G90(N)*F2
!95 H(N)=H85(N)*F1+H90(N)*F2
! GOTO 300
!C
!C
!C INTERPOLATE BETWEEN 1990 - 1995
!C
!100 F2=(year-1990.D0)/5.
! F1=1.D0-F2
! DO 105 N=1,66
! G(N)=G90(N)*F1+G95(N)*F2
!105 H(N)=H90(N)*F1+H95(N)*F2
! GOTO 300
!C
!C
!C INTERPOLATE BETWEEN 1995 - 2000
!C
!110 F2=(year-1995.D0)/5.D0
! F1=1.D0-F2
! DO 115 N=1,66
! G(N)=G95(N)*F1+G00(N)*F2
!115 H(N)=H95(N)*F1+H00(N)*F2
! GOTO 300
!C
!C INTERPOLATE BETWEEN 2000 - 2005
!C
!120 F2=(year-2000.D0)/5.D0
! F1=1.D0-F2
! DO 125 N=1,66
! G(N)=G00(N)*F1+G05(N)*F2
!125 H(N)=H00(N)*F1+H05(N)*F2
! GOTO 300
!C
!C INTERPOLATE BETWEEN 2005 - 2010
!C
!130 F2=(year-2005.D0)/5.D0
! F1=1.D0-F2
! DO 135 N=1,66
! G(N)=G05(N)*F1+G2010(N)*F2
!135 H(N)=H05(N)*F1+H2010(N)*F2
! GOTO 300
!C
!C
!300 continue
!c
!C
!999 format(//1x,
! * '*** WARNING -- Input year = ',F7.2,/
! * ' is out of valid range 1965-2015. Using nearest ***'//)
!c
! return
! end
!C
C---------------------------------------------------------------------
C
c
subroutine get_terms(g,h,alpha,beta,x0,y0,z0,b0)
c
implicit none
c
REAL*8 g(66),h(66),alpha,beta,x0,y0,z0
REAL*8 g10,g11,h11,g20,g21,g22,h21,h22
REAL*8 b02,b0,sq3,l0,l1,l2,e
c
g10 = g(2)
g11 = g(3)
h11 = h(3)
g20 = g(4)
g21 = g(5)
g22 = g(6)
h21 = h(5)
h22 = h(6)
c
b02 = g10*g10 + g11*g11 + h11*h11
b0 = sqrt(b02)
c write(6,*)b0
alpha = acos(-g10/b0)
beta = atan(h11/g11)
c
sq3 = sqrt(3.D0)
l0 = 2.D0*g10*g20 + sq3*(g11*g21 + h11*h21)
l1 = -g11*g20 + sq3*(g10*g21 + g11*g22 + h11*h22)
l2 = -h11*g20 + sq3*(g10*h21 - h11*g22 + g11*h22)
e = (l0*g10 + l1*g11 + l2*h11)/(4.*b02)
c
z0 = (l0 - g10*e)/(3.D0*b02)
x0 = (l1 - g11*e)/(3.D0*b02)
y0 = (l2 - h11*e)/(3.D0*b02)
! write(6,*)alpha*180.d0/3.14d0,beta*180.d0/3.14d0
c
return
end
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE calcul_GH1
C
IMPLICIT NONE
C
INTEGER*4 I,J,L,M,N
REAL*8 g(66),h(66)
REAL*8 GH1(144),GH(144)
REAL*8 X,F,F0
C
COMMON /MODEL/GH1
COMMON /dgrf/G,H
C
L = 0
M = 0
DO I = 1,10
DO J = 0,I
M = M + 1
L = L+1
GH(L) = G(M+1)
IF(J.NE.0) THEN
L = L+1
GH(L) = H(M+1)
ENDIF
ENDDO
ENDDO
C
GH1(1) = 0.0D0
I=2
F0= -1.D0
DO 9 N=1,10
X = N
F0 = F0 * X * X / (4.D0 * X - 2.D0)
F0 = F0 * (2.D0 * X - 1.D0) / X
F = F0 * 0.5D0 * SQRT(2.D0)
GH1(I) = GH(I-1) * F0
I = I+1
DO 9 M=1,N
F = F * (X + M) / (X - M + 1.D0)
F = F * SQRT((X - M + 1.D0) / (X + M))
GH1(I) = GH(I-1) * F
GH1(I+1) = GH(I) * F
I=I+2
9 CONTINUE
C
RETURN
END
C
C*********************************************************************
SUBROUTINE RLL_GDZ (rr,lati,longi,alti)
C rr adimensionne, alti en km, longi en degres, lati en degres
C
IMPLICIT NONE
C
REAL*8 lati,longi
REAL*8 alti,rr
REAL*8 CT,ST,CP,SP
REAL*8 D
REAL*8 bb,cc
REAL*8 ERA,AQUAD,BQUAD
C
COMMON/GENER/ERA,AQUAD,BQUAD
REAL*8 pi,rad
common /rconst/rad,pi
C
CALL INITIZE
C
CT = DSIN(lati*rad)
ST = DCOS(lati*rad)
D = DSQRT(AQUAD-(AQUAD-BQUAD)*CT*CT)
CP = COS(longi*rad)
SP = SIN(longi*rad)
bb = (BQUAD*CT*CT+AQUAD*ST*ST)/D
cc = (BQUAD*BQUAD*CT*CT+AQUAD*AQUAD*ST*ST)/D/D-rr*rr*ERA*ERA
alti = -bb+DSQRT(bb*bb-cc)
C
RETURN
END
C
C*********************************************************************
C alti en km
C lati,longi en degres
C x,y,z GEO en adimensionne
C
SUBROUTINE GDZ_GEO(lati,longi,alti,xx,yy,zz)
C
IMPLICIT NONE
C
REAL*8 lati,longi
REAL*8 alti,xx,yy,zz
REAL*8 CT,ST,CP,SP
REAL*8 D,RHO
REAL*8 ERA,AQUAD,BQUAD
C
COMMON/GENER/ERA,AQUAD,BQUAD
REAL*8 pi,rad
common /rconst/rad,pi
C
CALL INITIZE
C
CT = SIN(lati*rad)
ST = COS(lati*rad)
D = SQRT(AQUAD-(AQUAD-BQUAD)*CT*CT)
CP = COS(longi*rad)
SP = SIN(longi*rad)
zz = (alti+BQUAD/D)*CT/ERA
RHO= (alti+AQUAD/D)*ST/ERA
xx = RHO*CP
yy = RHO*SP
C
RETURN
END
C
C*********************************************************************
C x,y,z en adimensionne
C lati,longi en degres
C alti en km
C
SUBROUTINE GEO_GDZ(xx,yy,zz,lati,longi,alti)
C
IMPLICIT NONE
C
REAL*8 precision
PARAMETER (precision = 1.D-5) ! increase precision to 1E-5
C
REAL*8 lati,longi,lat0
REAL*8 alti,xx,yy,zz,alt0
REAL*8 D,RHO
REAL*8 ERA,AQUAD,BQUAD
integer*4 i
C
COMMON/GENER/ERA,AQUAD,BQUAD
REAL*8 pi,rad
common /rconst/rad,pi
C
CALL INITIZE
longi = ATAN2(yy,xx)/rad
RHO = SQRT(xx*xx+yy*yy)
lati = ATAN2(zz,RHO)
D = COS(lati)
IF (D.LT.1.D-15)THEN
lati = lati/rad
alti = (zz-1.D0)*SQRT(BQUAD)
ELSE
alti = RHO/D - 1.D0
C
i=0
10 CONTINUE
alt0 = alti
lat0 = lati
D = SQRT(AQUAD-(AQUAD-BQUAD)*sin(lat0)*sin(lat0))
lati = ATAN2(zz*(alt0+AQUAD/D/ERA),RHO*(alt0+BQUAD/D/ERA))
alti = RHO/cos(lati) - AQUAD/D/ERA
i=i+1
if (i .gt. 1000) then
alti=0.
lati=0.
!write(6,*)'geo2gdz ',i,alti
return
endif
IF (ABS(alti - alt0).GT. precision .OR. ABS(lati -
& lat0) .GT. precision) GOTO 10 ! keep looping until BOTH are within precision
alti = alti*ERA
lati = lati/rad
ENDIF
!write(6,*)'geo2gdz ',i,alti
RETURN
END
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
SUBROUTINE INITIZE
C----------------------------------------------------------------
C Initializes the parameters in COMMON/GENER/
C
C UMR = ATAN(1.0)*4./180. *UMR=
C ERA EARTH RADIUS FOR NORMALIZATION OF CARTESIAN
C COORDINATES (6371.2 KM)
C EREQU MAJOR HALF AXIS FOR EARTH ELLIPSOID (6378.160 KM)
C ERPOL MINOR HALF AXIS FOR EARTH ELLIPSOID (6356.775 KM)
C AQUAD SQUARE OF MAJOR HALF AXIS FOR EARTH ELLIPSOID
C BQUAD SQUARE OF MINOR HALF AXIS FOR EARTH ELLIPSOID
C
C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL
C ASTRONOMICAL UNION .
C-----------------------------------------------------------------
IMPLICIT NONE
REAL*8 ERA,AQUAD,BQUAD,EREQU,ERPOL
COMMON/GENER/ERA,AQUAD,BQUAD
real*8 rad,pi
common/rconst/rad,pi
ERA=6371.2D0
c WGS84 World Geodetic System 84 (GPS)
EREQU=6378.137D0
ERPOL=6356.752314D0
c Older World Geodetic System
c EREQU=6378.16D0
c ERPOL=6356.775D0
c
AQUAD=EREQU*EREQU
BQUAD=ERPOL*ERPOL
pi = 4.D0*ATAN(1.D0)
rad = pi/180.D0
RETURN
END