-
Notifications
You must be signed in to change notification settings - Fork 110
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Remove goto statements from chgres_cube (#775)
The WAM option was using the NRLMSISE-00 library, which contains numerous 'goto' statements. That library was replaced with the latest NRLMSIS 2.1 library, which does not use 'goto' statements and is also a welcome improvement in initializing the upper atmosphere. The new library included a unit test, which was added to run under Github actions. Fixes #759.
- Loading branch information
Showing
22 changed files
with
3,466 additions
and
2,808 deletions.
There are no files selected for viewing
Binary file not shown.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,26 @@ | ||
| include(FetchContent) | ||
| include(CheckFortranCompilerFlag) | ||
|
|
||
| add_library(msis2 | ||
| msis_utils.F90 | ||
| msis_constants.F90 | ||
| msis_init.F90 | ||
| msis_gfn.F90 | ||
| msis_tfn.F90 | ||
| msis_dfn.F90 | ||
| msis_calc.F90 | ||
| msis_gtd8d.F90) | ||
|
|
||
| set_target_properties(msis2 PROPERTIES Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/include) | ||
|
|
||
| target_include_directories(msis2 INTERFACE ${CMAKE_CURRENT_BINARY_DIR}/include) | ||
|
|
||
| if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU) | ||
| # msis_calc:bspline has argument mismatch on nodes variable | ||
| target_compile_options(msis2 PRIVATE -std=legacy) | ||
| endif() | ||
|
|
||
| target_link_libraries( | ||
| msis2 | ||
| PUBLIC | ||
| MPI::MPI_Fortran) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,223 @@ | ||
| !####################################################################### | ||
| ! 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 | ||
| !!! =========================================================================== | ||
| !!! | ||
| !!! MSISCALC: Interface with re-ordered input arguments and output arrays. | ||
| ! | ||
| ! PREREQUISITES: | ||
| ! Must first run MSISINIT to load parameters and set switches. The | ||
| ! MSISCALC subroutine checks for initialization and does a default | ||
| ! initialization if necessary. This self-initialization will be removed | ||
| ! in future versions. | ||
| ! | ||
| ! CALLING SEQUENCE: | ||
| ! CALL MSISCALC(DAY, UTSEC, Z, LAT, LON, SFLUXAVG, SFLUX, AP, TN, DN, [TEX]) | ||
| ! | ||
| ! INPUT VARIABLES: | ||
| ! DAY Day of year (1.0 to 365.0 or 366.0) | ||
| ! UTSEC Universal time (seconds) | ||
| ! Z Geodetic altitude (km) (default) or Geopotential height (km) | ||
| ! LAT Geodetic latitude (deg) | ||
| ! LON Geodetic longitude (deg) | ||
| ! SFLUXAVG 81 day average, centered on input time, of F10.7 solar | ||
| ! activity index | ||
| ! SFLUX Daily F10.7 for previous day | ||
| ! AP Geomagnetic activity index array: | ||
| ! (1) Daily Ap | ||
| ! (2) 3 hr ap index for current time | ||
| ! (3) 3 hr ap index for 3 hrs before current time | ||
| ! (4) 3 hr ap index for 6 hrs before current time | ||
| ! (5) 3 hr ap index for 9 hrs before current time | ||
| ! (6) Average of eight 3 hr ap indices from 12 to 33 hrs | ||
| ! prior to current time | ||
| ! (7) Average of eight 3 hr ap indices from 36 to 57 hrs | ||
| ! prior to current time | ||
| ! AP(2:7) are only used when switch_legacy(9) = -1.0 in MSISINIT | ||
| ! | ||
| ! NOTES ON INPUT VARIABLES: | ||
| ! - The day-of-year dependence of the model only uses the DAY argument. If | ||
| ! a continuous day-of-year dependence is desired, this argument should | ||
| ! include the fractional day (e.g., DAY = <day of year> + UTSEC/86400.0 | ||
| ! - If lzalt_type = .true. (default) in the MSISINIT call, then Z is | ||
| ! treated as geodetic altitude. | ||
| ! If lzalt_type = .false., then Z is treated as geopotential height. | ||
| ! - F107 and F107A values are the 10.7 cm radio flux at the Sun-Earth | ||
| ! distance, not the radio flux at 1 AU. | ||
| ! | ||
| ! OUTPUT VARIABLES: | ||
| ! TN Temperature at altitude (K) | ||
| ! DN(1) Total mass density (kg/m3) | ||
| ! DN(2) N2 number density (m-3) | ||
| ! DN(3) O2 number density (m-3) | ||
| ! DN(4) O number density (m-3) | ||
| ! DN(5) He number density (m-3) | ||
| ! DN(6) H number density (m-3) | ||
| ! DN(7) Ar number density (m-3) | ||
| ! DN(8) N number density (m-3) | ||
| ! DN(9) Anomalous oxygen number density (m-3) | ||
| ! DN(10) NO number density (m-3) | ||
| ! TEX Exospheric temperature (K) (optional argument) | ||
| ! | ||
| ! NOTES ON OUTPUT VARIABLES: | ||
| ! - Missing density values are returned as 9.999e-38 | ||
| ! - Species included in mass density calculation are set in MSISINIT | ||
| ! | ||
| !!! ========================================================================= | ||
|
|
||
| !************************************************************************************************** | ||
| ! MSIS_CALC Module: Contains main MSIS entry point | ||
| !************************************************************************************************** | ||
| module msis_calc | ||
|
|
||
| contains | ||
|
|
||
| !================================================================================================== | ||
| ! MSISCALC: The main MSIS subroutine entry point | ||
| !================================================================================================== | ||
| subroutine msiscalc(day,utsec,z,lat,lon,sfluxavg,sflux,ap,tn,dn,tex) | ||
|
|
||
| use msis_constants, only : rp, dmissing, lnp0, Mbarg0divkB, kB, nspec, nodesTN, nd, zetaF, zetaB, & | ||
| Hgamma, zetagamma, maxnbf | ||
| use msis_init, only : msisinit, initflag, zaltflag, specflag, massflag, masswgt, etaTN | ||
| use msis_gfn, only : globe | ||
| use msis_tfn, only : tnparm, tfnparm, tfnx | ||
| use msis_dfn, only : dnparm, dfnparm, dfnx | ||
| use msis_utils, only : alt2gph, bspline, dilog | ||
|
|
||
| implicit none | ||
|
|
||
| real(kind=rp), intent(in) :: day | ||
| real(kind=rp), intent(in) :: utsec | ||
| real(kind=rp), intent(in) :: z | ||
| real(kind=rp), intent(in) :: lat | ||
| real(kind=rp), intent(in) :: lon | ||
| real(kind=rp), intent(in) :: sfluxavg,sflux,ap(1:7) | ||
| real(kind=rp), intent(out) :: tn, dn(1:10) | ||
| real(kind=rp), intent(out), optional :: tex | ||
|
|
||
| real(kind=rp), save :: lastday = -9999.0 | ||
| real(kind=rp), save :: lastutsec = -9999.0 | ||
| real(kind=rp), save :: lastlat = -9999.0 | ||
| real(kind=rp), save :: lastlon = -9999.0 | ||
| real(kind=rp), save :: lastz = -9999.0 | ||
| real(kind=rp), save :: lastsflux = -9999.0 | ||
| real(kind=rp), save :: lastsfluxavg = -9999.0 | ||
| real(kind=rp), save :: lastap(1:7) = -9999.0 | ||
| real(kind=rp), save :: gf(0:maxnbf-1) | ||
| real(kind=rp), save :: Sz(-5:0,2:6) | ||
| integer, save :: iz | ||
| type(tnparm), save :: tpro | ||
| type(dnparm), save :: dpro(1:nspec-1) | ||
|
|
||
| real(8) :: zaltd, latd | ||
| real(kind=rp) :: zeta, lndtotz, Vz, Wz, HRfact, lnPz, delz | ||
| integer :: i, j, kmax, ispec | ||
|
|
||
| ! Check if model has been initialized; if not, perform default initialization | ||
| if (.not. initflag) call msisinit() | ||
|
|
||
| ! Calculate geopotential height, if necessary | ||
| if(zaltflag) then | ||
| zaltd = dble(z) | ||
| latd = dble(lat) | ||
| zeta = alt2gph(latd,zaltd) | ||
| else | ||
| zeta = z | ||
| endif | ||
|
|
||
| ! If only altitude changes then update the local spline weights | ||
| if (zeta .lt. zetaB) then | ||
| if (zeta .ne. lastz) then | ||
| if (zeta .lt. zetaF) then | ||
| kmax = 5 | ||
| else | ||
| kmax = 6 | ||
| endif | ||
| call bspline(zeta,nodesTN,nd+2,kmax,etaTN,Sz,iz) | ||
| lastz = zeta | ||
| endif | ||
| endif | ||
|
|
||
| ! If location, time, or solar/geomagnetic conditions change then recompute the profile parameters | ||
| if ((day .ne. lastday) .or. (utsec .ne. lastutsec) .or. & | ||
| (lat .ne. lastlat) .or. (lon .ne. lastlon) .or. & | ||
| (sflux .ne. lastsflux) .or. (sfluxavg .ne. lastsfluxavg) .or. & | ||
| any(ap .ne. lastap)) then | ||
| call globe(day,utsec,lat,lon,sfluxavg,sflux,ap,gf) | ||
| call tfnparm(gf,tpro) | ||
| do ispec = 2, nspec-1 | ||
| if (specflag(ispec)) call dfnparm(ispec,gf,tpro,dpro(ispec)) | ||
| enddo | ||
| lastday = day | ||
| lastutsec = utsec | ||
| lastlat = lat | ||
| lastlon = lon | ||
| lastsflux = sflux | ||
| lastsfluxavg = sfluxavg | ||
| lastap = ap | ||
| endif | ||
|
|
||
| ! Exospheric temperature | ||
| if (present(tex)) then | ||
| tex = tpro%tex | ||
| endif | ||
|
|
||
| ! Temperature at altitude | ||
| tn = tfnx(zeta,iz,Sz(-3:0,4),tpro) | ||
|
|
||
| ! Temperature integration terms at altitude, total number density | ||
| delz = zeta - zetaB | ||
| if (zeta .lt. zetaF) then | ||
| i = max(iz-4,0) | ||
| if (iz .lt. 4) then | ||
| j = -iz | ||
| else | ||
| j = -4 | ||
| endif | ||
| Vz = dot_product(tpro%beta(i:iz),Sz(j:0,5)) + tpro%cVS | ||
| Wz = 0.0_rp | ||
| lnPz = lnP0 - Mbarg0divkB*(Vz - tpro%Vzeta0) | ||
| lndtotz = lnPz - log(kB*tn) | ||
| else | ||
| if (zeta .lt. zetaB) then | ||
| Vz = dot_product(tpro%beta(iz-4:iz),Sz(-4:0,5)) + tpro%cVS | ||
| Wz = dot_product(tpro%gamma(iz-5:iz),Sz(-5:0,6)) + tpro%cVS*delz + tpro%cWS | ||
| else | ||
| Vz = (delz + log(tn/tpro%tex)/tpro%sigma)/tpro%tex + tpro%cVB | ||
| Wz = (0.5_rp*delz*delz + dilog(tpro%b*exp(-tpro%sigma*delz))/tpro%sigmasq)/tpro%tex & | ||
| + tpro%cVB*delz + tpro%cWB | ||
| endif | ||
| endif | ||
|
|
||
| ! Species number densities at altitude | ||
| HRfact = 0.5_rp * (1.0_rp + tanh(Hgamma*(zeta - zetagamma))) !Reduction factor for chemical/dynamical correction scale height below zetagamma | ||
| do ispec = 2, nspec-1 | ||
| if (specflag(ispec)) then | ||
| dn(ispec) = dfnx(zeta,tn,lndtotz,Vz,Wz,HRfact,tpro,dpro(ispec)) | ||
| else | ||
| dn(ispec) = dmissing | ||
| endif | ||
| enddo | ||
|
|
||
| ! Mass density | ||
| if (specflag(1)) then | ||
| dn(1) = dot_product(dn,masswgt) | ||
| else | ||
| dn(1) = dmissing | ||
| endif | ||
|
|
||
| return | ||
|
|
||
| end subroutine msiscalc | ||
|
|
||
| end module msis_calc |
Oops, something went wrong.