/[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.2 by mocchiut, Fri Apr 27 10:35:35 2007 UTC revision 1.12 by mocchiut, Mon Jan 19 12:32:13 2015 UTC
# Line 1  Line 1 
         subroutine igrf_sub(xlat,xlong,year,height,  
      &          xl,icode,dip,dec)  
 c----------------------------------------------------------------  
 c   INPUT:  
 c       xlat    geodatic latitude in degrees  
 c       xlong   geodatic longitude in degrees  
 c       year    decimal year (year+month/12.0-0.5 or year+day-of-year/365  
 c               or 366 if leap year)  
 c       height  height in km  
 c   OUTPUT:  
 c       xl      L value  
 c       icode   =1  L is correct; =2  L is not correct;  
 c               =3  an approximation is used  
 c       dip     geomagnetic inclination in degrees  
 c       dec     geomagnetic declination in degress  
 c----------------------------------------------------------------  
   
       REAL              LATI,LONGI  
       COMMON/GENER/     UMR,ERA,AQUAD,BQUAD  
       SAVE /GENER/  
1  C  C
       CALL INITIZE  
         ibbb=0  
       ALOG2=ALOG(2.)  
       ISTART=1  
         lati=xlat  
         longi=xlong  
 c  
 C----------------CALCULATE PROFILES-----------------------------------  
 c  
         CALL FELDCOF(YEAR,DIMO)  
         CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS)  
         CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1)  
         DIP=ASIN(BDOWN/BABS)/UMR  
         DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR  
       RETURN  
       END  
 c  
 c  
2  C SHELLIG.FOR, Version 2.0, January 1992  C SHELLIG.FOR, Version 2.0, January 1992
3  C  C
4  C 11/01/91-DKB- SHELLG: lowest starting point for B0 search is 2    C 11/01/91-DKB- SHELLG: lowest starting point for B0 search is 2  
# Line 191  C                          APPROXIMATION Line 153  C                          APPROXIMATION
153  C          B0           MAGNETIC FIELD STRENGTH IN GAUSS  C          B0           MAGNETIC FIELD STRENGTH IN GAUSS
154  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
155        DIMENSION         V(3),U(3,3),P(8,100),SP(3)        DIMENSION         V(3),U(3,3),P(8,100),SP(3)
156        COMMON            X(3),H(144)        COMMON            X(3),H(196)
157        COMMON/FIDB0/     SP        COMMON/FIDB0/     SP
158        SAVE /FIDB0/        SAVE /FIDB0/
159        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
160        SAVE /GENER/        SAVE /GENER/
161          REAL FLS
162  C  C
163  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3  C-- RMIN, RMAX ARE BOUNDARIES FOR IDENTIFICATION OF ICODE=2 AND 3
164  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 166  C-- STEQ IS STEP SIZE FOR INTEGRATION
166  C  C
167        DATA RMIN,RMAX    /0.05,1.01/        DATA RMIN,RMAX    /0.05,1.01/
168        DATA STEP,STEQ    /0.20,0.03/        DATA STEP,STEQ    /0.20,0.03/
169          BEQU=1.E10        BEQU=1.E10
170          FLS = FL
171  C*****ENTRY POINT  SHELLG  TO BE USED WITH GEODETIC CO-ORDINATES  C*****ENTRY POINT  SHELLG  TO BE USED WITH GEODETIC CO-ORDINATES
172        RLAT=GLAT*UMR        RLAT=GLAT*UMR
173        CT=SIN(RLAT)                                                      CT=SIN(RLAT)                                              
# Line 304  C*****INNER LOOP (FOR QUADRATURE)       Line 268  C*****INNER LOOP (FOR QUADRATURE)      
268        HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)                                    HLI=0.5*(((C3*T+C2)*T+C1)*T+C0)                            
269        ZQ=Z*Z        ZQ=Z*Z
270        R=HLI+SQRT(HLI*HLI+ZQ)        R=HLI+SQRT(HLI*HLI+ZQ)
271          IF(R.NE.R)THEN
272             FL = FLS
273             RETURN
274          ENDIF
275        IF(R.LE.RMIN)GOTO30                                      IF(R.LE.RMIN)GOTO30                              
276        RQ=R*R        RQ=R*R
277        FF=SQRT(1.+3.*ZQ/RQ)                                      FF=SQRT(1.+3.*ZQ/RQ)                              
278        RADIK=B0-((D2*T+D1)*T+D0)*R*RQ*FF                        RADIK=B0-((D2*T+D1)*T+D0)*R*RQ*FF                
279        IF(R-RMAX)44,44,45                                        IF((R-RMAX).le.0.) goto 44                                
280          IF((R-RMAX).gt.0.) goto 45                                
281  45    ICODE=2                                            45    ICODE=2                                          
282        RADIK=RADIK-12.*(R-RMAX)**2                              RADIK=RADIK-12.*(R-RMAX)**2                      
283  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 310  C-- The minimal allowable value of FI wa
310  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.
311  C-- D. Bilitza, Nov 87.  C-- D. Bilitza, Nov 87.
312  C  C
313  11      FI=0.5*ABS(FI)/SQRT(B0)+1E-12                        11         FI=0.5*ABS(FI)/SQRT(B0)+1E-12                      
314  C  C
315  C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.    C*****COMPUTE L FROM B AND I.  SAME AS CARMEL IN INVAR.  
316  C  C
317  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.  C-- Correct dipole moment is used here. D. Bilitza, Nov 87.
318  C  C
319        DIMOB0=DIMO/B0        DIMOB0=DIMO/B0
320          arg1=alog(FI)        arg1=alog(FI)
321          arg2=alog(DIMOB0)        arg2=alog(DIMOB0)
322  c       arg = FI*FI*FI/DIMOB0  c    arg = FI*FI*FI/DIMOB0
323  c       if(abs(arg).gt.88.0) arg=88.0  c    if(abs(arg).gt.88.0) arg=88.0
324        XX=3*arg1-arg2        XX=3*arg1-arg2
325        IF(XX.GT.23.0) GOTO 776          IF(XX.GT.23.0) GOTO 776  
326        IF(XX.GT.11.7) GOTO 775          IF(XX.GT.11.7) GOTO 775  
327        IF(XX.GT.+3.0) GOTO 774            IF(XX.GT.+3.0) GOTO 774    
328        IF(XX.GT.-3.0) GOTO 773          IF(XX.GT.-3.0) GOTO 773  
329        IF(XX.GT.-22.) GOTO 772          IF(XX.GT.-22.) GOTO 772  
330    771 GG=3.33338E-1*XX+3.0062102E-1                                  c  771 GG=3.33338E-1*XX+3.0062102E-1                                
331          GG=3.33338E-1*XX+3.0062102E-1                                
332        GOTO777                                                                  GOTO777                                                          
333    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+  
334       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 362  C* SUBROUTINE USED FOR FIELD LINE TRACIN
362  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *  C* CALLS ENTRY POINT FELDI IN GEOMAGNETIC FIELD SUBROUTINE FELDG   *
363  C*******************************************************************  C*******************************************************************
364        DIMENSION         P(7),U(3,3)        DIMENSION         P(7),U(3,3)
365        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
366  C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES            C*****XM,YM,ZM  ARE GEOMAGNETIC CARTESIAN INVERSE CO-ORDINATES          
367        ZM=P(3)                                                                  ZM=P(3)                                                          
368        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 444  C                 TO THE LOCAL GEODETIC
444  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST  C                 POINTING IN THE TANGENTIAL PLANE TO THE NORTH, EAST
445  C                 AND DOWNWARD.    C                 AND DOWNWARD.  
446  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
447        DIMENSION         V(3),B(3),G(144)        DIMENSION         V(3),B(3),G(196)
448        CHARACTER*258      NAME        CHARACTER*258      NAME
449        INTEGER NMAX        INTEGER NMAX
450        REAL TIME        REAL TIME
451        COMMON            XI(3),H(144)        COMMON            XI(3),H(196)
452        COMMON/MODEL/     G,NMAX,TIME,NAME        COMMON/MODEL/     G,NMAX,TIME,NAME
453        SAVE/MODEL/        SAVE/MODEL/
454        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD        COMMON/GENER/     UMR,ERA,AQUAD,BQUAD
# Line 528  C*****ENTRY POINT  FELDI  USED FOR L COM Line 498  C*****ENTRY POINT  FELDI  USED FOR L COM
498        Y=XI(2)*F                                                                Y=XI(2)*F                                                        
499        Z=XI(3)*(F+F)                                                            Z=XI(3)*(F+F)                                                    
500        I=I-2                                                                    I=I-2                                                            
501        IF(I-1)5,4,2                                                        c      print *,' I ',I
502          IF((I-1).lt.0) goto 5
503          IF((I-1).eq.0) goto 4
504          IF((I-1).gt.0) goto 2
505  2     DO 3 M=3,I,2                                                        2     DO 3 M=3,I,2                                                      
506        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))          
507       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 533  C*****ENTRY POINT  FELDI  USED FOR L COM
533        END                                                                      END                                                              
534  C  C
535  C  C
536          SUBROUTINE FELDCOF(YEAR,DIMO)        SUBROUTINE FELDCOF(YEAR,DIMO)
537  C------------------------------------------------------------------------  C------------------------------------------------------------------------
538  C  DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS  C     DETERMINES COEFFICIENTS AND DIPOL MOMENT FROM IGRF MODELS
539  C  C    
540  C       INPUT:  YEAR    DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO  C     INPUT:  YEAR    DECIMAL YEAR FOR WHICH GEOMAGNETIC FIELD IS TO
541  C                       BE CALCULATED  C     BE CALCULATED
542  C       OUTPUT: DIMO    GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED  C     OUTPUT: DIMO    GEOMAGNETIC DIPOL MOMENT IN GAUSS (NORMALIZED
543  C                       TO EARTH'S RADIUS) AT THE TIME (YEAR)  C     TO EARTH'S RADIUS) AT THE TIME (YEAR)
544  C  D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771,  C     D. BILITZA, NSSDC, GSFC, CODE 633, GREENBELT, MD 20771,
545  C       (301)286-9536   NOV 1987.  C     (301)286-9536   NOV 1987.
546  C  ### updated to IGRF-2000 version -dkb- 5/31/2000  C     ### updated to IGRF-2000 version -dkb- 5/31/2000
547  C  ### updated to IGRF-2005 version -dkb- 3/24/2000  C     ### updated to IGRF-2005 version -dkb- 3/24/2000
548  C-----------------------------------------------------------------------  C-----------------------------------------------------------------------
549          CHARACTER*258    FIL1, FIL2                  CHARACTER*258    FIL1, FIL2          
550          CHARACTER*258    FILMOD        CHARACTER*258    FILMOD
551  C ### FILMOD, DTEMOD arrays +1        DIMENSION       GH1(196),GH2(196),GHA(196),FILMOD(2)
552  c        DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(14),DTEMOD(14)        DOUBLE PRECISION X,F0,F
553          DIMENSION       GH1(144),GH2(120),GHA(144),FILMOD(3),DTEMOD(3)        INTEGER L1,L2,I1
554          DOUBLE PRECISION X,F0,F        INTEGER NMAX
555          INTEGER L1,L2,L3        REAL TIME
556          INTEGER NMAX        CHARACTER *258 P1,P2
557          REAL TIME        COMMON/PPATH/ I1,L1,L2,P1,P2
558          CHARACTER *258 P1,P2,P3        SAVE/PPATH/
559          COMMON/PPATH/ L1,L2,L3,P1, P2, P3        COMMON/MODEL/   GH1,NMAX,TIME,FIL1
560          SAVE/PPATH/        SAVE/MODEL/
561          COMMON/MODEL/   GH1,NMAX,TIME,FIL1        COMMON/GENER/   UMR,ERAD,AQUAD,BQUAD
562          SAVE/MODEL/        SAVE/GENER/
563          COMMON/GENER/   UMR,ERAD,AQUAD,BQUAD  c      print *, "qui"
         SAVE/GENER/  
 C ### updated to 2005  
 C      CHARACTER COEFPATH*80, COEF1*80, COEF2*80, COEF3*80  
   
 c      COEFPATH  = 'OrbitalInfo/src/'  
 c      COEF1 = 'dgrf00.dat'  
 c      COEF2 = 'igrf05.dat'  
 c      COEF3 = 'igrf05s.dat'  
 c      COEF1 = COEFPATH(1:16)//COEF1  
 c      COEF2 = COEFPATH(1:16)//COEF2  
 c      COEF3 = COEFPATH(1:16)//COEF3  
 c      FILMOD(1) = COEF1  
 c      FILMOD(2) = COEF2  
 c      FILMOD(3) = COEF3  
564        FILMOD(1) = P1(1:L1)        FILMOD(1) = P1(1:L1)
565        FILMOD(2) = P2(1:L2)        FILMOD(2) = P2(1:L2)
566        FILMOD(3) = P3(1:L3)  c      print *, "qua"
567  c      FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat'  C    
568  c      FILMOD(2) = 'OrbitalInfo/src/igrf05.dat'  C     IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION
569  c      FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'  C     IU  IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS
570  c      WRITE(*,*) FILMOD(1)  C    
571  c      WRITE(*,*) FILMOD(2)        IU = 10
572  c      WRITE(*,*) FILMOD(3)        IS = 0
573  c      DATA   FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/  C--   DETERMINE IGRF-YEARS FOR INPUT-YEAR
574        DATA   DTEMOD / 2000., 2005., 2010./        TIME = YEAR
575  c        IYEA = INT(YEAR/5.)*5
576  c        DATA            FILMOD /'dgrf45.dat', 'dgrf50.dat',                    L = IYEA + 5
577  c     1                  'dgrf55.dat', 'dgrf60.dat', 'dgrf65.dat',        C
578  c     2                  'dgrf70.dat', 'dgrf75.dat', 'dgrf80.dat',              DTE1 = REAL(IYEA)  
579  c     3                  'dgrf85.dat', 'dgrf90.dat', 'dgrf95.dat',        FIL1 = FILMOD(1)  
580  c     4                  'dgrf00.dat','igrf05.dat','igrf05s.dat'/        DTE2 = REAL(L)
581  c        DATA   DTEMOD / 1945., 1950., 1955., 1960., 1965., 1970.,        FIL2 = FILMOD(2)
582  c     1       1975., 1980., 1985., 1990., 1995., 2000.,2005.,2010./        c      print *,'IYEA ',IYEA,' L ',L
583  C  c      WRITE(*,*) FIL1
584  C ### numye = numye + 1 ; is number of years represented by IGRF  c      WRITE(*,*) FIL2
585  C  c      print *, "que"
586  c        NUMYE=13  C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
587          NUMYE=2        CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
588          IF (IER .NE. 0) STOP                          
589  C  c      print *, "quessss"
590  C  IS=0 FOR SCHMIDT NORMALIZATION   IS=1 GAUSS NORMALIZATION        CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
591  C  IU  IS INPUT UNIT NUMBER FOR IGRF COEFFICIENT SETS        IF (IER .NE. 0) STOP                    
592  C  c      print *, "quj"
593          IU = 10  C--   DETERMINE IGRF COEFFICIENTS FOR YEAR
594          IS = 0        IF (I1.EQ.0) THEN      
595  C-- DETERMINE IGRF-YEARS FOR INPUT-YEAR           CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,
596          TIME = YEAR       1        NMAX2, GH2, NMAX, GHA)                        
597          IYEA = INT(YEAR/5.)*5        ELSE              
598  c        L = (IYEA - 1945)/5 + 1           CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,    
599          L = (IYEA - 2000)/5 + 1       1        GH2, NMAX, GHA)                                    
600          IF(L.LT.1) L=1        ENDIF
601          IF(L.GT.NUMYE) L=NUMYE          c      print *, "quw"
602          DTE1 = DTEMOD(L)    C--   DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G
603          FIL1 = FILMOD(L)          F0=0.D0
604          DTE2 = DTEMOD(L+1)        DO 1234 J=1,3
605          FIL2 = FILMOD(L+1)           F = GHA(J) * 1.D-5
606  C-- GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS           F0 = F0 + F * F
607          CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)     1234 CONTINUE
608              IF (IER .NE. 0) STOP                                  DIMO = REAL(DSQRT(F0))
609          CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)          
610              IF (IER .NE. 0) STOP                            GH1(1) =  0.0
611  C-- DETERMINE IGRF COEFFICIENTS FOR YEAR        I=2          
612          IF (L .LE. NUMYE-1) THEN                                F0=1.D-5                
613            CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,        IF(IS.EQ.0) F0=-F0
614       1          NMAX2, GH2, NMAX, GHA)                                SQRT2=SQRT(2.)      
         ELSE                
           CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,      
      1          GH2, NMAX, GHA)                                      
         ENDIF  
 C-- DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G  
         F0=0.D0  
         DO 1234 J=1,3  
            F = GHA(J) * 1.D-5  
            F0 = F0 + F * F  
 1234    CONTINUE  
         DIMO = DSQRT(F0)  
   
         GH1(1) =  0.0  
         I=2            
         F0=1.D-5                  
         IF(IS.EQ.0) F0=-F0  
         SQRT2=SQRT(2.)        
615    
616    c      print *, "quq"
617          
618        DO 9 N=1,NMAX                  DO 9 N=1,NMAX          
619          X = N           X = N
620          F0 = F0 * X * X / (4.D0 * X - 2.D0)                         F0 = F0 * X * X / (4.D0 * X - 2.D0)              
621          IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X
622          F = F0 * 0.5D0                                               F = F0 * 0.5D0                                    
623          IF(IS.EQ.0) F = F * SQRT2           IF(IS.EQ.0) F = F * SQRT2
624          GH1(I) = GHA(I-1) * F0           GH1(I) = GHA(I-1) * REAL(F0)
625          I = I+1                                                   I = I+1                                        
626        DO 9 M=1,N                                               DO 9 M=1,N                                    
627          F = F * (X + M) / (X - M + 1.D0)                              F = F * (X + M) / (X - M + 1.D0)                
628          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))
629          GH1(I) = GHA(I-1) * F              GH1(I) = GHA(I-1) * REAL(F)
630          GH1(I+1) = GHA(I) * F              GH1(I+1) = GHA(I) * REAL(F)
631          I=I+2              I=I+2
632  9     CONTINUE                                             9       CONTINUE                                          
633          RETURN        RETURN
634          END        END
635  C  C    
636  C  C    
637          SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)                  SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)          
638                                                                                                                                                                    
639  C ===============================================================                C ===============================================================              
640  C                                                                                C                                                                              
641  C       Version 1.01                                                  C     Version 1.01                                                
642  C                                                                                C    
643  C       Reads spherical harmonic coefficients from the specified      C     Reads spherical harmonic coefficients from the specified    
644  C       file into an array.                                            C     file into an array.                                          
645  C                                                                                C    
646  C       Input:                                                        C     Input:                                                      
647  C           IU    - Logical unit number                                C     IU    - Logical unit number                              
648  C           FSPEC - File specification                                C     FSPEC - File specification                              
649  C                                                                                C    
650  C       Output:                                                        C     Output:                                                      
651  C           NMAX  - Maximum degree and order of model                  C     NMAX  - Maximum degree and order of model                
652  C           ERAD  - Earth's radius associated with the spherical      C     ERAD  - Earth's radius associated with the spherical    
653  C                   harmonic coefficients, in the same units as        C     harmonic coefficients, in the same units as      
654  C                   elevation                                          C     elevation                                        
655  C           GH    - Schmidt quasi-normal internal spherical            C     GH    - Schmidt quasi-normal internal spherical          
656  C                   harmonic coefficients                              C     harmonic coefficients                            
657  C           IER   - Error number: =  0, no error                      C     IER   - Error number: =  0, no error                    
658  C                                 = -2, records out of order          C     = -2, records out of order        
659  C                                 = FORTRAN run-time error number      C     = FORTRAN run-time error number    
660  C                                                                      C    
661  C       A. Zunde                                                      C     A. Zunde                                                    
662  C       USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225      C     USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225    
663  C                                                                                C    
664  C ===============================================================                C     ===============================================================              
665                                                                                          
666          CHARACTER  FSPEC*(*), FOUT*258                                            CHARACTER  FSPEC*(*), FOUT*258
667          DIMENSION       GH(*)                                                DIMENSION       GH(*)                                        
668  C ---------------------------------------------------------------                C     ---------------------------------------------------------------              
669  C       Open coefficient file. Read past first header record.          C     Open coefficient file. Read past first header record.        
670  C       Read degree and order of model and Earth's radius.            C     Read degree and order of model and Earth's radius.          
671  C ---------------------------------------------------------------                C     ---------------------------------------------------------------              
672        WRITE(FOUT,667) FSPEC        WRITE(FOUT,667) FSPEC
673  c 667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)  c     667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)
674  667   FORMAT(A258)   667  FORMAT(A258)
675    c      print *," gui"
676        OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)            OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)    
677    c      print *," gua"
678        READ (IU, *, IOSTAT=IER, ERR=999)        READ (IU, *, IOSTAT=IER, ERR=999)
679    c      print *," gue"
680        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD
681    c      print *," guo"
682  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
683  C       Read the coefficient file, arranged as follows:                C       Read the coefficient file, arranged as follows:              
684  C                                                                                C                                                                              
# Line 752  C       coefficient.                     Line 700  C       coefficient.                    
700  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
701                                                                                                                                                                    
702          I = 0                                                                  I = 0                                                        
703          DO 2211 NN = 1, NMAX                                                        DO 2211 NN = 1, NMAX              
704              DO 2233 MM = 0, NN                                                          DO 2233 MM = 0, NN            
705                  READ (IU, *, IOSTAT=IER, ERR=999) N, M, G, H                          READ (IU, *, IOSTAT=IER, ERR=999) N, M, G, H        
706                  IF (NN .NE. N .OR. MM .NE. M) THEN                                    IF (NN .NE. N .OR. MM .NE. M) THEN                  
707                      IER = -2                                                              IER = -2                                        
# Line 767  C -------------------------------------- Line 715  C --------------------------------------
715                  ENDIF                                                                  ENDIF                                                
716  2233        CONTINUE                                                      2233        CONTINUE                                                    
717  2211    CONTINUE                                                          2211    CONTINUE                                                        
718    c      print *," guj"
719                                                                                                                                                                    
720  999     CLOSE (IU)                                                    999     CLOSE (IU)                                                  
721    c      print *," guw IER",IER
722          if ( IER .eq. -1 ) IER = 0 ! gfortran 4.1.2 bug workaround... hoping not to create problems with other versions
723                                                                                                                                                                    
724          RETURN                                                                RETURN                                                      
725          END                                                                    END                                                          
# Line 824  C -------------------------------------- Line 774  C --------------------------------------
774          ELSE IF (NMAX1 .GT. NMAX2) THEN                                        ELSE IF (NMAX1 .GT. NMAX2) THEN                              
775              K = NMAX2 * (NMAX2 + 2)                                                K = NMAX2 * (NMAX2 + 2)                                  
776              L = NMAX1 * (NMAX1 + 2)                                                L = NMAX1 * (NMAX1 + 2)                                  
777              DO 1122 I = K + 1, L                                                        DO 1122 I = K + 1, L    
778  1122            GH(I) = GH1(I) + FACTOR * (-GH1(I))                    1122            GH(I) = GH1(I) + FACTOR * (-GH1(I))                  
779              NMAX = NMAX1                                                          NMAX = NMAX1                                            
780          ELSE                                                                  ELSE                                                        
781              K = NMAX1 * (NMAX1 + 2)                                                K = NMAX1 * (NMAX1 + 2)                                  
782              L = NMAX2 * (NMAX2 + 2)                                                L = NMAX2 * (NMAX2 + 2)                                  
783              DO 1133 I = K + 1, L                                                        DO 1133 I = K + 1, L
784  1133            GH(I) = FACTOR * GH2(I)                                1133            GH(I) = FACTOR * GH2(I)                              
785              NMAX = NMAX2                                                          NMAX = NMAX2                                            
786          ENDIF                                                                  ENDIF                                                        
787                                                                                                                                                                    
788          DO 1144 I = 1, K                                                            DO 1144 I = 1, K  
789  1144        GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I))                1144        GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I))              
790                                                                                                                                                                    
791          RETURN                                                                RETURN                                                      
# Line 891  C -------------------------------------- Line 841  C --------------------------------------
841          ELSE IF (NMAX1 .GT. NMAX2) THEN                                          ELSE IF (NMAX1 .GT. NMAX2) THEN                                
842              K = NMAX2 * (NMAX2 + 2)                                                  K = NMAX2 * (NMAX2 + 2)                                    
843              L = NMAX1 * (NMAX1 + 2)                                                  L = NMAX1 * (NMAX1 + 2)                                    
844              DO 1155 I = K + 1, L                                                          DO 1155 I = K + 1, L    
845  1155            GH(I) = GH1(I)                                          1155            GH(I) = GH1(I)                                        
846              NMAX = NMAX1                                                            NMAX = NMAX1                                              
847          ELSE                                                                    ELSE                                                          
848              K = NMAX1 * (NMAX1 + 2)                                                  K = NMAX1 * (NMAX1 + 2)                                    
849              L = NMAX2 * (NMAX2 + 2)                                                  L = NMAX2 * (NMAX2 + 2)                                    
850              DO 1166 I = K + 1, L                                                          DO 1166 I = K + 1, L  
851  1166            GH(I) = FACTOR * GH2(I)                                  1166            GH(I) = FACTOR * GH2(I)                                
852              NMAX = NMAX2                                                            NMAX = NMAX2                                              
853          ENDIF                                                                    ENDIF                                                          
854                                                                                                                                                                    
855          DO 1177 I = 1, K                                                              DO 1177 I = 1, K
856  1177        GH(I) = GH1(I) + FACTOR * GH2(I)                            1177        GH(I) = GH1(I) + FACTOR * GH2(I)                          
857                                                                                                                                                                    
858          RETURN                                                                  RETURN                                                        
859          END                                                                      END                                                            
860  C  C
861  C  C
862          SUBROUTINE INITIZE(TP1,TL1,TP2,TL2,TP3,TL3)          SUBROUTINE INITIZE(ISSEC,TP1,TL1,TP2,TL2)
863  C----------------------------------------------------------------  C----------------------------------------------------------------
864  C Initializes the parameters in COMMON/GENER/  C Initializes the parameters in COMMON/GENER/
865  C  C
# Line 924  C Line 874  C
874  C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL  C ERA, EREQU and ERPOL as recommended by the INTERNATIONAL
875  C ASTRONOMICAL UNION .  C ASTRONOMICAL UNION .
876  C-----------------------------------------------------------------  C-----------------------------------------------------------------
877          INTEGER TL1,TL2,TL3          INTEGER TL1,TL2,ISSEC
878          CHARACTER *258 TP1,TP2,TP3          CHARACTER (len=*) :: TP1,TP2
879          INTEGER L1,L2,L3          INTEGER L1,L2
880          CHARACTER *258 P1,P2,P3          CHARACTER *258 P1,P2
881          COMMON/PPATH/ L1,L2,L3,P1, P2, P3          COMMON/PPATH/ I1,L1,L2,P1,P2
882          SAVE/PPATH/          SAVE/PPATH/
883    
884          COMMON/GENER/UMR,ERA,AQUAD,BQUAD          COMMON/GENER/UMR,ERA,AQUAD,BQUAD
885          SAVE/GENER/          SAVE/GENER/
886    
887            I1 = ISSEC
888          L1 = TL1          L1 = TL1
889          L2 = TL2          L2 = TL2
890          L3 = TL3  
891          P1 = TP1(1:L1)          P1 = TP1(1:L1)
892          P2 = TP2(1:L2)          P2 = TP2(1:L2)
         P3 = TP3(1:L3)  
893    
894          ERA=6371.2          ERA=6371.2
895          EREQU=6378.16          EREQU=6378.16
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.12
  ViewVC Help
Powered by ViewVC 1.1.23