2、linux平台代码
! ap8ae8.f90
!
! FUNCTIONS/SUBROUTINES exported from ap8ae8.dll:
! ap8ae8 - subroutine
!
subroutine ap8ae8(OutputFileName,lenOfOutputFileName,
1coordType,
1L_FROMFILE_Start,L_FROMFILE_End, L_FROMFILE_Gap,
1BB0_FROMFILE_Start,BB0_FROMFILE_End,BB0_FROMFILE_Gap,xlatStart,
1xlatEnd,xlatGap,xlongStart,xlongEnd,xlongGap,yearStart,yearEnd,
1yearGap,heightStart,heightEnd,heightGap, MODEL, FLUXTYPE, E1, E2)
! Expose function apae20151108flux to users of this DLL
!
!DEC$ ATTRIBUTES DLLEXPORT::ap8ae8
! Variables
! Body of apae20151108flux
real L_FROMFILE_Start,L_FROMFILE_End, L_FROMFILE_Gap,
1BB0_FROMFILE_Start,BB0_FROMFILE_End,BB0_FROMFILE_Gap,
1xlatStart,xlatEnd,xlatGap,xlongStart,xlongEnd,xlongGap,
1yearStart,yearEnd,yearGap,heightStart,
1heightEnd,heightGap
real(4) E1, E2
integer coordType,MODEL, FLUXTYPE
CHARACTER, dimension(*), intent(in) :: OutputFileName
! CHARACTER(*) :: OutputFileName
! OutputFileName(1)='D'
! OutputFileName(2)=':'
! OutputFileName(3)='\'
! OutputFileName(4)='\'
! OutputFileName(5)='T'
! OutputFileName(6)='.'
! OutputFileName(7)='t'
! OutputFileName(8)='x'
! OutputFileName(9)='t'
CHARACTER(len=lenOfOutputFileName):: FNM
INTEGER :: lenFNM
do i=1,lenOfOutputFileName
FNM(i:i)=OutputFileName(i)
end do
! OPEN(UNIT=21,FILE="argument.txt", status="old" ,access='append',
! 1FORM='FORMATTED')
OPEN(UNIT=21,FILE="argument.txt",status="replace",access='append',
1FORM='FORMATTED')
WRITE(21,91113)
! WRITE(*,91113)
1coordType,
1L_FROMFILE_Start,L_FROMFILE_End, L_FROMFILE_Gap,
1BB0_FROMFILE_Start,BB0_FROMFILE_End,BB0_FROMFILE_Gap,xlatStart,
1xlatEnd,xlatGap,xlongStart,xlongEnd,xlongGap,yearStart,yearEnd,
1yearGap,heightStart,heightEnd,heightGap, MODEL, FLUXTYPE, E1, E2
91113 FORMAT(I8,',',F15.7,',',F15.7,',',F15.7,',',F15.7,
1',',
1F15.7,',',F15.7,',',F8.3,',',F8.3,',',F8.3,',',F8.3,',',F8.3,',',
1F8.3,',',F8.3,',',F8.3,',',F8.3,',',F10.3,',',F10.3,',',
1F10.3,',',
1I8,',',I8,',',F15.7,',',F15.7)
close(21,status='delete')
call ap8ae8main(OutputFileName,lenOfOutputFileName,coordType,
1L_FROMFILE_Start,L_FROMFILE_End, L_FROMFILE_Gap,
1BB0_FROMFILE_Start,BB0_FROMFILE_End,BB0_FROMFILE_Gap,xlatStart,
1xlatEnd,xlatGap,xlongStart,xlongEnd,xlongGap,yearStart,yearEnd,
1yearGap,heightStart,heightEnd,heightGap, MODEL, FLUXTYPE, E1, E2)
write(*,*)"hello!"
end subroutine
subroutine ap8ae8main(OutputFileName,lenOfOutputFileName,
1coordType,L_FROMFILE_Start,L_FROMFILE_End, L_FROMFILE_Gap,
1BB0_FROMFILE_Start,BB0_FROMFILE_End,BB0_FROMFILE_Gap,xlatStart,
1xlatEnd,xlatGap,xlongStart,xlongEnd,xlongGap,yearStart,yearEnd,
1yearGap,heightStart,heightEnd,heightGap, MODEL, FLUXTYPE, E1, E2)
real L_FROMFILE_Start,L_FROMFILE_End, L_FROMFILE_Gap,
1BB0_FROMFILE_Start,BB0_FROMFILE_End,BB0_FROMFILE_Gap,xlatStart,
1xlatEnd,xlatGap,xlongStart,xlongEnd,xlongGap,yearStart,yearEnd,
1yearGap,heightStart,heightEnd,heightGap
real(4) L_FROMFILE, BB0_FROMFILE, E1, E2,FLUX1,YEAR,LATI,LONGI,
1HEIGHT
integer MODEL, FLUXTYPE,coordType
LOGICAL(4) ALIVE
integer pp
common pp
CHARACTER, dimension(*), intent(in) :: OutputFileName
CHARACTER(len=lenOfOutputFileName):: FNM
INTEGER :: lenFNM
do i=1,lenOfOutputFileName
FNM(i:i)=OutputFileName(i)
end do
IF(coordType.EQ.1) THEN
OPEN(UNIT=10,FILE=FNM, status="replace",access='append',
1FORM='FORMATTED')
WRITE(10,90000)
90000 FORMAT(' L B/B0 MODEL FLUXTYPE
1E1 E2 FLUX1')
close(10)
do L_FROMFILE=L_FROMFILE_Start,L_FROMFILE_End, L_FROMFILE_Gap
do BB0_FROMFILE=BB0_FROMFILE_Start,BB0_FROMFILE_End,
1BB0_FROMFILE_Gap
call ap8ae8Computing(L_FROMFILE, BB0_FROMFILE, MODEL,
1 FLUXTYPE, E1, E2,FLUX1)
OPEN(UNIT=10,FILE=FNM, status="old" ,access='append',
1FORM='FORMATTED')
WRITE(10,90001) L_FROMFILE, BB0_FROMFILE, MODEL, FLUXTYPE,
1E1, E2,FLUX1
90001 FORMAT(F15.7,',',F15.7,',',I8,',',I8,',',F15.7,',',
1F15.7,',',E31.8)
close(10)
end do
end do
end if
IF(coordType.EQ.2) THEN
INQUIRE(FILE=FNM,EXIST=ALIVE)
! IF(ALIVE.EQ.0.0)THEN
IF(ALIVE)THEN
! goto 11234
goto 12134
ELSE
! goto 12134
goto 11234
ENDIF
11234 OPEN(UNIT=10,FILE=FNM, status="replace" ,access='append',
1FORM='FORMATTED')
WRITE(10,90002)
90002 FORMAT(' YEAR, LONG, LAT, HEIGHT, MODEL, FLUXTYPE,
1E1, E2, FLUX1')
close(10)
C OPEN(UNIT=20,FILE="no_cal.txt", status="old" ,access='append')
C WRITE(20,*)"hello"
C close(20)
12134 do LATI=xlatStart,xlatEnd,xlatGap
do LONGI=xlongStart,xlongEnd,xlongGap
do year=yearStart,yearEnd,yearGap
do height=heightStart,heightEnd,heightGap
pp=0
call getLBBO(YEAR,LATI,LONGI,HEIGHT,L_FROMFILE,BB0_FROMFILE)
call ap8ae8Computing(L_FROMFILE, BB0_FROMFILE, MODEL, FLUXTYPE,
1E1, E2,FLUX1)
OPEN(UNIT=10,FILE=FNM, status="old" ,access='append',
1FORM='FORMATTED')
WRITE(10,90003)YEAR,LONGI,LATI,HEIGHT,MODEL, FLUXTYPE, E1,
1E2,FLUX1
90003 FORMAT(F8.3,',',F8.3,',',F8.3,',',F10.3,',',I8,',',
1I8,',',F15.7,',',F15.7,',',E31.8)
pp=1
close(10)
end do
end do
end do
end do
ENDIF
IF(pp.EQ.1)THEN
! OPEN(UNIT=20,FILE="no_cal.txt", status="old" ,access='append',
! 1FORM='FORMATTED')
OPEN(UNIT=20,FILE="no_cal.txt", status="replace" ,access='append',
1FORM='FORMATTED')
WRITE(20,90113)
! WRITE(*,90113)
1YEAR,LONGI,LATI,HEIGHT,MODEL, FLUXTYPE, E1,
1E2
90113 FORMAT(F8.3,',',F8.3,',',F8.3,',',F10.3,',',I8,',',
1I8,',',F15.7,',',F15.7)
close(20,status='delete')
endif
return
end
subroutine ap8ae8Computing(L_FROMFILE, BB0_FROMFILE, MODEL,
1FLUXTYPE, E1, E2,FLUX1)
integer EI
INTEGER INFILE,OUTFILE
INTEGER MODEL, FLUXTYPE
REAL(4) L_FROMFILE, BB0_FROMFILE, E1, E2
9 DIMENSION E(7), IHEAD(8),MAP(20000),EDA(3)
REAL XL(2,10),FLUX(7)
DOUBLE PRECISION AF(7,10,10),DF(5,10,10)
DIMENSION OUTPUT_ARR(6)
DOUBLE PRECISION OUTPUT_ARR7
CHARACTER*1 ITE(5,10,10)
CHARACTER*4 BLTEX,LBTEX
CHARACTER*6 NAME,MNAME(8)
CHARACTER*10 FNAME
! CHARACTER*12 OUTFILENAME
! OUTFILENAME="ap8ae8output"
!soil CHARACTER*11 INFILENAME
!soil输出文件名赋值与处理
!CHARACTER, dimension(*), intent(in) :: JOUTFILENAME
!CHARACTER(len=lenOUTFILENAME):: OUTFILENAME
!INTEGER :: lenOUTFILENAME
DATA MNAME /'AP8MAX','AP8MIN','AE4MAX','AE4MIN',
1'AEI7HI', 'AEI7LO','AE8MAX','AE8MIN'/
DATA E,XL,AF,DF /1227*0.0/,EDA/.05,.1,.2/
!
DO 1567 I=1,5
DO 1567 K=1,10
DO 1567 L=1,10
1567 ITE(I,K,L)=' '
!
! I/O UNIT NUMBERS
!
EGNR=5 ! INPUT
MONITO=6 ! MONITOR
IUAEAP=15 ! MODEL COEFFICIENTS INPUT
OUTFILE=16
INFILE=17
!
! INPUT OF PARAMETERS
!
!soil WRITE(MONITO, *) 'TYPE INPUT FILE NAME: '
!soil READ '(13A)', INFILENAME
!soil OPEN(UNIT=INFILE,FILE=INFILENAME,STATUS='OLD', FORM='FORMATTED',ERR=8000)
!soil WRITE(MONITO, *) 'TYPE RESULT FILE NAME: '
!soil READ '(13A)', OUTFILENAME
!soil OPEN(UNIT=OUTFILE,FILE=OUTFILENAME, STATUS='NEW',FORM='FORMATTED')
! OPEN(UNIT=OUTFILE,FILE=OUTFILENAME, STATUS = 'REPLACE',FORM='FORMATTED')
! WRITE(OUTFILE,3910)
!3910 FORMAT('L B/B0 MODEL FLUXTYPE E1 E2 FLUX1')
!soil 1000 READ(INFILE, *, END=9999) L_FROMFILE, BB0_FROMFILE, MODEL, FLUXTYPE, E1, E2
IF (FLUXTYPE.EQ.1) THEN
NE = 1
E(1)=E1
endif
if (FLUXTYPE.EQ.2) THEN
NE = 2
E(1)=E1
E(2)=E2
ENDIF
!---------------4. window: model type-------------------------------
MTYPE = MODEL
NAME=MNAME(MTYPE)
WRITE(FNAME,1129) NAME
! ASCII
! Using the ASCII coefficient files instead of the binary
1129 FORMAT(A6,'.ASC')
OPEN(IUAEAP,FILE=FNAME,STATUS='OLD',FORM='FORMATTED')
READ(IUAEAP,1301) IHEAD
NMAP=IHEAD(8)
READ(IUAEAP,1301) (MAP(I),I=1,NMAP)
1301 FORMAT(1X,12I6)
CLOSE(IUAEAP)
!-------------------7. window: number of L-values--------------------
7012 NL=1
!----------------8. window (b): L-values-------------------------
2798 XL(1,1) = L_FROMFILE
!-----------------9. window: number of B/B0 values------------------
7013 NB=1
!------------10. window (a): B/B0 grid--------------------------------
3812 XL(2,1) = BB0_FROMFILE
!----------------THE L-VALUE LOOP-------------------------------
L=1
FL=XL(1,L)
!----------------THE B LOOP-------------------------------------
I=1
BB0 = XL(2,I)
CALL TRARA1(IHEAD,MAP,FL,BB0,E,FLUX,NE)
!----------------THE ENERGY LOOP--------------------------------
DO 1110 K=1,NE
AF(K,L,I)=0.0
IF(FLUX(K).GT.0.0) AF(K,L,I)=10.**FLUX(K)
IF(K.EQ.1) GOTO 1110
DF(K-1,L,I)=ABS(AF(K,L,I)-AF(K-1,L,I))/ (E(K)-E(K-1))
IF(AF(K,L,I).LE.0.0) DF(K-1,L,I)=0.0
1110 CONTINUE
!-----------------TESTS VALIDITY OF DIFFERENTIAL FLUX------------
DO 7788 K=1,NE-1
EI=E(K+1)
ITEST=0
EDIFF=EI-E(K)
!-----------------IS ENERGY INTERVALL LARGE ENOUGH ?-------------
IEI=1
IF(EI.GT.0.10) IEI=2
IF(EI.GT.0.25) IEI=3
IF(EDIFF.LT.EDA(IEI)) ITEST=1
DO 7788 L=1,NL
ITT=0
IF(XL(1,L).LT.1.2) ITT=1
DO 7788 I=1,NB
ITE(K,L,I)=' '
IF(ITEST+ITT.NE.0) ITE(K,L,I)='?'
7788 CONTINUE
!------------------TABLE OUTPUT------------------------------------
!soil WRITE(MONITO,9101)
!soil 9101 FORMAT(1X/)
! integral flux: FLUXTYPE=1
if(FLUXTYPE.EQ.1) then
JTAB = 1
end if
! differential flux: FLUXTYPE=2
if(FLUXTYPE.EQ.2) then
JTAB = 3
end if
IBL=2
BLTEX='B/B0'
LBTEX=' L '
N=NL
NN=NB
IBLTAB=1
IF(IBL.EQ.1) IBLTAB=2
INDE = 1
OUTPUT_ARR(1) = XL(1,1) !L
OUTPUT_ARR(2) = XL(2,1) !B/B0
OUTPUT_ARR(3) = MODEL !MODEL NUMBER
OUTPUT_ARR(4) = FLUXTYPE !FLUXTUPE
OUTPUT_ARR(5) = E1 !E1
OUTPUT_ARR(6) = E2 !E2
BLV=XL(IBL,INDE)
DO 1094 I=1,N
IL=INDE
IB=INDE
IF(IBL.EQ.1) THEN
IB=I
XX=XL(2,I)
ELSE
IL=I
XX=XL(1,I)
ENDIF
IF(JTAB.LT.3) THEN
OUTPUT_ARR7 = AF(1,IL,IB)
ELSE
OUTPUT_ARR7 = DF(1,IL,IB)
ENDIF
FLUX1=OUTPUT_ARR7
!WRITE(OUTFILE, *) OUTPUT_ARR
!WRITE(OUTFILE, *) OUTPUT_ARR7
! WRITE(OUTFILE,'(F15.7)') OUTPUT_ARR(1)
! WRITE(OUTFILE,'(F15.7)') OUTPUT_ARR(2)
! WRITE(OUTFILE,'(F15.7)') OUTPUT_ARR(3)
! WRITE(OUTFILE,'(F15.7)') OUTPUT_ARR(4)
! WRITE(OUTFILE,'(F15.7)') OUTPUT_ARR(5)
! WRITE(OUTFILE,'(F15.7)') OUTPUT_ARR(6)
! WRITE(OUTFILE,'(E31.8)') OUTPUT_ARR7
!soil WRITE(MONITO, *) OUTPUT_ARR
! close(OUTFILE)
1094 CONTINUE
!soil goto 1000
!soil 8000 write(MONITO, *) 'input file not found'
!soil 9999 STOP
return
!!
end
!!!!!trmfun.for
SUBROUTINE TRARA1(DESCR,MAP,FL,BB0,E,F,N)
!***********************************************************************
!*** TRARA1 FINDS PARTICLE FLUXES FOR GIVEN ENERGIES, MAGNETIC FIELD ***
!*** STRENGTH AND L-VALUE. FUNCTION TRARA2 IS USED TO INTERPOLATE IN ***
!*** B-L-SPACE. ***
!*** INPUT: DESCR(8) HEADER OF SPECIFIED TRAPPED RADITION MODEL ***
!*** MAP(...) MAP OF TRAPPED RADITION MODEL ***
!*** (DESCR AND MAP ARE EXPLAINED AT THE BEGIN ***
!*** OF THE MAIN PROGRAM MODEL) ***
!*** N NUMBER OF ENERGIES ***
!*** E(N) ARRAY OF ENERGIES IN MEV ***
!*** FL L-VALUE ***
!*** BB0 =B/B0 MAGNETIC FIELD STRENGTH NORMALIZED ***
!*** TO FIELD STRENGTH AT MAGNETIC EQUATOR ***
!*** OUTPUT: F(N) DECADIC LOGARITHM OF INTEGRAL FLUXES IN ***
!*** PARTICLES/(CM*CM*SEC) ***
!***********************************************************************
LOGICAL S0,S1,S2
DIMENSION E(N),F(N),MAP(1)
INTEGER DESCR(8)
COMMON/TRA2/FISTEP
DATA F1,F2/1.001,1.002/
!
FISTEP=DESCR(7)/DESCR(2)
ESCALE=DESCR(4)
FSCALE=DESCR(7)
XNL=AMIN1(15.6,ABS(FL))
NL=XNL*DESCR(5)
IF(BB0.LT.1.) BB0=1.
NB=(BB0-1.)*DESCR(6)
!
! I2 IS THE NUMBER OF ELEMENTS IN THE FLUX MAP FOR THE FIRST ENERGY.
! I3 IS THE INDEX OF THE LAST ELEMENT OF THE SECOND ENERGY MAP.
! L3 IS THE LENGTH OF THE MAP FOR THE THIRD ENERGY.
! E1 IS THE ENERGY OF THE FIRST ENERGY MAP (UNSCALED)
! E2 IS THE ENERGY OF THE SECOND ENERGY MAP (UNSCALED)
!
I1=0
I2=MAP(1)
I3=I2+MAP(I2+1)
L3=MAP(I3+1)
E1=MAP(I1+2)/ESCALE
E2=MAP(I2+2)/ESCALE
!
! S0, S1, S2 ARE LOGICAL VARIABLES WHICH INDICATE WHETHER THE FLUX FOR
! A PARTICULAR E, B, L POINT HAS ALREADY BEEN FOUND IN A PREVIOUS CALL
! TO FUNCTION TRARA2. IF NOT, S.. =.TRUE.
!
S1=.TRUE.
S2=.TRUE.
!
! ENERGY LOOP
!
DO 3 IE=1,N
!
! FOR EACH ENERGY E(I) FIND THE SUCCESSIVE ENERGIES E0,E1,E2 IN
! MODEL MAP, WHICH OBEY E0 < E1 < E(I) < E2 .
!
1 IF((E(IE).LE.E2).OR.(L3.EQ.0)) GOTO 2
I0=I1
I1=I2
I2=I3
I3=I3+L3
L3=MAP(I3+1)
E0=E1
E1=E2
E2=MAP(I2+2)/ESCALE
S0=S1
S1=S2
S2=.TRUE.
F0=F1
F1=F2
GOTO 1
!
! CALL TRARA2 TO INTERPOLATE THE FLUX-MAPS FOR E1,E2 IN L-B/B0-
! SPACE TO FIND FLUXES F1,F2 [IF THEY HAVE NOT ALREADY BEEN
! CALCULATED FOR A PREVIOUS E(I)].
!
2 IF(S1) F1=TRARA2(MAP(I1+3),NL,NB)/FSCALE
IF(S2) F2=TRARA2(MAP(I2+3),NL,NB)/FSCALE
S1=.FALSE.
S2=.FALSE.
!
! FINALLY, INTERPOLATE IN ENERGY.
!
F(IE)=F1+(F2-F1)*(E(IE)-E1)/(E2-E1)
IF(F2.GT.0.0) GOTO 3
IF(I1.EQ.0) GOTO 3
!
! --------- SPECIAL INTERPOLATION ---------------------------------
! IF THE FLUX FOR THE SECOND ENERGY CANNOT BE FOUND (I.E. F2=0.0),
! AND THE ZEROTH ENERGY MAP HAS BEEN DEFINED (I.E. I1 NOT EQUAL 0),
! THEN INTERPOLATE USING THE FLUX MAPS FOR THE ZEROTH AND FIRST
! ENERGY AND CHOOSE THE MINIMUM OF THIS INTERPOLATIONS AND THE
! INTERPOLATION THAT WAS DONE WITH F2=0.
!
IF(S0) F0=TRARA2(MAP(I0+3),NL,NB)/FSCALE
S0=.FALSE.
F(IE)=AMIN1(F(IE),F0+(F1-F0)*(E(IE)-E0)/(E1-E0))
!
! THE LOGARITHMIC FLUX IS ALWAYS KEPT GREATER OR EQUAL ZERO.
!
3 F(IE)=AMAX1(F(IE),0.)
RETURN
END
!
FUNCTION TRARA2(MAP,IL,IB)
!*****************************************************************
!*** TRARA2 INTERPOLATES LINEARLY IN L-B/B0-MAP TO OBTAIN ***
!*** THE LOGARITHM OF INTEGRAL FLUX AT GIVEN L AND B/B0. ***
!*** INPUT: MAP(..) IS SUB-MAP (FOR SPECIFIC ENERGY) OF ***
!*** TRAPPED RADIATION MODEL MAP ***
!*** IL SCALED L-VALUE ***
!*** IB SCALED B/B0-1 ***
!*** OUTPUT: TRARA2 SCALED LOGARITHM OF PARTICLE FLUX ***
!*****************************************************************
!*** SEE MAIN PROGRAM 'MODEL' FOR EXPLANATION OF MAP FORMAT ***
!*** SCALING FACTORS. ***
!*** THE STEPSIZE FOR THE PARAMETERIZATION OF THE LOGARITHM ***
!*** OF FLUX IS OBTAINED FROM 'COMMON/TRA2/'. ***
!*****************************************************************
DIMENSION MAP(1)
COMMON/TRA2/FISTEP
FNL=IL
FNB=IB
ITIME=0
I2=0
!
! FIND CONSECUTIVE SUB-SUB-MAPS FOR SCALED L-VALUES LS1,LS2,
! WITH IL LESS OR EQUAL LS2. L1,L2 ARE LENGTHS OF SUB-SUB-MAPS.
! I1,I2 ARE INDECES OF FIRST ELEMENTS MINUS 1.
!
1 L2=MAP(I2+1)
IF(MAP(I2+2).GT.IL) GOTO 2
I1=I2
L1=L2
I2=I2+L2
GOTO 1
2 CONTINUE
!
! IF SUB-SUB-MAPS ARE EMPTY, I. E. LENGTH LESS 4, THAN TRARA2=0
!
IF((L1.LT.4).AND.(L2.LT.4)) GOTO 50
!
! IF FLOG2 LESS FLOG1, THAN LS2 FIRST MAP AND LS1 SECOND MAP
!
IF(MAP(I2+3).GT.MAP(I1+3)) GOTO 10
5 KT=I1
I1=I2
I2=KT
KT=L1
L1=L2
L2=KT
!
! DETERMINE INTERPOLATE IN SCALED L-VALUE
!
10 FLL1=MAP(I1+2)
FLL2=MAP(I2+2)
DFL=(FNL-FLL1)/(FLL2-FLL1)
FLOG1=MAP(I1+3)
FLOG2=MAP(I2+3)
FKB1=0.
FKB2=0.
IF(L1.LT.4) GOTO 32
!
! B/B0 LOOP
!
DO 17 J2=4,L2
FINCR2=MAP(I2+J2)
IF(FKB2+FINCR2.GT.FNB) GOTO 23
FKB2=FKB2+FINCR2
17 FLOG2=FLOG2-FISTEP
ITIME=ITIME+1
IF(ITIME.EQ.1)GO TO 5
GO TO 50
23 IF(ITIME.EQ.1)GO TO 30
IF(J2.EQ.4)GO TO 28
SL2=FLOG2/FKB2
DO 27 J1=4,L1
FINCR1=MAP(I1+J1)
FKB1=FKB1+FINCR1
FLOG1=FLOG1-FISTEP
FKBJ1=((FLOG1/FISTEP)*FINCR1+FKB1)/((FINCR1/FISTEP)*SL2+1.)
IF(FKBJ1.LE.FKB1) GOTO 31
27 CONTINUE
IF(FKBJ1.LE.FKB2) GOTO 50
31 IF(FKBJ1.LE.FKB2) GOTO 29
FKB1=0.
30 FKB2=0.
32 J2=4
FINCR2=MAP(I2+J2)
FLOG2=MAP(I2+3)
FLOG1=MAP(I1+3)
28 FLOGM=FLOG1+(FLOG2-FLOG1)*DFL
FKBM=0.
FKB2=FKB2+FINCR2
FLOG2=FLOG2-FISTEP
SL2=FLOG2/FKB2
IF(L1.LT.4) GOTO 35
J1=4
FINCR1=MAP(I1+J1)
FKB1=FKB1+FINCR1
FLOG1=FLOG1-FISTEP
SL1=FLOG1/FKB1
GOTO 15
29 FKBM=FKBJ1+(FKB2-FKBJ1)*DFL
FLOGM=FKBM*SL2
FLOG2=FLOG2-FISTEP
FKB2=FKB2+FINCR2
SL1=FLOG1/FKB1
SL2=FLOG2/FKB2
15 IF(SL1.LT.SL2) GOTO 20
FKBJ2=((FLOG2/FISTEP)*FINCR2+FKB2)/((FINCR2/FISTEP)*SL1+1.)
FKB=FKB1+(FKBJ2-FKB1)*DFL
FLOG=FKB*SL1
IF(FKB.GE.FNB) GOTO 60
FKBM=FKB
FLOGM=FLOG
IF(J1.GE.L1) GOTO 50
J1=J1+1
FINCR1=MAP(I1+J1)
FLOG1=FLOG1-FISTEP
FKB1=FKB1+FINCR1
SL1=FLOG1/FKB1
GOTO 15
20 FKBJ1=((FLOG1/FISTEP)*FINCR1+FKB1)/((FINCR1/FISTEP)*SL2+1.)
FKB=FKBJ1+(FKB2-FKBJ1)*DFL
FLOG=FKB*SL2
IF(FKB.GE.FNB) GOTO 60
FKBM=FKB
FLOGM=FLOG
IF(J2.GE.L2) GOTO 50
J2=J2+1
FINCR2=MAP(I2+J2)
FLOG2=FLOG2-FISTEP
FKB2=FKB2+FINCR2
SL2=FLOG2/FKB2
GOTO 15
35 FINCR1=0.
SL1=-900000.
GOTO 20
60 IF(FKB.LT.FKBM+1.E-10) GOTO 50
TRARA2=FLOGM+(FLOG-FLOGM)*((FNB-FKBM)/(FKB-FKBM))
TRARA2=AMAX1(TRARA2,0.)
RETURN
50 TRARA2=0.
RETURN
END
!****************************************************************************
!
! PROGRAM: getLBBO
!
! PURPOSE: Entry point for the console application.
!
!****************************************************************************
SUBROUTINE getLBBO(YEAR,LATI,LONGI,HEIGHT,L_FROMFILE,BB0_FROMFILE)
! NASAIgrfDll.f90
!
! FUNCTIONS/SUBROUTINES exported from NASAIgrfDll.dll:
! NASAIgrfDll - subroutine
!
! subroutine igrf_sub(xlat,xlong,year,height,LONGI,LATI,DIMO,dip,dec,BABS,BDOWN,BNORTH,BEQU,BEAST,xl,ICODE)
!----------------------------------------------------------------
! INPUT:
! xlat geodatic latitude in degrees
! xlong geodatic longitude in degrees
! year decimal year (year+month/12.0-0.5 or year+day-of-year/365
! or 366 if leap year)
! height height in km
! OUTPUT:
! xl L value
! icode =1 L is correct; =2 L is not correct;
! =3 an approximation is used
! dip geomagnetic inclination in degrees
! dec geomagnetic declination in degress
!----------------------------------------------------------------
INTEGER EGNR,AGNR,OGNR
REAL LATI,LONGI,year,height,bbo,L_FROMFILE,BB0_FROMFILE
COMMON/GENER/ UMR,ERA,AQUAD,BQUAD
CALL INITIZE
ibbb=0
ALOG2=ALOG(2.)
ISTART=1
! lati=xlat
! longi=xlong
!
!----------------CALCULATE PROFILES-----------------------------------
!
CALL FELDCOF(YEAR,DIMO)
CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS)
CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1,BEQU)
DIP=ASIN(BDOWN/BABS)/UMR
DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR
bbo=BABS/BAB1
L_FROMFILE=XL
BB0_FROMFILE=bbo
return
end
!
!
! SHELLIG.FOR
!
! 11/01/91 SHELLG: lowest starting point for B0 search is 2
! 1/27/92 Adopted to IGRF-91 coeffcients model
! 2/05/92 Reduce variable-names: INTER(P)SHC,EXTRA(P)SHC,INITI(ALI)ZE
! 8/08/95 Updated to IGRF-45-95; new coeff. DGRF90, IGRF95, IGRF95S
! 5/31/00 Updated to IGRF-45-00; new coeff.: IGRF00, IGRF00s
! 3/24/05 Updated to IGRF-45-10; new coeff.: IGRF05, IGRF05s
! 4/25/05 ENTRY FELDI(XI,H) and DO 1111 I=1,7 [Alexey Petrov]
! 7/22/09 SHELLG: NMAX=13 for DGRF00 and IGRF05; H/G-arrays(195)
! 2/26/10 FELDCOF: Updated IGRF45-15; new coeff: DGRF05, IGRF10, IGRF10S
! 4/29/10 H/H-arrays(196); FELDCOF: corrected IGRF00 and ..00S
! 4/29/10 Change to new files dgrf%%%%.asc; new GETSHC; char*12 to 13
!*********************************************************************
! SUBROUTINES FINDB0, SHELLG, STOER, FELDG, FELDCOF, GETSHC, *
! INTERSHC, EXTRASHC, INITIZE *
!*********************************************************************
!*********************************************************************
!
!
SUBROUTINE FINDB0(STPS,BDEL,VALUE,BEQU,RR0)
!--------------------------------------------------------------------
! FINDS SMALLEST MAGNETIC FIELD STRENGTH ON FIELD LINE
!
!INPUT: STPS STEP SIZE FOR FIELD LINE TRACING
! COMMON/FIDB0/
! SP DIPOLE ORIENTED COORDINATES FORM SHELLG; P(1,*),
! P(2,*), P(3,*) CLOSEST TO MAGNETIC EQUATOR
! BDEL REQUIRED ACCURACY = [ B(LAST) - BEQU ] / BEQU
! B(LAST) IS FIELD STRENGTH BEFORE BEQU
!
! OUTPUT: VALUE =.FALSE., IF BEQU IS NOT MINIMAL VALUE ON FIELD LINE
! BEQU MAGNETIC FIELD STRENGTH AT MAGNETIC EQUATOR
! RR0 EQUATORIAL RADIUS NORMALIZED TO EARTH RADIUS
! BDEL FINAL ACHIEVED ACCURACY
!--------------------------------------------------------------------
DIMENSION P(8,4),SP(3)
LOGICAL VALUE
COMMON/FIDB0/ SP
!
STEP=STPS
IRUN=0
7777 IRUN=IRUN+1
IF(IRUN.GT.5) THEN
VALUE=.FALSE.
GOTO 8888
ENDIF
!*********************FIRST THREE POINTS
P(1,2)=SP(1)
P(2,2)=SP(2)
P(3,2)=SP(3)
STEP=-SIGN(STEP,P(3,2))
CALL STOER(P(1,2),BQ2,R2)
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)
!******************INVERT SENSE IF REQUIRED
IF(BQ3.LE.BQ1) GOTO 2
STEP=-STEP
R3=R1
BQ3=BQ1
DO 1 I=1,5
ZZ=P(I,1)
P(I,1)=P(I,3)
1 P(I,3)=ZZ
!******************INITIALIZATION
2 STEP12=STEP/12.
VALUE=.TRUE.
BMIN=1.E4
BOLD=1.E4
!******************CORRECTOR (FIELD LINE TRACING)
N=0
5555 P(1,3)=P(1,2)+STEP12*(5.*P(4,3)+8.*P(4,2)-P(4,1))
N=N+1
P(2,3)=P(2,2)+STEP12*(5.*P(5,3)+8.*P(5,2)-P(5,1))
!******************PREDICTOR (FIELD LINE TRACING)
P(1,4)=P(1,3)+STEP12*(23.*P(4,3)-16.*P(4,2)+5.*P(4,1))
P(2,4)=P(2,3)+STEP12*(23.*P(5,3)-16.*P(5,2)+5.*P(5,1))
P(3,4)=P(3,3)+STEP
CALL STOER(P(1,4),BQ3,R3)
DO 1111 J=1,3
! DO 1111 I=1,8
DO 1111 I=1,7
1111 P(I,J)=P(I,J+1)
B=SQRT(BQ3)
IF(B.LT.BMIN) BMIN=B
IF(B.LE.BOLD) THEN
BOLD=B
ROLD=1./R3
SP(1)=P(1,4)
SP(2)=P(2,4)
SP(3)=P(3,4)
GOTO 5555
ENDIF
IF(BOLD.NE.BMIN) THEN
VALUE=.FALSE.
ENDIF
BDELTA=(B-BOLD)/BOLD
IF(BDELTA.GT.BDEL) THEN
STEP=STEP/10.
GOTO 7777
ENDIF
8888 RR0=ROLD
BEQU=BOLD
BDEL=BDELTA
RETURN
END
!
!
SUBROUTINE SHELLG(GLAT,GLON,ALT,DIMO,FL,ICODE,B0,BEQU)
!--------------------------------------------------------------------
! CALCULATES L-VALUE FOR SPECIFIED GEODAETIC COORDINATES, ALTITUDE
! AND GEMAGNETIC FIELD MODEL.
! REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTER, INTERNAL NOTE
! NO. 67, 1970.
! G. KLUGE, COMPUTER PHYSICS COMMUNICATIONS 3, 31-35, 1972
!--------------------------------------------------------------------
! CHANGES (D. BILITZA, NOV 87):
! - USING CORRECT DIPOL MOMENT I.E.,DIFFERENT COMMON/MODEL/
! 09/07/22 NMAX=13 for DGRF00 and IGRF05; H/G-arrays(195)
!--------------------------------------------------------------------
! INPUT: ENTRY POINT SHELLG
! GLAT GEODETIC LATITUDE IN DEGREES (NORTH)
! GLON GEODETIC LONGITUDE IN DEGREES (EAST)
! ALT ALTITUDE IN KM ABOVE SEA LEVEL
!
! ENTRY POINT SHELLC
! V(3) CARTESIAN COORDINATES IN EARTH RADII (6371.2 KM)
! X-AXIS POINTING TO EQUATOR AT 0 LONGITUDE
! Y-AXIS POINTING TO EQUATOR AT 90 LONG.
! Z-AXIS POINTING TO NORTH POLE
!
! DIMO DIPOL MOMENT IN GAUSS (NORMALIZED TO EARTH RADIUS)
!
! COMMON
! X(3) NOT USED
! H(144) FIELD MODEL COEFFICIENTS ADJUSTED FOR SHELLG
!-----------------------------------------------------------------------
! OUTPUT: FL L-VALUE
! ICODE =1 NORMAL COMPLETION
! =2 UNPHYSICAL CONJUGATE POINT (FL MEANINGLESS)
! =3 SHELL PARAMETER GREATER THAN LIMIT UP TO
! WHICH ACCURATE CALCULATION IS REQUIRED;
! APPROXIMATION IS USED.
! B0 MAGNETIC FIELD STRENGTH IN GAUSS
!-----------------------------------------------------------------------
DIMENSION V(3),U(3,3),P(8,100),SP(3)
COMMON/IGRF2/ X(3),H(196)
COMMON/FIDB0/ SP
COMMON/GENER/ UMR,ERA,AQUAD,BQUAD
!
!-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3
!-- STEP IS STEP SIZE FOR FIELD LINE TRACING
!-- STEQ IS STEP SIZE FOR INTEGRATION
!
DATA RMIN,RMAX /0.05,1.01/
DATA STEP,STEQ /0.20,0.03/
BEQU=1.E10
!*****ENTRY POINT SHELLG TO BE USED WITH GEODETIC CO-ORDINATES
RLAT=GLAT*UMR
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)
GOTO 9
ENTRY SHELLC(V,FL,B0)
!*****ENTRY POINT SHELLC TO BE USED WITH CARTESIAN CO-ORDINATES
X(1)=V(1)
X(2)=V(2)
X(3)=V(3)
!*****CONVERT TO DIPOL-ORIENTED CO-ORDINATES
DATA U/ +0.3511737,-0.9148385,-0.1993679,
1+0.9335804,+0.3583680,+0.0000000, +0.0714471,-0.1861260,
1+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
!*****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)
!*****INVERT SENSE IF REQUIRED
IF(BQ3.LE.BQ1)GOTO 2
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
!*****SEARCH FOR LOWEST MAGNETIC FIELD STRENGTH
2 IF(BQ1.LT.BEQU) THEN
BEQU=BQ1
IEQU=1
ENDIF
IF(BQ2.LT.BEQU) THEN
BEQU=BQ2
IEQU=2
ENDIF
IF(BQ3.LT.BEQU) THEN
BEQU=BQ3
IEQU=3
ENDIF
!*****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
STP=STP/0.75
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))
!*****MAIN LOOP (FIELD LINE TRACING)
DO 3 N=3,3333
!*****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))
!*****PREPARE EXPANSION COEFFICIENTS FOR INTERPOLATION
!*****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)
D1=(P(6,N)-P(6,N-2))*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
!*****INNER LOOP (FOR QUADRATURE)
4 T=(Z-P(3,N-1))/STEP
IF(T.GT.1.)GOTO 5
HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)
ZQ=Z*Z
R=HLI+SQRT(HLI*HLI+ZQ)
IF(R.LE.RMIN) GOTO 30
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
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
GOTO 4
!*****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)
!*****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
SP(1)=P(1,IEQU-1)
SP(2)=P(2,IEQU-1)
SP(3)=P(3,IEQU-1)
IF(ORADIK.LT.1E-15) GOTO 11
FI=FI+STP/0.75*OTERM*ORADIK/(ORADIK-RADIK)
!
!-- The minimal allowable value of FI was changed from 1E-15 to 1E-12,
!-- because 1E-38 is the minimal allowable arg. for ALOG in our envir.
!-- D. Bilitza, Nov 87.
!
11 FI=0.5*ABS(FI)/SQRT(B0)+1E-12
!
!*****COMPUTE L FROM B AND I. SAME AS CARMEL IN INVAR.
!
!-- Correct dipole moment is used here. D. Bilitza, Nov 87.
!
DIMOB0=DIMO/B0
arg1=alog(FI)
arg2=alog(DIMOB0)
! arg = FI*FI*FI/DIMOB0
! 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
GOTO 777
772 GG=((((((((-8.1537735E-14*XX+8.3232531E-13)*XX+1.0066362E-9)*XX
1+8.1048663E-8)*XX+3.2916354E-6)*XX+8.2711096E-5)*XX+1.3714667E-3)
1*XX+1.5017245E-2)*XX+4.3432642E-1)*XX+6.2337691E-1
GOTO 777
773 GG=((((((((2.6047023E-10*XX+2.3028767E-9)*XX-2.1997983E-8)*XX-
15.3977642E-7)*XX-3.3408822E-6)*XX+3.8379917E-5)*XX+1.1784234E-3)*
1 XX+1.4492441E-2)*XX+4.3352788E-1)*XX+6.228644E-1
GOTO 777
774 GG=((((((((6.3271665E-10*XX-3.958306E-8)*XX+9.9766148E-07)*XX-
1 1.2531932E-5)*XX+7.9451313E-5)*XX-3.2077032E-4)*XX+2.1680398E-3)*
1 XX+1.2817956E-2)*XX+4.3510529E-1)*XX+6.222355E-1
GOTO 777
775 GG=(((((2.8212095E-8*XX-3.8049276E-6)*XX+2.170224E-4)*XX
1-6.7310339E-3)*XX+1.2038224E-1)*XX-1.8461796E-1)*XX+2.0007187E0
GOTO 777
776 GG=XX-3.0460681E0
777 FL=EXP(ALOG((1.+EXP(GG))*DIMOB0)/3.0)
RETURN
!*****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
!
!
SUBROUTINE STOER(P,BQ,R)
!*******************************************************************
!* SUBROUTINE USED FOR FIELD LINE TRACING IN SHELLG *
!* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG *
!
! 09/07/22 NMAX=13 for DGRF00 and IGRF05; H/G-arrays(195)
!*******************************************************************
DIMENSION P(7),U(3,3)
COMMON/IGRF2/ XI(3),H(196)
!*****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
!*****TRANSFORM TO GEOGRAPHIC CO-ORDINATE SYSTEM
DATA U/ +0.3511737,-0.9148385,-0.1993679,
1+0.9335804,+0.3583680,+0.0000000, +0.0714471,-0.1861260,
1+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)
!*****COMPUTE DERIVATIVES
! 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)
!*****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
!*****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
!
SUBROUTINE FELDG(GLAT,GLON,ALT,BNORTH,BEAST,BDOWN,BABS)
! Body of igrfExtraOut
!-------------------------------------------------------------------
! CALCULATES EARTH MAGNETIC FIELD FROM SPHERICAL HARMONICS MODEL
! REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTRE, INTERNAL NOTE 61,
! 1970.
!--------------------------------------------------------------------
! CHANGES (D. BILITZA, NOV 87):
! - FIELD COEFFICIENTS IN BINARY DATA FILES INSTEAD OF BLOCK DATA
! - CALCULATES DIPOL MOMENT
! 09/07/22 NMAX=13 for DGRF00 and IGRF05; H/G-arrays(195)
!--------------------------------------------------------------------
! INPUT: ENTRY POINT FELDG
! GLAT GEODETIC LATITUDE IN DEGREES (NORTH)
! GLON GEODETIC LONGITUDE IN DEGREES (EAST)
! ALT ALTITUDE IN KM ABOVE SEA LEVEL
!
! ENTRY POINT FELDC
! V(3) CARTESIAN COORDINATES IN EARTH RADII (6371.2 KM)
! X-AXIS POINTING TO EQUATOR AT 0 LONGITUDE
! Y-AXIS POINTING TO EQUATOR AT 90 LONG.
! Z-AXIS POINTING TO NORTH POLE
!
! COMMON BLANK AND ENTRY POINT FELDI ARE NEEDED WHEN USED
! IN CONNECTION WITH L-CALCULATION PROGRAM SHELLG.
!
! COMMON /MODEL/ AND /GENER/
! UMR = ATAN(1.0)*4./180. <DEGREE>*UMR=<RADIANT>
! ERA EARTH RADIUS FOR NORMALIZATION OF CARTESIAN
! COORDINATES (6371.2 KM)
! AQUAD, BQUAD SQUARE OF MAJOR AND MINOR HALF AXIS FOR
! EARTH ELLIPSOID AS RECOMMENDED BY INTERNATIONAL
! ASTRONOMICAL UNION (6378.160, 6356.775 KM).
! NMAX MAXIMUM ORDER OF SPHERICAL HARMONICS
! TIME YEAR (DECIMAL: 1973.5) FOR WHICH MAGNETIC
! FIELD IS TO BE CALCULATED
! G(M) NORMALIZED FIELD COEFFICIENTS (SEE FELDCOF)
! M=NMAX*(NMAX+2)
!------------------------------------------------------------------------
! OUTPUT: BABS MAGNETIC FIELD STRENGTH IN GAUSS
! BNORTH, BEAST, BDOWN COMPONENTS OF THE FIELD WITH RESPECT
! TO THE LOCAL GEODETIC COORDINATE SYSTEM, WITH AXIS
! POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
! AND DOWNWARD.
!-----------------------------------------------------------------------
DIMENSION V(3),B(3)
CHARACTER*13 NAME
COMMON/IGRF2/ XI(3),H(196)
COMMON/MODEL/ NAME,NMAX,TIME,G(196)
COMMON/GENER/ UMR,ERA,AQUAD,BQUAD
!
!-- IS RECORDS ENTRY POINT
!
!*****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)
RLON=GLON*UMR
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
ENTRY FELDC(V,B)
!*****ENTRY POINT FELDC TO BE USED WITH CARTESIAN CO-ORDINATES
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
ENTRY FELDI
!*****ENTRY POINT FELDI USED FOR L COMPUTATION
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)-H(IH+M-1))
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)-H(IH+M-2))
1+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
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
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
!
!
SUBROUTINE FELDCOF(YEAR,DIMO )
!------------------------------------------------------------------------
! DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS
!
! INPUT: YEAR DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO
! BE CALCULATED
! OUTPUT: DIMO GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED
! TO EARTH'S RADIUS) AT THE TIME (YEAR)
! D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771,
! (301)286-9536 NOV 1987.
! 05/31/2000 updated to IGRF-2000 version (###)
! 03/24/2000 updated to IGRF-2005 version (###)
!07/22/2009 NMAX=13 for DGRF00 and IGRF05; H/G-arrays(195)
! 02/26/2010 updated to IGRF-2010 version (###)
!-----------------------------------------------------------------------
CHARACTER*13 FILMOD, FIL1, FIL2
! ### FILMOD, DTEMOD array-size is number of IGRF maps
DIMENSION GH1(196),GH2(196),GHA(196),FILMOD(24),DTEMOD(24)
DOUBLE PRECISION X,F0,F
COMMON/MODEL/ FIL1,NMAX,TIME,GH1
COMMON/GENER/ UMR,ERAD,AQUAD,BQUAD
! ### updated coefficient file names and corresponding years
DATA FILMOD /'igrf1900.dat', 'igrf1905.dat', 'igrf1910.dat',
1 'igrf1915.dat', 'igrf1920.dat', 'igrf1925.dat', 'igrf1930.dat',
1 'igrf1935.dat', 'igrf1940.dat', 'dgrf1945.dat', 'dgrf1950.dat',
1 'dgrf1955.dat', 'dgrf1960.dat', 'dgrf1965.dat', 'dgrf1970.dat',
1 'dgrf1975.dat', 'dgrf1980.dat', 'dgrf1985.dat', 'dgrf1990.dat',
1'dgrf1995.dat', 'dgrf2000.dat', 'dgrf2005.dat', 'igrf2010.dat',
1 'igrf2010s.dat'/
DATA DTEMOD / 1900., 1905., 1910., 1915., 1920., 1925., 1930.,
1 1935., 1940., 1945., 1950., 1955., 1960., 1965., 1970., 1975.,
11980., 1985., 1990., 1995., 2000., 2005., 2010., 2015./
!
! ### numye is number of IGRF coefficient files minus 1
! ### istye is start year of IGRF coefficient files
!
NUMYE=23
ISTYE=1900
!
! IS=0 FOR SCHMIDT NORMALIZATION IS=1 GAUSS NORMALIZATION
! IU IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS
!
IU = 10
IS = 0
!-- DETERMINE IGRF-YEARS FOR INPUT-YEAR
TIME = YEAR
IYEA = INT(YEAR/5.)*5
L = (IYEA - ISTYE)/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)
!-- 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
!-- DETERMINE IGRF COEFFICIENTS FOR YEAR
IF (L .LE. NUMYE-1) THEN
CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2, NMAX2, GH2,
1 NMAX, GHA)
ELSE
CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2, GH2,
1 NMAX, GHA)
ENDIF
!-- DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G
F0=0.D0
DO 1234 J=1,3
F = GHA(J) * 1.D-5
F0 = F0 + F * F
1234 CONTINUE
DIMO = DSQRT(F0)
GH1(1) = 0.0
I=2
F0=1.D-5
IF(IS.EQ.0) F0=-F0
SQRT2=SQRT(2.)
DO 9 N=1,NMAX
X = N
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
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))
GH1(I) = GHA(I-1) * F
GH1(I+1) = GHA(I) * F
I=I+2
9 CONTINUE
RETURN
END
!
!
SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)
! ===============================================================
! Reads spherical harmonic coefficients from the specified
! file into an array.
! Input:
! IU - Logical unit number
! FSPEC - File specification
! Output:
! NMAX - Maximum degree and order of model
! ERAD - Earth's radius associated with the spherical
! harmonic coefficients, in the same units as
! elevation
! GH - Schmidt quasi-normal internal spherical
! harmonic coefficients
! IER - Error number: = 0, no error
! = -2, records out of order
! = FORTRAN run-time error number
! ===============================================================
CHARACTER FSPEC*(*), FOUT*80
DIMENSION GH(196)
do 1 j=1,196
1 GH(j)=0.0
! ---------------------------------------------------------------
! Open coefficient file. Read past first header record.
! Read degree and order of model and Earth's radius.
! ---------------------------------------------------------------
WRITE(FOUT,667) FSPEC
667 FORMAT(A13)
! 667 FORMAT('/var/www/omniweb/cgi/vitmo/IRI/',A13)
OPEN (IU, FILE=FOUT, 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 write(monito,100) FOUT
100 FORMAT('Error while reading ',A13)
888 CLOSE (IU)
RETURN
END
!
!
SUBROUTINE INTERSHC (DATE, DTE1, NMAX1, GH1, DTE2, NMAX2,
1GH2, NMAX, GH)
! ===============================================================
!
! Version 1.01
!
! Interpolates linearly, in time, between two spherical
! harmonic models.
!
! Input:
! DATE - Date of resulting model (in decimal year)
! DTE1 - Date of earlier model
! NMAX1 - Maximum degree and order of earlier model
! GH1 - Schmidt quasi-normal internal spherical
! harmonic coefficients of earlier model
! DTE2 - Date of later model
! NMAX2 - Maximum degree and order of later model
! GH2 - Schmidt quasi-normal internal spherical
! harmonic coefficients of later model
!
! Output:
! GH - Coefficients of resulting model
! NMAX - Maximum degree and order of resulting model
!
! A. Zunde
! USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225
!
! ===============================================================
DIMENSION GH1(*), GH2(*), GH(*)
! ---------------------------------------------------------------
! The coefficients (GH) of the resulting model, at date
! DATE, are computed by linearly interpolating between the
! coefficients of the earlier model (GH1), at date DTE1,
! and those of the later model (GH2), at date DTE2. If one
! model is smaller than the other, the interpolation is
! performed with the missing coefficients assumed to be 0.
! ---------------------------------------------------------------
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
!
!
SUBROUTINE EXTRASHC (DATE, DTE1, NMAX1, GH1, NMAX2,
1 GH2, NMAX, GH)
! ===============================================================
!
! Version 1.01
!
! Extrapolates linearly a spherical harmonic model with a
! rate-of-change model.
!
! Input:
! DATE - Date of resulting model (in decimal year)
! DTE1 - Date of base model
! NMAX1 - Maximum degree and order of base model
! GH1 - Schmidt quasi-normal internal spherical
! harmonic coefficients of base model
! NMAX2 - Maximum degree and order of rate-of-change
! model
! GH2 - Schmidt quasi-normal internal spherical
! harmonic coefficients of rate-of-change model
!
! Output:
! GH - Coefficients of resulting model
! NMAX - Maximum degree and order of resulting model
!
! A. Zunde
! USGS, MS 964, Box 25046 Federal Center, Denver, CO 80225
!
! ===============================================================
DIMENSION GH1(*), GH2(*), GH(*)
! ---------------------------------------------------------------
! The coefficients (GH) of the resulting model, at date
! DATE, are computed by linearly extrapolating the coef-
! ficients of the base model (GH1), at date DTE1, using
! those of the rate-of-change model (GH2), at date DTE2. If
! one model is smaller than the other, the extrapolation is
! performed with the missing coefficients assumed to be 0.
! ---------------------------------------------------------------
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
!
!
SUBROUTINE INITIZE
!----------------------------------------------------------------
! Initializes the parameters in COMMON/GENER/
!
! UMR = ATAN(1.0)*4./180. <DEGREE>*UMR=<RADIANT>
! ERA EARTH RADIUS FOR NORMALIZATION OF CARTESIAN
! COORDINATES (6371.2 KM)
! EREQU MAJOR HALF AXIS FOR EARTH ELLIPSOID (6378.160 KM)
! ERPOL MINOR HALF AXIS FOR EARTH ELLIPSOID (6356.775 KM)
! AQUAD SQUARE OF MAJOR HALF AXIS FOR EARTH ELLIPSOID
! BQUAD SQUARE OF MINOR HALF AXIS FOR EARTH ELLIPSOID
!
! ERA, EREQU and ERPOL as recommended by the INTERNATIONAL
! ASTRONOMICAL UNION .
!-----------------------------------------------------------------
COMMON/GENER/ UMR,ERA,AQUAD,BQUAD
ERA=6371.2
EREQU=6378.16
ERPOL=6356.775
AQUAD=EREQU*EREQU
BQUAD=ERPOL*ERPOL
UMR=ATAN(1.0)*4./180.
RETURN
END
以上是两个平台下分别正确的代码，标红部分不一样。
linux标红部分如果和windows平台下一样，用gfortran编译器编译时，则报错：
gfortran -shared -fPIC -o ap8ae8.so ap8ae8.for