diff --git a/.gitignore b/.gitignore index fa5040d..8ff29dd 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,5 @@ +baseIRI2016 +testiri2016 bin/ # Byte-compiled / optimized / DLL files diff --git a/.travis.yml b/.travis.yml index 657506e..6433790 100644 --- a/.travis.yml +++ b/.travis.yml @@ -21,8 +21,11 @@ before_install: # temp for pip < 9.1 - pip -q install numpy -install: pip install -e .[tests] +install: pip -q install -e .[tests] -script: coverage run tests/test_all.py -v +script: + - pytest -v -after_success: coveralls +after_success: + - coverage run tests/test_all.py + - coveralls diff --git a/AltitudeProfile.py b/AltitudeProfile.py new file mode 100755 index 0000000..7e976d6 --- /dev/null +++ b/AltitudeProfile.py @@ -0,0 +1,38 @@ +#!/usr/bin/env python +""" Height Profile Example """ +import pyiri2016 +# +import numpy as np +from matplotlib.pyplot import figure, show + +glat, glon = -11.95, -76.77 + +alt_km = np.arange(80,1000,20.) + +iri = pyiri2016.IRI('2012-08-21T12', alt_km, glat, glon) + +fig = figure(figsize=(16,6)) +axs = fig.subplots(1,2) + +fig.suptitle(f'{str(iri.time[0].values)[:-13]}\n Glat, Glon: {iri.glat}, {iri.attrs["glon"]}') + + +pn = axs[0] +pn.plot(iri['ne'].squeeze(), iri.alt_km, label='N$_e$') +#pn.set_title(iri2016Obj.title1) +pn.set_xlabel('Density (m$^{-3}$)') +pn.set_ylabel('Altitude (km)') +pn.set_xscale('log') +pn.legend(loc='best') +pn.grid(True) + +pn = axs[1] +pn.plot(iri['Ti'].squeeze(), iri.alt_km, label='T$_i$') +pn.plot(iri['Te'].squeeze(), iri.alt_km, label='T$_e$') +#pn.set_title(iri2016Obj.title2) +pn.set_xlabel('Temperature (K)') +pn.set_ylabel('Altitude (km)') +pn.legend(loc='best') +pn.grid(True) + +show() \ No newline at end of file diff --git a/LatitudeProfile.py b/LatitudeProfile.py new file mode 100755 index 0000000..3ef3fe8 --- /dev/null +++ b/LatitudeProfile.py @@ -0,0 +1,45 @@ +#!/usr/bin/env python +import pyiri2016 +from numpy import arange +from matplotlib.pyplot import figure, show + +""" Geog. Latitude Profile Example """ + +latlim = [-60, 60] +latstp = 2. +sim = pyiri2016.geoprofile(altkm=600, latlim=latlim, dlat=latstp, \ +glon=-76.77, time='2004-01-01T17') + +latbins = arange(latlim[0], latlim[1], latstp) + + +fig = figure(figsize=(8,12)) +axs = fig.subplots(2,1, sharex=True) + +pn = axs[0] + +pn.plot(latbins, sim['NmF2'].squeeze(), label='N$_m$F$_2$') +pn.plot(latbins, sim['NmF1'].squeeze(), label='N$_m$F$_1$') +pn.plot(latbins, sim['NmE'].squeeze(), label='N$_m$E') +pn.set_title(str(sim.time[0].values)[:-13] + ' latitude'+str(latlim)) +pn.set_xlim(latbins[[0, -1]]) +pn.set_xlabel('Geog. Lat. ($^\circ$)') +pn.set_ylabel('(m$^{-3}$)') +pn.set_yscale('log') + + +pn = axs[1] +pn.plot(latbins, sim['hmF2'].squeeze(), label='h$_m$F$_2$') +pn.plot(latbins, sim['hmF1'].squeeze(), label='h$_m$F$_1$') +pn.plot(latbins, sim['hmE'].squeeze(), label='h$_m$E') +pn.set_xlim(latbins[[0, -1]]) +pn.set_title(str(sim.time[0].values)[:-13] + ' latitude'+str(latlim)) +pn.set_xlabel('Geog. Lat. ($^\circ$)') +pn.set_ylabel('(km)') + +for a in axs: + a.legend(loc='best') + a.grid(True) + + +show() diff --git a/README.rst b/README.rst index 80ba998..4c6106b 100644 --- a/README.rst +++ b/README.rst @@ -31,20 +31,20 @@ Usage Height-profile --------------- -`plot density and temperatures vs height `_ +`plot density and temperatures vs height `_ .. image:: figures/iri1DExample01.png Latitudinal profile ------------------- -`plot densities and height at the peak of F2, F2, and E regions vs geographic latitude `_ +`plot densities and height at the peak of F2, F2, and E regions vs geographic latitude `_ .. image:: figures/iri1DExample02.png GMT profile ----------- -`plot densities and height at the peak of F2, F2, and E regions vs universal time `_ +`plot densities and height at the peak of F2, F2, and E regions vs universal time `_ .. image:: figures/iri1DExample08.png @@ -74,7 +74,7 @@ Fortran compile make - ./testiri2016 + make test f2py compile diff --git a/TimeProfile.py b/TimeProfile.py new file mode 100755 index 0000000..26729a9 --- /dev/null +++ b/TimeProfile.py @@ -0,0 +1,104 @@ +#!/usr/bin/env python +""" Time Profile: IRI2016 """ +import numpy as np +from datetime import timedelta +from matplotlib.pyplot import figure, show +try: + import ephem +except ImportError: + ephem = None +# +import pyiri2016 + + +# %% user parameters +#lat = -11.95; lon = -76.77 +#glat, glon = 0,0 +glat,glon = 65,-148 +alt_km = np.arange(120, 180, 20) +# %% ru +sim = pyiri2016.timeprofile(('2012-08-21','2012-08-22'),timedelta(hours=0.25), + alt_km,glat,glon) + +# %% Plots +Nplot=3 + +if Nplot>2: + fig = figure(figsize=(16,12)) + axs = fig.subplots(3,1, sharex=True).ravel() +else: + fig = figure(figsize=(16,6)) + axs = fig.subplots(1,2).ravel() + +fig.suptitle(f'{str(sim.time[0].values)[:-13]} to {str(sim.time[-1].values)[:-13]}\n Glat, Glon: {sim.glat}, {sim.glon}') + +ax = axs[0] +#NmF1 = pyiri2016.IRI2016()._RmNeg(sim.b[2, :]) +#NmE = sim.loc[4, :] +ax.plot(sim.time, sim['NmF2'].squeeze(), label='N$_m$F$_2$') +#ax.plot(sim.time, NmF1, label='N$_m$F$_1$') +#ax.plot(sim.time, NmE, label='N$_m$E') +ax.set_title('Maximum number densities vs. ionospheric layer') +ax.set_xlabel('Hour (UT)') +ax.set_ylabel('(m$^{-3}$)') +ax.set_yscale('log') +ax.legend(loc='best') + +ax = axs[1] +#hmF1 = pyiri2016.IRI2016()._RmNeg(sim.b[3, :]) +#hmE = sim.b[5, :] +ax.plot(sim.time, sim['hmF2'].squeeze(), label='h$_m$F$_2$') +#ax.plot(sim.time, hmF1, label='h$_m$F$_1$') +#ax.plot(sim.time, hmE, label='h$_m$E') +ax.set_title('Height of maximum density vs. ionospheric layer') +ax.set_xlabel('Hour (UT)') +ax.set_ylabel('(km)') +ax.legend(loc='best') +# %% +if Nplot > 2: + ax = axs[2] + + for a in sim.alt_km: + ax.plot(sim.time, sim['ne'].squeeze(), marker='.', label=f'{a.item()} km') + ax.set_xlabel('time UTC (hours)') + ax.set_ylabel('[m$^{-3}$]') + ax.set_title(f'$N_e$ vs. altitude and time') + ax.set_yscale('log') + ax.legend(loc='best') +# %% +if Nplot > 4: + ax = axs[4] + tec = sim.b[36, :] + ax.plot(sim.time, tec, label=r'TEC') + ax.set_xlabel('Hour (UT)') + ax.set_ylabel('(m$^{-2}$)') + #ax.set_yscale('log') + ax.legend(loc='best') + + ax = axs[5] + vy = sim.b[43, :] + ax.plot(sim.time, vy, label=r'V$_y$') + ax.set_xlabel('Hour (UT)') + ax.set_ylabel('(m/s)') + ax.legend(loc='best') + +for a in axs.ravel(): + a.grid(True) + + +if __name__ == '__main__': + from argparse import ArgumentParser + p = ArgumentParser(description='IRI2016 time profile plot') + p.add_argument('-t','--trange',help='START STOP STEP (hours) time [UTC]',nargs=3, + default=('2012-08-21','2012-08-22',0.25)) + p.add_argument('--alt',help='START STOP STEP altitude [km]',type=float, nargs=3,default=(120,180,20)) + p.add_argument('-c','--latlon',help='geodetic coordinates of simulation', + type=float,default=(65,-147.5)) + p.add_argument('--f107',type=float,default=200.) + p.add_argument('--f107a', type=float,default=200.) + p.add_argument('--ap', type=int, default=4) + p.add_argument('--species',help='species to plot',nargs='+',default=('ne')) + p = p.parse_args() + + +show() diff --git a/data/index/ig_rz.dat b/data/index/ig_rz.dat index b436d71..da91558 100644 --- a/data/index/ig_rz.dat +++ b/data/index/ig_rz.dat @@ -1,4 +1,4 @@ -10,17,2016, +11,2,2017, 1,1958,12,2018, @@ -61,10 +61,10 @@ 75.5, 75.8, 75.1, 74.5, 75.4, 78.0, 79.4, 81.6, 86.6, 90.8, 92.1, 91.8, 92.1, 92.5, 93.2, 94.8, 95.6, 96.0, 96.8, 96.8, 94.4, 92.0, 91.4, 91.2, 90.1, 88.1, 85.0, 81.0, 76.4, 71.3, 66.5, 62.0, 58.5, 55.5, 53.0, 50.1, - 41.0, 38.5, 37.0, 35.7, 35.3, 35.3, 34.1, 33.5, 33.2, 33.9, 34.4, 33.0, - 35.6, 30.3, 31.0, 29.6, 30.7, 28.9, 27.0, 23.6, 22.2, 15.7, 17.4, 11.8, - 9.0, 8.3, 6.9, 3.0, 4.1, 8.8, 6.2, 8.5, 11.5, 13.5, 15.2, 14.3, - 18.3, + 41.0, 38.5, 37.0, 35.2, 32.3, 28.6, 24.7, 21.2, 18.2, 15.9, 14.2, 12.6, + 11.0, 9.2, 8.0, 8.3, 8.7, 9.3, 9.8, 8.6, 9.8, 9.7, 9.7, 9.5, + 8.1, 9.5, 8.0, 9.6, 8.8, 10.0, 7.7, 5.9, 6.6, 5.7, 7.2, 5.2, + 2.5, 200.1, 199.0,200.9,201.3,196.8,191.4,186.8,185.2,184.9,183.8,182.2,180.7,180.5, @@ -124,9 +124,10 @@ 65.5, 66.9, 66.8, 64.6, 61.7, 58.9, 57.8, 58.2, 58.1, 58.6, 59.7, 59.6, 58.7, 58.4, 57.5, 57.8, 59.8, 62.6, 65.4, 68.9, 73.1, 75.0, 75.4, 76.0, 109.3,110.5,114.3,116.4,115.0,114.1,112.6,108.3,101.9, 97.3, 94.7, 92.2, - 89.3, 86.1, 82.2, 78.9, 76.1, 72.1, 68.3, 66.4, 65.9, 64.3, 61.3 57.8, - 54.5, 52.6, 50.6, 48.2, 46.3, 44.5, 42.7, 40.9, 39.2, 38.2, 37.1, 36.7, - 37.0, 35.9, 34.3, 33.2, 32.0, 30.7, 29.4, 28.3, 27.4, 26.7, 26.1, 25.3, - 24.5, 23.6, 22.6, 21.8, 21.2, 20.3, 19.5, 18.6, 17.4, 16.5, 15.9, 15.5, - 15.3, + 89.3, 86.1, 82.2, 78.9, 76.1, 72.1, 68.3, 66.4, 65.9, 64.3, 61.2, 57.8, + 54.4, 52.5, 50.4, 47.8, 44.8, 41.5, 38.5, 36.0, 33.2, 31.5, 29.9, 28.5, + 27.8, 26.5, 25.8, 24.9, 24.1, 24.1, 23.9, 23.4, 23.2, 22.6, 21.9, 21.7, + 21.5, 20.6, 18.6, 17.4, 17.0, 16.1, 15.4, 14.8, 14.3, 14.0, 13.7, 13.2, + 12.7 + diff --git a/examples/TimeProfile.py b/examples/TimeProfile.py deleted file mode 100755 index 50f2eba..0000000 --- a/examples/TimeProfile.py +++ /dev/null @@ -1,94 +0,0 @@ -#!/usr/bin/env python -from pyiri2016 import IRI2016,IRI2016Profile -from numpy import arange -from matplotlib.pyplot import figure, show - -""" Time Profile Example """ - -hrlim = [0, 24] # [0,24] does whole day(s) -hrstp = 0.1 # time step [decimal hours] -#lat = -11.95; lon = -76.77 -lat = 65; lon = -147.5 -alts = arange(100, 200, 10) -sim = IRI2016Profile(hrlim=hrlim, hrstp=hrstp, lat=lat, lon=lon, alt=150., - option='time', verbose=False, time='2017-08-21') - -hrbins = arange(hrlim[0], hrlim[1] + hrstp, hrstp) - -nhr = hrbins.size -index = range(nhr) - -Nplot=4 - - -if Nplot>2: - fig = figure(figsize=(16,12)) - axs = fig.subplots(2,2, sharex=True).ravel() -else: - fig = figure(figsize=(16,6)) - axs = fig.subplots(1,2).ravel() - -pn = axs[0] -NmF2 = sim.b[0, index] -NmF1 = IRI2016()._RmNeg(sim.b[2, index]) -NmE = sim.b[4, index] -pn.plot(hrbins, NmF2, label='N$_m$F$_2$') -pn.plot(hrbins, NmF1, label='N$_m$F$_1$') -pn.plot(hrbins, NmE, label='N$_m$E') -pn.set_title(sim.title1) -pn.set_xlim(hrbins[[0, -1]]) -pn.set_xlabel('Hour (UT)') -pn.set_ylabel('(m$^{-3}$)') -pn.set_yscale('log') -pn.legend(loc='best') - -pn = axs[1] -hmF2 = sim.b[1, index] -hmF1 = IRI2016()._RmNeg(sim.b[3, index]) -hmE = sim.b[5, index] -pn.plot(hrbins, hmF2, label='h$_m$F$_2$') -pn.plot(hrbins, hmF1, label='h$_m$F$_1$') -pn.plot(hrbins, hmE, label='h$_m$E') -pn.set_xlim(hrbins[[0, -1]]) -pn.set_title(sim.title2) -pn.set_xlabel('Hour (UT)') -pn.set_ylabel('(km)') -pn.legend(loc='best') - -if Nplot > 2: - pn = axs[2] - - for alt in alts: - sim = IRI2016Profile(hrlim=hrlim, hrstp=hrstp, lat=lat, lon=lon, alt=alt, - option='time', verbose=False, time='2017-08-10') - - pn.plot(sim.out.time, sim.out.loc[:,'ne'], marker='.', label=alt) - pn.set_xlabel('time UTC (hours)') - pn.set_ylabel('[m$^{-3}$]') - pn.set_title(f'$N_e$ vs. altitude and time') - pn.legend(loc='best') - -if Nplot > 4: - pn = axs[4] - tec = sim.b[36, index] - pn.plot(hrbins, tec, label=r'TEC') - pn.set_xlim(hrbins[[0, -1]]) - pn.set_xlabel('Hour (UT)') - pn.set_ylabel('(m$^{-2}$)') - #pn.set_yscale('log') - pn.legend(loc='best') - - pn = axs[5] - vy = sim.b[43, index] - pn.plot(hrbins, vy, label=r'V$_y$') - pn.set_xlim(hrbins[[0, -1]]) - pn.set_xlabel('Hour (UT)') - pn.set_ylabel('(m/s)') - pn.legend(loc='best') - -for a in axs.ravel(): - a.grid(True) - - - -show() diff --git a/examples/example01.py b/examples/example01.py index caff867..630eb35 100755 --- a/examples/example01.py +++ b/examples/example01.py @@ -1,9 +1,8 @@ #!/usr/bin/env python -from pyiri2016 import IRI2016 +from pyiri2016 import IRI -sim = IRI2016() -IRIData, IRIDATAAdd = sim.IRI() +IRIData, IRIDATAAdd = IRI() print('Ne {:.3e}'.format(IRIData['ne'])) print('NmF2 {:.3e}'.format(IRIDATAAdd['NmF2'])) print('hmF2 {:.3e}'.format(IRIDATAAdd['hmF2'])) diff --git a/examples/iri1DExample02.py b/examples/iri1DExample02.py deleted file mode 100755 index 5843c7e..0000000 --- a/examples/iri1DExample02.py +++ /dev/null @@ -1,52 +0,0 @@ -#!/usr/bin/env python -from pyiri2016 import IRI2016,IRI2016Profile -from numpy import arange -from matplotlib.pyplot import figure, show - -""" Geog. Latitude Profile Example """ - -latlim = [-60, 60] -latstp = 2. -iri2016Obj = IRI2016Profile(alt=600, latlim=latlim, latstp=latstp, \ -lon=-76.77, option='lat', time='2004-01-01T17') - -latbins = arange(latlim[0], latlim[1] + latstp, latstp) - - -index = range(latbins.size) - -fig = figure(figsize=(8,12)) -axs = fig.subplots(2,1, sharex=True) - -pn = axs[0] -NmF2 = iri2016Obj.b[0, index] -NmF1 = IRI2016()._RmNeg(iri2016Obj.b[2, index]) -NmE = iri2016Obj.b[4, index] -pn.plot(latbins, NmF2, label='N$_m$F$_2$') -pn.plot(latbins, NmF1, label='N$_m$F$_1$') -pn.plot(latbins, NmE, label='N$_m$E') -pn.set_title(iri2016Obj.title1) -pn.set_xlim(latbins[[0, -1]]) -pn.set_xlabel('Geog. Lat. ($^\circ$)') -pn.set_ylabel('(m$^{-3}$)') -pn.set_yscale('log') - - -pn = axs[1] -hmF2 = iri2016Obj.b[1, index] -hmF1 = IRI2016()._RmNeg(iri2016Obj.b[3, index]) -hmE = iri2016Obj.b[5, index] -pn.plot(latbins, hmF2, label='h$_m$F$_2$') -pn.plot(latbins, hmF1, label='h$_m$F$_1$') -pn.plot(latbins, hmE, label='h$_m$E') -pn.set_xlim(latbins[[0, -1]]) -pn.set_title(iri2016Obj.title2) -pn.set_xlabel('Geog. Lat. ($^\circ$)') -pn.set_ylabel('(km)') - -for a in axs: - a.legend(loc='best') - a.grid(True) - - -show() diff --git a/examples/iri_heightprofile.py b/examples/iri_heightprofile.py deleted file mode 100755 index 2ae0a53..0000000 --- a/examples/iri_heightprofile.py +++ /dev/null @@ -1,44 +0,0 @@ -#!/usr/bin/env python -""" Height Profile Example """ -from pyiri2016 import IRI2016Profile -# -from numpy import arange -from matplotlib.pyplot import figure, show - - -altlim = [100., 1000.] -altstp = 5. -lat, lon = -11.95, -76.77 - -iri2016Obj = IRI2016Profile(altlim=altlim, altstp=altstp, lat=lat, \ - lon=lon, time='2003-11-21', option='vertical', verbose=False) - -altbins = arange(altlim[0], altlim[1] + altstp, altstp) - -index = range(altbins.size) - -fig = figure(figsize=(16,6)) -axs = fig.subplots(1,2) - -pn = axs[0] -ne = iri2016Obj.a[0, index] -pn.plot(ne, altbins, label='N$_e$') -pn.set_title(iri2016Obj.title1) -pn.set_xlabel('Density (m$^{-3}$)') -pn.set_ylabel('Altitude (km)') -pn.set_xscale('log') -pn.legend(loc='best') -pn.grid(True) - -pn = axs[1] -ti = iri2016Obj.a[2, index] -te = iri2016Obj.a[3, index] -pn.plot(ti, altbins, label='T$_i$') -pn.plot(te, altbins, label='T$_e$') -pn.set_title(iri2016Obj.title2) -pn.set_xlabel('Temperature (K)') -pn.set_ylabel('Altitude (km)') -pn.legend(loc='best') -pn.grid(True) - -show() \ No newline at end of file diff --git a/scripts/iri2DExample01.py b/examples/scripts/iri2DExample01.py similarity index 100% rename from scripts/iri2DExample01.py rename to examples/scripts/iri2DExample01.py diff --git a/scripts/iri2DExample02.py b/examples/scripts/iri2DExample02.py similarity index 100% rename from scripts/iri2DExample02.py rename to examples/scripts/iri2DExample02.py diff --git a/fortran/CMakeLists.txt b/fortran/CMakeLists.txt index ca92925..cf4aff8 100644 --- a/fortran/CMakeLists.txt +++ b/fortran/CMakeLists.txt @@ -1,6 +1,11 @@ cmake_minimum_required(VERSION 2.8.12) project( iri2016 Fortran ) -add_compile_options(-mtune=native) +set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ..) -add_executable(testiri2016 iriwebg.for irisub.for irifun.for iritec.for iridreg.for igrf.for cira.for iriflip.for iritest.for) +add_compile_options(-mtune=native -g -fbacktrace) + +# not used for now +# add_executable(testiri2016 iriwebg.for irisub.for irifun.for iritec.for iridreg.for igrf.for cira.for iriflip.for iritest.for) + +add_executable(baseIRI2016 irisub.for irifun.for iritec.for iridreg.for igrf.for cira.for iriflip.for test.f90) diff --git a/fortran/igrf.for b/fortran/igrf.for index e903e12..3fa58b8 100644 --- a/fortran/igrf.for +++ b/fortran/igrf.for @@ -1,28 +1,28 @@ c igrf.for, version number can be found at the end of this comment. -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- C -C Subroutines to compute IGRF parameters for IRI and all functions and +C Subroutines to compute IGRF parameters for IRI and all functions and C subroutines required for this computation, including: -C IGRF_SUB, IGRF_DIP, FINDB0, SHELLG, STOER, FELDG, FELDCOF, GETSHC, +C IGRF_SUB, IGRF_DIP, FINDB0, SHELLG, STOER, FELDG, FELDCOF, GETSHC, C INTERSHC, EXTRASHC, GEODIP, fmodip C -C CGM coordinates : GEOCGM01, OVL_ANG, CGMGLA, CGMGLO, DFR1DR, -C AZM_ANG, MLTUT, MFC, FTPRNT, GEOLOW, CORGEO, GEOCOR, SHAG, RIGHT, +C CGM coordinates : GEOCGM01, OVL_ANG, CGMGLA, CGMGLO, DFR1DR, +C AZM_ANG, MLTUT, MFC, FTPRNT, GEOLOW, CORGEO, GEOCOR, SHAG, RIGHT, C IGRF, RECALC, SPHCAR, BSPCAR, GEOMAG, MAGSM, SMGSM C C MLT: CLCMLT, DPMTRX c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -c Required i/o units: +c Required i/o units: c KONSOL= 6 Program messages (used when jf(12)=.true. -> konsol) c KONSOL=11 Program messages (used when jf(12)=.false. -> MESSAGES.TXT) c -c COMMON/iounit/konsol,mess is used to pass the value of KONSOL from +c COMMON/iounit/konsol,mess is used to pass the value of KONSOL from c IRISUB to IRIFUN and IGRF. If mess=false then messages are turned off. -c +c c UNIT=14 IGRF/GETSHC: IGRF coeff. (DGRF%%%%.DAT or IGRF%%%%.DAT, %%%%=year) c- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C Corrections: -C 11/01/91 SHELLG: lowest starting point for B0 search is 2 +C 11/01/91 SHELLG: lowest starting point for B0 search is 2 C 1/27/92 Adopted to IGRF-91 coeffcients model C 2/05/92 Reduce variable names: INTER(P)SHC,EXTRA(P)SHC,INITI(ALI)ZE C 8/08/95 Updated to IGRF-45-95; new coeff. DGRF90, IGRF95, IGRF95S @@ -35,27 +35,27 @@ c 2000.01 10/28/02 replace TAB/6 blanks, enforce 72/line (D. Simpson) C 2000.02 11/08/02 change unit for coefficients to 14 C 2000.03 06/05/03 correct DIPL computation (V. Truhlik) C 2005.00 04/25/05 CALL FELDI and DO 1111 I=1,7 (Alexey Petrov) -C 2005.01 11/10/05 added igrf_dip and geodip (MLAT) +C 2005.01 11/10/05 added igrf_dip and geodip (MLAT) C 2005.02 11/10/05 FELDCOF: updated to IGRF-10 version C 2005.03 12/21/06 GH2(120) -> GH2(144) C 2007.00 05/18/07 Release of IRI-2007 -C 2007.08 07/30/09 SHELLG,STOER,FELDG,FELDCOF: NMAX=13; H/G-arrays(195) +C 2007.08 07/30/09 SHELLG,STOER,FELDG,FELDCOF: NMAX=13; H/G-arrays(195) C 2007.10 02/26/10 FELDCOF: updated to IGRF-11; DGRF05, IGRF10, IGRF10S C 2007.11 04/27/10 RECALC: updated to IGRF-11 -C 2007.11 04/27/10 Make all arrays(195) to arrays(196) +C 2007.11 04/27/10 Make all arrays(195) to arrays(196) C 2007.11 04/27/10 FELDCOF: corrected Filmod and also IGRF10.DAT C 2007.11 04/29/10 New files dgrf%%%%.asc; new GETSHC; char*12 to 13 C C 2012.00 10/05/11 IRI-2012: bottomside B0 B1 model (SHAMDB0D, SHAB1D), C 2012.00 10/05/11 bottomside Ni model (iriflip.for), auroral foE C 2012.00 10/05/11 storm model (storme_ap), Te with PF10.7 (elteik), -C 2012.00 10/05/11 oval kp model (auroral_boundary), IGRF-11(igrf.for), +C 2012.00 10/05/11 oval kp model (auroral_boundary), IGRF-11(igrf.for), C 2012.00 10/05/11 NRLMSIS00 (cira.for), CGM coordinates, F10.7 daily C 2012.00 10/05/11 81-day 365-day indices (apf107.dat), ap->kp (ckp), C 2012.00 10/05/11 array size change jf(50) outf(20,1000), oarr(100). C 2012.02 12/17/12 igrf_dip: Add magnetic declination as output parameter C 2014.01 07/20/14 igrf_dip,FTPRNT,RECALC: ASIN(x): abs(x)>1.0 x=sign(1.,x) -C 2014.02 07/24/14 COMMON/iounit: added 'mess' +C 2014.02 07/24/14 COMMON/iounit: added 'mess' C 2015.01 02/10/15 Updating to IGRF-12 (2015) C 2015.01 07/12/15 use mess,konsol in IGRF and RECALC C 2015.02 08/23/15 initialization of Earth constants moved to IRI_SUB @@ -64,16 +64,16 @@ C 2015.03 10/14/15 RECALC: update with IGRF-12 until 2020 C 2015.03 10/14/15 IGRF_SUB,_DIP: move CALL FELDCOF to IRISUB.FOR C 2015.03 10/14/15 FELDCOF,SHELLG: DIMO to COMMON/IGRF1/ C 2016.01 02/17/16 GEODIP: add PI to CONST -c----------------------------------------------------------------------- -C +c----------------------------------------------------------------------- +C subroutine igrf_sub(xlat,xlong,year,height, & xl,icode,dipl,babs) -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- c INPUT: c xlat geodatic latitude in degrees c xlong geodatic longitude in degrees -c year decimal year (year+(month-0.5)/12.0-0.5 or -c year+day-of-year/365 or ../366 if leap year) +c year decimal year (year+(month-0.5)/12.0-0.5 or +c year+day-of-year/365 or ../366 if leap year) c height height in km c OUTPUT: c xl L value @@ -81,11 +81,11 @@ c icode =1 L is correct; =2 L is not correct; c =3 an approximation is used c dipl dip latitude in degrees c babs magnetic field strength in Gauss -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- REAL LATI,LONGI COMMON /CONST/UMR,PI - + lati=xlat longi=xlong c CALL FELDCOF(YEAR,DIMO) @@ -97,19 +97,19 @@ c CALL FELDCOF(YEAR,DIMO) c c subroutine igrf_dip(xlat,xlong,year,height,dec,dip,dipl,ymodip) -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- c INPUT: c xlat geodatic latitude in degrees c xlong geodatic longitude in degrees -c year decimal year (year+month/12.0-0.5 or -c year+day-of-year/365 or ../366 if leap year) +c year decimal year (year+month/12.0-0.5 or +c year+day-of-year/365 or ../366 if leap year) c height height in km c OUTPUT: c dec magnetic declination in degrees c dip magnetic inclination (dip) in degrees c dipl dip latitude in degrees -c ymodip modified dip latitude = asin{dip/sqrt[dip^2+cos(LATI)]} -c----------------------------------------------------------------------- +c ymodip modified dip latitude = asin{dip/sqrt[dip^2+cos(LATI)]} +c----------------------------------------------------------------------- COMMON /CONST/UMR,PI @@ -126,10 +126,10 @@ c CALL FELDCOF(YEAR,DIMO) DIP=ASIN(BDBA) dipdiv=DIP/SQRT(DIP*DIP+cos(XLATI*UMR)) IF(ABS(dipdiv).GT.1.) dipdiv=SIGN(1.,dipdiv) - SMODIP=ASIN(dipdiv) + SMODIP=ASIN(dipdiv) c DIPL1=ATAN(0.5*TAN(DIP))/UMR DIPL=ATAN(BDOWN/2.0/sqrt(BNORTH*BNORTH+BEAST*BEAST))/umr - YMODIP=SMODIP/UMR + YMODIP=SMODIP/UMR DEC=DEC/UMR DIP=DIP/UMR RETURN @@ -138,7 +138,7 @@ c c C SHELLIG.FOR C -C 11/01/91 SHELLG: lowest starting point for B0 search is 2 +C 11/01/91 SHELLG: lowest starting point for B0 search is 2 C 1/27/92 Adopted to IGRF-91 coeffcients model C 2/05/92 Reduce variable-names: INTER(P)SHC,EXTRA(P)SHC,INITI(ALI)ZE C 8/08/95 Updated to IGRF-45-95; new coeff. DGRF90, IGRF95, IGRF95S @@ -160,18 +160,18 @@ C C SUBROUTINE SHELLG(GLAT,GLON,ALT,FL,ICODE,B0) c SUBROUTINE SHELLG(GLAT,GLON,ALT,DIMO,FL,ICODE,B0) -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- C CALCULATES L-VALUE FOR SPECIFIED GEODAETIC COORDINATES, ALTITUDE C AND GEMAGNETIC FIELD MODEL. -C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTER, INTERNAL NOTE +C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTER, INTERNAL NOTE C NO. 67, 1970. C G. KLUGE, COMPUTER PHYSICS COMMUNICATIONS 3, 31-35, 1972 -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- C CHANGES (D. BILITZA, NOV 87): C - USING CORRECT DIPOL MOMENT I.E.,DIFFERENT COMMON/MODEL/ C - USING IGRF EARTH MAGNETIC FIELD MODELS FROM 1945 TO 1990 C 09/07/22 NMAX=13 for DGRF00 and IGRF05; H/G-arrays(195) -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- C INPUT: ENTRY POINT SHELLG C GLAT GEODETIC LATITUDE IN DEGREES (NORTH) C GLON GEODETIC LONGITUDE IN DEGREES (EAST) @@ -183,12 +183,12 @@ C X-AXIS POINTING TO EQUATOR AT 0 LONGITUDE C Y-AXIS POINTING TO EQUATOR AT 90 LONG. C Z-AXIS POINTING TO NORTH POLE C -C DIMO DIPOL MOMENT IN GAUSS (NORMALIZED TO EARTH RADIUS) +C DIMO DIPOL MOMENT IN GAUSS (NORMALIZED TO EARTH RADIUS) C -C COMMON +C COMMON C X(3) NOT USED C H(144) FIELD MODEL COEFFICIENTS ADJUSTED FOR SHELLG -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- C OUTPUT: FL L-VALUE C ICODE =1 NORMAL COMPLETION C =2 UNPHYSICAL CONJUGATE POINT (FL MEANINGLESS) @@ -196,69 +196,69 @@ C =3 SHELL PARAMETER GREATER THAN LIMIT UP TO C WHICH ACCURATE CALCULATION IS REQUIRED; C APPROXIMATION IS USED. C B0 MAGNETIC FIELD STRENGTH IN GAUSS -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- DIMENSION V(3),U(3,3),P(8,100),SP(3) COMMON/IGRF2/ X(3),H(196) - COMMON/FIDB0/ SP /CONST/UMR,PI + COMMON/FIDB0/ SP /CONST/UMR,PI COMMON/IGRF1/ ERA,AQUAD,BQUAD,DIMO C C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3 C-- STEP IS STEP SIZE FOR FIELD LINE TRACING C-- STEQ IS STEP SIZE FOR INTEGRATION -C +C DATA RMIN,RMAX /0.05,1.01/ DATA STEP,STEQ /0.20,0.03/ BEQU=1.E10 C*****ENTRY POINT SHELLG TO BE USED WITH GEODETIC CO-ORDINATES RLAT=GLAT*UMR - CT=SIN(RLAT) - ST=COS(RLAT) + CT=SIN(RLAT) + ST=COS(RLAT) D=SQRT(AQUAD-(AQUAD-BQUAD)*CT*CT) X(1)=(ALT+AQUAD/D)*ST/ERA X(3)=(ALT+BQUAD/D)*CT/ERA RLON=GLON*UMR - X(2)=X(1)*SIN(RLON) - X(1)=X(1)*COS(RLON) - GOTO9 - ENTRY SHELLC(V,FL,B0) + X(2)=X(1)*SIN(RLON) + X(1)=X(1)*COS(RLON) + GOTO9 + ENTRY SHELLC(V,FL,B0) C*****ENTRY POINT SHELLC TO BE USED WITH CARTESIAN CO-ORDINATES - X(1)=V(1) - X(2)=V(2) - X(3)=V(3) -C*****CONVERT TO DIPOL-ORIENTED CO-ORDINATES - DATA U/ +0.3511737,-0.9148385,-0.1993679, - A +0.9335804,+0.3583680,+0.0000000, - B +0.0714471,-0.1861260,+0.9799247/ + X(1)=V(1) + X(2)=V(2) + X(3)=V(3) +C*****CONVERT TO DIPOL-ORIENTED CO-ORDINATES + DATA U/ +0.3511737,-0.9148385,-0.1993679, + A +0.9335804,+0.3583680,+0.0000000, + B +0.0714471,-0.1861260,+0.9799247/ 9 RQ=1./(X(1)*X(1)+X(2)*X(2)+X(3)*X(3)) - R3H=SQRT(RQ*SQRT(RQ)) - P(1,2)=(X(1)*U(1,1)+X(2)*U(2,1)+X(3)*U(3,1))*R3H - P(2,2)=(X(1)*U(1,2)+X(2)*U(2,2) )*R3H - P(3,2)=(X(1)*U(1,3)+X(2)*U(2,3)+X(3)*U(3,3))*RQ -C*****FIRST THREE POINTS OF FIELD LINE - STEP=-SIGN(STEP,P(3,2)) - CALL STOER(P(1,2),BQ2,R2) - B0=SQRT(BQ2) - P(1,3)=P(1,2)+0.5*STEP*P(4,2) - P(2,3)=P(2,2)+0.5*STEP*P(5,2) - P(3,3)=P(3,2)+0.5*STEP - CALL STOER(P(1,3),BQ3,R3) - P(1,1)=P(1,2)-STEP*(2.*P(4,2)-P(4,3)) - P(2,1)=P(2,2)-STEP*(2.*P(5,2)-P(5,3)) - P(3,1)=P(3,2)-STEP - CALL STOER(P(1,1),BQ1,R1) - P(1,3)=P(1,2)+STEP*(20.*P(4,3)-3.*P(4,2)+P(4,1))/18. - P(2,3)=P(2,2)+STEP*(20.*P(5,3)-3.*P(5,2)+P(5,1))/18. - P(3,3)=P(3,2)+STEP - CALL STOER(P(1,3),BQ3,R3) -C*****INVERT SENSE IF REQUIRED - IF(BQ3.LE.BQ1)GOTO2 - STEP=-STEP - R3=R1 - BQ3=BQ1 - DO 1 I=1,7 - ZZ=P(I,1) - P(I,1)=P(I,3) -1 P(I,3)=ZZ + R3H=SQRT(RQ*SQRT(RQ)) + P(1,2)=(X(1)*U(1,1)+X(2)*U(2,1)+X(3)*U(3,1))*R3H + P(2,2)=(X(1)*U(1,2)+X(2)*U(2,2) )*R3H + P(3,2)=(X(1)*U(1,3)+X(2)*U(2,3)+X(3)*U(3,3))*RQ +C*****FIRST THREE POINTS OF FIELD LINE + STEP=-SIGN(STEP,P(3,2)) + CALL STOER(P(1,2),BQ2,R2) + B0=SQRT(BQ2) + P(1,3)=P(1,2)+0.5*STEP*P(4,2) + P(2,3)=P(2,2)+0.5*STEP*P(5,2) + P(3,3)=P(3,2)+0.5*STEP + CALL STOER(P(1,3),BQ3,R3) + P(1,1)=P(1,2)-STEP*(2.*P(4,2)-P(4,3)) + P(2,1)=P(2,2)-STEP*(2.*P(5,2)-P(5,3)) + P(3,1)=P(3,2)-STEP + CALL STOER(P(1,1),BQ1,R1) + P(1,3)=P(1,2)+STEP*(20.*P(4,3)-3.*P(4,2)+P(4,1))/18. + P(2,3)=P(2,2)+STEP*(20.*P(5,3)-3.*P(5,2)+P(5,1))/18. + P(3,3)=P(3,2)+STEP + CALL STOER(P(1,3),BQ3,R3) +C*****INVERT SENSE IF REQUIRED + IF(BQ3.LE.BQ1)GOTO2 + STEP=-STEP + R3=R1 + BQ3=BQ1 + DO 1 I=1,7 + ZZ=P(I,1) + P(I,1)=P(I,3) +1 P(I,3)=ZZ C*****SEARCH FOR LOWEST MAGNETIC FIELD STRENGTH 2 IF(BQ1.LT.BEQU) THEN BEQU=BQ1 @@ -272,83 +272,83 @@ C*****SEARCH FOR LOWEST MAGNETIC FIELD STRENGTH BEQU=BQ3 IEQU=3 ENDIF -C*****INITIALIZATION OF INTEGRATION LOOPS +C*****INITIALIZATION OF INTEGRATION LOOPS STEP12=STEP/12. - STEP2=STEP+STEP - STEQ=SIGN(STEQ,STEP) - FI=0. - ICODE=1 - ORADIK=0. - OTERM=0. - STP=R2*STEQ - Z=P(3,2)+STP + STEP2=STEP+STEP + STEQ=SIGN(STEQ,STEP) + FI=0. + ICODE=1 + ORADIK=0. + OTERM=0. + STP=R2*STEQ + Z=P(3,2)+STP STP=STP/0.75 - P(8,1)=STEP2*(P(1,1)*P(4,1)+P(2,1)*P(5,1)) + P(8,1)=STEP2*(P(1,1)*P(4,1)+P(2,1)*P(5,1)) P(8,2)=STEP2*(P(1,2)*P(4,2)+P(2,2)*P(5,2)) -C*****MAIN LOOP (FIELD LINE TRACING) - DO 3 N=3,3333 -C*****CORRECTOR (FIELD LINE TRACING) - P(1,N)=P(1,N-1)+STEP12*(5.*P(4,N)+8.*P(4,N-1)-P(4,N-2)) - P(2,N)=P(2,N-1)+STEP12*(5.*P(5,N)+8.*P(5,N-1)-P(5,N-2)) -C*****PREPARE EXPANSION COEFFICIENTS FOR INTERPOLATION -C*****OF SLOWLY VARYING QUANTITIES - P(8,N)=STEP2*(P(1,N)*P(4,N)+P(2,N)*P(5,N)) - C0=P(1,N-1)**2+P(2,N-1)**2 - C1=P(8,N-1) - C2=(P(8,N)-P(8,N-2))*0.25 +C*****MAIN LOOP (FIELD LINE TRACING) + DO 3 N=3,3333 +C*****CORRECTOR (FIELD LINE TRACING) + P(1,N)=P(1,N-1)+STEP12*(5.*P(4,N)+8.*P(4,N-1)-P(4,N-2)) + P(2,N)=P(2,N-1)+STEP12*(5.*P(5,N)+8.*P(5,N-1)-P(5,N-2)) +C*****PREPARE EXPANSION COEFFICIENTS FOR INTERPOLATION +C*****OF SLOWLY VARYING QUANTITIES + P(8,N)=STEP2*(P(1,N)*P(4,N)+P(2,N)*P(5,N)) + C0=P(1,N-1)**2+P(2,N-1)**2 + C1=P(8,N-1) + C2=(P(8,N)-P(8,N-2))*0.25 C3=(P(8,N)+P(8,N-2)-C1-C1)/6.0 - D0=P(6,N-1) + D0=P(6,N-1) D1=(P(6,N)-P(6,N-2))*0.5 - D2=(P(6,N)+P(6,N-2)-D0-D0)*0.5 + D2=(P(6,N)+P(6,N-2)-D0-D0)*0.5 E0=P(7,N-1) - E1=(P(7,N)-P(7,N-2))*0.5 - E2=(P(7,N)+P(7,N-2)-E0-E0)*0.5 -C*****INNER LOOP (FOR QUADRATURE) -4 T=(Z-P(3,N-1))/STEP - IF(T.GT.1.)GOTO5 - HLI=0.5*(((C3*T+C2)*T+C1)*T+C0) + E1=(P(7,N)-P(7,N-2))*0.5 + E2=(P(7,N)+P(7,N-2)-E0-E0)*0.5 +C*****INNER LOOP (FOR QUADRATURE) +4 T=(Z-P(3,N-1))/STEP + IF(T.GT.1.)GOTO5 + HLI=0.5*(((C3*T+C2)*T+C1)*T+C0) ZQ=Z*Z R=HLI+SQRT(HLI*HLI+ZQ) - IF(R.LE.RMIN)GOTO30 + IF(R.LE.RMIN)GOTO30 RQ=R*R - FF=SQRT(1.+3.*ZQ/RQ) - RADIK=B0-((D2*T+D1)*T+D0)*R*RQ*FF - IF(R-RMAX)44,44,45 -45 ICODE=2 - RADIK=RADIK-12.*(R-RMAX)**2 + FF=SQRT(1.+3.*ZQ/RQ) + RADIK=B0-((D2*T+D1)*T+D0)*R*RQ*FF + IF(R-RMAX)44,44,45 +45 ICODE=2 + RADIK=RADIK-12.*(R-RMAX)**2 44 IF(RADIK+RADIK.LE.ORADIK) GOTO 10 - TERM=SQRT(RADIK)*FF*((E2*T+E1)*T+E0)/(RQ+ZQ) - FI=FI+STP*(OTERM+TERM) - ORADIK=RADIK - OTERM=TERM - STP=R*STEQ - Z=Z+STP - GOTO4 -C*****PREDICTOR (FIELD LINE TRACING) -5 P(1,N+1)=P(1,N)+STEP12*(23.*P(4,N)-16.*P(4,N-1)+5.*P(4,N-2)) - P(2,N+1)=P(2,N)+STEP12*(23.*P(5,N)-16.*P(5,N-1)+5.*P(5,N-2)) - P(3,N+1)=P(3,N)+STEP - CALL STOER(P(1,N+1),BQ3,R3) + TERM=SQRT(RADIK)*FF*((E2*T+E1)*T+E0)/(RQ+ZQ) + FI=FI+STP*(OTERM+TERM) + ORADIK=RADIK + OTERM=TERM + STP=R*STEQ + Z=Z+STP + GOTO4 +C*****PREDICTOR (FIELD LINE TRACING) +5 P(1,N+1)=P(1,N)+STEP12*(23.*P(4,N)-16.*P(4,N-1)+5.*P(4,N-2)) + P(2,N+1)=P(2,N)+STEP12*(23.*P(5,N)-16.*P(5,N-1)+5.*P(5,N-2)) + P(3,N+1)=P(3,N)+STEP + CALL STOER(P(1,N+1),BQ3,R3) C*****SEARCH FOR LOWEST MAGNETIC FIELD STRENGTH IF(BQ3.LT.BEQU) THEN IEQU=N+1 BEQU=BQ3 ENDIF 3 CONTINUE -10 IF(IEQU.lt.2) IEQU=2 +10 IF(IEQU.lt.2) IEQU=2 SP(1)=P(1,IEQU-1) SP(2)=P(2,IEQU-1) SP(3)=P(3,IEQU-1) - IF(ORADIK.LT.1E-15)GOTO11 - FI=FI+STP/0.75*OTERM*ORADIK/(ORADIK-RADIK) + IF(ORADIK.LT.1E-15)GOTO11 + FI=FI+STP/0.75*OTERM*ORADIK/(ORADIK-RADIK) C C-- The minimal allowable value of FI was changed from 1E-15 to 1E-12, C-- because 1E-38 is the minimal allowable arg. for ALOG in our envir. C-- D. Bilitza, Nov 87. C -11 FI=0.5*ABS(FI)/SQRT(B0)+1E-12 +11 FI=0.5*ABS(FI)/SQRT(B0)+1E-12 C -C*****COMPUTE L FROM B AND I. SAME AS CARMEL IN INVAR. +C*****COMPUTE L FROM B AND I. SAME AS CARMEL IN INVAR. C C-- Correct dipole moment is used here. D. Bilitza, Nov 87. C @@ -358,44 +358,44 @@ C c arg = FI*FI*FI/DIMOB0 c if(abs(arg).gt.88.0) arg=88.0 XX=3*arg1-arg2 - IF(XX.GT.23.0) GOTO 776 - IF(XX.GT.11.7) GOTO 775 - IF(XX.GT.+3.0) GOTO 774 - IF(XX.GT.-3.0) GOTO 773 - IF(XX.GT.-22.) GOTO 772 - 771 GG=3.33338E-1*XX+3.0062102E-1 - GOTO777 - 772 GG=((((((((-8.1537735E-14*XX+8.3232531E-13)*XX+1.0066362E-9)*XX+ - & 8.1048663E-8)*XX+3.2916354E-6)*XX+8.2711096E-5)*XX+ + IF(XX.GT.23.0) GOTO 776 + IF(XX.GT.11.7) GOTO 775 + IF(XX.GT.+3.0) GOTO 774 + IF(XX.GT.-3.0) GOTO 773 + IF(XX.GT.-22.) GOTO 772 + 771 GG=3.33338E-1*XX+3.0062102E-1 + GOTO777 + 772 GG=((((((((-8.1537735E-14*XX+8.3232531E-13)*XX+1.0066362E-9)*XX+ + & 8.1048663E-8)*XX+3.2916354E-6)*XX+8.2711096E-5)*XX+ & 1.3714667E-3)*XX+1.5017245E-2)*XX+4.3432642E-1)*XX+ - & 6.2337691E-1 - GOTO777 - 773 GG=((((((((2.6047023E-10*XX+2.3028767E-9)*XX-2.1997983E-8)*XX- - & 5.3977642E-7)*XX-3.3408822E-6)*XX+3.8379917E-5)*XX+ + & 6.2337691E-1 + GOTO777 + 773 GG=((((((((2.6047023E-10*XX+2.3028767E-9)*XX-2.1997983E-8)*XX- + & 5.3977642E-7)*XX-3.3408822E-6)*XX+3.8379917E-5)*XX+ & 1.1784234E-3)*XX+1.4492441E-2)*XX+4.3352788E-1)*XX+ - & 6.228644E-1 - GOTO777 - 774 GG=((((((((6.3271665E-10*XX-3.958306E-8)*XX+9.9766148E-07)*XX- - & 1.2531932E-5)*XX+7.9451313E-5)*XX-3.2077032E-4)*XX+ + & 6.228644E-1 + GOTO777 + 774 GG=((((((((6.3271665E-10*XX-3.958306E-8)*XX+9.9766148E-07)*XX- + & 1.2531932E-5)*XX+7.9451313E-5)*XX-3.2077032E-4)*XX+ & 2.1680398E-3)*XX+1.2817956E-2)*XX+4.3510529E-1)*XX+ - & 6.222355E-1 - GOTO777 + & 6.222355E-1 + GOTO777 775 GG=(((((2.8212095E-8*XX-3.8049276E-6)*XX+2.170224E-4)*XX- & 6.7310339E-3)*XX+1.2038224E-1)*XX-1.8461796E-1)*XX+ - & 2.0007187E0 - GOTO777 - 776 GG=XX-3.0460681E0 + & 2.0007187E0 + GOTO777 + 776 GG=XX-3.0460681E0 777 FL=EXP(ALOG((1.+EXP(GG))*DIMOB0)/3.0) - RETURN -C*****APPROXIMATION FOR HIGH VALUES OF L. -30 ICODE=3 - T=-P(3,N-1)/STEP - FL=1./(ABS(((C3*T+C2)*T+C1)*T+C0)+1E-15) - RETURN - END + RETURN +C*****APPROXIMATION FOR HIGH VALUES OF L. +30 ICODE=3 + T=-P(3,N-1)/STEP + FL=1./(ABS(((C3*T+C2)*T+C1)*T+C0)+1E-15) + RETURN + END C C - SUBROUTINE STOER(P,BQ,R) + SUBROUTINE STOER(P,BQ,R) C******************************************************************* C* SUBROUTINE USED FOR FIELD LINE TRACING IN SHELLG * C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG * @@ -404,55 +404,55 @@ C 09/07/22 NMAX=13 for DGRF00 and IGRF05; H/G-arrays(195) C******************************************************************* DIMENSION P(7),U(3,3) COMMON/IGRF2/ XI(3),H(196) -C*****XM,YM,ZM ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES - ZM=P(3) +C*****XM,YM,ZM ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES + ZM=P(3) FLI=P(1)*P(1)+P(2)*P(2)+1E-15 R=0.5*(FLI+SQRT(FLI*FLI+(ZM+ZM)**2)) RQ=R*R - WR=SQRT(R) - XM=P(1)*WR - YM=P(2)*WR -C*****TRANSFORM TO GEOGRAPHIC CO-ORDINATE SYSTEM - DATA U/ +0.3511737,-0.9148385,-0.1993679, - A +0.9335804,+0.3583680,+0.0000000, - B +0.0714471,-0.1861260,+0.9799247/ - XI(1)=XM*U(1,1)+YM*U(1,2)+ZM*U(1,3) - XI(2)=XM*U(2,1)+YM*U(2,2)+ZM*U(2,3) - XI(3)=XM*U(3,1) +ZM*U(3,3) -C*****COMPUTE DERIVATIVES -c CALL FELDI(XI,H) + WR=SQRT(R) + XM=P(1)*WR + YM=P(2)*WR +C*****TRANSFORM TO GEOGRAPHIC CO-ORDINATE SYSTEM + DATA U/ +0.3511737,-0.9148385,-0.1993679, + A +0.9335804,+0.3583680,+0.0000000, + B +0.0714471,-0.1861260,+0.9799247/ + XI(1)=XM*U(1,1)+YM*U(1,2)+ZM*U(1,3) + XI(2)=XM*U(2,1)+YM*U(2,2)+ZM*U(2,3) + XI(3)=XM*U(3,1) +ZM*U(3,3) +C*****COMPUTE DERIVATIVES +c CALL FELDI(XI,H) CALL FELDI - Q=H(1)/RQ - DX=H(3)+H(3)+Q*XI(1) - DY=H(4)+H(4)+Q*XI(2) - DZ=H(2)+H(2)+Q*XI(3) -C*****TRANSFORM BACK TO GEOMAGNETIC CO-ORDINATE SYSTEM - DXM=U(1,1)*DX+U(2,1)*DY+U(3,1)*DZ - DYM=U(1,2)*DX+U(2,2)*DY - DZM=U(1,3)*DX+U(2,3)*DY+U(3,3)*DZ - DR=(XM*DXM+YM*DYM+ZM*DZM)/R -C*****FORM SLOWLY VARYING EXPRESSIONS - P(4)=(WR*DXM-0.5*P(1)*DR)/(R*DZM) - P(5)=(WR*DYM-0.5*P(2)*DR)/(R*DZM) + Q=H(1)/RQ + DX=H(3)+H(3)+Q*XI(1) + DY=H(4)+H(4)+Q*XI(2) + DZ=H(2)+H(2)+Q*XI(3) +C*****TRANSFORM BACK TO GEOMAGNETIC CO-ORDINATE SYSTEM + DXM=U(1,1)*DX+U(2,1)*DY+U(3,1)*DZ + DYM=U(1,2)*DX+U(2,2)*DY + DZM=U(1,3)*DX+U(2,3)*DY+U(3,3)*DZ + DR=(XM*DXM+YM*DYM+ZM*DZM)/R +C*****FORM SLOWLY VARYING EXPRESSIONS + P(4)=(WR*DXM-0.5*P(1)*DR)/(R*DZM) + P(5)=(WR*DYM-0.5*P(2)*DR)/(R*DZM) DSQ=RQ*(DXM*DXM+DYM*DYM+DZM*DZM) BQ=DSQ*RQ*RQ - P(6)=SQRT(DSQ/(RQ+3.*ZM*ZM)) - P(7)=P(6)*(RQ+ZM*ZM)/(RQ*DZM) - RETURN - END + P(6)=SQRT(DSQ/(RQ+3.*ZM*ZM)) + P(7)=P(6)*(RQ+ZM*ZM)/(RQ*DZM) + RETURN + END C C - SUBROUTINE FELDG(GLAT,GLON,ALT,BNORTH,BEAST,BDOWN,BABS) -c----------------------------------------------------------------------- + SUBROUTINE FELDG(GLAT,GLON,ALT,BNORTH,BEAST,BDOWN,BABS) +c----------------------------------------------------------------------- C CALCULATES EARTH MAGNETIC FIELD FROM SPHERICAL HARMONICS MODEL -C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTRE, INTERNAL NOTE 61, +C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTRE, INTERNAL NOTE 61, C 1970. -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- C CHANGES (D. BILITZA, NOV 87): C - FIELD COEFFICIENTS IN BINARY DATA FILES INSTEAD OF BLOCK DATA C - CALCULATES DIPOL MOMENT C 09/07/22 NMAX=13 for DGRF00 and IGRF05; H/G-arrays(195) -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- C INPUT: ENTRY POINT FELDG C GLAT GEODETIC LATITUDE IN DEGREES (NORTH) C GLON GEODETIC LONGITUDE IN DEGREES (EAST) @@ -466,119 +466,119 @@ C Z-AXIS POINTING TO NORTH POLE C C COMMON BLANK AND ENTRY POINT FELDI ARE NEEDED WHEN USED C IN CONNECTION WITH L-CALCULATION PROGRAM SHELLG. -C +C C COMMON /MODEL/ AND /IGRF1/ C UMR = ATAN(1.0)*4./180. *UMR= -C ERA EARTH RADIUS FOR NORMALIZATION OF CARTESIAN +C ERA EARTH RADIUS FOR NORMALIZATION OF CARTESIAN C COORDINATES (6371.2 KM) -C AQUAD, BQUAD SQUARE OF MAJOR AND MINOR HALF AXIS OF -C EARTH ELLIPSOID AS RECOMMENDED BY INTERNAT. +C AQUAD, BQUAD SQUARE OF MAJOR AND MINOR HALF AXIS OF +C EARTH ELLIPSOID AS RECOMMENDED BY INTERNAT. C ASTRONOMICAL UNION (6378.160, 6356.775 KM). C NMAX MAXIMUM ORDER OF SPHERICAL HARMONICS -C TIME YEAR (DECIMAL: 1973.5) FOR WHICH MAGNETIC +C TIME YEAR (DECIMAL: 1973.5) FOR WHICH MAGNETIC C FIELD IS TO BE CALCULATED C G(M) NORMALIZED FIELD COEFFICIENTS (SEE FELDCOF) C M=NMAX*(NMAX+2) -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- C OUTPUT: BABS MAGNETIC FIELD STRENGTH IN GAUSS C BNORTH, BEAST, BDOWN COMPONENTS OF THE FIELD WITH RESPECT C TO THE LOCAL GEODETIC COORDINATE SYSTEM, WITH AXIS C POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST -C AND DOWNWARD. +C AND DOWNWARD. C----------------------------------------------------------------------- - DIMENSION V(3),B(3) + DIMENSION V(3),B(3) CHARACTER*13 NAME COMMON/IGRF2/XI(3),H(196) - COMMON/MODEL/NMAX,TIME,G(196),NAME + COMMON/MODEL/NMAX,TIME,G(196),NAME COMMON/IGRF1/ERA,AQUAD,BQUAD,DIMO /CONST/UMR,PI C C-- IS RECORDS ENTRY POINT C -C*****ENTRY POINT FELDG TO BE USED WITH GEODETIC CO-ORDINATES - IS=1 +C*****ENTRY POINT FELDG TO BE USED WITH GEODETIC CO-ORDINATES + IS=1 RLAT=GLAT*UMR - CT=SIN(RLAT) - ST=COS(RLAT) - D=SQRT(AQUAD-(AQUAD-BQUAD)*CT*CT) + CT=SIN(RLAT) + ST=COS(RLAT) + D=SQRT(AQUAD-(AQUAD-BQUAD)*CT*CT) RLON=GLON*UMR - CP=COS(RLON) - SP=SIN(RLON) + CP=COS(RLON) + SP=SIN(RLON) ZZZ=(ALT+BQUAD/D)*CT/ERA RHO=(ALT+AQUAD/D)*ST/ERA - XXX=RHO*CP - YYY=RHO*SP - GOTO 10 - -C*****ENTRY POINT FELDC TO BE USED WITH CARTESIAN CO-ORDINATES - ENTRY FELDC(V,B) - IS=2 - XXX=V(1) - YYY=V(2) - ZZZ=V(3) -10 RQ=1./(XXX*XXX+YYY*YYY+ZZZ*ZZZ) - XI(1)=XXX*RQ - XI(2)=YYY*RQ - XI(3)=ZZZ*RQ - GOTO 20 - -C*****ENTRY POINT FELDI USED FOR L COMPUTATION - ENTRY FELDI - IS=3 -20 IHMAX=NMAX*NMAX+1 - LAST=IHMAX+NMAX+NMAX - IMAX=NMAX+NMAX-1 - DO 8 I=IHMAX,LAST -8 H(I)=G(I) - DO 6 K=1,3,2 - I=IMAX - IH=IHMAX -1 IL=IH-I - F=2./FLOAT(I-K+2) - X=XI(1)*F - Y=XI(2)*F - Z=XI(3)*(F+F) - I=I-2 + XXX=RHO*CP + YYY=RHO*SP + GOTO 10 + +C*****ENTRY POINT FELDC TO BE USED WITH CARTESIAN CO-ORDINATES + ENTRY FELDC(V,B) + IS=2 + XXX=V(1) + YYY=V(2) + ZZZ=V(3) +10 RQ=1./(XXX*XXX+YYY*YYY+ZZZ*ZZZ) + XI(1)=XXX*RQ + XI(2)=YYY*RQ + XI(3)=ZZZ*RQ + GOTO 20 + +C*****ENTRY POINT FELDI USED FOR L COMPUTATION + ENTRY FELDI + IS=3 +20 IHMAX=NMAX*NMAX+1 + LAST=IHMAX+NMAX+NMAX + IMAX=NMAX+NMAX-1 + DO 8 I=IHMAX,LAST +8 H(I)=G(I) + DO 6 K=1,3,2 + I=IMAX + IH=IHMAX +1 IL=IH-I + F=2./FLOAT(I-K+2) + X=XI(1)*F + Y=XI(2)*F + Z=XI(3)*(F+F) + I=I-2 IF(I-1) 5,4,2 - -2 DO 3 M=3,I,2 - H(IL+M+1)=G(IL+M+1)+Z*H(IH+M+1)+X*(H(IH+M+3)- - A H(IH+M-1))-Y*(H(IH+M+2)+H(IH+M-2)) -3 H(IL+M)=G(IL+M)+Z*H(IH+M)+X*(H(IH+M+2)- + +2 DO 3 M=3,I,2 + H(IL+M+1)=G(IL+M+1)+Z*H(IH+M+1)+X*(H(IH+M+3)- + A H(IH+M-1))-Y*(H(IH+M+2)+H(IH+M-2)) +3 H(IL+M)=G(IL+M)+Z*H(IH+M)+X*(H(IH+M+2)- A H(IH+M-2))+Y*(H(IH+M+3)+H(IH+M-1)) - -4 H(IL+2)=G(IL+2)+Z*H(IH+2)+X*H(IH+4)-Y*(H(IH+3)+H(IH)) - H(IL+1)=G(IL+1)+Z*H(IH+1)+Y*H(IH+4)+X*(H(IH+3)-H(IH)) -5 H(IL)=G(IL)+Z*H(IH)+2.*(X*H(IH+1)+Y*H(IH+2)) - IH=IL - IF(I.GE.K) GOTO 1 + +4 H(IL+2)=G(IL+2)+Z*H(IH+2)+X*H(IH+4)-Y*(H(IH+3)+H(IH)) + H(IL+1)=G(IL+1)+Z*H(IH+1)+Y*H(IH+4)+X*(H(IH+3)-H(IH)) +5 H(IL)=G(IL)+Z*H(IH)+2.*(X*H(IH+1)+Y*H(IH+2)) + IH=IL + IF(I.GE.K) GOTO 1 6 CONTINUE - - IF(IS.EQ.3) RETURN - - S=.5*H(1)+2.*(H(2)*XI(3)+H(3)*XI(1)+H(4)*XI(2)) - T=(RQ+RQ)*SQRT(RQ) - BXXX=T*(H(3)-S*XXX) - BYYY=T*(H(4)-S*YYY) - BZZZ=T*(H(2)-S*ZZZ) - - IF(IS.EQ.2) GOTO 7 + + IF(IS.EQ.3) RETURN + + S=.5*H(1)+2.*(H(2)*XI(3)+H(3)*XI(1)+H(4)*XI(2)) + T=(RQ+RQ)*SQRT(RQ) + BXXX=T*(H(3)-S*XXX) + BYYY=T*(H(4)-S*YYY) + BZZZ=T*(H(2)-S*ZZZ) + + IF(IS.EQ.2) GOTO 7 BABS=SQRT(BXXX*BXXX+BYYY*BYYY+BZZZ*BZZZ) - BEAST=BYYY*CP-BXXX*SP - BRHO=BYYY*SP+BXXX*CP - BNORTH=BZZZ*ST-BRHO*CT - BDOWN=-BZZZ*CT-BRHO*ST - RETURN - -7 B(1)=BXXX - B(2)=BYYY - B(3)=BZZZ - RETURN - END + BEAST=BYYY*CP-BXXX*SP + BRHO=BYYY*SP+BXXX*CP + BNORTH=BZZZ*ST-BRHO*CT + BDOWN=-BZZZ*CT-BRHO*ST + RETURN + +7 B(1)=BXXX + B(2)=BYYY + B(3)=BZZZ + RETURN + END C C SUBROUTINE FELDCOF(YEAR) -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- C DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS C C INPUT: YEAR DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO @@ -587,66 +587,65 @@ C COMMON/IGRF1/ERAD,AQUAD,BQUAD,DIMO /CONST/UMR,PI C OUTPUT: COMMON/MODEL/NMAX,TIME,GH1,FIL1 C COMMON/DIPOL/GHI1,GHI2,GHI3 C -C THE GEOMAGNETIC DIPOL MOMENT (DIMO) IN GAUSS (NORMALIZED TO EARTH'S +C THE GEOMAGNETIC DIPOL MOMENT (DIMO) IN GAUSS (NORMALIZED TO EARTH'S C RADIUS) AT THE TIME (YEAR) IS COMPUTED BUT NOT USED. C -C 05/31/2000 updated to IGRF-2000 version (###) -C 03/24/2000 updated to IGRF-2005 version (###) +C 05/31/2000 updated to IGRF-2000 version (###) +C 03/24/2000 updated to IGRF-2005 version (###) C 07/22/2009 NMAX=13 for DGRF00 and IGRF05; H/G-arrays(195) -C 02/26/2010 update to IGRF-11 (2010) (###) +C 02/26/2010 update to IGRF-11 (2010) (###) C 10/05/2011 added COMMON/DIPOL/ for MLT computation in DPMTRX (IRIFUN) C 02/10/2015 update to IGRF-12 (2015) (###) -c----------------------------------------------------------------------- - CHARACTER*13 FILMOD, FIL1, FIL2 +c----------------------------------------------------------------------- + real, intent(in) :: year + CHARACTER(13) FILMOD, FIL1, FIL2 C ### FILMOD, DTEMOD array-size is number of IGRF maps DIMENSION GH1(196),GH2(196),GHA(196),FILMOD(16) DIMENSION DTEMOD(16) - DOUBLE PRECISION X,F0,F + DOUBLE PRECISION X,F0,F COMMON/MODEL/ NMAX,TIME,GH1,FIL1 COMMON/IGRF1/ ERAD,AQUAD,BQUAD,DIMO /CONST/UMR,PI COMMON/DIPOL/ GHI1,GHI2,GHI3 C ### updated coefficient file names and corresponding years - DATA FILMOD / 'dgrf1945.dat','dgrf1950.dat','dgrf1955.dat', + DATA FILMOD / 'dgrf1945.dat','dgrf1950.dat','dgrf1955.dat', 1 'dgrf1960.dat','dgrf1965.dat','dgrf1970.dat','dgrf1975.dat', 2 'dgrf1980.dat','dgrf1985.dat','dgrf1990.dat','dgrf1995.dat', 3 'dgrf2000.dat','dgrf2005.dat','dgrf2010.dat','igrf2015.dat', 4 'igrf2015s.dat'/ - DATA DTEMOD / 1945., 1950., 1955., 1960., 1965., + DATA DTEMOD / 1945., 1950., 1955., 1960., 1965., 1 1970., 1975., 1980., 1985., 1990., 1995., 2000.,2005., - 2 2010., 2015., 2020./ + 2 2010., 2015., 2020./ C C ### numye is number of IGRF coefficient files minus 1 C NUMYE=15 C C IS=0 FOR SCHMIDT NORMALIZATION IS=1 GAUSS NORMALIZATION -C IU IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS -C - IU = 14 + IS = 0 C-- DETERMINE IGRF-YEARS FOR INPUT-YEAR TIME = YEAR IYEA = INT(YEAR/5.)*5 L = (IYEA - 1945)/5 + 1 IF(L.LT.1) L=1 - IF(L.GT.NUMYE) L=NUMYE - DTE1 = DTEMOD(L) - FIL1 = FILMOD(L) - DTE2 = DTEMOD(L+1) - FIL2 = FILMOD(L+1) + IF(L.GT.NUMYE) L=NUMYE + DTE1 = DTEMOD(L) + FIL1 = FILMOD(L) + DTE2 = DTEMOD(L+1) + FIL2 = FILMOD(L+1) C-- GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS - CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER) - IF (IER .NE. 0) STOP - CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER) - IF (IER .NE. 0) STOP + CALL GETSHC (FIL1, NMAX1, ERAD, GH1) + + CALL GETSHC (FIL2, NMAX2, ERAD, GH2) + C-- DETERMINE IGRF COEFFICIENTS FOR YEAR - IF (L .LE. NUMYE-1) THEN - CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2, - 1 NMAX2, GH2, NMAX, GHA) - ELSE - CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2, - 1 GH2, NMAX, GHA) - ENDIF + IF (L .LE. NUMYE-1) THEN + CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2, + 1 NMAX2, GH2, NMAX, GHA) + ELSE + CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2, + 1 GH2, NMAX, GHA) + ENDIF C-- DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G F0=0.D0 DO 1234 J=1,3 @@ -654,219 +653,235 @@ C-- DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G F0 = F0 + F * F 1234 CONTINUE DIMO = DSQRT(F0) - GHI1=GHA(1) - GHI2=GHA(2) + GHI1=GHA(1) + GHI2=GHA(2) GHI3=GHA(3) GH1(1) = 0.0 - I=2 - F0=1.D-5 - IF(IS.EQ.0) F0=-F0 - SQRT2=SQRT(2.) + I=2 + F0=1.D-5 + IF(IS.EQ.0) F0=-F0 + SQRT2=SQRT(2.) - DO 9 N=1,NMAX + DO 9 N=1,NMAX X = N - F0 = F0 * X * X / (4.D0 * X - 2.D0) + F0 = F0 * X * X / (4.D0 * X - 2.D0) IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X - F = F0 * 0.5D0 + F = F0 * 0.5D0 IF(IS.EQ.0) F = F * SQRT2 GH1(I) = GHA(I-1) * F0 - I = I+1 - DO 9 M=1,N - F = F * (X + M) / (X - M + 1.D0) - IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M)) + I = I+1 + DO 9 M=1,N + F = F * (X + M) / (X - M + 1.D0) + IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M)) GH1(I) = GHA(I-1) * F GH1(I+1) = GHA(I) * F I=I+2 -9 CONTINUE +9 CONTINUE RETURN END C C - SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER) -C =============================================================== -C Reads spherical harmonic coefficients from the specified -C file into an array. -C Input: -C IU - Logical unit number -C FSPEC - File specification -C Output: -C NMAX - Maximum degree and order of model -C ERAD - Earth's radius associated with the spherical -C harmonic coefficients, in the same units as -C elevation -C GH - Schmidt quasi-normal internal spherical -C harmonic coefficients -C IER - Error number: = 0, no error -C = -2, records out of order -C = FORTRAN run-time error number -C =============================================================== - - CHARACTER FSPEC*(*), FOUT*80 - character(256) dirdata1, filename + SUBROUTINE GETSHC (FSPEC, NMAX, ERAD, GH) +C =============================================================== +C Reads spherical harmonic coefficients from the specified +C file into an array. +C Input: +C IU - Logical unit number +C FSPEC - File specification +C Output: +C NMAX - Maximum degree and order of model +C ERAD - Earth's radius associated with the spherical +C harmonic coefficients, in the same units as +C elevation +C GH - Schmidt quasi-normal internal spherical +C harmonic coefficients +C IER - Error number: = 0, no error +C = -2, records out of order +C = FORTRAN run-time error number +C =============================================================== + CHARACTER(*), intent(in) :: FSPEC + integer, intent(out) :: nmax + real, intent(out) :: gh + + integer ier + character(80) :: FOUT + character(256) :: dirdata1, filename DIMENSION GH(196) - LOGICAL mess - COMMON/iounit/konsol,mess + LOGICAL mess + integer IU + + COMMON/iounit/konsol,mess common /folders/ dirdata1 - do 1 j=1,196 + do 1 j=1,196 1 GH(j)=0.0 -C --------------------------------------------------------------- -C Open coefficient file. Read past first header record. -C Read degree and order of model and Earth's radius. -C --------------------------------------------------------------- +C --------------------------------------------------------------- +C Open coefficient file. Read past first header record. +C Read degree and order of model and Earth's radius. +C --------------------------------------------------------------- WRITE(FOUT,667) FSPEC 667 FORMAT(A13) c-web-for webversion c 667 FORMAT('/var/www/omniweb/cgi/vitmo/IRI/',A13) - filename = trim(trim(trim(dirdata1) // '/igrf/') // trim(FOUT)) - OPEN (IU, FILE=filename, STATUS='OLD', IOSTAT=IER, ERR=999) - READ (IU, *, IOSTAT=IER, ERR=999) - READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD, XMYEAR - nm=nmax*(nmax+2) - READ (IU, *, IOSTAT=IER, ERR=999) (GH(i),i=1,nm) - goto 888 - -999 if (mess) write(konsol,100) filename -100 FORMAT('Error while reading ',A13) - -888 CLOSE (IU) - RETURN - END -C -C - SUBROUTINE INTERSHC (DATE, DTE1, NMAX1, GH1, DTE2, - 1 NMAX2, GH2, NMAX, GH) - -C =============================================================== -C -C Version 1.01 -C -C Interpolates linearly, in time, between two spherical -C harmonic models. -C -C Input: -C DATE - Date of resulting model (in decimal year) -C DTE1 - Date of earlier model -C NMAX1 - Maximum degree and order of earlier model -C GH1 - Schmidt quasi-normal internal spherical -C harmonic coefficients of earlier model -C DTE2 - Date of later model -C NMAX2 - Maximum degree and order of later model -C GH2 - Schmidt quasi-normal internal spherical -C harmonic coefficients of later model -C -C Output: -C GH - Coefficients of resulting model -C NMAX - Maximum degree and order of resulting model -C -C A. Zunde -C USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 -C -C =============================================================== - - DIMENSION GH1(*), GH2(*), GH(*) - -C --------------------------------------------------------------- -C The coefficients (GH) of the resulting model, at date -C DATE, are computed by linearly interpolating between the -C coefficients of the earlier model (GH1), at date DTE1, -C and those of the later model (GH2), at date DTE2. If one -C model is smaller than the other, the interpolation is -C performed with the missing coefficients assumed to be 0. -C --------------------------------------------------------------- - - FACTOR = (DATE - DTE1) / (DTE2 - DTE1) - - IF (NMAX1 .EQ. NMAX2) THEN - K = NMAX1 * (NMAX1 + 2) - NMAX = NMAX1 - ELSE IF (NMAX1 .GT. NMAX2) THEN - K = NMAX2 * (NMAX2 + 2) - L = NMAX1 * (NMAX1 + 2) - DO 1122 I = K + 1, L -1122 GH(I) = GH1(I) + FACTOR * (-GH1(I)) - NMAX = NMAX1 - ELSE - K = NMAX1 * (NMAX1 + 2) - L = NMAX2 * (NMAX2 + 2) - DO 1133 I = K + 1, L -1133 GH(I) = FACTOR * GH2(I) - NMAX = NMAX2 - ENDIF - - DO 1144 I = 1, K -1144 GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I)) - - RETURN - END -C -C - SUBROUTINE EXTRASHC (DATE, DTE1, NMAX1, GH1, NMAX2, - 1 GH2, NMAX, GH) - -C =============================================================== -C -C Version 1.01 -C -C Extrapolates linearly a spherical harmonic model with a -C rate-of-change model. -C -C Input: -C DATE - Date of resulting model (in decimal year) -C DTE1 - Date of base model -C NMAX1 - Maximum degree and order of base model -C GH1 - Schmidt quasi-normal internal spherical -C harmonic coefficients of base model -C NMAX2 - Maximum degree and order of rate-of-change -C model -C GH2 - Schmidt quasi-normal internal spherical -C harmonic coefficients of rate-of-change model -C -C Output: -C GH - Coefficients of resulting model -C NMAX - Maximum degree and order of resulting model -C -C A. Zunde -C USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 -C -C =============================================================== - - DIMENSION GH1(*), GH2(*), GH(*) - -C --------------------------------------------------------------- -C The coefficients (GH) of the resulting model, at date -C DATE, are computed by linearly extrapolating the coef- -C ficients of the base model (GH1), at date DTE1, using -C those of the rate-of-change model (GH2), at date DTE2. If -C one model is smaller than the other, the extrapolation is -C performed with the missing coefficients assumed to be 0. -C --------------------------------------------------------------- - - FACTOR = (DATE - DTE1) - - IF (NMAX1 .EQ. NMAX2) THEN - K = NMAX1 * (NMAX1 + 2) - NMAX = NMAX1 - ELSE IF (NMAX1 .GT. NMAX2) THEN - K = NMAX2 * (NMAX2 + 2) - L = NMAX1 * (NMAX1 + 2) - DO 1155 I = K + 1, L -1155 GH(I) = GH1(I) - NMAX = NMAX1 - ELSE - K = NMAX1 * (NMAX1 + 2) - L = NMAX2 * (NMAX2 + 2) - DO 1166 I = K + 1, L -1166 GH(I) = FACTOR * GH2(I) - NMAX = NMAX2 - ENDIF - - DO 1177 I = 1, K -1177 GH(I) = GH1(I) + FACTOR * GH2(I) - - RETURN - END + filename = trim(trim(dirdata1) // '/igrf/' // FOUT) + OPEN(newunit=IU, FILE=filename, STATUS='OLD', IOSTAT=IER) + if(ier/=0) then + write(konsol,*) 'Error opening ',filename, + & ' in ',dirdata1,' with ',fout + error stop + endif + + READ (IU, *, IOSTAT=IER) + READ (IU, *, IOSTAT=IER) NMAX, ERAD, XMYEAR + + nm = nmax*(nmax+2) + + READ (IU, *, IOSTAT=IER) (GH(i),i=1,nm) + + + if (ier/=0) then + write(konsol,'(A,A13)') 'Error while reading ', filename + error stop + endif + + close(iu) + + END SUBROUTINE GETSHC + + + SUBROUTINE INTERSHC (DATE, DTE1, NMAX1, GH1, DTE2, + 1 NMAX2, GH2, NMAX, GH) + +C =============================================================== +C +C Version 1.01 +C +C Interpolates linearly, in time, between two spherical +C harmonic models. +C +C Input: +C DATE - Date of resulting model (in decimal year) +C DTE1 - Date of earlier model +C NMAX1 - Maximum degree and order of earlier model +C GH1 - Schmidt quasi-normal internal spherical +C harmonic coefficients of earlier model +C DTE2 - Date of later model +C NMAX2 - Maximum degree and order of later model +C GH2 - Schmidt quasi-normal internal spherical +C harmonic coefficients of later model +C +C Output: +C GH - Coefficients of resulting model +C NMAX - Maximum degree and order of resulting model +C +C A. Zunde +C USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 +C +C =============================================================== + + DIMENSION GH1(*), GH2(*), GH(*) + +C --------------------------------------------------------------- +C The coefficients (GH) of the resulting model, at date +C DATE, are computed by linearly interpolating between the +C coefficients of the earlier model (GH1), at date DTE1, +C and those of the later model (GH2), at date DTE2. If one +C model is smaller than the other, the interpolation is +C performed with the missing coefficients assumed to be 0. +C --------------------------------------------------------------- + + FACTOR = (DATE - DTE1) / (DTE2 - DTE1) + + IF (NMAX1 .EQ. NMAX2) THEN + K = NMAX1 * (NMAX1 + 2) + NMAX = NMAX1 + ELSE IF (NMAX1 .GT. NMAX2) THEN + K = NMAX2 * (NMAX2 + 2) + L = NMAX1 * (NMAX1 + 2) + DO 1122 I = K + 1, L +1122 GH(I) = GH1(I) + FACTOR * (-GH1(I)) + NMAX = NMAX1 + ELSE + K = NMAX1 * (NMAX1 + 2) + L = NMAX2 * (NMAX2 + 2) + DO 1133 I = K + 1, L +1133 GH(I) = FACTOR * GH2(I) + NMAX = NMAX2 + ENDIF + + DO 1144 I = 1, K +1144 GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I)) + + RETURN + END +C +C + SUBROUTINE EXTRASHC (DATE, DTE1, NMAX1, GH1, NMAX2, + 1 GH2, NMAX, GH) + +C =============================================================== +C +C Version 1.01 +C +C Extrapolates linearly a spherical harmonic model with a +C rate-of-change model. +C +C Input: +C DATE - Date of resulting model (in decimal year) +C DTE1 - Date of base model +C NMAX1 - Maximum degree and order of base model +C GH1 - Schmidt quasi-normal internal spherical +C harmonic coefficients of base model +C NMAX2 - Maximum degree and order of rate-of-change +C model +C GH2 - Schmidt quasi-normal internal spherical +C harmonic coefficients of rate-of-change model +C +C Output: +C GH - Coefficients of resulting model +C NMAX - Maximum degree and order of resulting model +C +C A. Zunde +C USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225 +C +C =============================================================== + + DIMENSION GH1(*), GH2(*), GH(*) + +C --------------------------------------------------------------- +C The coefficients (GH) of the resulting model, at date +C DATE, are computed by linearly extrapolating the coef- +C ficients of the base model (GH1), at date DTE1, using +C those of the rate-of-change model (GH2), at date DTE2. If +C one model is smaller than the other, the extrapolation is +C performed with the missing coefficients assumed to be 0. +C --------------------------------------------------------------- + + FACTOR = (DATE - DTE1) + + IF (NMAX1 .EQ. NMAX2) THEN + K = NMAX1 * (NMAX1 + 2) + NMAX = NMAX1 + ELSE IF (NMAX1 .GT. NMAX2) THEN + K = NMAX2 * (NMAX2 + 2) + L = NMAX1 * (NMAX1 + 2) + DO 1155 I = K + 1, L +1155 GH(I) = GH1(I) + NMAX = NMAX1 + ELSE + K = NMAX1 * (NMAX1 + 2) + L = NMAX2 * (NMAX2 + 2) + DO 1166 I = K + 1, L +1166 GH(I) = FACTOR * GH2(I) + NMAX = NMAX2 + ENDIF + + DO 1177 I = 1, K +1177 GH(I) = GH1(I) + FACTOR * GH2(I) + + RETURN + END C C SUBROUTINE GEODIP(IYR,SLA,SLO,DLA,DLO,J) @@ -880,9 +895,9 @@ C OUTPUT: DLA,DLO SLA,SLO C Last revision: November 2005 (Vladimir Papitashvili) C The code is modifed from GEOCOR written by V.Popov and V.Papitashvili -C in mid-1980s. +C in mid-1980s. - COMMON /CONST/UMR,PI + COMMON /CONST/UMR,PI C Earth's radius (km) RE = 6371.2 @@ -891,7 +906,7 @@ C RH = (RE + HI)/RE R = 1. if(j.gt.0) goto 1234 - + COL = (90.- SLA)*UMR RLO = SLO*UMR CALL SPHCAR(R,COL,RLO,X,Y,Z,1) @@ -902,7 +917,7 @@ C RH = (RE + HI)/RE DCO = TH/UMR DLA = 90.- DCO RETURN - + 1234 continue COL = (90.- DLA)*UMR RLO = DLO*UMR @@ -916,14 +931,14 @@ C RH = (RE + HI)/RE END SUBROUTINE GEODIP -C -C +C +C function fmodip(xlat) real, intent(in) :: xlat - + common/findRLAT/xlong,year - + call igrf_dip(xlat,xlong,year,300.,dec,dip,dipl,ymodip) fmodip=ymodip @@ -934,7 +949,7 @@ C C ********************************************************************* C Version 2011 for GEO-CGM.FOR (good through 2015) January 2011 C Version 2005 for GEO-CGM.FOR (good through 2010) November 2005 -C Nov 11, 2005 IGRF and RECALC are is modified to the IGRF-10 model +C Nov 11, 2005 IGRF and RECALC are is modified to the IGRF-10 model C and extended back to 1900 using the DGRF coeffcients C Apr 11, 2001 GEOLOW is modified to account for interpolation of C CGM meridians near equator across the 360/0 boundary @@ -944,7 +959,7 @@ C NASA/Goddard Space Flight Center, Greenbelt, Maryland) C Vladimir O. Papitashvili (IZMIRAN, Moscow, Russia, now at SPRL, C University of Michigan, Ann Arbor) C Conributions from Boris A. Belov and Vladimir A. Popov (both at -C IZMIRAN), Therese Moretto (DMI, DSRI, now at NSF), Freddy +C IZMIRAN), Therese Moretto (DMI, DSRI, now at NSF), Freddy C Christiansen (DMI, DSRI), and Scott Boardsen (NASA/GSFC). C The original version of this code is described in the brochure by @@ -1346,7 +1361,7 @@ C longitude. If cr360 is true, geolon+360 deg is returned when geolon C is less than 90 deg. If cr0 is true, geolon-360 deg is returned C when geolon is larger than 270 degrees. C ********************************************************************* - real clon + real clon logical cr360,cr0 common/cgmgeo/cclat,cr360,cr0,rh @@ -1387,7 +1402,7 @@ C ********************************************************************** ! f2py does not understand functions calling functions on the fly ! example: f2py -m igrf -c irifun.for igrf.for skip: dfridr implicit none - real, external :: CGMfunc + real, external :: CGMfunc REAL, intent(in) :: h,x real :: err @@ -1396,13 +1411,13 @@ C ********************************************************************** LOGICAL mess Real,PARAMETER :: CON=1.4,CON2=CON*CON,BIG=1.E30,SAFE=2. Integer, Parameter :: NTAB=10 - - COMMON/iounit/konsol,mess + + COMMON/iounit/konsol,mess INTEGER i,j REAL errt,fac,hh,a(NTAB,NTAB) if(h.eq.0.) then - if (mess) write(konsol,100) + if (mess) write(konsol,100) 100 FORMAT('h must be nonzero in dfridr') return endif @@ -1478,10 +1493,9 @@ C ********************************************************************* alfa = atan2(sb,st) AZM_ANG = alfa/rad - RETURN END -C -C + + SUBROUTINE MLTUT(SLA,SLO,CLA,PLA,PLO,UT) C ********************************************************************* C Calculates the MLT midnight in UT hours @@ -2418,9 +2432,9 @@ C ********************************************************************* * H1990(66),H1995(66),H2000(66),H2005(66),H2010(66) logical mess - + common /iounit/konsol,mess - + DATA G1900/ * 0., -31543., -2298., -677., 2905., 924., 1022., * -1469., 1256., 572., 876., 628., 660., -361., @@ -2978,11 +2992,11 @@ C ********************************************************************* * 72.8, 68.6, 76.0, -141.4, -22.9, 13.1, -77.9, * 80.4, -75.0, -4.7, 45.3, 14.0, 10.4, 1.6, * 4.9, 24.3, 8.2, -14.5, -5.7, -19.3, 11.6, - * 10.9, -14.1, -3.7, 5.4, 9.4, 3.4, -5.3, + * 10.9, -14.1, -3.7, 5.4, 9.4, 3.4, -5.3, * 3.1, -12.4, -0.8, 8.4, -8.4, -10.1, -2.0, * -6.3, 0.9, -1.1, -0.2, 2.5, -0.3, 2.2, * 3.1, -1.0, -2.8/ - + DATA H2010/ * 0.0, 0.0, 4945.1, 0.0, -2707.7, -575.4, 0.0, @@ -3001,7 +3015,7 @@ C ********************************************************************* * 0.0, 11.4, 16.7, -11.3, -3.9, 2.7, 1.3, * -3.9, -2.9, -8.1, -1.4, 2.0, -8.9, 4.4, * -2.3, -0.5, 0.5, -1.5, -0.7, 1.3, 1.4, - * -0.3, -0.3, -0.3, 1.9, -1.6, -0.2, 1.8, + * -0.3, -0.3, -0.3, 1.9, -1.6, -0.2, 1.8, * 0.2, -0.1, -0.6, 1.4, 0.3, 0.1, -0.8, * 0.4, -0.1, 0.1, -0.5, 0.3, -0.3, 0.3, * 0.2, -0.5, 0.2/ @@ -3073,7 +3087,7 @@ C 40 CONTINUE GOTO 300 -C INTERPOLATE BETWEEN YEARS +C INTERPOLATE BETWEEN YEARS C INTERPOLATE BETWEEN 1900 - 1905: 1900 F2=(IYR-1900)/5. @@ -3083,7 +3097,7 @@ C INTERPOLATE BETWEEN 1900 - 1905: H(N)=H1900(N)*F1+H1905(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1905 - 1910: 1905 F2=(IYR-1905)/5. F1=1.-F2 @@ -3092,7 +3106,7 @@ C INTERPOLATE BETWEEN 1905 - 1910: H(N)=H1905(N)*F1+H1910(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1910 - 1915: 1910 F2=(IYR-1910)/5. F1=1.-F2 @@ -3101,7 +3115,7 @@ C INTERPOLATE BETWEEN 1910 - 1915: H(N)=H1910(N)*F1+H1915(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1915 - 1920: 1915 F2=(IYR-1915)/5. F1=1.-F2 @@ -3110,7 +3124,7 @@ C INTERPOLATE BETWEEN 1915 - 1920: H(N)=H1915(N)*F1+H1920(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1920 - 1925: 1920 F2=(IYR-1920)/5. F1=1.-F2 @@ -3119,7 +3133,7 @@ C INTERPOLATE BETWEEN 1920 - 1925: H(N)=H1920(N)*F1+H1925(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1925 - 1930: 1925 F2=(IYR-1925)/5. F1=1.-F2 @@ -3128,7 +3142,7 @@ C INTERPOLATE BETWEEN 1925 - 1930: H(N)=H1925(N)*F1+H1930(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1930 - 1935: 1930 F2=(IYR-1930)/5. F1=1.-F2 @@ -3137,7 +3151,7 @@ C INTERPOLATE BETWEEN 1930 - 1935: H(N)=H1930(N)*F1+H1935(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1935 - 1940: 1935 F2=(IYR-1935)/5. F1=1.-F2 @@ -3146,7 +3160,7 @@ C INTERPOLATE BETWEEN 1935 - 1940: H(N)=H1935(N)*F1+H1940(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1940 - 1945: 1940 F2=(IYR-1940)/5. F1=1.-F2 @@ -3155,7 +3169,7 @@ C INTERPOLATE BETWEEN 1940 - 1945: H(N)=H1940(N)*F1+H1945(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1945 - 1950: 1945 F2=(IYR-1945)/5. F1=1.-F2 @@ -3164,7 +3178,7 @@ C INTERPOLATE BETWEEN 1945 - 1950: H(N)=H1945(N)*F1+H1950(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1950 - 1955: 1950 F2=(IYR-1950)/5. F1=1.-F2 @@ -3173,7 +3187,7 @@ C INTERPOLATE BETWEEN 1950 - 1955: H(N)=H1950(N)*F1+H1955(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1955 - 1960: 1955 F2=(IYR-1955)/5. F1=1.-F2 @@ -3182,7 +3196,7 @@ C INTERPOLATE BETWEEN 1955 - 1960: H(N)=H1955(N)*F1+H1960(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1960 - 1965: 1960 F2=(IYR-1960)/5. F1=1.-F2 @@ -3191,7 +3205,7 @@ C INTERPOLATE BETWEEN 1960 - 1965: H(N)=H1960(N)*F1+H1965(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1965 - 1970: 1965 F2=(IYR-1965)/5. F1=1.-F2 @@ -3200,7 +3214,7 @@ C INTERPOLATE BETWEEN 1965 - 1970: H(N)=H1965(N)*F1+H1970(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1970 - 1975: 1970 F2=(IYR-1970)/5. F1=1.-F2 @@ -3209,7 +3223,7 @@ C INTERPOLATE BETWEEN 1970 - 1975: H(N)=H1970(N)*F1+H1975(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1975 - 1980: 1975 F2=(IYR-1975)/5. F1=1.-F2 @@ -3218,7 +3232,7 @@ C INTERPOLATE BETWEEN 1975 - 1980: H(N)=H1975(N)*F1+H1980(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1980 - 1985: 1980 F2=(IYR-1980)/5. F1=1.-F2 @@ -3227,7 +3241,7 @@ C INTERPOLATE BETWEEN 1980 - 1985: H(N)=H1980(N)*F1+H1985(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1985 - 1990: 1985 F2=(IYR-1985)/5. F1=1.-F2 @@ -3236,7 +3250,7 @@ C INTERPOLATE BETWEEN 1985 - 1990: H(N)=H1985(N)*F1+H1990(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1990 - 1995: 1990 F2=(IYR-1990)/5. F1=1.-F2 @@ -3245,7 +3259,7 @@ C INTERPOLATE BETWEEN 1990 - 1995: H(N)=H1990(N)*F1+H1995(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 1995 - 2000: 1995 F2=(IYR-1995)/5. F1=1.-F2 @@ -3254,7 +3268,7 @@ C INTERPOLATE BETWEEN 1995 - 2000: H(N)=H1995(N)*F1+H2000(N)*F2 ENDDO GOTO 300 - + C INTERPOLATE BETWEEN 2000 - 2005: 2000 F2=(IYR-2000)/5. F1=1.-F2 @@ -3397,7 +3411,7 @@ C from 1945 (updated by V. Papitashvili, February 1995) C C Modified to accept years up to 2005 (V. Papitashvili, January 2001) C -C Modified to accept years from 1900 through 2010 using the DGRF & +C Modified to accept years from 1900 through 2010 using the DGRF & C IGRF-10 coeficients (updated by V. Papitashvili, November 2005) C C Modified to accept years up to 2015 (V. Papitashvili, January 2011) @@ -3428,8 +3442,8 @@ C ********************************************************************* COMMON/C1/ ST0,CT0,SL0,CL0,CTCL,STCL,CTSL,STSL,SFI,CFI,SPS,CPS, * SHI,CHI,HI,PSI,XMUT,A11,A21,A31,A12,A22,A32,A13,A23,A33,DS3, * K,IY,BA - - common/iounit/konsol,mess + + common/iounit/konsol,mess DATA IYE,IDE/2*0/ IF (IYR.EQ.IYE.AND.IDAY.EQ.IDE) GOTO 5 @@ -3949,19 +3963,19 @@ C-------------------------------------------------------------------- INTEGER IHOUR,MIN,ISEC REAL GST,SLONG,SRASN,SDEC REAL BE,CAL,SA(3),S,C,SG(3),SM(3) - REAL LAM,LAMS,DELLAM + REAL LAM,LAMS,DELLAM COMMON /CONST/DTOR,PI - + XG=COS(GLAT*DTOR)*COS(GLON*DTOR) YG=COS(GLAT*DTOR)*SIN(GLON*DTOR) ZG=SIN(GLAT*DTOR) CALL DPMTRX(IYYYY,DDD,XXM,YYM,ZZM) - + C transform XM=XXM(1)*XG+XXM(2)*YG+XXM(3)*ZG YM=YYM(1)*XG+YYM(2)*YG+YYM(3)*ZG ZM=ZZM(1)*XG+ZZM(2)*YG+ZZM(3)*ZG -C +C IHOUR=INT(UTHR) MIN=INT((UTHR-IHOUR)*60) ISEC=INT((UTHR-IHOUR-MIN/60.0)*3600) @@ -3976,12 +3990,12 @@ C C=COS(BE) SG(1)=C*SA(1)+S*SA(2) SG(2)=C*SA(2)-S*SA(1) - SG(3)=SA(3) + SG(3)=SA(3) C transform SM(1)=XXM(1)*SG(1)+XXM(2)*SG(2)+XXM(3)*SG(3) SM(2)=YYM(1)*SG(1)+YYM(2)*SG(2)+YYM(3)*SG(3) SM(3)=ZZM(1)*SG(1)+ZZM(2)*SG(2)+ZZM(3)*SG(3) -C +C LAM=ATAN2(YM,XM) LAMS=ATAN2(SM(2),SM(1)) DELLAM=LAM-LAMS @@ -3993,7 +4007,7 @@ C C SUBROUTINE DPMTRX(IYYYY,DDD,XM,YM,ZM) C-------------------------------------------------------------------------- -C calculates othonormal matrix (columns XM,YM,ZM) for transformation +C calculates othonormal matrix (columns XM,YM,ZM) for transformation C from geographic to magnetic coordinates C Inputs: C IYYYY..year @@ -4001,7 +4015,7 @@ C DDD..day of year (1.1 = 0) C Outputs: C XM,YM,ZM..colums of the matrix C Notes: -C MX(N),MY(N),MZ(N)..coordinates of the B vector in geographic system +C MX(N),MY(N),MZ(N)..coordinates of the B vector in geographic system C for years stored in YR(N) C N..number of elements of arrays MX,MY,MZ and YR C-------------------------------------------------------------------------- diff --git a/fortran/irifun.for b/fortran/irifun.for index 2dd8b8d..51a93c2 100644 --- a/fortran/irifun.for +++ b/fortran/irifun.for @@ -5312,7 +5312,7 @@ C SYNTHESIZES THE VALUE OF HMF2 FROM THE MODEL C ********************************************************** CALL SCHNEVPDH (RZ,RLAT,FLON,dum,T,L,dum,dum,HMF2) RETURN -9999 STOP +9999 error STOP END SUBROUTINE SHAMDHMF2 C C @@ -5537,7 +5537,7 @@ c change to intercept form of fourier series. BE = Y BV = Z RETURN - 999 STOP + 999 error STOP END SUBROUTINE SCHNEVPDH C C @@ -5727,7 +5727,7 @@ c .. array arguments .. c .. local scalars .. integer coeff_month_read(1:12) character(256) filedata - integer i, j + integer i, j,u c .. local arrays .. double precision coeff_month_all(0:148,0:47,1:12) save coeff_month_all @@ -5738,12 +5738,12 @@ c C if (coeff_month_read(month) .eq. 0) then write(filedata, 10) month+10 - filename=trim(trim(trim(dirdata1)//'/mcsat/')//trim(filedata)) - open(10, File=filename, status='old') + filename=trim(trim(dirdata1)//'/mcsat/'//filedata) + open(newunit=u, File=filename, status='old') do j=0,47 - read(10,20) (coeff_month_all(i,j,month),i=0,148) + read(u,20) (coeff_month_all(i,j,month),i=0,148) end do - close(10) + close(u) coeff_month_read(month) = 1 end if c @@ -6719,7 +6719,7 @@ C SYNTHESIZES THE VALUE OF B0 FROM THE MODEL C ********************************************************** CALL SCHNEVPD(RZ,RLAT,FLON,dum,T,L,dum,dum,B) RETURN -9999 STOP +9999 error STOP END C C @@ -6876,7 +6876,7 @@ C ********************************************************** CALL SCHNEVPD(RZ,FLAT,FLON,dum,T,L,dum,dum,B) C RETURN -9999 STOP +9999 error STOP END C C @@ -7100,7 +7100,7 @@ c change to intercept form of fourier series. BE = Y BV = Z RETURN - 999 STOP + 999 error STOP END C C @@ -7216,7 +7216,7 @@ C NUMERICAL ERROR UNACCEPTABLY LARGE DUE TO ADDING OF C LARGE AND SMALL NUMBERS. if (mess) WRITE(konsol,108) M,FN,CONST,J,A,B 108 FORMAT (//12H ** ERROR **/1X,I5,F10.5,E15.7,I5,2D15.7) - STOP + error STOP C SERIES TRUNCATED SUCCESSFULLY. 110 PS = PNM DPS = DPNM @@ -8369,7 +8369,7 @@ C END c c - subroutine read_ig_rz + subroutine read_ig_rz() c---------------------------------------------------------------- c Reads the Rz12 and IG12 indices file IG_RZ.DAT from I/O UNIT=12 c and stores the indices in COMMON: @@ -8418,11 +8418,13 @@ c---------------------------------------------------------------- integer iyst,iyend,iymst,iupd,iupm,iupy,imst,imend real aig(806),arz(806) character(256) dirdata1,filename + integer u common /igrz/aig,arz,iymst,iymend common /folders/ dirdata1 filename = trim(trim(dirdata1) // '/index/ig_rz.dat' ) - open(unit=12,file=filename,FORM='FORMATTED',status='old') + print *,'reading ',filename + open(newunit=u,file=filename,FORM='FORMATTED',status='old') c-web- special for web version c open(unit=12,file= @@ -8432,8 +8434,8 @@ c * FORM='FORMATTED',status='old') c Read the update date, the start date and the end date (mm,yyyy), and c get number of data points to read. - read(12,*) iupd,iupm,iupy - read(12,*) imst,iyst,imend,iyend + read(u,*) iupd,iupm,iupy + read(u,*) imst,iyst,imend,iyend iymst=iyst*100+imst iymend=iyend*100+imend @@ -8443,8 +8445,8 @@ c 1st year \ full years \last y\ before & after inum_vals= 3-imst+(iyend-iyst)*12 +imend c read all the IG12 (ionoindx) and Rz12 (indrz) values - read(12,*) (aig(i),i=1,inum_vals) - read(12,*) (arz(i),i=1,inum_vals) + read(u,*) (aig(i),i=1,inum_vals) + read(u,*) (arz(i),i=1,inum_vals) c do 1 jj=1,inum_vals c rrr=arz(jj) c ggg=aig(jj) @@ -8462,13 +8464,12 @@ c arz(jj)=rrr c aig(jj)=ggg c1 continue - close(unit=12) + close(u) - return - end + end subroutine read_ig_rz c c - subroutine tcon(yr,mm,day,idn,rz,ig,rsn,nmonth) + subroutine tcon(yr,mm,day,idn,rz,ig,rsn,nmonth) c---------------------------------------------------------------- c input: yr,mm,day year(yyyy),month(mm),day(dd) c idn day of year(ddd) @@ -8487,25 +8488,32 @@ c month. The indices for the given day are obtained by linear c interpolation and are stored in rz(3) and ig(3). c---------------------------------------------------------------- - integer yr,mm,day,iyst,iyend,iymst - integer imst,iymend - real ionoindx(806),indrz(806) - real ig(3),rz(3) - logical mess + integer,intent(in) :: yr,mm,day + real, intent(out) :: rsn, ig(3),rz(3) + integer, intent(out) :: nmonth + + integer iyst,iyend,iymst + integer imst,iymend + real ionoindx(806),indrz(806) - common /iounit/konsol,mess - common /igrz/ionoindx,indrz,iymst,iymend + logical mess + + common /iounit/konsol,mess + common /igrz/ionoindx,indrz,iymst,iymend iytmp=yr*100+mm - if (iytmp.lt.iymst.or.iytmp.gt.iymend) then - if(mess) write(konsol,8000) iytmp,iymst,iymend -8000 format(1x,I10,'** OUT OF RANGE **'/,5x, + if (iytmp < iymst.or.iytmp > iymend) then + print '(A5,I4,A1,I2,A1,I2)','date ',yr,'-',mm,'-',day + write(konsol,8000) iytmp,iymst,iymend + +8000 format(1x,I10,'** OUT OF RANGE **'/,5x, & 'The file IG_RZ.DAT which contains the indices Rz12', & ' and IG12'/5x,'currently only covers the time period', & ' (yymm) : ',I6,'-',I6) - nmonth=-1 - return - endif + + nmonth=-1 + error stop + endif iyst=iymst/100 imst=iymst-iyst*100 @@ -8557,8 +8565,8 @@ c if((yr/4*4.eq.yr).and.(yr/100*100.ne.yr)) idd2=381 ig(3)=ig(2)+(ig(1)-ig(2))*rsn 1927 nmonth=imm2 - return - end + + end subroutine tcon C C subroutine readapf107 diff --git a/fortran/irisub.for b/fortran/irisub.for index 8c026b9..782ea98 100644 --- a/fortran/irisub.for +++ b/fortran/irisub.for @@ -1,14 +1,14 @@ c irisub.for, version number can be found at the end of this comment. -c----------------------------------------------------------------------- -C Includes subroutines IRI_SUB and IRI_WEB to compute IRI parameters -C for specified location, date, time, and altitude range and subroutine -C IRI_WEB to computes IRI parameters for specified location, date, time -C and variable range; variable can be altitude, latitude, longitude, -C year, month, day of month, day of year, or hour (UT or LT). -C IRI_WEB requires IRI_SUB. Both subroutines require linking with the -c following library files IRIFUN.FOR, IRITEC.FOR, IRIDREG.FOR, +c----------------------------------------------------------------------- +C Includes subroutines IRI_SUB and IRI_WEB to compute IRI parameters +C for specified location, date, time, and altitude range and subroutine +C IRI_WEB to computes IRI parameters for specified location, date, time +C and variable range; variable can be altitude, latitude, longitude, +C year, month, day of month, day of year, or hour (UT or LT). +C IRI_WEB requires IRI_SUB. Both subroutines require linking with the +c following library files IRIFUN.FOR, IRITEC.FOR, IRIDREG.FOR, c CIRA.FOR, IGRF.FOR -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- c Programs using subroutine IRI_SUB need to include (see IRITEST.FOR): c c call read_ig_rz @@ -20,18 +20,18 @@ c do i=1,100 c oar(i,1)=-1.0 c enddo c -c -c----------------------------------------------------------------------- -c Required i/o units: +c +c----------------------------------------------------------------------- +c Required i/o units: c KONSOL= 6 IRISUB: Program messages (used when jf(12)=.true. -> konsol) c IUCCIR=10 IRISUB: CCIR and URSI coefficients (CCIR%%.ASC, %%=month+10) c KONSOL=11 IRISUB: Program messages (used when jf(12)=.false. -> MESSAGES.TXT) -c KONSOL=6/11 is also used in IRIFUN and IGRF. COMMON/iounit/konsol,mess +c KONSOL=6/11 is also used in IRIFUN and IGRF. COMMON/iounit/konsol,mess c is used to pass the value of KONSOL. If mess=false messages are turned off. -c UNIT=12 IRIFUN/TCON: Solar/ionospheric indices IG12, R12 (IG_RZ.DAT) -c UNIT=13 IRIFUN/APF..: Magnetic indices and F10.7 (APF107.DAT +c UNIT=12 IRIFUN/TCON: Solar/ionospheric indices IG12, R12 (IG_RZ.DAT) +c UNIT=13 IRIFUN/APF..: Magnetic indices and F10.7 (APF107.DAT c UNIT=14 IGRF/GETSHC: IGRF coeff. (DGRF%%%%.DAT or IGRF%%%%.DAT, %%%%=year) -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- C CHANGES FROM IRIS11.FOR TO IRIS12.FOR: C - CIRA-1986 INSTEAD OF CIRA-1972 FOR NEUTRAL TEMPERATURE C - 10/30/91 VNER FOR NIGHTTIME LAY-VERSION: ABS(..) @@ -75,7 +75,7 @@ C - 5/29/99 the limit for IG comp. from Rz12-input is 174 not 274 C - 11/08/99 jf(18)=t simple UT to LT conversion, otherwise UT_LT C - 11/09/99 added COMMON/const1/humr,dumr also for CIRA86 C CHANGES FROM IRIS13.FOR TO IRISUB.FOR: -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- C-Version-MM/DD/YY-Description (person reporting correction) C 2000.01 05/09/00 B0_98 replaces B0_TAB and B1: 1.9/day to 2.6/night C 2000.02 06/11/00 including new F1 and indermediate region @@ -84,7 +84,7 @@ C 2000.04 10/29/00 include special option for D region models C 2000.05 12/07/00 change name IRIS13 to IRISUB C 2000.06 12/14/00 jf(30),outf(20,100),oarr(50) C 2000.07 03/17/01 include Truhlik-Triskova Te model and IGRF -C 2000.08 05/07/01 include Fuller-Rowell-Condrescu storm model +C 2000.08 05/07/01 include Fuller-Rowell-Condrescu storm model C 2000.09 07/09/01 LATI instead of LAT1 in F00 call -------- M. Torkar C 2000.10 07/09/01 sdte instead of dte in ELTEIK call --- P. Wilkinson C 2000.11 09/18/01 correct computation of foF2 for Rz12 user input @@ -98,27 +98,27 @@ C 2000.18 02/06/03 jf(27) for IG12 user input; all F1 prob in oar C 2000.19 07/14/04 covsat<188 instead of covsat= outf(500), numhei=numstp=500 +C 2007.02 10/31/08 outf(100) -> outf(500), numhei=numstp=500 C 2007.03 02/12/09 Jf(24)=.false.-> outf(1,60-140km)=FIRI- M. Friedrich C 2007.04 03/14/09 SOCO(70->80;500->300km) --------------- R. Davidson C 2007.05 03/26/09 call for APF_ONLY includes F107M C 2007.09 08/17/09 STROM off if input; fof2in, fof1in,foein corr C 2007.10 02/03/10 F10.7D = F10.7M = COV if EOF -C 2007.11 04/19/10 Corrections in irifun.for, cira.for -C 2007.12 11/23/10 FNIGHT computed twice at 8334 --------- C. Vasly +C 2007.11 04/19/10 Corrections in irifun.for, cira.for +C 2007.12 11/23/10 FNIGHT computed twice at 8334 --------- C. Vasly C C 2012.00 10/05/11 IRI-2012: bottomside B0 B1 model (SHAMDB0D, SHAB1D), C 2012.00 10/05/11 bottomside Ni model (iriflip.for), auroral foE C 2012.00 10/05/11 storm model (storme_ap), Te with PF10.7 (elteik), -C 2012.00 10/05/11 oval kp model (auroral_boundary),IGRF-11(igrf.for), +C 2012.00 10/05/11 oval kp model (auroral_boundary),IGRF-11(igrf.for), C 2012.00 10/05/11 NRLMSIS00 (cira.for), CGM coordinates, F10.7 daily C 2012.00 10/05/11 81-day 365-day indices (apf107.dat), ap->kp (ckp), C 2012.00 10/05/11 array size change jf(50) outf(20,1000), oarr(100). @@ -134,7 +134,7 @@ C 2012.03 02/13/13 Move B1 before B0 for Gulyaeva-1987 C 2012.03 02/20/13 Use foot-point for CGM to be closer to AACGM C 2012.03 02/20/13 DAT(11,*) is UT time of MLT=0 C 2012.04 09/12/13 Replace HOUR with HOURUT in APFMSIS ---- P. Coisson -C 2014.00 01/22/14 TMAXN in GTD7 SEC->SECNI HOUR->0.0 +C 2014.00 01/22/14 TMAXN in GTD7 SEC->SECNI HOUR->0.0 C 2014.01 07/17/14 Change estromcor to estormcor -------- A.Shabanloui C 2014.02 07/24/14 COMMON/iounit/: added 'mess' C 2014.03 09/18/14 JF(18): FIELDG not UT_LT ............... A.Mazzella @@ -173,7 +173,10 @@ C***************************************************************** C C SUBROUTINE IRI_SUB(JF,JMAG,ALATI,ALONG,IYYYY,MMDD,DHOUR, - & HEIBEG,HEIEND,HEISTP,OUTF,OARR) + & altkm,Nalt,datadir, + & OUTF,OARR) + + use, intrinsic:: iso_fortran_env, only: error_unit C----------------------------------------------------------------- C C INPUT: JF(1:50) true/false switches for several options @@ -181,7 +184,7 @@ C JMAG =0 geographic = 1 geomagnetic coordinates C ALATI,ALONG LATITUDE NORTH AND LONGITUDE EAST IN DEGREES C IYYYY Year as YYYY, e.g. 1985 C MMDD (-DDD) DATE (OR DAY OF YEAR AS A NEGATIVE NUMBER) -C DHOUR LOCAL TIME (OR UNIVERSAL TIME + 25) IN DECIMAL +C DHOUR LOCAL TIME (OR UNIVERSAL TIME + 25) IN DECIMAL C HOURS C HEIBEG, HEIGHT RANGE IN KM; maximal 100 heights, i.e. C HEIEND,HEISTP int((heiend-heibeg)/heistp)+1.le.100 @@ -196,7 +199,7 @@ C 3 Ne & Ni computed Ni not computed t C 4 B0,B1 - Bil-2000 B0,B1 - other models jf(31) false C 5 foF2 - CCIR foF2 - URSI false C 6 Ni - DS-1995 & DY-1985 Ni - RBV-2010 & TTS-2005 false -C 7 Ne - Tops: f10.7<188 f10.7 unlimited t +C 7 Ne - Tops: f10.7<188 f10.7 unlimited t C 8 foF2 from model foF2 or NmF2 - user input t C 9 hmF2 from model hmF2 or M3000F2 - user input t C 10 Te - Standard Te - Using Te/Ne correlation t @@ -219,9 +222,9 @@ C 26 foF2 storm model no storm updating t C 27 IG12 from file IG12 - user t C 28 spread-F probability not computed false C 29 IRI01-topside new options as def. by JF(30) false -C 30 IRI01-topside corr. NeQuick topside model false +C 30 IRI01-topside corr. NeQuick topside model false C (29,30) = (t,t) IRIold, (f,t) IRIcor, (f,f) NeQuick -C 31 B0,B1 ABT-2009 B0 Gulyaeva-1987 h0.5 t +C 31 B0,B1 ABT-2009 B0 Gulyaeva-1987 h0.5 t C (4,31) = (t,t) Bil-00, (f,t) ABT-09, (f,f) Gul-87, (t,f) not used C 32 F10.7_81 from file F10.7_81 - user input (oarr(46)) t C 33 Auroral boundary model on/off true/false false @@ -236,10 +239,10 @@ C 41 Use COV=F10.7_365 COV=f(IG12) (IRI before Oct 2015) t C 42 Te with PF10.7 dep. w/o PF10.7 dependance t C 43 B0 from model B0 user input t C .... -C 50 +C 50 C ------------------------------------------------------------------ C -C Depending on the jf() settings additional INPUT parameters may +C Depending on the jf() settings additional INPUT parameters may c be required: C C Setting INPUT parameter @@ -247,11 +250,11 @@ C ----------------------------------------------------------------- C jf(8) =.false. OARR(1)=user input for foF2/MHz or NmF2/m-3 C jf(9) =.false. OARR(2)=user input for hmF2/km or M(3000)F2 C jf(10 )=.false. OARR(15),OARR(16)=user input for Ne(300km), -C Ne(400km)/m-3. Use OARR()=-1 if one of these values is not +C Ne(400km)/m-3. Use OARR()=-1 if one of these values is not C available. If jf(23)=.false. then Ne(300km), Ne(550km)/m-3. -C jf(13) =.false. OARR(3)=user input for foF1/MHz or NmF1/m-3 +C jf(13) =.false. OARR(3)=user input for foF1/MHz or NmF1/m-3 C jf(14) =.false. OARR(4)=user input for hmF1/km -C jf(15) =.false. OARR(5)=user input for foE/MHz or NmE/m-3 +C jf(15) =.false. OARR(5)=user input for foE/MHz or NmE/m-3 C jf(16) =.false. OARR(6)=user input for hmE/km C jf(17) =.flase. OARR(33)=user input for Rz12 C jf(25) =.false. OARR(41)=user input for daily F10.7 index @@ -264,7 +267,7 @@ C OUTF(1,*) ELECTRON DENSITY/M-3 C OUTF(2,*) NEUTRAL TEMPERATURE/K C OUTF(3,*) ION TEMPERATURE/K C OUTF(4,*) ELECTRON TEMPERATURE/K -C OUTF(5,*) O+ ION DENSITY/% or /M-3 if jf(22)=f +C OUTF(5,*) O+ ION DENSITY/% or /M-3 if jf(22)=f C OUTF(6,*) H+ ION DENSITY/% or /M-3 if jf(22)=f C OUTF(7,*) HE+ ION DENSITY/% or /M-3 if jf(22)=f C OUTF(8,*) O2+ ION DENSITY/% or /M-3 if jf(22)=f @@ -272,20 +275,20 @@ C OUTF(9,*) NO+ ION DENSITY/% or /M-3 if jf(22)=f C AND, IF JF(6)=.FALSE.: C OUTF(10,*) CLUSTER IONS DEN/% or /M-3 if jf(22)=f C OUTF(11,*) N+ ION DENSITY/% or /M-3 if jf(22)=f -C OUTF(12,*) -C OUTF(13,*) -C if(jf(24) OUTF(14,1:11) standard IRI-Ne for 60,65,..,110km -C =.false.) 12:22) Friedrich (FIRI) model at these heights -C 23:33) standard Danilov (SW=0, WA=0) -C 34:44) for minor Stratospheric Warming (SW=0.5) -C 45:55) for major Stratospheric Warming (SW=1) +C OUTF(12,*) +C OUTF(13,*) +C if(jf(24) OUTF(14,1:11) standard IRI-Ne for 60,65,..,110km +C =.false.) 12:22) Friedrich (FIRI) model at these heights +C 23:33) standard Danilov (SW=0, WA=0) +C 34:44) for minor Stratospheric Warming (SW=0.5) +C 45:55) for major Stratospheric Warming (SW=1) C 56:66) weak Winter Anomaly (WA=0.5) conditions C 67:77) strong Winter Anomaly (WA=1) conditions C OUTF(15-19,*) free C OUTF(20,*) GEOMAGNETIC FIELD MAGNITUDE C OUTF(21-30,*) NEUTRAL ATMOSPHERE DENSITIES & TEMPERATURES c -C OARR(1:100) ADDITIONAL OUTPUT PARAMETERS +C OARR(1:100) ADDITIONAL OUTPUT PARAMETERS C C #OARR(1) = NMF2/M-3 #OARR(2) = HMF2/KM C #OARR(3) = NMF1/M-3 #OARR(4) = HMF1/KM @@ -306,13 +309,13 @@ C OARR(31) = ISEASON (1=spring) OARR(32) = Geographic longitude C #OARR(33) = Rz12 OARR(34) = Covington Index C OARR(35) = B1 OARR(36) = M(3000)F2 C $OARR(37) = TEC/m-2 $OARR(38) = TEC_top/TEC*100. -C #OARR(39) = gind (IG12) OARR(40) = F1 probability +C #OARR(39) = gind (IG12) OARR(40) = F1 probability C #OARR(41) = F10.7 daily OARR(42) = c1 (F1 shape) -C OARR(43) = daynr OARR(44) = equatorial vertical +C OARR(43) = daynr OARR(44) = equatorial vertical C OARR(45) = foF2_storm/foF2_quiet ion drift in m/s -C #OARR(46) = F10.7_81 OARR(47) = foE_storm/foE_quiet -C OARR(48) = spread-F probability -C OARR(49) = Geomag. latitude OARR(50) = Geomag. longitude +C #OARR(46) = F10.7_81 OARR(47) = foE_storm/foE_quiet +C OARR(48) = spread-F probability +C OARR(49) = Geomag. latitude OARR(50) = Geomag. longitude C OARR(51) = ap at current time OARR(52) = daily ap C OARR(53) = invdip/degree OARR(54) = MLT-Te C OARR(55) = CGM-latitude OARR(56) = CGM-longitude @@ -329,11 +332,11 @@ C OARR(75) = CGM-lati(MLT=16) OARR(76) = CGM-lati for MLT=17 C OARR(77) = CGM-lati(MLT=18) OARR(78) = CGM-lati for MLT=19 C OARR(79) = CGM-lati(MLT=20) OARR(80) = CGM-lati for MLT=21 C OARR(81) = CGM-lati(MLT=22) OARR(82) = CGM-lati for MLT=23 -C OARR(83) = Kp at current time OARR(84) = magnetic declination -C OARR(85) = L-value OARR(86) = dipole moment +C OARR(83) = Kp at current time OARR(84) = magnetic declination +C OARR(85) = L-value OARR(86) = dipole moment C # INPUT as well as OUTPUT parameter C $ special for IRIWeb (only place-holders) -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- C***************************************************************** C*** THE ALTITUDE LIMITS ARE: LOWER (DAY/NIGHT) UPPER *** C*** ELECTRON DENSITY 60/80 KM 1500 KM *** @@ -350,14 +353,23 @@ C********* ALL TIMES ARE IN DECIMAL HOURS ************** C***************************************************************** C***************************************************************** C***************************************************************** - real, intent(inout) :: OARR(100) - real, intent(out) :: OUTF(30,1000) + real, intent(out) :: OARR(100) + real, intent(out) :: OUTF(30,Nalt) + +! leave JF intent off for f2py, due to setting of JF(26). Simple workaround. + logical :: JF(50) + integer, intent(in) :: JMAG,nalt, IYYYY,MMDD + real, intent(in) :: ALATI,ALONG,DHOUR,altkm(nalt) + + real HEIBEG,HEIEND,HEISTP INTEGER DAYNR,DDO,DO2,SEASON,SEADAY - REAL LATI,LONGI,MO2,MO,MODIP,NMF2,MAGBR,INVDIP,IAPO, + REAL LATI,LONGI,MO2,MO,MODIP,NMF2,MAGBR,INVDIP,IAPO, & NMF1,NME,NMD,MM,MLAT,MLONG,NMF2S,NMES,INVDPC - CHARACTER FILNAM*12 + CHARACTER(12) FILNAM + character(*), intent(in) :: datadir character(256) dirdata1,filename + integer iuccir, i c-web-for webversion c CHARACTER FILNAM*53 @@ -373,21 +385,21 @@ c CHARACTER FILNAM*53 LOGICAL EXT,SCHALT,TECON(2),sam_mon,sam_yea,sam_ut,sam_date, & F1REG,FOF2IN,HMF2IN,URSIF2,LAYVER,DY,DREG,rzino,FOF1IN, & HMF1IN,FOEIN,HMEIN,RZIN,sam_doy,F1_OCPRO,F1_L_COND,NODEN, - & NOTEM,NOION,TENEOP,OLD79,JF(50),URSIFO,igin,igino,mess, + & NOTEM,NOION,TENEOP,OLD79,URSIFO,igin,igino,mess, & dnight,enight,fnight,fstorm_on,estorm_on,B0IN COMMON /CONST/UMR,PI /const1/humr,dumr /ARGEXP/ARGMAX c & /const2/icalls,montho,nmono,iyearo,idaynro,ursifo,rzino, c & igino,ut0 & /IGRF1/ERA,AQUAD,BQUAD,DIMO - & /BLOCK1/HMF2,NMF2S,HMF1,F1REG /BLOCK2/B0,B1,C1 - & /BLOCK3/HZ,T,HST /BLOCK4/HME,NMES,HEF + & /BLOCK1/HMF2,NMF2S,HMF1,F1REG /BLOCK2/B0,B1,C1 + & /BLOCK3/HZ,T,HST /BLOCK4/HME,NMES,HEF & /BLOCK5/ENIGHT,E /BLOCK6/HMD,NMD,HDX & /BLOCK7/D1,XKK,FP30,FP3U,FP1,FP2 - & /BLOCK8/HS,TNHS,XSM,MM,DTI,MXSM + & /BLOCK8/HS,TNHS,XSM,MM,DTI,MXSM & /BLOTE/AHH,ATE1,STTE,DTE - & /BLO10/BETA,ETA,DELTA,ZETA /findRLAT/FLON,RYEAR - & /BLO11/B2TOP,TC3,itopn,alg10,hcor1 + & /BLO10/BETA,ETA,DELTA,ZETA /findRLAT/FLON,RYEAR + & /BLO11/B2TOP,TC3,itopn,alg10,hcor1 & /iounit/konsol,mess /CSW/SW(25),ISW,SWC(25) & /QTOP/Y05,H05TOP,QF,XNETOP,XM3000,HHALF,TAU common /folders/ dirdata1 @@ -396,21 +408,25 @@ c & igino,ut0 DATA icalls/0/ save - + + dirdata1 = datadir + call read_ig_rz + call readapf107 + mess=jf(34) - -c set switches for NRLMSIS00 + +c set switches for NRLMSIS00 ISW=0 do 6492 KI=1,25 6492 SWMI(KI)=1. - nummax=1000 + nummax=Nalt DO 7397 KI=1,30 do 7397 kk=1,nummax 7397 OUTF(KI,kk)=-1. C C oarr(1:6,15,16,33,39:41) is used for inputs -C +C oarr(7)=-1. oarr(8)=-1. oarr(9)=-1. @@ -427,44 +443,53 @@ C C C PROGRAM CONSTANTS AND INITIALIZATION C - if(icalls.lt.1) then - ARGMAX=88.0 - pi=ATAN(1.0)*4. - UMR=pi/180. - humr=pi/12. - dumr=pi/182.5 - ALOG2=ALOG(2.) - ALG10=ALOG(10.) - ALG100=ALOG(100.) - montho=-1 - nmono=-1 - iyearo=-1 - idaynro=-1 - rzino=.true. - igino=.true. - ut0=-1 - ursifo=.true. + if(icalls < 1) then + ARGMAX=88.0 + pi=ATAN(1.0)*4. + UMR=pi/180. + humr=pi/12. + dumr=pi/182.5 + + ALG10=LOG(10.) + + montho=-1 + nmono=-1 + iyearo=-1 + idaynro=-1 + rzino=.true. + igino=.true. + ut0=-1 + ursifo=.true. C Initialize parameters for COMMON/IGRF1/ -C ERA EARTH RADIUS (WGS-84: 6371.137 KM) +C ERA EARTH RADIUS (WGS-84: 6371.137 KM) C EREQU MAJOR HALF AXIS FOR EARTH ELLIPSOID (6378.160 KM) C ERPOL MINOR HALF AXIS FOR EARTH ELLIPSOID (6356.775 KM) C AQUAD SQUARE OF MAJOR HALF AXIS FOR EARTH ELLIPSOID C BQUAD SQUARE OF MINOR HALF AXIS FOR EARTH ELLIPSOID C EEXC Eccentricity of Earth's orbit -C DIMO Earth's dipole moment in Gauss -C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL +C DIMO Earth's dipole moment in Gauss +C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL C ASTRONOMICAL UNION . - ERA=6371.2 - EREQU=6378.16 - ERPOL=6356.775 - AQUAD=EREQU*EREQU - BQUAD=ERPOL*ERPOL - EEXC=0.01675 - dimo=0.311653 - endif - - numhei=int(abs(heiend-heibeg)/abs(heistp))+1 - if(numhei.gt.nummax) numhei=nummax + ERA=6371.2 + EREQU=6378.16 + ERPOL=6356.775 + AQUAD=EREQU*EREQU + BQUAD=ERPOL*ERPOL + EEXC=0.01675 + dimo=0.311653 + endif + + !numhei=int(abs(heiend-heibeg)/abs(heistp))+1 + heibeg = altkm(1) + heiend = altkm(Nalt) + heistp = altkm(2) - altkm(1) ! assumes uniform alt grid + if (.not.all(altkm(2:)-altkm(1:Nalt-1) == altkm(2)-altkm(1))) then + error stop 'uniform altitude grid required' + endif + + + numhei = Nalt + !if(numhei > nummax) numhei=nummax C C NEW-GUL------------------------------ c Y05=.6931473 @@ -506,15 +531,15 @@ C FIRST SPECIFY YOUR COMPUTERS CHANNEL NUMBERS .................... C AGNR=OUTPUT (OUTPUT IS DISPLAYED OR STORED IN FILE OUTPUT.IRI)... C IUCCIR=UNIT NUMBER FOR CCIR COEFFICIENTS ........................ c - IUCCIR=10 + !IUCCIR=10 c-web- special for web version -c-web- messages should be turned off with mess=jf(34)=.false. +c-web- messages should be turned off with mess=jf(34)=.false. - KONSOL=6 if(.not.jf(12).and.mess) then - konsol=11 - open(11,file='messages.txt') - endif + open(newunit=konsol,file='messages.txt') + else + KONSOL = error_unit + endif c c selection of density, temperature and ion composition options ...... c @@ -538,7 +563,7 @@ c else oarr(33)=-1. ENDIF - + IGIN=(.not.jf(27)) IF(IGIN) THEN AIGIN=OARR(39) @@ -561,12 +586,12 @@ c c Topside density .................................................... c if(jf(29)) then -c if (jf(30)) then +c if (jf(30)) then itopn=0 c else c itopn=3 ! Gulyaeva topside option c endif - else + else if (jf(30)) then itopn=1 else @@ -659,8 +684,8 @@ C IF(TENEOP) THEN DO 8154 JXNAR=1,2 XNAR(JXNAR)=OARR(JXNAR+14) - TECON(JXNAR)=.FALSE. -8154 IF(XNAR(JXNAR).GT.0.) TECON(JXNAR)=.TRUE. + TECON(JXNAR)=.FALSE. +8154 IF(XNAR(JXNAR).GT.0.) TECON(JXNAR)=.TRUE. else oarr(15)=-1. oarr(16)=-1. @@ -670,65 +695,65 @@ c lists the selected options before starting the table c if(icalls.ge.1.or.(.not.mess)) goto 8201 - write(konsol,2911) + write(konsol,2911) if(NODEN) goto 2889 - if(LAYVER) write(konsol,9012) - if(OLD79) write(konsol,9014) + if(LAYVER) write(konsol,9012) + if(OLD79) write(konsol,9014) if (itopn.eq.0) write(konsol,9207) if (itopn.eq.1) write(konsol,9204) if (itopn.eq.2) write(konsol,9205) c if (itopn.eq.3) write(konsol,9206) if(FOF2IN) then - write(konsol,9015) + write(konsol,9015) goto 2889 endif if(URSIF2) then - write(konsol,9016) + write(konsol,9016) else - write(konsol,9017) + write(konsol,9017) endif - if(HMF2IN) write(konsol,9018) - if(fof1in) write(konsol,9019) - if(HMF1IN.and.LAYVER) write(konsol,9021) + if(HMF2IN) write(konsol,9018) + if(fof1in) write(konsol,9019) + if(HMF1IN.and.LAYVER) write(konsol,9021) if(jf(4)) then write(konsol,9214) - else + else if (jf(31)) then write(konsol,9216) else write(konsol,9215) endif endif - if(foein) write(konsol,9022) - if(HMEIN) write(konsol,9023) - if(B0IN) write(konsol,9923) - if(F1_OCPRO) write(konsol,9024) - if(F1_L_COND) write(konsol,9025) - if(DREG) then - write(konsol,9026) + if(foein) write(konsol,9022) + if(HMEIN) write(konsol,9023) + if(B0IN) write(konsol,9923) + if(F1_OCPRO) write(konsol,9024) + if(F1_L_COND) write(konsol,9025) + if(DREG) then + write(konsol,9026) else - write(konsol,9027) + write(konsol,9027) endif if(jf(26)) then - if(fof2in) then - write(konsol,9028) + if(fof2in) then + write(konsol,9028) jf(26)=.false. else - write(konsol,9029) - endif + write(konsol,9029) endif + endif 2889 continue - if((.not.NOION).and.(DY)) write(konsol,9031) - if((.not.NOION).and.(.not.DY)) write(konsol,9039) + if((.not.NOION).and.(DY)) write(konsol,9031) + if((.not.NOION).and.(.not.DY)) write(konsol,9039) if(NOTEM) goto 8201 - if(TENEOP) write(konsol,9032) - if(jf(23)) then - write(konsol,9033) + if(TENEOP) write(konsol,9032) + if(jf(23)) then + write(konsol,9033) else - write(konsol,9034) + write(konsol,9034) endif if(jf(33)) then @@ -784,14 +809,14 @@ c if (itopn.eq.3) write(konsol,9206) 8201 continue C -C CALCULATION OF DAY OF YEAR OR MONTH/DAY AND DECIMAL YEAR -c NRDAYM is the number of days in the current month +C CALCULATION OF DAY OF YEAR OR MONTH/DAY AND DECIMAL YEAR +c NRDAYM is the number of days in the current month c IDAYY is the number of days in the current year c c leap year rule: years evenly divisible by 4 are leap years, except -c years also evenly divisible by 100 are not leap years, except years -c also evenly divisible by 400 are leap years. The year 2000 is a 100 -c and 400 year exception and therefore it is a normal leap year. +c years also evenly divisible by 100 are not leap years, except years +c also evenly divisible by 400 are leap years. The year 2000 is a 100 +c and 400 year exception and therefore it is a normal leap year. c The next 100 year exception will be in the year 2100! c @@ -820,25 +845,24 @@ C calculate center height for CGM computation C height_center=(HEIBEG+HEIEND)/2. - + C -C CALCULATION OF GEODETIC/GEOMAGNETIC COORDINATES (LATI, LONGI AND -C MLAT, MLONG), MAGNETIC INCLINATION (DIP), DIP LATITUDE (MAGBR) +C CALCULATION OF GEODETIC/GEOMAGNETIC COORDINATES (LATI, LONGI AND +C MLAT, MLONG), MAGNETIC INCLINATION (DIP), DIP LATITUDE (MAGBR) C AND MODIFIED DIP (MODIP), ALL IN DEGREES C - if(along.lt.0.) along = along + 360. ! -180/180 to 0-360 - IF(JMAG.GT.0) THEN MLAT=ALATI - MLONG=ALONG + MLONG= modulo(ALONG, 360.) ! -180/180 to 0-360 ELSE LATI=ALATI - LONGI=ALONG + LONGI= modulo(ALONG, 360.) ENDIF CALL GEODIP(IYEAR,LATI,LONGI,MLAT,MLONG,JMAG) +! print *,'data directory ' ,dirdata1 CALL FELDCOF(RYEAR) if(jf(18)) then @@ -866,53 +890,53 @@ c c C CALCULATION OF UT/LT and XMLT ............... c - IF(DHOUR.le.24.0) then + IF(DHOUR <= 24.) then HOUR=DHOUR ! dhour =< 24 is LT hourut=hour-longi/15. - if(hourut.lt.0) hourut=hourut+24. + if(hourut < 0) hourut=hourut+24. else hourut=DHOUR-25. ! dhour>24 is UT+25 hour=hourut+longi/15. ! hour becomes LT - if(hour.gt.24.) hour=hour-24. + if(hour > 24.) hour=hour-24. endif - + CALL CLCMLT(IYEAR,DAYNR,HOURUT,LATI,LONGI,XMLT) c -c SEASON assumes equal length seasons (92 days) with spring -c (SEASON=1) starting at day-of-year=45; for lati < 0 adjustment +c SEASON assumes equal length seasons (92 days) with spring +c (SEASON=1) starting at day-of-year=45; for lati < 0 adjustment c for southern hemisphere is made. Some models require the c seasonal month (ISEAMON) or the seasonal day-of year (SEADAY) c ZMONTH is decimal month (Jan 1 = 1.0 and Dec 31 = 12.97) c SDAY is the day number reduced to a 360 day year (TOPH05) -c NRDAYM is the number of days in the current month +c NRDAYM is the number of days in the current month c IDAYY is the number of days in the current year -c - +c + SEASON=INT((DAYNR+45.0)/92.0) IF(SEASON.LT.1) SEASON=4 NSEASN=SEASON ! Northern hemisphere season zmonth = month + (iday-1)*1./nrdaym C NEW-GUL------------------------------ - sday=daynr/idayy*360. + sday=daynr/idayy*360. C NEW-GUL------------------------------ seaday=daynr iseamon=month IF(LATI.GE.0.0) GOTO 5592 - SEASON=SEASON-2 + SEASON=SEASON-2 IF(SEASON.LT.1) SEASON=SEASON+4 iseamon=month+6 if(iseamon.gt.12) iseamon=iseamon-12 seaday=daynr+idayy/2. if(seaday.gt.idayy) seaday=seaday-idayy C NEW-GUL------------------------------ - sday=sday+180. - if (sday.gt.360.) sday=sday-360. + sday=sday+180. + if (sday.gt.360.) sday=sday-360. C NEW-GUL------------------------------ C -C 12-month running mean sunspot number (rssn) and Ionospheric Global -C index (gind), daily F10.7 cm solar radio flux (f107d) and monthly -C F10.7 (cov) index +C 12-month running mean sunspot number (rssn) and Ionospheric Global +C index (gind), daily F10.7 cm solar radio flux (f107d) and monthly +C F10.7 (cov) index C 5592 continue @@ -921,12 +945,12 @@ C sam_doy=(daynr.eq.idaynro) sam_date=(sam_yea.and.sam_doy) sam_ut=(hourut.eq.ut0) - + if(sam_date.and..not.rzino.and..not.rzin. & and..not.igin.and..not.igino) goto 2910 - + call tcon(iyear,month,iday,daynr,rzar,arig,ttt,nmonth) - if(nmonth.lt.0) goto 3330 ! jump to end of program + if(nmonth.lt.0) error stop 'invalid month' ! jump to end of program if(RZIN) then rrr = arzin @@ -951,11 +975,11 @@ c arig(3) = zi c rlimit=gind c COVSAT=63.75+rlimit*(0.728+rlimit*0.00089) -C Getting F10.7 index: daily (f107d), previous day (f107y; +C Getting F10.7 index: daily (f107d), previous day (f107y; C required by MSIS), 81-day average (f10781), 365-day average -C (f107365), and PF10.7=(F10.7_daily + F10.7_81_day)/2. -C F10.7 should be adjusted (to top of atmosphere) value not -C observed (at the ground) value. +C (f107365), and PF10.7=(F10.7_daily + F10.7_81_day)/2. +C F10.7 should be adjusted (to top of atmosphere) value not +C observed (at the ground) value. f107d=cov f107y=cov f10781=cov @@ -971,28 +995,28 @@ C observed (at the ground) value. endif endif if(.not.jf(25)) then - f107d=f107din ! user input: F10.7 daily - f107y=f107din ! user input: F10.7 previous day + f107d=f107din ! user input: F10.7 daily + f107y=f107din ! user input: F10.7 previous day endif if(.not.jf(32)) then f10781=f10781in ! F10.7_81 user input - endif + endif pf107 = (f107d+f10781)/2. c Correcting F10.7 adjusted flux from APF107.DAT to flux observed at -c Earth that is expected by NRLMSIS00 (GTD7) and CHEMION +c Earth that is expected by NRLMSIS00 (GTD7) and CHEMION - f_adj=radj*radj + f_adj=radj*radj f107yo=f107y/f_adj f10781o=f10781/f_adj - - if(jf(41)) cov=f107365 + + if(jf(41)) cov=f107365 COVSAT=cov if(covsat.gt.188.) covsat=188 C -C CALCULATION OF SOLAR ZENITH ANGLE (XHI/DEG), SUN DECLINATION ANGLE -C (SUNDEC),SOLAR ZENITH ANGLE AT NOON (XHINON) AND TIME OF LOCAL +C CALCULATION OF SOLAR ZENITH ANGLE (XHI/DEG), SUN DECLINATION ANGLE +C (SUNDEC),SOLAR ZENITH ANGLE AT NOON (XHINON) AND TIME OF LOCAL C SUNRISE/SUNSET (SAX, SUX; dec. hours) AT 70 KM (D-REGION), 110 KM C (E-REGION), 200 KM (F1-REGION), AND 500 KM (TE, TI). C @@ -1037,7 +1061,7 @@ C C CALCULATION OF ELECTRON DENSITY PARAMETERS................ C lower height boundary (HNEA), upper boundary (HNEE) C - + 1334 continue HNEA=65. IF(DNIGHT) HNEA=80. @@ -1065,7 +1089,7 @@ c c c F2 peak critical frequency foF2, density NmF2, and height hmF2 c -C READ CCIR AND URSI COEFFICIENT SET FOR CHOSEN MONTH +C READ CCIR AND URSI COEFFICIENT SET FOR CHOSEN MONTH C IF((FOF2IN).AND.(HMF2IN).and.(itopn.ne.2)) GOTO 501 IF(URSIF2.NEQV.URSIFO) GOTO 7797 @@ -1079,8 +1103,8 @@ C 104 FORMAT('ccir',I2,'.asc') c-web-for webversion c104 FORMAT('/var/www/omniweb/cgi/vitmo/IRI/ccir',I2,'.asc') - filename=trim(trim(trim(dirdata1) // '/ccir/') // trim(FILNAM)) - OPEN(IUCCIR,FILE=filename,STATUS='OLD',ERR=8448, + filename=trim(trim(dirdata1) // '/ccir/' // FILNAM) + OPEN(newunit=IUCCIR,FILE=filename,STATUS='OLD',ERR=8448, & FORM='FORMATTED') READ(IUCCIR,4689) F2,FM3 4689 FORMAT(1X,4E15.8) @@ -1093,16 +1117,16 @@ C 1144 FORMAT('ursi',I2,'.asc') c-web-for webversion c1144 FORMAT('/var/www/omniweb/cgi/vitmo/IRI/ursi',I2,'.asc') - filename=trim(trim(trim(dirdata1)//'/ursi/')//trim(FILNAM)) - OPEN(IUCCIR,FILE=filename,STATUS='OLD',ERR=8448, + filename=trim(trim(dirdata1)//'/ursi/'// FILNAM) + OPEN(newunit=IUCCIR,FILE=filename,STATUS='OLD',ERR=8448, & FORM='FORMATTED') READ(IUCCIR,4689) F2 CLOSE(IUCCIR) endif C -C READ CCIR AND URSI COEFFICIENT SET FOR NMONTH, i.e. previous -c month if day is less than 15 and following month otherwise +C READ CCIR AND URSI COEFFICIENT SET FOR NMONTH, i.e. previous +c month if day is less than 15 and following month otherwise C 4293 continue @@ -1112,8 +1136,8 @@ c first CCIR .............................................. c WRITE(FILNAM,104) NMONTH+10 - filename=trim(trim(trim(dirdata1) // '/ccir/') // trim(FILNAM)) - OPEN(IUCCIR,FILE=filename,STATUS='OLD',ERR=8448, + filename=trim(trim(dirdata1) // '/ccir/'// FILNAM) + OPEN(newunit=IUCCIR,FILE=filename,STATUS='OLD',ERR=8448, & FORM='FORMATTED') READ(IUCCIR,4689) F2N,FM3N CLOSE(IUCCIR) @@ -1123,19 +1147,19 @@ C then URSI if chosen ..................................... C if(URSIF2) then WRITE(FILNAM,1144) NMONTH+10 - filename=trim(trim(trim(dirdata1)//'/ursi/')//trim(FILNAM)) - OPEN(IUCCIR,FILE=filename,STATUS='OLD',ERR=8448, + filename=trim(trim(dirdata1)//'/ursi/'// FILNAM) + OPEN(newunit=IUCCIR,FILE=filename,STATUS='OLD',ERR=8448, & FORM='FORMATTED') READ(IUCCIR,4689) F2N CLOSE(IUCCIR) endif GOTO 4291 - + 8448 WRITE(konsol,8449) filename 8449 FORMAT(1X////, & ' The file ',A30,'is not in your directory.') - GOTO 3330 + error stop C C LINEAR INTERPOLATION IN SOLAR ACTIVITY. IG12 used for foF2 C @@ -1174,7 +1198,7 @@ C yfof2 = zfof2 + ttt * (fof2n-zfof2) xm3000= zm3000+ ttt * (xm300n-zm3000) endif - + 501 IF(FOF2IN) THEN FOF2=AFOF2 NMF2=ANMF2 @@ -1195,14 +1219,14 @@ c fstorm_on=jf(26).and.jf(8) estorm_on=jf(35).and.jf(15) if(fstorm_on.or.jf(33).or.estorm_on) then -c if(.not.sam_date.or..not.sam_ut) then +c if(.not.sam_date.or..not.sam_ut) then call apf(isdate,hourut,indap) endif - + c -c stormtime updating for foF2 (foF2s, NmF2s) +c stormtime updating for foF2 (foF2s, NmF2s) c - if(fstorm_on.and.(indap(1).gt.-1)) then + if(fstorm_on.and.(indap(1).gt.-1)) then icoord=1 kut=int(hourut) call STORM(indap,lati,longi,icoord,cglat,kut, @@ -1213,7 +1237,7 @@ c c c stormtime updating for foE (foEs, NmEs) c - if(estorm_on.and.(indap(1).gt.-1)) then + if(estorm_on.and.(indap(1).gt.-1)) then estormcor=STORME_AP(DAYNR,MLAT,INDAP(13)*1.0) if(estormcor.gt.-2.0) foes=foe*estormcor NMES=1.24E10*FOES*FOES @@ -1224,10 +1248,10 @@ c if(jf(33)) then if(indap(1).gt.-1) then xkp=ckp(indap(13)) - else + else xkp=3.0 - endif -c Corrected magnetic latitude CGM of equatorward boundary, + endif +c Corrected magnetic latitude CGM of equatorward boundary, c ab_mlat(48), for MLT=0.0,0.5,1.0 ... 23.5 h and kp=xkp call auroral_boundary(xkp,-1.0,cgmlat,ab_mlat) DAT(1,1)=lati @@ -1243,7 +1267,7 @@ c cgm_mlt00_ut=DAT(11,1) if(cgm_mlt.lt.0.) cgm_mlt=24.+hourut-cgm_mlt00_ut c print*,cgm_mlt c cgm_mlt_ut=DAT(11,1) -c cgm_mlt=cgm_mlt_ut+cgm_lon/15. +c cgm_mlt=cgm_mlt_ut+cgm_lon/15. c if(cgm_mlt.gt.24.) cgm_mlt=cgm_mlt-24. c CGM latitude of boundary (cgmlat) for present MLT value @@ -1251,7 +1275,7 @@ C 2012.02 12/17/12 Add magnetic declination as oarr(84) output zmlt=xmlt C zmlt=cgm_mlt cgmlat=100.0 - if(zmlt.ge.0.0.and.zmlt.le.24.0) + if(zmlt.ge.0.0.and.zmlt.le.24.0) & call auroral_boundary(xkp,zmlt,cgmlat,ab_mlat) endif IF(((.not.JF(4)).and.JF(31)).or.((.not.JF(39)).and.JF(40))) THEN @@ -1269,22 +1293,22 @@ C zmlt=cgm_mlt endif RLAT=XRLAT ENDIF - + IF(HMF2IN) THEN IF(AHMF2.LT.50.0) THEN XM3000=AHMF2 HMF2=HMF2ED(MAGBR,RSSN,FOF2/FOE,XM3000) - ELSE + ELSE HMF2=AHMF2 c XM3000=XM3000HM(MAGBR,RSSN,FOF2/FOE,HMF2) ENDIF ELSE IF(JF(39)) THEN ratf=fof2/foe -c set jf(36)=false to use foF2_storm in hmF2 formula +c set jf(36)=false to use foF2_storm in hmF2 formula if(.not.jf(36)) ratf=fof2s/foe HMF2=HMF2ED(MAGBR,RSSN,RATF,XM3000) ELSE IF(JF(40)) THEN -c AMTB digisonde model +c AMTB digisonde model CALL SHAMDHMF2(RLAT,FLON,ZMONTH,RSSN,HMF2) ELSE c SHUBIN-COSMIC model @@ -1316,10 +1340,10 @@ c EX1=EX+1 EPIN=4.*EX/(EX1*EX1) ETA1=-0.02*EPIN - ETA = 0.058798 + ETA1 - - & FLU * (0.014065 - 0.0069724 * COS2) + + ETA = 0.058798 + ETA1 - + & FLU * (0.014065 - 0.0069724 * COS2) + & FO1* (0.0024287 + 0.0042810 * COS2 - 0.0001528 * FO1) - ZETA = 0.078922 - 0.0046702 * COS2 - + ZETA = 0.078922 - 0.0046702 * COS2 - & FLU * (0.019132 - 0.0076545 * COS2) + & FO1* (0.0032513 + 0.0060290 * COS2 - 0.00020872 * FO1) BETA=-128.03 + 20.253 * COS2 - @@ -1339,10 +1363,10 @@ c zmp111 = zmp1 / (zmp11 * zmp11) zmp2 = exp(modip / 19.) zmp22 = 1. + zmp2 - zmp222 = zmp2 / (zmp22 * zmp22) - r2n = -0.84 - 1.6 * zmp111 - r2d = -0.84 - 0.64 * zmp111 - x1n = 230. - 700. * zmp222 + zmp222 = zmp2 / (zmp22 * zmp22) + r2n = -0.84 - 1.6 * zmp111 + r2d = -0.84 - 0.64 * zmp111 + x1n = 230. - 700. * zmp222 x1d = 550. - 1900. * zmp222 r2 = HPOL(HOUR,r2d,r2n,SAX300,SUX300,1.,1.) x1 = HPOL(HOUR,x1d,x1n,SAX300,SUX300,1.,1.) @@ -1379,7 +1403,7 @@ c endif c c Bottomside thickness parameter B0 and shape parameters B1 -c +c if(jf(4)) then B0=B0_98(HOUR,SAX200,SUX200,NSEASN,RSSN,LONGI,MODIP) B1=HPOL(HOUR,1.9,2.6,SAX200,SUX200,1.,1.) @@ -1410,7 +1434,7 @@ c F1 layer thickness parameter c1 c c1 = f1_c1(modip,hour,sux200,sax200) c -c F1 occurrence probability with Scotto et al. 1997 or Ducharme et al. +c F1 occurrence probability with Scotto et al. 1997 or Ducharme et al. c if jf(19)=f1_ocpro=.true. or .false. c If .not.jf(20)=f1_l_cond=.true. then Scotto model with L-condition c @@ -1420,17 +1444,17 @@ c if(f1_l_cond) f1pb = f1pbl f1reg=.false. if((fof1in).or.(f1pb.ge.0.5)) f1reg=.true. - else - f1pb = 0.0 - if(fof1in.or.((.not.fnight).and.(fof1.gt.0.0))) f1pb=1. + else + f1pb = 0.0 + if(fof1in.or.((.not.fnight).and.(fof1.gt.0.0))) f1pb=1. f1reg=.false. if(f1pb.gt.0.0) f1reg=.true. endif - + c -c E-valley: DEPTH=(NmE-N_deepest)/NmE*100, WIDTH=HEF-HmE, -c distance of deepest value point above E-peak(HDEEP), -c derivative at valley top divided by NmE (DLNDH), +c E-valley: DEPTH=(NmE-N_deepest)/NmE*100, WIDTH=HEF-HmE, +c distance of deepest value point above E-peak(HDEEP), +c derivative at valley top divided by NmE (DLNDH), c and height of valley top (HEF) c XDEL=XDELS(SEASON)/DELA @@ -1474,16 +1498,16 @@ c XDX=NMD*EXP(X*(FP1+X*(FP2+X*FP30))) DXDX=XDX*(FP1+X*(2.0*FP2+X*3.0*FP30)) X=HME-HDX - XKK=-DXDX*X/(XDX*ALOG(XDX/NMES)) + XKK=-DXDX*X/(XDX*LOG(XDX/NMES)) c -c if exponent xkk is larger than xkkmax, then xkk will be set to -c xkkmax and d1 will be determined such that the point hdx/xdx is +c if exponent xkk is larger than xkkmax, then xkk will be set to +c xkkmax and d1 will be determined such that the point hdx/xdx is c reached; derivative is no longer continuous. c xkkmax=5. if(xkk.gt.xkkmax) then xkk=xkkmax - d1=-alog(xdx/nmes)/(x**xkk) + d1=-log(xdx/nmes)/(x**xkk) else D1=DXDX/(XDX*XKK*X**(XKK-1.0)) endif @@ -1519,7 +1543,7 @@ c do ii=1,11 ddens(4,ii)=-1. if(ii.lt.8) ddens(4,ii)=10**(elg(ii)+6) - enddo + enddo f5sw=0. f6wa=1. call DRegion(xhi,month,f107d,vKp,f5SW,f6WA,elg) @@ -1550,7 +1574,7 @@ c omit F1 feature if nmf1*0.9 is smaller than nme goto 9427 endif goto 9245 - endif + endif CALL REGFA1(HEF,HMF2,XE2H,NMF2S,0.001,NMF1,XE2,SCHALT,HMF1) IF(.not.SCHALT) GOTO 3801 @@ -1558,7 +1582,7 @@ c c omit F1 feature .................................................... c -9427 if(mess) WRITE(KONSOL,11) +9427 if(mess) WRITE(KONSOL,11) 11 FORMAT(1X,'*NE* HMF1 IS NOT EVALUATED BY THE FUNCTION XE2'/ & 1X,'CORR.: NO F1 REGION, B1=3, C1=0.0') HMF1=0. @@ -1599,7 +1623,7 @@ C hf2=hef xf2=xe3_1(hf2) if(xf2.gt.nmes) goto 3885 - + CALL REGFA1(hf1,HF2,XF1,XF2,0.001,NMES,XE3_1,SCHALT,HST) if(schalt) goto 3885 @@ -1657,7 +1681,7 @@ C IAPO(1)=0. else SWMI(9)=-1.0 - endif + endif CALL TSELEC(SWMI) CALL GTD7(IYD,SEC,HEQUI,LATI,LONGI,HOUR,F10781o,F107Yo,IAPO,0, & D_MSIS,T_MSIS) @@ -1673,7 +1697,7 @@ C secni=utni*3600. C endif CALL GTD7(IYD,SECNI,HEQUI,LATI,LONGI,0.0,F10781o,F107Yo,IAPO,0, & D_MSIS,T_MSIS) - TN1NI=T_MSIS(2) + TN1NI=T_MSIS(2) ELSE TN1NI=T_MSIS(2) ENDIF @@ -1689,7 +1713,7 @@ c Te(120km) = Tn(120km) AHH(1)=120. ATE(1)=TN120 -C Te-MAXIMUM based on JICAMARCA and ARECIBO data +C Te-MAXIMUM based on JICAMARCA and ARECIBO data HMAXD=60.*EXP(-(MLAT/22.41)**2)+210. HMAXN=150. @@ -1700,13 +1724,13 @@ C Te-MAXIMUM based on JICAMARCA and ARECIBO data TMAXN=T_MSIS(2) ATE(2)=HPOL(HOUR,TMAXD,TMAXN,SAX200,SUX200,1.,1.) -c Te(300km), Te(400km) from AE-C, Te(1400km), Te(3000km) from +c Te(300km), Te(400km) from AE-C, Te(1400km), Te(3000km) from c ISIS, Brace and Theis DIPLAT=MAGBR CALL TEBA(DIPLAT,HOUR,NSEASN,TEA) - icd=0 + icd=0 if(jf(23)) then c Te at fixed heights taken from Brace and Theis @@ -1747,7 +1771,7 @@ c ise=0 ! season correction off do ijk=3,7 c call igrf_sub(lati,longi,ryear,ahh(ijk), c & xl,icode,dipl,babs) -c if(xl.gt.10.) xl=10. +c if(xl.gt.10.) xl=10. c call elteik(1,isa,invdip,xl6,dimo,babs6,dipl6, c & xmlt,ahh(ijk),daynr,pf107,teh2,sdte) call elteik(isa,invdip,xmlt,ahh(ijk),daynr,pf107, @@ -1776,7 +1800,7 @@ c Te corrected and Te > Tn enforced TNAHHI=T_MSIS(2) IF(ATE(I+1).LT.TNAHHI) ATE(I+1)=TNAHHI STTE2=(ATE(I+1)-ATE(I))/(AHH(I+1)-AHH(I)) - ATE(I)=ATE(I)-(STTE2-STTE1)*DTE(I-1)*ALOG2 + ATE(I)=ATE(I)-(STTE2-STTE1)*DTE(I-1)*LOG(2.) 1901 STTE1=STTE2 c Te gradients STTE are computed for each segment @@ -1808,7 +1832,7 @@ c Ti(430km) duirng nighttime from AEROS data c Ti(430km) for specified time using HPOL - TI1=TIN1 + TI1=TIN1 IF(TID1.GT.TIN1) TI1=HPOL(HOUR,TID1,TIN1,SAX300,SUX300,1.,1.) c Tn < Ti < Te enforced @@ -1844,7 +1868,7 @@ c XTETI is altitude where Te=Ti IF(XTTS.GT.0.1) GOTO 2390 XTETI=X+XTTS*5. -c Ti=Te above XTETI +c Ti=Te above XTETI MXSM=3 MM(3)=STTE(6) @@ -1884,7 +1908,7 @@ C kk=1 300 IF(NODEN) GOTO 330 - IF((HEIGHT.GT.HNEE).OR.(HEIGHT.LT.HNEA)) GOTO 330 + IF((HEIGHT > HNEE).OR.(HEIGHT < HNEA)) GOTO 330 IF(LAYVER) THEN ELEDE=-9. IF(IIQU.LT.2) ELEDE= @@ -1905,7 +1929,7 @@ c plasma temperatures c 330 IF(NOTEM) GOTO 7108 - IF((HEIGHT.GT.HTE).OR.(HEIGHT.LT.HTA)) GOTO 7108 + IF((HEIGHT > HTE).OR.(HEIGHT < HTA)) GOTO 7108 CALL GTD7(IYD,SEC,HEIGHT,LATI,LONGI,HOUR,F10781o,F107Yo, & IAPO,48,D_MSIS,T_MSIS) TNH=T_MSIS(2) @@ -1915,7 +1939,7 @@ c if(TIH.lt.TNH) TIH=TNH endif TEH=TNH - if(HEIGHT.GT.HEQUI) then + if(HEIGHT.GT.HEQUI) then TEH=ELTE(HEIGHT) if(TEH.lt.TIH) TEH=TIH endif @@ -1931,7 +1955,7 @@ C C C NEUTRAL DENSITIES (cm-3) AS OUTPUT -C +C OUTF(21,kk)=D_MSIS(1) ! He OUTF(22,kk)=D_MSIS(2) ! O OUTF(23,kk)=D_MSIS(3) ! N2 @@ -1983,7 +2007,7 @@ c Richards-Bilitza-Voglozin-2010 IDC model CALL CHEMION(jprint,height,F107Yo,F10781o,TEH,TIH,TNH, & D_MSIS(2),D_MSIS(4),D_MSIS(3),D_MSIS(1),D_MSIS(7), & -1.0,XN4S,EDENS,-1.0,xhi,ro,ro2,rno,rn2,rn,Den_NO, - & Den_N2D,INEWT) + & Den_N2D,INEWT) if(INEWT.gt.0) then sumion = edens/100. rox=ro/sumion @@ -2011,7 +2035,7 @@ c c ion densities are given in percent of total electron density; c - if(jf(22)) then + if(jf(22)) then xnorm=1 else xnorm=elede/100. @@ -2040,32 +2064,32 @@ c if(kk.le.numhei) goto 300 C -C END OF PARAMETER COMPUTATION LOOP +C END OF PARAMETER COMPUTATION LOOP C c c D region special: densities for 11 heights (60,65,70,..,110km) -c outf(14,1:11)=IRI-07, outf(14,12:22)=FIRI, -c outf(14,23:33)= Danilov et al.(1995) with SW=0,WA=0 -c outf(14,34:44)= with SW=0.5,WA=0, -c outf(14,45:55)= with SW=1,WA=0, -c outf(14,56:66)= with SW=0,WA=0.5, -c outf(14,67:77)= with SW=0,WA=1, +c outf(14,1:11)=IRI-07, outf(14,12:22)=FIRI, +c outf(14,23:33)= Danilov et al.(1995) with SW=0,WA=0 +c outf(14,34:44)= with SW=0.5,WA=0, +c outf(14,45:55)= with SW=1,WA=0, +c outf(14,56:66)= with SW=0,WA=0.5, +c outf(14,67:77)= with SW=0,WA=1, c if(.not.dreg) then do ii=1,11 - Htemp=55+ii*5 - outf(14,ii)=-1. - if(Htemp.ge.65.) outf(14,ii)=XE6(Htemp) + Htemp=55+ii*5 + outf(14,ii)=-1. + if(Htemp.ge.65.) outf(14,ii)=XE6(Htemp) outf(14,11+ii)=-1. call F00(Htemp,LATI,DAYNR,XHI,F107D,EDENS,IERROR) if(ierror.eq.0.or.ierror.eq.2) outf(14,11+ii)=edens - outf(14,22+ii)=ddens(1,ii) - outf(14,33+ii)=ddens(2,ii) - outf(14,44+ii)=ddens(3,ii) - outf(14,55+ii)=ddens(4,ii) - outf(14,66+ii)=ddens(5,ii) + outf(14,22+ii)=ddens(1,ii) + outf(14,33+ii)=ddens(2,ii) + outf(14,44+ii)=ddens(3,ii) + outf(14,55+ii)=ddens(4,ii) + outf(14,66+ii)=ddens(5,ii) enddo endif @@ -2130,11 +2154,11 @@ C OARR(25)=DIP OARR(26)=MAGBR OARR(27)=MODIP - OARR(28)=LATI + OARR(28)=LATI OARR(29)=SAX200 OARR(30)=SUX200 OARR(31)=SEASON - OARR(32)=LONGI + OARR(32)=LONGI OARR(33)=rssn OARR(34)=COV OARR(35)=B1 @@ -2163,7 +2187,7 @@ C OARR(37) used for TEC and 38 for TEC-top c include only every second auroral boundary point (MLT=0,1,2..23) jjj=58 do iii=1,47,2 - jjj=jjj+1 + jjj=jjj+1 oarr(jjj)=ab_mlat(iii) enddo OARR(83)=xkp @@ -2174,24 +2198,24 @@ c include only every second auroral boundary point (MLT=0,1,2..23) c output of solar indices used c write(6,10201) iyyyy,rssn,gind,cov,covsat,f107d,f10781, -c & f107365,pf107,cov-f10781,cov-f107365,cov-pf107 +c & f107365,pf107,cov-f10781,cov-f107365,cov-pf107 c10201 format(I5,11F6.1) icalls=icalls+1 - RETURN - END -c -c - subroutine iri_web(jmag,jf,alati,along,iyyyy,mmdd,iut,dhour, - & height,h_tec_max,ivar,vbeg,vend,vstp,a,b) -c----------------------------------------------------------------------- + + END subroutine iri_sub + + + subroutine iri_web(jmag,jf,alati,along,iyyyy,mmdd,iut,dhour, + & height,Nalt,h_tec_max,ivar,a,b) +c----------------------------------------------------------------------- c changes: c 11/16/99 jf(30) instead of jf(17) c 10/31/08 outf, a, b (100 -> 500) C 11/05/16 array size change a(30,1000), outf(30,1000) c -c----------------------------------------------------------------------- +c----------------------------------------------------------------------- c input: jmag,alati,along,iyyyy,mmdd,dhour see IRI_SUB c height height in km c h_tec_max =0 no TEC otherwise upper boundary for integral @@ -2206,35 +2230,44 @@ c output: a similar to outf in IRI_SUB c b similar to oarr in IRI_SUB c c numstp number of steps; maximal 1000 -c----------------------------------------------------------------------- - dimension outf(30,1000),oar(100),oarr(100),a(30,1000) - dimension xvar(8),b(100,1000) +c----------------------------------------------------------------------- + real outf(30,1000),oar(100),oarr(100) + real, intent(out) :: a(30,1000),b(100,1000) + dimension xvar(8) logical jf(50) + real :: height(:) + character(256) dirdata1 - nummax=1000 - numstp=int((vend-vbeg)/vstp)+1 - if(numstp.gt.nummax) numstp=nummax + common /folders/ dirdata1 + + nummax=1000 + !numstp=int((vend-vbeg)/vstp)+1 + !if(numstp.gt.nummax) numstp=nummax do 6249 i=1,100 -6249 oar(i)=b(i,1) +6249 oar(i)=b(i,1) if(ivar.eq.1) then do 1249 i=1,100 -1249 oarr(i)=oar(i) +1249 oarr(i)=oar(i) xhour=dhour+iut*25. +! TODO: need to pass vector of altitudes call IRI_SUB(JF,JMAG,ALATI,ALONG,IYYYY,MMDD,XHOUR, - & VBEG,VEND,VSTP,a,OARR) - if(h_tec_max.gt.50.) then + & height,Nalt,dirdata1,a,OARR) + if(h_tec_max.gt.50.) then call iri_tec (50.,h_tec_max,2,tec,tect,tecb) oarr(37)=tec oarr(38)=tect - endif + endif do 1111 i=1,100 1111 b(i,1)=oarr(i) return - endif + endif + + where (height <= 0.0) + height=100. + endwhere - if(height.le.0.0) height=100 xvar(2)=alati xvar(3)=along xvar(4)=iyyyy @@ -2248,18 +2281,18 @@ c----------------------------------------------------------------------- alati=xvar(2) along=xvar(3) iyyyy=int(xvar(4)) - if(ivar.eq.7) then + if(ivar == 7) then mmdd=-int(vbeg) else mmdd=int(xvar(5)*100+xvar(6)) endif dhour=xvar(8)+iut*25. - do 1 i=1,numstp + do i=1,numstp do 1349 iii=1,100 -1349 oarr(iii)=b(iii,i) +1349 oarr(iii)=b(iii,i) call IRI_SUB(JF,JMAG,ALATI,ALONG,IYYYY,MMDD,DHOUR, - & height,height,1.,OUTF,OARR) + & height,1,dirdata1,OUTF,OARR) if(h_tec_max.gt.50.) then call iri_tec (50.,h_tec_max,2,tec,tect,tecb) oarr(37)=tec @@ -2280,8 +2313,6 @@ c----------------------------------------------------------------------- mmdd=int(xvar(5)*100+xvar(6)) endif dhour=xvar(8)+iut*25. -1 continue - - return - end + enddo + end subroutine iri_web diff --git a/fortran/iritec.for b/fortran/iritec.for index 2b56932..2fae315 100644 --- a/fortran/iritec.for +++ b/fortran/iritec.for @@ -79,8 +79,8 @@ c call iri_tec (hbeg,hend,2,tec,tect,tecb) end subroutine IRIT13 -c -c + + pure real function ioncorr(tec,f) c----------------------------------------------------------------------- c computes ionospheric correction IONCORR (in m) for given vertical @@ -91,8 +91,8 @@ c----------------------------------------------------------------------- ioncorr = 40.3 * tec / (f*f) end function ioncorr -c -c + + subroutine iri_tec (hstart,hend,istep,tectot,tectop,tecbot) c----------------------------------------------------------------------- C subroutine to compute the total ionospheric content diff --git a/fortran/iritest.for b/fortran/iritest.for index 7450618..3fc5a89 100644 --- a/fortran/iritest.for +++ b/fortran/iritest.for @@ -135,7 +135,7 @@ c read(5,*) jchoice do i=1,50 jf(i)=.true. - enddo + enddo if(piktab.eq.4) jf(24)=.false. if(jchoice.eq.0) then c defaults for jf(1:50) diff --git a/fortran/iriwebg.for b/fortran/iriwebg.for index 122df88..6bdaf96 100644 --- a/fortran/iriwebg.for +++ b/fortran/iriwebg.for @@ -1,5 +1,5 @@ - subroutine iriwebg(inJMAG,inJF,inALATI,inALONG,inIYYYY, + subroutine iriwebg(inJMAG,inJF,inALATI,inALONG,inIYYYY, & inMMDD,inIUT,inDHOUR,inHEIGHT,inH_TEC_MAX, & inIVAR,inVBEG,inVEND,inVSTP,inADDINP,dirdata,outA,outB) @@ -18,11 +18,11 @@ Cf2py intent(in) inJMAG, inJF, inALATI, inALONG, inIYYYY, inMMDD Cf2py intent(in) inIUT, inDHOUR, inHEIGHT, inH_TEC_MAX, inIVAR Cf2py intent(in) inVBEG, inVEND, inVSTP, inADDINP, dirdata - + common /folders/ dirdata1 dirdata1 = trim(dirdata) - + call read_ig_rz call readapf107 @@ -59,21 +59,21 @@ C hmF2 or M(3000)F2 jf(9) = .false. b(2,1) = addinp(2) endif -C Ne(300km) +C Ne(300km) if(addinp(3).ne.-1) then jf(10) = .false. do i = 1, 1000 b(15,i) = addinp(3) end do endif -C Ne(400km) +C Ne(400km) if(addinp(4).ne.-1) then jf(10) = .false. do i = 1, 1000 b(16,i) = addinp(4) end do endif -C Ne(550km) +C Ne(550km) if(addinp(5).ne.-1) then jf(10) = .false. jf(23) = .false. @@ -81,34 +81,34 @@ C Ne(550km) b(16,i) = addinp(5) end do endif -C foF1 or NmF1 +C foF1 or NmF1 if(addinp(6).ne.-1) then jf(13) = .false. b(3,1) = addinp(6) endif -C hmF1 +C hmF1 if(addinp(7).ne.-1) then jf(14) = .false. b(4,1) = addinp(7) endif -C foE or NmE +C foE or NmE if(addinp(8).ne.-1) then jf(15) = .false. b(5,1) = addinp(8) endif -C hmE +C hmE if(addinp(9).ne.-1) then jf(16) = .false. b(6,1) = addinp(9) endif -C Rz12 +C Rz12 if(addinp(10).ne.-1) then jf(17) = .false. do i = 1, 1000 b(33,i) = addinp(10) end do endif -C F10.7 +C F10.7 if(addinp(11).ne.-1) then jf(21) = .true. jf(23) = .false. @@ -119,7 +119,7 @@ C F10.7 b(41,i) = addinp(11) end do endif -C IG12 +C IG12 if(addinp(12).ne.-1) then jf(27) = .false. do i = 1, 1000 @@ -152,7 +152,7 @@ C IG12 integer jmag,iyyyy,mmdd real alati,along,dhour,heibeg,heiend,heistp character*256 dirdata,dirdata1 - + Cf2py intent(in) jf,jmag,iyyyy,mmdd,alati,along,dhour Cf2py intent(in) heibeg,heiend,heistp,dirdata @@ -197,13 +197,13 @@ Cf2py integer intent(hide),depend(coordl) :: lenl=shape(coordl,0) do i=1,lenl - along = real(coordl(i,1),kind(along)) + along = real(coordl(i,1),kind(along)) alati = real(coordl(i,3),kind(alati)) heibeg = real(coordl(i,2),kind(heibeg)) heiend = heibeg + 1.0 heistp = 1.0 - + call iri_sub(jf,jag,alati,along,iyyyy,mmdd,dhour, & heibeg,heiend,heistp,outf,oarr) @@ -212,8 +212,8 @@ Cf2py integer intent(hide),depend(coordl) :: lenl=shape(coordl,0) end do do j=1,100 - oarr1(j,i) = oarr(j) - end do + oarr1(j,i) = oarr(j) + end do end do @@ -227,7 +227,7 @@ Cf2py integer intent(hide),depend(coordl) :: lenl=shape(coordl,0) integer yyyy,ddd,lenl,i real coordl(lenl,3) - real uhour + real uhour character*256 dirdata,dirdata1 integer mm,dd,nrdaymo,nmonth @@ -237,21 +237,21 @@ Cf2py integer intent(hide),depend(coordl) :: lenl=shape(coordl,0) Cf2py intent(in) yyyy,ddd,uhour,coordl,dirdata Cf2py integer intent(hide),depend(coordl) :: lenl=shape(coordl,0) - + common /folders/ dirdata1 - dirdata1 = trim(dirdata) + dirdata1 = trim(dirdata) call initialize call read_ig_rz - call moda(1,yyyy,mm,dd,ddd,nrdaymo) + call moda(1,yyyy,mm,dd,ddd,nrdaymo) call tcon(yyyy,mm,dd,ddd,rz,ig,rsn,nmonth) f107d = 63.75 + rz(3) * (0.728 + rz(3) * 0.00089) do i=1,lenl glon = real(coordl(i,1),kind(glon)) - hei = real(coordl(i,2),kind(hei)) + hei = real(coordl(i,2),kind(hei)) glat = real(coordl(i,3),kind(glat)) call ut_lt(0,uhour,lhour,glon,yyyy,ddd) @@ -260,7 +260,7 @@ Cf2py integer intent(hide),depend(coordl) :: lenl=shape(coordl,0) call f00(hei,glat,ddd,xhi,f107d1,edens,ierr) edens1(i) = edens ierr1(i) = ierr - + end do end subroutine firisubl diff --git a/fortran/test.f90 b/fortran/test.f90 new file mode 100644 index 0000000..27fda47 --- /dev/null +++ b/fortran/test.f90 @@ -0,0 +1,45 @@ +program basictest +implicit none + +logical, parameter :: jf(50) = .true. +integer, parameter :: jmag = 1, iyyyy=1980,mmdd=0321,dhour=12 +real, parameter :: glat=0., glon=0. +real,parameter :: altkm(*) = [130., 140., 150.] +integer,parameter :: Nalt = size(altkm) +character(*), parameter :: datadir = 'data/' + +real :: oarr(100), outf(30,Nalt) +integer :: i + +do i=1,2 + call IRI_SUB(JF,JMAG,glat,glon,IYYYY,MMDD,DHOUR+25, altkm,Nalt,datadir,OUTF,OARR) + + print '(A,ES10.3,A,F5.1,A)','NmF2 ',oarr(1),' [m^-3] hmF2 ',oarr(2),' [km] ' + print '(A,F10.3,A,I3,A,F10.3)','F10.7 ',oarr(41), ' Ap ',int(oarr(51)),' B0 ',oarr(10) + print *,'Ne ',outf(1,:) +enddo + +end program + + +! logical, parameter :: jf(50) = .true. +! integer, parameter :: jmag = 1, iyyyy=1980,mmdd=0321,dhour=12 +! real, parameter :: glat=0., glon=0. +! real,parameter :: altkm(*) = [130., 140., 150.] + +! ** IRI parameters are being calculated *** +!Ne: IRI-2001 for Topside +!Ne, foF2: CCIR model is used. +!Ne: B0,B1 Bil-2000 +!Ne, foF1: probability function used. +!Ne, D: IRI1990 +!Ne, foF2: storm model included +!Ion Com.: DS-95 & DY-85 +!Te: Aeros/AE/ISIS model +!Auroral boundary model on +!Ne, foE: storm model on +! +! NmF2 1.362E+12 [m^-3] hmF2 327.5 [km] +! F10.7 162.400 Ap 18 B0 88.054 +! Ne 1.05585536E+09 815260608. 1.58453235E+09 + diff --git a/pyiri2016/__init__.py b/pyiri2016/__init__.py index 9fd15f8..ef61869 100644 --- a/pyiri2016/__init__.py +++ b/pyiri2016/__init__.py @@ -1,244 +1,234 @@ from pathlib import Path -import datetime +from datetime import datetime, timedelta from dateutil.parser import parse import xarray -from iri2016 import iriwebg -from numpy import arange, nan, ones, squeeze - -class IRI2016(object): - - def __init__(self): - self.iriDataFolder = Path(__file__).parent / 'data' - - - def Switches(self): - """ - IRI switches to turn on/off several options - """ - - jf = ones(50) - - # i 1 0 standard version - # ------------------------------------------------------------------------ - jf[ 1 - 1] = 1; # 1 Ne computed Ne not computed 1 - jf[ 2 - 1] = 1; # 2 Te, Ti computed Te, Ti not computed 1 - jf[ 3 - 1] = 1; # 3 Ne & Ni computed Ni not computed 1 - jf[ 4 - 1] = 0; # 4 B0 - Table option B0 - other models jf(31) 0 - jf[ 5 - 1] = 0; # 5 foF2 - CCIR foF2 - URSI 0 - jf[ 6 - 1] = 0; # 6 Ni - DS-95 & DY-85 Ni - RBV-10 & TTS-03 0 - jf[ 7 - 1] = 1; # 7 Ne - Tops: f10.7<188 f10.7 unlimited 1 - jf[ 8 - 1] = 1; # 8 foF2 from model foF2 or NmF2 - user input 1 - jf[ 9 - 1] = 1; # 9 hmF2 from model hmF2 or M3000F2 - user input 1 - jf[10 - 1] = 1; # 10 Te - Standard Te - Using Te/Ne correlation 1 - jf[11 - 1] = 1; # 11 Ne - Standard Profile Ne - Lay-function formalism 1 - jf[12 - 1] = 1; # 12 Messages to unit 6 to meesages.text on unit 11 1 - jf[13 - 1] = 1; # 13 foF1 from model foF1 or NmF1 - user input 1 - jf[14 - 1] = 1; # 14 hmF1 from model hmF1 - user input (only Lay version)1 - jf[15 - 1] = 1; # 15 foE from model foE or NmE - user input 1 - jf[16 - 1] = 1; # 16 hmE from model hmE - user input 1 - jf[17 - 1] = 1; # 17 Rz12 from file Rz12 - user input 1 - jf[18 - 1] = 1; # 18 IGRF dip, magbr, modip old FIELDG using POGO68/10 for 1973 1 - jf[19 - 1] = 1; # 19 F1 probability model critical solar zenith angle (old) 1 - jf[20 - 1] = 1; # 20 standard F1 standard F1 plus L condition 1 - jf[21 - 1] = 1; # 21 ion drift computed ion drift not computed 0 - jf[22 - 1] = 0; # 22 ion densities in % ion densities in m-3 1 - jf[23 - 1] = 0; # 23 Te_tops (Aeros,ISIS) Te_topside (TBT-2011) 0 - jf[24 - 1] = 0; # 24 D-region: IRI-95 Special: 3 D-region models (FIRI) 1 - jf[25 - 1] = 1; # 25 F107D from APF107.DAT F107D user input (oarr(41)) 1 - jf[26 - 1] = 0; # 26 foF2 storm model no storm updating 1 - jf[27 - 1] = 1; # 27 IG12 from file IG12 - user 1 - jf[28 - 1] = 0; # 28 spread-F probability not computed 0 - jf[29 - 1] = 0; # 29 IRI01-topside new options as def. by JF(30) 0 - jf[30 - 1] = 0; # 30 IRI01-topside corr. NeQuick topside model 0 - # (29,30) = (1,1) IRIold, (0,1) IRIcor, (0,0) NeQuick, (1,0) Gulyaeva - jf[31 - 1] = 1; # 31 B0,B1 ABT-2009 B0 Gulyaeva h0.5 1 - jf[32 - 1] = 1; # 32 F10.7_81 from file PF10.7_81 - user input (oarr(46)) 1 - jf[33 - 1] = 0; # 33 Auroral boundary model on/off true/false 0 - jf[34 - 1] = 0; # 34 Messages on Messages off 1 - jf[35 - 1] = 0; # 35 foE storm model no foE storm updating 0 - # .. .... - # 50 .... - # ------------------------------------------------------------------ - - return jf +import numpy as np +# +import iri2016 # fortran + +proot = Path(__file__).parents[1] +simout = ['ne','Tn','Ti','Te','nO+','nH+','nHe+','nO2+','nNO+'] + +def datetimerange(start:datetime, end:datetime, step:timedelta) -> list: + """like range() for datetime!""" + if isinstance(start,str): + start = parse(start) + + if isinstance(end,str): + end = parse(end) + + assert isinstance(start, datetime) + assert isinstance(end, datetime) + assert isinstance(step, timedelta) + + return [start + i*step for i in range((end-start) // step)] + + +def Switches(): + """ + IRI switches to turn on/off several options + """ + + jf = np.ones(50,dtype=bool) + + # i 1 0 standard version + # ------------------------------------------------------------------------ + jf[ 1 - 1] = 1; # 1 Ne computed Ne not computed 1 + jf[ 2 - 1] = 1; # 2 Te, Ti computed Te, Ti not computed 1 + jf[ 3 - 1] = 1; # 3 Ne & Ni computed Ni not computed 1 + jf[ 4 - 1] = 0; # 4 B0 - Table option B0 - other models jf(31) 0 + jf[ 5 - 1] = 0; # 5 foF2 - CCIR foF2 - URSI 0 + jf[ 6 - 1] = 0; # 6 Ni - DS-95 & DY-85 Ni - RBV-10 & TTS-03 0 + jf[ 7 - 1] = 1; # 7 Ne - Tops: f10.7<188 f10.7 unlimited 1 + jf[ 8 - 1] = 1; # 8 foF2 from model foF2 or NmF2 - user input 1 + jf[ 9 - 1] = 1; # 9 hmF2 from model hmF2 or M3000F2 - user input 1 + jf[10 - 1] = 1; # 10 Te - Standard Te - Using Te/Ne correlation 1 + jf[11 - 1] = 1; # 11 Ne - Standard Profile Ne - Lay-function formalism 1 + jf[12 - 1] = 1; # 12 Messages to unit 6 to meesages.text on unit 11 1 + jf[13 - 1] = 1; # 13 foF1 from model foF1 or NmF1 - user input 1 + jf[14 - 1] = 1; # 14 hmF1 from model hmF1 - user input (only Lay version)1 + jf[15 - 1] = 1; # 15 foE from model foE or NmE - user input 1 + jf[16 - 1] = 1; # 16 hmE from model hmE - user input 1 + jf[17 - 1] = 1; # 17 Rz12 from file Rz12 - user input 1 + jf[18 - 1] = 1; # 18 IGRF dip, magbr, modip old FIELDG using POGO68/10 for 1973 1 + jf[19 - 1] = 1; # 19 F1 probability model critical solar zenith angle (old) 1 + jf[20 - 1] = 1; # 20 standard F1 standard F1 plus L condition 1 + jf[21 - 1] = 1; # 21 ion drift computed ion drift not computed 0 + jf[22 - 1] = 0; # 22 ion densities in % ion densities in m-3 1 + jf[23 - 1] = 0; # 23 Te_tops (Aeros,ISIS) Te_topside (TBT-2011) 0 + jf[24 - 1] = 1; # 24 D-region: IRI-95 Special: 3 D-region models (FIRI) 1 + jf[25 - 1] = 1; # 25 F107D from APF107.DAT F107D user input (oarr(41)) 1 + jf[26 - 1] = 0; # 26 foF2 storm model no storm updating 1 + jf[27 - 1] = 1; # 27 IG12 from file IG12 - user 1 + jf[28 - 1] = 0; # 28 spread-F probability not computed 0 + jf[29 - 1] = 0; # 29 IRI01-topside new options as def. by JF(30) 0 + jf[30 - 1] = 0; # 30 IRI01-topside corr. NeQuick topside model 0 + # (29,30) = (1,1) IRIold, (0,1) IRIcor, (0,0) NeQuick, (1,0) Gulyaeva + jf[31 - 1] = 1; # 31 B0,B1 ABT-2009 B0 Gulyaeva h0.5 1 + jf[32 - 1] = 1; # 32 F10.7_81 from file PF10.7_81 - user input (oarr(46)) 1 + jf[33 - 1] = 0; # 33 Auroral boundary model on/off true/false 0 + jf[34 - 1] = 0; # 34 Messages on Messages off 1 + jf[35 - 1] = 0; # 35 foE storm model no foE storm updating 0 + # .. .... + # 50 .... + # ------------------------------------------------------------------ + + return jf #%% +def IRI( time, altkm, glat, glon, ap=None, f107=None, ssn=None, var=None): - def IRI(self, ap=5, f107=150, glat=0., glon=0., time=datetime.datetime.now(), - ssn=150, var=1, vbeg=130., vend=130.+1., vstp=1.): + if isinstance(time, str): + time = parse(time) - if isinstance(time, str): - time = parse(time) -# doy = squeeze(TimeUtilities().CalcDOY(year, month, dom)) - - # IRI options - jf = self.Switches() - - # additional "input parameters" (necessary to scale the empirical model results - # to measurements) - addinp = -ones(12) - - #------------------------------------------------------------------------------ - # - if time.year < 1958: - - addinp[10 - 1] = ssn # RZ12 (This switches 'jf[17 - 1]' to '0' and - # uses correlation function to estimate IG12) - - jf[25 - 1] = 1 # 25 F107D from APF107.DAT F107D user input (oarr(41)) 1 - jf[27 - 1] = 1 # 27 IG12 from file IG12 - user 1 - jf[32 - 1] = 1 # 32 F10.7_81 from file PF10.7_81 - user input (oarr(46)) 1 + altkm = np.atleast_1d(altkm) - else: # case for solar and geomagnetic indices from files +# doy = squeeze(TimeUtilities().CalcDOY(year, month, dom)) - jf[17 - 1] = 1 # 17 Rz12 from file Rz12 - user input 1 - jf[26 - 1] = 1 # 26 foF2 storm model no storm updating 1 - jf[35 - 1] = 1 # 35 foE storm model no foE storm updating 0 - # - #------------------------------------------------------------------------------ + # IRI options + jf = Switches() + + # additional "input parameters" (necessary to scale the empirical model results + # to measurements) +# addinp = -np.ones(12) + + #------------------------------------------------------------------------------ + # +# if time.year < 1958: +# +# addinp[10 - 1] = ssn # RZ12 (This switches 'jf[17 - 1]' to '0' and +# # uses correlation function to estimate IG12) +# +# jf[25 - 1] = 1 # 25 F107D from APF107.DAT F107D user input (oarr(41)) 1 +# jf[27 - 1] = 1 # 27 IG12 from file IG12 - user 1 +# jf[32 - 1] = 1 # 32 F10.7_81 from file PF10.7_81 - user input (oarr(46)) 1 +# +# else: # case for solar and geomagnetic indices from files +# +# jf[17 - 1] = 1 # 17 Rz12 from file Rz12 - user input 1 +# jf[26 - 1] = 1 # 26 foF2 storm model no storm updating 1 +# jf[35 - 1] = 1 # 35 foE storm model no foE storm updating 0 +# # +# #------------------------------------------------------------------------------ + + mmdd = 100*time.month + time.day # month and dom (MMDD) + # hour + 25 denotes UTC time + dhour = (time.hour + 25) + time.minute/60. +# %% more inputs + jmag = 0 # 0: geographic; 1: geomagnetic + # iut = 0 # 0: for LT; 1: for UT + # height = 300. # in km + # h_tec_max = 2000 # 0: no TEC; otherwise: upper boundary for integral + # ivar = var # 1: altitude; 2: latitude; 3: longitude; ... - mmdd = int(1e2 * time.month) + time.day # month and dom (MMDD) + # Ionosphere (IRI) +# a, b = iriwebg(jmag, jf, glat, glon, int(time.year), mmdd, iut, time.hour, +# height, h_tec_max, ivar, ivbeg, ivend, ivstp, addinp, self.iriDataFolder) -# %% more inputs - jmag = 0 # 0: geographic; 1: geomagnetic - iut = 0 # 0: for LT; 1: for UT - height = 300. # in km - h_tec_max = 2000 # 0: no TEC; otherwise: upper boundary for integral - ivar = var # 1: altitude; 2: latitude; 3: longitude; ... - ivbeg = vbeg - ivend = vend - ivstp = vstp - # Ionosphere (IRI) - a, b = iriwebg(jmag, jf, glat, glon, int(time.year), mmdd, iut, time.hour, - height, h_tec_max, ivar, ivbeg, ivend, ivstp, addinp, self.iriDataFolder) + outf,oarr = iri2016.iri_sub(jf, jmag, glat, glon, + time.year, mmdd, dhour, altkm, + proot/'data/') - bins = arange(ivbeg, ivend + ivstp * 0., ivstp) - a = a[:, arange(len(bins))] - b = b[:, arange(len(bins))] +# %% collect output + dsf = {k: (('time','alt_km','lat','lon'),np.atleast_2d(v[None,:,None,None])) for (k,v) in zip(simout, outf[:9,:])} -# %% - # IRI Standard Ne (in m-3) - neIRI = squeeze(self._RmZeros(self._RmNeg(a[0, :]))[0]) + dsf.update({'NmF2':(('time','lat','lon'),np.atleast_3d(oarr[0]))}) + dsf.update({'hmF2':(('time','lat','lon'),np.atleast_3d(oarr[1]))}) + dsf.update({'NmF1':(('time','lat','lon'),np.atleast_3d(oarr[2]))}) + dsf.update({'hmF1':(('time','lat','lon'),np.atleast_3d(oarr[3]))}) + dsf.update({'NmE':(('time','lat','lon'),np.atleast_3d(oarr[4]))}) + dsf.update({'hmE':(('time','lat','lon'),np.atleast_3d(oarr[5]))}) + dsf.update({'B0':(('time','lat','lon'),np.atleast_3d(oarr[9]))}) - # IRI Temperature (in K) - teIRI = squeeze(a[4 - 1, :][0]) - tiIRI = squeeze(a[3 - 1, :][0]) + iri = xarray.Dataset(dsf, + coords={'time':[time],'alt_km':altkm,'lat':[glat],'lon':[glon]}, + attrs={'f107':oarr[40], 'ap':oarr[50], + 'glat':glat,'glon':glon,'time':time, + }) - # FIRI Ne (in m-3) - iri_ne_firi = squeeze(self._RmNeg(a[13 - 1, :])[0]) +# FIRI Ne (in m-3) +# iri_ne_firi = self._RmNeg(a[13 - 1, :])[0] - ######### Ionic density (NO+, O2+, O+, H+, He+, N+, Cluster Ions) - # Ionic density (O+, O2+, NO+) - oplusIRI = squeeze(self._RmZeros(a[5 - 1, :])[0]) / 100. * neIRI # in m-3 - o2plusIRI = squeeze(self._RmZeros(a[8 - 1, :])[0]) / 100. * neIRI # in m-3 - noplusIRI = squeeze(self._RmZeros(a[9 - 1, :])[0]) / 100. * neIRI # in m-3 +## Ionic density (NO+, O2+, O+, H+, He+, N+, Cluster Ions) - # more ionic densities (H+, He+, N+) - hplusIRI = squeeze(self._RmZeros(a[6 - 1, :])[0]) / 100. * neIRI # in m-3 - heplusIRI = squeeze(self._RmZeros(a[7 - 1, :])[0]) / 100. * neIRI # in m-3 - nplusIRI = squeeze(self._RmZeros(a[11 - 1, :])[0]) / 100. * neIRI # in m-3 +# iri = {'ne' : neIRI, 'te' : teIRI, 'ti' : tiIRI, 'neFIRI' : iri_ne_firi, +# 'oplus' : oplusIRI, 'o2plus' : o2plusIRI, 'noplus' : noplusIRI, +# 'hplus' : hplusIRI, 'heplus' : heplusIRI, 'nplus' : nplusIRI} - iri = {'ne' : neIRI, 'te' : teIRI, 'ti' : tiIRI, 'neFIRI' : iri_ne_firi, - 'oplus' : oplusIRI, 'o2plus' : o2plusIRI, 'noplus' : noplusIRI, - 'hplus' : hplusIRI, 'heplus' : heplusIRI, 'nplus' : nplusIRI} +# iriadd = { 'NmF2' : b[1 - 1, :][0], 'hmF2' : b[2 - 1, :][0], +# 'B0' : b[10 - 1, :][0] } - iriadd = { 'NmF2' : b[1 - 1, :][0], 'hmF2' : b[2 - 1, :][0], - 'B0' : b[10 - 1, :][0] } + return iri - return iri, iriadd +def timeprofile(tlim:tuple, dt:timedelta, + altkm:np.ndarray, glat:float, glon:float) -> xarray.Dataset: + """compute IRI90 altitude profile over time range for fixed lat/lon + """ - def _RmZeros(self, inputs): - """ Replace "zero" values with 'NaN' - """ - inputs[inputs == 0.0] = nan + T = datetimerange(tlim[0], tlim[1], dt) - return inputs + altkm = np.atleast_1d(altkm) + iono = None - def _RmNeg(self, inputs): - """ Replace negative values with 'NaN' """ - inputs[inputs < 0.0] = nan + f107 =[]; ap=[] + for t in T: + iri = IRI(t, altkm, glat, glon) + if iono is None: + iono = iri + else: + iono = xarray.merge((iono,iri)) - return inputs + f107.append(iri.f107) + ap.append(iri.ap) -class IRI2016Profile(IRI2016): + iono.attrs = iri.attrs + iono.attrs['f107'] = f107 + iono.attrs['ap'] = ap - def __init__(self, alt=None, altlim=[90.,150.], altstp=2., htecmax=0, - time=datetime.datetime.now(), hrlim=[0., 24.], hrstp=None, - iut=1, jmag=0, - lat=0., latlim=[-90, 90], latstp=10., - lon=0., lonlim=[-180,180], lonstp=20., - option='vertical', verbose=False): + return iono - if isinstance(time,str): - time = parse(time) - self.iriDataFolder = Path(__file__).parent / 'data' +def geoprofile(latlim:tuple, dlat:float, glon:float, + altkm:np.ndarray, time:datetime) -> xarray.Dataset: + """compute IRI90 altitude profiles at time, over lat or lon range + """ - self.jf = self.Switches() + glat = np.arange(*latlim,dlat) - self.addinp = list(map(lambda x : -1, range(12))) + altkm = np.atleast_1d(altkm) - self.option = option + iono = None - if option == 'vertical': # Height Profile - self.simtype = 1 - self.vbeg = altlim[0] - self.vend = altlim[1] - self.vstp = altstp - elif option == 'lat': # Latitude Profile - self.simtype = 2 - self.vbeg = latlim[0] - self.vend = latlim[1] - self.vstp = latstp - elif option == 'lon': # Longitude Profile - self.simtype = 3 - self.vbeg = lonlim[0] - self.vend = lonlim[1] - self.vstp = lonstp - elif option == 'time': # Local Time Profile - self.simtype = 8 - self.vbeg = hrlim[0] - self.vend = hrlim[1] - self.vstp = hrstp + f107 =[]; ap=[] + for l in glat: + iri = IRI(time, altkm, l, glon) + if iono is None: + iono = iri else: - raise ValueError(f'Invalid option {option}') + iono = xarray.merge((iono,iri)) + f107.append(iri.f107) + ap.append(iri.ap) - self.htecmax = htecmax - self.jmag = jmag - self.lat = lat - self.lon = lon - self.month = time.month - self.dom = time.day - self.mmdd = self.month * 100 + self.dom - self.year = time.year - self.iut = iut - self.hour = time.hour - self.alt = alt - self.verbose = verbose - self.numstp = int((self.vend - self.vbeg) / self.vstp) + 1 + iono.attrs = iri.attrs + iono.attrs['f107'] = f107 + iono.attrs['ap'] = ap - if option == 'vertical': - self.HeiProfile() - elif option == 'lat': - self.LatProfile() - elif option == 'lon': - self.LonProfile() - elif option == 'time': - self.HrProfile() + return iono - def _CallIRI(self): - self.a, self.b = iriwebg(self.jmag, self.jf, self.lat, self.lon, self.year, self.mmdd, - self.iut, self.hour, self.alt, self.htecmax, self.simtype, self.vbeg, - self.vend, self.vstp, self.addinp, self.iriDataFolder) +# def _CallIRI(self): +# +# self.a, self.b = iriwebg(self.jmag, self.jf, self.lat, self.lon, self.year, self.mmdd, +# self.iut, self.hour, self.alt, self.htecmax, self.simtype, self.vbeg, +# self.vend, self.vstp, self.addinp, self.iriDataFolder) def _Hr2HHMMSS(self): @@ -298,53 +288,3 @@ def HeiProfile(self): print('%8.3f %10.3f %8.3f %8.3f %8.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %6.3f %7.3f %7.3f %7.3f %7.3f' % (varval, edens, edratio, a[2,i], a[3,i], a[4,i], a[5,i], a[10,i], a[6,i], a[7,i], a[8,i], a[9,i], b[32,i], b[38,i], b[40,i], b[45,i], b[50,i], b[51,i])) - - - def LatProfile(self): - - self._CallIRI() - - self._GetTitle() - - if self.verbose: - - lat = arange(self.vbeg, self.numstp*self.vstp, self.vstp) - - print('\tGLON\tGLAT\tNmF2\t\thmF2\tB0') - for j in range(lat.size): - print('%8.3f %8.3f %8.3e %8.3f %8.3f' % (self.lon, lat[j], self.b[0, j], self.b[1, j], self.b[9, j])) - - - def LonProfile(self): - - self._CallIRI() - - self._GetTitle() - - if self.verbose: - - lon = arange(self.vbeg, self.numstp*self.vstp, self.vstp) - - print('\tGLON\tGLAT\tNmF2\t\thmF2\tB0') - for j in range(lon.size): - print('%8.3f %8.3f %8.3e %8.3f %8.3f' % (lon[j], self.lat, self.b[0, j], self.b[1, j], self.b[9, j])) - - - def HrProfile(self): - - self._CallIRI() - - self._GetTitle() - - t = arange(self.vbeg,self.numstp*self.vstp,self.vstp) - - self.out = xarray.DataArray(self.a[:9,:t.size].T, - coords={'time':t, - 'sim':['ne','Tn','Ti','Te','nO+','nH+','nHe+','nO2+','nNO+']}, - dims=['time','sim']) - - if self.verbose: - - print(' GLON GLAT\tHR\tNmF2\thmF2\tB0') - for j in range(t.size): - print('%8.3f %8.3f %8.3f %8.3e %8.3f %8.3f' % (self.lon, self.lat, t[j], self.b[0, j], self.b[1, j], self.b[9, j])) diff --git a/pyiri2016/iri2016prof2D.py b/pyiri2016/iri2016prof2D.py index 4720520..3bcfed6 100644 --- a/pyiri2016/iri2016prof2D.py +++ b/pyiri2016/iri2016prof2D.py @@ -3,7 +3,7 @@ from numpy import arange, array, ceil, empty, floor, isnan, linspace, \ log10, meshgrid, nan, tile, transpose, where, hypot from numpy.ma import masked_where -from matplotlib.pyplot import close, cm, colorbar, figure, savefig +from matplotlib.pyplot import close, cm, colorbar, figure try: from mpl_toolkits.basemap import Basemap except ImportError: @@ -351,8 +351,8 @@ def PlotLatVsFLFIRI(self, save=False, verbose=False): figname = gpath / 'firi-{:02d}{:02d}.jpg'.format(self.time[0], self.time[1]) if verbose: - print('Saving at: {:s}'.format(figname)) - savefig(str(figname), bbox_inches='tight', format='jpg', dpi=100) + print('Saving',figname) + fg.savefig(figname, bbox_inches='tight', format='jpg', dpi=100) close(fg) @@ -363,11 +363,11 @@ def Plot2D(self, save=False): logging.error('TODO: this needs to be updated to cartopy') return - f = figure(figsize=(24, 6)) + fg = figure(figsize=(24, 6)) if self.option == 1: - pn = f.add_subplot(131) + pn = fg.add_subplot(131) X, Y = meshgrid(self.data2D['hour'], self.data2D['alt']) ipc = pn.pcolor(X, Y, transpose(log10(self.data2D['Ne'])), cmap=cm.jet, vmax=13, vmin=9) pn.set_title(self.data2D['title1']) @@ -376,7 +376,7 @@ def Plot2D(self, save=False): cp1 = colorbar(ipc) cp1.set_label('Log$_{10}$N$_e$(m$^{-3}$)') - pn = f.add_subplot(132) + pn = fg.add_subplot(132) ipc = pn.pcolor(X, Y, transpose(self.data2D['Te']), cmap=cm.jet, vmax=4000, vmin=100) pn.set_title(self.data2D['title2']) pn.set_xlabel('Hour (UT)') @@ -384,7 +384,7 @@ def Plot2D(self, save=False): cp1 = colorbar(ipc) cp1.set_label('T$_e$ ($^\circ$)') - pn = f.add_subplot(133) + pn = fg.add_subplot(133) ipc = pn.pcolor(X, Y, transpose(self.data2D['Ti']), cmap=cm.jet, vmax=4000, vmin=100) pn.set_xlabel('Hour (UT)') pn.set_ylabel('Altitude (km)') @@ -394,7 +394,7 @@ def Plot2D(self, save=False): elif self.option == 2: - pn1 = f.add_subplot(111) + pn1 = fg.add_subplot(111) m = Basemap(llcrnrlon=self.data2D['lon'][0], llcrnrlat=self.data2D['lat'][0], \ urcrnrlon=self.data2D['lon'][-1], urcrnrlat=self.data2D['lat'][-1], \ @@ -426,7 +426,7 @@ def Plot2D(self, save=False): gpath.mkdir(parents=True, exist_ok=True) figname = gpath / 'iri-{:02d}{:02d}.jpg'.format(self.HH, self.MM) - savefig(str(figname), bbox_inches='tight', format='jpg', dpi=100) + fg.savefig(figname, bbox_inches='tight', format='jpg', dpi=100) # convert -resize 50% -delay 20 -loop 0 *.jpg myimage.gif diff --git a/setup.py b/setup.py index d473d2d..d162c40 100644 --- a/setup.py +++ b/setup.py @@ -1,6 +1,6 @@ #!/usr/bin/env python install_requires = ['python-dateutil','numpy','xarray'] -tests_require = ['nose','coveralls'] +tests_require = ['pytest','nose','coveralls'] name = 'pyiri2016' # %% from setuptools import find_packages @@ -9,7 +9,8 @@ from os.path import join -src = ['iriwebg.for', 'irisub.for', 'irifun.for', +src = [#'iriwebg.for', + 'irisub.for', 'irifun.for', 'iritec.for', 'iridreg.for', 'igrf.for', 'cira.for', 'iriflip.for'] src = [join('fortran', s) for s in src] @@ -37,17 +38,17 @@ setup(name=name, packages=find_packages(), - version='1.3.1', + version='1.4.0', author=['Michael Hirsch, Ph.D.','Ronald Ilma'], url = 'https://github.com/scivision/pyIRI2016', description='IRI2016 International Reference Ionosphere from Python', long_description=open('README.rst').read(), classifiers=[ - 'Intended Audience :: Science/Research', 'Development Status :: 5 - Production/Stable', + 'Intended Audience :: Science/Research', 'License :: OSI Approved :: MIT License', + 'Programming Language :: Python :: 3.6', 'Topic :: Scientific/Engineering :: Atmospheric Science', - 'Programming Language :: Python :: 3', ], ext_modules = [ext], data_files = iriDataFiles, @@ -56,4 +57,6 @@ 'tests':tests_require,}, tests_require = tests_require, python_requires='>=3.6', + scripts=['AltitudeProfile.py','TimeProfile.py','LatitudeProfile.py'], + include_package_data=True, ) diff --git a/tests/test_all.py b/tests/test_all.py index a794ba0..e2c91ab 100755 --- a/tests/test_all.py +++ b/tests/test_all.py @@ -1,15 +1,16 @@ #!/usr/bin/env python import numpy as np -import unittest -from pyiri2016 import IRI2016 +import pyiri2016 -class BasicTests(unittest.TestCase): - def test_main1(self): +def test_main1(): + + iri = pyiri2016.IRI('1980-03-21T12', 130., 0., 0.) + + np.testing.assert_allclose((iri['ne'].item(), iri.NmF2, iri.hmF2), + (267285184512.0, 2580958937088.0, 438.78643798828125)) + + print('assert passed') - Obj = IRI2016() - IRIData, IRIDATAAdd = Obj.IRI(time='1980-03-21T12') - np.testing.assert_allclose((IRIData['ne'], IRIDATAAdd['NmF2'], IRIDATAAdd['hmF2']), - (267285184512.0, 2580958937088.0, 438.78643798828125)) if __name__ == '__main__': - unittest.main() + np.testing.run_module_suite()