|
XBeach
|
00001 module morphevolution 00002 implicit none 00003 save 00004 00005 #ifdef HAVE_CONFIG_H 00006 #include "config.h" 00007 #endif 00008 contains 00009 subroutine transus(s,par) 00010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00011 ! Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University ! 00012 ! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski, ! 00013 ! Jaap van Thiel de Vries, Robert McCall ! 00014 ! ! 00015 ! d.roelvink@unesco-ihe.org ! 00016 ! UNESCO-IHE Institute for Water Education ! 00017 ! P.O. Box 3015 ! 00018 ! 2601 DA Delft ! 00019 ! The Netherlands ! 00020 ! ! 00021 ! This library is free software; you can redistribute it and/or ! 00022 ! modify it under the terms of the GNU Lesser General Public ! 00023 ! License as published by the Free Software Foundation; either ! 00024 ! version 2.1 of the License, or (at your option) any later version. ! 00025 ! ! 00026 ! This library is distributed in the hope that it will be useful, ! 00027 ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! 00028 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! 00029 ! Lesser General Public License for more details. ! 00030 ! ! 00031 ! You should have received a copy of the GNU Lesser General Public ! 00032 ! License along with this library; if not, write to the Free Software ! 00033 ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 ! 00034 ! USA ! 00035 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00036 use params 00037 use spaceparams 00038 use xmpi_module 00039 use interp 00040 use paramsconst 00041 ! use vsmumod 00042 00043 implicit none 00044 00045 type(spacepars),target :: s 00046 type(parameters) :: par 00047 00048 integer :: i,isig 00049 integer :: j,jg 00050 real*8 :: exp_ero 00051 00052 real*8,dimension(:),allocatable,save :: chain,cumchain 00053 real*8,dimension(:,:),allocatable,save :: vmag2,uau,uav,um,vm,ueu_sed,uev_sed,veu_sed,vev_sed 00054 real*8,dimension(:,:),allocatable,save :: ccvt,dcdz,dsigt,aref 00055 real*8,dimension(:,:),allocatable,save :: cc,ccb,cu,cv,Sus,Svs 00056 real*8,dimension(:,:),allocatable,save :: cub,cvb,Sub,Svb,pbbedu,pbbedv 00057 real*8,dimension(:,:),allocatable,save :: suq3d,svq3d,eswmax,eswbed,sigs,deltas 00058 real*8,dimension(:,:,:),allocatable,save :: dsig,ccv,sdif,cuq3d,cvq3d,fac 00059 logical,dimension(:,:),allocatable,save :: bermslopeindexbed,bermslopeindexsus,bermslopeindex 00060 00061 real*8,dimension(:,:),allocatable,save :: sinthm,costhm 00062 00063 real*8 :: delta,delta_x,shields,ftheta,psi_x,Sbtot ! Lodewijk: for direction of sediment transport (bed slope effect) 00064 real*8 :: Ssmtot, dzbds, Sbmtot ! Lodewijk: for magnitude of sediment transport (bed slope effect) 00065 00066 !include 's.ind' 00067 !include 's.inp' 00068 00069 if (.not. allocated(vmag2)) then 00070 allocate(vmag2 (s%nx+1,s%ny+1)) 00071 allocate(uau (s%nx+1,s%ny+1)) 00072 allocate(uav (s%nx+1,s%ny+1)) 00073 allocate(ueu_sed (s%nx+1,s%ny+1)) 00074 allocate(uev_sed (s%nx+1,s%ny+1)) 00075 allocate(veu_sed (s%nx+1,s%ny+1)) 00076 allocate(vev_sed (s%nx+1,s%ny+1)) 00077 allocate(cu (s%nx+1,s%ny+1)) 00078 allocate(cv (s%nx+1,s%ny+1)) 00079 allocate(cc (s%nx+1,s%ny+1)) 00080 allocate(ccb (s%nx+1,s%ny+1)) 00081 allocate(fac (s%nx+1,s%ny+1,par%ngd)) 00082 allocate(Sus (s%nx+1,s%ny+1)) 00083 allocate(Svs (s%nx+1,s%ny+1)) 00084 allocate(cub (s%nx+1,s%ny+1)) 00085 allocate(cvb (s%nx+1,s%ny+1)) 00086 allocate(Sub (s%nx+1,s%ny+1)) 00087 allocate(Svb (s%nx+1,s%ny+1)) 00088 allocate(pbbedu (s%nx+1,s%ny+1)) 00089 allocate(pbbedv (s%nx+1,s%ny+1)) 00090 allocate(ccvt (s%nx+1,s%ny+1)) 00091 allocate(dcdz (s%nx+1,s%ny+1)) 00092 allocate(dsigt (s%nx+1,s%ny+1)) 00093 allocate(dsig(s%nx+1,s%ny+1,par%kmax)) 00094 allocate(ccv(s%nx+1,s%ny+1,par%kmax)) 00095 allocate(sdif(s%nx+1,s%ny+1,par%kmax)) 00096 allocate(um (s%nx+1,s%ny+1)) 00097 allocate(vm (s%nx+1,s%ny+1)) 00098 allocate(deltas(s%nx+1,s%ny+1)) 00099 allocate(sigs(s%nx+1,s%ny+1)) 00100 allocate(eswmax(s%nx+1,s%ny+1)) 00101 allocate(eswbed(s%nx+1,s%ny+1)) 00102 allocate(suq3d(s%nx+1,s%ny+1)) 00103 allocate(svq3d(s%nx+1,s%ny+1)) 00104 allocate(cuq3d(s%nx+1,s%ny+1,par%kmax)) 00105 allocate(cvq3d(s%nx+1,s%ny+1,par%kmax)) 00106 allocate(aref(s%nx+1,s%ny+1)) 00107 allocate(chain(par%kmax)) 00108 allocate(cumchain(par%kmax)) 00109 allocate (sinthm(s%nx+1,s%ny+1)) 00110 allocate (costhm(s%nx+1,s%ny+1)) 00111 allocate(bermslopeindexbed(s%nx+1,s%ny+1)) 00112 allocate(bermslopeindexsus(s%nx+1,s%ny+1)) 00113 allocate(bermslopeindex(s%nx+1,s%ny+1)) 00114 00115 delta_x = 0.d0 ! Lodewijk 00116 shields = 0.d0 ! Lodewijk 00117 ftheta = 0.d0 ! Lodewijk 00118 psi_x = 0.d0 ! Lodewijk 00119 Sbtot = 0.d0 ! Lodewijk 00120 delta = 0.d0 ! Lodewijk 00121 uau = 0.d0 00122 uav = 0.d0 00123 um = 0.d0 00124 vm = 0.d0 00125 chain = 0.0d0 00126 cumchain = 0.0d0 00127 fac = 1.d0 00128 exp_ero = 0.d0 00129 ! generate sigma grid shape... 00130 do isig=2,par%kmax 00131 chain(isig) = par%sigfac**(isig-1) 00132 cumchain(isig) = cumchain(isig-1)+chain(isig) 00133 enddo 00134 bermslopeindex = .false. ! turned off unless needed 00135 bermslopeindexbed = .false. 00136 bermslopeindexsus = .false. 00137 endif 00138 00139 ! use eulerian velocities 00140 vmag2 = s%ue**2+s%ve**2 00141 cu = 0.0d0 00142 cv = 0.0d0 00143 cub = 0.0d0 00144 cvb = 0.0d0 00145 s%dcsdx = 0.0d0 00146 s%dcsdy = 0.0d0 00147 00148 sinthm = sin(s%thetamean-s%alfaz) 00149 costhm = cos(s%thetamean-s%alfaz) 00150 00151 ! short wave runup 00152 if (par%swrunup==1 .and. par%struct==1) then 00153 call hybrid(s,par) 00154 endif 00155 00156 ! compute turbulence due to wave breaking 00157 if (par%lwt==1 .or. par%turb == TURB_BORE_AVERAGED .or. par%turb == TURB_WAVE_AVERAGED) then 00158 call waveturb(s,par) 00159 endif 00160 00161 if (par%swave==1) then 00162 ! include wave skewness and assymetry in sediment advection velocity 00163 if (par%waveform==WAVEFORM_RUESSINK_VANRIJN)then 00164 call RvR(s,par) 00165 elseif (par%waveform==WAVEFORM_VANTHIEL) then 00166 call vT(s,par) 00167 endif 00168 00169 endif 00170 00171 ! calculate equilibrium concentration/sediment source 00172 select case (par%form) 00173 case (FORM_SOULSBY_VANRIJN,FORM_VANTHIEL_VANRIJN,FORM_VANRIJN1993) 00174 ! Soulsby van Rijn and Van Thiel de Vries & Reniers 2008 formulations 00175 call sedtransform(s,par) 00176 case (FORM_NIELSEN2006) 00177 call Nielsen2006(s,par) 00178 if (par%bulk==0) then 00179 return 00180 endif 00181 case (FORM_MCCALL_VANRIJN) 00182 call mccall_vanrijn(s,par) 00183 if (par%bulk==0) then 00184 return 00185 endif 00186 end select 00187 00188 ! compute reduction factor for sediment sources due to presence of hard layers 00189 do jg = 1,par%ngd 00190 do j=1,s%ny+1 00191 do i=1,s%nx+1 00192 exp_ero = par%morfac*par%dt/(1.d0-par%por)*s%hh(i,j)*(s%ceqsg(i,j,jg)*s%pbbed(i,j,1,jg)/s%Tsg(i,j,jg) & 00193 + s%ceqbg(i,j,jg)*s%pbbed(i,j,1,jg)/par%dt) 00194 ! limit erosion to available sediment on top of hard layer wwvv changed tiny into epsilon 00195 fac(i,j,jg) = min(1.d0,s%structdepth(i,j)*s%pbbed(i,j,1,jg)/max(epsilon(0.d0),exp_ero) ) 00196 !if (fac(i,j,jg)*exp_ero > dzbed(i,j,1)*pbbed(i,j,1,jg)) then 00197 ! limit erosion to available sand in top layer 00198 ! fac(i,j,jg) = min(fac(i,j,jg),dzbed(i,j,1)*pbbed(i,j,1,jg)/max(tiny(0.d0),exp_ero) ) 00199 ! write(*,*)'WARNING: expected erosion from top layer is larger than available sand in top layer' 00200 !endif 00201 enddo 00202 enddo 00203 enddo 00204 00205 ! compute diffusion coefficient 00206 s%Dc = par%facDc*s%nuh 00207 !Robert: removed this which is not grid-independent s%Dc = (par%nuh+par%nuhfac*s%hh*(s%DR/par%rho)**(1.d0/3.d0)) 00208 00209 do jg = 1,par%ngd 00210 cc = s%ccg(:,:,jg) 00211 if (s%D50(jg)>0.002d0) then 00212 ! RJ: set ceqsg to zero for gravel. 00213 ! Dano: try without this fix cc = 0.d0 ! Can be used to test total transport mode 00214 endif 00215 ! 00216 ! X-direction 00217 ! 00218 ! Get ua in u points and split out in u and v direction 00219 uau(1:s%nx,:) = 0.5*(s%ua(1:s%nx,:)*costhm(1:s%nx,:)+s%ua(2:s%nx+1,:)*costhm(1:s%nx,:)) 00220 uav(1:s%nx,:) = 0.5*(s%ua(1:s%nx,:)*sinthm(1:s%nx,:)+s%ua(2:s%nx+1,:)*sinthm(1:s%nx,:)) 00221 if (par%nz>1) then 00222 ueu_sed(1:s%nx,:) = 0.5*(s%ue_sed(1:s%nx,:)+s%ue_sed(2:s%nx+1,:)) 00223 else 00224 ueu_sed=s%ueu 00225 endif 00226 veu_sed(1:s%nx,:) = 0.5*(s%ve_sed(1:s%nx,:)+s%ve_sed(2:s%nx+1,:)) 00227 if (xmpi_isbot) then 00228 veu_sed(s%nx+1,:) = veu_sed(s%nx,:) 00229 endif 00230 00231 ! Compute vmagu including ua 00232 ! s%vmagu = sqrt((s%uu+uau)**2+(s%vu+uav)**2) 00233 s%vmagu = sqrt((ueu_sed+uau)**2+(veu_sed+uav)**2) 00234 ! sediment advection velocity for suspended load and bed load respectively 00235 ! REMARK: when vreps does not equal vv; no mass conservation 00236 ! s%ureps = s%ueu+uau 00237 ! s%urepb = s%ueu+uau ! RJ maybe reduce this velocity? 00238 s%ureps = ueu_sed+uau 00239 s%urepb = ueu_sed+uau 00240 ! 00241 do j=1,s%ny+1 00242 do i=1,s%nx 00243 if(s%ureps(i,j)>0.d0) then 00244 ! test cu(i,j)=cc(i,j) 00245 cu(i,j)=par%thetanum*cc(i,j)+(1.d0-par%thetanum)*cc(min(i+1,s%nx),j) 00246 cub(i,j)=par%thetanum*s%pbbed(i,j,1,jg)*s%ceqbg(i,j,jg)+(1.d0-par%thetanum)& 00247 *s%pbbed(min(i+1,s%nx),j,1,jg)*s%ceqbg(min(i+1,s%nx),j,jg) 00248 !cub(i,j)=par%thetanum*ccb(i,j)+(1.d0-par%thetanum)*ccb(min(i+1,nx),j) 00249 elseif(s%ureps(i,j)<0.d0) then 00250 cu(i,j)=par%thetanum*cc(i+1,j)+(1.d0-par%thetanum)*cc(max(i,2),j) 00251 cub(i,j)=par%thetanum*s%pbbed(i+1,j,1,jg)*s%ceqbg(i+1,j,jg)+(1.d0-par%thetanum)& 00252 *s%pbbed(max(i,2),j,1,jg)*s%ceqbg(max(i,2),j,jg) 00253 !cub(i,j)=par%thetanum*ccb(i+1,j)+(1.d0-par%thetanum)*ccb(max(i,2),j) 00254 else 00255 cu(i,j)=0.5d0*(cc(i,j)+cc(i+1,j)) 00256 cub(i,j)=0.5d0*(s%pbbed(i,j,1,jg)*s%ceqbg(i,j,jg)+s%pbbed(i+1,j,1,jg)*s%ceqbg(i+1,j,jg)) 00257 !cub(i,j)=0.5d0*(ccb(i,j)+ccb(i+1,j)) 00258 endif 00259 s%dcsdx(i,j)=(cc(i+1,j)-cc(i,j))/s%dsu(i,j) 00260 enddo 00261 enddo 00262 ! wwvv dcdx(nx:1,:) is still untouched, correct this ofr the parallel case 00263 cu(s%nx+1,:) = cc(s%nx+1,:) !Robert 00264 ! wwvv fix this in parallel case 00265 ! wwvv in parallel version, there will be a discrepancy between the values 00266 ! of dzbdx(nx+1,:). 00267 !wwvv so fix that 00268 ! 00269 Sus = 0.d0 00270 Sub = 0.d0 00271 ! 00272 ! suspended load, Lodewijk: no bed slope effect (yet) 00273 Sus=par%sus*(cu*s%ureps*s%hu-s%Dc*s%hu*s%dcsdx)*s%wetu 00274 ! bed load, Lodewijk: no bed slope effect (yet) 00275 Sub=par%bed*(cub*s%urepb*s%hu)*s%wetu 00276 ! 00277 ! 00278 ! Y-direction 00279 ! 00280 ! Jaap: get ua in v points and split out in u and v direction 00281 if (s%ny>0) then 00282 uau(:,1:s%ny) = 0.5*(s%ua(:,1:s%ny)*costhm(:,1:s%ny)+s%ua(:,2:s%ny+1)*costhm(:,1:s%ny)) 00283 uav(:,1:s%ny) = 0.5*(s%ua(:,1:s%ny)*sinthm(:,1:s%ny)+s%ua(:,2:s%ny+1)*sinthm(:,1:s%ny)) 00284 uau(:,s%ny+1) = uau(:,s%ny) ! Jaap 00285 uav(:,s%ny+1) = uav(:,s%ny) ! Jaap 00286 if (par%nz>1) then 00287 vev_sed(:,1:s%ny) = 0.5*(s%ve_sed(:,1:s%ny)+s%ve_sed(:,2:s%ny+1)) 00288 else 00289 vev_sed=s%vev 00290 endif 00291 uev_sed(:,1:s%ny) = 0.5*(s%ue_sed(:,1:s%ny)+s%ue_sed(:,2:s%ny+1)) 00292 if (xmpi_isright) then 00293 uev_sed(:,1:s%ny+1) = uev_sed(:,1:s%ny) 00294 endif 00295 else 00296 uau=s%ua*costhm 00297 uav=s%ua*sinthm 00298 uev_sed=s%ue_sed 00299 vev_sed=s%ve_sed 00300 endif 00301 ! Jaap: compute vmagv including ua 00302 ! s%vmagv = sqrt((s%uv+uau)**2+(s%vv+uav)**2) 00303 s%vmagv = sqrt((uev_sed+uau)**2+(vev_sed+uav)**2) 00304 ! sediment advection velocity for suspended load and bed load respectively 00305 ! REMARK: when vreps does not equal vv; no mass conservation 00306 ! s%vreps = s%vev+uav 00307 ! s%vrepb = s%vev+uav ! RJ maybe reduce this velocity? Should be s%vv instead of s%vev? 00308 s%vreps = vev_sed+uav 00309 s%vrepb = vev_sed+uav ! RJ maybe reduce this velocity? Should be s%vv instead of s%vev? 00310 ! 00311 if (s%ny>0) then 00312 do j=1,s%ny 00313 do i=1,s%nx+1 00314 if(s%vreps(i,j)>0) then 00315 cv(i,j)=par%thetanum*cc(i,j)+(1.d0-par%thetanum)*cc(i,min(j+1,s%ny)) 00316 cvb(i,j)=par%thetanum*s%pbbed(i,j,1,jg)*s%ceqbg(i,j,jg)+(1.d0-par%thetanum)& 00317 *s%pbbed(i,min(j+1,s%ny),1,jg)*s%ceqbg(i,min(j+1,s%ny),jg) 00318 !cvb(i,j)=par%thetanum*ccb(i,j)+(1.d0-par%thetanum)*ccb(i,min(j+1,ny)) 00319 elseif(s%vreps(i,j)<0) then 00320 cv(i,j)=par%thetanum*cc(i,j+1)+(1.d0-par%thetanum)*cc(i,max(j,2)) 00321 cvb(i,j)=par%thetanum*s%pbbed(i,j+1,1,jg)*s%ceqbg(i,j+1,jg)+(1.d0-par%thetanum)& 00322 *s%pbbed(i,max(j,2),1,jg)*s%ceqbg(i,max(j,2),jg) 00323 !cvb(i,j)=par%thetanum*ccb(i,j+1)+(1.d0-par%thetanum)*ccb(i,max(j,2)) 00324 else 00325 cv(i,j)=0.5d0*(cc(i,j)+cc(i,j+1)) !Jaap: cc instead of cv 00326 cvb(i,j)=0.5d0*(s%pbbed(i,j,1,jg)*s%ceqbg(i,j,jg)+s%pbbed(i,j+1,1,jg)*s%ceqbg(i,j+1,jg)) 00327 !cvb(i,j)=0.5d0*(ccb(i,j)+ccb(i,j+1)) 00328 end if 00329 s%dcsdy(i,j)=(cc(i,j+1)-cc(i,j))/s%dnv(i,j) !Jaap 00330 00331 end do 00332 end do 00333 ! wwvv dcdy(:,ny+1) is not filled in, so in parallel case: 00334 cv(:,s%ny+1) = cc(:,s%ny+1) !Robert 00335 ! wwvv in parallel version, there will be a discrepancy between the values 00336 ! of dzbdy(:,ny+1). 00337 ! wwvv so fix that 00338 else 00339 cv = cc 00340 cvb = s%ceqbg(:,:,jg) 00341 endif ! s%ny>0 00342 ! 00343 ! Compute sedimnent transport in v-direction 00344 ! 00345 Svs = 0.d0 00346 Svb = 0.d0 00347 ! Suspended load 00348 Svs=par%sus*(cv*s%vreps*s%hv-s%Dc*s%hv*s%dcsdy)*s%wetv 00349 ! Bed load 00350 Svb=par%bed*(cvb*s%vrepb*s%hv)*s%wetv 00351 ! 00352 ! 00353 ! Bed slope effects and bermslope model 00354 ! 00355 ! Where bermslope model should run 00356 if (par%bermslopetransport==1) then ! only update from false if bermslope model is used 00357 if (par%wavemodel == WAVEMODEL_SURFBEAT) then 00358 where (s%H/s%hu>par%bermslopegamma .or. s%hu<par%bermslopedepth) 00359 bermslopeindex = .true. 00360 elsewhere 00361 bermslopeindex = .false. 00362 endwhere 00363 else 00364 where (s%hu<par%bermslopedepth) 00365 bermslopeindex = .true. 00366 elsewhere 00367 bermslopeindex = .false. 00368 endwhere 00369 endif 00370 if(par%bermslopebed==1) then ! only update from false if used 00371 bermslopeindexbed = bermslopeindex 00372 endif 00373 if(par%bermslopesus==1) then ! only update from false if used 00374 bermslopeindexsus = bermslopeindex 00375 endif 00376 endif 00377 ! 00378 ! 00379 ! Do bermslope transports in bermslope locations [separated from previous to reduce block size] 00380 where(bermslopeindexbed) 00381 Sub = Sub-par%bed*(par%bermslopefac*cub*s%vmagu*s%hu*(s%dzbdx-par%bermslope))*s%wetu 00382 endwhere 00383 where(bermslopeindexsus) 00384 Sus = Sus-par%sus*(par%bermslopefac* cu*s%vmagu*s%hu*(s%dzbdx-par%bermslope))*s%wetu 00385 endwhere 00386 ! 00387 ! 00388 ! Do regular bed slope effects in other locations 00389 if (par%bdslpeffmag == BDSLPEFFMAG_ROELV_TOTAL) then 00390 where(.not. bermslopeindexbed) 00391 Sub = Sub-par%bed*(par%facsl*cub*s%vmagu*s%hu*s%dzbdx)*s%wetu 00392 endwhere 00393 where (.not. bermslopeindexsus) 00394 Sus = Sus-par%sus*(par%facsl*cu*s%vmagu*s%hu*s%dzbdx)*s%wetu 00395 endwhere 00396 Svb = Svb-par%bed*(par%facsl*cvb*s%vmagv*s%hv*s%dzbdy)*s%wetv 00397 Svs = Svs-par%sus*(par%facsl*cv*s%vmagv*s%hv*s%dzbdy)*s%wetv 00398 elseif (par%bdslpeffmag == BDSLPEFFMAG_ROELV_BED) then 00399 where(.not. bermslopeindexbed) 00400 Sub = Sub-par%bed*(par%facsl*cub*s%vmagu*s%hu*s%dzbdx)*s%wetu 00401 endwhere 00402 Svb = Svb-par%bed*(par%facsl*cvb*s%vmagv*s%hv*s%dzbdy)*s%wetv 00403 elseif (par%bdslpeffmag == BDSLPEFFMAG_SOULS_TOTAL .or. par%bdslpeffmag == BDSLPEFFMAG_SOULS_BED) then 00404 ! 00405 ! Bed slope magnitude effect (as Souslby intended) and change direction transport (see Van Rijn 1993 (section 7.2.6)) 00406 ! 00407 do j=1,s%ny+1 00408 do i=1,s%nx+1 00409 if ((dabs(Sub(i,j)) > 0.000001d0) .or. (dabs(Svb(i,j)) > 0.000001d0) .and. (.not. bermslopeindexbed(i,j)) ) then 00410 Sbmtot = dsqrt( Sub(i,j)**2.d0 + Svb(i,j)**2.d0 ) 00411 dzbds = s%dzbdx(i,j)*Sub(i,j)/Sbmtot + s%dzbdy(i,j)*Svb(i,j)/Sbmtot 00412 ! dzbdn = s%dzbdx*Svb/Sbtot + s%dzbdy*Sub/Sbtot 00413 Sub(i,j) = Sub(i,j)*(1.d0 - par%facsl*dzbds) 00414 Svb(i,j) = Svb(i,j)*(1.d0 - par%facsl*dzbds) 00415 ! 00416 endif 00417 if (((dabs(Sus(i,j)) > 0.000001d0) .or. (dabs(Svs(i,j)) > 0.000001d0)) .and. (.not. bermslopeindexsus(i,j)) & 00418 .and. par%bdslpeffmag == BDSLPEFFMAG_SOULS_TOTAL) then 00419 Ssmtot = dsqrt( Sus(i,j)**2.d0 + Svs(i,j)**2.d0 ) 00420 dzbds = s%dzbdx(i,j)*Sus(i,j)/Ssmtot + s%dzbdy(i,j)*Svs(i,j)/Ssmtot 00421 ! dzbdn = s%dzbdx*Svb/Sbtot + s%dzbdy*Sub/Sbtot 00422 Sus(i,j) = Sus(i,j)*(1.d0 - par%facsl*dzbds) 00423 Svs(i,j) = Svs(i,j)*(1.d0 - par%facsl*dzbds) 00424 endif 00425 enddo 00426 enddo 00427 endif 00428 ! 00429 ! Lodewijk: modify the direction of the bed load transport based on the bed slope, see Van Rijn 1993 (section 7.2.6) 00430 if (par%bdslpeffdir == BDSLPEFFDIR_TALMON) then 00431 do j=1,s%ny+1 00432 do i=1,s%nx+1 00433 if (((dabs(s%urepb(i,j)) > 0.0001d0) .or. (dabs(s%vrepb(i,j)) > 0.0001d0)) & 00434 .and. ((dabs(s%taubx(i,j)) > 0.0001d0) .or. (dabs(s%tauby(i,j)) > 0.0001d0))) then 00435 if (s%urepb(i,j) < 0.d0) then 00436 delta_x = datan(s%vrepb(i,j)/s%urepb(i,j))+par%px ! Angle between fluid velocity vector and the s%x-axis 00437 else 00438 delta_x = datan(s%vrepb(i,j)/s%urepb(i,j)) ! Angle between fluid velocity vector and the s%x-axis 00439 endif 00440 delta = (par%rhos-par%rho)/par%rho 00441 shields = sqrt(s%taubx(i,j)**2 + s%tauby(i,j)**2)/(delta*par%rho*par%g*s%D50(jg)) 00442 ! shields = (urepb(i,j)**2.d0+vrepb(i,j)**2.d0)*s%cf(i,j)/(par%g*D50(jg)*delta) 00443 ftheta = 1.d0/(9.d0*(s%D50(jg)/s%hh(i,j))**0.3d0*sqrt(shields)) ! Talmon 00444 psi_x = datan( (dsin(delta_x)-ftheta*s%dzbdy(i,j)) / (dcos(delta_x)-ftheta*s%dzbdx(i,j)) ) 00445 psi_x = par%bdslpeffdirfac*(psi_x - delta_x)+delta_x 00446 Sbtot = dsqrt( Sub(i,j)**2.d0 + Svb(i,j)**2.d0 ) ! Magnitude of sediment transport without direction modifcation 00447 ! Decompose the sediment transport again, know with the knowledge of the direction of the sediment transport vector 00448 Sub(i,j) = Sbtot * dcos(psi_x) 00449 Svb(i,j) = Sbtot * dsin(psi_x) 00450 else 00451 Sub(i,j) = 0.d0 00452 Svb(i,j) = 0.d0 00453 endif 00454 enddo 00455 enddo 00456 endif 00457 ! 00458 ! 00459 ! 00460 do j=1,s%ny+1 00461 do i=1,s%nx 00462 if(Sub(i,j)>0.d0) then 00463 pbbedu(i,j) = s%pbbed(i,j,1,jg) 00464 elseif(Sub(i,j)<0.d0) then 00465 pbbedu(i,j)= s%pbbed(i+1,j,1,jg) 00466 else 00467 pbbedu(i,j)=0.5d0*(s%pbbed(i,j,1,jg)+s%pbbed(i+1,j,1,jg)) 00468 endif 00469 enddo 00470 enddo 00471 ! 00472 Sub = pbbedu*Sub 00473 ! 00474 do j=1,s%ny 00475 do i=1,s%nx+1 00476 if(Svb(i,j)>0) then 00477 pbbedv(i,j)=s%pbbed(i,j,1,jg) 00478 else if(Svb(i,j)<0) then 00479 pbbedv(i,j)=s%pbbed(i,j+1,1,jg) 00480 else 00481 pbbedv(i,j)=0.5d0*(s%pbbed(i,j,1,jg)+s%pbbed(i,j+1,1,jg)) 00482 end if 00483 end do 00484 end do 00485 ! 00486 Svb = pbbedv*Svb 00487 ! 00488 ! BRJ: implicit concentration update (compute sources first, sink must be computed after updating actual sed.conc.) 00489 ! 00490 if (s%ny>0) then 00491 do j=jmin_zs,jmax_zs 00492 do i=imin_zs,imax_zs 00493 ! Changed to hh from hold by RJ (13072009) !**2/max(hh(i,j),par%hmin) 00494 s%ero(i,j,jg) = fac(i,j,jg)*s%hh(i,j)*s%ceqsg(i,j,jg)*s%pbbed(i,j,1,jg)/s%Tsg(i,j,jg) 00495 ! depo_ex(i,j,jg) = max(hold(i,j),0.01d0)*cc(i,j)/Tsg(i,j,jg) 00496 ! BRJ: the volume in the water column is updated and not the volume concentration. 00497 cc(i,j) = (par%dt*s%Tsg(i,j,jg))/(par%dt+s%Tsg(i,j,jg))* & 00498 (s%hold(i,j)*cc(i,j)/par%dt -((Sus(i,j)*s%dnu(i,j)-Sus(i-1,j)*s%dnu(i-1,j)+& 00499 Svs(i,j)*s%dsv(i,j)-Svs(i,j-1)*s%dsv(i,j-1))*s%dsdnzi(i,j)-& 00500 s%ero(i,j,jg))) 00501 00502 cc(i,j)=max(cc(i,j),0.0d0) ! Jaap: negative cc's are possible... 00503 cc(i,j)=min(cc(i,j),par%cmax*s%hh(i,j)) 00504 s%depo_ex(i,j,jg) = cc(i,j)/s%Tsg(i,j,jg) 00505 enddo 00506 enddo 00507 00508 else 00509 j=1 00510 do i=imin_zs,imax_zs 00511 ! Changed to hh from hold by RJ (13072009) !**2/max(hh(i,j),par%hmin) 00512 s%ero(i,j,jg) = fac(i,j,jg)*s%hh(i,j)*s%ceqsg(i,j,jg)*s%pbbed(i,j,1,jg)/s%Tsg(i,j,jg) 00513 ! depo_ex(i,j,jg) = max(hold(i,j),0.01d0)*cc(i,j)/Tsg(i,j,jg) 00514 ! BRJ: the volume in the water column is updated and not the volume concentration. 00515 cc(i,j) = (par%dt*s%Tsg(i,j,jg))/(par%dt+s%Tsg(i,j,jg))* & 00516 (s%hold(i,j)*cc(i,j)/par%dt -((Sus(i,j)*s%dnu(i,j)-Sus(i-1,j)*s%dnu(i-1,j))*s%dsdnzi(i,j)-& 00517 s%ero(i,j,jg))) 00518 00519 cc(i,j)=max(cc(i,j),0.0d0) ! Jaap: negative cc's are possible... 00520 cc(i,j)=min(cc(i,j),par%cmax*s%hh(i,j)) 00521 00522 s%depo_ex(i,j,jg) = cc(i,j)/s%Tsg(i,j,jg) 00523 enddo 00524 endif 00525 00526 00527 cc(imin_zs:imax_zs,jmin_zs:jmax_zs) = cc(imin_zs:imax_zs,jmin_zs:jmax_zs)/s%hh(imin_zs:imax_zs,jmin_zs:jmax_zs) 00528 00529 ! do lateral boundaries... 00530 if(xmpi_istop)then 00531 cc(1,:)=cc(2,:) 00532 s%ero(1,:,jg)=s%ero(2,:,jg) 00533 s%depo_ex(1,:,jg)=s%depo_ex(2,:,jg) 00534 endif 00535 if(xmpi_isleft .and. s%ny>0)then 00536 cc(:,1)=cc(:,2) 00537 s%ero(:,1,jg)=s%ero(:,2,jg) 00538 s%depo_ex(:,1,jg)=s%depo_ex(:,2,jg) 00539 endif 00540 if(xmpi_isbot)then 00541 cc(s%nx+1,:)=cc(s%nx,:) 00542 s%ero(s%nx+1,:,jg)=s%ero(s%nx,:,jg) 00543 s%depo_ex(s%nx+1,:,jg)=s%depo_ex(s%nx,:,jg) 00544 endif 00545 if(xmpi_isright .and. s%ny>0)then 00546 cc(:,s%ny+1)=cc(:,s%ny) 00547 s%ero(:,s%ny+1,jg)=s%ero(:,s%ny,jg) 00548 s%depo_ex(:,s%ny+1,jg)=s%depo_ex(:,s%ny,jg) 00549 endif 00550 ! Dano: fix nasty numerics 00551 where (cc<1.d-10) 00552 cc=0.d0 00553 endwhere 00554 ! wwvv fix the first and last rows and columns of cc in parallel case 00555 #ifdef USEMPI 00556 call xmpi_shift_ee(cc) 00557 #endif 00558 ! wwvv border columns and rows of ccg Svg and Sug have to be communicated 00559 ! 00560 ! An additional sensitivity to D50 (called par%alfaD50) 00561 ! By default par%alfaD50 = 0, so no additional sensitivity 00562 if (par%alfaD50 > 0) then 00563 do j=1,s%ny+1 00564 do i=1,s%nx+1 00565 Svs(i,j) = Svs(i,j) * (0.000225/par%D50(jg))**par%alfaD50 00566 Sus(i,j) = Sus(i,j) * (0.000225/par%D50(jg))**par%alfaD50 00567 Svb(i,j) = Svb(i,j) * (0.000225/par%D50(jg))**par%alfaD50 00568 Sub(i,j) = Sub(i,j) * (0.000225/par%D50(jg))**par%alfaD50 00569 enddo 00570 enddo 00571 endif 00572 00573 ! Save concentrations and sediment transport in s 00574 s%ccg(:,:,jg) = cc 00575 s%Svsg(:,:,jg) = Svs 00576 s%Susg(:,:,jg) = Sus 00577 s%Svbg(:,:,jg) = Svb 00578 s%Subg(:,:,jg) = Sub 00579 00580 end do ! number of sediment fractions 00581 00582 s%vmag=sqrt(max(vmag2,par%umin)) 00583 00584 end subroutine transus 00585 00586 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00587 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00588 00589 subroutine bed_update(s,par) 00590 use params, only: parameters 00591 use spaceparamsdef 00592 use spaceparams 00593 use xmpi_module 00594 00595 implicit none 00596 00597 type(spacepars),target :: s 00598 type(parameters) :: par 00599 00600 integer :: i,j,j1,jg,ii,ie,id,je,jd,jdz,ndz, hinterland 00601 integer , dimension(:,:,:),allocatable,save :: indSus,indSub,indSvs,indSvb 00602 real*8 :: dzb,dzmax,dzt,dzleft,sdz,dzavt,fac,Savailable,dAfac 00603 real*8 , dimension(:,:),allocatable,save :: Sout,hav,tempexchange 00604 real*8 , dimension(par%ngd) :: edg,edg1,edg2,dzg 00605 real*8 , dimension(:),pointer :: dz 00606 real*8 , dimension(:,:),pointer :: pb 00607 integer :: n_aval,im1,ip1,jm1,jp1 00608 real*8,save :: delta 00609 real*8 :: totslp 00610 00611 !include 's.ind' 00612 !include 's.inp' 00613 00614 if (.not. allocated(Sout)) then 00615 allocate(Sout(s%nx+1,s%ny+1)) 00616 allocate(hav(s%nx+1,s%ny+1)) 00617 allocate(indSus(s%nx+1,s%ny+1,par%ngd)) 00618 allocate(indSub(s%nx+1,s%ny+1,par%ngd)) 00619 allocate(indSvs(s%nx+1,s%ny+1,par%ngd)) 00620 allocate(indSvb(s%nx+1,s%ny+1,par%ngd)) 00621 if (par%gwflow==1) then 00622 allocate(tempexchange(s%nx+1,s%ny+1)) 00623 endif 00624 delta = (par%rhos-par%rho)/par%rho 00625 endif 00626 00627 ! Super fast 1D 00628 if (s%ny==0) then 00629 j1 = 1 00630 else 00631 j1 = 2 00632 endif 00633 s%dzbnow = 0.d0 00634 dzb = 0.d0 00635 if (par%t>=par%morstart .and. par%t < par%morstop .and. par%morfac > .999d0) then 00636 ! 00637 ! bed_predict 00638 ! 00639 ! reduce sediment transports when hard layer comes to surface 00640 ! this step is mainly necessary at the transition from hard layers to sand 00641 if (par%struct == 1) then 00642 00643 do jg = 1,par%ngd 00644 indSus = 0 00645 indSub = 0 00646 indSvs = 0 00647 indSvb = 0 00648 Sout = 0.d0 00649 do j=j1,s%ny+1 00650 do i=2,s%nx+1 00651 ! fluxes at i,j 00652 if (s%Subg(i,j,jg) > 0.d0) then ! bed load s%u-direction 00653 indSub(i,j,jg) = 1 00654 Sout(i,j) = Sout(i,j) + s%Subg(i,j,jg)*s%dnu(i,j) 00655 endif 00656 ! fluxes at i-1,j 00657 if (s%Subg(i-1,j,jg) < 0.d0 ) then ! bed load s%u-direction 00658 Sout(i,j) = Sout(i,j) - s%Subg(i-1,j,jg)*s%dnu(i-1,j) 00659 endif 00660 if (par%sourcesink==0) then 00661 ! fluxes at i,j 00662 if (s%Susg(i,j,jg) > 0.d0 ) then ! suspended load s%u-direction 00663 indSus(i,j,jg) = 1 00664 Sout(i,j) = Sout(i,j) + s%Susg(i,j,jg)*s%dnu(i,j) 00665 endif 00666 ! fluxes at i-1,j 00667 if (s%Susg(i-1,j,jg) < 0.d0 ) then ! suspended load s%u-direction 00668 Sout(i,j) = Sout(i,j) - s%Susg(i-1,j,jg)*s%dnu(i-1,j) 00669 endif 00670 endif 00671 enddo 00672 enddo 00673 if (s%ny>0) then 00674 do j=j1,s%ny+1 00675 do i=2,s%nx+1 00676 if (s%Svbg(i,j,jg) > 0.d0 ) then ! bed load s%v-direction 00677 indSvb(i,j,jg) = 1 00678 Sout(i,j) = Sout(i,j) + s%Svbg(i,j,jg)*s%dsv(i,j) 00679 endif 00680 ! fluxes at i,j-1 00681 if (s%Svbg(i,j-1,jg) < 0.d0 ) then ! bed load s%v-direction 00682 Sout(i,j) = Sout(i,j) - s%Svbg(i,j-1,jg)*s%dsv(i,j-1) 00683 endif 00684 if (par%sourcesink==0) then 00685 if (s%Svsg(i,j,jg) > 0.d0 ) then ! suspended load s%v-direction 00686 indSvs(i,j,jg) = 1 00687 Sout(i,j) = Sout(i,j) + s%Svsg(i,j,jg)*s%dsv(i,j) 00688 endif 00689 ! fluxes at i,j-1 00690 if (s%Svsg(i,j-1,jg) < 0.d0 ) then ! suspended load s%v-direction 00691 Sout(i,j) = Sout(i,j) - s%Svsg(i,j-1,jg)*s%dsv(i,j-1) 00692 endif 00693 endif ! sourcesink = 0 00694 enddo !s%nx+1 00695 enddo !s%ny+1 00696 endif !s%ny>0 00697 ! 00698 do j=j1,s%ny+1 00699 do i=2,s%nx+1 00700 Savailable = s%structdepth(i,j)*s%pbbed(i,j,1,jg)/par%morfac/par%dt*(1.d0-par%por)/s%dsdnzi(i,j) 00701 ! reduction factor for cell outgoing sediment transports wwvv changed tiny into epsilon 00702 fac = min(1.d0,Savailable/max(Sout(i,j),epsilon(0.d0)) ) 00703 ! fix sediment transports for the presence of a hard layer; remind indSus etc are 1 in cases of cell outgoing transports 00704 ! updated S oell outgoing transports cell incoming transports 00705 if (fac < 1.d0)then 00706 s%Subg(i,j,jg) = fac*indSub(i,j,jg)*s%Subg(i,j,jg) + (1-indSub(i,j,jg))*s%Subg(i,j,jg) 00707 s%Subg(i-1,j,jg) = fac*(1-indSub(i-1,j,jg))*s%Subg(i-1,j,jg) + indSub(i-1,j,jg)*s%Subg(i-1,j,jg) 00708 if (s%ny>0) then 00709 s%Svbg(i,j,jg) = fac*indSvb(i,j,jg)*s%Svbg(i,j,jg) + (1-indSvb(i,j,jg))*s%Svbg(i,j,jg) 00710 s%Svbg(i,j-1,jg) = fac*(1-indSvb(i,j-1,jg))*s%Svbg(i,j-1,jg) + indSvb(i,j-1,jg)*s%Svbg(i,j-1,jg) 00711 endif 00712 if (par%sourcesink==0) then 00713 s%Susg(i,j,jg) = fac*indSus(i,j,jg)*s%Susg(i,j,jg) + (1-indSus(i,j,jg))*s%Susg(i,j,jg) 00714 s%Susg(i-1,j,jg) = fac*(1-indSus(i-1,j,jg))*s%Susg(i-1,j,jg) + indSus(i-1,j,jg)*s%Susg(i-1,j,jg) 00715 if (s%ny>0) then 00716 s%Svsg(i,j,jg) = fac*indSvs(i,j,jg)*s%Svsg(i,j,jg) + (1-indSvs(i,j,jg))*s%Svsg(i,j,jg) 00717 s%Svsg(i,j-1,jg) = fac*(1-indSvs(i,j-1,jg))*s%Svsg(i,j-1,jg) + indSvs(i,j-1,jg)*s%Svsg(i,j-1,jg) 00718 endif !s%ny = 0 00719 endif ! sourcesink = 0 00720 endif !fac<1.d0 00721 enddo ! s%nx+1 00722 enddo !s%ny + 1 00723 enddo !par%ngd 00724 endif !struct == 1 00725 00726 if (s%ny>0) then 00727 do j=jmin_zs,jmax_zs 00728 do i=imin_zs,imax_zs 00729 00730 ! bed level changes per fraction in this morphological time step in meters sand including pores 00731 ! positive in case of erosion 00732 if (par%sourcesink==0) then 00733 dzg=par%morfac*par%dt/(1.d0-par%por)*( & ! Dano, dz from sus transport gradients 00734 ( s%Susg(i,j,:)*s%dnu(i,j)-s%Susg(i-1,j,:)*s%dnu(i-1,j) +& 00735 s%Svsg(i,j,:)*s%dsv(i,j)-s%Svsg(i,j-1,:)*s%dsv(i,j-1) +& 00736 ! dz from bed load transport gradients 00737 s%Subg(i,j,:)*s%dnu(i,j)-s%Subg(i-1,j,:)*s%dnu(i-1,j)+& 00738 s%Svbg(i,j,:)*s%dsv(i,j)-s%Svbg(i,j-1,:)*s%dsv(i,j-1) )*s%dsdnzi(i,j) ) 00739 elseif (par%sourcesink==1) then 00740 dzg=par%morfac*par%dt/(1.d0-par%por)*( & 00741 s%ero(i,j,:)-s%depo_ex(i,j,:) +& 00742 ( s%Subg(i,j,:)*s%dnu(i,j)-s%Subg(i-1,j,:)*s%dnu(i-1,j)+& 00743 s%Svbg(i,j,:)*s%dsv(i,j)-s%Svbg(i,j-1,:)*s%dsv(i,j-1) )*s%dsdnzi(i,j) ) 00744 endif 00745 00746 if (par%ngd==1) then ! Simple bed update in case one fraction 00747 00748 s%zb(i,j) = s%zb(i,j)-sum(dzg) 00749 s%dzbnow(i,j) = s%dzbnow(i,j)-sum(dzg) ! naamgeveing? 00750 s%dzbdt(i,j) = s%dzbnow(i,j)/par%dt 00751 s%sedero(i,j) = s%sedero(i,j)-sum(dzg) 00752 s%structdepth(i,j) = max(0.d0,s%structdepth(i,j)-sum(dzg)) 00753 00754 else 00755 ! erosion/deposition rate of sand mass (m/s) 00756 ! positive in case of erosion 00757 edg = dzg*(1.d0-par%por)/par%dt 00758 00759 00760 dz=>s%dzbed(i,j,:) 00761 pb=>s%pbbed(i,j,:,:) 00762 00763 call update_fractions(par,s,i,j,dz,pb,edg,sum(dzg)) 00764 00765 endif 00766 00767 enddo ! s%nx+1 00768 enddo ! s%ny+1 00769 else 00770 j=1 00771 do i=imin_zs,imax_zs 00772 ! bed level changes per fraction in this morphological time step in meters sand including pores 00773 ! positive in case of erosion 00774 if (par%sourcesink==0) then 00775 dzg=par%morfac*par%dt/(1.d0-par%por)*( & ! Dano, dz from sus transport gradients 00776 ( s%Susg(i,j,:)*s%dnu(i,j)-s%Susg(i-1,j,:)*s%dnu(i-1,j) +& 00777 s%Subg(i,j,:)*s%dnu(i,j)-s%Subg(i-1,j,:)*s%dnu(i-1,j) )*s%dsdnzi(i,j) +& 00778 (s%Svsg(i,j,:)+s%Svbg(i,j,:))*par%lsgrad) 00779 elseif (par%sourcesink==1) then 00780 dzg=par%morfac*par%dt/(1.d0-par%por)*( & 00781 s%ero(i,j,:)-s%depo_ex(i,j,:) +& 00782 ( s%Subg(i,j,:)*s%dnu(i,j)-s%Subg(i-1,j,:)*s%dnu(i-1,j) )*s%dsdnzi(i,j) +& 00783 (s%Svsg(i,j,:)+s%Svbg(i,j,:))*par%lsgrad) 00784 endif 00785 00786 if (par%ngd==1) then ! Simple bed update in case one fraction 00787 00788 s%zb(i,j) = s%zb(i,j)-sum(dzg) 00789 s%dzbnow(i,j) = s%dzbnow(i,j)-sum(dzg) 00790 s%dzbdt(i,j) = s%dzbnow(i,j)/par%dt 00791 s%sedero(i,j) = s%sedero(i,j)-sum(dzg) 00792 s%structdepth(i,j) = max(0.d0,s%structdepth(i,j)-sum(dzg)) 00793 00794 else ! multiple fractions... 00795 ! erosion/deposition rate of sand mass (m/s) 00796 ! positive in case of erosion 00797 edg = dzg*(1.d0-par%por)/par%dt 00798 00799 dz=>s%dzbed(i,j,:) 00800 pb=>s%pbbed(i,j,:,:) 00801 00802 call update_fractions(par,s,i,j,dz,pb,edg,sum(dzg)) 00803 00804 endif 00805 00806 enddo ! s%nx+1 00807 endif !s%ny>0 00808 #ifdef USEMPI 00809 call xmpi_shift_ee(s%zb) 00810 call xmpi_shift_ee(s%structdepth) 00811 #endif 00812 ! In case of groundwater, we need to update groundwater levels and surface water 00813 if(par%gwflow==1) then 00814 where (s%wetz==1 .and. s%dzbnow>0.d0) ! accretion in wet areas 00815 s%zs = s%zs + s%dzbnow*(1.d0-par%por) ! zs = zs + dzbnow - dzbnow*par%por 00816 s%gwlevel = s%gwlevel+s%dzbnow 00817 s%infil = s%infil + s%dzbnow*par%por/par%dt 00818 elsewhere (s%wetz==1 .and. s%dzbnow<0.d0) ! erosion in wet areas 00819 ! maximum water leaving groundwater = gwlevel(not updated) - zb(updated) 00820 ! note that exfiltration is negative infil 00821 tempexchange = min((s%zb-s%gwlevel)*par%por,0.d0) 00822 s%zs = s%zs + s%dzbnow - tempexchange 00823 s%gwlevel = s%gwlevel + tempexchange/par%por 00824 s%infil = s%infil + tempexchange 00825 endwhere 00826 #ifdef USEMPI 00827 call xmpi_shift_ee(s%gwlevel) 00828 call xmpi_shift_ee(s%infil) 00829 call xmpi_shift_ee(s%zs) 00830 #endif 00831 else ! only need to update water levels 00832 where (s%wetz(imin_zs:imax_zs,jmin_zs:jmax_zs)==1) 00833 s%zs(imin_zs:imax_zs,jmin_zs:jmax_zs) = s%zs(imin_zs:imax_zs,jmin_zs:jmax_zs)+s%dzbnow(imin_zs:imax_zs,jmin_zs:jmax_zs) 00834 endwhere 00835 #ifdef USEMPI 00836 call xmpi_shift_ee(s%zs) 00837 #endif 00838 endif 00839 00840 00841 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00842 ! 00843 ! Avalanching 00844 ! 00845 00846 00847 if (par%avalanching==1) then 00848 00849 if(par%fixedavaltime==0) then 00850 par%avaltime = par%nTrepavaltime*par%Trep 00851 endif 00852 00853 do ii=1,nint(par%morfac) 00854 !aval=.false. 00855 n_aval=0 !Dano fix for MPI 00856 do j=1,s%ny+1 00857 do i=1,s%nx 00858 s%dzbdx(i,j)=(s%zb(i+1,j)-s%zb(i,j))/s%dsu(i,j) 00859 enddo 00860 s%dzbdx(s%nx+1,j)=s%dzbdx(s%nx,j) 00861 enddo 00862 do j=1,s%ny 00863 do i=1,s%nx+1 00864 s%dzbdy(i,j)=(s%zb(i,j+1)-s%zb(i,j))/s%dnv(i,j) 00865 enddo 00866 enddo 00867 if (s%ny>0) then 00868 do i=1,s%nx+1 00869 s%dzbdy(i,s%ny+1)=s%dzbdy(i,s%ny) 00870 enddo 00871 endif 00872 ! 00873 hav = s%hh 00874 ! Fix Hav for short wave runup: 00875 if (par%swrunup == 1) then 00876 do j = 1,s%ny+1 00877 hinterland = 0 00878 do i = 1, s%nx+1 00879 if (hinterland == 0 .and. s%runup(j)+s%zs(nint(s%iwl(j)),j) < s%zb(i,j) + par%eps) then 00880 hinterland = 1; 00881 endif 00882 if (s%wetz(i,j) == 1) then 00883 hav(i,j) = max(par%eps,(s%hh(i,j) + s%runup(j))) 00884 elseif (hinterland == 0) then 00885 hav(i,j) = max(par%eps,s%runup(j)+s%zs(nint(s%iwl(j)),j)-s%zb(i,j) ) 00886 else 00887 hav(i,j) = par%eps; 00888 endif 00889 enddo 00890 enddo 00891 endif 00892 ! 00893 do i=2,s%nx !-1 Jaap -1 gives issues for bed updating at mpi boundaries 00894 do j=1,s%ny+1 00895 !if (max( max(hh(i,j),par%delta*H(i,j)), max(hh(i+1,j),par%delta*H(i+1,j)) )>par%hswitch+par%eps) then 00896 ip1=min(i+1,s%nx+1) 00897 jm1=max(j-1,1) 00898 totslp=sqrt(s%dzbdx(i,j)**2+(.25*(s%dzbdy(i,j)+s%dzbdy(ip1,j)+s%dzbdy(i,jm1)+s%dzbdy(ip1,jm1)))**2) 00899 if(max(hav(i,j),hav(i+1,j))>par%hswitch+par%eps) then ! Jaap instead of s%hh 00900 dzmax=par%wetslp*abs(s%dzbdx(i,j))/totslp 00901 ! tricks: seaward of istruct (transition from sand to structure) wetslope is set to 0.03; 00902 if (i>nint(s%istruct(j))) then 00903 !dzmax = 0.03d0 00904 dzmax = max(0.03d0,abs(s%dzbdx(i,j))*0.99d0) 00905 endif 00906 else 00907 dzmax=par%dryslp*abs(s%dzbdx(i,j))/totslp 00908 end if 00909 00910 if(abs(s%dzbdx(i,j))>dzmax+1.d-10 .and. s%structdepth(i+nint(max(0.d0,sign(1.d0,s%dzbdx(i,j)))),j)>par%eps) then 00911 n_aval=n_aval+1 00912 dzb=sign(1.0d0,s%dzbdx(i,j))*(abs(s%dzbdx(i,j))-dzmax)*s%dsu(i,j) 00913 00914 if (dzb >= 0.d0) then 00915 ie = i+1 ! index erosion point 00916 id = i ! index deposition point 00917 dAfac = s%dsdnzi(i,j)/s%dsdnzi(i+1,j) ! take into account varying grid sizes 00918 !dzb=min(dzb,par%dzmax*par%dt/s%dsu(i,j)) ! make sure dzb is not in conflict with par%dzmax 00919 dzb=dzb/par%avaltime*par%dt ! Dano more transparent formulation, removes direction dependence 00920 dzb=min(dzb,s%structdepth(i+1,j)) ! make sure dzb is not larger than sediment layer thickness 00921 else 00922 ie = i ! index erosion point 00923 id = i+1 ! index deposition point 00924 dAfac = s%dsdnzi(i+1,j)/s%dsdnzi(i,j) ! take into account varying grid sizes 00925 !dzb=max(dzb,-par%dzmax*par%dt/s%dsu(i,j)) 00926 dzb=dzb/par%avaltime*par%dt ! Dano more transparent formulation, removes direction dependence 00927 dzb=max(dzb,-s%structdepth(i,j)) 00928 endif 00929 00930 00931 if (par%ngd == 1) then ! Simple bed update in case one fraction 00932 00933 dzleft = abs(dzb) 00934 00935 s%zb(id,j) = s%zb(id,j)+dzleft*dAfac 00936 s%zb(ie,j) = s%zb(ie,j)-dzleft 00937 s%dzbnow(id,j) = s%dzbnow(id,j)+dzleft*dAfac ! naamgeveing? 00938 s%dzbnow(ie,j) = s%dzbnow(ie,j)-dzleft 00939 s%sedero(id,j) = s%sedero(id,j)+dzleft*dAfac 00940 s%sedero(ie,j) = s%sedero(ie,j)-dzleft 00941 s%structdepth(id,j) = max(0.d0,s%structdepth(id,j)+dzleft*dAfac) 00942 s%structdepth(ie,j) = max(0.d0,s%structdepth(ie,j)-dzleft) 00943 00944 s%zs(id,j) = s%zs(id,j)+dzleft*dAfac 00945 s%zs(ie,j) = s%zs(ie,j)-dzleft 00946 s%dzav(id,j)= s%dzav(id,j)+dzleft*dAfac 00947 s%dzav(ie,j)= s%dzav(ie,j)-dzleft 00948 00949 else ! multiple fractions... 00950 00951 ! now fix fractions.... 00952 dz => s%dzbed(ie,j,:) 00953 pb => s%pbbed(ie,j,:,:) 00954 00955 ! figure out how many depth layers (ndz) are eroded in point iii 00956 sdz = 0 00957 ndz = 0 00958 do while (sdz<abs(dzb)) 00959 ndz = ndz+1 00960 sdz = sdz+dz(ndz) 00961 enddo 00962 00963 ! now update bed and fractions by stepping through each layer seperately 00964 dzleft = abs(dzb) 00965 dzavt = 0.d0 00966 00967 do jdz=1,ndz 00968 00969 dzt = min(dz(jdz),dzleft) 00970 dzleft = dzleft-dzt; 00971 00972 ! erosion deposition per fraction (upwind or downwind); edg is positive in case of erosion 00973 do jg=1,par%ngd 00974 edg2(jg) = s%sedcal(jg)*dzt*pb(jdz,jg)*(1.d0-par%por)/par%dt ! erosion (dzt always > 0 ) 00975 edg1(jg) = -s%sedcal(jg)*edg2(jg)*dAfac ! deposition (dzt always < 0 ) 00976 enddo 00977 00978 dzavt = dzavt + sum(edg2)*par%dt/(1.d0-par%por) 00979 00980 call update_fractions(par,s,ie,j,s%dzbed(ie,j,:),s%pbbed(ie,j,:,:),edg2,dzavt) ! update bed in eroding point 00981 00982 call update_fractions(par,s,id,j,s%dzbed(id,j,:),s%pbbed(id,j,:,:),edg1,-dzavt*dAfac) ! update bed in deposition point 00983 00984 enddo 00985 00986 ! update water levels and dzav 00987 s%zs(ie,j) = s%zs(ie,j)-dzavt 00988 s%dzav(ie,j)= s%dzav(ie,j)-dzavt 00989 00990 s%zs(id,j) = s%zs(id,j)+dzavt*dAfac 00991 s%dzav(id,j)= s%dzav(id,j)+dzavt*dAfac 00992 00993 end if ! yes/no multiple fractions 00994 end if 00995 end do 00996 end do 00997 !Dano: computed earlier 00998 !do j=1,s%ny 00999 ! do i=1,s%nx+1 01000 ! s%dzbdy(i,j)=(s%zb(i,j+1)-s%zb(i,j))/s%dnv(i,j) 01001 ! enddo 01002 !enddo 01003 01004 do j=2,s%ny !-1 Jaap -1 gives issues for bed updating at mpi boundaries 01005 do i=1,s%nx+1 01006 jp1=min(j+1,s%ny+1) 01007 im1=max(i-1,1) 01008 totslp=sqrt(s%dzbdy(i,j)**2+(.25*(s%dzbdx(i,j)+s%dzbdx(i,jp1)+s%dzbdx(im1,j)+s%dzbdx(im1,jp1)))**2) 01009 01010 if(max(s%hh(i,j),s%hh(i,j+1))>par%hswitch+par%eps) then 01011 dzmax=par%wetslp*abs(s%dzbdy(i,j))/totslp 01012 else 01013 dzmax=par%dryslp*abs(s%dzbdy(i,j))/totslp 01014 end if 01015 if(abs(s%dzbdy(i,j))>dzmax .and. s%structdepth(i,j+nint(max(0.d0,sign(1.d0,s%dzbdy(i,j)))))>par%eps) then ! Jaap 01016 n_aval=n_aval+1 01017 dzb=sign(1.0d0,s%dzbdy(i,j))*(abs(s%dzbdy(i,j))-dzmax)*s%dnv(i,j) 01018 ! 01019 if (dzb >= 0.d0) then 01020 je = j+1 ! index erosion point 01021 jd = j ! index deposition point 01022 dAfac = s%dsdnzi(i,j)/s%dsdnzi(i,j+1) ! take into account varying grid sizes 01023 !dzb=min(dzb,par%dzmax*par%dt/s%dnv(i,j)) 01024 dzb=dzb/par%avaltime*par%dt ! Dano 01025 dzb=min(dzb,s%structdepth(i,j+1)) 01026 else 01027 je = j ! index erosion point 01028 jd = j+1 ! index deposition point 01029 dAfac = s%dsdnzi(i,j+1)/s%dsdnzi(i,j) ! take into account varying grid sizes 01030 !dzb=max(dzb,-par%dzmax*par%dt/s%dnv(i,j)) 01031 dzb=dzb/par%avaltime*par%dt ! Dano 01032 dzb=max(dzb,-s%structdepth(i,j)) 01033 endif 01034 01035 if (par%ngd == 1) then ! Simple bed update in case one fraction 01036 01037 dzleft = abs(dzb) 01038 01039 s%zb(i,jd) = s%zb(i,jd)+dzleft*dAfac 01040 s%zb(i,je) = s%zb(i,je)-dzleft 01041 s%dzbnow(i,jd) = s%dzbnow(i,jd)+dzleft*dAfac ! naamgeveing? 01042 s%dzbnow(i,je) = s%dzbnow(i,je)-dzleft 01043 s%sedero(i,jd) = s%sedero(i,jd)+dzleft*dAfac 01044 s%sedero(i,je) = s%sedero(i,je)-dzleft 01045 s%structdepth(i,jd) = max(0.d0,s%structdepth(i,jd)+dzleft*dAfac) 01046 s%structdepth(i,je) = max(0.d0,s%structdepth(i,je)-dzleft) 01047 01048 s%zs(i,jd) = s%zs(i,jd)+dzleft*dAfac 01049 s%zs(i,je) = s%zs(i,je)-dzleft 01050 s%dzav(i,jd)= s%dzav(i,jd)+dzleft*dAfac 01051 s%dzav(i,je)= s%dzav(i,je)-dzleft 01052 01053 else ! multiple fractions... 01054 01055 dz => s%dzbed(i,je,:) 01056 pb => s%pbbed(i,je,:,:) 01057 01058 ! figure out how many depth layers (ndz) are affected 01059 sdz = 0 01060 ndz = 0 01061 do while (sdz<abs(dzb)) 01062 ndz = ndz+1 01063 sdz = sdz+dz(ndz) 01064 enddo 01065 01066 ! now update bed and fractions by stepping through each layer seperately 01067 dzleft = abs(dzb) 01068 dzavt = 0.d0 01069 01070 do jdz=1,ndz 01071 dzt = min(dz(jdz),dzleft) 01072 dzleft = dzleft-dzt; 01073 01074 ! erosion deposition per fraction (upwind or downwind); edg is positive in case of erosion 01075 do jg=1,par%ngd 01076 edg2(jg) = s%sedcal(jg)*dzt*pb(jdz,jg)*(1.d0-par%por)/par%dt ! erosion (dzt always > 0 ) 01077 edg1(jg) = -s%sedcal(jg)*edg2(jg)*dAfac ! deposition (dzt always < 0 ) 01078 enddo 01079 01080 dzavt = dzavt + sum(edg2)*par%dt/(1.d0-par%por) 01081 01082 call update_fractions(par,s,i,je,s%dzbed(i,je,:),s%pbbed(i,je,:,:),edg2,dzavt) ! upwind point 01083 01084 call update_fractions(par,s,i,jd,s%dzbed(i,jd,:),s%pbbed(i,jd,:,:),edg1,-dzavt*dAfac) ! downwind point 01085 01086 enddo 01087 01088 ! update water levels and dzav 01089 s%zs(i,je) = s%zs(i,je)-dzavt 01090 s%dzav(i,je)= s%dzav(i,je)-dzavt 01091 01092 s%zs(i,jd) = s%zs(i,jd)+dzavt*dAfac 01093 s%dzav(i,jd)= s%dzav(i,jd)+dzavt*dAfac 01094 01095 endif !yes/no multiple fractions 01096 end if 01097 end do 01098 end do 01099 ! write(*,*)'t ',par%t,'ii ',ii,' n_aval ',n_aval,' rank ',xmpi_rank 01100 ! Dano: in parallel version bed update must take place AFTER EACH ITERATION 01101 #ifdef USEMPI 01102 call xmpi_allreduce(n_aval ,MPI_SUM) 01103 #endif 01104 if (n_aval>0) then 01105 #ifdef USEMPI 01106 call xmpi_shift_ee(s%zb) 01107 call xmpi_shift_ee(s%structdepth) 01108 #endif 01109 else 01110 exit 01111 endif 01112 01113 end do 01114 else 01115 s%dzbdx=0.d0 01116 s%dzbdy=0.d0 01117 do j=1,s%ny+1 01118 do i=1,s%nx 01119 s%dzbdx(i,j)=(s%zb(i+1,j)-s%zb(i,j))/s%dsu(i,j) 01120 enddo 01121 enddo 01122 do j=1,s%ny 01123 do i=1,s%nx+1 01124 s%dzbdy(i,j)=(s%zb(i,j+1)-s%zb(i,j))/s%dnv(i,j) 01125 enddo 01126 enddo 01127 end if 01128 ! 01129 ! bed boundary conditions 01130 ! 01131 ! Dano: the following Neumann boundaries introduce a mass error equal to 01132 ! the bed level change due to the b.c. times the cell area. 01133 if(xmpi_isleft .and. s%ny>0) then 01134 s%zb(:,1) = s%zb(:,2) 01135 s%dzbdt(:,1) = s%dzbdt(:,2) 01136 s%dzbnow(:,1) = s%dzbnow(:,2) 01137 s%sedero(:,1) = s%sedero(:,2) 01138 s%structdepth(:,1) = s%structdepth(:,2) 01139 s%pbbed(:,1,:,:)=s%pbbed(:,2,:,:) 01140 s%z0bed(:,1)=s%z0bed(:,2) 01141 s%dzbed(:,1,:)=s%dzbed(:,2,:) 01142 endif 01143 01144 if(xmpi_isright .and. s%ny>0) then 01145 s%zb(:,s%ny+1) = s%zb(:,s%ny) 01146 s%dzbdt(:,s%ny+1) = s%dzbdt(:,s%ny) 01147 s%dzbnow(:,s%ny+1) = s%dzbnow(:,s%ny) 01148 s%sedero(:,s%ny+1) = s%sedero(:,s%ny) 01149 s%structdepth(:,s%ny+1) = s%structdepth(:,s%ny) 01150 s%pbbed(:,s%ny+1,:,:)=s%pbbed(:,s%ny,:,:) 01151 s%z0bed(:,s%ny+1)=s%z0bed(:,s%ny) 01152 s%dzbed(:,s%ny+1,:)=s%dzbed(:,s%ny,:) 01153 endif 01154 01155 ! Update representative sed.diameter at the bed for flow friction and output 01156 if (par%ngd>1) then 01157 do j=j1,max(s%ny,1) 01158 do i=2,s%nx 01159 s%D50top = sum(s%pbbed(i,j,1,:)*s%D50) 01160 s%D90top = sum(s%pbbed(i,j,1,:)*s%D90) 01161 enddo 01162 enddo 01163 endif 01164 01165 endif ! if par%t>par%morstart 01166 01167 end subroutine bed_update 01168 01169 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01170 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01171 01172 subroutine update_fractions(par,s,i,j,dz,pb,edg,dzb) 01173 01174 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01175 ! Copyright (C) 2009 Technische Universiteit Delft 01176 ! Bram van Prooijen 01177 ! b.c.vanprooijen@tudelft.nl 01178 ! +31(0)15 2784070 01179 ! Faculty of Civil Engineering and Geosciences 01180 ! department of Hydraulic Engineering 01181 ! PO Box 5048 01182 ! 2600 GA Delft 01183 ! The Netherlands 01184 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01185 01186 use params, only: parameters 01187 use spaceparams 01188 use xmpi_module 01189 01190 01191 implicit none 01192 01193 type(spacepars),target :: s 01194 type(parameters) :: par 01195 01196 integer :: i,j,jg,jd 01197 real*8 :: ED,zbold,dzbt,fac,dzb_loc 01198 real*8, intent(in) :: dzb 01199 real*8 , dimension(par%ngd),intent(in) :: edg 01200 real*8 , dimension(s%nd(i,j)),intent(inout) :: dz 01201 real*8 , dimension(s%nd(i,j),par%ngd),intent(inout) :: pb 01202 01203 real*8 , dimension(:),allocatable,save :: Ap,b 01204 real*8 , dimension(:,:),allocatable,save :: Sm,A 01205 01206 if (.not. allocated(Ap)) then 01207 allocate(Ap (s%nd(1,1))) 01208 allocate(b (s%nd(1,1))) 01209 allocate(Sm (s%nd(1,1),par%ngd)) 01210 allocate(A (s%nd(1,1),3)) 01211 endif 01212 !TODO, dzb_loc is not initialized can be Nan, leading to infinite loop 01213 dzb_loc = dzb 01214 01215 !!!initialize Sm 01216 01217 ED = sum(edg) 01218 ! do t_sub=1,nt_sub !loop over subtimesteps 01219 do while (abs(dzb_loc) .gt. 0.d0) 01220 dzbt = min(dzb_loc,dz(par%nd_var)) ! make sure erosion (dzg is positive) is limited to thickness of variable layer 01221 dzbt = max(dzbt,-par%frac_dz*dz(par%nd_var+1)) ! make sure deposition (dzg is negative) is limited to thickness of first layer below variable layer 01222 01223 fac = dzbt/dzb_loc ! factor over mass change in cells to limit erosion and deposition 01224 01225 dzb_loc = dzb_loc-dzbt ! update dzb 01226 01227 do jg=1,par%ngd 01228 01229 A=0. 01230 Ap=0. 01231 b=0. 01232 01233 Sm(:,jg)=pb(:,jg)*dz*(1.d0-par%por) 01234 01235 !!!build matrix A 01236 select case(par%nd_var) 01237 case(1) 01238 !in this case: A=0 01239 case(2) 01240 A (1,1:3)= (/0.d0 , min(ED,0.d0) , max(ED,0.d0) /) 01241 A (2,1:3)= (/-min(ED,0.d0) , -max(ED,0.d0) , 0.d0 /) 01242 case(3:1000) 01243 A (1,1:3)= (/0.d0 , min(ED,0.d0) , max(ED,0.d0) /) 01244 01245 A (2:par%nd_var-1,1)=-min(ED,0.d0) 01246 A (2:par%nd_var-1,2)=-abs(ED) 01247 A (2:par%nd_var-1,3)=max(ED,0.d0) 01248 01249 A (par%nd_var,1:3)= (/-min(ED,0.d0) , -max(ED,0.d0) , 0.d0 /) 01250 end select 01251 01252 !!!determine RHS 01253 01254 ! Ap makes sure that layer nd_var varies in thickness in time 01255 ! Ap = 0 with single fraction (in case(1)) 01256 Ap(1) = sum(A(1,2:3)*pb(1:2,jg)) 01257 do jd = 2,par%nd_var 01258 Ap(jd) = sum(A (jd,:)*pb(jd-1:jd+1,jg)) 01259 enddo 01260 Ap(par%nd_var+1) = sum(A(par%nd_var+1,1:2)*pb(par%nd_var:par%nd_var+1,jg)) 01261 01262 ! b represents the actual erosion and deposition in the top layer. 01263 ! However, the thickness of the top layer remains constant and instead layer nd_var will breath in thickness 01264 b(1) = -edg(jg) 01265 01266 !!!update Sm 01267 01268 ! Sm is the sediment mass per fraction per layer 01269 ! Sm(1:par%nd_var+1,jg) = Sm(1:par%nd_var+1,jg) + par%dt/nt_sub*(Ap(1:par%nd_var+1)+b(1:par%nd_var+1)) 01270 Sm(1:par%nd_var+1,jg) = Sm(1:par%nd_var+1,jg) + par%dt*fac*(Ap(1:par%nd_var+1)+b(1:par%nd_var+1)) 01271 01272 enddo !fractions 01273 01274 ! From Sm we can compute the new fraction ratios per layer and the layer thickness... 01275 do jd=1,par%nd_var+1 01276 if (sum(Sm(jd,:))>0) then 01277 pb(jd,:) = Sm(jd,:)/sum(Sm(jd,:)) 01278 else 01279 pb(jd,:) = pb(jd,:) 01280 endif 01281 dz(jd) = sum(Sm(jd,:))/(1.d0-par%por) 01282 enddo 01283 01284 !!! modify grid 01285 01286 !merge two upper layers in case of erosion 01287 if (dz(par%nd_var) .lt. par%merge*dz(par%nd_var+1)) then 01288 forall (jg=1:par%ngd) 01289 pb(par%nd_var,jg) = (dz(par%nd_var)*pb(par%nd_var,jg) + dz(par%nd_var+1)* & 01290 pb(par%nd_var+1,jg))/(dz(par%nd_var)+dz(par%nd_var+1)) 01291 pb(par%nd_var+1:s%nd(i,j)-1,jg) = pb(par%nd_var+2:s%nd(i,j),jg) 01292 pb(s%nd(i,j),jg) = pb(s%nd(i,j),jg) 01293 endforall 01294 s%z0bed(i,j) = s%z0bed(i,j)-dz(par%nd_var+1) 01295 dz(par%nd_var) = dz(par%nd_var+1)+dz(par%nd_var) 01296 endif 01297 !split upper layer in case of sedimentation 01298 if (dz(par%nd_var)>par%split*dz(par%nd_var+1)) then 01299 pb(par%nd_var+1:s%nd(i,j),:) = pb(par%nd_var:s%nd(i,j)-1,:) 01300 s%z0bed(i,j) = s%z0bed(i,j)+dz(par%nd_var+1) 01301 dz(par%nd_var) = dz(par%nd_var)-dz(par%nd_var+1) 01302 endif 01303 enddo ! nt_sub 01304 01305 pb = max(0.d0,min(pb,1.d0)) 01306 01307 zbold = s%zb(i,j) 01308 s%zb(i,j) = s%z0bed(i,j)+sum(dz) 01309 s%dzbnow(i,j) = s%dzbnow(i,j)+(s%zb(i,j)-zbold) 01310 s%sedero(i,j) = s%sedero(i,j)+(s%zb(i,j)-zbold) 01311 s%structdepth(i,j) = max(0.d0,s%structdepth(i,j)+(s%zb(i,j)-zbold)) 01312 01313 end subroutine update_fractions 01314 01315 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01316 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01317 01318 subroutine sedtransform(s,par) 01319 use params, only: parameters 01320 use spaceparams 01321 use readkey_module 01322 use xmpi_module 01323 use paramsconst 01324 01325 implicit none 01326 01327 type(spacepars),target :: s 01328 type(parameters) :: par 01329 01330 integer :: i,j,jg, ii 01331 real*8 :: z0,dcf,dcfin,ML 01332 real*8 :: Te,Sster,cc1,cc2,wster,Ass 01333 real*8 :: kl,alpha,alpha1,alpha2,beta,psi 01334 real*8 , save :: delta,kvis,onethird,twothird,phi 01335 real*8 , dimension(:),allocatable ,save :: dster,ws0, shieldscrit, sigz, ceqssteps, hhsteps 01336 real*8 , dimension(:,:),allocatable ,save :: vmg,Cd,Asb,dhdx,dhdy,Ts,hfac 01337 real*8 , dimension(:,:),allocatable ,save :: urms2,Ucr,Ucrc,Ucrw,term1,B2,srfTotal,srfRhee,vero,Ucrb,Ucrs 01338 real*8 , dimension(:,:),allocatable ,save :: uandv,b,fslope,hloc,ceqs,ceqb,fallvelredfac 01339 real*8 , dimension(:,:,:),allocatable,save :: w 01340 real*8 , dimension(:,:),allocatable ,save :: uorb, A, ksw, fw, uw, tauwav, muw, fc, f1c, tauc, muc, taubcw, taucr, used, ue 01341 real*8 :: ca, refA, E, Me, M, Eswm, Eswb, Esw, Esc, uster, deltar, delw, ta, Escm, ceqexplicit, ceqimplicit, hhtmp, hhtmp2, dhhtmp, sigs, ceq_tmp 01342 01343 !include 's.ind' 01344 !include 's.inp' 01345 01346 if (.not. allocated(vmg)) then 01347 allocate (vmg (s%nx+1,s%ny+1)) 01348 allocate (term1 (s%nx+1,s%ny+1)) 01349 allocate (B2 (s%nx+1,s%ny+1)) 01350 allocate (Cd (s%nx+1,s%ny+1)) 01351 allocate (Asb (s%nx+1,s%ny+1)) 01352 allocate (Ucr (s%nx+1,s%ny+1)) 01353 allocate (Ucrc (s%nx+1,s%ny+1)) 01354 allocate (Ucrw (s%nx+1,s%ny+1)) 01355 allocate (urms2 (s%nx+1,s%ny+1)) 01356 allocate (hloc (s%nx+1,s%ny+1)) 01357 allocate (Ts (s%nx+1,s%ny+1)) 01358 allocate (ceqs (s%nx+1,s%ny+1)) 01359 allocate (ceqb (s%nx+1,s%ny+1)) 01360 allocate (srfTotal(s%nx+1,s%ny+1)) ! Lodewijk 01361 allocate (Ucrb (s%nx+1,s%ny+1)) ! Lodewijk 01362 allocate (Ucrs (s%nx+1,s%ny+1)) ! Lodewijk 01363 allocate (srfRhee(s%nx+1,s%ny+1)) ! Lodewijk 01364 allocate (vero (s%nx+1,s%ny+1)) ! Lodewijk 01365 allocate (fallvelredfac(s%nx+1,s%ny+1)) ! Lodewijk 01366 allocate (w (s%nx+1,s%ny+1,par%ngd)) ! Lodewijk 01367 allocate (dster (par%ngd)) 01368 allocate (ws0 (par%ngd)) 01369 allocate (shieldscrit(par%ngd)) 01370 allocate (dhdx (s%nx+1,s%ny+1)) 01371 allocate (dhdy (s%nx+1,s%ny+1)) 01372 allocate (uandv (s%nx+1,s%ny+1)) 01373 allocate (b (s%nx+1,s%ny+1)) 01374 allocate (fslope(s%nx+1,s%ny+1)) 01375 allocate (hfac (s%nx+1,s%ny+1)) 01376 ! Van Rijn (1993) by Kees 01377 allocate (used (s%nx+1,s%ny+1)) 01378 allocate (ue (s%nx+1,s%ny+1)) 01379 allocate (uorb (s%nx+1,s%ny+1)) 01380 allocate (A (s%nx+1,s%ny+1)) 01381 allocate (ksw (s%nx+1,s%ny+1)) 01382 allocate (fw (s%nx+1,s%ny+1)) 01383 allocate (uw(s%nx+1,s%ny+1)) 01384 allocate (tauwav (s%nx+1,s%ny+1)) 01385 allocate (muw (s%nx+1,s%ny+1)) 01386 allocate (fc (s%nx+1,s%ny+1)) 01387 allocate (f1c (s%nx+1,s%ny+1)) 01388 allocate (muc (s%nx+1,s%ny+1)) 01389 allocate (tauc (s%nx+1,s%ny+1)) 01390 allocate (taubcw(s%nx+1,s%ny+1)) 01391 allocate (taucr (s%nx+1,s%ny+1)) 01392 allocate(sigz(par%kmax)) 01393 allocate(hhsteps(par%kmax+1)) 01394 allocate(ceqssteps(par%kmax+1)) 01395 01396 vmg = 0.d0 01397 onethird=1.d0/3.d0 01398 twothird=2.d0/3.d0 01399 phi = par%reposeangle/180*par%px ! Angle of internal friction 01400 ! Robert: do only once, not necessary every time step 01401 do jg=1,par%ngd 01402 ! cjaap: compute fall velocity with simple expression from Ahrens (2000) 01403 Te = 20.d0 01404 kvis = 4.d0/(20.d0+Te)*1d-5 ! Van rijn, 1993 01405 Sster = s%D50(jg)/(4*kvis)*sqrt((par%rhos/par%rho-1)*par%g*s%D50(jg)) 01406 cc1 = 1.06d0*tanh(0.064d0*Sster*exp(-7.5d0/Sster**2)) 01407 cc2 = 0.22d0*tanh(2.34d0*Sster**(-1.18d0)*exp(-0.0064d0*Sster**2)) 01408 wster = cc1+cc2*Sster 01409 ws0(jg) = wster*sqrt((par%rhos/par%rho-1.d0)*par%g*s%D50(jg)) 01410 ! RJ: for modeling gravel 01411 delta = (par%rhos-par%rho)/par%rho 01412 dster(jg)=(delta*par%g/1.d-12)**onethird*s%D50(jg) 01413 shieldscrit(jg) = 0.3d0/(1+1.2d0*dster(jg))+0.055d0*(1-exp(-0.02d0*dster(jg))) ! Soulsby, 1997 (eq 76) 01414 if (par%fallvelred==0) then 01415 w(:,:,jg) = ws0(jg) 01416 endif 01417 enddo 01418 endif 01419 01420 01421 ! Lodewijk: calculate fall velocity reduction coefficient based on concentration of previous time step 01422 ! Rowe (1987) made an estimation of the exponent alpha by fitting a logarithmic function on a dataset of Richardson and Zaki (1954). 01423 if (par%fallvelred==1) then 01424 do jg=1,par%ngd 01425 alpha = 2.35d0*(2.d0+0.175d0*(ws0(jg)*s%D50(jg)/kvis)**0.75d0)/(1.d0+0.175d0*(ws0(jg)*s%D50(jg)/kvis)**0.75d0) 01426 fallvelredfac = (1.d0-s%ccg(:,:,jg))**(alpha) 01427 w(:,:,jg) = ws0(jg)*fallvelredfac 01428 enddo 01429 endif 01430 ! 01431 ! hloc = max(hh,0.01d0) ! Jaap 01432 hloc = max(s%hh,0.01) 01433 ! Compute mean fall velocity 01434 par%ws=0.d0 01435 do jg=1,par%ngd 01436 par%ws=par%ws+w(1,1,jg) 01437 enddo 01438 par%ws=par%ws/par%ngd 01439 01440 ! 01441 ! compute near bed turbulence 01442 ! 01443 ! due to short waves 01444 01445 if (par%swave==1) then 01446 ! wave breaking induced turbulence due to short waves 01447 do j=1,s%ny+1 01448 do i=1,s%nx+1 01449 ! compute mixing length 01450 ! ML = 2*R(i,j)*par%Trep/(par%rho*c(i,j)*max(H(i,j),par%eps)) 01451 ML = dsqrt(2*s%R(i,j)*par%Trep/(par%rho*max(s%c(i,j),1d-10))) 01452 ! ML = 0.9d0*H(i,j) 01453 ML = min(ML,hloc(i,j)); 01454 ! exponential decay turbulence over depth 01455 dcfin = exp(min(100.d0,hloc(i,j)/max(ML,0.00000000000000001d0))) 01456 dcf = min(1.d0,1.d0/(dcfin-1.d0)) 01457 ! Jaap: new approach: compute kb based on waveturb result 01458 s%kb(i,j) = s%kturb(i,j)*dcf 01459 ! 01460 if (par%turb == TURB_BORE_AVERAGED) then 01461 s%kb(i,j) = s%kb(i,j)*par%Trep/s%Tbore(i,j) 01462 endif 01463 01464 enddo 01465 ! Jaap: rundown jet creating additional turbulence 01466 if (par%swrunup==1)then 01467 s%kb(nint(s%istruct(j)),j) = s%kb(nint(s%istruct(j)),j) + par%jetfac* & 01468 (s%E(nint(s%istruct(j)),j)*s%strucslope(j)*sqrt(par%g/s%hh(nint(s%istruct(j)),j)))**twothird 01469 endif 01470 enddo 01471 elseif (par%wavemodel==WAVEMODEL_NONH) then 01472 do j=1,s%ny+1 01473 do i=1,s%nx+1 01474 !ML=max(s%rolthick(i,j),0.07d0*max(s%hh(i,j),par%hmin)) 01475 !s%kb(i,j) = s%kturb(i,j)*ML/max(s%hh(i,j),par%hmin) ! Simpler expression 01476 s%kb(i,j) = s%kturb(i,j) ! even simpler expression :-) 01477 enddo 01478 enddo 01479 endif !par%swave == 1 01480 01481 ! switch to include long wave stirring 01482 if (par%lws==1) then 01483 vmg = dsqrt(s%ue**2+s%ve**2) 01484 elseif (par%lws==0) then 01485 ! vmg lags on actual mean flow; but long wave contribution to mean flow is included... 01486 vmg = (1.d0-1.d0/par%cats/par%Trep*par%dt)*vmg + (1.d0/par%cats/par%Trep*par%dt)*dsqrt(s%ue**2+s%ve**2) 01487 endif 01488 01489 urms2 = s%urms**2.d0+1.45d0*s%kb 01490 01491 do jg = 1,par%ngd 01492 01493 Ts = par%tsfac*hloc/w(:,:,jg) 01494 if (par%oldTsmin==1) then 01495 s%Tsg(:,:,jg) = max(Ts,par%Tsmin) 01496 else 01497 s%Tsg(:,:,jg) = max(Ts,par%dtlimTs*par%dt) 01498 endif 01499 ! 01500 ! calculate treshold velocity Ucr 01501 ! 01502 if (par%form==FORM_SOULSBY_VANRIJN) then ! Soulsby van Rijn 01503 if(s%D50(jg)<=0.0005d0) then 01504 Ucr=0.19d0*s%D50(jg)**0.1d0*log10(4.d0*hloc/s%D90(jg)) 01505 else !Dano see what happens with coarse material 01506 Ucr=8.5d0*s%D50(jg)**0.6d0*log10(4.d0*hloc/s%D90(jg)) 01507 end if 01508 elseif ((par%form==FORM_VANTHIEL_VANRIJN) .or. (par%form==FORM_VANRIJN1993)) then 01509 ! needed for full Van Thiel de Vries & Reniers 2008 or bed load Van Rijn 1993 01510 if(s%D50(jg)<=0.0005) then 01511 Ucrc=0.19d0*s%D50(jg)**0.1d0*log10(4.d0*hloc/s%D90(jg)) !Shields 01512 Ucrw=0.24d0*(delta*par%g)**0.66d0*s%D50(jg)**0.33d0*par%Trep**0.33d0 !Komar and Miller (1975) 01513 else if(s%D50(jg)<=0.002) then 01514 Ucrc=8.5d0*s%D50(jg)**0.6d0*log10(4.d0*hloc/s%D90(jg)) !Shields 01515 Ucrw=0.95d0*(delta*par%g)**0.57d0*s%D50(jg)**0.43*par%Trep**0.14 !Komar and Miller (1975) 01516 else if(s%D50(jg)>0.002) then 01517 Ucrc=1.3d0*sqrt(delta*par%g*s%D50(jg))*(hloc/s%D50(jg))**(0.5d0*onethird) !Maynord (1978) --> also Neill (1968) where 1.3d0 = 1.4d0 01518 Ucrw=0.95d0*(delta*par%g)**0.57d0*s%D50(jg)**0.43*par%Trep**0.14 !Komar and Miller (1975) 01519 end if 01520 B2 = vmg/max(vmg+dsqrt(urms2),par%eps) 01521 Ucr = B2*Ucrc + (1-B2)*Ucrw !Van Rijn 2007 (Bed load transport paper) 01522 end if 01523 01524 ! Lodewijk: implementation Rhee (2010), reduction sediment transport due dilatancy 01525 srfRhee(:,:) = 0.d0 01526 srfTotal(:,:) = 1.d0 01527 if (par%dilatancy == 1) then 01528 vero(:,:) = max(0.d0,-1.d0*s%dzbdt) ! Erosion velocity, for now asume it equal to -s%dzbdt of the previous time step 01529 kl = par%g/(160.d0*kvis)*(s%D15(jg)**2.d0)*((par%por**3.d0)/(1.d0-par%por)**2.d0) ! Permeability, Adel 1987, which is based on the Ergun equation 01530 ! bed porosity, estimated here as maximum porosity, to be added to user input 01531 ! A=3/4 for single particles and A=1/(1-n0) for a continuum 01532 ! Reduction factor on the critical Shields parameter by dilatancy (Van Rhee, 2010) 01533 srfRhee(:,:) = vero(:,:)/kl*(par%pormax-par%por)/(1.d0-par%pormax)*par%rheeA/delta 01534 endif 01535 ! Lodewijk: implementation bed slope reduction on critical flow velocity 01536 if (par%bdslpeffini == BDSLPEFFINI_NONE) then 01537 srfTotal(:,:) = 1.d0 + srfRhee(:,:) 01538 elseif (par%bdslpeffini == BDSLPEFFINI_TOTAL .or. par%bdslpeffini == BDSLPEFFINI_BED) then 01539 do j=1,s%ny+1 01540 do i=1,s%nx+1 01541 ! Prevent NaN values if too small values 01542 if ((dabs(s%ue(i,j)) > 0.000001d0 .or. dabs(s%ve(i,j)) > 0.000001d0) .and. & 01543 (dabs(s%dzbdx(i,j))>0.000001d0 .or. dabs(s%dzbdy(i,j))>0.000001d0)) then 01544 ! Angle between the x-axis and the flow velocity 01545 ! REMARK: also include waves in the velocity? 01546 if (s%ue(i,j) < 0.d0) then 01547 alpha1 = datan(s%ve(i,j)/s%ue(i,j)) + par%px 01548 else 01549 alpha1 = datan(s%ve(i,j)/s%ue(i,j)) 01550 endif 01551 ! Angle between the x-axis and the bed slope vector directed in down-slope direction, derived in thesis Lodewijk 01552 if (s%dzbdy(i,j) >= 0.d0) then 01553 alpha2 = -datan(s%dzbdx(i,j)/s%dzbdy(i,j))+1.5d0*par%px 01554 else 01555 alpha2 = -datan(s%dzbdx(i,j)/s%dzbdy(i,j))+0.5d0*par%px 01556 endif 01557 psi = alpha1-(alpha2-par%px) ! Angle between the flow direction and the up-slope directed vector 01558 if (dabs(s%dzbdx(i,j))<0.000001d0) then ! A smaller slope could result in a NaN for beta 01559 ! Beta purely based on dzbdy 01560 beta = datan(dabs(s%dzbdy(i,j))) 01561 else 01562 beta = datan(dabs(s%dzbdx(i,j)/dsin(datan(s%dzbdx(i,j)/s%dzbdy(i,j))))) ! Maximum absolute bed slope angle, derived in thesis Lodewijk 01563 endif 01564 beta = min(beta,phi) ! Take min to exclude NaN's 01565 if (par%dilatancy == 1) then 01566 srfTotal(i,j) = (dcos(psi)*dsin(beta)+dsqrt( & 01567 ((srfRhee(i,j))**2+2*srfRhee(i,j)*dcos(beta)+dcos(beta)**2 & 01568 )*dtan(phi)**2-dsin(psi)**2*dsin(beta)**2 & 01569 ))/dtan(phi) ! Soulsby (1997), modified by Lodewijk (see Thesis) 01570 else 01571 srfTotal(i,j) = (dcos(psi)*dsin(beta)+ & 01572 dsqrt(dcos(beta)**2*dtan(phi)**2-dsin(psi)**2*dsin(beta)**2))/dtan(phi) ! Soulsby (1997) 01573 endif 01574 endif 01575 end do 01576 end do 01577 endif 01578 ! Calculate the new critical velocity based on the modification factors on the Shields parameter, bed slope only on bed load 01579 Ucrb(:,:) = Ucr(:,:)*sqrt(srfTotal) ! Lodewijk 01580 if (par%bdslpeffini == BDSLPEFFINI_TOTAL) then 01581 Ucrs = Ucrb 01582 else 01583 Ucrs = Ucr*(1.d0+sqrt(srfRhee)) ! Lodewijk, no effect on suspended load by bed slope 01584 endif 01585 01586 01587 if (par%form==FORM_SOULSBY_VANRIJN) then ! Soulsby van Rijn 01588 ! 01589 ! drag coefficient 01590 z0 = par%z0 01591 Cd=(0.40d0/(log(max(hloc,10.d0*z0)/z0)-1.0d0))**2 01592 ! 01593 ! transport parameters 01594 Asb=0.005d0*hloc*(s%D50(jg)/hloc/(delta*par%g*s%D50(jg)))**1.2d0 ! bed load coefficent 01595 Ass=0.012d0*s%D50(jg)*dster(jg)**(-0.6d0)/(delta*par%g*s%D50(jg))**1.2d0 ! suspended load coeffient 01596 ! 01597 term1=(vmg**2+0.018/Cd*par%sws*urms2) ! Make 0.018/Cd is always smaller than the flow friction coefficient 01598 ! 01599 ! reduce sediment suspensions for (inundation) overwash conditions with critical flow velocitties 01600 ! vmag2=min(vmag2,par%smax*par%C**2*D50(jg)*delta) 01601 ! vmag2=min(vmag2,par%smax*par%g/par%cf*D50(jg)*delta) ! In terms of cf 01602 ! term1=sqrt(vmag2+0.018d0/Cd*urms2) ! nearbed-velocity 01603 ! the two above lines are comment out and replaced by a limit on total velocity u2+urms2, robert 1/9 and ap 28/11 01604 ! 01605 term1=min(term1,par%smax*par%g/s%cf*s%D50(jg)*delta) 01606 term1=sqrt(term1) 01607 ! 01608 ceqb = 0.d0 !initialize ceqb 01609 ceqs = 0.d0 !initialize ceqs 01610 do j=1,s%ny+1 01611 do i=1,s%nx 01612 ! Lodewijk: separate bed load from suspended load since Ucr not anymore the same 01613 if(term1(i,j)>Ucrb(i,j) .and. hloc(i,j)>par%eps) then 01614 ceqb(i,j)=Asb(i,j)*(term1(i,j)-Ucrb(i,j))**2.4d0 01615 end if 01616 if(term1(i,j)>Ucrs(i,j) .and. hloc(i,j)>par%eps) then 01617 ceqs(i,j)=Ass*(term1(i,j)-Ucrs(i,j))**2.4d0 01618 end if 01619 end do 01620 end do 01621 ! 01622 elseif (par%form==FORM_VANTHIEL_VANRIJN) then ! Van Thiel de Vries & Reniers 2008 01623 ! 01624 ! transport parameters 01625 Asb=0.015d0*hloc*(s%D50(jg)/hloc)**1.2d0/(delta*par%g*s%D50(jg))**0.75d0 !bed load coefficent 01626 Ass=0.012d0*s%D50(jg)*dster(jg)**(-0.6d0)/(delta*par%g*s%D50(jg))**1.2d0 !suspended load coeffient 01627 ! 01628 ! Jaap: par%sws to set short wave stirring to zero 01629 ! Jaap: Van Rijn use Peak orbital flow velocity --> 0.64 corresponds to 0.4 coefficient regular waves Van Rijn (2007) 01630 term1=vmg**2+0.64d0*par%sws*urms2 01631 ! reduce sediment suspensions for (inundation) overwash conditions with critical flow velocitties 01632 term1=min(term1,par%smax*par%g/max(s%cf,1d-10)*s%D50(jg)*delta) 01633 term1=sqrt(term1) 01634 ! 01635 ceqb = 0.d0 !initialize ceqb 01636 ceqs = 0.d0 !initialize ceqs 01637 do j=1,s%ny+1 01638 do i=1,s%nx 01639 ! Lodewijk: sepperate bed load from suspended load since Ucr not anymore the same 01640 if(term1(i,j)>Ucrb(i,j) .and. hloc(i,j)>par%eps) then 01641 ceqb(i,j)=Asb(i,j)*(term1(i,j)-Ucrb(i,j))**1.5 01642 end if 01643 if(term1(i,j)>Ucrs(i,j) .and. hloc(i,j)>par%eps) then 01644 ceqs(i,j)=Ass*(term1(i,j)-Ucrs(i,j))**2.4 01645 end if 01646 end do 01647 end do 01648 ! 01649 elseif (par%form==FORM_VANRIJN1993) then ! Van Rijn (1993) a la Delft3D 01650 ! 01651 ! General parameters 01652 ue = dsqrt(vmg**2) ! NOT similar to Van Thiel -> velocity magnitude for stirring sediment 01653 used = vmg + s%ua ! velocity used for sediment transport (reduced with asymmetry) 01654 ! 01655 ! Bed load transport (Delft3D-FLOW manual; eq 11.82) 01656 ! 01657 do j=1,s%ny+1 01658 do i=1,s%nx+1 01659 if (ue(i,j)>Ucrb(i,j) .and. hloc(i,j)>par%eps) then ! Ucrs is weighted current (including waves) 01660 Me = (used(i,j)-Ucrb(i,j))**2 / ((delta*par%g*s%D50(jg))) ! Effective mobility parameter 01661 M = (used(i,j)-0)**2 / ((delta*par%g*s%D50(jg))) ! "normal" mobility parameter 01662 Asb(i,j) = 0.006d0 * s%D50(jg) * ws0(jg) ! bed load coefficent (without ps) 01663 ceqb(i,j) = Asb(i,j) * (M**0.5)* (Me**0.7) ! calculate the main part 01664 ceqb(i,j) = ceqb(i,j) / used(i,j) ! goal is sediment concentration: divide by velocity 01665 else 01666 ceqb(i,j) = 0 01667 end if 01668 end do 01669 end do 01670 ! 01671 ! Suspended transport (Van Rijn, 1993: Delft3D-FLOW manual) 01672 ! 01673 deltar = par%deltar ! ripple height input (eq. 11.78) 01674 ksw = par%rwave*deltar ! wave-related roughness -> based on ripple height (eq. 11.78) 01675 ksw = max(min(0.1, ksw), 0.01) ! numerical limits: not larger or smaller 01676 ! orbital velocityand peak orbital excursion 01677 uorb = par%px*s%H/par%Trep/sinh(min(max(s%k,0.01d0)*max(s%hh,par%delta*s%H),10.0d0)) ! orbital velocity (linear wave theory) 01678 A = uorb*par%Trep / (2.0*par%px) ! peak orbital excursion (below eq. 11.80) 01679 ! numerical limits for A, otherwise infinity problems 01680 where (A<0.1d0) 01681 A = 0.1d0 01682 endwhere 01683 ! 01684 ! Calculate bed-shear stress due to waves 01685 fw = min(exp( - 6.0 + 5.2*(A/ksw)**( - 0.19)), 0.3) ! Van Rijn & Kroon (1992); eq 13 below 01686 muw = 0.6/dster(jg) ! efficiency waves 01687 tauwav = 0.25 * par%rho * fw * uorb**2 ! bed shear stress due to waves (eq. 11.76) 01688 ! 01689 ! Calculate bed-shear stress due to currents 01690 f1c = 0.24*log10(12.0*s%hh/(3.0*3*s%D90(jg)))**( - 2) ! total-current relation factor (eq. 11.73) 01691 fc = 0.24*log10(12.0*s%hh/deltar)**( - 2) ! grain-related factor (eq. 11.72) 01692 muc = f1c/fc ! efficiency current 01693 tauc = 0.125d0* par%rho * fc * (s%ue**2+s%ve**2) ! Van Rijn & Kroon (1992); eq. 13 01694 ! 01695 ! Calculate bed shear stress ratio: combined waves and currents 01696 taubcw = muc*tauc + muw*tauwav ! combined shear stress (eq. 11.70) 01697 taucr = shieldscrit(jg) * par%g * (par%rhos-par%rho) * s%D50(jg) ! critical bed shear stress (Soulsby, 1997 at Eq. 74; page 104) 01698 ! 01699 ! Loop over every grid cell where critical bed shear stress < combined shear stress due to waves and currents 01700 ! clean variables 01701 s%ccz(:,:,:) = 0 01702 s%refA(:,:) = 0 01703 s%ca(:,:) = 0 01704 do j=1,s%ny+1 01705 do i=1,s%nx+1 01706 if (taubcw(i,j)>taucr(i,j) .and. hloc(i,j)>par%eps) then 01707 ! 01708 ! A) Reference level 01709 refA = 0.5*deltar ! reference level is half the ripple height 01710 refA = max(0.01, refA) ! minimum reference level 01711 refA = min(0.2*s%hh(i,j), refA) ! maximum reference level 20% of the water depth -> same Delft3D 01712 ! 01713 ! B) Compute wave boundary laver thickness 01714 delw = 0.072*A(i,j)*(A(i,j)/ksw(i,j))**(-0.25) 01715 ! 01716 ! C) Calculate Van Rijn's dimensionless bed-shear stress for reference concentration at z=a 01717 ta = max(0.0, (muc(i,j)*tauc(i,j) + muw(i,j)*tauwav(i,j))/taucr(i,j) - 1.0) 01718 ! 01719 ! D) Calculate Van Rijn's reference concentration 01720 ca = 0.015*s%d50(jg)*ta**1.5/(refA*dster(jg)**0.3) ! Van Rijn, 1984; no rho in there) 01721 ca = min(ca, 0.05) ! maximum ca value 01722 ! 01723 ! E) Depth-averaged mixing due to waves 01724 ! 1) Van Rijn & Kroon (1992); eq. 9 (maximum) 01725 Eswm = 0.035 * s%H(i,j)*s%hh(i,j) / par%Trep 01726 ! 2) Van Rijn & Kroon (1992); eq. 8 (at the bed) with 3 deltar = thicknes layer 01727 Eswb = 0.004*dster(jg)*uorb(i,j)*(3.0*delw) 01728 ! 01729 ! F) Depth-averaged mixing due to currents 01730 uster = ((s%taubx(i,j)**2 + s%tauby(i,j)**2)**0.5 / par%rho)**0.5 ! calculate bed shear velocity 01731 Escm = 0.25 * par%vonkar * uster * s%hh(i,j) ! Van Rijn & Kroon (1992); eq. 10 (maximum) 01732 ! 01733 ! G) Numerically integrate the concentration profile (no analyical expression with varying mixing exists!) 01734 ! one could assume a constant mixing coefficient: in that case exponential profile with analyical epxression 01735 ! 1) Create the grid on which we approximate the concentration profile 01736 ! Grid goes from reference level to water depth 01737 hhsteps(1) = refA 01738 do ii=1,(par%kmax) 01739 sigs = 0.01 ! starting percentile; 1% now 01740 sigz(ii) = sigs*(1/sigs)**(float(ii-1)/float(par%kmax-1)) ! could place outside the loop for efficiency 01741 hhsteps(ii+1) = refA + (s%hh(i,j)-refA) * sigz(ii) 01742 enddo 01743 ! 01744 ! 2) Now numerical integrate 01745 ! Start concentration is the reference concentration of Van Rijn, 1984 01746 ceqssteps(1) = ca 01747 do ii = 1, (par%kmax) 01748 ! 01749 ! i) Determine water depth 01750 hhtmp = hhsteps(ii) ! current water depth 01751 hhtmp2 = hhsteps(ii+1) ! water depth in next step 01752 dhhtmp = hhtmp2-hhtmp ! change in water depth in this step -> related to sigma layer 01753 ! 01754 ! ii) Determine e,waves 01755 if (hhtmp < (3*delw)) then 01756 Esw = Eswb 01757 elseif (hhtmp > (0.5*s%hh(i,j))) then 01758 Esw = Eswm 01759 else 01760 Esw = Eswb + (Eswm-Eswb)* (hhtmp - (3*delw)) / (0.5*s%hh(i,j) - 3*delw) 01761 endif 01762 ! 01763 ! iii) Determine e,currents 01764 if (hhtmp > (0.5*s%hh(i,j))) then 01765 Esc = Escm ! Van Rijn & Kroon (1992); eq. 10 (maximum) 01766 else 01767 Esc = par%vonkar * 1* uster * hhtmp * (1- hhtmp/s%hh(i,j)) ! Van Rijn & Kroon (1992); eq. 10 (at certain depth) 01768 endif 01769 ! 01770 ! iv) Combined mixing coefficient due to currents and waves 01771 E = (Esw**2 + Esc**2)**0.5 ! Van Rijn & Kroon (1992); eq. 5 01772 ! 01773 if (E > 0) then 01774 ! v) Calculcate new concentration 01775 ! a) Explicit first step 01776 ceqexplicit = ceqssteps(ii) - ws0(jg) * ceqssteps(ii)/E * dhhtmp 01777 if (ceqexplicit< 0) then 01778 ceqexplicit = 0 ! numerical limit; overcome overshooting of numerical approximation 01779 endif 01780 ! b) Implicit second step 01781 ceqimplicit = ceqssteps(ii) * (E / (dhhtmp * (ws0(jg)+E/dhhtmp))) ! no limit needed 01782 ! c) Combined method with theta = 0.5 -> second order 01783 ceqssteps(ii+1) = 0.5*ceqexplicit + 0.5*ceqimplicit 01784 if (ceqssteps(ii+1) < 0) then 01785 ceqssteps(ii+1) = 0 ! numerical limit; overcome overshooting of numerical approximation 01786 endif 01787 else 01788 ! No mixing; so no equilibrium in the vertical 01789 ceqssteps = 0 01790 endif 01791 enddo 01792 ! 01793 ! E) Save concentration profile, reference concentraton and reference level 01794 s%refA(i,j) = refA ! reference level 01795 s%ca(i,j) = ca ! reference concentration 01796 do ii = 1, (par%kmax) 01797 s%ccz(i,j,ii) = ceqssteps(ii)*0.5 + ceqssteps(ii+1)*0.5 01798 enddo 01799 ! 01800 ! F) Now determine the depth-averaged concentration (sum divided by total length) 01801 if (s%hh(i,j) < refA) then 01802 ! assume that the depth-averaged concentration is similar to reference concentration 01803 ! not needed anymore; 20% water depth constrain 01804 ceqs(i,j) = ca 01805 do ii = 1, (par%kmax) 01806 s%ccz(i,j,ii) = ca 01807 enddo 01808 else 01809 ! average of the 'active part of the water colum' + reduction factor based on reference layer thickness 01810 ceq_tmp = 0 01811 ceq_tmp = sigz(1) * (ceqssteps(2)*0.5 + ceqssteps(1)*0.5) 01812 do ii = 2, (par%kmax) 01813 ceq_tmp = ceq_tmp + (sigz(ii)-sigz(ii-1)) * (ceqssteps(ii+1)*0.5 + ceqssteps(ii)*0.5) 01814 enddo 01815 ceqs(i,j) =ceq_tmp 01816 endif 01817 else 01818 ceqs(i,j) = 0 01819 endif 01820 enddo 01821 enddo 01822 end if 01823 ! 01824 ceqb = min(ceqb/hloc,par%cmax/2) ! maximum equilibrium bed concentration 01825 s%ceqbg(:,:,jg) = (1-par%bulk)*ceqb*s%sedcal(jg)*s%wetz 01826 ceqs = min(ceqs/hloc,par%cmax/2) ! maximum equilibrium suspended concentration 01827 s%ceqsg(:,:,jg) = (ceqs+par%bulk*ceqb)*s%sedcal(jg)*s%wetz 01828 ! Jaap: old brute method to prevent strong coastline erosion 01829 ! where (hloc<=par%hmin) ceqsg(:,:,jg) = 0.d0 01830 enddo ! end of grain size classes 01831 ! end of grain size classes 01832 01833 end subroutine sedtransform 01834 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01835 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01836 01837 subroutine Nielsen2006(s,par) 01838 01839 use params, only: parameters 01840 use spaceparams 01841 use xmpi_module 01842 use math_tools 01843 use paramsconst 01844 01845 IMPLICIT NONE 01846 01847 type(spacepars),target :: s 01848 type(parameters) :: par 01849 01850 real*8 :: Tsmooth,factime,omegap,dstar 01851 real*8,save :: phirad,delta,khlim,reposerad,reposedzdx 01852 01853 integer :: i,jg 01854 01855 real*8,dimension(:,:),allocatable,save :: dudtsmooth 01856 real*8,dimension(:,:),allocatable,save :: fsed 01857 real*8,dimension(:,:),allocatable,save :: Arms,umeanupd,uvarupd 01858 real*8,dimension(:,:),allocatable,save :: umeanupdphi,uvarupdphi 01859 real*8,dimension(:,:),allocatable,save :: shields,qsedu 01860 real*8,dimension(:,:),allocatable,save :: blphi,facbl,facrw,facslp 01861 real*8,dimension(:,:),allocatable,save :: ulocal,ulocalold,philocal 01862 real*8,dimension(:,:),allocatable,save :: dcfinl,dcfl,fe,cffac,ustar 01863 real*8,dimension(:) ,allocatable,save :: shieldscrit 01864 01865 if (.not. allocated(dudtsmooth)) then 01866 allocate(dudtsmooth(s%nx+1,s%ny+1)) 01867 allocate(fsed(s%nx+1,s%ny+1)) 01868 allocate(shields(s%nx+1,s%ny+1)) 01869 allocate(qsedu(s%nx+1,s%ny+1)) 01870 allocate (blphi (s%nx+1,s%ny+1)) 01871 allocate (facbl (s%nx+1,s%ny+1)) 01872 allocate (facrw (s%nx+1,s%ny+1)) 01873 allocate (facslp (s%nx+1,s%ny+1)) 01874 allocate (ulocal (s%nx+1,s%ny+1)) 01875 allocate (ulocalold (s%nx+1,s%ny+1)) 01876 allocate(philocal(s%nx+1,s%ny+1)) 01877 allocate(umeanupdphi(s%nx+1,s%ny+1)) 01878 allocate(uvarupdphi(s%nx+1,s%ny+1)) 01879 allocate(dcfinl(s%nx+1,s%ny+1)) 01880 allocate(dcfl(s%nx+1,s%ny+1)) 01881 allocate(cffac(s%nx+1,s%ny+1)) 01882 cffac = 0.d0 01883 allocate(Arms(s%nx+1,s%ny+1)) 01884 allocate(umeanupd(s%nx+1,s%ny+1)) 01885 allocate(uvarupd(s%nx+1,s%ny+1)) 01886 allocate(ustar(s%nx+1,s%ny+1)) 01887 allocate(shieldscrit(par%ngd)) 01888 umeanupd = 0.d0 01889 uvarupd = 0.d0 01890 01891 if(par%streaming==1) then 01892 allocate(fe(s%nx+1,s%ny+1)) 01893 endif 01894 fe = 0.d0 01895 fsed = 0.d0 01896 phirad = par%phit/180*par%px 01897 delta = (par%rhos-par%rho)/par%rho 01898 ulocalold = s%uu 01899 ulocal = s%uu 01900 philocal = phirad 01901 umeanupdphi = 0.d0 01902 uvarupdphi = 0.d0 01903 ustar = 0.d0 01904 ! 01905 ! find water depth for which kh(based on Trep) ~ 0.1 01906 khlim = 0.01d0*par%Trep**2*par%g/4/par%px**2 01907 ! 01908 ! Angle of repose 01909 reposerad = par%reposeangle*par%px/180 01910 reposedzdx = tan(reposerad) 01911 ! 01912 ! 01913 do jg=1,par%ngd 01914 if(par%thetcr<0.d0) then 01915 dstar = (par%g*(par%rhos/par%rho-1)/1.d-6/1.d-6)**(1.d0/3)*par%D50(jg) 01916 shieldscrit(jg) = 0.3d0/(1+1.2d0*dstar)+0.055d0*(1-exp(-0.02d0*dstar)) 01917 else 01918 shieldscrit(jg) = par%thetcr 01919 endif 01920 enddo 01921 endif 01922 01923 ! radial velocity of representative wave 01924 omegap = 2*par%px/par%Trep 01925 01926 ! compute local filtered velocity and acceleration term 01927 Tsmooth = par%Trep/20 01928 factime = min(par%dt/Tsmooth,1.d0) 01929 ulocal = (1-factime)*ulocal + factime*s%uu 01930 !ulocal = uu 01931 where(s%wetu==1) 01932 dudtsmooth = (ulocal-ulocalold)/par%dt 01933 elsewhere 01934 dudtsmooth = 0.d0 01935 endwhere 01936 ulocalold = ulocal 01937 ! 01938 ! calculate vertical distribution for turbulence kturb (less turbulence at bottom for larger waterdepth) 01939 if(par%lwt==1 .and. par%yturb>0.d0) then 01940 dcfinl = exp(min(100.d0,(s%hh)/max(par%facthr*s%rolthick,0.01d0))) 01941 dcfl = min(1.d0,1.d0/(dcfinl-1.d0)) 01942 dcfl = max(0.d0,dcfl) 01943 endif 01944 01945 Tsmooth = par%Trep*5 01946 factime = min(par%dt/Tsmooth,1.d0) 01947 01948 umeanupd = (1.d0-factime)*umeanupd + factime * ulocal 01949 uvarupd = (1.d0-factime)*uvarupd + factime * (ulocal-umeanupd)**2 01950 01951 if (par%sedfricfac==SEDFRICFAC_NIELSEN) then 01952 Arms = sqrt(2.d0)/omegap*sqrt(uvarupd) 01953 endif 01954 01955 do jg = 1,par%ngd 01956 ! compute sediment friction factor term 01957 select case (par%sedfricfac) 01958 case (SEDFRICFAC_NIELSEN) 01959 where(s%wetu==1) 01960 fsed = exp(5.5d0*(2.5d0*s%D50(jg)/Arms)**0.2d0-6.3d0) 01961 elsewhere 01962 fsed=0.d0 01963 endwhere 01964 fsed = min(fsed,0.05) 01965 fsed = max(fsed,s%cf) 01966 case (SEDFRICFAC_SWART) 01967 where (s%wetu==1) 01968 fsed = exp(5.213d0*(2.5d0/par%Arms*s%D50(jg))**(0.194)-5.977) 01969 elsewhere 01970 fsed = 0.d0 01971 endwhere 01972 case (SEDFRICFAC_FLOWFRIC) 01973 fsed = s%cf 01974 case(SEDFRICFAC_WILSON) 01975 fsed = 0.114d0*(par%Arms/(par%g*delta*par%Trep**2))*0.4d0 01976 case(SEDFRICFAC_CONSTANT) 01977 fsed = par%fsed 01978 end select 01979 ! 01980 ! compute boundary layer velocity 01981 if (par%phaselag==1) then 01982 where(s%wetu==1) 01983 ! add turbulent term to ustar 01984 ustar = sqrt(fsed/2)*(cos(philocal)*ulocal + & 01985 sin(philocal)/omegap*dudtsmooth) 01986 elsewhere 01987 ustar = 0.d0 01988 endwhere 01989 else 01990 where(s%wetu==1) 01991 ! add turbulent term to ustar 01992 ustar = sqrt(fsed/2)*ulocal 01993 elsewhere 01994 ustar = 0.d0 01995 endwhere 01996 endif 01997 ! 01998 ! Add turbulent stress 01999 if(par%lwt==1 .and. par%yturb>0.d0) then 02000 where(s%wetu==1) 02001 ustar = ustar + sqrt(fsed/2)*(2*sqrt(s%kturb)*dcfl*par%yturb)*sign(1.d0,ustar) 02002 endwhere 02003 endif 02004 ! 02005 ! compute shields parameter 02006 where(s%wetu==1) 02007 shields = ustar**2/(delta*par%g*s%D50(jg)) 02008 elsewhere 02009 shields = 0.d0 02010 endwhere 02011 ! 02012 ! compute bed slope effect 02013 if (par%slopecorr==SLOPECORR_NIELSEN) then 02014 where(s%wetu==1) 02015 where (s%dzbdx>0.d0) 02016 shields = max(shields/cos(min(atan(s%dzbdx),reposerad))-min(s%dzbdx,reposedzdx)*sign(1.d0,ustar),0.d0) 02017 elsewhere(s%dzbdx<0.d0) 02018 shields = max(shields/cos(max(atan(s%dzbdx),-reposerad))-max(s%dzbdx,-reposedzdx)*sign(1.d0,ustar),0.d0) 02019 endwhere 02020 endwhere 02021 elseif (par%slopecorr==SLOPECORR_FREDSOE_DEIGAARD) then 02022 where(s%wetu==1) 02023 where (s%dzbdx>0.d0) 02024 shields = shields*cos(atan(min(s%dzbdx,reposerad)))* & 02025 max(1.d0-sign(1.d0,ustar)*min(s%dzbdx/reposedzdx,1.d0),0.d0) 02026 elsewhere(s%dzbdx<0.d0) 02027 shields = shields*cos(atan(max(s%dzbdx,-reposerad)))* & 02028 max(1.d0-sign(1.d0,ustar)*max(s%dzbdx/reposedzdx,-1.d0),0.d0) 02029 02030 endwhere 02031 endwhere 02032 endif 02033 02034 ! compute streaming effect 02035 if(par%streaming==1) then 02036 ! Note assumes wave propagation is in direction of positive s. Modify to negative 02037 ! contribution in case wave direction is in negative s direction 02038 where(s%wetu==1) 02039 where(uvarupd>1d0*abs(umeanupd)) 02040 fe = exp(5.5d0*(min(170*sqrt(max(shields-0.05d0,0.d0))*s%D50(jg),s%hh)/Arms)**0.2-6.3d0) 02041 elsewhere(uvarupd<0.25d0*abs(umeanupd)) 02042 fe = 0.d0 02043 elsewhere 02044 fe = ((uvarupd/abs(umeanupd))-0.25d0)/0.75d0 * & 02045 exp(5.5d0*(min(170*sqrt(max(shields-0.05d0,0.d0))*s%D50(jg),s%hh)/Arms)**0.2-6.3d0) 02046 endwhere 02047 02048 shields = shields + (2d0/3/par%px*fe*Arms**3*omegap**3/(par%g*s%hh)) / & 02049 (delta*par%g*s%D50(jg)) * & 02050 sign(1.d0,ustar) 02051 shields = max(shields,0.d0) 02052 endwhere 02053 endif 02054 02055 02056 ! compute sediment transport 02057 where(s%wetu==1 .and. shields>shieldscrit(jg)) 02058 qsedu = (par%Ctrans*(shields-shieldscrit(jg))*sqrt(shields))* & 02059 sqrt(delta*par%g*s%D50(jg)**3)*sign(1.d0,ustar) 02060 elsewhere 02061 qsedu = 0.d0 02062 endwhere 02063 02064 !qsedu = min(qsedu,ustar*100*D50(jg)*(1.d0-par%por)) 02065 ! minimum porosity ~0.60 for fluidized bed 02066 where(qsedu>0.d0) 02067 qsedu = min(qsedu,(1.d0-0.60d0)*(abs(s%qx))) 02068 elsewhere(qsedu<0.d0) 02069 qsedu = max(qsedu,(1.d0-0.60d0)*(-abs(s%qx))) 02070 endwhere 02071 02072 02073 if(par%thetanum<1.d0) then 02074 do i=2,s%nx 02075 s%Subg(i,:,jg) = (1.d0-par%thetanum)/2*(qsedu(i-1,:)+qsedu(i+1,:)) + & 02076 par%thetanum * qsedu(i,:) 02077 enddo 02078 if (xmpi_istop) then 02079 s%Subg(1,:,jg) = qsedu(i,:) 02080 endif 02081 if (xmpi_isbot) then 02082 s%Subg(s%nx+1,:,jg) = qsedu(s%nx+1,:) 02083 endif 02084 else 02085 s%Subg(:,:,jg) = qsedu 02086 endif 02087 02088 enddo 02089 02090 end subroutine Nielsen2006 02091 02092 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 02093 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 02094 subroutine mccall_vanrijn(s,par) 02095 02096 use params, only: parameters 02097 use spaceparams 02098 use xmpi_module 02099 use math_tools 02100 use paramsconst 02101 02102 IMPLICIT NONE 02103 02104 type(spacepars),target :: s 02105 type(parameters) :: par 02106 02107 real*8 :: omegap 02108 real*8,save :: phirad,delta,khlim,reposerad,reposedzdx 02109 02110 integer :: i,j,jg 02111 02112 real*8,dimension(:,:),allocatable,save :: dudtsmooth 02113 real*8,dimension(:,:),allocatable,save :: fsed 02114 real*8,dimension(:,:),allocatable,save :: Arms,umeanupd,uvarupd 02115 real*8,dimension(:,:),allocatable,save :: umeanupdphi,uvarupdphi 02116 real*8,dimension(:,:),allocatable,save :: shields,qsedu,qsedutemp,dist 02117 real*8,dimension(:,:),allocatable,save :: blphi,facbl,facrw,facslp,facrwf 02118 real*8,dimension(:,:),allocatable,save :: ulocal,ulocalold,philocal 02119 real*8,dimension(:,:),allocatable,save :: dcfinl,dcfl,cffac 02120 real*8,dimension(:) ,allocatable,save :: shieldscrit,dstar 02121 real*8,dimension(:,:),allocatable,save :: signShields,thetacrlocal 02122 real*8,dimension(:,:),allocatable,save :: phishields,nEF 02123 real*8,dimension(:,:),allocatable,save :: dzbdxf 02124 02125 real*8 :: Te,kvis,Sster,cc1,cc2,wster 02126 real*8 , dimension(:),allocatable,save :: w 02127 02128 if (.not. allocated(dudtsmooth)) then 02129 allocate(dudtsmooth(s%nx+1,s%ny+1)) 02130 allocate(fsed(s%nx+1,s%ny+1)) 02131 allocate(shields(s%nx+1,s%ny+1)) 02132 allocate(qsedu(s%nx+1,s%ny+1)) 02133 allocate(qsedutemp(s%nx+1,s%ny+1)) 02134 allocate(dist(s%nx+1,s%nx+1)) 02135 allocate (blphi (s%nx+1,s%ny+1)) 02136 allocate (facbl (s%nx+1,s%ny+1)) 02137 allocate (facrw (s%nx+1,s%ny+1)) 02138 allocate (facrwf(s%nx+1,s%ny+1)) 02139 allocate (facslp (s%nx+1,s%ny+1)) 02140 allocate (ulocal (s%nx+1,s%ny+1)) 02141 allocate (ulocalold (s%nx+1,s%ny+1)) 02142 allocate(philocal(s%nx+1,s%ny+1)) 02143 allocate(umeanupdphi(s%nx+1,s%ny+1)) 02144 allocate(uvarupdphi(s%nx+1,s%ny+1)) 02145 allocate(dcfinl(s%nx+1,s%ny+1)) 02146 allocate(dcfl(s%nx+1,s%ny+1)) 02147 allocate(cffac(s%nx+1,s%ny+1)) 02148 cffac = 0.d0 02149 allocate(Arms(s%nx+1,s%ny+1)) 02150 allocate(umeanupd(s%nx+1,s%ny+1)) 02151 allocate(uvarupd(s%nx+1,s%ny+1)) 02152 umeanupd = 0.d0 02153 uvarupd = 0.d0 02154 allocate(signShields(s%nx+1,s%ny+1)) 02155 signShields = 0.d0 02156 allocate(thetacrlocal(s%nx+1,s%ny+1)) 02157 thetacrlocal = par%thetcr 02158 allocate(shieldscrit(par%ngd)) 02159 allocate(dstar(par%ngd)) 02160 allocate(dzbdxf(s%nx+1,s%ny+1)) 02161 02162 if (par%form==FORM_WILCOCK_CROW) then 02163 allocate(phishields(s%nx+1,s%ny+1)) 02164 elseif (par%form==FORM_ENGELUND_FREDSOE) then 02165 allocate(nEF(s%nx+1,s%ny+1)) 02166 endif 02167 02168 if(par%bulk==1) then 02169 Te = 20.d0 02170 kvis = 4.d0/(20.d0+Te)*1d-5 ! Van rijn, 1993 02171 allocate(w(par%ngd)) 02172 do jg=1,par%ngd 02173 Sster = s%D50(jg)/(4*kvis)*sqrt((par%rhos/par%rho-1)*par%g*s%D50(jg)) 02174 cc1 = 1.06d0*tanh(0.064d0*Sster*exp(-7.5d0/Sster**2)) 02175 cc2 = 0.22d0*tanh(2.34d0*Sster**(-1.18d0)*exp(-0.0064d0*Sster**2)) 02176 wster = cc1+cc2*Sster 02177 w(jg) = wster*sqrt((par%rhos/par%rho-1.d0)*par%g*s%D50(jg)) 02178 enddo 02179 endif 02180 02181 fsed = 0.d0 02182 phirad = par%phit/180*par%px 02183 delta = (par%rhos-par%rho)/par%rho 02184 ulocalold = s%uu 02185 ulocal = s%uu 02186 philocal = phirad 02187 umeanupdphi = 0.d0 02188 uvarupdphi = 0.d0 02189 ! 02190 ! find water depth for which kh(based on Trep) ~ 0.1 02191 khlim = 0.01d0*par%Trep**2*par%g/4/par%px**2 02192 ! 02193 ! Angle of repose 02194 reposerad = par%reposeangle*par%px/180 02195 reposedzdx = tan(reposerad) 02196 do i=1,s%nx+1 02197 dist(:,i) = abs(s%xu(i,1)-s%xu(:,1)) 02198 dist(:,i) = dist(:,i)/0.25d0 02199 dist(:,i) = max(dist(:,i),1.d0) 02200 dist(:,i) = 1.d0-dist(:,i) 02201 enddo 02202 02203 do jg=1,par%ngd 02204 dstar(jg) = (par%g*(par%rhos/par%rho-1)/1.d-6/1.d-6)**(1.d0/3)*par%D50(jg) 02205 if(par%thetcr<0.d0) then 02206 shieldscrit(jg) = 0.3d0/(1+1.2d0*dstar(jg))+0.055d0*(1-exp(-0.02d0*dstar(jg))) 02207 else 02208 shieldscrit(jg) = par%thetcr 02209 endif 02210 enddo 02211 endif 02212 02213 02214 !do j=1,s%ny+1 02215 ! do i=1,nx 02216 ! Lloc = par%Trep*sqrt(par%g*max(0.1d0,hh(i,j)))/8 02217 ! irange1 = max(1, i-nint(Lloc/dsu(i,j)))! staggered grid means lower should be at least zb(i) 02218 ! irange2 = min(s%nx+1,i+nint(Lloc/dsu(i,j))) 02219 ! irange2 = max(i+1,irange2) ! staggered grid means upper should be at least zb(i+1) 02220 ! dzbdxf(i,j) = (zb(irange2,j)-zb(irange1,j))/(sdist(irange2,j)-sdist(irange1,j)) 02221 ! enddo 02222 !enddo 02223 !dzbdxf(s%nx+1,:) = dzbdxf(nx,:) 02224 dzbdxf = s%dzbdx 02225 02226 do jg = 1,par%ngd 02227 02228 ! compute shields parameter 02229 02230 if (par%gwflow==1 .and. par%inclrelweight==1) then 02231 !facrw = 0.5d0*max(infil/Kz,-1.d0) 02232 facrw = 0.5d0*s%infil/s%Kzinf 02233 else 02234 facrw = 0.d0 02235 endif 02236 02237 where(s%wetu==1) 02238 ! try adding (groundwater) head gradient 02239 shields = (abs(s%taubx-par%incldzdx*s%dzsdx*par%rho*par%g*s%D50(jg))) /(par%rho*par%g*s%D50(jg)*max(delta+facrw,1e-6)) 02240 signShields = sign(1.d0,s%taubx-par%incldzdx*s%dzsdx*par%rho*par%g*s%D50(jg)) 02241 elsewhere 02242 shields = 0.d0 02243 signShields = 0.d0 02244 endwhere 02245 ! 02246 ! compute bed slope effect 02247 if (par%slopecorr==SLOPECORR_NIELSEN) then 02248 where(s%wetu==1) 02249 where (s%dzbdx>0.d0) 02250 shields = max(shields/cos(min(atan(dzbdxf),reposerad))-min(dzbdxf,reposedzdx)*signShields,0.d0) 02251 elsewhere(s%dzbdx<0.d0) 02252 shields = max(shields/cos(max(atan(dzbdxf),-reposerad))-max(dzbdxf,-reposedzdx)*signShields,0.d0) 02253 endwhere 02254 endwhere 02255 elseif (par%slopecorr==SLOPECORR_FREDSOE_DEIGAARD) then 02256 where(s%wetu==1) 02257 where (s%dzbdx>0.d0) 02258 shields = shields*cos(atan(min(dzbdxf,reposerad)))* & 02259 max(1.d0-signShields*min(dzbdxf/reposedzdx,1.d0),0.d0) 02260 elsewhere(s%dzbdx<0.d0) 02261 shields = shields*cos(atan(max(dzbdxf,-reposerad)))* & 02262 max(1.d0-signShields*max(dzbdxf/reposedzdx,-1.d0),0.d0) 02263 endwhere 02264 endwhere 02265 endif 02266 02267 !! compute streaming effect 02268 !if(par%streaming==1) then 02269 ! ! Note assumes wave propagation is in direction of positive s. Modify to negative 02270 ! ! contribution in case wave direction is in negative s direction 02271 ! Arms = sqrt(2.d0)/omegap*sqrt(uvarupd) 02272 ! where(s%wetu==1) 02273 ! where(uvarupd>1d0*abs(umeanupd)) 02274 ! fe = exp(5.5d0*(min(170*sqrt(max(shields-0.05d0,0.d0))*s%D50(jg),s%hh)/Arms)**0.2-6.3d0) 02275 ! elsewhere(uvarupd<0.25d0*abs(umeanupd)) 02276 ! fe = 0.d0 02277 ! elsewhere 02278 ! fe = ((uvarupd/abs(umeanupd))-0.25d0)/0.75d0 * & 02279 ! exp(5.5d0*(min(170*sqrt(max(shields-0.05d0,0.d0))*s%D50(jg),s%hh)/Arms)**0.2-6.3d0) 02280 ! endwhere 02281 ! 02282 ! shields = shields + (2d0/3/par%px*fe*Arms**3*omegap**3/(par%g*s%hh)) / & 02283 ! ((delta+facrw)*par%g*s%D50(jg)) * & 02284 ! signShields 02285 ! shields = max(shields,0.d0) 02286 ! endwhere 02287 !endif 02288 02289 thetacrlocal = shieldscrit(jg) 02290 02291 ! compute sediment transport 02292 select case (par%form) 02293 case (FORM_MPM) 02294 par%Ctrans = 8.d0 02295 where(s%wetu==1 .and. shields>thetacrlocal) 02296 qsedu = (par%Ctrans*(shields-thetacrlocal)**1.5d0)* & 02297 sqrt(delta*par%g*s%D50(jg)**3)*signShields 02298 elsewhere 02299 qsedu = 0.d0 02300 endwhere 02301 case (FORM_WONG_PARKER) 02302 par%Ctrans = 3.97d0 02303 where(s%wetu==1 .and. shields>thetacrlocal) 02304 qsedu = (par%Ctrans*(shields-thetacrlocal)**1.5d0)* & 02305 sqrt(delta*par%g*s%D50(jg)**3)*signShields 02306 elsewhere 02307 qsedu = 0.d0 02308 endwhere 02309 case (FORM_FL_VB) ! Fernandez Luque, R., and van Beek, R.(1976). Erosion and Transport of Bed-Load 02310 ! Sediment. J. Hydraul, Res., 14(2), 127-144 02311 par%Ctrans = 5.7d0 02312 where(s%wetu==1 .and. shields>thetacrlocal) 02313 qsedu = (par%Ctrans*(shields-thetacrlocal)**1.5d0)* & 02314 sqrt(delta*par%g*s%D50(jg)**3)*signShields 02315 elsewhere 02316 qsedu = 0.d0 02317 endwhere 02318 case (FORM_WILCOCK_CROW) 02319 where(s%wetu==1) 02320 phishields = shields/thetacrlocal 02321 where(phishields>=1.35d0) 02322 qsedu = 14*(1.d0-0.894d0/sqrt(phishields))**4.5d0*shields**1.5*sqrt(delta*par%g)*s%D50(jg)**1.5*signShields 02323 elsewhere 02324 qsedu = 0.002d0*phishields**7.5d0*shields**1.5*sqrt(delta*par%g)*s%D50(jg)**1.5*signShields 02325 endwhere 02326 elsewhere 02327 qsedu = 0.d0 02328 endwhere 02329 case (FORM_ENGELUND_FREDSOE) 02330 where(s%wetu==1 .and. shields>thetacrlocal) 02331 ! note: pi/6 ~ 0.5236, dynamic friction ~ 0.5-0.65 (Fredsoe & Deigaard) 02332 nEF = (1.d0+(0.5236d0*0.5d0/(shields-thetacrlocal))**4)**-0.25d0 02333 qsedu = 5*nEF*(sqrt(shields)-0.7d0*sqrt(thetacrlocal))*sqrt(delta*par%g*s%D50(jg)**3)*signShields 02334 elsewhere 02335 qsedu = 0.d0 02336 endwhere 02337 case (FORM_FREDSOE_DEIGAARD) ! Equation 7.59 02338 where(s%wetu==1 .and. shields>thetacrlocal) 02339 ! note: 30/pi = 9.5493, mu = 0.65 -> 30/pi/mu = 14.6912 02340 qsedu = 14.6912d0*(shields-thetacrlocal)*(sqrt(shields)-0.7d0*sqrt(thetacrlocal)) * & 02341 sqrt(delta*par%g*s%D50(jg)**3)*signShields 02342 elsewhere 02343 qsedu = 0.d0 02344 endwhere 02345 case(FORM_MCCALL_VANRIJN) ! Eq 10. note: devision by rhos to get m2/s transport rate 02346 par%Ctrans = 0.5d0 02347 where(s%wetu==1 .and. shields>thetacrlocal) 02348 qsedu = par%Ctrans*s%D50(jg)*dstar(jg)**(-0.3d0)*sqrt(abs(s%taubx)/par%rho)* & 02349 ((shields-thetacrlocal)/thetacrlocal)*signShields 02350 elsewhere 02351 qsedu = 0.d0 02352 endwhere 02353 end select 02354 02355 02356 where(qsedu>0.d0) 02357 qsedu = qsedu * par%uprushfac 02358 elsewhere 02359 qsedu = qsedu * par%backwashfac 02360 endwhere 02361 02362 !qsedu = min(qsedu,ustar*100*D50(jg)*(1.d0-par%por)) 02363 ! minimum porosity ~0.60 for fluidized bed 02364 where(qsedu>0.d0) 02365 qsedu = min(qsedu,(1.d0-0.60d0)*(abs(s%qx))) 02366 elsewhere(qsedu<0.d0) 02367 qsedu = max(qsedu,(1.d0-0.60d0)*(-abs(s%qx))) 02368 endwhere 02369 02370 if (par%bulk==0) then 02371 if(par%thetanum<1.d0) then 02372 do i=2,s%nx 02373 s%Subg(i,:,jg) = (1.d0-par%thetanum)/2*(qsedu(i-1,:)+qsedu(i+1,:)) + & 02374 par%thetanum * qsedu(i,:) 02375 enddo 02376 if(xmpi_istop) s%Subg(1,:,jg) = qsedu(i,:) 02377 if(xmpi_isbot) s%Subg(s%nx+1,:,jg) = qsedu(s%nx+1,:) 02378 else 02379 s%Subg(:,:,jg) = qsedu 02380 endif 02381 ! To Do: fix this properly for MPI 02382 do j=1,s%ny+1 02383 do i=1,s%nx+1 02384 if(s%Subg(i,j,jg)>0.d0) then 02385 s%Subg(i,j,jg) = s%Subg(i,j,jg)*s%sedcal(jg)*s%pbbed(i,j,1,jg) 02386 else 02387 s%Subg(i,j,jg) = s%Subg(i,j,jg)*s%sedcal(jg)*s%pbbed(min(i+1,s%nx+1),j,1,jg) 02388 endif 02389 enddo 02390 enddo 02391 else 02392 if (par%oldTsmin==1) then 02393 s%Tsg(:,:,jg) = max(par%tsfac*s%hu/w(jg),par%Tsmin) 02394 else 02395 s%Tsg(:,:,jg) = max(par%tsfac*s%hu/w(jg),par%dtlimTs*par%dt) 02396 endif 02397 !s%Tsg(:,:,jg) = max(par%tsfac*s%hu/w(jg),par%Tsmin) 02398 s%ceqsg(:,:,jg) = qsedu/s%hu/s%uu 02399 s%ceqbg(:,:,jg) = 0.d0 02400 endif 02401 enddo 02402 02403 end subroutine mccall_vanrijn 02404 02405 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 02406 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 02407 02408 02409 subroutine waveturb(s,par) 02410 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 02411 ! Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University ! 02412 ! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski, ! 02413 ! Jaap van Thiel de Vries, Robert McCall ! 02414 ! ! 02415 ! d.roelvink@unesco-ihe.org ! 02416 ! UNESCO-IHE Institute for Water Education ! 02417 ! P.O. Box 3015 ! 02418 ! 2601 DA Delft ! 02419 ! The Netherlands ! 02420 ! ! 02421 ! This library is free software; you can redistribute it and/or ! 02422 ! modify it under the terms of the GNU Lesser General Public ! 02423 ! License as published by the Free Software Foundation; either ! 02424 ! version 2.1 of the License, or (at your option) any later version. ! 02425 ! ! 02426 ! This library is distributed in the hope that it will be useful, ! 02427 ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! 02428 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! 02429 ! Lesser General Public License for more details. ! 02430 ! ! 02431 ! You should have received a copy of the GNU Lesser General Public ! 02432 ! License along with this library; if not, write to the Free Software ! 02433 ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 ! 02434 ! USA ! 02435 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 02436 use params, only: parameters 02437 use spaceparams 02438 use xmpi_module 02439 use paramsconst, only: TURBADV_EULERIAN, TURBADV_LAGRANGIAN, TURBADV_NONE 02440 02441 implicit none 02442 02443 type(spacepars),target :: s 02444 type(parameters) :: par 02445 02446 integer :: i 02447 integer :: j 02448 real*8 :: ML, disturb 02449 real*8, save :: twothird 02450 real*8,dimension(:,:),allocatable,save :: ksource, kturbu,kturbv,Sturbu,Sturbv,dzsdt_cr 02451 02452 !include 's.ind' 02453 !include 's.inp' 02454 02455 if (.not. allocated(kturbu)) then 02456 allocate(ksource (s%nx+1,s%ny+1)) 02457 allocate(kturbu (s%nx+1,s%ny+1)) 02458 allocate(kturbv (s%nx+1,s%ny+1)) 02459 allocate(Sturbu (s%nx+1,s%ny+1)) 02460 allocate(Sturbv (s%nx+1,s%ny+1)) 02461 allocate(dzsdt_cr (s%nx+1,s%ny+1)) 02462 twothird=2.d0/3.d0 02463 endif 02464 ! use lagrangian velocities 02465 kturbu = 0.0d0 !Jaap 02466 kturbv = 0.0d0 !Jaap 02467 ! dzsdt_cr=par%beta*s%c 02468 dzsdt_cr=par%beta*sqrt(par%g*s%hh) 02469 ! Update roller thickness 02470 !s%rolthick=s%rolthick+par%dt*(abs(s%dzsdt)-dzsdt_cr) Dano: not abs! 02471 s%rolthick=s%rolthick+par%dt*(s%dzsdt-dzsdt_cr) 02472 s%rolthick=max(s%rolthick,0.d0) 02473 s%rolthick=min(s%rolthick,s%hh) 02474 where (s%wetz==0) 02475 s%rolthick=0.d0 02476 endwhere 02477 02478 ! Jaap compute sources and sinks 02479 ! long wave source 02480 ksource = 0 02481 if (par%lwt == 1) then 02482 ksource = par%g*s%rolthick*par%beta*sqrt(par%g*s%hh) !only important in shallow water, where s%c=sqrt(gh) 02483 endif 02484 if (par%turbadv == TURBADV_NONE) then 02485 s%kturb = (s%DR/par%rho)**twothird ! See Battjes, 1975 / 1985 02486 else 02487 ksource = ksource + s%DR/par%rho 02488 02489 02490 ! Jaap do long wave turb approach for both short waves and long waves 02491 ! 02492 ! Turbulence in uu-points 02493 ! 02494 do j=1,s%ny+1 02495 do i=1,s%nx 02496 if(s%uu(i,j)>0.d0) then 02497 kturbu(i,j)=par%thetanum*s%kturb(i,j)+(1.d0-par%thetanum)*s%kturb(min(i+1,s%nx),j) 02498 elseif(s%uu(i,j)<0.d0) then 02499 kturbu(i,j)=par%thetanum*s%kturb(i+1,j)+(1.d0-par%thetanum)*s%kturb(max(i,2),j) 02500 else 02501 kturbu(i,j)=0.5d0*(s%kturb(i,j)+s%kturb(i+1,j)) 02502 endif 02503 enddo 02504 enddo 02505 if (xmpi_isbot) kturbu(s%nx+1,:) = s%kturb(s%nx+1,:) 02506 ! 02507 ! Turbulence in vv-points 02508 ! 02509 do j=1,s%ny 02510 do i=1,s%nx+1 02511 if(s%vv(i,j)>0) then 02512 kturbv(i,j)=par%thetanum*s%kturb(i,j)+(1.d0-par%thetanum)*s%kturb(i,min(j+1,s%ny)) 02513 else if(s%vv(i,j)<0) then 02514 kturbv(i,j)=par%thetanum*s%kturb(i,j+1)+(1.d0-par%thetanum)*s%kturb(i,max(j,2)) 02515 else 02516 kturbv(i,j)=0.5d0*(kturbv(i,j)+kturbv(i,j+1)) 02517 endif 02518 enddo 02519 enddo 02520 kturbv(:,s%ny+1) = s%kturb(:,s%ny+1) !Robert 02521 ! 02522 ! Turbulence advection in X and Y direction 02523 ! 02524 if (par%turbadv == TURBADV_LAGRANGIAN) then 02525 Sturbu=kturbu*s%uu*s%hu*s%wetu 02526 Sturbv=kturbv*s%vv*s%hv*s%wetv 02527 elseif (par%turbadv == TURBADV_EULERIAN) then 02528 Sturbu=kturbu*s%ueu*s%hu*s%wetu 02529 Sturbv=kturbv*s%vev*s%hv*s%wetv 02530 endif 02531 ! 02532 ! Update turbulence 02533 ! 02534 if (s%ny>0) then 02535 do j=2,s%ny+1 02536 do i=2,s%nx+1 02537 if (par%betad>0) then 02538 disturb=par%betad*s%kturb(i,j)**1.5d0 02539 else 02540 !ML=max(s%rolthick(i,j),0.07d0*s%hstokes(i,j)) 02541 !disturb=0.08d0*s%hstokes(i,j)/ML*s%kturb(i,j)**1.5d0 02542 endif 02543 s%kturb(i,j) = (s%hold(i,j)*s%kturb(i,j)-par%dt*( & 02544 (Sturbu(i,j)*s%dnu(i,j)-Sturbu(i-1,j)*s%dnu(i-1,j)+& 02545 Sturbv(i,j)*s%dsv(i,j)-Sturbv(i,j-1)*s%dsv(i,j-1))*s%dsdnzi(i,j)-& 02546 (ksource(i,j)-disturb)))/s%hh(i,j) 02547 s%kturb(i,j)=max(s%kturb(i,j),0.0d0) 02548 02549 enddo 02550 enddo 02551 else 02552 j=1 02553 do i=2,s%nx+1 02554 ! Robert: perhaps better to re-write to implicit sink term 02555 s%kturb(i,j) = (s%hold(i,j)*s%kturb(i,j)-par%dt*( & 02556 (Sturbu(i,j)*s%dnu(i,j)-Sturbu(i-1,j)*s%dnu(i-1,j))*s%dsdnzi(i,j)-& 02557 (ksource(i,j)-par%betad*s%kturb(i,j)**1.5d0)))/s%hh(i,j) 02558 s%kturb(i,j)=max(s%kturb(i,j),0.0d0) 02559 02560 enddo 02561 endif 02562 02563 !Robert: cannot divide by hh here again 02564 !s%kturb = s%kturb/s%hh !max(s%hh,0.01d0) : 02565 02566 ! Jaap only required for advection mode? 02567 if (xmpi_istop) s%kturb(1,:)=s%kturb(2,:) 02568 if (xmpi_isbot) s%kturb(s%nx+1,:)=s%kturb(s%nx+1-1,:) 02569 if (s%ny>0) then 02570 if (xmpi_isleft) s%kturb(:,1)=s%kturb(:,2) 02571 if (xmpi_isright) s%kturb(:,s%ny+1)=s%kturb(:,s%ny+1-1) 02572 endif 02573 02574 #ifdef USEMPI 02575 call xmpi_shift_ee(s%kturb) 02576 #endif 02577 02578 endif ! turbadv == 'none' 02579 02580 end subroutine waveturb 02581 02582 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 02583 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 02584 02585 02586 subroutine RvR(s,par) 02587 02588 use params, only: parameters 02589 use spaceparams 02590 use xmpi_module 02591 02592 implicit none 02593 02594 type(spacepars),target :: s 02595 type(parameters) :: par 02596 02597 real*8 , save :: m1,m2,m3,m4,m5,m6,alpha,beta 02598 02599 real*8 , dimension(:,:),allocatable,save :: Urs,Bm,B1 02600 02601 !include 's.ind' 02602 !include 's.inp' 02603 02604 ! only in first timestep.. 02605 if (.not. allocated(Urs)) then 02606 02607 allocate (Urs (s%nx+1,s%ny+1)) 02608 allocate (Bm (s%nx+1,s%ny+1)) 02609 allocate (B1 (s%nx+1,s%ny+1)) 02610 02611 m1 = 0; ! a = 0 02612 m2 = 0.7939; ! b = 0.79 +/- 0.023 02613 m3 = -0.6065; ! c = -0.61 +/- 0.041 02614 m4 = 0.3539; ! d = -0.35 +/- 0.032 02615 m5 = 0.6373; ! e = 0.64 +/- 0.025 02616 m6 = 0.5995; ! f = 0.60 +/- 0.043 02617 alpha = -log10(exp(1.d0))/m4 02618 beta = exp(m3/m4) 02619 02620 endif 02621 02622 Urs = 3.d0/8.d0*sqrt(2.d0)*s%H*s%k/(s%k*s%hh)**3 !Ursell number 02623 Urs = max(Urs,0.000000000001d0) 02624 Bm = m1 + (m2-m1)/(1.d0+beta*Urs**alpha) !Boltzmann sigmoid (eq 6) 02625 B1 = (-90.d0+90.d0*tanh(m5/Urs**m6))*par%px/180.d0 02626 s%Sk = Bm*cos(B1) !Skewness (eq 8) 02627 s%As = Bm*sin(B1) !Asymmetry(eq 9) 02628 s%ua = par%sws*(par%facSk*s%Sk-par%facAs*s%As)*s%urms 02629 02630 ! multiply Sk and As with wetz to get zeros at h = 0 for output 02631 s%Sk = s%Sk*s%wetz 02632 s%As = s%As*s%wetz 02633 02634 end subroutine RvR 02635 02636 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 02637 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 02638 02639 subroutine vT(s,par) 02640 02641 use params, only: parameters 02642 use spaceparams 02643 use readkey_module 02644 use xmpi_module 02645 02646 implicit none 02647 02648 type(spacepars),target :: s 02649 type(parameters) :: par 02650 02651 integer :: i,j 02652 integer , save :: nh,nt 02653 integer :: ih0,it0,ih1,it1 02654 real*8 :: p,q,f0,f1,f2,f3 02655 real*8 :: t0fac,siguref,duddtmax,dudtmax,duddtmean,dudtmean,detadxmean 02656 real*8 , save :: dh,dt 02657 02658 real*8 , dimension(:,:),allocatable ,save :: h0,t0,detadxmax 02659 ! Robert: RF table now included in source code, rather than read from file 02660 ! Rienecker Fenton table with amongst others amplitudes non-linear components obtained with stream function theory 02661 include 'RF.inc' 02662 ! Robert: 'RF.inc' contains definition of RF as real*8(5,33,40) with "parameter" attribute 02663 ! so RF values may not be modified! To save memory, only rows 13,14,15,16 and 18 of the 02664 ! original matrix are stored. So new row 1 corresponds with old row 13, etc. 02665 02666 !include 's.ind' 02667 !include 's.inp' 02668 02669 02670 ! only in first timestep.. 02671 if (.not. allocated(h0)) then 02672 allocate (h0 (s%nx+1,s%ny+1)) 02673 allocate (t0 (s%nx+1,s%ny+1)) 02674 allocate (detadxmax (s%nx+1,s%ny+1)) 02675 dh = 0.03d0 02676 dt = 1.25d0 02677 nh = floor(0.99d0/dh); 02678 nt = floor(50.d0/dt); 02679 endif 02680 02681 ! non-linearity of short waves is listed in table as function of dimensionless wave height h0 and dimensionless wave period t0 02682 02683 ! compute dimensionless wave height and wave period in each grid point.. 02684 h0 = min(nh*dh,max(dh,min(s%H,s%hh)/s%hh)) 02685 t0 = min(nt*dt,max(dt,par%Trep*sqrt(par%g/s%hh))) 02686 02687 ! estimate Sk, As and ua by interpolating table values 02688 do j=1,s%ny+1 02689 do i=1,s%nx+1 02690 ! interpolate table values.... 02691 ih0=floor(h0(i,j)/dh); 02692 it0=floor(T0(i,j)/dt); 02693 ih1=min(ih0+1,nh); 02694 it1=min(it0+1,nt); 02695 p=(h0(i,j)-ih0*dh)/dh; 02696 q=(T0(i,j)-it0*dt)/dt; 02697 02698 f0=(1-p)*(1-q); 02699 f1=p*(1-q); 02700 f2=q*(1-p); 02701 f3=p*q; 02702 02703 ! Skewness and assymetry 02704 s%Sk(i,j) = f0*RF(1,ih0,it0)+f1*RF(1,ih1,it0)+ f2*RF(1,ih0,it1)+f3*RF(1,ih1,it1) 02705 s%As(i,j) = f0*RF(2,ih0,it0)+f1*RF(2,ih1,it0)+ f2*RF(2,ih0,it1)+f3*RF(2,ih1,it1) 02706 02707 ! Sediment advection velocity from Skewness and Assymetry 02708 ! ua(i,j) = par%sws*par%facua*(Sk(i,j)-As(i,j))*urms(i,j) 02709 s%ua(i,j) = par%sws*(par%facSk*s%Sk(i,j)-par%facAs*s%As(i,j))*s%urms(i,j) 02710 02711 ! Estimate bore period Tbore and mean slope bore front to feeded back in roller energy balance 02712 02713 ! correct slope in case 1.25>T0>50 02714 if (t0(i,j)==50.d0) then 02715 t0fac = 50.d0/max((par%Trep*sqrt(par%g/s%hh(i,j))),50.d0) 02716 elseif (t0(i,j)==1.25)then 02717 t0fac = 1.25d0/min((par%Trep*sqrt(par%g/s%hh(i,j))),1.25d0) 02718 else 02719 t0fac = 1.d0 02720 endif 02721 02722 ! detadxmax for Tbore... 02723 ! dimnesionless maximum acceleration under bore front 02724 duddtmax = f0*RF(3,ih0,it0)+f1*RF(3,ih1,it0)+ f2*RF(3,ih0,it1)+f3*RF(3,ih1,it1) 02725 siguref = f0*RF(4,ih0,it0)+f1*RF(4,ih1,it0)+ f2*RF(4,ih0,it1)+f3*RF(4,ih1,it1) 02726 ! translate dimensionless duddtmax to real world dudtmax 02727 ! /scale with variance and go from [-] to [m/s^2] /tableb./dimensionless dudtmax 02728 dudtmax = s%urms(i,j)/max(par%eps,siguref)*sqrt(par%g/s%hh(i,j))*t0fac*duddtmax 02729 detadxmax(i,j) = dudtmax*sinh(s%k(i,j)*s%hh(i,j))/max(max(s%c(i,j),sqrt(s%H(i,j)*par%g)),1d-10)/s%sigm(i,j) 02730 02731 ! detadxmean for roller energy balance dissipation... 02732 if (par%rfb==1) then 02733 duddtmean = f0*RF(5,ih0,it0)+f1*RF(5,ih1,it0)+ f2*RF(5,ih0,it1)+f3*RF(5,ih1,it1) 02734 dudtmean = s%urms(i,j)/max(par%eps,siguref)*sqrt(par%g/s%hh(i,j))*t0fac*duddtmean 02735 detadxmean = dudtmean*sinh(s%k(i,j)*s%hh(i,j))/max(s%c(i,j),sqrt(s%H(i,j)*par%g))/s%sigm(i,j) 02736 s%BR(i,j) = par%BRfac*sin(atan(detadxmean)) 02737 endif 02738 02739 enddo 02740 enddo 02741 02742 s%Tbore = max(par%Trep/25.d0,min(par%Trep/4.d0,s%H/(max(max(s%c,sqrt(s%H*par%g)),1d-10)*max(detadxmax,par%eps)))) 02743 s%Tbore = par%Tbfac*s%Tbore 02744 02745 ! multiply Sk and As with wetz to get zeros at h = 0 for output 02746 s%Sk = s%Sk*s%wetz 02747 s%As = s%As*s%wetz 02748 02749 end subroutine vT 02750 02751 subroutine setbathy_update(s, par) 02752 02753 use params, only: parameters 02754 use spaceparams 02755 use interp, only: linear_interp 02756 02757 implicit none 02758 02759 type(spacepars) :: s 02760 type(parameters) :: par 02761 02762 integer :: i,j,dummy 02763 real*8,dimension(s%nx+1,s%ny+1) :: zbnew 02764 02765 ! interpolate from file 02766 do j=1,s%ny+1 02767 do i=1,s%nx+1 02768 call LINEAR_INTERP(s%tsetbathy,s%setbathy(i,j,:),par%nsetbathy, & 02769 par%t,zbnew(i,j),dummy) 02770 enddo 02771 enddo 02772 ! update water level 02773 s%zs = s%zs+zbnew-s%zb 02774 ! update bed level 02775 s%zb = zbnew 02776 s%hh=max(s%zs-s%zb,par%eps) 02777 ! update wet and dry cells 02778 where (s%hh>par%eps) 02779 s%wetz=1 02780 elsewhere 02781 s%wetz=0 02782 s%zs = s%zb+par%eps 02783 s%hh = par%eps 02784 endwhere 02785 02786 end subroutine setbathy_update 02787 02788 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 02789 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 02790 02791 subroutine hybrid(s,par) 02792 02793 use params, only: parameters 02794 use interp, only: linear_interp 02795 use spaceparams 02796 use xmpi_module 02797 02798 implicit none 02799 02800 type(spacepars),target :: s 02801 type(parameters) :: par 02802 02803 integer :: i,j,j1,indx,first, nIter, maxIter 02804 integer , dimension(:), allocatable,save :: slopeind 02805 real*8 , dimension(:), allocatable,save :: hav1d 02806 real*8 :: irrb,runup_old 02807 02808 !include 's.ind' 02809 !include 's.inp' 02810 02811 if (.not. allocated(hav1d)) then 02812 allocate(hav1d (s%nx+1)) 02813 allocate(slopeind (s%nx+1)) 02814 endif 02815 02816 do j=1,s%ny+1 02817 indx = s%nx+1 02818 first = 0 02819 do i=1,s%nx 02820 if (s%wetz(i,j)-s%wetz(max(i-1,1),j)==-1 .and. first==0) then ! transition from wet to dry 02821 ! only consider first dry point 02822 first = 1 02823 ! find wave height for runup at L1 meter from water line 02824 call linear_interp(s%xz(:,j),s%H(:,j),s%nx+1,s%xz(i-1,j)-s%L1(i-1,j),s%Hrunup(j),indx) 02825 ! Find toe of runup slope if present (dzbdx > 0.15). 02826 ! If not present Hrunup will converge to H at the water line (where H = 0 per definition) 02827 do j1=indx,i-1 02828 ! TODO: Strange condition. In case of no toe, or no steep slope, 02829 ! there will still be extra turbulence at L1 meter from the water line... 02830 ! cross shore location structure toe 02831 if (s%dzbdx(j1,j)<0.15d0 .or. s%structdepth(j1,j)>0.1d0) then 02832 indx = j1 02833 endif 02834 enddo 02835 ! update Hrunup and runup x-location 02836 s%Hrunup(j) = s%H(indx,j) ! short wave height at revetment toe 02837 s%xHrunup(j) = s%xz(indx,j) ! cross-shore location revetment toe 02838 s%istruct(j) = indx*1.d0 ! cross-shore index revetment toe 02839 s%iwl(j) = (i-1)*1.d0 ! cross-shore location waterline (inlcuding lw-s%runup) 02840 02841 ! now iteratively compute runup 02842 hav1d = s%hh(:,j) 02843 runup_old = huge(0.d0) 02844 s%runup(j) = 0; 02845 nIter = 0; 02846 maxIter = 50; 02847 do while (abs(s%runup(j)-runup_old)>0.01d0 .and. nIter < maxIter) 02848 nIter = nIter +1; 02849 runup_old = s%runup(j) 02850 slopeind = 0 02851 where (hav1d>par%eps .and. s%dzbdx(:,j)>0.15d0) 02852 slopeind = 1 02853 endwhere 02854 s%strucslope(j) = sum(s%dzbdx(indx:s%nx,j)*s%dsu(indx:s%nx,j)*slopeind(indx:s%nx))/ & 02855 max(par%eps,sum(s%dsu(indx:s%nx,j)*slopeind(indx:s%nx))) 02856 if (s%strucslope(j) > 0.d0) then 02857 irrb = s%strucslope(j)/sqrt(2*par%px*max(s%Hrunup(j),par%eps)/par%g/par%Trep**2) 02858 s%runup(j) = par%facrun*min(irrb,2.3d0)*s%Hrunup(j)*cos(2*par%px/par%Trep*par%t) 02859 else 02860 s%runup(j) = 0.d0; 02861 endif 02862 ! This triggers s%runup to be calculated for the complete hinterland (also behind a dune). 02863 ! Only values calculated for the first dune front are used in the avalanching algorithm (morphevolution) 02864 hav1d = s%wetz(:,j)*max(par%eps,(s%hh(:,j) + s%runup(j))) + & 02865 (1.d0-s%wetz(:,j))*max(par%eps,s%runup(j)+s%zs(i-1,j)-s%zb(:,j) ) 02866 enddo 02867 endif 02868 enddo 02869 enddo 02870 02871 02872 02873 end subroutine hybrid 02874 02875 end module morphevolution
1.7.6.1
