|
XBeach
|
00001 module wave_stationary_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 contains 00032 00033 subroutine wave_stationary_directions(s,par,callType) 00034 use params 00035 use spaceparams 00036 use roelvink_module 00037 use wave_functions_module 00038 use xmpi_module 00039 use logging_module 00040 use paramsconst 00041 00042 IMPLICIT NONE 00043 00044 type(spacepars), target :: s 00045 type(parameters) :: par 00046 integer :: callType 00047 00048 integer,parameter :: callTypeStationary = 0 00049 integer,parameter :: callTypeDirections = 1 00050 00051 integer :: i,imax,i1 00052 integer :: j 00053 integer :: itheta,iter,itermax 00054 integer :: ntheta_local 00055 integer,save :: scheme_local,scheme_local_yadvec 00056 real*8 , dimension(:,:) ,allocatable,save :: dhdx,dhdy,dudx,dudy,dvdx,dvdy 00057 real*8 , dimension(:,:),allocatable,save :: uorb 00058 real*8 , dimension(:,:) ,allocatable,save :: sinh2kh 00059 real*8 , dimension(:,:) ,allocatable,save :: hhwlocal 00060 real*8 , dimension(:,:,:),allocatable,save :: xadvec,yadvec,thetaadvec,dd,drr,dder 00061 real*8 , dimension(:,:,:),allocatable,save :: xradvec,yradvec,thetaradvec 00062 real*8 , dimension(:,:,:),allocatable,save :: ee_local 00063 real*8 , dimension(:),allocatable,save :: Hprev,thetaprev,gammafac_i,RH 00064 real*8 :: Herr,thetaerr,thetaerr_cyc,dtw 00065 real*8 :: dtheta_local 00066 logical :: stopiterate 00067 ! 00068 ! set dimension sizes 00069 select case (callType) 00070 case (callTypeStationary) 00071 ntheta_local = s%ntheta 00072 case (callTypeDirections) 00073 ntheta_local = s%ntheta_s 00074 end select 00075 ! 00076 ! Allocate 00077 if (.not. allocated(xadvec)) then 00078 allocate(xadvec (s%nx+1,s%ny+1,ntheta_local)) 00079 allocate(yadvec (s%nx+1,s%ny+1,ntheta_local)) 00080 allocate(thetaadvec(s%nx+1,s%ny+1,ntheta_local)) 00081 allocate(dd (s%nx+1,s%ny+1,ntheta_local)) 00082 allocate(dder (s%nx+1,s%ny+1,ntheta_local)) 00083 allocate(ee_local (s%nx+1,s%ny+1,ntheta_local)) 00084 allocate(dhdx(s%nx+1,s%ny+1)) 00085 allocate(dhdy(s%nx+1,s%ny+1)) 00086 allocate(dudx(s%nx+1,s%ny+1)) 00087 allocate(dudy(s%nx+1,s%ny+1)) 00088 allocate(dvdx(s%nx+1,s%ny+1)) 00089 allocate(dvdy(s%nx+1,s%ny+1)) 00090 allocate(uorb(s%nx+1,s%ny+1)) 00091 allocate(sinh2kh(s%nx+1,s%ny+1)) 00092 allocate(hhwlocal(s%nx+1,s%ny+1)) 00093 allocate(gammafac_i(s%ny+1)) 00094 allocate(RH(s%ny+1)) 00095 allocate(Hprev(s%ny+1)) 00096 allocate(thetaprev(s%ny+1)) 00097 ! 00098 if(par%scheme==SCHEME_WARMBEAM) then 00099 scheme_local = SCHEME_UPWIND_2 ! note, this is already stated as warning in log file during read of params.txt 00100 else 00101 scheme_local = par%scheme 00102 endif 00103 if(scheme_local .ne. SCHEME_UPWIND_1) then 00104 scheme_local_yadvec = SCHEME_UPWIND_1 00105 call writelog('lws','(a)', 'Warning: y-advection in stationary wave solver set to upwind_1 for convergence') 00106 else 00107 scheme_local_yadvec = scheme_local 00108 endif 00109 endif 00110 ! 00111 ! Allocate additional variables for stationary wave computation 00112 if ((.not. allocated(xradvec)) .and. (callType==callTypeStationary)) then 00113 allocate(xradvec(s%nx+1,s%ny+1,ntheta_local)) 00114 allocate(yradvec(s%nx+1,s%ny+1,ntheta_local)) 00115 allocate(thetaradvec(s%nx+1,s%ny+1,ntheta_local)) 00116 allocate(drr (s%nx+1,s%ny+1,ntheta_local)) 00117 endif 00118 ! 00119 ! Initialize 00120 xadvec = 0.0d0 00121 yadvec = 0.0d0 00122 thetaadvec = 0.0d0 00123 dd = 0.0d0 00124 dder = 0.0d0 00125 dhdx = 0.0d0 00126 dhdy = 0.0d0 00127 dudx = 0.0d0 00128 dudy = 0.0d0 00129 dvdx = 0.0d0 00130 dvdy = 0.0d0 00131 sinh2kh = 0.0d0 00132 uorb = 0.0d0 00133 ! 00134 if (callType==callTypeStationary) then 00135 xradvec = 0.0d0 00136 yradvec = 0.0d0 00137 thetaradvec = 0.0d0 00138 drr = 0.0d0 00139 endif 00140 ! 00141 ! Set switch for pointer variables (not strictly defined as Fortran pointers) 00142 select case (callType) 00143 case (callTypeStationary) 00144 hhwlocal = s%hhw 00145 ee_local = s%ee 00146 dtheta_local = s%dtheta 00147 case (callTypeDirections) 00148 hhwlocal = s%hhws 00149 ee_local = s%ee_s 00150 dtheta_local = s%dtheta_s 00151 end select 00152 ! 00153 ! Slopes of water depth 00154 call slope2D(hhwlocal,s%nx,s%ny,s%dsu,s%dnv,dhdx,dhdy,s%wete) 00155 ! 00156 ! Dano: limit slopes used in refraction to avoid unrealistic refraction speeds 00157 dhdx=sign(1.d0,dhdx)*min(abs(dhdx),0.1d0) 00158 dhdy=sign(1.d0,dhdy)*min(abs(dhdy),0.1d0) 00159 ! 00160 ! slopes of flow (in case of WCI) 00161 if (par%wci==1) then 00162 select case (callType) 00163 case(callTypeStationary) 00164 call slope2D(s%u,s%nx,s%ny,s%dsu,s%dnv,dudx,dudy,s%wete) 00165 call slope2D(s%v,s%nx,s%ny,s%dsu,s%dnv,dvdx,dvdy,s%wete) 00166 case (callTypeDirections) 00167 call slope2D(s%uws,s%nx,s%ny,s%dsu,s%dnv,dudx,dudy,s%wete) 00168 call slope2D(s%vws,s%nx,s%ny,s%dsu,s%dnv,dvdx,dvdy,s%wete) 00169 end select 00170 else 00171 dudx = 0.d0 00172 dudy = 0.d0 00173 dvdx = 0.d0 00174 dvdy = 0.d0 00175 endif 00176 ! wwvv these slope routines are in wave_timestep, and are 00177 ! MPI-aware 00178 ! 00179 ! Calculate once sinh(2kh) 00180 where(s%wete==1 .and. 2*hhwlocal*s%k<=3000.d0) 00181 sinh2kh=sinh(min(2*s%k*hhwlocal,10.0d0)) 00182 elsewhere 00183 sinh2kh = 3000.d0 00184 endwhere 00185 ! 00186 ! all dry cells have zero energy 00187 do itheta=1,ntheta_local 00188 where(s%wete==0) 00189 ee_local(:,:,itheta) = 0.d0 00190 endwhere 00191 enddo 00192 where(s%wete==0) 00193 s%E=0.d0 00194 s%H=0.d0 00195 endwhere 00196 ! 00197 ! wave directions 00198 select case (callType) 00199 case(callTypeStationary) 00200 if (par%snells==0) then 00201 forall (i=1:s%nx+1,j=1:s%ny+1,s%wete(i,j)==1) 00202 s%thetamean(i,j)=(sum(ee_local(i,j,:)*s%thet(i,j,:))/ntheta_local) / & 00203 (max(sum(ee_local(i,j,:)),0.00001d0)/ntheta_local) 00204 endforall 00205 else 00206 s%thetamean=asin(sin(s%theta0-s%alfaz(1,1))*s%c/s%c(1,1))+s%alfaz(1,1) 00207 s%costh(:,:,1)=cos(s%thetamean-s%alfaz) 00208 s%sinth(:,:,1)=sin(s%thetamean-s%alfaz) 00209 endif 00210 case(callTypeDirections) 00211 forall (i=1:s%nx+1,j=1:s%ny+1,s%wete(i,j)==1) 00212 s%thetamean(i,j)=(sum(ee_local(i,j,:)*s%thet_s(i,j,:))/ntheta_local) / & 00213 (max(sum(ee_local(i,j,:)),0.00001d0)/ntheta_local) 00214 endforall 00215 end select 00216 ! 00217 ! split wave velocities in wave grid directions theta 00218 select case (callType) 00219 case(callTypeStationary) 00220 call compute_wave_direction_velocities(s,par,0,dhdx,dhdy,dudx,dudy,dvdx,dvdy,sinh2kh) 00221 case(callTypeDirections) 00222 call compute_wave_direction_velocities(s,par,2,dhdx,dhdy,dudx,dudy,dvdx,dvdy,sinh2kh) 00223 end select 00224 ! 00225 ! Initialize at the boundary 00226 s%E(1,:)=max(sum(ee_local(1,:,:),2)*dtheta_local,0.d0) 00227 s%H(1,:)=sqrt(s%E(1,:)/par%rhog8) 00228 if (callType==callTypeStationary) s%R(1,:)=max(sum(s%rr(1,:,:),2)*dtheta_local,0.0d0) 00229 ! 00230 ! 00231 ! Initialize loop over grid 00232 ! 00233 ! write to screen that waves are updated 00234 if(xmaster) then 00235 select case (callType) 00236 case(callTypeStationary) 00237 call writelog('ls','(a,f0.2,a)','Computing wave transformation at t = ',par%t,' s') 00238 case(callTypeDirections) 00239 call writelog('ls','(a,f0.2,a)','Computing wave directions at t = ',par%t,' s') 00240 end select 00241 endif 00242 call progress_indicator(.true.,0.d0,5.d0,2.d0) 00243 imax=s%nx 00244 itermax = 0 00245 ! 00246 ! 00247 ! Start of loop of the grid 00248 do i=2,imax 00249 call progress_indicator(.false.,dble(i)/imax*100,5.d0,2.d0) 00250 ! 00251 ! determine internal time step 00252 select case (callType) 00253 case(callTypeStationary) 00254 dtw= .5d0*minval(s%dsu(i-1:i+1,jmin_ee:jmax_ee))/max(1.0d-6,maxval(abs(s%cgx(i-1:i+1,jmin_ee:jmax_ee,:)))) 00255 dtw=min(dtw,.5d0*minval(s%dnv(i,jmin_ee:jmax_ee)) /max(1.0d-6,maxval(abs(s%cgy(i,jmin_ee:jmax_ee,:) )))) 00256 dtw=min(dtw,.5d0*dtheta_local /max(1.0d-6,maxval(abs(s%ctheta(i,jmin_ee:jmax_ee,:) )))) 00257 case(callTypeDirections) 00258 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,:)))) 00259 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,:) )))) 00260 dtw=min(dtw,.5d0*dtheta_local /max(1.0d-6,maxval(abs(s%ctheta_s(i,jmin_ee:jmax_ee,:) )))) 00261 end select 00262 ! 00263 !Dano: need to make sure all processes use the same dtw, min of all processes 00264 #ifdef USEMPI 00265 call xmpi_allreduce(dtw,MPI_MIN) 00266 #endif 00267 ! 00268 ! Reduced internal time step in case of WCI 00269 if(par%wci==1) then 00270 dtw = min(dtw,par%dt) 00271 endif 00272 ! 00273 ! Initialize iteration loop per gridline 00274 Herr=huge(1.d0) 00275 thetaerr = 2*par%px 00276 iter=0 00277 stopiterate=.false. 00278 ! 00279 ! Start of iteration loop per gridline 00280 do while (stopiterate .eqv. .false.) 00281 iter=iter+1 00282 Hprev=s%H(i,:) 00283 thetaprev = s%thetamean(i,:) 00284 ! 00285 ! transform energy to wave action 00286 i1=max(i-2,1) 00287 do itheta=1,ntheta_local 00288 where(s%wete(i1:i+1,:)==1) 00289 ee_local(i1:i+1,:,itheta) = ee_local(i1:i+1,:,itheta)/s%sigm(i1:i+1,:) 00290 endwhere 00291 enddo 00292 ! 00293 ! Upwind Euler timestep propagation 00294 ! 00295 ! M-Direction 00296 select case (callType) 00297 case(callTypeStationary) 00298 if (i==2) then 00299 call advecxho(ee_local(i-1:i+1,:,:),s%cgx(i-1:i+1,:,:),xadvec(i-1:i+1,:,:), & 00300 2,s%ny,ntheta_local,s%dnu(i-1:i+1,:),s%dsu(i-1:i+1,:),s%dsdnzi(i-1:i+1,:),SCHEME_UPWIND_1, & 00301 s%wete(i-1:i+1,:),dtw,s%dsz) 00302 else 00303 ! note: at the moment Warm-Beam is currently not a valid combination (automatic change to upwind_2 in params.F90) 00304 call advecxho(ee_local(i-2:i+1,:,:),s%cgx(i-2:i+1,:,:),xadvec(i-2:i+1,:,:), & 00305 3,s%ny,ntheta_local,s%dnu(i-2:i+1,:),s%dsu(i-2:i+1,:),s%dsdnzi(i-2:i+1,:),scheme_local, & 00306 s%wete(i-2:i+1,:),dtw,s%dsz) 00307 endif 00308 case(callTypeDirections) 00309 if (i==2) then 00310 call advecxho(ee_local(i-1:i+1,:,:),s%cgx_s(i-1:i+1,:,:),xadvec(i-1:i+1,:,:), & 00311 2,s%ny,ntheta_local,s%dnu(i-1:i+1,:),s%dsu(i-1:i+1,:),s%dsdnzi(i-1:i+1,:),SCHEME_UPWIND_1, & 00312 s%wete(i-1:i+1,:),dtw,s%dsz) 00313 else 00314 ! note: at the moment Warm-Beam is currently not a valid combination (automatic change to upwind_2 in params.F90) 00315 call advecxho(ee_local(i-2:i+1,:,:),s%cgx_s(i-2:i+1,:,:),xadvec(i-2:i+1,:,:), & 00316 3,s%ny,ntheta_local,s%dnu(i-2:i+1,:),s%dsu(i-2:i+1,:),s%dsdnzi(i-2:i+1,:),scheme_local, & 00317 s%wete(i-2:i+1,:),dtw,s%dsz) 00318 endif 00319 end select 00320 ! 00321 ! N-Direction 00322 if (s%ny>0) then 00323 select case (callType) 00324 case(callTypeStationary) 00325 call advecyho(ee_local(i,:,:),s%cgy(i,:,:),yadvec(i,:,:), & 00326 0,s%ny,ntheta_local,s%dsv(i,:),s%dnv(i,:),s%dsdnzi(i,:),scheme_local_yadvec,s%wete(i,:),dtw,s%dnz) 00327 case(callTypeDirections) 00328 call advecyho(ee_local(i,:,:),s%cgy_s(i,:,:),yadvec(i,:,:), & 00329 0,s%ny,ntheta_local,s%dsv(i,:),s%dnv(i,:),s%dsdnzi(i,:),scheme_local_yadvec,s%wete(i,:),dtw,s%dnz) 00330 end select 00331 endif 00332 ! 00333 ! Theta-Direction 00334 select case (callType) 00335 case(callTypeStationary) 00336 call advecthetaho(ee_local(i,:,:),s%ctheta(i,:,:),thetaadvec(i,:,:),0,s%ny,ntheta_local, & 00337 dtheta_local,scheme_local,s%wete(i,:)) 00338 case(callTypeDirections) 00339 call advecthetaho(ee_local(i,:,:),s%ctheta_s(i,:,:),thetaadvec(i,:,:),0,s%ny,ntheta_local, & 00340 dtheta_local,scheme_local,s%wete(i,:)) 00341 end select 00342 ! 00343 ! update wave action 00344 ee_local(i,:,:)=ee_local(i,:,:)-dtw*(xadvec(i,:,:) + yadvec(i,:,:) + thetaadvec(i,:,:)) 00345 ! 00346 ! transform back to wave energy 00347 do itheta=1,ntheta_local 00348 where(s%wete(i1:i+1,:)==1) 00349 ee_local(i1:i+1,:,itheta) = ee_local(i1:i+1,:,itheta)*s%sigm(i1:i+1,:) 00350 ee_local(i1:i+1,:,itheta)=max(ee_local(i1:i+1,:,itheta),0.0d0) 00351 endwhere 00352 enddo 00353 #ifdef USEMPI 00354 call xmpi_shift(ee_local(i-1:i+1,:,:),SHIFT_Y_R,1,2) 00355 call xmpi_shift(ee_local(i-1:i+1,:,:),SHIFT_Y_L,3,4) 00356 #endif 00357 ! 00358 ! Energy integrated over wave directions,Hrms 00359 where(s%wete(i,:)==1) 00360 s%E(i,:)=sum(ee_local(i,:,:),2)*dtheta_local 00361 s%H(i,:)=sqrt(s%E(i,:)/par%rhog8) 00362 endwhere 00363 ! 00364 ! adjust wave energy and height for gammax 00365 gammafac_i = (par%gammax*hhwlocal(i,:))/s%H(i,:) ! inverse of gammafac = s%H(i,:)/(par%gammax*hhwlocal(i,:)) 00366 do j=1,s%ny+1 00367 if (s%wete(i,j)==1 .and. gammafac_i(j)<1.d0) then 00368 do itheta=1,ntheta_local 00369 ee_local(i,j,itheta)=ee_local(i,j,itheta)*gammafac_i(j)**2 00370 enddo 00371 s%H(i,j)=s%H(i,j)*gammafac_i(j) 00372 s%E(i,j)=s%E(i,j)*gammafac_i(j)**2 00373 endif 00374 enddo 00375 ! 00376 ! redo wave directions (Dano: not for Snellius; Robert: also always for callTypeDirections) 00377 select case (callType) 00378 case(callTypeStationary) 00379 if (par%snells==0) then 00380 where(s%wete(i,:)==1) 00381 s%thetamean(i,:) = (sum(ee_local(i,:,:)*s%thet(i,:,:),2)/ntheta_local)/ & 00382 (max(sum(ee_local(i,:,:),2),0.000010d0)/ntheta_local) 00383 endwhere 00384 endif 00385 case(callTypeDirections) 00386 where(s%wete(i,:)==1) 00387 s%thetamean(i,:) = (sum(ee_local(i,:,:)*s%thet_s(i,:,:),2)/ntheta_local)/ & 00388 (max(sum(ee_local(i,:,:),2),0.000010d0)/ntheta_local) 00389 endwhere 00390 end select 00391 ! 00392 ! 00393 ! Compute wave dissipation 00394 ! 00395 ! Dissipation by breaking 00396 select case(par%break) 00397 case(BREAK_ROELVINK1,BREAK_ROELVINK2) 00398 call roelvink (par,s,i) 00399 case(BREAK_BALDOCK) 00400 call baldock (par,s,i) 00401 case(BREAK_JANSSEN) 00402 call janssen_battjes(par,s,i) 00403 end select 00404 ! 00405 ! Dissipation by bed friction 00406 where(s%fw(i,:)>0.d0 .and. s%wete(i,:)==1 .and. hhwlocal(i,:)<=par%fwcutoff) 00407 uorb(i,:)=par%px*s%H(i,:)/par%Trep/sinh(min(max(s%k(i,:),0.01d0)*hhwlocal(i,:),10.0d0)) 00408 s%Df(i,:)=0.28d0*par%rho*s%fw(i,:)*uorb(i,:)**3 ! Robert: note in wave-instationary the coefficient is 2/(3pi)?? 00409 elsewhere 00410 s%Df(i,:) = 0.d0 00411 end where 00412 ! 00413 ! Dissipation by vegetation computed in vegetation module 00414 ! 00415 ! 00416 ! Distribution of dissipation over directions and frequencies 00417 do itheta=1,ntheta_local 00418 where(s%wete(i,:)==1) 00419 ! just dissipation due to breaking, used for roller energy balance (in callTypeStationary only) 00420 dder(i,:,itheta)=ee_local(i,:,itheta)*s%D(i,:)/max(s%E(i,:),0.00001d0) 00421 ! Then all short wave energy dissipation, including bed friction and vegetation, used to reduce wave height 00422 ! Robert: note Dveg needs update for callTypeDirections, because based on only instantaneous H 00423 dd(i,:,itheta)=dder(i,:,itheta) + ee_local(i,:,itheta)*(s%Df(i,:)+s%Dveg(i,:))/max(s%E(i,:),0.00001d0) 00424 elsewhere 00425 dder(i,:,itheta)=0.d0 00426 dd(i,:,itheta)=0.d0 00427 endwhere 00428 end do 00429 ! 00430 ! 00431 ! Calculate roller energy balance (only needed for callTypeStationary) 00432 if (callType==callTypeStationary .and. par%roller==1) then 00433 ! 00434 ! M-direction 00435 if (i==2) then 00436 call advecxho(s%rr(i-1:i+1,:,:),s%cx(i-1:i+1,:,:),xradvec(i-1:i+1,:,:), & 00437 2,s%ny,ntheta_local,s%dnu(i-1:i+1,:),s%dsu(i-1:i+1,:),s%dsdnzi(i-1:i+1,:),SCHEME_UPWIND_1, & 00438 s%wete(i-1:i+1,:),dtw,s%dsz) 00439 else 00440 ! note: at the moment Warm-Beam is currently not a valid combination (automatic change to upwind_2 in params.F90) 00441 call advecxho(s%rr(i-2:i+1,:,:),s%cx(i-2:i+1,:,:),xradvec(i-2:i+1,:,:), & 00442 3,s%ny,ntheta_local,s%dnu(i-2:i+1,:),s%dsu(i-2:i+1,:),s%dsdnzi(i-2:i+1,:),scheme_local, & 00443 s%wete(i-2:i+1,:),dtw,s%dsz) 00444 endif 00445 ! 00446 ! N-Direction 00447 if (s%ny>0) then 00448 call advecyho(s%rr(i,:,:),s%cy(i,:,:),yradvec(i,:,:), & 00449 0,s%ny,ntheta_local,s%dsv(i,:),s%dnv(i,:),s%dsdnzi(i,:),scheme_local_yadvec,s%wete(i,:),dtw,s%dnz) 00450 endif 00451 ! 00452 ! Theta-Direction 00453 call advecthetaho(s%rr(i,:,:),s%ctheta(i,:,:),thetaradvec(i,:,:), & 00454 0,s%ny,ntheta_local,dtheta_local,scheme_local,s%wete(i,:)) 00455 00456 s%rr(i,:,:)=s%rr(i,:,:)-dtw*(xradvec(i,:,:)+yradvec(i,:,:)+thetaradvec(i,:,:)) 00457 s%rr(i,:,:)=max(s%rr(i,:,:),0.0d0) 00458 #ifdef USEMPI 00459 call xmpi_shift(s%rr(i-1:i,:,:),SHIFT_Y_R,1,2) 00460 call xmpi_shift(s%rr(i-1:i,:,:),SHIFT_Y_L,3,4) 00461 #endif 00462 endif 00463 ! 00464 ! 00465 ! Compute energy source and sink function 00466 ! 00467 ! update internal timestep 00468 do j=jmin_ee,jmax_ee 00469 if (s%wete(i,j)==1) then 00470 do itheta=1,ntheta_local 00471 if(ee_local(i,j,itheta)>0.d0 .and. dd(i,j,itheta)>0.d0) then 00472 dtw=min(dtw,.5*ee_local(i,j,itheta)/dd(i,j,itheta)) 00473 endif 00474 enddo 00475 endif 00476 enddo 00477 ! Dano: need to make sure all processes use the same dtw, min of all processes 00478 #ifdef USEMPI 00479 call xmpi_allreduce(dtw,MPI_MIN) 00480 #endif 00481 ! 00482 ! Timestep update of source and sink 00483 do j=jmin_ee,jmax_ee 00484 if(s%wete(i,j)==1) then 00485 do itheta=1,ntheta_local 00486 ee_local(i,j,itheta)=ee_local(i,j,itheta)-dtw*dd(i,j,itheta) 00487 ee_local(i,j,itheta)=max(ee_local(i,j,itheta),0.0d0) 00488 if(callType==callTypeStationary) then 00489 if (par%roller==1) then !Christophe 00490 drr(i,j,itheta) = 2*par%g*par%beta*max(s%rr(i,j,itheta),0.0d0)/& 00491 sqrt(s%cx(i,j,itheta)**2 +s%cy(i,j,itheta)**2) 00492 s%rr(i,j,itheta) = s%rr(i,j,itheta)+dtw*dder(i,j,itheta)-dtw*drr(i,j,itheta) 00493 s%rr(i,j,itheta) = max(s%rr(i,j,itheta),0.0d0) 00494 else 00495 drr(i,j,itheta)= 0.0d0 00496 s%rr(i,j,itheta)= 0.0d0 00497 end if 00498 endif 00499 enddo 00500 else 00501 ee_local(i,j,:)=0.0d0 00502 if(callType==callTypeStationary) then 00503 s%rr(i,j,:)=0.0d0 00504 drr(i,j,:)=0.0d0 00505 endif 00506 endif 00507 enddo 00508 ! 00509 ! Lateral boundary conditions 00510 if (xmpi_isleft .and. s%ny>0) then 00511 do itheta=1,ntheta_local 00512 ee_local(i,1,itheta)=ee_local(i,2,itheta) 00513 if(callType==callTypeStationary) s%rr(i,1,itheta)=s%rr(i,2,itheta) 00514 enddo 00515 s%k(:,1)=s%k(:,2) 00516 s%sigm(:,1)=s%sigm(:,2) 00517 endif 00518 if (xmpi_isright .and. s%ny>0) then 00519 do itheta=1,ntheta_local 00520 ee_local(i,s%ny+1,itheta)=ee_local(i,s%ny,itheta) 00521 if(callType==callTypeStationary) s%rr(i,s%ny+1,itheta)=s%rr(i,s%ny,itheta) 00522 end do 00523 s%k(:,s%ny+1)=s%k(:,s%ny) 00524 s%sigm(:,s%ny+1)=s%sigm(:,s%ny) 00525 endif 00526 ! 00527 ! 00528 ! Redo wave directions (Dano: not for Snellius; Robert: also always for callTypeDirections) 00529 select case (callType) 00530 case(callTypeStationary) 00531 if (par%snells==0) then 00532 where(s%wete(i,:)==1) 00533 s%thetamean(i,:) = (sum(ee_local(i,:,:)*s%thet(i,:,:),2)/ntheta_local)/ & 00534 (max(sum(ee_local(i,:,:),2),0.000010d0)/ntheta_local) 00535 endwhere 00536 endif 00537 case(callTypeDirections) 00538 where(s%wete(i,:)==1) 00539 s%thetamean(i,:) = (sum(ee_local(i,:,:)*s%thet_s(i,:,:),2)/ntheta_local)/ & 00540 (max(sum(ee_local(i,:,:),2),0.000010d0)/ntheta_local) 00541 endwhere 00542 end select 00543 ! 00544 ! forward-copy wave directions on dry cells in case they flood before the next wave update 00545 where(s%wete(i,:)==0) 00546 s%thetamean(i,:) = s%thetamean(i-1,:) 00547 endwhere 00548 ! 00549 ! 00550 ! Energy integrated over wave directions,Hrms 00551 where(s%wete(i,:)==1) 00552 s%E(i,:)=sum(ee_local(i,:,:),2)*dtheta_local 00553 s%H(i,:)=sqrt(s%E(i,:)/par%rhog8) 00554 elsewhere 00555 s%E(i,:)=0.d0 00556 s%H(i,:)=0.d0 00557 endwhere 00558 ! 00559 ! Roller energy 00560 if (callType==callTypeStationary) then 00561 where(s%wete(i,:)==1) 00562 s%R(i,:)=sum(s%rr(i,:,:),2)*dtheta_local 00563 s%DR(i,:)=sum(drr(i,:,:),2)*dtheta_local 00564 elsewhere 00565 s%R(i,:)=0.d0 00566 s%DR(i,:)=0.d0 00567 endwhere 00568 endif 00569 ! 00570 ! Correct roller for gammax 00571 ! adjust wave energy and height for gammax 00572 if(callType==callTypeStationary .and. par%rollergammax==1) then 00573 RH = sqrt(s%R(i,:)/par%rhog8) 00574 gammafac_i = (par%gammax*hhwlocal(i,:))/RH ! inverse of gammafac = s%H(i,:)/(par%gammax*hhwlocal(i,:)) 00575 do j=1,s%ny+1 00576 if (s%wete(i,j)==1 .and. gammafac_i(j)<1.d0) then 00577 do itheta=1,ntheta_local 00578 s%rr(i,j,itheta)=s%rr(i,j,itheta)*gammafac_i(j)**2 00579 enddo 00580 s%R(i,j)=RH(j)*gammafac_i(j) 00581 endif 00582 enddo 00583 endif 00584 ! 00585 ! 00586 ! Wave height and direction difference with last iteration 00587 Herr=maxval(abs(Hprev(jmin_ee:jmax_ee)-s%H(i,jmin_ee:jmax_ee)),mask=s%wete(i,jmin_ee:jmax_ee)==1) 00588 thetaerr = 0.d0 00589 do j=jmin_ee,jmax_ee 00590 thetaerr_cyc = minval( (/ abs(thetaprev(j)-(s%thetamean(i,j)-2*par%px)) , & 00591 abs(thetaprev(j)-(s%thetamean(i,j) )) , & 00592 abs(thetaprev(j)-(s%thetamean(i,j)+2*par%px)) /) ) 00593 thetaerr = max(thetaerr,thetaerr_cyc) 00594 enddo 00595 ! 00596 ! communicate error across MPI domains 00597 #ifdef USEMPI 00598 call xmpi_allreduce(Herr,MPI_MAX) 00599 call xmpi_allreduce(thetaerr,MPI_MAX) 00600 #endif 00601 ! 00602 ! 00603 ! Stopping criteria 00604 if (iter<par%maxiter) then 00605 if (Herr<par%maxerror .and. thetaerr<par%maxerror_angle) then 00606 stopiterate=.true. 00607 endif 00608 else 00609 stopiterate=.true. 00610 if(xmaster) then 00611 if (Herr>=par%maxerror) then 00612 call writelog('lsw','(a,i4,a,i4,a,f5.4)','Wave propagation row ',i,', iteration ',iter,', H error: ',Herr) 00613 endif 00614 if (thetaerr>=par%maxerror_angle) then 00615 call writelog('lsw','(a,i4,a,i4,a,f5.4)','Wave propagation row ',i,', iteration ', & 00616 iter,', dir error: ',thetaerr/par%px*180) 00617 endif 00618 endif 00619 endif 00620 ! 00621 enddo ! end interation loop per gridline 00622 ! 00623 ! Store maximum iterations update 00624 itermax = max(itermax,iter) 00625 enddo ! end loop over grid 2,nx 00626 ! 00627 ! 00628 ! write summary to log 00629 call writelog('ls','(a,i4)','Maximum number of iterations: ',itermax) 00630 ! 00631 ! 00632 ! store local energy "pointer" variable to s structure 00633 select case (callType) 00634 case(callTypeStationary) 00635 s%ee= ee_local 00636 case(callTypeDirections) 00637 s%ee_s = ee_local 00638 end select 00639 ! 00640 ! 00641 ! Boundary conditions at nx+1 00642 s%E(s%nx+1,:) = s%E(s%nx,:) 00643 s%H(s%nx+1,:) = s%H(s%nx,:) 00644 s%k(s%nx+1,:) = s%k(s%nx,:) 00645 s%sigm(s%nx+1,:) = s%sigm(s%nx,:) 00646 s%cg(s%nx+1,:) = s%cg(s%nx,:) 00647 s%c(s%nx+1,:) = s%c(s%nx,:) 00648 select case (callType) 00649 case(callTypeStationary) 00650 s%ee(s%nx+1,:,:) = s%ee(s%nx,:,:) 00651 s%rr(s%nx+1,:,:) = s%rr(s%nx,:,:) 00652 s%R(s%nx+1,:) = s%R(s%nx,:) 00653 s%DR(s%nx+1,:) = s%DR(s%nx,:) 00654 s%thet(s%nx+1,:,:) = s%thet(s%nx,:,:) 00655 case(callTypeDirections) 00656 s%ee_s(s%nx+1,:,:) = s%ee_s(s%nx,:,:) 00657 s%thet_s(s%nx+1,:,:) = s%thet_s(s%nx,:,:) 00658 end select 00659 ! 00660 ! 00661 ! Compute radiation stress and forcing terms (callType=callTypeStationary only) 00662 if (callType==callTypeStationary) then 00663 ! 00664 ! wave forces 00665 call compute_wave_forces(s) 00666 ! 00667 ! orbital velocity 00668 s%urms=uorb/sqrt(2.d0) 00669 ! 00670 ! Stokes drift and orbital velocities 00671 call compute_stokes_drift(s,par) 00672 endif 00673 ! 00674 end subroutine wave_stationary_directions 00675 end module wave_stationary_directions_module
1.7.6.1
