module msise00_glow

C module to avoid link conflicts

use, intrinsic :: iso_fortran_env, only : stderr=>error_unit

use msise00_data_glow, only : parm7, ptm, pdm, pavgm, imr

private

public :: gtd7, tselec, meters

contains

SUBROUTINE GTD7(IYD,SEC,ALT,GLAT,GLONG,STL,F107A,F107,AP,MASS,D,T)

C

C NRLMSISE-00

C -----------

C Neutral Atmosphere Empirical Model from the surface to lower

C exosphere

C

C NEW FEATURES:

C *Extensive satellite drag database used in model generation

C *Revised O2 (and O) in lower thermosphere

C *Additional nonlinear solar activity term

C *"ANOMALOUS OXYGEN" NUMBER DENSITY, OUTPUT D(9)

C At high altitudes (> 500 km), hot atomic oxygen or ionized

C oxygen can become appreciable for some ranges of subroutine

C inputs, thereby affecting drag on satellites and debris. We

C group these species under the term "anomalous oxygen," since

C their individual variations are not presently separable with

C the drag data used to define this model component.

C

C SUBROUTINES FOR SPECIAL OUTPUTS:

C

C HIGH ALTITUDE DRAG: EFFECTIVE TOTAL MASS DENSITY

C (SUBROUTINE GTD7D, OUTPUT D(6))

C For atmospheric drag calculations at altitudes above 500 km,

C call SUBROUTINE GTD7D to compute the "effective total mass

C density" by including contributions from "anomalous oxygen."

C See "NOTES ON OUTPUT VARIABLES" below on D(6).

C

C PRESSURE GRID (SUBROUTINE GHP7)

C See subroutine GHP7 to specify outputs at a pressure level

C rather than at an altitude.

C

C OUTPUT IN M-3 and KG/M3: CALL METERS(.TRUE.)

C

C INPUT VARIABLES:

C IYD - YEAR AND DAY AS YYDDD (day of year from 1 to 365 (or 366))

C (Year ignored in current model)

C SEC - UT(SEC)

C ALT - ALTITUDE(KM)

C GLAT - GEODETIC LATITUDE(DEG)

C GLONG - GEODETIC LONGITUDE(DEG)

C STL - LOCAL APPARENT SOLAR TIME(HRS; see Note below)

C F107A - 81 day AVERAGE OF F10.7 FLUX (centered on day DDD)

C F107 - DAILY F10.7 FLUX FOR PREVIOUS DAY

C AP - MAGNETIC INDEX(DAILY) OR WHEN SW(9)=-1. :

C - ARRAY CONTAINING:

C (1) DAILY AP

C (2) 3 HR AP INDEX FOR CURRENT TIME

C (3) 3 HR AP INDEX FOR 3 HRS BEFORE CURRENT TIME

C (4) 3 HR AP INDEX FOR 6 HRS BEFORE CURRENT TIME

C (5) 3 HR AP INDEX FOR 9 HRS BEFORE CURRENT TIME

C (6) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 12 TO 33 HRS PRIOR

C TO CURRENT TIME

C (7) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 36 TO 57 HRS PRIOR

C TO CURRENT TIME

C MASS - MASS NUMBER (ONLY DENSITY FOR SELECTED GAS IS

C CALCULATED. MASS 0 IS TEMPERATURE. MASS 48 FOR ALL.

C MASS 17 IS Anomalous O ONLY.)

C

C NOTES ON INPUT VARIABLES:

C UT, Local Time, and Longitude are used independently in the

C model and are not of equal importance for every situation.

C For the most physically realistic calculation these three

C variables should be consistent (STL=SEC/3600+GLONG/15).

C The Equation of Time departures from the above formula

C for apparent local time can be included if available but

C are of minor importance.

c

C F107 and F107A values used to generate the model correspond

C to the 10.7 cm radio flux at the actual distance of the Earth

C from the Sun rather than the radio flux at 1 AU. The following

C site provides both classes of values:

C ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/

C

C F107, F107A, and AP effects are neither large nor well

C established below 80 km and these parameters should be set to

C 150., 150., and 4. respectively.

C

C OUTPUT VARIABLES:

C D(1) - HE NUMBER DENSITY(CM-3)

C D(2) - O NUMBER DENSITY(CM-3)

C D(3) - N2 NUMBER DENSITY(CM-3)

C D(4) - O2 NUMBER DENSITY(CM-3)

C D(5) - AR NUMBER DENSITY(CM-3)

C D(6) - TOTAL MASS DENSITY(GM/CM3)

C D(7) - H NUMBER DENSITY(CM-3)

C D(8) - N NUMBER DENSITY(CM-3)

C D(9) - Anomalous oxygen NUMBER DENSITY(CM-3)

C T(1) - EXOSPHERIC TEMPERATURE

C T(2) - TEMPERATURE AT ALT

C

C NOTES ON OUTPUT VARIABLES:

C TO GET OUTPUT IN M-3 and KG/M3: CALL METERS(.TRUE.)

C

C O, H, and N are set to zero below 72.5 km

C

C T(1), Exospheric temperature, is set to global average for

C altitudes below 120 km. The 120 km gradient is left at global

C average value for altitudes below 72 km.

C

C D(6), TOTAL MASS DENSITY, is NOT the same for subroutines GTD7

C and GTD7D

C

C SUBROUTINE GTD7 -- D(6) is the sum of the mass densities of the

C species labeled by indices 1-5 and 7-8 in output variable D.

C This includes He, O, N2, O2, Ar, H, and N but does NOT include

C anomalous oxygen (species index 9).

C

C SUBROUTINE GTD7D -- D(6) is the "effective total mass density

C for drag" and is the sum of the mass densities of all species

C in this model, INCLUDING anomalous oxygen.

C

C SWITCHES: The following is for test and special purposes:

C

C TO TURN ON AND OFF PARTICULAR VARIATIONS CALL TSELEC(SW),

C WHERE SW IS A 25 ELEMENT ARRAY CONTAINING 0. FOR OFF, 1.

C FOR ON, OR 2. FOR MAIN EFFECTS OFF BUT CROSS TERMS ON

C FOR THE FOLLOWING VARIATIONS

C 1 - F10.7 EFFECT ON MEAN 2 - TIME INDEPENDENT

C 3 - SYMMETRICAL ANNUAL 4 - SYMMETRICAL SEMIANNUAL

C 5 - ASYMMETRICAL ANNUAL 6 - ASYMMETRICAL SEMIANNUAL

C 7 - DIURNAL 8 - SEMIDIURNAL

C 9 - DAILY AP 10 - ALL UT/LONG EFFECTS

C 11 - LONGITUDINAL 12 - UT AND MIXED UT/LONG

C 13 - MIXED AP/UT/LONG 14 - TERDIURNAL

C 15 - DEPARTURES FROM DIFFUSIVE EQUILIBRIUM

C 16 - ALL TINF VAR 17 - ALL TLB VAR

C 18 - ALL TN1 VAR 19 - ALL S VAR

C 20 - ALL TN2 VAR 21 - ALL NLB VAR

C 22 - ALL TN3 VAR 23 - TURBO SCALE HEIGHT VAR

C

C To get current values of SW: CALL TRETRV(SW)

C

Real,Intent(OUT):: D(9), T(2)

Real,Intent(In) :: SEC,ALT,GLAT,GLONG,STL,F107A,F107,AP(7)

Integer,Intent(IN)::IYD,MASS

Real DS(9),TS(2)

real ZN3(5),ZN2(4),SV(25)

COMMON/GTS3C/TLB,S,DB04,DB16,DB28,DB32,DB40,DB48,DB01,ZA,T0,Z0

& ,G0,RL,DD,DB14,TR12

COMMON/MESO7/TN1(5),TN2(4),TN3(5),TGN1(2),TGN2(2),TGN3(2)

COMMON/PARM7/PT(150),PD(150,9),PS(150),PDL(25,2),PTL(100,4),

$ PMA(100,10),SAM(100)

COMMON/CSW/SW(25),ISW,SWC(25)

COMMON/DMIX/DM04,DM16,DM28,DM32,DM40,DM01,DM14

COMMON/PARMB/GSURF,RE

COMMON/METSEL/IMR

SAVE

DATA MN3/5/,ZN3/32.5,20.,15.,10.,0./

DATA MN2/4/,ZN2/72.5,55.,45.,32.5/

DATA ZMIX/62.5/,ALAST/99999./,MSSL/-999/

DATA SV/25*1./

IF(ISW.NE.64999) CALL TSELEC(SV)

C

C Test for changed input

V1=VTST7(IYD,SEC,GLAT,GLONG,STL,F107A,F107,AP,1)

C Latitude variation of gravity (none for SW(2)=0)

XLAT=GLAT

IF(SW(2).EQ.0) XLAT=45.

CALL GLATF(XLAT,GSURF,RE)

C

XMM=PDM(5,3)

C

C THERMOSPHERE/MESOSPHERE (above ZN2(1))

ALTT=MAX(ALT,ZN2(1))

MSS=MASS

C Only calculate N2 in thermosphere if alt in mixed region

IF(ALT.LT.ZMIX.AND.MASS.GT.0) MSS=28

C Only calculate thermosphere if input parameters changed

C or altitude above ZN2(1) in mesosphere

IF(V1.EQ.1..OR.ALT.GT.ZN2(1).OR.ALAST.GT.ZN2(1).OR.MSS.NE.MSSL)

$ THEN

CALL GTS7(IYD,SEC,ALTT,GLAT,GLONG,STL,F107A,F107,AP,MSS,DS,TS)

DM28M=DM28

C metric adjustment

IF(IMR.EQ.1) DM28M=DM28*1.E6

MSSL=MSS

ENDIF

T(1)=TS(1)

T(2)=TS(2)

IF(ALT.GE.ZN2(1)) THEN

DO 5 J=1,9

D(J)=DS(J)

5 CONTINUE

GOTO 10

ENDIF

C

C LOWER MESOSPHERE/UPPER STRATOSPHERE [between ZN3(1) and ZN2(1)]

C Temperature at nodes and gradients at end nodes

C Inverse temperature a linear function of spherical harmonics

C Only calculate nodes if input changed

IF(V1.EQ.1..OR.ALAST.GE.ZN2(1)) THEN

TGN2(1)=TGN1(2)

TN2(1)=TN1(5)

TN2(2)=PMA(1,1)*PAVGM(1)/(1.-SW(20)*GLOB7S(PMA(1,1)))

TN2(3)=PMA(1,2)*PAVGM(2)/(1.-SW(20)*GLOB7S(PMA(1,2)))

TN2(4)=PMA(1,3)*PAVGM(3)/(1.-SW(20)*SW(22)*GLOB7S(PMA(1,3)))

TGN2(2)=PAVGM(9)*PMA(1,10)*(1.+SW(20)*SW(22)*GLOB7S(PMA(1,10)))

$ *TN2(4)*TN2(4)/(PMA(1,3)*PAVGM(3))**2

TN3(1)=TN2(4)

ENDIF

C Including ZN3(1) in the jump condition creates a model coverage gap at that exact altitude

C IF(ALT.GE.ZN3(1)) GOTO 6

IF(ALT.GT.ZN3(1)) GOTO 6

C

C LOWER STRATOSPHERE AND TROPOSPHERE [below ZN3(1)]

C Temperature at nodes and gradients at end nodes

C Inverse temperature a linear function of spherical harmonics

C Only calculate nodes if input changed

IF(V1.EQ.1..OR.ALAST.GE.ZN3(1)) THEN

TGN3(1)=TGN2(2)

TN3(2)=PMA(1,4)*PAVGM(4)/(1.-SW(22)*GLOB7S(PMA(1,4)))

TN3(3)=PMA(1,5)*PAVGM(5)/(1.-SW(22)*GLOB7S(PMA(1,5)))

TN3(4)=PMA(1,6)*PAVGM(6)/(1.-SW(22)*GLOB7S(PMA(1,6)))

TN3(5)=PMA(1,7)*PAVGM(7)/(1.-SW(22)*GLOB7S(PMA(1,7)))

TGN3(2)=PMA(1,8)*PAVGM(8)*(1.+SW(22)*GLOB7S(PMA(1,8)))

$ *TN3(5)*TN3(5)/(PMA(1,7)*PAVGM(7))**2

ENDIF

6 CONTINUE

IF(MASS.EQ.0) GOTO 50

C LINEAR TRANSITION TO FULL MIXING BELOW ZN2(1)

DMC=0

IF(ALT.GT.ZMIX) DMC=1.-(ZN2(1)-ALT)/(ZN2(1)-ZMIX)

DZ28=DS(3)

C ***** N2 DENSITY ****

DMR=DS(3)/DM28M-1.

D(3)=DENSM(ALT,DM28M,XMM,TZ,MN3,ZN3,TN3,TGN3,MN2,ZN2,TN2,TGN2)

D(3)=D(3)*(1.+DMR*DMC)

C ***** HE DENSITY ****

D(1)=0

IF(MASS.NE.4.AND.MASS.NE.48) GOTO 204

DMR=DS(1)/(DZ28*PDM(2,1))-1.

D(1)=D(3)*PDM(2,1)*(1.+DMR*DMC)

204 CONTINUE

C **** O DENSITY ****

D(2)=0

D(9)=0

216 CONTINUE

C ***** O2 DENSITY ****

D(4)=0

IF(MASS.NE.32.AND.MASS.NE.48) GOTO 232

DMR=DS(4)/(DZ28*PDM(2,4))-1.

D(4)=D(3)*PDM(2,4)*(1.+DMR*DMC)

232 CONTINUE

C ***** AR DENSITY ****

D(5)=0

IF(MASS.NE.40.AND.MASS.NE.48) GOTO 240

DMR=DS(5)/(DZ28*PDM(2,5))-1.

D(5)=D(3)*PDM(2,5)*(1.+DMR*DMC)

240 CONTINUE

C ***** HYDROGEN DENSITY ****

D(7)=0

C ***** ATOMIC NITROGEN DENSITY ****

D(8)=0

C

C TOTAL MASS DENSITY

C

IF(MASS.EQ.48) THEN

D(6) = 1.66E-24*(4.*D(1)+16.*D(2)+28.*D(3)+32.*D(4)+40.*D(5)+

& D(7)+14.*D(8))

IF(IMR.EQ.1) D(6)=D(6)/1000.

ENDIF

T(2)=TZ

10 CONTINUE

GOTO 90

50 CONTINUE

DD=DENSM(ALT,1.,0.,TZ,MN3,ZN3,TN3,TGN3,MN2,ZN2,TN2,TGN2)

T(2)=TZ

90 CONTINUE

ALAST=ALT

RETURN

END

C-----------------------------------------------------------------------

SUBROUTINE GTD7D(IYD,SEC,ALT,GLAT,GLONG,STL,F107A,F107,AP,MASS,

$ D,T)

C

C NRLMSISE-00

C -----------

C This subroutine provides Effective Total Mass Density for

C output D(6) which includes contributions from "anomalous

C oxygen" which can affect satellite drag above 500 km. This

C subroutine is part of the distribution package for the

C Neutral Atmosphere Empirical Model from the surface to lower

C exosphere. See subroutine GTD7 for more extensive comments.

C

C INPUT VARIABLES:

C IYD - YEAR AND DAY AS YYDDD (day of year from 1 to 365 (or 366))

C (Year ignored in current model)

C SEC - UT(SEC)

C ALT - ALTITUDE(KM)

C GLAT - GEODETIC LATITUDE(DEG)

C GLONG - GEODETIC LONGITUDE(DEG)

C STL - LOCAL APPARENT SOLAR TIME(HRS; see Note below)

C F107A - 81 day AVERAGE OF F10.7 FLUX (centered on day DDD)

C F107 - DAILY F10.7 FLUX FOR PREVIOUS DAY

C AP - MAGNETIC INDEX(DAILY) OR WHEN SW(9)=-1. :

C - ARRAY CONTAINING:

C (1) DAILY AP

C (2) 3 HR AP INDEX FOR CURRENT TIME

C (3) 3 HR AP INDEX FOR 3 HRS BEFORE CURRENT TIME

C (4) 3 HR AP INDEX FOR 6 HRS BEFORE CURRENT TIME

C (5) 3 HR AP INDEX FOR 9 HRS BEFORE CURRENT TIME

C (6) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 12 TO 33 HRS PRIOR

C TO CURRENT TIME

C (7) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 36 TO 57 HRS PRIOR

C TO CURRENT TIME

C MASS - MASS NUMBER (ONLY DENSITY FOR SELECTED GAS IS

C CALCULATED. MASS 0 IS TEMPERATURE. MASS 48 FOR ALL.

C MASS 17 IS Anomalous O ONLY.)

C

C NOTES ON INPUT VARIABLES:

C UT, Local Time, and Longitude are used independently in the

C model and are not of equal importance for every situation.

C For the most physically realistic calculation these three

C variables should be consistent (STL=SEC/3600+GLONG/15).

C The Equation of Time departures from the above formula

C for apparent local time can be included if available but

C are of minor importance.

c

C F107 and F107A values used to generate the model correspond

C to the 10.7 cm radio flux at the actual distance of the Earth

C from the Sun rather than the radio flux at 1 AU.

C

C OUTPUT VARIABLES:

C D(1) - HE NUMBER DENSITY(CM-3)

C D(2) - O NUMBER DENSITY(CM-3)

C D(3) - N2 NUMBER DENSITY(CM-3)

C D(4) - O2 NUMBER DENSITY(CM-3)

C D(5) - AR NUMBER DENSITY(CM-3)

C D(6) - TOTAL MASS DENSITY(GM/CM3) [includes anomalous oxygen]

C D(7) - H NUMBER DENSITY(CM-3)

C D(8) - N NUMBER DENSITY(CM-3)

C D(9) - Anomalous oxygen NUMBER DENSITY(CM-3)

C T(1) - EXOSPHERIC TEMPERATURE

C T(2) - TEMPERATURE AT ALT

C

DIMENSION D(9),T(2),AP(7),DS(9),TS(2)

COMMON/METSEL/IMR

CALL GTD7(IYD,SEC,ALT,GLAT,GLONG,STL,F107A,F107,AP,MASS,D,T)

C TOTAL MASS DENSITY

C

IF(MASS.EQ.48) THEN

D(6) = 1.66E-24*(4.*D(1)+16.*D(2)+28.*D(3)+32.*D(4)+40.*D(5)+

& D(7)+14.*D(8)+16.*D(9))

IF(IMR.EQ.1) D(6)=D(6)/1000.

ENDIF

END SUBROUTINE GTD7D

C-----------------------------------------------------------------------

SUBROUTINE GHP7(IYD,SEC,ALT,GLAT,GLONG,STL,F107A,F107,AP,

$ D,T,PRESS)

C FIND ALTITUDE OF PRESSURE SURFACE (PRESS) FROM GTD7

C INPUT:

C IYD - YEAR AND DAY AS YYDDD

C SEC - UT(SEC)

C GLAT - GEODETIC LATITUDE(DEG)

C GLONG - GEODETIC LONGITUDE(DEG)

C STL - LOCAL APPARENT SOLAR TIME(HRS)

C F107A - 3 MONTH AVERAGE OF F10.7 FLUX

C F107 - DAILY F10.7 FLUX FOR PREVIOUS DAY

C AP - MAGNETIC INDEX(DAILY) OR WHEN SW(9)=-1. :

C - ARRAY CONTAINING:

C (1) DAILY AP

C (2) 3 HR AP INDEX FOR CURRENT TIME

C (3) 3 HR AP INDEX FOR 3 HRS BEFORE CURRENT TIME

C (4) 3 HR AP INDEX FOR 6 HRS BEFORE CURRENT TIME

C (5) 3 HR AP INDEX FOR 9 HRS BEFORE CURRENT TIME

C (6) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 12 TO 33 HRS PRIOR

C TO CURRENT TIME

C (7) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 36 TO 59 HRS PRIOR

C TO CURRENT TIME

C PRESS - PRESSURE LEVEL(MB)

C OUTPUT:

C ALT - ALTITUDE(KM)

C D(1) - HE NUMBER DENSITY(CM-3)

C D(2) - O NUMBER DENSITY(CM-3)

C D(3) - N2 NUMBER DENSITY(CM-3)

C D(4) - O2 NUMBER DENSITY(CM-3)

C D(5) - AR NUMBER DENSITY(CM-3)

C D(6) - TOTAL MASS DENSITY(GM/CM3)

C D(7) - H NUMBER DENSITY(CM-3)

C D(8) - N NUMBER DENSITY(CM-3)

C D(9) - HOT O NUMBER DENSITY(CM-3)

C T(1) - EXOSPHERIC TEMPERATURE

C T(2) - TEMPERATURE AT ALT

C

COMMON/PARMB/GSURF,RE

COMMON/METSEL/IMR

DIMENSION D(9),T(2),AP(7)

SAVE

DATA BM/1.3806E-19/,RGAS/831.4/

DATA TEST/.00043/,LTEST/12/

PL=ALOG10(PRESS)

C Initial altitude estimate

IF(PL.GE.-5.) THEN

IF(PL.GT.2.5) ZI=18.06*(3.00-PL)

IF(PL.GT..75.AND.PL.LE.2.5) ZI=14.98*(3.08-PL)

IF(PL.GT.-1..AND.PL.LE..75) ZI=17.8*(2.72-PL)

IF(PL.GT.-2..AND.PL.LE.-1.) ZI=14.28*(3.64-PL)

IF(PL.GT.-4..AND.PL.LE.-2.) ZI=12.72*(4.32-PL)

IF(PL.LE.-4.) ZI=25.3*(.11-PL)

IDAY=MOD(IYD,1000)

CL=GLAT/90.

CL2=CL*CL

IF(IDAY.LT.182) CD=1.-IDAY/91.25

IF(IDAY.GE.182) CD=IDAY/91.25-3.

CA=0

IF(PL.GT.-1.11.AND.PL.LE.-.23) CA=1.0

IF(PL.GT.-.23) CA=(2.79-PL)/(2.79+.23)

IF(PL.LE.-1.11.AND.PL.GT.-3.) CA=(-2.93-PL)/(-2.93+1.11)

Z=ZI-4.87*CL*CD*CA-1.64*CL2*CA+.31*CA*CL

ENDIF

IF(PL.LT.-5.) Z=22.*(PL+4.)**2+110

C ITERATION LOOP

L=0

10 CONTINUE

L=L+1

CALL GTD7(IYD,SEC,Z,GLAT,GLONG,STL,F107A,F107,AP,48,D,T)

XN=D(1)+D(2)+D(3)+D(4)+D(5)+D(7)+D(8)

P=BM*XN*T(2)

IF(IMR.EQ.1) P=P*1.E-6

DIFF=PL-ALOG10(P)

IF(ABS(DIFF).LT.TEST .OR. L.EQ.LTEST) GOTO 20

XM=D(6)/XN/1.66E-24

IF(IMR.EQ.1) XM = XM*1.E3

G=GSURF/(1.+Z/RE)**2

SH=RGAS*T(2)/(XM*G)

C New altitude estimate using scale height

IF(L.LT.6) THEN

Z=Z-SH*DIFF*2.302

ELSE

Z=Z-SH*DIFF

ENDIF

GOTO 10

20 CONTINUE

IF(L.EQ.LTEST) write(stderr,100) PRESS,DIFF

100 FORMAT(1X,29HGHP7 NOT CONVERGING FOR PRESS, 1PE12.2,E12.2)

ALT=Z

END SUBROUTINE GHP7

C-----------------------------------------------------------------------

SUBROUTINE GLATF(LAT,GV,REFF)

C CALCULATE LATITUDE VARIABLE GRAVITY (GV) AND EFFECTIVE

C RADIUS (REFF)

REAL LAT

SAVE

DATA DGTR/1.74533E-2/

C2 = COS(2.*DGTR*LAT)

GV = 980.616*(1.-.0026373*C2)

REFF = 2.*GV/(3.085462E-6 + 2.27E-9*C2)*1.E-5

END SUBROUTINE GLATF

C-----------------------------------------------------------------------

FUNCTION VTST7(IYD,SEC,GLAT,GLONG,STL,F107A,F107,AP,IC)

C Test if geophysical variables or switches changed and save

C Return 0 if unchanged and 1 if changed

DIMENSION AP(7),IYDL(2),SECL(2),GLATL(2),GLL(2),STLL(2)

DIMENSION FAL(2),FL(2),APL(7,2),SWL(25,2),SWCL(25,2)

COMMON/CSW/SW(25),ISW,SWC(25)

SAVE

DATA IYDL/2*-999/,SECL/2*-999./,GLATL/2*-999./,GLL/2*-999./

DATA STLL/2*-999./,FAL/2*-999./,FL/2*-999./,APL/14*-999./

DATA SWL/50*-999./,SWCL/50*-999./

VTST7=0

IF(IYD.NE.IYDL(IC)) GOTO 10

IF(SEC.NE.SECL(IC)) GOTO 10

IF(GLAT.NE.GLATL(IC)) GOTO 10

IF(GLONG.NE.GLL(IC)) GOTO 10

IF(STL.NE.STLL(IC)) GOTO 10

IF(F107A.NE.FAL(IC)) GOTO 10

IF(F107.NE.FL(IC)) GOTO 10

DO 5 I=1,7

IF(AP(I).NE.APL(I,IC)) GOTO 10

5 CONTINUE

DO 7 I=1,25

IF(SW(I).NE.SWL(I,IC)) GOTO 10

IF(SWC(I).NE.SWCL(I,IC)) GOTO 10

7 CONTINUE

GOTO 20

10 CONTINUE

VTST7=1

IYDL(IC)=IYD

SECL(IC)=SEC

GLATL(IC)=GLAT

GLL(IC)=GLONG

STLL(IC)=STL

FAL(IC)=F107A

FL(IC)=F107

DO 15 I=1,7

APL(I,IC)=AP(I)

15 CONTINUE

DO 16 I=1,25

SWL(I,IC)=SW(I)

SWCL(I,IC)=SWC(I)

16 CONTINUE

20 CONTINUE

END FUNCTION VTST7

C-----------------------------------------------------------------------

SUBROUTINE GTS7(IYD,SEC,ALT,GLAT,GLONG,STL,F107A,F107,AP,MASS,D,T)

C

C Thermospheric portion of NRLMSISE-00

C See GTD7 for more extensive comments

C

C OUTPUT IN M-3 and KG/M3: CALL METERS(.TRUE.)

C

C INPUT VARIABLES:

C IYD - YEAR AND DAY AS YYDDD (day of year from 1 to 365 (or 366))

C (Year ignored in current model)

C SEC - UT(SEC)

C ALT - ALTITUDE(KM) (>72.5 km)

C GLAT - GEODETIC LATITUDE(DEG)

C GLONG - GEODETIC LONGITUDE(DEG)

C STL - LOCAL APPARENT SOLAR TIME(HRS; see Note below)

C F107A - 81 day AVERAGE OF F10.7 FLUX (centered on day DDD)

C F107 - DAILY F10.7 FLUX FOR PREVIOUS DAY

C AP - MAGNETIC INDEX(DAILY) OR WHEN SW(9)=-1. :

C - ARRAY CONTAINING:

C (1) DAILY AP

C (2) 3 HR AP INDEX FOR CURRENT TIME

C (3) 3 HR AP INDEX FOR 3 HRS BEFORE CURRENT TIME

C (4) 3 HR AP INDEX FOR 6 HRS BEFORE CURRENT TIME

C (5) 3 HR AP INDEX FOR 9 HRS BEFORE CURRENT TIME

C (6) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 12 TO 33 HRS PRIOR

C TO CURRENT TIME

C (7) AVERAGE OF EIGHT 3 HR AP INDICIES FROM 36 TO 57 HRS PRIOR

C TO CURRENT TIME

C MASS - MASS NUMBER (ONLY DENSITY FOR SELECTED GAS IS

C CALCULATED. MASS 0 IS TEMPERATURE. MASS 48 FOR ALL.

C MASS 17 IS Anomalous O ONLY.)

C

C NOTES ON INPUT VARIABLES:

C UT, Local Time, and Longitude are used independently in the

C model and are not of equal importance for every situation.

C For the most physically realistic calculation these three

C variables should be consistent (STL=SEC/3600+GLONG/15).

C The Equation of Time departures from the above formula

C for apparent local time can be included if available but

C are of minor importance.

c

C F107 and F107A values used to generate the model correspond

C to the 10.7 cm radio flux at the actual distance of the Earth

C from the Sun rather than the radio flux at 1 AU. The following

C site provides both classes of values:

C ftp://ftp.ngdc.noaa.gov/STP/SOLAR_DATA/SOLAR_RADIO/FLUX/

C

C F107, F107A, and AP effects are neither large nor well

C established below 80 km and these parameters should be set to

C 150., 150., and 4. respectively.

C

C OUTPUT VARIABLES:

C D(1) - HE NUMBER DENSITY(CM-3)

C D(2) - O NUMBER DENSITY(CM-3)

C D(3) - N2 NUMBER DENSITY(CM-3)

C D(4) - O2 NUMBER DENSITY(CM-3)

C D(5) - AR NUMBER DENSITY(CM-3)

C D(6) - TOTAL MASS DENSITY(GM/CM3) [Anomalous O NOT included]

C D(7) - H NUMBER DENSITY(CM-3)

C D(8) - N NUMBER DENSITY(CM-3)

C D(9) - Anomalous oxygen NUMBER DENSITY(CM-3)

C T(1) - EXOSPHERIC TEMPERATURE

C T(2) - TEMPERATURE AT ALT

C

DIMENSION ZN1(5),ALPHA(9)

COMMON/GTS3C/TLB,S,DB04,DB16,DB28,DB32,DB40,DB48,DB01,ZA,T0,Z0

& ,G0,RL,DD,DB14,TR12

COMMON/MESO7/TN1(5),TN2(4),TN3(5),TGN1(2),TGN2(2),TGN3(2)

DIMENSION D(9),T(2),MT(11),AP(7),ALTL(8)

COMMON/PARM7/PT(150),PD(150,9),PS(150),PDL(25,2),PTL(100,4),

$ PMA(100,10),SAM(100)

COMMON/CSW/SW(25),ISW,SWC(25)

COMMON/TTEST/TINFG,GB,ROUT,TT(15)

COMMON/DMIX/DM04,DM16,DM28,DM32,DM40,DM01,DM14

COMMON/METSEL/IMR

SAVE

DATA MT/48,0,4,16,28,32,40,1,49,14,17/

DATA ALTL/200.,300.,160.,250.,240.,450.,320.,450./

DATA MN1/5/,ZN1/120.,110.,100.,90.,72.5/

DATA DGTR/1.74533E-2/,DR/1.72142E-2/,ALAST/-999./

DATA ALPHA/-0.38,0.,0.,0.,0.17,0.,-0.38,0.,0./

C Test for changed input

V2=VTST7(IYD,SEC,GLAT,GLONG,STL,F107A,F107,AP,2)

C

YRD=IYD

ZA=PDL(16,2)

ZN1(1)=ZA

DO 2 J=1,9

D(J)=0.

2 CONTINUE

C TINF VARIATIONS NOT IMPORTANT BELOW ZA OR ZN1(1)

IF(ALT.GT.ZN1(1)) THEN

IF(V2.EQ.1..OR.ALAST.LE.ZN1(1)) TINF=PTM(1)*PT(1)

$ *(1.+SW(16)*GLOBE7(YRD,SEC,GLAT,GLONG,STL,F107A,F107,AP,PT))

ELSE

TINF=PTM(1)*PT(1)

ENDIF

T(1)=TINF

C GRADIENT VARIATIONS NOT IMPORTANT BELOW ZN1(5)

IF(ALT.GT.ZN1(5)) THEN

IF(V2.EQ.1.OR.ALAST.LE.ZN1(5)) G0=PTM(4)*PS(1)

$ *(1.+SW(19)*GLOBE7(YRD,SEC,GLAT,GLONG,STL,F107A,F107,AP,PS))

ELSE

G0=PTM(4)*PS(1)

ENDIF

C Calculate these temperatures only if input changed

IF(V2.EQ.1. .OR. ALT.LT.300.)

$ TLB=PTM(2)*(1.+SW(17)*GLOBE7(YRD,SEC,GLAT,GLONG,STL,

$ F107A,F107,AP,PD(1,4)))*PD(1,4)

S=G0/(TINF-TLB)

C Lower thermosphere temp variations not significant for

C density above 300 km

IF(ALT.LT.300.) THEN

IF(V2.EQ.1..OR.ALAST.GE.300.) THEN

TN1(2)=PTM(7)*PTL(1,1)/(1.-SW(18)*GLOB7S(PTL(1,1)))

TN1(3)=PTM(3)*PTL(1,2)/(1.-SW(18)*GLOB7S(PTL(1,2)))

TN1(4)=PTM(8)*PTL(1,3)/(1.-SW(18)*GLOB7S(PTL(1,3)))

TN1(5)=PTM(5)*PTL(1,4)/(1.-SW(18)*SW(20)*GLOB7S(PTL(1,4)))

TGN1(2)=PTM(9)*PMA(1,9)*(1.+SW(18)*SW(20)*GLOB7S(PMA(1,9)))

$ *TN1(5)*TN1(5)/(PTM(5)*PTL(1,4))**2

ENDIF

ELSE

TN1(2)=PTM(7)*PTL(1,1)

TN1(3)=PTM(3)*PTL(1,2)

TN1(4)=PTM(8)*PTL(1,3)

TN1(5)=PTM(5)*PTL(1,4)

TGN1(2)=PTM(9)*PMA(1,9)

$ *TN1(5)*TN1(5)/(PTM(5)*PTL(1,4))**2

ENDIF

C

Z0=ZN1(4)

T0=TN1(4)

TR12=1.

C

IF(MASS.EQ.0) GO TO 50

C N2 variation factor at Zlb

G28=SW(21)*GLOBE7(YRD,SEC,GLAT,GLONG,STL,F107A,F107,

& AP,PD(1,3))

DAY=AMOD(YRD,1000.)

C VARIATION OF TURBOPAUSE HEIGHT

ZHF=PDL(25,2)

$ *(1.+SW(5)*PDL(25,1)*SIN(DGTR*GLAT)*COS(DR*(DAY-PT(14))))

YRD=IYD

T(1)=TINF

XMM=PDM(5,3)

Z=ALT

C

DO 10 J = 1,11

IF(MASS.EQ.MT(J)) GO TO 15

10 CONTINUE

WRITE(6,100) MASS

GO TO 90

15 IF(Z.GT.ALTL(6).AND.MASS.NE.28.AND.MASS.NE.48) GO TO 17

C

C **** N2 DENSITY ****

C

C Diffusive density at Zlb

DB28 = PDM(1,3)*EXP(G28)*PD(1,3)

C Diffusive density at Alt

D(3)=DENSU(Z,DB28,TINF,TLB, 28.,ALPHA(3),T(2),PTM(6),S,MN1,ZN1,

& TN1,TGN1)

DD=D(3)

C Turbopause

ZH28=PDM(3,3)*ZHF

ZHM28=PDM(4,3)*PDL(6,2)

XMD=28.-XMM

C Mixed density at Zlb

B28=DENSU(ZH28,DB28,TINF,TLB,XMD,ALPHA(3)-1.,TZ,PTM(6),S,MN1,

& ZN1,TN1,TGN1)

IF(Z.GT.ALTL(3).OR.SW(15).EQ.0.) GO TO 17

C Mixed density at Alt

DM28=DENSU(Z,B28,TINF,TLB,XMM,ALPHA(3),TZ,PTM(6),S,MN1,

& ZN1,TN1,TGN1)

C Net density at Alt

D(3)=DNET(D(3),DM28,ZHM28,XMM,28.)

17 CONTINUE

GO TO (20,50,20,25,90,35,40,45,25,48,46), J

20 CONTINUE

C

C **** HE DENSITY ****

C

C Density variation factor at Zlb

G4 = SW(21)*GLOBE7(YRD,SEC,GLAT,GLONG,STL,F107A,F107,AP,PD(1,1))

C Diffusive density at Zlb

DB04 = PDM(1,1)*EXP(G4)*PD(1,1)

C Diffusive density at Alt

D(1)=DENSU(Z,DB04,TINF,TLB, 4.,ALPHA(1),T(2),PTM(6),S,MN1,ZN1,

& TN1,TGN1)

DD=D(1)

IF(Z.GT.ALTL(1).OR.SW(15).EQ.0.) GO TO 24

C Turbopause

ZH04=PDM(3,1)

C Mixed density at Zlb

B04=DENSU(ZH04,DB04,TINF,TLB,4.-XMM,ALPHA(1)-1.,

$ T(2),PTM(6),S,MN1,ZN1,TN1,TGN1)

C Mixed density at Alt

DM04=DENSU(Z,B04,TINF,TLB,XMM,0.,T(2),PTM(6),S,MN1,ZN1,TN1,TGN1)

ZHM04=ZHM28

C Net density at Alt

D(1)=DNET(D(1),DM04,ZHM04,XMM,4.)

C Correction to specified mixing ratio at ground

RL=ALOG(B28*PDM(2,1)/B04)

ZC04=PDM(5,1)*PDL(1,2)

HC04=PDM(6,1)*PDL(2,2)

C Net density corrected at Alt

D(1)=D(1)*CCOR(Z,RL,HC04,ZC04)

24 CONTINUE

IF(MASS.NE.48) GO TO 90

25 CONTINUE

C

C **** O DENSITY ****

C

C Density variation factor at Zlb

G16= SW(21)*GLOBE7(YRD,SEC,GLAT,GLONG,STL,F107A,F107,AP,PD(1,2))

C Diffusive density at Zlb

DB16 = PDM(1,2)*EXP(G16)*PD(1,2)

C Diffusive density at Alt

D(2)=DENSU(Z,DB16,TINF,TLB, 16.,ALPHA(2),T(2),PTM(6),S,MN1,

$ ZN1,TN1,TGN1)

DD=D(2)

IF(Z.GT.ALTL(2).OR.SW(15).EQ.0.) GO TO 34

C Corrected from PDM(3,1) to PDM(3,2) 12/2/85

C Turbopause

ZH16=PDM(3,2)

C Mixed density at Zlb

B16=DENSU(ZH16,DB16,TINF,TLB,16-XMM,ALPHA(2)-1.,

$ T(2),PTM(6),S,MN1,ZN1,TN1,TGN1)

C Mixed density at Alt

DM16=DENSU(Z,B16,TINF,TLB,XMM,0.,T(2),PTM(6),S,MN1,ZN1,TN1,TGN1)

ZHM16=ZHM28

C Net density at Alt

D(2)=DNET(D(2),DM16,ZHM16,XMM,16.)

C 3/16/99 Change form to match O2 departure from diff equil near 150

C km and add dependence on F10.7

C RL=ALOG(B28*PDM(2,2)*ABS(PDL(17,2))/B16)

RL=PDM(2,2)*PDL(17,2)*(1.+SW(1)*PDL(24,1)*(F107A-150.))

HC16=PDM(6,2)*PDL(4,2)

ZC16=PDM(5,2)*PDL(3,2)

HC216=PDM(6,2)*PDL(5,2)

D(2)=D(2)*CCOR2(Z,RL,HC16,ZC16,HC216)

C Chemistry correction

HCC16=PDM(8,2)*PDL(14,2)

ZCC16=PDM(7,2)*PDL(13,2)

RC16=PDM(4,2)*PDL(15,2)

C Net density corrected at Alt

D(2)=D(2)*CCOR(Z,RC16,HCC16,ZCC16)

34 CONTINUE

IF(MASS.NE.48.AND.MASS.NE.49) GO TO 90

35 CONTINUE

C

C **** O2 DENSITY ****

C

C Density variation factor at Zlb

G32= SW(21)*GLOBE7(YRD,SEC,GLAT,GLONG,STL,F107A,F107,AP,PD(1,5))

C Diffusive density at Zlb

DB32 = PDM(1,4)*EXP(G32)*PD(1,5)

C Diffusive density at Alt

D(4)=DENSU(Z,DB32,TINF,TLB, 32.,ALPHA(4),T(2),PTM(6),S,MN1,

$ ZN1,TN1,TGN1)

IF(MASS.EQ.49) THEN

DD=DD+2.*D(4)

ELSE

DD=D(4)

ENDIF

IF(SW(15).EQ.0.) GO TO 39

IF(Z.GT.ALTL(4)) GO TO 38

C Turbopause

ZH32=PDM(3,4)

C Mixed density at Zlb

B32=DENSU(ZH32,DB32,TINF,TLB,32.-XMM,ALPHA(4)-1.,

$ T(2),PTM(6),S,MN1,ZN1,TN1,TGN1)

C Mixed density at Alt

DM32=DENSU(Z,B32,TINF,TLB,XMM,0.,T(2),PTM(6),S,MN1,ZN1,TN1,TGN1)

ZHM32=ZHM28

C Net density at Alt

D(4)=DNET(D(4),DM32,ZHM32,XMM,32.)

C Correction to specified mixing ratio at ground

RL=ALOG(B28*PDM(2,4)/B32)

HC32=PDM(6,4)*PDL(8,2)

ZC32=PDM(5,4)*PDL(7,2)

D(4)=D(4)*CCOR(Z,RL,HC32,ZC32)

38 CONTINUE

C Correction for general departure from diffusive equilibrium above Zlb

HCC32=PDM(8,4)*PDL(23,2)

HCC232=PDM(8,4)*PDL(23,1)

ZCC32=PDM(7,4)*PDL(22,2)

RC32=PDM(4,4)*PDL(24,2)*(1.+SW(1)*PDL(24,1)*(F107A-150.))

C Net density corrected at Alt

D(4)=D(4)*CCOR2(Z,RC32,HCC32,ZCC32,HCC232)

39 CONTINUE

IF(MASS.NE.48) GO TO 90

40 CONTINUE

C

C **** AR DENSITY ****

C

C Density variation factor at Zlb

G40= SW(21)*GLOBE7(YRD,SEC,GLAT,GLONG,STL,F107A,F107,AP,PD(1,6))

C Diffusive density at Zlb

DB40 = PDM(1,5)*EXP(G40)*PD(1,6)

C Diffusive density at Alt

D(5)=DENSU(Z,DB40,TINF,TLB, 40.,ALPHA(5),T(2),PTM(6),S,MN1,

$ ZN1,TN1,TGN1)

DD=D(5)

IF(Z.GT.ALTL(5).OR.SW(15).EQ.0.) GO TO 44

C Turbopause

ZH40=PDM(3,5)

C Mixed density at Zlb

B40=DENSU(ZH40,DB40,TINF,TLB,40.-XMM,ALPHA(5)-1.,

$ T(2),PTM(6),S,MN1,ZN1,TN1,TGN1)

C Mixed density at Alt

DM40=DENSU(Z,B40,TINF,TLB,XMM,0.,T(2),PTM(6),S,MN1,ZN1,TN1,TGN1)

ZHM40=ZHM28

C Net density at Alt

D(5)=DNET(D(5),DM40,ZHM40,XMM,40.)

C Correction to specified mixing ratio at ground

RL=ALOG(B28*PDM(2,5)/B40)

HC40=PDM(6,5)*PDL(10,2)

ZC40=PDM(5,5)*PDL(9,2)

C Net density corrected at Alt

D(5)=D(5)*CCOR(Z,RL,HC40,ZC40)

44 CONTINUE

IF(MASS.NE.48) GO TO 90

45 CONTINUE

C

C **** HYDROGEN DENSITY ****

C

C Density variation factor at Zlb

G1 = SW(21)*GLOBE7(YRD,SEC,GLAT,GLONG,STL,F107A,F107,AP,PD(1,7))

C Diffusive density at Zlb

DB01 = PDM(1,6)*EXP(G1)*PD(1,7)

C Diffusive density at Alt

D(7)=DENSU(Z,DB01,TINF,TLB,1.,ALPHA(7),T(2),PTM(6),S,MN1,

$ ZN1,TN1,TGN1)

DD=D(7)

IF(Z.GT.ALTL(7).OR.SW(15).EQ.0.) GO TO 47

C Turbopause

ZH01=PDM(3,6)

C Mixed density at Zlb

B01=DENSU(ZH01,DB01,TINF,TLB,1.-XMM,ALPHA(7)-1.,

$ T(2),PTM(6),S,MN1,ZN1,TN1,TGN1)

C Mixed density at Alt

DM01=DENSU(Z,B01,TINF,TLB,XMM,0.,T(2),PTM(6),S,MN1,ZN1,TN1,TGN1)

ZHM01=ZHM28

C Net density at Alt

D(7)=DNET(D(7),DM01,ZHM01,XMM,1.)

C Correction to specified mixing ratio at ground

RL=ALOG(B28*PDM(2,6)*ABS(PDL(18,2))/B01)

HC01=PDM(6,6)*PDL(12,2)

ZC01=PDM(5,6)*PDL(11,2)

D(7)=D(7)*CCOR(Z,RL,HC01,ZC01)

C Chemistry correction

HCC01=PDM(8,6)*PDL(20,2)

ZCC01=PDM(7,6)*PDL(19,2)

RC01=PDM(4,6)*PDL(21,2)

C Net density corrected at Alt

D(7)=D(7)*CCOR(Z,RC01,HCC01,ZCC01)

47 CONTINUE

IF(MASS.NE.48) GO TO 90

48 CONTINUE

C

C **** ATOMIC NITROGEN DENSITY ****

C

C Density variation factor at Zlb

G14 = SW(21)*GLOBE7(YRD,SEC,GLAT,GLONG,STL,F107A,F107,AP,PD(1,8))

C Diffusive density at Zlb

DB14 = PDM(1,7)*EXP(G14)*PD(1,8)

C Diffusive density at Alt

D(8)=DENSU(Z,DB14,TINF,TLB,14.,ALPHA(8),T(2),PTM(6),S,MN1,

$ ZN1,TN1,TGN1)

DD=D(8)

IF(Z.GT.ALTL(8).OR.SW(15).EQ.0.) GO TO 49

C Turbopause

ZH14=PDM(3,7)

C Mixed density at Zlb

B14=DENSU(ZH14,DB14,TINF,TLB,14.-XMM,ALPHA(8)-1.,

$ T(2),PTM(6),S,MN1,ZN1,TN1,TGN1)

C Mixed density at Alt

DM14=DENSU(Z,B14,TINF,TLB,XMM,0.,T(2),PTM(6),S,MN1,ZN1,TN1,TGN1)

ZHM14=ZHM28

C Net density at Alt

D(8)=DNET(D(8),DM14,ZHM14,XMM,14.)

C Correction to specified mixing ratio at ground

RL=ALOG(B28*PDM(2,7)*ABS(PDL(3,1))/B14)

HC14=PDM(6,7)*PDL(2,1)

ZC14=PDM(5,7)*PDL(1,1)

D(8)=D(8)*CCOR(Z,RL,HC14,ZC14)

C Chemistry correction

HCC14=PDM(8,7)*PDL(5,1)

ZCC14=PDM(7,7)*PDL(4,1)

RC14=PDM(4,7)*PDL(6,1)

C Net density corrected at Alt

D(8)=D(8)*CCOR(Z,RC14,HCC14,ZCC14)

49 CONTINUE

IF(MASS.NE.48) GO TO 90

46 CONTINUE

C

C **** Anomalous OXYGEN DENSITY ****

C

G16H = SW(21)*GLOBE7(YRD,SEC,GLAT,GLONG,STL,F107A,F107,AP,PD(1,9))

DB16H = PDM(1,8)*EXP(G16H)*PD(1,9)

THO=PDM(10,8)*PDL(7,1)

DD=DENSU(Z,DB16H,THO,THO,16.,ALPHA(9),T2,PTM(6),S,MN1,

$ ZN1,TN1,TGN1)

ZSHT=PDM(6,8)

ZMHO=PDM(5,8)

ZSHO=SCALH(ZMHO,16.,THO)

D(9)=DD*EXP(-ZSHT/ZSHO*(EXP(-(Z-ZMHO)/ZSHT)-1.))

IF(MASS.NE.48) GO TO 90

C

C TOTAL MASS DENSITY

C

D(6) = 1.66E-24*(4.*D(1)+16.*D(2)+28.*D(3)+32.*D(4)+40.*D(5)+

& D(7)+14.*D(8))

DB48=1.66E-24*(4.*DB04+16.*DB16+28.*DB28+32.*DB32+40.*DB40+DB01+

& 14.*DB14)

GO TO 90

C TEMPERATURE AT ALTITUDE

50 CONTINUE

Z=ABS(ALT)

DDUM = DENSU(Z,1., TINF,TLB,0.,0.,T(2),PTM(6),S,MN1,ZN1,TN1,TGN1)

90 CONTINUE

C ADJUST DENSITIES FROM CGS TO KGM

IF(IMR.EQ.1) THEN

DO 95 I=1,9

D(I)=D(I)*1.E6

95 CONTINUE

D(6)=D(6)/1000.

ENDIF

ALAST=ALT

RETURN

100 FORMAT(1X,'MASS', I5, ' NOT VALID')

END Subroutine GTS7

C-----------------------------------------------------------------------

SUBROUTINE METERS(METER)

implicit none (external)

C Convert outputs to Kg & Meters if METER true

logical,Intent(In) :: METER

Integer IMR

COMMON/METSEL/IMR

SAVE

IMR=0

IF(METER) IMR=1

END SUBROUTINE METERS

C-----------------------------------------------------------------------

FUNCTION SCALH(ALT,XM,TEMP)

C Calculate scale height (km)

COMMON/PARMB/GSURF,RE

SAVE

DATA RGAS/831.4/

G=GSURF/(1.+ALT/RE)**2

SCALH=RGAS*TEMP/(G*XM)

END FUNCTION SCALH

C-----------------------------------------------------------------------

FUNCTION GLOBE7(YRD,SEC,LAT,LONG,TLOC,F107A,F107,AP,P)

C CALCULATE G(L) FUNCTION

C Upper Thermosphere Parameters

REAL LAT, LONG

DIMENSION P(150),SV(25),AP(7)

COMMON/TTEST/TINF,GB,ROUT,T(15)

COMMON/CSW/SW(25),ISW,SWC(25)

COMMON/LPOLY/PLG(9,4),CTLOC,STLOC,C2TLOC,S2TLOC,C3TLOC,S3TLOC,

$ IYR,DAY,DF,DFA,APD,APDF,APT(4),XLONG

SAVE

DATA DGTR/1.74533E-2/,DR/1.72142E-2/, XL/1000./,TLL/1000./

DATA SW9/1./,DAYL/-1./,P14/-1000./,P18/-1000./,P32/-1000./

DATA HR/.2618/,SR/7.2722E-5/,SV/25*1./,NSW/14/,P39/-1000./

C 3hr Magnetic activity functions

C Eq. A24d

G0(A)=(A-4.+(P(26)-1.)*(A-4.+(EXP(-ABS(P(25))*(A-4.))-1.)/ABS(P(25

*))))

C Eq. A24c

SUMEX(EX)=1.+(1.-EX**19)/(1.-EX)*EX**(.5)

C Eq. A24a

SG0(EX)=(G0(AP(2))+(G0(AP(3))*EX+G0(AP(4))*EX*EX+G0(AP(5))*EX**3

$ +(G0(AP(6))*EX**4+G0(AP(7))*EX**12)*(1.-EX**8)/(1.-EX))

$ )/SUMEX(EX)

IF(ISW.NE.64999) CALL TSELEC(SV)

DO 10 J=1,14

T(J)=0

10 CONTINUE

IF(SW(9).GT.0) SW9=1.

IF(SW(9).LT.0) SW9=-1.

IYR = YRD/1000.

DAY = YRD - IYR*1000.

XLONG=LONG

C Eq. A22 (remainder of code)

IF(XL.EQ.LAT) GO TO 15

C CALCULATE LEGENDRE POLYNOMIALS