diff --git a/.gitattributes b/.gitattributes new file mode 100644 index 0000000..44a9715 --- /dev/null +++ b/.gitattributes @@ -0,0 +1,20 @@ +.gitattributes text eol=lf +.gitignore text eol=lf +Makefile text eol=lf +*.yml text eol=lf +LICENSE text eol=lf +*.ipynb text eol=lf +*.txt text eol=lf +*.py text eol=lf +*.sh text eol=lf +*.c text eol=lf +*.cpp text eol=lf +*.f text eol=lf +*.f90 text eol=lf +*.for text eol=lf +*.md text eol=lf +*.rst text eol=lf +*.csv text eol=lf +*.m text eol=lf +*.grc text eol=lf +*.pas text eol=lf diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..e28ece6 --- /dev/null +++ b/.gitignore @@ -0,0 +1,17 @@ +*.out +fort.12 +driver +*.a +*.o +*.exe +*.mod +*.so +CMakeFiles/ +CMakeCache.txt +cmake_install.cmake +Makefile +build/ +*.egg-info/ +basic +__pycache__/ + diff --git a/.travis.yml b/.travis.yml new file mode 100644 index 0000000..90a0d4b --- /dev/null +++ b/.travis.yml @@ -0,0 +1,55 @@ +language: generic + +os: + - linux + - osx + +dist: trusty +group: edge + +env: TRAVIS_PYTHON_VERSION=3.6 FC=gfortran + +notifications: + email: false + +git: + depth: 3 + +addons: + apt: + packages: + - gfortran + +before_install: + - if [[ $TRAVIS_OS_NAME == osx ]]; then + wget https://repo.continuum.io/miniconda/Miniconda3-latest-MacOSX-x86_64.sh -O miniconda.sh; + brew update; + brew install gcc; + else + wget https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh -O miniconda.sh; + fi + + - bash miniconda.sh -b -p $HOME/miniconda + - export PATH=$HOME/miniconda/bin:$PATH + - hash -r + + - conda config --set always_yes yes --set changeps1 no + - conda update -q conda + - conda create -n test python=$TRAVIS_PYTHON_VERSION + - source activate test + - pip -q install coveralls + +install: + - python setup.py develop + +script: + # test fortran-only "basic" program + - cd build + - cmake .. + - make + - cd .. + # separately (previous 4 lines not needed) test python access to Glow + - coverage run tests/test.py -v + +after_success: coveralls + diff --git a/CMakeLists.txt b/CMakeLists.txt new file mode 100644 index 0000000..4cd46cc --- /dev/null +++ b/CMakeLists.txt @@ -0,0 +1,45 @@ +cmake_minimum_required (VERSION 2.8.12) +project(glow Fortran) + +set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ..) +set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ..) + +add_compile_options(-O3 -g) + +if(${CMAKE_Fortran_COMPILER_ID} STREQUAL "Intel") + add_compile_options(-debug full -traceback -fpe-all=0 -check all) +else() + add_compile_options(-mtune=native -Warray-bounds -fexceptions -fbacktrace) +endif() +# --------- GLOW library ------------------- +add_library(glow cglow.f90 glow.f90 bands.f90 conduct.f90 egrid.f90 ephoto.f90 etrans.f90 exsect.f fieldm.f gchem.f90 geomag.f90 maxt.f90 mzgrid.f90 qback.f90 rcolum.f90 rout.f90 snoem.f90 snoemint.f90 solzen.f90 ssflux.f90 iri90.f nrlmsise00.f) +#------------- basic example ------------ +add_executable(basic glowbasic.f90) +target_link_libraries(basic glow) +#include_directories(../build) # for .mod +#---------- NetCDF (optional) ------------ +find_package(netCDF) +if(netCDF_DIR) + include_directories(${netCDF_INCLUDE_DIR}) + + add_library(ionetcdf readtgcm.f90 output.f90 ) + target_link_libraries(ionetcdf netcdff glow) # yes double "ff" +else() + message(WARNING "NetCDF failed to be included " ${netCDF_INCLUDE_DIR}) +endif() +#---------- MPI (optional) ------------- +find_package(MPI) +if(MPI_Fortran_COMPILER) + # wrapper compiler mpif90 + #set(CMAKE_Fortran_COMPILER ${MPI_Fortran_COMPILER}) + + add_definitions(${MPI_Fortran_COMPILE_FLAGS}) + include_directories(${MPI_Fortran_INCLUDE_PATH}) + add_compile_options(-pthread) + + add_executable(driver glowdriver.f90 tzgrid.f90) + target_link_libraries(driver glow ionetcdf ${MPI_Fortran_LIBRARIES}) +else() + message(WARNING "MPI failed to be included " ${MPI_Fortran_LIBRARIES}) +endif() + diff --git a/Makefile.glowbasic b/Makefile.glowbasic index bd131c2..443f884 100644 --- a/Makefile.glowbasic +++ b/Makefile.glowbasic @@ -1,5 +1,5 @@ # -FC = ifort +#FC = ifort FFLAGS = -O3 #FFLAGS = -g $(DBGFLAGS) diff --git a/Makefile.glowdriver b/Makefile.glowdriver index 6206866..6750026 100644 --- a/Makefile.glowdriver +++ b/Makefile.glowdriver @@ -1,7 +1,7 @@ # FC = mpif90 - FFLAGS = -fc=ifort -O3 -I $(INC_NETCDF) + FFLAGS = -O3 -I $(INC_NETCDF) #FFLAGS = -fc=ifort -g $(DBGFLAGS) -I $(INC_NETCDF) DBGFLAGS = -debug full -traceback @@ -17,17 +17,17 @@ DBGFLAGS += -fpe-all=0 # this traps all floating point exceptions # # Sources (in order of dependency): # -SOURCES = cglow.f90 readtgcm.f90 output.f90 glowdriver.f90 glow.f90 bands.f90 conduct.f90 egrid.f90 ephoto.f90 etrans.f90 exsect.f fieldm.f gchem.f90 geomag.f90 maxt.f90 mzgrid.f90, qback.f90 rcolum.f90 rout.f90 snoem.f90 snoemint.f90 solzen.f90 ssflux.f90 tzgrid.f90, iri90.f nrlmsise00.f +SOURCES = cglow.f90 mpicdf/readtgcm.f90 mpicdf/output.f90 mpicdf/glowdriver.f90 glow.f90 bands.f90 conduct.f90 egrid.f90 ephoto.f90 etrans.f90 exsect.f fieldm.f gchem.f90 geomag.f90 maxt.f90 mzgrid.f90, qback.f90 rcolum.f90 rout.f90 snoem.f90 snoemint.f90 solzen.f90 ssflux.f90 mpicdf/tzgrid.f90, iri90.f nrlmsise00.f OBJS := $(addsuffix .o, $(basename $(SOURCES))) EXEC = glow.exe $(EXEC): $(OBJS) - $(FC) -o -fc=ifort $@ $(OBJS) $(LIBS) $(LDFLAGS) + $(FC) -o $@ $(OBJS) $(LIBS) $(LDFLAGS) LIB_NETCDF = /home/tgcm/intel/netcdf-4.1.1/lib INC_NETCDF = /home/tgcm/intel/netcdf-4.1.1/include -LIBS = -L /usr/lib64 -lcurl -L$(LIB_NETCDF) -lnetcdf +LIBS = -L /usr/lib64 -lcurl -L$(LIB_NETCDF) -lnetcdff clean: rm -f *.o *.mod $(EXEC) diff --git a/README.md b/README.md index 25887f9..de260dd 100644 --- a/README.md +++ b/README.md @@ -1,10 +1,90 @@ +[![Build Status](https://travis-ci.org/scivision/GLOW.svg?branch=master)](https://travis-ci.org/scivision/GLOW) + # GLOW The GLobal airglOW Model This directory contains: - Fortran-90 source code files, - Makefiles, - Example input and output files, - Example job script, - Subdirectory data/ contains input data files, - Subdirectory data/iri90 contains IRI input data files + + * Modern Fortran source code files, + * Makefile + * CMakeLists.txt + * Example input and output files + * Example job script + * `data/` contains input data files, + * `data/iri90/` contains IRI input data files + * [Release Notes](ReleaseNotes.rst) + + +## Python install +This command automatically compiles the Fortran code to access from Python on any platform. + + python setup.py develop + +You can then run the self-tests with + + python tests/test.py -v + +## Fortran-only mode +While many users use the Python interface to Glow, users on HPCC may want to use MPI directly in Fortran using the Makefiles or CMake. Here's how to compile the ``basic`` example with CMake. + + cd build + cmake .. + make + cd .. + +This creates: + + +executable | description +------------|-------------- +basic | basic GLOW +driver | MPI/NetCDF GLOW + +### Fortran examples +With regard to precision, at first try I occasionally see the least digit of precision in the text output files differ. +This can be related to Stan using Intel Fortran compiler and I'm using Gfortran. +E.g. I get 4.34e-9 and Stan got 4.33e-9. + +I also get the message:: + + Note: The following floating-point exceptions are signalling: IEEE_UNDERFLOW_FLAG IEEE_DENORMAL + +Which is a warning sign that different compilers and platforms may yield different results of unpredictable variance. + + + +#### Aurora example (night) +Compare your results vs. Stan's results with `meld`: + + ./basic < in.basic.aur > aur.basic.out + + meld out.basic.aur aur.basic.out + + +#### Aurora Example (night) + + ./basic < in.basic.day > day.basic.out + + meld out.basic.day day.basic.out + +### MPI Prereq +Allows parallel execution of GLOW Fortran code for HPCC. + + sudo apt install libopenmpi-dev + +### NetCDF prereq +allows reading/writing data in NetCDF (optional). +You need both of these. + + sudo apt install libnetcdf-dev libnetcdff-dev + +### Select Fortran compiler +Simply use the variable `FC`. Example + + cd build/ + rm -r * + + FC=ifort cmake .. + make + +NOTE: Using the Intel compiler requires that you have [built NetCDF using Intel Fortran](https://software.intel.com/en-us/articles/performance-tools-for-software-developers-building-netcdf-with-the-intel-compilers/)--this is an issue ANYTIME you use the Intel Compiler and NetCDF. diff --git a/ReleaseNotes.rst b/ReleaseNotes.rst new file mode 100644 index 0000000..c5c4b84 --- /dev/null +++ b/ReleaseNotes.rst @@ -0,0 +1,138 @@ +GLOW version 0.98 release notes + +:author: Stan Solomon +:date: 3/2017 + +.. contents:: + +See Quickstart.txt for very basic instructions +See Glow.txt for more details +See Glowlicense.txt for the Open Source Academic Research License Agreement + +This is now Modern Fortran code, with MPI and netCDF options. +No more common blocks. +No more header file. +Use-associated variables are defined in module ``cglow.f90`` +Most code converted to lower-case free-form style + +Example driver program glowbasic provides functionality in previous versions +Can run one profile at a time or loop +Single-processor, no MPI required +No netCDF required either +You can easily modify this program to suit particular applications + +New driver program glowdriver uses TIE-GCM, TIME-GCM or empirical model inputs +Ultimately this will include WACCM-X as well +Fully MPI code with user-specified number of processors +Namelist input from file in.namelist +netCDF output to a file specified in the namelist input + +Can still run with empirical models: +Default empirical neutral model updated to MSIS-2K +NOEM empirical nitric oxide model, based on SNOE data, used for empirical runs +IRI-90 still used for empirical ionosphere + +New features in v. 0.98 +======================= + +* Most of the new things in v. 0.98 are computational, not scientific. +* Move to MPI parallel processing forced the elimination of common blocks. +* Reads netCDF input files from GCMs, and writes global netCDF output files. +* All of this is optional - you can still run the model the "old" way. +(But old driver programs will have to be updated to use shared variables.) +* Empirical model calls are consolidated in subroutine MZGRID +* Calculates conductivities using Ryan McGranaghan's CONDUCT routine (adapted). +* N2 LBH, O 1356, N 1493, O 1304 added to airglow output (ZETA array). +(Note 1493 is very approximate since branching rations are poorly known.) +* Radiative recombination and mutual neutralization included in 1356. +* Subroutine BANDS added to allocate LBH among v' levels. +(Really just a stub now, will ultimately handle vibrational distributions.) + +Note that several bugfixes and parameter adjustments were made following the + "beta" release of v. 0.98 in early January 2017, so do not use that version. + +Other significant 0.98 changes +------------------------------- + +The quartic equation solver formerly used by GCHEM to calculate electron density has been abandoned because it produced minor errors and could not be verified. +Instead, a simple iterative method is used, which turns out to be more accurate and equally fast. +Thanks to Michael Hirsch, Ph.D. for identifying this problem. + +The high-energy electron precipitation example (hexexample) is no longer provided. +This functionality is built-in to the codes, but if you want to do a high-energy electron run you will have to increase NBINS in cglow.f90 (e.g., to 343 for energies up to 100 MeV) and extend the altitude grid to lower altitude by increasing JMAX in the driver program, and adding levels in mxgrid.f90. + +Conjugate hemisphere photoelectrons are no longer provided in the example drivers, since they are very minor, triple the run time, and weren't widely used. +They can easily be re-installed in a driver program, but note that the magnetic field needs to be updated. +The calculation is fairly approximate anyway because it doesn't consider exosphere/plasmasphere attenuation of the electron flux. + +The background ("night time") ionization rate calculated in QBACK was updated to conform to the forumlation now used in the TIE-GCM, TIME-GCM, WACCM, and WACCM-X. +It is still pretty approximate, but produces an night ionosphere in the E and F1 regions that is in reasonable agreement with observations. + +Several unused artifacts were eliminated, but arrays PIA and SESPEC were retained in anticipation of including proton aurora ionization in the future. + +Known issues v.0.98 +------------------- + (most of these carried over from v. 0.97): + +* X-rays shortward of 18 A need to be re-examined and updated. +* Magnetic field (GEOMAG, FIELDM) is out of date. GEOMAG is only used by the NOEM model, and isn't too different from the coordinate transform used to generate SNOE data on a magnetic grid. FIELDM only calculates the dip angle, which has a minimual effect on results. But these should be updated to time-varying IGRF/Apex coordinates. +* Temporary Y2K fix to SUNCOR. Also should update to J2000 epoch. Should be fine for 1950-2050, but really need to make if valid for either yyddd to yyyyddd date format, which will extend the range of validity to 1900-2100. (This only affects SZA.) +* O(1S) needs to be re-evaluated (still). Nightglow recombination source still not included. Also need to include O2 A-bands. +* O(1S) from O2 dissociation (BSO2) is currently hardwired - only works with LMAX=123 (although changes shortward of 800 A are OK). +* IRI should be updated. Usually, IRI is Only used for electron density above 200 km, and IRI hasn't changed much in the F-region, so it's OK for now. +* Cascade contributions to 7774 and subsequent cascade to 1356 is questionable. I reduced the effective 7774 cross section, so 7774 and 1356 are now self-consistent, and in reasonable agreement with GUVI data. +* 1493 sources and branching ratios are speculative. For now, I am presuming that it is produced during photodissociative ionization and electron impact dissociative ionization of N2. + +Version 0.973 release notes +=========================== +Stan Solomon, 3/2015 + +Version 0.973 is an incremental release of GLOW, mostly just the example drivers + +* updates example drivers to MSIS-2K +* adds NOEM empirical nitric oxide model, based on SNOE data, to example drivers +* some cleanup of example drivers +* fixs ssflux so it only reads file on first call (or if ISCALE changes). + + + +Version 0.97 release notes +========================== +Stan Solomon, 4/2005 + +New features: + +* Relativistic correction to electron impact cross sections included +* Maxwellian or monoenergetic fluxes generated by MAXT +* Possible to use any solar grid by changing only input files +* Photoabsorption and photoionization cross sections supplied in files +* SSFLUX completely re-written: + * Model parameters supplied in files + * Default is ~1 nm grid (5 nm in FUV) + * Hinteregger model still there (ISCALE=0) + * EUVAC also available (ISCALE=1) + * User grid and input supported by changing input file and LMAX (ISCALE=2) +* Common block CGLOW is unchanged + * Should facilitate upgrade path for existing programs + * But there are now several obsolete artifacts + +Issues addressed in v. 0.97 +--------------------------- + +* Fixed two problems with Auger electron production +* Fixed some small bugs in O(1S) calculation +* Adjusted N(2D)+O rate coefficient to Fell et al. value (6.9e-13) +* Adjusted C III on N2 cross section to fix O2 ionization rate problem +* Removed various artifacts, including EAURI +* Removed unnecessary N(2D) initial guess (now just set to zero) +* Now use standard energy and altitude grid in both day and aurora examples + +Known issues +============ + +* X-rays shortward of 18 A need to be re-examined and updated. +* Magnetic field (GEOMAG, FIELDM) is out of date, really need to update to IGRF, but at the resolutions typical here should be OK for now. +* Temporary Y2K fix to SUNCOR. Should be fine for 1950-2050, but really need to change from yyddd to yyyyddd date format, which will enable range of validity to extend from 1900-2100. (This only affects SZA.) +* O(1S) needs to be re-evaluated (still). +* O(1S) from O2 dissociation (BSO2) is currently hardwired - only works with LMAX=123 (although changes shortward of 800 A are OK). +* IRI should be updated. diff --git a/Releasenotes.txt b/Releasenotes.txt deleted file mode 100644 index 58c92d9..0000000 --- a/Releasenotes.txt +++ /dev/null @@ -1,146 +0,0 @@ -GLOW version 0.98 release notes, Stan Solomon, 3/2017 - -See Quickstart.txt for very basic instructions -See Glow.txt for more details -See Glowlicense.txt for the Open Source Academic Research License Agreement - -This is now a Fortran-90 code, with MPI and netCDF options - No more common blocks. - No more header file. - Use-associated variables are defined in module cglow.f90 - Most code converted to lower-case free-form style -Example driver program glowbasic provides functionality in previous versions - Can run one profile at a time or loop - Single-processor, no MPI required - No netCDF required either - You can easily modify this program to suit particular applications -New driver program glowdriver uses TIE-GCM, TIME-GCM or empirical model inputs - Ultimately this will include WACCM-X as well - Fully MPI code with user-specified number of processors - Namelist input from file in.namelist - netCDF output to a file specified in the namelist input -Can still run with empirical models: - Default empirical neutral model updated to MSIS-2K - NOEM empirical nitric oxide model, based on SNOE data, used for empirical runs - IRI-90 still used for empirical ionosphere - -New features in v. 0.98: - Most of the new things in v. 0.98 are computational, not scientific. - Move to MPI parallel processing forced the elimination of common blocks. - Reads netCDF input files from GCMs, and writes global netCDF output files. - All of this is optional - you can still run the model the "old" way. - (But old driver programs will have to be updated to use shared variables.) - Empirical model calls are consolidated in subroutine MZGRID - Calculates conductivities using Ryan McGranaghan's CONDUCT routine (adapted). - N2 LBH, O 1356, N 1493, O 1304 added to airglow output (ZETA array). - (Note 1493 is very approximate since branching rations are poorly known.) - Radiative recombination and mutual neutralization included in 1356. - Subroutine BANDS added to allocate LBH among v' levels. - (Really just a stub now, will ultimately handle vibrational distributions.) - -Note that several bugfixes and parameter adjustments were made following the - "beta" release of v. 0.98 in early January 2017, so do not use that version. - -Other significant changes: - - The quartic equation solver formerly used by GCHEM to calculate electron -density has been abandoned because it produced minor errors and could not -be verified. Instead, a simple iterative method is used, which turns out -to be more accurate and equally fast. Thanks to Matt Hirsch for identifying -this problem. - - The high-energy electron precipitation example (hexexample) is no longer -provided. This functionality is built-in to the codes, but if you want to -do a high-energy electron run you will have to increase NBINS in cglow.f90 -(e.g., to 343 for energies up to 100 MeV) and extend the altitude grid to -lower altitude by increasing JMAX in the driver program, and adding levels -in mxgrid.f90. - - Conjugate hemisphere photoelectrons are no longer provided in the example -drivers, since they are very minor, triple the run time, and weren't widely -used. They can easily be re-installed in a driver program, but note that -the magnetic field needs to be updated. The calculation is fairly approximate -anyway because it doesn't consider exosphere/plasmasphere attenuation of -the electron flux. - - The background ("night time") ionization rate calculated in QBACK was updated -to conform to the forumlation now used in the TIE-GCM, TIME-GCM, WACCM, and -WACCM-X. It is still pretty approximate, but produces an night ionosphere in -the E and F1 regions that is in reasonable agreement with observations. - - Several unused artifacts were eliminated, but arrays PIA and SESPEC were -retained in anticipation of including proton aurora ionization in the future. - -Known issues (most of these carried over from v. 0.97): - X-rays shortward of 18 A need to be re-examined and updated. - Magnetic field (GEOMAG, FIELDM) is out of date. GEOMAG is only used by - the NOEM model, and isn't too different from the coordinate transform - used to generate SNOE data on a magnetic grid. FIELDM only calculates - the dip angle, which has a minimual effect on results. But these should - be updated to time-varying IGRF/Apex coordinates. - Temporary Y2K fix to SUNCOR. Also should update to J2000 epoch. - Should be fine for 1950-2050, but really need to make if valid for either - yyddd to yyyyddd date format, which will extend the range of validity to - 1900-2100. (This only affects SZA.) - O(1S) needs to be re-evaluated (still). Nightglow recombination source - still not included. Also need to include O2 A-bands. - O(1S) from O2 dissociation (BSO2) is currently hardwired - only works with - LMAX=123 (although changes shortward of 800 A are OK). - IRI should be updated. Usually, IRI is Only used for electron density above - 200 km, and IRI hasn't changed much in the F-region, so it's OK for now. - Cascade contributions to 7774 and subsequent cascade to 1356 is questionable. - I reduced the effective 7774 cross section, so 7774 and 1356 are now - self-consistent, and in reasonable agreement with GUVI data. - 1493 sources and branching ratios are speculative. For now, I am presuming - that it is produced during photodissociative ionization and electron - impact dissociative ionization of N2. - -------------------------------------------------- - -Version 0.973 release notes, Stan Solomon, 3/2015 - -Version 0.973 is an incremental release of GLOW, mostly just the example drivers - updates example drivers to MSIS-2K - adds NOEM empirical nitric oxide model, based on SNOE data, to example drivers - some cleanup of example drivers - fixs ssflux so it only reads file on first call (or if ISCALE changes). - -------------------------------------------------- - -Version 0.97 release notes, Stan Solomon, 4/2005 - -New features: - Relativistic correction to electron impact cross sections included - Maxwellian or monoenergetic fluxes generated by MAXT - Possible to use any solar grid by changing only input files - Photoabsorption and photoionization cross sections supplied in files - SSFLUX completely re-written: - Model parameters supplied in files - Default is ~1 nm grid (5 nm in FUV) - Hinteregger model still there (ISCALE=0) - EUVAC also available (ISCALE=1) - User grid and input supported by changing input file and LMAX (ISCALE=2) - Common block CGLOW is unchanged - Should facilitate upgrade path for existing programs - But there are now several obsolete artifacts - -Issues addressed in v. 0.97: - Fixed two problems with Auger electron production - Fixed some small bugs in O(1S) calculation - Adjusted N(2D)+O rate coefficient to Fell et al. value (6.9e-13) - Adjusted C III on N2 cross section to fix O2 ionization rate problem - Removed various artifacts, including EAURI - Removed unnecessary N(2D) initial guess (now just set to zero) - Now use standard energy and altitude grid in both day and aurora examples - -Known issues: - X-rays shortward of 18 A need to be re-examined and updated. - Magnetic field (GEOMAG, FIELDM) is out of date, really need to update to - IGRF, but at the resolutions typical here should be OK for now. - Temporary Y2K fix to SUNCOR. Should be fine for 1950-2050, but really - need to change from yyddd to yyyyddd date format, which will enable - range of validity to extend from 1900-2100. (This only affects SZA.) - O(1S) needs to be re-evaluated (still). - O(1S) from O2 dissociation (BSO2) is currently hardwired - only works with - LMAX=123 (although changes shortward of 800 A are OK). - IRI should be updated. diff --git a/bin/.ignore b/bin/.ignore new file mode 100644 index 0000000..e69de29 diff --git a/build/.ignore b/build/.ignore new file mode 100644 index 0000000..e69de29 diff --git a/cglow.f90 b/cglow.f90 index c4126b8..2e279ae 100644 --- a/cglow.f90 +++ b/cglow.f90 @@ -50,6 +50,8 @@ module cglow implicit none save + real, parameter :: pi = atan(1.)*4. + ! Array dimensions, configurable: integer :: jmax ! number of vertical levels @@ -188,6 +190,24 @@ subroutine cglow_init end subroutine cglow_init + + Elemental Real Function sind(ang) + ! this is a non-standard function so include it explicity for cross-platform compatibility + real, intent(in) :: ang + + sind = sin(ang*pi/180.) + + End Function sind + + + Elemental Real Function cosd(ang) + ! this is a non-standard function so include it explicity for cross-platform compatibility + real, intent(in) :: ang + + cosd = cos(ang*pi/180.) + + End Function cosd + !----------------------------------------------------------------------- end module cglow diff --git a/ephoto.f90 b/ephoto.f90 index 5e6e202..143ff9a 100644 --- a/ephoto.f90 +++ b/ephoto.f90 @@ -115,7 +115,7 @@ subroutine ephoto ifirst = 0 filepath = trim(data_dir)//'ephoto_xn2.dat' - open(unit=1,file=filepath,status='old',readonly) + open(unit=1,file=filepath,status='old',action='read') read(1,*) read(1,*) read(1,*) @@ -126,7 +126,7 @@ subroutine ephoto close(1) filepath = trim(data_dir)//'ephoto_xo2.dat' - open(unit=1,file=filepath,status='old',readonly) + open(unit=1,file=filepath,status='old',action='read') read(1,*) read(1,*) read(1,*) @@ -137,7 +137,7 @@ subroutine ephoto close(1) filepath = trim(data_dir)//'ephoto_xo.dat' - open(unit=1,file=filepath,status='old',readonly) + open(unit=1,file=filepath,status='old',action='read') read(1,*) read(1,*) read(1,*) diff --git a/glow.f90 b/glow.f90 index ffa886d..47198b7 100644 --- a/glow.f90 +++ b/glow.f90 @@ -220,6 +220,4 @@ subroutine glow call bands - return - end subroutine glow diff --git a/glowbasic.f90 b/glowbasic.f90 index e062237..efae913 100644 --- a/glowbasic.f90 +++ b/glowbasic.f90 @@ -140,10 +140,12 @@ program glowbasic ! Output section: ! write(6,"(1x,i7,9f8.1)") idate,ut,glat,glong,f107a,f107,f107p,ap,ef,ec - write(6,"(' Z Tn O N2 NO Ne(in) Ne(out) Ionrate O+ O2+ NO+ N(2D) Pederson Hall')") + write(6,"(' Z Tn O N2 NO Ne(in) & + Ne(out) Ionrate O+ O2+ NO+ N(2D) Pederson Hall')") write(6,"(1x,0p,f5.1,f6.0,1p,12e10.2)") (z(j),ztn(j),zo(j),zn2(j),zno(j),ze(j), & ecalc(j),tir(j),zxden(3,j),zxden(6,j),zxden(7,j),zxden(10,j),pedcond(j),hallcond(j),j=1,jmax) - write(6,"(' Z 3371 4278 5200 5577 6300 7320 10400 3644 7774 8446 3726 LBH 1356 1493 1304')") + write(6,"(' Z 3371 4278 5200 5577 6300 7320 & + 10400 3644 7774 8446 3726 LBH 1356 1493 1304')") write(6,"(1x,f5.1,15f8.2)")(z(j),(zeta(ii,j),ii=1,15),j=1,jmax) enddo diff --git a/glowdriver.f90 b/glowdriver.f90 index 0792488..5110d3c 100644 --- a/glowdriver.f90 +++ b/glowdriver.f90 @@ -425,7 +425,7 @@ program glowdriver ! Create and define a new netCDF output file for each time (root task only): ! if (itask == 0) then - write (ifile,"('.',i3.3,'.nc')"),itime + write (ifile,"('.',i3.3,'.nc')") itime glow_ncfileit = trim(glow_ncfile) // ifile call create_ncfile(glow_ncfileit,tgcm_ncfile) call write_ncfile(glow_ncfileit) diff --git a/glowiono/__init__.py b/glowiono/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/nrlmsise00.f b/nrlmsise00.f index 7040a10..4cc07da 100644 --- a/nrlmsise00.f +++ b/nrlmsise00.f @@ -1657,6 +1657,7 @@ FUNCTION CCOR2(ALT, R,H1,ZH,H2) C----------------------------------------------------------------------- BLOCK DATA GTD7BK C MSISE-00 01-FEB-02 + CHARACTER(len=4) :: ISDATE,ISTIME,NAME COMMON/PARM7/PT1(50),PT2(50),PT3(50),PA1(50),PA2(50),PA3(50), $ PB1(50),PB2(50),PB3(50),PC1(50),PC2(50),PC3(50), $ PD1(50),PD2(50),PD3(50),PE1(50),PE2(50),PE3(50), diff --git a/output.f90 b/output.f90 index 1118dbf..0d9264d 100644 --- a/output.f90 +++ b/output.f90 @@ -14,7 +14,7 @@ module output ! any routines in this module are called. They are currently ! set by the main driver (glowdriver.f90) and module cglow.f90. - + use, intrinsic :: iso_fortran_env, only: error_unit use netcdf use cglow,only: nbins,lmax,nmaj,nei,nex,nw,nc,nst @@ -110,10 +110,10 @@ subroutine create_ncfile(ncfile,tgcm_ncfile) istat = nf90_create(ncfile,NF90_CLOBBER,ncid) if (istat /= NF90_NOERR) then - write(6,"('>>> Error creating netcdf output file ',a)") trim(ncfile) + write(error_unit,"('>>> Error creating netcdf output file ',a)") trim(ncfile) call handle_ncerr(istat,'Error from nf90_create',1) endif - write(6,"('Created netcdf dataset ',a)") trim(ncfile) + write(error_unit,"('Created netcdf dataset ',a)") trim(ncfile) ! ! Define dimensions: ! @@ -326,7 +326,7 @@ subroutine write_ncfile(ncfile) istat = nf90_open(ncfile,NF90_WRITE,ncid) if (istat /= NF90_NOERR) then - write(6,"('>>> Error opening file for writing: ncfile=',a)") trim(ncfile) + write(error_unit,"('>>> Error opening file for writing: ncfile=',a)") trim(ncfile) call handle_ncerr(istat,'Error from nf90_open',1) endif @@ -397,9 +397,9 @@ subroutine datetime(curdate,curtime) curtime = ' ' call date_and_time(date,time,zone,values) ! -! write(6,"('datetime: date=',a,' time=',a,' zone=',a)") +! write(error_unit,"('datetime: date=',a,' time=',a,' zone=',a)") ! | date,time,zone -! write(6,"('datetime: values=',8i8)") values +! write(error_unit,"('datetime: values=',8i8)") values ! curdate(1:2) = date(5:6) curdate(3:3) = '/' @@ -424,13 +424,14 @@ subroutine handle_ncerr(istat,msg,ifatal) integer,intent(in) :: istat,ifatal character(len=*),intent(in) :: msg ! - write(6,"(/72('-'))") - write(6,"('>>> Error from netcdf library:')") - write(6,"(a)") trim(msg) - write(6,"('istat=',i5)") istat - write(6,"(a)") nf90_strerror(istat) - write(6,"(72('-')/)") - if (ifatal > 0) stop('Fatal netcdf error') + write(error_unit, "(/72('-'))") + write(error_unit, "('>>> Error from netcdf library:')") + write(error_unit, "(a)") trim(msg) + write(error_unit, "('istat=',i5)") istat + write(error_unit, "(a)") nf90_strerror(istat) + write(error_unit, "(72('-')/)") + if (ifatal > 0) error stop 'Fatal netcdf error' + end subroutine handle_ncerr !----------------------------------------------------------------------- diff --git a/readtgcm.f90 b/readtgcm.f90 index f5e21f4..8545504 100644 --- a/readtgcm.f90 +++ b/readtgcm.f90 @@ -66,7 +66,7 @@ module readtgcm ! iday Model day, ddd ! ut Model universal time, hours - + use iso_fortran_env, only: error_unit use netcdf implicit none @@ -140,7 +140,7 @@ subroutine read_tgcm(ncfile,itime) call handle_ncerr(istat,trim(msg),1) else if (iprint > 0) & - write(6,"(/,'Opened file ',a)") trim(ncfile) + write(error_unit,"(/,'Opened file ',a)") trim(ncfile) endif ! ! Allocate arrays to read 2d and 3d spatial variables: @@ -151,9 +151,9 @@ subroutine read_tgcm(ncfile,itime) ! Read itime history (model time mtime(:,itime)): ! if (itime <= 0.or.itime > ntime) then - write(6,"('>>> read_tgcm: bad itime=',i4,' (ntime=',i4,')')") & + write(error_unit,"('>>> read_tgcm: bad itime=',i4,' (ntime=',i4,')')") & itime,ntime - stop 'itime' + error stop 'itime' endif istat = nf90_inq_varid(ncid,'year',id) @@ -171,7 +171,7 @@ subroutine read_tgcm(ncfile,itime) istat = nf90_inq_varid(ncid,'mtime',id) istat = nf90_get_var(ncid,id,mtime,(/1,itime/),(/3,1/)) - write(6,"(/,'read_tgcm: itime=',i4,' ntime=',i4,' iyear=',i5,' iday=',i4,' ut=',f8.2,' mtime=',3i4)") & + write(error_unit,"(/,'read_tgcm: itime=',i4,' ntime=',i4,' iyear=',i5,' iday=',i4,' ut=',f8.2,' mtime=',3i4)") & itime,ntime,iyear,iday,ut,mtime ! ! Read vars at current history/time: @@ -188,11 +188,11 @@ subroutine read_tgcm(ncfile,itime) f3d = 0. istat = nf90_get_var(ncid,id,f3d,(/1,1,1,itime/),(/nlon,nlat,nlev,1/)) if (istat /= NF90_NOERR) then - write(6,"('>>> Error reading 3d var ')") fnames(n) + write(error_unit,"('>>> Error reading 3d var ')") fnames(n) call handle_ncerr(istat,'Error from nf90_get_var',1) endif if (iprint > 0) & - write(6,"('Read 3d var: itime=',i4,' n=',i4,' fld ',a,' id=',i4,' min,max=',2es12.4)") & + write(error_unit,"('Read 3d var: itime=',i4,' n=',i4,' fld ',a,' id=',i4,' min,max=',2es12.4)") & itime,n,trim(fnames(n)),id,minval(f3d(:,:,1:nlev-1)),maxval(f3d(:,:,1:nlev-1)) ! ! 2d var+time: assume (nlon,nlat): @@ -201,11 +201,11 @@ subroutine read_tgcm(ncfile,itime) f2d = 0. istat = nf90_get_var(ncid,id,f2d,(/1,1,itime/),(/nlon,nlat,1/)) if (istat /= NF90_NOERR) then - write(6,"('>>> Error reading 2d var ')") fnames(n) + write(error_unit,"('>>> Error reading 2d var ')") fnames(n) call handle_ncerr(istat,'Error from nf90_get_var',1) endif if (iprint > 0) & - write(6,"('Read 2d var: itime=',i4,' n=',i4,' fld ',a,' id=',i4,' min,max=',2es12.4)") & + write(error_unit,"('Read 2d var: itime=',i4,' n=',i4,' fld ',a,' id=',i4,' min,max=',2es12.4)") & itime,n,trim(fnames(n)),id,minval(f2d(:,:)),maxval(f2d(:,:)) endif ! 2d or 3d var ! @@ -218,7 +218,7 @@ subroutine read_tgcm(ncfile,itime) tn = f3d tn(:,:,nlev) = tn(:,:,nlev-1) if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(tn),maxval(tn) found(n) = .true. case('UN') @@ -226,7 +226,7 @@ subroutine read_tgcm(ncfile,itime) un = f3d un(:,:,nlev) = un(:,:,nlev-1) if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(un),maxval(un) found(n) = .true. case('VN') @@ -234,42 +234,42 @@ subroutine read_tgcm(ncfile,itime) vn = f3d vn(:,:,nlev) = vn(:,:,nlev-1) if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(vn),maxval(vn) found(n) = .true. case('O2') if (.not.allocated(o2)) allocate(o2(nlon,nlat,nlev)) o2 = f3d if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(o2),maxval(o2) found(n) = .true. case('O1') if (.not.allocated(o1)) allocate(o1(nlon,nlat,nlev)) o1 = f3d if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(o1),maxval(o1) found(n) = .true. case('N2') if (.not.allocated(n2)) allocate(n2(nlon,nlat,nlev)) n2 = f3d if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(n2),maxval(n2) found(n) = .true. case('HE') if (.not.allocated(he)) allocate(he(nlon,nlat,nlev)) he = f3d if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(he),maxval(he) found(n) = .true. case('NO') if (.not.allocated(no)) allocate(no(nlon,nlat,nlev)) no = f3d if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(no),maxval(no) found(n) = .true. case('TI') @@ -277,7 +277,7 @@ subroutine read_tgcm(ncfile,itime) ti = f3d ti(:,:,nlev) = ti(:,:,nlev-1) if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(ti),maxval(ti) found(n) = .true. case('TE') @@ -285,42 +285,42 @@ subroutine read_tgcm(ncfile,itime) te = f3d te(:,:,nlev) = te(:,:,nlev-1) if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(te),maxval(te) found(n) = .true. case('NE') if (.not.allocated(ne)) allocate(ne(nlon,nlat,nlev)) ne = f3d if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(ne),maxval(ne) found(n) = .true. case('Z') if (.not.allocated(z)) allocate(z(nlon,nlat,nlev)) z = f3d if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(z),maxval(z) found(n) = .true. case('ZG') if (.not.allocated(zg)) allocate(zg(nlon,nlat,nlev)) zg = f3d*1.e-5 ! cm->km if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(zg),maxval(zg) found(n) = .true. case('N2D') if (.not.allocated(n2d)) allocate(n2d(nlon,nlat,nlev)) n2d = f3d if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(n2d),maxval(n2d) found(n) = .true. case('N4S') if (.not.allocated(n4s)) allocate(n4s(nlon,nlat,nlev)) n4s = f3d if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(n4s),maxval(n4s) found(n) = .true. ! @@ -330,46 +330,46 @@ subroutine read_tgcm(ncfile,itime) if (.not.allocated(cusp)) allocate(cusp(nlon,nlat)) cusp = f2d if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(cusp),maxval(cusp) found(n) = .true. case('DRIZZLE') if (.not.allocated(drizzle)) allocate(drizzle(nlon,nlat)) drizzle = f2d if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(drizzle),maxval(drizzle) found(n) = .true. case('ALFA') if (.not.allocated(alfa)) allocate(alfa(nlon,nlat)) alfa = f2d if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(alfa),maxval(alfa) found(n) = .true. case('NFLUX') if (.not.allocated(nflux)) allocate(nflux(nlon,nlat)) nflux = f2d if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(nflux),maxval(nflux) found(n) = .true. case('EFLUX') if (.not.allocated(eflux)) allocate(eflux(nlon,nlat)) eflux = f2d if (iprint > 0) & - write(6,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field ',a,' min,max=',2es12.4)") & fnames(n),minval(eflux),maxval(eflux) found(n) = .true. ! ! Unknown field (this should not happen): ! case default - write(6,"('>>> read_tgcm: unknown field: ',a)") trim(fnames(n)) - stop 'read tgcm 3d fields' + write(error_unit,"('>>> read_tgcm: unknown field: ',a)") trim(fnames(n)) + error stop 'read tgcm 3d fields' end select ! fnames(n) else - write(6,"('>>> WARNING: did not find field ',a,' on the tgcm file')") trim(fnames(n)) + write(error_unit,"('>>> WARNING: did not find field ',a,' on the tgcm file')") trim(fnames(n)) endif ! bottom of found variable conditional, scs 11/16 enddo ! n=1,nf ! @@ -380,12 +380,12 @@ subroutine read_tgcm(ncfile,itime) if (found(findx('O1')).and.found(findx('O2')).and.found(findx('HE'))) then n2 = 1.-o2-o1-he if (iprint > 0) & - write(6,"('read_tgcm: Field N2 (1-O2-O-HE) min,max=',2es12.4)") & + write(error_unit,"('read_tgcm: Field N2 (1-O2-O-HE) min,max=',2es12.4)") & minval(n2),maxval(n2) else - write(6,"('>>> FATAL: N2 not found on the file, and at least one of')") - write(6,"(' o2, o, he also not found on the file.')") - stop 'read_tgcm' + write(error_unit,"('>>> FATAL: N2 not found on the file, and at least one of')") + write(error_unit,"(' o2, o, he also not found on the file.')") + error stop 'read_tgcm' endif found(findx('N2'))=.true. endif ! N2 not found @@ -396,15 +396,15 @@ subroutine read_tgcm(ncfile,itime) if (.not.found(findx('TN')).or..not.found(findx('O2')).or.& .not.found(findx('O1')).or..not.found(findx('N2')).or.& .not.found(findx('Z' )).or..not.found(findx('HE'))) then - write(6,"('>>> FATAL: ZG not found on the file, and at least one of')") - write(6,"(' tn,o2,o,n2,he,z also not found on the file.')") - stop 'read_tgcm' + write(error_unit,"('>>> FATAL: ZG not found on the file, and at least one of')") + write(error_unit,"(' tn,o2,o,n2,he,z also not found on the file.')") + error stop 'read_tgcm' endif if (.not.allocated(zg)) allocate(zg(nlon,nlat,nlev)) call calczg(tn,o2,o1,n2,he,z,zg,glat,nlon,nlat,nlev,dlev) zg = zg*1.e-5 ! cm to km if (iprint > 0) & - write(6,"('Field ZG (from calczg) min,max (km)=',2es12.4)") & + write(error_unit,"('Field ZG (from calczg) min,max (km)=',2es12.4)") & minval(zg),maxval(zg) found(findx('ZG'))=.true. endif ! zg not on the file @@ -414,7 +414,7 @@ subroutine read_tgcm(ncfile,itime) if (.not.found(findx('N2D'))) then if (.not.allocated(n2d)) allocate(n2d(nlon,nlat,nlev)) n2d = 0. - write(6,"('read_tgcm: N2D not found on tgcm history, so set it to zero.')") + write(error_unit,"('read_tgcm: N2D not found on tgcm history, so set it to zero.')") found(findx('N2D'))=.true. endif ! @@ -423,8 +423,8 @@ subroutine read_tgcm(ncfile,itime) if (.not.found(findx('HE'))) then if (.not.allocated(he)) allocate(he(nlon,nlat,nlev)) he = 0. - write(6,"('read_tgcm: HE not found on tgcm history, so set it to zero.')") - write(6,"('>>> WARNING: This may cause erroneous O density at high altitude')") + write(error_unit,"('read_tgcm: HE not found on tgcm history, so set it to zero.')") + write(error_unit,"('>>> WARNING: This may cause erroneous O density at high altitude')") found(findx('HE'))=.true. endif ! @@ -433,31 +433,31 @@ subroutine read_tgcm(ncfile,itime) if (.not.found(findx('CUSP'))) then if (.not.allocated(cusp)) allocate(cusp(nlon,nlat)) cusp = 0. - write(6,"('read_tgcm: CUSP not found on tgcm history, so set it to zero.')") + write(error_unit,"('read_tgcm: CUSP not found on tgcm history, so set it to zero.')") found(findx('CUSP'))=.true. endif if (.not.found(findx('DRIZZLE'))) then if (.not.allocated(drizzle)) allocate(drizzle(nlon,nlat)) drizzle = 0. - write(6,"('read_tgcm: DRIZZLE not found on tgcm history, so set it to zero.')") + write(error_unit,"('read_tgcm: DRIZZLE not found on tgcm history, so set it to zero.')") found(findx('DRIZZLE'))=.true. endif if (.not.found(findx('ALFA'))) then if (.not.allocated(alfa)) allocate(alfa(nlon,nlat)) alfa = 0. - write(6,"('read_tgcm: ALFA not found on tgcm history, so set it to zero.')") + write(error_unit,"('read_tgcm: ALFA not found on tgcm history, so set it to zero.')") found(findx('ALFA'))=.true. endif if (.not.found(findx('NFLUX'))) then if (.not.allocated(nflux)) allocate(nflux(nlon,nlat)) nflux = 0. - write(6,"('read_tgcm: NFLUX not found on tgcm history, so set it to zero.')") + write(error_unit,"('read_tgcm: NFLUX not found on tgcm history, so set it to zero.')") found(findx('NFLUX'))=.true. endif if (.not.found(findx('EFLUX'))) then if (.not.allocated(eflux)) allocate(eflux(nlon,nlat)) eflux = 0. - write(6,"('read_tgcm: EFLUX not found on tgcm history, so set it to zero.')") + write(error_unit,"('read_tgcm: EFLUX not found on tgcm history, so set it to zero.')") found(findx('EFLUX'))=.true. endif ! @@ -465,9 +465,9 @@ subroutine read_tgcm(ncfile,itime) ! do n=1,nf if (.not.found(n)) then - write(6,"(/,'>>> FATAL: field ',a,' was not found on tgcm file ',a)") & + write(error_unit,"(/,'>>> FATAL: field ',a,' was not found on tgcm file ',a)") & fnames(n),trim(ncfile) - stop 'read_tgcm' + error stop 'read_tgcm' endif enddo ! @@ -491,13 +491,13 @@ subroutine read_tgcm(ncfile,itime) call denconv(n4s,14.,tn,o2mmr,o1mmr,n2mmr,hemmr,zlev,nlon,nlat,nlev) if (iprint > 0) then - write(6,"('After denconv: O2 (cm3) min,max=',2es12.4)") minval(o2),maxval(o2) - write(6,"('After denconv: O1 (cm3) min,max=',2es12.4)") minval(o1),maxval(o1) - write(6,"('After denconv: N2 (cm3) min,max=',2es12.4)") minval(n2),maxval(n2) - write(6,"('After denconv: HE (cm3) min,max=',2es12.4)") minval(he),maxval(he) - write(6,"('After denconv: NO (cm3) min,max=',2es12.4)") minval(no),maxval(no) - write(6,"('After denconv: N2D(cm3) min,max=',2es12.4)") minval(n2d),maxval(n2d) - write(6,"('After denconv: N4S(cm3) min,max=',2es12.4)") minval(n4s),maxval(n4s) + write(error_unit,"('After denconv: O2 (cm3) min,max=',2es12.4)") minval(o2),maxval(o2) + write(error_unit,"('After denconv: O1 (cm3) min,max=',2es12.4)") minval(o1),maxval(o1) + write(error_unit,"('After denconv: N2 (cm3) min,max=',2es12.4)") minval(n2),maxval(n2) + write(error_unit,"('After denconv: HE (cm3) min,max=',2es12.4)") minval(he),maxval(he) + write(error_unit,"('After denconv: NO (cm3) min,max=',2es12.4)") minval(no),maxval(no) + write(error_unit,"('After denconv: N2D(cm3) min,max=',2es12.4)") minval(n2d),maxval(n2d) + write(error_unit,"('After denconv: N4S(cm3) min,max=',2es12.4)") minval(n4s),maxval(n4s) endif ! ! Close the file: @@ -617,8 +617,8 @@ integer function findx(name) endif enddo if (findx==0) then - write(6,"('>>> findx: could not find index to field ',a)") name - stop 'findx' + write(error_unit,"('>>> findx: could not find index to field ',a)") name + error stop 'findx' endif end function findx @@ -676,19 +676,19 @@ integer function find_mtimes(ncfile,start_mtime,stop_mtime,mtimes,itimes) mtime_file(2,i)==start_mtime(2).and. & mtime_file(3,i)==start_mtime(3)) then itime0 = i - write(6,"('Found requested start_mtime ',3i4,' (history ',i4,' on the file).')") & + write(error_unit,"('Found requested start_mtime ',3i4,' (history ',i4,' on the file).')") & start_mtime,itime0 exit endif enddo endif if (itime0==0) then - write(6,"('>>> Could not find tgcm start_mtime ',3i4,' mtimes on the file are as follows:')") & + write(error_unit,"('>>> Could not find tgcm start_mtime ',3i4,' mtimes on the file are as follows:')") & start_mtime do i=1,ntimes - write(6,"('i=',i4,' mtime(i)=',3i4)") i,mtime_file(:,i) + write(error_unit,"('i=',i4,' mtime(i)=',3i4)") i,mtime_file(:,i) enddo - stop 'start_mtime' + error stop 'start_mtime' endif ! ! Search for stop time: @@ -701,33 +701,33 @@ integer function find_mtimes(ncfile,start_mtime,stop_mtime,mtimes,itimes) mtime_file(2,i)==stop_mtime(2).and. & mtime_file(3,i)==stop_mtime(3)) then itime1 = i - write(6,"('Found requested stop_mtime ',3i4,' (history ',i4,' on the file).')") & + write(error_unit,"('Found requested stop_mtime ',3i4,' (history ',i4,' on the file).')") & stop_mtime,itime1 exit endif enddo endif if (itime1==0) then - write(6,"('>>> Could not find tgcm stop_mtime ',3i4,' mtimes on the file are as follows:')") & + write(error_unit,"('>>> Could not find tgcm stop_mtime ',3i4,' mtimes on the file are as follows:')") & stop_mtime do i=1,ntimes - write(6,"('i=',i4,' mtime(i)=',3i4)") i,mtime_file(:,i) + write(error_unit,"('i=',i4,' mtime(i)=',3i4)") i,mtime_file(:,i) enddo - stop 'start_mtime' + error stop 'start_mtime' endif ! ! This should not happen, but you never know... if (itime1 < itime0) then - write(6,"('>>> find_mtimes: bad itime0 must be <= itime1: itime0,1=',2i4)") & + write(error_unit,"('>>> find_mtimes: bad itime0 must be <= itime1: itime0,1=',2i4)") & itime0,itime1 - write(6,"('>>> mtime_file(itime0)=',3i4,' mtime_file(itime1)=',3i4)") & + write(error_unit,"('>>> mtime_file(itime0)=',3i4,' mtime_file(itime1)=',3i4)") & mtime_file(:,itime0),mtime_file(:,itime1) - stop 'itime0,1' + error stop 'itime0,1' endif ! ! Return values: find_mtimes = itime1-itime0+1 ! number of model times - write(6,"('Number of model times to read = ',i4)") find_mtimes + write(error_unit,"('Number of model times to read = ',i4)") find_mtimes do i=itime0,itime1 mtimes(:,i-itime0+1) = mtime_file(:,i) ! model times from start to stop itimes(i-itime0+1) = i ! file index of each model time @@ -757,14 +757,14 @@ subroutine read_tgcm_coords(ncfile) write(msg,"('Error opening file ',a)") trim(ncfile) call handle_ncerr(istat,trim(msg),1) else - write(6,"(/,'Opened file ',a)") trim(ncfile) + write(error_unit,"(/,'Opened file ',a)") trim(ncfile) endif ! ! Get number of times (histories) (length of unlimited variable): ! istat = nf90_inq_dimid(ncid,'time',idunlim) istat = nf90_inquire_dimension(ncid,idunlim,varname,ntime) -! write(6,"('read_tgcm_coords: ntime=',i4)") ntime +! write(error_unit,"('read_tgcm_coords: ntime=',i4)") ntime istat = nf90_inq_dimid(ncid,'lon',id) istat = nf90_inquire_dimension(ncid,id,varname,nlon) @@ -798,13 +798,13 @@ subroutine read_tgcm_coords(ncfile) istat = nf90_inq_varid(ncid,'time',idv_time) istat = nf90_get_var(ncid,idv_time,time,(/1/),(/ntime/)) -! write(6,"('read_tgcm_coords: nlon =',i4,' glon =',/,(8f10.3))") nlon,glon -! write(6,"('read_tgcm_coords: nlat =',i4,' glat =',/,(8f10.3))") nlat,glat -! write(6,"('read_tgcm_coords: nlev =',i4,' zlev =',/,(8f10.3))") nlev,zlev -! write(6,"('read_tgcm_coords: nilev=',i4,' zilev=',/,(8f10.3))") nilev,zilev -! write(6,"('read_tgcm_coords: nilev=',i4,' zilev(nilev)=',f9.3,' zilev(1)=',f9.3,' dlev=',f9.3)") & +! write(error_unit,"('read_tgcm_coords: nlon =',i4,' glon =',/,(8f10.3))") nlon,glon +! write(error_unit,"('read_tgcm_coords: nlat =',i4,' glat =',/,(8f10.3))") nlat,glat +! write(error_unit,"('read_tgcm_coords: nlev =',i4,' zlev =',/,(8f10.3))") nlev,zlev +! write(error_unit,"('read_tgcm_coords: nilev=',i4,' zilev=',/,(8f10.3))") nilev,zilev +! write(error_unit,"('read_tgcm_coords: nilev=',i4,' zilev(nilev)=',f9.3,' zilev(1)=',f9.3,' dlev=',f9.3)") & ! nilev,zilev(nilev),zilev(1),dlev - write(6,"('read_tgcm_coords: ntime=',i4,' nlon=',i4,' nlat=',i4,' nlev=',i4)") ntime,nlon,nlat,nlev + write(error_unit,"('read_tgcm_coords: ntime=',i4,' nlon=',i4,' nlat=',i4,' nlev=',i4)") ntime,nlon,nlat,nlev ! ! Allocate and read 1d time-dependent variables (or, this could be done ! on a per-history basis below, reading into scalars instead of arrays @@ -852,17 +852,17 @@ subroutine read_tgcm_coords(ncfile) istat = nf90_inq_varid(ncid,'ed',id) istat = nf90_get_var (ncid,id,ed) -! write(6,"('read_tgcm_coords: ntime=',i4,' f107d=',/,(8f9.2))") ntime,f107d -! write(6,"('read_tgcm_coords: ntime=',i4,' f107a=',/,(8f9.2))") ntime,f107a -! write(6,"('read_tgcm_coords: ntime=',i4,' hpower= ',/,(8f9.2))") ntime,hpower -! write(6,"('read_tgcm_coords: ntime=',i4,' e1= ',/,(8f9.2))") ntime,e1 -! write(6,"('read_tgcm_coords: ntime=',i4,' e2= ',/,(8f9.2))") ntime,e2 -! write(6,"('read_tgcm_coords: ntime=',i4,' h1= ',/,(8f9.2))") ntime,h1 -! write(6,"('read_tgcm_coords: ntime=',i4,' h2= ',/,(8f9.2))") ntime,h2 -! write(6,"('read_tgcm_coords: ntime=',i4,' alfac=',/,(8f9.2))") ntime,alfac -! write(6,"('read_tgcm_coords: ntime=',i4,' alfad=',/,(8f9.2))") ntime,alfad -! write(6,"('read_tgcm_coords: ntime=',i4,' ec= ',/,(8f9.2))") ntime,ec -! write(6,"('read_tgcm_coords: ntime=',i4,' ed= ',/,(8f9.2))") ntime,ed +! write(error_unit,"('read_tgcm_coords: ntime=',i4,' f107d=',/,(8f9.2))") ntime,f107d +! write(error_unit,"('read_tgcm_coords: ntime=',i4,' f107a=',/,(8f9.2))") ntime,f107a +! write(error_unit,"('read_tgcm_coords: ntime=',i4,' hpower= ',/,(8f9.2))") ntime,hpower +! write(error_unit,"('read_tgcm_coords: ntime=',i4,' e1= ',/,(8f9.2))") ntime,e1 +! write(error_unit,"('read_tgcm_coords: ntime=',i4,' e2= ',/,(8f9.2))") ntime,e2 +! write(error_unit,"('read_tgcm_coords: ntime=',i4,' h1= ',/,(8f9.2))") ntime,h1 +! write(error_unit,"('read_tgcm_coords: ntime=',i4,' h2= ',/,(8f9.2))") ntime,h2 +! write(error_unit,"('read_tgcm_coords: ntime=',i4,' alfac=',/,(8f9.2))") ntime,alfac +! write(error_unit,"('read_tgcm_coords: ntime=',i4,' alfad=',/,(8f9.2))") ntime,alfad +! write(error_unit,"('read_tgcm_coords: ntime=',i4,' ec= ',/,(8f9.2))") ntime,ec +! write(error_unit,"('read_tgcm_coords: ntime=',i4,' ed= ',/,(8f9.2))") ntime,ed istat = nf90_close(ncid) @@ -877,13 +877,13 @@ subroutine handle_ncerr(istat,msg,ifatal) integer,intent(in) :: istat,ifatal character(len=*),intent(in) :: msg ! - write(6,"(/72('-'))") - write(6,"('>>> Error from netcdf library:')") - write(6,"(a)") trim(msg) - write(6,"('istat=',i5)") istat - write(6,"(a)") nf90_strerror(istat) - write(6,"(72('-')/)") - if (ifatal > 0) stop('Fatal netcdf error') + write(error_unit,"(/72('-'))") + write(error_unit,"('>>> Error from netcdf library:')") + write(error_unit,"(a)") trim(msg) + write(error_unit,"('istat=',i5)") istat + write(error_unit,"(a)") nf90_strerror(istat) + write(error_unit,"(72('-')/)") + if (ifatal > 0) error stop 'Fatal netcdf error' end subroutine handle_ncerr !----------------------------------------------------------------------- diff --git a/rout.f90 b/rout.f90 index 72c9921..5756741 100644 --- a/rout.f90 +++ b/rout.f90 @@ -46,7 +46,9 @@ SUBROUTINE ROUT(ROFILE,LUN,EF,EZ,ITAIL,FRACO,FRACO2,FRACN2) write(lun,"(4f8.2)") f107,f107p,f107a,xuvfac write(lun,"(' Eflux ',' Ezero ',' Itail ',' FracO ',' FracO2 ',' FracN2 ')") write(lun,"(f8.2,f8.1,i8,3f8.2)") ef, ez, itail, fraco, fraco2, fracn2 - write(lun,"(' Alt Tn Ti Te O O2 N2 He N Ne O+ 1356 1304 1027 989 LBH')") + write(lun,"(' Alt Tn Ti Te O O2 N2 & + He N Ne O+ 1356 1304 1027 989& + LBH')") do j=1,jmax write(lun,"(0p,f6.1,3f6.0,1p,12e9.2)") & diff --git a/setup.py b/setup.py new file mode 100644 index 0000000..7c1b24f --- /dev/null +++ b/setup.py @@ -0,0 +1,36 @@ +#!/usr/bin/env python +req = ['nose','numpy','matplotlib','seaborn'] +pipreq=['sciencedates'] +# %% +import pip +try: + import conda.cli + conda.cli.main('install',*req) +except Exception: + pip.main(['install'] + req) +pip.main(['install'] + pipreq) +# %% +import setuptools # enables develop +from numpy.distutils.core import Extension, setup + +src = 'cglow.f90 glow.f90 bands.f90 conduct.f90 egrid.f90 ephoto.f90 etrans.f90 exsect.f fieldm.f gchem.f90 geomag.f90 maxt.f90 mzgrid.f90 qback.f90 rcolum.f90 rout.f90 snoem.f90 snoemint.f90 solzen.f90 ssflux.f90 iri90.f nrlmsise00.f'.split(' ') + +ext = Extension(name='glow', sources=src, + f2py_options=['--quiet']) + + +setup(name='glowiono', + packages=['glowiono'], + version='1.0', + author=['Stan Solomon','Liam Kilcommons','Michael Hirsch, Ph.D.'], + url = 'https://github.com/NCAR/GLOW', + description='GLobal airglOW model', + classifiers=[ + 'Intended Audience :: Science/Research', + 'Development Status :: 5 - Production/Stable', + 'Topic :: Scientific/Engineering :: Atmospheric Science', + 'Programming Language :: Python :: 3', + ], + ext_modules=[ ext ], + data_files=[], + ) diff --git a/snoem.f90 b/snoem.f90 index 2c85e01..86bd2e9 100644 --- a/snoem.f90 +++ b/snoem.f90 @@ -11,7 +11,7 @@ subroutine snoem(doy, kp, f107, z, mlat, nozm) - use cglow,only: data_dir + use cglow,only: data_dir,sind,cosd implicit none save @@ -34,7 +34,7 @@ subroutine snoem(doy, kp, f107, z, mlat, nozm) if (ifirst == 1) then ifirst = 0 filepath = trim(data_dir)//'snoem_eof.dat' - open(unit=1,file=filepath,status='old',readonly) + open(unit=1,file=filepath,status='old', action='read') read(1,*) (z(k),k=1,16) read(1,*) (mlat(j),j=1,33) read(1,*) ((no_mean(j,k),j=1,33),k=1,16) diff --git a/ssflux.f90 b/ssflux.f90 index 1e17b18..1d66c21 100644 --- a/ssflux.f90 +++ b/ssflux.f90 @@ -137,7 +137,7 @@ subroutine ssflux (iscale,f107,f107a,xuvfac,wave1,wave2,sflux) if (iscale == 0) then if (islast /= iscale) then filepath = trim(data_dir)//'ssflux_hint.dat' - open(unit=1,file=filepath,status='old',readonly) + open(unit=1,file=filepath,status='old',action='read') read(1,*) do l=lmax,1,-1 read(1,*) waves(l),wavel(l),rflux(l),scale1(l),scale2(l) @@ -161,7 +161,7 @@ subroutine ssflux (iscale,f107,f107a,xuvfac,wave1,wave2,sflux) if (iscale == 1) then if (islast /= iscale) then filepath = trim(data_dir)//'ssflux_euvac.dat' - open(unit=1,file=filepath,status='old',readonly) + open(unit=1,file=filepath,status='old',action='read') read(1,*) do l=lmax,1,-1 read(1,*) waves(l),wavel(l),rflux(l),a(l) @@ -184,7 +184,7 @@ subroutine ssflux (iscale,f107,f107a,xuvfac,wave1,wave2,sflux) if (iscale == 2) then if (islast /= iscale) then filepath = trim(data_dir)//'ssflux_user.dat' - open(unit=1,file=filepath,status='old',readonly) + open(unit=1,file=filepath,status='old',action='read') read(1,*) do l=lmax,1,-1 read(1,*) waves(l),wavel(l),uflux(l) diff --git a/tests/test.py b/tests/test.py new file mode 100755 index 0000000..88a82f6 --- /dev/null +++ b/tests/test.py @@ -0,0 +1,146 @@ +#!/usr/bin/env python +""" +Registration testing of GLOW +Michael Hirsch + +""" +import logging +from datetime import datetime +from itertools import chain +from numpy import array,zeros,float32,log,isclose,nan +from numpy.testing import assert_allclose,run_module_suite +# +#from glowiono import makeeigen +from sciencedates import datetime2yd,datetime2gtd +# +import glow +msise00=None +#%% test inputs +z = list(range(30,110+1,1)) +z += ( + [111.5,113.,114.5,116.] + + list(chain(range(118,150+2,2),range(153,168+3,3),range(172,180+4,4), + range(185,205+5,5),range(211,223+6,6),range(230,244+7,7), + range(252,300+8,8),range(309,345+9,9),range(355,395+10,10), + range(406,428+11,11))) + + [440,453,467,482,498,515,533,551] + + list(range(570,950+20,20)) + ) +z = array(z) + +nbins = 190; jmax=170 #glow.h + +eflux = 1. +e0 = 1e3 +maxind = 112 +glat = 70; glon=0 #like aurora.in +ap=4; f107=100; f107a=100 +nmaj=3; nst=6 +dtime = datetime(1999,12,21) +# +yd,utsec = datetime2yd(dtime)[:2] + +def test_egrid_maxt(): + ener,dE = glow.egrid(199) + assert_allclose(ener[[maxind,maxind+10,-1]],[1017.7124,1677.9241,75005.171875],rtol=1e-5) +#%% test of maxt + phi = glow.maxt(eflux,e0,ener, dE, itail=0, fmono=nan, emono=nan) + assert phi.argmax() == maxind + assert_allclose(phi[[maxind,maxind+10]],[ 114810.6,97814.438]) + +#%% test vquart (quartic root) KNOWN DEFECTIVE FORTRAN ALGORITHM +#Aquart = tile([-1,0,0,0,1],(jmax,1)) +#qroot = glow.vquart(Aquart,1) +#assert_allclose(qroot[0],roots(Aquart[0,-1])) +#Aquart = array([[-1,0,0,0,1], +# [-1,0,0,1,1]]) +#nq = Aquart.shape[0] +#Aquart = tile(Aquart,(jmax//nq,1)) +#qroot = glow.vquartmod.vquart(Aquart, nq) +#try: +# assert_allclose(qroot[:nq], +# [1,0.8191725133961643]) +#except AssertionError as e: +# print('this mismatch is in discussion with S. Solomon. {}'.format(e)) + +def test_solzen(): + sza = glow.solzen(yd,utsec,glat,glon) + assert isclose(sza, 133.43113708496094) + return sza + +#def test_snoem(): +# doy = datetime2gtd(dtime)[0] +# zno,maglat,nozm = glow.snoem(doy,1.75*log(0.4*ap),f107) + # assert_allclose((nozm[12,15],nozm[-2,-1]),(35077728.0, 1.118755e+08)) +# return nozm + +def test_snoemint(): + if msise00 is not None: + densd,tempd = rungtd1d(dtime,z,glat,glon,f107a,f107,[ap]*7) + # (nighttime background ionization) + print(tempd) + znoint = glow.snoemint(dtime.strftime('%Y%j'),glat,glon,f107,ap,z, + tempd.loc[:,'Tn']) + assert_allclose(znoint[[28,143]], (1.262170e+08, 3.029169e+01),rtol=1e-5) #arbitrary + return znoint + +#def test_fieldm(): +# xdip,ydip,zdip,totfield,dipang,decl,smodip = glow.fieldm(glat,glon%360,z[50]) +# assert isclose(xdip,0.1049523800611496) + # assert isclose(totfield,0.5043528079986572) +# assert isclose(dipang,77.72911071777344) + +#def test_ssflux(): +# iscale=1; hlybr=0.; hlya=0.; fexvir=0.; heiew=0.; xuvfac=3. +# wave1,wave2,sflux = glow.ssflux(iscale,f107,f107a,hlybr,fexvir,hlya,heiew,xuvfac) +# assert_allclose(sflux[[11,23]],(4.27225743e+11, 5.54400400e+07)) + +def test_rcolum_qback(): + if msise00 is not None: + densd,tempd = rungtd1d(dtime,z,glat,glon,f107a,f107,[ap]*7) + + """ VCD: Vertical Column Density """ + sza = test_solzen() + zcol,zvcd = glow.rcolum(sza,z*1e5, + densd.loc[:,['O','O2','N2']].values.T, + tempd.loc[:,'Tn']) + # FIXME these tests were numerically unstable (near infinity values) + assert isclose(zcol[0,0], 1e30) #see rcolum comments for sun below horizon 1e30 + assert isclose(zvcd[2,5],5.97157e+28,rtol=1e-2) #TODO changes a bit between python 2 / 3 + #%% skipping EPHOTO since we care about night time more for now + znoint = test_snoemint() + # zeros because nighttime + photoi = zeros((nst,nmaj,jmax),dtype=float32,order='F') + phono = zeros((nst,jmax),dtype=float32,order='F') + glow.qback(zmaj=densd.loc[:,['O','O2','N2']].values.T, + zno=znoint, + zvcd=zvcd, + photoi=photoi,phono=phono) + #arbitrary point check + assert isclose(photoi[0,0,77],1.38091e-18,rtol=1e-5) + assert isclose(phono[0,73],0.0,rtol=1e-5) + else: + logging.warning('skipped rcolum qback due to missing external msise00') + +#def test_glow(): + # electron precipitation + #First enact "glow" subroutine, which calls QBACK, ETRANS and GCHEM among others + +# glow.glow() #no args + + #%% ver and constituants +# """ +# currently using common block CGLOW, in future use module +# """ + # zceta = glow.cglow.zceta.T + # zeta = glow.cglow.zeta.T[:,:11] +# zcsum = zceta.sum(axis=-1)[:,:11] +# assert_allclose(zcsum,zeta,rtol=1e-6) + +#def test_eigen(): +# ener,dE = glow.egrid() +# ver,photIon,isr,phitop,zceta,sza,prates,lrates,tezs,sion=makeeigen(ener,ones_like(ener),dtime,(glat,glon)) + +if __name__ == '__main__': + #test_snoemint() + run_module_suite() diff --git a/tzgrid.f90 b/tzgrid.f90 index e183874..cc68869 100644 --- a/tzgrid.f90 +++ b/tzgrid.f90 @@ -35,6 +35,7 @@ subroutine tzgrid(i,l,jmax,z,zo,zo2,zn2,zns,znd,zno,ztn,zun,zvn,ze,zti,zte) + use, intrinsic:: iso_fortran_env,only: error_unit use readtgcm,only: nlev,zg,tn,un,vn,o2,o1,n2,n4s,n2d,no,ti,te,ne implicit none @@ -45,8 +46,8 @@ subroutine tzgrid(i,l,jmax,z,zo,zo2,zn2,zns,znd,zno,ztn,zun,zvn,ze,zti,zte) integer :: j,jj if (nlev/=29 .and. nlev/=57 .and. nlev/=49 .and. nlev/=97) then - write(6,"('zgrid: unknown NLEV = ',i5)") nlev - stop 'zgrid' + write(error_unit,"('zgrid: unknown NLEV = ',i5)") nlev + error stop 'zgrid' endif if (nlev == 29) then ! low-res TIE-GCM @@ -171,7 +172,7 @@ subroutine tzgrid(i,l,jmax,z,zo,zo2,zn2,zns,znd,zno,ztn,zun,zvn,ze,zti,zte) endif enddo z(jmax) = z(jmax-1) + (z(jmax-1)-z(jmax-2)) !patch top level - ze(jmax) = ze(jmax-1)**2 / + ze(jmax-2) + ze(jmax) = ze(jmax-1)**2 / (+ ze(jmax-2)) endif if (nlev == 97) then ! high-res TIME-GCM