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