|
XBeach
|
00001 module wave_instationary_module 00002 implicit none 00003 save 00004 private 00005 public:: wave_instationary 00006 00007 contains 00008 00009 subroutine wave_instationary(s,par) 00010 00011 use params, only: parameters 00012 use spaceparams 00013 use roelvink_module, only: roelvink, baldock, janssen_battjes 00014 use wave_functions_module,only: compute_wave_direction_velocities, slope2d,breakerdelay, advecqy, advecthetaho, & 00015 & advecxho, advecqx, advecyho, compute_wave_forces, compute_stokes_drift 00016 use xmpi_module 00017 use mnemmodule 00018 use interp, only: Linear_interp 00019 use paramsconst 00020 00021 implicit none 00022 00023 type(spacepars), target :: s 00024 type(parameters) :: par 00025 00026 integer :: i 00027 integer :: j 00028 integer :: itheta 00029 integer :: dummy 00030 00031 00032 00033 real*8 , dimension(:,:) ,allocatable,save :: dhdx,dhdy,dudx,dudy,dvdx,dvdy,ustw,Erfl 00034 real*8 , dimension(:,:) ,allocatable,save :: km,kmx,kmy,xwadvec,ywadvec,sinh2kh !,wm 00035 real*8 , dimension(:,:,:),allocatable,save :: xadvec,yadvec,thetaadvec,dd,drr,dder 00036 real*8 , dimension(:,:,:),allocatable,save :: xradvec,yradvec,thetaradvec 00037 real*8 , dimension(:,:) ,allocatable,save :: dkmxdx,dkmxdy,dkmydx,dkmydy,cgxm,cgym,arg,fac 00038 real*8 , dimension(:,:) ,allocatable,save :: uorb,hhwlocal,RH 00039 real*8 , dimension(:) ,allocatable,save :: wcrestpos 00040 logical, dimension(:,:) ,allocatable,save :: gammax_correct, mask_ee 00041 real*8 :: coffshore 00042 00043 00044 00045 if (.not. allocated(drr)) then 00046 allocate(drr (s%nx+1,s%ny+1,s%ntheta)) 00047 allocate(xadvec (s%nx+1,s%ny+1,s%ntheta)) 00048 allocate(yadvec (s%nx+1,s%ny+1,s%ntheta)) 00049 allocate(thetaadvec (s%nx+1,s%ny+1,s%ntheta)) 00050 allocate(xradvec (s%nx+1,s%ny+1,s%ntheta)) 00051 allocate(yradvec (s%nx+1,s%ny+1,s%ntheta)) 00052 allocate(thetaradvec (s%nx+1,s%ny+1,s%ntheta)) 00053 allocate(dd (s%nx+1,s%ny+1,s%ntheta)) 00054 allocate(dder (s%nx+1,s%ny+1,s%ntheta)) 00055 00056 allocate(dhdx (s%nx+1,s%ny+1)) 00057 allocate(dhdy (s%nx+1,s%ny+1)) 00058 allocate(dudx (s%nx+1,s%ny+1)) 00059 allocate(dudy (s%nx+1,s%ny+1)) 00060 allocate(dvdx (s%nx+1,s%ny+1)) 00061 allocate(dvdy (s%nx+1,s%ny+1)) 00062 allocate(km (s%nx+1,s%ny+1)) 00063 allocate(kmx (s%nx+1,s%ny+1)) 00064 allocate(kmy (s%nx+1,s%ny+1)) 00065 ! allocate(wm (nx+1,ny+1)) 00066 allocate(ustw (s%nx+1,s%ny+1)) 00067 allocate(Erfl (s%nx+1,s%ny+1)) ! wwvv not used 00068 allocate(xwadvec (s%nx+1,s%ny+1)) 00069 allocate(ywadvec (s%nx+1,s%ny+1)) 00070 allocate(sinh2kh (s%nx+1,s%ny+1)) 00071 allocate(dkmxdx (s%nx+1,s%ny+1)) 00072 allocate(dkmxdy (s%nx+1,s%ny+1)) 00073 allocate(dkmydx (s%nx+1,s%ny+1)) 00074 allocate(dkmydy (s%nx+1,s%ny+1)) 00075 allocate(cgxm (s%nx+1,s%ny+1)) 00076 allocate(cgym (s%nx+1,s%ny+1)) 00077 allocate(arg (s%nx+1,s%ny+1)) 00078 allocate(fac (s%nx+1,s%ny+1)) 00079 allocate(uorb (s%nx+1,s%ny+1)) 00080 allocate(hhwlocal (s%nx+1,s%ny+1)) 00081 allocate(RH (s%nx+1,s%ny+1)) 00082 allocate(wcrestpos (s%nx+1)) 00083 allocate(gammax_correct(s%nx+1,s%ny+1)) 00084 allocate(mask_ee (s%nx+1,s%ny+1)) 00085 00086 00087 ! wwvv todo: I think these iniailization are superfluous 00088 drr = 0.d0 00089 xadvec = 0.d0 00090 yadvec = 0.d0 00091 thetaadvec = 0.d0 00092 xradvec = 0.d0 00093 yradvec = 0.d0 00094 thetaradvec = 0.d0 00095 dd = 0.d0 00096 dder = 0.d0 00097 dhdx = 0.d0 00098 dhdy = 0.d0 00099 dudx = 0.d0 00100 dudy = 0.d0 00101 dvdx = 0.d0 00102 dvdy = 0.d0 00103 km = 0.d0 00104 kmx = 0.d0 00105 kmy = 0.d0 00106 ! wm = 0.d0 00107 ustw = 0.d0 00108 Erfl = 0.d0 00109 xwadvec = 0.d0 00110 ywadvec = 0.d0 00111 sinh2kh = 0.d0 00112 dkmxdx = 0.d0 00113 dkmxdy = 0.d0 00114 dkmydx = 0.d0 00115 dkmydy = 0.d0 00116 cgxm = 0.d0 00117 cgym = 0.d0 00118 arg = 0.d0 00119 fac = 0.d0 00120 uorb = 0.d0 00121 hhwlocal = s%hh 00122 RH = 0.d0 00123 gammax_correct = .false. 00124 mask_ee = .false. 00125 s%Fx = 0.d0 ! in spacepars 00126 s%Fy = 0.d0 ! in spacepars 00127 endif 00128 00129 mask_ee (imin_ee:imax_ee,jmin_ee:jmax_ee) = .true. ! Mask tu use in where loops, MPI ok 00130 00131 ! the local water depth to use in this subroutine generally depend on whether we're using wci or not 00132 ! both options include par%delta*s%H effect 00133 ! dispersion routine in wave_timestep already accounts for these differences 00134 if (par%wci==1) then 00135 hhwlocal = s%hhwcins 00136 else 00137 hhwlocal = s%hhw 00138 endif 00139 00140 if (par%single_dir==1) then 00141 s%costh(:,:,1)=cos(s%thetamean-s%alfaz) 00142 s%sinth(:,:,1)=sin(s%thetamean-s%alfaz) 00143 00144 elseif (par%snells==0) then 00145 s%thetamean=(sum(s%ee*s%thet,3)/s%ntheta)/(max(sum(s%ee,3),0.00001d0)/s%ntheta) 00146 elseif(par%snells==1) then !Dano: Snellius 00147 ! Check for borderline cases where critical c/c(1,1) is reached.... 00148 #ifdef USEMPI 00149 if (xmpi_istop .and. xmpi_isleft) then 00150 coffshore = s%c(1,1) 00151 else 00152 coffshore = -huge(0.d0) 00153 endif 00154 call xmpi_allreduce(coffshore,MPI_MAX) 00155 #else 00156 coffshore = s%c(1,1) 00157 #endif 00158 s%thetamean=asin(max(-1.0d0, min(1.0d0, sin(s%theta0-s%alfaz(1,1))*s%c/coffshore)))+s%alfaz(1,1) 00159 s%costh(:,:,1)=cos(s%thetamean-s%alfaz) 00160 s%sinth(:,:,1)=sin(s%thetamean-s%alfaz) 00161 endif 00162 00163 ! Slopes of water depth 00164 call slope2D(hhwlocal,s%nx,s%ny,s%dsu,s%dnv,dhdx,dhdy,s%wete) 00165 if (par%wci==1) then 00166 call slope2D(s%uwcins,s%nx,s%ny,s%dsu,s%dnv,dudx,dudy,s%wete) 00167 call slope2D(s%vwcins,s%nx,s%ny,s%dsu,s%dnv,dvdx,dvdy,s%wete) 00168 else 00169 dudx = 0.d0 00170 dudy = 0.d0 00171 dvdx = 0.d0 00172 dvdy = 0.d0 00173 endif 00174 ! 00175 ! Calculate once sinh(2kh) 00176 where(2*hhwlocal*s%k<=3000.d0) 00177 sinh2kh=sinh(min(2*s%k*hhwlocal,10.0d0)) 00178 elsewhere 00179 sinh2kh = 3000.d0 00180 endwhere 00181 00182 ! split wave velocities in wave grid directions theta 00183 call compute_wave_direction_velocities(s,par,1,dhdx,dhdy,dudx,dudy,dvdx,dvdy,sinh2kh) 00184 ! 00185 ! transform to wave action 00186 ! 00187 do itheta = 1,s%ntheta 00188 where(s%wete == 1) 00189 s%ee(:,:,itheta) = s%ee(:,:,itheta)/s%sigt(:,:,itheta) 00190 endwhere 00191 enddo 00192 ! 00193 ! Upwind Euler timestep propagation 00194 ! 00195 call advecxho(s%ee,s%cgx,xadvec,s%nx,s%ny,s%ntheta,s%dnu,s%dsu,s%dsdnzi,par%scheme,s%wete,par%dt,s%dsz) 00196 if (s%ny>0) then 00197 call advecyho(s%ee,s%cgy,yadvec,s%nx,s%ny,s%ntheta,s%dsv,s%dnv,s%dsdnzi,par%scheme,s%wete,par%dt,s%dnz) 00198 endif 00199 call advecthetaho(s%ee,s%ctheta,thetaadvec,s%nx,s%ny,s%ntheta,s%dtheta,par%scheme,s%wete) 00200 ! 00201 do itheta = 1,s%ntheta 00202 where(s%wete==1) 00203 s%ee(:,:,itheta)=s%ee(:,:,itheta)-par%dt*(xadvec(:,:,itheta)+yadvec(:,:,itheta)+thetaadvec(:,:,itheta)) 00204 endwhere 00205 enddo 00206 ! 00207 ! transform back to wave energy 00208 ! 00209 do itheta = 1,s%ntheta 00210 where(s%wete == 1) 00211 s%ee(:,:,itheta) = max(s%ee(:,:,itheta)*s%sigt(:,:,itheta),0.d0) 00212 elsewhere 00213 s%ee(:,:,itheta) = 0.d0 00214 endwhere 00215 enddo 00216 ! 00217 ! Energy integrated over wave directions,Hrms 00218 ! 00219 where(s%wete == 1) 00220 s%E=sum(s%ee,3)*s%dtheta 00221 elsewhere 00222 s%E = 0.d0 00223 endwhere 00224 s%H=sqrt(s%E/par%rhog8) 00225 ! 00226 ! Correct for gammax in these areas 00227 where(s%H>par%gammax*s%hh .and. s%wete==1 .and. mask_ee) 00228 gammax_correct = .true. 00229 elsewhere 00230 gammax_correct = .false. 00231 endwhere 00232 ! 00233 do itheta=1,s%ntheta 00234 ! note: here we do use instantaneous water depth (and excluding par%delta*s%H effect) 00235 where(gammax_correct) 00236 s%ee(:,:,itheta)=s%ee(:,:,itheta)/(s%H/(par%gammax*s%hh))**2 00237 endwhere 00238 enddo 00239 where(gammax_correct) 00240 s%H=min(s%H,par%gammax*s%hh) 00241 endwhere 00242 where(gammax_correct) 00243 s%E=par%rhog8*s%H**2 ! thread-safe computation of all H before recomputing E 00244 endwhere 00245 ! 00246 ! Total dissipation 00247 ! 00248 select case(par%break) 00249 case(BREAK_ROELVINK1,BREAK_ROELVINK2) 00250 call roelvink(par,s) 00251 case(BREAK_BALDOCK) 00252 call baldock(par,s) 00253 case(BREAK_JANSSEN) 00254 call janssen_battjes(par,s) 00255 case(BREAK_ROELVINK_DALY) ! Dano: probably not MPI ok 00256 cgxm = s%c*dcos(s%thetamean) 00257 cgym = s%c*dsin(s%thetamean) 00258 call advecqx(cgxm,s%Qb,xwadvec,s%nx,s%ny,s%dsu,s%wete) 00259 if (s%ny>0) then 00260 call advecqy(cgym,s%Qb,ywadvec,s%nx,s%ny,s%dnv,s%wete) 00261 s%Qb = s%Qb-par%dt*(xwadvec+ywadvec) 00262 else 00263 s%Qb = s%Qb-par%dt*xwadvec 00264 endif 00265 call roelvink(par,s) 00266 end select 00267 00268 ! Dissipation by bed friction 00269 uorb=par%px*s%H/par%Trep/sinh(min(max(s%k,0.01d0)*s%hhw,10.0d0)) 00270 s%Df= 2d0/(3d0*par%px)*par%rho*s%fw*uorb**3 00271 where (s%hh>par%fwcutoff) 00272 s%Df = 0.d0 00273 end where 00274 ! 00275 ! Distribution of dissipation over directions and frequencies 00276 ! 00277 do itheta=1,s%ntheta 00278 ! Only calculate for E>0 FB 00279 ! First just the dissipation that is fed to the roller 00280 dder(:,:,itheta)=s%ee(:,:,itheta)*s%D/max(s%E,0.00001d0) 00281 ! Then all short wave energy dissipation, including bed friction and vegetation 00282 dd(:,:,itheta)=dder(:,:,itheta) + s%ee(:,:,itheta)*(s%Df+s%Dveg)/max(s%E,0.00001d0) 00283 enddo 00284 ! 00285 ! Euler step dissipation 00286 ! 00287 ! calculate roller energy balance 00288 ! 00289 call advecxho(s%rr,s%cx,xradvec,s%nx,s%ny,s%ntheta,s%dnu,s%dsu,s%dsdnzi,par%scheme,s%wete,par%dt,s%dsz) 00290 if (s%ny>0) then 00291 call advecyho(s%rr,s%cy,yradvec,s%nx,s%ny,s%ntheta,s%dsv,s%dnv,s%dsdnzi,par%scheme,s%wete,par%dt,s%dnz) 00292 endif 00293 !call advectheta(rr*ctheta,thetaradvec,nx,ny,ntheta,dtheta) 00294 call advecthetaho(s%rr,s%ctheta,thetaradvec,s%nx,s%ny,s%ntheta,s%dtheta,par%scheme,s%wete) 00295 do itheta=1,s%ntheta 00296 where(s%wete==1 .and. mask_ee) 00297 s%rr(:,:,itheta)=s%rr(:,:,itheta) & 00298 -par%dt*(xradvec(:,:,itheta) & 00299 +yradvec(:,:,itheta) & 00300 +thetaradvec(:,:,itheta)) 00301 s%rr(:,:,itheta)=max(s%rr(:,:,itheta),0.0d0) 00302 elsewhere (s%wete==1) 00303 s%rr(:,:,itheta) = 0.d0 00304 endwhere 00305 enddo 00306 ! 00307 ! euler step roller energy dissipation (source and sink function) 00308 ! 00309 do itheta=1,s%ntheta 00310 do j=jmin_ee,jmax_ee 00311 do i=imin_ee,imax_ee 00312 if(s%wete(i,j)==1) then 00313 s%ee(i,j,itheta)=s%ee(i,j,itheta)-par%dt*dd(i,j,itheta) 00314 if(par%roller==1) then 00315 drr(i,j,itheta) = 2*par%g*s%BR(i,j)*max(s%rr(i,j,itheta),0.0d0)/ & 00316 sqrt(s%cx(i,j,itheta)**2 +s%cy(i,j,itheta)**2) 00317 s%rr(i,j,itheta)=s%rr(i,j,itheta)+par%dt*(dder(i,j,itheta) & ! Robert: changed from dd to dder 00318 -drr(i,j,itheta)) ! (only from wave breaking, 00319 else ! not vegetation or bed friction) 00320 s%rr(i,j,itheta)= 0.0d0 00321 drr(i,j,itheta)= 0.0d0 00322 endif 00323 s%ee(i,j,itheta)=max(s%ee(i,j,itheta),0.0d0) 00324 s%rr(i,j,itheta)=max(s%rr(i,j,itheta),0.0d0) 00325 else 00326 s%ee(i,j,itheta)=0.0d0 00327 s%rr(i,j,itheta)=0.0d0 00328 drr(i,j,itheta)=0.0d0 00329 endif 00330 enddo 00331 enddo 00332 enddo 00333 ! turn off roller energy in non-wet points 00334 do itheta=1,s%ntheta 00335 where(s%wetz==0) 00336 s%rr(:,:,itheta) = 0.d0 00337 drr(:,:,itheta) = 0.d0 00338 endwhere 00339 enddo 00340 ! 00341 ! Bay boundary Robert + Jaap 00342 ! 00343 ! wwvv 00344 ! this has consequences for the parallel version, 00345 ! but also, if we do nothing, there are discrepancies 00346 ! between ee and rr in the different processes. We need to 00347 ! get valid values for ee(nx+1,:,:) and rr(nx+1,:,:) from 00348 ! the neighbour below. We cannot postpone this until this 00349 ! subroutine ends, because ee and rr are used in this subroutine 00350 00351 if (xmpi_isbot) then 00352 s%ee(s%nx+1,:,:) =s%ee(s%nx,:,:) 00353 s%rr(s%nx+1,:,:) =s%rr(s%nx,:,:) 00354 endif 00355 00356 if (par%t>0.0d0) then 00357 if (s%ny>0) then 00358 if (xmpi_isleft)then ! Jaap 00359 if (par%lateralwave==LATERALWAVE_NEUMANN) then 00360 ! 00361 ! Lateral boundary at y=0; 00362 ! 00363 s%ee(1:s%nx+1,1,:)=s%ee(1:s%nx+1,2,:) 00364 s%rr(1:s%nx+1,1,:)=s%rr(1:s%nx+1,2,:) 00365 elseif (par%lateralwave==LATERALWAVE_WAVECREST) then 00366 ! wcrestpos=xz+tan(thetamean(:,2))*(yz(2)-yz(1)) 00367 wcrestpos=s%sdist(:,1)+tan(s%thetamean(:,2)-s%alfaz(:,2))*s%dnv(:,1) 00368 do itheta=1,s%ntheta 00369 do i=1,s%nx+1 00370 call Linear_interp((/-huge(0.d0) ,wcrestpos ,huge(0.d0)/),& 00371 (/s%ee(1,2,itheta),s%ee(:,2,itheta),s%ee(s%nx+1,2,itheta)/),& 00372 s%nx+1,s%sdist(i,1),s%ee(i,1,itheta),dummy) 00373 call Linear_interp((/-huge(0.d0) ,wcrestpos ,huge(0.d0)/),& 00374 (/s%rr(1,2,itheta),s%rr(:,2,itheta),s%rr(s%nx+1,2,itheta)/),& 00375 s%nx+1,s%sdist(i,1),s%rr(i,1,itheta),dummy) 00376 s%ee(i,1,itheta)=max(s%ee(i,1,itheta),0.d0) 00377 s%rr(i,1,itheta)=max(s%rr(i,1,itheta),0.d0) 00378 enddo 00379 enddo 00380 endif 00381 endif 00382 if (xmpi_isright)then 00383 if (par%lateralwave==LATERALWAVE_NEUMANN) then 00384 ! 00385 ! lateral; boundary at y=ny*dy 00386 ! 00387 s%ee(1:s%nx+1,s%ny+1,:)=s%ee(1:s%nx+1,s%ny,:) 00388 s%rr(1:s%nx+1,s%ny+1,:)=s%rr(1:s%nx+1,s%ny,:) 00389 elseif (par%lateralwave==LATERALWAVE_WAVECREST) then 00390 ! wcrestpos=xz-tan(thetamean(:,ny))*(yz(ny+1)-yz(ny)) 00391 wcrestpos=s%sdist(:,s%ny+1)-tan(s%thetamean(:,s%ny)-s%alfaz(:,s%ny))*s%dnv(:,s%ny) 00392 do itheta=1,s%ntheta 00393 do i=1,s%nx+1 00394 call Linear_interp((/-huge(0.d0) ,wcrestpos ,huge(0.d0)/),& 00395 (/s%ee(1,s%ny,itheta),s%ee(:,s%ny,itheta),s%ee(s%nx+1,s%ny,itheta)/),& 00396 s%nx+1,s%sdist(i,s%ny+1),s%ee(i,s%ny+1,itheta),dummy) 00397 call Linear_interp((/-huge(0.d0) ,wcrestpos ,huge(0.d0)/),& 00398 (/s%rr(1,s%ny,itheta),s%rr(:,s%ny,itheta),s%rr(s%nx+1,s%ny,itheta)/),& 00399 s%nx+1,s%sdist(i,s%ny+1),s%rr(i,s%ny+1,itheta),dummy) 00400 s%ee(i,s%ny+1,itheta)=max(s%ee(i,s%ny+1,itheta),0.d0) 00401 s%rr(i,s%ny+1,itheta)=max(s%rr(i,s%ny+1,itheta),0.d0) 00402 enddo 00403 enddo 00404 endif 00405 endif 00406 endif 00407 endif 00408 ! wwvv communicate ee(:,1,:) 00409 #ifdef USEMPI 00410 call xmpi_shift_ee(s%ee) 00411 call xmpi_shift_ee(s%rr) 00412 #endif 00413 00414 ! 00415 ! Energy integrated over wave directions,Hrms 00416 ! 00417 s%E = sum(s%ee,3)*s%dtheta 00418 s%R = sum(s%rr,3)*s%dtheta 00419 s%DR = sum(drr,3)*s%dtheta 00420 s%H = sqrt(s%E/par%rhog8) 00421 !Dano: recompute uorb as it is affected by the MPI shifts and time stepping 00422 uorb=par%px*s%H/par%Trep/sinh(min(max(s%k,0.01d0)*s%hhw,10.0d0)) 00423 ! 00424 ! Correct roller for gammax in these areas 00425 if (par%rollergammax==1) then 00426 RH = sqrt(s%R/par%rhog8) 00427 where(RH>par%gammax*s%hhw .and. s%wete==1 .and. mask_ee) 00428 gammax_correct = .true. 00429 elsewhere 00430 gammax_correct = .false. 00431 endwhere 00432 ! 00433 do itheta=1,s%ntheta 00434 ! note: here we do use instantaneous water depth (and excluding par%delta*s%H effect) 00435 where(gammax_correct) 00436 s%rr(:,:,itheta)=s%rr(:,:,itheta)/(RH/(par%gammax*s%hhw))**2 00437 endwhere 00438 enddo 00439 where(gammax_correct) 00440 s%R=min(s%R,par%gammax*s%hhw) 00441 endwhere 00442 endif 00443 ! 00444 ! Compute mean wave direction 00445 ! 00446 if (par%snells==0 .and. par%single_dir==0) then 00447 where(s%wete == 1) 00448 s%thetamean=(sum(s%ee*s%thet,3)/s%ntheta)/(max(sum(s%ee,3),0.00001d0)/s%ntheta) 00449 endwhere 00450 endif 00451 ! 00452 ! compute wave forces 00453 call compute_wave_forces(s) 00454 ! 00455 ! orbital velocity 00456 s%urms=uorb/sqrt(2.d0) 00457 ! 00458 ! Stokes drift and orbital velocities 00459 call compute_stokes_drift(s,par) 00460 00461 end subroutine wave_instationary 00462 00463 end module wave_instationary_module
1.7.6.1
