c     mr.f Calculates average and frequency of soundings.
c     This program considers both Capital and noncapital
c     letters in filenames
      parameter (n=70)          !  n = rows in file unit 1
      character*3 date
      character*4 time
      character*4 ext(2)        !extension file: WIN or win
      character*1 fore(2)       !fore part of the file: P or p
      character*12 filename(4)
      character*9 fileoutput
      character*2 month         
      character*2 amopm
      character*2 stn
      integer*1 go
      dimension u(n),v(n),hm(n)
      integer j(n),tmax
      data ext/'.win','.WIN'/
      data fore/'p','P'/
c
c    Initialize variables
      izero = 0
      tmax = 0
      ie = 0
      num = 0
      do 1 i=1,n
        j(i)=0
        hm(i)=0
        u(i)=0.
        v(i)=0.
1     continue
c
c    Enter information from keyboard
      write(*,*)'station number:'
          read(*,'(a2)')stn
      write(*,*)'input month (ex: 05 for May)'
          read(*,'(a2)')month
      write(*,*)'choose am, pm or total (am,pm or to ,respectively)'
          read(*,'(a2)')amopm
c
c    Use information from keyboard ...
      CALL days(lim,month)      !..to define "lim" (tot.number of days)
      fileoutput  = fore(1)//stn//month//amopm//'.s'
      write(*,*)'This month has ',int(lim/2.),'days'
      write(*,*)'The mean will be written in ',fileoutput
c
c    Opening fixed files
      open(2,file= fileoutput)              !output file
      write(2,11)                           !first row in output file
      open(3,file='fechas.d',status='old')  !reading file
      open(4,file='comments')
c
c    Main loop
      do 5 ii=1,lim
       CALL morni_after(date,num,amopm,go)
       if(go.eq.1)goto 5                !reject file
       filename(1) = fore(1)//stn//month//date//ext(1)
       filename(2) = fore(2)//stn//month//date//ext(2)
       filename(3) = fore(1)//stn//month//date//ext(2)
       filename(4) = fore(2)//stn//month//date//ext(1)
c      Open wind data file (either lower or upper case filename)
       CALL readwri(filename,n,hm,u,v,j,tmax,count,izero)
5     continue
c
c    Writing means
35    write(2,13)(lim-num)-izero,izero,(lim-num)
      write(2,14)tmax/float((lim-num)-izero)
      write(2,*)
      write(2,*)' height   mag    dir      u       v    freq '
      do 40 k=1,n
        if(j(k).ne.0)then
           umean=u(k)/float(j(k))
           vmean=v(k)/float(j(k))
           xmag=sqrt(umean**2+vmean**2)
           dir=asin(-umean/xmag)*180/3.1416
           if(vmean.gt.0)dir=180-dir
           if(umean.gt.0.and.vmean.lt.0)dir=360+dir
           write(2,12)hm(k),xmag,dir,umean,vmean,j(k)
        endif
40    continue
c
11    format(36h    file           No.obs   MaxH(m) )
12    format(x,f6.0,2x,f5.2,2x,f5.1,2x,f6.2,2x,f6.2,4x,i2)
13    format(/,10hobserved= ,i2,3x,9hmissing= ,i2,3x,7htotal= ,i2)
14    format(/,47hAveraged height (excluding files with 0 obs.)= ,f6.0)
      stop
      end
c
       subroutine days(lim,month)
       character*2 month         
       lim=0
       if(month.eq.'01'.or.month.eq.'03'.or.month.eq.'05'.or.
     .    month.eq.'07'.or.month.eq.'09'.or.month.eq.'11')then
          lim=62
       else
          lim=60                 !lim=twice days of the month       
       endif
       if(month.eq.'02')lim=56   !change in case of lap year
       return
       end
c
       subroutine morni_after(date,num,amopm,go)
       character*3 date
       character*2 amopm
       integer*1 go 
       go = 0
       read(3,11)date
        if(amopm.ne.'to')then
          read(date,'(i3)')idate
          if(amopm.eq.'pm')then
             if(mod(idate,2).ne.0)then       !afternoon
              num = num +1
              go = 1
             endif
          endif
          if(amopm.eq.'am')then
             if(mod(idate,2).eq.0)then       !morning
              num = num +1
              go = 1
             endif
          endif
        endif
11     format(1x,A3)
       return
       end
c
       subroutine readwri(filename,n,hm,u,v,j,tmax,count,izero)
       character*12 filename(4)
       character*5 label         !to differentiate radiosonde from pibal
       dimension u(n),v(n),hm(n)
       integer j(n),tmax
c
c      Open wind data file (either lower or upper case filename)
       open(1,file=filename(1),status='old',iostat=ierr)      !input file
       if(ierr.ne.0)then
        close(1)
        open(1,file=filename(2),status='old',iostat=ierr2)    !input file 2nd. option
       endif
       if(ierr2.ne.0)then
        close(1)
        open(1,file=filename(3),status='old',iostat=ierr3)    !input file 3rd. option
       endif
       if(ierr3.ne.0)then
        close(1)
        open(1,file=filename(4),status='old',iostat=ierr4)    !input file 4th. option
       endif
       if(ierr4.ne.0)then                                     !data file not found
         write(2,*)'   ',filename(1),'       0'
         write(4,*)'file  ',filename(1),'   not found'
         izero = izero+1
         return                                               !leave the subrutine
       endif
c
       read(1,*,end=20)
       read(1,'(a5)',end=20)label                             !second row to distinguish pibal 
c
       open(10,file='altura30')                               !open file with correct heights
       iban = 0
       do 10 k=1,n
         if(label.eq.' t[s]')then
           read(1,11,end=40,err=30)time,zt,x,x,x,x,x,u0,v0    !pibal data
         else
           read(1,12,end=40,err=30)time,zt,x,x,x,x,x,u0,v0    !radiosonde data
         endif
c
c        input correct height
         read(10,*,end=30)time2,zt2
         if(zt.eq.zt2.and.iban.eq.1)then
           iban=0
           write(4,*)'gotcha again!',zt,zt2
         endif
         if(zt.ne.zt2)then
           write(4,*)'error in Height ',filename(1)
           write(4,*)zt,zt2
           iban=1
           goto 10
         endif
         hm(k)=zt
         u(k)=u(k)+u0
         v(k)=v(k)+v0
         j(k)=j(k)+1
         maxz=int(zt)
         count=k
10     continue
c
c       Finally write!
40     tmax=tmax+maxz
       write(2,13)filename(1),int(count),maxz
       rewind(10)
       return
c
20     write(4,*)'No data found in file ',filename(1)
       rewind(10)
       return
c
30     write(*,*)'corrupted data (or file too long)',filename(1),' row '
     . ,int(count+1)
       write(*,*)'the program will compute data up to this height'
       goto 40
c
11     format(x,f4.1,4(3x,f7.1),3x,f7.2,3(3x,f7.3))
12     format(f4.1,3x,4(x,f7.1),x,f6.2,x,f6.3,x,2(x,f7.3))
13     format(4x,A12,5x,i3,5x,i5)
       end
