svd.for

svd气象中用的fortran程序

     *+ABS(SX(I+3))+ABS(SX(I+4))+ABS(SX(I+5))
   50 CONTINUE
   60 SASUM=STEMP
      RETURN
      END
CCCCC SAXPY
      SUBROUTINE SAXPY(N,SA,SX,INCX,SY,INCY)
      REAL SX(1),SY(1),SA
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
      IF(N.LE.0) RETURN
      IF(SA.EQ.0.0) RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1) GO TO 20
      IX=1
      IY=1
      IF(INCX.LT.0) IX=(-N+1)*INCX+1
      IF(INCY.LT.0) IY=(-N+1)*INCY+1
      DO 10 I=1,N
      SY(IY)=SY(IY)+SA*SX(IX)
      IX=IX+INCX
      IY=IY+INCY
   10 CONTINUE
      RETURN
   20 M=MOD(N,4)
      IF(M.EQ.0) GO TO 40
      DO 30 I=1,M
      SY(I)=SY(I)+SA*SX(I)
   30 CONTINUE
      IF(N.LT.4) RETURN
   40 MP1=M+1
      DO 50 I=MP1,N,4
      SY(I)=SY(I)+SA*SX(I)
      SY(I+1)=SY(I+1)+SA*SX(I+1)
      SY(I+2)=SY(I+2)+SA*SX(I+2)
      SY(I+3)=SY(I+3)+SA*SX(I+3)
   50 CONTINUE
      RETURN
      END
CCCCC SCOPY
      SUBROUTINE SCOPY(N,SX,INCX,SY,INCY)
      REAL SX(1),SY(1)
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
      IF(N.LE.0) RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1)GO TO 20
      IX=1
      IY=1
      IF(INCX.LT.0) IX=(-N+1)*INCX+1
      IF(INCY.LT.0) IY=(-N+1)*INCY+1
      DO 10 I=1,N
      SY(IY)=SX(IX)
      IX=IX+INCX
      IY=IY+INCY
   10 CONTINUE
      RETURN
   20 M=MOD(N,7)
      IF(M.EQ.0) GO TO 40
      DO 30 I=1,M
      SY(I)=SX(I)
   30 CONTINUE
      IF(N.LT.7) RETURN
   40 MP1=M+1
      DO 50 I=MP1,N,7
      SY(I)=SX(I)
      SY(I+1)=SX(I+1)
      SY(I+2)=SX(I+2)
      SY(I+3)=SX(I+3)
      SY(I+4)=SX(I+4)
      SY(I+5)=SX(I+5)
      SY(I+6)=SX(I+6)
   50 CONTINUE
      RETURN
      END
CCCCC SDOT
      REAL FUNCTION SDOT(N,SX,INCX,SY,INCY)
      REAL SX(1),SY(1),STEMP
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
      STEMP=0.0E0
      SDOT=0.0E0
      IF(N.LE.0) RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1) GO TO 20
      IX=1
      IY=1
      IF(INCX.LT.0) IX=(-N+1)*INCX+1
      IF(INCY.LT.0) IY=(-N+1)*INCY+1
      DO 10 I=1,N
      STEMP=STEMP+SX(IX)*SY(IY)
      IX=IX+INCX
      IY=IY+INCY
   10 CONTINUE
      SDOT=STEMP
      RETURN
   20 M=MOD(N,5)
      IF(M.EQ.0) GO TO 40
      DO 30 I=1,M
      STEMP=STEMP+SX(I)*SY(I)
   30 CONTINUE
      IF(N.LT.5) GO TO 60
   40 MP1=M+1
      DO 50 I=MP1,N,5
      STEMP=STEMP+SX(I)*SY(I)+SX(I+1)*SY(I+1)+
     *SX(I+2)*SY(I+2)+SX(I+3)*SY(I+3)+SX(I+4)*SY(I+4)
   50 CONTINUE
   60 SDOT=STEMP
      RETURN
      END
CCCCC SNRM2
      REAL FUNCTION SNRM2(N,SX,INCX)
      INTEGER NEXT
      REAL SX(1),CUTLO,CUTHI,HITEST,SUM,XMAX,ZERO,ONE
      DATA ZERO,ONE/0.0E0,1.0E0/
      DATA CUTLO,CUTHI/4.441E-16,1.304E19/
      IF(N.GT.0) GO TO 10
      SNRM2=ZERO
      GO TO 300
   10 ASSIGN 30 TO NEXT
      SUM=ZERO
      NN=N*INCX
      I=1
   20 GO TO NEXT,(30,50,70,110)
   30 IF(ABS(SX(I)).GT.CUTLO) GO TO 85
      ASSIGN 50 TO NEXT
      XMAX=ZERO
   50 IF(SX(I).EQ.ZERO) GO TO 200
      IF(ABS(SX(I)).GT.CUTLO) GO TO 85
      ASSIGN 70 TO NEXT
      GO TO 105
  100 I=J
      ASSIGN 110 TO NEXT
      SUM=(SUM/SX(I))/SX(I)
  105 XMAX=ABS(SX(I))
      GO TO 115
   70 IF(ABS(SX(I)).GT.CUTLO) GO TO 75
  110 IF(ABS(SX(I)).LE.XMAX) GO TO 115
      SUM=ONE+SUM*(XMAX/SX(I))**2
      XMAX=ABS(SX(I))
      GO TO 200
  115 SUM=SUM+(SX(I)/XMAX)**2
      GO TO 200
   75 SUM=(SUM*XMAX)*XMAX
   85 HITEST=CUTHI/FLOAT(N)
      DO 95 J=I,NN,INCX
      IF(ABS(SX(J)).GE.HITEST) GO TO 100
   95 SUM=SUM+SX(J)**2
      SNRM2=SQRT(SUM)
      GO TO 300
  200 CONTINUE
      I=I+INCX
      IF(I.LE.NN) GO TO 20
      SNRM2=XMAX*SQRT(SUM)
  300 CONTINUE
      RETURN
      END
CCCCC SROT
      SUBROUTINE SROT(N,SX,INCX,SY,INCY,C,S)
      REAL SX(1),SY(1),STEMP,C,S
      INTEGER I,INCX,INCY,IX,IY,N
      IF(N.LE.0) RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1) GO TO 20
      IX=1
      IY=1
      IF(INCX.LT.0) IX=(-N+1)*INCX+1
      IF(INCY.LT.0) IY=(-N+1)*INCY+1
      DO 10 I=1,N
      STEMP=C*SX(IX)+S*SY(IY)
      SY(IY)=C*SY(IY)-S*SX(IX)
      SX(IX)=STEMP
      IX=IX+INCX
      IY=IY+INCY
   10 CONTINUE
      RETURN
   20 DO 30 I=1,N
      STEMP=C*SX(I)+S*SY(I)
      SY(I)=C*SY(I)-S*SX(I)
      SX(I)=STEMP
   30 CONTINUE
      RETURN
      END
CCCCC SROTG
      SUBROUTINE SROTG(SA,SB,C,S)
      REAL SA,SB,C,S,ROE,SCALE,R,Z
      ROE=SB
      IF(ABS(SA).GT.ABS(SB)) ROE=SA
      SCALE=ABS(SA)+ABS(SB)
      IF(SCALE.NE.0.0) GO TO 10
      C=1.0
      S=0.0
      R=0.0
      GO TO 20
   10 R=SCALE*SQRT((SA/SCALE)**2+(SB/SCALE)**2)
      R=SIGN(1.0,ROE)*R
      C=SA/R
      S=SB/R
   20 Z=1.0
      IF(ABS(SA).GT.ABS(SB)) Z=S
      IF(ABS(SB).GE.ABS(SA).AND.C.NE.0.0)  Z=1.0/C
      SA=R
      SB=Z
      RETURN
      END
CCCCC SSCAL
      SUBROUTINE SSCAL(N,SA,SX,INCX)
      REAL SA,SX(1)
      INTEGER I,INCX,M,MP1,N,NINCX
      IF(N.LE.0) RETURN
      IF(INCX.EQ.1) GO TO 20
      NINCX=N*INCX
      DO 10 I=1,NINCX,INCX
      SX(I)=SA*SX(I)
   10 CONTINUE
      RETURN
   20 M=MOD(N,5)
      IF(M.EQ.0) GO TO 40
      DO 30 I=1,M
      SX(I)=SA*SX(I)
   30 CONTINUE
      IF(N.LT.5) RETURN
   40 MP1=M+1
      DO 50 I=MP1,N,5
      SX(I)=SA*SX(I)
      SX(I+1)=SA*SX(I+1)
      SX(I+2)=SA*SX(I+2)
      SX(I+3)=SA*SX(I+3)
      SX(I+4)=SA*SX(I+4)
   50 CONTINUE
      RETURN
      END
CCCCC SSWAP
      SUBROUTINE SSWAP(N,SX,INCX,SY,INCY)
      REAL SX(1),SY(1),STEMP
      INTEGER I,INCX,INCY,IX,IY,M,MP1,N
      IF(N.LE.0) RETURN
      IF(INCX.EQ.1.AND.INCY.EQ.1) GO TO 20
      IX=1
      IY=1
      IF(INCX.LT.0) IX=(-N+1)*INCX+1
      IF(INCY.LT.0) IY=(-N+1)*INCY+1
      DO 10 I=1,N
      STEMP=SX(IX)
      SX(IX)=SY(IY)
      SY(IY)=STEMP
      IX=IX+INCX
      IY=IY+INCY
   10 CONTINUE
      RETURN
   20 M=MOD(N,3)
      IF(M.EQ.0) GO TO 40
      DO 30 I=1,M
      STEMP=SX(I)
      SX(I)=SY(I)
      SY(I)=STEMP
   30 CONTINUE
      IF(N.LT.3) RETURN
   40 MP1=M+1
      DO 50 I=MP1,N,3
      STEMP=SX(I)
      SX(I)=SY(I)
      SY(I)=STEMP
      STEMP=SX(I+1)
      SX(I+1)=SY(I+1)
      SY(I+1)=STEMP
      STEMP=SX(I+2)
      SX(I+2)=SY(I+2)
      SY(I+2)=STEMP
   50 CONTINUE
      RETURN
      END
CCCCC TCFC
      SUBROUTINE TCFC(UV,M,TC,X,AM,N)
c                TCFC(U,N,TC1,SLP,AM1,IT)
      DIMENSION UV(M,M),TC(M,N),AM(M),X(M,N)
      DO 11 I=1,M
      DO 11 J=1,N
      C=0.0
      DO 12 K=1,M
   12 C=C+UV(K,I)*X(K,J)
   11 TC(I,J)=C
C     WRITE(6,'(5X,16HTIME COEFFICIENT)')
C     WRITE(6,33) ((TC(I,J),I=1,16),J=1,N)
C     WRITE(MM,33) ((TC(I,J),I=1,16),J=1,N)
C  33 FORMAT(1X,16F8.2)
C     DO 13 I=1,M
C     DO 13 J=1,N
C     TCC(J,I)=TC(I,J)
C  13 CONTINUE
      DO 14 I=1,M
      D=0.0
      DO 15 K=1,N
   15 D=D+TC(I,K)*TC(I,K)
   14 AM(I)=D
      E=0.0
      DO 16 II=1,M
      E=E+AM(II)
   16 CONTINUE
      DO 17 II=1,M
      AM(II)=AM(II)/E
   17 CONTINUE
C     WRITE(6,68)
C  68 FORMAT(1X,27HFANG CA AND FANG CA LEI JIA)
C     WRITE(6,44) (AM(II),II=1,16)
C  44 FORMAT(1X,16F8.3)
C     DO 18 I=2,N
C  18 AM(I)=AM(I-1)+AM(I)
C     WRITE(6,55) (AM(I),I=1,16)
C  55 FORMAT(1X,16F8.3)
      RETURN
      END	  

CCCCC TCHXG~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      SUBROUTINE TCHXG(R,TC1,TC2,M,IP,N,IT)
c                  TCHXG(R,TC1,TC2,M,IT)
      DIMENSION R(M),XC(16),XC1(16),XC2(16),
     *          T(16),TC1(N,IT),TC2(IP,IT)

      DO 19 I=1,M
      T(I)=0.0
      DO 20 J=1,it
      T(I)=T(I)+TC1(I,J)
   20 CONTINUE
      T(I)=T(I)/REAL(it)
   19 CONTINUE
      DO 21 I=1,M
      DO 21 J=1,it
      TC1(I,J)=TC1(I,J)-T(I)
   21 CONTINUE
      DO 22 I=1,M
      T(I)=0.0
      DO 23 J=1,it
      T(I)=T(I)+TC2(I,J)
   23 CONTINUE
      T(I)=T(I)/REAL(it)
   22 CONTINUE
      DO 24 I=1,M
      DO 24 J=1,it
      TC2(I,J)=TC2(I,J)-T(I)
   24 CONTINUE
      DO 25 I=1,M
      XC(I)=0.0
   25 CONTINUE
      DO 26 I=1,M
      DO 27 J=1,it
      XC(I)=XC(I)+TC1(I,J)*TC2(I,J)
   27 CONTINUE
   26 CONTINUE
      DO 28 I=1,M
      XC1(I)=0.0
      XC2(I)=0.0
   28 CONTINUE
      DO 29 I=1,M
      DO 30 J=1,it
      XC1(I)=XC1(I)+TC1(I,J)*TC1(I,J)
      XC2(I)=XC2(I)+TC2(I,J)*TC2(I,J)
   30 CONTINUE
   29 CONTINUE
      DO 31 I=1,M
      R(I)=XC(I)/SQRT(XC1(I)*XC2(I))
   31 CONTINUE
      WRITE(6,'(5X,2HRC)')
      WRITE(6,66) (R(I),I=1,16)
   66 FORMAT(1X,16F8.4)
      RETURN
      END
CCCCC CACOR
      SUBROUTINE CACOR(SSTA,H5A,X,N,IP,IT)
      DIMENSION SSTA(IP,IT),H5A(N,IT),X(N,IP)
      DO 100 I=1,N
      DO 100 K=1,IP
      C=0.0
      DO 95 J=1,IT
   95 C=C+H5A(I,J)*SSTA(K,J)
      X(I,K)=C/REAL(IT)
  100 CONTINUE
      RETURN
      END
      SUBROUTINE nomal(N,X,XX)
      DIMENSION X(N),XX(N)
      PX=0.
      DO 10 I=1,N
  10  PX=PX+X(I)
      PX=PX/FLOAT(N)
      S=0.
      DO 20 I=1,N
  20  S=S+(X(I)-PX)**2
      S=SQRT(S/FLOAT(N))
      DO 30 I=1,N
  30  XX(I)=(X(I)-PX)/S
      RETURN
      END