MODULE WINDWAVE

USE GLOBAL
IMPLICIT NONE

INTEGER(4) ,PARAMETER::UFET=214       !FETCH.OUT
INTEGER(4) ,PARAMETER::UWIN=215       !LIJXY.OUT
INTEGER(4) ,PARAMETER::UTAU=216       !TAUW.OUT
REAL(RKD)  ,PARAMETER::WHMI=1D-3      !MINIMUM WAVE HEIGHT
REAL(RKD)            ::ROTAT          !ANTICLOCKWISE ROTATION OF DOMAIN [0,360]
INTEGER(4),PRIVATE,PARAMETER:: NZONE=8              !NUMBER OF ZONES
REAL(RKD) ,PRIVATE,PARAMETER:: FETANG(NZONE) = (/0,45,90,135,180,225,270,315/)

TYPE WAVE
  REAL(4)   , POINTER::TWX(:)           !BED SHEAR STRESS BY WAVE
  REAL(4)   , POINTER::TWY(:)
  REAL(RKD) , POINTER::UDEL(:)          !ORBITAL VELOCITY
  REAL(RKD) , POINTER::RLS(:)           !WAVE LENGTH
  REAL(RKD)          ::FW               !WAVE FRICTION COEFFICIENT
END TYPE WAVE

TYPE(WAVE)           ::WV
PRIVATE              ::FETZONE

CONTAINS

  SUBROUTINE WINDWAVETUR 
    USE GLOBAL  
    INTEGER(4)   ::L,IWVRDC  
    REAL(4)      ::RA,CDTMP,AEXTMP,UWORBIT,VISMUDD,REYWAVE
    REAL,EXTERNAL::CSEDVIS 

    CALL WINDWAVECAL         !OUTPUT: WVWHA,UDEL,WDIR,WVFRQL,RLS
    
    !**  GENERATE WAVE TABLE: IT IS NOT USED BY WIND WAVE  
    DO L=2,LA  
      WVKHP(L)=1.  
      IF(WVWHA(L).GE.WHMI) WVKHP(L)=2*PI/WV%RLS(L)*HP(L)
    ENDDO  

    ! *** ADJUST WAVE HEIGHTS DEPENDING ON VEGETATION 
    IF (ISVEG>0) THEN
      DO L=2,LA  
        IF(MVEGL(L).NE.MVEGOW)THEN
          ! *** FORCE WAVE HEIGHT TO BE ZERO FOR VEGETATIVE CELLS
          WVWHA(L)=0.  
        ELSE
          IWVRDC=0  
          ! *** IF SURROUNDING CELLS VEGETATIVE, REDUCE WAVE HIEGHT
          IF(MVEGL(L-1).NE.MVEGOW) IWVRDC=1  
          IF(MVEGL(L+1).NE.MVEGOW) IWVRDC=1  
          IF(MVEGL(LSC(L)).NE.MVEGOW) IWVRDC=1  
          IF(MVEGL(LNC(L)).NE.MVEGOW) IWVRDC=1  
          IF(IWVRDC.GT.0) WVWHA(L)=0.5*WVWHA(L) 
        ENDIF 
      ENDDO  
    ENDIF
    
    !**  INITIALIZE WAVE-CURRENT BOUNDARY LAYER MODEL CALCULATING  
    !**  THE WAVE TURBULENT INTENSITY, QQWV  
    !**  AND SQUARED HORIZONTAL WAVE OBRITAL VELOCITY MAGNITUDE  

    DO L=2,LA  
      IF(WVWHA(L).GE.WHMI)THEN
        UWORBIT=WV%UDEL(L)         
        AEXTMP=MAX(KSW,UWORBIT/WVFRQL(L))        !TO CONTROL FW
        UWVSQ(L)=UWORBIT*UWORBIT  
        IF (UWVSQ(L)<1.E-6) UWVSQ(L)=0.          ! PMC

        IF(ZBR(L).LE.0.)THEN  
          !** TURBULENT SMOOTH WAVE BOUNDARY LAYER  
          VISMUDD=1.E-6  
          IF(ISMUD.GE.1) VISMUDD=CSEDVIS(SED(L,1,1))  
          REYWAVE=UWORBIT*AEXTMP/VISMUDD  
          CDTMP=0.012/(REYWAVE**0.123)  
          QQWV1(L)=CDTMP*UWORBIT*UWORBIT   
        ELSE  
          !** TURBULENT ROUGH WAVE BOUNDARY LAYER  
          RA= AEXTMP/KSW
          WV%FW = EXP(5.5*RA**(-0.2)-6.3)  ! *** Nielsen (1992) for all RA's
          CDTMP=0.5*WV%FW
          QQWV1(L)=CDTMP*UWORBIT*UWORBIT  
        ENDIF  
      ELSE
        QQWV1(L)=0.  
        UWVSQ(L)=0.
        WVWHA(L)=0.
        WVFRQL(L)=0.
      ENDIF
      
      WV%TWX(L)=RHO*QQWV1(L)*WV%TWX(L)
      WV%TWY(L)=RHO*QQWV1(L)*WV%TWY(L)
    ENDDO  

    IF (TIMEDAY>=SNAPSHOTS(NSNAPSHOTS)) THEN
      WRITE(UTAU) (WV%TWX(L),WV%TWY(L),L=2,LA)
	  ![ykchoi 10.04.26
      !FLUSH(UTAU)
	  CALL FLUSH(UTAU)
	  !ykchoi]
    ENDIF  

  END SUBROUTINE

  SUBROUTINE WINDWAVECAL
    !CALCULATING WAVE PARAMETERS FOR EVERY CELL
    !BASED ON COMPUTED WIND PARAMETERS FROM WSER.INP AND SHELTERING
    !INPUT:
    !WNDVELE(L),WNDVELN(L),HP(L)
    !OUTPUT:
    !WVWHA(L),WVFRQL(L),WACCWE(L),WV%UDEL(L)
    !  WVWHA(L)  - WAVE HEIGHT (M)
    !  WACCWE(L) - WAVE ANGLE (RADIANS)
    !  WVFRQL(L) - WAVE FREQENCY (SEC)
    !  WV%TWX(L),WV%TWY(L)
    INTEGER(4) ::L,ZONE
    REAL(RKD)  ::FW,WB,TAUW,TP
    REAL(RKD)  ::AVEDEP,WVEL2,FC1,FC2,FC3
    REAL(RKD)  ::WDIR           ! WIND DIRECTION IN DEG [0,360]
    REAL(RKD)  ::WINX,WINY      !IN CURVI-LINEAR SYS
    REAL(RKD)  ::WVEL           !INTERPOLATED WIND VELOCITY

    !CALCULATING WAVE HEIGHT,PERIOD,ORBITAL VELOCITY AND LENGTH
    AVEDEP=SUM(HP(2:LA))/FLOAT(LA-1)
    DO L=2,LA
      !WINX  = WNDVELE(L)  !X IS TRUE EAST, Y IS TRUE NORTH
      !WINY  = WNDVELN(L)  
      WINX =  CVN(L)*WNDVELE(L) - CVE(L)*WNDVELN(L)  !CURVI-LINEAR SYS
      WINY = -CUN(L)*WNDVELE(L) + CUE(L)*WNDVELN(L)
      
      WVEL2 = WINX**2+WINY**2
      WVEL  = SQRT(WVEL2)

      IF (HP(L)>HDRY.AND.WVEL>1D-6) THEN
        WV%TWX(L)=WINX/WVEL
        WV%TWY(L)=WINY/WVEL
        !AVEDEP=HP(L)
        IF(WINX>=0) THEN
          WDIR  = ACOS(WV%TWY(L))*180./PI     !DEG. (NORTH,WIND TO)
        ELSE
          WDIR  = 360-ACOS(WV%TWY(L))*180./PI
        ENDIF
        ZONE = FETZONE(WDIR)
     
        ! *** WAVE HEIGHT
        FC3 =TANH(0.530*(9.81*AVEDEP/WVEL2)**0.75)
        FC1=WVEL2/9.81*0.283*FC3  
        FC2=TANH(0.0125*(9.81*FWDIR(L,ZONE)/WVEL2)**0.42/FC3)   
        WVWHA(L)=MIN(0.75*HP(L),FC1*FC2)           !INCLUDING BREAKING WAVE

        ! *** WAVE FREQUENCY
        FC3 = TANH(0.833*(9.81*AVEDEP/WVEL2)**0.375)
        FC1=(WVEL/9.81)*7.54*FC3
        FC2=TANH(0.077*(9.81*FWDIR(L,ZONE)/WVEL2)**0.25/FC3)   
        TP=MAX(1.0E-6,FC1*FC2)                      ! PERIOD
        WVFRQL(L)=2.0*PI/TP                       ! FREQUENCY OMEGA

        ! *** ORBITAL VELOCITY
        FC1=(2.0*PI/TP)**2*HP(L)/9.8
        FC2=FC1+1.0/(1.0+0.6522*(FC1)+0.4622*(FC1)**2+0.0864*(FC1)**4+0.0675*(FC1)**5)
        WV%RLS(L)=TP*SQRT(9.81*HP(L)/FC2)   ! *** WAVE LENGTH
        IF (HP(L)/WV%RLS(L)<100) THEN
          WV%UDEL(L)=MAX(1.0E-2,PI*WVWHA(L)/(TP*SINH(HP(L)*2.0*PI/WV%RLS(L))))
        ELSE
          WV%UDEL(L)= 1D-6
        ENDIF
        
        ! *** WAVE DIRECTION (RADIANS) ANTICLOCKWISE (CELL-EAST AXIS,WAVE)
        WACCWE(L)=(90-WDIR-ROTAT)*PI/180._8
        
      ELSE
        WVWHA(L)  = 0
        WVFRQL(L) = 0
        WACCWE(L) = 0
        WV%RLS(L) = 0
        WV%UDEL(L)= 0   
      ENDIF
    ENDDO

  END SUBROUTINE

  FUNCTION FETZONE(WDIR) RESULT(ZONE)
    !DETERMINING FETCH ZONE AND FETCH MAIN ANGLE
    !BASED ON THE GIVEN WIND DIRECTION WDIR
    !WDIR     : INTERPOLATED WIND DIRECTION FROM WSER.INP
    !UNIT     : [0,360]
    !FORMATION: ANGLE BY (NORTH,WIND TO)IN CLOCKWISE DIRECTION
    !ZONE 1: NORTH       >337.5 OR <=22.5
    !ZONE 2: NORTH-EAST  >22.5  OR <=67.5
    !ZONE 3: EAST        >
    !ZONE 4: SOUTH-EAST
    !ZONE 5: SOUTH
    !ZONE 6: SOUTH-WEST
    !ZONE 7: WEST
    !ZONE 8: NORTH-WEST
    REAL(RKD) ,INTENT(IN )::WDIR   ![0,360]
    INTEGER(4)::ZONE

    IF     (WDIR>337.5 .OR. WDIR <= 22.5) THEN
     ZONE = 1
    ELSEIF (WDIR>22.5 .AND. WDIR <= 67.5) THEN
     ZONE = 2
    ELSEIF (WDIR>67.5 .AND. WDIR <= 112.5) THEN
     ZONE = 3
    ELSEIF (WDIR>112.5 .AND. WDIR <= 157.5) THEN
     ZONE = 4
    ELSEIF (WDIR>157.5 .AND. WDIR <= 202.5) THEN
     ZONE = 5
    ELSEIF (WDIR>202.5 .AND. WDIR <= 247.5) THEN
     ZONE = 6
    ELSEIF (WDIR>247.5 .AND. WDIR <= 292.5) THEN
     ZONE = 7
    ELSEIF (WDIR>292.5 .AND. WDIR <= 337.5) THEN
     ZONE = 8
    ENDIF
  END FUNCTION
 
  SUBROUTINE FETCH
    !DETERMINING THE FETCHES OF ALL CELLS:
    !OUTPUT: FWDIR(2:LA,1:NZONE) IN M
    USE DRIFTER,ONLY:INSIDECELL
    REAL(RKD)::AL(NZONE),RL,XM,YM,RL0
    INTEGER(4)::I,J,L,NZ,IM,JM,LM,STATUS,MUL

    OPEN(UFET,FILE='FETCH.OUT')
    FWDIR = 0
    AL = (180+90-FETANG-ROTAT)*PI/180._8   !ANTICLOCKWISE (X',WIND FR)
    RL0 = 0.25*MIN(MINVAL(DXP(2:LA)),MINVAL(DYP(2:LA)))
    DO L=2,LA
     DO NZ=1,NZONE
       RL=0
       IM=IL(L)
       JM=JL(L)
       LOOP1:DO WHILE(.TRUE.)
         STATUS=0
         RL=RL+RL0                        !SEARCH FETCH FOR EVERY 10M
         XM = XCOR(L,5)+RL*COS(AL(NZ))    !UPWIND DISTANCE=FETCH
         YM = YCOR(L,5)+RL*SIN(AL(NZ))
         LOOP2:DO J=JM-1,JM+1
           DO I=IM-1,IM+1
             LM=LIJ(I,J)
             IF (LM<2) CYCLE
             IF (INSIDECELL(LM,XM,YM)) THEN
               STATUS=1
               EXIT LOOP2
             ENDIF
           ENDDO
         ENDDO LOOP2
         IF (STATUS==0) EXIT LOOP1
         IM=IL(LM)
         JM=JL(LM)
       ENDDO LOOP1
       FWDIR(L,NZ)=RL
     ENDDO
     WRITE(UFET,'(2I6,8F15.4)') IL(L),JL(L),(FWDIR(L,NZ),NZ=1,NZONE)
    ENDDO
    CLOSE(UFET)
  END SUBROUTINE 
 
  SUBROUTINE WINDWAVEINIT   
    ! *** INITIALIZES WAVE VARIABLES AND GENERATES FETCH.OUT
    USE GLOBAL  
    INTEGER(4)   ::L,K

    ALLOCATE(WV%TWX(LA),WV%TWY(LA))
    ALLOCATE(WV%UDEL(LA),WV%RLS(LA))
    WV%TWX  = 0.0
    WV%TWY  = 0.0
    WV%UDEL = 0.0
    WV%RLS  = 0.0
    ROTAT   = 0.0

    KSW = MAX(1.0E-6,KSW)
    PRINT *,'COMPUTING FETCH'
    CALL FETCH
    OPEN(UWIN,FILE='LIJXY.OUT',ACTION='WRITE')
    DO L=2,LA
      WRITE(UWIN,'(3I10,2F15.5)') L,IL(L),JL(L),DLON(L),DLAT(L)
    ENDDO
    CLOSE(UWIN)
    OPEN(UTAU,FILE='TAUW.OUT',FORM='UNFORMATTED')

    JSWRPH=1  
    DO L=1,LC  
      HMPW(L)=0.  
      HMCW(L)=0.  
      HMUW(L)=0.  
      HMVW(L)=0.  
      WVWHA(L)=0.  
      WVKHP(L)=0.  
      WVKHC(L)=0.  
      WVKHU(L)=0.  
      WVKHV(L)=0.  
      WVTMP1(L)=0.  
      WVTMP2(L)=0.  
      WVTMP3(L)=0.  
      WVTMP4(L)=0.  
      UWVMAG(L)=0.  
      VWVMAG(L)=0.  
      WVENEP(L)=0.  
      UWVSQ(L)=0.  
      QQWC(L)=1.E-12  
      QQWCR(L)=1.E-12  
      QQWV1(L)=1.E-12  
      QQWV2(L)=0.
      QQWV3(L)=1.E-12  
      WACCWE(L)=0.  
    ENDDO  
    DO K=1,KC  
      DO L=1,LC  
       WVHUU(L,K)=0.  
       WVHVV(L,K)=0.  
       WVHUV(L,K)=0.  
       WVPP(L,K)=0.  
       WVPU(L,K)=0.  
       WVPV(L,K)=0.  
       WVDISP(L,K)=0.  
       FXWAVE(L,K)=0.  
       FYWAVE(L,K)=0.  
      ENDDO  
    ENDDO  
   
    ITWCBL1=1  
    ITWCBL2=0  
    ITWCBL3=0  
    ISWCBL =0  ! PMC ONLY FOR ISWAVE = 1 OR 3

  END SUBROUTINE
 
END MODULE