!**************************************
!
!     1 ¼Þ¹ÞÝ ¶¼®³ ͯÝÄÞ³ ¹²»Ý
!
!       ´Ý¼­³ ÓÝÀÞ² 19
!
!**************************************
!
      DIMENSION X(60),DX(60),H(60),Z(60),B(60),US(60),TS(60),ZZ(60),QSU(60),C(60),PTIM(6),OT(6)
      REAL IE(60)
!
      DATA NJ,Q,H0,SN,D,S/51,500.,12.,0.02,0.0001,1.65/
      DATA G,RAMDA,EPS/9.8,0.4,0.001/
      DATA IDT1,IDT2,ETIM/50,100,86400./
      DATA OT/0.,1.,4.,8.,12.,24./
!----------------------------------------------------------
      TIM=0.0
      ITM=1
      DO 100 I=2,10
         PTIM(I)=OT(I)*3600.
  100 CONTINUE
!----------------------------------------------------------
      DO 110 J=1,NJ
         X(J)=200.*(J-1)
         B(J)=100.
         Z(J)=1./700.*(J-1)*200.
         ZZ(J)=Z(J)
  110 CONTINUE
      DO 120 J=2,NJ
         DX(J)=X(J)-X(J-1)
!      WRITE(*,'(A,I5,A,F8.1,A,F8.3)')' J=',J,' DX=',DX(J),' Z=',Z(J)
  120 CONTINUE
!**********************************************************
!
      CALL CALWF(S,G,D,WF)
!
!----------------------------------------------------------
!
 1000 CONTINUE
!
!----------------------------------------------------------
!
      CALL CALH(H,Z,B,Q,SN,H0,G,DX,EPS,IE,NJ)
           DO 200 J=1,NJ
              TS(J)=H(J)*IE(J)/(S*D)
              US(J)=SQRT(G*H(J)*IE(J))
!      WRITE(*,'(A,I5,A,F8.3,A,F8.3)') ' J=',J,' H=',H(J),' TS=',TS(J)
  200      CONTINUE
      CALL CALIWA(S,G,D,USC,TSC)
      CALL CALQSU(S,G,D,NJ,US,TS,TSC,WF,QSU)
      CALL CALC(NJ,X,B,QSU,WF,C,Q,H,IDT1,IDT2)
      CALL POUT(Z,ZZ,H,TS,US,QSU,C,NJ,TIM,PTIM,ITM)
      CALL CALDZ(NJ,RAMDA,IDT2,WF,Z,QSU,C)
!
!**********************************************************
!
      TIM=TIM+IDT2
      IF (TIM.LE.ETIM) GO TO 1000
!
!----------------------------------------------------------
      STOP
      END
!     -----WF É ¹²»Ý-----
      SUBROUTINE CALWF(S,G,D,WF)
      DD=D*100.
      A=36.*0.01**2./(S*G*100.)
      DDD=DD**3.
      WF=(SQRT(2./3.+A/DDD)-SQRT(A/DDD))*SQRT(S*G*100.*DD)/100.
      RETURN
      END
!     -----USC,TSC É ¹²»Ý-----
      SUBROUTINE CALIWA(S,G,D,USC,TSC)
!     ***  IWAGAKI METHOD  ***
      DD=D*100.
      IF (DD.GE.0.303) USC2=80.9*DD
      IF ((0.188.LE.DD).AND.(DD.LT.0.303)) USC2=134.6*DD**(31./22.)
      IF ((0.0565.LE.DD).AND.(DD.LT.0.118)) USC2=55.*DD
      IF ((0.0065.LE.DD).AND.(DD.LT.0.0565)) USC2=8.41*DD**(11./32.)
      IF (DD.LT.0.0065) USC2=226.*DD
      USC=USC2**(1./2.)/100.
      TSC=USC**2./(S*G*D)
      RETURN
      END
!     ----- H É ¹²»Ý -----
      SUBROUTINE CALH(H,Z,B,Q,SN,H0,G,DX,EPS,IE,NJ)
      DIMENSION H(60),Z(60),B(60),DX(60)
      REAL IE(60)
      QQ=Q*Q
      SNN=SN*SN
      H(1)=H0-Z(1)
      DO 10 J=1,NJ
      JJ=J-1
      IF (J.EQ.1) GO TO 20
      HH=H(JJ)**2
      BB=B(JJ)**2
      H3=H(JJ)**(10./3.)
      FD=H(JJ)+Z(JJ)+QQ/(2.*G*BB*HH)+SNN*QQ*DX(J)/(2.*BB*H3)
      H(J)=H(JJ)
   30 CONTINUE
      HH=H(J)**2
      BB=B(J)**2
      H3=H(J)**(10./3.)
      FU=H(J)+Z(J)+QQ/(2.*G*BB*HH)-SNN*QQ*DX(J)/(2.*BB*H3)
      FH=FU-FD
      IF (ABS(FH).GT.EPS) THEN
      DFDH1=QQ/(G*BB*H(J)**3)
      DFDH2=(5./3.)*SNN*QQ*DX(J)/(BB*H(J)**(13./3.))
      DFDH=1.-DFDH1+DFDH2
      H(J)=H(J)-FH/DFDH
      GO TO 30
      END IF
   20 CONTINUE
      IE(J)=SNN*QQ/(B(J)*B(J)*H(J)**(10./3.))
   10 CONTINUE
      RETURN
      END
!     ----- QSU É ¹²»Ý -----
      SUBROUTINE CALQSU(S,G,D,NJ,US,TS,TSC,WF,QSU)
      DIMENSION US(60),QSU(60),TS(60)
      DATA A1,A2,A3/0.436184,-0.120168,0.937298/
      DATA AS,BS,PK/0.14,0.143,0.008/
      DATA RW,RS,ET/1.0,2.65,0.5/
      PI=3.141592
      SSS=(S*G*D)**0.5
      DO 10 J=1,NJ
      AD=BS/TS(J)-1./ET
      ADD=ABS(AD)
      T=1./(1+0.33267*2**0.5*ADD)
      D1=1./(2.*PI**0.5)*EXP(-ADD*ADD)
      D2=1./(2.*PI)**0.5*EXP(-ADD*ADD)*(A1*T+A2*T*T+A3*T**3.)
      IF (AD.LT.0.0) D2=1-D2
      OMEGA=TS(J)/BS*D1/D2+TS(J)/BS/ET-1.
      QSU(J)=PK*(AS*RW/RS*OMEGA*SSS/TS(J)**0.5-WF)
      IF (QSU(J).LT.0.) QSU(J)=0.0
      IF (US(J)/WF.LT.1.08) QSU(J)=0.0
      IF (TS(J).LT.TSC) QSU(J)=0.0
   10 CONTINUE
      RETURN
      END
!     ----- C É ¹²»Ý -----
      SUBROUTINE CALC(NJ,X,B,QSU,WF,C,Q,H,IDT1,IDT2)
      DIMENSION X(60),B(60),QSU(60),C(60),H(60),CC(60),CA(60)
      INTEGER TIM
      CC(NJ)=QSU(NJ)/WF
      DO 100 J=1,NJ
         C(J)=CC(J)
  100 CONTINUE
      DO 200 TIM=IDT1,IDT2,IDT1
        DO 300 J=1,NJ
           CA(J)=CC(J)
  300   CONTINUE
        DO 400 J=NJ-1,1,-1
           DXB=(X(J+1)-X(J))*B(J)
           CC(J)=CA(J)+IDT1/H(J)*(QSU(J)-WF*CA(J)-Q*(CA(J)-CA(J+1))/DXB)
  400  CONTINUE
  200 CONTINUE
      RETURN
      END
!     ----- PRINT OUT -----
      SUBROUTINE POUT(Z,ZZ,H,TS,US,QSU,C,NJ,TIM,PTIM,ITM)
      DIMENSION Z(60),ZZ(60),H(60),TS(60),US(60),QSU(60),C(60),HH(60),PTIM(6)
!      OPEN(6,FILE='PRN')
      IF (TIM.LT.PTIM(ITM)) GO TO 900
      ITM=ITM+1
      WRITE(6,630) TIM/3600.
  630 FORMAT(1H ,3X,'TIME=',F5.0,' HR'//)
      WRITE(6,600)
  600 FORMAT(1H ,3X,'NO',3X,'Z (M)',6X,'DZ (M)',4X,'HH (M)',4X,'H (M)',7X,'T*',4X,'U* (CM/S)',1X,'QSU (CM/S)',4X,'C'/)
      DO 700 J=1,NJ
         HH(J)=H(J)+Z(J)
         DZ=Z(J)-ZZ(J)
         WRITE(6,610) J,Z(J),DZ,HH(J),H(J),TS(J),US(J)*100., QSU(J)*100.,C(J)
  610    FORMAT(1H ,I5,8F10.5)
  700 CONTINUE
      WRITE(6,620)
  620 FORMAT(1H1)
  900 RETURN
      END
!     ----- DZ É ¹²»Ý -----
      SUBROUTINE CALDZ(NJ,RAMDA,IDT2,WF,Z,QSU,C)
      DIMENSION Z(60),QSU(60),C(60)
      DO 100 J=1,NJ
         DZ=1./(1-RAMDA)*(WF*C(J)-QSU(J))*IDT2
         Z(J)=Z(J)+DZ
  100 CONTINUE
      RETURN
      END