|
XBeach
|
00001 module timestep_module 00002 00003 00004 ! Time steps in XBeach 00005 ! 00006 ! Relevant time variables 00007 ! 00008 ! [Main time loop] 00009 ! tstop: stop time simulation (set in params, used in boundaryconditions groundwater output params readtide readwind timestep varoutput waveparams xbeach) 00010 ! 00011 ! [Output] 00012 ! tint: time interval output global values (set in params, used in nonhydrostatic output params timestep timestep_old, varoutput) 00013 ! 00014 ! tintg: time interval output global values (set in params, used in params timestep) 00015 ! tintm: time interval output mean global values (set in params timestep varoutput, used in params timestep varianceupdate varoutput) 00016 ! tintp: time interval output point values (set in params, used in drifters params timestep) 00017 ! 00018 ! tsglobal: file with list of output times for global values (set in timestep, used in timestep) 00019 ! tsmean: file with list of output times for meanglobal values (set in timestep, used in timestep) 00020 ! tspoints: file with list of output times for point values (set in timestep, used in timestep) 00021 ! 00022 ! tnext: (set in params timestep timestep_old varoutput, used in params timestep timestep_old varoutput) 00023 ! 00024 ! [Hydrodynamic] 00025 ! tstart: start time of simulation (or output?!) (set in params, used in nonhydrostatic output params timestep) 00026 ! CFL: maximum courant number (set in params, used in params timestep timestep_old) 00027 ! 00028 ! [Boundary condition] 00029 ! rt: the record length of the boundary condition (set in params, waveparams, used in boundaryconditions params timestep waveparams) 00030 ! dtbc: boundary condition file time step dtbc (set in params, waveparams, used in params timestep waveparams) 00031 ! taper: time to spin up wave boundary conditions (set in params, used in nonhydrostatic output params timestep timestep_old varoutput) 00032 ! 00033 ! [Morphology] 00034 ! morfac: (set in params, used in boundaryconditions morphevolution params readtide readwind timestep varoutput waveparams) 00035 ! morstart: (set in params, used in morphevolution params timestep) 00036 ! 00037 ! [Waves] 00038 ! wavint: interval between stationary wave module calls (set in params, used in params timestep xbeach) 00039 ! 00040 ! [Tides] 00041 ! tidelen: length of tidal record (set in params readtide, used in boundaryconditions params readtide timestep) 00042 ! tideloc: number of input tidal time series (set in boundaryconditions groundwater params, used in boundaryconditions groundwater initialize params readtide timestep) 00043 ! 00044 ! Relevant Time functions and subroutines 00045 ! outputtimes_update 00046 ! 00047 ! XBeach time loop 00048 ! 0 tstop[s] 00049 ! |------------------------------------------------------------------------------------------| 00050 ! 00051 ! Hydrodynamic time scale flow_timestep(input: par%dt, output:s%hh) 00052 ! |--------------------|------------------|------------------|---------------|--------|------| 00053 ! The hydrodynamic timeloop computes at varying intervals. 00054 ! The interval (input:s%hh, par%CFL, par%dt, par%xz, s%nx, s%xz, s%nx, output=par%dt) is computed 00055 ! in timestep module 00056 ! 00057 ! The morphodynamic time scale(input: par%dt, par%morstart, p%morfac) 00058 ! |--------|------------------|------------------|---------------|--------|------| 00059 ! This scale can start at a different point in time (morstart) and uses a multiplication factor 00060 ! (morfac [s/s] (>=1)). 00061 ! 00062 ! Wave time scale 00063 ! |----------------|----------------|----------------|----------------|----------------|-----| 00064 ! At different time scales using wavint. 00065 ! Wave dispersion 00066 ! | 00067 ! Only calculated if par%t == par%dt 00068 ! 00069 ! Global output time scale fixed interval (tintg) 00070 ! |---------------------------|---------------------------|---------------------------|------- 00071 ! 00072 ! Mean output time scale fixed interval (tintm) 00073 ! |------------------------|------------------------|------------------------|---------------- 00074 ! 00075 ! Point out time scale fixed inteval(tintp) 00076 ! |---------------------------------|---------------------------------|------------------------ 00077 ! 00078 ! Global output time scale varying interval (tsglobal) 00079 ! |--------------------|-----------|---------------|-|-|------------|-------------|------------- 00080 ! 00081 ! Global output time scale varying interval (tsmean) 00082 ! |-----------------|-----------|-----------|-|-|------------|-------------|-------------------- 00083 ! 00084 ! Global output time scale varying interval (tspoint) 00085 ! |-------------------------|-----------|---------------|-|-|------------|-------------|-------- 00086 ! 00087 ! Boundary condition length [s] (rt) 00088 ! |-------------|-------------|-------------|-------------|-------------|-------------|---------- 00089 00090 implicit none 00091 save 00092 00093 private 00094 public timepars, timestep_init, timestep, outputtimes_update, compute_dt 00095 ! This structure contains the time variables which are computed during the simulation. 00096 ! Parameters which can be configured by users should go in params structure. 00097 type timepars 00098 real*8,dimension(:),allocatable :: tpg ! global output times 00099 real*8,dimension(:),allocatable :: tpp ! point oputput times 00100 real*8,dimension(:),allocatable :: tpm ! time average output times 00101 real*8,dimension(:),allocatable :: tph ! hotstart file output times 00102 real*8,dimension(:),allocatable :: tpw ! wave computation times 00103 real*8 :: tnext = -123 ! next time point for output 00104 00105 integer :: itg ! global output index (1 based) 00106 integer :: itp ! point output index (1 based) 00107 integer :: itm ! time average output index (1 based) 00108 integer :: itw ! wave computation index (1 based) 00109 integer :: ith ! hotstart computation index (1 based) 00110 00111 logical :: outputg ! output global variables? 00112 logical :: outputp ! output point variables? 00113 logical :: outputm ! output time average variables? 00114 logical :: outputw ! stop for wave computation 00115 logical :: outputh ! stop for hotstart file output 00116 logical :: output ! output any variable 00117 end type timepars 00118 00119 00120 00121 contains 00122 00123 subroutine timestep_init(par, tpar) 00124 use params, only: parameters 00125 use xmpi_module 00126 use readkey_module 00127 use logging_module, only: report_file_read_error, writelog 00128 use paramsconst 00129 implicit none 00130 00131 type(parameters),intent(inout) :: par 00132 type(timepars),intent(inout) :: tpar 00133 00134 real*8 :: iir 00135 integer :: i,ii,ier 00136 !character(80) :: fname 00137 00138 #ifdef USEMPI 00139 logical, parameter :: toall = .true. 00140 logical, parameter :: tocomp = .false. 00141 #endif 00142 00143 ! initialize the time to 0 00144 par%t = 0.0d0 00145 00146 tpar%tnext = 0.0d0 00147 ! initialize timestep counters 00148 tpar%itg = 0 00149 tpar%itp = 0 00150 tpar%itm = 0 00151 tpar%itw = 0 00152 ! initialize output, default no output 00153 tpar%outputg = .FALSE. 00154 tpar%outputp = .FALSE. 00155 tpar%outputm = .FALSE. 00156 tpar%outputw = .FALSE. 00157 tpar%output = .FALSE. 00158 00159 00160 !!!!! OUTPUT TIME POINTS !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00161 ! If working in instat 0 or instat 40 we want time to coincide at wavint timesteps 00162 if (par%wavemodel==WAVEMODEL_STATIONARY .or. par%single_dir==1) then 00163 if(xmaster) then 00164 ii=ceiling(par%tstop/par%wavint) 00165 allocate(tpar%tpw(ii)) 00166 do i=1,ii 00167 tpar%tpw(i)=real(i)*par%wavint 00168 enddo 00169 endif 00170 #ifdef USEMPI 00171 call xmpi_bcast(ii,toall) 00172 if (.not.xmaster) then 00173 allocate(tpar%tpw(ii)) 00174 endif 00175 call xmpi_bcast(tpar%tpw,toall) 00176 #endif 00177 else 00178 allocate(tpar%tpw(0)) 00179 endif 00180 00181 ! If we want global output then 00182 if (par%nglobalvar/=0) then 00183 if(xmaster) then 00184 if (par%tsglobal/=' ') then 00185 open(10,file=par%tsglobal) 00186 read(10,*,iostat=ier)iir 00187 if (ier .ne. 0) then 00188 call report_file_read_error(par%tsglobal) 00189 endif 00190 ii = nint(iir) 00191 allocate(tpar%tpg(ii)) 00192 do i=1,ii 00193 read(10,*,iostat=ier)tpar%tpg(i) 00194 if (ier .ne. 0) then 00195 call report_file_read_error(par%tsglobal) 00196 endif 00197 enddo 00198 tpar%tpg=tpar%tpg/max(par%morfac,1.d0) 00199 close(10) 00200 else 00201 ii=floor((par%tstop-par%tstart)/par%tintg)+1 00202 allocate(tpar%tpg(ii)) 00203 do i=1,ii 00204 tpar%tpg(i)=par%tstart+par%tintg*real(i-1) 00205 enddo 00206 endif 00207 endif !xmaster 00208 #ifdef USEMPI 00209 call xmpi_bcast(ii,toall) 00210 if (.not.xmaster) then 00211 allocate(tpar%tpg(ii)) 00212 endif 00213 call xmpi_bcast(tpar%tpg,toall) 00214 #endif 00215 else 00216 allocate(tpar%tpg(0)) 00217 endif ! nglobalvar /=0 00218 00219 ! If we want point output then 00220 if ((par%npoints+par%nrugauge+par%ndrifter)>0) then 00221 if (xmaster) then 00222 if (par%tspoints/=' ') then 00223 open(10,file=par%tspoints) 00224 read(10,*,iostat=ier)iir 00225 if (ier .ne. 0) then 00226 call report_file_read_error(par%tspoints) 00227 endif 00228 ii = nint(iir) 00229 allocate(tpar%tpp(ii)) 00230 do i=1,ii 00231 read(10,*,iostat=ier)tpar%tpp(i) 00232 if (ier .ne. 0) then 00233 call report_file_read_error(par%tspoints) 00234 endif 00235 enddo 00236 tpar%tpp=tpar%tpp/max(par%morfac,1.d0) 00237 close(10) 00238 else 00239 ii=floor((par%tstop-par%tstart)/par%tintp)+1 00240 allocate(tpar%tpp(ii)) 00241 do i=1,ii 00242 tpar%tpp(i)=par%tstart+par%tintp*real(i-1) 00243 enddo 00244 endif 00245 endif ! xmaster 00246 #ifdef USEMPI 00247 call xmpi_bcast(ii,toall) 00248 if (.not.xmaster) then 00249 allocate(tpar%tpp(ii)) 00250 endif 00251 call xmpi_bcast(tpar%tpp,toall) 00252 #endif 00253 else 00254 allocate(tpar%tpp(0)) 00255 endif ! (npoints+nrugauge)>0 00256 00257 ! If we want time-average output then 00258 if (par%nmeanvar>0) then 00259 if(xmaster) then 00260 if (par%tsmean/=' ') then 00261 open(10,file=par%tsmean) 00262 read(10,*,iostat=ier)iir 00263 if (ier .ne. 0) then 00264 call report_file_read_error(par%tsmean) 00265 endif 00266 ii = nint(iir) 00267 allocate(tpar%tpm(ii)) 00268 do i=1,ii 00269 read(10,*,iostat=ier)tpar%tpm(i) 00270 if (ier .ne. 0) then 00271 call report_file_read_error(par%tsmean) 00272 endif 00273 enddo 00274 tpar%tpm=tpar%tpm/max(par%morfac,1.d0) 00275 close(10) 00276 else 00277 ii=floor((par%tstop-par%tstart)/par%tintm)+1 00278 if (ii<=1) then 00279 call writelog('lse','','Tintm is larger than output simulation time ') 00280 call halt_program 00281 endif 00282 allocate(tpar%tpm(ii)) 00283 do i=1,ii 00284 tpar%tpm(i)=par%tstart+par%tintm*real(i-1) 00285 enddo 00286 endif 00287 endif 00288 #ifdef USEMPI 00289 call xmpi_bcast(ii,toall) 00290 if (.not.xmaster) then 00291 allocate(tpar%tpm(ii)) 00292 endif 00293 call xmpi_bcast(tpar%tpm,toall) 00294 #endif 00295 else 00296 allocate(tpar%tpm(0)) 00297 00298 endif ! nmeanvar > 0 00299 00300 ! If we want hotstart output then 00301 if (par%writehotstart==1) then 00302 if(xmaster) then 00303 if (par%tshotstart/=' ') then 00304 open(10,file=par%tshotstart) 00305 read(10,*,iostat=ier)iir 00306 if (ier .ne. 0) then 00307 call report_file_read_error(par%tshotstart) 00308 endif 00309 ii = nint(iir) 00310 allocate(tpar%tph(ii)) 00311 do i=1,ii 00312 read(10,*,iostat=ier)tpar%tph(i) 00313 if (ier .ne. 0) then 00314 call report_file_read_error(par%tshotstart) 00315 endif 00316 enddo 00317 tpar%tph=tpar%tph/max(par%morfac,1.d0) 00318 close(10) 00319 else 00320 ii=floor(par%tstop/par%tinth) ! don't want hotstart at t=0 00321 if (ii<=0) then 00322 call writelog('lse','','Tinth is larger than output simulation time ') 00323 call halt_program 00324 endif 00325 allocate(tpar%tph(ii)) 00326 do i=1,ii 00327 tpar%tph(i)=par%tinth*real(i) ! don't start at t=0 00328 enddo 00329 endif 00330 endif 00331 #ifdef USEMPI 00332 call xmpi_bcast(ii,toall) 00333 if (.not.xmaster) then 00334 allocate(tpar%tph(ii)) 00335 endif 00336 call xmpi_bcast(tpar%tph,toall) 00337 #endif 00338 else 00339 allocate(tpar%tph(0)) 00340 00341 endif ! wrtitehotstart == 1 00342 00343 end subroutine timestep_init 00344 00345 00346 subroutine outputtimes_update(par, tpar) 00347 use params, only: parameters 00348 use xmpi_module 00349 use logging_module, only: writelog 00350 implicit none 00351 type(parameters),intent(inout) :: par 00352 type(timepars), intent(inout) :: tpar 00353 real*8 :: t1,t2,t3,t4,t5 00354 00355 !if (tpar%outputm .or. tpar%itm .eq. 0) then 00356 ! par%tintm = tpar%tpm(min(tpar%itm+2,size(tpar%tpm)))-tpar%tpm(tpar%itm+1) 00357 ! par%tintm = max(par%tintm,tiny(0.d0)) 00358 !end if 00359 00360 ! Current timestep: 00361 ! check if we need output at the current timestep 00362 ! do we have anymore output time points, is it equal to the current? 00363 if (size(tpar%tpg) .gt. 0) then 00364 tpar%outputg = abs(tpar%tpg(min(tpar%itg+1,size(tpar%tpg)))-par%t) .le. 0.0000001d0 00365 ! if global output has last output before tstop, this minimum is necessary 00366 endif 00367 if (size(tpar%tpp) .gt. 0) then 00368 tpar%outputp = abs(tpar%tpp(min(tpar%itp+1,size(tpar%tpp)))-par%t) .le. 0.0000001d0 00369 endif 00370 if (size(tpar%tpm) .gt. 0) then 00371 tpar%outputm = abs(tpar%tpm(min(tpar%itm+1,size(tpar%tpm)))-par%t) .le. 0.0000001d0 00372 endif 00373 if (size(tpar%tpw) .gt. 0) then 00374 tpar%outputw = abs(tpar%tpw(min(tpar%itw+1,size(tpar%tpw)))-par%t) .le. 0.0000001d0 00375 endif 00376 if (size(tpar%tph) .gt. 0) then 00377 tpar%outputh = abs(tpar%tph(min(tpar%ith+1,size(tpar%tph)))-par%t) .le. 0.0000001d0 00378 endif 00379 00380 ! tpar%outputg = (size(tpar%tpg) .gt. tpar%itg) .and. (abs(tpar%tpg(tpar%itg+1)-par%t) .le. 0.0000001d0) 00381 ! tpar%outputp = (size(tpar%tpp) .gt. tpar%itp) .and. (abs(tpar%tpp(tpar%itp+1)-par%t) .le. 0.0000001d0) 00382 ! tpar%outputm = (size(tpar%tpm) .gt. tpar%itm) .and. (abs(tpar%tpm(tpar%itm+1)-par%t) .le. 0.0000001d0) 00383 ! tpar%outputc = (size(tpar%tpc) .gt. tpar%itc) .and. (abs(tpar%tpc(tpar%itc+1)-par%t) .le. 0.0000001d0) 00384 ! tpar%outputw = (size(tpar%tpw) .gt. tpar%itw) .and. (abs(tpar%tpw(tpar%itw+1)-par%t) .le. 0.0000001d0) 00385 tpar%output = (tpar%outputg .or. tpar%outputp .or. tpar%outputm .or. tpar%outputw .or. tpar%outputh) 00386 00387 ! if we are updating at the current timestep, increase the counter by 1 00388 ! 00389 ! update output counter 00390 if (tpar%outputg) tpar%itg = tpar%itg + 1 00391 if (tpar%outputp) tpar%itp = tpar%itp + 1 00392 if (tpar%outputm) tpar%itm = tpar%itm + 1 00393 if (tpar%outputw) tpar%itw = tpar%itw + 1 00394 if (tpar%outputh) tpar%ith = tpar%ith + 1 00395 00396 ! Next time step: 00397 ! calculate the output for the following time 00398 ! this routine computes the next timestep at which output should be generated 00399 ! Determine next time step 00400 t1=minval(tpar%tpg,MASK=tpar%tpg .gt. par%t+0.0000001d0) 00401 t2=minval(tpar%tpp,MASK=tpar%tpp .gt. par%t+0.0000001d0) 00402 t3=minval(tpar%tpm,MASK=tpar%tpm .gt. par%t+0.0000001d0) 00403 t4=minval(tpar%tpw,MASK=tpar%tpw .gt. par%t+0.0000001d0) 00404 t5=minval(tpar%tph,MASK=tpar%tph .gt. par%t+0.0000001d0) 00405 00406 tpar%tnext=min(t1,t2,t3,t4,t5,par%tstop) 00407 00408 if (tpar%tnext .eq. huge(tpar%tpg) .and. par%t .eq. 0) then 00409 call writelog('ls','', 'no output times found, setting tnext to tstop') 00410 tpar%tnext = par%tstop 00411 end if 00412 00413 end subroutine outputtimes_update 00414 00415 00416 subroutine compute_dt(s,par, tpar, it, ilim, jlim, dtref, dt, ierr, limtypeout) 00417 use params, only: parameters 00418 use spaceparams 00419 use xmpi_module 00420 use logging_module 00421 use paramsconst 00422 00423 IMPLICIT NONE 00424 00425 type(spacepars), intent(inout) :: s 00426 type(parameters), intent(inout) :: par 00427 type(timepars), intent(inout) :: tpar 00428 integer, intent(inout) :: it 00429 integer, intent(out), optional :: ierr,limtypeout 00430 real*8, intent(in), optional :: dt 00431 real*8, intent(inout) :: dtref 00432 integer, intent(inout) :: ilim, jlim 00433 integer :: i 00434 integer :: j,j1,j2 00435 integer :: n,limtype 00436 real*8 :: mdx,mdy,tny 00437 real*8,save :: dtold 00438 real*8 :: celeru 00439 00440 if(.not. xcompute) return 00441 00442 celeru = 0.d0 00443 limtype = 0 00444 00445 ! Super fast 1D 00446 if (s%ny==0) then 00447 j1 = 1 00448 else 00449 j1 = 2 00450 endif 00451 00452 ! Calculate dt based on the Courant number. 00453 ! Check when we need next output. 00454 ! Next time step will be, min(output times, t+dt) 00455 00456 if (present(ierr)) then 00457 ierr = 0 00458 endif 00459 tny = tiny(0.d0) 00460 00461 ! Robert new time step criterion 00462 if (par%t<=0.0d0 .and. par%hotstart==0) then ! conservative estimate in case of zeros initialization 00463 par%dt = par%CFL*min(minval(s%dsu(1:s%nx+1,1:s%ny+1)), & 00464 minval(s%dnv(1:s%nx+1,1:s%ny+1))) & 00465 /sqrt(maxval(s%hh)*par%g) 00466 dtref = 0.d0 00467 if (tpar%tnext > par%dt) then 00468 n = ceiling((tpar%tnext-par%t)/par%dt) 00469 par%dt = (tpar%tnext-par%t)/n 00470 end if 00471 elseif(par%t<=0.0d0 .and. par%hotstart==1) then ! conservative estimate in case of hotstart initialization 00472 par%dt = par%CFL*min(minval(s%dsu(1:s%nx+1,1:s%ny+1)), & 00473 minval(s%dnv(1:s%nx+1,1:s%ny+1))) & 00474 /maxval(sqrt(max(s%zs-s%zb,par%eps)*par%g)+sqrt(s%uu**2+s%vv**2)) 00475 dtref = 0.d0 00476 if (tpar%tnext > par%dt) then 00477 n = ceiling((tpar%tnext-par%t)/par%dt) 00478 par%dt = (tpar%tnext-par%t)/n 00479 end if 00480 else 00481 par%dt=huge(0.0d0) ! Seed dt 00482 dtold = par%dt 00483 limtype = 0 00484 if (s%ny>2) then 00485 do j=jmin_ee,jmax_ee 00486 do i=imin_ee,imax_ee 00487 if(s%wetu(i,j)==1) then 00488 ! u-points 00489 mdx=s%dsu(i+1,j) 00490 mdy=min(s%dnv(i,j),s%dnv(i,j-1)) 00491 ! x-component 00492 par%dt=min(par%dt,mdx/max(tny,max(sqrt(par%g*s%hu(i,j))+abs(s%uu(i,j)),abs(s%ueu(i,j))))) !Jaap: include sediment advection velocities 00493 ! y-component 00494 par%dt=min(par%dt,mdy/max(tny,(sqrt(par%g*s%hu(i,j))+abs(s%vu(i,j))))) 00495 if (par%dt<dtold) then 00496 ilim = i 00497 jlim = j 00498 limtype = 1 00499 dtold = par%dt 00500 endif 00501 endif 00502 if(s%wetv(i,j)==1) then 00503 ! v-points 00504 mdx=min(s%dsu(i,j),s%dsu(i-1,j)) 00505 mdy=s%dnv(i,j) 00506 ! x-component 00507 par%dt=min(par%dt,mdx/max(tny,(sqrt(par%g*s%hv(i,j))+abs(s%uv(i,j))))) 00508 ! y-component 00509 par%dt=min(par%dt,mdy/max(tny,max(sqrt(par%g*s%hv(i,j))+abs(s%vv(i,j)),abs(s%vev(i,j))))) !Jaap: include sediment advection velocities 00510 00511 if (par%dt<dtold) then 00512 ilim = i 00513 jlim = j 00514 limtype = 2 00515 dtold = par%dt 00516 endif 00517 endif 00518 if(s%wetz(i,j)==1) then 00519 mdx = min(s%dsu(i,j),s%dsz(i,j))**2 00520 mdy = min(s%dnv(i,j),s%dnz(i,j))**2 00521 00522 par%dt=min(par%dt,0.5d0*mdx*mdy/(mdx+mdy)/max(s%nuh(i,j),1e-6)) 00523 00524 if (par%dt<dtold) then 00525 ilim = i 00526 jlim = j 00527 limtype = 3 00528 dtold = par%dt 00529 endif 00530 00531 !Bas: the following criterion is not yet tested for 2D 00532 !par%dt=min(par%dt,0.5d0*mdx/max(s%Dc(i,j)*s%wetz(i,j),1e-6)) 00533 !par%dt=min(par%dt,0.5d0*mdy/max(s%Dc(i,j)*s%wetz(i,j),1e-6)) 00534 endif 00535 enddo 00536 enddo 00537 else 00538 j2=max(s%ny,1) ! Robert: hmmm, in sf1D this should be 1, in "old" 1d this should be 2 00539 do i=imin_ee,imax_ee 00540 if(s%wetu(i,j2)==1) then 00541 ! u-points 00542 mdx=s%dsu(i,j2) 00543 par%dt=min(par%dt,mdx/max(tny,max(sqrt(par%g*s%hu(i,j2))+abs(s%uu(i,j2)),abs(s%ueu(i,j2))))) !Jaap: include sediment advection velocities 00544 if (par%dt<dtold) then 00545 ilim = i 00546 jlim = j2 00547 limtype = 1 00548 dtold = par%dt 00549 endif 00550 ! u-points bed celerity ! Dano 00551 celeru=0.1d0*abs(sum(s%Subg(i,j2,:))+sum(s%Susg(i,j2,:)))*par%morfac/(1.d0-par%por) 00552 par%dt=min(par%dt,mdx/max(tny,celeru)) 00553 endif 00554 if(s%wetz(i,j2)==1) then 00555 ! v-points 00556 mdx=min(s%dsu(i,j2),s%dsu(i-1,j2)) 00557 par%dt=min(par%dt,mdx/max(tny,(sqrt(par%g*s%hv(i,j2))+abs(s%uv(i,j2))))) 00558 if (par%dt<dtold) then 00559 ilim = i 00560 jlim = j2 00561 limtype = 2 00562 dtold = par%dt 00563 endif 00564 00565 mdx = min(s%dsu(i,j2),s%dsz(i,j2))**2 00566 00567 if (par%dy > 0.d0) then 00568 mdy = par%dy 00569 else 00570 mdy = mdx 00571 endif 00572 00573 par%dt=min(par%dt,0.5d0*mdx*mdy/(mdx+mdy)/max(s%nuh(i,j2),1e-6)) 00574 if (par%dt<dtold) then 00575 ilim = i 00576 jlim = j2 00577 limtype = 3 00578 dtold = par%dt 00579 endif 00580 if (par%sedtrans==1 .and. par%facDc>0.d0) then 00581 par%dt=min(par%dt,0.5d0*mdx/max(s%Dc(i,j2),1e-6)) 00582 endif 00583 endif 00584 enddo 00585 endif 00586 00587 if (s%ny>0) then 00588 par%dt=par%dt*par%CFL*0.5d0 00589 else 00590 par%dt=par%dt*par%CFL*0.5d0 00591 endif 00592 00593 ! groundwater stability 00594 if (par%gwflow==1 .and. par%gwnonh==0) then 00595 j2=max(s%ny,1) 00596 do j=j1,j2 00597 do i=2,s%nx 00598 ! Vreugdenhil: 00599 ! 2*kx*gwheight*dt/dx2 < CFL -> dt = (CFL*dx2)/(2*kx*gwheight) 00600 ! par%dt=min(par%dt,par%CFL*s%dsc(i,j)**2/(2*par%kx*s%gwheight(i,j))) 00601 ! 00602 ! Wikipedia: http://en.wikipedia.org/wiki/Von_Neumann_stability_analysis for FTCS: 00603 ! 2*a*dt/dx2 < 1 , where dh/dt = a *d2h/dx2 -> a = k/por 00604 ! 2*k/por*dt/dx2 = CFL -> dt = CFL*dx2*por/2/k 00605 par%dt=min(par%dt,par%CFL*s%dsc(i,j)**2*par%por/(2*par%kx)) 00606 if (par%dt<dtold) then 00607 ilim = i 00608 jlim = j 00609 limtype = 4 00610 dtold = par%dt 00611 endif 00612 enddo 00613 enddo 00614 if (par%ny>2) then 00615 do j=j1,j2 00616 do i=2,s%nx 00617 ! 2*ky*gwheight*dt/dy2 < CFL -> dt = (CFL*dy2)/(2*ky*gwheight) 00618 ! par%dt=min(par%dt,par%CFL*s%dnc(i,j)**2/(2*par%ky*s%gwheight(i,j))) 00619 par%dt=min(par%dt,par%CFL*s%dnc(i,j)**2*par%por/(2*par%ky)) 00620 if (par%dt<dtold) then 00621 ilim = i 00622 jlim = j 00623 limtype = 4 00624 dtold = par%dt 00625 endif 00626 enddo 00627 enddo 00628 endif 00629 endif 00630 00631 if (par%swave==1) then 00632 #ifdef USEMPI 00633 ! Dano: no refraction if ntheta==1 so no need for this check 00634 if (s%ntheta>1) then 00635 par%dt=min(par%dt,par%CFL*s%dtheta/(maxval(maxval(abs(s%ctheta),3)*real(s%wetz))+tny)) 00636 endif 00637 #else 00638 ! if (par%instat(1:4)/='stat') then 00639 ! wwvv: hoping that this is a correct translation: 00640 select case (par%wavemodel) 00641 case(WAVEMODEL_STATIONARY) 00642 continue 00643 case default 00644 par%dt=min(par%dt,par%CFL*s%dtheta/(maxval(maxval(abs(s%ctheta),3)*real(s%wetz))+tny)) 00645 end select 00646 #endif 00647 endif 00648 !To avoid large timestep differences due to output, which can cause instabities 00649 !in the hanssen (leapfrog) scheme, we smooth the timestep. 00650 00651 if (tpar%tnext > par%t) then 00652 ! sometime dt too small and then integer overflow. Catch with check beforehand 00653 if ((tpar%tnext-par%t)/par%dt<dble(0.9*huge(0))) then 00654 n = ceiling((tpar%tnext-par%t)/par%dt) 00655 par%dt = (tpar%tnext-par%t)/n 00656 endif 00657 end if 00658 end if ! first, or other time step 00659 00660 if (dtref .eq. 0.0d0) then 00661 ! If dtref is not yet set. 00662 dtref = par%dt 00663 #ifdef USEMPI 00664 ! Use the same dtref everywhere. 00665 call xmpi_allreduce(dtref,MPI_MIN) 00666 #endif 00667 end if 00668 00669 ! wwvv: In the mpi version par%dt will be calculated different 00670 ! on different processes. So, we take the minimum of all dt's 00671 #ifdef USEMPI 00672 call xmpi_allreduce(par%dt,MPI_MIN) 00673 #endif 00674 00675 ! maximize dt with given dt 00676 if (present(dt)) then 00677 par%dt = min(dt, par%dt) 00678 end if 00679 00680 if (present(limtypeout)) then 00681 limtypeout = limtype 00682 endif 00683 00684 end subroutine compute_dt 00685 00686 subroutine timestep(s,par, tpar, it, dt, ierr) 00687 use params 00688 use spaceparams 00689 use xmpi_module 00690 use logging_module 00691 use paramsconst 00692 00693 IMPLICIT NONE 00694 00695 type(spacepars), intent(inout) :: s 00696 type(parameters), intent(inout) :: par 00697 type(timepars), intent(inout) :: tpar 00698 integer, intent(inout) :: it 00699 integer, intent(out), optional :: ierr 00700 real*8, intent(in), optional :: dt 00701 integer :: ilim 00702 integer :: jlim 00703 integer :: limtype 00704 real*8 :: fac 00705 real*8,save :: dtref 00706 00707 if(.not. xcompute) return 00708 00709 call compute_dt(s,par, tpar, it, ilim, jlim, dtref, dt=dt, ierr=ierr,limtypeout=limtype) 00710 00711 par%t=par%t+par%dt 00712 00713 if(par%t >= tpar%tnext) then 00714 par%dt = par%dt-(par%t-tpar%tnext) 00715 par%t = tpar%tnext 00716 it = it+1 00717 end if 00718 00719 ! Sometimes dtref is set incorrectly in intialization. We want to stop sudden changes 00720 ! in time step size, so allow dtref to change slowely to current situations. If simulation 00721 ! requires smaller timestep due to e.g. overwash, this should not stop the simulation. 00722 ! Errors occur in O(1-10) timesteps. Therefore filter on O(1000-10000) timesteps 00723 fac = 1.d0/5000 00724 dtref = (1.d0-fac)*dtref + fac*par%dt 00725 00726 if (((dtref/par%dt>par%maxdtfac .and. par%t<tpar%tnext) .or. par%dt/dtref>par%maxdtfac) & 00727 .and. par%defuse==1) then 00728 00729 call writelog('lswe','','Quit XBeach since computational time implodes/explodes') 00730 call writelog('lswe','','Please check output at the end of the simulation') 00731 select case (limtype) 00732 case (0) 00733 call writelog('lswe','(a)','Time constrain criterium exceeded in all cells') 00734 case (1) 00735 call writelog('lswe','(a,i0,a,i0,a)','U-velocities are too high in cell (',ilim,',',jlim,')(M,N)') 00736 case (2) 00737 call writelog('lswe','(a,i0,a,i0,a)','V-velocities are too high in cell (',ilim,',',jlim,')(M,N)') 00738 case (3) 00739 call writelog('lswe','(a,i0,a,i0,a)','Viscosity condition is too high in cell (',ilim,',',jlim,')(M,N)') 00740 case (4) 00741 call writelog('lswe','(a,i0,a,i0,a)','Groundwater condition is too high in cell (',ilim,',',jlim,')(M,N)') 00742 end select 00743 call writelog('lswe','','dtref: ',dtref) 00744 call writelog('lswe','','par%dt: ',par%dt) 00745 00746 ! Force output before exit in main time loop 00747 tpar%outputg = (size(tpar%tpg) .gt. 0) 00748 tpar%outputp = (size(tpar%tpp) .gt. 0) 00749 tpar%outputm = (size(tpar%tpm) .gt. 0) 00750 tpar%outputh = (size(tpar%tph) .gt. 0) 00751 ! set output time to current time 00752 if (tpar%outputg) tpar%tpg=min(tpar%tpg,par%t) 00753 if (tpar%outputp) tpar%tpp=min(tpar%tpp,par%t) 00754 if (tpar%outputm) tpar%tpm=min(tpar%tpm,par%t) 00755 if (tpar%outputh) tpar%tph=min(tpar%tph,par%t) 00756 tpar%output = (tpar%outputg .or. tpar%outputp .or. tpar%outputm .or. tpar%outputh) 00757 if (present(ierr)) then 00758 ierr = 1 00759 endif 00760 endif 00761 00762 end subroutine timestep 00763 end module timestep_module
1.7.6.1
