/[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.9 by mocchiut, Fri Jan 17 09:54:31 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
# Line 308  C*****INNER LOOP (FOR QUADRATURE)       Line 308  C*****INNER LOOP (FOR QUADRATURE)      
308        RQ=R*R        RQ=R*R
309        FF=SQRT(1.+3.*ZQ/RQ)                                      FF=SQRT(1.+3.*ZQ/RQ)                              
310        RADIK=B0-((D2*T+D1)*T+D0)*R*RQ*FF                        RADIK=B0-((D2*T+D1)*T+D0)*R*RQ*FF                
311        IF(R-RMAX)44,44,45                                        IF((R-RMAX).le.0.) goto 44                                
312          IF((R-RMAX).gt.0.) goto 45                                
313  45    ICODE=2                                            45    ICODE=2                                          
314        RADIK=RADIK-12.*(R-RMAX)**2                              RADIK=RADIK-12.*(R-RMAX)**2                      
315  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 342  C-- The minimal allowable value of FI wa
342  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.
343  C-- D. Bilitza, Nov 87.  C-- D. Bilitza, Nov 87.
344  C  C
345  11      FI=0.5*ABS(FI)/SQRT(B0)+1E-12                        11         FI=0.5*ABS(FI)/SQRT(B0)+1E-12                      
346  C  C
347  C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.    C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.  
348  C  C
349  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.
350  C  C
351        DIMOB0=DIMO/B0        DIMOB0=DIMO/B0
352          arg1=alog(FI)        arg1=alog(FI)
353          arg2=alog(DIMOB0)        arg2=alog(DIMOB0)
354  c       arg = FI*FI*FI/DIMOB0  c    arg = FI*FI*FI/DIMOB0
355  c       if(abs(arg).gt.88.0) arg=88.0  c    if(abs(arg).gt.88.0) arg=88.0
356        XX=3*arg1-arg2        XX=3*arg1-arg2
357        IF(XX.GT.23.0) GOTO 776          IF(XX.GT.23.0) GOTO 776  
358        IF(XX.GT.11.7) GOTO 775          IF(XX.GT.11.7) GOTO 775  
359        IF(XX.GT.+3.0) GOTO 774            IF(XX.GT.+3.0) GOTO 774    
360        IF(XX.GT.-3.0) GOTO 773          IF(XX.GT.-3.0) GOTO 773  
361        IF(XX.GT.-22.) GOTO 772          IF(XX.GT.-22.) GOTO 772  
362    771 GG=3.33338E-1*XX+3.0062102E-1                                  c  771 GG=3.33338E-1*XX+3.0062102E-1                                
363          GG=3.33338E-1*XX+3.0062102E-1                                
364        GOTO777                                                                  GOTO777                                                          
365    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+  
366       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 394  C* SUBROUTINE USED FOR FIELD LINE TRACIN
394  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *
395  C*******************************************************************  C*******************************************************************
396        DIMENSION         P(7),U(3,3)        DIMENSION         P(7),U(3,3)
397        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
398  C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES            C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES          
399        ZM=P(3)                                                                  ZM=P(3)                                                          
400        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 476  C                 TO THE LOCAL GEODETIC
476  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
477  C                 AND DOWNWARD.    C                 AND DOWNWARD.  
478  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
479        DIMENSION         V(3),B(3),G(144)        DIMENSION         V(3),B(3),G(196)
480        CHARACTER*258      NAME        CHARACTER*258      NAME
481        INTEGER NMAX        INTEGER NMAX
482        REAL TIME        REAL TIME
483        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
484        COMMON/MODEL/     G,NMAX,TIME,NAME        COMMON/MODEL/     G,NMAX,TIME,NAME
485        SAVE/MODEL/        SAVE/MODEL/
486        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
# Line 528  C*****ENTRY POINT  FELDI  USED FOR L COM Line 530  C*****ENTRY POINT  FELDI  USED FOR L COM
530        Y=XI(2)*F                                                                Y=XI(2)*F                                                        
531        Z=XI(3)*(F+F)                                                            Z=XI(3)*(F+F)                                                    
532        I=I-2                                                                    I=I-2                                                            
533        IF(I-1)5,4,2                                                        c      print *,' I ',I
534          IF((I-1).lt.0) goto 5
535          IF((I-1).eq.0) goto 4
536          IF((I-1).gt.0) goto 2
537  2     DO 3 M=3,I,2                                                        2     DO 3 M=3,I,2                                                      
538        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))          
539       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 565  C*****ENTRY POINT  FELDI  USED FOR L COM
565        END                                                                      END                                                              
566  C  C
567  C  C
568          SUBROUTINE FELDCOF(YEAR,DIMO)        SUBROUTINE FELDCOF(YEAR,DIMO)
569  C------------------------------------------------------------------------  C------------------------------------------------------------------------
570  C  DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS  C     DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS
571  C  C    
572  C       INPUT:  YEAR    DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO  C     INPUT:  YEAR    DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO
573  C                       BE CALCULATED  C     BE CALCULATED
574  C       OUTPUT: DIMO    GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED  C     OUTPUT: DIMO    GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED
575  C                       TO EARTH'S RADIUS) AT THE TIME (YEAR)  C     TO EARTH'S RADIUS) AT THE TIME (YEAR)
576  C  D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771,  C     D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771,
577  C       (301)286-9536   NOV 1987.  C     (301)286-9536   NOV 1987.
578  C  ### updated to IGRF-2000 version -dkb- 5/31/2000  C     ### updated to IGRF-2000 version -dkb- 5/31/2000
579  C  ### updated to IGRF-2005 version -dkb- 3/24/2000  C     ### updated to IGRF-2005 version -dkb- 3/24/2000
580  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
581          CHARACTER*258    FIL1, FIL2                  CHARACTER*258    FIL1, FIL2          
582          CHARACTER*258    FILMOD        CHARACTER*258    FILMOD
583  C ### FILMOD, DTEMOD arrays +1  C     ### FILMOD, DTEMOD arrays +1
584  c        DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)  c     DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)
585          DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(3),DTEMOD(3)        DIMENSION       GH1(196),GH2(196),GHA(196),FILMOD(3),DTEMOD(3)
586          DOUBLE PRECISION X,F0,F        DOUBLE PRECISION X,F0,F
587          INTEGER L1,L2,L3        INTEGER L1,L2,L3
588          INTEGER NMAX        INTEGER NMAX
589          REAL TIME        REAL TIME
590          CHARACTER *258 P1,P2,P3        CHARACTER *258 P1,P2,P3
591          COMMON/PPATH/ L1,L2,L3,P1, P2, P3        COMMON/PPATH/ L1,L2,L3,P1, P2, P3
592          SAVE/PPATH/        SAVE/PPATH/
593          COMMON/MODEL/   GH1,NMAX,TIME,FIL1        COMMON/MODEL/   GH1,NMAX,TIME,FIL1
594          SAVE/MODEL/        SAVE/MODEL/
595          COMMON/GENER/   UMR,ERAD,AQUAD,BQUAD        COMMON/GENER/   UMR,ERAD,AQUAD,BQUAD
596          SAVE/GENER/        SAVE/GENER/
597  C ### updated to 2005  C     ### updated to 2005
598  C      CHARACTER COEFPATH*80, COEF1*80, COEF2*80, COEF3*80  C     CHARACTER COEFPATH*80, COEF1*80, COEF2*80, COEF3*80
599          
600  c      COEFPATH  = 'OrbitalInfo/src/'  c     COEFPATH  = 'OrbitalInfo/src/'
601  c      COEF1 = 'dgrf00.dat'  c     COEF1 = 'dgrf00.dat'
602  c      COEF2 = 'igrf05.dat'  c     COEF2 = 'igrf05.dat'
603  c      COEF3 = 'igrf05s.dat'  c     COEF3 = 'igrf05s.dat'
604  c      COEF1 = COEFPATH(1:16)//COEF1  c     COEF1 = COEFPATH(1:16)//COEF1
605  c      COEF2 = COEFPATH(1:16)//COEF2  c     COEF2 = COEFPATH(1:16)//COEF2
606  c      COEF3 = COEFPATH(1:16)//COEF3  c     COEF3 = COEFPATH(1:16)//COEF3
607  c      FILMOD(1) = COEF1  c     FILMOD(1) = COEF1
608  c      FILMOD(2) = COEF2  c     FILMOD(2) = COEF2
609  c      FILMOD(3) = COEF3  c     FILMOD(3) = COEF3
610    c      print *, "qui"
611        FILMOD(1) = P1(1:L1)        FILMOD(1) = P1(1:L1)
612        FILMOD(2) = P2(1:L2)        FILMOD(2) = P2(1:L2)
613        FILMOD(3) = P3(1:L3)        FILMOD(3) = P3(1:L3)
614  c      FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'  c      print *, "qua"
615  c      FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'  c     FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'
616  c      FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'  c     FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'
617    c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'
618  c      WRITE(*,*) FILMOD(1)  c      WRITE(*,*) FILMOD(1)
619  c      WRITE(*,*) FILMOD(2)  c      WRITE(*,*) FILMOD(2)
620  c      WRITE(*,*) FILMOD(3)  c      WRITE(*,*) FILMOD(3)
621  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/
622        DATA   DTEMOD / 2000., 2005., 2010./        DATA   DTEMOD / 2005., 2010., 2015./
623  c  c    
624  c        DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',              c     DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',            
625  c     1                  'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',        c     1                  'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',      
626  c     2                  'dgrf70.dat', 'dgrf75.dat', 'dgrf80.dat',        c     2                  'dgrf70.dat', 'dgrf75.dat', 'dgrf80.dat',      
627  c     3                  'dgrf85.dat', 'dgrf90.dat', 'dgrf95.dat',  c     3                  'dgrf85.dat', 'dgrf90.dat', 'dgrf95.dat',
628  c     4                  'dgrf00.dat','igrf05.dat','igrf05s.dat'/  c     4                  'dgrf00.dat','igrf05.dat','igrf05s.dat'/
629  c        DATA   DTEMOD / 1945., 1950., 1955., 1960., 1965., 1970.,  c     DATA   DTEMOD / 1945., 1950., 1955., 1960., 1965., 1970.,
630  c     1       1975., 1980., 1985., 1990., 1995., 2000.,2005.,2010./        c     1       1975., 1980., 1985., 1990., 1995., 2000.,2005.,2010./      
631  C  C    
632  C ### numye = numye + 1 ; is number of years represented by IGRF  C     ### numye = numye + 1 ; is number of years represented by IGRF
633  C  C    
634  c        NUMYE=13  c     NUMYE=13
635          NUMYE=2        NUMYE=2
636    c      print *, "quo"
637  C        
638  C  IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION  C    
639  C  IU  IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS  C     IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION
640  C  C     IU  IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS
641          IU = 10  C    
642          IS = 0        IU = 10
643  C-- DETERMINE IGRF-YEARS FOR INPUT-YEAR        IS = 0
644          TIME = YEAR  C--   DETERMINE IGRF-YEARS FOR INPUT-YEAR
645          IYEA = INT(YEAR/5.)*5        TIME = YEAR
646  c        L = (IYEA - 1945)/5 + 1        IYEA = INT(YEAR/5.)*5
647          L = (IYEA - 2000)/5 + 1  c     L = (IYEA - 1945)/5 + 1
648          IF(L.LT.1) L=1        L = (IYEA - 2000)/5 + 1
649          IF(L.GT.NUMYE) L=NUMYE                IF(L.LT.1) L=1
650          DTE1 = DTEMOD(L)          IF(L.GT.NUMYE) L=NUMYE        
651          FIL1 = FILMOD(L)          DTE1 = DTEMOD(L)  
652          DTE2 = DTEMOD(L+1)        FIL1 = FILMOD(L)  
653          FIL2 = FILMOD(L+1)        DTE2 = DTEMOD(L+1)
654  C-- GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS        FIL2 = FILMOD(L+1)
655          CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)    c      WRITE(*,*) FIL1
656              IF (IER .NE. 0) STOP                            c      WRITE(*,*) FIL2
657          CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)    c      print *, "que"
658              IF (IER .NE. 0) STOP                      C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
659  C-- DETERMINE IGRF COEFFICIENTS FOR YEAR        CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
660          IF (L .LE. NUMYE-1) THEN                                IF (IER .NE. 0) STOP                          
661            CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,  c      print *, "quessss"
662       1          NMAX2, GH2, NMAX, GHA)                                CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
663          ELSE                      IF (IER .NE. 0) STOP                    
664            CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,      c      print *, "quj"
665       1          GH2, NMAX, GHA)                                      C--   DETERMINE IGRF COEFFICIENTS FOR YEAR
666          ENDIF        IF (L .LE. NUMYE-1) THEN                        
667  C-- DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G           CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,
668          F0=0.D0       1        NMAX2, GH2, NMAX, GHA)                        
669          DO 1234 J=1,3        ELSE              
670             F = GHA(J) * 1.D-5           CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,    
671             F0 = F0 + F * F       1        GH2, NMAX, GHA)                                    
672  1234    CONTINUE        ENDIF
673          DIMO = DSQRT(F0)  c      print *, "quw"
674    C--   DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G
675          GH1(1) =  0.0        F0=0.D0
676          I=2                  DO 1234 J=1,3
677          F0=1.D-5                           F = GHA(J) * 1.D-5
678          IF(IS.EQ.0) F0=-F0           F0 = F0 + F * F
679          SQRT2=SQRT(2.)         1234 CONTINUE
680          DIMO = REAL(DSQRT(F0))
681          
682          GH1(1) =  0.0
683          I=2          
684          F0=1.D-5                
685          IF(IS.EQ.0) F0=-F0
686          SQRT2=SQRT(2.)      
687    
688    c      print *, "quq"
689          
690        DO 9 N=1,NMAX                  DO 9 N=1,NMAX          
691          X = N           X = N
692          F0 = F0 * X * X / (4.D0 * X - 2.D0)                         F0 = F0 * X * X / (4.D0 * X - 2.D0)              
693          IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X
694          F = F0 * 0.5D0                                               F = F0 * 0.5D0                                    
695          IF(IS.EQ.0) F = F * SQRT2           IF(IS.EQ.0) F = F * SQRT2
696          GH1(I) = GHA(I-1) * F0           GH1(I) = GHA(I-1) * REAL(F0)
697          I = I+1                                                   I = I+1                                        
698        DO 9 M=1,N                                               DO 9 M=1,N                                    
699          F = F * (X + M) / (X - M + 1.D0)                              F = F * (X + M) / (X - M + 1.D0)                
700          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))
701          GH1(I) = GHA(I-1) * F              GH1(I) = GHA(I-1) * REAL(F)
702          GH1(I+1) = GHA(I) * F              GH1(I+1) = GHA(I) * REAL(F)
703          I=I+2              I=I+2
704  9     CONTINUE                                             9       CONTINUE                                          
705          RETURN        RETURN
706          END        END
707  C  C    
708  C  C    
709          SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)                  SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)          
710                                                                                                                                                                    
711  C ===============================================================                C ===============================================================              
712  C                                                                                C                                                                              
713  C       Version 1.01                                                  C     Version 1.01                                                
714  C                                                                                C    
715  C       Reads spherical harmonic coefficients from the specified      C     Reads spherical harmonic coefficients from the specified    
716  C       file into an array.                                            C     file into an array.                                          
717  C                                                                                C    
718  C       Input:                                                        C     Input:                                                      
719  C           IU    - Logical unit number                                C     IU    - Logical unit number                              
720  C           FSPEC - File specification                                C     FSPEC - File specification                              
721  C                                                                                C    
722  C       Output:                                                        C     Output:                                                      
723  C           NMAX  - Maximum degree and order of model                  C     NMAX  - Maximum degree and order of model                
724  C           ERAD  - Earth's radius associated with the spherical      C     ERAD  - Earth's radius associated with the spherical    
725  C                   harmonic coefficients, in the same units as        C     harmonic coefficients, in the same units as      
726  C                   elevation                                          C     elevation                                        
727  C           GH    - Schmidt quasi-normal internal spherical            C     GH    - Schmidt quasi-normal internal spherical          
728  C                   harmonic coefficients                              C     harmonic coefficients                            
729  C           IER   - Error number: =  0, no error                      C     IER   - Error number: =  0, no error                    
730  C                                 = -2, records out of order          C     = -2, records out of order        
731  C                                 = FORTRAN run-time error number      C     = FORTRAN run-time error number    
732  C                                                                      C    
733  C       A. Zunde                                                      C     A. Zunde                                                    
734  C       USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225      C     USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225    
735  C                                                                                C    
736  C ===============================================================                C     ===============================================================              
737                                                                                          
738          CHARACTER  FSPEC*(*), FOUT*258                                            CHARACTER  FSPEC*(*), FOUT*258
739          DIMENSION       GH(*)                                                DIMENSION       GH(*)                                        
740  C ---------------------------------------------------------------                C     ---------------------------------------------------------------              
741  C       Open coefficient file. Read past first header record.          C     Open coefficient file. Read past first header record.        
742  C       Read degree and order of model and Earth's radius.            C     Read degree and order of model and Earth's radius.          
743  C ---------------------------------------------------------------                C     ---------------------------------------------------------------              
744        WRITE(FOUT,667) FSPEC        WRITE(FOUT,667) FSPEC
745  c 667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)  c     667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)
746  667   FORMAT(A258)   667  FORMAT(A258)
747    c      print *," gui"
748        OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)            OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)    
749    c      print *," gua"
750        READ (IU, *, IOSTAT=IER, ERR=999)        READ (IU, *, IOSTAT=IER, ERR=999)
751    c      print *," gue"
752        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD
753    c      print *," guo"
754  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
755  C       Read the coefficient file, arranged as follows:                C       Read the coefficient file, arranged as follows:              
756  C                                                                                C                                                                              
# Line 752  C       coefficient.                     Line 772  C       coefficient.                    
772  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
773                                                                                                                                                                    
774          I = 0                                                                  I = 0                                                        
775          DO 2211 NN = 1, NMAX                                                        DO 2211 NN = 1, NMAX              
776              DO 2233 MM = 0, NN                                                          DO 2233 MM = 0, NN            
777                  READ (IU, *, IOSTAT=IER, ERR=999) N, M, G, H                          READ (IU, *, IOSTAT=IER, ERR=999) N, M, G, H        
778                  IF (NN .NE. N .OR. MM .NE. M) THEN                                    IF (NN .NE. N .OR. MM .NE. M) THEN                  
779                      IER = -2                                                              IER = -2                                        
# Line 767  C -------------------------------------- Line 787  C --------------------------------------
787                  ENDIF                                                                  ENDIF                                                
788  2233        CONTINUE                                                      2233        CONTINUE                                                    
789  2211    CONTINUE                                                          2211    CONTINUE                                                        
790    c      print *," guj"
791                                                                                                                                                                    
792  999     CLOSE (IU)                                                    999     CLOSE (IU)                                                  
793    c      print *," guw IER",IER
794          if ( IER .eq. -1 ) IER = 0 ! gfortran 4.1.2 bug workaround... hoping not to create problems with other versions
795                                                                                                                                                                    
796          RETURN                                                                RETURN                                                      
797          END                                                                    END                                                          
# Line 824  C -------------------------------------- Line 846  C --------------------------------------
846          ELSE IF (NMAX1 .GT. NMAX2) THEN                                        ELSE IF (NMAX1 .GT. NMAX2) THEN                              
847              K = NMAX2 * (NMAX2 + 2)                                                K = NMAX2 * (NMAX2 + 2)                                  
848              L = NMAX1 * (NMAX1 + 2)                                                L = NMAX1 * (NMAX1 + 2)                                  
849              DO 1122 I = K + 1, L                                                        DO 1122 I = K + 1, L    
850  1122            GH(I) = GH1(I) + FACTOR * (-GH1(I))                    1122            GH(I) = GH1(I) + FACTOR * (-GH1(I))                  
851              NMAX = NMAX1                                                          NMAX = NMAX1                                            
852          ELSE                                                                  ELSE                                                        
853              K = NMAX1 * (NMAX1 + 2)                                                K = NMAX1 * (NMAX1 + 2)                                  
854              L = NMAX2 * (NMAX2 + 2)                                                L = NMAX2 * (NMAX2 + 2)                                  
855              DO 1133 I = K + 1, L                                                        DO 1133 I = K + 1, L
856  1133            GH(I) = FACTOR * GH2(I)                                1133            GH(I) = FACTOR * GH2(I)                              
857              NMAX = NMAX2                                                          NMAX = NMAX2                                            
858          ENDIF                                                                  ENDIF                                                        
859                                                                                                                                                                    
860          DO 1144 I = 1, K                                                            DO 1144 I = 1, K  
861  1144        GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I))                1144        GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I))              
862                                                                                                                                                                    
863          RETURN                                                                RETURN                                                      
# Line 891  C -------------------------------------- Line 913  C --------------------------------------
913          ELSE IF (NMAX1 .GT. NMAX2) THEN                                          ELSE IF (NMAX1 .GT. NMAX2) THEN                                
914              K = NMAX2 * (NMAX2 + 2)                                                  K = NMAX2 * (NMAX2 + 2)                                    
915              L = NMAX1 * (NMAX1 + 2)                                                  L = NMAX1 * (NMAX1 + 2)                                    
916              DO 1155 I = K + 1, L                                                          DO 1155 I = K + 1, L    
917  1155            GH(I) = GH1(I)                                          1155            GH(I) = GH1(I)                                        
918              NMAX = NMAX1                                                            NMAX = NMAX1                                              
919          ELSE                                                                    ELSE                                                          
920              K = NMAX1 * (NMAX1 + 2)                                                  K = NMAX1 * (NMAX1 + 2)                                    
921              L = NMAX2 * (NMAX2 + 2)                                                  L = NMAX2 * (NMAX2 + 2)                                    
922              DO 1166 I = K + 1, L                                                          DO 1166 I = K + 1, L  
923  1166            GH(I) = FACTOR * GH2(I)                                  1166            GH(I) = FACTOR * GH2(I)                                
924              NMAX = NMAX2                                                            NMAX = NMAX2                                              
925          ENDIF                                                                    ENDIF                                                          
926                                                                                                                                                                    
927          DO 1177 I = 1, K                                                              DO 1177 I = 1, K
928  1177        GH(I) = GH1(I) + FACTOR * GH2(I)                            1177        GH(I) = GH1(I) + FACTOR * GH2(I)                          
929                                                                                                                                                                    
930          RETURN                                                                  RETURN                                                        
# Line 925  C ERA, EREQU and ERPOL as recommended by Line 947  C ERA, EREQU and ERPOL as recommended by
947  C ASTRONOMICAL UNION .  C ASTRONOMICAL UNION .
948  C-----------------------------------------------------------------  C-----------------------------------------------------------------
949          INTEGER TL1,TL2,TL3          INTEGER TL1,TL2,TL3
950          CHARACTER *258 TP1,TP2,TP3          CHARACTER (len=258) TP1,TP2,TP3
951          INTEGER L1,L2,L3          INTEGER L1,L2,L3
952          CHARACTER *258 P1,P2,P3          CHARACTER *258 P1,P2,P3
953          COMMON/PPATH/ L1,L2,L3,P1, P2, P3          COMMON/PPATH/ L1,L2,L3,P1, P2, P3
954          SAVE/PPATH/          SAVE/PPATH/
955    
956          COMMON/GENER/   UMR,ERA,AQUAD,BQUAD          COMMON/GENER/UMR,ERA,AQUAD,BQUAD
957          SAVE/GENER/          SAVE/GENER/
958    
         P1 = TP1  
         P2 = TP2  
         P3 = TP3  
959          L1 = TL1          L1 = TL1
960          L2 = TL2          L2 = TL2
961          L3 = TL3          L3 = TL3
962            P1 = TP1(1:L1)
963            P2 = TP2(1:L2)
964            P3 = TP3(1:L3)
965    
966          ERA=6371.2          ERA=6371.2
967          EREQU=6378.16          EREQU=6378.16
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.9
  ViewVC Help
Powered by ViewVC 1.1.23