/[PAMELA software]/DarthVader/OrbitalInfo/src/igrf_sub.for
ViewVC logotype

Diff of /DarthVader/OrbitalInfo/src/igrf_sub.for

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.1 by mocchiut, Tue Jan 23 17:03:39 2007 UTC revision 1.8 by mocchiut, Thu Jan 16 15:29:26 2014 UTC
# Line 1  Line 1 
1          subroutine igrf_sub(xlat,xlong,year,height,  c        subroutine igrf_sub(xlat,xlong,year,height,
2       &          xl,icode,dip,dec)  c     &          xl,icode,dip,dec)
3  c----------------------------------------------------------------  c----------------------------------------------------------------
4  c   INPUT:  c   INPUT:
5  c       xlat    geodatic latitude in degrees  c     xlat     geodatic latitude in degrees
6  c       xlong   geodatic longitude in degrees  c     xlong    geodatic longitude in degrees
7  c       year    decimal year (year+month/12.0-0.5 or year+day-of-year/365  c     year     decimal year (year+month/12.0-0.5 or year+day-of-year/365
8  c               or 366 if leap year)  c          or 366 if leap year)
9  c       height  height in km  c     height     height in km
10  c   OUTPUT:  c   OUTPUT:
11  c       xl      L value  c     xl     L value
12  c       icode   =1  L is correct; =2  L is not correct;  c     icode     =1  L is correct; =2  L is not correct;
13  c               =3  an approximation is used  c          =3  an approximation is used
14  c       dip     geomagnetic inclination in degrees  c     dip     geomagnetic inclination in degrees
15  c       dec     geomagnetic declination in degress  c     dec     geomagnetic declination in degress
16  c----------------------------------------------------------------  c----------------------------------------------------------------
17    c
18        REAL              LATI,LONGI  c      REAL              LATI,LONGI
19        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD  c      COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
20        SAVE /GENER/  c      SAVE /GENER/
21  C  C
22        CALL INITIZE  c      CALL INITIZE
23          ibbb=0  c      ibbb=0
24        ALOG2=ALOG(2.)  c      ALOG2=ALOG(2.)
25        ISTART=1  c      ISTART=1
26          lati=xlat  c        lati=xlat
27          longi=xlong  c        longi=xlong
28  c  c
29  C----------------CALCULATE PROFILES-----------------------------------  C----------------CALCULATE PROFILES-----------------------------------
30  c  c
31          CALL FELDCOF(YEAR,DIMO)  c        CALL FELDCOF(YEAR,DIMO)
32          CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS)  c        CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS)
33          CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1)  c        CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1)
34          DIP=ASIN(BDOWN/BABS)/UMR  c        DIP=ASIN(BDOWN/BABS)/UMR
35          DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR  c      DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR
36        RETURN  c      RETURN
37        END  c      END
38  c  c
39  c  c
40  C SHELLIG.FOR, Version 2.0, January 1992  C SHELLIG.FOR, Version 2.0, January 1992
# Line 55  C*************************************** Line 55  C***************************************
55  C  C
56  C  C
57        SUBROUTINE FINDB0(STPS,BDEL,VALUE,BEQU,RR0)        SUBROUTINE FINDB0(STPS,BDEL,VALUE,BEQU,RR0)
58          IMPLICIT REAL(8)(A-H)
59          IMPLICIT REAL(8)(O-Z)
60  C--------------------------------------------------------------------  C--------------------------------------------------------------------
61  C FINDS SMALLEST MAGNETIC FIELD STRENGTH ON FIELD LINE  C FINDS SMALLEST MAGNETIC FIELD STRENGTH ON FIELD LINE
62  C  C
# Line 70  C          BEQU   MAGNETIC FIELD STRENGT Line 72  C          BEQU   MAGNETIC FIELD STRENGT
72  C          RR0    EQUATORIAL RADIUS NORMALIZED TO EARTH RADIUS  C          RR0    EQUATORIAL RADIUS NORMALIZED TO EARTH RADIUS
73  C          BDEL   FINAL ACHIEVED ACCURACY  C          BDEL   FINAL ACHIEVED ACCURACY
74  C--------------------------------------------------------------------  C--------------------------------------------------------------------
75          REAL(8) P
76        DIMENSION         P(8,4),SP(3)        DIMENSION         P(8,4),SP(3)
77        LOGICAL           VALUE        LOGICAL           VALUE
78        COMMON/FIDB0/     SP        COMMON/FIDB0/     SP
# Line 154  C        DO 1111 I=1,8 Line 157  C        DO 1111 I=1,8
157  C  C
158  C  C
159        SUBROUTINE SHELLG(GLAT,GLON,ALT,DIMO,FL,ICODE,B0)        SUBROUTINE SHELLG(GLAT,GLON,ALT,DIMO,FL,ICODE,B0)
160          IMPLICIT REAL(8)(A-H)
161          IMPLICIT REAL(8)(O-Z)
162  C--------------------------------------------------------------------  C--------------------------------------------------------------------
163  C CALCULATES L-VALUE FOR SPECIFIED GEODAETIC COORDINATES, ALTITUDE  C CALCULATES L-VALUE FOR SPECIFIED GEODAETIC COORDINATES, ALTITUDE
164  C AND GEMAGNETIC FIELD MODEL.  C AND GEMAGNETIC FIELD MODEL.
# Line 190  C                          WHICH ACCURAT Line 195  C                          WHICH ACCURAT
195  C                          APPROXIMATION IS USED.  C                          APPROXIMATION IS USED.
196  C          B0           MAGNETIC FIELD STRENGTH IN GAUSS  C          B0           MAGNETIC FIELD STRENGTH IN GAUSS
197  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
198          REAL(8) AQUAD,BQUAD,ERA
199          REAL(8) CT,ST,D,ALT,X,RQ,P,U,STEP
200        DIMENSION         V(3),U(3,3),P(8,100),SP(3)        DIMENSION         V(3),U(3,3),P(8,100),SP(3)
201        COMMON            X(3),H(144)        COMMON            X(3),H(196)
202        COMMON/FIDB0/     SP        COMMON/FIDB0/     SP
203        SAVE /FIDB0/        SAVE /FIDB0/
204        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     ERA,AQUAD,BQUAD,UMR
205        SAVE /GENER/        SAVE /GENER/
206  C  C
207  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3
# Line 308  C*****INNER LOOP (FOR QUADRATURE)       Line 315  C*****INNER LOOP (FOR QUADRATURE)      
315        RQ=R*R        RQ=R*R
316        FF=SQRT(1.+3.*ZQ/RQ)                                      FF=SQRT(1.+3.*ZQ/RQ)                              
317        RADIK=B0-((D2*T+D1)*T+D0)*R*RQ*FF                        RADIK=B0-((D2*T+D1)*T+D0)*R*RQ*FF                
318        IF(R-RMAX)44,44,45                                        IF((R-RMAX).le.0.) goto 44                                
319          IF((R-RMAX).gt.0.) goto 45                                
320  45    ICODE=2                                            45    ICODE=2                                          
321        RADIK=RADIK-12.*(R-RMAX)**2                              RADIK=RADIK-12.*(R-RMAX)**2                      
322  44    IF(RADIK+RADIK.LE.ORADIK) GOTO 10  44    IF(RADIK+RADIK.LE.ORADIK) GOTO 10
# Line 341  C-- The minimal allowable value of FI wa Line 349  C-- The minimal allowable value of FI wa
349  C-- because 1E-38 is the minimal allowable arg. for ALOG in our envir.  C-- because 1E-38 is the minimal allowable arg. for ALOG in our envir.
350  C-- D. Bilitza, Nov 87.  C-- D. Bilitza, Nov 87.
351  C  C
352  11      FI=0.5*ABS(FI)/SQRT(B0)+1E-12                        11         FI=0.5*ABS(FI)/SQRT(B0)+1E-12                      
353  C  C
354  C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.    C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.  
355  C  C
356  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.
357  C  C
358        DIMOB0=DIMO/B0        DIMOB0=DIMO/B0
359          arg1=alog(FI)        arg1=dlog(FI)
360          arg2=alog(DIMOB0)        arg2=dlog(DIMOB0)
361  c       arg = FI*FI*FI/DIMOB0  c    arg = FI*FI*FI/DIMOB0
362  c       if(abs(arg).gt.88.0) arg=88.0  c    if(abs(arg).gt.88.0) arg=88.0
363        XX=3*arg1-arg2        XX=3*arg1-arg2
364        IF(XX.GT.23.0) GOTO 776          IF(XX.GT.23.0) GOTO 776  
365        IF(XX.GT.11.7) GOTO 775          IF(XX.GT.11.7) GOTO 775  
366        IF(XX.GT.+3.0) GOTO 774            IF(XX.GT.+3.0) GOTO 774    
367        IF(XX.GT.-3.0) GOTO 773          IF(XX.GT.-3.0) GOTO 773  
368        IF(XX.GT.-22.) GOTO 772          IF(XX.GT.-22.) GOTO 772  
369    771 GG=3.33338E-1*XX+3.0062102E-1                                  c  771 GG=3.33338E-1*XX+3.0062102E-1                                
370          GG=3.33338E-1*XX+3.0062102E-1                                
371        GOTO777                                                                  GOTO777                                                          
372    772 GG=((((((((-8.1537735E-14*XX+8.3232531E-13)*XX+1.0066362E-9)*XX+      772 GG=((((((((-8.1537735E-14*XX+8.3232531E-13)*XX+1.0066362E-9)*XX+  
373       18.1048663E-8)*XX+3.2916354E-6)*XX+8.2711096E-5)*XX+1.3714667E-3)*       18.1048663E-8)*XX+3.2916354E-6)*XX+8.2711096E-5)*XX+1.3714667E-3)*
# Line 376  c      if(abs(arg).gt.88.0) arg=88.0 Line 385  c      if(abs(arg).gt.88.0) arg=88.0
385       1E-3)*XX+1.2038224E-1)*XX-1.8461796E-1)*XX+2.0007187E0                   1E-3)*XX+1.2038224E-1)*XX-1.8461796E-1)*XX+2.0007187E0            
386        GOTO777                                                                  GOTO777                                                          
387    776 GG=XX-3.0460681E0                                                    776 GG=XX-3.0460681E0                                                
388    777 FL=EXP(ALOG((1.+EXP(GG))*DIMOB0)/3.0)    777 FL=EXP(dLOG((1.+EXP(GG))*DIMOB0)/3.0)
389        RETURN                                                                    RETURN                                                            
390  C*****APPROXIMATION FOR HIGH VALUES OF L.                                C*****APPROXIMATION FOR HIGH VALUES OF L.                              
391  30    ICODE=3                                                            30    ICODE=3                                                          
# Line 387  C*****APPROXIMATION FOR HIGH VALUES OF L Line 396  C*****APPROXIMATION FOR HIGH VALUES OF L
396  C  C
397  C  C
398        SUBROUTINE STOER(P,BQ,R)                                                  SUBROUTINE STOER(P,BQ,R)                                          
399          IMPLICIT REAL(8)(A-H)
400          IMPLICIT REAL(8)(O-Z)
401    
402  C*******************************************************************  C*******************************************************************
403  C* SUBROUTINE USED FOR FIELD LINE TRACING IN SHELLG                *  C* SUBROUTINE USED FOR FIELD LINE TRACING IN SHELLG                *
404  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *
405  C*******************************************************************  C*******************************************************************
406          REAL(8) P,ZM,FLI,WR,XM,R,YM,XI,DR
407        DIMENSION         P(7),U(3,3)        DIMENSION         P(7),U(3,3)
408        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
409  C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES            C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES          
410        ZM=P(3)                                                                  ZM=P(3)                                                          
411        FLI=P(1)*P(1)+P(2)*P(2)+1E-15        FLI=P(1)*P(1)+P(2)*P(2)+1E-15
# Line 432  C*****FORM SLOWLY VARYING EXPRESSIONS   Line 445  C*****FORM SLOWLY VARYING EXPRESSIONS  
445  C  C
446  C  C
447        SUBROUTINE FELDG(GLAT,GLON,ALT,BNORTH,BEAST,BDOWN,BABS)                  SUBROUTINE FELDG(GLAT,GLON,ALT,BNORTH,BEAST,BDOWN,BABS)          
448          IMPLICIT REAL(8)(A-H)
449          IMPLICIT REAL(8)(O-Z)
450  C-------------------------------------------------------------------  C-------------------------------------------------------------------
451  C CALCULATES EARTH MAGNETIC FIELD FROM SPHERICAL HARMONICS MODEL  C CALCULATES EARTH MAGNETIC FIELD FROM SPHERICAL HARMONICS MODEL
452  C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTRE, INTERNAL NOTE 61,  C REF: G. KLUGE, EUROPEAN SPACE OPERATIONS CENTRE, INTERNAL NOTE 61,
# Line 474  C                 TO THE LOCAL GEODETIC Line 489  C                 TO THE LOCAL GEODETIC
489  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
490  C                 AND DOWNWARD.    C                 AND DOWNWARD.  
491  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
492        DIMENSION         V(3),B(3),G(144)        REAL(8) G,H,X,Z,F,S,T,XI,Y,XXX,BXXX,BYYY,BZZZ,B
493          REAL(8) YYY,ZZZ,BABS,BEAST,BRHO,BNORTH,BDOWN,ST,CT,CP,SP
494          REAL(8) ERA, AQUAD, BQUAD,ALT,RHO,RQ,D
495          DIMENSION         V(3),B(3),G(196)
496        CHARACTER*258      NAME        CHARACTER*258      NAME
497        INTEGER NMAX        INTEGER NMAX
498        REAL TIME        REAL TIME
499        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
500          
501        COMMON/MODEL/     G,NMAX,TIME,NAME        COMMON/MODEL/     G,NMAX,TIME,NAME
502        SAVE/MODEL/        SAVE/MODEL/
503        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     ERA,AQUAD,BQUAD,UMR
504        SAVE/GENER/        SAVE/GENER/
505  C  C
506  C-- IS RECORDS ENTRY POINT  C-- IS RECORDS ENTRY POINT
# Line 528  C*****ENTRY POINT  FELDI  USED FOR L COM Line 547  C*****ENTRY POINT  FELDI  USED FOR L COM
547        Y=XI(2)*F                                                                Y=XI(2)*F                                                        
548        Z=XI(3)*(F+F)                                                            Z=XI(3)*(F+F)                                                    
549        I=I-2                                                                    I=I-2                                                            
550        IF(I-1)5,4,2                                                        c      print *,' I ',I
551          IF((I-1).lt.0) goto 5
552          IF((I-1).eq.0) goto 4
553          IF((I-1).gt.0) goto 2
554  2     DO 3 M=3,I,2                                                        2     DO 3 M=3,I,2                                                      
555        H(IL+M+1)=G(IL+M+1)+Z*H(IH+M+1)+X*(H(IH+M+3)-H(IH+M-1))                  H(IL+M+1)=G(IL+M+1)+Z*H(IH+M+1)+X*(H(IH+M+3)-H(IH+M-1))          
556       A                               -Y*(H(IH+M+2)+H(IH+M-2))                 A                               -Y*(H(IH+M+2)+H(IH+M-2))          
# Line 560  C*****ENTRY POINT  FELDI  USED FOR L COM Line 582  C*****ENTRY POINT  FELDI  USED FOR L COM
582        END                                                                      END                                                              
583  C  C
584  C  C
585          SUBROUTINE FELDCOF(YEAR,DIMO)        SUBROUTINE FELDCOF(YEAR,DIMO)
586          IMPLICIT REAL(8)(A-H)
587          IMPLICIT REAL(8)(O-Z)
588  C------------------------------------------------------------------------  C------------------------------------------------------------------------
589  C  DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS  C     DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS
590  C  C    
591  C       INPUT:  YEAR    DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO  C     INPUT:  YEAR    DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO
592  C                       BE CALCULATED  C     BE CALCULATED
593  C       OUTPUT: DIMO    GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED  C     OUTPUT: DIMO    GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED
594  C                       TO EARTH'S RADIUS) AT THE TIME (YEAR)  C     TO EARTH'S RADIUS) AT THE TIME (YEAR)
595  C  D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771,  C     D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771,
596  C       (301)286-9536   NOV 1987.  C     (301)286-9536   NOV 1987.
597  C  ### updated to IGRF-2000 version -dkb- 5/31/2000  C     ### updated to IGRF-2000 version -dkb- 5/31/2000
598  C  ### updated to IGRF-2005 version -dkb- 3/24/2000  C     ### updated to IGRF-2005 version -dkb- 3/24/2000
599  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
600          CHARACTER*258    FIL1, FIL2                  CHARACTER*258    FIL1, FIL2          
601          CHARACTER*258    FILMOD        CHARACTER*258    FILMOD
602  C ### FILMOD, DTEMOD arrays +1  C     ### FILMOD, DTEMOD arrays +1
603  c        DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)  c     DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)
604          DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(3),DTEMOD(3)        REAL(8) GH1, GH2, GHA
605          DOUBLE PRECISION X,F0,F        DIMENSION       GH1(196),GH2(196),GHA(196),FILMOD(3),DTEMOD(3)
606          INTEGER L1,L2,L3        DOUBLE PRECISION X,F0,F
607          INTEGER NMAX        DOUBLE PRECISION DIMO
608          REAL TIME        INTEGER L1,L2,L3
609          CHARACTER *258 P1,P2,P3        INTEGER NMAX
610          COMMON/PPATH/ L1,L2,L3,P1, P2, P3        REAL YEAR
611          SAVE/PPATH/        REAL TIME
612          COMMON/MODEL/   GH1,NMAX,TIME,FIL1        CHARACTER *258 P1,P2,P3
613          SAVE/MODEL/        REAL(8) AQUAD,BQUAD,ERAD
614          COMMON/GENER/   UMR,ERAD,AQUAD,BQUAD        COMMON/PPATH/ L1,L2,L3,P1, P2, P3
615          SAVE/GENER/        SAVE/PPATH/
616  C ### updated to 2005        COMMON/MODEL/   GH1,NMAX,TIME,FIL1
617  C      CHARACTER COEFPATH*80, COEF1*80, COEF2*80, COEF3*80        SAVE/MODEL/
618          COMMON/GENER/   ERAD,AQUAD,BQUAD,UMR
619  c      COEFPATH  = 'OrbitalInfo/src/'        SAVE/GENER/
620  c      COEF1 = 'dgrf00.dat'  C     ### updated to 2005
621  c      COEF2 = 'igrf05.dat'  C     CHARACTER COEFPATH*80, COEF1*80, COEF2*80, COEF3*80
622  c      COEF3 = 'igrf05s.dat'        
623  c      COEF1 = COEFPATH(1:16)//COEF1  c     COEFPATH  = 'OrbitalInfo/src/'
624  c      COEF2 = COEFPATH(1:16)//COEF2  c     COEF1 = 'dgrf00.dat'
625  c      COEF3 = COEFPATH(1:16)//COEF3  c     COEF2 = 'igrf05.dat'
626  c      FILMOD(1) = COEF1  c     COEF3 = 'igrf05s.dat'
627  c      FILMOD(2) = COEF2  c     COEF1 = COEFPATH(1:16)//COEF1
628  c      FILMOD(3) = COEF3  c     COEF2 = COEFPATH(1:16)//COEF2
629    c     COEF3 = COEFPATH(1:16)//COEF3
630    c     FILMOD(1) = COEF1
631    c     FILMOD(2) = COEF2
632    c     FILMOD(3) = COEF3
633    c      print *, "qui"
634        FILMOD(1) = P1(1:L1)        FILMOD(1) = P1(1:L1)
635        FILMOD(2) = P2(1:L2)        FILMOD(2) = P2(1:L2)
636        FILMOD(3) = P3(1:L3)        FILMOD(3) = P3(1:L3)
637  c      FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'  c      print *, "qua"
638  c      FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'
639  c      FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'
640    c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'
641  c      WRITE(*,*) FILMOD(1)  c      WRITE(*,*) FILMOD(1)
642  c      WRITE(*,*) FILMOD(2)  c      WRITE(*,*) FILMOD(2)
643  c      WRITE(*,*) FILMOD(3)  c      WRITE(*,*) FILMOD(3)
644  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/
645        DATA   DTEMOD / 2000., 2005., 2010./        DATA   DTEMOD / 2005., 2010., 2015./
646  c  c    
647  c        DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',              c     DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',            
648  c     1                  'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',        c     1                  'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',      
649  c     2                  'dgrf70.dat', 'dgrf75.dat', 'dgrf80.dat',        c     2                  'dgrf70.dat', 'dgrf75.dat', 'dgrf80.dat',      
650  c     3                  'dgrf85.dat', 'dgrf90.dat', 'dgrf95.dat',  c     3                  'dgrf85.dat', 'dgrf90.dat', 'dgrf95.dat',
651  c     4                  'dgrf00.dat','igrf05.dat','igrf05s.dat'/  c     4                  'dgrf00.dat','igrf05.dat','igrf05s.dat'/
652  c        DATA   DTEMOD / 1945., 1950., 1955., 1960., 1965., 1970.,  c     DATA   DTEMOD / 1945., 1950., 1955., 1960., 1965., 1970.,
653  c     1       1975., 1980., 1985., 1990., 1995., 2000.,2005.,2010./        c     1       1975., 1980., 1985., 1990., 1995., 2000.,2005.,2010./      
654  C  C    
655  C ### numye = numye + 1 ; is number of years represented by IGRF  C     ### numye = numye + 1 ; is number of years represented by IGRF
656  C  C    
657  c        NUMYE=13  c     NUMYE=13
658          NUMYE=2        NUMYE=2
659    c      print *, "quo"
660  C        
661  C  IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION  C    
662  C  IU  IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS  C     IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION
663  C  C     IU  IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS
664          IU = 10  C    
665          IS = 0        IU = 10
666  C-- DETERMINE IGRF-YEARS FOR INPUT-YEAR        IS = 0
667          TIME = YEAR  C--   DETERMINE IGRF-YEARS FOR INPUT-YEAR
668          IYEA = INT(YEAR/5.)*5        TIME = YEAR
669  c        L = (IYEA - 1945)/5 + 1        IYEA = INT(YEAR/5.)*5
670          L = (IYEA - 2000)/5 + 1  c     L = (IYEA - 1945)/5 + 1
671          IF(L.LT.1) L=1        L = (IYEA - 2000)/5 + 1
672          IF(L.GT.NUMYE) L=NUMYE                IF(L.LT.1) L=1
673          DTE1 = DTEMOD(L)          IF(L.GT.NUMYE) L=NUMYE        
674          FIL1 = FILMOD(L)          DTE1 = DTEMOD(L)  
675          DTE2 = DTEMOD(L+1)        FIL1 = FILMOD(L)  
676          FIL2 = FILMOD(L+1)        DTE2 = DTEMOD(L+1)
677  C-- GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS        FIL2 = FILMOD(L+1)
678          CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)    c      WRITE(*,*) FIL1
679              IF (IER .NE. 0) STOP                            c      WRITE(*,*) FIL2
680          CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)    c      print *, "que"
681              IF (IER .NE. 0) STOP                      C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
682  C-- DETERMINE IGRF COEFFICIENTS FOR YEAR        CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
683          IF (L .LE. NUMYE-1) THEN                                IF (IER .NE. 0) STOP                          
684            CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,  c      print *, "quessss"
685       1          NMAX2, GH2, NMAX, GHA)                                CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
686          ELSE                      IF (IER .NE. 0) STOP                    
687            CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,      c      print *, "quj"
688       1          GH2, NMAX, GHA)                                      C--   DETERMINE IGRF COEFFICIENTS FOR YEAR
689          ENDIF        IF (L .LE. NUMYE-1) THEN                        
690  C-- DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G           CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,
691          F0=0.D0       1        NMAX2, GH2, NMAX, GHA)                        
692          DO 1234 J=1,3        ELSE              
693             F = GHA(J) * 1.D-5           CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,    
694             F0 = F0 + F * F       1        GH2, NMAX, GHA)                                    
695  1234    CONTINUE        ENDIF
696          DIMO = DSQRT(F0)  c      print *, "quw"
697    C--   DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G
698          GH1(1) =  0.0        F0=0.D0
699          I=2                  DO 1234 J=1,3
700          F0=1.D-5                           F = GHA(J) * 1.D-5
701          IF(IS.EQ.0) F0=-F0           F0 = F0 + F * F
702          SQRT2=SQRT(2.)         1234 CONTINUE
703          DIMO = DSQRT(F0)
704          
705          GH1(1) =  0.0
706          I=2          
707          F0=1.D-5                
708          IF(IS.EQ.0) F0=-F0
709          SQRT2=SQRT(2.)      
710    
711    c      print *, "quq"
712          
713        DO 9 N=1,NMAX                  DO 9 N=1,NMAX          
714          X = N           X = N
715          F0 = F0 * X * X / (4.D0 * X - 2.D0)                         F0 = F0 * X * X / (4.D0 * X - 2.D0)              
716          IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X
717          F = F0 * 0.5D0                                               F = F0 * 0.5D0                                    
718          IF(IS.EQ.0) F = F * SQRT2           IF(IS.EQ.0) F = F * SQRT2
719          GH1(I) = GHA(I-1) * F0           GH1(I) = GHA(I-1) * F0
720          I = I+1                                                   I = I+1                                        
721        DO 9 M=1,N                                               DO 9 M=1,N                                    
722          F = F * (X + M) / (X - M + 1.D0)                              F = F * (X + M) / (X - M + 1.D0)                
723          IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M))                          IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M))
724          GH1(I) = GHA(I-1) * F              GH1(I) = GHA(I-1) * F
725          GH1(I+1) = GHA(I) * F              GH1(I+1) = GHA(I) * F
726          I=I+2              I=I+2
727  9     CONTINUE                                             9       CONTINUE                                          
728          RETURN        RETURN
729          END        END
730  C  C    
731  C  C    
732          SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)                  SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)          
733                                                                                          IMPLICIT REAL(8)(A-H)
734          IMPLICIT REAL(8)(O-Z)
735  C ===============================================================                C ===============================================================              
736  C                                                                                C                                                                              
737  C       Version 1.01                                                  C     Version 1.01                                                
738  C                                                                                C    
739  C       Reads spherical harmonic coefficients from the specified      C     Reads spherical harmonic coefficients from the specified    
740  C       file into an array.                                            C     file into an array.                                          
741  C                                                                                C    
742  C       Input:                                                        C     Input:                                                      
743  C           IU    - Logical unit number                                C     IU    - Logical unit number                              
744  C           FSPEC - File specification                                C     FSPEC - File specification                              
745  C                                                                                C    
746  C       Output:                                                        C     Output:                                                      
747  C           NMAX  - Maximum degree and order of model                  C     NMAX  - Maximum degree and order of model                
748  C           ERAD  - Earth's radius associated with the spherical      C     ERAD  - Earth's radius associated with the spherical    
749  C                   harmonic coefficients, in the same units as        C     harmonic coefficients, in the same units as      
750  C                   elevation                                          C     elevation                                        
751  C           GH    - Schmidt quasi-normal internal spherical            C     GH    - Schmidt quasi-normal internal spherical          
752  C                   harmonic coefficients                              C     harmonic coefficients                            
753  C           IER   - Error number: =  0, no error                      C     IER   - Error number: =  0, no error                    
754  C                                 = -2, records out of order          C     = -2, records out of order        
755  C                                 = FORTRAN run-time error number      C     = FORTRAN run-time error number    
756  C                                                                      C    
757  C       A. Zunde                                                      C     A. Zunde                                                    
758  C       USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225      C     USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225    
759  C                                                                                C    
760  C ===============================================================                C     ===============================================================              
761                                                                                          
762          CHARACTER  FSPEC*(*), FOUT*258                                            CHARACTER  FSPEC*(*), FOUT*258
763          DIMENSION       GH(*)                                                REAL(8) GH,ERAD
764  C ---------------------------------------------------------------                      DIMENSION       GH(*)                                        
765  C       Open coefficient file. Read past first header record.          C     ---------------------------------------------------------------              
766  C       Read degree and order of model and Earth's radius.            C     Open coefficient file. Read past first header record.        
767  C ---------------------------------------------------------------                C     Read degree and order of model and Earth's radius.          
768    C     ---------------------------------------------------------------              
769        WRITE(FOUT,667) FSPEC        WRITE(FOUT,667) FSPEC
770  c 667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)  c     667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)
771  667   FORMAT(A258)   667  FORMAT(A258)
772    c      print *," gui"
773        OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)            OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)    
774    c      print *," gua"
775        READ (IU, *, IOSTAT=IER, ERR=999)        READ (IU, *, IOSTAT=IER, ERR=999)
776    c      print *," gue"
777        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD
778    c      print *," guo"
779  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
780  C       Read the coefficient file, arranged as follows:                C       Read the coefficient file, arranged as follows:              
781  C                                                                                C                                                                              
# Line 752  C       coefficient.                     Line 797  C       coefficient.                    
797  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
798                                                                                                                                                                    
799          I = 0                                                                  I = 0                                                        
800          DO 2211 NN = 1, NMAX                                                        DO 2211 NN = 1, NMAX              
801              DO 2233 MM = 0, NN                                                          DO 2233 MM = 0, NN            
802                  READ (IU, *, IOSTAT=IER, ERR=999) N, M, G, H                          READ (IU, *, IOSTAT=IER, ERR=999) N, M, G, H        
803                  IF (NN .NE. N .OR. MM .NE. M) THEN                                    IF (NN .NE. N .OR. MM .NE. M) THEN                  
804                      IER = -2                                                              IER = -2                                        
# Line 767  C -------------------------------------- Line 812  C --------------------------------------
812                  ENDIF                                                                  ENDIF                                                
813  2233        CONTINUE                                                      2233        CONTINUE                                                    
814  2211    CONTINUE                                                          2211    CONTINUE                                                        
815    c      print *," guj"
816                                                                                                                                                                    
817  999     CLOSE (IU)                                                    999     CLOSE (IU)                                                  
818    c      print *," guw IER",IER
819          if ( IER .eq. -1 ) IER = 0 ! gfortran 4.1.2 bug workaround... hoping not to create problems with other versions
820                                                                                                                                                                    
821          RETURN                                                                RETURN                                                      
822          END                                                                    END                                                          
# Line 777  C Line 824  C
824  C  C
825          SUBROUTINE INTERSHC (DATE, DTE1, NMAX1, GH1, DTE2,                    SUBROUTINE INTERSHC (DATE, DTE1, NMAX1, GH1, DTE2,          
826       1                        NMAX2, GH2, NMAX, GH)                         1                        NMAX2, GH2, NMAX, GH)                  
827                                                                                            IMPLICIT REAL(8)(A-H)
828            IMPLICIT REAL(8)(O-Z)
829            REAL DATE
830  C ===============================================================                C ===============================================================              
831  C                                                                                C                                                                              
832  C       Version 1.01                                                  C       Version 1.01                                                
# Line 805  C       USGS, MS 964, Box 25046 Federal Line 854  C       USGS, MS 964, Box 25046 Federal
854  C                                                                                C                                                                              
855  C ===============================================================                C ===============================================================              
856                                                                                                                                                                    
857            REAL(8) GH1, GH2, GH
858          DIMENSION       GH1(*), GH2(*), GH(*)                                  DIMENSION       GH1(*), GH2(*), GH(*)                        
859                                                                                                                                                                    
860  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
# Line 824  C -------------------------------------- Line 874  C --------------------------------------
874          ELSE IF (NMAX1 .GT. NMAX2) THEN                                        ELSE IF (NMAX1 .GT. NMAX2) THEN                              
875              K = NMAX2 * (NMAX2 + 2)                                                K = NMAX2 * (NMAX2 + 2)                                  
876              L = NMAX1 * (NMAX1 + 2)                                                L = NMAX1 * (NMAX1 + 2)                                  
877              DO 1122 I = K + 1, L                                                        DO 1122 I = K + 1, L    
878  1122            GH(I) = GH1(I) + FACTOR * (-GH1(I))                    1122            GH(I) = GH1(I) + FACTOR * (-GH1(I))                  
879              NMAX = NMAX1                                                          NMAX = NMAX1                                            
880          ELSE                                                                  ELSE                                                        
881              K = NMAX1 * (NMAX1 + 2)                                                K = NMAX1 * (NMAX1 + 2)                                  
882              L = NMAX2 * (NMAX2 + 2)                                                L = NMAX2 * (NMAX2 + 2)                                  
883              DO 1133 I = K + 1, L                                                        DO 1133 I = K + 1, L
884  1133            GH(I) = FACTOR * GH2(I)                                1133            GH(I) = FACTOR * GH2(I)                              
885              NMAX = NMAX2                                                          NMAX = NMAX2                                            
886          ENDIF                                                                  ENDIF                                                        
887                                                                                                                                                                    
888          DO 1144 I = 1, K                                                            DO 1144 I = 1, K  
889  1144        GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I))                1144        GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I))              
890                                                                                                                                                                    
891          RETURN                                                                RETURN                                                      
# Line 844  C Line 894  C
894  C  C
895          SUBROUTINE EXTRASHC (DATE, DTE1, NMAX1, GH1, NMAX2,                    SUBROUTINE EXTRASHC (DATE, DTE1, NMAX1, GH1, NMAX2,          
896       1                        GH2, NMAX, GH)                                 1                        GH2, NMAX, GH)                          
897            IMPLICIT REAL(8)(A-H)
898            IMPLICIT REAL(8)(O-Z)
899            REAL DATE
900                                                                                                                                                                    
901  C ===============================================================                C ===============================================================              
902  C                                                                                C                                                                              
# Line 872  C       USGS, MS 964, Box 25046 Federal Line 925  C       USGS, MS 964, Box 25046 Federal
925  C                                                                                C                                                                              
926  C ===============================================================                C ===============================================================              
927                                                                                                                                                                    
928            REAL(8) GH1, GH2, GH
929          DIMENSION       GH1(*), GH2(*), GH(*)                                    DIMENSION       GH1(*), GH2(*), GH(*)                          
930                                                                                                                                                                    
931  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
# Line 891  C -------------------------------------- Line 945  C --------------------------------------
945          ELSE IF (NMAX1 .GT. NMAX2) THEN                                          ELSE IF (NMAX1 .GT. NMAX2) THEN                                
946              K = NMAX2 * (NMAX2 + 2)                                                  K = NMAX2 * (NMAX2 + 2)                                    
947              L = NMAX1 * (NMAX1 + 2)                                                  L = NMAX1 * (NMAX1 + 2)                                    
948              DO 1155 I = K + 1, L                                                          DO 1155 I = K + 1, L    
949  1155            GH(I) = GH1(I)                                          1155            GH(I) = GH1(I)                                        
950              NMAX = NMAX1                                                            NMAX = NMAX1                                              
951          ELSE                                                                    ELSE                                                          
952              K = NMAX1 * (NMAX1 + 2)                                                  K = NMAX1 * (NMAX1 + 2)                                    
953              L = NMAX2 * (NMAX2 + 2)                                                  L = NMAX2 * (NMAX2 + 2)                                    
954              DO 1166 I = K + 1, L                                                          DO 1166 I = K + 1, L  
955  1166            GH(I) = FACTOR * GH2(I)                                  1166            GH(I) = FACTOR * GH2(I)                                
956              NMAX = NMAX2                                                            NMAX = NMAX2                                              
957          ENDIF                                                                    ENDIF                                                          
958                                                                                                                                                                    
959          DO 1177 I = 1, K                                                              DO 1177 I = 1, K
960  1177        GH(I) = GH1(I) + FACTOR * GH2(I)                            1177        GH(I) = GH1(I) + FACTOR * GH2(I)                          
961                                                                                                                                                                    
962          RETURN                                                                  RETURN                                                        
# Line 910  C -------------------------------------- Line 964  C --------------------------------------
964  C  C
965  C  C
966          SUBROUTINE INITIZE(TP1,TL1,TP2,TL2,TP3,TL3)          SUBROUTINE INITIZE(TP1,TL1,TP2,TL2,TP3,TL3)
967            IMPLICIT REAL(8)(A-H)
968            IMPLICIT REAL(8)(O-Z)
969  C----------------------------------------------------------------  C----------------------------------------------------------------
970  C Initializes the parameters in COMMON/GENER/  C Initializes the parameters in COMMON/GENER/
971  C  C
# Line 925  C ERA, EREQU and ERPOL as recommended by Line 981  C ERA, EREQU and ERPOL as recommended by
981  C ASTRONOMICAL UNION .  C ASTRONOMICAL UNION .
982  C-----------------------------------------------------------------  C-----------------------------------------------------------------
983          INTEGER TL1,TL2,TL3          INTEGER TL1,TL2,TL3
984          CHARACTER *258 TP1,TP2,TP3          CHARACTER (len=258) TP1,TP2,TP3
985          INTEGER L1,L2,L3          INTEGER L1,L2,L3
986          CHARACTER *258 P1,P2,P3          CHARACTER *258 P1,P2,P3
987            REAL(8) AQUAD,BQUAD,ERA
988    
989          COMMON/PPATH/ L1,L2,L3,P1, P2, P3          COMMON/PPATH/ L1,L2,L3,P1, P2, P3
990          SAVE/PPATH/          SAVE/PPATH/
991    
992          COMMON/GENER/   UMR,ERA,AQUAD,BQUAD          COMMON/GENER/ERA,AQUAD,BQUAD,UMR
993          SAVE/GENER/          SAVE/GENER/
994    
         P1 = TP1  
         P2 = TP2  
         P3 = TP3  
995          L1 = TL1          L1 = TL1
996          L2 = TL2          L2 = TL2
997          L3 = TL3          L3 = TL3
998            P1 = TP1(1:L1)
999            P2 = TP2(1:L2)
1000            P3 = TP3(1:L3)
1001    
1002          ERA=6371.2          ERA=6371.2
1003          EREQU=6378.16          EREQU=6378.16
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.8
  ViewVC Help
Powered by ViewVC 1.1.23