From 5be6f8a5b0f72876a3384a8740418cf5c084c774 Mon Sep 17 00:00:00 2001 From: Antoine Brunet Date: Fri, 27 Aug 2021 16:41:38 +0200 Subject: [PATCH 1/5] Changed ntime_max to ntime on most routines to allow variable-length arrays. --- source/AE8_AP8.f | 12 ++++----- source/AFRL_CRRES_models.f | 20 +++++++-------- source/CoordTrans.f | 4 +-- source/LAndI2Lstar.f | 46 +++++++++++++++++------------------ source/get_bderivs.f | 50 +++++++++++++++++++------------------- source/get_hemi.f | 14 +++++------ source/onera_desp_lib.f | 50 +++++++++++++++++++------------------- 7 files changed, 98 insertions(+), 98 deletions(-) diff --git a/source/AE8_AP8.f b/source/AE8_AP8.f index b8ffc9db..99c16b56 100644 --- a/source/AE8_AP8.f +++ b/source/AE8_AP8.f @@ -51,7 +51,7 @@ ! 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 usefull if imput 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]) @@ -82,10 +82,10 @@ SUBROUTINE fly_in_nasa_aeap1(ntime,sysaxes,whichm,whatf,nene, INTEGER*4 nene_max PARAMETER (nene_max=25) INTEGER*4 ntime,sysaxes,whichm,whatf,Nene - INTEGER*4 iyear(ntime_max),idoy(ntime_max) + INTEGER*4 iyear(ntime),idoy(ntime) REAL*8 energy(2,nene_max) - REAL*8 UT(ntime_max) - real*8 xIN1(ntime_max),xIN2(ntime_max),xIN3(ntime_max) + 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 @@ -94,8 +94,8 @@ 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) diff --git a/source/AFRL_CRRES_models.f b/source/AFRL_CRRES_models.f index f58a4791..b0305e74 100644 --- a/source/AFRL_CRRES_models.f +++ b/source/AFRL_CRRES_models.f @@ -41,7 +41,7 @@ ! 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 usefull if imput 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]) @@ -65,10 +65,10 @@ SUBROUTINE fly_in_afrl_crres1(ntime,sysaxes,whichm,whatf,nene, INTEGER*4 nene_max PARAMETER (nene_max=25) INTEGER*4 ntime,sysaxes,whichm,whatf,nene - INTEGER*4 iyear(ntime_max),idoy(ntime_max) + INTEGER*4 iyear(ntime),idoy(ntime) REAL*8 energy(2,nene_max) - REAL*8 UT(ntime_max) - real*8 xIN1(ntime_max),xIN2(ntime_max),xIN3(ntime_max) + 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 @@ -78,9 +78,9 @@ 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) @@ -186,8 +186,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 @@ -209,7 +209,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 diff --git a/source/CoordTrans.f b/source/CoordTrans.f index 64d9752d..dd965ece 100644 --- a/source/CoordTrans.f +++ b/source/CoordTrans.f @@ -2720,8 +2720,8 @@ SUBROUTINE coord_trans_vec1(ntime,sysaxesIN,sysaxesOUT, 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 diff --git a/source/LAndI2Lstar.f b/source/LAndI2Lstar.f index 9ecac9d8..a0d20b59 100644 --- a/source/LAndI2Lstar.f +++ b/source/LAndI2Lstar.f @@ -20448,11 +20448,11 @@ SUBROUTINE LAndI2Lstar1(ntime,kext,options,sysaxes,iyearsat, 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 @@ -20468,9 +20468,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 @@ -20532,12 +20532,12 @@ SUBROUTINE LAndI2Lstar_shell_splitting1(ntime,Nipa,kext,options, INTEGER*4 kext,k_ext,k_l,options(5),Nalp,Nipa PARAMETER (Nalp=25) 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) + INTEGER*4 iyearsat(ntime) + integer*4 idoy(ntime) + real*8 UT(ntime) + real*8 xIN1(ntime),xIN2(ntime),xIN3(ntime) real*8 alpha(Nalp) - real*8 maginput(25,ntime_max) + real*8 maginput(25,ntime) c 1: Kp c 2: Dst c 3: dens @@ -20552,16 +20552,16 @@ 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 c c Declare output variables - REAL*8 BLOCAL(ntime_max,Nalp),BMIN(ntime_max) - REAL*8 XJ(ntime_max,Nalp),MLT(ntime_max) + REAL*8 BLOCAL(ntime_max,Nalp),BMIN(ntime) + REAL*8 XJ(ntime_max,Nalp),MLT(ntime) REAL*8 Lm(ntime_max,Nalp),Lstar(ntime_max,Nalp) C c This method to compute L* is only available for IGRF + Olson-Pfitzer quiet @@ -20675,9 +20675,9 @@ SUBROUTINE EmpiricalLstar1(ntime,kext,options,iyearsat,idoy, 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 @@ -20689,7 +20689,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 @@ -20708,7 +20708,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 diff --git a/source/get_bderivs.f b/source/get_bderivs.f index 5b2cad50..7cc6f48d 100644 --- a/source/get_bderivs.f +++ b/source/get_bderivs.f @@ -23,10 +23,10 @@ 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 @@ -39,17 +39,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 @@ -139,20 +139,20 @@ SUBROUTINE compute_grad_curv_curl(ntime,Bgeo,Bmag,gradBmag,diffB, 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 diff --git a/source/get_hemi.f b/source/get_hemi.f index ba4a1a30..0b2af9d9 100644 --- a/source/get_hemi.f +++ b/source/get_hemi.f @@ -91,7 +91,7 @@ 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 @@ -99,14 +99,14 @@ SUBROUTINE GET_HEMI_MULTI(ntime,kext,options,sysaxes,iyearsat, 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 diff --git a/source/onera_desp_lib.f b/source/onera_desp_lib.f index 08a968d2..8094cefc 100644 --- a/source/onera_desp_lib.f +++ b/source/onera_desp_lib.f @@ -56,11 +56,11 @@ SUBROUTINE make_lstar1(ntime,kext,options,sysaxes,iyearsat, 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 c Declare internal variables INTEGER*4 isat,iyear,kint,ifail @@ -70,9 +70,9 @@ SUBROUTINE make_lstar1(ntime,kext,options,sysaxes,iyearsat, real*8 alti,lati,longi 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 COMMON /magmod/k_ext,k_l,kint COMMON /flag_L/Ilflag @@ -173,12 +173,12 @@ SUBROUTINE make_lstar_shell_splitting1(ntime,Nipa,kext,options, INTEGER*4 kext,k_ext,k_l,options(5),Nalp,Nipa PARAMETER (Nalp=25) 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) + INTEGER*4 iyearsat(ntime) + integer*4 idoy(ntime) + real*8 UT(ntime) + real*8 xIN1(ntime),xIN2(ntime),xIN3(ntime) real*8 alpha(Nalp) - real*8 maginput(25,ntime_max) + real*8 maginput(25,ntime) c c c Declare internal variables @@ -190,8 +190,8 @@ SUBROUTINE make_lstar_shell_splitting1(ntime,Nipa,kext,options, real*8 alti,lati,longi c c Declare output variables - REAL*8 BLOCAL(ntime_max,Nalp),BMIN(ntime_max) - REAL*8 XJ(ntime_max,Nalp),MLT(ntime_max) + REAL*8 BLOCAL(ntime_max,Nalp),BMIN(ntime) + REAL*8 XJ(ntime_max,Nalp),MLT(ntime) REAL*8 Lm(ntime_max,Nalp),Lstar(ntime_max,Nalp) C COMMON /magmod/k_ext,k_l,kint @@ -290,15 +290,15 @@ SUBROUTINE Lstar_Phi1(ntime,whichinv,options,iyearsat, c declare inputs INTEGER*4 whichinv,options(5) INTEGER*4 ntime - INTEGER*4 iyearsat(ntime_max) - integer*4 idoy(ntime_max) + INTEGER*4 iyearsat(ntime) + integer*4 idoy(ntime) c c Declare internal variables INTEGER*4 isat,kint REAL*8 Bo,xc,yc,zc,ct,st,cp,sp c c Declare output variables - REAL*8 Phi(ntime_max),Lstar(ntime_max) + REAL*8 Phi(ntime),Lstar(ntime) C COMMON /dipigrf/Bo,xc,yc,zc,ct,st,cp,sp REAL*8 pi,rad @@ -748,7 +748,7 @@ SUBROUTINE GET_FIELD1(kext,options,sysaxes,iyearsat,idoy,UT, SUBROUTINE GET_FIELD_MULTI(ntime,kext,options,sysaxes,iyearsat, & idoy,UT,xIN1,xIN2,xIN3,maginput,BxGEO,Bl) -C Call get_field1 many times (ntime, in fact, up to ntime = ntime_max) +C Call get_field1 many times (ntime, in fact, up to ntime = ntime) c IMPLICIT NONE INCLUDE 'variables.inc' @@ -758,14 +758,14 @@ SUBROUTINE GET_FIELD_MULTI(ntime,kext,options,sysaxes,iyearsat, INTEGER*4 ntime INTEGER*4 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 output variables - REAL*8 BxGEO(3,ntime_max),Bl(ntime_max) + REAL*8 BxGEO(3,ntime),Bl(ntime) C c Declare internal variables integer*4 isat From d94a586539ddd2fb565d371164975fa10c245b63 Mon Sep 17 00:00:00 2001 From: Antoine Brunet Date: Fri, 27 Aug 2021 17:07:01 +0200 Subject: [PATCH 2/5] Removed unnecessary 'ntime_max.inc' includes. Put back 'ntime_max' array sizes on drift_shell routines inputs to be coherent with outputs. --- source/LAndI2Lstar.f | 15 +++++++-------- source/onera_desp_lib.f | 17 +++++++---------- 2 files changed, 14 insertions(+), 18 deletions(-) diff --git a/source/LAndI2Lstar.f b/source/LAndI2Lstar.f index a0d20b59..39d06e44 100644 --- a/source/LAndI2Lstar.f +++ b/source/LAndI2Lstar.f @@ -20443,7 +20443,6 @@ 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) @@ -20532,12 +20531,12 @@ SUBROUTINE LAndI2Lstar_shell_splitting1(ntime,Nipa,kext,options, INTEGER*4 kext,k_ext,k_l,options(5),Nalp,Nipa PARAMETER (Nalp=25) INTEGER*4 ntime,sysaxes - INTEGER*4 iyearsat(ntime) - integer*4 idoy(ntime) - real*8 UT(ntime) - real*8 xIN1(ntime),xIN2(ntime),xIN3(ntime) + 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 alpha(Nalp) - real*8 maginput(25,ntime) + real*8 maginput(25,ntime_max) c 1: Kp c 2: Dst c 3: dens @@ -20560,8 +20559,8 @@ SUBROUTINE LAndI2Lstar_shell_splitting1(ntime,Nipa,kext,options, INTEGER*4 imagin c c Declare output variables - REAL*8 BLOCAL(ntime_max,Nalp),BMIN(ntime) - REAL*8 XJ(ntime_max,Nalp),MLT(ntime) + REAL*8 BLOCAL(ntime_max,Nalp),BMIN(ntime_max) + REAL*8 XJ(ntime_max,Nalp),MLT(ntime_max) REAL*8 Lm(ntime_max,Nalp),Lstar(ntime_max,Nalp) C c This method to compute L* is only available for IGRF + Olson-Pfitzer quiet diff --git a/source/onera_desp_lib.f b/source/onera_desp_lib.f index 8094cefc..e1238bf2 100644 --- a/source/onera_desp_lib.f +++ b/source/onera_desp_lib.f @@ -51,7 +51,6 @@ SUBROUTINE make_lstar1(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) @@ -173,12 +172,12 @@ SUBROUTINE make_lstar_shell_splitting1(ntime,Nipa,kext,options, INTEGER*4 kext,k_ext,k_l,options(5),Nalp,Nipa PARAMETER (Nalp=25) INTEGER*4 ntime,sysaxes - INTEGER*4 iyearsat(ntime) - integer*4 idoy(ntime) - real*8 UT(ntime) - real*8 xIN1(ntime),xIN2(ntime),xIN3(ntime) + 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 alpha(Nalp) - real*8 maginput(25,ntime) + real*8 maginput(25,ntime_max) c c c Declare internal variables @@ -190,8 +189,8 @@ SUBROUTINE make_lstar_shell_splitting1(ntime,Nipa,kext,options, real*8 alti,lati,longi c c Declare output variables - REAL*8 BLOCAL(ntime_max,Nalp),BMIN(ntime) - REAL*8 XJ(ntime_max,Nalp),MLT(ntime) + REAL*8 BLOCAL(ntime_max,Nalp),BMIN(ntime_max) + REAL*8 XJ(ntime_max,Nalp),MLT(ntime_max) REAL*8 Lm(ntime_max,Nalp),Lstar(ntime_max,Nalp) C COMMON /magmod/k_ext,k_l,kint @@ -285,7 +284,6 @@ SUBROUTINE Lstar_Phi1(ntime,whichinv,options,iyearsat, c IMPLICIT NONE INCLUDE 'variables.inc' - INCLUDE 'ntime_max.inc' C c declare inputs INTEGER*4 whichinv,options(5) @@ -752,7 +750,6 @@ SUBROUTINE GET_FIELD_MULTI(ntime,kext,options,sysaxes,iyearsat, c IMPLICIT NONE INCLUDE 'variables.inc' - INCLUDE 'ntime_max.inc' ! include file created by make, defines ntime_max C c declare inputs INTEGER*4 ntime From 069eb73ab171645a5d5a6dbe3cecee1d784695fb Mon Sep 17 00:00:00 2001 From: Antoine Brunet Date: Thu, 11 Apr 2024 15:32:30 +0200 Subject: [PATCH 3/5] Correction of typos in comments --- source/AE8_AP8.f | 2 +- source/AFRL_CRRES_models.f | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/source/AE8_AP8.f b/source/AE8_AP8.f index 99c16b56..c2afabae 100644 --- a/source/AE8_AP8.f +++ b/source/AE8_AP8.f @@ -51,7 +51,7 @@ ! 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), long array(ntime), double array(ntime)) +! 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]) diff --git a/source/AFRL_CRRES_models.f b/source/AFRL_CRRES_models.f index b0305e74..66d6e969 100644 --- a/source/AFRL_CRRES_models.f +++ b/source/AFRL_CRRES_models.f @@ -41,7 +41,7 @@ ! 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), long array(ntime), double array(ntime)) +! 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]) From 10137c75664af43a861ae12253ada3e31c83fddd Mon Sep 17 00:00:00 2001 From: Antoine Brunet Date: Tue, 16 Apr 2024 15:25:09 +0200 Subject: [PATCH 4/5] Updated documentation to reflect the possibility to use variable length arrays --- docs/frames/frames.html | 161 +++++++++++++++++++--------------------- source/CoordTrans.f | 6 +- source/LAndI2Lstar.f | 1 - source/get_bderivs.f | 2 - source/get_hemi.f | 1 - source/onera_desp_lib.f | 38 +++++----- 6 files changed, 95 insertions(+), 114 deletions(-) diff --git a/docs/frames/frames.html b/docs/frames/frames.html index 65d8c280..5758b0b6 100755 --- a/docs/frames/frames.html +++ b/docs/frames/frames.html @@ -34,14 +34,12 @@

Common Argument Definiti format. This section provides the definitions for those arguments. The particular definitions given here are for the functions that perform a particular calculation for multiple - points (ntime, up to NTIME_MAX, defined in - IRBEM_GET_NTIME_MAX()). For a few functions, there is also a + points (ntime). For a few functions, there is also a pitch-angle dimension given by Nipa, which is up to NALPHA_MAX = 25. For functions that can compute for only a single point or angle the NTIME_MAX argument can be ignored.

- ntime: long integer number of time in arrays - (max allowed is NTIME_MAX) + ntime: long integer number of time in arrays

nipa: long integer number of pitch angles in arrays @@ -237,27 +235,27 @@

Common Argument Definiti

- iyear, iyr or iyearsat: array(NTIME_MAX) of long integer year when + iyear, iyr or iyearsat: array(ntime) of long integer year when measurements are given.
- idoy: array(NTIME_MAX) of long integer doy when + idoy: array(ntime) of long integer doy when measurements are given
- UT or secs: array(NTIME_MAX) of double, UT in seconds + UT or secs: array(ntime) of double, UT in seconds
Note: In matlab, year/doy/UT arguments are replaced by the single matlabd argument, which is a Matlab date numbers (a double precision day counter created with datenum). They can be generated, if needed, using the helper function [iyear,idoy,UT] = onera_desp_lib_matlabd2yds(matlabd).

- x1: array(NTIME_MAX) of double, first coordinate + x1: array(ntime) of double, first coordinate according to sysaxes. If sysaxes is 0 then altitude has to be in km otherwise use dimensionless variables (in Re)

- x2: array(NTIME_MAX) of double, second coordinate + x2: array(ntime) of double, second coordinate according to sysaxes. If sysaxes is 0 then latitude has to be in degrees otherwise use dimensionless variables (in Re)

- x3: array(NTIME_MAX) of double, third coordinate + x3: array(ntime) of double, third coordinate according to sysaxes. If sysaxes is 0 then longitude has to be in degrees otherwise use dimensionless variables (in Re).

@@ -270,7 +268,7 @@

Common Argument Definiti distance along the drift orbit (use R0 < 1) for the drift loss cone.

- maginput: array (25,NTIME_MAX) of double to specify + maginput: array (25,ntime) of double to specify magnetic fields inputs such as:
@@ -368,7 +366,7 @@

Common Argument Definiti

- Lm: L McIlwain (array(NTIME_MAX,NALPHA_MAX) of double) + Lm: L McIlwain (array(ntime) of double, unless specified)

@@ -379,21 +377,18 @@

Common Argument Definiti

Lstar or Φ : L Roederer  or Φ=2π*Bo*/Lstar - [nT Re2] (array(NTIME_MAX,NALPHA_MAX) of double) + [nT Re2] (array(ntime) of double, unless specified)
- Blocal: magnitude of magnetic field at point (array(NTIME_MAX,NALPHA_MAX) - of double) - [nT] + Blocal: magnitude of magnetic field at point (array(ntime) of double, unless specified) - [nT]
- Bmir or Bmirror: magnitude of magnetic field at mirror point (array(NTIME_MAX,NALPHA_MAX) - of double) - [nT] + Bmir or Bmirror: magnitude of magnetic field at mirror point (array(ntime) of double, unless specified) - [nT]
Bmin: magnitude of magnetic field at equator - (array(NTIME_MAX,NALPHA_MAX) of double) - [nT] + (array(ntime) of double, unless specified) - [nT]
XJ: I, related to second adiabatic invariant - (array(NTIME_MAX,NALPHA_MAX) of double) - [Re]
- MLT: magnetic local time in hours, (array(NTIME_MAX,NALPHA_MAX) - of double) - [hour] + (array(ntime) of double, unless specified) - [Re]
+ MLT: magnetic local time in hours, (array(ntime) of double, unless specified) - [hour]

@@ -549,13 +544,6 @@

DESCRIPTION:
- - - - - - -

INPUTS:

@@ -605,11 +593,13 @@

DESCRIPTION: This function allows one to compute the magnetic coordinates at any s/c position and pitch-angle, i.e. L, L*, Bmirror, Bequator. A set of internal/external - field can be selected. + field can be selected. + All the output arrays have shape (NTIME_MAX, Nipa). +

INPUTS:

- ntime, Nipa, kext, options, sysaxes, iyear, idoy, UT, x1, x2, x3, alpha, maginput

+ ntime (up to NTIME_MAX), Nipa, kext, options, sysaxes, iyear, idoy, UT, x1, x2, x3, alpha, maginput

see Common Argument Definitions @@ -654,14 +644,16 @@

DESCRIPTION: + The errors in L* values are less than 2%. + All the output arrays have shape (NTIME_MAX, Nipa). +

INPUTS:

- ntime, Nipa, kext, options, sysaxes, iyear, idoy, UT, x1, x2, x3, alpha, maginput

+ ntime (up to NTIME_MAX), Nipa, kext, options, sysaxes, iyear, idoy, UT, x1, x2, x3, alpha, maginput

see Common Argument Definitions @@ -1041,9 +1033,9 @@

INPUTS:

OUTPUTS:

Bgeo: BxGEO,ByGEO and BzGEO of the magnetic field - (nT), (array(3,NTIME_MAX) of double) + (nT), (array(3,ntime) of double)
- Bl: magnitude of magnetic field (array(NTIME_MAX) of double) (nT) + Bl: magnitude of magnetic field (array(ntime) of double) (nT)
 

CALLING SEQUENCE FROM MATLAB: @@ -1104,13 +1096,13 @@

INPUTS:

OUTPUTS:

Bgeo: BxGEO,ByGEO and BzGEO of the magnetic field - (nT), (array(3,NTIME_MAX) of double) + (nT), (array(3,ntime) of double)
- Bmag: magnitude of magnetic field (array(NTIME_MAX) of double) (nT) + Bmag: magnitude of magnetic field (array(ntime) of double) (nT)
- gradBmag: gradient (GEO) of magnitude of magnetic field (array(3,NTIME_MAX) of double) (nT) + gradBmag: gradient (GEO) of magnitude of magnetic field (array(3,ntime) of double) (nT)
- diffB: derivatives of magnetic field vector (array(3,3,NTIME_MAX) of double) (nT)
+ diffB: derivatives of magnetic field vector (array(3,3,ntime) of double) (nT)
diffB(i,j,t) = dB_i/dx_j at t'th location. GEO coordinates.
  @@ -1149,14 +1141,14 @@

INPUTS:

OUTPUTS:

- grad_par: gradient of Bmag along B nT/RE, (array(NTIME_MAX) of double)
- grad_perp: gradient of Bmag perpendicular to B nT/RE (array(3,NTIME_MAX) of double)
- grad_drift: (bhat x grad_perp)/B, 1/RE (part of gradient drift velocity) (array(3,NTIME_MAX) of double)
- curvature: (bhat dot grad)bhat, 1/RE (part of curvature force) (array(3,NTIME_MAX) of double)
- Rcurv: 1/|curvature| RE (radius of curvature) (array(3,NTIME_MAX) of double)
- curv_drift: (bhat x curvature), 1/RE (part of curvature drift) (array(NTIME_MAX) of double)
- curlB: curl of B (nT/RE) (part of electrostatic current term) (array(3,NTIME_MAX) of double)
- divB: divergence of B (nT/RE) (should be zero!) (array(NTIME_MAX) of double)
+ grad_par: gradient of Bmag along B nT/RE, (array(ntime) of double)
+ grad_perp: gradient of Bmag perpendicular to B nT/RE (array(3,ntime) of double)
+ grad_drift: (bhat x grad_perp)/B, 1/RE (part of gradient drift velocity) (array(3,ntime) of double)
+ curvature: (bhat dot grad)bhat, 1/RE (part of curvature force) (array(3,ntime) of double)
+ Rcurv: 1/|curvature| RE (radius of curvature) (array(3,ntime) of double)
+ curv_drift: (bhat x curvature), 1/RE (part of curvature drift) (array(ntime) of double)
+ curlB: curl of B (nT/RE) (part of electrostatic current term) (array(3,ntime) of double)
+ divB: divergence of B (nT/RE) (should be zero!) (array(ntime) of double)

@@ -1442,7 +1434,7 @@

INPUTS or OUTPUTS (accor

Phi: Φ=2π*Bo/Lstar [nT Re2] - (array(GET_IRBEM_NTIME_MAX()) of double)  + (array(ntime) of double) 

CALLING SEQUENCE FROM MATLAB: @@ -1474,7 +1466,7 @@

COORD_TRANS_VECDESCRIPTION:

Generic coordinates transformation from one geophysical or - heliospheric coordinate system to another. Handles up to GET_IRBEM_NTIME_MAX() times at once. + heliospheric coordinate system to another.
 

INPUTS: @@ -1491,7 +1483,7 @@

INPUTS: sysaxesOUT : long integer indicating output coordinates system
- xIN : array (3,NTIME_MAX) of double with position + xIN : array (3,ntime) of double with position in input coordinates system. @@ -2149,7 +2141,7 @@

INPUTS:

OUTPUTS:

- xOUT : array (3,NTIME_MAX) of double with position + xOUT : array (3,ntime) of double with position in output coordinates system. (systems defined by sysaxesOUT)
  @@ -4320,8 +4312,7 @@

DESCRIPTION:INPUTS:

- Ntime: long integer number of time in arrays - (max allowed is GET_IRBEM_NTIME_MAX()) + Ntime: long integer number of time in arrays

WhichAp: long integer to select the kind @@ -4346,35 +4337,35 @@

INPUTS:

- DOY: array(GET_IRBEM_NTIME_MAX()) of long integer, Number of + DOY: array(ntime) of long integer, Number of day of year (DOY=1 is for January 1st)

- UT: array(GET_IRBEM_NTIME_MAX()) of double, UT in + UT: array(ntime) of double, UT in seconds when positions are given..

- Alt: array(GET_IRBEM_NTIME_MAX()) of double, Altitude + Alt: array(ntime) of double, Altitude in km (greater than 85 km).

- Lat: array(GET_IRBEM_NTIME_MAX()) of double, Geodetic + Lat: array(ntime) of double, Geodetic latitude (Deg.)

- Long: array(GET_IRBEM_NTIME_MAX()) of double, Geodetic + Long: array(ntime) of double, Geodetic longitude (Deg.)

- F107A: array(GET_IRBEM_NTIME_MAX()) of double, 3 month + F107A: array(ntime) of double, 3 month average of F10.7 flux.

- F107: array(GET_IRBEM_NTIME_MAX()) of double, daily + F107: array(ntime) of double, daily F10.7 flux for previous day.

- Ap: array(7,GET_IRBEM_NTIME_MAX()) of double where:
+ Ap: array(7,ntime) of double where:

1st element=Daily Ap
@@ -4393,7 +4384,7 @@

INPUTS:

OUTPUTS:

- Dens: array(8,GET_IRBEM_NTIME_MAX()) of double where:
+ Dens: array(8,ntime) of double where:

Dens(1st element,*)=He number density (cm-3) @@ -4407,7 +4398,7 @@

OUTPUTS: Dens(8th element,*)=N number density (cm-3)

- Temp: array(2,GET_IRBEM_NTIME_MAX()) of double where:
+ Temp: array(2,ntime) of double where:

Temp(1st element,*)=Exospheric temperature (K) @@ -4463,8 +4454,7 @@

DESCRIPTION:INPUTS:

- Ntime: long integer number of time in arrays - (max allowed is GET_IRBEM_NTIME_MAX()) + Ntime: long integer number of time in arrays

WhichAp: long integer to select the kind @@ -4485,35 +4475,35 @@

INPUTS:

- DOY: array(GET_IRBEM_NTIME_MAX()) of long integer, Number of + DOY: array(ntime) of long integer, Number of day of year (DOY=1 is for January 1st)

- UT: array(GET_IRBEM_NTIME_MAX()) of double, UT in + UT: array(ntime) of double, UT in seconds when positions are given..

- Alt: array(GET_IRBEM_NTIME_MAX()) of double, Altitude + Alt: array(ntime) of double, Altitude (km).

- Lat: array(GET_IRBEM_NTIME_MAX()) of double, Geodetic + Lat: array(ntime) of double, Geodetic latitude (Deg.)

- Long: array(GET_IRBEM_NTIME_MAX()) of double, Geodetic + Long: array(ntime) of double, Geodetic longitude (Deg.)

- F107A: array(GET_IRBEM_NTIME_MAX()) of double, 3 month + F107A: array(ntime) of double, 3 month average of F10.7 flux.

- F107: array(GET_IRBEM_NTIME_MAX()) of double, daily + F107: array(ntime) of double, daily F10.7 flux for previous day.

- Ap: array(7,GET_IRBEM_NTIME_MAX()) of double where:
+ Ap: array(7,ntime) of double where:

1st element=Daily Ap
@@ -4538,7 +4528,7 @@

INPUTS:

OUTPUTS:

- Dens: array(8,GET_IRBEM_NTIME_MAX()) of double where:
+ Dens: array(8,ntime) of double where:

Dens(1st element,*)=He number density (cm-3) @@ -4552,7 +4542,7 @@

INPUTS: Dens(8th element,*)=N number density (cm-3)

- Temp: array(2,GET_IRBEM_NTIME_MAX()) of double where:
+ Temp: array(2,ntime) of double where:

Temp(1st element,*)=Exospheric temperature (K) @@ -4612,8 +4602,7 @@

DESCRIPTION:INPUTS:

- Ntime: long integer number of time in arrays - (max allowed is GET_IRBEM_NTIME_MAX()) + Ntime: long integer number of time in arrays

WhichAp: long integer to select the kind @@ -4640,35 +4629,35 @@

INPUTS:

- DOY: array(GET_IRBEM_NTIME_MAX()) of long integer, Number of + DOY: array(ntime) of long integer, Number of day of year (DOY=1 is for January 1st)

- UT: array(GET_IRBEM_NTIME_MAX()) of double, UT in + UT: array(ntime) of double, UT in seconds when positions are given..

- Alt: array(GET_IRBEM_NTIME_MAX()) of double, Altitude + Alt: array(ntime) of double, Altitude in km (greater than 85 km).

- Lat: array(GET_IRBEM_NTIME_MAX()) of double, Geodetic + Lat: array(ntime) of double, Geodetic latitude (Deg.)

- Long: array(GET_IRBEM_NTIME_MAX()) of double, Geodetic + Long: array(ntime) of double, Geodetic longitude (Deg.)

- F107A: array(GET_IRBEM_NTIME_MAX()) of double, 3 month + F107A: array(ntime) of double, 3 month average of F10.7 flux.

- F107: array(GET_IRBEM_NTIME_MAX()) of double, daily + F107: array(ntime) of double, daily F10.7 flux for previous day.

- Ap: array(7,GET_IRBEM_NTIME_MAX()) of double where:
+ Ap: array(7,ntime) of double where:

1st element=Daily Ap
@@ -4692,7 +4681,7 @@

INPUTS:

OUTPUTS:

- Dens: array(9,GET_IRBEM_NTIME_MAX()) of double where:
+ Dens: array(9,ntime) of double where:

Dens(1st element,*)=He number density (cm-3) @@ -4708,7 +4697,7 @@

INPUTS: Dens(9th element,*)=Anomalous oxygen number density (cm-3)

- Temp: array(2,GET_IRBEM_NTIME_MAX()) of double where:
+ Temp: array(2,ntime) of double where:

Temp(1st element,*)=Exospheric temperature (K) diff --git a/source/CoordTrans.f b/source/CoordTrans.f index dd965ece..212f7f35 100644 --- a/source/CoordTrans.f +++ b/source/CoordTrans.f @@ -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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ @@ -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 @@ -2716,7 +2713,6 @@ 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 diff --git a/source/LAndI2Lstar.f b/source/LAndI2Lstar.f index 39d06e44..eec23286 100644 --- a/source/LAndI2Lstar.f +++ b/source/LAndI2Lstar.f @@ -20669,7 +20669,6 @@ 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) diff --git a/source/get_bderivs.f b/source/get_bderivs.f index 7cc6f48d..185f7257 100644 --- a/source/get_bderivs.f +++ b/source/get_bderivs.f @@ -30,7 +30,6 @@ SUBROUTINE GET_Bderivs(ntime,kext,options,sysaxes,dX, 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 @@ -133,7 +132,6 @@ 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 diff --git a/source/get_hemi.f b/source/get_hemi.f index 0b2af9d9..b8f4d33d 100644 --- a/source/get_hemi.f +++ b/source/get_hemi.f @@ -94,7 +94,6 @@ SUBROUTINE GET_HEMI_MULTI(ntime,kext,options,sysaxes,iyearsat, 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) diff --git a/source/onera_desp_lib.f b/source/onera_desp_lib.f index e1238bf2..7250c077 100644 --- a/source/onera_desp_lib.f +++ b/source/onera_desp_lib.f @@ -172,12 +172,12 @@ SUBROUTINE make_lstar_shell_splitting1(ntime,Nipa,kext,options, INTEGER*4 kext,k_ext,k_l,options(5),Nalp,Nipa PARAMETER (Nalp=25) 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) + INTEGER*4 iyearsat(ntime) + integer*4 idoy(ntime) + real*8 UT(ntime) + real*8 xIN1(ntime),xIN2(ntime),xIN3(ntime) real*8 alpha(Nalp) - real*8 maginput(25,ntime_max) + real*8 maginput(25,ntime) c c c Declare internal variables @@ -190,8 +190,8 @@ SUBROUTINE make_lstar_shell_splitting1(ntime,Nipa,kext,options, c c Declare output variables REAL*8 BLOCAL(ntime_max,Nalp),BMIN(ntime_max) - REAL*8 XJ(ntime_max,Nalp),MLT(ntime_max) - REAL*8 Lm(ntime_max,Nalp),Lstar(ntime_max,Nalp) + REAL*8 XJ(ntime_max,Nipa),MLT(ntime_max) + REAL*8 Lm(ntime_max,Nipa),Lstar(ntime_max,Nipa) C COMMON /magmod/k_ext,k_l,kint COMMON /flag_L/Ilflag @@ -1561,11 +1561,11 @@ subroutine msis86(ntime,whichAp,DOY,UT,ALT,LAT,LONG,F107A, &F107,AP,Dens,Temp) c IMPLICIT NONE - INTEGER*4 ISW,ntime,whichAp,DOY(100000),I,J + INTEGER*4 ISW,ntime,whichAp,DOY(ntime),I,J REAL*8 SV(25),STL - REAL*8 UT(100000),ALT(100000),LAT(100000),LONG(100000) - REAL*8 F107A(100000),F107(100000),AP(7,100000) - REAL*8 Dens(8,100000),Temp(2,100000),APin(7),D(8),T(2) + REAL*8 UT(ntime),ALT(ntime),LAT(ntime),LONG(ntime) + REAL*8 F107A(ntime),F107(ntime),AP(7,ntime) + REAL*8 Dens(8,ntime),Temp(2,ntime),APin(7),D(8),T(2) c COMMON/CSWI/ISW DO I=1,25 @@ -1605,11 +1605,11 @@ subroutine msise90(ntime,whichAp,DOY,UT,ALT,LAT,LONG,F107A, &F107,AP,Dens,Temp) c IMPLICIT NONE - INTEGER*4 ISW,ntime,whichAp,DOY(100000),I,J + INTEGER*4 ISW,ntime,whichAp,DOY(ntime),I,J REAL*8 SV(25),STL - REAL*8 UT(100000),ALT(100000),LAT(100000),LONG(100000) - REAL*8 F107A(100000),F107(100000),AP(7,100000) - REAL*8 Dens(8,100000),Temp(2,100000),APin(7),D(8),T(2) + REAL*8 UT(ntime),ALT(ntime),LAT(ntime),LONG(ntime) + REAL*8 F107A(ntime),F107(ntime),AP(7,ntime) + REAL*8 Dens(8,ntime),Temp(2,ntime),APin(7),D(8),T(2) c COMMON/CSWI/ISW DO I=1,25 @@ -1639,11 +1639,11 @@ subroutine nrlmsise00(ntime,whichAp,DOY,UT,ALT,LAT,LONG,F107A, &F107,AP,Dens,Temp) c IMPLICIT NONE - INTEGER*4 ISW,ntime,whichAp,DOY(100000),I,J + INTEGER*4 ISW,ntime,whichAp,DOY(ntime),I,J REAL*8 SV(25),STL - REAL*8 UT(100000),ALT(100000),LAT(100000),LONG(100000) - REAL*8 F107A(100000),F107(100000),AP(7,100000) - REAL*8 Dens(9,100000),Temp(2,100000),APin(7),D(8),T(2) + REAL*8 UT(ntime),ALT(ntime),LAT(ntime),LONG(ntime) + REAL*8 F107A(ntime),F107(ntime),AP(7,ntime) + REAL*8 Dens(9,ntime),Temp(2,ntime),APin(7),D(8),T(2) c COMMON/CSWI/ISW DO I=1,25 From 3400ec8e70cf37f88bea65cacfc1e2157fbdc852 Mon Sep 17 00:00:00 2001 From: Antoine Brunet Date: Tue, 16 Apr 2024 15:35:40 +0200 Subject: [PATCH 5/5] Uniform treatment of Nene in radiation models, using variable array dimension everywhere, and updated documentation --- docs/frames/frames.html | 24 ++++++++++++------------ source/AE8_AP8.f | 18 ++++++++---------- source/AFRL_CRRES_models.f | 10 ++++------ 3 files changed, 24 insertions(+), 28 deletions(-) diff --git a/docs/frames/frames.html b/docs/frames/frames.html index 5758b0b6..2f9cc78c 100755 --- a/docs/frames/frames.html +++ b/docs/frames/frames.html @@ -3546,10 +3546,10 @@

INPUTS:

Nene: long integer number of energies - in array energy (max allowed is 25) + in array energy

- energy: array(2,25) of double. If whatf=1 or + energy: array(2,Nene) of double. If whatf=1 or 3 then energy(2,*) is not considered.
       Note: if energy is in MeV then differential flux will be per MeV. @@ -3588,7 +3588,7 @@

INPUTS:

OUTPUTS:

flux: flux according to selection of whatf for - all times and energies (array(GET_IRBEM_NTIME_MAX(),25) of double) + all times and energies (array(GET_IRBEM_NTIME_MAX(),Nene) of double)

CALLING SEQUENCE FROM MATLAB:

@@ -3661,10 +3661,10 @@

INPUTS:

Nene: long integer number of energies - in array energy (max allowed is 25) + in array energy

- energy: array(2,25) of double. If whatf=1 or + energy: array(2,Nene) of double. If whatf=1 or 3 then energy(2,*) is not considered.
       Note: if energy is in MeV then differential flux will be per MeV. @@ -3691,7 +3691,7 @@

OUTPUTS:

flux: flux according to selection of whatf for - all times and energies (array(GET_IRBEM_NTIME_MAX(),25) of double) + all times and energies (array(GET_IRBEM_NTIME_MAX(),Nene) of double)

CALLING SEQUENCE FROM MATLAB: @@ -3793,11 +3793,11 @@

INPUTS:

Nene: long integer number of energies in - array energy (max allowed is 25) + array energy

- energy: array(2,25) of double (MeV). If whatf=1 + energy: array(2,Nene) of double (MeV). If whatf=1 or 3 then energy(2,*) is not considered.

@@ -3850,7 +3850,7 @@

INPUTS:

OUTPUTS:

flux: flux according to selection of whatf for - all times and energies (array(GET_IRBEM_NTIME_MAX(),25) of double) + all times and energies (array(GET_IRBEM_NTIME_MAX(),Nene) of double)

CALLING SEQUENCE FROM MATLAB:

@@ -3926,12 +3926,12 @@

INPUTS:

Nene: long integer number of energies in - array energy (max allowed is 25) + array energy (max allowed is Nene)

- energy: array(2,25) of double (MeV). If whatf=1 + energy: array(2,Nene) of double (MeV). If whatf=1 or 3 then energy(2,*) is not considered.

@@ -3970,7 +3970,7 @@

OUTPUTS:

flux: flux according to selection of whatf for - all times and energies (array(GET_IRBEM_NTIME_MAX(),25) of double) + all times and energies (array(GET_IRBEM_NTIME_MAX(),Nene) of double)

CALLING SEQUENCE FROM MATLAB: diff --git a/source/AE8_AP8.f b/source/AE8_AP8.f index c2afabae..25b367c1 100644 --- a/source/AE8_AP8.f +++ b/source/AE8_AP8.f @@ -57,7 +57,7 @@ ! 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 @@ -79,11 +79,9 @@ 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),idoy(ntime) - REAL*8 energy(2,nene_max) + REAL*8 energy(2,nene) REAL*8 UT(ntime) real*8 xIN1(ntime),xIN2(ntime),xIN3(ntime) c Declare internal variables @@ -98,7 +96,7 @@ SUBROUTINE fly_in_nasa_aeap1(ntime,sysaxes,whichm,whatf,nene, 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 @@ -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 @@ -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 @@ -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 diff --git a/source/AFRL_CRRES_models.f b/source/AFRL_CRRES_models.f index 66d6e969..33f352ce 100644 --- a/source/AFRL_CRRES_models.f +++ b/source/AFRL_CRRES_models.f @@ -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]) +! 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) !--------------------------------------------------------------------------------------------------- @@ -62,11 +62,9 @@ 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),idoy(ntime) - REAL*8 energy(2,nene_max) + REAL*8 energy(2,nene) REAL*8 UT(ntime) real*8 xIN1(ntime),xIN2(ntime),xIN3(ntime) c Declare internal variables @@ -83,7 +81,7 @@ SUBROUTINE fly_in_afrl_crres1(ntime,sysaxes,whichm,whatf,nene, 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