|
XBeach
|
00001 module wave_functions_module 00002 00003 00004 implicit none 00005 save 00006 private 00007 public:: dispersion, slope2d, breakerdelay, & 00008 & advecthetaho, advecxho, advecyho, iteratedispersion, advecqy, advecwy, & 00009 & advecwx, advecqx, wave_dispersion, update_means_wave_flow, compute_wave_direction_velocities, compute_stokes_drift, compute_wave_forces 00010 00011 contains 00012 00013 subroutine update_means_wave_flow(s,par) 00014 use params, only: parameters 00015 use spaceparams 00016 use paramsconst 00017 implicit none 00018 00019 type(spacepars), target :: s 00020 type(parameters) :: par 00021 real*8 :: factime 00022 00023 00024 ! For the "stationary" part: 00025 if (par%wavemodel == WAVEMODEL_SURFBEAT .and. par%single_dir == 1) then 00026 ! we need smoothed "stationary" values for the wave directions 00027 factime = 1.d0/(par%wavint*2)*par%dt 00028 s%hhws = max(factime*s%hhw + (1-factime)*s%hhws,par%eps) 00029 if (par%wci==1) then 00030 s%uws = factime*s%u + (1-factime)*s%uws 00031 s%vws = factime*s%v + (1-factime)*s%vws 00032 endif 00033 endif 00034 00035 ! for the instationary part 00036 if (par%wavemodel == WAVEMODEL_SURFBEAT .and. par%wci==1) then 00037 ! we need smoothed water depth and velocities for wci in wave instationary 00038 if (par%single_dir == 1) then 00039 ! maintain consistency with stationary wave directions model 00040 factime = 1.d0/(par%wavint*2)*par%dt 00041 else 00042 ! maintain consistency with boundary conditions smoothing 00043 factime = 1.d0/par%cats/par%Trep*par%dt 00044 endif 00045 s%hhwcins = max(factime*s%hhw + (1-factime)*s%hhwcins,par%eps) 00046 s%uwcins = factime*s%u + (1-factime)*s%uwcins 00047 s%vwcins = factime*s%v + (1-factime)*s%vwcins 00048 endif 00049 00050 end subroutine update_means_wave_flow 00051 00052 subroutine slope2D(h,nx,ny,dsu,dnv,dhdx,dhdy,wete) 00053 use xmpi_module 00054 implicit none 00055 00056 integer :: i,j,nx,ny 00057 real*8, dimension(nx+1,ny+1) :: h,dhdx,dhdy 00058 real*8, dimension(nx+1,ny+1) :: dsu 00059 real*8, dimension(nx+1,ny+1) :: dnv 00060 integer, dimension(nx+1,ny+1),intent(in) :: wete 00061 00062 00063 ! wwvv dhdx(2:nx,:) is computed, dhdx(1,:) and dhdx(nx+1,:) 00064 ! get boundary values, so in the parallel case, we need 00065 ! to do something about that: get the boundaries from 00066 ! upper and lower neighbours 00067 00068 ! u-gradients 00069 if(nx+1>=2) then 00070 forall(i=2:nx,j=1:ny+1,wete(i,j)==1) 00071 dhdx(i,j)=(h(i+1,j)-h(i-1,j))/(dsu(i,j)+dsu(i-1,j)) 00072 endforall 00073 forall(j=1:ny+1,wete(1,j)==1) 00074 dhdx(1,j)=(h(2,j)-h(1,j))/dsu(1,j) 00075 endforall 00076 forall(j=1:ny+1,wete(nx+1,j)==1) 00077 dhdx(nx+1,j)=(h(nx+1,j)-h(nx,j))/dsu(nx,j) 00078 endforall 00079 where(wete==0) 00080 dhdx=0.d0 00081 endwhere 00082 else 00083 dhdx=0.d0 00084 endif 00085 !#ifdef USEMPI 00086 ! call xmpi_shift(dhdx,SHIFT_X_U,3,4) 00087 ! call xmpi_shift(dhdx,SHIFT_X_D,1,2) 00088 !#endif 00089 ! 00090 ! v-gradients 00091 if(ny+1>=2) then 00092 forall(i=1:nx+1,j=2:ny,wete(i,j)==1) 00093 dhdy(i,j)=(h(i,j+1)-h(i,j-1))/(dnv(i,j)+dnv(i,j-1)) 00094 endforall 00095 forall(i=1:nx+1,wete(i,1)==1) 00096 dhdy(i,1)=(h(i,2)-h(i,1))/dnv(i,1) 00097 endforall 00098 forall(i=1:nx+1,wete(i,ny+1)==1) 00099 dhdy(i,ny+1)=(h(i,ny+1)-h(i,ny))/dnv(i,ny) 00100 endforall 00101 where(wete==0) 00102 dhdy=0.d0 00103 endwhere 00104 else 00105 dhdy=0.d0 00106 endif 00107 !#ifdef USEMPI 00108 ! call xmpi_shift(dhdy,SHIFT_Y_R,1,2) 00109 ! call xmpi_shift(dhdy,SHIFT_Y_L,3,4) 00110 !#endif 00111 00112 end subroutine slope2D 00113 00114 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00115 00116 subroutine advecxho(ee,cgx,xadvec,nx,ny,ntheta,dnu,dsu,dsdnzi,scheme,wete,dt,dsz) 00117 use spaceparams 00118 !use xmpi_module 00119 use paramsconst 00120 00121 implicit none 00122 00123 integer :: i,j,nx,ny,ntheta 00124 integer, intent(in) :: scheme 00125 integer, dimension(nx+1,ny+1),intent(in) :: wete 00126 real*8, intent(in) :: dt 00127 integer :: itheta 00128 real*8 , dimension(nx+1,ny+1) :: dnu,dsu,dsz,dsdnzi,fluxx 00129 real*8 , dimension(nx+1,ny+1,ntheta) :: xadvec,ee,cgx 00130 real*8 :: cgxu,eupw 00131 00132 integer :: scheme_now 00133 integer :: imin_local,imax_local 00134 00135 if(nx>3) then 00136 imin_local = imin_ee 00137 imax_local = imax_ee 00138 else 00139 ! Need nx here in case of stationary wave solver (nx=3 given in subroutine call; arrays from (i-2:i+1), need to fill 00140 ! 3rd place in the arrays (i) ) 00141 imin_local = nx 00142 imax_local = nx 00143 endif 00144 00145 xadvec = 0.d0 00146 fluxx = 0.d0 00147 00148 00149 ! split into schemes first, less split loops -> more efficiency 00150 scheme_now=scheme 00151 select case(scheme_now) 00152 case(SCHEME_UPWIND_1) 00153 do itheta=1,ntheta 00154 do j=1,ny+1 00155 do i=1,nx ! Whole domain 00156 if(wete(i,j)==1) then 00157 cgxu=.5*(cgx(i+1,j,itheta)+cgx(i,j,itheta)) 00158 if (cgxu>0) then 00159 fluxx(i,j)=ee(i,j,itheta)*cgxu*dnu(i,j) 00160 else 00161 fluxx(i,j)=ee(i+1,j,itheta)*cgxu*dnu(i,j) 00162 endif 00163 endif 00164 enddo 00165 enddo 00166 !do j=1,ny+1 ! 00167 do j=jmin_ee,jmax_ee 00168 do i=imin_local,imax_local 00169 if(wete(i,j)==1) then 00170 xadvec(i,j,itheta)=(fluxx(i,j)-fluxx(i-1,j))*dsdnzi(i,j) 00171 endif 00172 enddo 00173 enddo 00174 enddo 00175 case(SCHEME_UPWIND_2,SCHEME_WARMBEAM) 00176 do itheta=1,ntheta 00177 do j=1,ny+1 00178 do i=2,nx-1 00179 if(wete(i,j)==1) then 00180 cgxu=.5*(cgx(i+1,j,itheta)+cgx(i,j,itheta)) 00181 if (cgxu>0) then 00182 ! eupw=((dsu(i,j)+.5*dsu(i-1,j))*ee(i,j,itheta)-.5*dsu(i-1,j)*ee(i-1,j,itheta))/dsu(i-1,j) 00183 eupw=((dsu(i-1,j)+.5*dsu(i,j))*ee(i,j,itheta)-.5*dsu(i,j)*ee(i-1,j,itheta))/dsu(i-1,j) 00184 if (eupw<0.d0) eupw=ee(i,j,itheta) 00185 fluxx(i,j)=eupw*cgxu*dnu(i,j) 00186 else 00187 ! eupw=((dsu(i+1,j)+.5*dsu(i+2,j))*ee(i+1,j,itheta)-.5*dsu(i+2,j)*ee(i+2,j,itheta))/dsu(i+1,j) 00188 eupw=((dsu(i+1,j)+.5*dsu(i,j))*ee(i+1,j,itheta)-.5*dsu(i,j)*ee(i+2,j,itheta))/dsu(i+1,j) 00189 if (eupw<0.d0) eupw=ee(i+1,j,itheta) 00190 fluxx(i,j)=eupw*cgxu*dnu(i,j) 00191 endif 00192 endif 00193 enddo 00194 if (xmpi_istop) then 00195 i=1 ! only compute for i==1 00196 if(wete(i,j)==1) then 00197 cgxu=.5*(cgx(i+1,j,itheta)+cgx(i,j,itheta)) 00198 if (cgxu>0) then 00199 fluxx(i,j)=ee(i,j,itheta)*cgxu*dnu(i,j) 00200 else 00201 ! eupw=((dsu(i+1,j)+.5*dsu(i+2,j))*ee(i+1,j,itheta)-.5*dsu(i+2,j)*ee(i+2,j,itheta))/dsu(i+1,j) 00202 eupw=((dsu(i+1,j)+.5*dsu(i,j))*ee(i+1,j,itheta)-.5*dsu(i,j)*ee(i+2,j,itheta))/dsu(i+1,j) 00203 if (eupw<0.d0) eupw=ee(i+1,j,itheta) 00204 fluxx(i,j)=eupw*cgxu*dnu(i,j) 00205 endif 00206 endif 00207 endif 00208 if (xmpi_isbot) then 00209 i=nx ! only compute for i==nx0 00210 if(wete(i,j)==1) then 00211 cgxu=.5*(cgx(i+1,j,itheta)+cgx(i,j,itheta)) 00212 if (cgxu>0) then 00213 ! eupw=((dsu(i,j)+.5*dsu(i-1,j))*ee(i,j,itheta)-.5*dsu(i-1,j)*ee(i-1,j,itheta))/dsu(i-1,j) 00214 eupw=((dsu(i-1,j)+.5*dsu(i,j))*ee(i,j,itheta)-.5*dsu(i,j)*ee(i-1,j,itheta))/dsu(i-1,j) 00215 if (eupw<0.d0) eupw=ee(i,j,itheta) 00216 fluxx(i,j)=eupw*cgxu*dnu(i,j) 00217 else 00218 fluxx(i,j)=ee(i+1,j,itheta)*cgxu*dnu(i,j) 00219 endif 00220 endif 00221 endif 00222 enddo 00223 do j=jmin_ee,jmax_ee 00224 !do i=2,nx 00225 do i=imin_local,imax_local 00226 if(wete(i,j)==1) then 00227 xadvec(i,j,itheta)=(fluxx(i,j)-fluxx(i-1,j))*dsdnzi(i,j) 00228 endif 00229 enddo 00230 enddo 00231 enddo 00232 end select 00233 if (scheme_now==SCHEME_WARMBEAM) then 00234 do itheta=1,ntheta 00235 do j=jmin_ee,jmax_ee 00236 !do i=2,nx 00237 do i=imin_local,imax_local 00238 if(wete(i,j)==1) then 00239 xadvec(i,j,itheta)= xadvec(i,j,itheta) & 00240 -((ee(i+1,j,itheta)-ee(i ,j,itheta))/dsu(i ,j) & 00241 -(ee(i ,j,itheta)-ee(i-1,j,itheta))/dsu(i-1,j))/ & 00242 dsz(i,j)*dt/2*cgx(i,j,itheta)**2 00243 endif 00244 enddo 00245 enddo 00246 enddo 00247 endif 00248 00249 end subroutine advecxho 00250 00251 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00252 00253 subroutine advecthetaho(ee,ctheta,thetaadvec,nx,ny,ntheta,dtheta,scheme,wete) 00254 use spaceparams 00255 !use xmpi_module 00256 use paramsconst 00257 00258 implicit none 00259 00260 integer :: i,j,nx,ny,ntheta 00261 integer, intent(in) :: scheme 00262 integer, dimension(nx+1,ny+1),intent(in) :: wete 00263 integer :: itheta 00264 real*8 , dimension(ntheta) :: fluxtheta 00265 real*8 , dimension(nx+1,ny+1,ntheta) :: thetaadvec,ee,ctheta 00266 real*8 :: dtheta,ctheta_between,eupw 00267 00268 integer :: scheme_now 00269 integer :: imin_local,imax_local 00270 00271 if(nx>0) then 00272 imin_local = imin_ee 00273 imax_local = imax_ee 00274 else 00275 ! Need from 1 here in case of stationary wave solver (nx=0 given in subroutine call) 00276 imin_local = 1 00277 imax_local = 1 00278 endif 00279 00280 thetaadvec = 0.d0 00281 fluxtheta = 0.d0 00282 00283 ! No refraction caan take place if ntheta==1 00284 if (ntheta>1) then 00285 00286 ! split into schemes first, less split loops -> more efficiency 00287 scheme_now=scheme 00288 select case(scheme_now) 00289 case(SCHEME_UPWIND_1) 00290 do j=jmin_ee,jmax_ee 00291 do i=imin_local,imax_local 00292 if(wete(i,j)==1) then 00293 do itheta=1,ntheta-1 00294 ctheta_between=.5*(ctheta(i,j,itheta+1)+ctheta(i,j,itheta)) 00295 if (ctheta_between>0) then 00296 fluxtheta(itheta)=ee(i,j,itheta)*ctheta_between 00297 else 00298 fluxtheta(itheta)=ee(i,j,itheta+1)*ctheta_between 00299 endif 00300 enddo 00301 thetaadvec(i,j,1)=(fluxtheta(1)-0.d0)/dtheta ! No flux across lower boundary theta grid 00302 do itheta=2,ntheta-1 00303 thetaadvec(i,j,itheta)=(fluxtheta(itheta)-fluxtheta(itheta-1))/dtheta 00304 enddo 00305 thetaadvec(i,j,ntheta)=(0.d0-fluxtheta(ntheta-1))/dtheta ! No flux across upper boundary theta grid 00306 endif 00307 enddo 00308 enddo 00309 case(SCHEME_UPWIND_2,SCHEME_WARMBEAM) 00310 do j=jmin_ee,jmax_ee 00311 do i=imin_local,imax_local 00312 if(wete(i,j)==1) then 00313 do itheta=2,ntheta-2 00314 ctheta_between=.5*(ctheta(i,j,itheta+1)+ctheta(i,j,itheta)) 00315 if (ctheta_between>0) then 00316 eupw=1.5d0*ee(i,j,itheta)-.5*ee(i,j,itheta-1) 00317 if (eupw<0.d0) eupw=ee(i,j,itheta) 00318 fluxtheta(itheta)=eupw*ctheta_between 00319 else 00320 eupw=1.5d0*ee(i,j,itheta+1)-.5*ee(i,j,itheta+2) 00321 if (eupw<0.d0) eupw=ee(i,j,itheta+1) 00322 fluxtheta(itheta)=eupw*ctheta_between 00323 endif 00324 enddo 00325 itheta=1 ! only compute for itheta==1 00326 ctheta_between=.5*(ctheta(i,j,itheta+1)+ctheta(i,j,itheta)) 00327 if (ctheta_between>0) then 00328 fluxtheta(itheta)=ee(i,j,itheta)*ctheta_between 00329 else 00330 eupw=1.5d0*ee(i,j,itheta+1)-.5*ee(i,j,itheta+2) 00331 if (eupw<0.d0) eupw=ee(i,j,itheta+1) 00332 fluxtheta(itheta)=eupw*ctheta_between 00333 endif 00334 itheta=ntheta-1 ! only compute for itheta==ntheta-1 00335 ctheta_between=.5*(ctheta(i,j,itheta+1)+ctheta(i,j,itheta)) 00336 if (ctheta_between>0) then 00337 eupw=1.5d0*ee(i,j,itheta)-.5*ee(i,j,itheta-1) 00338 if (eupw<0.d0) eupw=ee(i,j,itheta) 00339 fluxtheta(itheta)=eupw*ctheta_between 00340 else 00341 eupw=ee(i,j,itheta+1) 00342 fluxtheta(itheta)=eupw*ctheta_between 00343 endif 00344 thetaadvec(i,j,1)=(fluxtheta(1)-0.d0)/dtheta ! No flux across lower boundary theta grid 00345 do itheta=2,ntheta-1 00346 thetaadvec(i,j,itheta)=(fluxtheta(itheta)-fluxtheta(itheta-1))/dtheta 00347 enddo 00348 thetaadvec(i,j,ntheta)=(0.d0-fluxtheta(ntheta-1))/dtheta ! No flux across upper boundary theta grid 00349 endif 00350 enddo 00351 enddo 00352 end select 00353 !if (scheme_now==SCHEME_WARMBEAM) then 00354 ! do itheta=2,ntheta-1 00355 ! do j=jmin_ee,jmax_ee 00356 ! do i=imin_local,imax_local 00357 ! if(wete(i,j)==1) then 00358 ! thetaadvec(i,j,itheta)= thetaadvec(i,j,itheta) & 00359 ! - (ee(i,j,itheta+1)-2*ee(i,j,itheta)+ee(i,j,itheta-1))/ & 00360 ! dtheta**2*dt/2*ctheta(i,j,itheta)**2 00361 ! endif 00362 ! enddo 00363 ! enddo 00364 ! enddo 00365 !endif 00366 00367 00368 endif !ntheta>1 00369 end subroutine advecthetaho 00370 00371 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00372 00373 subroutine advecyho(ee,cgy,yadvec,nx,ny,ntheta,dsv,dnv,dsdnzi,scheme,wete,dt,dnz) 00374 00375 use spaceparams 00376 use paramsconst 00377 00378 implicit none 00379 00380 integer :: i,j,nx,ny,ntheta 00381 integer, intent(in) :: scheme 00382 integer, dimension(nx+1,ny+1),intent(in) :: wete 00383 real*8, intent(in) :: dt 00384 integer :: itheta 00385 real*8 , dimension(nx+1,ny+1) :: dsv,dnv,dnz,dsdnzi,fluxy 00386 real*8 , dimension(nx+1,ny+1,ntheta) :: yadvec,ee,cgy 00387 real*8 :: cgyv,eupw 00388 00389 integer :: scheme_now 00390 integer :: imin_local,imax_local 00391 00392 if(nx>0) then 00393 imin_local = imin_ee 00394 imax_local = imax_ee 00395 else 00396 ! Need from 1 here in case of stationary wave solver (nx=0 given in subroutine call) 00397 imin_local = 1 00398 imax_local = 1 00399 endif 00400 00401 yadvec = 0.d0 00402 fluxy = 0.d0 00403 00404 ! split into schemes first, less split loops -> more efficiency 00405 scheme_now=scheme 00406 select case(scheme_now) 00407 case(SCHEME_UPWIND_1) 00408 do itheta=1,ntheta 00409 do j=1,ny 00410 do i=1,nx+1 ! Whole domain 00411 if(wete(i,j)==1) then 00412 cgyv=.5*(cgy(i,j+1,itheta)+cgy(i,j,itheta)) 00413 if (cgyv>0) then 00414 fluxy(i,j)=ee(i,j,itheta)*cgyv*dsv(i,j) 00415 else 00416 fluxy(i,j)=ee(i,j+1,itheta)*cgyv*dsv(i,j) 00417 endif 00418 endif 00419 enddo 00420 enddo 00421 do j=jmin_ee,jmax_ee 00422 do i=imin_local,imax_local 00423 if(wete(i,j)==1) then 00424 yadvec(i,j,itheta)=(fluxy(i,j)-fluxy(i,j-1))*dsdnzi(i,j) 00425 endif 00426 enddo 00427 enddo 00428 enddo 00429 case(SCHEME_UPWIND_2,SCHEME_WARMBEAM) 00430 do itheta=1,ntheta 00431 do j=2,ny-1 00432 do i=1,nx+1 00433 if(wete(i,j)==1) then 00434 cgyv=.5*(cgy(i,j+1,itheta)+cgy(i,j,itheta)) 00435 if (cgyv>0) then 00436 ! eupw=((dnv(i,j)+.5*dnv(i,j-1))*ee(i,j,itheta)-.5*dnv(i,j-1)*ee(i,j-1,itheta))/dnv(i,j-1) 00437 eupw=((dnv(i,j-1)+.5*dnv(i,j))*ee(i,j,itheta)-.5*dnv(i,j)*ee(i,j-1,itheta))/dnv(i,j-1) 00438 if (eupw<0.d0) eupw=ee(i,j,itheta) 00439 fluxy(i,j)=eupw*cgyv*dsv(i,j) 00440 else 00441 ! eupw=((dnv(i,j+1)+.5*dnv(i,j+2))*ee(i,j+1,itheta)-.5*dnv(i,j+2)*ee(i,j+2,itheta))/dnv(i,j+1) 00442 eupw=((dnv(i,j+1)+.5*dnv(i,j))*ee(i,j+1,itheta)-.5*dnv(i,j)*ee(i,j+2,itheta))/dnv(i,j+1) 00443 if (eupw<0.d0) eupw=ee(i,j+1,itheta) 00444 fluxy(i,j)=eupw*cgyv*dsv(i,j) 00445 endif 00446 endif 00447 enddo 00448 enddo 00449 if (xmpi_isleft) then 00450 j=1 ! only compute for j==1 00451 do i=1,nx+1 00452 if(wete(i,j)==1) then 00453 cgyv=.5*(cgy(i,j+1,itheta)+cgy(i,j,itheta)) 00454 if (cgyv>0) then 00455 fluxy(i,j)=ee(i,j,itheta)*cgyv*dsv(i,j) 00456 else 00457 ! eupw=((dnv(i,j+1)+.5*dnv(i,j+2))*ee(i,j+1,itheta)-.5*dnv(i,j+2)*ee(i,j+2,itheta))/dnv(i,j+1) 00458 eupw=((dnv(i,j+1)+.5*dnv(i,j))*ee(i,j+1,itheta)-.5*dnv(i,j)*ee(i,j+2,itheta))/dnv(i,j+1) 00459 if (eupw<0.d0) eupw=ee(i,j+1,itheta) 00460 fluxy(i,j)=eupw*cgyv*dsv(i,j) 00461 endif 00462 endif 00463 enddo 00464 endif 00465 if (xmpi_isright) then 00466 j=ny ! only compute for j==ny 00467 do i=1,nx+1 00468 if(wete(i,j)==1) then 00469 cgyv=.5*(cgy(i,j+1,itheta)+cgy(i,j,itheta)) 00470 if (cgyv>0) then 00471 ! eupw=((dnv(i,j)+.5*dnv(i,j-1))*ee(i,j,itheta)-.5*dnv(i,j-1)*ee(i,j-1,itheta))/dnv(i,j-1) 00472 eupw=((dnv(i,j-1)+.5*dnv(i,j))*ee(i,j,itheta)-.5*dnv(i,j)*ee(i,j-1,itheta))/dnv(i,j-1) 00473 if (eupw<0.d0) eupw=ee(i,j,itheta) 00474 fluxy(i,j)=eupw*cgyv*dsv(i,j) 00475 else 00476 fluxy(i,j)=ee(i,j+1,itheta)*cgyv*dsv(i,j) 00477 endif 00478 endif 00479 enddo 00480 endif 00481 do j=2,ny 00482 do i=imin_local,imax_local 00483 if(wete(i,j)==1) then 00484 yadvec(i,j,itheta)=(fluxy(i,j)-fluxy(i,j-1))*dsdnzi(i,j) 00485 endif 00486 enddo 00487 enddo 00488 enddo 00489 end select 00490 if (scheme_now==SCHEME_WARMBEAM) then 00491 do itheta=1,ntheta 00492 do j=jmin_ee,jmax_ee 00493 do i=imin_local,imax_local 00494 if(wete(i,j)==1) then 00495 yadvec(i,j,itheta) = yadvec(i,j,itheta) & 00496 -((ee(i,j+1,itheta)-ee(i,j ,itheta))/dnv(i,j ) & 00497 -(ee(i,j ,itheta)-ee(i,j-1,itheta))/dnv(i,j-1))/ & 00498 dnz(i,j)*dt/2*cgy(i,j,itheta)**2 00499 endif 00500 enddo 00501 enddo 00502 enddo 00503 endif 00504 end subroutine advecyho 00505 00506 00507 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00508 00509 subroutine advecwx(arrin2d,xwadvec,kmx,nx,ny,dsu,wete) 00510 use xmpi_module 00511 00512 implicit none 00513 00514 integer :: i,nx,ny 00515 integer :: j 00516 real*8 , dimension(nx+1,ny+1) :: dsu 00517 real*8 , dimension(nx+1,ny+1) :: xwadvec,arrin2d,kmx 00518 integer, dimension(nx+1,ny+1),intent(in) :: wete 00519 00520 xwadvec = 0.d0 00521 00522 do j=2,ny 00523 do i=2,nx 00524 if(wete(i,j)==1) then 00525 if (kmx(i,j)>0) then 00526 xwadvec(i,j)=(arrin2d(i,j)-arrin2d(i-1,j))/dsu(i-1,j) 00527 elseif (kmx(i,j)<0) then 00528 xwadvec(i,j)=(arrin2d(i+1,j)-arrin2d(i,j))/dsu(i,j) 00529 else 00530 xwadvec(i,j)=(arrin2d(i+1,j)-arrin2d(i-1,j))/(dsu(i,j)+dsu(i-1,j)) 00531 endif 00532 endif 00533 end do 00534 end do 00535 00536 ! wwvv here we miss the computations of the first and last columns and rows, 00537 ! in the parallel case we shift these in form neighbours 00538 00539 end subroutine advecwx 00540 00541 00542 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1 00543 00544 subroutine advecwy(arrin2d,ywadvec,kmy,nx,ny,dnv,wete) 00545 use xmpi_module 00546 use xmpi_module 00547 implicit none 00548 00549 integer :: i,nx,ny 00550 integer :: j 00551 real*8 , dimension(nx+1,ny+1) :: dnv 00552 real*8 , dimension(nx+1,ny+1) :: ywadvec,arrin2d,kmy 00553 integer, dimension(nx+1,ny+1),intent(in) :: wete 00554 00555 ywadvec = 0.d0 00556 00557 do j=2,ny 00558 do i=2,nx 00559 if(wete(i,j)==1) then 00560 if (kmy(i,j)>0) then 00561 ywadvec(i,j)=(arrin2d(i,j)-arrin2d(i,j-1))/dnv(i,j-1) 00562 elseif (kmy(i,j)<0) then 00563 ywadvec(i,j)=(arrin2d(i,j+1)-arrin2d(i,j))/dnv(i,j) 00564 else 00565 ywadvec(i,j)=(arrin2d(i,j+1)-arrin2d(i,j-1))/(dnv(i,j)+dnv(i,j-1)) 00566 endif 00567 endif 00568 end do 00569 end do 00570 00571 if(ny>0) then 00572 where(wete(:,1)==1) 00573 ywadvec(:,1)= ywadvec(:,2) !Ap 00574 endwhere 00575 where(wete(:,ny+1)==1) 00576 ywadvec(:,ny+1) = ywadvec(:,ny) !Ap 00577 endwhere 00578 endif 00579 00580 00581 end subroutine advecwy 00582 00583 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00584 00585 subroutine advecqx(c,arrin2d,xwadvec,nx,ny,dsu,wete) 00586 use xmpi_module 00587 00588 implicit none 00589 00590 integer :: i,nx,ny 00591 integer :: j 00592 real*8 , dimension(nx+1,ny+1) :: dsu 00593 real*8 , dimension(nx+1,ny+1) :: xwadvec,arrin2d,c 00594 integer, dimension(nx+1,ny+1),intent(in) :: wete 00595 00596 xwadvec = 0.d0 00597 00598 do j=2,ny 00599 do i=2,nx 00600 if(wete(i,j)==1) then 00601 if (c(i,j)>0) then 00602 xwadvec(i,j)=c(i,j)*(arrin2d(i,j)-arrin2d(i-1,j))/dsu(i-1,j) 00603 elseif (c(i,j)<0) then 00604 xwadvec(i,j)=c(i,j)*(arrin2d(i+1,j)-arrin2d(i,j))/dsu(i,j) 00605 else 00606 xwadvec(i,j)=c(i,j)*(arrin2d(i+1,j)-arrin2d(i-1,j))/(dsu(i,j)+dsu(i-1,j)) 00607 endif 00608 endif 00609 end do 00610 end do 00611 00612 if(ny>0) then 00613 where(wete(:,1)==1) 00614 xwadvec(:,1)= xwadvec(:,2) !Ap 00615 endwhere 00616 where(wete(:,ny+1)==1) 00617 xwadvec(:,ny+1) = xwadvec(:,ny) !Ap 00618 endwhere 00619 endif 00620 00621 00622 end subroutine advecqx 00623 00624 00625 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1 00626 00627 subroutine advecqy(c,arrin2d,ywadvec,nx,ny,dnv,wete) 00628 use xmpi_module 00629 use xmpi_module 00630 implicit none 00631 00632 integer :: i,nx,ny 00633 integer :: j 00634 real*8 , dimension(nx+1,ny+1) :: dnv 00635 real*8 , dimension(nx+1,ny+1) :: ywadvec,arrin2d,c 00636 integer, dimension(nx+1,ny+1),intent(in) :: wete 00637 00638 ywadvec = 0.d0 00639 00640 do j=2,ny 00641 do i=2,nx 00642 if(wete(i,j)==1) then 00643 if (c(i,j)>0) then 00644 ywadvec(i,j)=c(i,j)*(arrin2d(i,j)-arrin2d(i,j-1))/dnv(i,j-1) 00645 elseif (c(i,j)<0) then 00646 ywadvec(i,j)=c(i,j)*(arrin2d(i,j+1)-arrin2d(i,j))/dnv(i,j) 00647 else 00648 ywadvec(i,j)=c(i,j)*(arrin2d(i,j+1)-arrin2d(i,j-1))/(dnv(i,j)+dnv(i,j-1)) 00649 endif 00650 endif 00651 end do 00652 end do 00653 00654 end subroutine advecqy 00655 00656 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00657 00658 subroutine dispersion(par,s,h) 00659 use params, only: parameters 00660 use spaceparams 00661 use logging_module, only: writelog 00662 00663 ! Robert: iteration along L=L0tanh(2pih/L) 00664 00665 implicit none 00666 00667 type(spacepars) :: s 00668 type(parameters) :: par 00669 real*8, dimension(1:s%nx+1,1:s%ny+1),intent(in) :: h ! water depth for dispersion can vary for stationary, single_dir and 00670 ! instationary computations, with and without wci 00671 00672 ! Robert: for some reason, MPI version does not like the "allocatable" version 00673 ! of this subroutine. 00674 real*8, dimension(1:s%nx+1,1:s%ny+1) :: L0,kh,Ltemp 00675 integer :: i,j,j1,j2 00676 real*8 :: backdis,disfac 00677 integer :: index 00678 00679 if (s%ny==0) then 00680 j1=1 00681 j2=1 00682 else 00683 j1=1 00684 j2=s%ny+1 00685 endif 00686 ! 00687 ! 00688 where(s%wete==1) 00689 L0 = par%g*par%Trep**2/(2*par%px) 00690 elsewhere 00691 L0=par%eps 00692 endwhere 00693 00694 if (.not. associated(s%L1)) then 00695 allocate(s%L1(s%nx+1,s%ny+1)) 00696 s%L1=L0 00697 endif 00698 Ltemp = L0 00699 00700 do j = j1,j2 00701 do i = 1,s%nx+1 00702 if(s%wete(i,j)==1) then 00703 Ltemp(i,j) = iteratedispersion(L0(i,j),Ltemp(i,j),par%px,h(i,j)) 00704 if (Ltemp(i,j)<0.d0) then ! this is an error from iteratedispersion 00705 Ltemp(i,j) = -Ltemp(i,j) 00706 call writelog('lws','','Warning: no convergence in dispersion relation iteration at t = ', & 00707 par%t*max(par%morfac*par%morfacopt,1.d0)) 00708 endif 00709 endif 00710 end do 00711 end do 00712 if (par%shoaldelay==1) then 00713 ! find Lmod looking back over distance par%facsd*L1 00714 ! presumes sigma direction is shore normal 00715 s%L1 = 0.d0 ! modified wave length, initially set to L1 00716 do j = j1,j2 00717 do i = 2,s%nx+1 00718 if(s%wete(i,j)==1) then 00719 index = i ! start index 00720 backdis = 0.d0 ! relative distance backward 00721 do while (backdis<1.d0) 00722 ! use average wavelength over distance dsc 00723 disfac = s%dsc(index,j)/(par%facsd*0.5d0*(Ltemp(index,j)+Ltemp(max(index-1,1),j))) 00724 disfac = min(disfac,1.d0-backdis) 00725 00726 00727 s%L1(i,j) = s%L1(i,j)+disfac*0.5d0*(Ltemp(index,j)+Ltemp(max(index-1,1),j)) 00728 backdis = backdis+disfac 00729 00730 index = max(index-1,1) 00731 enddo 00732 endif 00733 enddo 00734 enddo 00735 where(s%wete(1,:)==1) 00736 s%L1(1,:) = Ltemp(1,:) 00737 endwhere 00738 else 00739 where(s%wete==1) 00740 s%L1 = Ltemp 00741 endwhere 00742 endif 00743 00744 ! boundary copies for non superfast 1D 00745 if (s%ny>0) then 00746 if (xmpi_isleft) then 00747 where(s%wete(:,1)==1) 00748 s%L1(:,1)=s%L1(:,2) 00749 endwhere 00750 endif 00751 if (xmpi_isright) then 00752 where(s%wete(:,s%ny+1)==1) 00753 s%L1(:,s%ny+1)=s%L1(:,s%ny) 00754 endwhere 00755 endif 00756 endif 00757 where(s%wete==1) 00758 s%k = 2*par%px/s%L1 00759 s%c = s%sigm/s%k 00760 kh = min(s%k*h,10.0d0) 00761 s%n=0.5d0+kh/sinh(2*kh) 00762 s%cg=s%c*s%n 00763 elsewhere 00764 s%k = 25.d0 00765 s%c = sqrt(par%g*par%eps) 00766 s%n = 1.d0 00767 s%cg= sqrt(par%g*par%eps) 00768 endwhere 00769 00770 00771 end subroutine dispersion 00772 00773 elemental function iteratedispersion(L0,Lestimate,px,h) result(L) 00774 00775 implicit none 00776 ! input 00777 real*8,intent(in) :: L0 00778 real*8,intent(in) :: Lestimate 00779 real*8,intent(in) :: px 00780 real*8,intent(in) :: h 00781 ! output 00782 real*8 :: L 00783 ! internal 00784 real*8 :: L1,L2 00785 integer :: iter 00786 real*8 :: err 00787 real*8,parameter :: aphi = 1.d0/(((1.0d0 + sqrt(5.0d0))/2)+1) 00788 real*8,parameter :: bphi = ((1.0d0 + sqrt(5.0d0))/2)/(((1.0d0 + sqrt(5.0d0))/2)+1) 00789 integer,parameter :: itermax = 150 00790 real*8,parameter :: errmax = 0.00001d0 00791 00792 00793 err = huge(0.0d0) 00794 iter = 0 00795 L1 = Lestimate 00796 do while (err > errmax .and. iter < itermax) 00797 iter = iter+1 00798 L2 = L0*tanh(2*px*h/L1) 00799 L1 = (L1*aphi + L2*bphi) ! Golden ratio 00800 err = abs(L2 - L1) 00801 end do 00802 00803 if (iter<=itermax) then 00804 L = L1 00805 else 00806 ! signal this went wrong 00807 L = -L1 00808 endif 00809 00810 end function iteratedispersion 00811 00812 subroutine breakerdelay(par,s) 00813 !!! Dano: this routine is not MPI compliant and assumes xz to be in cross-shore; TBD 00814 !!! -------------------------------------------------------------------------------- 00815 use params, only: parameters 00816 use spaceparams 00817 use xmpi_module 00818 00819 implicit none 00820 00821 type(spacepars),target :: s 00822 type(parameters) :: par 00823 00824 real*8 :: Lbr 00825 real*8, dimension(s%nx+1) :: utemp 00826 integer :: jx,jy,i,j1,nbr,tempxid 00827 integer, dimension(s%nx+1) :: ibr 00828 00829 !include 's.ind' 00830 !include 's.inp' 00831 00832 ! Superfast 1D 00833 if (s%ny>0) then 00834 j1 = 2 00835 else 00836 j1 = 1 00837 endif 00838 00839 do jy = j1,max(1,s%ny) 00840 s%usd(1,jy) = s%ustr(1,jy) 00841 00842 do jx = 2,s%nx+1 00843 if(s%wete(jx,jy)==1) then 00844 nbr = 0 00845 Lbr = sqrt(par%g*s%hhw(jx,jy))*par%Trep*par%breakerdelay 00846 i = jx-1 00847 do while (abs(s%xz(i,jy)-s%xz(jx,jy))<=Lbr .and. i>1) 00848 nbr = nbr+1 00849 i = i-1 00850 end do 00851 00852 if(nbr.gt.1) then 00853 do i = 1,nbr+1 00854 ibr(i) = i 00855 tempxid = jx-nbr+i-1 00856 utemp(i) = s%ustr(tempxid,jy) 00857 enddo 00858 00859 s%usd(jx,jy) = sum(ibr(1:nbr+1)*utemp(1:nbr+1))/sum(ibr(1:nbr+1)) 00860 else 00861 s%usd(jx,jy) = s%ustr(jx,jy) 00862 end if 00863 endif 00864 end do 00865 end do 00866 00867 ! lateral boundaries 00868 if (xmpi_istop ) then 00869 where(s%wete(1,:)==1) 00870 s%usd(1,:) = s%usd(2,:) 00871 endwhere 00872 endif 00873 if (xmpi_isbot ) then 00874 where(s%wete(s%nx+1,:)==1) 00875 s%usd(s%nx+1,:) = s%usd(s%nx,:) 00876 endwhere 00877 endif 00878 if (xmpi_isleft .and. s%ny>0) then 00879 where(s%wete(:,1)==1) 00880 s%usd(:,1) = s%usd(:,2) 00881 endwhere 00882 endif 00883 if (xmpi_isright .and. s%ny>0) then 00884 where(s%wete(:,s%ny+1)==1) 00885 s%usd(:,s%ny+1) = s%usd(:,s%ny) 00886 endwhere 00887 endif 00888 00889 ! wwvv for the parallel version, shift in the columns and rows 00890 00891 end subroutine breakerdelay 00892 00893 subroutine wave_dispersion(s,par,useAverageDepthSwitch) 00894 00895 use params, only: parameters 00896 use spaceparams 00897 00898 implicit none 00899 00900 type(spacepars), target :: s 00901 type(parameters) :: par 00902 integer,intent(in),optional :: useAverageDepthSwitch 00903 00904 real*8,dimension(:,:),allocatable,save :: km,kmx,kmy 00905 real*8,dimension(:,:),allocatable,save :: hhlocal,ulocal,vlocal,relangle 00906 real*8,dimension(:,:),allocatable,save :: arg,fac 00907 real*8,dimension(:,:),allocatable,save :: cgym,cgxm 00908 real*8,dimension(:,:),allocatable,save :: dkmxdx,dkmxdy,dkmydx,dkmydy 00909 real*8,dimension(:,:),allocatable,save :: xwadvec,ywadvec 00910 real*8,dimension(:),allocatable,save :: L0,L1 00911 real*8,save :: Trepold 00912 integer :: itheta 00913 integer :: luseAverageDepthSwitch 00914 00915 if(.not.allocated(km)) then 00916 allocate(km(s%nx+1,s%ny+1)) 00917 allocate(kmx(s%nx+1,s%ny+1)) 00918 allocate(kmy(s%nx+1,s%ny+1)) 00919 allocate(arg(s%nx+1,s%ny+1)) 00920 allocate(fac(s%nx+1,s%ny+1)) 00921 allocate(cgym(s%nx+1,s%ny+1)) 00922 allocate(cgxm(s%nx+1,s%ny+1)) 00923 allocate(dkmxdx(s%nx+1,s%ny+1)) 00924 allocate(dkmxdy(s%nx+1,s%ny+1)) 00925 allocate(dkmydx(s%nx+1,s%ny+1)) 00926 allocate(dkmydy(s%nx+1,s%ny+1)) 00927 allocate(xwadvec(s%nx+1,s%ny+1)) 00928 allocate(ywadvec(s%nx+1,s%ny+1)) 00929 allocate(hhlocal(s%nx+1,s%ny+1)) 00930 if (par%wci==1) then 00931 allocate(ulocal(s%nx+1,s%ny+1)) 00932 allocate(vlocal(s%nx+1,s%ny+1)) 00933 allocate(relangle(s%nx+1,s%ny+1)) 00934 allocate(L0(s%ny+1)) 00935 allocate(L1(s%ny+1)) 00936 L0 = 0.d0 00937 L1 = -huge(0.d0) 00938 endif 00939 Trepold = 0.d0 00940 call dispersion(par,s,s%hhw) ! at initialisation, water depth is always s%hhw 00941 km=s%k 00942 endif 00943 00944 ! water depth and velocities to use depends on wave mode and presence of wci 00945 00946 ! default to use instantaneous water depth and velocities 00947 if (present(useAverageDepthSwitch)) then 00948 luseAverageDepthSwitch = useAverageDepthSwitch 00949 else 00950 luseAverageDepthSwitch = 0 00951 endif 00952 00953 if (luseAverageDepthSwitch==0) then ! default: use instantaneous depth and velocities 00954 hhlocal = s%hhw 00955 if (par%wci==1) then 00956 ulocal = s%u 00957 vlocal = s%v 00958 endif 00959 elseif (luseAverageDepthSwitch==1) then ! used averaged depths for single_dir computation 00960 hhlocal = s%hhws 00961 if (par%wci==1) then 00962 ulocal = s%uws 00963 vlocal = s%vws 00964 endif 00965 elseif (luseAverageDepthSwitch==2) then ! used averaged depths for instationary computation with wci 00966 hhlocal = s%hhwcins 00967 ulocal = s%uwcins 00968 vlocal = s%vwcins 00969 else 00970 ! this should never occur 00971 endif 00972 00973 00974 if(par%wci==1) then 00975 ! Dano NEED TO CHECK THIS FOR CURVI 00976 !if (xmpi_istop) then 00977 ! ! requires boundary condition at the offshore boundary, where we assume zero flow 00978 ! L0 = par%g*par%Trep**2/(2*par%px) 00979 ! if (any(L1<=0.d0)) then 00980 ! L1 = L0 00981 ! endif 00982 ! do j=1,s%ny+1 00983 ! L1(j) = iteratedispersion(L0(j),L1(j),par%px,s%hhw(1,j)) 00984 ! enddo 00985 ! km(1,:) = 2*par%px/max(L1,0.00001d0) 00986 !endif 00987 arg = min(100.0d0,km*hhlocal) 00988 fac = ( 1.d0 + ((km*s%H/2.d0)**2)) ! use deep water correction 00989 s%sigm(1,:) = sqrt( par%g*km(1,:)*tanh(arg(1,:))) ! *( 1.d0+ ((km(1,:)*s%H(1,:)/2.d0)**2))) 00990 ! calculate change in intrinsic frequency 00991 relangle = s%thetamean-s%alfaz 00992 kmx = km*dcos(relangle) 00993 kmy = km*dsin(relangle) 00994 s%wm = s%sigm+kmx*ulocal*min(& 00995 min(hhlocal/par%hwci,1.d0), & 00996 min(1.d0,(1.d0-hhlocal/par%hwcimax)) & 00997 )+ & 00998 kmy*vlocal*min( & 00999 min(hhlocal/par%hwci,1.d0), & 01000 min(1.d0,(1.d0-hhlocal/par%hwcimax)) & 01001 ) 01002 01003 where(km>0.01d0) 01004 s%c = s%sigm/km 01005 ! cg = c*(0.5d0+arg/sinh(2.0d0*arg)) ! Linear 01006 s%cg = s%c*(0.5d0+arg/sinh(2*arg))*sqrt(fac) ! & to include more 01007 ! + km*(H/2)**2*sqrt(max(par%g*km*tanh(arg),0.001d0))/sqrt(max(fac,0.001d0)) ! include wave steepness 01008 s%n=0.5d0+km*hhlocal/sinh(2*max(km,0.00001d0)*hhlocal) 01009 elsewhere 01010 s%c = 0.01d0 01011 s%cg = 0.01d0 01012 s%n = 1.d0 01013 endwhere 01014 01015 cgym = s%cg*dsin(relangle)+vlocal*min(min(hhlocal/par%hwci,1.d0),min(1.d0,(1.d0-hhlocal/par%hwcimax))) 01016 cgxm = s%cg*dcos(relangle)+ulocal*min(min(hhlocal/par%hwci,1.d0),min(1.d0,(1.d0-hhlocal/par%hwcimax))) 01017 01018 call slope2D(kmx,s%nx,s%ny,s%dsu,s%dnv,dkmxdx,dkmxdy,s%wete) 01019 call slope2D(kmy,s%nx,s%ny,s%dsu,s%dnv,dkmydx,dkmydy,s%wete) 01020 call advecwx(s%wm,xwadvec,kmx,s%nx,s%ny,s%dsu,s%wete) ! cjaap: s%xz or s%xu? kmx = kmx -par%dt*xwadvec -par%dt*cgym*(dkmydx-dkmxdy) 01021 kmx = kmx -par%dt*xwadvec -par%dt*cgym*(dkmydx-dkmxdy) 01022 if (s%ny>0) then 01023 if (xmpi_isright) kmx(:,s%ny+1) = kmx(:,s%ny) ! lateral bc 01024 if (xmpi_isleft) kmx(:,1) = kmx(:,2) ! lateral bc 01025 endif 01026 01027 call advecwy(s%wm,ywadvec,kmy,s%nx,s%ny,s%dnv,s%wete) ! cjaap: s%yz or s%yv? 01028 kmy = kmy-par%dt*ywadvec + par%dt*cgxm*(dkmydx-dkmxdy) 01029 if (s%ny>0) then 01030 if (xmpi_isright) kmy(:,s%ny+1) = kmy(:,s%ny) ! lateral bc 01031 if (xmpi_isleft) kmy(:,1) = kmy(:,2) ! lateral bc 01032 endif 01033 01034 #ifdef USEMPI 01035 call xmpi_shift_ee(kmx) 01036 call xmpi_shift_ee(kmy) 01037 #endif 01038 01039 ! update km 01040 km = sqrt(kmx**2+kmy**2) 01041 km = min(km,25.d0) ! limit to gravity waves 01042 ! non-linear dispersion 01043 arg = min(100.0d0,km*hhlocal) 01044 arg = max(arg,0.0001) 01045 ! fac = ( 1.d0 + ((km*H/2.d0)**2)*( (8.d0+(cosh(min(4.d0*arg,10.0d0)))**1.d0-2.d0*(tanh(arg))**2.d0 ) /(8.d0*(sinh(arg))**4.d0) ) ) 01046 fac = ( 1.d0 + ((km*s%H/2.d0)**2)) ! use deep water correction instead of expression above (waves are short near blocking point anyway) 01047 ! fac = 1.d0 ! Linear 01048 s%sigm = sqrt( par%g*km*tanh(arg)*fac) 01049 s%sigm = max(s%sigm,0.010d0) 01050 ! update intrinsic frequency 01051 do itheta=1,s%ntheta 01052 s%sigt(:,:,itheta) = s%sigm 01053 enddo 01054 ! 01055 ! update k 01056 s%k = km 01057 else 01058 ! check if we need to recompute sigm and sigt 01059 if (abs(par%Trep-Trepold)/par%Trep > 0.0001d0) then 01060 Trepold = par%Trep 01061 do itheta=1,s%ntheta 01062 s%sigt(:,:,itheta) = 2*par%px/par%Trep 01063 end do 01064 s%sigm = 2*par%px/par%Trep 01065 endif 01066 call dispersion(par,s,hhlocal) 01067 endif ! end wave current interaction 01068 01069 01070 end subroutine wave_dispersion 01071 01072 subroutine compute_wave_direction_velocities(s,par,flavour,dhdx,dhdy,dudx,dudy,dvdx,dvdy,sinh2kh) 01073 01074 use params, only: parameters 01075 use spaceparams 01076 01077 implicit none 01078 01079 type(spacepars), target :: s 01080 type(parameters) :: par 01081 integer,intent(in) :: flavour ! 0 for stationary, 1 for instationary, 2 for directions part of single_dir 01082 real*8,dimension(s%nx+1,s%ny+1),intent(in) :: dhdx,dhdy,dudx,dudy,dvdx,dvdy,sinh2kh 01083 01084 real*8,dimension(:,:,:),allocatable,save :: cgx,cgy,cx,cy,ctheta 01085 real*8,dimension(:,:),allocatable,save :: uwci,vwci 01086 01087 integer :: nthetalocal,nthetamax 01088 integer :: itheta,j,i 01089 real*8 :: cs,sn 01090 01091 if(.not. allocated(cgx)) then 01092 if (par%single_dir==1) then 01093 nthetamax = max(s%ntheta,s%ntheta_s) 01094 else 01095 nthetamax = s%ntheta 01096 endif 01097 allocate(cgx(s%nx+1,s%ny+1,nthetamax)) 01098 allocate(cgy(s%nx+1,s%ny+1,nthetamax)) 01099 allocate(cx(s%nx+1,s%ny+1,nthetamax)) 01100 allocate(cy(s%nx+1,s%ny+1,nthetamax)) 01101 allocate(ctheta(s%nx+1,s%ny+1,nthetamax)) 01102 if (par%wci==1) then 01103 allocate(uwci(s%nx+1,s%ny+1)) 01104 allocate(vwci(s%nx+1,s%ny+1)) 01105 endif 01106 endif 01107 01108 ! set local variables according to flavour 01109 select case (flavour) 01110 case (0) 01111 nthetalocal = s%ntheta 01112 if (par%wci==1) then 01113 uwci = s%u 01114 vwci = s%v 01115 endif 01116 case (1) 01117 nthetalocal = s%ntheta 01118 if (par%wci==1) then 01119 uwci = s%uwcins 01120 vwci = s%vwcins 01121 endif 01122 case (2) 01123 nthetalocal = s%ntheta_s 01124 if (par%wci==1) then 01125 uwci = s%uws 01126 vwci = s%vws 01127 endif 01128 end select 01129 ! 01130 ! split wave velocities in wave grid directions theta 01131 do itheta=1,nthetalocal 01132 do j=1,s%ny+1 01133 do i=1,s%nx+1 01134 if (s%wete(i,j) == 1) then 01135 ! wave grid directions theta with respect to spatial grid s,n 01136 if (flavour == 2) then 01137 cs=s%costh_s(i,j,itheta) 01138 sn=s%sinth_s(i,j,itheta) 01139 else 01140 cs=s%costh(i,j,itheta) 01141 sn=s%sinth(i,j,itheta) 01142 endif 01143 01144 ! split wave velocities over theta bins 01145 ! note: cg is already correct for each flavour, as long as wave_dispersion is called before this subroutine 01146 if (par%wci==1) then 01147 cgx(i,j,itheta)= s%cg(i,j)*cs+uwci(i,j) 01148 cgy(i,j,itheta)= s%cg(i,j)*sn+vwci(i,j) 01149 cx(i,j,itheta) = s%c(i,j)*cs+uwci(i,j) 01150 cy(i,j,itheta) = s%c(i,j)*sn+vwci(i,j) 01151 else 01152 cgx(i,j,itheta)= s%cg(i,j)*cs 01153 cgy(i,j,itheta)= s%cg(i,j)*sn 01154 cx(i,j,itheta) = s%c(i,j)*cs 01155 cy(i,j,itheta) = s%c(i,j)*sn 01156 endif 01157 01158 ! compute refraction velocity 01159 ! note: sigm is already correct for each flavour, as long as wave_dispersion is called before this subroutine 01160 if (par%wci==1) then 01161 ctheta(i,j,itheta)= s%sigm(i,j)/sinh2kh(i,j)*(dhdx(i,j)*sn-dhdy(i,j)*cs) + & 01162 (cs*(sn*dudx(i,j)-cs*dudy(i,j)) + sn*(sn*dvdx(i,j)-cs*dvdy(i,j))) 01163 else 01164 ctheta(i,j,itheta)= s%sigm(i,j)/sinh2kh(i,j)*(dhdx(i,j)*sn-dhdy(i,j)*cs) 01165 endif 01166 else 01167 cgx(i,j,itheta) = 0.d0 01168 cgy(i,j,itheta) = 0.d0 01169 cx(i,j,itheta) = 0.d0 01170 cy(i,j,itheta) = 0.d0 01171 ctheta(i,j,itheta) = 0.d0 01172 endif 01173 enddo 01174 enddo 01175 enddo 01176 01177 ! Dano Limit unrealistic refraction speed to 1/2 pi per wave period 01178 ctheta=sign(1.d0,ctheta)*min(abs(ctheta),.5*par%px/par%Trep) 01179 01180 ! return to the right parts in s 01181 if (flavour == 2) then 01182 s%cgx_s(:,:,1:s%ntheta_s) = cgx(:,:,1:s%ntheta_s) 01183 s%cgy_s(:,:,1:s%ntheta_s) = cgy(:,:,1:s%ntheta_s) 01184 s%ctheta_s(:,:,1:s%ntheta_s) = ctheta(:,:,1:s%ntheta_s) 01185 ! no cx and cy needed for roller energy 01186 else 01187 s%cgx(:,:,1:s%ntheta) = cgx(:,:,1:s%ntheta) 01188 s%cgy(:,:,1:s%ntheta) = cgy(:,:,1:s%ntheta) 01189 s%ctheta(:,:,1:s%ntheta) = ctheta(:,:,1:s%ntheta) 01190 s%cx(:,:,1:s%ntheta) = cx(:,:,1:s%ntheta) 01191 s%cy(:,:,1:s%ntheta) = cy(:,:,1:s%ntheta) 01192 endif 01193 01194 end subroutine compute_wave_direction_velocities 01195 01196 subroutine compute_wave_forces(s) 01197 use spaceparams 01198 01199 implicit none 01200 01201 type(spacepars), target :: s 01202 01203 integer :: i,j 01204 01205 ! 01206 ! wave radiation stress 01207 s%Sxx=(s%n*sum((1.d0+s%costh**2)*s%ee,3)-.5d0*sum(s%ee,3))*s%dtheta 01208 s%Syy=(s%n*sum((1.d0+s%sinth**2)*s%ee,3)-.5d0*sum(s%ee,3))*s%dtheta 01209 s%Sxy=s%n*sum(s%sinth*s%costh*s%ee,3)*s%dtheta 01210 ! 01211 ! 01212 ! add roller contribution 01213 s%Sxx = s%Sxx + sum((s%costh**2)*s%rr,3)*s%dtheta 01214 s%Syy = s%Syy + sum((s%sinth**2)*s%rr,3)*s%dtheta 01215 s%Sxy = s%Sxy + sum(s%sinth*s%costh*s%rr,3)*s%dtheta 01216 ! 01217 ! 01218 ! Wave forces 01219 if (s%ny>0) then 01220 do j=2,s%ny 01221 do i=1,s%nx 01222 s%Fx(i,j)=-(s%Sxx(i+1,j)-s%Sxx(i,j))/s%dsu(i,j) & 01223 -(s%Sxy(i,j+1)+s%Sxy(i+1,j+1)-s%Sxy(i,j-1)-s%Sxy(i+1,j-1))/ & 01224 (s%dnv(i,j-1)+s%dnv(i,j)+s%dnv(i+1,j-1)+s%dnv(i+1,j)) 01225 enddo 01226 enddo 01227 01228 do j=1,s%ny 01229 do i=2,s%nx 01230 s%Fy(i,j)=-(s%Syy(i,j+1)-s%Syy(i,j))/s%dnv(i,j) & 01231 -(s%Sxy(i+1,j)+s%Sxy(i+1,j+1)-s%Sxy(i-1,j)-s%Sxy(i-1,j+1))/ & 01232 (s%dsu(i-1,j)+s%dsu(i,j)+s%dsu(i-1,j+1)+s%dsu(i,j+1)) 01233 enddo 01234 enddo 01235 else 01236 j=1 01237 do i=1,s%nx 01238 s%Fx(i,j)=-(s%Sxx(i+1,j)-s%Sxx(i,j))/s%dsu(i,j) 01239 enddo 01240 do i=2,s%nx 01241 s%Fy(i,j)=-(s%Sxy(i+1,j)-s%Sxy(i-1,j))/ (s%dsu(i-1,j)+s%dsu(i,j)) 01242 enddo 01243 endif 01244 ! wwvv in the previous, Fx and Fy are computed. The missing elements 01245 ! elements are Fx(:,1), Fx(nx+1,:), Fx(:,ny+1) 01246 ! Fy(1,:), Fy(nx+1,:), Fy(:,ny+1) 01247 01248 ! wwvv so, Fx(:ny+1) and Fy(:ny+1) are left zero and Fx(nx+1,:) and Fy(nx+1,:) 01249 ! are made zero. In the parallel case, Fx(:,1) and Fy(1,:) don't get a 01250 ! value if the submatrices are not on suitable border. 01251 ! I guess that it is necessary to communicate with neighbours the values of these elements 01252 01253 ! wwvv todo the following has consequences for // version 01254 if(xmpi_istop) then 01255 s%Fy(1,:)=s%Fy(2,:) 01256 endif 01257 if(xmpi_isbot) then 01258 s%Fx(s%nx+1,:) = 0.0d0 01259 s%Fy(s%nx+1,:) = 0.0d0 01260 endif 01261 if(xmpi_isleft .and. s%ny>0) then 01262 s%Fx(:,1)=s%Fx(:,2) 01263 ! Fy(:,1)=Fy(:,2) 01264 ! ! where (Fy(:,1)>0.d0) Fy(:,1)=0.d0; !Jaap + Bas: Don't do this before calling Dano :-) 01265 endif 01266 if (xmpi_isright .and. s%ny>0) then 01267 ! Fy(:,ny+1)=Fy(:,ny) ! only a dummy point in non mpi 01268 s%Fx(:,s%ny+1)=s%Fx(:,s%ny) ! only a dummy point in non mpi 01269 ! ! where (Fy(:,ny+1)<0.d0) Fy(:,ny+1)=0.d0; !Jaap + Bas: Don't do this before calling Dano :-) 01270 endif 01271 01272 end subroutine compute_wave_forces 01273 01274 subroutine compute_stokes_drift(s,par) 01275 use params 01276 use spaceparams 01277 01278 implicit none 01279 01280 type(spacepars), target :: s 01281 type(parameters) :: par 01282 01283 real*8 , dimension(:,:) ,allocatable,save :: ustw 01284 01285 if (.not. allocated(ustw)) then 01286 allocate(ustw (s%nx+1,s%ny+1)) 01287 endif 01288 01289 where(s%wete==1) 01290 ! wave induced mass flux 01291 ustw=s%E/s%c/par%rho/s%hstokes ! Jaap (Robert: used in wave instationary) 01292 ! ustw=s%E*s%k/s%sigm/par%rho/max(s%hh,.001d0) ! Robert: old version used in wave stationary. Superceded? 01293 s%uwf = ustw*cos(s%thetamean-s%alfaz) 01294 s%vwf = ustw*sin(s%thetamean-s%alfaz) 01295 ! 01296 ! roller contribution 01297 s%ustr=2*s%R/s%c/par%rho/s%hstokes ! Jaap (Robert: used in wave instationary) 01298 !s%ustr = 0.d0 01299 !s%ustr=2*s%R*s%k/s%sigm/par%rho/max(s%hh,.001d0) ! Robert: old version used in wave stationary. Superceded? 01300 elsewhere 01301 ustw = 0.d0 01302 s%uwf = 0.d0 01303 s%vwf = 0.d0 01304 s%ustr = 0.d0 01305 endwhere 01306 ! 01307 ! introduce breaker delay 01308 if (par%breakerdelay == 1) then 01309 call breakerdelay(par,s) 01310 s%ust = ustw+s%usd 01311 else 01312 s%ust = ustw+s%ustr 01313 endif 01314 ! 01315 ! boundary conditions 01316 if(xmpi_istop) then 01317 s%ust(1,:) = s%ust(2,:) 01318 endif 01319 if(xmpi_isleft.and.s%ny>0) then 01320 s%ust(:,1) = s%ust(:,2) 01321 endif 01322 if(xmpi_isright.and. s%ny>0) then 01323 s%ust(:,s%ny+1) = s%ust(:,s%ny) 01324 endif 01325 01326 end subroutine compute_stokes_drift 01327 01328 end module wave_functions_module
1.7.6.1
