/[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.11 by mocchiut, Fri Oct 31 14:02:53 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 191  C                          APPROXIMATION Line 191  C                          APPROXIMATION
191  C          B0           MAGNETIC FIELD STRENGTH IN GAUSS  C          B0           MAGNETIC FIELD STRENGTH IN GAUSS
192  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
193        DIMENSION         V(3),U(3,3),P(8,100),SP(3)        DIMENSION         V(3),U(3,3),P(8,100),SP(3)
194        COMMON            X(3),H(144)        COMMON            X(3),H(196)
195        COMMON/FIDB0/     SP        COMMON/FIDB0/     SP
196        SAVE /FIDB0/        SAVE /FIDB0/
197        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
198        SAVE /GENER/        SAVE /GENER/
199          REAL FLS
200  C  C
201  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3
202  C-- STEP IS STEP SIZE FOR FIELD LINE TRACING  C-- STEP IS STEP SIZE FOR FIELD LINE TRACING
# Line 203  C-- STEQ IS STEP SIZE FOR INTEGRATION Line 204  C-- STEQ IS STEP SIZE FOR INTEGRATION
204  C  C
205        DATA RMIN,RMAX    /0.05,1.01/        DATA RMIN,RMAX    /0.05,1.01/
206        DATA STEP,STEQ    /0.20,0.03/        DATA STEP,STEQ    /0.20,0.03/
207          BEQU=1.E10        BEQU=1.E10
208          FLS = FL
209  C*****ENTRY POINT  SHELLG  TO BE USED WITH GEODETIC CO-ORDINATES  C*****ENTRY POINT  SHELLG  TO BE USED WITH GEODETIC CO-ORDINATES
210        RLAT=GLAT*UMR        RLAT=GLAT*UMR
211        CT=SIN(RLAT)                                                      CT=SIN(RLAT)                                              
# Line 304  C*****INNER LOOP (FOR QUADRATURE)       Line 306  C*****INNER LOOP (FOR QUADRATURE)      
306        HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)                                    HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)                            
307        ZQ=Z*Z        ZQ=Z*Z
308        R=HLI+SQRT(HLI*HLI+ZQ)        R=HLI+SQRT(HLI*HLI+ZQ)
309          IF(R.NE.R)THEN
310             FL = FLS
311             RETURN
312          ENDIF
313        IF(R.LE.RMIN)GOTO30                                      IF(R.LE.RMIN)GOTO30                              
314        RQ=R*R        RQ=R*R
315        FF=SQRT(1.+3.*ZQ/RQ)                                      FF=SQRT(1.+3.*ZQ/RQ)                              
316        RADIK=B0-((D2*T+D1)*T+D0)*R*RQ*FF                        RADIK=B0-((D2*T+D1)*T+D0)*R*RQ*FF                
317        IF(R-RMAX)44,44,45                                        IF((R-RMAX).le.0.) goto 44                                
318          IF((R-RMAX).gt.0.) goto 45                                
319  45    ICODE=2                                            45    ICODE=2                                          
320        RADIK=RADIK-12.*(R-RMAX)**2                              RADIK=RADIK-12.*(R-RMAX)**2                      
321  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 348  C-- The minimal allowable value of FI wa
348  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.
349  C-- D. Bilitza, Nov 87.  C-- D. Bilitza, Nov 87.
350  C  C
351  11      FI=0.5*ABS(FI)/SQRT(B0)+1E-12                        11         FI=0.5*ABS(FI)/SQRT(B0)+1E-12                      
352  C  C
353  C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.    C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.  
354  C  C
355  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.
356  C  C
357        DIMOB0=DIMO/B0        DIMOB0=DIMO/B0
358          arg1=alog(FI)        arg1=alog(FI)
359          arg2=alog(DIMOB0)        arg2=alog(DIMOB0)
360  c       arg = FI*FI*FI/DIMOB0  c    arg = FI*FI*FI/DIMOB0
361  c       if(abs(arg).gt.88.0) arg=88.0  c    if(abs(arg).gt.88.0) arg=88.0
362        XX=3*arg1-arg2        XX=3*arg1-arg2
363        IF(XX.GT.23.0) GOTO 776          IF(XX.GT.23.0) GOTO 776  
364        IF(XX.GT.11.7) GOTO 775          IF(XX.GT.11.7) GOTO 775  
365        IF(XX.GT.+3.0) GOTO 774            IF(XX.GT.+3.0) GOTO 774    
366        IF(XX.GT.-3.0) GOTO 773          IF(XX.GT.-3.0) GOTO 773  
367        IF(XX.GT.-22.) GOTO 772          IF(XX.GT.-22.) GOTO 772  
368    771 GG=3.33338E-1*XX+3.0062102E-1                                  c  771 GG=3.33338E-1*XX+3.0062102E-1                                
369          GG=3.33338E-1*XX+3.0062102E-1                                
370        GOTO777                                                                  GOTO777                                                          
371    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+  
372       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 392  C* SUBROUTINE USED FOR FIELD LINE TRACIN Line 400  C* SUBROUTINE USED FOR FIELD LINE TRACIN
400  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *
401  C*******************************************************************  C*******************************************************************
402        DIMENSION         P(7),U(3,3)        DIMENSION         P(7),U(3,3)
403        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
404  C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES            C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES          
405        ZM=P(3)                                                                  ZM=P(3)                                                          
406        FLI=P(1)*P(1)+P(2)*P(2)+1E-15        FLI=P(1)*P(1)+P(2)*P(2)+1E-15
# Line 474  C                 TO THE LOCAL GEODETIC Line 482  C                 TO THE LOCAL GEODETIC
482  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
483  C                 AND DOWNWARD.    C                 AND DOWNWARD.  
484  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
485        DIMENSION         V(3),B(3),G(144)        DIMENSION         V(3),B(3),G(196)
486        CHARACTER*258      NAME        CHARACTER*258      NAME
487        INTEGER NMAX        INTEGER NMAX
488        REAL TIME        REAL TIME
489        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
490        COMMON/MODEL/     G,NMAX,TIME,NAME        COMMON/MODEL/     G,NMAX,TIME,NAME
491        SAVE/MODEL/        SAVE/MODEL/
492        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
# Line 528  C*****ENTRY POINT  FELDI  USED FOR L COM Line 536  C*****ENTRY POINT  FELDI  USED FOR L COM
536        Y=XI(2)*F                                                                Y=XI(2)*F                                                        
537        Z=XI(3)*(F+F)                                                            Z=XI(3)*(F+F)                                                    
538        I=I-2                                                                    I=I-2                                                            
539        IF(I-1)5,4,2                                                        c      print *,' I ',I
540          IF((I-1).lt.0) goto 5
541          IF((I-1).eq.0) goto 4
542          IF((I-1).gt.0) goto 2
543  2     DO 3 M=3,I,2                                                        2     DO 3 M=3,I,2                                                      
544        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))          
545       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 571  C*****ENTRY POINT  FELDI  USED FOR L COM
571        END                                                                      END                                                              
572  C  C
573  C  C
574          SUBROUTINE FELDCOF(YEAR,DIMO)        SUBROUTINE FELDCOF(YEAR,DIMO)
575  C------------------------------------------------------------------------  C------------------------------------------------------------------------
576  C  DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS  C     DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS
577  C  C    
578  C       INPUT:  YEAR    DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO  C     INPUT:  YEAR    DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO
579  C                       BE CALCULATED  C     BE CALCULATED
580  C       OUTPUT: DIMO    GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED  C     OUTPUT: DIMO    GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED
581  C                       TO EARTH'S RADIUS) AT THE TIME (YEAR)  C     TO EARTH'S RADIUS) AT THE TIME (YEAR)
582  C  D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771,  C     D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771,
583  C       (301)286-9536   NOV 1987.  C     (301)286-9536   NOV 1987.
584  C  ### updated to IGRF-2000 version -dkb- 5/31/2000  C     ### updated to IGRF-2000 version -dkb- 5/31/2000
585  C  ### updated to IGRF-2005 version -dkb- 3/24/2000  C     ### updated to IGRF-2005 version -dkb- 3/24/2000
586  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
587          CHARACTER*258    FIL1, FIL2                  CHARACTER*258    FIL1, FIL2          
588          CHARACTER*258    FILMOD        CHARACTER*258    FILMOD
589  C ### FILMOD, DTEMOD arrays +1  C     ### FILMOD, DTEMOD arrays +1
590  c        DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)  c     DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)
591          DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(3),DTEMOD(3)        DIMENSION       GH1(196),GH2(196),GHA(196),FILMOD(3),DTEMOD(3)
592          DOUBLE PRECISION X,F0,F        DOUBLE PRECISION X,F0,F
593          INTEGER L1,L2,L3        INTEGER L1,L2,L3
594          INTEGER NMAX        INTEGER NMAX
595          REAL TIME        REAL TIME
596          CHARACTER *258 P1,P2,P3        CHARACTER *258 P1,P2,P3
597          COMMON/PPATH/ L1,L2,L3,P1, P2, P3        COMMON/PPATH/ L1,L2,L3,P1, P2, P3
598          SAVE/PPATH/        SAVE/PPATH/
599          COMMON/MODEL/   GH1,NMAX,TIME,FIL1        COMMON/MODEL/   GH1,NMAX,TIME,FIL1
600          SAVE/MODEL/        SAVE/MODEL/
601          COMMON/GENER/   UMR,ERAD,AQUAD,BQUAD        COMMON/GENER/   UMR,ERAD,AQUAD,BQUAD
602          SAVE/GENER/        SAVE/GENER/
603  C ### updated to 2005  C     ### updated to 2005
604  C      CHARACTER COEFPATH*80, COEF1*80, COEF2*80, COEF3*80  C     CHARACTER COEFPATH*80, COEF1*80, COEF2*80, COEF3*80
605          
606  c      COEFPATH  = 'OrbitalInfo/src/'  c     COEFPATH  = 'OrbitalInfo/src/'
607  c      COEF1 = 'dgrf00.dat'  c     COEF1 = 'dgrf00.dat'
608  c      COEF2 = 'igrf05.dat'  c     COEF2 = 'igrf05.dat'
609  c      COEF3 = 'igrf05s.dat'  c     COEF3 = 'igrf05s.dat'
610  c      COEF1 = COEFPATH(1:16)//COEF1  c     COEF1 = COEFPATH(1:16)//COEF1
611  c      COEF2 = COEFPATH(1:16)//COEF2  c     COEF2 = COEFPATH(1:16)//COEF2
612  c      COEF3 = COEFPATH(1:16)//COEF3  c     COEF3 = COEFPATH(1:16)//COEF3
613  c      FILMOD(1) = COEF1  c     FILMOD(1) = COEF1
614  c      FILMOD(2) = COEF2  c     FILMOD(2) = COEF2
615  c      FILMOD(3) = COEF3  c     FILMOD(3) = COEF3
616    c      print *, "qui"
617        FILMOD(1) = P1(1:L1)        FILMOD(1) = P1(1:L1)
618        FILMOD(2) = P2(1:L2)        FILMOD(2) = P2(1:L2)
619        FILMOD(3) = P3(1:L3)        FILMOD(3) = P3(1:L3)
620  c      FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'  c      print *, "qua"
621  c      FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'
622  c      FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'
623    c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'
624  c      WRITE(*,*) FILMOD(1)  c      WRITE(*,*) FILMOD(1)
625  c      WRITE(*,*) FILMOD(2)  c      WRITE(*,*) FILMOD(2)
626  c      WRITE(*,*) FILMOD(3)  c      WRITE(*,*) FILMOD(3)
627  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/
628        DATA   DTEMOD / 2000., 2005., 2010./        DATA   DTEMOD / 2005., 2010., 2015./
629  c  c    
630  c        DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',              c     DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',            
631  c     1                  'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',        c     1                  'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',      
632  c     2                  'dgrf70.dat', 'dgrf75.dat', 'dgrf80.dat',        c     2                  'dgrf70.dat', 'dgrf75.dat', 'dgrf80.dat',      
633  c     3                  'dgrf85.dat', 'dgrf90.dat', 'dgrf95.dat',  c     3                  'dgrf85.dat', 'dgrf90.dat', 'dgrf95.dat',
634  c     4                  'dgrf00.dat','igrf05.dat','igrf05s.dat'/  c     4                  'dgrf00.dat','igrf05.dat','igrf05s.dat'/
635  c        DATA   DTEMOD / 1945., 1950., 1955., 1960., 1965., 1970.,  c     DATA   DTEMOD / 1945., 1950., 1955., 1960., 1965., 1970.,
636  c     1       1975., 1980., 1985., 1990., 1995., 2000.,2005.,2010./        c     1       1975., 1980., 1985., 1990., 1995., 2000.,2005.,2010./      
637  C  C    
638  C ### numye = numye + 1 ; is number of years represented by IGRF  C     ### numye = numye + 1 ; is number of years represented by IGRF
639  C  C    
640  c        NUMYE=13  c     NUMYE=13
641          NUMYE=2        NUMYE=2
642    c      print *, "quo"
643  C        
644  C  IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION  C    
645  C  IU  IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS  C     IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION
646  C  C     IU  IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS
647          IU = 10  C    
648          IS = 0        IU = 10
649  C-- DETERMINE IGRF-YEARS FOR INPUT-YEAR        IS = 0
650          TIME = YEAR  C--   DETERMINE IGRF-YEARS FOR INPUT-YEAR
651          IYEA = INT(YEAR/5.)*5        TIME = YEAR
652  c        L = (IYEA - 1945)/5 + 1        IYEA = INT(YEAR/5.)*5
653          L = (IYEA - 2000)/5 + 1  c     L = (IYEA - 1945)/5 + 1
654          IF(L.LT.1) L=1        L = (IYEA - 2000)/5 + 1
655          IF(L.GT.NUMYE) L=NUMYE                IF(L.LT.1) L=1
656          DTE1 = DTEMOD(L)          IF(L.GT.NUMYE) L=NUMYE        
657          FIL1 = FILMOD(L)          DTE1 = DTEMOD(L)  
658          DTE2 = DTEMOD(L+1)        FIL1 = FILMOD(L)  
659          FIL2 = FILMOD(L+1)        DTE2 = DTEMOD(L+1)
660  C-- GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS        FIL2 = FILMOD(L+1)
661          CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)    c      WRITE(*,*) FIL1
662              IF (IER .NE. 0) STOP                            c      WRITE(*,*) FIL2
663          CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)    c      print *, "que"
664              IF (IER .NE. 0) STOP                      C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
665  C-- DETERMINE IGRF COEFFICIENTS FOR YEAR        CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
666          IF (L .LE. NUMYE-1) THEN                                IF (IER .NE. 0) STOP                          
667            CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,  c      print *, "quessss"
668       1          NMAX2, GH2, NMAX, GHA)                                CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
669          ELSE                      IF (IER .NE. 0) STOP                    
670            CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,      c      print *, "quj"
671       1          GH2, NMAX, GHA)                                      C--   DETERMINE IGRF COEFFICIENTS FOR YEAR
672          ENDIF        IF (L .LE. NUMYE-1) THEN                        
673  C-- DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G           CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,
674          F0=0.D0       1        NMAX2, GH2, NMAX, GHA)                        
675          DO 1234 J=1,3        ELSE              
676             F = GHA(J) * 1.D-5           CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,    
677             F0 = F0 + F * F       1        GH2, NMAX, GHA)                                    
678  1234    CONTINUE        ENDIF
679          DIMO = DSQRT(F0)  c      print *, "quw"
680    C--   DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G
681          GH1(1) =  0.0        F0=0.D0
682          I=2                  DO 1234 J=1,3
683          F0=1.D-5                           F = GHA(J) * 1.D-5
684          IF(IS.EQ.0) F0=-F0           F0 = F0 + F * F
685          SQRT2=SQRT(2.)         1234 CONTINUE
686          DIMO = REAL(DSQRT(F0))
687          
688          GH1(1) =  0.0
689          I=2          
690          F0=1.D-5                
691          IF(IS.EQ.0) F0=-F0
692          SQRT2=SQRT(2.)      
693    
694    c      print *, "quq"
695          
696        DO 9 N=1,NMAX                  DO 9 N=1,NMAX          
697          X = N           X = N
698          F0 = F0 * X * X / (4.D0 * X - 2.D0)                         F0 = F0 * X * X / (4.D0 * X - 2.D0)              
699          IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X
700          F = F0 * 0.5D0                                               F = F0 * 0.5D0                                    
701          IF(IS.EQ.0) F = F * SQRT2           IF(IS.EQ.0) F = F * SQRT2
702          GH1(I) = GHA(I-1) * F0           GH1(I) = GHA(I-1) * REAL(F0)
703          I = I+1                                                   I = I+1                                        
704        DO 9 M=1,N                                               DO 9 M=1,N                                    
705          F = F * (X + M) / (X - M + 1.D0)                              F = F * (X + M) / (X - M + 1.D0)                
706          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))
707          GH1(I) = GHA(I-1) * F              GH1(I) = GHA(I-1) * REAL(F)
708          GH1(I+1) = GHA(I) * F              GH1(I+1) = GHA(I) * REAL(F)
709          I=I+2              I=I+2
710  9     CONTINUE                                             9       CONTINUE                                          
711          RETURN        RETURN
712          END        END
713  C  C    
714  C  C    
715          SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)                  SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)          
716                                                                                                                                                                    
717  C ===============================================================                C ===============================================================              
718  C                                                                                C                                                                              
719  C       Version 1.01                                                  C     Version 1.01                                                
720  C                                                                                C    
721  C       Reads spherical harmonic coefficients from the specified      C     Reads spherical harmonic coefficients from the specified    
722  C       file into an array.                                            C     file into an array.                                          
723  C                                                                                C    
724  C       Input:                                                        C     Input:                                                      
725  C           IU    - Logical unit number                                C     IU    - Logical unit number                              
726  C           FSPEC - File specification                                C     FSPEC - File specification                              
727  C                                                                                C    
728  C       Output:                                                        C     Output:                                                      
729  C           NMAX  - Maximum degree and order of model                  C     NMAX  - Maximum degree and order of model                
730  C           ERAD  - Earth's radius associated with the spherical      C     ERAD  - Earth's radius associated with the spherical    
731  C                   harmonic coefficients, in the same units as        C     harmonic coefficients, in the same units as      
732  C                   elevation                                          C     elevation                                        
733  C           GH    - Schmidt quasi-normal internal spherical            C     GH    - Schmidt quasi-normal internal spherical          
734  C                   harmonic coefficients                              C     harmonic coefficients                            
735  C           IER   - Error number: =  0, no error                      C     IER   - Error number: =  0, no error                    
736  C                                 = -2, records out of order          C     = -2, records out of order        
737  C                                 = FORTRAN run-time error number      C     = FORTRAN run-time error number    
738  C                                                                      C    
739  C       A. Zunde                                                      C     A. Zunde                                                    
740  C       USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225      C     USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225    
741  C                                                                                C    
742  C ===============================================================                C     ===============================================================              
743                                                                                          
744          CHARACTER  FSPEC*(*), FOUT*258                                            CHARACTER  FSPEC*(*), FOUT*258
745          DIMENSION       GH(*)                                                DIMENSION       GH(*)                                        
746  C ---------------------------------------------------------------                C     ---------------------------------------------------------------              
747  C       Open coefficient file. Read past first header record.          C     Open coefficient file. Read past first header record.        
748  C       Read degree and order of model and Earth's radius.            C     Read degree and order of model and Earth's radius.          
749  C ---------------------------------------------------------------                C     ---------------------------------------------------------------              
750        WRITE(FOUT,667) FSPEC        WRITE(FOUT,667) FSPEC
751  c 667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)  c     667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)
752  667   FORMAT(A258)   667  FORMAT(A258)
753    c      print *," gui"
754        OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)            OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)    
755    c      print *," gua"
756        READ (IU, *, IOSTAT=IER, ERR=999)        READ (IU, *, IOSTAT=IER, ERR=999)
757    c      print *," gue"
758        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD
759    c      print *," guo"
760  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
761  C       Read the coefficient file, arranged as follows:                C       Read the coefficient file, arranged as follows:              
762  C                                                                                C                                                                              
# Line 752  C       coefficient.                     Line 778  C       coefficient.                    
778  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
779                                                                                                                                                                    
780          I = 0                                                                  I = 0                                                        
781          DO 2211 NN = 1, NMAX                                                        DO 2211 NN = 1, NMAX              
782              DO 2233 MM = 0, NN                                                          DO 2233 MM = 0, NN            
783                  READ (IU, *, IOSTAT=IER, ERR=999) N, M, G, H                          READ (IU, *, IOSTAT=IER, ERR=999) N, M, G, H        
784                  IF (NN .NE. N .OR. MM .NE. M) THEN                                    IF (NN .NE. N .OR. MM .NE. M) THEN                  
785                      IER = -2                                                              IER = -2                                        
# Line 767  C -------------------------------------- Line 793  C --------------------------------------
793                  ENDIF                                                                  ENDIF                                                
794  2233        CONTINUE                                                      2233        CONTINUE                                                    
795  2211    CONTINUE                                                          2211    CONTINUE                                                        
796    c      print *," guj"
797                                                                                                                                                                    
798  999     CLOSE (IU)                                                    999     CLOSE (IU)                                                  
799    c      print *," guw IER",IER
800          if ( IER .eq. -1 ) IER = 0 ! gfortran 4.1.2 bug workaround... hoping not to create problems with other versions
801                                                                                                                                                                    
802          RETURN                                                                RETURN                                                      
803          END                                                                    END                                                          
# Line 824  C -------------------------------------- Line 852  C --------------------------------------
852          ELSE IF (NMAX1 .GT. NMAX2) THEN                                        ELSE IF (NMAX1 .GT. NMAX2) THEN                              
853              K = NMAX2 * (NMAX2 + 2)                                                K = NMAX2 * (NMAX2 + 2)                                  
854              L = NMAX1 * (NMAX1 + 2)                                                L = NMAX1 * (NMAX1 + 2)                                  
855              DO 1122 I = K + 1, L                                                        DO 1122 I = K + 1, L    
856  1122            GH(I) = GH1(I) + FACTOR * (-GH1(I))                    1122            GH(I) = GH1(I) + FACTOR * (-GH1(I))                  
857              NMAX = NMAX1                                                          NMAX = NMAX1                                            
858          ELSE                                                                  ELSE                                                        
859              K = NMAX1 * (NMAX1 + 2)                                                K = NMAX1 * (NMAX1 + 2)                                  
860              L = NMAX2 * (NMAX2 + 2)                                                L = NMAX2 * (NMAX2 + 2)                                  
861              DO 1133 I = K + 1, L                                                        DO 1133 I = K + 1, L
862  1133            GH(I) = FACTOR * GH2(I)                                1133            GH(I) = FACTOR * GH2(I)                              
863              NMAX = NMAX2                                                          NMAX = NMAX2                                            
864          ENDIF                                                                  ENDIF                                                        
865                                                                                                                                                                    
866          DO 1144 I = 1, K                                                            DO 1144 I = 1, K  
867  1144        GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I))                1144        GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I))              
868                                                                                                                                                                    
869          RETURN                                                                RETURN                                                      
# Line 891  C -------------------------------------- Line 919  C --------------------------------------
919          ELSE IF (NMAX1 .GT. NMAX2) THEN                                          ELSE IF (NMAX1 .GT. NMAX2) THEN                                
920              K = NMAX2 * (NMAX2 + 2)                                                  K = NMAX2 * (NMAX2 + 2)                                    
921              L = NMAX1 * (NMAX1 + 2)                                                  L = NMAX1 * (NMAX1 + 2)                                    
922              DO 1155 I = K + 1, L                                                          DO 1155 I = K + 1, L    
923  1155            GH(I) = GH1(I)                                          1155            GH(I) = GH1(I)                                        
924              NMAX = NMAX1                                                            NMAX = NMAX1                                              
925          ELSE                                                                    ELSE                                                          
926              K = NMAX1 * (NMAX1 + 2)                                                  K = NMAX1 * (NMAX1 + 2)                                    
927              L = NMAX2 * (NMAX2 + 2)                                                  L = NMAX2 * (NMAX2 + 2)                                    
928              DO 1166 I = K + 1, L                                                          DO 1166 I = K + 1, L  
929  1166            GH(I) = FACTOR * GH2(I)                                  1166            GH(I) = FACTOR * GH2(I)                                
930              NMAX = NMAX2                                                            NMAX = NMAX2                                              
931          ENDIF                                                                    ENDIF                                                          
932                                                                                                                                                                    
933          DO 1177 I = 1, K                                                              DO 1177 I = 1, K
934  1177        GH(I) = GH1(I) + FACTOR * GH2(I)                            1177        GH(I) = GH1(I) + FACTOR * GH2(I)                          
935                                                                                                                                                                    
936          RETURN                                                                  RETURN                                                        
# Line 925  C ERA, EREQU and ERPOL as recommended by Line 953  C ERA, EREQU and ERPOL as recommended by
953  C ASTRONOMICAL UNION .  C ASTRONOMICAL UNION .
954  C-----------------------------------------------------------------  C-----------------------------------------------------------------
955          INTEGER TL1,TL2,TL3          INTEGER TL1,TL2,TL3
956          CHARACTER *258 TP1,TP2,TP3          CHARACTER (len=*) :: TP1,TP2,TP3
957          INTEGER L1,L2,L3          INTEGER L1,L2,L3
958          CHARACTER *258 P1,P2,P3          CHARACTER *258 P1,P2,P3
959          COMMON/PPATH/ L1,L2,L3,P1, P2, P3          COMMON/PPATH/ L1,L2,L3,P1, P2, P3
960          SAVE/PPATH/          SAVE/PPATH/
961    
962          COMMON/GENER/   UMR,ERA,AQUAD,BQUAD          COMMON/GENER/UMR,ERA,AQUAD,BQUAD
963          SAVE/GENER/          SAVE/GENER/
964    
         P1 = TP1  
         P2 = TP2  
         P3 = TP3  
965          L1 = TL1          L1 = TL1
966          L2 = TL2          L2 = TL2
967          L3 = TL3          L3 = TL3
968            P1 = TP1(1:L1)
969            P2 = TP2(1:L2)
970            P3 = TP3(1:L3)
971    
972          ERA=6371.2          ERA=6371.2
973          EREQU=6378.16          EREQU=6378.16
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.11
  ViewVC Help
Powered by ViewVC 1.1.23