|
XBeach
|
00001 module wave_stationary_module 00002 implicit none 00003 save 00004 private 00005 public:: wave_stationary 00006 contains 00007 subroutine wave_stationary(s,par) 00008 use params, only: parameters 00009 use spaceparams 00010 use roelvink_module, only: roelvink, baldock, janssen_battjes 00011 use wave_functions_module, only: dispersion, slope2d, breakerdelay, & 00012 & advecthetaho, advecxho, advecyho, compute_wave_direction_velocities 00013 use xmpi_module 00014 use logging_module, only: writelog, progress_indicator 00015 use paramsconst 00016 00017 ! wwvv in my testcase, this routine was not called, so it is not 00018 ! tested. Nevertheless, I put in code for the parallel version. 00019 00020 IMPLICIT NONE 00021 00022 type(spacepars), target :: s 00023 type(parameters) :: par 00024 00025 integer :: i,imax,i1 00026 integer :: j 00027 integer :: itheta,iter,itermax 00028 integer,save :: scheme_local_yadvec 00029 real*8 , dimension(:,:) ,allocatable,save :: dhdx,dhdy,dudx,dudy,dvdx,dvdy,ustw 00030 real*8 , dimension(:,:) ,allocatable,save :: sinh2kh ! ,wm 00031 real*8 , dimension(:,:,:),allocatable,save :: xadvec,yadvec,thetaadvec,dd,drr,dder 00032 real*8 , dimension(:,:,:),allocatable,save :: xradvec,yradvec,thetaradvec 00033 real*8 , dimension(:),allocatable,save :: Hprev,thetaprev 00034 real*8 :: Herr,thetaerr,thetaerr_cyc,dtw 00035 real*8 , dimension(:,:),allocatable,save :: uorb 00036 logical :: stopiterate 00037 00038 !include 's.ind' 00039 !include 's.inp' 00040 00041 if (.not. allocated(xadvec)) then 00042 allocate(xadvec (s%nx+1,s%ny+1,s%ntheta)) 00043 allocate(yadvec (s%nx+1,s%ny+1,s%ntheta)) 00044 allocate(thetaadvec(s%nx+1,s%ny+1,s%ntheta)) 00045 allocate(xradvec (s%nx+1,s%ny+1,s%ntheta)) 00046 allocate(yradvec (s%nx+1,s%ny+1,s%ntheta)) 00047 allocate(thetaradvec(s%nx+1,s%ny+1,s%ntheta)) 00048 allocate(dd (s%nx+1,s%ny+1,s%ntheta)) 00049 allocate(drr (s%nx+1,s%ny+1,s%ntheta)) 00050 allocate(dder (s%nx+1,s%ny+1,s%ntheta)) 00051 00052 allocate(dhdx(s%nx+1,s%ny+1)) 00053 allocate(dhdy(s%nx+1,s%ny+1)) 00054 allocate(dudx(s%nx+1,s%ny+1)) 00055 allocate(dudy(s%nx+1,s%ny+1)) 00056 allocate(dvdx(s%nx+1,s%ny+1)) 00057 allocate(dvdy(s%nx+1,s%ny+1)) 00058 allocate(ustw(s%nx+1,s%ny+1)) 00059 allocate(sinh2kh(s%nx+1,s%ny+1)) 00060 allocate(Hprev(s%ny+1)) 00061 allocate(thetaprev(s%ny+1)) 00062 00063 allocate(uorb (s%nx+1,s%ny+1)) 00064 00065 00066 if(par%cyclic==1 .and. par%scheme .ne. SCHEME_UPWIND_1) then 00067 scheme_local_yadvec = SCHEME_UPWIND_1 00068 call writelog('lws','(a)', 'Warning: y-advection in stationary wave solver set to upwind_1 for convergence') 00069 else 00070 scheme_local_yadvec = par%scheme 00071 endif 00072 endif 00073 00074 xadvec = 0.0d0 00075 yadvec = 0.0d0 00076 thetaadvec = 0.0d0 00077 xradvec = 0.0d0 00078 yradvec = 0.0d0 00079 thetaradvec = 0.0d0 00080 dd = 0.0d0 00081 drr = 0.0d0 00082 dder = 0.0d0 00083 dhdx = 0.0d0 00084 dhdy = 0.0d0 00085 dudx = 0.0d0 00086 dudy = 0.0d0 00087 dvdx = 0.0d0 00088 dvdy = 0.0d0 00089 ustw = 0.0d0 00090 sinh2kh = 0.0d0 00091 00092 uorb = 0.0d0 00093 00094 ! Slopes of water depth 00095 call slope2D(s%hhw,s%nx,s%ny,s%dsu,s%dnv,dhdx,dhdy,s%wete) 00096 ! Dano limit slopes used in refraction to avoid unrealistic refraction speeds 00097 dhdx=sign(1.d0,dhdx)*min(abs(dhdx),0.1d0) 00098 dhdy=sign(1.d0,dhdy)*min(abs(dhdy),0.1d0) 00099 if (par%wci==1) then 00100 call slope2D(s%u,s%nx,s%ny,s%dsu,s%dnv,dudx,dudy,s%wete) 00101 call slope2D(s%v,s%nx,s%ny,s%dsu,s%dnv,dvdx,dvdy,s%wete) 00102 else 00103 dudx = 0.d0 00104 dudy = 0.d0 00105 dvdx = 0.d0 00106 dvdy = 0.d0 00107 endif 00108 ! 00109 ! wwvv these slope routines are in wave_timestep, and are 00110 ! MPI-aware 00111 ! 00112 ! Calculate once sinh(2kh) 00113 where(2*s%hhw*s%k<=3000.d0) 00114 sinh2kh=sinh(min(2*s%k*s%hhw,10.0d0)) 00115 elsewhere 00116 sinh2kh = 3000.d0 00117 endwhere 00118 ! 00119 ! 00120 if (par%snells==0) then 00121 s%thetamean=(sum(s%ee*s%thet,3)/s%ntheta)/(max(sum(s%ee,3),0.00001d0)/s%ntheta) 00122 else !Dano: Snellius 00123 s%thetamean=asin(sin(s%theta0-s%alfaz(1,1))*s%c/s%c(1,1))+s%alfaz(1,1) 00124 s%costh(:,:,1)=cos(s%thetamean-s%alfaz) 00125 s%sinth(:,:,1)=sin(s%thetamean-s%alfaz) 00126 endif 00127 ! 00128 ! split wave velocities in wave grid directions theta 00129 call compute_wave_direction_velocities(s,par,0,dhdx,dhdy,dudx,dudy,dvdx,dvdy,sinh2kh) 00130 ! 00131 ! 00132 s%E(1,:)=max(sum(s%ee(1,:,:),2)*s%dtheta,0.0d0) 00133 s%R(1,:)=max(sum(s%rr(1,:,:),2)*s%dtheta,0.0d0) 00134 s%H(1,:)=sqrt(s%E(1,:)/par%rhog8) 00135 00136 imax=s%nx 00137 !Dano This is ok, since we will set mpiboundary to y in stationary mode 00138 ! 00139 ! write to screen that waves are updated 00140 if(xmaster) call writelog('ls','(a,f0.2,a)','Computing wave transformation at t = ',par%t,' s') 00141 call progress_indicator(.true.,0.d0,5.d0,2.d0) 00142 itermax = 0 00143 do i=2,imax 00144 call progress_indicator(.false.,dble(i)/imax*100,5.d0,2.d0) 00145 dtw=.5*minval(s%dsu(i:i+1,jmin_ee:jmax_ee))/max(maxval(s%cgx(i-1:i+1,jmin_ee:jmax_ee,:)),1d-10) 00146 dtw=min(dtw,.5*minval(s%dnv(i,jmin_ee:jmax_ee))/max(maxval(abs(s%cgy(i,jmin_ee:jmax_ee,:))),1d-10)) 00147 dtw=min(dtw,.5*s%dtheta/max(1.0d-6,maxval(abs(s%ctheta(i,jmin_ee:jmax_ee,:))))) 00148 !Dano: need to make sure all processes use the same dtw, min of all processes 00149 #ifdef USEMPI 00150 call xmpi_allreduce(dtw,MPI_MIN) 00151 #endif 00152 Herr=1. 00153 thetaerr = 2*par%px 00154 iter=0 00155 stopiterate=.false. 00156 do while (stopiterate .eqv. .false.) 00157 iter=iter+1 00158 Hprev=s%H(i,:) 00159 thetaprev = s%thetamean(i,:) 00160 ! 00161 ! transform to wave action 00162 ! 00163 i1=max(i-2,1) 00164 s%ee(i1:i+1,:,:) = s%ee(i1:i+1,:,:)/s%sigt(i1:i+1,:,:) 00165 ! 00166 ! Upwind Euler timestep propagation 00167 ! 00168 if (i>2 .and. (par%scheme==SCHEME_UPWIND_2 .or. par%scheme==SCHEME_WARMBEAM)) then 00169 ! note: at the moment Warm-Beam is currently not a valid combination (automatic change to upwind_2 in params.F90) 00170 call advecxho(s%ee(i-2:i+1,:,:),s%cgx(i-2:i+1,:,:),xadvec(i-2:i+1,:,:), & 00171 3,s%ny,s%ntheta,s%dnu(i-2:i+1,:),s%dsu(i-2:i+1,:),s%dsdnzi(i-2:i+1,:),par%scheme, & 00172 s%wete(i-2:i+1,:),dtw,s%dsz) 00173 else 00174 call advecxho(s%ee(i-1:i+1,:,:),s%cgx(i-1:i+1,:,:),xadvec(i-1:i+1,:,:), & 00175 2,s%ny,s%ntheta,s%dnu(i-1:i+1,:),s%dsu(i-1:i+1,:),s%dsdnzi(i-1:i+1,:),SCHEME_UPWIND_1, & 00176 s%wete(i-1:i+1,:),dtw,s%dsz) 00177 endif 00178 if (s%ny>0) then 00179 call advecyho(s%ee(i,:,:),s%cgy(i,:,:),yadvec(i,:,:), & ! RMc: why was this always set to upwind_1? 00180 0,s%ny,s%ntheta,s%dsv(i,:),s%dnv(i,:),s%dsdnzi(i,:),scheme_local_yadvec,s%wete(i,:),dtw,s%dnz) 00181 endif 00182 00183 call advecthetaho(s%ee(i,:,:),s%ctheta(i,:,:),thetaadvec(i,:,:),0,s%ny,s%ntheta,s%dtheta,par%scheme,s%wete(i,:)) 00184 00185 s%ee(i,:,:)=s%ee(i,:,:)-dtw*(xadvec(i,:,:) + yadvec(i,:,:) & 00186 + thetaadvec(i,:,:)) 00187 #ifdef USEMPI 00188 call xmpi_shift(s%ee(i-1:i,:,:),SHIFT_Y_R,1,2) 00189 call xmpi_shift(s%ee(i-1:i,:,:),SHIFT_Y_L,3,4) 00190 #endif 00191 ! 00192 ! transform back to wave energy 00193 ! 00194 s%ee(i1:i+1,:,:) = s%ee(i1:i+1,:,:)*s%sigt(i1:i+1,:,:) 00195 s%ee(i,:,:)=max(s%ee(i,:,:),0.0d0) 00196 00197 00198 ! 00199 ! Energy integrated over wave directions,Hrms 00200 ! 00201 s%E(i,:)=sum(s%ee(i,:,:),2)*s%dtheta 00202 s%H(i,:)=sqrt(s%E(i,:)/par%rhog8) 00203 do itheta=1,s%ntheta 00204 s%ee(i,:,itheta)=s%ee(i,:,itheta)/max(1.0d0,(s%H(i,:)/(par%gammax*s%hh(i,:)))**2) 00205 enddo 00206 s%H(i,:)=min(s%H(i,:),par%gammax*s%hh(i,:)) 00207 s%E(i,:)=par%rhog8*s%H(i,:)**2 00208 00209 if (par%snells==0) then !Dano not for SNellius 00210 s%thetamean(i,:) = (sum(s%ee(i,:,:)*s%thet(i,:,:),2)/s%ntheta)/(max(sum(s%ee(i,:,:),2),0.000010d0)/s%ntheta) 00211 endif 00212 ! 00213 ! Total dissipation 00214 00215 select case(par%break) 00216 case(BREAK_ROELVINK1,BREAK_ROELVINK2) 00217 call roelvink (par,s,i) 00218 case(BREAK_BALDOCK) 00219 call baldock (par,s,i) 00220 case(BREAK_JANSSEN) 00221 call janssen_battjes(par,s,i) 00222 end select 00223 00224 ! Dissipation by bed friction 00225 uorb(i,:)=par%px*s%H(i,:)/par%Trep/sinh(min(max(s%k(i,:),0.01d0)*s%hhw(i,:),10.0d0)) 00226 s%Df(i,:)=0.28d0*par%rho*s%fw(i,:)*uorb(i,:)**3 00227 where (s%hh>par%fwcutoff) 00228 s%Df = 0.d0 00229 end where 00230 ! 00231 ! Distribution of dissipation over directions and frequencies 00232 ! 00233 do itheta=1,s%ntheta 00234 dder(i,:,itheta)=s%ee(i,:,itheta)*s%D(i,:)/max(s%E(i,:),0.00001d0) 00235 ! Then all short wave energy dissipation, including bed friction and vegetation 00236 dd(i,:,itheta)=dder(i,:,itheta) + s%ee(i,:,itheta)*(s%Df(i,:)+s%Dveg(i,:))/max(s%E(i,:),0.00001d0) 00237 end do 00238 ! 00239 ! Euler step dissipation 00240 ! 00241 ! calculate roller energy balance 00242 ! 00243 ! RMc: changed to match ee balance 00244 if (i>2 .and. (par%scheme==SCHEME_UPWIND_2 .or. par%scheme==SCHEME_WARMBEAM)) then 00245 ! note: at the moment Warm-Beam is currently not a valid combination (automatic change to upwind_2 in params.F90) 00246 call advecxho(s%rr(i-2:i+1,:,:),s%cx(i-2:i+1,:,:),xradvec(i-2:i+1,:,:), & 00247 3,s%ny,s%ntheta,s%dnu(i-2:i+1,:),s%dsu(i-2:i+1,:),s%dsdnzi(i-2:i+1,:),par%scheme, & 00248 s%wete(i-2:i+1,:),dtw,s%dsz) 00249 else 00250 call advecxho(s%rr(i-1:i+1,:,:),s%cx(i-1:i+1,:,:),xradvec(i-1:i+1,:,:), & 00251 2,s%ny,s%ntheta,s%dnu(i-1:i+1,:),s%dsu(i-1:i+1,:),s%dsdnzi(i-1:i+1,:),SCHEME_UPWIND_1, & 00252 s%wete(i-1:i+1,:),dtw,s%dsz) 00253 endif 00254 if (s%ny>0) then 00255 call advecyho(s%rr(i,:,:),s%cy(i,:,:),yradvec(i,:,:), & 00256 0,s%ny,s%ntheta,s%dsv(i,:),s%dnv(i,:),s%dsdnzi(i,:),scheme_local_yadvec,s%wete(i,:),dtw,s%dnz) 00257 endif 00258 call advecthetaho(s%rr(i,:,:),s%ctheta(i,:,:),thetaradvec(i,:,:),0,s%ny,s%ntheta,s%dtheta,par%scheme,s%wete(i,:)) 00259 00260 s%rr(i,:,:)=s%rr(i,:,:)-dtw*(xradvec(i,:,:)+yradvec(i,:,:)+thetaradvec(i,:,:)) 00261 s%rr(i,:,:)=max(s%rr(i,:,:),0.0d0) 00262 #ifdef USEMPI 00263 call xmpi_shift(s%rr(i-1:i,:,:),SHIFT_Y_R,1,2) 00264 call xmpi_shift(s%rr(i-1:i,:,:),SHIFT_Y_L,3,4) 00265 #endif 00266 ! 00267 ! euler step roller energy dissipation (source and sink function) 00268 do j=jmin_ee,jmax_ee 00269 do itheta=1,s%ntheta 00270 if (dtw*dd(i,j,itheta)>s%ee(i,j,itheta)) then 00271 dtw=min(dtw,.5*s%ee(i,j,itheta)/dd(i,j,itheta)) 00272 endif 00273 enddo 00274 enddo 00275 !Dano: need to make sure all processes use the same dtw, min of all processes 00276 #ifdef USEMPI 00277 call xmpi_allreduce(dtw,MPI_MIN) 00278 #endif 00279 do j=1,s%ny+1 00280 do itheta=1,s%ntheta 00281 if(s%wete(i,j)==1) then 00282 s%ee(i,j,itheta)=s%ee(i,j,itheta)-dtw*dd(i,j,itheta) 00283 if (par%roller==1) then !Christophe 00284 s%rr(i,j,itheta)=s%rr(i,j,itheta)+dtw*dder(i,j,itheta)& 00285 -dtw*2.*par%g*par%beta*s%rr(i,j,itheta)& 00286 /sqrt(s%cx(i,j,itheta)**2+s%cy(i,j,itheta)**2) 00287 drr(i,j,itheta) = 2*par%g*par%beta*max(s%rr(i,j,itheta),0.0d0)/ & 00288 sqrt(s%cx(i,j,itheta)**2 +s%cy(i,j,itheta)**2) 00289 else 00290 s%rr(i,j,itheta)= 0.0d0 00291 drr(i,j,itheta)= 0.0d0 00292 end if 00293 s%ee(i,j,itheta)=max(s%ee(i,j,itheta),0.0d0) 00294 s%rr(i,j,itheta)=max(s%rr(i,j,itheta),0.0d0) 00295 else 00296 s%ee(i,j,itheta)=0.0d0 00297 s%rr(i,j,itheta)=0.0d0 00298 drr(i,j,itheta)=0.0d0 00299 end if 00300 end do 00301 end do 00302 ! Lateral boundary condition 00303 if (xmpi_isleft .and. s%ny>0) then 00304 do itheta=1,s%ntheta 00305 !if (s%sinth(i,1,itheta)>=0.) then 00306 s%ee(i,1,itheta)=s%ee(i,2,itheta) 00307 s%rr(i,1,itheta)=s%rr(i,2,itheta) 00308 !endif 00309 enddo 00310 s%k(:,1)=s%k(:,2) 00311 s%sigm(:,1)=s%sigm(:,2) 00312 endif 00313 if (xmpi_isright .and. s%ny>0) then 00314 do itheta=1,s%ntheta 00315 !if (s%sinth(i,s%ny+1,itheta)<=0.) then 00316 s%ee(i,s%ny+1,itheta)=s%ee(i,s%ny,itheta) 00317 s%rr(i,s%ny+1,itheta)=s%rr(i,s%ny,itheta) 00318 !endif 00319 end do 00320 s%k(:,s%ny+1)=s%k(:,s%ny) 00321 s%sigm(:,s%ny+1)=s%sigm(:,s%ny) 00322 endif 00323 ! 00324 ! Compute mean wave direction 00325 ! 00326 if (par%snells==0) then 00327 s%thetamean(i,:)=(sum(s%ee(i,:,:)*s%thet(i,:,:),2)/size(s%ee(i,:,:),2)) & 00328 /(max(sum(s%ee(i,:,:),2),0.000010d0) /size(s%ee(i,:,:),2)) 00329 endif 00330 ! 00331 ! Energy integrated over wave directions,Hrms 00332 ! 00333 s%E(i,:)=sum(s%ee(i,:,:),2)*s%dtheta 00334 s%R(i,:)=sum(s%rr(i,:,:),2)*s%dtheta 00335 s%DR(i,:)=sum(drr(i,:,:),2)*s%dtheta 00336 s%H(i,:)=sqrt(s%E(i,:)/par%rhog8) 00337 Herr=maxval(abs(Hprev(jmin_ee:jmax_ee)-s%H(i,jmin_ee:jmax_ee))) 00338 00339 thetaerr = 0.d0 00340 do j=jmin_ee,jmax_ee 00341 thetaerr_cyc = minval( (/ abs(thetaprev(j)-(s%thetamean(i,j)-2*par%px)) , & 00342 abs(thetaprev(j)-(s%thetamean(i,j) )) , & 00343 abs(thetaprev(j)-(s%thetamean(i,j)+2*par%px)) /) ) 00344 thetaerr = max(thetaerr,thetaerr_cyc) 00345 enddo 00346 00347 #ifdef USEMPI 00348 call xmpi_allreduce(Herr,MPI_MAX) 00349 call xmpi_allreduce(thetaerr,MPI_MAX) 00350 #endif 00351 ! Stopping criteria 00352 if (iter<par%maxiter) then 00353 if (Herr<par%maxerror .and. thetaerr< par%maxerror_angle) then 00354 stopiterate=.true. 00355 !if(xmaster) call writelog('ls','(a,i4,a,i4)','Wave propagation row ',i,', iteration ',iter) 00356 endif 00357 else 00358 stopiterate=.true. 00359 if(xmaster) call writelog('ls','(a,i4,a,i4,a,f5.4)','Wave propagation row ',i,', iteration ',iter,', error: ',Herr) 00360 endif 00361 enddo ! End while loop 00362 itermax = max(itermax,iter) 00363 enddo ! End do i=2:s%nx loop 00364 call writelog('ls','(a,i4)','Maximum number of iterations: ',itermax) 00365 00366 00367 s%ee(s%nx+1,:,:) = s%ee(s%nx,:,:) 00368 s%rr(s%nx+1,:,:) = s%rr(s%nx,:,:) 00369 s%E(s%nx+1,:) = s%E(s%nx,:) 00370 s%R(s%nx+1,:) = s%R(s%nx,:) 00371 s%DR(s%nx+1,:) = s%DR(s%nx,:) 00372 s%H(s%nx+1,:) = s%H(s%nx,:) 00373 s%k(s%nx+1,:) = s%k(s%nx,:) 00374 s%sigm(s%nx+1,:) = s%sigm(s%nx,:) 00375 s%cg(s%nx+1,:) = s%cg(s%nx,:) 00376 s%c(s%nx+1,:) = s%c(s%nx,:) 00377 s%thet(s%nx+1,:,:) = s%thet(s%nx,:,:) 00378 00379 ! 00380 ! Radiation stresses and forcing terms 00381 ! 00382 s%n=s%cg/max(s%c,1d-10) 00383 s%Sxx=(s%n*sum((1.d0+s%costh**2.d0)*s%ee,3)-.5d0*sum(s%ee,3))*s%dtheta 00384 s%Syy=(s%n*sum((1.d0+s%sinth**2.d0)*s%ee,3)-.5d0*sum(s%ee,3))*s%dtheta 00385 s%Sxy=s%n*sum(s%sinth*s%costh*s%ee,3)*s%dtheta 00386 00387 ! add roller contribution 00388 00389 s%Sxx = s%Sxx + sum((s%costh**2)*s%rr,3)*s%dtheta 00390 s%Syy = s%Syy + sum((s%sinth**2)*s%rr,3)*s%dtheta 00391 s%Sxy = s%Sxy + sum(s%sinth*s%costh*s%rr,3)*s%dtheta 00392 00393 if (s%ny>0) then 00394 do j=2,s%ny 00395 do i=1,s%nx 00396 s%Fx(i,j)=-(s%Sxx(i+1,j)-s%Sxx(i,j))/s%dsu(i,j) & 00397 -(s%Sxy(i,j+1)+s%Sxy(i+1,j+1)-s%Sxy(i,j-1)-s%Sxy(i+1,j-1))/ & 00398 (s%dnv(i,j-1)+s%dnv(i,j)+s%dnv(i+1,j-1)+s%dnv(i+1,j)) 00399 enddo 00400 enddo 00401 00402 do j=1,s%ny 00403 do i=2,s%nx 00404 s%Fy(i,j)=-(s%Syy(i,j+1)-s%Syy(i,j))/s%dnv(i,j) & 00405 -(s%Sxy(i+1,j)+s%Sxy(i+1,j+1)-s%Sxy(i-1,j)-s%Sxy(i-1,j+1))/ & 00406 (s%dsu(i-1,j)+s%dsu(i,j)+s%dsu(i-1,j+1)+s%dsu(i,j+1)) 00407 enddo 00408 enddo 00409 else 00410 j=1 00411 do i=1,s%nx 00412 s%Fx(i,j)=-(s%Sxx(i+1,j)-s%Sxx(i,j))/s%dsu(i,j) 00413 enddo 00414 do i=2,s%nx 00415 s%Fy(i,j)=-(s%Sxy(i+1,j)-s%Sxy(i-1,j))/ (s%dsu(i-1,j)+s%dsu(i,j)) 00416 enddo 00417 endif 00418 00419 if (s%ny>0) then 00420 s%Fx(:,1)=s%Fx(:,2) 00421 ! Fy(:,1)=Fy(:,2) 00422 s%Fx(:,s%ny+1)=s%Fx(:,s%ny) 00423 ! Fy(:,ny+1)=Fy(:,ny) 00424 endif 00425 s%Fx(1,:)=s%Fx(2,:) 00426 ! Fy(1,:)=Fy(2,:) 00427 00428 s%urms=uorb/sqrt(2.d0) 00429 00430 ! wave induced mass flux 00431 ustw=s%E*s%k/s%sigm/par%rho/max(s%hh,.001d0) ! Robert: why different from wave_instationary ??? 00432 s%uwf = ustw*dcos(s%thetamean-s%alfaz) 00433 s%vwf = ustw*dsin(s%thetamean-s%alfaz) 00434 00435 ! roller contribution 00436 s%ustr=2.*s%R*s%k/s%sigm/par%rho/max(s%hh,.001d0) ! Robert: why different from wave_instationary ??? 00437 00438 !ust = usd 00439 s%ust=s%usd+ustw ! Robert: why different from wave_instationary ??? 00440 !lateral boundaries 00441 s%ust(1,:) = s%ust(2,:) 00442 if (s%ny>0) then 00443 s%ust(:,1) = s%ust(:,2) 00444 s%ust(:,s%ny+1) = s%ust(:,s%ny) 00445 endif 00446 00447 end subroutine wave_stationary 00448 end module wave_stationary_module
1.7.6.1
