|
XBeach
|
00001 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00002 ! Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University ! 00003 ! Robert McCall, Dano Roelvink, Ap van Dongeren ! 00004 ! ! 00005 ! ! 00006 ! d.roelvink@unesco-ihe.org ! 00007 ! UNESCO-IHE Institute for Water Education ! 00008 ! P.O. Box 3015 ! 00009 ! 2601 DA Delft ! 00010 ! The Netherlands ! 00011 ! ! 00012 ! This library is free software; you can redistribute it and/or ! 00013 ! modify it under the terms of the GNU Lesser General Public ! 00014 ! License as published by the Free Software Foundation; either ! 00015 ! version 2.1 of the License, or (at your option) any later version. ! 00016 ! ! 00017 ! This library is distributed in the hope that it will be useful, ! 00018 ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! 00019 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! 00020 ! Lesser General Public License for more details. ! 00021 ! ! 00022 ! You should have received a copy of the GNU Lesser General Public ! 00023 ! License along with this library; if not, write to the Free Software ! 00024 ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 ! 00025 ! USA ! 00026 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00027 module compute_tide_module 00028 00029 00030 implicit none 00031 save 00032 private 00033 00034 real*8,dimension(4),private :: ndistcorners, sdistcorners 00035 00036 public tide_init ! routine called on global xmaster for full grid to initialise tide 00037 public tide_boundary_timestep ! routine called by all processes to update tide boundary during the run 00038 private timeinterp_tide ! routine called by any processes to interpolate tide time series to current timestep 00039 private boundaryinterp_tide ! routine called by any processes to interpolate current tide level at model corners to current 00040 ! level along the model boundary 00041 private boundaryinterp_tide_complex ! routine called by any processes to compute tide levels on complex boundaries 00042 ! (interpolation, splitting, MPI, etc.) 00043 private fill_tide_grid ! routine called only by xmaster process to fill an initial field with zs0 from zs0 time series 00044 ! and values generated at the boundary. Copies method of boundary_tide_complex to initialise in the 00045 ! centre of the domain 00046 00047 contains 00048 00049 subroutine tide_init(s,par) 00050 ! This routine is called by xmaster on full grid 00051 use params, only: parameters 00052 use spaceparamsdef, only: spacepars 00053 use filefunctions, only: create_new_fid 00054 use logging_module, only: report_file_read_error 00055 00056 implicit none 00057 save 00058 00059 type(parameters),intent(in) :: par 00060 type(spacepars),target :: s 00061 integer :: i,j,fid,ier 00062 logical :: exists 00063 00064 ! Set global domain corner distances 00065 ndistcorners(1) = s%ndist(1,1) 00066 ndistcorners(4) = s%ndist(s%nx+1,1) 00067 ndistcorners(3) = s%ndist(s%nx+1,s%ny+1) 00068 ndistcorners(2) = s%ndist(1,s%ny+1) 00069 sdistcorners(1) = s%sdist(1,1) 00070 sdistcorners(4) = s%sdist(s%nx+1,1) 00071 sdistcorners(3) = s%sdist(s%nx+1,s%ny+1) 00072 sdistcorners(2) = s%sdist(1,s%ny+1) 00073 00074 ! if there is a file with initial zs0, then use that 00075 inquire(file=par%zsinitfile,exist=exists) 00076 if (exists) then 00077 fid = create_new_fid() 00078 open(fid,file=par%zsinitfile) 00079 do j=1,s%ny+1 00080 read(fid,*,iostat=ier)(s%zs0(i,j),i=1,s%nx+1) 00081 if (ier .ne. 0) then 00082 call report_file_read_error(par%zsinitfile) 00083 endif 00084 enddo 00085 close(fid) 00086 else ! no initial zs0 file exists 00087 ! Call tidal water level at four corners of the model at current time (returns zs01-zs04) 00088 call timeinterp_tide(s,par,0.d0) 00089 00090 ! call tide levels at the boundaries 00091 call boundaryinterp_tide(s,par,globalonly=.true.,correctforwetz=.false.) ! in this case we want all to be carried out on 00092 ! global grid, not local grids, and not correct 00093 ! for wetz, as this is not yet initialised 00094 00095 ! call zs0 water levels at intermediate points 00096 call fill_tide_grid(s) 00097 endif 00098 00099 end subroutine tide_init 00100 00101 subroutine tide_boundary_timestep(s,par) 00102 use params, only: parameters 00103 use spaceparamsdef, only: spacepars 00104 #ifdef USEMPI 00105 use xmpi_module 00106 #endif 00107 use paramsconst 00108 implicit none 00109 save 00110 00111 type(parameters),intent(in) :: par 00112 type(spacepars),target :: s 00113 logical,save :: initialised = .false. 00114 integer :: i 00115 00116 ! First we need common data from xmaster regarding corner points (initialisation done only on xmaster) 00117 #ifdef USEMPI 00118 if (.not. initialised) then 00119 do i=1,4 00120 call xmpi_bcast(ndistcorners(i)) 00121 call xmpi_bcast(sdistcorners(i)) 00122 enddo 00123 initialised = .true. 00124 endif 00125 #endif 00126 00127 ! Call tidal water level at four corners of the model at current time (returns zs01-zs04) 00128 call timeinterp_tide(s,par,par%t) 00129 00130 ! call tide levels at the boundaries 00131 call boundaryinterp_tide(s,par) ! in this case we want all to be carried out on all local grids 00132 00133 if(par%tidetype==TIDETYPE_INSTANT) then 00134 call fill_tide_grid(s) 00135 endif 00136 00137 00138 end subroutine tide_boundary_timestep 00139 00140 subroutine timeinterp_tide(s,par,timenow) 00141 ! This routine may be called by all processes and interpolates tide time series to current time step at all four corners of 00142 ! the global model grid 00143 use params, only: parameters 00144 use spaceparamsdef, only: spacepars 00145 use interp, only: linear_interp 00146 use paramsconst 00147 00148 implicit none 00149 save 00150 00151 type(parameters),intent(in) :: par 00152 type(spacepars),target :: s 00153 real*8, intent(in) :: timenow 00154 integer :: indt 00155 00156 00157 select case (par%tideloc) 00158 case (0) 00159 s%zs01 = par%zs0 00160 s%zs02 = par%zs0 00161 s%zs03 = par%zs0 00162 s%zs04 = par%zs0 00163 case (1) 00164 ! Need to interpolate to the correct moment in time. First point in tidal 00165 ! record not necessarily == 0.0 00166 call LINEAR_INTERP(s%tideinpt,s%tideinpz(:,1),s%tidelen,timenow, s%zs01, indt) 00167 s%zs02 = s%zs01 00168 s%zs03 = s%zs01 00169 s%zs04 = s%zs01 00170 case (2) 00171 select case (par%paulrevere) 00172 case(PAULREVERE_LAND) 00173 ! One time series on bay and one on sea 00174 call LINEAR_INTERP(s%tideinpt,s%tideinpz(:,1),s%tidelen,timenow, s%zs01, indt) 00175 s%zs02 = s%zs01 00176 call LINEAR_INTERP(s%tideinpt,s%tideinpz(:,2),s%tidelen,timenow, s%zs03, indt) 00177 s%zs04 = s%zs03 00178 case(PAULREVERE_SEA) 00179 ! One time series for each sea corner (copied also to bay corners) 00180 call LINEAR_INTERP(s%tideinpt,s%tideinpz(:,1),s%tidelen,timenow, s%zs01, indt) 00181 s%zs04 = s%zs01 00182 call LINEAR_INTERP(s%tideinpt,s%tideinpz(:,2),s%tidelen,timenow, s%zs02, indt) 00183 s%zs03 = s%zs02 00184 end select 00185 case (4) 00186 ! One time series per corner 00187 call LINEAR_INTERP(s%tideinpt,s%tideinpz(:,1),s%tidelen,timenow, s%zs01, indt) 00188 call LINEAR_INTERP(s%tideinpt,s%tideinpz(:,2),s%tidelen,timenow, s%zs02, indt) 00189 call LINEAR_INTERP(s%tideinpt,s%tideinpz(:,3),s%tidelen,timenow, s%zs03, indt) 00190 call LINEAR_INTERP(s%tideinpt,s%tideinpz(:,4),s%tidelen,timenow, s%zs04, indt) 00191 end select 00192 00193 end subroutine timeinterp_tide 00194 00195 subroutine boundaryinterp_tide(s,par,globalonly,correctforwetz) 00196 ! This routine interpolates in space the current tide at the four corners of the model to the (local MPI subdomain) grid 00197 ! boundary points 00198 use params, only: parameters 00199 use spaceparams 00200 use paramsconst 00201 00202 implicit none 00203 save 00204 00205 type(parameters),intent(in) :: par 00206 type(spacepars),target :: s 00207 logical,intent(in),optional :: globalonly,correctforwetz 00208 logical :: globalonly_local,correctforwetz_local 00209 integer :: i,j 00210 real*8 :: tidegradient 00211 00212 if(present(globalonly)) then 00213 globalonly_local = globalonly 00214 else 00215 globalonly_local = .false. 00216 endif 00217 00218 if(present(correctforwetz)) then 00219 correctforwetz_local = correctforwetz 00220 else 00221 correctforwetz_local = .true. 00222 endif 00223 00224 select case (par%tideloc) 00225 case (0,1) 00226 ! No spatially-varying boundary conditions 00227 if (xmpi_istop .or. globalonly_local) s%zs0(1:2,:) = s%zs01 ! carry out also if running in global grid 00228 if (xmpi_isbot .or. globalonly_local) s%zs0(s%nx:s%nx+1,:) = s%zs01 00229 if (xmpi_isleft .or. globalonly_local) s%zs0(:,1) = s%zs01 00230 if (par%ny>0) then ! stop unnecessary additional work 00231 if (xmpi_isright .or. globalonly_local) s%zs0(:,s%ny+1) = s%zs01 00232 endif 00233 case (2) 00234 select case (par%paulrevere) 00235 case(PAULREVERE_LAND) 00236 ! One time series on bay and one on sea, easy for onshore and offshore boundaries 00237 if (xmpi_istop .or. globalonly_local) s%zs0(1:2,:) = s%zs01 00238 if (xmpi_isbot .or. globalonly_local) s%zs0(s%nx:s%nx+1,:) = s%zs03 00239 call boundaryinterp_tide_complex(s,globalonly_local,3) ! left 00240 if (par%ny>0) then 00241 call boundaryinterp_tide_complex(s,globalonly_local,4) ! right 00242 endif 00243 case(PAULREVERE_SEA) 00244 ! One time series for each of the two sea points (copied to land points) 00245 if (xmpi_isleft .or. globalonly_local) s%zs0(:,1) = s%zs01 00246 if (par%ny>0) then 00247 if (xmpi_isright .or. globalonly_local) s%zs0(:,s%ny+1) = s%zs02 00248 endif 00249 ! add horizontal tide gradient at boundaries 00250 if(s%ny==0) then 00251 tidegradient = (s%zs02-s%zs01)/s%dnz(1,1) 00252 else 00253 tidegradient = (s%zs02-s%zs01)/(ndistcorners(2)-ndistcorners(1)) 00254 endif 00255 if (xmpi_isleft .or. globalonly_local) s%dzs0dn(:,1) = tidegradient*s%wetz(:,1) 00256 if (par%ny>0) then 00257 if (xmpi_isright .or. globalonly_local) s%dzs0dn(:,s%ny+1) = tidegradient*s%wetz(:,s%ny+1) 00258 endif 00259 ! 00260 ! 00261 call boundaryinterp_tide_complex(s,globalonly_local,1) ! front 00262 call boundaryinterp_tide_complex(s,globalonly_local,2) ! back 00263 if (xmpi_istop .or. globalonly_local) s%zs0(2,:) = s%zs0(1,:) 00264 if (xmpi_isbot .or. globalonly_local) s%zs0(s%nx,:) = s%zs0(s%nx+1,:) 00265 end select 00266 case(4) 00267 ! All four corners can differ 00268 call boundaryinterp_tide_complex(s,globalonly_local,1) ! front 00269 call boundaryinterp_tide_complex(s,globalonly_local,2) ! back 00270 call boundaryinterp_tide_complex(s,globalonly_local,3) ! left 00271 if (par%ny>0) then 00272 call boundaryinterp_tide_complex(s,globalonly_local,4) ! right 00273 endif 00274 if (xmpi_istop .or. globalonly_local) s%zs0(2,:) = s%zs0(1,:) 00275 if (xmpi_isbot .or. globalonly_local) s%zs0(s%nx,:) = s%zs0(s%nx+1,:) 00276 end select 00277 ! 00278 ! now set zs0 to zb level in all cells where the boundary neighbour is dry 00279 ! wetz not initialised until after tide initialised, so only apply this step after the first time step 00280 if (correctforwetz_local) then 00281 ! Do whole domain 00282 where(s%wetz==0) 00283 s%zs0 = s%zb 00284 endwhere 00285 ! Fix 2D cases 00286 if (par%ny>0) then 00287 do i = 1,s%nx+1 00288 if (s%wetz(i,2)==0) then 00289 s%zs0(i,1) = min(s%zb(i,1),s%zb(i,2)) 00290 endif 00291 if (s%wetz(i,s%ny)==0) then 00292 s%zs0(i,s%ny+1) = min(s%zb(i,s%ny+1),s%zb(i,s%ny)) 00293 endif 00294 enddo 00295 endif 00296 ! Fix front and back 00297 do j=1,s%ny+1 00298 ! front 00299 if(s%wetz(2,j)==0) then 00300 s%zs0(1,j) = min(s%zb(1,j),s%zb(2,j)) 00301 endif 00302 ! back 00303 if(s%wetz(s%nx,j)==0) then 00304 s%zs0(s%nx+1,j) = min(s%zb(s%nx,j),s%zb(s%nx+1,j)) 00305 endif 00306 enddo 00307 endif ! correctforwetz_local 00308 end subroutine boundaryinterp_tide 00309 00310 subroutine boundaryinterp_tide_complex(s,globalonly_local,boundaryID) 00311 use spaceparams 00312 use xmpi_module 00313 00314 implicit none 00315 save 00316 00317 type(spacepars),target :: s 00318 logical,intent(in) :: globalonly_local 00319 integer,intent(in) :: boundaryID 00320 logical :: mpidiscretisation 00321 integer :: i,j 00322 real*8 :: maxzbg,maxzbl 00323 real*8 :: distmaxzbg,distmaxzbl 00324 integer :: relevantxmpisize 00325 integer :: imin,imax,jmin,jmax 00326 logical :: relevantxmpidomain 00327 real*8 :: tide1,tide2,disttide1,disttide2 00328 real*8 :: tidegradient1,tidegradient2 00329 real*8 :: fac,ddist 00330 00331 00332 ! Based on boundaryID (1 = front, 2 = back, 3 = left, 4 = right), we need to set some common switches 00333 #ifdef USEMPI 00334 select case (boundaryID) 00335 case (1,2) 00336 relevantxmpisize = xmpi_n 00337 case (3,4) 00338 relevantxmpisize = xmpi_m 00339 end select 00340 #else 00341 relevantxmpisize = 0 00342 #endif 00343 select case (boundaryID) 00344 case (1) ! front 00345 imin = 1 00346 imax = 1 00347 jmin = 1 00348 jmax = s%ny+1 00349 relevantxmpidomain = xmpi_istop 00350 tide1 = s%zs01 00351 tide2 = s%zs02 00352 disttide1 = ndistcorners(1) 00353 disttide2 = ndistcorners(2) 00354 ! we don't use tidal gradient on front and back boundaries 00355 tidegradient1 = 0.d0 ! (s%zs04-s%zs01)/(sdistcorners(4)-sdistcorners(1)) 00356 tidegradient2 = 0.d0 !(s%zs03-s%zs02)/(sdistcorners(3)-sdistcorners(2)) 00357 case (2) ! back 00358 imin = s%nx+1 00359 imax = s%nx+1 00360 jmin = 1 00361 jmax = s%ny+1 00362 relevantxmpidomain = xmpi_isbot 00363 tide1 = s%zs04 00364 tide2 = s%zs03 00365 disttide1 = ndistcorners(4) 00366 disttide2 = ndistcorners(3) 00367 tidegradient1 = 0.d0 ! (s%zs04-s%zs01)/(sdistcorners(4)-sdistcorners(1)) 00368 tidegradient2 = 0.d0 ! (s%zs03-s%zs02)/(sdistcorners(3)-sdistcorners(2)) 00369 case (3) ! left (flow = right mpi) 00370 imin = 1 00371 imax = s%nx+1 00372 jmin = s%ny+1 00373 jmax = s%ny+1 00374 relevantxmpidomain = xmpi_isright 00375 tide1 = s%zs02 00376 tide2 = s%zs03 00377 disttide1 = sdistcorners(2) 00378 disttide2 = sdistcorners(3) 00379 if(s%ny==0) then 00380 tidegradient1 = (s%zs02-s%zs01)/s%dnz(1,1) 00381 tidegradient2 = (s%zs03-s%zs04)/s%dnz(s%nx+1,1) 00382 else 00383 tidegradient1 = (s%zs02-s%zs01)/(ndistcorners(2)-ndistcorners(1)) 00384 tidegradient2 = (s%zs03-s%zs04)/(ndistcorners(3)-ndistcorners(4)) 00385 endif 00386 case (4) ! right (flow = left mpi) 00387 imin = 1 00388 imax = s%nx+1 00389 jmin = 1 00390 jmax = 1 00391 relevantxmpidomain = xmpi_isleft 00392 tide1 = s%zs01 00393 tide2 = s%zs04 00394 disttide1 = sdistcorners(1) 00395 disttide2 = sdistcorners(4) 00396 if(s%ny==0) then 00397 tidegradient1 = (s%zs02-s%zs01)/s%dnz(1,1) 00398 tidegradient2 = (s%zs03-s%zs04)/s%dnz(s%nx+1,1) 00399 else 00400 tidegradient1 = (s%zs02-s%zs01)/(ndistcorners(2)-ndistcorners(1)) 00401 tidegradient2 = (s%zs03-s%zs04)/(ndistcorners(3)-ndistcorners(4)) 00402 endif 00403 end select 00404 ! 00405 ! We decide if there is MPI subdomain discretization along the boundary. This is not the case if run 00406 ! without MPI, or if MPI is passing global field (optional input "globalonly"), or if the MPI division is 00407 ! such that there is discretisation along the boundary (i.e., xmpi_m = 1, or xmpi_n = 1) 00408 ! 00409 #ifdef USEMPI 00410 if (globalonly_local .or. relevantxmpisize==1) then ! only runnnig on global s, or domain has no MPI discretization 00411 ! in cross-shore direction 00412 mpidiscretisation = .false. 00413 else 00414 mpidiscretisation = .true. 00415 endif 00416 #else 00417 mpidiscretisation = .false. 00418 #endif 00419 ! 00420 ! find maximum bed level on boundary of each (or single) domain 00421 maxzbl = maxval(s%zb(imin:imax,jmin:jmax)) 00422 ! find sdist at which maximum bed level occurs 00423 select case (boundaryID) 00424 case(1,2) 00425 distmaxzbl = s%ndist(imin,maxval(maxloc(s%zb(imin,:)))) 00426 case(3,4) 00427 distmaxzbl = s%sdist(maxval(maxloc(s%zb(:,jmin))),jmin) 00428 end select 00429 ! 00430 ! if not running in single domain, then we need to communicate maximum bed levels and location thereof 00431 if (mpidiscretisation) then 00432 #ifdef USEMPI 00433 ! remove data from domains not on the boundary 00434 if(.not. relevantxmpidomain) maxzbl = -huge(0.d0) 00435 ! store results in one array to be communicated 00436 maxzbg = maxzbl 00437 ! find the maximum of local maxima (in precompile statement to allow to compile) 00438 call xmpi_allreduce(maxzbg,MPI_MAX) 00439 ! now find the corresponding distance value if the maximum value is in your domain 00440 if(abs(maxzbg-maxzbl)<1.d-8) then 00441 distmaxzbg = distmaxzbl 00442 else 00443 distmaxzbg = - huge(0.d0) 00444 endif 00445 ! and communicate location of global maximum 00446 call xmpi_allreduce(distmaxzbg,MPI_MAX) 00447 #endif 00448 else 00449 maxzbg = maxzbl 00450 distmaxzbg = distmaxzbl 00451 endif 00452 00453 ! Now 4 cases for imposing water levels on boundary: 00454 ! 1) Water level offshore (OWL) and bay (BWL) both lower than maximum bed level (zbmax), then OWB from offshore 00455 ! to smaxzbg, and BWL from smaxzbg to bay 00456 ! 2) OWL and BWL greater than maxzb, then linear interpolation along boundary 00457 ! 3) OWL greater than maxzb, but BWL not, then OWL from offshore to smaxzbg and interpolation from smaxzbg to 00458 ! bay 00459 ! 4) BWL greater than maxzb, but OWL not, then BWL from bay to smaxzbg and interpolation from smaxzbg to 00460 ! offshore 00461 ! 00462 00463 ! Only carry this out on the relevant mpi subdomains or if running everything on the global grid (xmaster is not 00464 ! necessarily also top, bottom, left and right) 00465 if (relevantxmpidomain .or. globalonly_local) then 00466 ! (1) 00467 if (tide1<=maxzbg .and. tide2<=maxzbg) then 00468 select case (boundaryID) 00469 case(1,2) ! top, bottom 00470 where(s%ndist(imin,:)<=distmaxzbg) 00471 s%zs0(imin,:) = tide1 00472 elsewhere 00473 s%zs0(imin,:) = tide2 00474 endwhere 00475 case(3,4) ! left, right 00476 where(s%sdist(:,jmin)<=distmaxzbg) 00477 s%zs0(:,jmin) = tide1 00478 s%dzs0dn(:,jmin) = tidegradient1*s%wetz(:,jmin) 00479 elsewhere 00480 s%zs0(:,jmin) = tide2 00481 s%dzs0dn(:,jmin) = tidegradient2*s%wetz(:,jmin) 00482 endwhere 00483 end select 00484 ! (2) 00485 elseif (tide1>maxzbg .and. tide2>maxzbg) then 00486 select case (boundaryID) 00487 case(1,2) ! top, bottom 00488 ! don't want to call this using linear_interp function, because that requires having the whole array of ndist 00489 ddist = 1.d0/(disttide2-disttide1) 00490 do j = jmin,jmax 00491 fac = (s%ndist(imin,j)-disttide1)*ddist 00492 s%zs0(imin,j) = (1.d0-fac)*tide1+fac*tide2 00493 enddo 00494 case(3,4) ! left, right 00495 ddist = 1.d0/(disttide2-disttide1) 00496 do i = imin,imax 00497 fac = (s%sdist(i,jmin)-disttide1)*ddist 00498 s%zs0(i,jmin) = (1.d0-fac)*tide1+fac*tide2 00499 s%dzs0dn(i,jmin) = ((1.d0-fac)*tidegradient1+fac*tidegradient2)*s%wetz(i,jmin) 00500 enddo 00501 end select 00502 ! (3) 00503 elseif (tide1>maxzbg .and. tide2<=maxzbg) then 00504 select case (boundaryID) 00505 case(1,2) ! top, bottom 00506 ! don't want to call this using linear_interp function, because that requires having the whole array of ndist 00507 ddist = 1.d0/(disttide2-distmaxzbg) 00508 do j = jmin,jmax 00509 if (s%ndist(imin,j)<=distmaxzbg) then 00510 s%zs0(imin,j) = tide1 00511 else 00512 fac = (s%ndist(imin,j)-distmaxzbg)*ddist 00513 s%zs0(imin,j) = (1.d0-fac)*tide1+fac*tide2 00514 endif 00515 enddo 00516 case (3,4) ! left,right 00517 ddist = 1.d0/(disttide2-distmaxzbg) 00518 do i = imin,imax 00519 if (s%sdist(i,jmin)<=distmaxzbg) then 00520 s%zs0(i,jmin) = tide1 00521 s%dzs0dn(i,jmin) = tidegradient1*s%wetz(i,jmin) 00522 else 00523 fac = (s%sdist(i,jmin)-distmaxzbg)*ddist 00524 s%zs0(i,jmin) = (1.d0-fac)*tide1+fac*tide2 00525 s%dzs0dn(i,jmin) = ((1.d0-fac)*tidegradient1+fac*tidegradient2)*s%wetz(i,jmin) 00526 endif 00527 enddo 00528 end select 00529 ! (4) 00530 elseif(tide1<=maxzbg .and. tide2>maxzbg) then 00531 select case (boundaryID) 00532 case(1,2) ! top, bottom 00533 ! don't want to call this using linear_interp function, because that requires having the whole array of ndist 00534 ddist = 1.d0/(distmaxzbg-disttide1) 00535 do j = jmin,jmax 00536 if (s%ndist(imin,j)<=distmaxzbg) then 00537 fac = (s%ndist(imin,j)-disttide1)*ddist 00538 s%zs0(imin,j) = (1.d0-fac)*tide1+fac*tide2 00539 else 00540 s%zs0(imin,j) = tide2 00541 endif 00542 enddo 00543 case (3,4) ! left,right 00544 ddist = 1.d0/(distmaxzbg-disttide1) 00545 do i = imin,imax 00546 if (s%sdist(i,jmin)<=distmaxzbg) then 00547 fac = (s%sdist(i,jmin)-disttide1)*ddist 00548 s%zs0(i,jmin) = (1.d0-fac)*tide1+fac*tide2 00549 s%dzs0dn(i,jmin) = ((1.d0-fac)*tidegradient1+fac*tidegradient2)*s%wetz(i,jmin) 00550 else 00551 s%zs0(i,jmin) = tide2 00552 s%dzs0dn(i,jmin) = tidegradient1*s%wetz(i,jmin) 00553 endif 00554 enddo 00555 end select 00556 endif ! selection of tide combinations 00557 ! 00558 ! make sure that zs0 is not below zb 00559 select case (boundaryID) 00560 case(1,2) ! top, bottom 00561 s%zs0(imin,:) = max(s%zs0(imin,:),s%zb(imin,:)) 00562 case(3,4) ! left, right 00563 s%zs0(:,jmin) = max(s%zs0(:,jmin),s%zb(:,jmin)) 00564 end select 00565 endif ! relevant mpi subdomain 00566 00567 end subroutine boundaryinterp_tide_complex 00568 00569 subroutine fill_tide_grid(s) 00570 ! This subroutine fills an initial zs0 field, done by xmaster. Not called during the remainder of the simulation. 00571 use spaceparamsdef, only: spacepars 00572 00573 implicit none 00574 save 00575 00576 type(spacepars),target :: s 00577 integer :: i,j 00578 real*8 :: maxzbg 00579 real*8 :: distmaxzbg 00580 real*8 :: tide1,tide2 00581 real*8 :: fac,ddist 00582 00583 ! Run down each cross-shore grid line. Note that this will automatically be skipped if ny=0. This is fine, zs0 at ny=1 00584 ! is already set by boundary condition call 00585 do j=2,s%ny 00586 ! find maximum bed level in each grid line 00587 maxzbg = maxval(s%zb(:,j)) 00588 ! find sdist at which maximum bed level occurs 00589 distmaxzbg = s%sdist(maxloc(s%zb(:,j),dim=1),j) 00590 ! find offshore (tide1) and bay (tide2) tide level 00591 tide1 = s%zs0(1,j) 00592 tide2 = s%zs0(s%nx+1,j) 00593 00594 ! Now 4 cases for imposing water levels on boundary: 00595 ! 1) Water level offshore (OWL) and bay (BWL) both lower than maximum bed level (zbmax), then OWB from offshore 00596 ! to smaxzbg, and BWL from smaxzbg to bay 00597 ! 2) OWL and BWL greater than maxzb, then linear interpolation along boundary 00598 ! 3) OWL greater than maxzb, but BWL not, then OWL from offshore to smaxzbg and interpolation from smaxzbg to 00599 ! bay 00600 ! 4) BWL greater than maxzb, but OWL not, then BWL from bay to smaxzbg and interpolation from smaxzbg to 00601 ! offshore 00602 ! 00603 ! 00604 ! Note that we only fill in i = 3:nx-1, for some reason i=2 is currently always copied from i=1 and i=nx from i=nx+1. 00605 ! Best not disturb that. 00606 ! (1) 00607 if (tide1<=maxzbg .and. tide2<=maxzbg) then 00608 where(s%sdist(3:s%nx-1,j)<=distmaxzbg) 00609 s%zs0(3:s%nx-1,j) = tide1 00610 elsewhere 00611 s%zs0(3:s%nx-1,j) = tide2 00612 endwhere 00613 ! (2) 00614 elseif (tide1>maxzbg .and. tide2>maxzbg) then 00615 ddist = 1.d0/(s%sdist(s%nx+1,j)-s%sdist(1,j)) 00616 do i = 3,s%nx-1 00617 fac = (s%sdist(i,j)-s%sdist(1,j))*ddist 00618 s%zs0(i,j) = (1.d0-fac)*tide1+fac*tide2 00619 enddo 00620 ! (3) 00621 elseif (tide1>maxzbg .and. tide2<=maxzbg) then 00622 ddist = 1.d0/(s%sdist(s%nx+1,j)-distmaxzbg) 00623 do i = 3,s%nx-1 00624 if (s%sdist(i,j)<=distmaxzbg) then 00625 s%zs0(i,j) = tide1 00626 else 00627 fac = (s%sdist(i,j)-distmaxzbg)*ddist 00628 s%zs0(i,j) = (1.d0-fac)*tide1+fac*tide2 00629 endif 00630 enddo 00631 ! (4) 00632 elseif(tide1<=maxzbg .and. tide2>maxzbg) then 00633 ddist = 1.d0/(distmaxzbg-s%sdist(1,j)) 00634 do i = 3,s%nx-1 00635 if (s%sdist(i,j)<=distmaxzbg) then 00636 fac = (s%sdist(i,j)-s%sdist(1,j))*ddist 00637 s%zs0(i,j) = (1.d0-fac)*tide1+fac*tide2 00638 else 00639 s%zs0(i,j) = tide2 00640 endif 00641 enddo 00642 endif ! selection of tide combinations 00643 s%zs0(3:s%nx-1,j) = max(s%zs0(3:s%nx-1,j),s%zb(3:s%nx-1,j)) 00644 enddo ! j=2,ny 00645 00646 end subroutine fill_tide_grid 00647 00648 end module
1.7.6.1
