--- DarthVader/OrbitalInfo/src/igrf_sub.for 2012/05/15 14:31:14 1.6 +++ DarthVader/OrbitalInfo/src/igrf_sub.for 2014/08/29 07:06:41 1.10 @@ -1,5 +1,5 @@ - subroutine igrf_sub(xlat,xlong,year,height, - & xl,icode,dip,dec) +c subroutine igrf_sub(xlat,xlong,year,height, +c & xl,icode,dip,dec) c---------------------------------------------------------------- c INPUT: c xlat geodatic latitude in degrees @@ -14,27 +14,27 @@ c dip geomagnetic inclination in degrees c dec geomagnetic declination in degress c---------------------------------------------------------------- - - REAL LATI,LONGI - COMMON/GENER/ UMR,ERA,AQUAD,BQUAD - SAVE /GENER/ -C - CALL INITIZE - ibbb=0 - ALOG2=ALOG(2.) - ISTART=1 - lati=xlat - longi=xlong +c +c REAL LATI,LONGI +c COMMON/GENER/ UMR,ERA,AQUAD,BQUAD +c SAVE /GENER/ +C +c CALL INITIZE +c ibbb=0 +c ALOG2=ALOG(2.) +c ISTART=1 +c lati=xlat +c 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 CALL FELDCOF(YEAR,DIMO) +c CALL FELDG(LATI,LONGI,HEIGHT,BNORTH,BEAST,BDOWN,BABS) +c CALL SHELLG(LATI,LONGI,HEIGHT,DIMO,XL,ICODE,BAB1) +c DIP=ASIN(BDOWN/BABS)/UMR +c DEC=ASIN(BEAST/SQRT(BEAST*BEAST+BNORTH*BNORTH))/UMR +c RETURN +c END c c C SHELLIG.FOR, Version 2.0, January 1992 @@ -607,17 +607,17 @@ c FILMOD(1) = COEF1 c FILMOD(2) = COEF2 c FILMOD(3) = COEF3 - print *, "qui" +c print *, "qui" FILMOD(1) = P1(1:L1) FILMOD(2) = P2(1:L2) FILMOD(3) = P3(1:L3) - print *, "qua" +c print *, "qua" c FILMOD(1) = 'OrbitalInfo/src/dgrf00.dat' c FILMOD(2) = 'OrbitalInfo/src/igrf05.dat' c FILMOD(3) = 'OrbitalInfo/src/igrf05s.dat' - WRITE(*,*) FILMOD(1) - WRITE(*,*) FILMOD(2) - WRITE(*,*) FILMOD(3) +c WRITE(*,*) FILMOD(1) +c WRITE(*,*) FILMOD(2) +c WRITE(*,*) FILMOD(3) c DATA FILMOD / 'dgrf00.dat', 'igrf05.dat', 'igrf05s.dat'/ DATA DTEMOD / 2005., 2010., 2015./ c @@ -633,7 +633,7 @@ C c NUMYE=13 NUMYE=2 - print *, "quo" +c print *, "quo" C C IS=0 FOR SCHMIDT NORMALIZATION IS=1 GAUSS NORMALIZATION @@ -652,16 +652,16 @@ FIL1 = FILMOD(L) DTE2 = DTEMOD(L+1) FIL2 = FILMOD(L+1) - WRITE(*,*) FIL1 - WRITE(*,*) FIL2 - print *, "que" +c WRITE(*,*) FIL1 +c WRITE(*,*) FIL2 +c print *, "que" C-- GET IGRF COEFFICIENTS FOR THE BOUNDARY YEARS CALL GETSHC (IU, FIL1, NMAX1, ERAD, GH1, IER) IF (IER .NE. 0) STOP - print *, "quessss" +c print *, "quessss" CALL GETSHC (IU, FIL2, NMAX2, ERAD, GH2, IER) IF (IER .NE. 0) STOP - print *, "quj" +c print *, "quj" C-- DETERMINE IGRF COEFFICIENTS FOR YEAR IF (L .LE. NUMYE-1) THEN CALL INTERSHC (YEAR, DTE1, NMAX1, GH1, DTE2, @@ -670,14 +670,14 @@ CALL EXTRASHC (YEAR, DTE1, NMAX1, GH1, NMAX2, 1 GH2, NMAX, GHA) ENDIF - print *, "quw" +c print *, "quw" 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) + DIMO = REAL(DSQRT(F0)) GH1(1) = 0.0 I=2 @@ -693,13 +693,13 @@ IF(IS.EQ.0) F0 = F0 * (2.D0 * X - 1.D0) / X F = F0 * 0.5D0 IF(IS.EQ.0) F = F * SQRT2 - GH1(I) = GHA(I-1) * F0 + GH1(I) = GHA(I-1) * REAL(F0) I = I+1 DO 9 M=1,N F = F * (X + M) / (X - M + 1.D0) IF(IS.EQ.0) F = F * DSQRT((X - M + 1.D0) / (X + M)) - GH1(I) = GHA(I-1) * F - GH1(I+1) = GHA(I) * F + GH1(I) = GHA(I-1) * REAL(F) + GH1(I+1) = GHA(I) * REAL(F) I=I+2 9 CONTINUE RETURN @@ -744,13 +744,13 @@ WRITE(FOUT,667) FSPEC c 667 FORMAT('/usr/local/etc/httpd/cgi-bin/natasha/IRI/',A12) 667 FORMAT(A258) - print *," gui" +c print *," gui" OPEN (IU, FILE=FOUT, STATUS='OLD', IOSTAT=IER, ERR=999) - print *," gua" +c print *," gua" READ (IU, *, IOSTAT=IER, ERR=999) - print *," gue" +c print *," gue" READ (IU, *, IOSTAT=IER, ERR=999) NMAX, ERAD - print *," guo" +c print *," guo" C --------------------------------------------------------------- C Read the coefficient file, arranged as follows: C @@ -787,10 +787,10 @@ ENDIF 2233 CONTINUE 2211 CONTINUE - print *," guj" +c print *," guj" 999 CLOSE (IU) - print *," guw IER",IER +c print *," guw IER",IER if ( IER .eq. -1 ) IER = 0 ! gfortran 4.1.2 bug workaround... hoping not to create problems with other versions RETURN @@ -947,7 +947,7 @@ C ASTRONOMICAL UNION . C----------------------------------------------------------------- INTEGER TL1,TL2,TL3 - CHARACTER (len=258) TP1,TP2,TP3 + CHARACTER (len=*) :: TP1,TP2,TP3 INTEGER L1,L2,L3 CHARACTER *258 P1,P2,P3 COMMON/PPATH/ L1,L2,L3,P1, P2, P3