SUBROUTINE SBOT (ABRBOT ,DEP2 ,ECOS ,ESIN ,KWAVE , & SPCSIG ,UBOT ,UX2 ,UY2 ,IDCMIN , & IDCMAX ,ISSTOP ,VARFR ,FRCOEF ,IG , & IMATRA ) ! (This subroutine has been tested only for case IBOT = 1) ! !**************************************************************** ! ! Computation of the source terms due to bottom friction ! ! Method ! ! In SWAN several bottom dissipation models are computed, i.e.: ! ! IBOT = 1 Jonswap bottom friction model ! IBOT = 2 Collins bottom friction model ! IBOT = 3 Madsen bottom friction model (see Tolman) ! ! Both methods are implemented in SWAN and the user has to make ! a choice in the input file. ! ! 1. Jonswap model: ! ----------------- ! ! The bottom interaction term SEbf(s,d) is supposed to take the ! Jonswap form (Hasselman et al. 1973): ! 2 ! sigma E(s,d) ! SEbf = -GAMMA ---------------- ! 2 2 ! g sinh (kD) ! 2 -3 ! where GAMMA is the decay parameter ,(default GAMMA = 0.067 m s ). ! In the Jonswap form the current velocities are not taken into ! account. ! ! 2. COLLINS model: ! ----------------- ! ! The energy dissipation due to bottom friction is modelled ! according the quadratic friction law: ! 2 ! SE = Tau * |U| ! ! which for a spectrum can be written as: ! 2 ! sigma ! SE(s,d)= - ---------------- * (Cfw.Ub + Cfc.Uc) * E(s,d) ! 2 ! g sinh (K(s) * D) ! ! Ub is the velocity due to the wave at the bottom ! ! The current velocity is Uc ! ! 2. MADSEN formulation: ! ---------------------- ! ! The bottom dissipation applying Madsen formulation is as ! follows: ! ! fw [n - 1/2] UBR E(s,d) ! [1] Sdsb(s,d) = - ------------------------ ! D ! ! in which : ! 2 ! s * D ! [1a] (n - 1/2) = ------------- ! 2 ! 2 g sinh (kD) ! ! UBOT(IX,IY) is computed in the subroutine SINTGRL. The friction ! factor fw is estimated using the formulation of Jonsson (1963, ! 1966a): ! ! 1 1 Ab,r ! [2] -------- + log { ---------- } = mf + log { ----- } ! 4 sqrt(fw) 10 4 sqrt(fw) 10 Kn ! ! with: ! ! 2 // 1 ! [3] Ab,r = 2 * // -------------- E(s,d) ds dd ! // 2 ! sinh (kD) ! ! with: Ab,r is the representative near bottom excursion ! amplitude ! Kn equivalent roughness ! mf constant ( mf = -0.08) (determined by Jonsson ! and Carlssen 1976 ) ! ! [2] is only valid for Ab,r/Kn larger than approximately 1. ! For smaller values a constant value of fw is used (fw = 0.3 ! for Ab,r/Kn < 1.57 ) ! USE SWCOMM3 USE SWCOMM4 USE OCPCOMM4 ! USE ALL_VARS, ONLY : MT,AC2 USE VARS_WAVE, ONLY : MT,AC2 IMPLICIT NONE REAL SPCSIG(MSC) INTEGER ID,IS,ISSTOP,J,IDDUM,IG REAL :: XDUM,KD,SBOTEO,CFBOT,FACB,CFW,FW,CURR,UC,ABRBOT,ADUM,CDUM,DDUM LOGICAL VARFR REAL :: DEP2(MT),ECOS(MDC),ESIN(MDC),IMATDA(MDC,MSC),IMATRA(MDC,MSC), & KWAVE(MSC,ICMAX),PLBTFR(MDC,MSC,NPTST),UBOT(MT), & UX2(MT),UY2(MT),DISSC1(MDC,MSC),FRCOEF(MT) INTEGER :: IDCMIN(MSC),IDCMAX(MSC) REAL :: AKN IF(IBOT >= 1 .AND. DEP2(IG) > 0.)THEN IF(IBOT == 1)THEN ! ! *** Jonswap model *** ! ! PBOT(3) = GAMMA (a) in the Jonswap equation, ! CFBOT = PBOT(3) / GRAV_W**2 ! CFBOT = CFBOT*10.0 !JQI added ELSE IF (IBOT == 2)THEN ! ! *** Collins model *** ! ! PBOT(2) = [cfw] ! IF(VARFR)THEN CFW = FRCOEF(IG) ELSE CFW = PBOT(2) ENDIF CFBOT = CFW * UBOT(IG) / GRAV_W ELSE IF(IBOT == 3)THEN ! ! *** Madsen model *** ! IF(VARFR)THEN AKN = FRCOEF(IG) ELSE AKN = PBOT(5) ENDIF ! ! *** PBOT(4) = Mf *** ! *** AKN = PBOT(5) = [kn] (roughness) *** ! IF((ABRBOT/AKN) > 1.57)THEN XDUM = PBOT(4) + LOG10 ( ABRBOT / AKN ) ! ! *** solving the implicit equation using a Newton *** ! *** Rapshon iteration proces : a + log a = b *** ! *** the start value for ADUM = 0.3 because 0.3626 *** ! *** is the minimum value of ADUM with b=-0.08. *** ! ADUM = 0.3 DO J = 1, 50 CDUM = ADUM DDUM = ( ADUM + LOG10(ADUM) - XDUM ) / ( 1.+ ( 1. / ADUM) ) ADUM = ADUM - DDUM IF(ABS(CDUM - ADUM) < 1.E-4) EXIT END DO IF(J .eq. 51 ) WRITE(*,*) ' error in iteration fw: Madsen formulation' ! 1 1 ! *** computation of FW --> A = ----- --> FW = ----- ! 4 uFW 16 A**2 FW = 1. / (16. * ADUM**2) ELSE FW = 0.3 ENDIF CFBOT = UBOT(IG) * FW / (SQRT(2.) * GRAV_W) ENDIF DO IS = 1, ISSTOP KD = KWAVE(IS,1) * DEP2(IG) IF(KD < 10.)THEN FACB = CFBOT * (SPCSIG(IS) / SINH(KD)) **2 ! DO IDDUM = IDCMIN(IS) , IDCMAX(IS) ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1 ! SBOTEO = FACB IF(IBOT == 2 .AND. ICUR == 1 .AND. PBOT(1) > 0.)THEN ! additional dissipation due to current, seldom used CURR = UX2(IG)*ECOS(ID) + UY2(IG)*ESIN(ID) UC = ABS(CURR) ! PBOT(1) = [cfc] SBOTEO = FACB + PBOT(1) * UC * (SPCSIG(IS) / SINH(KD)) **2 END IF ! ! *** store the results in the array IMATDA *** ! IMATRA(ID,IS) = IMATRA(ID,IS) - SBOTEO*AC2(ID,IS,IG) IMATDA(ID,IS) = IMATDA(ID,IS) + SBOTEO ! IF (TESTFL) PLBTFR(ID,IS,IPTST) = -1.* SBOTEO ! DISSC1(ID,IS) = DISSC1(ID,IS) + SBOTEO END DO !IDDUM END IF END DO !IS ! ENDIF ! RETURN END SUBROUTINE SBOT ! !**************************************************************** ! SUBROUTINE FRABRE ( HM, ETOT, QBLOC ) ! !**************************************************************** ! ! to compute the fraction of breaking waves in point ix,iy ! of the computational grid ! ! Method (updated...) ! ! The fraction of breaking waves in a point ix,iy is given by ! the implicit relation: ! ! 1 - Qb ETOT ! ------ = -8 * ----- ! ln Qb HM**2 ! ! from which Qb can be found by solving the equation: ! ! ETOT ! F = 1 - Qb + 8 * ---- * ln(Qb) = 0. ! 2 ! HM ! ! The following appproximation is applied: ! ! 2 ! (1)| B = sqrt( 8 ETOT/HM ), i.e. B = Hrms/HM ! ! ! | Qo = 0. B <= 0.5 ! (2)| 2 ! | Qo = ( 2B -1 ) 0.5 < B <= 1 ! ! ! applying the Newton-Raphson procedure (for 0.2= 1.0 ! | ! USE SWCOMM4 USE OCPCOMM4 ! IMPLICIT NONE REAL :: ETOT, HM, QBLOC INTEGER :: IENT REAL :: B, B2, QO, Z IF((HM > 0.) .AND. (ETOT >= 0.))THEN B = SQRT(8. * ETOT / (HM*HM) ) ELSE B = 0.0 END IF ! IF(B <= 0.5)THEN QO = 0. ELSE IF(B <= 1.0)THEN QO = (2.*B - 1.)**2 END IF ! IF(B <= 0.2)THEN QBLOC = 0.0 ELSE IF(B < 1.0)THEN ! ! *** second iteration to find Qb *** ! B2 = B*B Z = EXP((QO-1.)/B2) QBLOC = QO - B2 * (QO-Z)/(B2-Z) ELSE QBLOC = 1.0 END IF ! RETURN END SUBROUTINE FRABRE ! !**************************************************************** ! SUBROUTINE SSURF (ETOT ,HM ,QB ,SMEBRK , & IMATRA ,IMATDA ,IDCMIN ,IDCMAX , & PLWBRK ,ISSTOP ,IG ) ! !**************************************************************** ! Computation of the source term due to wave breaking. ! White capping is not taken into account ! ! Method ! ! The source term for surf breaking is implemented following ! the approach of Battjes/Janssen (1978) for the energy dissipation: ! ! Alpha - 2 - SMEBRK ! Dtot = ---- Qb * f * Hm with f = ------ ! 4 2 * Pi ! ! Now the source term is: ! ! SIGMA * AC2(ID,IS,IX,IY) ! Sbr = - Dtot * ------------------------ = ! Etot ! ! ! Alpha * SMEBRK * Qb * Hm * Hm SIGMA * AC2(ID,IS,IX,IY) ! = - ------------------------------ * ------------------------- ! 8 * Pi Etot ! ! ! = WS * SIGMA * AC2(ID,IS,IX,IY) = WS * E ! ! ! with ! ! Alpha = PSURF(1) ; ! ! SMEBRK Qb ! WS = Alpha ------ -- ; ! Pi BB ! 2 ! BB = 8 Etot / Hm = - (1 - Qb) / ln (Qb) ; ! ! ! The local maximum wave height Hm and mean frequency SMEBRK are computed ! in subroutine SINTGRL. ! The fraction of breaking waves Qb is calculated in the subroutine FRABRE ! ! The new value for the dissipation is computed implicitly using ! the last computed value for the action density Nold (at the spatial ! gridpoint under consideration). ! ! Sbr = WS * N ! ! = Sbr_new + (d Sbr/d N) (Nnew - Nold) ! ! = WS * Nnew + SbrD * (Nnew - Nold) ! ! = (WS + SbrD)* Nnew - SbrD * Nold ! ! = SURFA1 * Nnew - SURFA0 * Nold ! ! In order to do this we need the derivative ! of the source term Sbr to the action density N ! ! d Sbr d WS ! SbrD = ----- = ---- * N + WS ! d N d N ! ! Since BB and SMEBRK * N are proportional, we have ! ! d Sbr d WS SMEBRK (d Qb/ d BB) *BB - Qb ! ----- = ---- * BB + WS = Alpha ------ --------------------- * BB + WS = ! d N d BB Pi sqr(BB) ! ! ! SMEBRK d Qb ! Alpha ------ ---- ! Pi d BB ! ! With: ! ! d Qb 1 ! ---- = ------------- ; ! d BB (d BB / d Qb) ! ! 2 ! d Qb ln (Qb) ! ---- = --------------------------- ! d BB ln (Qb) + (1 - Qb) (1 / Qb) ! ! Qb (1 - Qb) ! = ------------ ; ! BB (BB - Qb) ! ! ------------------------------------------------------------ ! Get HM, QB and ETOT from the subroutine SINTGRL ! For spectral direction IS and ID do, ! get the mean energy frequency average over the full spectrum ! If ETOT > 0 then ! compute source term for energy dissipation SURFA0 and SURFA1 ! Else ! source term for wave breaking is 0. ! End if ! ---------------------------------------------------------- ! Compute source terms for energy averaged frequency ! Store results in the arrays IMATDA and IMATRA ! ------------------------------------------------------------ USE SWCOMM3 USE SWCOMM4 USE OCPCOMM4 ! USE ALL_VARS, ONLY : MT,AC2 USE VARS_WAVE, ONLY : MT,AC2 ! IMPLICIT NONE INTEGER :: ISSTOP,IDCMIN(MSC),IDCMAX(MSC) REAL :: DISSC0(MDC,MSC),DISSC1(MDC,MSC), & IMATDA(MDC,MSC),IMATRA(MDC,MSC),PLWBRK(MDC,MSC,NPTST) REAL :: ETOT,HM,QB,SMEBRK INTEGER :: ID,IDDUM,IENT,IS,IG DOUBLE PRECISION BB,DIS0,SbrD,SURFA0,SURFA1,WS ! ! ALFA = PSURF(1) ! BB = 8. * DBLE(ETOT) / ( DBLE(HM)**2 ) ! SURFA0 = 0. ! SURFA1 = 0. ! IF(REAL(BB) > 0. .AND. REAL(ABS(BB - DBLE(QB))) > 0.)THEN ! IF(BB < 1.)THEN ! WS = ( DBLE(PSURF(1)) / DBLE(PI_W)) * DBLE(QB) * DBLE(SMEBRK) / BB ! SbrD = WS * (1. - DBLE(QB)) / (BB - DBLE(QB)) ! ELSE ! WS = ( DBLE(PSURF(1)) / DBLE(PI_W)) * DBLE(SMEBRK) ! SbrD = 0. ! END IF ! SURFA0 = SbrD ! SURFA1 = WS + SbrD ! ELSE ! SURFA0 = 0. ! SURFA1 = 0. ! END IF IF(REAL(BB) > 0. .AND. REAL(ABS(BB - DBLE(QB))) > 0.)THEN WS = ( DBLE(PSURF(1)) / DBLE(PI_W)) * DBLE(QB) * DBLE(SMEBRK) / BB ELSE WS = 0. ENDIF ! ! *** store the results for surf wave breaking *** ! *** in the matrices IMATDA and IMATRA *** ! DO IS = 1, ISSTOP DO IDDUM = IDCMIN(IS), IDCMAX(IS) ID = MOD ( IDDUM - 1 + MDC , MDC ) + 1 IMATDA(ID,IS) = IMATDA(ID,IS) + REAL(WS) DIS0 = WS * DBLE(AC2(ID,IS,IG)) IMATRA(ID,IS) = IMATRA(ID,IS) - REAL(DIS0) END DO END DO RETURN END SUBROUTINE SSURF ! !**************************************************************** ! SUBROUTINE SWCAP (SPCDIR ,SPCSIG ,KWAVE ,IDCMIN , & IDCMAX ,ISSTOP ,ETOT ,IMATDA , & IMATRA ,PLWCAP ,CGO ,UFRIC , & DEP2 ,DISSIP ,DISIMP ,IG ) ! (This subroutine has not been fully tested. Only tested for ! case IWCAP = 2) ! !**************************************************************** ! ! Calculates the dissipation due to whitecapping ! ! Method ! ! Whitecapping dissipation is formulated as follows: ! ! S_wc(sig,th) = - C_wc E(sig,th) ! ! where the coefficient C_wc has four basic forms: ! ! C_wc1 = C_K sig~ (k/k~): According to Komen (generalised) ! C_wc2 = C_BJ sig~ (k/k~): According to Battjes-Janssen (modified) ! C_wc3 = C_LH : According to Longuet Higgins (1969) ! C_wc4 = C_CSM : According to cumulative steepness ! method (Ris et al. (1999) and Donelan (1999)) ! ! In these formulations C_K is defined as (Komen; 1994 p. 145): ! ! n1 n2 ! C_K = C1 [(1-delta) + delta (k/k~) ] (S~/S_PM) ! ! where C1, delta, n1 and n2 can be varied ! ! ! C_BJ is defined as: ! ! 2 ! alpha Hm Qb ! C_BJ = ----------- ! 8 Pi m0 ! ! where alpha can be varied. ! ! for Hrms > Hm the formulation changes in a limit to (Hrms->Hm; Qb->1): ! ! alpha ! C_BJ = ----- ! Pi ! ! and C_LH is defined as (Komen; 1994 p. 151-152): ! ! 4 2 2 ! C_LH = C3 Sqrt[(m0 sig0 )/g ] exp(A) sig0 (sig/sig0) ! ! where ! 2 2 4 2 2 ! A = -1/8 [1-eps ] [g /(m0 sig0 )] with [1-eps ] = [m2 ] / [m0 m4] ! ! and C3 can be varied ! ! ! Cumulative steepness method. Basic idea is that dissipation depends ! on the steepness of the wave spectrum at and below a particular ! frequency. Details can be found in "SWAN fysica plus" (2002), ! report H3937/A832, implemented by Delf Hydraulics and Alkyon by ! order of RIKZ/RWS as a part of the project HR-Ontwikkeling. ! 40.41 ! Recent modifications are described in: 40.41 ! Hurdle, D.P., and G.Ph. van Vledder, 2004: Improved spectral wave 40.41 ! modelling of white-capping dissipation in swell sea systems. Proc. 40.41 ! Offshore Mechanics, Arctic Engineering Conference, June 2004, 40.41 ! Vancouver, Canada, paper OMAE2004-51562. 40.41 ! 40.41 ! Default values: C_st = 4 40.41 ! m = 2 40.41 ! 40.41 ! C_CSM = A*C_st S(sig,th) 40.41 ! ! where ! ! C_st is a tuneable coefficient and ! ! S(sig,th) = integrals from 0 to sig and from 0 to 2pi of ! ! k^2 |cos(th-th')|^m E(sig,th) dsig dth' ! ! where coefficient m controls the directional dependence in case ! of straining mechanism (if dominant then m = 2 else m > 10) and ! th is the actual direction where straining mechanism takes ! place. ! 40.41 ! A is a normalisation coefficient computed 40.41 ! 40.41 ! 1 GAMMA(0.5*m+1.) 40.41 ! A = -------- ---------------- 40.41 ! SQRT(PI) GAMMA(0.5*m+0.5) 40.41 ! ! In these equations the variables have the following meaning: ! ! Hm : Maximum wave height ! Hrms : Root mean square of the wave heights ! eps^2: Measure for the spectral bandwidth ! m0 : Total wave energy density (=ETOT) ! m2 : Second moment of the variance spectrum (=ETOT2) ! m4 : Fourth moment of the variance spectrum (=ETOT4) ! k : Wave number (=KWAVE(IS,1)) ! k~ : Mean wave number ! Qb : Fraction of breaking waves ! sig : Frequency (=SPCSIG(IS)) ! sig0 : Average zero crossing frequency ! sig~ : Mean frequency ! S~ : Overall steepness (STP_OV) ! S_PM : Overall steepness for a Pierson-Moskowitz spectrum ! th : direction theta (=SPCDIR(ID)) ! ! ------------------------------------------------------------ ! Calculate needed parameters ! ------------------------------------------------------------ ! If IWCAP = 1, 2, or 5; Calculate C_K ! If IWCAP = 4 or 5; Calcualte C_BJ ! If IWCAP = 3; Calculate C_LH ! If IWCAP = 6; Calculate cumulative steepness spectrum ! IF IWCAP = 7; Calculate with Alves and Banner (2003) ! ------------------------------------------------------------ ! For frequency dependent part of the spectrum ! Calculate dissipation term due to whitecapping ! ------------------------------------------------------------ ! For the whole frequency domain ! Fill the matrices and PLWCAP-array ! ------------------------------------------------------------ ! End of SWCAP ! ------------------------------------------------------------ USE SWCOMM3 USE SWCOMM4 USE OCPCOMM4 USE M_WCAP ! USE ALL_VARS, ONLY : MT,AC2 USE VARS_WAVE, ONLY : MT,AC2 IMPLICIT NONE INTEGER, INTENT(IN) :: ISSTOP, IDCMIN(MSC), IDCMAX(MSC) REAL, INTENT(IN) :: DEP2(MT) REAL, INTENT(IN) :: ETOT REAL, INTENT(IN) :: KWAVE(MSC,MICMAX) REAL, INTENT(IN) :: SPCDIR(MDC,6), SPCSIG(MSC) REAL, INTENT(OUT) :: PLWCAP(MDC,MSC,NPTST) REAL, INTENT(IN OUT) :: IMATDA(MDC,MSC), IMATRA(MDC,MSC) REAL, INTENT(IN OUT) :: DISSIP(MDC,MSC), DISIMP(MDC,MSC) REAL, INTENT(IN) :: UFRIC REAL, INTENT(IN) :: CGO(MSC,MICMAX) INTEGER, SAVE :: IENT = 0 INTEGER :: ID, IDDUM, IS, ID1, ID2, IF, IL, MXWCP,IG REAL :: A, C_BJ, HM, HRMS, N1, N2 REAL :: QB_WC, SIG0, STP_OV, STP_PM REAL :: DDIF, DELTA, DSTEEP, EBIN, XFAC REAL :: CPOW, CTOT, GAMMA REAL :: BINSIZE REAL, ALLOCATABLE :: C_K(:), C_LH(:), WCAP(:), WCIMPL(:) REAL, ALLOCATABLE :: CUMSTP(:,:), DSIGMA(:) REAL, ALLOCATABLE :: FCOS(:) REAL :: B, P REAL :: EF(MSC) CHARACTER*20 NUMSTR, CHARS CHARACTER*80 MSGSTR MXWCP = 7 IF(IWCAP > MXWCP)THEN ! ! Error message ! CHARS = NUMSTR(MXWCP+1,RNAN,'(I1)') CALL TXPBLA(CHARS,IF,IL) MSGSTR = 'Value for IWCAP should be less than '// CHARS(IF:IL) CALL MSGERR ( 4, MSGSTR ) RETURN END IF ! ! Initialisation ! IF(ETOT <= 0.) RETURN IF(ETOT2 <= 0.) RETURN IF(ETOT4 <= 0.) RETURN IF(ACTOT <= 0.) RETURN IF(EDRKTOT <= 0.) RETURN ! ALLOCATE (C_K(MSC), C_LH(MSC), WCAP(MSC), WCIMPL(MSC)) ALLOCATE (CUMSTP(0:MSC,MDC), DSIGMA(1:MSC)) ALLOCATE (FCOS(1:MDC)) WCIMPL(1:MSC) = 0. ! ! Calculate coefficients ! IF((IWCAP == 1) .OR. (IWCAP == 2) .OR. (IWCAP == 5))THEN ! ! Calculate C_K ! STP_OV = KM_WAM * SQRT(ETOT) STP_PM = SQRT(PWCAP(2)) N1 = PWCAP(11) N2 = 2. * PWCAP(9) C_K(:) = PWCAP(1) * (1. - PWCAP(10) + & PWCAP(10) * (KWAVE(:,1) / KM_WAM)**N1) * & (STP_OV / STP_PM)**N2 ! ENDIF ! IF((IWCAP == 4) .OR. (IWCAP == 5))THEN ! ! Calculate values for Hm and Qb ! HRMS = SQRT(8. * ETOT) IF (IWCAP.EQ.4) HM = PWCAP(6) / KM01 IF (IWCAP.EQ.5) HM = PWCAP(6) / (PWCAP(8) * KM_WAM) CALL FRABRE(HM, ETOT, QB_WC) ! ! Calculate C_BJ ! IF(HRMS >= HM)THEN C_BJ = PWCAP(7) / PI_W ELSE IF(HRMS > 0.)THEN C_BJ = (PWCAP(7) * HM**2 * QB_WC) / (PI_W * HRMS**2) ELSE C_BJ = 0. END IF ENDIF ! IF(IWCAP == 3)THEN ! ! Calculate C_LH ! SIG0 = SQRT(ETOT2 / ETOT) ! ! A = -(1./8.)*(ETOT2**2/(ETOT*ETOT4))*(GRAV_W**2/(ETOT*SIG0**4)) ! rewrite to prevent underflow ! A = -(1./8.) * GRAV_W**2 / ETOT4 DO IS=1, ISSTOP ! C_LH(IS) = PWCAP(5) * SQRT((ETOT * SIG0**4) / GRAV_W**2) * ! & EXP(A) * SIG0 * (SPCSIG(IS) / SIG0)**2 ! rewrite to prevent underflow: ! C_LH(IS) = PWCAP(5) * EXP(A) * SQRT(ETOT2) * SPCSIG(IS)**2 / GRAV_W END DO END IF ! ! In case of cumulative steepness method, calculate its spectrum ! IF(IWCAP == 6)THEN XFAC = (SPCSIG(3)-SPCSIG(1))/(2.*SPCSIG(2)) DSIGMA(:) = XFAC*SPCSIG(:) DELTA = SPCDIR(2,1)-SPCDIR(1,1) ! ! --- precompute cos dependence ! DO ID1 = 1, MDC DDIF = REAL(ID1-1)*DELTA FCOS(ID1) = (ABS(COS(DDIF)))**PWCAP(13) END DO ! ! --- compute normalisation coefficient ! CPOW = PWCAP(13) IF(CPOW > 10)THEN CTOT = SQRT(CPOW/(2.*PI_W))*(1.+0.25/CPOW) ELSE CTOT = 1./SQRT(PI_W)*GAMMA(0.5*CPOW+1.)/GAMMA(0.5*CPOW+0.5) END IF ! CUMSTP = 0. DO IS = 1,ISSTOP BINSIZE = SPCSIG(IS)*DELTA*DSIGMA(IS) DO ID1 = 1,MDC CUMSTP(IS,ID1) = CUMSTP(IS-1,ID1) DO ID2 = 1,MDC EBIN= AC2(ID2,IS,IG)*BINSIZE ! EBIN= AC2(ID2,IS,kcgrd(1))*BINSIZE DSTEEP = KWAVE(IS,1)**2*EBIN*FCOS(ABS(ID1-ID2)+1) CUMSTP(IS,ID1) = CUMSTP(IS,ID1) + DSTEEP END DO END DO END DO CUMSTP = PWCAP(12)*CUMSTP CUMSTP = CUMSTP*CTOT END IF ! ! Calculate dissipation according to Alves & Banner (2003) ! IF(IWCAP == 7)THEN ! ! Loop to calculate B(k) ! DO IS = 1, ISSTOP ! ! Calculate E(f) ! EF(IS) = 0. DO ID = 1,MDC EF(IS) = EF(IS) + AC2(ID,IS,IG)*SPCSIG(IS)*PI2_W*DDIR END DO ! ! Calculate saturation spectrum B(k) from E(f) ! B = (1./PI2_W) * CGO(IS,1) * KWAVE(IS,1)**3 * EF(IS) ! ! Calculate exponent P of the relative saturation (B/Br) ! PWCAP(10)= 3. + TANH(25.76*(UFRIC*KWAVE(IS,1)/SPCSIG(IS)-0.1)) P = 0.5*PWCAP(10)*(1. + TANH( 10.*( (B/PWCAP(12))**0.5 - 1.))) ! ! Calculate WCAP(IS) from B(k) and P ! STP_OV = KM_WAM * SQRT(ETOT) WCAP(IS) = PWCAP(1)*(B/PWCAP(12))**(P/2.) * & STP_OV**PWCAP(9) * (KWAVE(IS,1)/KM_WAM)**PWCAP(11) * & (GRAV_W**(0.5)*KWAVE(IS,1)**(0.5)/SPCSIG(IS))**(PWCAP(10)/2-1) * & GRAV_W**(0.5)*KWAVE(IS,1)**(0.5) END DO END IF ! IF(IWCAP < 6)THEN ! ! Calculate the whitecapping source term WCAP(IS) ! DO IS=1, ISSTOP IF((IWCAP == 1) .OR. (IWCAP == 2) .OR. & ((IWCAP == 5) .AND. (C_BJ <= C_K(IS))))THEN WCAP(IS) = C_K(IS) * SIGM_10 * (KWAVE(IS,1) / KM_WAM) ELSE IF(IWCAP == 3)THEN WCAP(IS) = C_LH(IS) ELSE IF((IWCAP == 4) .OR. ((IWCAP == 5) .AND. (C_BJ >= C_K(IS))))THEN IF(IWCAP == 4) WCAP(IS) = C_BJ*SIGM01 *(KWAVE(IS,1)/KM01 ) IF(IWCAP == 5) WCAP(IS) = C_BJ*SIGM_10*(KWAVE(IS,1)/KM_WAM) ! ! Calculate a term that is added to both sides of the equation to compensate ! for the strong non-linearity in the fraction of breaking waves Qb ! IF(HRMS < HM)THEN WCIMPL(IS)=WCAP(IS) * ((1.-QB_WC)/((HRMS**2/HM**2)-QB_WC)) WCAP(IS) =WCAP(IS) + WCIMPL(IS) END IF ELSE CALL MSGERR(2,'Whitecapping is inactive') WRITE (PRINTF,*) 'Occurs in gridpoint: ', IG END IF END DO END IF ! ! Fill the diagonal of the matrix and the PLWCAP-array ! IF(IWCAP /= 6)THEN DO IS=1, ISSTOP ! ! Only fill the values for the current sweep ! DO IDDUM = IDCMIN(IS), IDCMAX(IS) ID = MOD(IDDUM - 1 + MDC, MDC) + 1 IMATDA(ID,IS) = IMATDA(ID,IS) + WCAP(IS) IMATRA(ID,IS) = IMATRA(ID,IS) - WCAP(IS) * AC2(ID,IS,IG) DISSIP(ID,IS) = DISSIP(ID,IS) + WCAP(IS) IF (TESTFL) PLWCAP(ID,IS,IPTST) = -1.*(WCAP(IS)-WCIMPL(IS)) END DO END DO ELSE DO IS=1, ISSTOP ! ! Only fill the values for the current sweep ! DO IDDUM = IDCMIN(IS), IDCMAX(IS) ID = MOD(IDDUM - 1 + MDC, MDC) + 1 IMATDA(ID,IS) = IMATDA(ID,IS) + CUMSTP(IS,ID) DISSIP(ID,IS) = DISSIP(ID,IS) + CUMSTP(IS,ID) IF (TESTFL) PLWCAP(ID,IS,IPTST) = -1.*CUMSTP(IS,ID) END DO END DO END IF ! ! Add the implicit part to the right-hand side ! IF((IWCAP == 4) .OR. (IWCAP == 5))THEN DO IS=1, ISSTOP ! ! Only fill the values for the current sweep ! DO IDDUM = IDCMIN(IS), IDCMAX(IS) ID = MOD(IDDUM - 1 + MDC, MDC) + 1 IMATRA(ID,IS) = IMATRA(ID,IS) + WCIMPL(IS) * AC2(ID,IS,IG) DISIMP(ID,IS) = DISIMP(ID,IS) + WCIMPL(IS) * AC2(ID,IS,IG) END DO END DO END IF ! DEALLOCATE (C_K, C_LH, WCAP, WCIMPL) DEALLOCATE (CUMSTP, DSIGMA, FCOS) ! RETURN END SUBROUTINE SWCAP !**************************************************************** ! !!$ SUBROUTINE BRKPAR(BRCOEF,ECOS,ESIN,AC2,SPCSIG,DEP2,RDX,RDY,IG) SUBROUTINE BRKPAR(BRCOEF,ECOS,ESIN,AC2,SPCSIG,DEP2,IG) ! (This subroutine has not been fully tested yet) ! !**************************************************************** ! USE SWCOMM3 USE SWCOMM4 USE OCPCOMM4 USE MOD_PREC USE ALL_VARS, ONLY : MT,NTSN,NBSN,VX,VY,ART2 # if defined (SPHERICAL) USE MOD_SPHERICAL, ONLY : TPI,DEG2RAD # endif ! IMPLICIT NONE ! ! 2. Purpose ! ! Determine the bottom slope in upwave direction and calculate ! the slope dependent breaking parameter according to Nelson (1987) ! Note that Nelson (1987) is used here since in Nelson (1994a,1994b) ! an error is present in the equation. ! ! 3. Method ! ! The breaker parameter is given by: ! ! Hm / d = 0.55 + 0.88 exp ( -0.012 * cot beta) ! ! with beta the angle the bed makes with the horizontal. This ! above equation is only valid for positive slopes (Negative ! slopes were not considered by Nelson. For very steep slopes ! (>0.05 say) a very large breaker parameter is obtained (>>1). ! ! To ensure wave breaking in laboratory cases (with very steep ! slopes an upper limit of 0.81 (which corresponds to a bottom ! slope of 0.01) is imposed on the model of Nelson. ! ! For negative bottom slopes (not considered by Nelson) a value ! op 0.73 is imposed (which is the average value in Table 2 of ! Battjes and Janssen (1978). ! REAL, INTENT(OUT) :: BRCOEF ! variable breaker coefficient REAL, INTENT(IN) :: SPCSIG(MSC) REAL, INTENT(IN) :: AC2(MDC,MSC,0:MT) ! action densities REAL, INTENT(IN) :: ECOS(MDC), ESIN(MDC) ! Cos and Sin of Theta REAL, INTENT(IN) :: DEP2(MT) ! depths at grid points !!$ REAL, INTENT(IN) :: RDX(10), RDY(10) INTEGER, INTENT(IN) :: IG ! 9. STRUCTURE ! ! ------------------------------------------------------------ ! Calculate total energy per direction for all frequencies ! determine action density per direction weighted with cos/sin ! determine mean propagation direction of energy ! ------------------------------------------------------------ ! calculate the depth derivative in the mean wave direction ! according to dd/ds (see also subroutine SPROSD) ! calculate the slope dependend breaking coefficient ! ------------------------------------------------------------ ! End of NELSON ! ------------------------------------------------------------- ! !************************************************************************ ! INTEGER :: ID ,IS ! counters INTEGER :: J,I1,I2 REAL(SP) :: F1 # if defined (SPHERICAL) REAL(DP) XTMP,XTMP1 # endif ! ! REAL :: ETOTS,EEX,EEY,EAD,SIGMA1,COSDIR,SINDIR,DDDX,DDDY,DDDS,DETOT ! ! *** determine the average wave direction *** ! EEX = 0. EEY = 0. ETOTS = 0. DO ID = 1, MDC EAD = 0. DO IS = 1, MSC SIGMA1 = SPCSIG(IS) DETOT = SIGMA1**2 * AC2(ID,IS,KCGRD(1)) EAD = EAD + DETOT END DO ETOTS = ETOTS + EAD EEX = EEX + EAD * ECOS(ID) EEY = EEY + EAD * ESIN(ID) END DO ! IF(ETOTS > 0.)THEN COSDIR = EEX / ETOTS SINDIR = EEY / ETOTS ELSE COSDIR = 1. SINDIR = 0. END IF ! ! *** Determine bottom slope in average wave propagation direction *** ! ! DDDX = RDX(1) * (DEP2(KCGRD(1)) - DEP2(KCGRD(2))) & ! + RDX(2) * (DEP2(KCGRD(1)) - DEP2(KCGRD(3))) ! DDDY = RDY(1) * (DEP2(KCGRD(1)) - DEP2(KCGRD(2))) & ! + RDY(2) * (DEP2(KCGRD(1)) - DEP2(KCGRD(3))) DDDX = 0.0_SP DDDY = 0.0_SP DO J=1,NTSN(IG)-1 I1 = NBSN(IG,J) I2 = NBSN(IG,J+1) F1 = 0.50_SP*(DEP2(I1)+DEP2(I2)) # if defined (SPHERICAL) DDDX=DDDX+F1*(VY(I1)-VY(I2))*TPI XTMP = VX(I2)*TPI-VX(I1)*TPI XTMP1 = VX(I2)-VX(I1) IF(XTMP1 > 180.0_SP)THEN XTMP = -360.0_SP*TPI+XTMP ELSE IF(XTMP1 < -180.0_SP)THEN XTMP = 360.0_SP*TPI+XTMP END IF DDDY=DDDY+F1*XTMP*COS(DEG2RAD*VY(IG)) # else DDDX = DDDX + F1*(VY(I1)-VY(I2)) DDDY = DDDY + F1*(VX(I2)-VX(I1)) # endif END DO DDDX = DDDX/ART2(IG) DDDY = DDDY/ART2(IG) ! DDDS = -1. * ( DDDX * COSDIR + DDDY * SINDIR ) ! ! *** calculate breaking coefficient according to Nelson (1987) *** ! IF(DDDS >= 0.)THEN DDDS = MAX ( 1.E-6 , DDDS) BRCOEF = PSURF(4) + PSURF(7) * EXP ( -PSURF(8) / DDDS ) BRCOEF = MIN ( PSURF(5) , BRCOEF ) ELSE BRCOEF = PSURF(6) END IF ! ! PSURF(2) = BRKVAR deleted ! PSURF(2) is no longer used to transmit br. coefficient RETURN END SUBROUTINE BRKPAR