C       THIS PROGRAM IS CALLED PREPON4FOR. IT IS LIKE PREPON3.FOR EXCEPT
C       IT HAS  BEEN ALTERED TO be a subroutine

C PROGRAM FOR ASSESSING PREVAILING WIND EROSION DIRECTION AND PREPONDERANCE
C OF WIND EROSION FORCES IN THE PREVAILING WIND EROSION DIRECTION
C DEFINITION:
C (1) U(I) IS THE MEAN WINDSPEED WITHIN THE ITH GOUP
C (2) F(I,J) IS A DURATION FACTOR WHICH IS EXPRESSED AS THE PERCENTAGE
C      OF THE TOTAL OBSERVATIONS THAT OCCUR IN THE JTH DIRECTION WITHIN
C     THE ITH SPEED GROUP
C (5) R(J) IS THE MAGNITUDE OF A WIND EROSION FORCE IN THE JTH DIRECTION
C (6) SUM IS THE SUM OF MAGNITUDE OF WIND EROSION FORCES
C (7) DELTR1 IS THE ANGLE OF ORIENTATION FOR OBTAINING THE MAXIMUM
C     PREPONDERNCE-PREVAILING WIND EROSION DIRECTION
C (8) CMAX IS THE MAXIMUM PREPONDERANCE
C (9) FPP IS THE POSITIVE PARALELL FORCE
C (10) FPN IS THE NEGATIVE PARALELL FORCE
C (11) M IS THE NUMBER OF DATA GROUPS, THIS MAY BE MONTHS OF THE YEAR
C (12) N IS THE NUMBER OF WIND SPEED GROUPS
C (13) CARD(I) GIVES ADDRESS DATA AND LOCATION OF DATA SOURCE

C*******THE VALUE OF N SHOULD BE SUBSTITUTED FOR THE NUMBER OF
C       INTEGER CONSTANTS TO BE TRANSMITTED IN FORMAT STATEMENT
C       NUMBERS 301 AND 200. *******************************************

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                   DATA ARRANGEMENT                            C

C (1) THE FIRST CARD: ONE INTEGER M INDICATES THE NUMBER OF DATA
C     GROUPS WHICH WE USE. (FORMAT I3) E.G. SALINA JAN, FEB, & MAR-
C     3 DATA GROUPS. 1 M CARD PER RUN.

C (2) THE FIRST CARD OF EACH DATA GROUP: ONE INTEGER N INDICATES THE
C     NUMBER OF WIND SPEED GROUPS IN THIS DATA GROUP. (FORMAT I3).

C (3) THE 2ND CARD OF EACH DATA GROUP: IS THE ADDRESS CARD - LOCATION
C     DATE AND RECORD.

C (4) THE 3TH CARD IS FIRST DATA CARD:  TOTAL OF 16 DATA CARDS. ONE
C     FOR EACH CARDINAL DIRECTION. (FORMAT NF5.1). DATA ARE IN PERCENT
C       WITH DECIMAL FORMAT.

C (5) AFTER THE 16 DATA CARDS HAVE A CARD THAT HAS THE MEAN WINDSPEED
C     IN MPH FOR EACH OF THE WINDSPEED GROUPS (FORMAT NF6.3).

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      subroutine prepon4(kk,f,ht,pwed,prepon,pospar)

C DECLARATION
      REAL F(16,19),U(20),UF(360),DF(360),R(16),
     1DR(360),AT,ACO,ASI,SUM,CR(360),DELTR1,
     2CMAX,D1,U1,U2,AT1,CO1,ACO1,FPP,FPN,
     *FPPA,FPNA,
     3DDR,CAT,ACAT,SAT,ASAT
     4,PWED(12),PREPON(12),POSPAR(12),HT
       integer kk,m,n,i,j,j1,k

c      LOGICAL FEXIST
c      CHARACTER RUNFIL*25, OUTFIL*25
C I/O FORMAT SEQUENCE
   30 FORMAT(I3)
   10 FORMAT(15A4)
   20 FORMAT(1X,15A4///)
   40 FORMAT(I3,F6.1)
   50 FORMAT(4I5)
   55 format(' prep',i3,5f7.1)
  100 FORMAT(I2)
  111 FORMAT(1H1)
  112 FORMAT(' ','DF(',I3,') IS LESS THAN OF EQUAL TO 0.0 .')
  200 FORMAT(4F6.3)
  300 FORMAT(4I5)
  301 FORMAT(16F5.1)
  400 FORMAT(1X,'J=1:N J=2:NNE J=3:NE J=4:ENE J=5:E'/1X,
     1'J=6:ESE J=7:SE J=8:SSE J=9:S J=10:SSW'/1X,
     2'J=11:SW J=12:WSW J=13:W J=14:WNW J=15:NW'/1X,
     3'J=16:NNW'///)
  500 FORMAT(1X,'MAGNITUDE OF WIND EROSION FORCES'/1X,
     1'FROM J=1 TO J=16'/)
  600 FORMAT(1X,'R(',I2,')= ',F10.3)
  700 FORMAT(1X,'SUM OF MAGNITUDE OF ALL WIND EROSION FORCES=',F6.1)
  800 FORMAT(1X,'ANGLE OF ORIENTATION = ',F5.1/1X,
     1'THE PREPONDERANCE R= ',F5.2/1X,
     2'POSITIVE PARALELL FORCE= ',F6.3/1X,
     3'NEGATIVE PARALELL FORCE= ',F6.3/////)
  900 FORMAT(1X,'OBTAINING THE MAXIMUM PREPONDERANCE')
  910 FORMAT(2X,'PARAMETER  JAN  FEB  MAR  APR  MAY  JUN  JUL  AUG  SEP
     1 OCT  NOV  DEC'/)
  950 FORMAT(F5.1)
C
C               OPEN FILES
C     OPEN(1,FILE='LUBBOCK.OUT',ACCESS='SEQUENTIAL',STATUS='OLD')
C     OPEN(2,FILE='PREPON.OUT',ACCESS='SEQUENTIAL',STATUS='UNKNOWN')
C     OPEN(3,FILE='PWED2.OUT',ACCESS='SEQUENTIAL',STATUS='UNKNOWN')
C
C
C                INPUT
C

C             ASSIGN APPROPRIATE VALUES FOR M AND N
C
      M=12
      N=19

C             ASSIGN APPROPRIATE WIND SPEED TO EACH OF THE N WIND SPEED GROUPS
C
      DO 60 I=1,14
      U(I)=I
  60  CONTINUE
      U(15)=17
      U(16)=22
      U(17)=27
      U(18)=32
      U(19)=37
C
c      DO 10000 KK=1,M
C
C READ IN THE WIND SPEED/FREQUENCY DATA AND ANEMOMETER HEIGHT
C
c      READ(5,950)HT

C CALCULATE THE WIND EROSION FORCES AND THEIR SUM

      SUM=0.0
      DO 1 J=1,16
      R(J)=0.0
      DO 2 I=1,N
    2 R(J)=R(J)+((U(I)*((10/HT)**.1429))**3)*F(J,I)
      SUM=SUM+R(J)
    1 CONTINUE
C CALCULATE THE PREPONDERANCE:
      DDR=ASIN(1.0)/90.0
      DO 3 K=1,360
      UF(K)=0.
      DF(K)=0.
      DR(K)=DDR*(K-1)
      DO 4 J=1,16
      J1=J-1
      AT=(J1*22.5*DDR)-DR(K)
      CAT=COS(AT)
      ACAT=ABS(CAT)
      SAT=SIN(AT)
      ASAT=ABS(SAT)
      ACO=ACAT*R(J)
      ASI=ASAT*R(J)
      UF(K)=UF(K)+ACO
      DF(K)=DF(K)+ASI
    4 CONTINUE
      IF (DF(K).EQ. 0.0)
     1WRITE(*,112) K
      IF (DF(K).EQ. 0.0)
     1GO TO 3
      CR(K)=UF(K)/DF(K)
    3 CONTINUE
C FIND THE ANGLE OF ORIENTATION, POSITIVE PARALELL FORCE AND NEGATIVE
C PARALELL FORCE FOR OBTAINING THE MAXIMUM PREPONDERANCE:
      CMAX=CR(1)
      DELTR1=DR(1)
      DO 6 K=1,360
      IF(CR(K).LT.CMAX) GO TO 6
      CMAX=CR(K)
      PREPON(KK)=CMAX
      DELTR1=DR(K)
    6 CONTINUE
      if (prepon(kk) .gt. 99.9) prepon(kk) = 99.9
      U1=0.
      U2=0.
      D1=0.
      DO 7 J=1,16
      J1=J-1
      AT1=J1*22.5*DDR-DELTR1
      CO1=R(J)*COS(AT1)
      ACO1=ABS(CO1)
      D1=D1+ACO1
      IF(CO1.GE.0.) GO TO 11
      U2=U2+CO1
      GO TO 7
   11 U1=U1+CO1
    7 CONTINUE
      FPP=U1/D1
      POSPAR(KK)=FPP
      FPN=U2/D1
      FPNA=ABS(FPN)
      FPPA=FPP
      DELTR1=DELTR1/DDR
      IF(FPNA.LE.FPPA) GO TO 1978
      FPP=FPNA
      POSPAR(KK)=FPP
      FPN=-FPPA
      DELTR1=DELTR1-180.0
      IF(DELTR1.GE.0.) GO TO 1978
      DELTR1 = DELTR1 + 360.
C PRINT OUT THE RESULTS:
 1978 PWED(KK)=DELTR1

      return
      END