//WGENPAR  JOB (M461,54-D,S40,10,R1),'WGEN PARAM'
//*TAMU P=5
//STEP  EXEC WATFIV,REGION=256K
//FT02F001  DD  DSN=USR.M461.R1.AUST02,DISP=SHARE
//SYSIN  DD  DATA
//$OPTIONS
C     ******************************************************************
C     *                                                                *
C     *                 DEE ALLEN WRIGHT                               *
C     *                 COMPUTER PROGRAMMER                            *
C     *                 FEBRUARY 25, 1983                              *
C     *                                                                *
C     ******************************************************************
      DIMENSION TMAX(20,365),TMIN(20,365),RAIN(20,365),RAD(20,365)
      DIMENSION AA(20)
      DIMENSION RC(365)
      DIMENSION XDATA(30),YDATA(4,12)
      DO1I=1,30
      XDATA(I)=0.0
1     CONTINUE
      DO2I=1,4
      DO2J=1,12
      YDATA(I,J)=0.0
2     CONTINUE
C*******************************************************************************
C
C-----INPUT #1  TITLE OF DATA SET (20A4)
C
C*******************************************************************************
      READ(2,100)(AA(I),I=1,20)
100   FORMAT(20A4)
      WRITE(6,101) (AA(I),I=1,20)
101   FORMAT('1',//,30X,'DATA IS -----',20A4,//)
C*******************************************************************************
C
C-----INPUT #2  NUMBER OF YEARS AND LATITUDE
C
C*******************************************************************************
      READ(2,102) NYRS,ALAT
102   FORMAT(I10,F10.0)
      XYRS=NYRS
      XLAT=ALAT*6.2832/360.
      DO 6 I = 1,365
      XI = I
      SD = 0.4102*SIN(0.0172*(XI-80.25))
      CH = -TAN(XLAT)*TAN(SD)
      IF (CH .GT. 1.0) H = 0.
      IF (CH .GT. 1.0) GO TO 5
      IF(CH .LT. -1.0) H = 3.1416
      IF(CH .LT. -1.0) GO TO 5
      H = ARCOS(CH)
5     DD=1.0+0.0335*SIN(0.0172 *(XI+88.2))
      RC(I)=889.2305*DD*((H*SIN(XLAT)*SIN(SD))+(COS(XLAT)*COS(SD)*SIN(H)
     *))
      RC(I)=RC(I)*0.80
6     CONTINUE
      DO 7 I = 1,NYRS
      DO 7 J = 1,365
C*******************************************************************************
C
C----- INPUT #3  DAY, MO, YEAR, RADIATION, MAX TEMP, MIN TEMP, RAINFALL
C
C*******************************************************************************
8     READ(2,103) IYR,IMO,IDA,V1,V2,V3,V4
      IF(IMO .EQ. 2 .AND. IDA .EQ. 29) GO TO 8
      TMAX(I,J) = V1
      TMIN(I,J) = V2
      RAIN(I,J) = V3
      RAD(I,J) = V4
      IF(RAD(I,J) .GT. RC(J)) RAD(I,J)=RC(J)
 7    CONTINUE
 103  FORMAT(10X,4F10.2)
C103  FORMAT(2X,I4,2I3,F8.2,2F8.0,8X,F8.0)
      WRITE(6,104)
104   FORMAT('1',///,5X,'MAXIMUM TEMPERATURE',/)
      CALL MSD(NYRS,TMAX,RAIN,1,XDATA)
      WRITE(6,105)
105   FORMAT('1',///5X,'MINIMUM TEMPERATURE',/)
      CALL MSD(NYRS,TMIN,RAIN,2,XDATA)
      WRITE(6,106)
106   FORMAT('1',///,5X,'SOLAR RADIATION ',/)
C     CALL MSD(NYRS,RAD,RAIN,3,XDATA)
      WRITE(6,107)
107   FORMAT('1',///,5X,'PRECIPITATION')
      CALL PPRAIN(RAIN,NYRS,YDATA)
C--PHASE ANGLE SHIFTED BY 180 DEGREES
C--AND SIGNS CHANGED ON AMPLITUDES
      XDATA(02)=XDATA(02)*(-1.0)
      XDATA(04)=XDATA(04)*(-1.0)
      XDATA(10)=XDATA(10)*(-1.0)
      XDATA(12)=XDATA(12)*(-1.0)
      XDATA(14)=XDATA(14)*(-1.0)
      WRITE(6,101) (AA(I),I=1,20)
      WRITE(6,707)
707   FORMAT(///,20X,'INPUT CARDS FOR THE WEATHER GENERATOR ARE AS FOLLO
     *WS ----------',////)
      WRITE(6,403)
403   FORMAT(///)
      WRITE(6,513) (YDATA(1,J),J=1,12)
      WRITE(6,514) (YDATA(2,J),J=1,12)
      WRITE(6,515) (YDATA(3,J),J=1,12)
      WRITE(6,516) (YDATA(4,J),J=1,12)
513   FORMAT(5X,'INPUT #  3 ---- P(W/W) ---- ',12F6.3)
514   FORMAT(5X,'INPUT #  4 ---- P(W/D) ---- ',12F6.3)
515   FORMAT(5X,'INPUT #  5 ---- ALPHA  ---- ',12F6.3)
516   FORMAT(5X,'INPUT #  6 ---- BETA   ---- ',12F6.3)
      WRITE(6,400)
400   FORMAT(////)
      WRITE(6,701)
      WRITE(6,501) XDATA(01)
      WRITE(6,502) XDATA(02)
      WRITE(6,503) XDATA(03)
      WRITE(6,504) XDATA(04)
701   FORMAT(5X,'INPUT #  7 ----',/)
501   FORMAT(15X,' 1    TXMD    ----  ',F10.3)
502   FORMAT(15X,' 2    ATX     ----  ',F10.3)
503   FORMAT(15X,' 3    CVTX    ----  ',F10.3)
504   FORMAT(15X,' 4    ACVTX   ----  ',F10.3)
      WRITE(6,702)
702   FORMAT(//,5X,'INPUT #  8 ----',/)
      WRITE(6,505) XDATA(05)
505   FORMAT(15X,' 5    TXMW    ----  ',F10.3)
      WRITE(6,703)
703   FORMAT(//,5X,'INPUT #  9 ----',/)
      WRITE(6,506) XDATA(09)
      WRITE(6,507) XDATA(10)
      WRITE(6,508) XDATA(11)
      WRITE(6,509) XDATA(12)
506   FORMAT(15X,' 6    TN      ----  ',F10.3)
507   FORMAT(15X,' 7    ATN     ----  ',F10.3)
508   FORMAT(15X,' 8    CVTN    ----  ',F10.3)
509   FORMAT(15X,' 9    ACVTN   ----  ',F10.3)
      WRITE(6,704)
704   FORMAT(//,5X,'INPUT # 10 ----',/)
      WRITE(6,510) XDATA(13)
      WRITE(6,511) XDATA(14)
510   FORMAT(15X,'10    RMD     ----  ',F10.3)
511   FORMAT(15X,'11    AR      ----  ',F10.3)
      WRITE(6,705)
705   FORMAT(//,5X,'INPUT # 11 ----',/)
      WRITE(6,512) XDATA(17)
512   FORMAT(15X,'12    RMW     ----  ',F10.3)
      WRITE(6,600)
600   FORMAT('1')
      STOP
      END
      SUBROUTINE FOUR(XM,SD,CV,XDATA)
      DIMENSION XM(13),SD(13)
      DIMENSION CV(13)
      DIMENSION XDATA(30)
      DATA JCT /0/
      S = 0.
      S1 = 0.
      S2 = 0.
      WRITE(6,200)
200   FORMAT(//,33X,'    PERIOD   MEAN    STD DEV     CV')
      DO 10 I = 1,13
      WRITE(6,201)I,XM(I),SD(I),CV(I)
201   FORMAT(30X,I10,3F10.2)
      S = S + XM(I)
      S1 = S1 + SD(I)
      S2 = S2 + CV(I)
10    CONTINUE
      XBAR=S/13.
      XBAR1=S1/13.
      XBAR2=S2/13.
      SUMA = 0.
      SUMB = 0.
      SUMA1 = 0.
      SUMB1 = 0.
      SUMA2=0.
      SUMB2 = 0.
      DO 15 K = 1,13
      XK = K
      SUMA=SUMA+(XM(K)-XBAR)*COS(6.2832*XK/13.)
      SUMA1=SUMA1+(SD(K)-XBAR1)*COS(6.2832*XK/13.)
      SUMA2=SUMA2+(CV(K)-XBAR2)*COS(6.2832*XK/13.)
      SUMB=SUMB+(XM(K)-XBAR)*SIN(6.2832*XK/13.)
      SUMB1=SUMB1+(SD(K)-XBAR1)*SIN(6.2832*XK/13.)
      SUMB2=SUMB2+(CV(K)-XBAR2)*SIN(6.2832*XK/13.)
15    CONTINUE
      A = SUMA*(2./13.)
      A1= SUMA1*(2./13.)
      A2 = SUMA2*(2./13.)
      B = SUMB*(2./13.)
      B1 = SUMB1*(2./13.)
      B2 = SUMB2*(2./13.)
      T = ATAN(-B/A)
      T1= ATAN(-B1/A1)
      T2=ATAN(-B2/A2)
      C = A/COS(T)
      C1 = A1/COS(T1)
      C2 = A2/COS(T2)
      WRITE(6,100)
100   FORMAT(/,15X,'FOURIER COEFFICIENTS--MEAN')
      WRITE(6,101) XBAR,C,T
101   FORMAT(15X,'MEAN =',F10.4,5X,'AMPLITUDE =',F10.4,5X,'PHASE =',
     *F10.4)
      JCT=JCT+1
      XDATA(JCT)=XBAR
      JCT=JCT+1
      XDATA(JCT)=C
      WRITE(6,102)
102   FORMAT(/,15X,'FOURIER COEFFICIENTS--STD. DEV.')
      WRITE(6,101)XBAR1,C1,T1
      WRITE(6,103)
103   FORMAT(/,15X,'FOURIER COEEFICIENTS--CV')
      WRITE(6,101) XBAR2,C2,T2
      JCT=JCT+1
      XDATA(JCT)=XBAR2
      JCT=JCT+1
      XDATA(JCT)=C2
      RETURN
      END
      SUBROUTINE MSD(NYRS,W,RAIN,ID,XDATA)
      DIMENSION W(20,365),RAIN(20,365),XM(13),XM1(13),SD(13),SD1(13)
      DIMENSION CX(13),CX1(13)
      DIMENSION XDATA(30)
      DO 20 I = 1,13
      NF=I*28
      NI=NF-27
      XN = 0.
      XN1 = 0.
      SUM = 0.
      SUM1 = 0.
      SS = 0.
      SS1 = 0.
      DO 15 JD=NI,NF
      DO 15 JY = 1,NYRS
      IF(ID .EQ. 2) GO TO 11
      IF(RAIN(JY,JD))11,11,12
11    CONTINUE
      XN = XN + 1.
      SUM = SUM+W(JY,JD)
      SS = SS + (W(JY,JD)*W(JY,JD))
      GO TO 15
12     CONTINUE
      XN1 = XN1 + 1.
      SUM1 = SUM1 + W(JY,JD)
      SS1=SS1+(W(JY,JD)*W(JY,JD))
15    CONTINUE
      IF(XN .LE. 2.) XM(I) = 0.
      IF(XN .LE. 2.) SD(I) = 0.
      IF(XN .LE. 2.) CX(I) = 0.
      IF(XN .LE. 2.) GO TO 400
      XM(I) = SUM / XN
      SD(I) = SQRT((SS-SUM*SUM/XN)/(XN-1.))
      IF(XM(I) .LT. 0.001) XM(I) = 0.001
      CX(I) = SD(I) / XM(I)
400   CONTINUE
      IF(ID .EQ. 2) GO TO 20
      IF(XN1 .LE. 2.) XM1(I)=0.
      IF(XN1 .LE. 2.) SD1(I) = 0.
      IF(XN1 .LE. 2.) CX1(I) = 0.
      IF(XN1 .LE. 2.) GO TO 500
      XM1(I) = SUM1 / XN1
      SD1(I)=SQRT((SS1-SUM1*SUM1/XN1)/(XN1-1.))
      IF(XM1(I) .LT. 0.001) XM1(I) = 0.001
      CX1(I)=SD1(I)/XM1(I)
500   CONTINUE
20    CONTINUE
      IF(ID .EQ. 2) GO TO 25
      WRITE(6,100)
100   FORMAT(10X,'DRY DAYS')
      CALL FOUR(XM,SD,CX,XDATA)
      WRITE(6,101)
101   FORMAT(/,10X,'WET DAYS')
      CALL FOUR(XM1,SD1,CX1,XDATA)
      GO TO 30
25    WRITE(6,102)
102   FORMAT(10X,'WET AND DRY DAYS')
      CALL FOUR(XM,SD,CX,XDATA)
30    CONTINUE
      RETURN
      END
      SUBROUTINE PPRAIN(XRAIN,NYR,YDATA)
      DIMENSION XRAIN(20,365)
      DIMENSION NWD(12),NDD(12),NDW(12),NWW(12)
      DIMENSION SUM(12),SUM2(12),SUM3(12)
      DIMENSION SL(12),PWW(12),PWD(12),RBAR(12)
      DIMENSION ALPHA(12),BETA(12)
      DIMENSION NW(12),IC(12),SUML(12)
      DIMENSION RLBAR(12),AL2(12),BE2(12),DATE(12)
      DIMENSION PPPW(12),ND(12)
      DIMENSION YDATA(4,12)
      CHARACTER *36 A(2)
      DATA DATE /'JAN.','FEB.','MAR.','APR.','MAY ','JUNE',
     *           'JULY','AUG.','SEP.','OCT.','NOV.','DEC.'/
      DATA A(1) /'                                    '/
      DATA A(2) /'NOT ENOUGH DATA TO DEFINE PARAMETERS'/
      DO 10 I =1,12
      ND(I) = 0
      PPPW(I) = 0.
      NWD(I) = 0
      NWW(I) = 0
      NDD(I) = 0
      NDW(I) =0
      NW(I) = 0
      SL(I) = 0.
      SUML(I) = 0.
      SUM(I) = 0.
      SUM2(I) = 0.
      PWW(I) = 0.
      PWD(I) = 0.
      ALPHA(I) = 0.
      BETA(I) = 0.
10    SUM3(I) = 0.
      XYR=NYR
      RIM1 = 0.
      DO 20 J = 1,NYR
      DO 30 K = 1,365
      IF(K .GE. 001 .AND. K .LE. 031) MO = 1
      IF(K .GE. 032 .AND. K .LE. 059) MO = 2
      IF(K .GE. 060 .AND. K .LE. 090) MO = 3
      IF(K .GE. 091 .AND. K .LE. 120) MO = 4
      IF(K .GE. 121 .AND. K .LE. 151) MO = 5
      IF(K .GE. 152 .AND. K .LE. 181) MO = 6
      IF(K .GE. 182 .AND. K .LE. 212) MO = 7
      IF(K .GE. 213 .AND. K .LE. 243) MO = 8
      IF(K .GE. 244 .AND. K .LE. 273) MO = 9
      IF(K .GE. 274 .AND. K .LE. 304) MO = 10
      IF(K .GE. 305 .AND. K .LE. 334) MO = 11
      IF(K .GE. 335 .AND. K .LE. 365) MO = 12
      RAIN=XRAIN(J,K)
      IF(RAIN .GT. 0.00) NW(MO)=NW(MO)+1
      ND(MO)=ND(MO)+1
      IF(RAIN) 5,5,3
  3   IF(RIM1)2,2,4
2     NWD(MO)=NWD(MO)+1
      GO TO 6
4     NWW(MO)=NWW(MO)+1
6     CONTINUE
      SUML(MO)=SUML(MO)+ALOG(RAIN)
      SUM(MO)=SUM(MO)+RAIN
      SUM2(MO)=SUM2(MO) + RAIN * RAIN
      SUM3(MO)=SUM3(MO)+RAIN*RAIN*RAIN
      SL(MO) = SL(MO)+ALOG(RAIN)
      GO TO 9
5     IF(RIM1) 7,7,8
 7    NDD(MO)=NDD(MO)+1
      GO TO 9
8     NDW(MO)=NDW(MO)+1
9     RIM1 = RAIN
30    CONTINUE
20    CONTINUE
      DO 120 I = 1,12
      XXND=ND(I)
      YYNW=NW(I)
      PPPW(I) = YYNW/XXND
      III=1
      IF(NW(I) .LT. 3) III=2
      IC(I) = III
      IF(NW(I) .LT. 3) GO TO 120
      XNWW=NWW(I)
      XNWD=NWD(I)
      XXNW=NWW(I)+NDW(I)
      XND=NDD(I)+NWD(I)
      XNW=NW(I)
      PWW(I)=XNWW/XXNW
      PWD(I)=XNWD/XND
      RBAR(I)=SUM(I)/XNW
      RLBAR(I)=SUML(I)/XNW
      Y=ALOG(RBAR(I))-RLBAR(I)
      ANUM=8.898919+9.05995*Y+0.9775373*Y*Y
      ADOM=Y*(17.79728+11.968477*Y+Y*Y)
      ALPHA2=ANUM/ADOM
      IF(ALPHA2 .GE. 1.0) ALPHA2=0.998
      BETA2=RBAR(I)/ALPHA2
      ALPHA(I)=ALPHA2
      BETA(I)=BETA2
 120  CONTINUE
      WRITE(6,201)
201   FORMAT(///,8X,'--MARKOV CHAIN--                -GAMMA DIST-',/,
     *1X,' MONTH  P(W/W)    P(W/D)                ALPHA   BETA',/)
      DO 130 I = 1,12
      WRITE(6,202)DATE(I),PWW(I),PWD(I),ALPHA(I),BETA(I),A(IC(I))
202   FORMAT(1X,A4,F8.3,F10.3,11X,F11.3,F7.3,5X,A36)
130   CONTINUE
      DO400J=1,12
      YDATA(1,J)=PWW(J)
      YDATA(2,J)=PWD(J)
      YDATA(3,J)=ALPHA(J)
      YDATA(4,J)=BETA(J)
400   CONTINUE
      RETURN
      END
//$DATA
