*"    等間隔データのリニアトレンドの除去 ( 回帰直線を引き去る )
*
*"    93/11/16 horinout
*"      ( SAISHE -- hasegawa )
*"    電脳 MATH1 使用（欠損値処理）
***********************************************************************
      SUBROUTINE SUBTRD(X,K,LD,L,M)
*
      IMPLICIT  REAL (A-H,O-Z)
*  [MODIFTY]
      DIMENSION  X(K,LD,M) !" 第２次元整合寸法 LD, 計算には L 個のみ使う
*
      LOGICAL    LFIRST
      DATA       LFIRST/ .TRUE. /
      SAVE       LFIRST
      REAL*4     RMISS
*
      IF(LFIRST) THEN
        LFIRST = .FALSE.
        CALL GLPGET('RMISS',RMISS)
      ENDIF
*
      DO 10 J=1,M
      DO 10 I=1,K
        CALL SAISHE(X(I,1,J),L,K,A,B,RMISS)
*
        DO 30 N=1,L
          IF( X(K,N,J).NE.RMISS ) THEN
            XTRD = A*N + B
            X(K,N,J) = X(K,N,J) - XTRD    !" トレンドを引く
          ENDIF
   30   CONTINUE
   10 CONTINUE
*
      RETURN
      END
C######################################################
      SUBROUTINE SAISHE(Y,ND,KI,A,B,RMISS)
C######################################################

C     saishou 2jouhou de liner trend wo motomeru
C     Y :data
C     ND:data number, KI:data interval
C     A ------- katamuki
C     B ------- seppen  

C     originary made by Hasegawa,F. reformed by horinout

      IMPLICIT  REAL (A-H,O-Z)
      DIMENSION Y(*)
      REAL*4    RMISS
      SX1=0
      SX2=0
      SY1=0
      SXY=0
      N=0
      
      DO 10 I=1,ND*KI,KI
        IF (Y(I).NE.RMISS) THEN
          X = I
          SX1=SX1+X
          SX2=SX2+X**2
          SY1=SY1+Y(I)
          SXY=SXY+X*Y(I)
          N=N+1
          ENDIF
   10 CONTINUE

      IF (N.LE.1) THEN
        A=RMISS
        B=RMISS
      ELSE 
        A=( N*SXY-SX1*SY1 ) / ( N*SX2-SX1**2 )
        B=( SX2*SY1-SX1*SXY ) / ( N*SX2-SX1**2 )
      ENDIF

      RETURN
      END
