|
XBeach
|
00001 module wave_directions_module 00002 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00003 ! Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University ! 00004 ! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski, ! 00005 ! Jaap van Thiel de Vries, Robert McCall ! 00006 ! ! 00007 ! d.roelvink@unesco-ihe.org ! 00008 ! UNESCO-IHE Institute for Water Education ! 00009 ! P.O. Box 3015 ! 00010 ! 2601 DA Delft ! 00011 ! The Netherlands ! 00012 ! ! 00013 ! This library is free software; you can redistribute it and/or ! 00014 ! modify it under the terms of the GNU Lesser General Public ! 00015 ! License as published by the Free Software Foundation; either ! 00016 ! version 2.1 of the License, or (at your option) any later version. ! 00017 ! ! 00018 ! This library is distributed in the hope that it will be useful, ! 00019 ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! 00020 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! 00021 ! Lesser General Public License for more details. ! 00022 ! ! 00023 ! You should have received a copy of the GNU Lesser General Public ! 00024 ! License along with this library; if not, write to the Free Software ! 00025 ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 ! 00026 ! USA ! 00027 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00028 implicit none 00029 save 00030 00031 private 00032 public wave_directions 00033 contains 00034 00035 subroutine wave_directions(s,par) 00036 use params, only: parameters 00037 use spaceparams 00038 use roelvink_module, only: roelvink, baldock, janssen_battjes 00039 use wave_functions_module, only: slope2d, advecthetaho, advecxho, advecyho, compute_wave_direction_velocities 00040 use xmpi_module 00041 use logging_module, only: writelog, progress_indicator 00042 use paramsconst 00043 00044 ! wwvv in my testcase, this routine was not called, so it is not 00045 ! tested. Nevertheless, I put in code for the parallel version. 00046 00047 IMPLICIT NONE 00048 00049 type(spacepars), target :: s 00050 type(parameters) :: par 00051 00052 integer :: i,imax,i1 00053 integer :: j 00054 integer :: itheta,iter,itermax 00055 integer,save :: scheme_local,scheme_local_yadvec 00056 real*8,dimension(:),allocatable,save :: dist,factor,e01 00057 real*8 , dimension(:,:) ,allocatable,save :: dhdx,dhdy,dudx,dudy,dvdx,dvdy 00058 real*8 , dimension(:,:) ,allocatable,save :: uorb 00059 real*8 , dimension(:,:) ,allocatable,save :: sinh2kh ! ,wm 00060 real*8 , dimension(:,:,:),allocatable,save :: xadvec,yadvec,thetaadvec,dd 00061 real*8 , dimension(:),allocatable,save :: Hprev,thetaprev 00062 real*8 :: Herr,thetaerr,thetaerr_cyc,dtw 00063 logical :: stopiterate 00064 00065 if (.not. allocated(e01)) then 00066 allocate(e01(1:s%ntheta_s)) 00067 allocate(dist(1:s%ntheta_s)) 00068 allocate(factor(1:s%ntheta_s)) 00069 allocate(xadvec (s%nx+1,s%ny+1,s%ntheta_s)) 00070 allocate(yadvec (s%nx+1,s%ny+1,s%ntheta_s)) 00071 allocate(thetaadvec(s%nx+1,s%ny+1,s%ntheta_s)) 00072 allocate(dd (s%nx+1,s%ny+1,s%ntheta_s)) 00073 00074 00075 allocate(dhdx(s%nx+1,s%ny+1)) 00076 allocate(dhdy(s%nx+1,s%ny+1)) 00077 allocate(dudx(s%nx+1,s%ny+1)) 00078 allocate(dudy(s%nx+1,s%ny+1)) 00079 allocate(dvdx(s%nx+1,s%ny+1)) 00080 allocate(dvdy(s%nx+1,s%ny+1)) 00081 allocate(uorb(s%nx+1,s%ny+1)) 00082 allocate(sinh2kh(s%nx+1,s%ny+1)) 00083 allocate(Hprev(s%ny+1)) 00084 allocate(thetaprev(s%ny+1)) 00085 00086 00087 if(par%scheme==SCHEME_WARMBEAM) then 00088 scheme_local = SCHEME_UPWIND_2 ! note, this is already stated as warning in log file during read of params.txt 00089 else 00090 scheme_local = par%scheme 00091 endif 00092 if(par%cyclic==1 .and. scheme_local .ne. SCHEME_UPWIND_1) then 00093 scheme_local_yadvec = SCHEME_UPWIND_1 00094 call writelog('lws','(a)', 'Warning: y-advection in stationary wave solver set to upwind_1 for convergence') 00095 else 00096 scheme_local_yadvec = scheme_local 00097 endif 00098 endif 00099 00100 00101 xadvec = 0.0d0 00102 yadvec = 0.0d0 00103 thetaadvec = 0.0d0 00104 dd = 0.0d0 00105 dhdx = 0.0d0 00106 dhdy = 0.0d0 00107 dudx = 0.0d0 00108 dudy = 0.0d0 00109 dvdx = 0.0d0 00110 dvdy = 0.0d0 00111 sinh2kh = 0.0d0 00112 00113 00114 ! Slopes of water depth 00115 call slope2D(s%hhws,s%nx,s%ny,s%dsu,s%dnv,dhdx,dhdy,s%wete) 00116 ! Dano limit slopes used in refraction to avoid unrealistic refraction speeds 00117 dhdx=sign(1.d0,dhdx)*min(abs(dhdx),0.1d0) 00118 dhdy=sign(1.d0,dhdy)*min(abs(dhdy),0.1d0) 00119 if (par%wci==1) then 00120 call slope2D(s%uws,s%nx,s%ny,s%dsu,s%dnv,dudx,dudy,s%wete) 00121 call slope2D(s%vws,s%nx,s%ny,s%dsu,s%dnv,dvdx,dvdy,s%wete) 00122 else 00123 dudx = 0.d0 00124 dudy = 0.d0 00125 dvdx = 0.d0 00126 dvdy = 0.d0 00127 endif 00128 ! wwvv these slope routines are in wave_timestep, and are 00129 ! MPI-aware 00130 ! 00131 ! Calculate once sinh(2kh) 00132 00133 where(s%wete==1 .and. 2*s%hhws*s%k<=3000.d0) 00134 sinh2kh=sinh(min(2*s%k*s%hhws,10.0d0)) 00135 elsewhere 00136 sinh2kh = 3000.d0 00137 endwhere 00138 00139 ! all dry cells have zero energy 00140 do itheta=1,s%ntheta_s 00141 where(s%wete==0) 00142 s%ee_s(:,:,itheta) = 0.d0 00143 endwhere 00144 enddo 00145 where(s%wete==0) 00146 s%E=0.d0 00147 s%H=0.d0 00148 endwhere 00149 00150 forall (i=1:s%nx+1,j=1:s%ny+1,s%wete(i,j)==1) 00151 s%thetamean(i,j)=(sum(s%ee_s(i,j,:)*s%thet_s(i,j,:))/s%ntheta_s)/(max(sum(s%ee_s(i,j,:)),0.00001d0)/s%ntheta_s) 00152 endforall 00153 00154 ! 00155 ! split wave velocities in wave grid directions theta 00156 call compute_wave_direction_velocities(s,par,2,dhdx,dhdy,dudx,dudy,dvdx,dvdy,sinh2kh) 00157 00158 !dtw=.9*minval(xz(2:nx+1)-xz(1:nx))/maxval(cgx_s) 00159 00160 ! Initialize at the boundary 00161 s%E(1,:)=max(sum(s%ee_s(1,:,:),2)*s%dtheta_s,0.d0) 00162 s%H(1,:)=sqrt(s%E(1,:)/par%rhog8) 00163 00164 imax=s%nx 00165 !Dano This is ok, since we will set mpiboundary to y in stationary mode 00166 ! 00167 ! write to screen that waves are updated 00168 if(xmaster) call writelog('ls','(a,f0.2,a)','Computing wave directions at t = ',par%t,' s') 00169 call progress_indicator(.true.,0.d0,5.d0,2.d0) 00170 itermax = 0 00171 do i=2,imax 00172 call progress_indicator(.false.,dble(i)/imax*100,5.d0,2.d0) 00173 00174 if(par%wci==0) then 00175 dtw= .5d0*minval(s%dsu(i-1:i+1,jmin_ee:jmax_ee))/max(1.0d-6,maxval(abs(s%cgx_s(i-1:i+1,jmin_ee:jmax_ee,:)))) 00176 dtw=min(dtw,.5d0*minval(s%dnv(i,jmin_ee:jmax_ee)) /max(1.0d-6,maxval(abs(s%cgy_s(i,jmin_ee:jmax_ee,:) ))) ) 00177 dtw=min(dtw,.5d0*s%dtheta_s /max(1.0d-6,maxval(abs(s%ctheta_s(i,jmin_ee:jmax_ee,:) ))) ) 00178 !Dano: need to make sure all processes use the same dtw, min of all processes 00179 #ifdef USEMPI 00180 call xmpi_allreduce(dtw,MPI_MIN) 00181 #endif 00182 else 00183 dtw = par%dt 00184 endif 00185 00186 Herr=huge(1.d0) 00187 thetaerr = 2*par%px 00188 iter=0 00189 stopiterate=.false. 00190 do while (stopiterate .eqv. .false.) 00191 iter=iter+1 00192 Hprev=s%H(i,:) 00193 thetaprev = s%thetamean(i,:) 00194 ! 00195 ! transform to wave action 00196 ! 00197 i1=max(i-2,1) 00198 do itheta=1,s%ntheta_s 00199 where(s%wete(i1:i+1,:)==1) 00200 s%ee_s(i1:i+1,:,itheta) = s%ee_s(i1:i+1,:,itheta)/s%sigm(i1:i+1,:) 00201 endwhere 00202 enddo 00203 ! 00204 ! Upwind Euler timestep propagation 00205 ! 00206 ! Robert: can we also use WARMBEAM here? [Robert: answer is no!] 00207 if (i>2.and. scheme_local==SCHEME_UPWIND_2) then 00208 call advecxho(s%ee_s(i-2:i+1,:,:),s%cgx_s(i-2:i+1,:,:),xadvec(i-2:i+1,:,:), & 00209 3,s%ny,s%ntheta_s,s%dnu(i-2:i+1,:),s%dsu(i-2:i+1,:),s%dsdnzi(i-2:i+1,:),SCHEME_UPWIND_2, & 00210 s%wete(i-2:i+1,:),dtw,s%dsz) 00211 else 00212 call advecxho(s%ee_s(i-1:i+1,:,:),s%cgx_s(i-1:i+1,:,:),xadvec(i-1:i+1,:,:), & 00213 2,s%ny,s%ntheta_s,s%dnu(i-1:i+1,:),s%dsu(i-1:i+1,:),s%dsdnzi(i-1:i+1,:),SCHEME_UPWIND_1,& 00214 s%wete(i-1:i+1,:),dtw,s%dsz) 00215 endif 00216 ! Robert: could this also be warm-beam? 00217 if (s%ny>0) then 00218 call advecyho(s%ee_s(i,:,:),s%cgy_s(i,:,:),yadvec(i,:,:), & 00219 0,s%ny,s%ntheta_s,s%dsv(i,:),s%dnv(i,:),s%dsdnzi(i,:),scheme_local_yadvec,s%wete(i,:),dtw,s%dnz) 00220 endif 00221 ! Robert: could this also be warm-beam? 00222 call advecthetaho(s%ee_s(i,:,:),s%ctheta_s(i,:,:),thetaadvec(i,:,:),0,s%ny,s%ntheta_s, & 00223 s%dtheta_s,scheme_local,s%wete(i,:)) 00224 00225 s%ee_s(i,:,:)=s%ee_s(i,:,:)-dtw*(xadvec(i,:,:) + yadvec(i,:,:) + thetaadvec(i,:,:)) 00226 ! 00227 ! transform back to wave energy 00228 ! 00229 do itheta=1,s%ntheta_s 00230 where(s%wete(i1:i+1,:)==1) 00231 s%ee_s(i1:i+1,:,itheta) = s%ee_s(i1:i+1,:,itheta)*s%sigm(i1:i+1,:) 00232 s%ee_s(i1:i+1,:,itheta)=max(s%ee_s(i1:i+1,:,itheta),0.0d0) 00233 endwhere 00234 enddo 00235 #ifdef USEMPI 00236 call xmpi_shift(s%ee_s(i-1:i+1,:,:),SHIFT_Y_R,1,2) 00237 call xmpi_shift(s%ee_s(i-1:i+1,:,:),SHIFT_Y_L,3,4) 00238 #endif 00239 ! 00240 ! Energy integrated over wave directions,Hrms 00241 ! 00242 forall(j=1:s%ny+1,s%wete(i,j)==1) 00243 s%E(i,j)=sum(s%ee_s(i,j,:))*s%dtheta_s 00244 endforall 00245 where(s%wete(i,:)==1) 00246 s%H(i,:)=sqrt(s%E(i,:)/par%rhog8) 00247 endwhere 00248 ! check with gammax 00249 do j=1,s%ny+1 00250 if (s%wete(i,j)==1 .and. s%H(i,j)>par%gammax*s%hhws(i,j)) then 00251 do itheta=1,s%ntheta_s 00252 s%ee_s(i,j,itheta)=s%ee_s(i,j,itheta)/(s%H(i,j)/(par%gammax*s%hhws(i,j)))**2 00253 enddo 00254 s%H(i,j)=par%gammax*s%hhws(i,j) 00255 s%E(i,j)=par%rhog8*s%H(i,j)**2 00256 endif 00257 enddo 00258 ! 00259 if (par%snells==0) then !Dano not for SNellius 00260 where(s%wete(i,:)==1) 00261 s%thetamean(i,:) = (sum(s%ee_s(i,:,:)*s%thet_s(i,:,:),2)/s%ntheta_s)/ & 00262 (max(sum(s%ee_s(i,:,:),2),0.000010d0)/s%ntheta_s) 00263 endwhere 00264 endif 00265 ! 00266 ! Total dissipation 00267 00268 select case(par%break) 00269 case(BREAK_ROELVINK1,BREAK_ROELVINK2) 00270 call roelvink (par,s,i) 00271 case(BREAK_BALDOCK) 00272 call baldock (par,s,i) 00273 case(BREAK_JANSSEN) 00274 call janssen_battjes(par,s,i) 00275 end select 00276 00277 00278 ! Dissipation by bed friction 00279 where(s%fw(i,:)>0.d0 .and. s%wete(i,:)==1 .and. s%hhws(i,:)<=par%fwcutoff) 00280 uorb(i,:)=par%px*s%H(i,:)/par%Trep/sinh(min(max(s%k(i,:),0.01d0)*s%hhws(i,:),10.0d0)) 00281 s%Df(i,:)=0.6666666d0/par%px*par%rho*s%fw(i,:)*uorb(i,:)**3 00282 elsewhere 00283 s%Df(i,:) = 0.d0 00284 end where 00285 ! 00286 ! Distribution of dissipation over directions and frequencies 00287 ! 00288 do itheta=1,s%ntheta_s 00289 where(s%wete(i,:)==1) 00290 dd(i,:,itheta)=s%ee_s(i,:,itheta)*(s%D(i,:)+s%Df(i,:))/max(s%E(i,:),0.00001d0) ! Robert: missing Dveg here! 00291 elsewhere 00292 dd(i,:,itheta)=0.d0 00293 endwhere 00294 end do 00295 ! 00296 ! Euler step dissipation 00297 ! 00298 do j=jmin_ee,jmax_ee 00299 do itheta=1,s%ntheta_s 00300 if (dtw*dd(i,j,itheta)>s%ee_s(i,j,itheta) .and. s%wete(i,j)==1) then 00301 dtw=min(dtw,.5*s%ee_s(i,j,itheta)/dd(i,j,itheta)) 00302 endif 00303 enddo 00304 enddo 00305 !Dano: need to make sure all processes use the same dtw, min of all processes 00306 #ifdef USEMPI 00307 call xmpi_allreduce(dtw,MPI_MIN) 00308 #endif 00309 do j=1,s%ny+1 00310 do itheta=1,s%ntheta_s 00311 if(s%wete(i,j)==1) then 00312 s%ee_s(i,j,itheta)=s%ee_s(i,j,itheta)-dtw*dd(i,j,itheta) 00313 s%ee_s(i,j,itheta)=max(s%ee_s(i,j,itheta),0.0d0) 00314 else 00315 s%ee_s(i,j,itheta)=0.0d0 00316 end if 00317 end do 00318 end do 00319 ! Lateral boundary condition 00320 if (xmpi_isleft .and. s%ny>0) then 00321 do itheta=1,s%ntheta_s 00322 !if (s%sinth_s(i,1,itheta)>=0.) then ! Robert: this is different to wave_stationary 00323 s%ee_s(i,1,itheta)=s%ee_s(i,2,itheta) 00324 !endif 00325 enddo 00326 s%k(:,1)=s%k(:,2) 00327 s%sigm(:,1)=s%sigm(:,2) 00328 endif 00329 if (xmpi_isright .and. s%ny>0) then 00330 do itheta=1,s%ntheta_s 00331 !if (s%sinth_s(i,s%ny+1,itheta)<=0.) then 00332 s%ee_s(i,s%ny+1,itheta)=s%ee_s(i,s%ny,itheta) 00333 !endif 00334 end do 00335 s%k(:,s%ny+1)=s%k(:,s%ny) 00336 s%sigm(:,s%ny+1)=s%sigm(:,s%ny) 00337 endif 00338 ! 00339 ! Compute mean wave direction 00340 ! 00341 if (par%snells==0) then 00342 where(s%wete(i,:)==1) 00343 s%thetamean(i,:)=(sum(s%ee_s(i,:,:)*s%thet_s(i,:,:),2)/size(s%ee_s(i,:,:),2)) & 00344 /(max(sum(s%ee_s(i,:,:),2),0.000010d0) /size(s%ee_s(i,:,:),2)) 00345 elsewhere 00346 ! forward-copy wave directions on dry cells in case they flood before the next wave update 00347 s%thetamean(i,:) = s%thetamean(i-1,:) 00348 endwhere 00349 endif 00350 ! 00351 ! Energy integrated over wave directions,Hrms 00352 ! 00353 where(s%wete(i,:)==1) 00354 s%E(i,:)=sum(s%ee_s(i,:,:),2)*s%dtheta_s 00355 s%H(i,:)=sqrt(s%E(i,:)/par%rhog8) 00356 elsewhere 00357 s%E(i,:)=0.d0 00358 s%H(i,:)=0.d0 00359 endwhere 00360 Herr=maxval(abs(Hprev(jmin_ee:jmax_ee)-s%H(i,jmin_ee:jmax_ee)),mask=s%wete(i,:)==1) 00361 00362 thetaerr = 0.d0 00363 do j=jmin_ee,jmax_ee 00364 thetaerr_cyc = minval( (/ abs(thetaprev(j)-(s%thetamean(i,j)-2*par%px)) , & 00365 abs(thetaprev(j)-(s%thetamean(i,j) )) , & 00366 abs(thetaprev(j)-(s%thetamean(i,j)+2*par%px)) /) ) 00367 thetaerr = max(thetaerr,thetaerr_cyc) 00368 enddo 00369 00370 #ifdef USEMPI 00371 call xmpi_allreduce(Herr,MPI_MAX) 00372 call xmpi_allreduce(thetaerr,MPI_MAX) 00373 #endif 00374 00375 ! Stopping criteria 00376 if (iter<par%maxiter) then 00377 if (Herr<par%maxerror .and. thetaerr<par%maxerror_angle) then 00378 stopiterate=.true. 00379 !if(xmaster) call writelog('ls','(a,i4,a,i4)','Wave propagation row ',i,', iteration ',iter) 00380 endif 00381 else 00382 stopiterate=.true. 00383 if(xmaster) call writelog('ls','(a,i4,a,i4,a,f5.4)','Wave propagation row ',i,', iteration ',iter,', error: ',Herr) 00384 endif 00385 enddo ! End while loop 00386 itermax = max(itermax,iter) 00387 enddo ! End do i=2:s%nx loop 00388 call writelog('ls','(a,i4)','Maximum number of iterations: ',itermax) 00389 00390 s%ee_s(s%nx+1,:,:) = s%ee_s(s%nx,:,:) 00391 s%E(s%nx+1,:) = s%E(s%nx,:) 00392 s%H(s%nx+1,:) = s%H(s%nx,:) 00393 s%k(s%nx+1,:) = s%k(s%nx,:) 00394 s%sigm(s%nx+1,:) = s%sigm(s%nx,:) 00395 s%cg(s%nx+1,:) = s%cg(s%nx,:) 00396 s%c(s%nx+1,:) = s%c(s%nx,:) 00397 s%thet_s(s%nx+1,:,:) = s%thet_s(s%nx,:,:) 00398 00399 end subroutine wave_directions 00400 end module wave_directions_module
1.7.6.1
