|
| 1 | +!####################################################################### |
| 2 | +! MSIS® (NRL-SOF-014-1) SOFTWARE |
| 3 | +! NRLMSIS® empirical atmospheric model software. Use is governed by the |
| 4 | +! Open Source Academic Research License Agreement contained in the file |
| 5 | +! nrlmsis2.1_license.txt, which is part of this software package. BY |
| 6 | +! USING OR MODIFYING THIS SOFTWARE, YOU ARE AGREEING TO THE TERMS AND |
| 7 | +! CONDITIONS OF THE LICENSE. |
| 8 | +!####################################################################### |
| 9 | + |
| 10 | +!!! =========================================================================== |
| 11 | +!!! NRLMSIS 2.1: |
| 12 | +!!! Neutral atmosphere empirical model from the surface to lower exosphere |
| 13 | +!!! =========================================================================== |
| 14 | +!!! |
| 15 | +!!! MSISCALC: Interface with re-ordered input arguments and output arrays. |
| 16 | +! |
| 17 | +! PREREQUISITES: |
| 18 | +! Must first run MSISINIT to load parameters and set switches. The |
| 19 | +! MSISCALC subroutine checks for initialization and does a default |
| 20 | +! initialization if necessary. This self-initialization will be removed |
| 21 | +! in future versions. |
| 22 | +! |
| 23 | +! CALLING SEQUENCE: |
| 24 | +! CALL MSISCALC(DAY, UTSEC, Z, LAT, LON, SFLUXAVG, SFLUX, AP, TN, DN, [TEX]) |
| 25 | +! |
| 26 | +! INPUT VARIABLES: |
| 27 | +! DAY Day of year (1.0 to 365.0 or 366.0) |
| 28 | +! UTSEC Universal time (seconds) |
| 29 | +! Z Geodetic altitude (km) (default) or Geopotential height (km) |
| 30 | +! LAT Geodetic latitude (deg) |
| 31 | +! LON Geodetic longitude (deg) |
| 32 | +! SFLUXAVG 81 day average, centered on input time, of F10.7 solar |
| 33 | +! activity index |
| 34 | +! SFLUX Daily F10.7 for previous day |
| 35 | +! AP Geomagnetic activity index array: |
| 36 | +! (1) Daily Ap |
| 37 | +! (2) 3 hr ap index for current time |
| 38 | +! (3) 3 hr ap index for 3 hrs before current time |
| 39 | +! (4) 3 hr ap index for 6 hrs before current time |
| 40 | +! (5) 3 hr ap index for 9 hrs before current time |
| 41 | +! (6) Average of eight 3 hr ap indices from 12 to 33 hrs |
| 42 | +! prior to current time |
| 43 | +! (7) Average of eight 3 hr ap indices from 36 to 57 hrs |
| 44 | +! prior to current time |
| 45 | +! AP(2:7) are only used when switch_legacy(9) = -1.0 in MSISINIT |
| 46 | +! |
| 47 | +! NOTES ON INPUT VARIABLES: |
| 48 | +! - The day-of-year dependence of the model only uses the DAY argument. If |
| 49 | +! a continuous day-of-year dependence is desired, this argument should |
| 50 | +! include the fractional day (e.g., DAY = <day of year> + UTSEC/86400.0 |
| 51 | +! - If lzalt_type = .true. (default) in the MSISINIT call, then Z is |
| 52 | +! treated as geodetic altitude. |
| 53 | +! If lzalt_type = .false., then Z is treated as geopotential height. |
| 54 | +! - F107 and F107A values are the 10.7 cm radio flux at the Sun-Earth |
| 55 | +! distance, not the radio flux at 1 AU. |
| 56 | +! |
| 57 | +! OUTPUT VARIABLES: |
| 58 | +! TN Temperature at altitude (K) |
| 59 | +! DN(1) Total mass density (kg/m3) |
| 60 | +! DN(2) N2 number density (m-3) |
| 61 | +! DN(3) O2 number density (m-3) |
| 62 | +! DN(4) O number density (m-3) |
| 63 | +! DN(5) He number density (m-3) |
| 64 | +! DN(6) H number density (m-3) |
| 65 | +! DN(7) Ar number density (m-3) |
| 66 | +! DN(8) N number density (m-3) |
| 67 | +! DN(9) Anomalous oxygen number density (m-3) |
| 68 | +! DN(10) NO number density (m-3) |
| 69 | +! TEX Exospheric temperature (K) (optional argument) |
| 70 | +! |
| 71 | +! NOTES ON OUTPUT VARIABLES: |
| 72 | +! - Missing density values are returned as 9.999e-38 |
| 73 | +! - Species included in mass density calculation are set in MSISINIT |
| 74 | +! |
| 75 | +!!! ========================================================================= |
| 76 | + |
| 77 | +!************************************************************************************************** |
| 78 | +! MSIS_CALC Module: Contains main MSIS entry point |
| 79 | +!************************************************************************************************** |
| 80 | +module msis_calc |
| 81 | + |
| 82 | +contains |
| 83 | + |
| 84 | + !================================================================================================== |
| 85 | + ! MSISCALC: The main MSIS subroutine entry point |
| 86 | + !================================================================================================== |
| 87 | + subroutine msiscalc(day,utsec,z,lat,lon,sfluxavg,sflux,ap,tn,dn,tex) |
| 88 | + |
| 89 | + use msis_constants, only : rp, dmissing, lnp0, Mbarg0divkB, kB, nspec, nodesTN, nd, zetaF, zetaB, & |
| 90 | + Hgamma, zetagamma, maxnbf |
| 91 | + use msis_init, only : msisinit, initflag, zaltflag, specflag, massflag, masswgt, etaTN |
| 92 | + use msis_gfn, only : globe |
| 93 | + use msis_tfn, only : tnparm, tfnparm, tfnx |
| 94 | + use msis_dfn, only : dnparm, dfnparm, dfnx |
| 95 | + use msis_utils, only : alt2gph, bspline, dilog |
| 96 | + |
| 97 | + implicit none |
| 98 | + |
| 99 | + real(kind=rp), intent(in) :: day |
| 100 | + real(kind=rp), intent(in) :: utsec |
| 101 | + real(kind=rp), intent(in) :: z |
| 102 | + real(kind=rp), intent(in) :: lat |
| 103 | + real(kind=rp), intent(in) :: lon |
| 104 | + real(kind=rp), intent(in) :: sfluxavg,sflux,ap(1:7) |
| 105 | + real(kind=rp), intent(out) :: tn, dn(1:10) |
| 106 | + real(kind=rp), intent(out), optional :: tex |
| 107 | + |
| 108 | + real(kind=rp), save :: lastday = -9999.0 |
| 109 | + real(kind=rp), save :: lastutsec = -9999.0 |
| 110 | + real(kind=rp), save :: lastlat = -9999.0 |
| 111 | + real(kind=rp), save :: lastlon = -9999.0 |
| 112 | + real(kind=rp), save :: lastz = -9999.0 |
| 113 | + real(kind=rp), save :: lastsflux = -9999.0 |
| 114 | + real(kind=rp), save :: lastsfluxavg = -9999.0 |
| 115 | + real(kind=rp), save :: lastap(1:7) = -9999.0 |
| 116 | + real(kind=rp), save :: gf(0:maxnbf-1) |
| 117 | + real(kind=rp), save :: Sz(-5:0,2:6) |
| 118 | + integer, save :: iz |
| 119 | + type(tnparm), save :: tpro |
| 120 | + type(dnparm), save :: dpro(1:nspec-1) |
| 121 | + |
| 122 | + real(8) :: zaltd, latd |
| 123 | + real(kind=rp) :: zeta, lndtotz, Vz, Wz, HRfact, lnPz, delz |
| 124 | + integer :: i, j, kmax, ispec |
| 125 | + |
| 126 | + ! Check if model has been initialized; if not, perform default initialization |
| 127 | + if (.not. initflag) call msisinit() |
| 128 | + |
| 129 | + ! Calculate geopotential height, if necessary |
| 130 | + if(zaltflag) then |
| 131 | + zaltd = dble(z) |
| 132 | + latd = dble(lat) |
| 133 | + zeta = alt2gph(latd,zaltd) |
| 134 | + else |
| 135 | + zeta = z |
| 136 | + endif |
| 137 | + |
| 138 | + ! If only altitude changes then update the local spline weights |
| 139 | + if (zeta .lt. zetaB) then |
| 140 | + if (zeta .ne. lastz) then |
| 141 | + if (zeta .lt. zetaF) then |
| 142 | + kmax = 5 |
| 143 | + else |
| 144 | + kmax = 6 |
| 145 | + endif |
| 146 | + call bspline(zeta,nodesTN,nd+2,kmax,etaTN,Sz,iz) |
| 147 | + lastz = zeta |
| 148 | + endif |
| 149 | + endif |
| 150 | + |
| 151 | + ! If location, time, or solar/geomagnetic conditions change then recompute the profile parameters |
| 152 | + if ((day .ne. lastday) .or. (utsec .ne. lastutsec) .or. & |
| 153 | + (lat .ne. lastlat) .or. (lon .ne. lastlon) .or. & |
| 154 | + (sflux .ne. lastsflux) .or. (sfluxavg .ne. lastsfluxavg) .or. & |
| 155 | + any(ap .ne. lastap)) then |
| 156 | + call globe(day,utsec,lat,lon,sfluxavg,sflux,ap,gf) |
| 157 | + call tfnparm(gf,tpro) |
| 158 | + do ispec = 2, nspec-1 |
| 159 | + if (specflag(ispec)) call dfnparm(ispec,gf,tpro,dpro(ispec)) |
| 160 | + enddo |
| 161 | + lastday = day |
| 162 | + lastutsec = utsec |
| 163 | + lastlat = lat |
| 164 | + lastlon = lon |
| 165 | + lastsflux = sflux |
| 166 | + lastsfluxavg = sfluxavg |
| 167 | + lastap = ap |
| 168 | + endif |
| 169 | + |
| 170 | + ! Exospheric temperature |
| 171 | + if (present(tex)) then |
| 172 | + tex = tpro%tex |
| 173 | + endif |
| 174 | + |
| 175 | + ! Temperature at altitude |
| 176 | + tn = tfnx(zeta,iz,Sz(-3:0,4),tpro) |
| 177 | + |
| 178 | + ! Temperature integration terms at altitude, total number density |
| 179 | + delz = zeta - zetaB |
| 180 | + if (zeta .lt. zetaF) then |
| 181 | + i = max(iz-4,0) |
| 182 | + if (iz .lt. 4) then |
| 183 | + j = -iz |
| 184 | + else |
| 185 | + j = -4 |
| 186 | + endif |
| 187 | + Vz = dot_product(tpro%beta(i:iz),Sz(j:0,5)) + tpro%cVS |
| 188 | + Wz = 0.0_rp |
| 189 | + lnPz = lnP0 - Mbarg0divkB*(Vz - tpro%Vzeta0) |
| 190 | + lndtotz = lnPz - log(kB*tn) |
| 191 | + else |
| 192 | + if (zeta .lt. zetaB) then |
| 193 | + Vz = dot_product(tpro%beta(iz-4:iz),Sz(-4:0,5)) + tpro%cVS |
| 194 | + Wz = dot_product(tpro%gamma(iz-5:iz),Sz(-5:0,6)) + tpro%cVS*delz + tpro%cWS |
| 195 | + else |
| 196 | + Vz = (delz + log(tn/tpro%tex)/tpro%sigma)/tpro%tex + tpro%cVB |
| 197 | + Wz = (0.5_rp*delz*delz + dilog(tpro%b*exp(-tpro%sigma*delz))/tpro%sigmasq)/tpro%tex & |
| 198 | + + tpro%cVB*delz + tpro%cWB |
| 199 | + endif |
| 200 | + endif |
| 201 | + |
| 202 | + ! Species number densities at altitude |
| 203 | + HRfact = 0.5_rp * (1.0_rp + tanh(Hgamma*(zeta - zetagamma))) !Reduction factor for chemical/dynamical correction scale height below zetagamma |
| 204 | + do ispec = 2, nspec-1 |
| 205 | + if (specflag(ispec)) then |
| 206 | + dn(ispec) = dfnx(zeta,tn,lndtotz,Vz,Wz,HRfact,tpro,dpro(ispec)) |
| 207 | + else |
| 208 | + dn(ispec) = dmissing |
| 209 | + endif |
| 210 | + enddo |
| 211 | + |
| 212 | + ! Mass density |
| 213 | + if (specflag(1)) then |
| 214 | + dn(1) = dot_product(dn,masswgt) |
| 215 | + else |
| 216 | + dn(1) = dmissing |
| 217 | + endif |
| 218 | + |
| 219 | + return |
| 220 | + |
| 221 | + end subroutine msiscalc |
| 222 | + |
| 223 | +end module msis_calc |
0 commit comments