c********************************************************************
C
C  PROGRAMM zur Auswertung der Disdrometerdaten
C
c     (c) by DLR 1999
c
c     Author: Martin Hagen
c             DLR, Oberpfaffenhofen, D-82234 Wessling, Germany
c     Contact: martin.hagen@dlr.de
c
c     Groessenklassen um 16 % vermindert, nach Vergleich mit 
c     den Regensummen der Automatischen Station (ASTA) am OTL
c     (nach pers. Mitteilung von Juerg Joss)
c
c     5-min-Mittelung der Tropfenspektren und integralen
c     Groessen berichtigt
c
c     Letzte Aenderung: 
c      10.11.99   Ulrich Blahak: Fehler und 16% Verminderung
c      09.06.00   mh: angepasst fuer alte Daten
c
c**********************************************************************

      DIMENSION dg(21), dm(20), rdwvel(20), spekt(20,1440),
     &          dm3(20), dm6(20), smk(20,20), iread(30), 
     &          nspek(20), nkspek(20), vel(20), dwidth(20), 
     &          nktropf(1440), ntropf(1440), spektj(20),
     &          ombro(1440), temp(1440), wind(1440)

      character indsn*80, outdsnt*80, text*120, zeile*120, cmon*3, ck*2,
     &          ctmon(12)*3, outdsnd*80, outdsns*80, c0d*1, cext*4

      data dg /0.40,0.56,0.72,0.80, 0.90,0.97,1.11,1.31,
     &         1.50,1.67,1.81,2.10, 2.42,2.68,3.03,3.28,
     &         3.69,4.11,4.47,4.92, 5.52/   
C                 ! letzter Wert nur damit die Klasse 20 eine Breite hat
      data ctmon /'Jan','Feb','Mae','Apr','Mai','Jun','Jul','Aug','Sep',
     &            'Okt','Nov','Dez'/
      data spekt /28800*0.0/, ntropf /1440*0/

c------ Vorbereitungen -------------------------------------------------

C --- Groessenklassen um 16% vermindern:
      do i=1,20
         dg(i) = 0.84 * dg(i)
      end do

C  Integrationsintervall, Messflaeche
      tint = 60.                               ! Messintervall (Sek.)
      area = 0.005                             ! Flaeche (m^2)
      rfakt = 0.0018849 / tint / area
      r60 = 1./60.
      c0d = char(13)

C  Mittleren Durchmesser, Klassenbreite und Fallgeschwindigk. bestimmen
      do 10 i=1,20
        dm(i) = (dg(i)+dg(i+1)) * 0.5          ! mittlerer Durchmesser
        dm3(i) = dm(i)**3
        dm6(i) = dm3(i) * dm3(i)
        dwidth(i) = dg(i+1) - dg(i)            ! Klassenbreite
        vel(i) = 9.65 - 10.3*exp(-0.6*dm(i))   ! Fallgeschwindigkeit
        rdwvel(i) = 1. / (dwidth(i) * vel(i) * area * tint) 
 10   continue

C  Korrektur Matrix fuer Totzeit (nach Waldvogel, Steiner, Henrich)
      do 13 i=1,20
        tmk = 0.85 * dm(i)
        do 14 j=1,20
          if (dm(j) .lt. tmk) then
            smk(i,j) = 0.0
          else
            smk(i,j) = log(dm(j) / (0.85*(dm(i)-0.25)))
          endif
 14     continue
 13   continue

C------ Daten lesen ----------------------------------------------------

C  Eingabedatei
      cext = '.dsd'                             ! Filename Extension
      print *, 'Eingabedatei ?'
      read(*,'(a)') indsn
      if (index(indsn,cext) .eq. 0) indsn(index(indsn,' '):) = cext
      write(*,*) indsn(1:index(indsn,' '))
      open(unit=10,file=indsn,access='sequential',status='old')

C  Header mit Datum lesen
      read(10,'(//A/)',end=29,err=29) text
      read(text,'(38x,i2,2x,a3,2x,i4)') itag, cmon, jahr
      do 24 i=1,12
        if (cmon .eq. ctmon(i)) monat = i
 24   continue
      write(*,*) 'Disdrometer data from ', itag, monat, jahr

C  Zeile fuer Zeile lesen
      iline = 0
 27   continue
      read(10,'(A)',end=29) text
      if (text(1:1) .eq. char(0)) goto 27
      i0d = index(text,c0d)-1                ! fuer DOS Dateien
      if (i0d .lt. 0) i0d = len(text)-2      ! fuer UNIX Dateien
      zeile(1:i0d) = text(1:i0d)
      zeile(i0d+1:i0d+2) = ' ;'

c  Entschluesseln der eingelesenen Zeile
      call zsplit(zeile,iread,nr)
      if (nr .eq. 0) goto 27
      if (nr .ge. 1) iline = iline + 1

C  Bearbeitung von unterschiedlichen Datensaetzen:
C   MAP Daten sind immer mit Uhrzeit, andere nur alle 30 Minuten
C   MAP Daten nur wenn Niederschlag (bzw. alle 30 Min), 
C                                      andere kontinuierlich
C   5 oder 26 Werte                  3, 5, 24, 26 
C  Datenzaetze ohne Uhrzeit werden ergaenzt
      if (nr.eq.3 .or. nr.eq.24) then
        do i=nr,1,-1
          iread(i+2) = iread(i)
        enddo
        nr = nr + 2
        izeit = izeit + 1
        iread(1) = izeit / 60
        iread(2) = mod(izeit,60)
      endif

C  Zeit
      istd = iread(1)
      imin = iread(2)
      izeit = istd*60 + imin
cx      print *, iline, nr, iread(1), iread(2), izeit

C  Messungen von Ombrometer, Thermometer und Anemometer speichern
      wind(izeit+1) = iread(3) * 0.1
      temp(izeit+1) = iread(4) * 0.1
      ombro(izeit+1) = iread(5) * 0.1
C   falsche Eichung Ombrometer korrigieren
      if (ombro(izeit+1).gt.0.0 .and. ombro(izeit+1).lt.9.0) then
        itropf = nint((-4.296e-3 + sqrt(1.8456e-5 -1.7712e-5*
     &       (0.0171-(ombro(izeit+1)/60.)))) / 8.856e-6)
        ombro(izeit+1) = 0.3 * max(0., itropf)
      endif

C  Gesamtzahl der Tropfen (Spektren mit weniger als 10 werden ignoriert)
      if (nr .ge. 7) then                    ! Disdro Daten vorhanden
        ntot = 0                             ! Gesamtzahl der Tropfen
        do 31 i=1,20
          ntot = ntot + iread(i+6)
 31     continue
        if (ntot .lt. 10) then
          nr = 0
cxcx      print *, iread(1), iread(2), ' too few drops ', ntot
        endif
      endif

C  Testsignal (Kanal 6-8 sehr viel mehr Tropfen als im Rest)
      if (nr .ge. 7) then                    ! Disdro Daten vorhanden
        itest = iread(6+6) + iread(6+7) + iread(6+8)
        if (itest .gt. ntot*0.7) then
          nr = 0
          print *, iread(1), iread(2), ' test signal ', itest, ntot
        endif
      endif

C  Auswertung der Spektren
      if (nr .ge. 7) then                    ! Disdro Daten vorhanden

        aktiv = -iread(6) * 0.01

C  Totzeitkorrektur (nach Waldvogel, Steiner, Henrich)
        nktot = 0                            ! korrigierte Gesamtzahl
        do 35 i=1,20
          nspek(i) = iread(i+6)              ! unkorrigierte Anzahl
          if (nspek(i) .gt. 1) then
            tkorr = 0.0
            do 36 j=1,20
              tkorr = tkorr + nspek(j)*smk(i,j)
 36         continue
            tkorr = exp((0.035/tint)*tkorr)
            nkspek(i) = nspek(i) * tkorr     ! korrigierte Anzahl
            nktot = nktot + nkspek(i)
          else
            nkspek(i) = nspek(i)
          endif
 35     continue
      
C  Spektren (aus korrigierter Anzahl)
        do 41 i=1,20
          spekt(i,izeit+1) = aktiv*nkspek(i)*rdwvel(i) ! korrig. Spekt.
 41     continue
        ntropf(izeit+1) = ntot    ! Gesamtzahl
        nktropf(izeit+1) = nktot  ! korrigierte Gesamtzahl
      
      endif

      goto 27

 29   continue

C  Alle Spektren sind eingelesen, zuerst werden die 1 Minuten Spektren
C  bearbeitet und ausgegeben. Danach laeuft das ganze noch mal fuer
C  die 5 Minuten Werte durch. 

      do 50 k=1,5,4

        write(ck,'(I2.2)') k
C  Ausgabedatei fuer Spektrumparameter (Regenrate, Reflektivitaet, ...)
        outdsnt = indsn
        outdsnt(index(indsn,cext):) = '_c.'//ck//'t'
        open(unit=10+k,file=outdsnt,access='sequential')
        write(10+k,'(2a)')'Calculated with 16% diminished', 
     &                   ' drop diameters'
        write(10+k,811)
 811    format('   Date  ',' Time ','Time (d)',
     &         ' N_tot',' Nktot','  Rainr','  Refl','    LWC',
     &         '    N_0 ','   D_0',' R_sum',' RR_Omb',' R_s_O')

C  Ausgabedatei fuer Spektren
        outdsnd = indsn
        outdsnd(index(indsn,cext):) = '_c.'//ck//'d'
        open(unit=11+k,file=outdsnd,access='sequential')
        write(11+k,821) dm
        write(11+k,'(2a)')'Calculated with 16% diminished', 
     &                   ' drop diameters'
 821    format(' Drop size distribution N(D) for mean diameters'/
     &         '   Date  ',' Time ','Time (d)', 20f5.2)

C  Ausgabedatei fuer Zusammenfassung
        outdsns = indsn
        outdsns(index(indsn,cext):) = '_c.'//ck//'s'
        open(unit=12+k,file=outdsns,access='sequential')

        reratag = 0.0           ! Gesamtregen
        reramax = 0.0           ! Max. Regenrate
        nregen  = 0             ! Zahl der Regeninterval.
        rsumomb = 0.0           ! Gesamtregen Ombrometer

        sx = 0.0                ! Groessen fuer Z-R Regr.
        sy = 0.0                !
        sxx = 0.0               !
        sxy = 0.0               !
        syy = 0.0               !
        nreg = 0                !

        do 51 j=1,1440,k

C  Die Zeiten geben die Zeit an, zu der das angegebene Spektrum
C  beendet wurde: 
C    z.B. bei 1 Min. Daten:  12:00 =>  11:59:00 bis 11:59:59.999
C    z.B. bei 5 Min. Daten:  12:00 =>  11:55:00 bis 11:59:59.999
          jk = j+k-1
          if (k .eq. 1) jk = jk - 1
          istd = jk / 60
          imin = mod(jk,60)
          tzeit = istd + imin/60.

C  Spektrum fuer das Zeitintervall aufsummieren
          ntropfj = 0
          nktropfj = 0
          do 52 i=1,20
            spektj(i) = 0.0
 52       continue
          do 53 jk=1,k
            ntropfj = ntropfj + ntropf(j+jk-1)
            nktropfj = nktropfj + nktropf(j+jk-1)
            do 54 i=1,20
              spektj(i) = spektj(i) + spekt(i,j+jk-1) / k
 54         continue
 53       continue

C  Ombrometermessungen
          rromb = 0.0
          do jk=1,k
            rromb = rromb + ombro(j+jk-1) / k
          enddo
          rsumomb = rsumomb + rromb/60.*k

          if (ntropfj .gt. 0) then
C  Regenrate, Reflektivitaet, LWC (aus korrigierter Anzahl)
            rlwc = 0.0          ! Regenwassergehalt D^3
            rera = 0.0          ! Regenrate
            reflin = 0.0        ! Reflektivitaet D^6

            do 43 i=1,20
              rlwc = rlwc + spektj(i) * dwidth(i) * dm3(i) * 0.5236
              rera = rera + spektj(i)*dwidth(i)*vel(i)*dm3(i)*0.0018849
              reflin = reflin + spektj(i) * dwidth(i) * dm6(i)
 43         continue
c
c     Fehler: Mittelung wird schon bei den Spektren gemacht!
c            rlwc = rlwc/k
c            rera = rera/k
c            reflin = reflin/k

            refdbz = 10. * log10(reflin) ! Reflektivitaet in dBZ

C  Z - R Beziehung (Werte nur fuer Regenraten ueber 1 mm/h)
            if (rera .ge. 1.0) then
              nreg = nreg + 1
              sx = sx + log10(rera)
              sy = sy + refdbz * 0.1
              sxx = sxx + log10(rera)**2
              sxy = sxy + log10(rera)*refdbz*0.1
              syy = syy + refdbz*0.1*refdbz*0.1
            endif

C  Regenstatistik
            nregen = nregen + 1
            reratag = reratag + rera * tint * k / 3600.
            reramax = max(rera,reramax)
            

C  N0 und D0 ueber Momente (Waldvogel, 1974)
            rnnull = 446. * rlwc * (rlwc/reflin)**1.333333 ! N null
            dnull = 0.5997 * (rlwc/reflin)**(-0.333333) ! D null

C  Ausgabe in Datei
            write(10+k,812) jahr, monat, itag, istd, imin, tzeit,
     &                      ntropfj, nktropfj, rera, refdbz, rlwc, 
     &                      rnnull, dnull, reratag, rromb, rsumomb
 812        format(i5,i2.2,i2.2,1x,i2.2,':',i2.2,f8.4,
     &             i6,i6,f7.2,f6.1,f7.1,
     &             f8.0,f6.2,f6.1,f7.2,f6.1)

            write(11+k,822) jahr, monat, itag, istd, imin, tzeit, 
     &                      (spektj(i),i=1,20)
 822        format(i5,i2.2,i2.2,1x,i2.2,':',i2.2,f8.4,
     &             20f8.1)
          endif

 51     continue
C  Berechnung a und b fure Z - R Beziehung
        if (nreg .eq. 0) then
          regb = 0.0
          rega = 0.0
        else
          regb = (nreg*sxy - sx*sy) / (nreg*sxx - sx*sx)
          rega = sy/real(nreg) - regb*sx/real(nreg)
        endif

C     Ausgabe Summary
        write(   *,*)
        write(   *,831)nregen,k,reratag,reramax,nint(10**rega),regb,nreg
        write(12+k,831)nregen,k,reratag,reramax,nint(10**rega),regb,nreg
 831    format(' Number of precip. intervalls:',i5,' (',i2,
     &         ' min. interv.)'/
     &         ' Precipitation sum:',f7.1,' mm'/
     &         ' Maximum rain rate:',f7.1,' mm/h'/
     &         ' Z - R relation: Z =',i4,' R^',f4.1,' (',i5,' interv.)')
        write(12+k,'(2a)')'Calculated with 16% diminished', 
     &                   ' drop diameters'

 50   continue

      STOP
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      subroutine zsplit(zeile,iread,nr)
c  zerlegt ZEILE in einzelne Integer Werte
      dimension iread(*)
      character zeile*(*), cb*1, cl*1, cj*1
      data cb/' '/

      lenz = index(zeile,';') - 1
c  Uebergang von Blank auf Zeichen zaehlen
      cl = zeile(1:1)
      if (cl .ne. cb) then
        nr = 1
      else
        nr = 0
      endif
      do 10 i=2,lenz
        cj = zeile(i:i)
        if (cl.eq.cb .and. cj.ne.cb) nr = nr + 1
        cl = cj
 10   continue

c  Werte lesen
      if (nr .gt. 0) then
        read(zeile,*) (iread(i),i=1,nr)
      else
        iread(1) = 0
      endif

      return
      end 
