JB2006

Fortran Source Code
SET home
[Home]
[Intro] [Indices] [Code]
[Publications]
[Contacts]
[Figures]
[License]
[SET/Spacewx]

Instructions:

REV A - Oct 2006

Smoothing of density corrections and scale height through different altitude bands in the latitude-local time correction subroutine DTSUB, with dTx correction replaced with dTc correction


The FORTRAN source code is composed of two files: A main "Test Driver" and a collection of subroutines. 

Because of variations in preferences and client settings in retrieving source code, serveral methods are provided here to retrieve the files:
     (You may also wish to separate the main driver from the subroutines for your specific applications.)
C
C TEST CASE FOR RUNNING JB2006
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION SUN(2),SAT(3),GEO(3),TEMP(2)
C
C
C INPUT TO JB2006:
C
C AMJD : Date and Time, in modified Julian Days
C and Fraction (MJD = JD-2400000.5)
C SUN(1) : Right Ascension of Sun (radians)
C SUN(2) : Declination of Sun (radians)
C SAT(1) : Right Ascension of Position (radians)
C SAT(2) : Geocentric Latitude of Position (radians)
C SAT(3) : Height of Position (km)
C GEO(1) : 10.7-cm Solar Flux (1.0E-22*Watt/(M**2*Hertz))
C (Tabular time 1.0 day earlier)
C GEO(2) : 10.7-cm Solar Flux, ave.
C 81-day centered on the input time
C GEO(3) : Geomagnetic Planetary 3-hour Index
C A-sub-p for a tabular time 0.279 days earlier
C (6.7 hours earlier)
C S10 : EUV index (26-34 nm) scaled to F10
C (Tabular time 1.0 day earlier)
C S10B : EUV 81-day ave. centered index
C XM10 : MG2 index scaled to F10
C XM10B : MG2 81-day ave. centered index
C (Tabular time 5.0 days earlier)
C
C OUTPUT FROM JB2006:
C
C TEMP(1): Exospheric Temperature above Input Position (deg K)
C TEMP(2): Temperature at Input Position (deg K)
C RHO : Total Mass-Desnity at Input Position (kg/m**3)
C
C
DATA PI/3.1415927D0/
C
C
C SET SOLAR INDICES
C
C USE 1 DAY LAG FOR EUV AND F10 INFLUENCE
S10 = 140
S10B = 100
F10 = 135
F10B = 95
C
C USE 5 DAY LAG FOR MG FUV INFLUENCE
XM10 = 130
XM10B = 95
C USE 6.7 HR LAG FOR ap INFLUENCE
AP = 30
C
GEO(1) = F10
GEO(2) = F10B
GEO(3) = AP
C
C SET POINT OF INTEREST LOCATION (RADIANS AND KM)
C
SAT(1) = 90.*PI/180.
SAT(2) = 45.*PI/180.
SAT(3) = 400.
C
C SET SUN LOCATION (RADIANS)
C
SUN(1) = 90.*PI/180.
SUN(2) = 20.*PI/180.
C
C SET TIME OF INTEREST
C
IYR = 01
IDAY = 200
IF (IYR.LT.50) IYR = IYR + 100
IYY = (IYR-50)*365 + ((IYR-1)/4-12)
ID1950 = IYY + IDAY
D1950 = ID1950
AMJD = D1950 + 33281.D0
C
C
C COMPUTE DENSITY KG/M3 RHO
C
CALL JB2006 (AMJD,SUN,SAT,GEO,S10,S10B,XM10,XM10B,TEMP,RHO)
C
C
WRITE(*,100) D1950,S10,S10B,F10,F10B,XM10,XM10B,AP,TEMP,RHO
100 FORMAT(F8.0,2X,7F5.0,2F10.1,D14.4)
C
C
C OUTPUT RESULTS:
C
C 18828. 140. 100. 135. 95. 130. 95. 30. 1145.8 1137.7 0.4066D-11
C
C
STOP
END
C
SUBROUTINE JB2006 (AMJD,SUN,SAT,GEO,S10,S10B,XM10,XM10B,TEMP,RHO)
C
C Jacchia-Bowman 2006 Model Atmosphere
C
C
C***********************************************************************
C
C This is the CIRA "Integration Form" of a Jacchia Model.
C There are no tabular values of density. Instead, the barometric
C equation and diffusion equation are integrated numerically using
C the Newton-Coates method to produce the density profile up to the
C input position.
C
C INPUT:
C
C AMJD : Date and Time, in modified Julian Days
C and Fraction (MJD = JD-2400000.5)
C SUN(1) : Right Ascension of Sun (radians)
C SUN(2) : Declination of Sun (radians)
C SAT(1) : Right Ascension of Position (radians)
C SAT(2) : Geocentric Latitude of Position (radians)
C SAT(3) : Height of Position (km)
C GEO(1) : 10.7-cm Solar Flux (1.0E-22*Watt/(M**2*Hertz))
C (Tabular time 1.0 day earlier)
C GEO(2) : 10.7-cm Solar Flux, ave.
C 81-day centered on the input time
C GEO(3) : Geomagnetic Planetary 3-hour Index
C A-sub-p for a tabular time 0.279 days earlier
C (6.7 hours earlier)
C S10 : EUV index (26-34 nm) scaled to F10
C (Tabular time 1.0 day earlier)
C S10B : EUV 81-day ave. centered index
C XM10 : MG2 index scaled to F10
C XM10B : MG2 81-day ave. centered index
C (Tabular time 5.0 days earlier)
C
C OUTPUT:
C
C TEMP(1): Exospheric Temperature above Input Position (deg K)
C TEMP(2): Temperature at Input Position (deg K)
C RHO : Total Mass-Desnity at Input Position (kg/m**3)
C
C
C JB2006 Model Development: (Ref. 1)
C
C 1. Start with the CIRA72 model (Jacchia 71).
C
C 2. Convert to a CIRA70 model replacing equations from Jacchia 70
C model (Ref. 5)
C
C 3. Replace Tc equation using new solar indices (Ref. 1 and 2)
C
C 4. Replace semiannual equation with new global model based
C on F10B (Ref. 1 and 3)
C
C 5. Add correction for local solar time and latitude errors
C (Ref. 1)
C Added smooth transition between altitude bands
C
C 6. Add high altitude ( z > 1500 km ) correction
C (Ref. 1 and 4)
C
C
C All equation references below refer to the original
C Jacchia model papers.
C
C
C References:
C
C 1. Bowman, Bruce R., etc. : "A New Empirical Thermospheric
C Density Model JB2006 Using New Solar Indices",
C AIAA/AAS Astrodynamics Specialists Conference, Keystone, CO,
C 21-24 Aug 2006, (Paper AIAA 2006-6166).
C
C 2. Bowman, Bruce R., etc. : "Improvements in Modeling
C Thermospheric Densities Using New EUV and FUV Solar Indices",
C AAS/AIAA Space Flight Mechanics Meeting, Tampa, FL,
C 23-26 Jan 2006, (Paper AAS 06-237).
C
C 3. Bowman, Bruce R.: "The Semiannual Thermospheric Density
C Variation From 1970 to 2002 Between 200-1100 km",
C AAS/AIAA Space Flight Mechanics Meeting, Maui, HI,
C 8-12 Feb 2004, (Paper AAS 04-174).
C
C 4. Bowman, Bruce R.; "Atmospheric Density Variations at
C 1500 km to 4000 km Height Determined from Long Term
C Orbit Perturbation Analysis", AAS/AIAA Space Flight
C Mechanics Meeting, Santa Barbara, CA, 11-14 Feb 2001,
C (Paper AAS 01-132).
C
C 5. Jacchia, Luigi G.; "New Static Models of the
C Thermosphere and Exosphere with Empirical Temperature
C Profiles", (Smithsonian Astrophysical Observatory
C Special Report 313), 6 May 1970.
C
C
C Written by: Bruce R Bowman (HQ AFSPC, Space Analysis Division),
C Feb 2006
C
C REV A - Oct 2006
C Smoothing of density corrections and scale height
C through different altitude bands in the latitude-
C local time correction subroutine DTSUB
C dTx correction replaced with dTc correction
C
C
C***********************************************************************
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION SUN(2), SAT(3), GEO(3), TEMP(2), AL10N(6), CHT(4)
DIMENSION ALN(6), ALPHA(5), AMW(6), FRAC(4), TC(4), WT(5)
C The alpha are the thermal diffusion coefficients in Eq. (6)
DATA ALPHA(1) /0.D0/
DATA ALPHA(2) /0.D0/
DATA ALPHA(3) /0.D0/
DATA ALPHA(4) /0.D0/
DATA ALPHA(5) /-0.38D0/
C AL10 is DLOG(10.0)
DATA AL10 /2.3025851D0/
C The AMW are the molecular weights in order: N2, O2, O, Ar, He & H
DATA AMW(1) /28.0134D0/
DATA AMW(2) /31.9988D0/
DATA AMW(3) /15.9994D0/
DATA AMW(4) /39.9480D0/
DATA AMW(5) /4.0026D0/
DATA AMW(6) /1.00797D0/
C AVOGAD is Avogadro's number in mks units (molecules/kmol)
DATA AVOGAD /6.02257D26/
DATA FOURPI /12.566371D0/
DATA TWOPI /6.2831853D0/
DATA PI /3.1415927D0/
DATA PIOV2 /1.5707963D0/
DATA PIOV4 /0.78539816D0/
C The FRAC are the assumed sea-level volume fractions in order:
C N2, O2, Ar, and He
DATA FRAC(1) /0.78110D0/
DATA FRAC(2) /0.20955D0/
DATA FRAC(3) /9.3400D-3/
DATA FRAC(4) /1.2890D-5/
C RSTAR is the universal gas-constant in mks units (joules/K/kmol)
DATA RSTAR /8314.32D0/
C The R# are values used to establish height step sizes in
C the regimes 90km to 105km, 105km to 500km and 500km upward.
DATA R1 /0.010D0/
DATA R2 /0.025D0/
DATA R3 /0.075D0/
C The WT are weights for the Newton-Cotes Five-Point Quad. formula
DATA WT(1) /0.311111111111111D0/
DATA WT(2) /1.422222222222222D0/
DATA WT(3) /0.533333333333333D0/
DATA WT(4) /1.422222222222222D0/
DATA WT(5) /0.311111111111111D0/
C The CHT are coefficients for high altitude density correction
DATA CHT /0.22D0,-0.20D-02,0.115D-02,-0.211D-05/
C
DEGRAD = PI / 180.D0
C
IDEBUG = 0
C
IF (IDEBUG.EQ.1) THEN
WRITE(8,101) (SAT(I),I=1,3),(GEO(I),I=1,3)
101 FORMAT(' JB2006 - SAT,GEO ',3F20.6,3F8.1)
ENDIF
C
F10B = GEO(2)
AP = GEO(3)
C
C Equation (14)
C
TSUBC = 379.D0 + 3.353D0*GEO(2) + 0.358D0*(GEO(1)-GEO(2))
* + 2.094D0*(S10-S10B) + 0.343D0*(XM10-XM10B)
C Equation (15)
ETA = 0.5D0 * DABS(SAT(2) - SUN(2))
THETA = 0.5D0 * DABS(SAT(2) + SUN(2))
C Equation (16)
H = SAT(1) - SUN(1)
TAU = H - 0.64577182D0 + 0.10471976D0 * DSIN(H + 0.75049158D0)
C
GLAT = SAT(2)
ZHT = SAT(3)
Z200 = ZHT/200.D0
SOLAT = SUN(2)
GLST = H + PI
GLSTHR = (GLST/DEGRAD)*(24.D0/360.D0)
IF (GLSTHR.GE.24.D0) GLSTHR = GLSTHR - 24.D0
IF (GLSTHR.LT. 0.D0) GLSTHR = GLSTHR + 24.D0
C
IF (IDEBUG.EQ.1) THEN
WRITE(8,131) ETA,THETA,TAU
131 FORMAT(' JB2006 - ETA,THETA,TAU',3D25.12)
ENDIF
C
C Equation (17)
C = DCOS(ETA)**2.5
S = DSIN(THETA)**2.5
C
IF (IDEBUG.EQ.1) THEN
WRITE(8,132) C,S
132 FORMAT(' JB2006 - C,S',2D25.12)
ENDIF
C
DF = S + (C - S) * DABS(DCOS(0.5 * TAU))**3
TSUBL = TSUBC * (1.D0 + 0.31D0 * DF)
C Equation (18)
EXPAP = DEXP(-0.08D0 * AP)
DTG = AP + 100.D0 * (1.D0 - EXPAP)
C
C Compute correction to dTc for local solar time and lat correction
C
F10 = GEO(1)
CALL DTSUB (F10,GLSTHR,GLAT,ZHT,DTCLST)
C
IF (IDEBUG.EQ.2) THEN
WRITE(8,133) SAT(3),DTCLST
133 FORMAT(' HT,DTCLST ',2F10.1)
ENDIF
C
C Compute the local exospheric temperature.
TINF = TSUBL + DTG + DTCLST
TEMP(1) = TINF
C Equation (9)
TSUBX = 444.3807D0 + 0.02385D0 * TINF - 392.8292D0 *
> DEXP(-0.0021357D0 * TINF)
C
IF (IDEBUG.EQ.1) THEN
WRITE(8,142) TINF,TSUBX,CAPPHI
142 FORMAT(' JB2006 - TINF,TSUBX,CAPPHI',3D20.8)
ENDIF
C
C Equation (11)
GSUBX = 0.054285714D0 * (TSUBX - 183.D0)
C The TC array will be an argument in the call to
C XLOCAL, which evaluates Equation (10) or Equation (13)
TC(1) = TSUBX
TC(2) = GSUBX
C A AND GSUBX/A OF Equation (13)
TC(3) = (TINF - TSUBX)/PIOV2
TC(4) = GSUBX/TC(3)
C
IF (IDEBUG.EQ.1) THEN
WRITE(8,123) TC
123 FORMAT(' JB2006 - TC',4F20.6)
ENDIF
C
C Equation (5)
Z1 = 90.D0
Z2 = DMIN1(SAT(3),105.D0)
AL = DLOG(Z2/Z1)
N = IDINT(AL/R1) + 1
ZR = DEXP(AL/DFLOAT(N))
AMBAR1 = XAMBAR(Z1)
TLOC1 = XLOCAL(Z1,TC)
ZEND = Z1
SUM2 = 0.D0
AIN = AMBAR1 * XGRAV(Z1)/TLOC1
C
IF (IDEBUG.EQ.1) THEN
WRITE(8,163) ZR,AMBAR1,TLOC1,AIN
163 FORMAT(' JB2006 - ZR-AIN',4D20.8)
ENDIF
C
DO 2 I = 1,N
Z = ZEND
ZEND = ZR * Z
DZ = 0.25D0 * (ZEND-Z)
SUM1 = WT(1)*AIN
DO 1 J = 2,5
Z = Z + DZ
AMBAR2 = XAMBAR(Z)
TLOC2 = XLOCAL(Z,TC)
GRAVL = XGRAV(Z)
AIN = AMBAR2 * GRAVL/TLOC2
C
C IF (IDEBUG.EQ.1) THEN
C WRITE(8,168) I,J,SUM1,TLOC2
C168 FORMAT(' JB2006 - SUM1,TLOC2',2I4,2D20.8)
C ENDIF
C
1 SUM1 = SUM1 + WT(J) * AIN
2 SUM2 = SUM2 + DZ * SUM1
FACT1 = 1000.0/RSTAR
RHO = 3.46D-6 * AMBAR2 * TLOC1 * DEXP(-FACT1*SUM2) /AMBAR1 /TLOC2
C Equation (2)
ANM = AVOGAD * RHO
AN = ANM/AMBAR2
C Equation (3)
FACT2 = ANM/28.960D0
ALN(1) = DLOG(FRAC(1)*FACT2)
ALN(4) = DLOG(FRAC(3)*FACT2)
ALN(5) = DLOG(FRAC(4)*FACT2)
C Equation (4)
ALN(2) = DLOG(FACT2 * (1.D0 + FRAC(2)) - AN)
ALN(3) = DLOG(2.D0 * (AN - FACT2))
C
IF (IDEBUG.EQ.1) THEN
WRITE(8,143) ALN
143 FORMAT(' JB2006 - ALN',6D20.8)
ENDIF
C
IF (SAT(3) .GT. 105.D0) GO TO 3
TEMP(2) = TLOC2
C Put in negligible hydrogen for use in DO-LOOP 13
ALN(6) = ALN(5) - 25.D0
GO TO 11
C Equation (6)
3 Z3 = DMIN1(SAT(3),500.D0)
AL = DLOG(Z3/Z)
N = IDINT(AL/R2) + 1
ZR = DEXP(AL/DFLOAT(N))
SUM2 = 0.D0
AIN = GRAVL/TLOC2
DO 5 I = 1,N
Z = ZEND
ZEND = ZR * Z
DZ = 0.25D0 * (ZEND - Z)
SUM1 = WT(1) * AIN
DO 4 J = 2,5
Z = Z + DZ
TLOC3 = XLOCAL(Z,TC)
GRAVL = XGRAV(Z)
AIN = GRAVL/TLOC3
4 SUM1 = SUM1 + WT(J) * AIN
5 SUM2 = SUM2 + DZ * SUM1
C
IF (IDEBUG.EQ.1) THEN
WRITE(8,124) SUM2,TLOC1,TLOC2,TLOC3
124 FORMAT(' JB2006 - SUM2,TLOC..',4D20.6)
ENDIF
C
Z4 = DMAX1(SAT(3),500.D0)
AL = DLOG(Z4/Z)
R = R2
IF (SAT(3) .GT. 500.D0) R = R3
N = IDINT(AL/R) + 1
ZR = DEXP(AL/DFLOAT(N))
SUM3 = 0.D0
DO 7 I=1,N
Z = ZEND
ZEND = ZR * Z
DZ = 0.25D0 * (ZEND - Z)
SUM1 = WT(1) * AIN
DO 6 J = 2,5
Z = Z + DZ
TLOC4 = XLOCAL(Z,TC)
GRAVL = XGRAV(Z)
AIN = GRAVL/TLOC4
6 SUM1 = SUM1 + WT(J) * AIN
7 SUM3 = SUM3 + DZ * SUM1
C
IF(IDEBUG.EQ.1) THEN
WRITE(8,106) SUM3,TLOC4
106 FORMAT(' SUM3,TLOC4 ',D20.9,F16.5)
ENDIF
C
IF (SAT(3) .GT. 500.D0) GO TO 8
T500 = TLOC4
TEMP(2) = TLOC3
ALTR = DLOG(TLOC3/TLOC2)
FACT2 = FACT1 * SUM2
HSIGN = 1.D0
GO TO 9
8 T500 = TLOC3
TEMP(2) = TLOC4
ALTR = DLOG(TLOC4/TLOC2)
FACT2 = FACT1 * (SUM2 + SUM3)
HSIGN = -1.D0
9 DO 10 I = 1,5
10 ALN(I) = ALN(I) - (1.0 + ALPHA(I)) * ALTR - FACT2 * AMW(I)
C Equation (7) - Note that in CIRA72, AL10T5 = DLOG10(T500)
C
IF(IDEBUG.EQ.1) THEN
WRITE(8,149) ALTR,FACT2
149 FORMAT(' JB2006 ALTR,FACT2 ',2D20.8)
ENDIF
C
C
AL10T5 = DLOG10(TINF)
ALNH5 = (5.5D0 * AL10T5 - 39.40D0) * AL10T5 + 73.13D0
ALN(6) = AL10 * (ALNH5 + 6.D0) + HSIGN * (DLOG(TLOC4/TLOC3)
> + FACT1 * SUM3 * AMW(6))
C
IF(IDEBUG.EQ.1) THEN
WRITE(8,146) ALN
146 FORMAT(' JB2006 ALN ',6D20.8)
ENDIF
C
11 CONTINUE
C Equation (24) - J70 Seasonal-Latitudinal Variation
C
TRASH = (AMJD - 36204.D0) / 365.2422D0
CAPPHI = DMOD(TRASH,1.D0)
C
DLRSL = 0.02D0 * (SAT(3) - 90.D0)
> * DEXP(-0.045D0 * (SAT(3) - 90.D0))
> * DSIGN(1.D0,SAT(2)) * DSIN(TWOPI * CAPPHI+ 1.72D0)
> * DSIN(SAT(2))**2
C
C Equation (23) - Computes the semiannual variation
C
IF (Z.LT.2000.) THEN
D1950 = AMJD - 33281.D0
CALL TMOUTD (D1950,IYR,YRDAY)
C Use new semiannual model
CALL SEMIAN (YRDAY,ZHT,F10B,FZZ,GTZ,DLRSA)
IF (FZZ.LT.0.0D0) DLRSA = 0.D0
ENDIF
C
C
IF(IDEBUG.EQ.1) THEN
WRITE(8,147) DLRSL, DLRSA
147 FORMAT(' JB2006 DLRSL',2D20.8)
ENDIF
C
C
C Sum the delta-log-rhos and apply to the number densities.
C In CIRA72 the following equation contains an actual sum,
C namely DLR = AL10 * (DLRGM + DLRSA + DLRSL)
C However, for Jacchia 70, there is no DLRGM or DLRSA.
DLR = AL10 * (DLRSL + DLRSA)
DO 12 I = 1,6
12 ALN(I) = ALN(I) + DLR
C Compute mass-density and mean-molecular-weight and
C convert number density logs from natural to common.
SUMN = 0.D0
SUMNM = 0.D0
DO 13 I = 1,6
AN = DEXP(ALN(I))
SUMN = SUMN + AN
SUMNM = SUMNM + AN*AMW(I)
AL10N(I) = ALN(I)/AL10
IF(IDEBUG.EQ.1) THEN
XVAL = 10.D0**AL10N(I)
WRITE (8,100) I,SUMNM,XVAL,ALN(I)
100 FORMAT(' JB2006 I,SUMN,XVAL,ALN ',I5,3D20.9)
ENDIF
13 CONTINUE
RHO = SUMNM/AVOGAD
C Compute the high altitude exospheric density correction factor
FEX = 1.D0
C
IF ((ZHT.GE.1000.D0).AND.(ZHT.LT.1500.D0)) THEN
ZETA = (ZHT - 1000.D0) * 0.002D0
ZETA2 = ZETA * ZETA
ZETA3 = ZETA * ZETA2
F15C = CHT(1) + CHT(2)*F10B + CHT(3)*1500.D0
* + CHT(4)*F10B*1500.D0
F15C_ZETA = (CHT(3) + CHT(4)*F10B) * 500.D0
FEX2 = 3.D0 * F15C - F15C_ZETA - 3.D0
FEX3 = F15C_ZETA - 2.D0 * F15C + 2.D0
FEX = 1.D0 + FEX2 * ZETA2 + FEX3 * ZETA3
ENDIF
IF (ZHT .GE. 1500.D0) THEN
FEX = CHT(1) + CHT(2)*F10B + CHT(3)*ZHT + CHT(4)*F10B*ZHT
ENDIF
C Apply the exospheric density correction factor.
RHO = FEX * RHO
C
IF(IDEBUG.EQ.1) THEN
WRITE (8,160) RHO
160 FORMAT(' JB2006 RHO ',D20.8,F10.4)
ENDIF
C
RETURN
END
C***********************************************************************
DOUBLE PRECISION FUNCTION XAMBAR(Z)
C Evaluates Equation (1)
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION C(7)
DATA C(1) /28.15204D0/
DATA C(2) /-8.5586D-2/
DATA C(3) /+1.2840D-4/
DATA C(4) /-1.0056D-5/
DATA C(5) /-1.0210D-5/
DATA C(6) /+1.5044D-6/
DATA C(7) /+9.9826D-8/
DZ = Z - 100.D0
AMB = C(7)
DO 1 I = 1,6
J = 7-I
1 AMB = DZ * AMB + C(J)
XAMBAR = AMB
RETURN
END
C***********************************************************************
DOUBLE PRECISION FUNCTION XGRAV(Z)
C Evaluates Equation (8)
IMPLICIT REAL*8(A-H,O-Z)
XGRAV = 9.80665D0/(1.D0 + Z/6356.766D0)**2
RETURN
END
C***********************************************************************
DOUBLE PRECISION FUNCTION XLOCAL(Z,TC)
C Evaluates Equation (10) or Equation (13), depending on Z
IMPLICIT REAL*8(A-H,O-Z)
DIMENSION TC(4)
DZ = Z - 125.D0
IF (DZ .GT. 0.D0) GO TO 1
XLOCAL = ((-9.8204695D-6 * DZ - 7.3039742D-4) * DZ**2 + 1.0)
> * DZ * TC(2) + TC(1)
RETURN
1 XLOCAL = TC(1) + TC(3) * DATAN(TC(4)*DZ*(1.D0 + 4.5D-6*DZ**2.5))
RETURN
END
C***********************************************************************
SUBROUTINE DTSUB (F10,XLST,XLAT,ZHT,DTC)
C
C COMPUTE dTc correction for Jacchia-Bowman model
C
C Calling Args:
C ------------
C F10 = (I) F10 FLUX
C XLST = (I) LOCAL SOLAR TIME (HOURS 0-23.999)
C XLAT = (I) XLAT = SAT LAT (RAD)
C ZHT = (I) ZHT = HEIGHT (KM)
C DTC = (O) dTc correction
C
IMPLICIT REAL*8(A-H,O-Z)
C
DIMENSION B(19),C(23)
C
DATA B /
C
1 -0.457512297D+01, -0.512114909D+01, -0.693003609D+02,
2 0.203716701D+03, 0.703316291D+03, -0.194349234D+04,
3 0.110651308D+04, -0.174378996D+03, 0.188594601D+04,
4 -0.709371517D+04, 0.922454523D+04, -0.384508073D+04,
5 -0.645841789D+01, 0.409703319D+02, -0.482006560D+03,
6 0.181870931D+04, -0.237389204D+04, 0.996703815D+03,
7 0.361416936D+02/
C
DATA C /
C
1 -0.155986211D+02, -0.512114909D+01, -0.693003609D+02,
2 0.203716701D+03, 0.703316291D+03, -0.194349234D+04,
3 0.110651308D+04, -0.220835117D+03, 0.143256989D+04,
4 -0.318481844D+04, 0.328981513D+04, -0.135332119D+04,
5 0.199956489D+02, -0.127093998D+02, 0.212825156D+02,
6 -0.275555432D+01, 0.110234982D+02, 0.148881951D+03,
7 -0.751640284D+03, 0.637876542D+03, 0.127093998D+02,
8 -0.212825156D+02, 0.275555432D+01/
C
C
DTC = 0.D0
C
tx = XLST/24.D0
ycs = DCOS(XLAT)
F = (F10 - 100.D0)/100.D0
C
C calculates dTc
C
IF (ZHT.GE.120.0.AND.ZHT.LE.200.0) THEN
H = (ZHT - 200.D0)/50.D0
DTC200 =
1 + C(17) + C(18)*tx*ycs + C(19)*tx**2*ycs
2 + C(20)*tx**3*ycs + C(21)*F*ycs + C(22)*tx*F*ycs
3 + C(23)*tx**2*F*ycs
sum = C(1) + B(2)*F + C(3)*tx*F + C(4)*tx**2*F
1 + C(5)*tx**3*F + C(6)*tx**4*F + C(7)*tx**5*F
2 + C(8)*tx*ycs + C(9)*tx**2*ycs + C(10)*tx**3*ycs
3 + C(11)*tx**4*ycs + C(12)*tx**5*ycs + C(13)*ycs
4 + C(14)*F*ycs + C(15)*tx*F*ycs + C(16)*tx**2*F*ycs
DTC200DZ = SUM
CC = 3.*DTC200 - DTC200DZ
DD = DTC200 - CC
ZP = (ZHT-120.)/80.D0
DTC = CC*ZP*ZP + DD*ZP*ZP*ZP
ENDIF
C
C
IF (ZHT.GT.200.0.AND.ZHT.LE.240.0) THEN
H = (ZHT - 200.D0)/50.D0
sum = C(1)*H + B(2)*F*H + C(3)*tx*F*H + C(4)*tx**2*F*H
1 + C(5)*tx**3*F*H + C(6)*tx**4*F*H + C(7)*tx**5*F*H
2 + C(8)*tx*ycs*H + C(9)*tx**2*ycs*H + C(10)*tx**3*ycs*H
3 + C(11)*tx**4*ycs*H + C(12)*tx**5*ycs*H + C(13)*ycs*H
4 + C(14)*F*ycs*H + C(15)*tx*F*ycs*H + C(16)*tx**2*F*ycs*H
5 + C(17) + C(18)*tx*ycs + C(19)*tx**2*ycs
6 + C(20)*tx**3*ycs + C(21)*F*ycs + C(22)*tx*F*ycs
7 + C(23)*tx**2*F*ycs
DTC = sum
ENDIF
C
C
IF (ZHT.GT.240.0.AND.ZHT.LE.300.0) THEN
H = (40.D0)/50.D0
sum = C(1)*H + B(2)*F*H + C(3)*tx*F*H + C(4)*tx**2*F*H
1 + C(5)*tx**3*F*H + C(6)*tx**4*F*H + C(7)*tx**5*F*H
2 + C(8)*tx*ycs*H + C(9)*tx**2*ycs*H + C(10)*tx**3*ycs*H
3 + C(11)*tx**4*ycs*H + C(12)*tx**5*ycs*H + C(13)*ycs*H
4 + C(14)*F*ycs*H + C(15)*tx*F*ycs*H + C(16)*tx**2*F*ycs*H
5 + C(17) + C(18)*tx*ycs + C(19)*tx**2*ycs
6 + C(20)*tx**3*ycs + C(21)*F*ycs + C(22)*tx*F*ycs
7 + C(23)*tx**2*F*ycs
AA = SUM
BB = C(1) + B(2)*F + C(3)*tx*F + C(4)*tx**2*F
1 + C(5)*tx**3*F + C(6)*tx**4*F + C(7)*tx**5*F
2 + C(8)*tx*ycs + C(9)*tx**2*ycs + C(10)*tx**3*ycs
3 + C(11)*tx**4*ycs + C(12)*tx**5*ycs + C(13)*ycs
4 + C(14)*F*ycs + C(15)*tx*F*ycs + C(16)*tx**2*F*ycs
H = 300./100.D0
sum = B(1) + B(2)*F + B(3)*tx*F + B(4)*tx**2*F
1 + B(5)*tx**3*F + B(6)*tx**4*F + B(7)*tx**5*F
2 + B(8)*tx*ycs + B(9)*tx**2*ycs + B(10)*tx**3*ycs
3 + B(11)*tx**4*ycs + B(12)*tx**5*ycs + B(13)*H*ycs
4 + B(14)*tx*H*ycs + B(15)*tx**2*H*ycs + B(16)*tx**3*H*ycs
5 + B(17)*tx**4*H*ycs + B(18)*tx**5*H*ycs + B(19)*ycs
DTC300 = SUM
sum = B(13)*ycs
1 + B(14)*tx*ycs + B(15)*tx**2*ycs + B(16)*tx**3*ycs
2 + B(17)*tx**4*ycs + B(18)*tx**5*ycs
DTC300DZ = SUM
CC = 3.*DTC300 - DTC300DZ - 3.*AA - 2.*BB
DD = DTC300 - AA - BB - CC
ZP = (ZHT-240.)/60.D0
DTC = AA + BB*ZP + CC*ZP*ZP + DD*ZP*ZP*ZP
ENDIF
C
C
IF (ZHT.GT.300.0.AND.ZHT.LE.600.0) THEN
H = ZHT/100.D0
sum = B(1) + B(2)*F + B(3)*tx*F + B(4)*tx**2*F
1 + B(5)*tx**3*F + B(6)*tx**4*F + B(7)*tx**5*F
2 + B(8)*tx*ycs + B(9)*tx**2*ycs + B(10)*tx**3*ycs
3 + B(11)*tx**4*ycs + B(12)*tx**5*ycs + B(13)*H*ycs
4 + B(14)*tx*H*ycs + B(15)*tx**2*H*ycs + B(16)*tx**3*H*ycs
5 + B(17)*tx**4*H*ycs + B(18)*tx**5*H*ycs + B(19)*ycs
DTC = sum
ENDIF
C
C
IF (ZHT.GT.600.0.AND.ZHT.LE.800.0) THEN
ZP = (ZHT - 600.D0)/100.
HP = 600.D0/100.D0
AA = B(1) + B(2)*F + B(3)*tx*F + B(4)*tx**2*F
1 + B(5)*tx**3*F + B(6)*tx**4*F + B(7)*tx**5*F
2 + B(8)*tx*ycs + B(9)*tx**2*ycs + B(10)*tx**3*ycs
3 + B(11)*tx**4*ycs + B(12)*tx**5*ycs + B(13)*HP*ycs
4 + B(14)*tx*HP*ycs + B(15)*tx**2*HP*ycs+ B(16)*tx**3*HP*ycs
5 + B(17)*tx**4*HP*ycs + B(18)*tx**5*HP*ycs + B(19)*ycs
BB = B(13)*ycs
1 + B(14)*tx*ycs + B(15)*tx**2*ycs + B(16)*tx**3*ycs
2 + B(17)*tx**4*ycs + B(18)*tx**5*ycs
CC = -(3.*AA+4.*BB)/4.
DD = (AA+BB)/4.
DTC = AA + BB*ZP + CC*ZP*ZP + DD*ZP*ZP*ZP
ENDIF
C
RETURN
END SUBROUTINE DTSUB
C
C******************************************************************
C
SUBROUTINE SEMIAN (DAY,HT,F10BAR,FZZ,GTZ,DRLOG)
C
C COMPUTE SEMIANNUAL VARIATION (DELTA LOG RHO)
C INPUT DAY, HEIGHT, F10BAR
C 025. 650. 150.
C OUTPUT FUNCTIONS FZ, GT, AND DEL LOG RHO VALUE
C
C DAY (I) DAY OF YEAR
C HT (I) HEIGHT (KM)
C F10BAR (I) AVE 81-DAY CENTERED F10
C FZZ (O) SEMIANNUAL AMPLITUDE
C GTZ (O) SEMIANNUAL PHASE FUNCTION
C DRLOG (O) DELTA LOG RHO
C
IMPLICIT REAL*8(A-H,O-Z)
C
REAL*8 FZM(6),GTM(23)
C
DATA TWOPI/6.2831853072D0/
C
C FZ GLOBAL MODEL VALUES
C 1978-2004 FIT:
DATA FZM /
1 0.111613D+00,-0.159000D-02, 0.126190D-01,
2 -0.100064D-01,-0.237509D-04, 0.260759D-04/
C
C GT GLOBAL MODEL VALUES
C 1978-2004 FIT:
DATA GTM /
1 -0.833646D+00,-0.265450D+00, 0.467603D+00,-0.299906D+00,
2 -0.105451D+00,-0.165537D-01,-0.380037D-01,-0.150991D-01,
3 -0.541280D-01, 0.119554D-01, 0.437544D-02,-0.369016D-02,
4 0.206763D-02,-0.142888D-02,-0.867124D-05, 0.189032D-04,
5 0.156988D-03, 0.491286D-03,-0.391484D-04,-0.126854D-04,
6 0.134078D-04,-0.614176D-05, 0.343423D-05/
C
C
F10B = F10BAR
FB2 = F10BAR**2
C
HTZ = HT/1000.D0
C
FZZ = FZM(1) + FZM(2)*F10B + FZM(3)*F10B*HTZ
* + FZM(4)*F10B*HTZ**2 + FZM(5)*F10B**2*HTZ
* + FZM(6)*F10B**2*HTZ**2
C
TAU = (DAY-1.D0)/365.D0
SIN1P = DSIN(TWOPI*TAU)
COS1P = DCOS(TWOPI*TAU)
SIN2P = DSIN(2.*TWOPI*TAU)
COS2P = DCOS(2.*TWOPI*TAU)
SIN3P = DSIN(3.*TWOPI*TAU)
COS3P = DCOS(3.*TWOPI*TAU)
SIN4P = DSIN(4.*TWOPI*TAU)
COS4P = DCOS(4.*TWOPI*TAU)
C
GTZ = GTM(1) + GTM(2)*SIN1P + GTM(3)*COS1P
1 + GTM(4)*SIN2P + GTM(5)*COS2P
2 + GTM(6)*SIN3P + GTM(7)*COS3P
3 + GTM(8)*SIN4P + GTM(9)*COS4P
4 + GTM(10)*F10B + GTM(11)*F10B*SIN1P + GTM(12)*F10B*COS1P
5 + GTM(13)*F10B*SIN2P + GTM(14)*F10B*COS2P
6 + GTM(15)*F10B*SIN3P + GTM(16)*F10B*COS3P
7 + GTM(17)*F10B*SIN4P + GTM(18)*F10B*COS4P
8 + GTM(19)*FB2 + GTM(20)*FB2 *SIN1P + GTM(21)*FB2 *COS1P
9 + GTM(22)*FB2 *SIN2P + GTM(23)*FB2 *COS2P
C
IF (FZZ.LT.1.D-6) FZZ = 1.D-6
C
DRLOG = FZZ*GTZ
C
RETURN
END SUBROUTINE SEMIAN
C
C******************************************************************
C
SUBROUTINE TMOUTD(D1950,IYR,DAY)
C
C COMPUTE DAY AND YEAR FROM TIME D1950 (DAYS SINCE 1950)
C
IMPLICIT REAL*8(A-H,O-Z)
C
IYDAY = D1950
FRAC = D1950 - IYDAY
IYDAY = IYDAY + 364
ITEMP = IYDAY/1461
IYDAY = IYDAY - ITEMP*1461
IYR = 1949 + 4*ITEMP
ITEMP = IYDAY/365
IF (ITEMP.GE.3) ITEMP = 3
IYR = IYR + ITEMP
IYDAY = IYDAY - 365*ITEMP + 1
IYR = IYR - 1900
DAY = IYDAY + FRAC
C
RETURN
END SUBROUTINE TMOUTD
C


[Home]
[Intro] [Indices] [Code]
[Publications]
[Contacts]
[Figures]
[License]
[SET/Spacewx]