/[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.4 by mocchiut, Tue Aug 11 12:58:22 2009 UTC
# Line 2  Line 2 
2       &          xl,icode,dip,dec)       &          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    
18        REAL              LATI,LONGI        REAL              LATI,LONGI
# Line 20  c--------------------------------------- Line 20  c---------------------------------------
20        SAVE /GENER/        SAVE /GENER/
21  C  C
22        CALL INITIZE        CALL INITIZE
23          ibbb=0        ibbb=0
24        ALOG2=ALOG(2.)        ALOG2=ALOG(2.)
25        ISTART=1        ISTART=1
26          lati=xlat          lati=xlat
# Line 32  c Line 32  c
32          CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS)          CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS)
33          CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1)          CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1)
34          DIP=ASIN(BDOWN/BABS)/UMR          DIP=ASIN(BDOWN/BABS)/UMR
35          DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR        DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR
36        RETURN        RETURN
37        END        END
38  c  c
# 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 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(144),GH2(120),GHA(144),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      WRITE(*,*) FILMOD(1)  c     FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat'
618  c      WRITE(*,*) FILMOD(2)  c     WRITE(*,*) FILMOD(1)
619  c      WRITE(*,*) FILMOD(3)  c     WRITE(*,*) FILMOD(2)
620    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 / 2000., 2005., 2010./
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      print *, "que"
656              IF (IER .NE. 0) STOP                            C--   GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS
657          CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)          CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER)  
658              IF (IER .NE. 0) STOP                            IF (IER .NE. 0) STOP                          
659  C-- DETERMINE IGRF COEFFICIENTS FOR YEAR        CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER)  
660          IF (L .LE. NUMYE-1) THEN                                IF (IER .NE. 0) STOP                    
661            CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,  c      print *, "quj"
662       1          NMAX2, GH2, NMAX, GHA)                          C--   DETERMINE IGRF COEFFICIENTS FOR YEAR
663          ELSE                      IF (L .LE. NUMYE-1) THEN                        
664            CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,               CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2,
665       1          GH2, NMAX, GHA)                                           1        NMAX2, GH2, NMAX, GHA)                        
666          ENDIF        ELSE              
667  C-- DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G           CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2,    
668          F0=0.D0       1        GH2, NMAX, GHA)                                    
669          DO 1234 J=1,3        ENDIF
670             F = GHA(J) * 1.D-5  c      print *, "quw"
671             F0 = F0 + F * F  C--   DETERMINE MAGNETIC DIPOL MOMENT AND COEFFIECIENTS G
672  1234    CONTINUE        F0=0.D0
673          DIMO = DSQRT(F0)        DO 1234 J=1,3
674             F = GHA(J) * 1.D-5
675          GH1(1) =  0.0           F0 = F0 + F * F
676          I=2             1234 CONTINUE
677          F0=1.D-5                        DIMO = DSQRT(F0)
678          IF(IS.EQ.0) F0=-F0        
679          SQRT2=SQRT(2.)              GH1(1) =  0.0
680          I=2          
681          F0=1.D-5                
682          IF(IS.EQ.0) F0=-F0
683          SQRT2=SQRT(2.)      
684    
685    c      print *, "quq"
686          
687        DO 9 N=1,NMAX                  DO 9 N=1,NMAX          
688          X = N           X = N
689          F0 = F0 * X * X / (4.D0 * X - 2.D0)                         F0 = F0 * X * X / (4.D0 * X - 2.D0)              
690          IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X           IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X
691          F = F0 * 0.5D0                                               F = F0 * 0.5D0                                    
692          IF(IS.EQ.0) F = F * SQRT2           IF(IS.EQ.0) F = F * SQRT2
693          GH1(I) = GHA(I-1) * F0           GH1(I) = GHA(I-1) * F0
694          I = I+1                                                   I = I+1                                        
695        DO 9 M=1,N                                               DO 9 M=1,N                                    
696          F = F * (X + M) / (X - M + 1.D0)                              F = F * (X + M) / (X - M + 1.D0)                
697          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))
698          GH1(I) = GHA(I-1) * F              GH1(I) = GHA(I-1) * F
699          GH1(I+1) = GHA(I) * F              GH1(I+1) = GHA(I) * F
700          I=I+2              I=I+2
701  9     CONTINUE                                             9       CONTINUE                                          
702          RETURN        RETURN
703          END        END
704  C  C    
705  C  C    
706          SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)                  SUBROUTINE GETSHC (IU, FSPEC, NMAX, ERAD, GH, IER)          
707                                                                                                                                                                    
708  C ===============================================================                C ===============================================================              
709  C                                                                                C                                                                              
710  C       Version 1.01                                                  C     Version 1.01                                                
711  C                                                                                C    
712  C       Reads spherical harmonic coefficients from the specified      C     Reads spherical harmonic coefficients from the specified    
713  C       file into an array.                                            C     file into an array.                                          
714  C                                                                                C    
715  C       Input:                                                        C     Input:                                                      
716  C           IU    - Logical unit number                                C     IU    - Logical unit number                              
717  C           FSPEC - File specification                                C     FSPEC - File specification                              
718  C                                                                                C    
719  C       Output:                                                        C     Output:                                                      
720  C           NMAX  - Maximum degree and order of model                  C     NMAX  - Maximum degree and order of model                
721  C           ERAD  - Earth's radius associated with the spherical      C     ERAD  - Earth's radius associated with the spherical    
722  C                   harmonic coefficients, in the same units as        C     harmonic coefficients, in the same units as      
723  C                   elevation                                          C     elevation                                        
724  C           GH    - Schmidt quasi-normal internal spherical            C     GH    - Schmidt quasi-normal internal spherical          
725  C                   harmonic coefficients                              C     harmonic coefficients                            
726  C           IER   - Error number: =  0, no error                      C     IER   - Error number: =  0, no error                    
727  C                                 = -2, records out of order          C     = -2, records out of order        
728  C                                 = FORTRAN run-time error number      C     = FORTRAN run-time error number    
729  C                                                                      C    
730  C       A. Zunde                                                      C     A. Zunde                                                    
731  C       USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225      C     USGS, MS 964, Box 25046 Federal Center, Denver, CO  80225    
732  C                                                                                C    
733  C ===============================================================                C     ===============================================================              
734                                                                                          
735          CHARACTER  FSPEC*(*), FOUT*258                                            CHARACTER  FSPEC*(*), FOUT*258
736          DIMENSION       GH(*)                                                DIMENSION       GH(*)                                        
737  C ---------------------------------------------------------------                C     ---------------------------------------------------------------              
738  C       Open coefficient file. Read past first header record.          C     Open coefficient file. Read past first header record.        
739  C       Read degree and order of model and Earth's radius.            C     Read degree and order of model and Earth's radius.          
740  C ---------------------------------------------------------------                C     ---------------------------------------------------------------              
741        WRITE(FOUT,667) FSPEC        WRITE(FOUT,667) FSPEC
742  c 667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)  c     667  FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12)
743  667   FORMAT(A258)   667  FORMAT(A258)
744    c      print *," gui"
745        OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)            OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999)    
746    c      print *," gua"
747        READ (IU, *, IOSTAT=IER, ERR=999)        READ (IU, *, IOSTAT=IER, ERR=999)
748    c      print *," gue"
749        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD        READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD
750    c      print *," guo"
751  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
752  C       Read the coefficient file, arranged as follows:                C       Read the coefficient file, arranged as follows:              
753  C                                                                                C                                                                              
# Line 752  C       coefficient.                     Line 769  C       coefficient.                    
769  C ---------------------------------------------------------------                C ---------------------------------------------------------------              
770                                                                                                                                                                    
771          I = 0                                                                  I = 0                                                        
772          DO 2211 NN = 1, NMAX                                                        DO 2211 NN = 1, NMAX              
773              DO 2233 MM = 0, NN                                                          DO 2233 MM = 0, NN            
774                  READ (IU, *, IOSTAT=IER, ERR=999) N, M, G, H                          READ (IU, *, IOSTAT=IER, ERR=999) N, M, G, H        
775                  IF (NN .NE. N .OR. MM .NE. M) THEN                                    IF (NN .NE. N .OR. MM .NE. M) THEN                  
776                      IER = -2                                                              IER = -2                                        
# Line 767  C -------------------------------------- Line 784  C --------------------------------------
784                  ENDIF                                                                  ENDIF                                                
785  2233        CONTINUE                                                      2233        CONTINUE                                                    
786  2211    CONTINUE                                                          2211    CONTINUE                                                        
787    c      print *," guj"
788                                                                                                                                                                    
789  999     CLOSE (IU)                                                    999     CLOSE (IU)                                                  
790    c      print *," guw IER",IER
791          if ( IER .eq. -1 ) IER = 0 ! gfortran 4.1.2 bug workaround... hoping not to create problems with other versions
792                                                                                                                                                                    
793          RETURN                                                                RETURN                                                      
794          END                                                                    END                                                          
# Line 824  C -------------------------------------- Line 843  C --------------------------------------
843          ELSE IF (NMAX1 .GT. NMAX2) THEN                                        ELSE IF (NMAX1 .GT. NMAX2) THEN                              
844              K = NMAX2 * (NMAX2 + 2)                                                K = NMAX2 * (NMAX2 + 2)                                  
845              L = NMAX1 * (NMAX1 + 2)                                                L = NMAX1 * (NMAX1 + 2)                                  
846              DO 1122 I = K + 1, L                                                        DO 1122 I = K + 1, L    
847  1122            GH(I) = GH1(I) + FACTOR * (-GH1(I))                    1122            GH(I) = GH1(I) + FACTOR * (-GH1(I))                  
848              NMAX = NMAX1                                                          NMAX = NMAX1                                            
849          ELSE                                                                  ELSE                                                        
850              K = NMAX1 * (NMAX1 + 2)                                                K = NMAX1 * (NMAX1 + 2)                                  
851              L = NMAX2 * (NMAX2 + 2)                                                L = NMAX2 * (NMAX2 + 2)                                  
852              DO 1133 I = K + 1, L                                                        DO 1133 I = K + 1, L
853  1133            GH(I) = FACTOR * GH2(I)                                1133            GH(I) = FACTOR * GH2(I)                              
854              NMAX = NMAX2                                                          NMAX = NMAX2                                            
855          ENDIF                                                                  ENDIF                                                        
856                                                                                                                                                                    
857          DO 1144 I = 1, K                                                            DO 1144 I = 1, K  
858  1144        GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I))                1144        GH(I) = GH1(I) + FACTOR * (GH2(I) - GH1(I))              
859                                                                                                                                                                    
860          RETURN                                                                RETURN                                                      
# Line 891  C -------------------------------------- Line 910  C --------------------------------------
910          ELSE IF (NMAX1 .GT. NMAX2) THEN                                          ELSE IF (NMAX1 .GT. NMAX2) THEN                                
911              K = NMAX2 * (NMAX2 + 2)                                                  K = NMAX2 * (NMAX2 + 2)                                    
912              L = NMAX1 * (NMAX1 + 2)                                                  L = NMAX1 * (NMAX1 + 2)                                    
913              DO 1155 I = K + 1, L                                                          DO 1155 I = K + 1, L    
914  1155            GH(I) = GH1(I)                                          1155            GH(I) = GH1(I)                                        
915              NMAX = NMAX1                                                            NMAX = NMAX1                                              
916          ELSE                                                                    ELSE                                                          
917              K = NMAX1 * (NMAX1 + 2)                                                  K = NMAX1 * (NMAX1 + 2)                                    
918              L = NMAX2 * (NMAX2 + 2)                                                  L = NMAX2 * (NMAX2 + 2)                                    
919              DO 1166 I = K + 1, L                                                          DO 1166 I = K + 1, L  
920  1166            GH(I) = FACTOR * GH2(I)                                  1166            GH(I) = FACTOR * GH2(I)                                
921              NMAX = NMAX2                                                            NMAX = NMAX2                                              
922          ENDIF                                                                    ENDIF                                                          
923                                                                                                                                                                    
924          DO 1177 I = 1, K                                                              DO 1177 I = 1, K
925  1177        GH(I) = GH1(I) + FACTOR * GH2(I)                            1177        GH(I) = GH1(I) + FACTOR * GH2(I)                          
926                                                                                                                                                                    
927          RETURN                                                                  RETURN                                                        
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.4
  ViewVC Help
Powered by ViewVC 1.1.23