subroutine windg
c
c     + + + PURPOSE + + +
c     Simulates daily average wind velocity using mean, standard deviation,
c     and skew coefficient of wind speed.  Wind direction is selected from
c     a uniform distribution.
c
c     + + + COMMON BLOCKS + + +
      include 'cbk1.inc'
c      read: dir,pi2
c      modify: wvl,wv,th
c     dir(i,j)    - Cumulative % time (fraction) from dir1, dir1+dir2, ...
c                          dim1: month
c                          dim2: compass direction
c     pi2         - Pi * 2; ie, a full rotation.
c     wvl(i,j,k)  - array of wind paramters where:
c                     i - ith direction (1 - north  - 16 nnw)
c                     j - parameters (1 - 4)
c                          1 - % time from direction i
c                          2 - mean speed from direction i
c                          3 - standard deviation of speed from direction i
c                          4 - skew coeficient of speed from direction i
c                     k - month (1=Jan, 2=Feb...)
c     wv          - Wind Velocity (m/sec)
c     th          - Wind Direction (radians from North)
c
      include 'cbk3.inc'
c -- XXX -- Note that 'j' is a counter in a loop! -- CRM 10/18/99
c      modify: j
c
      include 'cbk4.inc'
c      read: mo
c     mo          - The current month (1=Jan, 2=Feb...).
c
      include 'cbk7.inc'
c      read: k6,k8
c      modify: v9
c     k1 ... k9   - Seeds for random number generation.
c     v1 ... v11  - Random numbers used to generate various daily values.
c
CC    include 'ctap.inc'
CC    include 'ctap2.inc'
c       write: tap6
c
      include 'crandom.inc'
c        read: fxx
c
c     + + + FUNCTION DECLARATIONS + + +
      real dstn1
c
c     + + + LOCAL VARIABLES + + +
      real    tap6, xlv, r6, v10, xj1, g
      integer j1, ndflag
c     + + + END SPECIFICATIONS + + +
c
      fx=fxx(dax)
      j = 0
 1    continue
        j = j + 1
        j1=j-1
        ndflag = 0
        if (dir(mo,j).gt.fx) ndflag = 2
      if(ndflag.eq.0 .and. j.lt.16) goto 1
c
c       calm condition found
c
c
c **** L1 IF ****
      if(ndflag.eq.0) then
        wv=0.0
        th=0.0
c **** L1 ELSE ****
      else 
c
c       Wind direction calculations
c
        if (j.eq.1) then
          g=fx/dir(mo,j)
        else
          g=(fx-dir(mo,j1))/(dir(mo,j)-dir(mo,j1))
        endif
        xj1=j1
        th=pi2*(g+xj1-.5)/16.
        if (th.lt.0.) th=pi2+th
c
c         Wind speed calculations
c
        v10=v10x(dax)
        if(wvl(j,4,mo).eq.0.0) wvl(j,4,mo)=0.01
        r6=wvl(j,4,mo)/6.0
        xlv=(dstn1(v9,v10)-r6)*r6+1.0
        xlv=(xlv**3-1.0)*2.0/wvl(j,4,mo)
        v9=v10
	tap6=xlv
        wv=xlv*wvl(j,3,mo)+wvl(j,2,mo)
        if (wv.lt.0.) wv=.1
c **** L1 ENDIF ****
      endif
c
      return
      end