subroutine vwave( C Write only - 1 fWaveSurf, fUnsteadyVel, fFluidAcc, 2 fUnsteadyPress, fUnsteadyDPressDZ, C Read/Write - 1 ScaleSteady, ScaleUnsteady, C Read only - 1 kStep, kInc, 2 nblock, ndim, nprops, naquaconst, nwindconst, 3 nstatevar, nfieldvar, iElemType, iLoadType, sname, 4 lUpdFluidVar, Coord, Velocity, StateVar, FieldVar, 5 DirVec, AquaSteadyConstants, WindConstants, 6 fSteadyVel, Props, dt, timeTotal, timeStep) C include 'vaba_param.inc' C parameter ( j_upd_FreeSurf = 0, 1 j_upd_FluidVarBuoyancy = 1, 2 j_upd_FluidVarDrag = 2, 3 j_upd_FluidVarInertia = 3, C The types of distributed loads: 1 j_lcr_PB = 51, 2 j_lcr_DragFDD = 53, j_lcr_DragWDD = 54, 3 j_lcr_DragFDT = 55, j_lcr_DragFI = 56, 4 j_lcr_DragFD1 = 57, j_lcr_DragFD2 = 58, 5 j_lcr_DragWD1 = 59, j_lcr_DragWD2 = 60, 6 j_lcr_DragFI1 = 61, j_lcr_DragFI2 = 62, C The types of concentrated loads: 1 j_ccr_TSB = 1002, j_ccr_DragTFD = 1004, 2 j_ccr_DragTWD = 1005, j_ccr_DragTSI = 1006 ) C parameter (nMaxTrains = 10, nConstPerTrain = 5) parameter (pi=3.14159265358979d0, twopi = 2.0d0 * pi) parameter (zero = 0.d0, two = 2.d0, one = 1.d0, wltol = 5.0d-3) parameter (i_v3d_X = 1, i_v3d_Y = 2, i_v3d_Z = 3) parameter (i_v2d_X = 1, i_v2d_Y = 2) parameter (pointnine = 0.9d0, r_sys_MaxVal = 1.0d+38) C dimension Props(nProps), Coord(nblock,ndim), 1 Velocity(nblock,ndim), StateVar(nblock,nstatevar), 2 FieldVar(nblock,nfieldvar), DirVec(nblock,ndim), 3 fWaveSurf(nblock), fFluidAcc(nblock,ndim), 4 fUnsteadyVel(nblock,ndim), fSteadyVel(nblock,ndim), 5 fUnsteadyPress(nblock), fUnsteadyDPressDZ(nblock), 6 AquaSteadyConstants(naquaconst), WindConstants(nwindconst) C character*80 sname C C Local Variables dimension Amplitude(nMaxTrains), RLambda(nMaxTrains), * Phaserad(nMaxTrains), DCosX(nMaxTrains), DCosY(nMaxTrains), * AiryWaveK1(nMaxTrains), AiryWaveK2(nMaxTrains), Tau(nMaxTrains) C c The first constant indicates whether Wave Period is specified (1) c or Wave Length is specified (0). lWavePeriod = int(Props(1)) if (lWavePeriod .ne. 0 .AND. lWavePeriod .ne. 1) then write(*,*) ' VWAVE: Error in Wave Period Flag (0 OR 1)' write(*,*) ' VWAVE: lWavePeriod = ', lWavePeriod write(*,*) ' VWAVE: Props(1) = ', Props(1) call xplb_exit end if * * Currently limited to max. 10 Airy Wave Trains (nMaxTrains) * Reads 5 constants per train (nConstPerTrain) nairywaves = (nprops-1) / nConstPerTrain if (nairywaves * nConstPerTrain .lt. nprops-1) then write(*,*) ' VWAVE: Error in number of properties' write(*,*) ' VWAVE: nairywaves = ', nairywaves, * ' nConstPerTrain = ', nConstPerTrain, * ' nprops = ', nprops call xplb_exit end if c c Initialize required variables from Steady Aqua Data zb = AquaSteadyConstants(1) zs = AquaSteadyConstants(2) grav = AquaSteadyConstants(3) dens = AquaSteadyConstants(4) c AquaDepth = zs - zb rnegdensgrav = -grav * dens rnegtwopigrav = -twopi * grav rnegtwopidensgrav = twopi * rnegdensgrav c floatmax = pointnine * r_sys_MaxVal rlnfloatmax = log(floatmax) c c Initialize required variables from Airy Wave Data if (lWavePeriod .eq. 0) then do i=1,nairywaves iprev = (i-1)*nConstPerTrain + 1 Amplitude(i) = Props(iprev + 1) RLambda(i) = Props(iprev + 2) rk = twopi / RLambda(i) rc = sqrt(grav * tanh(rk * AquaDepth) / rk) Tau(i) = twopi / (rc * rk) PhaseRad(i) = Props(iprev + 3) DCosX(i) = Props(iprev + 4) DCosY(i) = Props(iprev + 5) AiryWaveK1(i) = twopi / RLambda(i) AiryWaveK2(i) = twopi / Tau(i) write(*,*) '-------Train Number = ', i,'---------' write(*,*) ' Amplitude = ', Amplitude(i) write(*,*) ' Lambda = ', RLambda(i) write(*,*) ' Tau = ', Tau(i) write(*,*) ' PhaseRad = ', PhaseRad(i) write(*,*) ' DCosX = ', DCosX(i) write(*,*) ' DCosY = ', DCosY(i) write(*,*) ' AiryWaveK1 = ', AiryWaveK1(i) write(*,*) ' AiryWaveK2 = ', AiryWaveK2(i) end do elseif (lWavePeriod .eq. 1) then do i=1,nairywaves iprev = (i-1)*nConstPerTrain + 1 Amplitude(i) = Props(iprev + 1) Tau(i) = Props(iprev + 2) tausq2p = Tau(i) * Tau(i) / twopi wlength = grav * tausq2p twopih = twopi * AquaDepth do knewton = 1,50 tgthyp = tanh (twopih / wlength) cothyp = one / tgthyp flambda = tausq2p * tgthyp * grav / wlength - one dflambda = -tausq2p *(tgthyp + twopih * cothyp / wlength) 1 * grav / (wlength * wlength) wlength = wlength - flambda / dflambda if (abs(flambda) .lt. wltol) go to 200 end do write(*,*) ' VWAVE: Calculation for computing the'// * 'wave length from the wave period did not converge'// * 'in 50 newton steps... hence exiting' call xplb_exit 200 continue RLambda(i) = wlength PhaseRad(i) = Props(iprev + 3) DCosX(i) = Props(iprev + 4) DCosY(i) = Props(iprev + 5) AiryWaveK1(i) = twopi / RLambda(i) AiryWaveK2(i) = twopi / Tau(i) write(*,*) '-------Train Number = ', i,'---------' write(*,*) ' Amplitude = ', Amplitude(i) write(*,*) ' Lambda = ', RLambda(i) write(*,*) ' Tau = ', Tau(i) write(*,*) ' PhaseRad = ', PhaseRad(i) write(*,*) ' DCosX = ', DCosX(i) write(*,*) ' DCosY = ', DCosY(i) write(*,*) ' AiryWaveK1 = ', AiryWaveK1(i) write(*,*) ' AiryWaveK2 = ', AiryWaveK2(i) end do else write(*,*) ' VWAVE: Error in Wave Period Flag (0 OR 1)' write(*,*) ' VWAVE: lWavePeriod = ', lWavePeriod call xplb_exit end if C C The following if-loop and do-loop structure illustrates C proper usage of this user-subroutine. C if (lUpdFluidVar .eq. j_upd_FreeSurf) then C This part gets executed at the first call. C if (ndim .eq. 3) then do kn = 1,nblock c user coding to update fWaveSurf FreeSurfAiry = zero x = Coord(kn,i_v3d_X) y = Coord(kn,i_v3d_Y) do kw=1,nairywaves sd = x * DCosX(kw) + y * DCosY(kw) * the time term (cos and sine) tterm = sd * AiryWaveK1(kw) - * timeTotal * AiryWaveK2(kw) + PhaseRad(kw) ttermc = cos(tterm) FreeSurfAiry = FreeSurfAiry + Amplitude(kw) * ttermc end do fWaveSurf(kn) = zs - FreeSurfAiry end do elseif (ndim .eq. 2) then do kn = 1,nblock c user coding to update fWaveSurf FreeSurfAiry = zero x = Coord(kn,i_v2d_X) do kw=1,nairywaves sd = sign(x,DCosX(kw)) * the time term (cos and sine) tterm = sd * AiryWaveK1(kw) - * timeTotal * AiryWaveK2(kw) + PhaseRad(kw) ttermc = cos(tterm) FreeSurfAiry = FreeSurfAiry + Amplitude(kw) * ttermc end do fWaveSurf(kn) = zs - FreeSurfAiry end do end if c optionally update StateVar C else C This part gets executed at the second call. C if (lUpdFluidVar .eq. j_upd_FluidVarBuoyancy) then C Update variables for Buoyancy Loads (PB, TSB): if (ndim .eq. 3) then do kn = 1,nblock c user coding to update fUnsteadyPress, fUnsteadyDPressDZ press = zero dpressdz = zero do kw=1,nairywaves sd = Coord(kn,i_v3d_X) * DCosX(kw) * + Coord(kn,i_v3d_Y) * DCosY(kw) z = (Coord(kn,i_v3d_Z) - zb) * AiryWaveK1(kw) zdep = AquaDepth * AiryWaveK1(kw) * the time term (cos) tterm = (sd * AiryWaveK1(kw)) * - (timeTotal * AiryWaveK2(kw)) * + PhaseRad(kw) ttermc = cos(tterm) * the depth terms z = max(zero,z) zdep = max(zero,zdep) rZminusZdep = z - zdep rZminusZdep = min(rlnfloatmax,rZminusZdep) rExpNegTwoZ = exp(-z - z) A2cosh = one + rExpNegTwoZ A2sinh = one - rExpNegTwoZ B2 = one + exp(-zdep - zdep) rterm = exp(rZminusZdep) rterm = rterm / B2 rCoshZbyCoshZ0 = rterm * A2cosh rSinhZbyCoshZ0 = rterm * A2Sinh * the pressure components presscon = Amplitude(kw) * ttermc press = press + presscon * rCoshZbyCoshZ0 dpressdz = dpressdz + (presscon * rSinhZbyCoshZ0 / * RLambda(kw)) end do fUnsteadyPress(kn) = rnegdensgrav * press fUnsteadyDPressDZ(kn) = rnegtwopidensgrav * dpressdz end do elseif (ndim .eq. 2) then do kn = 1,nblock c user coding to update fUnsteadyPress, fUnsteadyDPressDZ press = zero dpressdz = zero do kw=1,nairywaves sd = sign(Coord(kn,i_v2d_X),DCosX(kw)) z = (Coord(kn,i_v2d_Y) - zb) * AiryWaveK1(kw) zdep = AquaDepth * AiryWaveK1(kw) * the time term (cos) tterm = (sd * AiryWaveK1(kw)) * - (timeTotal * AiryWaveK2(kw)) * + PhaseRad(kw) ttermc = cos(tterm) * the depth terms z = max(zero,z) zdep = max(zero,zdep) rZminusZdep = z - zdep rZminusZdep = min(rlnfloatmax,rZminusZdep) rExpNegTwoZ = exp(-z - z) A2cosh = one + rExpNegTwoZ A2sinh = one - rExpNegTwoZ B2 = one + exp(-zdep - zdep) rterm = exp(rZminusZdep) rterm = rterm / B2 rCoshZbyCoshZ0 = rterm * A2cosh rSinhZbyCoshZ0 = rterm * A2Sinh * the pressure components presscon = Amplitude(kw) * ttermc press = press + presscon * rCoshZbyCoshZ0 dpressdz = dpressdz + (presscon * rSinhZbyCoshZ0 / * RLambda(kw)) end do fUnsteadyPress(kn) = rnegdensgrav * press fUnsteadyDPressDZ(kn) = rnegtwopidensgrav * dpressdz end do end if c optionally update multipliers ScaleSteady, ScaleUnsteady else if (lUpdFluidVar .eq. j_upd_FluidVarDrag) then C Update variables for Drag Loads (FDD, FDT, FD1, FD2, TFD): if (ndim .eq. 3) then do kn = 1, nblock c user coding to update fUnsteadyVel vx = zero vy = zero vz = zero do kw=1,nairywaves sd = Coord(kn,i_v3d_X) * DCosX(kw) * + Coord(kn,i_v3d_Y) * DCosY(kw) z = (Coord(kn,i_v3d_Z) - zb) * AiryWaveK1(kw) zdep = AquaDepth * AiryWaveK1(kw) * the time term (cos and sine) tterm = sd * AiryWaveK1(kw) * - timeTotal * AiryWaveK2(kw) * + PhaseRad(kw) ttermc = cos(tterm) tterms = sin(tterm) * the depth terms z = max(zero,z) zdep = max(zero,zdep) rZminusZdep = z - zdep rZminusZdep = min(rlnfloatmax,rZminusZdep) rExpNegTwoZ = exp(-z - z) A2cosh = one + rExpNegTwoZ A2sinh = one - rExpNegTwoZ B2 = one + exp(-zdep - zdep) rterm = exp(rZminusZdep) rterm = rterm / B2 rCoshZbyCoshZ0 = rterm * A2cosh rSinhZbyCoshZ0 = rterm * A2Sinh * * the velocity components velcon = Amplitude(kw) * Tau(kw) / RLambda(kw) velcon1 = velcon * rCoshZbyCoshZ0 * ttermc vx = vx + DCosX(kw) * velcon1 vy = vy + DCosY(kw) * velcon1 vz = vz + velcon * rSinhZbyCoshZ0 * tterms end do fUnsteadyVel(kn,i_v3d_X) = -grav * vx fUnsteadyVel(kn,i_v3d_Y) = -grav * vy fUnsteadyVel(kn,i_v3d_Z) = -grav * vz end do elseif (ndim .eq. 2) then do kn = 1, nblock c user coding to update fUnsteadyVel vx = zero vy = zero do kw=1,nairywaves sd = sign(Coord(kn,i_v3d_X),DCosX(kw)) z = (Coord(kn,i_v2d_Y) - zb) * AiryWaveK1(kw) zdep = AquaDepth * AiryWaveK1(kw) * the time term (cos and sine) tterm = sd * AiryWaveK1(kw) * - timeTotal * AiryWaveK2(kw) * + PhaseRad(kw) ttermc = cos(tterm) tterms = sin(tterm) * the depth terms z = max(zero,z) zdep = max(zero,zdep) rZminusZdep = z - zdep rZminusZdep = min(rlnfloatmax,rZminusZdep) rExpNegTwoZ = exp(-z - z) A2cosh = one + rExpNegTwoZ A2sinh = one - rExpNegTwoZ B2 = one + exp(-zdep - zdep) rterm = exp(rZminusZdep) rterm = rterm / B2 rCoshZbyCoshZ0 = rterm * A2cosh rSinhZbyCoshZ0 = rterm * A2Sinh * * the velocity components velcon = Amplitude(kw) * Tau(kw) / RLambda(kw) velcon1 = velcon * rCoshZbyCoshZ0 * ttermc vx = vx + sign(velcon1, DCosX(kw)) vy = vy + velcon * rSinhZbyCoshZ0 * tterms end do fUnsteadyVel(kn,i_v2d_X) = -grav * vx fUnsteadyVel(kn,i_v2d_Y) = -grav * vy end do end if c optionally update multipliers ScaleSteady, ScaleUnsteady else if (lUpdFluidVar .eq. j_upd_FluidVarInertia) then C Update variables for Inertia Loads (FI, FI1, FI2, TSI): if (ndim .eq. 3) then do kn = 1,nblock c user coding to update fFluidAcc ax = zero ay = zero az = zero do kw=1,nairywaves sd = Coord(kn,i_v3d_X) * DCosX(kw) * + Coord(kn,i_v3d_Y) * DCosY(kw) z = (Coord(kn,i_v3d_Z) - zb) * AiryWaveK1(kw) zdep = AquaDepth * AiryWaveK1(kw) * the time term (cos and sine) tterm = sd * AiryWaveK1(kw) * - timeTotal * AiryWaveK2(kw) * + PhaseRad(kw) ttermc = cos(tterm) tterms = sin(tterm) * the depth terms z = max(zero,z) zdep = max(zero,zdep) rZminusZdep = z - zdep rZminusZdep = min(rlnfloatmax,rZminusZdep) rExpNegTwoZ = exp(-z - z) A2cosh = one + rExpNegTwoZ A2sinh = one - rExpNegTwoZ B2 = one + exp(-zdep - zdep) rterm = exp(rZminusZdep) rterm = rterm / B2 rCoshZbyCoshZ0 = rterm * A2cosh rSinhZbyCoshZ0 = rterm * A2Sinh * * the acceleration components aclcon = Amplitude(kw) / RLambda(kw) aclcon1 = aclcon * rCoshZbyCoshZ0 * tterms ax = ax + DCosX(kw) * aclcon1 ay = ay + DCosY(kw) * aclcon1 az = az + aclcon * rSinhZbyCoshZ0 * ttermc end do fFluidAcc(kn,i_v3d_X) = rnegtwopigrav * ax fFluidAcc(kn,i_v3d_Y) = rnegtwopigrav * ay fFluidAcc(kn,i_v3d_Z) = rnegtwopigrav * az end do elseif (ndim .eq. 2) then do kn = 1,nblock c user coding to update fFluidAcc c user coding to update fFluidAcc ax = zero ay = zero do kw=1,nairywaves sd = sign(Coord(kn,i_v2d_X),DCosX(kw)) z = (Coord(kn,i_v2d_Y) - zb) * AiryWaveK1(kw) zdep = AquaDepth * AiryWaveK1(kw) * the time term (cos and sine) tterm = sd * AiryWaveK1(kw) * - timeTotal * AiryWaveK2(kw) * + PhaseRad(kw) ttermc = cos(tterm) tterms = sin(tterm) * the depth terms (cosh, sinh and denom) z = max(zero,z) zdep = max(zero,zdep) rZminusZdep = z - zdep rZminusZdep = min(rlnfloatmax,rZminusZdep) rExpNegTwoZ = exp(-z - z) A2cosh = one + rExpNegTwoZ A2sinh = one - rExpNegTwoZ B2 = one + exp(-zdep - zdep) rterm = exp(rZminusZdep) rterm = rterm / B2 rCoshZbyCoshZ0 = rterm * A2cosh rSinhZbyCoshZ0 = rterm * A2Sinh * the acceleration components aclcon = Amplitude(kw) / RLambda(kw) aclcon1 = aclcon * rCoshZbyCoshZ0 * tterms ax = ax + sign(aclcon1,DCosX(kw)) ay = ay + aclcon * rSinhZbyCoshZ0 * ttermc end do fFluidAcc(kn,i_v2d_X) = rnegtwopigrav * ax fFluidAcc(kn,i_v2d_Y) = rnegtwopigrav * ay end do end if c optionally update multipliers ScaleSteady, ScaleUnsteady end if C end if C return end