|
XBeach
|
00001 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00002 ! Copyright (C) 2011 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 ! 00028 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00029 ! VEGETATION MODULE XBEACH: ATTENUATION OF SHORT WAVES, IG WAVES, FLOW, AND WAVE SETUP 00030 ! 00031 ! Version 1.0: 00032 ! Attenuation of short waves and IG waves 00033 ! Jaap van Thiel de Vries, okt 2013, 00034 ! see Linh K. Phan, Jaap S.M. van Thiel de Vries, and Marcel J.F. Stive (2015) Coastal Mangrove Squeeze in the Mekong Delta. Journal of Coastal Research: Volume 31, Issue 2: pp. 233 – 243. 00035 ! 00036 ! Version 2.0: 00037 ! Attenuation of short waves, IG waves and nonlinear wave effects 00038 ! Arnold van Rooijen, okt 2015, 00039 ! see Van Rooijen, McCall, Van Thiel de Vries, Van Dongeren, Reniers and Roelvink (2016), Modeling the effect of wave-vegetation interaction on wave setup, JGR Oceans 121, pp 4341-4359. 00040 ! 00041 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00042 00043 module vegetation_module 00044 use typesandkinds, only: slen 00045 use constants, only: pi 00046 implicit none 00047 save 00048 00049 private 00050 public veggie 00051 type veggie 00052 character(slen) :: name ! Name of vegetation specification file 00053 integer :: nsec ! Number of sections used in vertical schematization of vegetation [-] 00054 real*8 , dimension(:) , allocatable :: ah ! Height of vertical sections used in vegetation schematization [m wrt zb_ini (zb0)] 00055 real*8 , dimension(:) , allocatable :: Cd ! Bulk drag coefficient [-] 00056 real*8 , dimension(:) , allocatable :: bv ! Width/diameter of individual vegetation stems [m] 00057 integer , dimension(:) , allocatable :: N ! Number of vegetation stems per unit horizontal area [m-2] 00058 end type veggie 00059 00060 public veggie_init 00061 public vegatt 00062 public porcanflow ! porous in-canopy model 00063 00064 contains 00065 00066 subroutine veggie_init(s,par) 00067 use params, only: parameters 00068 use xmpi_module 00069 use spaceparamsdef, only: spacepars 00070 use readkey_module, only: readkey_dblvec, readkey_int, count_lines 00071 use filefunctions, only: create_new_fid, check_file_exist 00072 use logging_module, only: writelog, report_file_read_error 00073 use interp 00074 use paramsconst 00075 00076 IMPLICIT NONE 00077 00078 type(parameters) :: par 00079 type(spacepars), target :: s 00080 00081 !character(1) :: ch 00082 type(veggie), dimension(:), allocatable :: veg 00083 integer :: i,j,fid,ier,is,m,ind 00084 00085 if (par%vegetation == 1) then 00086 ! INITIALIZATION OF VEGETATION 00087 ! Read files with vegetation properties: 00088 ! file 1: list of species 00089 ! file 2: vegetation properties per specie (could be multiple files) 00090 ! file 3: distribution of species over space 00091 call writelog('l','','--------------------------------') 00092 call writelog('l','','Initializing vegetation input settings ') 00093 00094 ! 1) Read veggiefile with veggie species 00095 par%nveg = count_lines(par%veggiefile) 00096 00097 if (xmaster) then 00098 allocate(veg(par%nveg)) 00099 fid=create_new_fid() 00100 call check_file_exist(par%veggiefile) 00101 open(fid,file=par%veggiefile) 00102 do i=1,par%nveg 00103 read(fid,'(a)',iostat=ier) veg(i)%name 00104 enddo 00105 close(fid) 00106 00107 allocate(s%vegtype(par%nx+1, par%ny+1)) 00108 allocate(s%Dveg(par%nx+1, par%ny+1)) 00109 allocate(s%Fvegu(par%nx+1, par%ny+1)) 00110 allocate(s%Fvegv(par%nx+1, par%ny+1)) 00111 allocate(s%ucan(par%nx+1, par%ny+1)) 00112 allocate(s%vcan(par%nx+1, par%ny+1)) 00113 00114 s%vegtype = 0 00115 s%Dveg = 0.d0 00116 s%Fvegu = 0.d0 00117 s%Fvegv = 0.d0 00118 s%ucan = 0.d0 00119 s%vcan = 0.d0 00120 s%nsecvegmax = 1 00121 00122 ! 2) Read spatial distribution of all vegetation species 00123 ! NB: vegtype = 1 corresponds to first vegetation specified in veggiefile etc. 00124 fid=create_new_fid() ! see filefunctions.F90 00125 call check_file_exist(par%veggiemapfile) 00126 00127 select case(par%gridform) 00128 case(GRIDFORM_XBEACH) 00129 open(fid,file=par%veggiemapfile) 00130 do j=1,s%ny+1 00131 read(fid,*,iostat=ier)(s%vegtype(i,j),i=1,s%nx+1) 00132 if (ier .ne. 0) then 00133 call report_file_read_error(par%veggiemapfile) 00134 endif 00135 enddo 00136 close(fid) 00137 case (GRIDFORM_DELFT3D) 00138 open(fid,file=par%veggiemapfile,status='old') 00139 do j=1,s%ny+1 00140 read(fid,*,iostat=ier)(s%vegtype(i,j),i=1,s%nx+1) 00141 if (ier .ne. 0) then 00142 call report_file_read_error(par%veggiemapfile) 00143 endif 00144 enddo 00145 close(fid) 00146 end select 00147 00148 ! 3) Allocate and read vegetation properties for every species 00149 do is=1,par%nveg ! for each species 00150 call check_file_exist(veg(is)%name) 00151 veg(is)%nsec = readkey_int(veg(is)%name,'nsec', 1, 1, 100, silent=.true., bcast=.false.) 00152 ! Number of vertical sections in vegetation schematization (max = 100) 00153 00154 allocate (veg(is)%ah(veg(is)%nsec)) 00155 allocate (veg(is)%Cd(veg(is)%nsec)) 00156 allocate (veg(is)%bv(veg(is)%nsec)) 00157 allocate (veg(is)%N(veg(is)%nsec)) 00158 00159 veg(is)%ah = readkey_dblvec(veg(is)%name,'ah',veg(is)%nsec,size(veg(is)%ah), 0.1d0, 0.01d0, 2.d0, bcast=.false. ) 00160 veg(is)%bv = readkey_dblvec(veg(is)%name,'bv',veg(is)%nsec,size(veg(is)%bv), 0.01d0, 0.001d0, 1.0d0, bcast=.false. ) 00161 veg(is)%N = nint(readkey_dblvec(veg(is)%name,'N', veg(is)%nsec,size(veg(is)%N) ,100.0d0, 1.0d0, 5000.d0, bcast=.false. )) 00162 veg(is)%Cd = readkey_dblvec(veg(is)%name,'Cd',veg(is)%nsec,size(veg(is)%Cd), 0.0d0, 0.0d0, 3d0, bcast=.false. ) 00163 00164 ! Get maximum number of vegetation sections within model domain - needed to set size of Cd, ah, bv and Nv matrix 00165 s%nsecvegmax = max(s%nsecvegmax, veg(is)%nsec) 00166 enddo 00167 00168 ! Create spatially varying nsec, ah, bv, N and Cd within s-structure 00169 allocate(s%nsecveg(par%nx+1, par%ny+1)) 00170 allocate(s%Cdveg(par%nx+1, par%ny+1, s%nsecvegmax)) 00171 allocate(s%ahveg(par%nx+1, par%ny+1, s%nsecvegmax)) 00172 allocate(s%bveg(par%nx+1, par%ny+1, s%nsecvegmax)) 00173 allocate(s%Nveg(par%nx+1, par%ny+1, s%nsecvegmax)) 00174 00175 do j = 1,s%ny+1 00176 do i = 1,s%nx+1 00177 ind = s%vegtype(i,j) 00178 if (ind > 0) then ! set vegetation properties at locations where vegetation is present 00179 s%nsecveg(i,j) = veg(ind)%nsec 00180 do m=1,s%nsecveg(i,j) 00181 s%Cdveg(i,j,m) = veg(ind)%Cd(m) 00182 s%ahveg(i,j,m) = veg(ind)%ah(m) 00183 s%bveg(i,j,m) = veg(ind)%bv(m) 00184 s%Nveg(i,j,m) = veg(ind)%N(m) 00185 enddo 00186 else ! set to zero at locations of no vegetation 00187 s%nsecveg(i,j) = 0 00188 s%Cdveg(i,j,:) = 0.d0 00189 s%ahveg(i,j,:) = 0.d0 00190 s%bveg(i,j,:) = 0.d0 00191 s%Nveg(i,j,:) = 0.d0 00192 endif 00193 enddo 00194 enddo 00195 deallocate(veg) 00196 endif 00197 00198 call writelog('l','','--------------------------------') 00199 call writelog('l','','Finished reading vegetation input... ') 00200 00201 else ! par%vegetation == 0 00202 if (xmaster) then 00203 ! just allocate address for memory, only on xmaster, rest is 00204 ! done automatically by call from libxbeach 00205 allocate(s%vegtype(par%nx+1, par%ny+1)) 00206 allocate(s%nsecveg(par%nx+1, par%ny+1)) 00207 allocate(s%Cdveg(par%nx+1, par%ny+1, 1)) 00208 allocate(s%ahveg(par%nx+1, par%ny+1, 1)) 00209 allocate(s%bveg(par%nx+1, par%ny+1, 1)) 00210 allocate(s%Nveg(par%nx+1, par%ny+1, 1)) 00211 allocate(s%Dveg(par%nx+1, par%ny+1)) 00212 allocate(s%Fvegu(par%nx+1, par%ny+1)) 00213 allocate(s%Fvegv(par%nx+1, par%ny+1)) 00214 allocate(s%ucan(par%nx+1, par%ny+1)) 00215 allocate(s%vcan(par%nx+1, par%ny+1)) 00216 s%vegtype = 0 00217 s%nsecveg = 0 00218 s%nsecvegmax = 1 00219 s%Cdveg = 0.d0 00220 s%ahveg = 0.d0 00221 s%bveg = 0.d0 00222 s%Nveg = 0.d0 00223 s%Dveg = 0.d0 00224 s%Fvegu = 0.d0 00225 s%Fvegv = 0.d0 00226 s%ucan = 0.d0 00227 s%vcan = 0.d0 00228 endif 00229 endif 00230 end subroutine veggie_init 00231 00232 subroutine vegatt(s,par) 00233 use params, only: parameters 00234 use spaceparams 00235 use readkey_module 00236 use xmpi_module 00237 00238 type(parameters) :: par 00239 type(spacepars) :: s 00240 00241 ! local variables 00242 integer :: i,j,m 00243 real*8 :: Cdterm 00244 00245 ! Skip in case of using porous in-canopy model 00246 if (par%porcanflow == 1) then 00247 call porcanflow(s,par) 00248 return 00249 endif 00250 00251 ! First compute drag coefficient (if not user-defined) 00252 do j=1,s%ny+1 00253 do i=1,s%nx+1 00254 if (s%nsecveg(i,j) > 0) then ! only in case vegetation is present 00255 do m=1,s%nsecveg(i,j) ! for each vertical vegetation section 00256 if (s%Cdveg(i,j,m) < 0.d0) then ! If Cd is not user specified: call subroutine of M. Bendoni 00257 call bulkdragcoeff(s,par,s%ahveg(i,j,m)+s%zb0(i,j)-s%zb(i,j),m,i,j,Cdterm) 00258 s%Cdveg(i,j,m) = Cdterm 00259 endif 00260 enddo 00261 endif 00262 enddo 00263 enddo 00264 00265 ! Attenuation by vegetation is computed in wave action balance (swvegatt) and the momentum balance (momeqveg); 00266 ! 00267 ! 1) Short wave dissipation by vegetation 00268 call swvegatt(s,par) 00269 00270 ! 2) Mom.Eq.: Long wave dissipation, mean flow dissipation, nonlinear short wave effects, effect of emerged vegetation 00271 call momeqveg(s,par) 00272 00273 end subroutine vegatt 00274 00275 subroutine swvegatt(s,par) 00276 use params, only: parameters 00277 use spaceparams 00278 use readkey_module 00279 00280 00281 00282 type(parameters) :: par 00283 type(spacepars), target :: s 00284 !type(veggie), dimension(:), pointer :: veg 00285 00286 ! local variables 00287 integer :: i,j,m ! indices of actual x,y point 00288 real*8 :: aht,hterm,htermold,Dvgt,ahtold 00289 real*8, dimension(s%nx+1,s%ny+1) :: Dvg,kmr 00290 00291 !include 's.ind' 00292 !include 's.inp' 00293 00294 kmr = min(max(s%k, 0.01d0), 100.d0) 00295 00296 ! Set dissipation in vegetation to zero everywhere for a start 00297 Dvg = 0.d0 00298 do j=1,s%ny+1 00299 do i=1,s%nx+1 00300 htermold = 0.d0 00301 ahtold = 0.d0 00302 if (s%nsecveg(i,j)>0) then ! only if vegetation is present at (i,j) 00303 do m=1,s%nsecveg(i,j) 00304 00305 ! Determine height of vegetation section (restricted to current bed level) 00306 !aht = veg(ind)%ah(m)+ahtold !+s%zb0(i,j)-s%zb(i,j)!(max(veg(ind)%zv(m)+s%zb0(i,j),s%zb(i,j))) 00307 aht = s%ahveg(i,j,m)+ahtold 00308 00309 ! restrict vegetation height to local water depth 00310 aht = min(aht,s%hh(i,j)) 00311 00312 ! compute hterm based on ah 00313 hterm = (sinh(kmr(i,j)*aht)**3+3*sinh(kmr(i,j)*aht))/(3.d0*kmr(i,j)*cosh(kmr(i,j)*s%hh(i,j))**3) 00314 00315 ! compute dissipation based on aht and correct for lower elevated dissipation layers (following Suzuki et al. 2012) 00316 Dvgt = 0.5d0/sqrt(par%px)*par%rho*s%Cdveg(i,j,m)*s%bveg(i,j,m)*s%Nveg(i,j,m)*(0.5d0*kmr(i,j)*par%g/s%sigm(i,j))**3*(hterm-htermold)*s%H(i,j)**3 00317 00318 ! save hterm to htermold to correct possibly in next vegetation section 00319 htermold = hterm 00320 ahtold = aht 00321 00322 ! add dissipation current vegetation section 00323 Dvg(i,j) = Dvg(i,j) + Dvgt 00324 enddo 00325 endif 00326 enddo 00327 enddo 00328 s%Dveg = Dvg 00329 00330 end subroutine swvegatt 00331 00332 subroutine momeqveg(s,par) 00333 use params, only: parameters 00334 use spaceparams 00335 use readkey_module 00336 00337 use paramsconst 00338 00339 type(parameters) :: par 00340 type(spacepars) :: s 00341 !type(veggie), dimension(:), pointer :: veg 00342 00343 ! local variables 00344 integer :: i,j,m ! indices of actual x,y point 00345 real*8 :: aht,ahtold,Fvgtu,Fvgtv,FvgStu,FvgStv,watr,wacr,uabsu,vabsv 00346 real*8 :: Fvgnlt,Fvgnlu,Fvgnlv,FvgCan,FvgCav,FvgCau,ucan,uabsunl !uabsunl,vabsvnl,hterm,htermold, 00347 real*8, dimension(s%nx+1,s%ny+1) :: Fvgu,Fvgv,kmr 00348 real*8, save :: totT 00349 real*8, dimension(s%nx+1,s%ny+1,50) :: unl0,etaw0 00350 real*8, save, allocatable, dimension(:,:,:) :: unl,etaw 00351 real*8, dimension(50) :: hvegeff,Fvgnlu0 00352 real*8, dimension(:,:), allocatable,save :: sinthm, costhm 00353 00354 !include 's.ind' 00355 !include 's.inp' 00356 00357 ! Compute one force related to vegetation present in the water column: 00358 ! 1) Long wave velocity (ul) 00359 ! 2) Stokes velocity (us) 00360 ! 3) Non linear short wave velocity residual (ua) 00361 ! 4) return flow / undertow (ue) 00362 ! 5) wave-induced in-canopy flow (?) 00363 00364 ! only allocate in 1st timestep 00365 if (.not. allocated(sinthm)) then 00366 allocate (sinthm(s%nx+1,s%ny+1)) 00367 allocate (costhm(s%nx+1,s%ny+1)) 00368 endif 00369 kmr = min(max(s%k, 0.01d0), 100.d0) 00370 00371 Fvgu = 0.d0 00372 Fvgv = 0.d0 00373 Fvgnlt = 0.d0 00374 Fvgnlu = 0.d0 00375 Fvgnlv = 0.d0 00376 FvgStu = 0.d0 00377 FvgStv = 0.d0 00378 ucan = 0.d0 00379 FvgCan = 0.d0 00380 FvgCav = 0.d0 00381 FvgCau = 0.d0 00382 uabsunl = 0.d0 00383 00384 costhm = cos(s%thetamean-s%alfaz) 00385 sinthm = sin(s%thetamean-s%alfaz) 00386 00387 ! initialize totT 00388 if (par%dt == par%t) then 00389 totT = par%Trep 00390 endif 00391 if (par%vegnonlin == 1 .and. par%wavemodel/=WAVEMODEL_NONH) then 00392 ! only compute new nonlinear velocity profile every Trep s 00393 if(totT >= par%Trep) then 00394 call swvegnonlin(s,par,unl0,etaw0) 00395 unl = unl0 00396 etaw = etaw0 00397 totT = 0.d0 00398 else 00399 totT = totT + par%dt 00400 endif 00401 endif 00402 00403 do j=1,s%ny+1 00404 do i=1,s%nx+1 00405 ahtold = 0.d0 00406 if (s%nsecveg(i,j)>0) then ! Only if vegetation is present 00407 00408 ! Compute uabsu for calculation of Fveg 00409 uabsu = 0.d0 00410 vabsv = 0.d0 00411 Fvgnlu0 = 0.d0 00412 00413 watr = 0d0 00414 wacr = 0d0 00415 do m=1,s%nsecveg(i,j) 00416 ! Determine height of vegetation section (restricted to current bed level) 00417 aht = s%ahveg(i,j,m)+s%zb0(i,j)-s%zb(i,j) 00418 00419 ! Determine which part of the vegetation is below the wave trough, and between trough and crest 00420 if (par%vegnonlin == 1 .and. par%wavemodel/=WAVEMODEL_NONH) then 00421 watr = minval(etaw(i,j,:)) 00422 watr = s%hh(i,j) + watr ! wave trough level 00423 wacr = maxval(etaw(i,j,:)) 00424 wacr = s%hh(i,j) + wacr ! wave crest level 00425 else 00426 watr = s%hh(i,j) 00427 wacr = s%hh(i,j) 00428 endif 00429 00430 if (ahtold > wacr) then ! if plant section is entirely above wave crest, then do nothing 00431 00432 ! mean and long wave flow (ue) 00433 Fvgtu = 0d0 00434 Fvgtv = 0d0 00435 00436 ! nonlinear waves 00437 Fvgnlu = 0.d0 00438 Fvgnlv = 0.d0 00439 00440 else ! vegetation section is located (partly) in between wave trough and crest level 00441 if (par%veguntow == 1) then 00442 ! mean and long wave flow (ue, ve) 00443 Fvgtu = max((min(aht,watr)-ahtold),0d0)*0.5d0*s%Cdveg(i,j,m)*s%bveg(i,j,m)*s%Nveg(i,j,m)*(s%ueu(i,j)*s%vmageu(i,j)) 00444 Fvgtv = max((min(aht,watr)-ahtold),0d0)*0.5d0*s%Cdveg(i,j,m)*s%bveg(i,j,m)*s%Nveg(i,j,m)*(s%vev(i,j)*s%vmageu(i,j)) 00445 else 00446 ! Only long wave velocity (assume undertow is diverted over vegetation) 00447 Fvgtu = max((min(aht,watr)-ahtold),0d0)*0.5d0*s%Cdveg(i,j,m)*s%bveg(i,j,m)*s%Nveg(i,j,m)*(s%uu(i,j)*s%vmagu(i,j)) 00448 Fvgtv = max((min(aht,watr)-ahtold),0d0)*0.5d0*s%Cdveg(i,j,m)*s%bveg(i,j,m)*s%Nveg(i,j,m)*(s%vv(i,j)*s%vmagu(i,j)) 00449 endif 00450 00451 ! nonlinear waves (including emerged vegetation effect) 00452 !etaw = 0.d0 00453 if (par%vegnonlin == 1 .and. par%wavemodel/=WAVEMODEL_NONH) then 00454 hvegeff = max(etaw(i,j,:) + s%hh(i,j)-ahtold,0.d0) ! effective vegetation height over a wave cycle 00455 Fvgnlt = trapz(((0.5d0*s%Cdveg(i,j,m)*s%bveg(i,j,m)*s%Nveg(i,j,m))*min(hvegeff,aht)*unl(i,j,:)*abs(unl(i,j,:))),par%Trep/50)/s%hh(i,j) 00456 00457 ! decompose in u and v-direction 00458 Fvgnlu = Fvgnlt*costhm(i,j) 00459 Fvgnlv = Fvgnlt*sinthm(i,j) 00460 endif 00461 00462 ! wave induced incanopy flow (Luhar et al., 2010) 00463 ucan = sqrt(4.d0*kmr(i,j)*par%Trep*s%urms(i,j)**3/(6.d0*par%px**2)) 00464 FvgCan = max((min(aht,watr)-ahtold),0d0)/s%hh(i,j)*0.5d0*s%Cdveg(i,j,m)*s%bveg(i,j,m)*s%Nveg(i,j,m)*ucan**2 00465 00466 ! decompose in u and v-direction 00467 FvgCau = FvgCan*costhm(i,j) 00468 FvgCav = FvgCan*sinthm(i,j) 00469 endif 00470 00471 ! save aht to ahtold to correct possibly in next vegetation section 00472 ahtold = aht 00473 00474 ! add Forcing current layer 00475 Fvgu(i,j) = Fvgu(i,j) + Fvgtu 00476 Fvgv(i,j) = Fvgv(i,j) + Fvgtv 00477 00478 if (par%vegnonlin == 1 .and. par%wavemodel/=WAVEMODEL_NONH) then ! add nonlin wave effect 00479 Fvgu(i,j) = Fvgu(i,j) + Fvgnlu 00480 Fvgv(i,j) = Fvgv(i,j) + Fvgnlv 00481 endif 00482 if (par%vegcanflo == 1) then ! add in canopy flow (Luhar et al., 2010) 00483 Fvgu(i,j) = Fvgu(i,j) + FvgCau 00484 Fvgv(i,j) = Fvgv(i,j) + FvgCav 00485 endif 00486 enddo 00487 endif 00488 enddo 00489 enddo 00490 00491 s%Fvegu = Fvgu*par%rho ! make sure units of drag force are consistent (N/m2) 00492 s%Fvegv = Fvgv*par%rho ! make sure units of drag force are consistent (N/m2) 00493 00494 end subroutine momeqveg 00495 00496 subroutine swvegnonlin(s,par,unl0,etaw0) 00497 use params, only: parameters 00498 use spaceparams 00499 00500 IMPLICIT NONE 00501 00502 type(parameters) :: par 00503 type(spacepars) :: s 00504 00505 integer :: i,j 00506 integer :: irf,ih0,it0,jrf,ih1,it1 !,m,ind,ih0,it0,ih1,it1,irf,jrf ! indices of actual x,y point 00507 integer, save :: nh,nt 00508 real*8 :: p,q,f0,f1,f2,f3 !,uabsunl,vabsvnl 00509 real*8, save :: dh,dt 00510 real*8, dimension(s%nx+1,s%ny+1) :: kmr,Urs,phi,w1,w2 00511 real*8, dimension(8),save :: urf0 00512 real*8, dimension(50),save :: urf2,urf !,urfueurfu 00513 real*8, dimension(50,8),save :: cs,sn,urf1 00514 real*8, dimension(:,:),save,allocatable :: h0,t0 00515 real*8, dimension(s%nx+1,s%ny+1,50),intent(out) :: unl0,etaw0 00516 00517 ! Subroutine to compute a net drag force due to wave skewness. Based on (matlab based) roller model with veggies by Ad. 00518 ! 00519 ! Background: 00520 ! The drag force (Fveg) is a function of u*abs(u), which is zero for linear waves. For non-linear, skewed waves the 00521 ! depth-averaged velocity integrated over the wave period is zero. However, due to the sharp peaks and flat troughs 00522 ! the integral of u*abs(u) is non-zero, and can significantly reduce wave setup, or even lead to set-down (e.g. Dean & Bender,2006). 00523 ! 00524 ! Here we use a method based on Rienecker & Fenton (1981), similar to the method used for onshore sediment transport due to wave asymmetry/ 00525 ! skewness (see also morphevolution.F90 + Van Thiel de Vries Phd thesis par 6.2.3). 00526 ! 00527 00528 ! load Ad's RF-table (update for depth averaged velocities?) 00529 include 'RFveg.inc' 00530 00531 ! Initialize/Prepare for interpolation of RF-value from RFveg-table 00532 if (.not. allocated(h0)) then 00533 allocate (h0(s%nx+1,s%ny+1)) 00534 allocate (t0(s%nx+1,s%ny+1)) 00535 00536 dh = 0.03d0 00537 dt = 1.25d0 00538 nh = floor(0.54d0/dh); 00539 nt = floor(25.d0/dt); 00540 !construct velocity profile based on cosine/sine functions / Fourier components 00541 do irf=1,8 00542 do jrf=1,50 00543 cs(jrf,irf) = cos((jrf*2*par%px/50)*irf) 00544 sn(jrf,irf) = sin((jrf*2*par%px/50)*irf) 00545 enddo 00546 enddo 00547 endif 00548 00549 h0 = min(nh*dh,max(dh,min(s%H,s%hh)/s%hh)) 00550 t0 = min(nt*dt,max(dt,par%Trep*sqrt(par%g/s%hh))) 00551 00552 ! Initialize 00553 urf0 = 0.d0 00554 urf1 = 0.d0 00555 urf2 = 0.d0 00556 urf = 0.d0 00557 w1 = 0.d0 00558 w2 = 0.d0 00559 phi = 0.d0 00560 Urs = 0.d0 00561 kmr = 0.d0 00562 00563 ! ! Now compute weight factors (w1,w2) for relative contribution of cosine and sine functions (for w1 = 1: only cosines -> 00564 ! fully skewed Stokes wave, for w2 = 1: only sines -> fully asymmetric wave) based on Ruessink. 00565 kmr = min(max(s%k, 0.01d0), 100.d0) 00566 Urs = s%H/kmr/kmr/(s%hh**3)! Ursell number 00567 00568 ! Compute phase and weight factors 00569 phi = par%px/2*(1-tanh(0.815/(Urs**0.672)))! according to Ruessink et al 2012 (eq 10): p5 = 0.815 ipv 0.64; ip6 = 0.672 ipv 0.6, Dano&Ad book: 0.64 and 0.6 00570 w1 = 1-phi/(par%px/2)!w1 = 1.d0 if fully skewed waves 00571 w2 = 1.d0-w1 00572 ! or use relation between w1 and phi as in Phd thesis Jaap (eq 6.13)?? 00573 00574 ! Interpolate RieneckerFenton velocity from RFveg table from Ad 00575 ! in ftab-dimension, only read 4:11 and sum later 00576 do j=1,s%ny+1 00577 do i=1,s%nx+1 00578 ! interpolate RF table values.... 00579 ih0=floor(h0(i,j)/dh) 00580 it0=floor(t0(i,j)/dt) 00581 ih1=min(ih0+1,nh) 00582 it1=min(it0+1,nt) 00583 p=(h0(i,j)-ih0*dh)/dh 00584 q=(t0(i,j)-it0*dt)/dt 00585 f0=(1-p)*(1-q) 00586 f1=p*(1-q) 00587 f2=q*(1-p) 00588 f3=p*q 00589 00590 ! Compute velocity amplitude per component 00591 do irf=1,8 00592 urf0(irf) = f0*RFveg(irf+3,ih0,it0)+f1*RFveg(irf+3,ih1,it0)+ f2*RFveg(irf+3,ih0,it1)+f3*RFveg(irf+3,ih1,it1) 00593 enddo 00594 00595 ! fill velocity amplitude matrix urf1([50 time points, 8 components]) 00596 do irf=1,8 00597 urf1(:,irf) = urf0(irf) 00598 enddo 00599 00600 ! Compute velocity profile matrix per component 00601 urf1 = urf1*(w1(i,j)*cs+w2(i,j)*sn) 00602 00603 ! Add velocity components 00604 urf2 = sum(urf1,2) 00605 00606 ! Scale the results to get velocity profile over wave period 00607 unl0(i,j,:) = urf2*sqrt(par%g*s%hh(i,j)) 00608 etaw0(i,j,:) = unl0(i,j,:)*sqrt(max(s%hh(i,j),0.d0)/par%g) 00609 enddo 00610 enddo 00611 00612 end subroutine swvegnonlin 00613 00614 function trapz(y,dx) result (value) 00615 implicit none 00616 real*8 :: integral,value,dx 00617 real*8, dimension(:) :: y 00618 integer :: i,n 00619 00620 integral = 0.d0 00621 n = size(y)-1.d0 00622 do i=1,n 00623 integral = integral+dx*(y(i+1)+y(i))/2 00624 end do 00625 value = integral 00626 00627 end function trapz 00628 00629 subroutine bulkdragcoeff(s,par,ahh,m,i,j,Cdterm) 00630 ! Michele Bendoni: subroutine to calculate bulk drag coefficient for short wave 00631 ! energy dissipation based on the Keulegan-Carpenter number 00632 ! Ozeren et al. (2013) or Mendez and Losada (2004) 00633 00634 use params 00635 use spaceparams 00636 00637 implicit none 00638 00639 !type(veggie), dimension(:), pointer :: veg 00640 00641 type(parameters) :: par 00642 type(spacepars) :: s 00643 real*8, intent(out) :: Cdterm 00644 real*8, intent(in) :: ahh ! [m] plant (total) height 00645 integer, intent(in) :: m,i,j 00646 00647 ! Local variables 00648 real*8 :: alfav ! [-] ratio between plant height and water depth 00649 real*8 :: um ! [m/s] typical velocity acting on the plant 00650 real*8 :: Tp ! [s] reference wave period 00651 real*8 :: KC ! [-] Keulegan-Carpenter number 00652 real*8 :: Q ! [-] modified Keulegan-Carpenter number 00653 integer :: myflag ! 1 => Ozeren et al. (2013); 2 => Mendez and Losada (2004) 00654 ! 00655 ! 00656 myflag = 2 00657 ! 00658 ! Representative wave period 00659 Tp = 2*par%px/s%sigm(i,j) 00660 ! 00661 ! Coefficient alfa 00662 if (ahh>=s%hh(i,j)) then 00663 alfav = 1.d0 00664 else 00665 alfav = ahh/s%hh(i,j) 00666 end if 00667 ! 00668 ! Representative orbital velocity 00669 ! (Could we also use urms here?) 00670 um = 0.5d0*s%H(i,j)*s%sigm(i,j)*cosh(s%k(i,j)*alfav*s%hh(i,j))/sinh(s%k(i,j)*s%hh(i,j)) 00671 ! 00672 ! Keulegan-Carpenter number 00673 KC = um*Tp/s%bveg(i,j,m) 00674 if (KC > 0d0) then 00675 KC = KC 00676 endif 00677 ! 00678 ! Bulk drag coefficient 00679 if (myflag == 1) then 00680 ! 00681 ! Approach from Ozeren et al. (2013), eq? 00682 ! 00683 if (KC>=10.d0) then 00684 Cdterm = 0.036d0+50.d0/(KC**0.926d0) 00685 else 00686 Cdterm = 0.036d0+50.d0/(10.d0**0.926d0) 00687 endif 00688 elseif (myflag == 2) then 00689 ! 00690 ! Approach from Mendez and Losada (2004), eq. 40 00691 ! Only applicable for Laminaria Hyperborea (kelp)??? 00692 ! 00693 Q = KC/(alfav**0.76d0) 00694 if (Q>=7) then 00695 Cdterm = exp(-0.0138*Q)/(Q**0.3d0) 00696 else 00697 Cdterm = exp(-0.0138*7)/(7**0.3d0) 00698 endif 00699 endif 00700 ! 00701 end subroutine bulkdragcoeff 00702 00703 subroutine porcanflow(s,par) 00704 ! porous in-canopy model. Computes the in-canopy flow and vegetation force. 00705 use params 00706 use spaceparams 00707 use readkey_module 00708 use xmpi_module 00709 use filefunctions 00710 use interp 00711 use logging_module 00712 use spaceparamsdef, only: spacepars 00713 00714 implicit none 00715 00716 type(parameters) :: par 00717 type(spacepars) :: s 00718 !type(veggie), dimension(:), pointer :: veg 00719 00720 ! local variables 00721 integer :: i,j,imax,j1,switch_drag 00722 real*8 :: p,mu,lamp,Kp,beta,hcan,Cf,Cm,A,Fcanu,Fcanv,U,V,ucan_old,vcan_old !rhs, 00723 00724 ! Initialization paramters 00725 mu = 10.d0**(-6) ! kinematic viscosity 00726 Kp = par%Kp ! permeability 00727 Cm = par%Cm ! inertia coefficient 00728 U = 0.d0 00729 V = 0.d0 00730 ucan_old=0.d0 00731 vcan_old=0.d0 00732 switch_drag=1 00733 lamp = 0 ! compiler is not sure that lamp always gets a value 00734 hcan = 0 ! idem 00735 beta = 0 ! idem 00736 00737 ! Superfast 1D 00738 if (s%ny==0) then 00739 j1 = 1 00740 else 00741 j1 = 2 00742 endif 00743 00744 ! In canopy momentum balance 00745 do j=j1,max(s%ny,1) 00746 imax = s%nx 00747 do i=2,imax 00748 ! Only compute ucan if vegeation is pressent 00749 if(s%vegtype(i,j)>0.d0) then 00750 ! vegetation type 00751 p = s%Nveg(i,j,1)/100.d0 ! porosity 00752 lamp = (1-p) ! lambda parameters (Britter and Hanna, 2003) 00753 hcan = s%ahveg(i,j,1) ! canopy height 00754 beta = s%Cdveg(i,j,1) ! Drag 00755 Cf = s%bveg(i,j,1) ! Friction 00756 00757 !Emergent case. hcan> h 00758 if (s%hu(i,j) < hcan) then 00759 hcan = s%hu(i,j) 00760 endif 00761 00762 !Implicit term momentum equation 00763 A = (1+Cm*lamp/(1-lamp))/par%dt 00764 00765 ! Select free stream velocity for top shear stress. 00766 ! if (par%nonhq3d == 1 .and. par%nhlay > 0.0d0 .and. par%switch_2dv==1) then 00767 ! ! hcan < 0.5 alpha * h 00768 ! if (hcan < 0.5d0 * par%nhlay*s%hh(i,j)) then 00769 ! !Free stream velocity is velocity layer 1 00770 ! U = s%u(i,j) + (1.d0 - par%nhlay) * s%du(i,j) 00771 ! V = s%v(i,j) + (1.d0 - par%nhlay) * s%dv(i,j) 00772 ! ! hcan > 0.5 alpha * h + alpha * h 00773 ! elseif (hcan > 0.5d0 * par%nhlay * s%hh(i,j) + par%nhlay*s%hh(i,j)) then 00774 ! !Free stream velocity is velocity layer 2 00775 ! U = s%u(i,j) - par%nhlay * s%du(i,j) 00776 ! V = s%v(i,j) - par%nhlay * s%dv(i,j) 00777 ! ! 0.5 alpha * h > hcan > 0.5 alpha * h + alpha * h 00778 ! else 00779 ! ! Velocity layer 1 00780 ! u11 = s%u(i,j) + (1.d0 - par%nhlay) * s%du(i,j) 00781 ! ! Velocity layer 2 00782 ! u22 = s%u(i,j) - par%nhlay * s%du(i,j) 00783 ! ! Interpolate between layers for smooth transition. 00784 ! U = (u22-u11)/(0.5d0*(1.d0 - par%nhlay) *s%hh(i,j) + 0.5d0 * par%nhlay * s%hh(i,j) ) * (hcan - 1.d0/2.d0 * par%nhlay * s%hh(i,j)) 00785 ! 00786 ! ! Velocity layer 1 00787 ! v11 = s%v(i,j) + (1.d0 - par%nhlay) * s%dv(i,j) 00788 ! ! Velocity layer 2 00789 ! v22 = s%v(i,j) - par%nhlay * s%dv(i,j) 00790 ! ! Interpolate between layers for smooth transition. 00791 ! V = (v22-v11)/(0.5d0*(1.d0 - par%nhlay) *s%hh(i,j) + 0.5d0 * par%nhlay * s%hh(i,j) ) * (hcan - 0.5d0 * par%nhlay * s%hh(i,j)) 00792 ! endif 00793 ! else 00794 ! ! Depth averaged velocity 00795 ! U = s%u(i,j) 00796 ! V = s%v(i,j) 00797 ! endif 00798 00799 ! free stream velocty 00800 U = s%u(i,j) 00801 V = s%v(i,j) 00802 00803 ! ucan previous time step 00804 ucan_old = s%ucan(i,j) 00805 vcan_old = s%vcan(i,j) 00806 00807 ! Compute u-incanopy velocity when cell is wet 00808 if (s%wetu(i,j)>0) then 00809 s%ucan(i,j) = (-1 * par%g*s%dzsdx(i,j) + 0.5d0*Cf/hcan*abs(U-s%ucan(i,j)) * (U-s%ucan(i,j)) + A * s%ucan(i,j))/(A + mu * (1-lamp)/Kp + beta * abs(s%ucan(i,j))) 00810 ! prevent high in-canopy for flooding 00811 !if ((s%ucan(i,j)-ucan_old)/par%dt > 0.2) then 00812 ! s%ucan(i,j) = ucan_old 00813 !endif 00814 ! Zero velocity if dry 00815 else 00816 s%ucan(i,j) = 0.d0 00817 endif 00818 00819 ! Compute v-incanopy velocity when cell is wet 00820 if (s%wetv(i,j)>0 .and. s%ny>1) then 00821 s%vcan(i,j) = (-1 * par%g*s%dzsdx(i,j) + 0.5d0*Cf/hcan*abs(V-s%vcan(i,j)) * (V-s%vcan(i,j)) + A * s%vcan(i,j))/(A + mu * (1-lamp)/Kp + beta * abs(s%vcan(i,j))) 00822 ! prevent high in-canopy for flooding 00823 !if (abs(s%vcan(i,j)) > abs(10 * vcan_old)) then 00824 ! s%vcan(i,j) = vcan_old 00825 !endif 00826 ! Zero velocity if dry 00827 else 00828 s%vcan(i,j) = 0.d0 00829 endif 00830 00831 ! old shear stress formulation. u|u| instead of (u-uc)|u-uc| 00832 !s%ucan(i,j) = (-1 * par%g*s%dzsdx(i,j) + abs(U)*U/(2.d0*hcan/Cf) + A * s%ucan(i,j))/(A + mu * (1-lamp)/Kp + beta * abs(s%ucan(i,j))) 00833 00834 !Zero velocity if no vegetation 00835 else 00836 s%ucan(i,j) = 0.0d0 00837 s%vcan(i,j) = 0.0d0 00838 endif 00839 00840 !Upper limit canopy flow??? 00841 ! not used, because there can be a phase difference between uc and u. 00842 !if ( abs(s%ucan(i,j))>abs(U)) then 00843 ! s%ucan(i,j) = U 00844 !endif 00845 00846 ! Compute canopy drag force 00847 if(s%nsecveg(i,j)>0.d0 .and. switch_drag==1) then 00848 ! Compute vegetation force 00849 if (s%wetu(i,j)>0) then 00850 Fcanu = abs(s%ucan(i,j))*ucan_old*beta + mu*(1-lamp)/Kp*ucan_old + Cm*lamp/(1-lamp) * (s%ucan(i,j)-ucan_old)/par%dt 00851 else 00852 Fcanu = 0.d0 00853 endif 00854 00855 00856 if (s%wetv(i,j)>0 .and. s%ny>1) then 00857 Fcanv = abs(s%vcan(i,j))*vcan_old*beta + mu*(1-lamp)/Kp*vcan_old + Cm*lamp/(1-lamp) * (s%vcan(i,j)-vcan_old)/par%dt 00858 else 00859 Fcanv = 0.d0 00860 endif 00861 00862 !dfu for XBeach-nh+. 00863 !if (hcan > par%nhlay*s%hh(i,j) .and. par%switch_2dv==1) then 00864 ! s%dFvegu(i,j) = Fcanu * par%nhlay*s%hh(i,j)* par%rho - Fcanu * (hcan - par%nhlay*s%hh(i,j))* par%rho 00865 !else if (hcan < par%nhlay*s%hh(i,j) .and. par%switch_2dv==1) then 00866 ! s%dFvegu(i,j) = Fcanu * hcan * par%rho 00867 !else 00868 ! s%dFvegu(i,j) = 0 00869 !endif 00870 00871 ! Force times height and rho (divide by rho in momentum eq). 00872 Fcanu = Fcanu * par%rho * hcan 00873 Fcanv = Fcanv * par%rho * hcan 00874 else 00875 Fcanu = 0. 00876 Fcanv = 0. 00877 endif 00878 s%Fvegu(i,j) = Fcanu 00879 s%Fvegv(i,j) = Fcanv 00880 end do 00881 end do 00882 00883 end subroutine porcanflow 00884 00885 end module vegetation_module
1.7.6.1
