Skip to content

Allow variable-length inputs and outputs #23

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
AntoineBrunet wants to merge 5 commits into
base: main
Choose a base branch
from
Open

Allow variable-length inputs and outputs #23

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
5be6f8a
Changed ntime_max to ntime on most routines to allow variable-length …
AntoineBrunet Aug 27, 2021
d94a586
Removed unnecessary 'ntime_max.inc' includes. Put back 'ntime_max' ar…
AntoineBrunet Aug 27, 2021
069eb73
Correction of typos in comments
AntoineBrunet Apr 11, 2024
10137c7
Updated documentation to reflect the possibility to use variable leng…
AntoineBrunet Apr 16, 2024
3400ec8
Uniform treatment of Nene in radiation models, using variable array d…
AntoineBrunet Apr 16, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
185 changes: 87 additions & 98 deletions docs/frames/frames.html
View file
Open in desktop

Large diffs are not rendered by default.

30 changes: 14 additions & 16 deletions source/AE8_AP8.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -51,13 +51,13 @@
! whatf -> which kind of flux, 1=differential 2=E range 3=integral (long integer)
! Nene -> Number of energy channels to compute
! energy -> energy (MeV) at which fluxes must be computed (double array [2,25])
! iyear,idoy,UT -> times when flux are to be computed (not usefull if imput position is not in GSE, GSM, SM,GEI) (respectively long array(ntime_max), long array(ntime_max), double array(ntime_max))
! iyear,idoy,UT -> times when flux are to be computed (not useful if input position is not in GSE, GSM, SM,GEI) (respectively long array(ntime), long array(ntime), double array(ntime))
! xIN1 -> first coordinate in the chosen system (double array [ntime_max])
! xIN2 -> second coordinate in the chosen system (double array [ntime_max])
! xIN3 -> third coordinate in the chosen system (double array [ntime_max])
!
! OUTPUTS:
! flux -> Computed fluxes (MeV-1 cm-2 s-1 or cm-2 s-1) (double array [ntime_max,25])
! flux -> Computed fluxes (MeV-1 cm-2 s-1 or cm-2 s-1) (double array [ntime_max,Nene])
!
! COMMON BLOCKS:
! COMMON/GENER/ERA,AQUAD,BQUAD
Expand All @@ -79,13 +79,11 @@ SUBROUTINE fly_in_nasa_aeap1(ntime,sysaxes,whichm,whatf,nene,
INCLUDE 'ntime_max.inc'
C
c declare inputs
INTEGER*4 nene_max
PARAMETER (nene_max=25)
INTEGER*4 ntime,sysaxes,whichm,whatf,Nene
INTEGER*4 iyear(ntime_max),idoy(ntime_max)
REAL*8 energy(2,nene_max)
REAL*8 UT(ntime_max)
real*8 xIN1(ntime_max),xIN2(ntime_max),xIN3(ntime_max)
INTEGER*4 iyear(ntime),idoy(ntime)
REAL*8 energy(2,nene)
REAL*8 UT(ntime)
real*8 xIN1(ntime),xIN2(ntime),xIN3(ntime)
c Declare internal variables
INTEGER*4 k_ext,k_l,isat,kint
INTEGER*4 t_resol,r_resol,Ilflag
Expand All @@ -94,11 +92,11 @@ SUBROUTINE fly_in_nasa_aeap1(ntime,sysaxes,whichm,whatf,nene,
REAL*8 xGSM(3),xSM(3),xGEI(3),xGSE(3)
real*8 alti,lati,longi,UT_dip,psi,tilt
REAL*8 ERA,AQUAD,BQUAD
REAL*8 BLOCAL(ntime_max),BMIN(ntime_max),XJ(ntime_max)
REAL*8 Lm(ntime_max),Lstar(ntime_max),BBo(ntime_max)
REAL*8 BLOCAL(ntime),BMIN(ntime),XJ(ntime)
REAL*8 Lm(ntime),Lstar(ntime),BBo(ntime)
c
c Declare output variables
REAL*8 flux(ntime_max,nene_max)
REAL*8 flux(ntime_max,Nene)
C
COMMON/GENER/ERA,AQUAD,BQUAD
COMMON /magmod/k_ext,k_l,kint
Expand Down Expand Up @@ -193,12 +191,12 @@ SUBROUTINE fly_in_nasa_aeap1(ntime,sysaxes,whichm,whatf,nene,
! whichm -> which model to use, 1=AE8min 2=AE8max 3=AP8min 4=AP8max (long integer)
! whatf -> which kind of flux, 1=differential 2=E range 3=integral (long integer)
! Nene -> Number of energy channels to compute
! energy -> energy (MeV) at which fluxes must be computed (double array [2,25])
! energy -> energy (MeV) at which fluxes must be computed (double array [2,Nene])
! BBo -> Blocal/Bequator (double array [ntime_max])
! L -> McIlwain L (double array [ntime_max])
!
! OUTPUTS:
! flux -> Computed fluxes (MeV-1 cm-2 s-1 or cm-2 s-1) (double array [ntime_max,25])
! flux -> Computed fluxes (MeV-1 cm-2 s-1 or cm-2 s-1) (double array [ntime_max,Nene])
!
! COMMON BLOCKS:
! COMMON /PROMIN/ IHEADPMIN, MAPPMIN
Expand Down Expand Up @@ -233,8 +231,8 @@ SUBROUTINE get_AE8_AP8_flux(ntmax,whichm,whatf,nene,
INTEGER*4 whatf
INTEGER*4 MAPPMIN(16582), MAPPMAX(16291),
& MAPEMIN(13168), MAPEMAX(13548)
REAL*8 energy(2,25)
REAL*8 flux(ntime_max,25),BBo(ntime_max),L(ntime_max)
REAL*8 energy(2,nene)
REAL*8 flux(ntime_max,nene),BBo(ntime_max),L(ntime_max)
REAL*8 FL,FL1
INTEGER*4 IHEADPMIN(8),IHEADPMAX(8),IHEADEMIN(8),IHEADEMAX(8)
c
Expand All @@ -245,7 +243,7 @@ SUBROUTINE get_AE8_AP8_flux(ntmax,whichm,whatf,nene,
c
c init
DO i=1,ntmax
do ieny=1,25
do ieny=1,Nene
Flux(i,ieny) = baddata
enddo
enddo
Expand Down
30 changes: 14 additions & 16 deletions source/AFRL_CRRES_models.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -40,14 +40,14 @@
! whichm -> which model to use, 1=pro quiet 2=pro active 3=ele ave 4=ele worst case 5=ele Ap15 (long integer)
! whatf -> which kind of flux, 1=differential 2=E range 3=integral (long integer)
! Nene -> Number of energy channels to compute
! energy -> energy (MeV) at which fluxes must be computed (double array [2,25])
! iyear,idoy,UT -> times when flux are to be computed (not usefull if imput position is not in GSE, GSM, SM,GEI) (respectively long array(ntime_max), long array(ntime_max), double array(ntime_max))
! energy -> energy (MeV) at which fluxes must be computed (double array [2,Nene])
! iyear,idoy,UT -> times when flux are to be computed (not useful if input position is not in GSE, GSM, SM,GEI) (respectively long array(ntime), long array(ntime), double array(ntime))
! xIN1 -> first coordinate in the chosen system (double array [ntime_max])
! xIN2 -> second coordinate in the chosen system (double array [ntime_max])
! xIN3 -> third coordinate in the chosen system (double array [ntime_max])
! Ap15 -> 15 previous days average of Ap index assuming a one day delay (double array [ntime_max])
!
! OUTPUT: flux -> Computed fluxes (MeV-1 cm-2 s-1 or cm-2 s-1) (double array [ntime_max,25])
! OUTPUT: flux -> Computed fluxes (MeV-1 cm-2 s-1 or cm-2 s-1) (double array [ntime_max,Nene])
!
! CALLING SEQUENCE: CALL fly_in_afrl_crres1(ntime,sysaxes,whichm,whatf,energy,xIN1,xIN2,xIN3,flux)
!---------------------------------------------------------------------------------------------------
Expand All @@ -62,13 +62,11 @@ SUBROUTINE fly_in_afrl_crres1(ntime,sysaxes,whichm,whatf,nene,
INTEGER*4 STRLEN
BYTE ascii_path(strlen)
c
INTEGER*4 nene_max
PARAMETER (nene_max=25)
INTEGER*4 ntime,sysaxes,whichm,whatf,nene
INTEGER*4 iyear(ntime_max),idoy(ntime_max)
REAL*8 energy(2,nene_max)
REAL*8 UT(ntime_max)
real*8 xIN1(ntime_max),xIN2(ntime_max),xIN3(ntime_max)
INTEGER*4 iyear(ntime),idoy(ntime)
REAL*8 energy(2,nene)
REAL*8 UT(ntime)
real*8 xIN1(ntime),xIN2(ntime),xIN3(ntime)
c Declare internal variables
INTEGER*4 k_ext,k_l,isat,kint
INTEGER*4 t_resol,r_resol,Ilflag
Expand All @@ -78,12 +76,12 @@ SUBROUTINE fly_in_afrl_crres1(ntime,sysaxes,whichm,whatf,nene,
REAL*8 xGSM(3),xSM(3),xGEI(3),xGSE(3)
real*8 alti,lati,longi,psi,tilt
REAL*8 ERA,AQUAD,BQUAD
REAL*8 BLOCAL(ntime_max),BMIN(ntime_max),XJ(ntime_max)
REAL*8 Lm(ntime_max),Lstar(ntime_max),BBo(ntime_max)
REAL*8 Ap15(ntime_max)
REAL*8 BLOCAL(ntime),BMIN(ntime),XJ(ntime)
REAL*8 Lm(ntime),Lstar(ntime),BBo(ntime)
REAL*8 Ap15(ntime)
c
c Declare output variables
REAL*8 flux(ntime_max,nene_max)
REAL*8 flux(ntime_max,nene)
c
CHARACTER*(500) afrl_crres_path
C
Expand Down Expand Up @@ -186,8 +184,8 @@ SUBROUTINE get_crres_flux(ntmax,whichm,whatf,nene,
c
INTEGER*4 Ne,Nl,Nbb0,ind
c
REAL*8 energy(2,25)
REAL*8 flux(ntime_max,25),BBo(ntime_max),L(ntime_max)
REAL*8 energy(2,nene)
REAL*8 flux(ntime_max,nene),BBo(ntime_max),L(ntime_max)
REAL*8 Ap15(ntime_max)
REAL*8 pente,cste,Flux1,Flux2
c
Expand All @@ -209,7 +207,7 @@ SUBROUTINE get_crres_flux(ntmax,whichm,whatf,nene,
c
c init
DO i=1,ntmax
do ieny=1,25
do ieny=1,nene
Flux(i,ieny) = baddata
enddo
enddo
Expand Down
10 changes: 3 additions & 7 deletions source/CoordTrans.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,8 +23,7 @@
! SUBROUTINE coord_trans1: Generic coordinate transformation from one Earth or Heliospheric coordinate
! system to another one
! SUBROUTINE coord_trans_vec1: Generic coordinate transformation from one Earth or Heliospheric coordinate
! system to another one (handle up to
! ntime_max positions)
! system to another one
!
!***************************************************************************************************
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
Expand Down Expand Up @@ -2706,8 +2705,6 @@ SUBROUTINE coord_trans1(sysaxesIN,sysaxesOUT,iyr,idoy,
C REAL*8 xINV(3,nmax), xOUTV(3,nmax)
C INTEGER*4 numpoints
C
C As with all (most?) onera library calls, the maximum array size
C is limited to ntime_max elements
C Contributed by Timothy Guild, 2.2.07
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
Expand All @@ -2716,12 +2713,11 @@ SUBROUTINE coord_trans_vec1(ntime,sysaxesIN,sysaxesOUT,
& iyear,idoy,secs,xINV,xOUTV)

IMPLICIT NONE
INCLUDE 'ntime_max.inc'
INCLUDE 'variables.inc'

INTEGER*4 nmax,i,ntime, sysaxesIN, sysaxesOUT
INTEGER*4 iyear(ntime_max),idoy(ntime_max),y,d
REAL*8 secs(ntime_max),xINV(3,ntime_max),xOUTV(3,ntime_max)
INTEGER*4 iyear(ntime),idoy(ntime),y,d
REAL*8 secs(ntime),xINV(3,ntime),xOUTV(3,ntime)
! local vars
REAL*8 xIN(3),xOUT(3),s

Expand Down
34 changes: 16 additions & 18 deletions source/LAndI2Lstar.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20443,16 +20443,15 @@ SUBROUTINE LAndI2Lstar1(ntime,kext,options,sysaxes,iyearsat,
c
IMPLICIT NONE
INCLUDE 'variables.inc'
INCLUDE 'ntime_max.inc'
c
c declare inputs
INTEGER*4 kext,k_ext,k_l,options(5)
INTEGER*4 ntime,sysaxes
INTEGER*4 iyearsat(ntime_max)
integer*4 idoy(ntime_max)
real*8 UT(ntime_max)
real*8 xIN1(ntime_max),xIN2(ntime_max),xIN3(ntime_max)
real*8 maginput(25,ntime_max)
INTEGER*4 iyearsat(ntime)
integer*4 idoy(ntime)
real*8 UT(ntime)
real*8 xIN1(ntime),xIN2(ntime),xIN3(ntime)
real*8 maginput(25,ntime)
c 1: Kp
c 2: Dst
c 3: dens
Expand All @@ -20468,9 +20467,9 @@ SUBROUTINE LAndI2Lstar1(ntime,kext,options,sysaxes,iyearsat,
INTEGER*4 option1,isat
c
c Declare output variables
REAL*8 BLOCAL(ntime_max),BMIN(ntime_max),XJ(ntime_max)
REAL*8 MLT(ntime_max)
REAL*8 Lm(ntime_max),Lstar(ntime_max)
REAL*8 BLOCAL(ntime),BMIN(ntime),XJ(ntime)
REAL*8 MLT(ntime)
REAL*8 Lm(ntime),Lstar(ntime)
C
c This method to compute L* is only available for IGRF + Olson-Pfitzer quiet
if (options(5) .ne. 0) options(5)=0 ! force internal field to be IGRF
Expand Down Expand Up @@ -20552,9 +20551,9 @@ SUBROUTINE LAndI2Lstar_shell_splitting1(ntime,Nipa,kext,options,
c Declare internal variables
INTEGER*4 option1,isat,IPA,ntime_tmp,sysaxesOUT,sysaxesIN
REAL*8 alti, lati, longi
REAL*8 BLOCAL_tmp(ntime_max),BMIN_tmp(ntime_max)
REAL*8 XJ_tmp(ntime_max),MLT_tmp(ntime_max)
REAL*8 Lm_tmp(ntime_max),Lstar_tmp(ntime_max)
REAL*8 BLOCAL_tmp(ntime),BMIN_tmp(ntime)
REAL*8 XJ_tmp(ntime),MLT_tmp(ntime)
REAL*8 Lm_tmp(ntime),Lstar_tmp(ntime)
REAL*8 xIN(3),xOUT(3),BL,BMIR,xGEO(3)
REAL*8 maginput_tmp(25)
INTEGER*4 imagin
Expand Down Expand Up @@ -20670,14 +20669,13 @@ SUBROUTINE EmpiricalLstar1(ntime,kext,options,iyearsat,idoy,
c
IMPLICIT NONE
INCLUDE 'variables.inc'
INCLUDE 'ntime_max.inc'
c
c declare inputs
INTEGER*4 kext,options(5)
INTEGER*4 ntime
INTEGER*4 iyearsat(ntime_max)
integer*4 idoy(ntime_max)
real*8 maginput(25,ntime_max)
INTEGER*4 iyearsat(ntime)
integer*4 idoy(ntime)
real*8 maginput(25,ntime)
c 1: Kp
c 2: Dst
c 3: dens
Expand All @@ -20689,7 +20687,7 @@ SUBROUTINE EmpiricalLstar1(ntime,kext,options,iyearsat,idoy,
c 9: G2
c 10: G3
c
REAL*8 XJ(ntime_max),Lm(ntime_max)
REAL*8 XJ(ntime),Lm(ntime)
c
c Declare internal variables
INTEGER*4 isat,DayIndexL,DayIndexR,i,iflag,kint
Expand All @@ -20708,7 +20706,7 @@ SUBROUTINE EmpiricalLstar1(ntime,kext,options,iyearsat,idoy,
common /rconst/rad,pi
c
c Declare output variables
REAL*8 Lstar(ntime_max)
REAL*8 Lstar(ntime)
C
COMMON/LAndI2LstarCom/Lmax,Imax,Lupper,Iupper,Lm4,A0,A1,A2,A3,A4,
&Lm5,A50,A51,A52,A53,A54,A55
Expand Down
52 changes: 25 additions & 27 deletions source/get_bderivs.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,14 +23,13 @@ SUBROUTINE GET_Bderivs(ntime,kext,options,sysaxes,dX,
C computes derivatives of B (vector and magnitude)
C inputs: ntime through maginput have the usual meaning, except dX
C REAL*8 dX is the step size, in RE for the numerical derivatives (recommend 1E-3?)
C real*8 Bgeo(3,ntime_max) - components of B in GEO, nT
C real*8 Bmag(ntime_max) - magnitude of B in nT
C real*8 gradBmag(3,ntime_max) - gradient of Bmag in GEO, nT/RE
C real*8 diffB(3,3,ntime_max) - derivatives of Bgeo in GEO, nT/RE
C real*8 Bgeo(3,ntime) - components of B in GEO, nT
C real*8 Bmag(ntime) - magnitude of B in nT
C real*8 gradBmag(3,ntime) - gradient of Bmag in GEO, nT/RE
C real*8 diffB(3,3,ntime) - derivatives of Bgeo in GEO, nT/RE
C diffB(i,j,t) = dB_i/dx_j for point t (t=1 to ntime)

IMPLICIT NONE
INCLUDE 'ntime_max.inc' ! include file created by make, defines ntime_max
INCLUDE 'variables.inc'
C
COMMON /magmod/k_ext,k_l,kint
Expand All @@ -39,17 +38,17 @@ SUBROUTINE GET_Bderivs(ntime,kext,options,sysaxes,dX,
INTEGER*4 ntime,kext,options(5)
INTEGER*4 sysaxes
REAL*8 dX
INTEGER*4 iyearsat(ntime_max)
integer*4 idoy(ntime_max)
real*8 UT(ntime_max)
real*8 xIN1(ntime_max),xIN2(ntime_max),xIN3(ntime_max)
real*8 maginput(25,ntime_max)
INTEGER*4 iyearsat(ntime)
integer*4 idoy(ntime)
real*8 UT(ntime)
real*8 xIN1(ntime),xIN2(ntime),xIN3(ntime)
real*8 maginput(25,ntime)

c declare outputs
real*8 Bgeo(3,ntime_max) ! components of B in GEO, nT
real*8 Bmag(ntime_max) ! magnitude of B in nT
real*8 gradBmag(3,ntime_max) ! gradient of Bmag in GEO, nT/RE
real*8 diffB(3,3,ntime_max) ! derivatives of Bgeo in GEO, nT/RE
real*8 Bgeo(3,ntime) ! components of B in GEO, nT
real*8 Bmag(ntime) ! magnitude of B in nT
real*8 gradBmag(3,ntime) ! gradient of Bmag in GEO, nT/RE
real*8 diffB(3,3,ntime) ! derivatives of Bgeo in GEO, nT/RE

c declare internal variables
integer*4 isat
Expand Down Expand Up @@ -133,26 +132,25 @@ SUBROUTINE compute_grad_curv_curl(ntime,Bgeo,Bmag,gradBmag,diffB,
C computes gradient factors, curvature factors, and curl of B

IMPLICIT NONE
INCLUDE 'ntime_max.inc' ! include file created by make, defines ntime_max
INCLUDE 'variables.inc'

C all coordinates are in GEO reference frame
C inputs:
integer*4 ntime ! number of points
real*8 Bgeo(3,ntime_max) ! components of B in GEO, nT
real*8 Bmag(ntime_max) ! magnitude of B in nT
real*8 gradBmag(3,ntime_max) ! gradient of Bmag in GEO, nT/RE
real*8 diffB(3,3,ntime_max) ! derivatives of Bgeo in GEO, nT/RE
real*8 Bgeo(3,ntime) ! components of B in GEO, nT
real*8 Bmag(ntime) ! magnitude of B in nT
real*8 gradBmag(3,ntime) ! gradient of Bmag in GEO, nT/RE
real*8 diffB(3,3,ntime) ! derivatives of Bgeo in GEO, nT/RE
c diffB(i,j,t) = dB_i/dx_j for point t (t=1 to ntime)
c outputs:
real*8 grad_par(ntime_max) ! gradient of Bmag along B nT/RE
real*8 grad_perp(3,ntime_max) ! gradient of Bmag perpendicular to B nT/RE
real*8 grad_drift(3,ntime_max) ! (bhat x grad_perp)/B, 1/RE (part of gradient drift velocity)
real*8 curvature(3,ntime_max) ! (bhat dot grad)bhat, 1/RE (part of curvature force)
real*8 Rcurv(ntime_max) ! 1/|curvature| RE (radius of curvature)
real*8 curv_drift(3,ntime_max) ! (bhat x curvature), 1/RE (part of curvature drift)
real*8 curlB(3,ntime_max) ! curl of B (nT/RE) (part of electrostatic current term)
real*8 divB(ntime_max) ! divergence of B (nT/RE) (should be zero!)
real*8 grad_par(ntime) ! gradient of Bmag along B nT/RE
real*8 grad_perp(3,ntime) ! gradient of Bmag perpendicular to B nT/RE
real*8 grad_drift(3,ntime) ! (bhat x grad_perp)/B, 1/RE (part of gradient drift velocity)
real*8 curvature(3,ntime) ! (bhat dot grad)bhat, 1/RE (part of curvature force)
real*8 Rcurv(ntime) ! 1/|curvature| RE (radius of curvature)
real*8 curv_drift(3,ntime) ! (bhat x curvature), 1/RE (part of curvature drift)
real*8 curlB(3,ntime) ! curl of B (nT/RE) (part of electrostatic current term)
real*8 divB(ntime) ! divergence of B (nT/RE) (should be zero!)


c internal variables
Expand Down
15 changes: 7 additions & 8 deletions source/get_hemi.f
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -91,22 +91,21 @@ SUBROUTINE GET_HEMI1(kext,options,sysaxes,iyearsat,idoy,UT,

SUBROUTINE GET_HEMI_MULTI(ntime,kext,options,sysaxes,iyearsat,
& idoy,UT,xIN1,xIN2,xIN3,maginput,xHEMI)
c calls get_hemi1 multiple times (ntime, <= ntime_max)
c calls get_hemi1 multiple times (ntime, <= ntime)
IMPLICIT NONE
INCLUDE 'variables.inc'
INCLUDE 'ntime_max.inc' ! include file created by make, defines ntime_max

c declare inputs
INTEGER*4 ntime,kext,options(5)
INTEGER*4 sysaxes
INTEGER*4 iyearsat(ntime_max)
integer*4 idoy(ntime_max)
real*8 UT(ntime_max)
real*8 xIN1(ntime_max),xIN2(ntime_max),xIN3(ntime_max)
real*8 maginput(25,ntime_max)
INTEGER*4 iyearsat(ntime)
integer*4 idoy(ntime)
real*8 UT(ntime)
real*8 xIN1(ntime),xIN2(ntime),xIN3(ntime)
real*8 maginput(25,ntime)

c declare outputs
integer*4 xHEMI(ntime_max)
integer*4 xHEMI(ntime)

c declare internal variables
integer*4 isat
Expand Down
Loading