|
XBeach
|
00001 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00002 ! Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University ! 00003 ! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski, ! 00004 ! Jaap van Thiel de Vries, Robert McCall ! 00005 ! ! 00006 ! d.roelvink@unesco-ihe.org ! 00007 ! UNESCO-IHE Institute for Water Education ! 00008 ! P.O. Box 3015 ! 00009 ! 2601 DA Delft ! 00010 ! The Netherlands ! 00011 ! ! 00012 ! This library is free software; you can redistribute it and/or ! 00013 ! modify it under the terms of the GNU Lesser General Public ! 00014 ! License as published by the Free Software Foundation; either ! 00015 ! version 2.1 of the License, or (at your option) any later version. ! 00016 ! ! 00017 ! This library is distributed in the hope that it will be useful, ! 00018 ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! 00019 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! 00020 ! Lesser General Public License for more details. ! 00021 ! ! 00022 ! You should have received a copy of the GNU Lesser General Public ! 00023 ! License along with this library; if not, write to the Free Software ! 00024 ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 ! 00025 ! USA ! 00026 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00027 module bedroughness_module 00028 00029 implicit none 00030 save 00031 private 00032 real*8,dimension(:,:),allocatable,private :: kru,krv,kru50,krv50,kru90,krv90 00033 real*8,dimension(:,:),allocatable,private :: urms_upd,u2_upd 00034 real*8,dimension(:,:),allocatable,private :: facbl,blphi,infilb,Ubed,Ventilation 00035 real*8,dimension(:,:),allocatable,private :: ueuf,uevf,vevf,veuf 00036 real*8,dimension(:,:),allocatable,private :: ueuold,uevold,vevold,veuold 00037 real*8,dimension(:,:),allocatable,private :: dudtsmooth,dvdtsmooth 00038 real*8,dimension(:,:),allocatable,private :: shieldsu,shieldsv 00039 real*8,dimension(2), private :: dtold 00040 real*8, private :: delta,rhogdelta 00041 logical,dimension(:,:),allocatable,private :: isveggie 00042 double precision,dimension(:,:),allocatable,private :: nman0 00043 00044 double precision, private :: nmanmin = -123 00045 public bedroughness_init 00046 public bedroughness_update 00047 00048 00049 contains 00050 00051 subroutine bedroughness_init(s,par) 00052 use params, only: parameters 00053 use spaceparamsdef, only: spacepars 00054 use paramsconst 00055 00056 implicit none 00057 00058 type(parameters),intent(in) :: par 00059 type(spacepars),target :: s 00060 integer :: i,j 00061 00062 !include 's.ind' 00063 !include 's.inp' 00064 00065 s%taubx_add = 0.d0 00066 s%tauby_add = 0.d0 00067 00068 select case (par%bedfriction) 00069 00070 case(BEDFRICTION_CHEZY) 00071 s%cfu = par%g/s%bedfriccoef**2 00072 s%cfv = par%g/s%bedfriccoef**2 00073 case(BEDFRICTION_CF) 00074 s%cfu = s%bedfriccoef 00075 s%cfv = s%bedfriccoef 00076 case (BEDFRICTION_MANNING) 00077 00078 if ( par%dynamrough .eq. 1 ) then 00079 allocate (nman0(s%nx+1,s%ny+1)) 00080 allocate (isveggie(s%nx+1,s%ny+1)) 00081 nman0=s%bedfriccoef 00082 nmanmin=minval( s%bedfriccoef ) 00083 where (s%bedfriccoef>nmanmin) 00084 isveggie=.true. 00085 elsewhere 00086 isveggie=.false. 00087 endwhere 00088 endif 00089 00090 where(s%wetu==1) 00091 s%cfu = par%g*s%bedfriccoef**2/s%hu**(1.d0/3) 00092 s%cfu = min(s%cfu,par%maxcf) 00093 elsewhere 00094 s%cfu = par%maxcf 00095 endwhere 00096 where(s%wetv==1) 00097 s%cfv = par%g*s%bedfriccoef**2/s%hv**(1.d0/3) 00098 s%cfv= min(s%cfv,par%maxcf) 00099 elsewhere 00100 s%cfv = par%maxcf 00101 endwhere 00102 case (BEDFRICTION_WHITE_COLEBROOK_GRAINSIZE) 00103 allocate (kru(s%nx+1,s%ny+1)) 00104 allocate (krv(s%nx+1,s%ny+1)) 00105 ! compute initial D90 at cell interface 00106 if (par%ngd==1) then 00107 kru = 3.d0*s%D90top 00108 krv = 3.d0*s%D90top 00109 else 00110 ! Rougness height in u-points 00111 do j=1,s%ny+1 00112 do i=1,s%nx 00113 kru(i,j) = 1.5d0*(s%D90top(i,j)+s%D90top(i+1,j)) 00114 enddo 00115 enddo 00116 kru(s%nx+1,:) = kru(s%nx,:) 00117 ! Roughness height in v-points 00118 if (par%ny==0) then 00119 ! v-point is central on top of z-point 00120 krv = 3*s%D90top 00121 else 00122 do j=1,s%ny 00123 do i=1,s%nx+1 00124 krv(i,j) = 1.5d0*(s%D90top(i,j)+s%D90top(i,j+1)) 00125 enddo 00126 enddo 00127 krv(:,s%ny+1) = krv(:,s%ny) 00128 endif 00129 endif 00130 where(s%wetu==1) 00131 s%cfu = par%g/(18*log10(12*max(s%hu,kru)/kru))**2 00132 s%cfu = min(s%cfu,par%maxcf) 00133 elsewhere 00134 s%cfu = par%maxcf 00135 endwhere 00136 where(s%wetv==1) 00137 s%cfv = par%g/(18*log10(12*max(s%hv,krv)/krv))**2 00138 s%cfv = min(s%cfv,par%maxcf) 00139 elsewhere 00140 s%cfv = par%maxcf 00141 endwhere 00142 case (BEDFRICTION_WHITE_COLEBROOK) 00143 allocate (kru(s%nx+1,s%ny+1)) 00144 allocate (krv(s%nx+1,s%ny+1)) 00145 kru = s%bedfriccoef 00146 krv = s%bedfriccoef 00147 where(s%wetu==1) 00148 s%cfu = par%g/(18*log10(12*max(s%hu,kru)/kru))**2 00149 s%cfu = min(s%cfu,par%maxcf) 00150 elsewhere 00151 s%cfu = par%maxcf 00152 endwhere 00153 where(s%wetv==1) 00154 s%cfv = par%g/(18*log10(12*max(s%hv,krv)/krv))**2 00155 s%cfv = min(s%cfv,par%maxcf) 00156 elsewhere 00157 s%cfv = par%maxcf 00158 endwhere 00159 end select 00160 00161 00162 00163 end subroutine bedroughness_init 00164 00165 00166 subroutine bedroughness_update(s,par) 00167 00168 use params, only: parameters 00169 use spaceparamsdef, only: spacepars 00170 use paramsconst 00171 00172 implicit none 00173 00174 type(parameters),intent(in) :: par 00175 type(spacepars),target :: s 00176 integer :: i,j 00177 00178 00179 00180 select case (par%bedfriction) 00181 00182 case(BEDFRICTION_CHEZY) 00183 ! do nothing, this is constant bed roughness, already set in initialise 00184 case(BEDFRICTION_CF) 00185 ! do nothing, this is constant bed roughness, already set in initialise 00186 case (BEDFRICTION_MANNING) 00187 ! C = H**(1/6)/n 00188 ! cf = g/C**2 = g/(hu**(1/6)/n)**2 = g*n**2/hu**(1/3) 00189 00190 if (par%dynamrough .eq. 1) then 00191 ! ML 2 KN changed back to one where loop instead of 3 separate to ensure just 1 of the 3 statements is executed per cell. 00192 where ((isveggie) .and. (s%sedero .gt. 0d0)) ! linear function do to deposition ( sedero > 0 ) 00193 s%bedfriccoef = nmanmin + min( max( (par%dstem-s%sedero)/par%dstem , 0.d0), 1.0d0) * (nman0 - nmanmin) 00194 elsewhere ((isveggie) .and. (s%sedero .lt. (par%droot*-1) ) ) ! step function do to erosion larger than root than always minimum ( sedero < -droot ) 00195 s%bedfriccoef = nmanmin 00196 isveggie = .false. 00197 elsewhere ! linear function do to deposition ( -droot < sedero < 0 ) 00198 s%bedfriccoef = nmanmin + min( max( (par%droot+s%sedero)/par%droot , 0.d0), 1.0d0) * (nman0 - nmanmin) 00199 endwhere 00200 endif 00201 00202 where(s%wetu==1) 00203 s%cfu = par%g*s%bedfriccoef**2/s%hu**(1.d0/3) 00204 s%cfu = min(s%cfu,par%maxcf) 00205 elsewhere 00206 s%cfu = par%maxcf 00207 endwhere 00208 where(s%wetv==1) 00209 s%cfv = par%g*s%bedfriccoef**2/s%hv**(1.d0/3) 00210 s%cfv = min(s%cfv,par%maxcf) 00211 elsewhere 00212 s%cfv = par%maxcf 00213 endwhere 00214 case (BEDFRICTION_WHITE_COLEBROOK_GRAINSIZE) 00215 if (par%ngd>1) then 00216 ! Rougness height in u-points 00217 do j=1,s%ny+1 00218 do i=1,s%nx 00219 if (s%wetu(i,j)==1) then 00220 kru(i,j) = 1.5d0*(s%D90top(i,j)+s%D90top(i+1,j)) 00221 endif 00222 enddo 00223 enddo 00224 kru(s%nx+1,:) = kru(s%nx,:) 00225 ! Roughness height in v-points 00226 if (par%ny==0) then 00227 ! v-point is central on top of z-point 00228 where (s%wetv==1) 00229 krv = 3*s%D90top 00230 endwhere 00231 else 00232 do j=1,s%ny 00233 do i=1,s%nx+1 00234 if (s%wetv(i,j)==1) then 00235 krv(i,j) = 1.5d0*(s%D90top(i,j)+s%D90top(i,j+1)) 00236 endif 00237 enddo 00238 enddo 00239 krv(:,s%ny+1) = krv(:,s%ny) 00240 endif 00241 endif 00242 where(s%wetu==1) 00243 s%cfu = par%g/(18*log10(12*max(s%hu,kru)/kru))**2 00244 s%cfu = min(s%cfu,0.1d0) 00245 elsewhere 00246 s%cfu = par%maxcf 00247 endwhere 00248 where(s%wetv==1) 00249 s%cfv = par%g/(18*log10(12*max(s%hv,krv)/krv))**2 00250 s%cfv = min(s%cfv,par%maxcf) 00251 elsewhere 00252 s%cfv = par%maxcf 00253 endwhere 00254 case (BEDFRICTION_WHITE_COLEBROOK) 00255 where(s%wetu==1) 00256 s%cfu = par%g/(18*log10(12*max(s%hu,kru)/kru))**2 00257 s%cfu = min(s%cfu,par%maxcf) 00258 elsewhere 00259 s%cfu = par%maxcf 00260 endwhere 00261 where(s%wetv==1) 00262 s%cfv = par%g/(18*log10(12*max(s%hv,krv)/krv))**2 00263 s%cfv = min(s%cfv,par%maxcf) 00264 elsewhere 00265 s%cfv = par%maxcf 00266 endwhere 00267 end select 00268 00269 ! Additional bed friction terms due to acceleration 00270 s%taubx_add = 0.d0 00271 s%tauby_add = 0.d0 00272 00273 ! Increased bed shear due to acceleration terms (phase-shift type or Morison type) 00274 if (par%friction_acceleration==CF_ACC_NIELSEN) then 00275 call acceleration_boundary_layer_effect_Nielsen(s,par) 00276 elseif (par%friction_acceleration==CF_ACC_MCCALL) then 00277 call acceleration_boundary_layer_effect(s,par) 00278 endif 00279 00280 ! turbulence and acceleration add to shear stress, but do 00281 ! not necessarily also need enhanced cf from infiltration 00282 if (par%friction_turbulence==1) then 00283 call turbulence_boundary_layer_effect(s,par) 00284 endif 00285 00286 ! Infiltration enhances normal bed shear stress through 00287 ! increase of cfu and cfv 00288 if (par%friction_infiltration==1) then 00289 call infiltration_boundary_layer_effect(s,par) 00290 endif 00291 00292 00293 s%cfu = min(s%cfu,1.d0) ! check why needed 00294 s%cfv = min(s%cfv,1.d0) 00295 00296 00297 end subroutine bedroughness_update 00298 00299 subroutine infiltration_boundary_layer_effect(s,par) 00300 00301 use params, only: parameters 00302 use spaceparamsdef, only: spacepars 00303 00304 IMPLICIT NONE 00305 00306 type(parameters),intent(in) :: par 00307 type(spacepars) :: s 00308 integer :: i,j 00309 real*8,parameter :: epsVentilation = 1.d0 00310 real*8,parameter :: maxEnhancement1 = 3.d0 ! Pers. Corr. Conley 2014 00311 real*8,parameter :: maxEnhancement2 = 0.1d0 ! Pers. Corr. Conley 2014 00312 00313 00314 if (.not.allocated(facbl)) then 00315 ! arrays available in bed friction module only 00316 allocate (facbl(s%nx+1,s%ny+1)) 00317 allocate (blphi(s%nx+1,s%ny+1)) 00318 allocate (infilb(s%nx+1,s%ny+1)) 00319 allocate (Ubed(s%nx+1,s%ny+1)) 00320 allocate (Ventilation(s%nx+1,s%ny+1)) 00321 endif 00322 00323 00324 ! Include infiltration effect on boundary layer thinning 00325 ! Taken from eq 13 Butt, Russell, Turner 2001, based on 00326 ! Conley and Inman 1994 00327 ! In C&I derivation f = 2*cf and w is positive up, and b=0.9-2.0 so modified to: 00328 ! Tau_w/Tau_0 = C/(2*cf) * -infil/abs(u) 00329 ! where C is a constant of ~2 00330 ! 00331 ! Can re-write to taubx_total = taubx + taubx_add 00332 ! 00333 ! Tau ~ cf 00334 ! cf_w/cf_0 = Tau_w/Tau_0 00335 ! 00336 ! Therefore: taubx_total = cf_0*rho*u*abs(u) + (1-cf_w/cf_0)*cf*rho*u*abs(u) 00337 ! (here we neglect the urms part to the drag in taubx_add for surf-beat, as intrawave velocities and intrawave infiltration 00338 ! not resolved in surf-beat/stationary type approach) 00339 00340 ! U-points 00341 do j=1,s%ny+1 00342 do i=1,s%nx 00343 if(s%wetu(i,j)==1) then 00344 infilb(i,j) = 0.5d0*(s%infil(i,j)+s%infil(i+1,j)) 00345 else 00346 infilb(i,j) = 0.d0 00347 endif 00348 enddo 00349 enddo 00350 infilb(s%nx+1,:) = infilb(s%nx,:) 00351 where(s%wetu==1) 00352 !!!Ubed = sqrt(cfu)*abs(vmageu) older 00353 !Ubed = abs(s%uu) ! original XBeach-G 00354 Ubed = abs(s%vmageu) 00355 Ventilation = -infilb/max(Ubed,1.d-6) 00356 Ventilation = max(min(Ventilation,epsVentilation),-epsVentilation) 00357 blphi = 0.9d0/2/s%cfu*Ventilation 00358 ! need to account for devision by zero 00359 where(infilb>=0.d0) 00360 blphi = min(blphi,-1d-4) 00361 facbl = blphi/(exp(blphi)-1.d0) 00362 facbl = min(facbl,maxEnhancement1) 00363 elsewhere 00364 blphi = max(blphi,1d-4) 00365 facbl = blphi/(exp(blphi)-1.d0) 00366 facbl = max(facbl,maxEnhancement2) 00367 endwhere 00368 !s%cfu = s%cfu*facbl 00369 !s%cfu = min(s%cfu,1.d0) 00370 ! backward compatibility maximum value (almost never met) 00371 facbl = min(facbl,1.d0/s%cfu) ! maximizes cfu at 1.0 (=huge) 00372 s%taubx_add = s%taubx_add + (facbl-1.d0)*s%cfu*par%rho*s%ueu*s%vmageu 00373 endwhere 00374 00375 00376 ! V-points 00377 if(par%ny>0) then 00378 do j=1,s%ny 00379 do i=1,s%nx+1 00380 if(s%wetv(i,j)==1) then 00381 infilb(i,j) = 0.5d0*(s%infil(i,j)+s%infil(i,j+1)) 00382 else 00383 infilb(i,j) = 0.d0 00384 endif 00385 enddo 00386 enddo 00387 infilb(:,s%ny+1) = s%infil(:,s%ny) 00388 else 00389 ! v-points centered on cell centre 00390 infilb = s%infil 00391 endif 00392 where(s%wetv==1) 00393 !!!Ubed = sqrt(cfv)*abs(vmagev) ! older 00394 ! Ubed = abs(vv) ! orginal XBeach-G 00395 Ubed = abs(s%vmagev) 00396 Ventilation = -infilb/max(Ubed,1.d-6) 00397 Ventilation = max(min(Ventilation,epsVentilation),-epsVentilation) 00398 blphi = 0.9d0/2/s%cfv*Ventilation 00399 ! need to account for devision by zero 00400 where(infilb>=0.d0) 00401 blphi = min(blphi,-1d-4) 00402 facbl = blphi/(exp(blphi)-1.d0) 00403 facbl = min(facbl,maxEnhancement1) 00404 elsewhere 00405 blphi = max(blphi,1d-4) 00406 facbl = blphi/(exp(blphi)-1.d0) 00407 facbl = max(facbl,maxEnhancement2) 00408 endwhere 00409 !s%cfv = s%cfv*facbl 00410 !s%cfv = min(s%cfv,1.d0) 00411 facbl = min(facbl,1.d0/s%cfv) 00412 s%tauby_add = s%tauby_add + (facbl-1.d0)*s%cfv*par%rho*s%vev*s%vmagev 00413 endwhere 00414 00415 end subroutine infiltration_boundary_layer_effect 00416 00417 subroutine acceleration_boundary_layer_effect(s,par) 00418 00419 use params, only: parameters 00420 use spaceparamsdef, only: spacepars 00421 use paramsconst 00422 00423 IMPLICIT NONE 00424 00425 type(parameters),intent(in) :: par 00426 type(spacepars) :: s 00427 00428 real*8 :: Tsmooth,factime 00429 real*8,dimension(:,:),allocatable,save :: Fi 00430 00431 00432 if(.not.allocated(dudtsmooth)) then 00433 ! arrays available in bed friction module only 00434 allocate (dudtsmooth(s%nx+1,s%ny+1)) 00435 allocate (dvdtsmooth(s%nx+1,s%ny+1)) 00436 allocate (ueuold (s%nx+1,s%ny+1)) 00437 allocate (uevold (s%nx+1,s%ny+1)) 00438 allocate (vevold (s%nx+1,s%ny+1)) 00439 allocate (veuold (s%nx+1,s%ny+1)) 00440 allocate (ueuf (s%nx+1,s%ny+1)) 00441 allocate (uevf (s%nx+1,s%ny+1)) 00442 allocate (vevf (s%nx+1,s%ny+1)) 00443 allocate (veuf (s%nx+1,s%ny+1)) 00444 allocate (Fi (s%nx+1,s%ny+1)) 00445 ueuold = 0.d0 00446 uevold = 0.d0 00447 vevold = 0.d0 00448 veuold = 0.d0 00449 ueuf = 0.d0 00450 uevf = 0.d0 00451 vevf = 0.d0 00452 veuf = 0.d0 00453 dtold = par%dt 00454 endif 00455 00456 ! Problem with derivation of DUdt in U and V points 00457 ! ueu and veu are not known until end of flow timestep, so 00458 ! both are based on previous time step. The temporal gradient 00459 ! can only be found by going back one more time step. 00460 ! We need to keep track of par%dt so the correct time step for 00461 ! the gradient is used: 00462 ! 00463 ! i-2 i-1 i (now) i+1 (future) 00464 ! dtold(1) dtold(2) par%dt 00465 ! ueuold ueuf unknown 00466 ! 00467 ! Note that this assumes that the acceleration term does not change significantly 00468 ! over 2 timesteps and that the location mask wetu can still be used, 00469 ! even though it is not the wetu mask of 2 timesteps ago. 00470 Tsmooth = par%Trep/20 00471 factime = min(dtold(1)/Tsmooth,1.d0) 00472 00473 where(s%wetu==1) 00474 ueuf = (1-factime)*ueuf + factime*s%ueu ! note, because taub_add has linear relation with total acceleration (dvmag/dt) 00475 ! the x-component taubx_add can be computed from only the u-component of the 00476 ! acceleration 00477 dudtsmooth = (ueuf-ueuold)/dtold(1) 00478 !dudtsmooth = (ueuf-ueuold)/dtold(1)+ududx*wetu_ududx ! Following Baldock attempt to correct du/dt to dp/dx 00479 ! cut off rediculous values 00480 where (dudtsmooth>100*par%g) 00481 dudtsmooth = 100*par%g 00482 elsewhere (dudtsmooth<-100*par%g) 00483 dudtsmooth = -100*par%g 00484 endwhere 00485 s%taubx_add = s%taubx_add + par%ci*par%rho*min(s%D50top,s%hu)*dudtsmooth 00486 elsewhere 00487 ueuf = 0.d0 00488 veuf = 0.d0 00489 endwhere 00490 00491 ! now set the old variable to the i-1 variables 00492 ueuold = ueuf 00493 veuold = veuf 00494 00495 ! in some cases we don't need to compute v-acceleration 00496 ! - 1D non-h has only u components 00497 ! - 1D surfbeat / stationary has no v-components if boundary set to wall 00498 if (par%ny==0 .and. ((par%left==LR_WALL .and. par%right==LR_WALL) .or. par%swave==0)) then 00499 ! do nothing (don't change tauby_add) 00500 else 00501 where(s%wetv==1) 00502 vevf = (1-factime)*vevf + factime*s%vev 00503 dvdtsmooth = (vevf-vevold)/dtold(1) 00504 where (dvdtsmooth>100*par%g) 00505 dvdtsmooth = 100*par%g 00506 elsewhere (dvdtsmooth<-100*par%g) 00507 dvdtsmooth = -100*par%g 00508 endwhere 00509 s%tauby_add = s%tauby_add + par%ci*par%rho*min(s%D50top,s%hv)*dvdtsmooth 00510 elsewhere 00511 vevf = 0.d0 00512 uevf = 0.d0 00513 endwhere 00514 00515 ! now set the old variable to the i-1 variables 00516 uevold = uevf 00517 vevold = vevf 00518 endif 00519 00520 ! update timekeeping 00521 dtold(1) = dtold(2) 00522 dtold(2) = par%dt 00523 00524 end subroutine acceleration_boundary_layer_effect 00525 00526 subroutine acceleration_boundary_layer_effect_Nielsen(s,par) 00527 00528 use params, only: parameters 00529 use spaceparamsdef, only: spacepars 00530 use constants, only: pi 00531 use paramsconst 00532 00533 IMPLICIT NONE 00534 00535 type(parameters),intent(in) :: par 00536 type(spacepars) :: s 00537 real*8 :: omegap,Tsmooth,factime,iomegap 00538 real*8,save :: phirad 00539 00540 00541 if(.not.allocated(dudtsmooth)) then 00542 allocate (dudtsmooth(s%nx+1,s%ny+1)) 00543 allocate (dvdtsmooth(s%nx+1,s%ny+1)) 00544 allocate (ueuold (s%nx+1,s%ny+1)) 00545 !allocate (uevold (s%nx+1,s%ny+1)) 00546 allocate (vevold (s%nx+1,s%ny+1)) 00547 !allocate (veuold (s%nx+1,s%ny+1)) 00548 allocate (ueuf (s%nx+1,s%ny+1)) 00549 !allocate (uevf (s%nx+1,s%ny+1)) 00550 allocate (vevf (s%nx+1,s%ny+1)) 00551 !allocate (veuf (s%nx+1,s%ny+1)) 00552 ueuold = 0.d0 00553 !uevold = 0.d0 00554 vevold = 0.d0 00555 !veuold = 0.d0 00556 ueuf = 0.d0 00557 !uevf = 0.d0 00558 vevf = 0.d0 00559 !veuf = 0.d0 00560 dtold = par%dt 00561 phirad = par%phit/180*par%px 00562 endif 00563 00564 omegap = 2*par%px/par%Trep 00565 iomegap = 1/omegap 00566 ! Problem with derivation of DUdt in U and V points 00567 ! ueu and veu are not known until end of flow timestep, so 00568 ! both are based on previous time step. The temporal gradient 00569 ! can only be found by going back one more time step. 00570 ! We need to keep track of par%dt so the correct time step for 00571 ! the gradient is used: 00572 ! 00573 ! i-2 i-1 i (now) i+1 (future) 00574 ! dtold(1) dtold(2) par%dt 00575 ! ueuold ueuf unknown 00576 ! 00577 ! Note that this assumes that the acceleration term does not change significantly 00578 ! over 2 timesteps and that the location mask wetu can still be used, 00579 ! even though it is not the wetu mask of 2 timesteps ago. 00580 Tsmooth = par%Trep/20 00581 factime = min(dtold(1)/Tsmooth,1.d0) 00582 where(s%wetu==1) 00583 !ueuf = (1-factime)*ueuf + factime*s%ueu 00584 !veuf = (1-factime)*veuf + factime*s%veu 00585 !dudtsmooth = sqrt(((ueuf-ueuold)/dtold(1))**2 + & 00586 ! ((veuf-veuold)/dtold(1))**2) 00587 ! 00588 ! Robert: stick to 1D approach, similar to acceleration_boundary_layer_effect subroutine 00589 ueuf = (1-factime)*ueuf + factime*s%ueu 00590 dudtsmooth = (ueuf-ueuold)/dtold(1) 00591 elsewhere 00592 ueuf = 0.d0 00593 !veuf = 0.d0 00594 dudtsmooth = 0.d0 00595 endwhere 00596 00597 ! in some cases we don't need to compute v-acceleration 00598 ! - 1D cases without short wave forcing (e.g., non-h) has only u components 00599 ! - 1D surfbeat / stationary with short wave forcing has no v-components if boundary set to wall 00600 if (par%ny==0 .and. ((par%left==LR_WALL .and. par%right==LR_WALL) .or. par%swave==0)) then 00601 dvdtsmooth = 0.d0 00602 else 00603 where(s%wetv==1) 00604 vevf = (1-factime)*vevf + factime*s%vev 00605 dvdtsmooth = (vevf-vevold)/dtold(1) 00606 !uevf = (1-factime)*uevf + factime*uev 00607 !dvdtsmooth = sqrt(((vevf-vevold)/dtold(1))**2 + & 00608 ! ((uevf-uevold)/dtold(1))**2) 00609 elsewhere 00610 vevf = 0.d0 00611 !uevf = 0.d0 00612 dvdtsmooth = 0.d0 00613 endwhere 00614 endif 00615 ! now set the old variable to the i-1 variables 00616 ueuold = ueuf 00617 !uevold = uevf 00618 vevold = vevf 00619 !veuold = veuf 00620 dtold(1) = dtold(2) 00621 dtold(2) = par%dt 00622 ! New approach: 00623 ! we want to use Nielsen to compute bed shear stress (taubx_total): 00624 ! taubx_total = cf*rho*(cos(phi)*ue + 1/w*sin(phi)*dudt)^2*sign(u) 00625 ! 00626 ! in flow timestep the "regular" taubx is computed 00627 ! (note, here we again ignore urms component, Nielsen not meant for wave-averaged approach): 00628 ! taubx = cf*rho*ueu*umageu 00629 ! 00630 ! therefore: 00631 ! taubx_add = taubx_total-taubx = cf*rho*((cos(phi)*ue + 1/w*sin(phi)*dudt)^2*sign(u)-ueu*umageu) 00632 ! 00633 ! there is inconsistency in definitions of ue and umageu (v-component), which is not resolved in Nielsen paper, but if we 00634 ! rewrite to maintain consistency for stationary case (phi = 0 and dudt = 0) then: 00635 ! 00636 ! taubx_total = cf*rho* (cos^2(phi)*ueu*vmageu + 2*cos(phi)*ueu*sin(phi)*1/w*dudt + (sin(phi)*1/w*dudt)^2*sign(ueu)) 00637 ! and 00638 ! taubx_add = cf*rho * ((cos^2(phi)-1)*ueu*vmageu + 2*cos(phi)*ueu*sin(phi)*1/w*dudt + (sin(phi)*1/w*dudt)^2*sign(ueu)) 00639 ! (ignoring urms component in taubx) 00640 where(s%wetu==1) 00641 s%taubx_add = s%taubx_add + s%cfu*par%rho*( ((cos(phirad))**2-1)*s%ueu*s%vmageu & 00642 + 2*cos(phirad)*s%ueu*sin(phirad)*iomegap*dudtsmooth & 00643 + (sin(phirad)*iomegap*dudtsmooth)**2*sign(1.d0,s%ueu) & 00644 ) 00645 endwhere 00646 if (par%ny==0 .and. ((par%left==LR_WALL .and. par%right==LR_WALL) .or. par%swave==0)) then 00647 ! do nothing, leave tauby_add as is 00648 else 00649 where(s%wetv==1) 00650 s%tauby_add = s%tauby_add + s%cfv*par%rho*( ((cos(phirad))**2-1)*s%vev*s%vmagev & 00651 + 2*cos(phirad)*s%vev*sin(phirad)*iomegap*dvdtsmooth & 00652 + (sin(phirad)*iomegap*dvdtsmooth)**2*sign(1.d0,s%vev) & 00653 ) 00654 endwhere 00655 endif 00656 00657 00658 end subroutine acceleration_boundary_layer_effect_Nielsen 00659 00660 subroutine turbulence_boundary_layer_effect(s,par) 00661 00662 use params, only: parameters 00663 use spaceparamsdef, only: spacepars 00664 00665 IMPLICIT NONE 00666 00667 type(parameters),intent(in) :: par 00668 type(spacepars) :: s 00669 real*8,dimension(:,:),allocatable,save :: kbl 00670 integer :: i,j 00671 00672 ! Following Reniers et al (2004) 00673 ! ubed_turb = sqrt(ubed^2+gamma*kb) :: gamma ~ 1 00674 ! 00675 ! New way: compute taubx_add and tauby_add 00676 ! taubx_total = cf*rho*(ubed^2+gamma*kb) 00677 ! taubx = cf*rho*(ubed^2) 00678 ! taubx_add = cf*rho*(gamma*kb) 00679 00680 if (.not.allocated(kbl)) then 00681 allocate (kbl(s%nx+1,s%ny+1)) 00682 endif 00683 00684 ! U-points 00685 do j=1,s%ny+1 00686 do i=1,s%nx 00687 if(s%wetu(i,j)==1) then 00688 kbl(i,j) = 0.5d0*(s%kb(i,j)+s%kb(i+1,j)) 00689 else 00690 kbl(i,j) = 0.d0 00691 endif 00692 enddo 00693 enddo 00694 kbl(s%nx+1,:) = kbl(s%nx,:) 00695 where(s%wetu==1 .and. kbl>0.d0) 00696 ! New way: 00697 s%taubx_add = s%taubx_add + s%cfu*par%rho*par%gamma_turb*kbl*sign(1.d0,s%ue) 00698 00699 ! Old way: 00700 !cfu = cfu*(1.d0+(kbl/max(uu**2,1.d-6))) 00701 !cfu = cfu*(1.d0+(2.d0**(1.25d0)*sqrt(kbl))/max(abs(uu),1.d-6)+(sqrt(2.d0)*kbl)/max(uu**2,1.d-6)) 00702 endwhere 00703 ! V-points 00704 if(par%ny>0) then 00705 do j=1,s%ny 00706 do i=1,s%nx+1 00707 if(s%wetv(i,j)==1) then 00708 kbl(i,j) = 0.5d0*(s%kb(i,j)+s%kb(i,j+1)) 00709 else 00710 kbl(i,j) = 0.d0 00711 endif 00712 enddo 00713 enddo 00714 kbl(:,s%ny+1) = kbl(:,s%ny) 00715 else 00716 ! v-points centered on cell centre 00717 kbl = s%kb 00718 endif 00719 where(s%wetv==1 .and. kbl>0.d0) 00720 ! New way: 00721 s%tauby_add = s%tauby_add + s%cfv*par%rho*par%gamma_turb*kbl*sign(1.d0,s%ve) 00722 00723 ! Old way: 00724 !cfv = cfv*(1.d0+(kbl/max(vv**2,1.d-6))) 00725 !cfv = cfv*(1.d0+(2.d0**(1.25d0)*sqrt(kbl))/max(abs(vv),1.d-6)+(sqrt(2.d0)*kbl)/max(vv**2,1.d-6)) 00726 endwhere 00727 00728 end subroutine turbulence_boundary_layer_effect 00729 00730 00731 end module bedroughness_module
1.7.6.1
