|
|
C Written by Shunrong Zhang
|
|
|
C Imported into Madrigal from modelRecovery/fullAnalytic/
|
|
|
C isrim.f on 10/21/2005 - Was revision 1.6
|
|
|
C
|
|
|
C Modified to change COEFDIR to Madrigal path, added comment
|
|
|
C on 4 = ion drift
|
|
|
C
|
|
|
C $Id: isrim.f 3375 2011-05-11 17:00:38Z brideout $
|
|
|
C
|
|
|
C
|
|
|
SUBROUTINE ISRIM(KINST,DOY,SLT,ALT,GLAT,F107,AP,IPAR,OUTPUT)
|
|
|
C
|
|
|
C INCOHERENT SCATTER RADAR IONOSPHERIC MODEL (ISRIM)
|
|
|
C
|
|
|
C ISRIMs have been developed from long-term datasets of seven
|
|
|
C incoherent scatter radars spanning invariant latitudes from
|
|
|
C 25 to 75 degrees in American, European and Asian longitudes
|
|
|
C at Svalbard, Tromso, Sondrestrom, Millstone Hill, St. Santin,
|
|
|
C Arecibo and Shigaraki. These models represent electron density,
|
|
|
C ion and electron temperatures, and ion drifts in the E and F
|
|
|
C regions, giving a comprehensive quantitative description of
|
|
|
C ionospheric properties.
|
|
|
C
|
|
|
C REFERENCES
|
|
|
C
|
|
|
C Zhang, S.-R., J. M. Holt, A. P. van Eyken, M. McCready, C.
|
|
|
C Amory-Mazaudier, S. Fukao, and M. Sulzer, Ionospheric local model
|
|
|
C and climatology from long-term databases of multiple incoherent
|
|
|
C scatter radars, Geophys. Res. Lett. 32, L20102, doi:10.1029/
|
|
|
C 2005GL023603, 2005.
|
|
|
C
|
|
|
C Zhang, S.-R., J. M. Holt, A. M. Zalucha, and C. Amory-Mazaudier,
|
|
|
C Mid-latitude ionospheric plasma temperature climatology and
|
|
|
C empirical model based on Saint Santin incoherent scatter radar data
|
|
|
C from 1966-1987, J. Geophys. Res., 109, A11311,
|
|
|
C doi:10.1029/2004JA010709, 2004.
|
|
|
C
|
|
|
C Holt, J. M., S.-R. Zhang, and M. J. Buonsanto, Regional and local
|
|
|
C ionospheric models based on Millstone Hill incoherent scatter radar
|
|
|
C data, Geophys. Res. Lett.,, 29(8), 10.1029/2002GL014678, 2002.
|
|
|
C
|
|
|
C
|
|
|
C ONLINE VERSION
|
|
|
C
|
|
|
C http://madrigal.haystack.mit.edu/models/
|
|
|
C
|
|
|
C
|
|
|
C CONTACT INFO
|
|
|
C
|
|
|
C Dr. Shun-Rong Zhang (shunrong@haystack.mit.edu)
|
|
|
C Dr. John M. Holt
|
|
|
C MIT Haystack Observatory
|
|
|
C Route 40,
|
|
|
C Westford, MA 01886
|
|
|
C USA
|
|
|
C PHONE: 781-981-5725; FAX: 781-981-5766
|
|
|
C
|
|
|
C
|
|
|
C REVISIONS
|
|
|
C
|
|
|
C Updated by SRZ, MAY 2007 --
|
|
|
C ALLOW FOR KINST = 30 FOR MH LOCAL MODEL
|
|
|
C
|
|
|
C Updated by SRZ, October 2005 --
|
|
|
C READ IN ALL COEFFICIENT FILES, INSTEAD OF READING IN
|
|
|
C AS NEEDED FOR A SPECIFIC PARAMETER. THIS CHANGE IS TO ACCOMMODATE
|
|
|
C THE MADRIGAL APPLICATIONS TO IMPROVE EFFICIENCY WHEN MULTIPLE
|
|
|
C PARAMETERS ARE COMPUTERED.
|
|
|
C
|
|
|
C Updated by SRZ, May 2005 ---
|
|
|
C ADD KINST (INSTRUMENT CODE); REMOVE THE FUNCTION OF OUTPUT(5)
|
|
|
C
|
|
|
C FIRST Version by SRZ, July 2002 ---
|
|
|
C MILLSTONE HILL IONOSPHERIC MODEL
|
|
|
C THIS SUBROUTINE RETURNS IONOSPHERIC PARAMETERS (NE, TE, AND TI)
|
|
|
C BASED ON MILLSTONE HILL INCOHERENT SCATTER RADAR OBSERVATION MADE
|
|
|
C DURING 1976-2001.
|
|
|
C
|
|
|
C
|
|
|
C INPUTS
|
|
|
C
|
|
|
C KINST - STATION/INSTRUMENT CODE
|
|
|
C 31 -> MILLSTONE HILL STEERABLE (regional model)
|
|
|
C 32 -> MILLSTONE HILL ZENITH (local model)
|
|
|
C 30 -> MILLSTONE HILL ZENITH (local model)
|
|
|
C !! 10 -> JICAMARCA - NOT AVAILABLE!!
|
|
|
C 20 -> ARECIBO
|
|
|
C 25 -> MU RADAR
|
|
|
C 40 -> ST. SANTIN
|
|
|
C 72 -> EISCAT TROMSO
|
|
|
C 80 -> SONDRESTROM (local)
|
|
|
C 95 -> EISCAT SVALBARD
|
|
|
C 0 -> AMERICA REGIONAL MODEL (GEODETC 15-70 DEGREES)
|
|
|
C 1 -> WESTERN eUROPEAN MODEL (GEODETIC 45-78 DEGREES)
|
|
|
C DOY: DAY NUMBER OF YEAR; DOUBLE PRECISION SCALAR
|
|
|
C [ 1 - 365 ] (always = an int value)
|
|
|
C SLT: LOCAL TIME (HOUR); DOUBLE PRECISION SCALAR
|
|
|
C [ 0 - 24]
|
|
|
C ALT: GEODETIC ALTITUDE (KM); DOUBLE PRECISION SCALAR
|
|
|
C LOCAL MODEL: [100 - 1000]; for Millstone Hill
|
|
|
C LOCAL MODEL: [200 - 500]; for Shigaraki
|
|
|
C LOCAL MODEL: [100 - 600]; for other sites
|
|
|
C REGIONAL MODEL: [100 - 600]: America/Europe
|
|
|
C REGIONAL MODEL: [200 - 600]: Millstone Hill
|
|
|
C GLAT: (ALSO AS OUTPUT IF LOCAL MODEL IS CALLED)
|
|
|
C GEODETIC LATITUDE DEGREE, DOUBLE PRECISION SCALAR
|
|
|
C NEEDED ONLY FOR REGIONAL MODELS
|
|
|
C F107: SOLAR 10.7 CM FLUX FOR THE PREVIOUS DAY;
|
|
|
C DOUBLE PRECISION SCALAR
|
|
|
C AP: 3-HOUR AP INDEX AT 3 HOURS BEFORE THE CURRENT TIME;
|
|
|
C DOUBLE PRECISION SCALAR
|
|
|
C IPAR: IONOSPHERIC PARAMETER TO BE RETURNED;
|
|
|
C INTEGER SCALAR
|
|
|
C = 1, FOR NE
|
|
|
C = 2, FOR TE
|
|
|
C = 3, FOR TI
|
|
|
C = 4, FOR PARALLEL ION DRIFT, + UPWARD
|
|
|
C OUTPUT(4): FLAG FOR HMAX AND NMAX OUTPUT;
|
|
|
C DOUBLE PRECISION SCALAR
|
|
|
C SET OUTPUT(4) TO ANY VALUE >0, THEN OUTPUT(2) AND
|
|
|
C OUTPUT(3) RETURN HMAX AND NMAX.
|
|
|
C
|
|
|
C OUTPUT
|
|
|
C
|
|
|
C OUTPUT: OUTPUT DATA; DOUBLE PRECISION ARRAY, 10 DIMS
|
|
|
C OUTPUT(1): OUTPUT DATA; NE 1/M^3 FOR IPAR=1
|
|
|
C TE, K FOR IPAR=2
|
|
|
C TI, K FOR IPAR=3
|
|
|
C OUTPUT(2): HMAX, KM FOR IPAR=1 AND OUTPUT(4)>0
|
|
|
C OUTPUT(3): NMAX, 1/M^3 FOR IPAR=1 AND OUTPUT(4)>0
|
|
|
C OUTPUT(4): MODEL DIAGNOSTIC STATUS
|
|
|
C = 0 NORMAL
|
|
|
C = -1 STATION NOT INCLUDED
|
|
|
C = -2 PARAMETER NOT INCLUDED
|
|
|
C = -3 STATION AND PARAMETER NOT MATCH
|
|
|
C GLAT: GEODETIC LATITUDE DEGREE, DOUBLE PRECISION SCALAR
|
|
|
C ONLY WHEN A LOCAL MODEL IS COMPUTED
|
|
|
C
|
|
|
C .. Scalar Arguments ..
|
|
|
DOUBLE PRECISION ALT,AP,DOY,F107,GLAT,SLT
|
|
|
INTEGER IPAR,KINST
|
|
|
C ..
|
|
|
C .. Array Arguments ..
|
|
|
DOUBLE PRECISION OUTPUT(10)
|
|
|
C ..
|
|
|
C .. Local Scalars ..
|
|
|
DOUBLE PRECISION DMAX,HMAX
|
|
|
INTEGER I,IAREA,IIPAR,IKINST,IKINST3,IKINST4,ISN,J,NCALLED,IND
|
|
|
INTEGER LIAREA
|
|
|
C LIAREA - last IAREA value
|
|
|
CHARACTER*2 CPAR
|
|
|
CHARACTER*128 COEFDIR,COEFILE,MADDIR
|
|
|
C ..
|
|
|
C .. Local Arrays ..
|
|
|
DOUBLE PRECISION GLATS(11),DWK1,DWK2
|
|
|
INTEGER STA(4),STATUS3(4),STATUS4(4)
|
|
|
CHARACTER*2 SNAME(11)
|
|
|
C ..
|
|
|
C .. External Subroutines ..
|
|
|
EXTERNAL PARAMETER3D,PARAMETER4D,READCOEF,SETKNOTS3D,
|
|
|
* SETKNOTS4D
|
|
|
C ..
|
|
|
C .. Save statement ..
|
|
|
SAVE NCALLED,IKINST,IKINST3,IKINST4,IIPAR,LIAREA
|
|
|
SAVE STATUS3,STATUS4
|
|
|
C .. IKINST is the kinst of the last instrument
|
|
|
C .. IKINST3 is the kinst of the last local instrument
|
|
|
C .. IKINST4 is the kinst of the last regional instrument
|
|
|
C .. Coeff from both regional and local are cached
|
|
|
C ..
|
|
|
C .. Intrinsic Functions ..
|
|
|
INTRINSIC INDEX
|
|
|
C ..
|
|
|
C .. Data statements ..
|
|
|
C
|
|
|
DATA SNAME/'AR','MU','MH','ST','TR','SO','ES','RM','EU','AM','M2'/
|
|
|
DATA GLATS/18D0,35D0,43D0,45D0,70D0,67D0,79D0,43D0,67D0,0D0,43D0/
|
|
|
DATA IKINST3/999/,IKINST4/999/,IIPAR/999/,LIAREA/999/
|
|
|
C ..
|
|
|
C .. DIRECTORY OF COEFFICIENT FILES ...............
|
|
|
C COEFDIR='coefiles/'
|
|
|
C this line modified by B. Rideout
|
|
|
C path should be automatically modified during installation
|
|
|
MADDIR = 'MADROOT'
|
|
|
IND = index(MADDIR, ' ') - 1
|
|
|
COEFDIR = MADDIR(1:IND) // '/source/madf/geolib/3D4DCoef/'
|
|
|
C COEFDIR='/opt/madrigal/empiricalModels/3D4DCoef/'
|
|
|
C ..................................................
|
|
|
C
|
|
|
C IAREA determines whether the local (2) or regional (1) used
|
|
|
IAREA = 2
|
|
|
C OUTPUT(4) = 0D0
|
|
|
IF (KINST.EQ.30) KINST = 32
|
|
|
C add other instruments near Millstone
|
|
|
IF (KINST.EQ.5340) KINST = 32
|
|
|
IF (KINST.EQ.5360) KINST = 32
|
|
|
IF (AP.GT.100) AP = 100D0
|
|
|
IF (KINST.EQ.31 .AND. OUTPUT(4) .LE. 0.0D0) THEN
|
|
|
IAREA = 1
|
|
|
ELSE IF (KINST.EQ.31 .AND. OUTPUT(4) .GT. 0.0D0) THEN
|
|
|
KINST = 32
|
|
|
END IF
|
|
|
IF (KINST.EQ.0) IAREA = 1
|
|
|
IF (KINST.EQ.1) IAREA = 1
|
|
|
C determine the cofficient file name to be read
|
|
|
CPAR = 'NO'
|
|
|
IF (IPAR.EQ.1) CPAR = 'Ne'
|
|
|
IF (IPAR.EQ.2) CPAR = 'Te'
|
|
|
IF (IPAR.EQ.3) CPAR = 'Ti'
|
|
|
IF (IPAR.EQ.4) CPAR = 'Vo'
|
|
|
ISN = 999
|
|
|
IF (KINST.EQ.20) ISN = 1
|
|
|
IF (KINST.EQ.25) ISN = 2
|
|
|
IF (KINST.EQ.32) ISN = 3
|
|
|
IF (KINST.EQ.40) ISN = 4
|
|
|
IF (KINST.EQ.72) ISN = 5
|
|
|
IF (KINST.EQ.80) ISN = 6
|
|
|
IF (KINST.EQ.95) ISN = 7
|
|
|
IF (KINST.EQ.31 .AND. OUTPUT(4) .LE. 0.0D0) ISN = 8
|
|
|
IF (KINST.EQ.31 .AND. OUTPUT(4) .GT. 0.0D0) ISN = 3
|
|
|
IF (KINST.EQ.1) ISN = 9
|
|
|
IF (KINST.EQ.0) ISN = 10
|
|
|
IF (KINST.EQ.33) ISN = 11
|
|
|
C .. invalid inputs ..
|
|
|
IF (ISN.EQ.999) OUTPUT(4) = -1D0
|
|
|
IF (CPAR.EQ.'NO') OUTPUT(4) = -2D0
|
|
|
IF (KINST.EQ.25 .AND. CPAR.EQ.'Vo') OUTPUT(4) = -3D0
|
|
|
IF (KINST.EQ.40 .AND. CPAR.EQ.'Vo') OUTPUT(4) = -3D0
|
|
|
IF (KINST.EQ.72 .AND. CPAR.EQ.'Vo') OUTPUT(4) = -3D0
|
|
|
IF (KINST.EQ.95 .AND. CPAR.EQ.'Vo') OUTPUT(4) = -3D0
|
|
|
IF (KINST.EQ.1 .AND. CPAR.EQ.'Vo') OUTPUT(4) = -3D0
|
|
|
IF (KINST.EQ.0 .AND. CPAR.EQ.'Vo') OUTPUT(4) = -3D0
|
|
|
IF (KINST.EQ.31 .AND. CPAR.EQ.'Vo') OUTPUT(4) = -3D0
|
|
|
IF (KINST.EQ.33 .AND. CPAR.EQ.'Vo') OUTPUT(4) = -3D0
|
|
|
C
|
|
|
C Altitude rules are unneeded if just want hmax, nmax
|
|
|
IF (OUTPUT(4) .LE. 0.0D0) THEN
|
|
|
IF (ALT.LT.100D0.OR.ALT.GT.1000D0) OUTPUT(4) = -4D0
|
|
|
C IF (ALT.LT.200D0.AND.IAREA.EQ.1) OUTPUT(4) = -4D0
|
|
|
IF (ALT.LT.200D0.AND.KINST.EQ.31) OUTPUT(4) = -4D0
|
|
|
IF (ALT.GT.600D0.AND.IAREA.EQ.1) OUTPUT(4) = -4D0
|
|
|
IF (ALT.LT.200D0.AND.KINST.EQ.25) OUTPUT(4) = -4D0
|
|
|
IF (ALT.GT.450D0.AND.KINST.EQ.25) OUTPUT(4) = -4D0
|
|
|
ELSE
|
|
|
ALT=300.0D0
|
|
|
END IF
|
|
|
|
|
|
IF (OUTPUT(4).LT.0D0) THEN
|
|
|
OUTPUT(1) = -999D0
|
|
|
RETURN
|
|
|
END IF
|
|
|
C
|
|
|
IF (IAREA.NE.1) GLAT = GLATS(ISN)
|
|
|
COEFILE = COEFDIR(1:INDEX(COEFDIR,' ')-1)//SNAME(ISN)
|
|
|
IF (IAREA.EQ.1) THEN
|
|
|
COEFILE = COEFDIR(1:INDEX(COEFDIR,' ')-1)//SNAME(ISN)
|
|
|
END IF
|
|
|
C to check if the coefficient file has been read into common block
|
|
|
IF (IAREA .NE. 1) THEN
|
|
|
IF (IKINST3.NE.KINST) NCALLED = 0
|
|
|
ELSE
|
|
|
IF (IKINST4.NE.KINST) NCALLED = 0
|
|
|
END IF
|
|
|
IF (NCALLED.NE.12345) THEN
|
|
|
IF (IAREA .NE. 1) THEN
|
|
|
CALL READCOEF(0,COEFILE,IAREA,STATUS3)
|
|
|
ELSE
|
|
|
CALL READCOEF(0,COEFILE,IAREA,STATUS4)
|
|
|
END IF
|
|
|
NCALLED = 12345
|
|
|
END IF
|
|
|
IF (IAREA .NE. 1) THEN
|
|
|
DO I = 1,4
|
|
|
IF (STATUS3(I).NE.0 .AND. IPAR.EQ.I) THEN
|
|
|
OUTPUT(1) = -999D0
|
|
|
OUTPUT(4) = -3D0
|
|
|
RETURN
|
|
|
END IF
|
|
|
END DO
|
|
|
ELSE
|
|
|
DO I = 1,4
|
|
|
IF (STATUS4(I).NE.0 .AND. IPAR.EQ.I) THEN
|
|
|
OUTPUT(1) = -999D0
|
|
|
OUTPUT(4) = -3D0
|
|
|
RETURN
|
|
|
END IF
|
|
|
END DO
|
|
|
END IF
|
|
|
IF (IAREA .NE. 1) THEN
|
|
|
IF (IIPAR.NE.IPAR .OR. KINST.NE.IKINST3 .OR.
|
|
|
* LIAREA.NE.IAREA) THEN
|
|
|
CALL READCOEF(IPAR,COEFILE,IAREA,STA)
|
|
|
CALL SETKNOTS3D
|
|
|
IKINST3 = KINST
|
|
|
END IF
|
|
|
ELSE
|
|
|
IF (IIPAR.NE.IPAR .OR. KINST.NE.IKINST4 .OR.
|
|
|
* LIAREA.NE.IAREA) THEN
|
|
|
CALL READCOEF(IPAR,COEFILE,IAREA,STA)
|
|
|
CALL SETKNOTS4D
|
|
|
IKINST4 = KINST
|
|
|
END IF
|
|
|
END IF
|
|
|
IKINST = KINST
|
|
|
LIAREA = IAREA
|
|
|
IIPAR = IPAR
|
|
|
|
|
|
IF (IAREA.EQ.2) CALL PARAMETER3D(DOY,SLT,GLAT,ALT,F107,AP,
|
|
|
* OUTPUT(1))
|
|
|
IF (IAREA.EQ.1) CALL PARAMETER4D(DOY,SLT,GLAT,ALT,F107,AP,
|
|
|
* OUTPUT(1))
|
|
|
IF (IPAR.EQ.1 .AND. OUTPUT(1).LE.0D0) OUTPUT(4) = -5D0
|
|
|
IF (IPAR.EQ.1 .AND. OUTPUT(1).LE.0D0) OUTPUT(1) = -OUTPUT(1)
|
|
|
|
|
|
IF (IPAR.EQ.1 .AND. KINST.EQ.31) OUTPUT(1) = OUTPUT(1)*1D12
|
|
|
IF (IPAR.EQ.1 .AND. KINST.NE.31) OUTPUT(1) = 10**OUTPUT(1)
|
|
|
C find hmax and Nmax at 1 km resolution - only use local
|
|
|
IF (IPAR.EQ.1 .AND. OUTPUT(4).GT.0.d0) THEN
|
|
|
DWK1 = 0.0D0
|
|
|
DWK2 = 0.0D0
|
|
|
HMAX = 0.0D0
|
|
|
DMAX = -10.0D6
|
|
|
IF (SLT.LT.5.OR.SLT.GT.22) THEN
|
|
|
IEND = 250
|
|
|
ELSE
|
|
|
IEND = 200
|
|
|
END IF
|
|
|
C first loop - rough search in 20 km steps - look for decrease
|
|
|
DO I = 500,IEND,-20
|
|
|
DWK1 = I*1.d0
|
|
|
CALL PARAMETER3D(DOY,SLT,GLAT,DWK1,F107,AP,DWK2)
|
|
|
IF (DWK2 .LT. DMAX .AND. DWK2 .GT. 0.0D0) THEN
|
|
|
HMAX = DWK1 + 20.0D0
|
|
|
GOTO 100
|
|
|
ELSE
|
|
|
DMAX = DWK2
|
|
|
END IF
|
|
|
END DO
|
|
|
C if we reach here no peak was found, set error values
|
|
|
OUTPUT(2) = 0.0D0
|
|
|
OUTPUT(3) = 0.0D0
|
|
|
OUTPUT(4) = -999.D0
|
|
|
RETURN
|
|
|
C second loop - fine search in 1 km steps
|
|
|
100 DO I = INT(HMAX)-15,INT(HMAX)+15,1
|
|
|
DWK1 = I*1.d0
|
|
|
CALL PARAMETER3D(DOY,SLT,GLAT,DWK1,F107,AP,DWK2)
|
|
|
IF (DWK2 .GT. DMAX) THEN
|
|
|
DMAX = DWK2
|
|
|
HMAX = DWK1
|
|
|
END IF
|
|
|
END DO
|
|
|
OUTPUT(2) = HMAX
|
|
|
OUTPUT(3) = 10**DMAX
|
|
|
OUTPUT(4) = -999.D0
|
|
|
END IF
|
|
|
RETURN
|
|
|
END
|
|
|
C
|
|
|
C -------------------------------------------------------------
|
|
|
C
|
|
|
SUBROUTINE READCOEF(SW,COEFILEPATH,MODEL,STATUS)
|
|
|
C
|
|
|
C INPUT: SW, COEFILEPATH,IAREA
|
|
|
C OUTPUT: STATUS(4)
|
|
|
C SW: SW=0 IF ONLY READ COEFFICIENTS FROM FILES
|
|
|
C SW.NE.0 IF ONLY LOAD COEFFICIENTS INTO COMMON BLOCKS
|
|
|
C FOR THE SPECIFIC PARAMETER; I.E.,
|
|
|
C SW=1, NE
|
|
|
C SW=2, TE
|
|
|
C SW=3, TI
|
|
|
C SW=4, VO
|
|
|
C MODEL - 1 if regional model, 2 if local model
|
|
|
C
|
|
|
C .. Scalar Arguments ..
|
|
|
INTEGER SW,MODEL
|
|
|
CHARACTER*(*) COEFILEPATH
|
|
|
C ..
|
|
|
C .. Local Scalars ..
|
|
|
INTEGER I,IAREA,IC,IPAR,J
|
|
|
CHARACTER*2 CPAR
|
|
|
CHARACTER*256 COEFILE
|
|
|
C ..
|
|
|
C .. Common blocks ..
|
|
|
|
|
|
COMMON /BASINFO1/KKS,NX,NY,NS,NG
|
|
|
COMMON /BASINFO13/KKS3,NX3,NY3,NS3,NG3
|
|
|
COMMON /BASINFO14/KKS4,NX4,NY4,NS4,NG4
|
|
|
COMMON /BASINFO2/XA,XB,YA,YB,SA,SB,GA,GB
|
|
|
COMMON /BASINFO23/XA3,XB3,YA3,YB3,SA3,SB3,GA3,GB3
|
|
|
COMMON /BASINFO24/XA4,XB4,YA4,YB4,SA4,SB4,GA4,GB4
|
|
|
COMMON /COEF/COEFS,N,M,ICC
|
|
|
COMMON /COEF3/COEFS3,N3,M3,ICC3
|
|
|
COMMON /COEF4/COEFS4,N4,M4,ICC4
|
|
|
COMMON /NODES/TT,WL,NBKL
|
|
|
COMMON /NODES3/TT3,WL3,NBKL3
|
|
|
COMMON /NODES4/TT4,WL4,NBKL4
|
|
|
DOUBLE PRECISION GA,GB,SA,SB,XA,XB,YA,YB
|
|
|
INTEGER ICC,KKS,M,N,NBKL,NG,NS,NX,NY
|
|
|
DOUBLE PRECISION COEFS(10,5000),COEFS4(4,10,5000),GA4(4),GB4(4),
|
|
|
* SA4(4),SB4(4),TT(100),TT4(4,100),WL(50),
|
|
|
* WL4(4,50),XA4(4),XB4(4),YA4(4),YB4(4)
|
|
|
DOUBLE PRECISION COEFS3(4,10,5000),GA3(4),GB3(4),
|
|
|
* SA3(4),SB3(4),TT3(4,100),
|
|
|
* WL3(4,50),XA3(4),XB3(4),YA3(4),YB3(4)
|
|
|
INTEGER ICC4(4),KKS4(4),M4(4),N4(4),NBKL4(4),NG4(4),NS4(4),NX4(4),
|
|
|
* NY4(4)
|
|
|
INTEGER ICC3(4),KKS3(4),M3(4),N3(4),NBKL3(4),NG3(4),NS3(4),NX3(4),
|
|
|
* NY3(4)
|
|
|
C ..
|
|
|
C .. Common blocks ..
|
|
|
C ..
|
|
|
|
|
|
C .. Array Arguments ..
|
|
|
INTEGER STATUS(4)
|
|
|
C ..
|
|
|
C .. Intrinsic Functions ..
|
|
|
INTRINSIC INDEX
|
|
|
C ..
|
|
|
IF (SW.EQ.0) THEN
|
|
|
DO 20 J = 1,4
|
|
|
IF (J.EQ.1) CPAR = 'Ne'
|
|
|
IF (J.EQ.2) CPAR = 'Te'
|
|
|
IF (J.EQ.3) CPAR = 'Ti'
|
|
|
IF (J.EQ.4) CPAR = 'Vo'
|
|
|
IF (MODEL.EQ.2) THEN
|
|
|
COEFILE = COEFILEPATH(1:INDEX(COEFILEPATH,' ')-1)//
|
|
|
* 'Coef3D'//CPAR//'.out'
|
|
|
ELSE
|
|
|
COEFILE = COEFILEPATH(1:INDEX(COEFILEPATH,' ')-1)//
|
|
|
* 'Coef4D'//CPAR//'.out'
|
|
|
END IF
|
|
|
|
|
|
OPEN (1,FILE=COEFILE(1:INDEX(COEFILE,' ')-1),STATUS='OLD',
|
|
|
* IOSTAT=STATUS(J))
|
|
|
C READ (1,FMT=9000) IPAR,IAREA,N,M,ICC,KKS,NX,NY,NS,NG,XA,XB,YA,YB,
|
|
|
C * SA,SB,GA,GB
|
|
|
IF (STATUS(J).NE.0) GO TO 20
|
|
|
IF (MODEL.EQ.1) THEN
|
|
|
C read in regional model
|
|
|
READ (1,FMT=*) IPAR,IAREA,N4(J),M4(J),ICC4(J),KKS4(J),
|
|
|
* NX4(J),NY4(J),NS4(J),NG4(J),XA4(J),XB4(J),YA4(J),YB4(J),
|
|
|
* SA4(J),SB4(J),GA4(J),GB4(J)
|
|
|
READ (1,FMT='(100F7.1)') (TT4(J,I),I=1,NY4(J)+KKS4(J))
|
|
|
DO I = 1,N4(J)
|
|
|
READ (1,FMT=9010) (COEFS4(J,IC,I),IC=1,ICC4(J))
|
|
|
END DO
|
|
|
NBKL4(J) = 0
|
|
|
WL4(J,1) = 0D0
|
|
|
READ (1,END=10,FMT='(I4)') NBKL4(J)
|
|
|
IF (NBKL4(J).LE.0) GO TO 10
|
|
|
READ (1,FMT='(100F7.1)') (WL4(J,I),I=1,
|
|
|
* NBKL4(J)+KKS4(J)+KKS4(J)-2)
|
|
|
|
|
|
ELSE
|
|
|
C read in local model
|
|
|
READ (1,FMT=*) IPAR,IAREA,N3(J),M3(J),ICC3(J),KKS3(J),
|
|
|
* NX3(J),NY3(J),NS3(J),NG3(J),XA3(J),XB3(J),YA3(J),YB3(J),
|
|
|
* SA3(J),SB3(J),GA3(J),GB3(J)
|
|
|
READ (1,FMT='(100F7.1)') (TT3(J,I),I=1,NY3(J)+KKS3(J))
|
|
|
DO I = 1,N3(J)
|
|
|
READ (1,FMT=9010) (COEFS3(J,IC,I),IC=1,ICC3(J))
|
|
|
END DO
|
|
|
NBKL3(J) = 0
|
|
|
WL3(J,1) = 0D0
|
|
|
READ (1,END=10,FMT='(I4)') NBKL3(J)
|
|
|
IF (NBKL3(J).LE.0) GO TO 10
|
|
|
READ (1,FMT='(100F7.1)') (WL3(J,I),I=1,
|
|
|
* NBKL3(J)+KKS3(J)+KKS3(J)-2)
|
|
|
|
|
|
END IF
|
|
|
|
|
|
10 CLOSE (1)
|
|
|
9000 FORMAT (2i3,I6,I7,6I3,8F8.2)
|
|
|
9010 FORMAT (4D15.7)
|
|
|
20 CONTINUE
|
|
|
END IF
|
|
|
|
|
|
IF (SW.NE.0) THEN
|
|
|
IF (MODEL.EQ.1) THEN
|
|
|
KKS = KKS4(SW)
|
|
|
NX = NX4(SW)
|
|
|
NY = NY4(SW)
|
|
|
NS = NS4(SW)
|
|
|
NG = NG4(SW)
|
|
|
XA = XA4(SW)
|
|
|
XB = XB4(SW)
|
|
|
YA = YA4(SW)
|
|
|
YB = YB4(SW)
|
|
|
SA = SA4(SW)
|
|
|
SB = SB4(SW)
|
|
|
GA = GA4(SW)
|
|
|
GB = GB4(SW)
|
|
|
N = N4(SW)
|
|
|
M = M4(SW)
|
|
|
ICC = ICC4(SW)
|
|
|
NBKL = NBKL4(SW)
|
|
|
DO I = 1,100
|
|
|
TT(I) = TT4(SW,I)
|
|
|
END DO
|
|
|
DO I = 1,50
|
|
|
WL(I) = WL4(SW,I)
|
|
|
END DO
|
|
|
DO I = 1,10
|
|
|
DO J = 1,5000
|
|
|
COEFS(I,J) = COEFS4(SW,I,J)
|
|
|
END DO
|
|
|
END DO
|
|
|
ELSE
|
|
|
KKS = KKS3(SW)
|
|
|
NX = NX3(SW)
|
|
|
NY = NY3(SW)
|
|
|
NS = NS3(SW)
|
|
|
NG = NG3(SW)
|
|
|
XA = XA3(SW)
|
|
|
XB = XB3(SW)
|
|
|
YA = YA3(SW)
|
|
|
YB = YB3(SW)
|
|
|
SA = SA3(SW)
|
|
|
SB = SB3(SW)
|
|
|
GA = GA3(SW)
|
|
|
GB = GB3(SW)
|
|
|
N = N3(SW)
|
|
|
M = M3(SW)
|
|
|
ICC = ICC3(SW)
|
|
|
NBKL = NBKL3(SW)
|
|
|
DO I = 1,100
|
|
|
TT(I) = TT3(SW,I)
|
|
|
END DO
|
|
|
DO I = 1,50
|
|
|
WL(I) = WL3(SW,I)
|
|
|
END DO
|
|
|
DO I = 1,10
|
|
|
DO J = 1,5000
|
|
|
COEFS(I,J) = COEFS3(SW,I,J)
|
|
|
END DO
|
|
|
END DO
|
|
|
|
|
|
END IF
|
|
|
END IF
|
|
|
|
|
|
RETURN
|
|
|
END
|
|
|
C
|
|
|
C -------------------------------------------------------------
|
|
|
C
|
|
|
SUBROUTINE SETKNOTS3D
|
|
|
INCLUDE 'basis3.h'
|
|
|
C .. Local Scalars ..
|
|
|
INTEGER I,KS
|
|
|
C ..
|
|
|
C .. External Subroutines ..
|
|
|
C ..
|
|
|
C .. Common blocks ..
|
|
|
COMMON /BASINFO1/KKS,NX,NY,NS,NG
|
|
|
COMMON /BASINFO2/XA,XB,YA,YB,SA,SB,GA,GB
|
|
|
COMMON /NODES/TT,WL,NBKL
|
|
|
DOUBLE PRECISION GA,GB,SA,SB,XA,XB,YA,YB
|
|
|
INTEGER KKS,NBKL,NG,NS,NX,NY
|
|
|
DOUBLE PRECISION TT(100),WL(50)
|
|
|
C ..
|
|
|
KS = KKS
|
|
|
DO I = 1,NY + KKS
|
|
|
T2(I) = TT(I)
|
|
|
END DO
|
|
|
XA = 0.d0
|
|
|
XB = 2.0D0*3.1415926535897D0
|
|
|
SA = 0.d0
|
|
|
SB = 2.0D0*3.1415926535897D0
|
|
|
N1 = NX
|
|
|
N2 = NY
|
|
|
N3 = NS
|
|
|
RETURN
|
|
|
END
|
|
|
C
|
|
|
C -------------------------------------------------------------
|
|
|
C
|
|
|
SUBROUTINE SETKNOTS4D
|
|
|
INCLUDE 'basis4.h'
|
|
|
C .. Local Scalars ..
|
|
|
INTEGER I,KS
|
|
|
C ..
|
|
|
C .. External Subroutines ..
|
|
|
C ..
|
|
|
C .. Common blocks ..
|
|
|
COMMON /BASINFO1/KKS,NX,NY,NS,NG
|
|
|
COMMON /BASINFO2/XA,XB,YA,YB,SA,SB,GA,GB
|
|
|
COMMON /NODES/TT,WL,NBKL
|
|
|
DOUBLE PRECISION GA,GB,SA,SB,XA,XB,YA,YB
|
|
|
INTEGER KKS,NBKL,NG,NS,NX,NY
|
|
|
DOUBLE PRECISION TT(100),WL(50)
|
|
|
C ..
|
|
|
KS = KKS
|
|
|
DO I = 1,NY + KKS
|
|
|
T2(I) = TT(I)
|
|
|
END DO
|
|
|
C ... IF LATITUDINAL VARIATION IS A B-SPLINE ..
|
|
|
C .. COMMENT OUT BELOW IF POLYNOMIA USED
|
|
|
IF (NBKL.LE.0) THEN
|
|
|
T4(1) = 12345D0
|
|
|
ELSE
|
|
|
DO I = 1,NBKL + KKS + KKS - 2
|
|
|
T4(I) = WL(I)
|
|
|
END DO
|
|
|
END IF
|
|
|
C ..
|
|
|
XA = 0.d0
|
|
|
XB = 2.0D0*3.1415926535897D0
|
|
|
SA = 0.d0
|
|
|
SB = 2.0D0*3.1415926535897D0
|
|
|
N1 = NX
|
|
|
N2 = NY
|
|
|
N3 = NS
|
|
|
N4 = NG
|
|
|
RETURN
|
|
|
END
|
|
|
C
|
|
|
C -------------------------------------------------------------
|
|
|
C
|
|
|
SUBROUTINE PARAMETER3D(DOY,SLT,LAT,ALT,F107,AP,V3D)
|
|
|
INCLUDE 'basis3.h'
|
|
|
|
|
|
C .. Scalar Arguments ..
|
|
|
DOUBLE PRECISION ALT,AP,DOY,F107,LAT,SLT,V3D
|
|
|
C ..
|
|
|
C .. Local Scalars ..
|
|
|
DOUBLE PRECISION APM,F107M,TWOPI
|
|
|
INTEGER I,IC
|
|
|
C ..
|
|
|
C .. Local Arrays ..
|
|
|
DOUBLE PRECISION BAS(30000),WK(10),WK2(10)
|
|
|
C ..
|
|
|
C .. External Subroutines ..
|
|
|
EXTERNAL SPLBAS3
|
|
|
C ..
|
|
|
C .. Common blocks ..
|
|
|
COMMON /COEF/COEFS,N,M,ICC
|
|
|
INTEGER ICC,M,N
|
|
|
DOUBLE PRECISION COEFS(10,5000)
|
|
|
C ..
|
|
|
TWOPI = 2.d0*3.1415926535897d0
|
|
|
IX1(1) = 1
|
|
|
IX2(1) = 1
|
|
|
IX3(1) = 1
|
|
|
X1(1) = SLT
|
|
|
X2(1) = ALT
|
|
|
X3(1) = DOY
|
|
|
C for peridiocal spline
|
|
|
X1(1) = X1(1)*TWOPI/24.d0
|
|
|
X3(1) = X3(1)*TWOPI/365.d0
|
|
|
C determine basis value, bas
|
|
|
CALL SPLBAS3(N,1,X1(1),BAS)
|
|
|
DO IC = 1,ICC
|
|
|
WK(IC) = 0.d0
|
|
|
DO I = 1,N
|
|
|
WK(IC) = WK(IC) + BAS(I)*COEFS(IC,I)
|
|
|
END DO
|
|
|
END DO
|
|
|
F107M = 135.D0
|
|
|
APM = 15.D0
|
|
|
DO IC = 1,ICC
|
|
|
WK2(IC) = 0.d0
|
|
|
END DO
|
|
|
WK2(1) = 1.d0
|
|
|
WK2(2) = (F107-F107M)/100.d0
|
|
|
WK2(3) = (AP-APM)/10.d0
|
|
|
WK2(4) = WK2(2)*WK2(3)
|
|
|
V3D = 0.d0
|
|
|
DO IC = 1,ICC
|
|
|
V3D = V3D + WK(IC)*WK2(IC)
|
|
|
END DO
|
|
|
RETURN
|
|
|
END
|
|
|
C
|
|
|
C -------------------------------------------------------------
|
|
|
C
|
|
|
SUBROUTINE PARAMETER4D(DOY,SLT,LAT,ALT,F107,AP,V4D)
|
|
|
INCLUDE 'basis4.h'
|
|
|
C .. Scalar Arguments ..
|
|
|
DOUBLE PRECISION ALT,AP,DOY,F107,LAT,SLT,V4D
|
|
|
C ..
|
|
|
C .. Local Scalars ..
|
|
|
DOUBLE PRECISION APM,F107M,TWOPI
|
|
|
INTEGER I,IC
|
|
|
C ..
|
|
|
C .. Local Arrays ..
|
|
|
DOUBLE PRECISION BAS(30000),WK(10),WK2(10)
|
|
|
C ..
|
|
|
C .. External Subroutines ..
|
|
|
EXTERNAL SPLBASP4
|
|
|
C ..
|
|
|
C .. Common blocks ..
|
|
|
COMMON /COEF/COEFS,N,M,ICC
|
|
|
INTEGER ICC,M,N
|
|
|
DOUBLE PRECISION COEFS(10,5000)
|
|
|
C ..
|
|
|
TWOPI = 2.d0*3.1415926535897d0
|
|
|
IX1(1) = 1
|
|
|
IX2(1) = 1
|
|
|
IX3(1) = 1
|
|
|
IX4(1) = 1
|
|
|
X1(1) = SLT
|
|
|
X2(1) = ALT
|
|
|
X3(1) = DOY
|
|
|
X4(1) = LAT
|
|
|
C for peridiocal spline
|
|
|
X1(1) = X1(1)*TWOPI/24.d0
|
|
|
X3(1) = X3(1)*TWOPI/365.d0
|
|
|
C determine basis value, bas
|
|
|
CALL SPLBASP4(N,1,X1(1),BAS)
|
|
|
DO IC = 1,ICC
|
|
|
WK(IC) = 0.d0
|
|
|
DO I = 1,N
|
|
|
WK(IC) = WK(IC) + BAS(I)*COEFS(IC,I)
|
|
|
END DO
|
|
|
END DO
|
|
|
C print*,t4(1),wk(1),coefs(1,1),coefs(2,1),coefs(3,1)
|
|
|
F107M = 135.D0
|
|
|
APM = 15.D0
|
|
|
DO IC = 1,ICC
|
|
|
WK2(IC) = 0.d0
|
|
|
END DO
|
|
|
WK2(1) = 1.d0
|
|
|
WK2(2) = (F107-F107M)/100.d0
|
|
|
WK2(3) = (AP-APM)/15.d0
|
|
|
WK2(4) = WK2(2)*WK2(3)
|
|
|
V4D = 0.d0
|
|
|
DO IC = 1,ICC
|
|
|
V4D = V4D + WK(IC)*WK2(IC)
|
|
|
END DO
|
|
|
C WRITE(18,'(F8.1,1000E15.7)')ALT,(BAS(I),I=1,N)
|
|
|
RETURN
|
|
|
END
|
|
|
|
