|
XBeach
|
00001 module vsm_u_xb_module 00002 implicit none 00003 save 00004 private 00005 public vsm_u_XB 00006 contains 00007 subroutine vsm_u_XB ( ue ,ve ,h ,ks , & 00008 rho ,tauwx ,tauwy ,theta , & 00009 urms ,omega ,k ,diss , & 00010 e ,er ,fcvisc,facdel,ws , & 00011 fw ,facdf ,swglm ,n ,sigz , & 00012 uz ,vz ,ustz ,nutz , & 00013 ue_sed,ve_sed ) 00014 00015 implicit none 00016 00017 ! Input/output parameters 00018 real*8, intent(in) :: ue ! depth-averaged eulerian velocity, ksi comp 00019 real*8, intent(in) :: ve ! depth-averaged eulerian velocity, eta comp 00020 real*8, intent(in) :: h ! water depth 00021 real*8, intent(in) :: ks ! nikuradse roughness 00022 real*8, intent(in) :: rho ! density 00023 real*8, intent(in) :: tauwx ! wind shear stress, ksi comp 00024 real*8, intent(in) :: tauwy ! wind shear stress, eta comp 00025 real*8, intent(in) :: theta ! wave angle w.r.t. grid ksi-direction 00026 real*8, intent(in) :: urms ! orbital velocity 00027 real*8, intent(in) :: omega ! angular peak frequency 00028 real*8, intent(in) :: k ! wave number 00029 real*8, intent(in) :: diss ! breaking dissipation 00030 real*8, intent(in) :: e ! wave energy 00031 real*8, intent(in) :: er ! roller energy 00032 real*8, intent(in) :: fcvisc ! calibration factor breaking-induced viscosity (~0.05) 00033 real*8, intent(in) :: facdel ! calibration factor bbl thickness (~20) 00034 real*8, intent(in) :: facdf ! calibration factor bottom dissipation (~1) 00035 real*8, intent(in) :: ws ! average fall velocity 00036 real*8, intent(in) :: fw ! wave friction coefficient 00037 integer, intent(in) :: swGLM ! switch to use GLM formulation (1) or not (0) 00038 integer, intent(in) :: n ! number of vertical layers 00039 real*8, intent(out) :: sigz(n) ! vertical sigma grid (0 at bottom, 1 at surface) 00040 real*8, intent(out) :: uz(n) ! vertically distributed velocity, ksi comp. 00041 real*8, intent(out) :: vz(n) ! vertically distributed velocity, eta comp. 00042 real*8, intent(out) :: ustz(n) ! vertically distributed stokes drift in wave direction 00043 real*8, intent(out) :: nutz(n) ! vertically distributed viscosity 00044 real*8, intent(out) :: ue_sed ! advection velocity sediment (ksi-dir) 00045 real*8, intent(out) :: ve_sed ! advection velocity sediment (eta-dir) 00046 ! Local variables 00047 real*8 :: kappa,nut1,nut2,nut3,g,uorb 00048 real*8 :: nutda,nusur,nutb,numin 00049 real*8 :: hrms,z0,cw,aorb,delta,phiw 00050 real*8 :: Df,Dfx,Dfy,tauwav,Qw,Qwe,Qwx,Qwy 00051 real*8 :: ustokes,tausx,tausy,tswi,cf,sigma0,uold 00052 real*8 :: tbnw,dhdy,dhdx,sigmas,phis,phinub,sigs1 00053 real*8 :: facA,facA1,facG,facG1,facH,facH1 00054 real*8 :: facCx,facCy,facBx,facBy,facC1x,facC1y,facB1x,facB1y 00055 real*8 :: uxmean,uymean,uabs,err,sigmat,facln1,facln2 00056 real*8 :: vut,vvt,vud,vvd,Qstok,difstok,sigma,hd 00057 real*8 :: ustar,ast,rousepar 00058 real*8 :: crel(n),dsig(n) 00059 integer :: order_st=2,itmax,iter,i 00060 !------------------------------------------------------------------- 00061 !- Modified version Dirk-Jan Walstra (08-06-2007 --> copied from V205) 00062 !- It is assumed that the mass flux has a vertical distribution 00063 !- as predicted by 2nd or 3rd order Stokes drift. 00064 !- The bottom shear stress is assumed to react only to the Stokes 00065 !- drift at the bottom. This shear stress is used to construct the 00066 !- undertow profile etc. 00067 !- At the end the stokes drift is subtracted from the computed profile. 00068 !- The roller contribution to the mass flux is included in a depth 00069 !- averaged manner. 00070 !------------------------------------------------------------------- 00071 00072 itmax = 30 00073 numin = 1.d-6 00074 kappa = 0.41 00075 g = 9.81 00076 hrms = sqrt(8.*e/rho/g) 00077 z0 = ks/33. 00078 uorb =urms*sqrt(2.d0) 00079 cw=omega/k 00080 if (omega.gt.0.) then 00081 aorb = uorb/omega 00082 delta = 0.09*(aorb/ks)**0.82*ks/h 00083 ! fw = min(1.39 * (uorb/(omega*z0))**(-.52),0.3) 00084 endif 00085 00086 delta = facdel*max(delta,exp(1.)*z0/h) 00087 delta = min(delta,0.5) 00088 phiw = 6./delta/delta 00089 00090 Df = facDf*.283*rho*fw*uorb**3 00091 Dfx = Df*cos(theta) 00092 Dfy = Df*sin(theta) 00093 00094 ! Separate terms in depth-averaged viscosity 00095 00096 if (cw.gt.0.) then 00097 tauwav = diss/cw 00098 if (SWGLM.gt.0) then 00099 Qwe = (e+2.0*er)/cw/rho 00100 ! Determine stokes velocity at the bottom 00101 call stokes(hrms,omega,k,h,0.d0,ustokes,swglm) 00102 Qw = ustokes*h 00103 else 00104 Qw = (e+2.0*er)/cw/rho 00105 endif 00106 else 00107 tauwav=0. 00108 Qw =0. 00109 Qwe=0. 00110 endif 00111 00112 tausx = tauwx + tauwav*cos(theta) 00113 tausy = tauwy + tauwav*sin(theta) 00114 tswi = sqrt(tauwx**2+tauwy**2) 00115 Qwx = Qw*cos(theta) 00116 Qwy = Qw*sin(theta) 00117 00118 nut2 = h*kappa/3.*sqrt(tswi/rho) 00119 nut3 = fcvisc*hrms*(Diss/rho)**(1./3.) 00120 00121 ! Viscosity at surface level 00122 00123 nusur = 1.5*sqrt(nut2*nut2+nut3*nut3) 00124 nusur = max(nusur,numin) 00125 00126 ! Viscosity at bottom 00127 00128 cf = fw/2. 00129 nutb = cf*cf*uorb*uorb/omega 00130 00131 sigma0=z0/h 00132 uold=0. 00133 00134 ! Iteration for u en dhdx 00135 00136 do iter=1,itmax 00137 00138 tbnw = abs(rho*g*h*sqrt(dhdy*dhdy+dhdx*dhdx)) 00139 00140 nut1 = h*kappa/6.*sqrt(tbnw/rho) 00141 00142 ! Total depth-averaged viscosity 00143 00144 nutda = sqrt(nut1*nut1+nut2*nut2+nut3*nut3) 00145 nutda=max(nutda,numin) 00146 00147 sigmas = (nutda-nusur/3.)/(nutda-nusur/2.) 00148 phis = nusur/nutda/(sigmas-1.) 00149 00150 facA = h/(rho*phis*nutda) 00151 phinub = phis*nutda + phiw*nutb 00152 sigs1 = (phis*nutda*sigmas + phiw*nutb*delta)/ & 00153 (phis*nutda + phiw*nutb ) 00154 facA1 = h/(rho*phinub) 00155 00156 facG = facA/sigmas* ( & 00157 log(1./delta) + & 00158 (sigmas-1.)*log((sigmas-1.)/(sigmas-delta)) ) 00159 00160 facG1= facA1/sigs1* ( & 00161 log(delta/sigma0) + & 00162 (sigs1-1.)*log((sigs1-delta)/(sigs1-sigma0)) ) 00163 00164 facH = facA* ( & 00165 (sigmas-1.)*log((sigmas-1.)/(sigmas-delta)) + & 00166 (1.-delta) ) 00167 00168 facH1= facA1* ( & 00169 (sigs1-1.)*log((sigs1-delta)/(sigs1-sigma0)) + & 00170 (delta-sigma0) ) 00171 00172 ! facCx = (-Qwx/h-(facG1+facG)*tausx-(facG1-facH1/delta)*Dfx/cw)/ 00173 facCx = ( ue -(facG1+facG)*tausx-(facG1-facH1/delta)*Dfx/cw)/ & 00174 (facH1+facH-facG1-facG) 00175 dhdx = facCx/(rho*g*h) 00176 ! facCy = rho*g*h*dhdy 00177 facCy = ( ve -(facG1+facG)*tausy-(facG1-facH1/delta)*Dfy/cw)/ & 00178 (facH1+facH-facG1-facG) 00179 dhdy = facCy/(rho*g*h) 00180 00181 facBx = tausx-facCx 00182 facBy = tausy-facCy 00183 00184 facC1x= facCx-Dfx/delta/cw 00185 facC1y= facCy-Dfy/delta/cw 00186 00187 facB1x= facBx+Dfx/cw 00188 facB1y= facBy+Dfy/cw 00189 00190 uxmean = facG*facBx+facH*facCx+facG1*facB1x+facH1*facC1x 00191 uymean = facG*facBy+facH*facCy+facG1*facB1y+facH1*facC1y 00192 uabs=sqrt(max(uxmean*uxmean+uymean*uymean,1.d-6)) 00193 00194 err = abs((uold-uabs)/uabs) 00195 00196 if (err.lt.1.e-3) exit 00197 uold=uabs 00198 00199 enddo 00200 00201 if (err.gt.1.e-3) then 00202 stop 'no convergence in vsm_u' 00203 endif 00204 00205 ! end iteration 00206 00207 ! generate vertical grid 00208 00209 hd=1.d0 00210 do i=1,n 00211 sigz(i)=0.01d0*(hd/0.01d0)**(float(i-1)/float(n-1))/hd 00212 enddo 00213 00214 00215 ! determine velocity at top boundary layer and at e*z0 00216 00217 sigmat = 2.71828*sigma0 00218 00219 facln1=log(sigmat/sigma0) 00220 facln2=log((sigs1-sigmat)/(sigs1-sigma0)) 00221 00222 vut = facA1/sigs1 * ( & 00223 facB1x * facln1 - & 00224 (facB1x+facC1x*sigs1)* facln2 ) 00225 vvt = facA1/sigs1 * ( & 00226 facB1y * facln1 - & 00227 (facB1y+facC1y*sigs1)* facln2 ) 00228 00229 facln1=log(delta/sigma0) 00230 facln2=log((sigs1-delta)/(sigs1-sigma0)) 00231 00232 vud = facA1/sigs1 * ( & 00233 facB1x * facln1 - & 00234 (facB1x+facC1x*sigs1)* facln2 ) 00235 vvd = facA1/sigs1 * ( & 00236 facB1y * facln1 - & 00237 (facB1y+facC1y*sigs1)* facln2 ) 00238 00239 if (SWGLM.gt.0) then 00240 !- Construct vertical profile of 2nd or 3rd order stokes drift 00241 do i=1,n 00242 call stokes(hrms,omega,k,h,sigz(i)*h,ustz(i),order_st) 00243 enddo 00244 !- Determine mass flux of stokes drift 00245 Qstok=.5*sigz(1)*h*ustz(1) 00246 do i=2,n 00247 Qstok=Qstok+.5*h*((sigz(i)-sigz(i-1))*(ustz(i)+ustz(i-1))) 00248 enddo 00249 !- Determine difference between stokes mass flux and total mass flux 00250 !- including roller contribution 00251 difstok= Qwe - Qstok 00252 !- difstok=0. 00253 endif 00254 00255 ! Construct velocity profile for bottom layer (both sigma < sigmat and 00256 ! sigma > sigmat) and middle layer 00257 00258 do 10 i=1,n 00259 00260 sigma=sigz(i) 00261 00262 if (sigma .lt. sigmat) then 00263 00264 uz(i) = (sigma/sigmat)*vut 00265 vz(i) = (sigma/sigmat)*vvt 00266 nutz(i) = nutb 00267 00268 elseif (sigma .lt. delta) then 00269 00270 facln1=log(sigma/sigma0) 00271 facln2=log((sigs1-sigma)/(sigs1-sigma0)) 00272 00273 uz(i) = facA1/sigs1 * ( & 00274 facB1x * facln1 - & 00275 (facB1x+facC1x*sigs1) * facln2 ) 00276 vz(i) = facA1/sigs1 * ( & 00277 facB1y * facln1 - & 00278 (facB1y+facC1y*sigs1) * facln2 ) 00279 nutz(i) = phinub*sigma*(sigs1-sigma) 00280 00281 else 00282 00283 facln1=log(sigma/delta) 00284 facln2=log((sigmas-sigma)/(sigmas-delta)) 00285 00286 uz(i) = vud + facA/sigmas * ( & 00287 facBx * facln1 - & 00288 (facBx+facCx*sigmas)* facln2 ) 00289 vz(i) = vvd + facA/sigmas * ( & 00290 facBy * facln1 - & 00291 (facBy+facCy*sigmas)* facln2 ) 00292 nutz(i) = phis*nutda*sigma* ( sigmas - sigma ) 00293 00294 endif 00295 !- construct eulerian velocity profile (stokes corrected and roller 00296 !- contribution included as a depth averaged flux) 00297 !- Ueuler=Uglm-Ustokes (i.e. Veloc(4,i)=uz(i)-Qroller/h+Ustokesbottom) 00298 !- NOTE that uz(i) now is the GLM velocities which is also used for 00299 !- Suspended Transport Computations !!! 00300 if (SWGLM.gt.0) then 00301 uz(i) = uz(i) - (ustz(i)+difstok/h)*cos(theta) + Qwx/h 00302 vz(i) = vz(i) - (ustz(i)+difstok/h)*sin(theta) + Qwy/h 00303 endif 00304 00305 00306 10 continue 00307 00308 ! Construct representative velocity for sediment transport using Rouse profile 00309 ustar=sqrt(0.5d0*fw)*max(urms,1.e-6) 00310 rousepar=ws/kappa/ustar 00311 ast=ks/h 00312 do i=1,n 00313 if (sigz(i)>ast) then 00314 crel(i)=(sigz(i)/ast*(ast-1.d0)/(sigz(i)-1.d0))**(-rousepar) 00315 else 00316 crel(i)=1.d0 00317 endif 00318 enddo 00319 dsig(1)=0.5d0*sigz(2) 00320 do i=2,n-1 00321 dsig(i)=(sigz(i+1)-sigz(i-1))*0.5d0 00322 enddo 00323 dsig(n)=0.5d0*(sigz(n)-sigz(n-1)) 00324 if (crel(1)>1.d-6) then 00325 ue_sed=sum(dsig*uz*crel)/sum(dsig*crel) 00326 ve_sed=sum(dsig*vz*crel)/sum(dsig*crel) 00327 else 00328 ue_sed=uz(1) 00329 ve_sed=vz(1) 00330 endif 00331 00332 end subroutine vsm_u_XB 00333 00334 subroutine stokes (hrms,omega,k,h,z,ustokes,order_stk) 00335 real*8 :: hrms,omega,k,h,z,ustokes 00336 integer :: order_stk 00337 00338 real*8 :: a,f3,s 00339 a=hrms/2. 00340 if (omega.gt.0.) then 00341 if (order_stk.eq.2) then 00342 !-------------------------------------------------------------------- 00343 !--Determine Second order Stokes drift 00344 !-------------------------------------------------------------------- 00345 ustokes=omega*k*a**2*(cosh(2*k*z)/(2.*sinh(k*h)*sinh(k*h))) 00346 else 00347 !-------------------------------------------------------------------- 00348 !--Determine Third order Stokes drift 00349 !-------------------------------------------------------------------- 00350 f3 = 3./16.*(8.*(cosh(k*h))**6. + 1)/(sinh(k*h)**6.) 00351 s = 108.*(hrms/(.5*k*k*f3)) + 12.*sqrt(12./((.25*k*k*f3)**3.) & 00352 + 81.*((hrms/(.5*k*k*f3))**2.)) 00353 a = s**(1./3.)/6. - 2.*1./(.25*k*k*f3*s**(1./3.)) 00354 ustokes=.5*(omega/k)*cosh(2*k*z)*(k*a/sinh(k*h))**2* & 00355 (1.-(k*a/sinh(k*h))*cosh(k*z)) 00356 endif 00357 else 00358 ustokes=0. 00359 endif 00360 end subroutine stokes 00361 00362 end module vsm_u_xb_module
1.7.6.1
