program wind_genf 

      include 'wind_db.fi'

c-------------------1. declear variable and format --------------------
      integer lnblnk, julday

c     for file units (read and write)
      integer i_unit	!input file unit (wind_gen database record)
      integer o_unit	!output file unit (wind_gen output file)
      integer e_unit	!error output file unit 

c     output filename
      character file_wdb*128, file_out*128, file_out1*128

c     desired wban station no.
      integer station_no

      character*1 answer

c     used for loop
      integer i, m, j

c     julian start and stop days
      integer jstart, jstop

c     random number
      real rnd1, rnd2, rand

c     max, min, and representative wind speed
      real umax, umin, urep

c     a term
      real a

      ! directional frequency for each month (computed)
      real dirf(0:NDIR,0:NMONTHS-1)

c     beginning year of simulation and number of yeras to simulate
      integer begin_year, no_years
	  integer y

c     the month for scale access, calendar year and day
      integer month, cal_year, cal_day

c     the direction of the wind
      integer idir      ! index for wind directions
      real    wdir      ! actual computed wind direction in degrees

c     beginning calendar year of simulation
      begin_year = 1

c     number of years to simulate
      no_years = 1

c     Set unit numbers for input and output file devices.
c     (stdin = 5, stdout = 6, stderr = 0)
      i_unit = 8	!stdin		!can't use stdin until station no. is obtained differently
      o_unit = 6	!stdout
      e_unit = 0	!stderr

c     for the database header
  100 format (a,1x,i5,1x,a,a,1x,a)
  110 format (i2,1x,i2,1x,a,1x,i3,1x,i2,1x,a,1x,i2,1x,a,1x,a,1x,a)
  200 format (12(1x,f4.1))  !for dir(i,m) 
  300 format (12(2x, f4.2)) !for scale(i,m) and shape(i,m) 
  400 format (12(2x,f3.1))  !for ratio(m) 
  500 format (12(i2, 3x))   !for hr_max(m) 

c     for the output header
  600 format (1x,'station: ',a,', ',a,2x,a)
  610 format (1x,'lat:  ', i3, 'd  ', i2,'m ', a,'  lon: ',
     +        i3,'d  ',i2,'m ',a)  
  620 format (1x, 'period:  ',a,'-',a,'  el: ',i4, 'm') 
  630 format (1x, 'day mo year  dir   umax  umin hrmax')
  631 format (1x, 'day mo year  dir   umax  umin hrmax',
     *    '  urep   a   shape scale idir rand1 rand2')
  640 format (1x, '             deg    m/s   m/s')
  650 format ('-------------------------------------')

c     for the output data
  700 format (i3,i3,1x,i5,f6.1,f6.1,f6.1,i5)

c     for long output format 
  800 format (i3,i3,1x,i5,f6.1,f6.1,f6.1,i5,
     *   f6.1,f7.3,f6.2,f6.2,i3,f6.3,f6.3)


c     Open input file (if changed from stdin).
	  if (i_unit .ne. 5) then
         file_wdb = 'wind_gen.wdb' ! default wind_gen db file name
         write(0,*) 'Enter the name of wind_gen database input file:'
         read*, file_wdb          
         open (i_unit, file=file_wdb)
      else
         file_wdb = 'stdin'
      endif

      write(0,*) 'Enter the station no:'
      read*, station_no          

      write(0,*) 'Do we want "long" output format? (y/n)'
      read*, answer

c     Open output file (if changed from stdout).
      if (o_unit .ne. 6) then
        print*, 'Enter the name of daily-wind output file:'
        read*, file_out1   
        open(o_unit, file = file_out1)
      endif

      call read_wdb(i_unit, station_no)

      if (answer .eq. 'y' .or. answer .eq. 'Y') then
          call disp_wdb(o_unit)		! part of long output format
          !call write_wdb(o_unit)	! to verify read_wdb()
      endif


c------------------ 4. do the daily wind speed simulation ---------------

! print output header
      write(o_unit,600)
     *      st_city(1:lnblnk(st_city)), st_state, st_country
      write(o_unit,610)
     *      lati_deg, lati_min, lati_hem,
     *      long_deg, long_min, long_hem
      write(o_unit,620) s_date, e_date, elevation
      if (answer .eq. 'y' .or. answer .eq. 'Y') then
          write(o_unit,631)
      else
          write(o_unit,630)
      endif
      write(o_unit,640)
      write(o_unit,650)

         
      jstart = julday(1,1,begin_year)
	  y = begin_year+no_years-1
      jstop  = julday(31,12,y)
      rnd1 = rand(1)      !initialize random number generator

      call get_dirf(dirf) !initialize dirf() array

      do 80 j=jstart, jstop  !for each day....

         call caldat(j, cal_day, month, cal_year)
         rnd1 = rand(0)     !get random no for computing wind direction
         rnd2 = rand(0)     !get random no for computing wind speed 


         call day_wind( rnd1, rnd2, month, dirf(0,0),
     +      idir, wdir, a,
     +      urep, umax, umin) 
       
         if (answer .eq. 'y' .or. answer .eq. 'Y') then
            if (idir .eq. 16) then
               write(o_unit,800) cal_day, month, cal_year,
     +           wdir, umax, umin, -1,
     +           urep, 0.0, -1.0, -1.0,
     +           idir,rnd1,rnd2 
            else
               write(o_unit,800) cal_day, month, cal_year,
     +           wdir, umax, umin, hr_max(month-1),
     +           urep,a,shape(idir,month-1),scale(idir,month-1),
     +           idir,rnd1,rnd2 
            endif
         else
           write(o_unit,700) cal_day, month, cal_year,
     +           wdir, umax, umin, hr_max(month-1) 
         endif

   80    continue         

      if (i_unit .ne. 5) then
          close(i_unit)
      endif
      if (o_unit .ne. 6) then
          close(o_unit)
      endif

      stop
      end