-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathmsis_constants.F90
196 lines (184 loc) · 15.4 KB
/
msis_constants.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
!#######################################################################
! MSIS® (NRL-SOF-014-1) SOFTWARE
! NRLMSIS® empirical atmospheric model software. Use is governed by the
! Open Source Academic Research License Agreement contained in the file
! nrlmsis2.1_license.txt, which is part of this software package. BY
! USING OR MODIFYING THIS SOFTWARE, YOU ARE AGREEING TO THE TERMS AND
! CONDITIONS OF THE LICENSE.
!#######################################################################
!!! ===========================================================================
!!! NRLMSIS 2.1:
!!! Neutral atmosphere empirical model from the surface to lower exosphere
!!! ===========================================================================
!**************************************************************************************************
! MSIS_CONSTANTS Module: Contains constants and hardwired parameters
!**************************************************************************************************
module msis_constants
implicit none
! Floating Point Precision
#ifdef DBLE
integer, parameter :: rp = 8
#else
integer, parameter :: rp = 4
#endif
! Missing density value
real(kind=rp),parameter :: dmissing = 9.999e-38_rp
! Trigonometric constants
real(kind=rp), parameter :: pi = 3.1415926535897932384626433832795_rp
real(kind=rp), parameter :: deg2rad = pi / 180.0_rp
real(kind=rp), parameter :: doy2rad = 2.0_rp*pi / 365.0_rp
real(kind=rp), parameter :: lst2rad = pi / 12.0_rp
!real(kind=rp), parameter :: tanh1 = 0.761594155955765485_rp ! tanh(1.0)
real(kind=rp), parameter :: tanh1 = tanh(1.0_rp)
! Thermodynamic constants
! Boltzmann constant (CODATA 2018) (J/kg)
real(kind=rp), parameter :: kB = 1.380649e-23_rp
! Avogadro constant (CODATA 2018)
real(kind=rp), parameter :: NA = 6.02214076e23_rp
! Reference gravity (CIMO Guide 2014) (m/s^2) (specified separately in alt2gph, in msis_utils.F90)
real(kind=rp), parameter :: g0 = 9.80665_rp
! Species molecular masses (kg/molecule) (CIPM 2007)
real(kind=rp), parameter :: specmass(1:10) = (/ 0.0_rp, & ! Mass density (dummy value)
28.0134_rp, & ! N2
31.9988_rp, & ! O2
31.9988_rp/2.0_rp, & ! O
4.0_rp, & ! He
1.0_rp, & ! H
39.948_rp, & ! Ar
28.0134_rp/2.0_rp, & ! N
31.9988_rp/2.0_rp, & ! Anomalous O
(28.0134_rp+31.9988_rp)/2.0_rp /) & ! NO
/ (1.0e3_rp * NA) ! Convert from g/mol to kg/molecule
! Dry air mean mass in fully mixed atmosphere (CIPM 2007) (includes CO2 and other trace species that are not yet in MSIS)
real(kind=rp), parameter :: Mbar = 28.96546_rp / (1.0e3_rp * NA) ! kg/molecule
! Dry air log volume mixing ratios (CIPM 2007)
real(kind=rp), parameter :: lnvmr(1:10) = log( (/ 1.0_rp, & ! Mass density (dummy value)
0.780848_rp, & ! N2
0.209390_rp, & ! O2
1.0_rp, & ! O (dummy value)
0.0000052_rp, & ! He
1.0_rp, & ! H (dummy value)
0.009332_rp, & ! Ar
1.0_rp, & ! N (dummy value)
1.0_rp, & ! Anomalous O (dummy value)
1.0_rp /) ) ! NO (dummy value)
! Natural log of global average surface pressure (Pa)
!real(kind=rp), parameter :: lnP0 = 11.5080482 !+ 0.00759597 After calibration with MERRA2
real(kind=rp), parameter :: lnP0 = 11.515614
! Derived constants
real(kind=rp), parameter :: g0divkB = g0/kB * 1.0e3_rp ! K/(kg km)
real(kind=rp), parameter :: Mbarg0divkB = Mbar*g0/kB * 1.0e3_rp ! K/km
! References:
! CODATA Internationally recommended 2018 values of the fundamental physical constants.
! https://pml.nist.gov/cuu/Constants/; https://pml.nist.gov/cuu/pdf/wallet_2018.pdf
! Picard, A., Davis, R. S., Glaeser, M., and Fujii, K. (2007). Revised formula for the density of
! air (CIPM 2007). Metrologia 45, 149–155. doi:10.1088/0026-1394/45/2/004
! World Meteorological Organization (2014). WMO guide to meteorological instruments and methods of observation
! (the CIMO Guide). Part I, Chapter 12. https://www.wmo.int/pages/prog/www/IMOP/CIMO-Guide.html
! Vertical profile parameters
integer, parameter :: nspec = 11 !Number of species including temperature
integer, parameter :: nd = 27 !Number of temperature profile nodes
integer, parameter :: p = 4 !Spline order
integer, parameter :: nl = nd - p !Last temperature profile level index
integer, parameter :: nls = 9 !Last parameter index for each species (excluding O, NO splines)
real(kind=rp), parameter :: bwalt = 122.5_rp ! Reference geopotential height for Bates Profile
real(kind=rp), parameter :: zetaF = 70.0_rp ! Fully mixed below this, uses constant mixing ratios
real(kind=rp), parameter :: zetaB = bwalt ! Bates Profile above this altitude
real(kind=rp), parameter :: zetaA = 85.0_rp ! Default reference height for active minor species
real(kind=rp), parameter :: zetagamma = 100.0_rp ! Reference height of tanh taper of chemical/dynamical correction scale height
real(kind=rp), parameter :: Hgamma = 1.0_rp/30.0_rp ! Inverse scale height of tanh taper of chemical/dynamical correction scale height
real(kind=rp), parameter :: nodesTN(0:nd+2) = & !Nodes for temperature profile splines
(/ -15., -10., -5., 0., 5., 10., 15., 20., 25., 30., 35., 40., 45., 50., &
55., 60., 65., 70., 75., 80., 85., 92.5, 102.5, 112.5, 122.5, 132.5, 142.5, &
152.5, 162.5, 172.5/)
integer, parameter :: izfmx = 13 ! fully mixed below this spline index
integer, parameter :: izfx = 14 ! Spline index at zetaF
integer, parameter :: izax = 17 ! Spline index at zetaA
integer, parameter :: itex = nl ! Index of Bates exospheric temperature
integer, parameter :: itgb0 = nl - 1 ! Index of Bates temperature gradient at lower boundary
integer, parameter :: itb0 = nl - 2 ! Index of Bates temperature at lower boundary
! O1 Spline parameters
integer, parameter :: ndO1 = 13
integer, parameter :: nsplO1 = ndO1-5 !Number of unconstrained spline parameters for O1 (there are 2 additional C1-constrained splines)
real(kind=rp), parameter :: nodesO1(0:ndO1) = & !Nodes for O1 splines (Domain 50-85 km)
(/ 35., 40., 45., 50., 55., 60., 65., 70., 75., 80., 85., 92.5, 102.5, 112.5/)
real(kind=rp), parameter :: zetarefO1 = zetaA !Joining height for O1 splines, and reference height for O1 density
! NO Spline parameters
integer, parameter :: ndNO = 13
integer, parameter :: nsplNO = ndNO-5 !Number of unconstrained spline parameters for NO (there are 2 additional C1-constrained splines)
real(kind=rp), parameter :: nodesNO(0:ndNO) = & !Nodes for NO splines (Domain 70-122.5 km)
(/ 47.5, 55., 62.5, 70., 77.5, 85., 92.5, 100., 107.5, 115., 122.5, 130., 137.5, 145./)
real(kind=rp), parameter :: zetarefNO = zetaB !Joining height for NO splines, and reference height for NO density
!C2 Continuity matrix for temperature; Last 3 splines are constrained (must be recomputed if nodes change)
real(kind=rp), parameter :: c2tn(3,3) = reshape((/1.0_rp, -10.0_rp, 33.333333333333336_rp, &
1.0_rp, 0.0_rp, -16.666666666666668_rp, &
1.0_rp, 10.0_rp, 33.333333333333336_rp/), &
(/3,3/))
!C1 Continuity for O1; Last 2 splines are constrained (must be recomputed if nodes change)
real(kind=rp), parameter :: c1o1(2,2) = reshape((/ 1.75_rp, -2.916666573405061_rp, &
-1.624999900076852_rp, 21.458332647194382_rp /), &
(/2,2/))
real(kind=rp), parameter :: c1o1adj(2) = (/0.257142857142857_rp, -0.102857142686844_rp/) !Weights for coefficents on 3rd to last spline; product to be subtracted from RHS of continuity equation
!C1 Continuity for NO; Last 2 splines are constrained (must be recomputed if nodes change)
real(kind=rp), parameter :: c1NO(2,2) = reshape((/ 1.5_rp, -3.75_rp, &
0.0_rp, 15.0_rp /), &
(/2,2/))
real(kind=rp), parameter :: c1NOadj(2) = (/0.166666666666667_rp, -0.066666666666667_rp/) !Weights for coefficents on 3rd to last spline; product to be subtracted from RHS of continuity equation
! Anomalous Oxygen parameters (legacy profile from NRLMSISE-00)
real(kind=rp),parameter :: zetarefOA = zetaB !Reference height for anomalous oxygen density
real(kind=rp),parameter :: TOA = 4000. !Temperature of anomalous oxygen density (K)
real(kind=rp),parameter :: HOA = (kB * TOA) / ( (16.0_rp/(1.0e3_rp*NA)) * g0 ) * 1.0e-3_rp !Hydrostatic scale height of anomalous oxygen density (km)
! Horizontal and time-dependent basis function (gfn) parameters
integer, parameter :: maxnbf = 512 ! Number of basis functions to be allocated
integer, parameter :: maxn = 6 ! Maximum latitude (Legendre) spectral degree
integer, parameter :: maxl = 3 ! Maximum local time (tidal) spectral order
integer, parameter :: maxm = 2 ! Maximum longitude (stationary planetary wave) order
integer, parameter :: maxs = 2 ! Maximimum day of year (intra-annual) Fourier order
integer, parameter :: amaxn = 6 ! Maximum Legendre degree used in time independent and intra-annual zonal mean terms
integer, parameter :: amaxs = 2 ! Maximum intra-annual order used in zonal mean terms
integer, parameter :: tmaxl = 3 ! Maximum tidal order used
integer, parameter :: tmaxn = 6 ! Maximum Legendre degree coupled with tides
integer, parameter :: tmaxs = 2 ! Maximum intra-annual order coupled with tides
integer, parameter :: pmaxm = 2 ! Maximum stationary planetary wave order used
integer, parameter :: pmaxn = 6 ! Maximum Legendre degree coupled with SPW
integer, parameter :: pmaxs = 2 ! Maximum intra-annual order coupled with SPW
integer, parameter :: nsfx = 5 ! Number of linear solar flux terms
integer, parameter :: nsfxmod = 5 ! Number of nonlinear modulating solar flux terms (legacy NRLMSISE-00 terms)
integer, parameter :: nmag = 54 ! Number of terms in NRLMSISE-00 legacy geomagnetic parameterization
integer, parameter :: nut = 12 ! Number of terms in NRLMSISE-00 legacy UT parameterization
integer, parameter :: ctimeind = 0 ! Starting index of time-independent terms
integer, parameter :: cintann = ctimeind + (amaxn+1) ! Starting index of zonal mean intra-annual terms
integer, parameter :: ctide = cintann + ((amaxn+1)*2*amaxs) ! Starting index of zonal mean intra-annual terms
integer, parameter :: cspw = ctide + (4*tmaxs+2)*(tmaxl*(tmaxn+1)-(tmaxl*(tmaxl+1))/2) ! Starting index of SPW terms
integer, parameter :: csfx = cspw + (4*pmaxs+2)*(pmaxm*(pmaxn+1)-(pmaxm*(pmaxm+1))/2) ! Starting index of linear solar flux terms
integer, parameter :: cextra = csfx + nsfx ! Starting index of time-independent terms
integer, parameter :: mbf = 383 ! Last index of linear terms
integer, parameter :: cnonlin = mbf + 1 ! Starting index of nonlinear terms
integer, parameter :: csfxmod = cnonlin ! Starting index of modulating solar flux terms
integer, parameter :: cmag = csfxmod + nsfxmod ! Starting index of daily geomagnetic terms
integer, parameter :: cut = cmag + nmag ! Starting index of UT terms
! Weights for calculation log pressure spline coefficients from temperature coefficients (must be recalcuated if nodes change)
real(kind=rp), parameter :: gwht(0:3) = (/ 5.0_rp/24.0_rp, 55.0_rp/24.0_rp, 55.0_rp/24.0_rp, 5.0_rp/24.0_rp /)
! Constants needed for analytical integration by parts of hydrostatic piecewise effective mass profile
real(kind=rp), parameter :: wbeta(0:nl) = (nodesTN(4:nd) - nodesTN(0:nl)) / 4.0_rp !Weights for 1st spline integration
real(kind=rp), parameter :: wgamma(0:nl) = (nodesTN(5:nd+1)- nodesTN(0:nl)) / 5.0_rp !Weights for 2nd spline integration
! Non-zero bspline values at zetaB (5th and 6th order) (must be recalcuated if nodes change)
real(kind=rp), parameter :: S5zetaB(0:3) = (/0.041666666666667_rp, 0.458333333333333_rp, 0.458333333333333_rp, &
0.041666666666667_rp/)
real(kind=rp), parameter :: S6zetaB(0:4) = (/0.008771929824561_rp, 0.216228070175439_rp, 0.550000000000000_rp, &
0.216666666666667_rp, 0.008333333333333_rp/)
!Weights for calculating temperature gradient at zetaA (must be recalcuated if nodes change)
real(kind=rp), parameter :: wghtAxdz(0:2) = (/-0.102857142857_rp, 0.0495238095238_rp, 0.053333333333_rp/)
!Non-zero bspline values at zetaA (4th, 5th and 6th order) (must be recalcuated if nodes change)
real(kind=rp), parameter :: S4zetaA(0:2) = (/0.257142857142857_rp, 0.653968253968254_rp, 0.088888888888889_rp/)
real(kind=rp), parameter :: S5zetaA(0:3) = (/0.085714285714286_rp, 0.587590187590188_rp, 0.313020313020313_rp, &
0.013675213675214_rp/)
real(kind=rp), parameter :: S6zetaA(0:4) = (/0.023376623376623_rp, 0.378732378732379_rp, 0.500743700743701_rp, &
0.095538448479625_rp, 0.001608848667672_rp/)
!Non-zero bspline values at zetaF (4th and 5th order) (must be recalcuated if nodes change)
real(kind=rp), parameter :: S4zetaF(0:2) = (/0.166666666666667_rp, 0.666666666666667_rp, 0.166666666666667_rp/)
real(kind=rp), parameter :: S5zetaF(0:3) = (/0.041666666666667_rp, 0.458333333333333_rp, 0.458333333333333_rp, &
0.041666666666667_rp/)
!Non-zero bspline values at zeta=0 (5th order) (must be recalcuated if nodes change)
real(kind=rp), parameter :: S5zeta0(0:2) = (/0.458333333333333_rp, 0.458333333333333_rp, 0.041666666666667_rp/)
end module msis_constants