|
XBeach
|
00001 module flow_timestep_module 00002 implicit none 00003 private 00004 public flow 00005 save 00006 contains 00007 subroutine flow(s,par) 00008 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00009 ! Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University ! 00010 ! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski, ! 00011 ! Jaap van Thiel de Vries, Robert McCall ! 00012 ! ! 00013 ! d.roelvink@unesco-ihe.org ! 00014 ! UNESCO-IHE Institute for Water Education ! 00015 ! P.O. Box 3015 ! 00016 ! 2601 DA Delft ! 00017 ! The Netherlands ! 00018 ! ! 00019 ! This library is free software; you can redistribute it and/or ! 00020 ! modify it under the terms of the GNU Lesser General Public ! 00021 ! License as published by the Free Software Foundation; either ! 00022 ! version 2.1 of the License, or (at your option) any later version. ! 00023 ! ! 00024 ! This library is distributed in the hope that it will be useful, ! 00025 ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! 00026 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! 00027 ! Lesser General Public License for more details. ! 00028 ! ! 00029 ! You should have received a copy of the GNU Lesser General Public ! 00030 ! License along with this library; if not, write to the Free Software ! 00031 ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 ! 00032 ! USA ! 00033 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00034 use params, only: parameters 00035 use spaceparamsdef 00036 use xmpi_module 00037 use boundaryconditions, only: flow_lat_bc, discharge_boundary_h, discharge_boundary_v 00038 use flow_secondorder_module, only: flow_secondorder_huhv, flow_secondorder_advuv 00039 use nonh_module, only: nonh_cor 00040 use bedroughness_module, only: bedroughness_init, bedroughness_update 00041 use vsm_u_xb_module, only: vsm_u_XB 00042 use paramsconst 00043 use rainfall_module 00044 00045 IMPLICIT NONE 00046 00047 type(spacepars),target :: s 00048 type(parameters) :: par 00049 00050 integer :: i 00051 integer :: j,j1,jp1,jm1 00052 real*8,dimension(:,:),allocatable,save :: vv_old !Velocity at previous timestep 00053 real*8,dimension(:,:),allocatable,save :: uu_old !Velocity at previous timestep 00054 real*8,dimension(:,:),allocatable,save :: zs_old 00055 real*8,dimension(:,:),allocatable,save :: vsu,usu,vsv,usv,veu,uev 00056 real*8,dimension(:,:),allocatable,save :: dudx,dvdy 00057 real*8,dimension(:,:),allocatable,save :: us,vs 00058 real*8 :: nuh1,nuh2 00059 real*8 :: dudx1,dudx2,dudy1,dudy2 00060 real*8 :: dvdy1,dvdy2,dvdx1,dvdx2 !Jaap 00061 real*8 :: dalfa !difference in grid angles 00062 real*8 :: uin,vin !s%u resp s%v-velocity corrected for grid angle change 00063 real*8 :: qin !specific discharge entering cell 00064 real*8 :: dzsdnavg !alongshore water level slope 00065 real*8,save :: fc 00066 real*8,dimension(:,:),allocatable,save :: sinthm,costhm 00067 real*8 :: fcvisc=0.1d0,facdel=5.d0,facdf=1.d0,ks 00068 real*8 :: tauw,tauwx,tauwy 00069 integer :: imax,jmax,jmin,swglm=0 00070 00071 00072 if (.not. allocated(vsu) ) then 00073 allocate ( vsu(s%nx+1,s%ny+1)) 00074 allocate ( usu(s%nx+1,s%ny+1)) 00075 allocate ( vsv(s%nx+1,s%ny+1)) 00076 allocate ( usv(s%nx+1,s%ny+1)) 00077 allocate ( veu(s%nx+1,s%ny+1)) 00078 allocate ( uev(s%nx+1,s%ny+1)) 00079 allocate ( dudx(s%nx+1,s%ny+1)) 00080 allocate ( dvdy(s%nx+1,s%ny+1)) 00081 allocate ( us(s%nx+1,s%ny+1)) 00082 allocate ( vs(s%nx+1,s%ny+1)) 00083 allocate (sinthm(s%nx+1,s%ny+1)) 00084 allocate (costhm(s%nx+1,s%ny+1)) 00085 ! 00086 if (par%secorder == 1 .or. par%wavemodel==WAVEMODEL_NONH) then 00087 allocate(vv_old(s%nx+1,s%ny+1)); vv_old = s%vv 00088 allocate(uu_old(s%nx+1,s%ny+1)); uu_old = s%uu 00089 allocate(zs_old(s%nx+1,s%ny+1)); zs_old = s%zs 00090 endif 00091 00092 s%vu =0.d0 00093 vsu =0.d0 00094 usu =0.d0 00095 vsv =0.d0 00096 usv =0.d0 00097 s%uv =0.d0 00098 veu =0.d0 00099 uev =0.d0 00100 s%ueu =0.d0 00101 s%vev =0.d0 00102 dudx =0.d0 00103 s%ududx =0.d0 00104 dvdy =0.d0 00105 s%vdvdy =0.d0 00106 s%udvdx =0.d0 00107 s%vdudy =0.d0 00108 s%viscu =0.d0 00109 s%viscv =0.d0 00110 us =0.d0 00111 vs =0.d0 00112 s%u =0.d0 00113 s%v =0.d0 00114 s%ue =0.d0 00115 s%ve =0.d0 00116 fc =2.d0*par%wearth*sin(par%lat) 00117 00118 call bedroughness_init(s,par) ! note, this is not yet designed for initialisation 00119 ! on sglobal, so don't call from initialize.F90 00120 endif 00121 00122 ! Super fast 1D 00123 if (s%ny==0) then 00124 j1 = 1 00125 else 00126 j1 = 2 00127 endif 00128 00129 ! Add vertical discharges 00130 call discharge_boundary_v(s,par) 00131 00132 ! 00133 ! s%zs=s%zs*s%wetz 00134 ! Water level slopes 00135 do j=1,s%ny+1 00136 do i=2,s%nx 00137 s%dzsdx(i,j)=(s%zs(i+1,j)+s%ph(i+1,j)-s%zs(i,j)-s%ph(i,j))/s%dsu(i,j) 00138 end do 00139 end do 00140 ! do j=2,ny 00141 do j=1,s%ny ! Dano need to get correct slope on boundary s%y=0 00142 do i=1,s%nx+1 00143 s%dzsdy(i,j)=(s%zs(i,j+1)+s%ph(i,j+1)-s%zs(i,j)-s%ph(i,j))/s%dnv(i,j) 00144 end do 00145 end do 00146 ! 00147 ! Compute velocity gradients for viscosity terms. 00148 ! Robert: Check whether should be same gradients as advection terms? 00149 do j=j1,max(s%ny,1) 00150 do i=2,s%nx+1 00151 dudx(i,j) = (s%uu(i,j)-s%uu(i-1,j))/s%dsz(i,j) 00152 enddo 00153 enddo 00154 ! wwvv: added: xmpi_istop 00155 if (xmpi_istop) then 00156 dudx(1,:) = 0.d0 ! Robert: by defintion of Neumann boundary 00157 endif 00158 if (s%ny>2) then 00159 do j=2,s%ny+1 00160 do i=1,s%nx+1 00161 dvdy(i,j) = (s%vv(i,j)-s%vv(i,j-1))/s%dnz(i,j) 00162 enddo 00163 enddo 00164 ! wwvv: added: xmpi_isleft 00165 if (xmpi_isleft) then 00166 dvdy(:,1) = 0.d0 ! Robert: by defintion of Neumann boundary 00167 endif 00168 else 00169 dvdy = 0.d0 ! Robert: by definition of 1D model 00170 endif 00171 00172 ! Update bed roughness coefficient 00173 call bedroughness_update(s,par) 00174 !cf = cfu : Robert: cf is not used anymore 00175 00176 ! 00177 ! Advection x-direction 00178 ! 00179 do j=j1,max(s%ny,1) 00180 do i=2,s%nx 00181 s%ududx(i,j) = 0.d0 00182 qin = .5d0*(s%qx(i,j)+s%qx(i-1,j)) 00183 if (qin>0) then 00184 dalfa = s%alfau(i,j)-s%alfau(i-1,j) 00185 uin = s%uu(i-1,j)*cos(dalfa) + s%vu(i-1,j)*sin(dalfa) 00186 if ((s%uu(i,j)-s%uu(i-1,j))>par%eps_sd) then 00187 ! Conservation of energy head 00188 s%ududx(i,j) = s%ududx(i,j) + 0.5d0*(s%uu(i-1,j)+s%uu(i,j))*(s%uu(i,j)-uin)*s%dnz(i,j)*s%dsdnui(i,j) 00189 else 00190 ! Conservation of momentum 00191 s%ududx(i,j) = s%ududx(i,j) + qin/s%hum(i,j) *(s%uu(i,j)-uin)*s%dnz(i,j)*s%dsdnui(i,j) 00192 endif 00193 endif 00194 qin = -.5d0*(s%qx(i,j)+s%qx(i+1,j)) 00195 if (qin>0) then 00196 dalfa = s%alfau(i,j)-s%alfau(i+1,j) 00197 uin = s%uu(i+1,j)*cos(dalfa) + s%vu(i+1,j)*sin(dalfa) 00198 if ((s%uu(i+1,j)-s%uu(i,j))>par%eps_sd) then 00199 ! Conservation of energy head 00200 s%ududx(i,j) = s%ududx(i,j) - 0.5d0*(s%uu(i+1,j)+s%uu(i,j))*(s%uu(i,j)-uin)*s%dnz(i+1,j)*s%dsdnui(i,j) 00201 else 00202 ! Conservation of momentum 00203 s%ududx(i,j) = s%ududx(i,j) + qin/s%hum(i,j) *(s%uu(i,j)-uin)*s%dnz(i+1,j)*s%dsdnui(i,j) 00204 endif 00205 endif 00206 end do 00207 end do 00208 do j=2,s%ny 00209 do i=1,s%nx 00210 s%vdudy(i,j) = 0.d0 00211 qin = .5d0*(s%qy(i,j-1)+s%qy(i+1,j-1)) 00212 if (qin>0) then 00213 dalfa = s%alfau(i,j)-s%alfau(i,j-1) 00214 uin = s%uu(i,j-1)*cos(dalfa) + s%vu(i,j-1)*sin(dalfa) 00215 if ((s%vv(i,j)-s%vv(i,j-1))>par%eps_sd) then 00216 ! Conservation of energy head 00217 s%vdudy(i,j) = s%vdudy(i,j) + 0.5d0*(s%vv(i,j-1)+s%vv(i+1,j-1))*(s%uu(i,j)-uin)*s%dsc(i,j-1)*s%dsdnui(i,j) 00218 else 00219 ! Conservation of momentum 00220 s%vdudy(i,j) = s%vdudy(i,j) + qin/s%hum(i,j) *(s%uu(i,j)-uin)*s%dsc(i,j-1)*s%dsdnui(i,j) 00221 endif 00222 endif 00223 qin = -.5d0*(s%qy(i,j)+s%qy(i+1,j)) 00224 if (qin>0) then 00225 dalfa = s%alfau(i,j)-s%alfau(i,j+1) 00226 uin = s%uu(i,j+1)*cos(dalfa) + s%vu(i,j+1)*sin(dalfa) 00227 if ((s%vv(i,j+1)-s%vv(i,j))>par%eps_sd) then 00228 ! Conservation of energy head 00229 s%vdudy(i,j) = s%vdudy(i,j) - 0.5d0*(s%vv(i,j)+s%vv(i+1,j))*(s%uu(i,j)-uin)*s%dsc(i,j)*s%dsdnui(i,j) 00230 else 00231 ! Conservation of momentum 00232 s%vdudy(i,j) = s%vdudy(i,j) + qin/s%hum(i,j) *(s%uu(i,j)-uin)*s%dsc(i,j)*s%dsdnui(i,j) 00233 endif 00234 endif 00235 end do 00236 end do 00237 ! 00238 ! Advection y-direction 00239 00240 if (xmpi_isright .and. s%ny>0) then ! no such condition needed for _isleft, because s%vdvdy(:,1) not needed in mpi subdomain 00241 jmax = s%ny-1 00242 elseif (xmpi_isright .and. s%ny==0) then 00243 jmax = 1 00244 else 00245 jmax = s%ny 00246 endif 00247 s%vdvdy=0.d0 00248 ! calculate true vdvdy up to ny in central domains and up to ny-1 on isright 00249 do j=2,jmax 00250 do i=2,s%nx 00251 qin = .5d0*(s%qy(i,j)+s%qy(i,j-1)) 00252 if (qin>0) then 00253 dalfa = s%alfav(i,j)-s%alfav(i,j-1) 00254 vin = s%vv(i,j-1)*cos(dalfa) - s%uv(i,j-1)*sin(dalfa) 00255 if ((s%vv(i,j)-s%vv(i,j-1))>par%eps_sd) then 00256 ! Conservation of energy head 00257 s%vdvdy(i,j) = s%vdvdy(i,j) + 0.5d0*(s%vv(i,j-1)+s%vv(i,j))*(s%vv(i,j)-vin)*s%dsz(i,j)*s%dsdnvi(i,j) 00258 else 00259 ! Conservation of momentum 00260 s%vdvdy(i,j) = s%vdvdy(i,j) + qin/s%hvm(i,j) *(s%vv(i,j)-vin)*s%dsz(i,j)*s%dsdnvi(i,j) 00261 endif 00262 endif 00263 qin = -.5d0*(s%qy(i,j)+s%qy(i,j+1)) 00264 if (qin>0) then 00265 dalfa = s%alfav(i,j)-s%alfav(i,j+1) 00266 vin = s%vv(i,j+1)*cos(dalfa) - s%uv(i,j+1)*sin(dalfa) 00267 if ((s%vv(i,j+1)-s%vv(i,j))>par%eps_sd) then 00268 ! Conservation of energy head 00269 s%vdvdy(i,j) = s%vdvdy(i,j) - 0.5d0*(s%vv(i,j+1)+s%vv(i,j))*(s%vv(i,j)-vin)*s%dsz(i,j+1)*s%dsdnvi(i,j) 00270 else 00271 ! Conservation of momentum 00272 s%vdvdy(i,j) = s%vdvdy(i,j) + qin/s%hvm(i,j) *(s%vv(i,j)-vin)*s%dsz(i,j+1)*s%dsdnvi(i,j) 00273 endif 00274 endif 00275 enddo 00276 enddo 00277 if (s%ny>0) then 00278 ! Global boundary conditions for vdvdy(:,1) and vdvdy(:,ny), global vdvdy(:,ny+1) not needed anywhere 00279 if (xmpi_isleft) then 00280 ! (vv(:,1)-vv(:,0))/dy == 0 so only second part of the vdvdy equation: 00281 do i=2,s%nx 00282 qin = -.5d0*(s%qy(i,1)+s%qy(i,2)) 00283 if (qin>0) then 00284 dalfa = s%alfav(i,1)-s%alfav(i,2) 00285 vin = s%vv(i,2)*cos(dalfa) - s%uv(i,2)*sin(dalfa) 00286 s%vdvdy(i,1) = s%vdvdy(i,1) + qin*(s%vv(i,1)-vin)*s%dsz(i,2)/s%hvm(i,1)*s%dsdnvi(i,1) 00287 endif 00288 enddo 00289 endif 00290 if (xmpi_isright) then 00291 ! (vv(:,ny+1)-vv(:,ny))/dy == 0 so only first part of the vdvdy equation: 00292 do i=2,s%nx 00293 qin = .5d0*(s%qy(i,s%ny)+s%qy(i,s%ny-1)) 00294 if (qin>0) then 00295 dalfa = s%alfav(i,s%ny)-s%alfav(i,s%ny-1) 00296 vin = s%vv(i,s%ny-1)*cos(dalfa) - s%uv(i,s%ny-1)*sin(dalfa) 00297 s%vdvdy(i,s%ny) = s%vdvdy(i,s%ny) + qin*(s%vv(i,s%ny)-vin)*s%dsz(i,s%ny)/s%hvm(i,s%ny)*s%dsdnvi(i,s%ny) 00298 endif 00299 enddo 00300 endif 00301 endif 00302 00303 s%udvdx=0.d0 00304 if (s%ny>0) then 00305 ! Robert: udvdx not usually needed at j = 1 00306 do j=1,s%ny !1,s%ny instead of 2,s%ny 00307 do i=2,s%nx 00308 qin = .5d0*(s%qx(i-1,j)+s%qx(i-1,j+1)) 00309 if (qin>0) then 00310 dalfa = s%alfav(i,j)-s%alfav(i-1,j) 00311 vin = s%vv(i-1,j)*cos(dalfa) - s%uv(i-1,j)*sin(dalfa) 00312 if ((s%uu(i,j)-s%uu(i-1,j))>par%eps_sd) then 00313 ! Conservation of energy head 00314 s%udvdx(i,j) = s%udvdx(i,j) + 0.5d0*(s%uu(i-1,j)+s%uu(i-1,j+1))*(s%vv(i,j)-vin)*s%dnc(i-1,j)*s%dsdnvi(i,j) 00315 else 00316 ! Conservation of momentum 00317 s%udvdx(i,j) = s%udvdx(i,j) + qin/s%hvm(i,j) *(s%vv(i,j)-vin)*s%dnc(i-1,j)*s%dsdnvi(i,j) 00318 endif 00319 endif 00320 qin = -.5d0*(s%qx(i,j)+s%qx(i,j+1)) 00321 if (qin>0) then 00322 dalfa = s%alfav(i,j)-s%alfav(i+1,j) 00323 vin = s%vv(i+1,j)*cos(dalfa) - s%uv(i+1,j)*sin(dalfa) 00324 if ((s%uu(i+1,j)-s%uu(i,j))>par%eps_sd) then 00325 ! Conservation of energy head 00326 s%udvdx(i,j) = s%udvdx(i,j) - 0.5d0*(s%uu(i,j)+s%uu(i,j+1))*(s%vv(i,j)-vin)*s%dnc(i,j)*s%dsdnvi(i,j) 00327 else 00328 ! Conservation of momentum 00329 s%udvdx(i,j) = s%udvdx(i,j) + qin/s%hvm(i,j) *(s%vv(i,j)-vin)*s%dnc(i,j)*s%dsdnvi(i,j) 00330 endif 00331 endif 00332 end do 00333 end do 00334 else 00335 do i=2,s%nx 00336 qin = s%qx(i-1,1) 00337 if (qin>0) then 00338 dalfa = s%alfav(i,1)-s%alfav(i-1,1) 00339 vin = s%vv(i-1,1)*cos(dalfa) - s%uv(i-1,1)*sin(dalfa) 00340 s%udvdx(i,1) = s%udvdx(i,1) + qin*(s%vv(i,1)-vin)*s%dnc(i-1,1)/s%hvm(i,1)*s%dsdnvi(i,1) 00341 endif 00342 qin = -s%qx(i,1) 00343 if (qin>0) then 00344 dalfa = s%alfav(i,1)-s%alfav(i+1,1) 00345 vin = s%vv(i+1,1)*cos(dalfa) - s%uv(i+1,1)*sin(dalfa) 00346 s%udvdx(i,1) = s%udvdx(i,1) + qin*(s%vv(i,1)-vin)*s%dnc(i,1)/s%hvm(i,1)*s%dsdnvi(i,1) 00347 endif 00348 enddo 00349 endif 00350 ! 00351 ! 00352 ! There is a possibility to turn viscosity off 00353 ! 00354 if (par%viscosity==0) then ! Dano: put it here so it actualy saves a lot of computation time 00355 s%nuh = 0.d0 00356 s%viscu = 0.d0 00357 else 00358 00359 ! Jaap: Slightly changes approach; 1) background viscosity is user defined or obtained from Smagorinsky, 2) nuh = max(nuh,roller induced viscosity) 00360 if (par%smag == 1) then 00361 !Use smagorinsky subgrid model 00362 call visc_smagorinsky(s,par) 00363 else 00364 s%nuh = par%nuh 00365 endif 00366 ! 00367 ! Add viscosity for wave breaking effects 00368 if (par%swave == 1) then 00369 do j=j1,max(s%ny,1) 00370 do i=2,s%nx 00371 s%nuh(i,j) = max(s%nuh(i,j),par%nuhfac*s%hh(i,j)*(s%DR(i,j)/par%rho)**(1.0d0/3.0d0)) ! Ad: change to max 00372 end do 00373 end do 00374 elseif (par%swave==0 .and. par%wavemodel==WAVEMODEL_NONH) then 00375 select case (par%nhbreaker) 00376 case (1) 00377 where (s%breaking/=0) 00378 s%nuh = par%breakviscfac*s%nuh 00379 endwhere 00380 case (2) 00381 where (s%breaking==1) 00382 s%nuh = s%nuh + (par%nuh*par%breakvisclen*s%hh)**2*sqrt(dudx**2+dvdy**2) 00383 endwhere 00384 case (3) 00385 ! Ad en Arnold: Battjes 1975 formulation to smoothen front of lf wave bore in the swash 00386 ! compute (long) wave turbulence due to breaking 00387 !nuh = max(par%nuh, par%avis*hloc*sqrt(kturb)) 00388 case default 00389 ! do nothing to increase viscosity 00390 end select 00391 endif 00392 00393 ! 00394 do j=jmin_uu,jmax_uu 00395 do i=imin_uu,imax_uu 00396 00397 dudx1 = s%nuh(i+1,j)*s%hh(i+1,j)*(s%uu(i+1,j)-s%uu(i,j))/s%dsz(i+1,j) 00398 dudx2 = s%nuh(i,j) *s%hh(i ,j)*(s%uu(i,j)-s%uu(i-1,j))/s%dsz(i,j) 00399 s%viscu(i,j) = (1.0d0/s%hum(i,j))*( 2*(dudx1-dudx2)/(s%dsz(i,j)+s%dsz(i+1,j)) ) 00400 end do 00401 end do 00402 00403 if (par%smag == 1) then 00404 ! 00405 ! For non constant eddy viscosity the stress terms read: 00406 ! 00407 ! d d d d d 00408 ! -- [ 2* mu -- (U) ] + --[ mu --(U) +mu --(V) ] 00409 ! dx dx dy dy dx 00410 ! 00411 ! d d 00412 ! Only when -- [mu] = --[mu] = 0 we have (for incompressible flow) 00413 ! dx dy 00414 ! 00415 ! 2 2 00416 ! d d 00417 ! mu --2 [U] + mu --2 [U] 00418 ! dx dy 00419 ! 00420 s%viscu = 2.0d0*s%viscu 00421 endif 00422 00423 do j=jmin_uu,jmax_uu 00424 jp1=min(j+1,1) 00425 jm1=max(j-1,1) 00426 do i=imin_uu,imax_uu 00427 !Nuh is defined at eta points, interpolate from four surrounding points 00428 nuh1 = .25d0*(s%nuh(i,j)+s%nuh(i+1,j)+s%nuh(i+1,jp1)+s%nuh(i,jp1)) 00429 nuh2 = .25d0*(s%nuh(i,j)+s%nuh(i+1,j)+s%nuh(i+1,jm1)+s%nuh(i,jm1)) 00430 00431 dudy1 = nuh1 *.5d0*(s%hvm(i,j )+s%hvm(i+1,j ))*(s%uu(i,jp1)-s%uu(i,j))/s%dnc(i,j) 00432 dudy2 = nuh2 *.5d0*(s%hvm(i,jm1)+s%hvm(i+1,jm1))*(s%uu(i,j)-s%uu(i,jm1))/s%dnc(i,jm1) 00433 s%viscu(i,j) = s%viscu(i,j) + (1.0d0/s%hum(i,j))* & 00434 ( 2.0d0*(dudy1-dudy2)/(s%dnc(i,j)+s%dnc(i,jm1)) )*s%wetu(i,jp1)*s%wetu(i,jm1) 00435 end do 00436 end do 00437 00438 if (par%smag == 1) then 00439 do j=jmin_uu,jmax_uu 00440 jp1=min(j+1,1) 00441 jm1=max(j-1,1) 00442 do i=imin_uu,imax_uu 00443 !Nuh is defined at eta points, interpolate from four surrounding points 00444 nuh1 = .25d0*(s%nuh(i,j)+s%nuh(i+1,j)+s%nuh(i+1,jp1)+s%nuh(i,jp1)) 00445 nuh2 = .25d0*(s%nuh(i,j)+s%nuh(i+1,j)+s%nuh(i+1,jm1)+s%nuh(i,jm1)) 00446 00447 dvdx1 = nuh1*.5d0*(s%hvm(i,j )+s%hvm(i+1,j ))*(s%vv(i+1,j )-s%vv(i,j ))/s%dsc(i,j) 00448 dvdx2 = nuh2*.5d0*(s%hvm(i,jm1)+s%hvm(i+1,jm1))*(s%vv(i+1,jm1)-s%vv(i,jm1))/s%dsc(i,jm1) 00449 s%viscu(i,j) = s%viscu(i,j) + (1.d0/s%hum(i,j))*(dvdx1-dvdx2)/s%dnz(i,j)*real(s%wetv(i+1,j) & * s%wetv(i,j)*s%wetv(i+1,jm1)*s%wetv(i,jm1),8) 00450 enddo 00451 enddo 00452 endif !smag ==1 and s%ny>0 00453 ! 00454 ! Viscosity y-direction 00455 ! 00456 s%viscv =0.d0 00457 do j=jmin_vv,jmax_vv 00458 jp1=min(j+1,1) 00459 jm1=max(j-1,1) 00460 do i=imin_vv,imax_vv 00461 dvdy1 = s%nuh(i,jp1)*s%hh(i,jp1)*(s%vv(i,jp1)-s%vv(i,j))/s%dnz(i,jp1) 00462 dvdy2 = s%nuh(i,j) *s%hh(i,j )*(s%vv(i,j)-s%vv(i,jm1))/s%dnz(i,j) 00463 s%viscv(i,j) = (1.0d0/s%hvm(i,j))* 2*(dvdy1-dvdy2)/(s%dnz(i,j)+s%dnz(i,jp1))*s%wetv(i,jp1)*s%wetv(i,jm1) 00464 end do 00465 end do 00466 ! Robert: global boundary at (:,1) edge 00467 if (s%ny>0) then 00468 if (xmpi_isleft) then 00469 s%viscv(:,1) = s%viscv(:,2) 00470 endif 00471 if (xmpi_isright) then 00472 s%viscv(:,s%ny) = s%viscv(:,s%ny-1) 00473 endif 00474 endif 00475 ! 00476 ! Viscosity 00477 if (par%smag == 1) then 00478 s%viscv = 2.0d0*s%viscv 00479 endif 00480 ! 00481 s%nuh = par%nuhv*s%nuh !Robert en Ap: increase s%nuh interaction in d2v/dx2 00482 do j=jmin_vv,jmax_vv 00483 jp1=min(j+1,1) 00484 jm1=max(j-1,1) 00485 do i=imin_vv,imax_vv 00486 !Nuh is defined at eta points, interpolate from four surrounding points 00487 nuh1 = .25d0*(s%nuh(i,j)+s%nuh(i+1,j)+s%nuh(i+1,jp1)+s%nuh(i,jp1)) 00488 nuh2 = .25d0*(s%nuh(i,j)+s%nuh(i-1,j)+s%nuh(i-1,jp1)+s%nuh(i,jp1)) 00489 00490 dvdx1 = nuh1*.5d0*(s%hum(i ,j)+s%hum(i ,jp1))*(s%vv(i+1,j)-s%vv(i,j))/s%dsc(i,j) 00491 dvdx2 = nuh2*.5d0*(s%hum(i-1,j)+s%hum(i-1,jp1))*(s%vv(i,j)-s%vv(i-1,j))/s%dsc(i-1,j) 00492 s%viscv(i,j) = s%viscv(i,j) + (1.0d0/s%hvm(i,j))*( 2*(dvdx1-dvdx2)/(s%dsc(i-1,j)+s%dsc(i,j)) ) & 00493 *s%wetv(i+1,j)*s%wetv(i-1,j) 00494 end do 00495 end do 00496 ! 00497 if (par%smag == 1) then 00498 do j=jmin_vv,jmax_vv 00499 jp1=min(j+1,1) 00500 jm1=max(j-1,1) 00501 do i=imin_vv,imax_vv 00502 !Nuh is defined at eta points, interpolate from four surrounding points 00503 nuh1 = .25d0*(s%nuh(i,j)+s%nuh(i+1,j)+s%nuh(i+1,jp1)+s%nuh(i,jp1)) 00504 nuh2 = .25d0*(s%nuh(i,j)+s%nuh(i-1,j)+s%nuh(i-1,jp1)+s%nuh(i,jp1)) 00505 00506 dudy1 = nuh1 *.5d0*(s%hum(i ,j)+s%hum(i ,jp1))*(s%uu(i,jp1 )-s%uu(i,j ))/s%dnc(i,j) 00507 dudy2 = nuh2 *.5d0*(s%hum(i-1,j)+s%hum(i-1,jp1))*(s%uu(i-1,jp1)-s%uu(i-1,j))/s%dnc(i-1,j) 00508 00509 s%viscv(i,j) = s%viscv(i,j) + (1.d0/s%hvm(i,j))*(dudy1-dudy2)/s%dsz(i,j) & 00510 * real(s%wetu(i,jp1)*s%wetu(i,j)*s%wetu(i-1,jp1)*s%wetv(i-1,j),8) 00511 enddo 00512 enddo 00513 endif 00514 00515 endif !par%viscosity==0 00516 ! 00517 ! Bed friction term x 00518 ! 00519 where (s%wetu==1) 00520 s%taubx=s%cfu*par%rho*s%ueu*sqrt((1.16d0*s%urms)**2+s%vmageu**2) !Ruessink et al, 2001 00521 s%taubx = s%taubx + s%taubx_add ! Account for infiltration etc. effects 00522 elsewhere 00523 s%taubx = 0.d0 00524 endwhere 00525 ! 00526 ! 00527 ! limit bed friction term to accelerate less than realistic value 00528 where (abs(s%taubx)>100*par%g*par%rho*s%hu) 00529 s%taubx = sign(1.d0,s%taubx)*100*par%g*par%rho*s%hu 00530 endwhere 00531 ! 00532 ! Bed friction term y 00533 ! 00534 where (s%wetv==1) 00535 s%tauby=s%cfv*par%rho*s%vev*sqrt((1.16d0*s%urms)**2+s%vmagev**2) !Ruessink et al, 2001 00536 s%tauby = s%tauby + s%tauby_add 00537 elsewhere 00538 s%tauby = 0.d0 00539 endwhere 00540 where (abs(s%tauby)>100*par%g*par%rho*s%hv) 00541 s%tauby = sign(1.d0,s%tauby)*100*par%g*par%rho*s%hv 00542 endwhere 00543 ! 00544 ! Explicit Euler step momentum u-direction 00545 ! 00546 do j=jmin_uu,jmax_uu 00547 do i=imin_uu,imax_uu 00548 if(s%wetu(i,j)==1) then 00549 s%uu(i,j)=s%uu(i,j)-par%dt*(s%ududx(i,j)+s%vdudy(i,j)-s%viscu(i,j) & !Ap,Robert,Jaap 00550 + par%g*s%dzsdx(i,j) & 00551 + s%taubx(i,j)/(par%rho*s%hu(i,j)) & ! Dano: changed s%hum to s%hu NOT s%cf volume approach 00552 + s%Fvegu(i,j)/(par%rho*s%hu(i,j)) & 00553 - par%lwave*s%Fx(i,j)/(par%rho*s%hum(i,j)) & 00554 - fc*s%vu(i,j) & 00555 - par%rhoa*par%Cd*s%windsu(i,j)*sqrt(s%windsu(i,j)**2+s%windnv(i,j)**2)/(par%rho*s%hum(i,j))) ! Kees: wind correction 00556 else 00557 s%uu(i,j)=0.0d0 00558 end if 00559 end do 00560 end do 00561 ! Lateral boundary conditions for uu 00562 if (s%ny>0) then 00563 if (xmpi_isleft) then !Dano/Robert only on outer boundary 00564 s%uu(1:s%nx+1,1)=s%uu(1:s%nx+1,2) ! RJ: can also be done after continuity but more appropriate here 00565 endif 00566 ! Lateral boundary at y=ny*dy 00567 if (xmpi_isright) then !Dano/Robert only at outer boundary 00568 s%uu(1:s%nx+1,s%ny+1)=s%uu(1:s%nx+1,s%ny) ! RJ: can also be done after continuity but more appropriate here 00569 endif 00570 endif 00571 ! 00572 ! Explicit Euler step momentum v-direction 00573 ! 00574 do j=jmin_vv,jmax_vv 00575 do i=imin_vv,imax_vv 00576 if(s%wetv(i,j)==1) then 00577 ! Robert: ensure taubx always has the same sign as uu (always decelerates) 00578 ! Dano: I don't agree 00579 s%vv(i,j)=s%vv(i,j)-par%dt*(s%udvdx(i,j)+s%vdvdy(i,j)-s%viscv(i,j)& !Ap,Robert,Jaap 00580 + par%g*s%dzsdy(i,j)& 00581 + s%tauby(i,j)/(par%rho*s%hv(i,j)) & ! Dano: s%hv instead of s%hvm, NOT cf volume approach 00582 + s%Fvegv(i,j)/(par%rho*s%hv(i,j)) & 00583 - par%lwave*s%Fy(i,j)/(par%rho*s%hvm(i,j)) & 00584 + fc*s%uv(i,j) & 00585 - par%rhoa*par%Cd*s%windnv(i,j)*sqrt(s%windsu(i,j)**2+s%windnv(i,j)**2)/(par%rho*s%hvm(i,j))) ! Kees: wind correction 00586 else 00587 s%vv(i,j)=0.0d0 00588 end if 00589 end do 00590 end do 00591 ! Robert: Boundary conditions along the global boundaries 00592 ! Function flow_lat_bc located in boundaryconditions.F90 00593 ! function call takes care of 1D vs 2D models and boundary condition types 00594 if (s%ny>0) then 00595 if (xmpi_isleft) then 00596 s%vv(:,1)=flow_lat_bc(s,par,par%right,1,2,s%udvdx(:,1),s%vdvdy(:,1),s%viscv(:,1)) 00597 endif 00598 if (xmpi_isright) then 00599 s%vv(:,s%ny)=flow_lat_bc(s,par,par%left,s%ny,s%ny-1,s%udvdx(:,s%ny),s%vdvdy(:,s%ny),s%viscv(:,s%ny)) 00600 endif 00601 endif 00602 00603 #ifdef USEMPI 00604 call xmpi_shift_ee(s%uu) 00605 call xmpi_shift_ee(s%vv) 00606 #endif 00607 00608 if (par%wavemodel==WAVEMODEL_NONH) then 00609 !Do explicit predictor step with pressure 00610 call nonh_cor(s,par,0,uu_old,vv_old) 00611 00612 #ifdef USEMPI 00613 call xmpi_shift_ee(s%uu) 00614 call xmpi_shift_ee(s%vv) 00615 #endif 00616 end if 00617 00618 if (par%secorder==1) then 00619 !Call second order correction to the advection 00620 call flow_secondorder_advUV(s,par,uu_old,vv_old) 00621 #ifdef USEMPI 00622 call xmpi_shift_ee(s%uu) 00623 call xmpi_shift_ee(s%vv) 00624 #endif 00625 end if 00626 00627 if (par%wavemodel==WAVEMODEL_NONH) then 00628 ! do non-hydrostatic pressure compensation to solve short waves 00629 call nonh_cor(s,par,1,uu_old,vv_old) 00630 ! note: MPI shift in subroutine nonh_cor 00631 end if 00632 00633 ! update hu en hv for continuity 00634 do j=1,s%ny+1 00635 do i=1,s%nx+1 00636 ! Water depth in u-points do continuity equation: upwind 00637 if (s%uu(i,j)>par%umin) then 00638 if (par%oldhu == 1) then 00639 s%hu(i,j)=s%hh(i,j) 00640 else 00641 s%hu(i,j)=s%zs(i,j)-max(s%zb(i,j),s%zb(min(s%nx,i)+1,j)) 00642 endif 00643 elseif (s%uu(i,j)<-par%umin) then 00644 if (par%oldhu == 1) then 00645 s%hu(i,j)=s%hh(min(s%nx,i)+1,j) 00646 else 00647 s%hu(i,j)=s%zs(min(s%nx,i)+1,j)-max(s%zb(i,j),s%zb(min(s%nx,i)+1,j)) 00648 endif 00649 else 00650 s%hu(i,j)=max(max(s%zs(i,j),s%zs(min(s%nx,i)+1,j))-max(s%zb(i,j),s%zb(min(s%nx,i)+1,j)),par%eps) 00651 end if 00652 end do 00653 end do 00654 00655 s%hu = max(s%hu,0.d0) 00656 00657 do j=1,s%ny+1 00658 do i=1,s%nx+1 00659 ! Water depth in v-points do continuity equation: upwind 00660 if (s%vv(i,j)>par%umin) then 00661 if (par%oldhu == 1) then 00662 s%hv(i,j)=s%hh(i,j) 00663 else 00664 s%hv(i,j)=s%zs(i,j)-max(s%zb(i,j),s%zb(i,min(s%ny,j)+1)) 00665 endif 00666 elseif (s%vv(i,j)<-par%umin) then 00667 if (par%oldhu == 1) then 00668 s%hv(i,j)=s%hh(i,min(s%ny,j)+1) 00669 else 00670 s%hv(i,j)=s%zs(i,min(s%ny,j)+1)-max(s%zb(i,j),s%zb(i,min(s%ny,j)+1)) 00671 endif 00672 else 00673 s%hv(i,j)=max(max(s%zs(i,j),s%zs(i,min(s%ny,j)+1))-max(s%zb(i,j),s%zb(i,min(s%ny,j)+1)),par%eps) 00674 end if 00675 end do 00676 end do 00677 s%hv = max(s%hv,0.d0) 00678 00679 if (par%secorder == 1) then 00680 ! 00681 ! Correct the waterdepths in U/V points using second-order limited expressions 00682 ! 00683 call flow_secondorder_huhv(s,par) 00684 ! 00685 end if 00686 00687 ! Flux in u-point 00688 s%qx=s%uu*s%hu 00689 ! Flux in v-points 00690 ! first column of qy is used later, and it is defined in the loop above 00691 ! no communication necessary at this point 00692 s%qy=s%vv*s%hv 00693 ! 00694 ! Add horizontal discharges 00695 ! 00696 call discharge_boundary_h(s,par) 00697 ! 00698 ! Add rainfall 00699 if (par%rainfall==1) then 00700 call rainfall_update(s,par) 00701 endif 00702 ! 00703 ! Update water level using continuity eq. 00704 ! 00705 if (s%ny>0) then 00706 do j=jmin_zs,jmax_zs 00707 do i=imin_zs,imax_zs 00708 s%dzsdt(i,j) = (-1.d0)*( s%qx(i,j)*s%dnu(i,j)-s%qx(i-1,j)*s%dnu(i-1,j) & 00709 + s%qy(i,j)*s%dsv(i,j)-s%qy(i,j-1)*s%dsv(i,j-1) )*s%dsdnzi(i,j) & 00710 - s%infil(i,j) + s%rainfallrate(i,j) 00711 end do 00712 end do 00713 s%zs(imin_zs:imax_zs,jmin_zs:jmax_zs) = s%zs(imin_zs:imax_zs,jmin_zs:jmax_zs) & 00714 +s%dzsdt(imin_zs:imax_zs,jmin_zs:jmax_zs)*par%dt 00715 else 00716 j=1 00717 do i=imin_zs,imax_zs 00718 s%dzsdt(i,j) = (-1.d0)*( s%qx(i,j)*s%dnu(i,j)-s%qx(i-1,j)*s%dnu(i-1,j) )*s%dsdnzi(i,j) & 00719 - s%infil(i,j) + s%rainfallrate(i,j) 00720 end do 00721 s%zs(imin_zs:imax_zs,1) = s%zs(imin_zs:imax_zs,1)+s%dzsdt(imin_zs:imax_zs,1)*par%dt 00722 endif !s%ny>0 00723 00724 if (par%secorder == 1) then 00725 ! 00726 ! P.B. Smit: The second order McCormack correction to the cont. eq. introduced 00727 ! minor damping. For this reason I removed it (Nov. 2014). 00728 ! call flow_secondorder_con(s,par,zs_old) 00729 endif 00730 00731 ! wwvv zs, uu, vv have to be communicated now, because they are used later on 00732 #ifdef USEMPI 00733 call xmpi_shift_ee(s%zs) 00734 #endif 00735 00736 00737 if (par%secorder == 1 .or. par%wavemodel==WAVEMODEL_NONH) then 00738 vv_old = s%vv 00739 uu_old = s%uu 00740 zs_old = s%zs 00741 endif 00742 00743 ! 00744 ! U and V in cell centre; do output and sediment stirring 00745 ! 00746 s%u(2:s%nx,:)=0.5d0*(s%uu(1:s%nx-1,:)+s%uu(2:s%nx,:)) 00747 ! offshore boundary 00748 if(xmpi_istop) then 00749 s%u(1,:)=s%uu(1,:) 00750 endif 00751 if(xmpi_isbot) then 00752 s%u(s%nx+1,:)=s%u(s%nx,:) 00753 endif 00754 00755 if (s%ny>0) then 00756 s%v(:,2:s%ny)=0.5d0*(s%vv(:,1:s%ny-1)+s%vv(:,2:s%ny)) 00757 if(xmpi_isleft) then 00758 s%v(:,1)=s%vv(:,1) 00759 endif 00760 if(xmpi_isright) then 00761 s%v(:,s%ny+1)=s%v(:,s%ny) ! bas: need this for calculation of s%ee in wci routine 00762 endif 00763 !Ap 00764 s%v(s%nx+1,:)=s%v(s%nx,:) 00765 else ! Dano 00766 s%v=s%vv 00767 endif !s%ny>0 00768 00769 ! Robert + Jaap: compute derivatives of u and v 00770 ! 00771 00772 sinthm = sin(s%thetamean-s%alfaz) 00773 costhm = cos(s%thetamean-s%alfaz) 00774 00775 ! V-velocities at u-points 00776 if (s%ny>0) then 00777 s%vu(1:s%nx,2:s%ny)= 0.25d0*(s%vv(1:s%nx,1:s%ny-1)+s%vv(1:s%nx,2:s%ny)+ & 00778 s%vv(2:s%nx+1,1:s%ny-1)+s%vv(2:s%nx+1,2:s%ny)) 00779 ! how about boundaries? 00780 if(xmpi_isleft) then 00781 s%vu(:,1) = s%vu(:,2) 00782 endif 00783 if(xmpi_isright) then 00784 s%vu(:,s%ny+1) = s%vu(:,s%ny) 00785 endif 00786 else 00787 s%vu(1:s%nx,1)= 0.5d0*(s%vv(1:s%nx,1)+s%vv(2:s%nx+1,1)) 00788 endif !s%ny>0 00789 ! wwvv fill in vu(:1) and vu(:ny+1) for non-left and non-right processes 00790 ! and vu(nx+1,:) 00791 s%vu=s%vu*s%wetu 00792 ! V-stokes velocities at U point 00793 if (s%ny>0) then 00794 vsu(1:s%nx,2:s%ny)=0.5d0*(s%ust(1:s%nx,2:s%ny)*sinthm(1:s%nx,2:s%ny)+ & 00795 s%ust(2:s%nx+1,2:s%ny)*sinthm(2:s%nx+1,2:s%ny)) 00796 if(xmpi_isleft) then 00797 vsu(:,1)=vsu(:,2) 00798 endif 00799 if(xmpi_isright) then 00800 vsu(:,s%ny+1) = vsu(:,s%ny) 00801 endif 00802 else 00803 vsu(1:s%nx,1)=0.5d0*(s%ust(1:s%nx,1)*sinthm(1:s%nx,1)+ & 00804 s%ust(2:s%nx+1,1)*sinthm(2:s%nx+1,1)) 00805 endif !s%ny>0 00806 ! wwvv same for vsu 00807 vsu = vsu*s%wetu 00808 ! U-stokes velocities at U point 00809 if (s%ny>0) then 00810 usu(1:s%nx,2:s%ny)=0.5d0*(s%ust(1:s%nx ,2:s%ny)*costhm(1:s%nx ,2:s%ny)+ & 00811 s%ust(2:s%nx+1,2:s%ny)*costhm(2:s%nx+1,2:s%ny)) 00812 if(xmpi_isleft) then 00813 usu(:,1)=usu(:,2) 00814 endif 00815 if(xmpi_isright) then 00816 usu(:,s%ny+1)=usu(:,s%ny) 00817 endif 00818 else 00819 usu(1:s%nx,1)=0.5d0*(s%ust(1:s%nx,1)*costhm(1:s%nx,1)+ & 00820 s%ust(2:s%nx+1,1)*costhm(2:s%nx+1,1)) 00821 endif !s%ny>0 00822 ! wwvv same for usu 00823 usu=usu*s%wetu 00824 00825 ! V-euler velocities at u-point 00826 veu = s%vu - vsu 00827 ! U-euler velocties at u-point 00828 s%ueu = s%uu - usu 00829 ! Velocity magnitude at u-points 00830 if (par%sedtrans == 0) then 00831 s%vmagu=sqrt(s%uu**2+s%vu**2) 00832 endif 00833 ! Eulerian velocity magnitude at u-points 00834 s%vmageu=sqrt(s%ueu**2+veu**2) 00835 00836 ! Ue and Ve in cell centre; do output and sediment stirring 00837 s%ue(2:s%nx,:)=0.5d0*(s%ueu(1:s%nx-1,:)+s%ueu(2:s%nx,:)) 00838 s%ue(1,:)=s%ueu(1,:) 00839 ! wwvv ue(nx+1,:) ? 00840 00841 if (s%ny>0) then 00842 s%ve(:,2:s%ny)=0.5d0*(s%vev(:,1:s%ny-1)+s%vev(:,2:s%ny)) !Jaap s%ny+1 00843 s%ve(:,1)=s%vev(:,1) 00844 else 00845 s%ve(:,1) = s%vev(:,1) 00846 endif !s%ny>0 00847 ! U-velocities at v-points 00848 if (s%ny>0) then 00849 s%uv(2:s%nx,1:s%ny)= .25d0*(s%uu(1:s%nx-1,1:s%ny)+s%uu(2:s%nx,1:s%ny)+ & 00850 s%uu(1:s%nx-1,2:s%ny+1)+s%uu(2:s%nx,2:s%ny+1)) 00851 ! boundaries? 00852 ! wwvv and what about uv(:,1) ? 00853 if(xmpi_isright) then 00854 s%uv(:,s%ny+1) = s%uv(:,s%ny) 00855 endif 00856 else 00857 s%uv(2:s%nx,1)= .5d0*(s%uu(1:s%nx-1,1)+s%uu(2:s%nx,1)) 00858 endif !s%ny>0 00859 ! wwvv fix uv(:,ny+1) for non-right processes 00860 ! uv(1,:) and uv(nx+1,:) need to be filled in for 00861 ! non-bot or top processes 00862 s%uv=s%uv*s%wetv 00863 ! V-stokes velocities at V point 00864 if (s%ny>0) then 00865 vsv(2:s%nx,1:s%ny)=0.5d0*(s%ust(2:s%nx,1:s%ny)*sinthm(2:s%nx,1:s%ny)+& 00866 s%ust(2:s%nx,2:s%ny+1)*sinthm(2:s%nx,2:s%ny+1)) 00867 if(xmpi_isleft) then 00868 vsv(:,1) = vsv(:,2) 00869 endif 00870 if(xmpi_isright) then 00871 vsv(:,s%ny+1) = vsv(:,s%ny) 00872 endif 00873 else 00874 vsv(2:s%nx,1)= s%ust(2:s%nx,1)*sinthm(2:s%nx,1) 00875 endif !s%ny>0 00876 ! wwvv fix vsv(:,1) and vsv(:,ny+1) and vsv(1,:) and vsv(nx+1,:) 00877 00878 vsv=vsv*s%wetv 00879 ! U-stokes velocities at V point 00880 if (s%ny>0) then 00881 usv(2:s%nx,1:s%ny)=0.5d0*(s%ust(2:s%nx,1:s%ny)*costhm(2:s%nx,1:s%ny)+& 00882 s%ust(2:s%nx,2:s%ny+1)*costhm(2:s%nx,2:s%ny+1)) 00883 if(xmpi_isleft) then 00884 usv(:,1) = usv(:,2) 00885 endif 00886 if(xmpi_isright) then 00887 usv(:,s%ny+1) = usv(:,s%ny) 00888 endif 00889 else 00890 usv(2:s%nx,1)=s%ust(2:s%nx,1)*costhm(2:s%nx,1) 00891 endif !s%ny>0 00892 ! wwvv fix usv(:,1) and usv(:,ny+1) and usv(1,:) and usv(nx+1,:) 00893 usv=usv*s%wetv 00894 00895 ! V-euler velocities at V-point 00896 s%vev = s%vv - vsv 00897 ! U-euler velocties at V-point 00898 uev = s%uv - usv 00899 ! Velocity magnitude at v-points 00900 if (par%sedtrans==0) then 00901 s%vmagv=sqrt(s%uv**2+s%vv**2) 00902 endif 00903 ! Eulerian velocity magnitude at v-points 00904 s%vmagev=sqrt(uev**2+s%vev**2) 00905 00906 ! wwvv vev(nx+1,:) ? 00907 ! 00908 s%hold =s%hh ! wwvv ? s%hold is never else used 00909 ! 00910 s%hh = max(s%zs-s%zb,par%eps) ! water depth 00911 s%zs1 = s%zs-s%zs0 ! water level minus tide 00912 00913 s%maxzs=max(s%zs,s%maxzs) 00914 s%minzs=min(s%zs,s%minzs) 00915 00916 ! XBeach uses same input for Van Rijn, 1993 (quasi 3d), but doesn't use this module 00917 if (par%nz>1 .and. par%form /= FORM_VANRIJN1993) then 00918 do j=1,s%ny+1 00919 do i=1,s%nx+1 00920 if (s%wetz(i,j) .eq. 1) then 00921 ks=12.d0*s%hh(i,j)/10.d0**(sqrt(par%g/s%cfu(i,j))/18.d0) 00922 tauw=par%rhoa*par%Cd*sqrt(s%windsu(i,j)**2+s%windnv(i,j)**2) 00923 tauwx=s%windsu(i,j)*tauw 00924 tauwy=s%windnv(i,j)*tauw 00925 00926 call vsm_u_XB ( s%ue(i,j) ,s%ve(i,j) ,s%hh(i,j) ,ks , & 00927 par%rho ,tauwx ,tauwy ,s%thetamean(i,j)-s%alfaz(i,j), & 00928 s%urms(i,j) ,s%sigm(i,j) ,s%k(i,j) ,s%Dr(i,j) , & 00929 s%E(i,j) ,s%R(i,j) ,fcvisc,facdel ,par%ws , & 00930 s%fw(i,j) ,facdf ,swglm ,par%nz ,s%sigz , & 00931 s%uz(i,j,:) ,s%vz(i,j,:) ,s%ustz(i,j,:) ,s%nutz(i,j,:), & 00932 s%ue_sed(i,j),s%ve_sed(i,j) ) 00933 else 00934 s%uz(i,j,:)=0.d0 00935 s%vz(i,j,:)=0.d0 00936 s%ustz(i,j,:)=0.d0 00937 s%nutz(i,j,:)=0.d0 00938 s%ue_sed(i,j)=0.d0 00939 s%ve_sed(i,j)=0.d0 00940 endif 00941 enddo 00942 enddo 00943 else 00944 s%ue_sed=s%ue 00945 s%ve_sed=s%ve 00946 endif 00947 00948 end subroutine flow 00949 00950 subroutine visc_smagorinsky(s,par) 00951 use params, only: parameters 00952 use spaceparamsdef 00953 use xmpi_module 00954 00955 IMPLICIT NONE 00956 00957 ! DATE AUTHOR CHANGES 00958 ! 00959 ! December 2010 Pieter Bart Smit New Subroutine 00960 ! March 2010 Pieter Bart Smit Changed formulation to standard smag. model 00961 00962 !------------------------------------------------------------------------------- 00963 ! DECLARATIONS 00964 !------------------------------------------------------------------------------- 00965 00966 !-------------------------- PURPOSE ---------------------------- 00967 ! 00968 ! Calculates the turbulent viscocity coefficient nuh according to the smagorinsky 00969 ! subgrid model. 00970 ! 00971 !-------------------------- METHOD ---------------------------- 00972 ! 00973 ! The turbulent viscocity is given as: 00974 ! 00975 ! nuh = C^2*dx*dy*Tau 00976 ! 00977 ! Tau =2^(1/2) * [ (du/dx)^2 + (dv/dy)^2 + 1/2 * (du/dy + dv/dx)^2 ] ^ (1/2) 00978 ! 00979 ! Where 00980 ! 00981 ! dx,dy : grid size 00982 ! C : Constant ~0.15 (set by par%nuh) 00983 ! Tau : Measure for the magnitude of the turbulent stresses 00984 ! 00985 !-------------------------- ARGUMENTS ---------------------------- 00986 00987 type(spacepars),target ,intent(inout) :: s 00988 type(parameters) ,intent(in) :: par 00989 00990 !-------------------------- LOCAL VARIABLES ---------------------------- 00991 00992 real*8 :: dudx !U Velocity gradient in x-dir 00993 real*8 :: dudy !U Velocity gradient in y-dir 00994 real*8 :: dvdx !V Velocity gradient in x-dir 00995 real*8 :: dvdy !V Velocity gradient in y-dir 00996 real*8 :: Tau !Measure for magnitude viscous stresses 00997 real*8 :: l !Local gridcell area 00998 integer :: i !Index variable 00999 integer :: j !Index variable 01000 01001 01002 !------------------------------------------------------------------------------- 01003 ! IMPLEMENTATION 01004 !------------------------------------------------------------------------------- 01005 01006 !MPI WARNING -> Check loop indices 01007 if (s%ny>2) then 01008 do j=2,s%ny 01009 do i=2,s%nx 01010 dudx = (s%uu(i,j)-s%uu(i-1,j))/s%dsz(i,j) 01011 dudy = .5d0*(s%uu(i,j+1) - s%uu(i,j-1) + s%uu(i-1,j+1) - s%uu(i-1,j-1))/(s%dnv(i,j)+s%dnv(i,j-1)) 01012 dvdx = .5d0*(s%vv(i+1,j) - s%vv(i-1,j) + s%vv(i+1,j-1) - s%vv(i-1,j-1))/(s%dsu(i,j)+s%dsu(i-1,j)) 01013 dvdy = (s%vv(i,j)-s%vv(i,j-1))/s%dnz(i,j) 01014 Tau = sqrt(2.0d0 * dudx**2+2.0d0 * dvdy**2 + (dvdx+dudy)**2) 01015 l = 1.d0/s%dsdnzi(i,j) 01016 ! Robert: probably need something like nuh = max(par%nuh,par%nuhfac*s%hh(i,j)*(s%DR(i,j)/par%rho)**(1.0d0/3.0d0)) 01017 ! s%nuh(i,j) = nuh**2 * l * Tau * real(s%wetu(i,j)*s%wetu(i-1,j)*s%wetv(i,j)*s%wetv(i,j-1),kind=8) 01018 s%nuh(i,j) = par%nuh**2 * l * Tau * real(s%wetu(i,j)*s%wetu(i-1,j)*s%wetv(i,j)*s%wetv(i,j-1),kind=8) 01019 enddo 01020 enddo 01021 01022 else 01023 01024 j = max(s%ny,1) 01025 01026 do i=2,s%nx 01027 dudx = (s%uu(i,j)-s%uu(i-1,j))/s%dsz(i,j) 01028 dvdx = (s%vv(i+1,j) - s%vv(i-1,j) )/(s%dsu(i,j)+s%dsu(i-1,j)) 01029 Tau = sqrt(2.0d0 * dudx**2 + dvdx**2) 01030 01031 if (par%dy > -1.d0) then 01032 l = s%dsz(i,j)*par%dy 01033 else 01034 l = s%dsz(i,j)**2 01035 endif 01036 01037 s%nuh(i,j) = par%nuh**2 * l * Tau * real(s%wetu(i,j)*s%wetu(i-1,j),kind=8) 01038 enddo 01039 01040 endif !s%ny>2 01041 01042 if (s%ny>0) then 01043 if (xmpi_isleft) s%nuh(:,1) = s%nuh(:,2) 01044 if (xmpi_isright) s%nuh(:,s%ny+1) = s%nuh(:,s%ny) 01045 endif 01046 01047 if (xmpi_istop) s%nuh(1,:) = s%nuh(2,:) 01048 if (xmpi_isbot) s%nuh(s%nx+1,:) = s%nuh(s%nx,:) 01049 01050 end subroutine visc_smagorinsky 01051 01052 end module flow_timestep_module 01053
1.7.6.1
