        program int_csv
! 17.11.1998 Bei Karte Diff: einen 2-ten parameter eingefuehrt, falls =1
!   gesetzt, wird in die Ausgabedatei der interpolierte Wert gesetzt. (sinnvoll
!   nur bei erstem parameter =0 (oder groesser)
! 21.4.98 Karte Wweg: eingefuegt
! Wweg: 6,1400,1600,16,8,9
! setzt die werte auf position 16 auf Default, die zwischen 8 und 9 liegen und
! wo gleichzeitig die Variable 6 zwischen 1400 und 1600 liegt 
! 10.3.98 parameter maxvar eingefuegt und entsprechend von 30 auf 50 erhoeht.
!  bei VAR= kann jetzt oxyg+nutri auch nicht angegeben werden, werden dann auf
!  1 gesetzt (eine fehlermoglichkeit, aber was solls)
!  Die Defaultbehandlung ist noch nicht ganz astrein, am bestens ist es wohl die
!  Krte Defa: -9. oder so ganz am Anfang zu setzen.
!  Keine Schwierigkeiten sollte es mit Defaulwert -9999. geben.
! neu 1.jun.96 label Z-Bereich
! beim Einlesen (Karte muss daher vor EING: angegeben werden, werden
! nur Zyklen in dem Angegebenen Bereich der Variable uebernohmen. (aehnlich
! D-Bereich, nur wird der ganze Zyklus beim einlesen weggeschmissen, schon
! eingelesen Zyklen sind davon nicht betroffen
!
! NEU 16.2.99 label-24= eich:
        use mk4lib
        parameter (max_jcsv=30000)
        parameter (maxvar=50)
        parameter (labmax=24)
!        character*4 label(labmax)/'A-Format','EINGabe','AUSGabe','ENDE',
!     *   'INTE','D=INterpoliert','STAT','DEPT','PARA',
!     *   'Moden','Range','Default','D-Bereich','Diff','START',
!     *   'VAR=','B-multi','M-AUsgabe','Info','MASKE','Z-Bereich'/
         character*4 label(labmax)
         data label/'A-Fo','EING','AUSG','ENDE',
     *   'INTE','D=IN','STAT','DEPT','PARA',
     *   'Mode','Rang','Defa','D-Be','Diff','STAR',
     *   'VAR=','B-mu','M-AU','Info','MASK','Z-Be','NORM','Wweg','eich'/
       logical m_ausgabe/.false./
! eingabe, angabe des eingabedatei
! Ausgabe, angabe der Ausgabedatei
! ende=schluss
! int1, zu interpolierende Variablen (nur mit originaldaten zu interpolieren)
! stat, Variable der Station und +- wert
! dept, Tiefenvariable und +- Wert und %-wert
! para, die Variablen die bei der interpolation benutzt werden sollen
        CHARACTER*132 TEXT
        character Variable*4
        integer i_csv,imoden      ! anzahl Werte, bzw. Moden
        integer j_csv,jmoden      ! anzahl der Zyklen, bzw. Tiefenstufen
        integer mod_ros,ivar(maxvar)
        real r_csv(max_jcsv,maxvar)
        real r_csv_aus(max_jcsv,maxvar)
        integer*1 qualt(max_jcsv,maxvar)
        real Moden(300,0:5)
        real*8 r8,stddif,meandif
        real xint(maxvar)
C in xint stehen die Werte des zu interpoliereden Zyklus drin
C (normierte Werte falls NORM: angegeben wurde 
        real b(5000),a(5000,maxvar),x(maxvar)
        real w(maxvar),v(maxvar,maxvar)
        integer int1(maxvar),para(maxvar)
        integer Schon_interpoliert/0/
        ! 0=noch nicht interpoliert
        ! 1=schon interpoliert
        ! 2=schon interpoliert aber neue Parameter ausgewaehlt
        real def/-9999./      ! Defaultwert
	logical L_Norm

        TYPE     S_D_Bereich 
          logical da 
          real min
          real max
          real stdabw
        END TYPE    
! da=definiert, min,max grenzen stdabw=genauigkeit,standardabweichung 
        TYPE   (s_d_bereich) d_bereich(maxvar)
        TYPE   (s_d_bereich) z_bereich(maxvar)

        TYPE     s_range 
         integer var
         real r
         real pro
        END TYPE     
        parameter (max_range=8)
        TYPE   (s_range)range(max_range)
        integer irange/0/
        real rang_min(max_range),rang_max(max_range)
        character*132 AUS_format
        data AUS_format/'##'/
        integer iw/0/      ! variable der unit zum wegschreiben der einzelstatistik
        real Bereichfaktor/0./,Max_faktor/100./
        integer info/0/
        integer maske_x/0/,maske_y/0/
        logical p_maske
        integer r_maske
        external r_maske,p_maske
        integer Vdiff(5)
        data Vdiff/1,3,6,4,5/
        integer jint1/0/,jpara/0/,n_interpol/0/,j_interpol/0/
        integer n/0/,idif/0/
        data imoden/0/,jmoden/0/,j_csv/0/,i_csv/0/
! ende Deklarationenen ============================================
!        do i=1,maxvar
!         d_bereich(i)%da=.false.
!         z_bereich(i)%da=.false.
!        enddo
         d_bereich%da=.false.
         z_bereich%da=.false.
        qualt=0
        i_csv=0
        j_csv=1
        errunit=0
        mod_ros=2
        i_Ausgabe=-1

    5   DO WHILE(.TRUE.)
           READ(5,'(A132)',END=50) TEXT
           write(*,'(1x,a)') text(1:60)
           ilab=label_such(text,label,labmax)
           IF(ilab.eq.0) then        ! Kommentarzeile
           elseIF(ilab.eq.1) then        ! format
              text=adjustl(text)
              AUS_Format=text
           elseIF(ilab.eq.2) then        ! Eingabe
              if(Schon_interpoliert.gt.0) Schon_interpoliert=2
! eingabe, angabe des eingabedatei
              i=rles(text,1)
              if(i.gt.0 .and. i_csv.eq.i) then      ! Zweite Eingabedatei
                 write(*,*) ' Neue Daten werden an fruehere Daten angehaengt'
                 j_csv=j_csv+1
              else
                 if(i_csv.gt.0) then
                    write(*,*) ' Fruehere Dateien werden ueberschrieben'
                 endif
                 j_csv=1
                 i_csv=abs(i)
              endif
              if(i_csv .LE. 0) stop ' FF: CSV-Eingabevariablen <=0 ?'
              if(i_csv .gt.maxvar) stop ' FF: CSV mehr als maxvar Variablen'

              open(11,file=adjustl(cles(text,2)),status='old')
              ros_def=rles(cles(text,3),1)
              if(ros_def.eq.0) ros_def=-9999.
              do while(.true.)
  810       continue             
  547            read(11,*,err=810,end=710) (r_csv(j_csv,i),i=1,i_csv)
                 do i=1,i_csv
                    if(z_bereich(i)%da .and.
     *               (r_csv(j_csv,i).lt. z_bereich(i)%min .or.
     *               r_csv(j_csv,i).gt. z_bereich(i)%max     )) goto 547
                 enddo
                 ! Rest des Zyklus mit def-fuellen
                 do i=i_csv+1,maxvar
                    r_csv(j_csv,i)=def
                    r_csv_aus(j_csv,i)=r_csv(j_csv,i)
                 enddo
                 ! rosettenwerte ins Ausgabearraykopieren
                 do i=1,i_csv
                    r_csv_aus(j_csv,i)=r_csv(j_csv,i)
                 enddo
                 j_csv=j_csv+1
                 if(j_csv.gt.max_jcsv) stop 'Rosettenarray zu klein'
              enddo
!  810         stop 'Fehler beim lesen CSV-daten'
  710         j_csv=j_csv-1
              write(*,'('' Gesamte Anzahl CSV-Werte:'',i8)') j_csv
              close(11)
              if(.not.m_ausgabe) m_csv=j_csv
           elseIF(ilab.eq.3) then        !
              ! Ausgabe, angabe der Ausgabedatei
              close(33)
              open(33,file=adjustl(text),status='new')
           elseIF(ilab.eq.4) then        !
              i_Ausgabe=-1
              ! ende=schluss
              goto 50
           elseIF(ilab.eq.5) then        !
              if(Schon_interpoliert.gt.0) Schon_interpoliert=2
              !C int1, zu interpolierende Variablen (nur mit originaldaten zu interpolieren)
              jint1=0
              do while(.true.)
                 jint1=jint1+1
                 int1(jint1)=rles(text,jint1)
                 if(int1(jint1).le.0) goto 678
              enddo
  678         jint1=jint1-1
           elseIF(ilab.eq.6) then        ! D=INterpoliert
              if(Schon_interpoliert.gt.0) Schon_interpoliert=2
              if(i_csv.le.0) stop ' D=IN, i_csv=0'
              if(m_csv.le.0) stop ' D=IN, m_csv=0'
              jj=rles(text,1)
              do i=1,i_csv
                 if(jj.eq.1 .or. jj.eq.2) then
                    do j=1,m_csv
                       r_csv(j,i)=r_csv_aus(j,i)
                    enddo
                 elseif(jj.eq.0 .or. jj.eq.-1) then
                    do j=1,m_csv
                       if(qualt(j,i).le.0) then
                          r_csv(j,i)=r_csv_aus(j,i)
                       else
                          r_csv(j,i)=-9                       
                       endif
                    enddo
                 else
                    write(*,*) text
                    write(*,*) jj
                    stop 'Falsche Eingabe von Label D=IN:'
                 endif
              enddo
              if(jj.eq.2 .or. jj.eq.-1) qualt=0

           elseIF(ilab.eq.7) then        !
              if(Schon_interpoliert.gt.0) Schon_interpoliert=2
              !C stat, Variable der Station und +- wert
              if(istat.eq.0) then
                 irange=irange+1
                 if(irange.gt.max_range) stop ' zuviele Range-einschraenkungen'
                 istat=irange
              endif
              range(istat)%var=rles(text,1)
              range(istat)%r  =rles(text,2)
              range(istat)%pro=0.

           elseIF(ilab.eq.8) then        !
              !C dept, Tiefenvariable und +- Wert und %-wert
              if(Schon_interpoliert.gt.0) Schon_interpoliert=2
              if(idept.eq.0) then
                 irange=irange+1
                 if(irange.gt.max_range) stop ' zuviele Range-einschraenkungen'
                 idept=irange
              endif
              range(idept)%var=rles(text,1)
              range(idept)%r  =rles(text,2)
              range(idept)%pro=rles(text,3)
              if (range(idept)%pro.gt.1)
     *         range(idept)%pro=range(idept)%pro/100.

           elseIF(ilab.eq.9) then        !
              !C para, die Variablen die bei der interpolation benutzt werden sollen
              if(Schon_interpoliert.gt.0) Schon_interpoliert=2
              jpara=0
              do while(.true.)
                 jpara=jpara+1
                 para(jpara)=rles(text,jpara)
                 if(para(jpara).eq.0) goto 676
              enddo
  676         jpara=jpara-1
           elseIF(ilab.eq.10) then        !
              if(Schon_interpoliert.gt.0) Schon_interpoliert=2
              !C ' Moden-daten'
              imoden=rles(text,1)
              if(imoden.le.0) then
                 write(*,*) ' keine Moden verwendet'
                 mod_ros=2
              elseif(imoden.gt.5) then
                 stop 'zuviele Moden'
              else
                 text=cles(text,2)
                 open(11,file=text,status='old')
                 read(11,'(a)',err=811) text
                 jmoden=1
                 do while(.true.)
                    read(11,*,err=811,end=711) (moden(jmoden,i),i=0,imoden)
                    jmoden=jmoden+1
                    if(jmoden.gt.300) stop 'Modalarray zu klein'
                 enddo
  811            stop 'Fehler beim lesen Modalzerlegung'
  711            jmoden=jmoden-1
                 mod_ros=0
                 close(11)
              endif

           elseIF(ilab.eq.11) then        ! Range-angabe
              if(Schon_interpoliert.gt.0) Schon_interpoliert=2
              irange=irange+1
              if(irange.gt.max_range) stop ' zuviele Range-einschraenkungen'
              range(irange)%var=rles(text,1)
              range(irange)%r  =rles(text,2)
              range(irange)%pro=rles(text,3)
              if (range(irange)%pro.gt.1)
     *         range(irange)%pro=range(irange)%pro/100.
              if(range(irange)%var.le.0) then
                 write(*,*) ' alle Range geloescht'
                 idept=0
                 istat=0
                 irange=0
              endif
           elseIF(ilab.eq.12) then        ! Defa: Default-angabe
              if(Schon_interpoliert.gt.0) Schon_interpoliert=2
              def=rles(text,1)
              if(def.eq.0) def=-9999.
           elseIF(ilab.eq.13) then        ! D-Be: definitions Bereichangabe
              if(Schon_interpoliert.gt.0) Schon_interpoliert=2
              i=rles(text,1)
              if(abs(i).gt.maxvar) then
                 stop ' Variablenposition zu gross'
              elseif(i.gt.0) then
                 D_Bereich(i)%da =.true.
                 D_Bereich(i)%min=rles(text,2)
                 D_Bereich(i)%max=rles(text,3)
                 D_Bereich(i)%stdabw=rles(text,4)
                 if(D_Bereich(i)%stdabw.le.0) D_Bereich(i)%stdabw=1e32
                 if(D_Bereich(i)%max.lt.D_Bereich(i)%min)
     *            stop 'Max<=min bei Definitionsbereich'
                 ! Aenderung 17.4.96, bei Angabe der Karte D-Bereich werden schon
                 ! eingelesene Zyklen auf den Angegebenen Bereich ueberprueft
                 ! und gegebenenfalls auf default gesetzt
                 do j=1,j_csv
                    if(r_csv(j,i).lt.d_bereich(i)%min .or.
     *                 r_csv(j,i).gt.d_bereich(i)%max) then
                       r_csv(j,i)=def
                       if(qualt(j,i).le.0) r_csv_aus(j,i)=def
                    endif
                 enddo
              elseif(i.lt.0) then
                 D_Bereich(abs(i))%da =.false.
              else
                 do i=1,maxvar
                    D_Bereich(i)%da =.false.
                 enddo
              endif

           elseIF(ilab.eq.14) then        ! Diff:
              i_diff=rles(text,1)
              j_diff=rles(text,2)
           elseIF(ilab.eq.15) then          ! start:
              i_Ausgabe=rles(text,1)
              goto 50
           elseIF(ilab.eq.16) then          ! VAR=:
              if(Schon_interpoliert.gt.0) Schon_interpoliert=2
              Variable=adjustl(cles(text,1))
              i=index(text,',')
              text=text(i+1:132)
              if(i_csv.le.0) stop ' VAR=, i_csv=0'
              if(j_csv.le.0) stop ' VAR=, j_csv=0'
              i1=rles(text,1)
              if(i1.le.0) then
                 i_csv=i_csv+1
                 i1=i_csv
                 if(i1.gt.maxvar) stop 'VAR zuviel Variablen'
              endif
              if(i1.gt.i_csv) stop 'VAR: VAR position zu gross'
              i2=rles(text,2)
              if(i2.gt.maxvar) stop 'VAR: PRES position zu gross'
              if(i2.le.0) stop 'VAR: PRES position <=0'
              if(i2.eq.i1) stop 'VAR: PRES position =  VAR position'
              i3=rles(text,3)
              if(abs(i3).gt.maxvar) stop 'VAR: TEMP position zu gross'
              if(abs(i3).le.0) stop 'VAR: TEMP position <=0'
              if(abs(i3).eq.i1) stop 'VAR: TEMP position =  VAR position'
              i4=rles(text,4)
              if(i4.le.0) stop 'VAR: S position <=0'
              if(i4.eq.i1) stop 'VAR: S position = VAR position'
              if(i4.gt.maxvar) stop 'VAR: S position zu gross'
              i5=rles(text,5)
              if(i5.lt.0) stop 'VAR: Oxyg position <=0'
              if(i5.eq.i1) stop 'VAR: Oxyg position =  VAR position'
              if(i5.gt.maxvar) stop 'VAR: Oxyg position zu gross'
              i6=rles(text,6)
              if(i6.lt.0) stop 'VAR: nutr. position <=0'
              if(i6.eq.i1) stop 'VAR: nutr. position =  VAR position'
              if(i6.gt.maxvar) stop 'VAR: nutr. position zu gross'
              rr=0
              rflag=0
              do j=1,j_csv
                 if(r_csv(j,i1).eq.def) then
                    if(i3.gt.0) then
                       rtemp=r_csv(j,i3)
                    elseif(r_csv(j,abs(i3)).ne.def) then
                       rtemp=r_csv(j,abs(i3))*1.00024        ! t90 in T68
                    else
                       rtemp=def
                    endif
                    r_csv(j,i1)=oz_variab(
     *                 r_csv(j,i2),rtemp,r_csv(j,i4),r_csv(j,i5),
     *                 r_csv(j,i6),rr,Variable,rflag)
                 endif
                 r_csv_aus(j,i1)=r_csv(j,i1)
              enddo
           elseIF(ilab.eq.17) then          ! B-multi:
              ! Vergrosserungsfaktor falls Zuwenig Daten im Bereich
              ! der Kernbereich (d.h. ohne prozent) wird bis ein
              ! maximum von Max_Faktor vergroessert
              ! Bereichfaktor muss >1 sein um effektiv zu sein
              ! Max_faktor muss grosser als Bereichfaktor sein
              ! sonst wird bis Faktor 100 vergroessert
              Bereichfaktor=rles(text,1)
              Max_faktor=rles(text,2)
              if(max_faktor.lt. BereichFaktor) Max_faktor=100.
           elseIF(ilab.eq.18) then          ! M-AUsgabe
              m_ausgabe=.true.
              m_csv=j_csv
           elseIF(ilab.eq.19) then          ! Info =0 keine Meldungen
              ! user zuwenig paramerte,etc
              ! geandert 17april96, frueher version mit info=1
              Info=rles(text,1)
           elseIF(ilab.eq.20) then          ! Maske, xpos,ypos,dateiname
              maske_x=rles(text,1)
              if(maske_x .gt. 0) then
                 maske_y=rles(text,2)
                 if(maske_y .le. 0) stop ' maskenfehler y<=0'
                 text=cles(text,3)
                 if(r_maske(text).ne.0) stop ' Maskenfehler'
              endif
           elseIF(ilab.eq.21) then        ! Z-Be: definitions Zyklenbereich
              i=rles(text,1)
              if(abs(i).gt.maxvar) then
                 stop ' Variablenposition zu gross'
              elseif(i.gt.0) then
                 Z_Bereich(i)%da =.true.
                 z_Bereich(i)%min=rles(text,2)
                 z_Bereich(i)%max=rles(text,3)
                 if(z_Bereich(i)%max.lt.z_Bereich(i)%min)
     *            stop 'Max<=min bei Zyklenbereich'
              elseif(i.lt.0) then
                 z_Bereich(abs(i))%da =.false.
              else
                 do i=1,maxvar
                    z_Bereich(i)%da =.false.
                 enddo
              endif
           elseIF(ilab.eq.22) then        ! NORM:
             if(L_Norm) then
               L_norm=.false.
             else
               L_norm=.true.
             endif
           elseIF(ilab.eq.23) then        ! Wweg:
              i1=rles(text,1)
              if(abs(i1).gt.maxvar) stop ' Variablenposition zu gross'
              r11=rles(text,2) ; r12=rles(text,3)
              i2=rles(text,4)
              if(abs(i2).gt.maxvar) stop ' Variablenposition zu gross'
              r21=rles(text,5) ; r22=rles(text,6)
              write(*,*) i1,r11,r12,i2,r21,r22
              do j=1,j_csv
                 if(r_csv(j,i1).gt.r11 .and . r_csv(j,i1).lt.r12 .and.
     &              r_csv(j,i2).gt.r21 .and . r_csv(j,i2).lt.r22) then
                    write(*,*) (r_csv(j,jj),jj=1,8)
                    r_csv(j,i2)=def
                 endif
              enddo
           elseIF(ilab.eq.24) then        ! eich:
              i1=rles(text,1)
              Bmin=rles(text,2)
              Bmax=rles(text,3)
              r0=rles(text,4)
              r1=rles(text,5)
              r2=rles(text,6)
              write(*,*) i1,bmin,bmax,r0,r1,r2,j_csv
              do j=1,j_csv
                 if(r_csv(j,i1).gt.Bmin .and . r_csv(j,i1).lt.Bmax) then
                    r_csv(j,i1)=r_csv(j,i1)+r0+r1*r_csv(j,i1)+r2*r_csv(j,i1)**2
                    if(r_csv(j,i1).lt.Bmin) r_csv(j,i1)=Bmin 
                    if(r_csv(j,i1).gt.Bmax) r_csv(j,i1)=Bmax
                    r_csv_aus(j,i1)=r_csv(j,i1)
                 endif
              enddo
           else
              write(*,*) ilab,text(1:40)
              stop ' nicht bekanntes label'
           endif
        enddo


   50   continue    ! ende, d.h. abarbeiten und fertig
        if(Schon_interpoliert.eq.1) goto 150
        Schon_interpoliert=1

        if(mod_ros.eq.0) then    ! moden auf rosettenwerten interpolieren
           if(imoden.le.0) stop ' Keine Modalzerlegung eingelesen'
           do i=1,imoden
              call jhlaint(moden(1,0),moden(1,i),jmoden,
     *         r_csv(1,idept),r_csv(1,i_csv+i),j_csv,2.,0.)
           enddo
           mod_ros=1
        endif

        do i=1,jint1
           if(int1(i).gt.i_csv) stop ' Nr. der int.Var. > max.Anzahl'
           if(int1(i).le.0) stop ' Nr. der int.Var. <=0'
        enddo
        do i=1,jpara
           if(para(i).gt.i_csv) stop ' Nr. der Par.Var. > max.Anzahl'
           if(para(i).le.0) stop ' Nr. der Par.Var. <=0'
        enddo

        do i=1,i_csv
           if(d_bereich(i)%da) then
              do j=1,j_csv
                 if(r_csv(j,i).lt.d_bereich(i)%min .or.
     *            r_csv(j,i).gt.d_bereich(i)%max) r_csv(j,i)=def
              enddo
           endif
        enddo
! interpoloation-1
        do ii=1,jint1          ! ueber anzahl variablen
           if(i_diff.eq.1 .or. i_diff.eq.2) then
              iw=iw+1
C              write(text,'(''Difout.v'',i2.2)') iw
              write(text,'(''Difout.v'',i2.2)') int1(ii)
              open(unit=20,file=text(1:10),status='new')
           endif
           meandif= 0
           stddif = 0
           idif=0
           difmax=-1e32
           difmin= 1e32
           n_interpol=0      !anzahl nicht interpolierter Werte
           j_interpol=0      !Anzahl interpolierter def.
           n_multi=0             !Anzahl der Vergroesserungsschritte

           do jj=1,m_csv      ! uber zu interpolierende Zyklen

              if(i_diff.eq.-1 .and. r_csv(jj,int1(ii)).ne.def) goto 888
              !Daten da, brauchen nicht interpoliert zu werden
              if(r_csv(jj,int1(ii)).eq.def) then
                 if(qualt(jj,int1(ii)).eq.0) then
                    qualt(jj,int1(ii))=-9
                 elseif(qualt(jj,int1(ii)).gt.0) then
                    qualt(jj,int1(ii))=qualt(jj,int1(ii))+2
                 endif
                 if(i_diff.eq. 4 ) goto 888
              endif
C i_diff=3 zum setzen der qualityflags bei allen Daten
C i_diff=4 ist zum setzen der qualityflags bei vorhandenen Daten
C flags bei vorhandenen Daten 0=gut, 
C pro interpolation (siehe D=in:) wird aufaddiert: +1=differenz>stdabw +2=differenz>3*stdabw
C oder +2 falls quality schon auf >0 (d.h. schon als fehlerhaft erkannt)
C flags bei fehlenden daten -9=fehlwert, -1=interpoliert, 
C   -2=interp.aber auf min gesetzt, -3=interp. aber auf max gesetzt

              ! suchen der moeglichen interpolationsvariablen, dieses haengt
              ! von den Variablen im zu interpolierenden Zyklus ab.
              l=0
              do i=1,jpara
                 if(int1(ii).ne.para(i) .and. r_csv(jj,para(i)).ne.def) then
                    l=l+1
                    ivar(l)=para(i)
                 endif
              enddo
              if(l.le.2) then
                 if(info.ne.0) then
                    write(*,1102) (r_csv(jj,kk),kk=1,5)
                    write(*,1030) int1(ii),(ivar(i),i=1,l)
                    write(*,1032)
                 endif
                 n_interpol=n_interpol+1
                 goto 888
              endif

              ! Bestimmen des Bereiches aus dem die Zyklen zur Bestimmung der MLR
              ! verwendet werden sollen
              rmulti=1.      ! Normaler Bereich
              k=1+l+imoden

  756         do kk=1,irange
                 if(r_csv(jj,range(kk)%var).ne.def) then
                    r= rmulti*range(kk)%r
     *               + r_csv(jj,range(kk)%var)*range(kk)%pro
                    rang_min(kk)=r_csv(jj,range(kk)%var) - r
                    rang_max(kk)=r_csv(jj,range(kk)%var) + r
                 else
                    rang_min(kk)=-1e32
                    rang_max(kk)=+1e32
                 endif
              enddo

              ! auffuellen des Vektors b und der Matrix a
              ! Die Auswahl der Zyklen erfolgt mit diesen Kriterien:
              ! a) Alle Regressionsvariablen mussen vorhanden sein
              ! b) Zyklus muss in dem Bereich liegen
              ! c) Die Verbindungsgerade zwischen dem zu interpolierenden
              !    Zyklus und dem aktuellen Punkt darf die in der Maske
              !    definierten Geradenstuecke nicht schneiden
              ! d) zu interpolierender Zyklus wird nicht verwendet

              n=1
              do j=1,j_csv
                 if(jj.eq.j) goto 777                  ! Kriterium D
                 if(r_csv(j,int1(ii)).eq.def) goto 777      ! Kriterium A
                 do i=1,l
                    if(r_csv(j,ivar(i)).eq.def) goto 777
                 enddo
                 do kk=1,irange                        ! Kriterium B
                    if( r_csv(j,range(kk)%var) .ne.def .and.
     *               (r_csv(j,range(kk)%var).gt.rang_max(kk) .or.
     *               r_csv(j,range(kk)%var).lt.rang_min(kk)) ) goto 777
                 enddo
                 ! Kriterium C;
                 if(maske_x.gt.0) then
                    ! falls Maskenvariablen=Fehlwert wird die Maske NICHT ueberprueft
                    if(r_csv(jj,maske_x) .ne. def .and. r_csv(jj,maske_y) .ne. def .and.
     *              r_csv( j,maske_x) .ne. def .and. r_csv( j,maske_y) .ne. def) then
                       if( P_maske(r_csv(jj,maske_x),r_csv(jj,maske_y),
     *                  r_csv( j,maske_x),r_csv( j,maske_y)) ) goto 777
                    endif
                 endif

                 b(n)=r_csv(j,int1(ii))
                 a(n,1)=1.
                 do i=1,l
                    a(n,i+1)=r_csv(j,ivar(i))
                 enddo
                 do i=1,imoden
                    a(n,l+1+i)=r_csv(j,i_csv+i)
                 enddo
                 ! Daten der gleichen Station werden hoeher gewichtet, falls LAbel STAT: angegeben
                 if(istat.gt.0. .and.
     *           r_csv(j,range(istat)%var).eq.r_csv(jj,range(istat)%var)) then
                    b(n)=b(n)*5.
                    do i=1,k
                       a(n,i)=a(n,i)*5.
                    enddo
                 endif
                 n=n+1
  777            continue
              enddo
              n=n-1

              if(n.le.3+k) then
                 if(Bereichfaktor.gt.1.) then
                    rmulti=rmulti*Bereichfaktor
                    if(max_faktor.gt.rmulti) then
                       n_multi=n_multi+1
                       !! write(*,1018) int1(ii),rmulti
                       goto 756
                    endif
                 endif
                 if(info.ne.0) then
                    write(*,1022)
                    write(*,1102) (r_csv(jj,kk),kk=1,5)
                    write(*,1020) int1(ii),n,k
                 endif
                 if(rmulti.gt.1.and.info.ne.0)
     *            write(*,1024) rmulti/Bereichfaktor
                 n_interpol=n_interpol+1
                 goto 888
              endif
		
C Normieren
              if(L_Norm) then
                do i=1,l
                 a(n+1,i+1)=r_csv(jj,ivar(i))
                 ! zu interpolierende Daten werden mit benutzt
                 call moment(a(1,i+1),n+1,ave,adev,sdev,varian,skew,curt)
                 fehler=d_bereich(ivar(i))%stdabw
                 if(fehler.le.0) fehler=1
                 do j=1,n+1
                  a(j,i+1)=(a(j,i+1)-ave)*(sdev/fehler)
                 enddo
                 xint(i+1)=a(n+1,i+1) ! auf xint stehen jetzt die normierten werte
                enddo 
              else
                ! die Werte des zu interpolierendem Zyklus werden auf xint kopiert
                 do i=1,l
                    xint(i+1)=r_csv(jj,ivar(i))
                 enddo
              endif
C Eigenwertzerlegung
              ! x bestimmen mit SVD
              call svdcmp(a,n,k,5000,maxvar,w,v)
              ! setzen von kleinen Eigenwerten=0
              wmax=0
              do i=1,k
                 wmax=max(wmax,w(i))
              enddo
              do i=1,k
                 if( w(i).le. (wmax*1.0e-6)) w(i)=0
              enddo
              call svbksb(a,w,v,n,k,5000,maxvar,b,x)
C den interpolierten Wert berechnen
              ! anwenden von x auf die Rosettendaten
              r8=x(1)
              do i=1,l
                 r8=r8+ x(i+1)*xint(i+1)
              enddo
              do i=1,imoden
                 r8=r8+ x(i+l+1)*r_csv(jj,i_csv+i)
              enddo
              if(i_diff.eq.2) write(20,'(1x,3i5,f12.3,20(i3,2f12.3))')
     &            n,k,l,r8,
     &            int1(ii),x(1),r_csv(jj,int1(ii)),
     &            (ivar(i),x(i+1),xint(i+1),i=1,l)
C     &            (ivar(i),x(i+1),r_csv(jj,ivar(i)),i=1,l)
 
              if(qualt(jj,int1(ii)).lt.0) qualt(jj,int1(ii))=-1
              ! Interpolierte werte ausser Definitionsbereich auf max. bzw. min.
              ! setzen
              if(d_bereich(int1(ii))%da) then
                 if(r8.lt.d_bereich(int1(ii))%min) then
                    r8=d_bereich(int1(ii))%min
                    qualt(jj,int1(ii))=-2
                 elseif(r8.gt.d_bereich(int1(ii))%max) then 
                    r8=d_bereich(int1(ii))%max
                    qualt(jj,int1(ii))=-3
                 endif
              endif

C zur ueberprufen ob hier nicht if(r_csv(jj,int1(ii)).eq.def) then  hinkommen sollte
              ! if(qualt(jj,int1(ii)).lt.0) then
              if(r_csv(jj,int1(ii)).eq.def) then
                 j_interpol=j_interpol+1
                 r_csv_aus(jj,int1(ii))=r8
              elseif(qualt(jj,int1(ii)).eq.0) then
                 idif=idif+1
                 dif=r_csv(jj,int1(ii))-r8
                 if(abs(dif).gt.  d_bereich(int1(ii))%stdabw) qualt(jj,int1(ii))=qualt(jj,int1(ii))+1
                 if(abs(dif).gt.3*d_bereich(int1(ii))%stdabw) qualt(jj,int1(ii))=qualt(jj,int1(ii))+2
                 meandif= meandif+ dif
                 stddif = stddif + dif*dif      ! annahme mittelwert=0
                 difmin=min(difmin,dif)
                 difmax=max(difmax,dif)
                 if(i_diff.eq.1) write(20,'(1x,8f10.3,3i5)')
     *            (r_csv(jj,Vdiff(kk)),kk=1,5),r_csv(jj,int1(ii)),r8,dif,n,k,l
C     *            (r_csv(jj,kk),kk=1,5),r_csv(jj,int1(ii)),r8,dif,n,k,l
              endif
              if(j_diff.eq.1) r_csv_aus(jj,int1(ii))=r8


  888         continue

           enddo

           write(*,1002) int1(ii),j_interpol
           if(n_multi.gt.1) write(*,1004) N_multi
           write(*,1006) n_interpol
           if(idif.gt.1) then
              write(*,1008)  idif
              write(*,1010) difmin,difmax
              write(*,1012) meandif/idif
              write(*,1014) sqrt(stddif/(idif-1))
           endif
           if(i_diff.ge.1) close(20)

        enddo
        if(i_ausgabe.eq.0) goto 5

 150    if(AUS_format(1:2).eq.'##') then
           write(AUS_format,5000) i_csv,jint1
 5000      format('(1x,',i2,'F15.3,',i2,'i3)')
           write(*,'(a)') aus_format
        endif
        if(i_diff.eq.3 .or. i_diff.eq.4) then
           do j=1,m_csv         ! wegzuschreibende Zyklen
              write(33,AUS_format) (r_csv_aus(j,i),i=1,i_csv),(qualt(j,int1(i)),i=1,jint1)
           enddo
        else
           do j=1,m_csv         ! wegzuschreibende Zyklen
              write(33,AUS_format) (r_csv_aus(j,i),i=1,i_csv)
           enddo
        endif
        if(i_ausgabe.lt.0) stop
        if(i_ausgabe.eq.1) then
           i_ausgabe=-1
           goto 5
        endif

        stop
 1002   format (' Variable:',i3,'  Interpolierte Fehlwerte:',i8)
 1004   format ('  Anzahl der Bereichsmultiplikationen:',i8)
 1006   format ('  nicht interpolierte Werte:',i8)
 1008   format (' Anzahl der Differenzwerte:',i8)
 1010   format (' Min,max  Differenz= ',2f10.3)
 1012   format (' Mittlere Differenz= ',f10.3)
 1014   format (' Standardabweichung (bei mittelwert=0)= ',f10.3)

 1018   format (' Var=',i3,' Bereich vergroessert um Faktor=',f10.3)
 1020   format (' Zuwenig Daten zur MLR; Var=',i3,'  Anzahl der Messwerte=',i8
     *   ,' Anzahl Par.',i3)
 1022   format (' nicht interpoliert, zuwenig Daten fuer die MLR')
 1024   format (' Bei maximaler Bereichsvergrosserung von ',f10.3)

 1030   format (' Var=',i3,' Parameter:',30i3)
 1032   format (' Zuwenig Parameter, nicht interpoliert ')

 1102   format (' werte(1-5)=',5f12.3)
        end                                                                                                                                 
!***********************************************************************                                                            
!***********************************************************************                                                            
!***********************************************************************                                                            
