CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C23456789012345678901234567890123456789012345678901234567890

      Program Wilgotnosc

      real ts, tw, p, t(1000), Ew(1000), El(1000)
      real t1, t2, EE1, EE2, EE, e, EEw, Es
      real a, q, D, td
      integer i, j, nrec
      character znak
      Write(*,*)
      Write(*,*) '----------------------------------------------'
      Write(*,*) ' PROGRAM WILG1'
      Write(*,*) ' Autor Krzysztof Markowicz'
      Write(*,*) 'e-mail: kmark@igf.fuw.edu.pl'
      Write(*,*) 'www.igf.fuw.edu.pl/meteo/stacja'
      Write(*,*)
      Write(*,*) 'Program oblicza wilgotnosci powietrza'
      Write(*,*) 'na podstawie temperatury z termometru suchego'
      Write(*,*) 'oraz zwilzonego'
      Write(*,*) '----------------------------------------------'
      Write(*,*)
      Write(*,*) 'Podaj temperature termometru suchego:'
      Read(*,*) ts
      Write(*,*) 'Podaj temperature termometru wilgotnego:'
      Read(*,*) tw
      Write(*,*) 'Podaj cisnienie powietrza w [hpa]:'
      Read(*,*) p
      Write(*,*) 'Podaj w dla wody oraz l dla lodu'
      Write(*,*) 'okreslajace stan batystu na termometrze wilgotnym'
      Read(*,*) znak

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC read data from file
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      open(21,file='cis_pary.dat',status='unknown')
      read(21,*) nrec
      do i=1, nrec
      read(21,*) t(i), Ew(i), El(i)
      end do
      close(21)
      if (znak.eq.'w') then
       do i=1,nrec
        if (tw.ge.t(i)) then
         EE2=Ew(i+1)
         EE1=Ew(i)
         t2=t(i+1)
         t1=t(i)
         EEw=EE1+(EE2-EE1)/(t2-t1)*(tw-t1)
         e=EEw-0.0007946*p*(ts-tw)
        end if
        if (ts.ge.t(i)) then
         EE2=Ew(i+1)
         EE1=Ew(i)
         t2=t(i+1)
         t1=t(i)
         Es=EE1+(EE2-EE1)/(t2-t1)*(ts-t1)
        end if
       end do
      end if
      if (znak.eq.'l') then
        do i=1,nrec
         if (tw.ge.t(i)) then
          EE2=El(i+1)
          EE1=El(i)
          t2=t(i+1)
          t1=t(i)
          EEw=EE1+(EE2-EE1)/(t2-t1)*(tw-t1)
          e=EEw-0.0007060*p*(ts-tw)
         end if
         if (ts.ge.t(i)) then
          EE2=Ew(i+1)
          EE1=Ew(i)
          t2=t(i+1)
          t1=t(i)
          Es=EE1+(EE2-EE1)/(t2-t1)*(ts-t1)
         end if
        end do
       end if
      f=e/Es*100
      a=217*e/(273.16+ts)
      q=622*e/p
      D=Es-e

C     Oblicza temperature punktu rosy.
      do i=1,nrec
       if (e.ge.Ew(i)) then
        t1=t(i)
        t2=t(i+1)
        E1=Ew(i)
        E2=Ew(i+1)
        td=t1+(t2-t1)/(E2-E1)*(e-E1)
       end if
      end do

      Write(*,*) '************************************************'
      Write(*,*) 'Wilgotnosc wzgledna w [%]:'
      Write(6,'(f5.1)') f
      Write(*,*) 'Wilgotnosc bezwzgledna w [g/m^3]'
      Write(6,'(f5.1)') a
      Write(*,*) 'Wilgotnosc wlasciwa w [g/kg]'
      Write(6,'(f5.1)') q
      Write(*,*) 'Cisnienie pary wodej w [hpa]'
      Write(6,'(f5.1)') e
      Write(*,*) 'Niedosyt wilgotnosci powietrza w [hpa]'
      Write(6,'(f5.1)') D
      Write(*,*) 'Temperatura punktu rosy'
      Write(6,'(f5.1)') td
      end