1) bug fix *) a fix to synthesizing surface elevation from a wave spectrum with sub- and super-harmonic transfers --- SwashBCStokeswave.ftn90 2024-05-07 10:23:35.000000000 +0200 +++ SwashBCStokeswave.ftn90 2024-08-13 15:30:36.838505771 +0200 @@ -73,8 +73,12 @@ real :: fac ! multiplication factor real :: kd ! dimensionless still water depth real :: kwav ! wave number of a Fourier component + real :: omeg ! angular frequency of a Fourier component real :: sph ! sine of phase real :: s2ph ! sine of two times phase + real :: wn ! last angular frequency + ! + logical :: spectrum ! indicates whether wave spectrum is dealt with or not ! ! Structure ! @@ -84,6 +88,16 @@ ! if (ltrace) call strace (ient,'SwashBCStokeswave') ! + ! is Fourier series a spectrum? + ! + if ( istok == 3 .and. nfreq > 2 ) then + spectrum = .true. + else + spectrum = .false. + endif + ! + wn = bcfour%omega(nfreq) + ! if ( kmax == 1 ) then ! do j = 1, nfreq @@ -116,10 +130,12 @@ ! do j = 1, nfreq ! - ampl = bcfour%ampl(j) + ampl = bcfour%ampl (j) + omeg = bcfour%omega(j) kwav = kwave(ibgrpt,j) ! if ( .not. ampl /= 0. ) cycle + if ( spectrum .and. 2.*omeg > wn ) exit ! kd = kwav * swd ! @@ -186,10 +202,12 @@ ! do j = 1, nfreq ! - ampl = bcfour%ampl(j) + ampl = bcfour%ampl (j) + omeg = bcfour%omega(j) kwav = kwave(ibgrpt,j) ! if ( .not. ampl /= 0. ) cycle + if ( spectrum .and. 2.*omeg > wn ) exit ! kd = kwav * swd ! @@ -268,10 +286,12 @@ ! do j = 1, nfreq ! - ampl = bcfour%ampl(j) + ampl = bcfour%ampl (j) + omeg = bcfour%omega(j) kwav = kwave(ibgrpt,j) ! if ( .not. ampl /= 0. ) cycle + if ( spectrum .and. 2.*omeg > wn ) exit ! kd = kwav * swd ! @@ -360,10 +380,12 @@ ! do j = 1, nfreq ! - ampl = bcfour%ampl(j) + ampl = bcfour%ampl (j) + omeg = bcfour%omega(j) kwav = kwave(ibgrpt,j) ! if ( .not. ampl /= 0. ) cycle + if ( spectrum .and. 2.*omeg > wn ) exit ! kd = kwav * swd ! --- SwashBCtransferfnc.ftn90 2024-05-07 10:23:35.000000000 +0200 +++ SwashBCtransferfnc.ftn90 2024-08-13 15:30:36.954505801 +0200 @@ -30,10 +30,12 @@ ! Authors ! ! 1.00: Panagiotis Vasarmidis +! 10.05: Panagiotis Vasarmidis ! ! Updates ! ! 1.00, November 2023: New subroutine +! 10.05, August 2024: bug fix: exclude very small omega3 ! ! Purpose ! @@ -81,8 +83,9 @@ ! ! Parameter variables ! - real, parameter :: khmin = 0.4 ! minimum dimensionless depth - real, parameter :: khmax = 10. ! maximum dimensionless depth + real, parameter :: khmin = 0.2 ! minimum dimensionless depth + real, parameter :: khmax = 5. ! maximum dimensionless depth + real, parameter :: odmin = 0.02 ! minimum difference angular frequency ! ! Local variables ! @@ -113,6 +116,8 @@ real :: theta2 ! wave direction of second primary wave component with respect to computational coordinates real :: theta3 ! wave direction of bound wave component with respect to problem coordinates real, dimension(kmax) :: vel ! local 2nd order transfer function of layer-dependent velocity + real :: w0 ! first angular frequency + real :: wn ! last angular frequency ! ! Structure ! @@ -140,6 +145,8 @@ endif ! dw = bcfour%omega(2) - bcfour%omega(1) + w0 = bcfour%omega(1) + wn = bcfour%omega(nfreq) ! floop: do j = 1, nfreq ! @@ -275,116 +282,138 @@ ! if ( kd2 > khmin .and. kd2 < khmax .and. ampl2 /= 0. ) then ! - ! calculate frequency index + ! calculate wave frequency of bound sub-harmonic component ! omega3 = omega2 - omega1 ! - l = nint(omega3/dw) - ! - ! calculate bound wave number - ! - kwav3 = real(sqrt( dble(kwav1)**2 + dble(kwav2)**2 - 2.*dble(kwav1)*dble(kwav2)*cos(theta1-theta2) )) - ! - ! calculate direction of bound sub-harmonic component - ! - theta3 = alpc + atan( ( kwav2*sin(theta2) - kwav1*sin(theta1) ) / ( kwav2*cos(theta2) - kwav1*cos(theta1) ) ) - ! - ! calculate phase of bound sub-harmonic component - ! - phase3 = phase2 - phase1 - ! - ! include phase shift related to wave direction and wave number - ! - phase3 = phase3 + kwav3 * ( cos(theta3)*xp + sin(theta3)*yp ) - ! - ! compute sub-harmonic transfer functions - ! - if ( kmax == 1 ) then - eta = etasub1() - vel(1) = velsb11() - else if ( kmax == 2 ) then - eta = etasub2() - vel(1) = velsb12() - vel(2) = velsb22() - else if ( kmax == 3 ) then - eta = etasub3() - vel(1) = velsb13() - vel(2) = velsb23() - vel(3) = velsb33() - else if ( kmax == 4 ) then - eta = etasub4() - vel(1) = velsb14() - vel(2) = velsb24() - vel(3) = velsb34() - vel(4) = velsb44() - endif - ! - ! store surface elevation and velocity of sub-harmonic component - ! - rval = ampl1 * ampl2 - ! - subzc(ibgrpt,l) = subzc(ibgrpt,l) + eta * rval * cos( phase3 ) - subzs(ibgrpt,l) = subzs(ibgrpt,l) + eta * rval * sin( phase3 ) - ! - if ( oned .or. optg == 5 ) then - subuc(ibgrpt,l,:) = subuc(ibgrpt,l,:) + vel(:) * rval * cos( phase3 ) - subus(ibgrpt,l,:) = subus(ibgrpt,l,:) + vel(:) * rval * sin( phase3 ) - else - subuc(ibgrpt,l,:) = subuc(ibgrpt,l,:) + vel(:) * rval * cos( theta3 ) * cos( phase3 ) - subus(ibgrpt,l,:) = subus(ibgrpt,l,:) + vel(:) * rval * cos( theta3 ) * sin( phase3 ) - subvc(ibgrpt,l,:) = subvc(ibgrpt,l,:) + vel(:) * rval * sin( theta3 ) * cos( phase3 ) - subvs(ibgrpt,l,:) = subvs(ibgrpt,l,:) + vel(:) * rval * sin( theta3 ) * sin( phase3 ) - endif - ! - ! calculate direction of bound super-harmonic component - ! - theta3 = alpc + atan( ( kwav1*sin(theta1) + kwav2*sin(theta2) ) / ( kwav1*cos(theta1) + kwav2*cos(theta2) ) ) - ! - ! calculate phase of bound super-harmonic component - ! - phase3 = phase1 + phase2 - ! - ! include phase shift related to wave direction and wave number - ! - phase3 = phase3 + kwav3 * ( cos(theta3)*xp + sin(theta3)*yp ) - ! - ! compute super-harmonic transfer functions - ! - if ( kmax == 1 ) then - eta = etasup1() - vel(1) = velsp11() - else if ( kmax == 2 ) then - eta = etasup2() - vel(1) = velsp12() - vel(2) = velsp22() - else if ( kmax == 3 ) then - eta = etasup3() - vel(1) = velsp13() - vel(2) = velsp23() - vel(3) = velsp33() - else if ( kmax == 4 ) then - eta = etasup4() - vel(1) = velsp14() - vel(2) = velsp24() - vel(3) = velsp34() - vel(4) = velsp44() + if ( nfreq == 2 .or. .not. omega3 < odmin ) then + ! + ! calculate frequency index + ! + l = nint(omega3/dw) + ! + ! calculate wave number of bound sub-harmonic component + ! + kwav3 = real(sqrt( dble(kwav1)**2 + dble(kwav2)**2 - 2.*dble(kwav1)*dble(kwav2)*cos(theta1-theta2) )) + ! + ! calculate direction of bound sub-harmonic component + ! + theta3 = alpc + atan( ( kwav2*sin(theta2) - kwav1*sin(theta1) ) / ( kwav2*cos(theta2) - kwav1*cos(theta1) ) ) + ! + ! calculate phase of bound sub-harmonic component + ! + phase3 = phase2 - phase1 + ! + ! include phase shift related to wave direction and wave number + ! + phase3 = phase3 + kwav3 * ( cos(theta3)*xp + sin(theta3)*yp ) + ! + ! compute sub-harmonic transfer functions + ! + if ( kmax == 1 ) then + eta = etasub1() + vel(1) = velsb11() + else if ( kmax == 2 ) then + eta = etasub2() + vel(1) = velsb12() + vel(2) = velsb22() + else if ( kmax == 3 ) then + eta = etasub3() + vel(1) = velsb13() + vel(2) = velsb23() + vel(3) = velsb33() + else if ( kmax == 4 ) then + eta = etasub4() + vel(1) = velsb14() + vel(2) = velsb24() + vel(3) = velsb34() + vel(4) = velsb44() + endif + ! + ! store surface elevation and velocity of sub-harmonic component + ! + rval = ampl1 * ampl2 + ! + subzc(ibgrpt,l) = subzc(ibgrpt,l) + eta * rval * cos( phase3 ) + subzs(ibgrpt,l) = subzs(ibgrpt,l) + eta * rval * sin( phase3 ) + ! + if ( oned .or. optg == 5 ) then + subuc(ibgrpt,l,:) = subuc(ibgrpt,l,:) + vel(:) * rval * cos( phase3 ) + subus(ibgrpt,l,:) = subus(ibgrpt,l,:) + vel(:) * rval * sin( phase3 ) + else + subuc(ibgrpt,l,:) = subuc(ibgrpt,l,:) + vel(:) * rval * cos( theta3 ) * cos( phase3 ) + subus(ibgrpt,l,:) = subus(ibgrpt,l,:) + vel(:) * rval * cos( theta3 ) * sin( phase3 ) + subvc(ibgrpt,l,:) = subvc(ibgrpt,l,:) + vel(:) * rval * sin( theta3 ) * cos( phase3 ) + subvs(ibgrpt,l,:) = subvs(ibgrpt,l,:) + vel(:) * rval * sin( theta3 ) * sin( phase3 ) + endif + ! endif ! - ! store surface elevation and velocity of super-harmonic component - ! - rval = ampl1 * ampl2 + ! calculate wave frequency of bound super-harmonic component ! - supzc(ibgrpt,l) = supzc(ibgrpt,l) + eta * rval * cos( phase3 ) - supzs(ibgrpt,l) = supzs(ibgrpt,l) + eta * rval * sin( phase3 ) + omega3 = omega1 + omega2 ! - if ( oned .or. optg == 5 ) then - supuc(ibgrpt,l,:) = supuc(ibgrpt,l,:) + vel(:) * rval * cos( phase3 ) - supus(ibgrpt,l,:) = supus(ibgrpt,l,:) + vel(:) * rval * sin( phase3 ) - else - supuc(ibgrpt,l,:) = supuc(ibgrpt,l,:) + vel(:) * rval * cos( theta3 ) * cos( phase3 ) - supus(ibgrpt,l,:) = supus(ibgrpt,l,:) + vel(:) * rval * cos( theta3 ) * sin( phase3 ) - supvc(ibgrpt,l,:) = supvc(ibgrpt,l,:) + vel(:) * rval * sin( theta3 ) * cos( phase3 ) - supvs(ibgrpt,l,:) = supvs(ibgrpt,l,:) + vel(:) * rval * sin( theta3 ) * sin( phase3 ) + if ( nfreq == 2 .or. .not. omega3 > wn ) then + ! + ! calculate frequency index + ! + l = nint((omega3-2.*w0)/dw) + ! + ! calculate wave number of bound super-harmonic component + ! + kwav3 = real(sqrt( dble(kwav1)**2 + dble(kwav2)**2 + 2.*dble(kwav1)*dble(kwav2)*cos(theta1-theta2) )) + ! + ! calculate direction of bound super-harmonic component + ! + theta3 = alpc + atan( ( kwav1*sin(theta1) + kwav2*sin(theta2) ) / ( kwav1*cos(theta1) + kwav2*cos(theta2) ) ) + ! + ! calculate phase of bound super-harmonic component + ! + phase3 = phase1 + phase2 + ! + ! include phase shift related to wave direction and wave number + ! + phase3 = phase3 + kwav3 * ( cos(theta3)*xp + sin(theta3)*yp ) + ! + ! compute super-harmonic transfer functions + ! + if ( kmax == 1 ) then + eta = etasup1() + vel(1) = velsp11() + else if ( kmax == 2 ) then + eta = etasup2() + vel(1) = velsp12() + vel(2) = velsp22() + else if ( kmax == 3 ) then + eta = etasup3() + vel(1) = velsp13() + vel(2) = velsp23() + vel(3) = velsp33() + else if ( kmax == 4 ) then + eta = etasup4() + vel(1) = velsp14() + vel(2) = velsp24() + vel(3) = velsp34() + vel(4) = velsp44() + endif + ! + ! store surface elevation and velocity of super-harmonic component + ! + rval = ampl1 * ampl2 + ! + supzc(ibgrpt,l) = supzc(ibgrpt,l) + eta * rval * cos( phase3 ) + supzs(ibgrpt,l) = supzs(ibgrpt,l) + eta * rval * sin( phase3 ) + ! + if ( oned .or. optg == 5 ) then + supuc(ibgrpt,l,:) = supuc(ibgrpt,l,:) + vel(:) * rval * cos( phase3 ) + supus(ibgrpt,l,:) = supus(ibgrpt,l,:) + vel(:) * rval * sin( phase3 ) + else + supuc(ibgrpt,l,:) = supuc(ibgrpt,l,:) + vel(:) * rval * cos( theta3 ) * cos( phase3 ) + supus(ibgrpt,l,:) = supus(ibgrpt,l,:) + vel(:) * rval * cos( theta3 ) * sin( phase3 ) + supvc(ibgrpt,l,:) = supvc(ibgrpt,l,:) + vel(:) * rval * sin( theta3 ) * cos( phase3 ) + supvs(ibgrpt,l,:) = supvs(ibgrpt,l,:) + vel(:) * rval * sin( theta3 ) * sin( phase3 ) + endif + ! endif ! endif --- SwashInit.ftn90 2024-05-07 10:23:41.000000000 +0200 +++ SwashInit.ftn90 2024-08-13 15:30:42.954507438 +0200 @@ -78,7 +78,7 @@ VERTXT = BLANK VERNUM = 10.05 write (VERTXT, '(f5.2)') VERNUM - !call BUGFIX ('A') + call BUGFIX ('A') !call BUGFIX ('B') ! call OCPINI ( 'swashinit', .true., inerr ) --- SwashUpdateData.ftn90 2024-05-07 10:23:44.000000000 +0200 +++ SwashUpdateData.ftn90 2024-08-13 15:30:45.742508199 +0200 @@ -807,8 +807,8 @@ ! do j = 1, nfreq-1 ! - omegsb = curbfs%omega(1) - curbfs%omega(j+1) - omegsp = curbfs%omega(1) + curbfs%omega(j+1) + omegsb = curbfs%omega(j+1) - curbfs%omega(1) + omegsp = curbfs%omega(j+1) + curbfs%omega(1) ! fvalu = fvalu + subuc(i,j,1) * cos( omegsb*timco ) + subus(i,j,1) * sin( omegsb*timco ) & + supuc(i,j,1) * cos( omegsp*timco ) + supus(i,j,1) * sin( omegsp*timco ) @@ -913,8 +913,8 @@ ! do j = 1, nfreq-1 ! - omegsb = curbfs%omega(1) - curbfs%omega(j+1) - omegsp = curbfs%omega(1) + curbfs%omega(j+1) + omegsb = curbfs%omega(j+1) - curbfs%omega(1) + omegsp = curbfs%omega(j+1) + curbfs%omega(1) ! fvalu = fvalu + subuc(i,j,1) * cos( omegsb*timco ) + subus(i,j,1) * sin( omegsb*timco ) & + supuc(i,j,1) * cos( omegsp*timco ) + supus(i,j,1) * sin( omegsp*timco ) @@ -1061,8 +1061,8 @@ ! do j = 1, nfreq-1 ! - omegsb = curbfs%omega(1) - curbfs%omega(j+1) - omegsp = curbfs%omega(1) + curbfs%omega(j+1) + omegsb = curbfs%omega(j+1) - curbfs%omega(1) + omegsp = curbfs%omega(j+1) + curbfs%omega(1) ! fvalu = fvalu + subuc(i,j,k) * cos( omegsb*timco ) + subus(i,j,k) * sin( omegsb*timco ) & + supuc(i,j,k) * cos( omegsp*timco ) + supus(i,j,k) * sin( omegsp*timco ) @@ -1210,8 +1210,8 @@ ! do j = 1, nfreq-1 ! - omegsb = curbfs%omega(1) - curbfs%omega(j+1) - omegsp = curbfs%omega(1) + curbfs%omega(j+1) + omegsb = curbfs%omega(j+1) - curbfs%omega(1) + omegsp = curbfs%omega(j+1) + curbfs%omega(1) ! fvalu = fvalu + subuc(i,j,k) * cos( omegsb*timco ) + subus(i,j,k) * sin( omegsb*timco ) & + supuc(i,j,k) * cos( omegsp*timco ) + supus(i,j,k) * sin( omegsp*timco ) @@ -1603,8 +1603,8 @@ ! do j = 1, nfreq-1 ! - omegsb = curbfs%omega(1) - curbfs%omega(j+1) - omegsp = curbfs%omega(1) + curbfs%omega(j+1) + omegsb = curbfs%omega(j+1) - curbfs%omega(1) + omegsp = curbfs%omega(j+1) + curbfs%omega(1) ! fvalu = fvalu + ( subuc(i,j,1) + fac*subzc(i,j) ) * cos( omegsb*timco ) + ( subus(i,j,1) + fac*subzs(i,j) ) * sin( omegsb*timco ) & + ( supuc(i,j,1) + fac*supzc(i,j) ) * cos( omegsp*timco ) + ( supus(i,j,1) + fac*supzs(i,j) ) * sin( omegsp*timco ) @@ -1716,8 +1716,8 @@ ! do j = 1, nfreq-1 ! - omegsb = curbfs%omega(1) - curbfs%omega(j+1) - omegsp = curbfs%omega(1) + curbfs%omega(j+1) + omegsb = curbfs%omega(j+1) - curbfs%omega(1) + omegsp = curbfs%omega(j+1) + curbfs%omega(1) ! fvalu = fvalu + ( subuc(i,j,1) + fac*subzc(i,j) * cosbnd ) * cos( omegsb*timco ) + ( subus(i,j,1) + fac*subzs(i,j) * cosbnd ) * sin( omegsb*timco ) & + ( supuc(i,j,1) + fac*supzc(i,j) * cosbnd ) * cos( omegsp*timco ) + ( supus(i,j,1) + fac*supzs(i,j) * cosbnd ) * sin( omegsp*timco ) @@ -1868,8 +1868,8 @@ ! do j = 1, nfreq-1 ! - omegsb = curbfs%omega(1) - curbfs%omega(j+1) - omegsp = curbfs%omega(1) + curbfs%omega(j+1) + omegsb = curbfs%omega(j+1) - curbfs%omega(1) + omegsp = curbfs%omega(j+1) + curbfs%omega(1) ! fvalu = fvalu + ( subuc(i,j,k) + fac*subzc(i,j) ) * cos( omegsb*timco ) + ( subus(i,j,k) + fac*subzs(i,j) ) * sin( omegsb*timco ) & + ( supuc(i,j,k) + fac*supzc(i,j) ) * cos( omegsp*timco ) + ( supus(i,j,k) + fac*supzs(i,j) ) * sin( omegsp*timco ) @@ -2024,8 +2024,8 @@ ! do j = 1, nfreq-1 ! - omegsb = curbfs%omega(1) - curbfs%omega(j+1) - omegsp = curbfs%omega(1) + curbfs%omega(j+1) + omegsb = curbfs%omega(j+1) - curbfs%omega(1) + omegsp = curbfs%omega(j+1) + curbfs%omega(1) ! fvalu = fvalu + ( subuc(i,j,k) + fac*subzc(i,j) * cosbnd ) * cos( omegsb*timco ) + ( subus(i,j,k) + fac*subzs(i,j) * cosbnd ) * sin( omegsb*timco ) & + ( supuc(i,j,k) + fac*supzc(i,j) * cosbnd ) * cos( omegsp*timco ) + ( supus(i,j,k) + fac*supzs(i,j) * cosbnd ) * sin( omegsp*timco ) --- SwashUpdateUData.ftn90 2024-05-07 10:23:44.000000000 +0200 +++ SwashUpdateUData.ftn90 2024-08-13 15:30:45.946508254 +0200 @@ -560,8 +560,8 @@ ! do j = 1, nfreq-1 ! - omegsb = curbfs%omega(1) - curbfs%omega(j+1) - omegsp = curbfs%omega(1) + curbfs%omega(j+1) + omegsb = curbfs%omega(j+1) - curbfs%omega(1) + omegsp = curbfs%omega(j+1) + curbfs%omega(1) ! fvalu = fvalu + subuc(i,j,1) * cos( omegsb*timco ) + subus(i,j,1) * sin( omegsb*timco ) & + supuc(i,j,1) * cos( omegsp*timco ) + supus(i,j,1) * sin( omegsp*timco ) @@ -694,8 +694,8 @@ ! do j = 1, nfreq-1 ! - omegsb = curbfs%omega(1) - curbfs%omega(j+1) - omegsp = curbfs%omega(1) + curbfs%omega(j+1) + omegsb = curbfs%omega(j+1) - curbfs%omega(1) + omegsp = curbfs%omega(j+1) + curbfs%omega(1) ! fvalu = fvalu + subuc(i,j,k) * cos( omegsb*timco ) + subus(i,j,k) * sin( omegsb*timco ) & + supuc(i,j,k) * cos( omegsp*timco ) + supus(i,j,k) * sin( omegsp*timco ) @@ -925,8 +925,8 @@ ! do j = 1, nfreq-1 ! - omegsb = curbfs%omega(1) - curbfs%omega(j+1) - omegsp = curbfs%omega(1) + curbfs%omega(j+1) + omegsb = curbfs%omega(j+1) - curbfs%omega(1) + omegsp = curbfs%omega(j+1) + curbfs%omega(1) ! fvalu = fvalu + ( subuc(i,j,1) + fac*subzc(i,j) ) * cos( omegsb*timco ) + ( subus(i,j,1) + fac*subzs(i,j) ) * sin( omegsb*timco ) & + ( supuc(i,j,1) + fac*supzc(i,j) ) * cos( omegsp*timco ) + ( supus(i,j,1) + fac*supzs(i,j) ) * sin( omegsp*timco ) @@ -1063,8 +1063,8 @@ ! do j = 1, nfreq-1 ! - omegsb = curbfs%omega(1) - curbfs%omega(j+1) - omegsp = curbfs%omega(1) + curbfs%omega(j+1) + omegsb = curbfs%omega(j+1) - curbfs%omega(1) + omegsp = curbfs%omega(j+1) + curbfs%omega(1) ! fvalu = fvalu + ( subuc(i,j,k) + fac*subzc(i,j) ) * cos( omegsb*timco ) + ( subus(i,j,k) + fac*subzs(i,j) ) * sin( omegsb*timco ) & + ( supuc(i,j,k) + fac*supzc(i,j) ) * cos( omegsp*timco ) + ( supus(i,j,k) + fac*supzs(i,j) ) * sin( omegsp*timco )