C***************************************************************
C*****  CALIBRATION FILTER SUBROUTINE LIBRARY BY A.LANGE *******
C*****  COPYRIGHT: THE FINNISH METEOROLOGICAL INSTITUTE 1982****
C***************************************************************
$LARGE: bjk,dbjk,cbjjk,mik,mjk      ! possibly over 64KB
$LARGE: GXJLK, CXGJLK, GYLK, XYJK   ! possibly over 64KB
c$LARGE: yik, xijk, gilk  ! needed only if RAO's MINQUE
      SUBROUTINE CALPRE(K,IMX,JMX,KMX,LMX,MI,MJ,MK,ML,
     1           MIK,MJK,CBJJK,CXGJLK,dbjk,CCALL,
     2           GL,GRYL,dsdkjj,
     3           GXJLK,GYLK,XYJK)
C    INPUT DATA ARRAYS
      DIMENSION MI(IMX),MJ(JMX),MK(KMX),ML(LMX)
C    OUTPUT DATA ARRAYS
      DIMENSION CBJJK(JMX,JMX,KMX),
     1          CXGJLK(JMX,LMX,KMX),dbjk(JMX,KMX),CCALL(LMX,LMX),
     2          GL(LMX),GRYL(LMX),dsdkjj(JMX,JMX),
     3          GXJLK(JMX,LMX,KMX),GYLK(LMX,KMX),XYJK(JMX,KMX),
     4          MIK(IMX,KMX),MJK(JMX,KMX)
C
C
      logical*2 MI,MJ,MK,ML,MIK,MJK
C
      double precision dsdkjj        ! 4 Nov 1997
      DOUBLE PRECISION CCALL
C
    1 DO 5 J=1,JMX
      DO 5 J2=1,JMX
      dsdkjj(J,J2)=0.0
    5 CONTINUE
C
    6 DO 8 L=1,LMX
      GL(L)=0.0
      GRYL(L)=0.0
      DO 8 L2=1,LMX
      CCALL(L,L2)=0.0
    8 CONTINUE
C
        DO 11 L=1,LMX
         GYLK(L,K)=0.0
   11   CONTINUE
        DO 20 J=1,JMX
          dbjk(J,K)=0.0
          XYJK(J,K)=0.0
          DO 12 J2=1,JMX
            CBJJK(J,J2,K)=0.0
   12     CONTINUE
          DO 15 L=1,LMX
            GXJLK(J,L,K)=0.0
            CXGJLK(J,L,K)=0.0
   15     CONTINUE
   20   CONTINUE
C
C  MAHDOTTOMAT KOMBINAATIOT POISTETAAN
C
        DO 991 I=1,IMX
          MIK(I,K)=.FALSE.
          IF(MI(I).OR.MK(K)) MIK(I,K)=.TRUE.
  991   CONTINUE
        DO 992 J=1,JMX
          MJK(J,K)=.FALSE.
          IF(MJ(J).OR.MK(K)) MJK(J,K)=.TRUE.
  992   CONTINUE
C
      RETURN
      END
C***************************************************************
C*****  CALIBRATION FILTER SUBROUTINE LIBRARY BY A.LANGE *******
C*****  COPYRIGHT: THE FINNISH METEOROLOGICAL INSTITUTE 1982****
C***************************************************************
$LARGE: bjk,dbjk,cbjjk,mik,mjk      ! possibly over 64KB
$LARGE: GXJLK, CXGJLK, GYLK, XYJK   ! possibly over 64KB
c$LARGE: yik, xijk, gilk  ! needed only if RAO's MINQUE
      SUBROUTINE CALINI(IMX,JMX,KMX,LMX,MI,MJ,MK,ML,
     1           MIK,MJK,CBJJK,CXGJLK,dbjk,CCALL,
     2           GL,GRYL,dsdkjj,
     3           GXJLK,GYLK,XYJK)
C    INPUT DATA ARRAYS
      DIMENSION MI(IMX),MJ(JMX),MK(KMX),ML(LMX)
C    OUTPUT DATA ARRAYS
      DIMENSION CBJJK(JMX,JMX,KMX),
     1          CXGJLK(JMX,LMX,KMX),dbjk(JMX,KMX),CCALL(LMX,LMX),
     2          GL(LMX),GRYL(LMX),dsdkjj(JMX,JMX),
     3          GXJLK(JMX,LMX,KMX),GYLK(LMX,KMX),XYJK(JMX,KMX),
     4          MIK(IMX,KMX),MJK(JMX,KMX)
C
C
      logical*2 MI,MJ,MK,ML,MIK,MJK
C
      double precision dsdkjj     ! 4 Nov 1977
      DOUBLE PRECISION CCALL
C
    1 DO 5 J=1,JMX
      DO 5 J2=1,JMX
      dsdkjj(J,J2)=0.0
    5 CONTINUE
C
    6 DO 8 L=1,LMX
      GL(L)=0.0
      GRYL(L)=0.0
      DO 8 L2=1,LMX
      CCALL(L,L2)=0.0
    8 CONTINUE
C
      DO 100 K=1,KMX
        DO 11 L=1,LMX
         GYLK(L,K)=0.0
   11   CONTINUE
        DO 20 J=1,JMX
          dbjk(J,K)=0.0
          XYJK(J,K)=0.0
          DO 12 J2=1,JMX
            CBJJK(J,J2,K)=0.0
   12     CONTINUE
          DO 15 L=1,LMX
            GXJLK(J,L,K)=0.0
            CXGJLK(J,L,K)=0.0
   15     CONTINUE
   20   CONTINUE
C
C  MAHDOTTOMAT KOMBINAATIOT POISTETAAN
C
        DO 991 I=1,IMX
          MIK(I,K)=.FALSE.
          IF(MI(I).OR.MK(K)) MIK(I,K)=.TRUE.
  991   CONTINUE
        DO 992 J=1,JMX
          MJK(J,K)=.FALSE.
          IF(MJ(J).OR.MK(K)) MJK(J,K)=.TRUE.
  992   CONTINUE
  100 CONTINUE
C
      RETURN
      END
C***************************************************************
C*****  CALIBRATION FILTER SUBROUTINE LIBRARY BY A.LANGE *******
C*****  COPYRIGHT: THE FINNISH METEOROLOGICAL INSTITUTE 1982****
C***************************************************************
C*      ACCumulation of cross products for CALibration
$LARGE: bjk,dbjk,cbjjk,mik,mjk      ! possibly over 64KB
$LARGE: GXJLK, CXGJLK, GYLK, XYJK   ! possibly over 64KB
c$LARGE: yik, xijk, gilk   ! needed only if RAO's MINQUE
      SUBROUTINE CALACC(K,IMX,JMX,KMX,LMX,MI,MJ,MK,ML,
     1           YKI,XKIJ,GKIL,MIK,MJK,XXJJK,ccall,
     2           GL,GXJLK,GYLK,XYJK)
C    INPUT DATA ARRAYS
      DIMENSION MI(IMX),MJ(JMX),MK(KMX),ML(LMX),
     1          YKI(IMX),XKIJ(IMX,JMX),GKIL(IMX,LMX),
     2          MIK(IMX,KMX),MJK(JMX,KMX)
C    OUTPUT DATA ARRAYS
      DIMENSION XXJJK(JMX,JMX,KMX),ccall(LMX,LMX),GL(LMX),
     1          GXJLK(JMX,LMX,KMX),GYLK(LMX,KMX),XYJK(JMX,KMX)
C
C
      logical*2 MI,MJ,MK,ML,MIK,MJK
C
      DOUBLE PRECISION ccall
C
C     Compute X'*Y and X'*X 
      DO 35 I=1,IMX
      IF(MIK(I,K)) GO TO 35
      DO 30 J=1,JMX
      IF(MJK(J,K)) GO TO 30
      XYJK(J,K)=XYJK(J,K) + XKIJ(I,J)*YKI(I)
      DO 25 J2=1,JMX
      IF(MJK(J2,K)) GO TO 25
      XXJJK(J,J2,K)=XXJJK(J,J2,K) + XKIJ(I,J)*XKIJ(I,J2)
   25 CONTINUE
   30 CONTINUE
   35 CONTINUE
C
C     Compute  GL=norms, GYLK=G'*Y and  CCALL=G'*G 
      DO 57 I=1,IMX
      IF(MIK(I,K)) GO TO 57
      DO 55 L=1,LMX
      IF(ML(L)) GO TO 55
      GL(L)=GL(L) + GKIL(I,L)**2  ! norms of the calibration predictors
      GYLK(L,K)=GYLK(L,K) + GKIL(I,L)*YKI(I)
      DO 50 L2=1,LMX
      IF(ML(L2)) GO TO 50
      ccall(L,L2)=ccall(L,L2) + GKIL(I,L)*GKIL(I,L2)
   50 CONTINUE
   55 CONTINUE
   57 CONTINUE
C
C     Compute GXJLK=G'*X
   58 DO 70 J=1,JMX
      IF(MJK(J,K)) GO TO 70
      DO 65 L=1,LMX
      IF(ML(L)) GO TO 65
      DO 60 I=1,IMX
      IF(MIK(I,K)) GO TO 60
      GXJLK(J,L,K)=GXJLK(J,L,K) + XKIJ(I,J)*GKIL(I,L)
   60 CONTINUE
   65 CONTINUE
   70 CONTINUE
C
      RETURN
      END
C***************************************************************
C*****  CALIBRATION FILTER SUBROUTINE LIBRARY BY A.LANGE *******
C*****  COPYRIGHT: THE FINNISH METEOROLOGICAL INSTITUTE 1982****
C***************************************************************
$LARGE: bjk,dbjk,cbjjk,mik,mjk      ! possibly over 64KB
$LARGE: GXJLK, CXGJLK, GYLK, XYJK   ! possibly over 64KB
c$LARGE: yik, xijk, gilk   ! needed only if RAO's MINQUE
      SUBROUTINE CALSTR(K,IMX,JMX,KMX,LMX,MI,MJ,MK,ML,
     1           TOLDET,JSTEP,JMIN,
     2           MIK,MJK,CBJJK,CXGJLK,dbjk,CCALL,GRYL,
     3           dsdkjj,GXJLK,GYLK,XYJK)
C
C    INPUT/OUTPUT DATA ARRAYS
C
      DIMENSION MI(IMX),MJ(JMX),MK(KMX),ML(LMX),
     1          MIK(IMX,KMX),MJK(JMX,KMX)
C
      DIMENSION CBJJK(JMX,JMX,KMX),
     1          CXGJLK(JMX,LMX,KMX),dbjk(JMX,KMX),CCALL(LMX,LMX),
     2          GRYL(LMX),dsdkjj(JMX,JMX),
     3          GXJLK(JMX,LMX,KMX),GYLK(LMX,KMX),XYJK(JMX,KMX)
C
C
      logical*2 MI,MJ,MK,ML,MIK,MJK
C
      double precision dsdkjj,toldet   
      DOUBLE PRECISION CCALL,det
C
    1 JDEL=0
C
    2 DO 15 J=1,JMX
      DO 10 J2=1,JMX
      dsdkjj(J,J2)=CBJJK(J,J2,K)
   10 CONTINUE
   15 CONTINUE
C
      call calinv(dsdkjj,jmx,jmx,det,*20)
ccccc      CALL GRAMIN(XXKJJ,JMX,JMX,DET,*20)
C
C*      LASKIMME MATRIISIN INVERSIN JA DETERMINANTIN          *
        IF (DET.EQ.1.0) GO TO 30
        IF (DET.GT.TOLDET) GO TO 50
C
C       OSITTEEN KOVARIANSSIMATRIISI ON SINGULAARINEN,
   20  WRITE(6,6600) K,DET
     1  ,(CBJJK(J,J,K),J=1,JMX)
 6600   FORMAT(' STRM:',I5,', DETERM:',F60.50/(' XXJJ:',10F25.20))
C       JOTEN HYLKAAMME PREDIKTOREITA HANTAPAASTA ASKELIN JSTEP
   21   JDEL=JDEL+JSTEP
        JBEG=JMX-JDEL+1
        IF(JBEG.LE.0) GO TO 30
C       POISTAMME VASTAAVAT TULOSUMMAT LASKUISTA (XY,XX JA XG):
        DO 25 JHYL=JBEG,JMX
   22   MJK(JHYL,K)=.TRUE.
        XYJK(JHYL,K)=0.0
        DO 23 J=1,JMX
        CBJJK(JHYL,J,K)=0.0
        CBJJK(J,JHYL,K)=0.0
   23   CONTINUE
        DO 24 L=1,LMX
        GXJLK(JHYL,L,K)=0.0
   24   CONTINUE
   25   CONTINUE
C
C       EMME YRITA VAHEMMALLA KUIN JMIN:LLA PREDIKTORILLA
        IF(JHYL.GT.JMIN) GO TO 2
C
C       KAIKKI PREDIKTORIT TULIVAT HYLATYIKSI
   30   MK(K)=.TRUE.
        DO 991 I=1,IMX
  991   MIK(I,K)=.TRUE.
        DO 992 J=1,JMX
        dbjk(J,K)=0.0
        MJK(J,K)=.TRUE.
  992   CONTINUE
        RETURN
C
C
C     OSITTEEN NORMAALI KASITTELY JATKUU
C
C     Compute CXGJLK=inv(XXKJJ)*GXJLK
   50 DO 85 J=1,JMX
      DO 80 J2=1,JMX
      DO 75 L=1,LMX
      CXGJLK(J,L,K)=CXGJLK(J,L,K)+dsdkjj(J,J2)*GXJLK(J2,L,K)
   75 CONTINUE
      CBJJK(J,J2,K)=dsdkjj(J,J2)
   80 CONTINUE
   85 CONTINUE
C
C
C     Compute dbjk=CBJJK*XYJK  
      DO 95 J=1,JMX
      dbjk(J,K)=0.0
      DO 90 J2=1,JMX
      dbjk(J,K)=dbjk(J,K)+dsdkjj(J,J2)*XYJK(J2,K)
   90 CONTINUE
   95 CONTINUE
c      WRITE(6,6611) K,DET,(dbjk(J,K),J=1,JMX)
c 6611 FORMAT(' STRM:',I5,', DETERM:',F99.30/(' BIJK:',20F6.1))
C
C      
C     Compute CCALL=+G'G-G'X*CX'G     
C     and      GRYL=+G'X*CX'Y-G'Y   summing over index k 
C      
      DO 130 L=1,LMX
      DO 110 L2=1,LMX
      DO 100 J=1,JMX
      CCALL(L,L2)=CCALL(L,L2)-GXJLK(J,L,K)*CXGJLK(J,L2,K)
  100 CONTINUE
  110 CONTINUE
      DO 120 J=1,JMX
      GRYL(L)=GRYL(L)+GXJLK(J,L,K)*dbjk(J,K)
  120 CONTINUE
      GRYL(L)=GRYL(L)-GYLK(L,K)
  130 CONTINUE
      RETURN
      END
C***************************************************************
C*****  CALIBRATION FILTER SUBROUTINE LIBRARY BY A.LANGE *******
C*****  COPYRIGHT: THE FINNISH METEOROLOGICAL INSTITUTE 1982****
C***************************************************************
c  kaantaa gram-matriisin paikallaan nollakonstilla 89-10-24/a.lange
      subroutine calinv(a,n,mx,det,*)
      dimension a(mx,mx)
      double precision a,det
      det=1.0
      do 60 j=1,n
      pvt=a(j,j)
      if(pvt.eq.0.0) go to 60
      det=det*pvt
c pienin luku (single pr) itseisarvoltaan msfortranissa on 1.1754944e-38
c ja vastaavasti (double pr) on 2.225073858507201d-308
      if(abs(det).lt.1.1754944d-38) return 1
      a(j,j)=1.0
      do 20 k=1,n
   20 a(j,k)=a(j,k)/pvt
      do 50 k=1,n
      if(k-j) 30,50,30
   30 t=a(k,j)
      if(t.eq.0.0) go to 50
      a(k,j)=0.0
      do 40 l=1,n
   40 a(k,l)=a(k,l)-a(j,l)*t
   50 continue
   60 continue
      return
      end
C***************************************************************
C*****  CALIBRATION FILTER SUBROUTINE LIBRARY BY A.LANGE *******
C*****  COPYRIGHT: THE FINNISH METEOROLOGICAL INSTITUTE 1982****
C***************************************************************
$LARGE: bjk,dbjk,cbjjk,mik,mjk      ! possibly over 64KB
$LARGE: GXJLK, CXGJLK, GYLK, XYJK   ! possibly over 64KB
c$LARGE: yik, xijk, gilk  ! needed only if RAO's MINQUE
      SUBROUTINE CALFKF(IMX,JMX,KMX,LMX,MI,MJ,MK,ML,
     1                  TOLXJK,TOLGL,TOLDET,JSTEP,JMIN,
     2                  MIK,MJK,
     3                  dbjk,CBJJK,dcal,CCALL,CXGJLK,
     4                  GL,GRYL,DSDKJJ,
     5                  GXJLK,GYLK,XYJK)
C
C    INPUT/OUTPUT DATA ARRAYS
      DIMENSION MI(IMX),MJ(JMX),MK(KMX),ML(LMX)
C
      DIMENSION MIK(IMX,KMX),MJK(JMX,KMX),
     1          dbjk(JMX,KMX),CBJJK(JMX,JMX,KMX),dcal(LMX),
     2          CCALL(LMX,LMX),CXGJLK(JMX,LMX,KMX)
C
C    WORK ARRAYS
      DIMENSION GL(LMX),GRYL(LMX),DSDKJJ(JMX,JMX),
     2 GXJLK(JMX,LMX,KMX),GYLK(LMX,KMX),XYJK(JMX,KMX)
C
C
      logical*2 MI,MJ,MK,ML,MIK,MJK
      double precision dsdkjj,toldet
      DOUBLE PRECISION CCALL,DETERM
C                      CCALL=SUMMAK(GK'GK-GK'XK*INV(X'X)*XK'GK)
C
C
CC    TULOSUMMAMATRIISIT ON KUMULOITU CALACC-RUTIINISSA
CC
C
C       G-PREDIKTORIN NELIO VAHINT. 0.0000001 SINGULAARISUUDEN VALTT.
    1 IF(TOLGL.LE.0.0) TOLGL=0.0000001
c      WRITE(6,6001) TOLGL
c 6001 FORMAT(' TOLGL =',F20.10)
C
C       X-PREDIKTORIN NELIO VAHINT. 0.0000001 SINGULAARISUUDEN VALTT.
    2 IF(TOLXJK.LE.0.0) TOLXJK=0.0000001
c      WRITE(6,6002) TOLXJK
c 6002 FORMAT(' TOLXJK =',F20.10)
C
C        POISTETAAN SINGULAARISET G-PREDIKTORIT (PITUUS)
      DO 3 L=1,LMX
      IF(ML(L)) GO TO 3
      IF(GL(L).GE.TOLGL) GO TO 3
      ML(L)=.TRUE.
      WRITE(6,6003) L,GL(L),TOLGL
 6003 FORMAT(' G-PREDICTOR WITHDRAWN L,GL,TOL:',I5,2F20.10)
    3 CONTINUE
C
C       POISTETAAN SINGULAARISET X-PREDICTORIT (PITUUS)
C
      DO 15 K=1,KMX
      IF(MK(K)) GO TO 15
      DO 14 J=1,JMX
      IF(MJ(J)) GO TO 14
      IF(CBJJK(J,J,K).GE.TOLXJK) GO TO 14
      MJK(J,K)=.TRUE.
      IF(K.EQ.100*(K/100)) WRITE(6,6004) J,K,CBJJK(J,J,K),TOLXJK
 6004 FORMAT(' XJK-PREDICTOR WITHDRAWN J,K,XXJK,TOL:',2I5,2F20.10)
   14 CONTINUE
   15 CONTINUE
C
C       TASSA ON K-LUUPPI YLI KAIKKIEN OSITTEIDEN
C
C  ALKUNOLLAUKSET HOIDETTIIN CALINI-RUTIINISSA ENNEN CALACC:UA
C
C  CALSTR-RUTIINI KAY LAPI OSITTEET JOHTAEN SISAISET REGRESSIOT
C
   17 DO 20 K=1,KMX
      IF(MK(K)) GO TO 20
      CALL CALSTR(K,IMX,JMX,KMX,LMX,MI,MJ,MK,ML,
     1 TOLDET,JSTEP,JMIN,
     2 MIK,MJK,CBJJK,CXGJLK,dbjk,CCALL,GRYL,
     3 DSDKJJ,GXJLK,GYLK,XYJK)
C
   20 CONTINUE
C
C
C**** KAANNETAAN MATRIISI CCALL PAIKALLAAN **************
C
cc      DO 25 L=1,LMX
cc   25 WRITE(6,6131) L,(CCALL(L,L2),L2=1,LMX)
cc 6131 FORMAT(' CCALL',I4,6F10.6/(10X,5F10.6))
C
C
      CALL calinv(CCALL,LMX,LMX,DETERM,*26)
C
       if(determ.gt.toldet) go to 30
   26  WRITE(6,6025) DETERM
 6025 FORMAT(' Singular, sorry, DETERM  =',D25.16)
C
c      IF(DETERM.ge.TOLDET) GO TO 31
C
C
C     SINGULAARINEN RATKAISU COMBINED REGRESSION SUHTEEN,
C     JOTEN POISTAMME COMBINED REGRESSION PREDIKTORIT KOKONAAN
C
      DO 28 L=1,LMX
      DO 27 L2=1,LMX
      CCALL(L,L2)=0.0
   27 CONTINUE
   28 CONTINUE
C
   30 continue
C
c   31 DO 32 L=1,LMX  
c   32 WRITE(6,6132) L,(CCALL(L,L2),L2=1,LMX)
c 6132 FORMAT(' COVAR',I4,6F10.6/(10X,5F10.6))
c      CALL calinv(CCALL,LMX,LMX,DETERM,*3333)
c 3333 WRITE(6,6025) DETERM
c      DO 33 L=1,LMX
c   33 WRITE(6,6131) L,(CCALL(L,L2),L2=1,LMX)
c      CALL calinv(CCALL,LMX,LMX,DETERM,*4444)
c 4444 WRITE(6,6025) DETERM
c      DO 34 L=1,LMX
c   34 WRITE(6,6132) L,(CCALL(L,L2),L2=1,LMX)
C
C
C
C     CCALL-MATRIISISSA ON NYT SYST.VIRHEIDEN KOVARIANSSIT JA
C     DIAGONAALILLA VARIANSSIT (JOS EI KAANTYNYT: NO COMBINED EFFECTS)
C
C
C***********************************************************
C
c      WRITE(6,6040) (GRYL(L),L=1,LMX)
c 6040 FORMAT(' GRYL     ',12F10.1/(10X,12F10.1))
C
C       LASKETAAN CALIBROINTIVIRHEET dcal=-CCALL*GRYL
      DO 50 L=1,LMX
      dcal(L)=0.0
      DO 40 L2=1,LMX
      dcal(L)=dcal(L)-CCALL(L,L2)*GRYL(L2)
   40 CONTINUE
   50 CONTINUE
cc      WRITE(6,6050) (dcal(L),L=1,LMX)
cc 6050 FORMAT(' Iteration steps for calibration parameters:',/,
cc     1         f5.0,5f5.1,10F5.0)
C
C
C       FILTTEROIDUT REGR. PARAMETRIT KULLEKIN OSITTEELLE K
C
      DO 100 K=1,KMX
          DO 90 L=1,LMX
               DO 80 J=1,JMX
               dbjk(J,K)=dbjk(J,K)-CXGJLK(J,L,K)*dcal(L)
   80          CONTINUE
   90     CONTINUE
  100 CONTINUE
      DO 150 K=1,KMX
CCCC  WRITE(6,6150) K,(CBJJK(J,J,K),J=1,JMX)
C6150 FORMAT(I4,10F12.7)
          DO 140 J=1,JMX
          DO 130 J2=1,JMX
          DSDKJJ(J,J2)=0.0
              DO 120 L=1,LMX
              DO 110 L2=1,LMX
              DSDKJJ(J,J2)=DSDKJJ(J,J2)+CXGJLK(J,L,K)
     1                                 *CCALL(L,L2)
     2                                 *CXGJLK(J2,L2,K)
  110         CONTINUE
  120         CONTINUE
          CBJJK(J,J2,K)=CBJJK(J,J2,K)+DSDKJJ(J,J2)
  130     CONTINUE
  140     CONTINUE
CCCC  WRITE(6,6150) K,(CBJJK(J,J,K),J=1,JMX)
  150 CONTINUE
  999 RETURN
      END