|
XBeach
|
00001 module boundaryconditions 00002 use typesandkinds, only: slen 00003 use constants, only: pi 00004 implicit none 00005 save 00006 private 00007 public flow_lat_bc, wave_bc, flow_bc, discharge_boundary_h, discharge_boundary_v 00008 contains 00009 subroutine wave_bc(sg,sl,par) 00010 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00011 ! Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University ! 00012 ! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski, ! 00013 ! Jaap van Thiel de Vries, Robert McCall ! 00014 ! ! 00015 ! d.roelvink@unesco-ihe.org ! 00016 ! UNESCO-IHE Institute for Water Education ! 00017 ! P.O. Box 3015 ! 00018 ! 2601 DA Delft ! 00019 ! The Netherlands ! 00020 ! ! 00021 ! This library is free software; you can redistribute it and/or ! 00022 ! modify it under the terms of the GNU Lesser General Public ! 00023 ! License as published by the Free Software Foundation; either ! 00024 ! version 2.1 of the License, or (at your option) any later version. ! 00025 ! ! 00026 ! This library is distributed in the hope that it will be useful, ! 00027 ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! 00028 ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! 00029 ! Lesser General Public License for more details. ! 00030 ! ! 00031 ! You should have received a copy of the GNU Lesser General Public ! 00032 ! License along with this library; if not, write to the Free Software ! 00033 ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 ! 00034 ! USA ! 00035 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00036 use params, only: parameters 00037 use spaceparamsdef 00038 use interp, only: linear_interp 00039 use wave_functions_module, only: iteratedispersion, dispersion 00040 use xmpi_module 00041 use logging_module, only: writelog, report_file_read_error 00042 use nonh_module, only: nonh_init_wcoef 00043 use spectral_wave_bc_module, only: spectral_wave_bc 00044 use filefunctions, only: create_new_fid 00045 use paramsconst 00046 #ifdef USEMPI 00047 use spaceparams, only: space_distribute 00048 #endif 00049 00050 00051 implicit none 00052 00053 type(spacepars), target :: sg, sl 00054 type(spacepars), pointer :: s 00055 type(parameters) :: par 00056 00057 integer, save :: nt 00058 integer :: i,new,reclen,wordsize 00059 integer, save :: old 00060 integer, save :: recpos, curline 00061 integer :: j 00062 integer :: itheta 00063 integer :: E_idx 00064 integer :: ier,ier2 00065 real*8 :: E1,ei,dum,Hm0, dum1, spreadpar, bcdur, dum2 00066 real*8, save :: dtbcfile,rt,bcendtime,bcstarttime 00067 real*8 :: em,tshifted,tnew 00068 real*8, save :: Emean,Llong 00069 real*8,dimension(:) ,allocatable,save :: e01,L0,L,Lest,kbw,wbw,tanhkhwb,kxmwt 00070 real*8,dimension(:) ,allocatable,save :: fac1,fac2 00071 real*8,dimension(:) ,allocatable,save :: tE,dataE,databi 00072 real*8,dimension(:,:) ,allocatable,save :: ht 00073 real*8,dimension(:,:) ,allocatable,save :: q1,q2,q 00074 real*8,dimension(:) ,allocatable,save :: wcrestpos 00075 real*8,dimension(:,:) ,allocatable,save :: ee1, ee2 00076 real*8,dimension(:,:) ,allocatable,save :: gq1,gq2,gq 00077 real*8,dimension(:,:) ,allocatable,save :: gee1, gee2 00078 character(len=1) :: bl 00079 character(slen),save :: ebcfname,qbcfname,nhbcfname,esbcfname 00080 real*8 :: E0,arms 00081 real*8,dimension(:),allocatable,save :: dist,factor 00082 logical :: startbcf 00083 logical :: isSet_U,isSet_V,isSet_Z,isSet_W,isSet_dU,isSet_dV,isSet_Q 00084 real*8,dimension(:) ,allocatable,save :: uig,vig,zig,wig,duig,dvig 00085 real*8,dimension(:) ,allocatable,save :: tempu,tempv,tempdu,tempdv 00086 00087 ! Initialize to false only once.... 00088 logical, save :: bccreated = .false. 00089 00090 s=>sl 00091 #ifndef USEMPI 00092 s=>sg 00093 #endif 00094 00095 if (.not. allocated(fac1)) then 00096 allocate(ht (2,s%ny+1)) 00097 allocate(fac1 (s%nx)) 00098 allocate(fac2 (s%nx)) 00099 allocate(e01(1:s%ntheta)) 00100 allocate(dist(1:s%ntheta)) 00101 allocate(factor(1:s%ntheta)) 00102 allocate(wcrestpos(s%nx+1)) 00103 allocate(L(s%ny+1)) 00104 allocate(L0(s%ny+1)) 00105 allocate(Lest(s%ny+1)) 00106 L0 = par%g*par%Trep**2/2/par%px 00107 L = L0 00108 allocate(kbw(s%ny+1)) 00109 allocate(tanhkhwb(s%ny+1)) 00110 allocate(kxmwt(s%ny+1)) 00111 allocate(wbw(s%ny+1)) 00112 wbw = 2*par%px/par%Trep 00113 endif 00114 ! 00115 ! BLOCK I: GENERATE AND READ-IN WAVE BOUNDARY CONDITIONS 00116 ! 00117 startbcf=.false. 00118 if(.not. bccreated ) then 00119 if (xmaster) then 00120 call writelog('ls','','Setting up boundary conditions') 00121 endif 00122 bccreated=.true. 00123 startbcf=.true. ! trigger read from bcf for instat 3,4,5,7 00124 bcendtime=huge(0.0d0) ! initial assumption for instat 3,4,5,7 00125 s%newstatbc=1 00126 ! Part A) read-in / generate different wbctypes (and waveforms) 00127 ! 1) TS_1 00128 if (par%wbctype==WBCTYPE_TS_1) then 00129 if(xmaster) then 00130 open( unit=7, file='bc/gen.ezs') 00131 endif 00132 if (xmaster) then 00133 5 continue 00134 read(7,'(a)',iostat=ier)bl 00135 if (ier .ne. 0) then 00136 call report_file_read_error('bc/gen.ezs') 00137 endif 00138 if(bl.eq.'*') goto 5 00139 read(7,*,iostat=ier)nt 00140 if (ier .ne. 0) then 00141 call report_file_read_error('bc/gen.ezs') 00142 endif 00143 endif 00144 #ifdef USEMPI 00145 call xmpi_bcast(nt) 00146 #endif 00147 allocate(dataE (nt)) 00148 allocate(tE (nt)) 00149 do i=1,nt 00150 if(xmaster) then 00151 read(7,*,iostat=ier) tE(i),dum,dataE(i) 00152 if (ier .ne. 0) then 00153 call report_file_read_error('bc/gen.ezs') 00154 endif 00155 endif 00156 #ifdef USEMPI 00157 call xmpi_bcast(tE(i)) 00158 call xmpi_bcast(dataE(i)) 00159 #endif 00160 end do 00161 if (xmaster) then 00162 close(7) 00163 endif 00164 Emean=sum(dataE)/nt 00165 call dispersion(par,s,s%hh) 00166 ! 2) TS_2 00167 elseif (par%wbctype==WBCTYPE_TS_2) then 00168 if (xmaster) then 00169 open( unit=7, file='bc/gen.ezs') 00170 6 continue 00171 read(7,'(a)',iostat=ier)bl 00172 if (ier .ne. 0) then 00173 call report_file_read_error('bc/gen.ezs') 00174 endif 00175 if(bl.eq.'*') goto 6 00176 read(7,*,iostat=ier)nt 00177 if (ier .ne. 0) then 00178 call report_file_read_error('bc/gen.ezs') 00179 endif 00180 endif 00181 #ifdef USEMPI 00182 call xmpi_bcast(nt) 00183 #endif 00184 allocate(dataE (nt)) 00185 allocate(databi (nt)) 00186 allocate(tE (nt)) 00187 do i=1,nt 00188 if(xmaster) then 00189 read(7,*,iostat=ier) tE(i),databi(i),dataE(i) 00190 if (ier .ne. 0) then 00191 call report_file_read_error('bc/gen.ezs') 00192 endif 00193 endif 00194 #ifdef USEMPI 00195 call xmpi_bcast(tE(i)) 00196 call xmpi_bcast(databi(i)) 00197 call xmpi_bcast(dataE(i)) 00198 #endif 00199 end do 00200 if (xmaster) then 00201 close(7) 00202 endif 00203 Emean=sum(dataE)/nt 00204 call dispersion(par,s,s%hh) 00205 ! 3) jons_table in combination with stationary 00206 elseif (par%wbctype==WBCTYPE_JONS_TABLE .and. par%wavemodel==WAVEMODEL_STATIONARY) then 00207 if (xmaster) then 00208 open( unit=7, file=par%bcfile) 00209 read(7,*,iostat=ier) Hm0, par%Trep,par%dir0, dum1, spreadpar, bcendtime, dum2 00210 if (ier .ne. 0) then 00211 call report_file_read_error(par%bcfile) 00212 endif 00213 par%Hrms = Hm0/sqrt(2.d0) 00214 par%m = int(2*spreadpar) 00215 if (par%morfacopt==1) bcendtime=bcendtime/max(par%morfac,1.d0) 00216 s%theta0=(1.5d0*par%px)-par%dir0*atan(1.d0)/45.d0 00217 do while(s%theta0<-par%px) 00218 s%theta0=s%theta0+2.d0*par%px 00219 enddo 00220 do while(s%theta0>par%px) 00221 s%theta0=s%theta0-2.d0*par%px 00222 enddo 00223 endif 00224 s%newstatbc=1 00225 #ifdef USEMPI 00226 call xmpi_bcast(bcendtime) 00227 call xmpi_bcast(par%Hrms) 00228 call xmpi_bcast(par%Trep) 00229 call xmpi_bcast(par%m) 00230 call xmpi_bcast(s%theta0) 00231 #endif 00232 do itheta=1,s%ntheta 00233 s%sigt(:,:,itheta) = 2.d0*par%px/par%Trep 00234 end do 00235 s%sigm = sum(s%sigt,3)/s%ntheta 00236 call dispersion(par,s,s%hh) 00237 00238 ! 4) jons_table in combination with surfbeat or nonh 00239 elseif ((par%wbctype==WBCTYPE_PARAMETRIC .or. par%wbctype==WBCTYPE_JONS_TABLE) .and. & 00240 (par%wavemodel/=WAVEMODEL_STATIONARY .and. xmaster)) then 00241 call spectral_wave_bc(sg,par,curline) 00242 ! 5) SWAN 00243 elseif (par%wbctype==WBCTYPE_SWAN.and.xmaster) then 00244 call spectral_wave_bc(sg,par,curline) 00245 ! 6) vardens 00246 elseif (par%wbctype==WBCTYPE_VARDENS.and.xmaster) then 00247 call spectral_wave_bc(sg,par,curline) 00248 ! 7) reuse 00249 elseif (par%wbctype==WBCTYPE_REUSE.and.xmaster) then 00250 curline = 1 00251 ! 8) ts_nonh 00252 elseif (par%wbctype==WBCTYPE_TS_NONH) then 00253 if (.not. allocated(uig)) then 00254 if (xmaster) then 00255 allocate(uig(sg%ny+1)) 00256 allocate(vig(sg%ny+1)) 00257 allocate(zig(sg%ny+1)) 00258 allocate(wig(sg%ny+1)) 00259 allocate(duig(sg%ny+1)) 00260 allocate(dvig(sg%ny+1)) 00261 else ! only need valid address for MPI-distribution 00262 allocate(uig(1)) 00263 allocate(vig(1)) 00264 allocate(zig(1)) 00265 allocate(wig(1)) 00266 allocate(duig(1)) 00267 allocate(dvig(1)) 00268 endif 00269 endif 00270 if (xmaster) then 00271 call velocity_Boundary('boun_U.bcf',uig,vig,zig,wig,duig,dvig,sg%ny,par%t, & 00272 isSet_U,isSet_V,isSet_Z,isSet_W,isSet_dU,isSet_dV,isSet_Q) 00273 endif 00274 endif ! done with the different wbctypes 00275 ! 00276 ! Part B) Directional distribution for non-spectral wbctypes 00277 if (par%wbctype==WBCTYPE_PARAMS .or. par%wbctype==WBCTYPE_TS_1 .or. par%wbctype==WBCTYPE_TS_2 .or. & 00278 (par%wbctype==WBCTYPE_JONS_TABLE .and. par%wavemodel==WAVEMODEL_STATIONARY) ) then 00279 dist=(cos(s%theta-s%theta0))**par%m 00280 do i=1,s%ntheta 00281 if(cos(s%theta(i)-s%theta0)<0.d0) then 00282 dist(i)=0.0d0 00283 end if 00284 end do 00285 if (par%wbctype==WBCTYPE_TS_1 .or. par%wbctype==WBCTYPE_TS_2) then 00286 par%Hrms=sqrt(8*Emean/(par%rho*par%g)) 00287 endif 00288 E0=par%rhog8*par%Hrms**2 00289 ! Energy density distribution 00290 if (sum(dist)>0.d0) then 00291 factor = (dist/sum(dist))/s%dtheta 00292 else 00293 factor=0.d0 00294 endif 00295 e01 = factor*E0; 00296 e01 = max(e01,0.0d0); 00297 Llong = par%Tlong*s%cg(1,1) 00298 endif 00299 if (xmaster) then 00300 call writelog('sl','','Boundary conditions complete, starting computation') 00301 endif 00302 end if ! if bc was not created 00303 ! 00304 ! BLOCK II: Recalculate bcf-file 00305 ! 00306 if (par%t .ge. bcendtime) then 00307 ! Go over the different wbctypes (and wavemodel) 00308 ! 1) jons_table in combination with stationary 00309 if (par%wbctype==WBCTYPE_JONS_TABLE .and. par%wavemodel==WAVEMODEL_STATIONARY) then 00310 if (xmaster) then 00311 call writelog('ls','','Reading new wave conditions') 00312 ! read new wave conditions until bcendtime exceeds current time step 00313 ! necessary in case of external time management 00314 do while (par%t .ge. bcendtime) 00315 read(7,*,iostat=ier) Hm0, par%Trep,par%dir0, dum1, spreadpar, bcdur, dum2 00316 if (ier .ne. 0) then 00317 call report_file_read_error(par%bcfile) 00318 endif 00319 par%Hrms = Hm0/sqrt(2.d0) 00320 ! set taper time to 1 second for new conditions (likewise in waveparams) 00321 par%taper = 0.d0 00322 par%m = int(2*spreadpar) 00323 if (par%morfacopt==1) then 00324 bcendtime=bcendtime+bcdur/max(par%morfac,1.d0) 00325 else 00326 bcendtime=bcendtime+bcdur 00327 endif 00328 s%theta0=(1.5d0*par%px)-par%dir0*atan(1.d0)/45.d0 00329 ! Make shore theta0 is also set between -px and px for new wave conditions 00330 do while(s%theta0<-par%px) 00331 s%theta0=s%theta0+2.d0*par%px 00332 enddo 00333 do while(s%theta0>par%px) 00334 s%theta0=s%theta0-2.d0*par%px 00335 enddo 00336 enddo 00337 endif 00338 s%newstatbc=1 00339 #ifdef USEMPI 00340 call xmpi_bcast(bcendtime) 00341 call xmpi_bcast(par%Hrms) 00342 call xmpi_bcast(par%Trep) 00343 call xmpi_bcast(par%m) 00344 call xmpi_bcast(s%theta0) 00345 #endif 00346 do itheta=1,s%ntheta 00347 s%sigt(:,:,itheta) = 2*par%px/par%Trep 00348 end do 00349 s%sigm = sum(s%sigt,3)/s%ntheta 00350 call dispersion(par,s,s%hh) 00351 ! Directional distribution 00352 dist=(cos(s%theta-s%theta0))**par%m 00353 do i=1,s%ntheta 00354 if(cos(s%theta(i)-s%theta0)<0.d0) then 00355 dist(i)=0 00356 end if 00357 end do 00358 E0=par%rhog8*par%Hrms**2 00359 ! Energy density distribution 00360 if (sum(dist)>0.d0) then 00361 factor = (dist/sum(dist))/s%dtheta 00362 else 00363 factor=0.d0 00364 endif 00365 e01 = factor*E0; 00366 e01 = max(e01,0.0d0); 00367 ! 2) JONSWAP (table) 00368 elseif ((par%wbctype==WBCTYPE_PARAMETRIC .or. par%wbctype==WBCTYPE_JONS_TABLE) .and. & 00369 (par%wavemodel/=WAVEMODEL_STATIONARY .and. xmaster)) then 00370 close(71) 00371 close(72) 00372 call spectral_wave_bc(sg,par,curline) 00373 startbcf=.true. 00374 ! 3) SWAN 00375 elseif (par%wbctype==WBCTYPE_SWAN.and.xmaster) then 00376 close(71) 00377 close(72) 00378 call spectral_wave_bc(sg,par,curline) 00379 startbcf=.true. 00380 ! 4) vardens 00381 elseif (par%wbctype==WBCTYPE_VARDENS.and.xmaster) then 00382 close(71) 00383 close(72) 00384 call spectral_wave_bc(sg,par,curline) 00385 startbcf=.true. 00386 ! 5) reuse 00387 elseif (par%wbctype==WBCTYPE_REUSE.and.xmaster) then 00388 close(71) 00389 close(72) 00390 startbcf=.true. 00391 if (par%t <= (par%tstop-par%dt)) then 00392 ! this is not compatible with external time management 00393 curline = curline + 1 00394 end if 00395 end if ! go over the different wbctypes 00396 #ifdef USEMPI 00397 call xmpi_bcast(startbcf) 00398 #endif 00399 end if ! done recalculating bcf-file 00400 ! 00401 ! BLOCK III: COMPUTE WAVE BOUNDARY CONDITIONS CURRENT TIMESTEP 00402 ! 00403 ! Part A) Stationary: go over the different waveforms 00404 if (par%wavemodel==WAVEMODEL_STATIONARY) then 00405 ! Go over all excepted wbctypes (all the same) 00406 if (par%wbctype == WBCTYPE_PARAMS .or. par%wbctype == WBCTYPE_JONS_TABLE) then 00407 ! Dano reset waves for new wave condition in stationary case, to avoid 'leftovers' 00408 if (s%newstatbc==1) s%ee=0.d0 00409 if(par%taper>tiny(0.d0)) then 00410 do j=1,s%ny+1 00411 s%ee(1,j,:) = e01*min(par%t/par%taper,1.0d0) 00412 enddo 00413 else 00414 do j=1,s%ny+1 00415 s%ee(1,j,:) = e01 00416 enddo 00417 endif 00418 s%bi(1) = 0.d0 00419 s%ui(1,:) = 0.0d0 00420 else 00421 endif ! Go over the different wbctypes for wavemodel = stationary 00422 ! 00423 ! Part B) Surfbeat: go over the different waveforms 00424 elseif (par%wavemodel == WAVEMODEL_SURFBEAT) then 00425 ! 1) Bichromatic waves 00426 if ( (par%wbctype == WBCTYPE_PARAMS) .and. (par%Tlong /= -123) .and. xmpi_istop ) then 00427 do j=1,s%ny+1 00428 s%ee(1,j,:)=e01*0.5d0 * & 00429 (1.d0+cos(2*par%px*(par%t/par%Tlong-( sin(s%theta0)*(s%yz(1,j)-s%yz(1,1)) & !jaap : -sum(s%alfaz(1,1:j))/j 00430 +cos(s%theta0)*(s%xz(1,j)-s%xz(1,1)) )/Llong))) * & !jaap: -sum(s%alfaz(1,1:j))/j 00431 min(par%t/par%taper,1.d0) 00432 em = (sum(0.5d0*e01))*s%dtheta *min(par%t/par%taper,1.d0) 00433 ei = sum(s%ee(1,j,1:s%ntheta))*s%dtheta 00434 s%bi(1) = -(2*s%cg(1,j)/s%c(1,j)-0.5d0)*(em-ei)/(s%cg(1,j)**2-par%g*s%hh(1,j))/par%rho 00435 ht=s%zs0(1:2,:)-s%zb(1:2,:) 00436 s%ui(1,j) = s%cg(1,j)*s%bi(1)/ht(1,j)*cos(s%theta0-s%alfaz(1,j)) 00437 end do 00438 ! 2) TS_1: bound long wave is calculate by the bound long wave theory of LG & Steward (1964) 00439 elseif (par%wbctype==WBCTYPE_TS_1 .and. xmpi_istop) then 00440 do j=1,s%ny+1 00441 if (abs(s%theta0-sum(s%alfaz(1,1:j))/j)<1e-3) then 00442 call linear_interp(tE,dataE,nt,par%t,E1,E_idx) 00443 else 00444 tshifted = max(par%t-(s%yz(1,j)-s%yz(1,1))*sin(s%theta0)/s%cg(1,1) & !-sum(s%alfaz(1,1:j))/j 00445 +(s%xz(1,j)-s%xz(1,1))*cos(s%theta0)/s%cg(1,1),0.d0) !-sum(s%alfaz(1,1:j))/j 00446 call linear_interp(tE,dataE,nt,tshifted,E1,E_idx) 00447 endif 00448 s%ee(1,j,:)=e01*E1/max(Emean,0.000001d0)*min(par%t/par%taper,1.d0) 00449 em = Emean *min(par%t/par%taper,1.d0) 00450 ei = sum(s%ee(1,j,:))*s%dtheta 00451 s%bi(1) = -(2*s%cg(1,j)/s%c(1,j)-0.5d0)*(em-ei)/(s%cg(1,j)**2-par%g*s%hh(1,j))/par%rho 00452 ht=s%zs0(1:2,:)-s%zb(1:2,:) 00453 s%ui(1,j) = s%cg(1,j)*s%bi(1)/ht(1,j)*cos(s%theta0-s%alfaz(1,j)) 00454 end do 00455 ! 3) TS_2: bound long wave is specificed by the user 00456 elseif (par%wbctype==WBCTYPE_TS_2 .and. xmpi_istop) then 00457 ht=s%zs0(1:2,:)-s%zb(1:2,:) 00458 do j=1,s%ny+1 00459 if (abs(s%theta0-sum(s%alfaz(1,1:j))/j)<1e-3) then 00460 call linear_interp(tE,dataE,nt,par%t,E1,E_idx) 00461 call linear_interp(tE,databi,nt,par%t,s%bi(1),E_idx) 00462 else 00463 tshifted = max(par%t-(s%yz(1,j)-s%yz(1,1))*sin(s%theta0)/s%cg(1,1) & !-sum(s%alfaz(1,1:j))/j 00464 +(s%xz(1,j)-s%xz(1,1))*cos(s%theta0)/s%cg(1,1),0.d0) !-sum(s%alfaz(1,1:j))/j 00465 call linear_interp(tE,dataE,nt,tshifted,E1,E_idx) 00466 call linear_interp(tE,databi,nt,tshifted,s%bi(1),E_idx) 00467 endif 00468 s%ee(1,j,:)=e01*E1/max(Emean,0.000001d0)*min(par%t/par%taper,1.d0) 00469 if (par%freewave==1) then 00470 s%ui(1,j) = sqrt(par%g/ht(1,j))*s%bi(1) 00471 else 00472 s%ui(1,j) = s%cg(1,j)*s%bi(1)/ht(1,j)*cos(s%theta0-s%alfaz(1,j))*min(par%t/par%taper,1.d0) 00473 endif 00474 end do 00475 ! 4) go over the different wbctypes (spectral) 00476 else 00477 if ((par%wbctype==WBCTYPE_PARAMETRIC) .or. (par%wbctype==WBCTYPE_JONS_TABLE) .or. & 00478 (par%wbctype==WBCTYPE_SWAN) .or. (par%wbctype==WBCTYPE_VARDENS) .or. (par%wbctype==WBCTYPE_REUSE)) then 00479 ! open file if first time 00480 if (startbcf) then 00481 if(xmaster) then 00482 open(53,file='ebcflist.bcf',form='formatted',position='rewind') 00483 open(54,file='qbcflist.bcf',form='formatted',position='rewind') 00484 if(par%single_dir==1 .and. par%wbctype==WBCTYPE_REUSE) then 00485 open(55,file='esbcflist.bcf',form='formatted',position='rewind') 00486 endif 00487 endif 00488 if (xmaster) then 00489 do i=1,curline 00490 read(53,*,iostat=ier)bcendtime,rt,dtbcfile,par%Trep,s%theta0,ebcfname 00491 if (ier .ne. 0) then 00492 call report_file_read_error('ebcflist.bcf') 00493 endif 00494 read(54,*,iostat=ier)bcendtime,rt,dtbcfile,par%Trep,s%theta0,qbcfname 00495 if (ier .ne. 0) then 00496 call report_file_read_error('qbcflist.bcf') 00497 endif 00498 if(par%single_dir==1 .and. par%wbctype==WBCTYPE_REUSE) then 00499 read(55,*,iostat=ier)bcendtime,rt,dtbcfile,par%Trep,s%theta0,esbcfname 00500 if (ier .ne. 0) then 00501 call report_file_read_error('esbcflist.bcf') 00502 endif 00503 endif 00504 enddo ! wwvv strange 00505 if(par%single_dir==1 .and. par%wbctype==WBCTYPE_REUSE) then 00506 inquire(iolength=wordsize) 1.d0 00507 reclen=wordsize*(sg%ny+1)*(sg%ntheta_s) 00508 open(73,file=esbcfname,status='old',form='unformatted',access='direct',recl=reclen) 00509 read(73,rec=1,iostat=ier)sg%ee_s(1,:,:) 00510 if (ier .ne. 0) then 00511 call report_file_read_error(esbcfname) 00512 endif 00513 close(73) 00514 endif 00515 endif 00516 #ifdef USEMPI 00517 call xmpi_bcast(bcendtime) 00518 call xmpi_bcast(rt) 00519 call xmpi_bcast(dtbcfile) 00520 call xmpi_bcast(par%Trep) 00521 call xmpi_bcast(s%theta0) 00522 call xmpi_bcast(ebcfname) 00523 if (par%single_dir==1) then 00524 call space_distribute("y",sl,sg%ee_s(1,:,:),s%ee_s(1,:,:)) 00525 endif 00526 #else 00527 if (par%single_dir==1) then 00528 s%ee_s=sg%ee_s 00529 endif 00530 #endif 00531 if (xmaster) then 00532 close(53) 00533 close(54) 00534 if(par%single_dir==1 .and. par%wbctype==WBCTYPE_REUSE) then 00535 close(55) 00536 endif 00537 endif 00538 ! Robert and Jaap : Initialize for new wave conditions 00539 do itheta=1,s%ntheta 00540 s%sigt(:,:,itheta) = 2*par%px/par%Trep 00541 end do 00542 s%sigm = sum(s%sigt,3)/s%ntheta 00543 call dispersion(par,s,s%hh) 00544 ! End initialize 00545 if (xmaster) then 00546 inquire(iolength=wordsize) 1.d0 00547 reclen=wordsize*(sg%ny+1)*(sg%ntheta) 00548 open(71,file=ebcfname,status='old',form='unformatted',access='direct',recl=reclen) 00549 reclen=wordsize*((sg%ny+1)*4) 00550 open(72,file=qbcfname,status='old',form='unformatted',access='direct',recl=reclen) 00551 endif 00552 ! wwvv note that we need the global value of ny here 00553 ! masterprocess reads and distributes 00554 if (xmaster) then 00555 if (.not. allocated(gq1) ) then 00556 allocate(gq1(sg%ny+1,4),gq2(sg%ny+1,4),gq(sg%ny+1,4)) 00557 allocate(gee1(sg%ny+1,s%ntheta),gee2(sg%ny+1,s%ntheta)) 00558 endif 00559 else 00560 if (.not. allocated(gq1) ) then ! to get valid addresses for 00561 ! gq1, gq2, gq, gee1, gee2 00562 allocate(gq1(1,4),gq2(1,4),gq(1,4)) 00563 allocate(gee1(1,s%ntheta),gee2(1,s%ntheta)) 00564 endif 00565 endif 00566 if (.not. allocated(q1) ) then 00567 allocate(q1(s%ny+1,4),q2(s%ny+1,4),q(s%ny+1,4)) 00568 allocate(ee1(s%ny+1,s%ntheta),ee2(s%ny+1,s%ntheta)) 00569 end if 00570 if (xmaster) then 00571 read(71,rec=1,iostat=ier)gee1 ! Earlier in time 00572 read(71,rec=2,iostat=ier2)gee2 ! Later in time 00573 if (ier+ier2 .ne. 0) then 00574 call report_file_read_error(ebcfname) 00575 endif 00576 read(72,rec=1,iostat=ier)gq1 ! Earlier in time 00577 read(72,rec=2,iostat=ier2)gq2 ! Later in time 00578 if (ier+ier2 .ne. 0) then 00579 call report_file_read_error(qbcfname) 00580 endif 00581 endif 00582 #ifdef USEMPI 00583 call space_distribute("y",sl,gee1,ee1) 00584 call space_distribute("y",sl,gee2,ee2) 00585 call space_distribute("y",sl,gq1,q1) 00586 call space_distribute("y",sl,gq2,q2) 00587 #else 00588 ee1=gee1 00589 ee2=gee2 00590 q1=gq1 00591 q2=gq2 00592 #endif 00593 old=floor((par%t/dtbcfile)+1) 00594 recpos=1 00595 end if 00596 new=floor((par%t/dtbcfile)+1) 00597 ! Check for next level in boundary condition file 00598 if (new/=old) then 00599 recpos=recpos+(new-old) 00600 ! Check for how many bcfile steps are jumped 00601 if (new-old>1) then ! Many steps further in the bc file 00602 if(xmaster) then 00603 read(72,rec=recpos+1,iostat=ier)gq2 00604 if (ier .ne. 0) then 00605 call report_file_read_error(qbcfname) 00606 endif 00607 read(71,rec=recpos+1,iostat=ier)gee2 00608 if (ier .ne. 0) then 00609 call report_file_read_error(ebcfname) 00610 endif 00611 read(72,rec=recpos,iostat=ier)gq1 00612 if (ier .ne. 0) then 00613 call report_file_read_error(qbcfname) 00614 endif 00615 read(71,rec=recpos,iostat=ier)gee1 00616 if (ier .ne. 0) then 00617 call report_file_read_error(ebcfname) 00618 endif 00619 endif 00620 #ifdef USEMPI 00621 call space_distribute("y",sl,gee2,ee2) 00622 call space_distribute("y",sl,gq2,q2) 00623 call space_distribute("y",sl,gee1,ee1) 00624 call space_distribute("y",sl,gq1,q1) 00625 #else 00626 ee1=gee1 00627 ee2=gee2 00628 q1=gq2 00629 q2=gq2 00630 #endif 00631 else ! Only one step further in the bc file 00632 ee1=ee2 00633 q1=q2 00634 if(xmaster) then 00635 read(72,rec=recpos+1,iostat=ier)gq2 00636 if (ier .ne. 0) then 00637 call report_file_read_error(qbcfname) 00638 endif 00639 read(71,rec=recpos+1,iostat=ier)gee2 00640 if (ier .ne. 0) then 00641 call report_file_read_error(ebcfname) 00642 endif 00643 endif 00644 #ifdef USEMPI 00645 call space_distribute("y",sl,gee2,ee2) 00646 call space_distribute("y",sl,gq2,q2) 00647 #else 00648 ee2=gee2 00649 q2=gq2 00650 #endif 00651 endif 00652 old=new 00653 end if 00654 ht=s%zs0(1:2,:)-s%zb(1:2,:) 00655 tnew = dble(new)*dtbcfile 00656 if (xmpi_istop) s%ee(1,:,:) = (dtbcfile-(tnew-par%t))/dtbcfile*ee2 + & !Jaap 00657 (tnew-par%t)/dtbcfile*ee1 00658 q = (dtbcfile-(tnew-par%t))/dtbcfile*q2 + & !Jaap 00659 (tnew-par%t)/dtbcfile*q1 00660 ! be aware: ui and vi are defined w.r.t. the grid, not w.r.t. the coordinate system! 00661 s%ui(1,:) = (q(:,1)*dcos(s%alfau(1,:)) + q(:,2)*dsin(s%alfau(1,:)))/ht(1,:)*min(par%t/par%taper,1.0d0) 00662 s%vi(1,:) = (-q(:,1)*dsin(s%alfau(1,:)) + q(:,2)*dcos(s%alfau(1,:)))/ht(1,:)*min(par%t/par%taper,1.0d0) 00663 if (xmpi_istop) s%ee(1,:,:)=s%ee(1,:,:)*min(par%t/par%taper,1.0d0) 00664 else 00665 endif 00666 endif ! Go over the different wbctypes for surfbeat 00667 ! Part C) Nonh: go over the different waveforms 00668 elseif (par%wavemodel == WAVEMODEL_NONH) then 00669 ! 1) monochromatic waves 00670 if (par%wbctype==WBCTYPE_PARAMS) then 00671 wbw = 2*par%px/par%Trep 00672 do j=1,s%ny+1 00673 Lest(j) = L(j) 00674 L(j) = iteratedispersion(L0(j),Lest(j),par%px,s%hh(1,j)) 00675 enddo 00676 kbw = 2*par%px/L 00677 arms = par%Hrms/2 00678 if(par%order==1) then 00679 s%zi(1,:) = arms*dsin(wbw*par%t-kbw*( dsin(s%theta0)*(s%yz(1,:)-s%yz(1,1)) & 00680 +dcos(s%theta0)*(s%xz(1,:)-s%xz(1,1)) & 00681 ) & 00682 ) 00683 s%ui(1,:) = 1.d0/s%hh(1,:)*wbw*arms/sinh(kbw*s%hh(1,:)) * & 00684 dsin(wbw*par%t & 00685 -kbw*( dsin(s%theta0)*(s%yz(1,:)-s%yz(1,1)) & 00686 +dcos(s%theta0)*(s%xz(1,:)-s%xz(1,1)) & 00687 ) & 00688 ) * & 00689 1.d0/kbw*sinh(kbw*s%hh(1,:)) *dcos(s%theta0-s%alfaz(1,:)) 00690 elseif(par%order==2) then 00691 ! 2nd order Stokes wave: To do more extensive testing 00692 tanhkhwb = tanh(kbw*s%hh(1,:)) 00693 kxmwt = kbw*(dsin(s%theta0)*(s%yz(1,:)-s%yz(1,1))+dcos(s%theta0)*(s%xz(1,:)-s%xz(1,1))) - wbw*par%t 00694 s%zi(1,:) = arms*(dcos(kxmwt) + kbw*arms*(3-tanhkhwb**2)/(4*tanhkhwb**3)*dcos(2*kxmwt)) 00695 s%ui(1,:) = (wbw/kbw)/s%hh(1,:)*s%zi(1,:) *dcos(s%theta0-s%alfaz(1,:)) 00696 elseif(par%order==3) then 00697 ! 3rd order Stokes wave: only deep water 00698 !s%zi(1,:) = arms*(dcos(wbw *par%t) + 0.5d0*kbw*arms*dcos(2*wbw*par%t) + 0.375d0*(kbw*arms)**2*dcos(3*wbw*par%t)) 00699 !s%ui(1,:) = (wbw/kbw)/s%hh(1,:)*s%zi(1,:) 00700 endif 00701 ! 2) spectral 00702 elseif ((par%wbctype==WBCTYPE_PARAMETRIC) .or. (par%wbctype==WBCTYPE_JONS_TABLE) .or. & 00703 (par%wbctype==WBCTYPE_SWAN) .or. (par%wbctype==WBCTYPE_VARDENS) .or. (par%wbctype==WBCTYPE_REUSE)) then 00704 if (startbcf) then 00705 if(xmaster) then 00706 open(53,file='nhbcflist.bcf',form='formatted',position='rewind') 00707 do i=1,curline 00708 read(53,*,iostat=ier)bcstarttime,bcendtime,par%Trep,nhbcfname 00709 if (ier .ne. 0) then 00710 call report_file_read_error('nhbcflist.bcf') 00711 endif 00712 enddo 00713 close(53) 00714 endif 00715 if (.not. allocated(uig)) then 00716 if (xmaster) then 00717 allocate(uig(sg%ny+1)) 00718 allocate(vig(sg%ny+1)) 00719 allocate(zig(sg%ny+1)) 00720 allocate(wig(sg%ny+1)) 00721 allocate(duig(sg%ny+1)) 00722 allocate(dvig(sg%ny+1)) 00723 else ! only need valid address for MPI-distribution 00724 allocate(uig(1)) 00725 allocate(vig(1)) 00726 allocate(zig(1)) 00727 allocate(wig(1)) 00728 allocate(duig(1)) 00729 allocate(dvig(1)) 00730 endif 00731 endif 00732 if (xmaster) then 00733 ! Robert: optimise? 00734 if (trim(nhbcfname)=='nh_reuse.bcf') then 00735 ! Reuse time series will start at t=0, so we need to add the offset time to 00736 ! the time vector in the boundary condition file using the optional bcst variable 00737 call velocity_Boundary(nhbcfname,uig,vig,zig,wig,duig,dvig,sg%ny,par%t, & 00738 isSet_U,isSet_V,isSet_Z,isSet_W,isSet_dU,isSet_dV,isSet_Q, & 00739 force_init=.true.,bcst=bcstarttime) 00740 else 00741 ! If individual spectra are read, then the time vector in the boundary condition 00742 ! file will already be corrected for the start of the spectrum condition. No need 00743 ! to add the optional bcst variable 00744 call velocity_Boundary(nhbcfname,uig,vig,zig,wig,duig,dvig,sg%ny,par%t, & 00745 isSet_U,isSet_V,isSet_Z,isSet_W,isSet_dU,isSet_dV,isSet_Q, & 00746 force_init=.true.) 00747 endif 00748 endif 00749 #ifdef USEMPI 00750 call space_distribute("y",sl,uig,s%ui(1,:)) 00751 call space_distribute("y",sl,vig,s%vi(1,:)) 00752 call space_distribute("y",sl,zig,s%zi(1,:)) 00753 call space_distribute("y",sl,wig,s%wi(1,:)) 00754 call space_distribute("y",sl,duig,s%dui(1,:)) 00755 call space_distribute("y",sl,dvig,s%dvi(1,:)) 00756 call xmpi_bcast(isSet_U) 00757 call xmpi_bcast(isSet_V) 00758 call xmpi_bcast(isSet_Z) 00759 call xmpi_bcast(isSet_W) 00760 call xmpi_bcast(isSet_dU) 00761 call xmpi_bcast(isSet_dV) 00762 call xmpi_bcast(bcendtime) 00763 #else 00764 s%ui(1,:) = uig 00765 s%vi(1,:) = vig 00766 s%zi(1,:) = zig 00767 s%wi(1,:) = wig 00768 s%dui(1,:) = duig 00769 s%dvi(1,:) = dvig 00770 #endif 00771 if (.not. isSet_U) s%ui(1,:) = 0.d0 00772 if (.not. isSet_Z) s%zi(1,:) = s%zs(2,:) 00773 if (.not. isSet_W) s%wi(1,:) = s%ws(2,:) 00774 if (.not. isSet_V) then 00775 s%vi(1,:) = 0.d0 00776 s%isSet_Vbc = 0 00777 else 00778 s%isSet_Vbc = 1 00779 endif 00780 if (par%nonhq3d == 1 .and. par%nhlay > 0.0d0) then 00781 ! Only if the reduced 2-layer model is initiated (P.B. Smit) 00782 if (.not. isSet_dU) s%dUi(1,:) = 0.d0 00783 if (.not. isSet_dV) s%dVi(1,:) = 0.d0 00784 endif 00785 call nonh_init_wcoef(s,par) 00786 else 00787 if (xmaster) then 00788 if (trim(nhbcfname)=='nh_reuse.bcf') then 00789 ! Reuse time series will start at t=0, so we need to add the offset time to 00790 ! the time vector in the boundary condition file using the optional bcst variable 00791 call velocity_Boundary(nhbcfname,uig,vig,zig,wig,duig,dvig,sg%ny,par%t,& 00792 isSet_U,isSet_V,isSet_Z,isSet_W,isSet_dU,isSet_dV,isSet_Q, & 00793 bcst=bcstarttime) 00794 else 00795 ! If individual spectra are read, then the time vector in the boundary condition 00796 ! file will already be corrected for the start of the spectrum condition. No need 00797 ! to add the optional bcst variable 00798 call velocity_Boundary(nhbcfname,uig,vig,zig,wig,duig,dvig,sg%ny,par%t, & 00799 isSet_U,isSet_V,isSet_Z,isSet_W,isSet_dU,isSet_dV,isSet_Q) 00800 endif 00801 endif 00802 #ifdef USEMPI 00803 call space_distribute("y",sl,uig,s%ui(1,:)) 00804 call space_distribute("y",sl,vig,s%vi(1,:)) 00805 call space_distribute("y",sl,zig,s%zi(1,:)) 00806 call space_distribute("y",sl,wig,s%wi(1,:)) 00807 call space_distribute("y",sl,duig,s%dui(1,:)) 00808 call space_distribute("y",sl,dvig,s%dvi(1,:)) 00809 call xmpi_bcast(isSet_U) 00810 call xmpi_bcast(isSet_V) 00811 call xmpi_bcast(isSet_Z) 00812 call xmpi_bcast(isSet_W) 00813 call xmpi_bcast(isSet_dU) 00814 call xmpi_bcast(isSet_dV) 00815 #else 00816 s%dui(1,:) = duig 00817 s%dvi(1,:) = dvig 00818 s%ui(1,:) = uig 00819 s%vi(1,:) = vig 00820 s%zi(1,:) = zig 00821 s%wi(1,:) = wig 00822 #endif 00823 if (.not. isSet_U) s%ui(1,:) = 0.d0 00824 if (.not. isSet_Z) s%zi(1,:) = s%zs(2,:) 00825 if (.not. isSet_W) s%wi(1,:) = s%ws(2,:) 00826 if (.not. isSet_V) then 00827 s%vi(1,:) = 0.d0 00828 s%isSet_Vbc = 0 00829 else 00830 s%isSet_Vbc = 1 00831 endif 00832 if (par%nonhq3d == 1 .and. par%nhlay > 0.0d0) then 00833 ! Only if the reduced 2-layer model is initiated (P.B. Smit) 00834 if (.not. isSet_dU) s%dUi(1,:) = 0.d0 00835 if (.not. isSet_dV) s%dVi(1,:) = 0.d0 00836 ! 00837 endif 00838 endif 00839 ! be aware: ui and vi are defined w.r.t. the grid, not w.r.t. the coordinate system! 00840 if (.not. allocated(tempu)) then 00841 ! allocate to local size grid 00842 allocate(tempu(s%ny+1)) 00843 allocate(tempv(s%ny+1)) 00844 allocate(tempdv(s%ny+1)) 00845 allocate(tempdu(s%ny+1)) 00846 endif 00847 tempu = s%ui(1,:) 00848 tempv = s%vi(1,:) 00849 tempdu = s%dUi(1,:) 00850 tempdv = s%dVi(1,:) 00851 ! rotate to grid directions 00852 s%ui(1,:) = tempu*dcos(s%alfau(1,:)) + tempv*dsin(s%alfau(1,:)) 00853 s%vi(1,:) = -tempu*dsin(s%alfau(1,:)) + tempv*dcos(s%alfau(1,:)) 00854 s%dUi(1,:) = tempdu*dcos(s%alfau(1,:)) + tempdv*dsin(s%alfau(1,:)) 00855 00856 s%ui(1,:) = s%ui(1,:)*min(par%t/par%taper,1.0d0) 00857 s%zi(1,:) = s%zi(1,:)*min(par%t/par%taper,1.0d0) 00858 s%wi(1,:) = s%wi(1,:)*min(par%t/par%taper,1.0d0) 00859 s%dUi(1,:) = s%dUi(1,:)*min(par%t/par%taper,1.0d0) 00860 ! 3) ts_nonh 00861 elseif (par%wbctype==WBCTYPE_TS_NONH) then 00862 if (xmaster) then 00863 call velocity_Boundary('boun_U.bcf',uig,vig,zig,wig,duig,dvig,sg%ny,& 00864 par%t,isSet_U,isSet_V,isSet_Z,isSet_W,isSet_dU,isSet_dV,isSet_Q) 00865 endif 00866 #ifdef USEMPI 00867 call space_distribute("y",sl,uig,s%ui(1,:)) 00868 call space_distribute("y",sl,vig,s%vi(1,:)) 00869 call space_distribute("y",sl,zig,s%zi(1,:)) 00870 call space_distribute("y",sl,wig,s%wi(1,:)) 00871 call space_distribute("y",sl,duig,s%dui(1,:)) 00872 call space_distribute("y",sl,dvig,s%dvi(1,:)) 00873 call xmpi_bcast(isSet_U) 00874 call xmpi_bcast(isSet_Z) 00875 call xmpi_bcast(isSet_W) 00876 call xmpi_bcast(isSet_dU) 00877 call xmpi_bcast(isSet_dV) 00878 #else 00879 if ( isSet_Q ) then 00880 if (par%nonhq3d == 1 .and. par%nhlay > 0.0d0) then 00881 ! Calculate the difference velocity from the difference discharges (here both duig and uig are still discharges) 00882 duig = ( (1.0d0 - 2.0d0*par%nhlay) * uig + duig ) / ( 2.0d0 * par%nhlay * s%hu(1,:) *( 1.0d0 - par%nhlay ) ) 00883 endif 00884 uig = uig/s%hu(1,:) 00885 endif 00886 s%ui(1,:) = uig 00887 s%vi(1,:) = vig 00888 s%zi(1,:) = zig 00889 s%wi(1,:) = wig 00890 s%dUi(1,:) = duig 00891 s%dVi(1,:) = dvig 00892 #endif 00893 if (.not. isSet_U) s%ui(1,:) = 0.d0 00894 if (.not. isSet_Z) s%zi(1,:) = s%zs(2,:) 00895 if (.not. isSet_W) s%wi(1,:) = s%ws(2,:) 00896 if (.not. isSet_V) then 00897 s%vi(1,:) = 0.d0 00898 s%isSet_Vbc = 0 00899 else 00900 s%isSet_Vbc = 1 00901 endif 00902 if (par%nonhq3d == 1 .and. par%nhlay > 0.0d0) then 00903 ! Only if the reduced 2-layer model is initiated (P.B. Smit) 00904 if (.not. isSet_dU) s%dUi(1,:) = 0.d0 00905 if (.not. isSet_dV) s%dVi(1,:) = 0.d0 00906 endif 00907 else 00908 endif ! Go over the different wbctypes for nonh 00909 endif ! Go over the different wavemodels 00910 ! 00911 ! Finalize: swave and lwave = 0 has influence on boundary 00912 s%ee(1,:,:) = par%swave*s%ee(1,:,:) 00913 if (par%nonhspectrum==0) then 00914 s%ui = par%lwave*(par%order-1.d0)*s%ui 00915 endif 00916 end subroutine wave_bc 00917 00918 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00919 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00920 ! 00921 ! FLOW BOUNDARY CONDITIONS 00922 ! 00923 subroutine flow_bc(s,par) 00924 use params, only: parameters 00925 use spaceparams 00926 use interp, only: linear_interp 00927 use compute_tide_module, only: tide_boundary_timestep 00928 #ifdef USEMPI 00929 use xmpi_module 00930 #endif 00931 use paramsconst 00932 00933 implicit none 00934 00935 type(spacepars), target :: s 00936 type(parameters) :: par 00937 00938 integer :: i,j,jj,j1,indt 00939 real*8 :: alphanew,vert,windxnow,windynow,factime 00940 real*8 :: udnusum,dnusum,vdnvsum,dnvsum 00941 !real*8 , dimension(2) :: xzs0,yzs0,szs0 00942 real*8 , dimension(:,:) ,allocatable,save :: zs0old,dzs0 00943 real*8 , dimension(:,:) ,allocatable,save :: zsmean 00944 real*8 , dimension(:,:) ,allocatable,save :: ht,beta,betanp1 00945 real*8 , dimension(:) ,allocatable,save :: bn,alpha2,thetai 00946 real*8 , dimension(:,:) ,allocatable,save :: dhdx,dhdy,dvdx,dvdy 00947 real*8 , dimension(:,:) ,allocatable,save :: dbetadx,dbetady,dvudy 00948 real*8 , dimension(:) ,allocatable,save :: inv_ht 00949 00950 00951 if (.not. allocated(ht)) then 00952 allocate(ht (2,s%ny+1)) 00953 allocate(dhdx (2,s%ny+1)) 00954 allocate(dhdy (2,s%ny+1)) ! wwvv not used 00955 allocate(dvdx (2,s%ny+1)) ! wwvv not used 00956 allocate(dvdy (2,s%ny+1)) 00957 allocate(dvudy (1,s%ny+1)) ! wwvv not used 00958 allocate(inv_ht (s%ny+1)) 00959 allocate(dbetadx (2,s%ny+1)) 00960 allocate(dbetady (2,s%ny+1)) 00961 allocate(beta (2,s%ny+1)) 00962 allocate(zsmean (2,s%ny+1)) 00963 allocate(bn (s%ny+1)) 00964 allocate(alpha2 (s%ny+1)) 00965 allocate(thetai (s%ny+1)) 00966 allocate(betanp1 (1,s%ny+1)) 00967 if (par%tidetype==TIDETYPE_INSTANT) then 00968 allocate(zs0old(s%nx+1,s%ny+1)) 00969 allocate(dzs0 (s%nx+1,s%ny+1)) 00970 zs0old = s%zs0 00971 dzs0 = 0.d0 00972 endif 00973 ! initialize zsmean 00974 zsmean(1,:) = s%zs(1,:) 00975 zsmean(2,:) = s%zs(s%nx,:) 00976 if (par%hotstart==0) then 00977 ! ToDo: also complete this for zsmean (would need to be added to s% structure) 00978 s%umean = s%uu 00979 s%vmean = s%vv 00980 endif 00981 thetai = 0.d0 00982 endif 00983 00984 ! Super fast 1D 00985 if (s%ny==0) then 00986 j1 = 1 00987 else 00988 j1 = 2 00989 endif 00990 00991 ! factime=1.d0/par%cats/par%Trep*par%dt !Jaap: not used anymore 00992 ! jaap: compute epsi based on offshore boundary conditions if epsi = -1 00993 if (par%epsi==-1.d0) then 00994 if (par%swave==1) then 00995 factime = 1.d0/par%cats/par%Trep*par%dt 00996 else 00997 ! need to do something here 00998 factime = 1.d0/60*par%dt 00999 endif 01000 else 01001 factime = par%epsi 01002 endif 01003 01004 !allocate(szs0(1:2)) ! wwvv changed this, now defined as an automatic array 01005 !allocate(yzs0(1:2)) 01006 ! 01007 ! Sea boundary at x=0; 01008 ! 01009 ! 01010 ! UPDATE TIDE AND SURGE 01011 if (par%tidetype==TIDETYPE_INSTANT) then 01012 zs0old = s%zs0 ! store zs0 of old time step 01013 call tide_boundary_timestep(s,par) ! also updates zs0 across domain 01014 where (s%wetz==1) 01015 dzs0 = s%zs0-zs0old 01016 s%zs = s%zs+dzs0 01017 elsewhere 01018 dzs0 = 0.d0 01019 endwhere 01020 else 01021 call tide_boundary_timestep(s,par) ! only updates zs0 on boundaries 01022 endif 01023 ! 01024 !Dano: compute umean, vmean only at this location, same for all options; MPI compliant 01025 if (par%tidetype==TIDETYPE_VELOCITY) then 01026 if (xmpi_istop) then 01027 ! make sure we have same umean along whole offshore boundary 01028 udnusum = sum(s%uu (1,jmin_uu:jmax_uu)*s%dnu(1,jmin_uu:jmax_uu)) 01029 dnusum = sum(s%dnu(1,jmin_uu:jmax_uu)) 01030 vdnvsum = sum(s%vv (1,jmin_vv:jmax_vv)*s%dnv(1,jmin_vv:jmax_vv)) 01031 dnvsum = sum(s%dnv(1,jmin_vv:jmax_vv)) 01032 else 01033 udnusum = 0.d0 01034 dnusum = 0.d0 01035 vdnvsum = 0.d0 01036 dnvsum = 0.d0 01037 endif 01038 #ifdef USEMPI 01039 call xmpi_allreduce(udnusum ,MPI_SUM) 01040 call xmpi_allreduce(dnusum ,MPI_SUM) 01041 call xmpi_allreduce(vdnvsum ,MPI_SUM) 01042 call xmpi_allreduce(dnvsum ,MPI_SUM) 01043 #endif 01044 s%umean(1,:)=factime*udnusum/dnusum+(1.0d0-factime)*s%umean(1,:) 01045 s%vmean(1,:)=factime*vdnvsum/dnvsum+(1.0d0-factime)*s%vmean(1,:) 01046 01047 else 01048 s%umean(1,:) = 0.d0 01049 s%vmean(1,:) = 0.d0 01050 endif 01051 ! 01052 ! UPDATE (LONG) WAVES 01053 ! 01054 #ifdef USEMPI 01055 call xmpi_shift(s%ui,SHIFT_Y_R,1,2) 01056 call xmpi_shift(s%vi,SHIFT_Y_R,1,2) 01057 call xmpi_shift(s%zi,SHIFT_Y_R,1,2) 01058 call xmpi_shift(s%wi,SHIFT_Y_R,1,2) 01059 01060 call xmpi_shift(s%ui,SHIFT_Y_L,3,4) 01061 call xmpi_shift(s%vi,SHIFT_Y_L,3,4) 01062 call xmpi_shift(s%zi,SHIFT_Y_L,3,4) 01063 call xmpi_shift(s%wi,SHIFT_Y_L,3,4) 01064 #endif 01065 if (par%wbctype/=WBCTYPE_OFF)then 01066 ! wwvv the following is probably only to do in the top processes, but take care for 01067 ! the mpi_shift calls in horizontal directions 01068 if(xmpi_istop) then 01069 if (par%front==FRONT_ABS_1D) then ! Ad's radiating boundary 01070 !ht(1:2,:)=max(zs0(1:2,:)-zb(1:2,:),par%eps) 01071 !uu(1,:)=2.d0*ui(1,:)-(sqrt(par%g/hh(1,:))*(zs(2,:)-zs0(2,:))) 01072 if (par%freewave==1) then 01073 s%uu(1,jmin_zs:jmax_zs)=2.0d0*s%ui(1,jmin_zs:jmax_zs)-(sqrt(par%g/s%hh(1,jmin_zs:jmax_zs))*(s%zs(2,jmin_zs:jmax_zs)-s%zs0(2,jmin_zs:jmax_zs))) + s%umean(1,jmin_zs:jmax_zs) 01074 else 01075 s%uu(1,jmin_zs:jmax_zs)=(1.0d0+sqrt(par%g*s%hh(1,jmin_zs:jmax_zs))/s%cg(1,jmin_zs:jmax_zs))*s%ui(1,jmin_zs:jmax_zs) - & 01076 (sqrt(par%g/s%hh(1,jmin_zs:jmax_zs))*(s%zs(2,jmin_zs:jmax_zs)-s%zs0(2,jmin_zs:jmax_zs))) + s%umean(1,jmin_zs:jmax_zs) 01077 endif 01078 s%vv(1,jmin_zs:jmax_zs)=s%vv(2,jmin_zs:jmax_zs) 01079 s%zs(1,jmin_zs:jmax_zs)=s%zs(2,jmin_zs:jmax_zs) 01080 if (par%wavemodel==WAVEMODEL_NONH) s%ws(1,jmin_zs:jmax_zs) = s%ws(2,jmin_zs:jmax_zs) 01081 elseif (par%front==FRONT_WAVEFLUME) then ! Ad's radiating boundary 01082 if (par%tidetype==TIDETYPE_VELOCITY) then 01083 zsmean(1,jmin_zs:jmax_zs) = factime*s%zs(1,jmin_zs:jmax_zs)+(1-factime)*zsmean(1,jmin_zs:jmax_zs) 01084 endif 01085 s%uu(1,jmin_zs:jmax_zs)=(1.0d0+sqrt(par%g*s%hh(1,jmin_zs:jmax_zs))/s%cg(1,jmin_zs:jmax_zs))*s%ui(1,jmin_zs:jmax_zs)-(sqrt(par%g/s%hh(1,jmin_zs:jmax_zs))*(s%zs(2,jmin_zs:jmax_zs)-zsmean(1,jmin_zs:jmax_zs))) 01086 s%vv(1,jmin_zs:jmax_zs)=s%vv(2,jmin_zs:jmax_zs) 01087 s%zs(1,jmin_zs:jmax_zs)=s%zs(2,jmin_zs:jmax_zs) 01088 if (par%wavemodel==WAVEMODEL_NONH) s%ws(1,jmin_zs:jmax_zs) = s%ws(2,jmin_zs:jmax_zs) 01089 elseif (par%front==FRONT_ABS_2D) then ! Van Dongeren (1997), weakly reflective boundary condition 01090 ht(1:2,:)=max(s%zs0(1:2,:)-s%zb(1:2,:),par%eps) 01091 beta=s%uu(1:2,:)-2.*dsqrt(par%g*s%hum(1:2,:)) 01092 01093 !do j=j1,max(s%ny,1) 01094 do j=jmin_zs,jmax_zs 01095 ! compute gradients in u-points 01096 dhdx (1,j) = ( ht (2,j) - ht (1,j) ) / s%dsu(1,j) 01097 dbetadx(1,j) = ( beta(2,j) - beta(1,j) ) / s%dsz(2,j) 01098 01099 if (s%ny>0) then 01100 dvdy (1,j) = ( s%vu (1,j+1) - s%vu (1,j-1) ) / ( s%dnu(1,j)+(s%dnu(1,j-1)+s%dnu(1,j+1))/2 ) 01101 dbetady(1,j) = ( beta(1,j+1) - beta(1,j-1) ) / ( s%dnu(1,j)+(s%dnu(1,j-1)+s%dnu(1,j+1))/2 ) 01102 else 01103 dvdy (1,j) = 0.d0 01104 dbetady(1,j) = 0.d0 01105 endif 01106 01107 inv_ht(j) = 1.d0/s%hum(1,j) 01108 01109 bn (j) = -( s%uu(1,j) - dsqrt(par%g*s%hum(1,j)) ) * dbetadx(1,j) & 01110 - s%vu(1,j) * dbetady(1,j) & 01111 + dsqrt(par%g*s%hum(1,j)) * dvdy (1,j) & 01112 + par%g * dhdx (1,j) & 01113 + s%Fx(1,j) * inv_ht(j) / par%rho & 01114 - s%cfu(1,j) * sqrt(s%uu(1,j)**2 + s%vu(1,j)**2) * s%uu(1,j) / s%hum(1,j) 01115 end do 01116 01117 ! Jaap: Compute angle of incominge wave 01118 thetai = datan(s%vi(1,:)/(s%ui(1,:)+1.d-16)) 01119 01120 !do j=j1,max(s%ny,1) 01121 do j=jmin_zs,jmax_zs 01122 01123 betanp1(1,j) = beta(1,j)+ bn(j)*par%dt 01124 alpha2(j)=-(s%theta0-s%alfaz(1,j)) ! Jaap: this is first estimate 01125 alphanew = 0.d0 01126 01127 do jj=1,50 01128 !---------- Lower order bound. cond. --- 01129 01130 if (par%freewave==1) then ! assuming incoming long wave propagates at sqrt(g*h) 01131 s%ur(1,j) = dcos(alpha2(j))/(dcos(alpha2(j))+1.d0)& 01132 *(betanp1(1,j)-s%umean(1,j)+2.d0*DSQRT(par%g*0.5d0*(ht(1,j)+ht(2,j))) & 01133 -s%ui(1,j)*(dcos(thetai(j))-1.d0)/dcos(thetai(j))) 01134 else ! assuming incoming long wave propagates at s%cg 01135 s%ur(1,j) = dcos(alpha2(j))/(dcos(alpha2(j))+1.d0)& 01136 *(betanp1(1,j)-s%umean(1,j)+2.d0*DSQRT(par%g*0.5d0*(ht(1,j)+ht(2,j)))& 01137 -s%ui(1,j)*(s%cg(1,j)*dcos(thetai(j))- DSQRT(par%g*0.5d0*(ht(1,j)+ht(2,j))) )/& 01138 (s%cg(1,j)*dcos(thetai(j)))) 01139 endif 01140 01141 !vert = velocity of the reflected wave = total-specified 01142 vert = s%vu(1,j)-s%vmean(1,j)-s%vi(1,j) ! Jaap use s%vi instead of s%ui(1,j)*tan(s%theta0) 01143 alphanew = datan(vert/(s%ur(1,j)+1.d-16)) ! wwvv can use atan2 here 01144 if (alphanew .gt. (par%px*0.5d0)) alphanew=alphanew-par%px 01145 if (alphanew .le. (-par%px*0.5d0)) alphanew=alphanew+par%px 01146 if(dabs(alphanew-alpha2(j)).lt.0.001d0) goto 1000 ! wwvv can use exit here 01147 alpha2(j) = alphanew 01148 end do 01149 1000 continue 01150 if (par%ARC==0) then 01151 s%uu(1,j) = (par%order-1)*s%ui(1,j) + s%umean(1,j) 01152 s%zs(1,j) = s%zs(2,j) 01153 else 01154 ! 01155 s%uu(1,j) = (par%order-1.d0)*s%ui(1,j) + s%ur(1,j) + s%umean(1,j) 01156 ! zs(1,:) is dummy variable 01157 ! with a taylor expansion to get to the zs point at index 1 from uu(1) and uu(2) 01158 s%zs(1,j) = 1.5d0*((betanp1(1,j)-s%uu(1,j))**2/4.d0/par%g+.5d0*(s%zb(1,j)+s%zb(2,j)))- & 01159 0.5d0*((beta(2,j) -s%uu(2,j))**2/4.d0/par%g+.5d0*(s%zb(2,j)+s%zb(3,j))) 01160 ! Ad + Jaap: zs does indeed influence hydrodynamics at boundary --> do higher order taylor expansions to check influence 01161 ! zs(1,j) = 13.d0/8.d0*((betanp1(1,j)-uu(1,j))**2.d0/4.d0/par%g+.5d0*(zb(1,j)+zb(2,j))) - & 01162 ! 0.75d0*((beta(2,j)-uu(2,j))**2.d0/4.d0/par%g+.5d0*(zb(2,j)+zb(3,j))) + & 01163 ! 0.125d0*0.5d0*(zs(3,j)+zs(4,j)) 01164 end if 01165 end do 01166 01167 if (par%wavemodel==WAVEMODEL_NONH.and.xmpi_istop) s%ws(1,jmin_zs:jmax_zs) = s%ws(2,jmin_zs:jmax_zs) 01168 s%vv(1,jmin_zs:jmax_zs)=s%vv(2,jmin_zs:jmax_zs) 01169 else if (par%front==FRONT_WALL) then 01170 ! uu(1,:)=0.d0 01171 ! zs(1,:)=max(zs(2,:),zb(1,:)) 01172 else if (par%front==FRONT_NONH_1D) then 01173 !Timeseries Boundary for nonh, only use in combination with instat=8 01174 if (par%ARC==0) then 01175 s%uu(1,jmin_zs:jmax_zs) = s%ui(1,jmin_zs:jmax_zs) 01176 s%vv(1,jmin_zs:jmax_zs) = s%vi(1,jmin_zs:jmax_zs) 01177 s%zs(1,jmin_zs:jmax_zs) = s%zi(1,jmin_zs:jmax_zs)+s%zs0(1,jmin_zs:jmax_zs) 01178 s%ws(1,jmin_zs:jmax_zs) = s%wi(1,jmin_zs:jmax_zs) 01179 if (par%nonhq3d == 1 .and. par%nhlay > 0.0d0) then 01180 ! Only if the reduced 2-layer model is initiated (P.B. Smit) 01181 s%dU(1,jmin_zs:jmax_zs) = s%dui(1,jmin_zs:jmax_zs) 01182 ! 01183 endif 01184 ! 01185 else 01186 !Radiating boundary for short waves: 01187 s%uu(1,jmin_zs:jmax_zs) = s%ui(1,jmin_zs:jmax_zs)-sqrt(par%g/s%hh(1,jmin_zs:jmax_zs))* & 01188 (s%zs(2,jmin_zs:jmax_zs)-s%zi(1,jmin_zs:jmax_zs)-s%zs0(2,jmin_zs:jmax_zs)) + & 01189 s%umean(1,jmin_zs:jmax_zs) 01190 ! 01191 ! uu(1,:)= 2*ui(1,:)-sqrt(par%g/hh(1,:))*(zs(2,:) -zs0(2,:)) 01192 if (par%nonhq3d == 1 .and. par%nhlay > 0.0d0) then 01193 ! 01194 01195 s%dU(1,jmin_zs:jmax_zs) = s%dui(1,jmin_zs:jmax_zs) 01196 ! write(*,*) s%dU(1,%) 01197 ! 01198 endif 01199 s%vv(1,jmin_zs:jmax_zs) = s%vi(1,jmin_zs:jmax_zs) 01200 s%zs(1,jmin_zs:jmax_zs) = s%zs(2,jmin_zs:jmax_zs) 01201 s%ws(1,jmin_zs:jmax_zs) = s%ws(2,jmin_zs:jmax_zs) 01202 endif 01203 if (s%isSet_Vbc == 0) then 01204 s%vv(1,jmin_zs:jmax_zs)=s%vv(2,jmin_zs:jmax_zs) 01205 endif 01206 endif ! par%front 01207 01208 if (s%ny>0) then 01209 if (xmpi_isleft) then 01210 s%uu(1,1)=s%uu(1,2) 01211 s%vv(1,1)=s%vv(1,2) 01212 s%zs(1,1)=max(s%zs(1,2) - s%dzs0dn(1,1)*s%dnv(1,1),s%zb(1,1)) 01213 endif 01214 if (xmpi_isright) then 01215 s%uu(1,s%ny+1)=s%uu(1,s%ny) 01216 s%vv(1,s%ny+1)=s%vv(1,s%ny) 01217 s%zs(1,s%ny+1)=max(s%zs(1,s%ny) + s%dzs0dn(1,s%ny+1)*s%dnv(1,s%ny),s%zb(1,s%ny+1)) 01218 endif 01219 endif 01220 endif ! xmpi_istop 01221 endif ! par%wbctype/=WBCTYPE_OFF 01222 ! 01223 ! Radiating boundary at x=nx*dx 01224 ! 01225 !Dano: compute umean, vmean only at this location, same for all options; MPI compliant 01226 if (par%tidetype==TIDETYPE_VELOCITY) then 01227 ! make sure we have same umean along whole offshore boundary 01228 if (xmpi_isbot) then 01229 udnusum = sum(s%uu (s%nx,jmin_uu:jmax_uu)*s%dnu(s%nx,jmin_uu:jmax_uu)) 01230 dnusum = sum(s%dnu(s%nx,jmin_uu:jmax_uu)) 01231 vdnvsum = sum(s%vv (s%nx,jmin_vv:jmax_vv)*s%dnv(s%nx,jmin_vv:jmax_vv)) 01232 dnvsum = sum(s%dnv(s%nx,jmin_vv:jmax_vv)) 01233 else 01234 udnusum = 0.d0 01235 dnusum = 0.d0 01236 vdnvsum = 0.d0 01237 dnvsum = 0.d0 01238 endif 01239 #ifdef USEMPI 01240 call xmpi_allreduce(udnusum ,MPI_SUM) 01241 call xmpi_allreduce(dnusum ,MPI_SUM) 01242 call xmpi_allreduce(vdnvsum ,MPI_SUM) 01243 call xmpi_allreduce(dnvsum ,MPI_SUM) 01244 #endif 01245 s%umean(s%nx,:)=factime*udnusum/dnusum+(1.0d0-factime)*s%umean(s%nx,:) 01246 s%vmean(s%nx,:)=factime*vdnvsum/dnvsum+(1.0d0-factime)*s%vmean(s%nx,:) 01247 01248 else 01249 s%umean(s%nx,:) = 0.d0 01250 s%vmean(s%nx,:) = 0.d0 01251 endif 01252 if (xmpi_isbot) then 01253 if (par%back==BACK_WALL) then 01254 s%uu(s%nx,jmin_zs:jmax_zs) = 0.d0 01255 s%zs(s%nx+1,jmin_zs:jmax_zs) = s%zs(s%nx,jmin_zs:jmax_zs) 01256 ! zs(nx+1,2:ny) = zs(nx+1,2:ny) + par%dt*hh(nx,2:ny)*uu(nx,2:ny)/(xu(nx+1)-xu(nx)) -par%dt*(hv(nx+1,2:ny+1)*vv(nx+1,2:ny+1)-hv(nx+1,1:ny)*vv(nx+1,1:ny))/(yv(2:ny+1)-yv(1:ny)) 01257 elseif (par%back==BACK_ABS_1D) then 01258 s%uu(s%nx,jmin_zs:jmax_zs)=sqrt(par%g/s%hh(s%nx,jmin_zs:jmax_zs))*(s%zs(s%nx,jmin_zs:jmax_zs)-max(s%zb(s%nx,jmin_zs:jmax_zs),s%zs0(s%nx,jmin_zs:jmax_zs)))+s%umean(s%nx,jmin_zs:jmax_zs) ! cjaap: make sure if the last cell is dry no radiating flow is computed... 01259 s%zs(s%nx+1,:)=s%zs(s%nx,:) 01260 elseif (par%back==BACK_ABS_2D) then 01261 ht(1:2,:)=max(s%zs0(s%nx:s%nx+1,:)-s%zb(s%nx:s%nx+1,:),par%eps) !cjaap; make sure ht is always larger than zero 01262 01263 beta=s%uu(s%nx-1:s%nx,:)+2.*dsqrt(par%g*s%hum(s%nx-1:s%nx,:)) !cjaap : replace s%hh with s%hum 01264 01265 do j=jmin_zs,jmax_zs 01266 if (s%wetu(s%nx,j)==1) then ! Robert: dry back boundary points 01267 ! Compute gradients in u-points.... 01268 dvdy(2,j)=(s%vu(s%nx,j+1)-s%vu(s%nx,j-1))/(2.d0*s%dnu(s%nx,j)) 01269 dhdx(2,j)=(ht(2,j)-ht(1,j))/s%dsu(s%nx,j) 01270 dbetadx(2,j)=(beta(2,j)-beta(1,j))/s%dsz(s%nx,j) 01271 dbetady(2,j)=(beta(1,j+1)-beta(1,j-1))/(2.0*s%dnu(s%nx,j)) 01272 01273 inv_ht(j) = 1.d0/s%hum(s%nx,j) 01274 01275 bn(j)=-(s%uu(s%nx,j)+dsqrt(par%g*s%hum(s%nx,j)))*dbetadx(2,j) & 01276 -s%vu(s%nx,j)*dbetady(2,j)& !Ap s%vu 01277 -dsqrt(par%g*s%hum(s%nx,j))*dvdy(2,j)& 01278 +s%Fx(s%nx,j)*inv_ht(j)/par%rho & 01279 -s%cfu(s%nx,j)*sqrt(s%uu(s%nx,j)**2+s%vu(s%nx,j)**2)*s%uu(s%nx,j)/s%hum(s%nx,j)& 01280 +par%g*dhdx(2,j) 01281 endif ! Robert: dry back boundary points 01282 enddo 01283 01284 do j=jmin_zs,jmax_zs 01285 if (s%wetu(s%nx,j)==1) then ! Robert: dry back boundary points 01286 betanp1(1,j) = beta(2,j)+ bn(j)*par%dt !Ap toch? 01287 alpha2(j)= s%theta0 01288 alphanew = 0.d0 01289 01290 do jj=1,50 01291 ! 01292 !---------- Lower order bound. cond. --- 01293 ! 01294 s%ur(s%nx,j) = dcos(alpha2(j))/(dcos(alpha2(j))+1.d0)& 01295 *((betanp1(1,j)-s%umean(s%nx,j)-2.d0*DSQRT(par%g*0.5*(ht(1,j)+ht(2,j))))) 01296 01297 !vert = velocity of the reflected wave = total-specified 01298 vert = s%vu(s%nx,j)-s%vmean(s%nx,j) !Jaap included s%vmean 01299 alphanew = datan(vert/(s%ur(s%nx,j)+1.d-16)) ! wwvv maybe better atan2 01300 if (alphanew .gt. (par%px*0.5d0)) alphanew=alphanew-par%px 01301 if (alphanew .le. (-par%px*0.5d0)) alphanew=alphanew+par%px 01302 if(dabs(alphanew-alpha2(j)).lt.0.001) exit ! wwvv can use exit here 01303 alpha2(j) = alphanew 01304 end do 01305 01306 s%uu(s%nx,j) = s%ur(s%nx,j) + s%umean(s%nx,j) 01307 ! Ap replaced zs with extrapolation. 01308 s%zs(s%nx+1,j) = 1.5*((betanp1(1,j)-s%uu(s%nx,j))**2.d0/4.d0/par%g+.5*(s%zb(s%nx,j)+s%zb(s%nx+1,j)))-& 01309 0.5*((beta(1,j)-s%uu(s%nx-1,j))**2.d0/4.d0/par%g+.5*(s%zb(s%nx-1,j)+s%zb(s%nx,j))) 01310 01311 endif ! Robert: dry back boundary points 01312 enddo 01313 endif !par%back 01314 s%vv(s%nx+1,jmin_zs:jmax_zs)=s%vv(s%nx,jmin_zs:jmax_zs) 01315 s%uu(s%nx+1,jmin_zs:jmax_zs)=s%uu(s%nx,jmin_zs:jmax_zs) 01316 if (s%ny>0) then 01317 if (xmpi_isleft) then 01318 s%uu(s%nx+1,1)=s%uu(s%nx+1,2) 01319 s%vv(s%nx+1,1)=s%vv(s%nx+1,2) 01320 s%zs(s%nx+1,1)=max(s%zs(s%nx+1,2) - s%dzs0dn(1,1)*s%dnv(1,1),s%zb(1,1)) 01321 endif 01322 if (xmpi_isright) then 01323 s%uu(s%nx+1,s%ny+1)=s%uu(s%nx+1,s%ny) 01324 s%vv(s%nx+1,s%ny+1)=s%vv(s%nx+1,s%ny) 01325 s%zs(s%nx+1,s%ny+1)=max(s%zs(s%nx+1,s%ny) + s%dzs0dn(1,s%ny+1)*s%dnv(1,s%ny),s%zb(1,s%ny+1)) 01326 endif 01327 endif 01328 endif !xmpi_isbot 01329 ! 01330 ! Lateral tide boundary conditions 01331 ! 01332 if (s%ny>0) then 01333 do i=imin_zs,imax_zs 01334 if (xmpi_isleft) then 01335 s%zs(i,1)=max(s%zs(i,2) - s%dzs0dn(i,1)*s%dnv(i,1),s%zb(i,1)) 01336 endif 01337 if (xmpi_isright) then 01338 s%zs(i,s%ny+1)=max(s%zs(i,s%ny) + s%dzs0dn(i,s%ny+1)*s%dnv(i,s%ny),s%zb(i,s%ny+1)) 01339 endif 01340 enddo 01341 endif !s%ny>0 01342 01343 01344 01345 !!! Wind boundary conditions 01346 if (s%windlen>1) then ! only if non-stationary wind, otherwise waste of computational time 01347 call LINEAR_INTERP(s%windinpt,s%windxts,s%windlen,par%t,windxnow,indt) 01348 call LINEAR_INTERP(s%windinpt,s%windyts,s%windlen,par%t,windynow,indt) 01349 s%windsu = windxnow*dcos(s%alfau) + windynow*dsin(s%alfau) 01350 ! Can we do this more efficiently? 01351 s%windnv = windynow*dcos(s%alfav-0.5d0*par%px) - windxnow*dsin(s%alfav-0.5d0*par%px) 01352 endif 01353 end subroutine flow_bc 01354 01355 function flow_lat_bc(s,par,bctype,jbc,jn,udvdxb,vdvdyb,viscvb) result(vbc) 01356 use params, only: parameters 01357 use spaceparams 01358 use paramsconst, only: LR_WALL, LR_NEUMANN_V, LR_NO_ADVEC, LR_NEUMANN, LR_ABS_1D 01359 01360 type(spacepars),target :: s 01361 type(parameters) :: par 01362 integer :: bctype 01363 integer :: jbc,jn ! j-index of boundary and j-index of neighbouring cell 01364 ! varies for xmpi_isleft and xmpi_isright 01365 real*8,dimension(s%nx+1) :: udvdxb,vdvdyb,viscvb ! advection terms from flow timestep 01366 real*8,dimension(s%nx+1) :: vbc ! result vector for boundary flow 01367 integer :: i ! internal variables 01368 real*8 :: advterm 01369 real*8,save :: fc 01370 01371 01372 if (par%t<=par%dt) fc =2.d0*par%wearth*sin(par%lat) 01373 01374 if (s%ny==0) then ! 1D models have special case 01375 if (bctype==LR_WALL) then 01376 vbc(:) = 0.d0 ! Only this bc scheme can be applied to s%ny==0 models 01377 else 01378 vbc(:) = s%vv(:,jbc) ! return whatever was already calculated in the flow routine for the cross shore row 01379 endif 01380 else 01381 select case (bctype) 01382 case (LR_WALL) 01383 ! No flow at boundary 01384 vbc(:) = 0.d0 01385 case (LR_NEUMANN_V) 01386 ! Calculate vv at boundary using identical forcing as (:,2), but with dzsdy determined from (:,1), which includes 01387 ! the tide-driven water level gradient as zs(:,1) = zs(:,2) + dzsdytide 01388 ! 01389 ! in flow timestep we calculate: 01390 ! vv(jn) = vvold(jn)-dt(gDZS(jn) + R(jn)) 01391 ! -> vvold(jn) = vv(jn) + dt(gDZS(jn) + R(jn)) 01392 ! 01393 ! we state in bc that: 01394 ! vv(jbc) = vvold(jn)-dt(gDZS(jbc) + R(jn)) 01395 ! -> vv(jbc) = ( vv(jn) + dt(gDZS(jn) + R(jn)) )-dt(gDZS(jbc) + R(jn)) 01396 ! -> vv(jbc) = vv(jn) + dt(gDZS(jn)) - dt(gDZS(jbc)) + dt(R(jn)) - dt(R(jn)) 01397 ! -> vv(jbc) = vv(jn) + dt(gDZS(jn)-gDZS(jbc)) 01398 ! -> vv(jbc) = vv(jn) + dtg(DZS(jn)-DZS(jbc)) 01399 ! vbc(:) = (vv(:,jn)+par%dt*par%g*(dzsdy(i,jn)-dzsdy(i,jbc))) * wetv(:,jbc) 01400 ! AANPASSING: Neumann zonder getij 01401 ! Dano/Jaap: Je kunt niet druk gradient a.g.v getij meenemen zonder ook andere forcerings termen (bijv bodem wrijving) 01402 ! uit te rekenen. In geval van getij gradient heb je een meer complexe neumann rvw die jij "free" noemt 01403 vbc(:) = s%vv(:,jn) 01404 case (LR_NO_ADVEC) 01405 ! We allow vv at the boundary to be calulated from NLSWE, but only include the advective terms if 01406 ! they decrease the flow velocity 01407 do i=2,s%nx 01408 if (s%wetv(i,jbc)==1) then 01409 advterm = udvdxb(i)+vdvdyb(i)-viscvb(i)+ fc*s%uv(i,jbc) 01410 if (s%vv(i,jbc)>0.d0) then 01411 advterm = max(advterm,0.d0) 01412 elseif (s%vv(i,jbc)<0.d0) then 01413 advterm = min(advterm,0.d0) 01414 else 01415 advterm = 0.d0 01416 endif 01417 vbc(i) = s%vv(i,jbc)- par%dt*(advterm & 01418 + par%g*s%dzsdy(i,jbc)& 01419 + s%tauby(i,jbc)/(par%rho*s%hvm(i,jbc)) & 01420 - par%lwave*s%Fy(i,jbc)/(par%rho*s%hvm(i,jbc)) & 01421 - par%rhoa*par%Cd*s%windnv(i,jbc)*sqrt(s%windsu(i,jbc)**2+s%windnv(i,jbc)**2)/(par%rho*s%hvm(i,jbc))) ! Kees: wind correction 01422 ! Robert: no Coriolis because dependecy u in this combination of wall / Neumann 01423 else 01424 vbc(i) = 0.d0 01425 endif 01426 enddo 01427 ! now fill i=1 and i=nx+1 01428 vbc(1) = vbc(2) * s%wetv(1,jbc) 01429 vbc(s%nx+1) = vbc(s%nx) * s%wetv(s%nx+1,jbc) 01430 case (LR_NEUMANN) 01431 ! Dano/Jaap: En dit is dan dus complexe vorm van Neumann (zoals in Delft3D) met getij 01432 ! We allow vv at the boundary to be calulated from NLSWE without any limitations 01433 do i=2,s%nx 01434 if (s%wetv(i,jbc)==1) then 01435 vbc(i) = s%vv(i,jbc)- par%dt*(udvdxb(i)+vdvdyb(i)-viscvb(i)& 01436 + par%g*s%dzsdy(i,jbc)& 01437 + s%tauby(i,jbc)/(par%rho*s%hvm(i,jbc)) & 01438 - par%lwave*s%Fy(i,jbc)/(par%rho*s%hvm(i,jbc)) & 01439 + fc*s%uv(i,jbc) & 01440 - par%rhoa*par%Cd*s%windnv(i,jbc)*sqrt(s%windsu(i,jbc)**2+s%windnv(i,jbc)**2)/(par%rho*s%hvm(i,jbc))) ! Kees: wind correction 01441 else 01442 vbc(i) = 0.d0 01443 endif 01444 enddo 01445 ! now fill i=1 and i=nx+1 01446 vbc(1) = vbc(2) * s%wetv(1,jbc) 01447 vbc(s%nx+1) = vbc(s%nx) * s%wetv(s%nx+1,jbc) 01448 case (LR_ABS_1D) 01449 ! Dano: weakly reflective boundary 01450 do i=2,s%nx 01451 if (s%wetv(i,jbc)==1) then 01452 vbc(i) = sqrt(par%g/s%hh(i,jn))*(s%zs(i,jn)-max(s%zb(i,jn),s%zs0(i,jn)))*(jbc-jn) 01453 s%zs(i,jbc)=s%zs(i,jn) 01454 else 01455 vbc(i) = 0.d0 01456 endif 01457 enddo 01458 ! now fill i=1 and i=nx+1 01459 vbc(1) = vbc(2) * s%wetv(1,jbc) 01460 vbc(s%nx+1) = vbc(s%nx) * s%wetv(s%nx+1,jbc) 01461 end select 01462 endif 01463 end function flow_lat_bc 01464 01465 01466 subroutine discharge_boundary_h(s,par) 01467 use params, only: parameters 01468 use spaceparamsdef 01469 use xmpi_module 01470 use interp, only: linear_interp 01471 01472 01473 implicit none 01474 01475 type(spacepars),target :: s 01476 type(parameters) :: par 01477 01478 integer :: i,indx 01479 integer :: m1,m2,n1,n2 01480 integer :: m1l,m2l,n1l,n2l 01481 real*8 :: qnow,A,CONST 01482 logical :: indomain,isborder 01483 integer :: ishift,jshift 01484 integer :: iminl,imaxl,jminl,jmaxl 01485 01486 #ifdef USEMPI 01487 ishift = s%is(xmpi_rank+1)-1 01488 jshift = s%js(xmpi_rank+1)-1 01489 #else 01490 ishift = 0 01491 jshift = 0 01492 #endif 01493 iminl = imin_zs 01494 imaxl = imax_zs 01495 jminl = jmin_zs 01496 jmaxl = jmax_zs 01497 if (xmpi_istop) then 01498 iminl = imin_zs-1 01499 endif 01500 if (xmpi_isbot) then 01501 imaxl = imax_zs+1 01502 endif 01503 if (xmpi_isleft) then 01504 jminl = jmin_zs-1 01505 endif 01506 if (xmpi_isright) then 01507 jmaxl = jmax_zs+1 01508 endif 01509 01510 ! loop through discharge location 01511 do i = 1,par%ndischarge 01512 if (s%pntdisch(i).eq.0) then 01513 01514 ! interpolate discharge timeseries to current time 01515 call linear_interp(s%tdisch,s%qdisch(:,i),par%ntdischarge,par%t,qnow,indx) 01516 01517 m1 = s%pdisch(i,1) 01518 n1 = s%pdisch(i,2) 01519 m2 = s%pdisch(i,3) 01520 n2 = s%pdisch(i,4) 01521 01522 ! convert to local coordinates 01523 m1l = m1-ishift 01524 m2l = m2-ishift 01525 n1l = n1-jshift 01526 n2l = n2-jshift 01527 01528 indomain = ( m1l >= iminl .and. m1l <= imaxl .and. & 01529 n1l >=jminl .and. n1l <= jmaxl ) .or. & 01530 ( m2l >=iminl .and. m2l <= imaxl .and. & 01531 n2l >=jminl .and. n2l <= jmaxl ) 01532 01533 m1l = min(max(m1l,1),s%nx) 01534 m2l = min(max(m2l,1),s%nx) 01535 n1l = min(max(n1l,1),s%ny) 01536 n2l = min(max(n2l,1),s%ny) 01537 01538 A = 0.d0 01539 01540 ! add momentum in either u- or v-direction 01541 if (n1.eq.n2) then 01542 01543 if (indomain) then 01544 isborder = (n1l.eq.1 .and. xmpi_isleft) .or. (n1l.eq.s%ny .and. xmpi_isright) 01545 m1l=max(m1l,iminl) 01546 m2l=min(m2l,imaxl) 01547 01548 ! make sure discharge is an inflow at border and water levels are defined 01549 if (isborder) then 01550 if (n1l.gt.1) then 01551 qnow = -1*qnow 01552 s%hv(m1l:m2l,n1l) = s%hh(m1l:m2l,n1l) 01553 elseif (n1l.lt.s%ny) then 01554 s%hv(m1l:m2l,n1l) = s%hh(m1l:m2l,n1l+1) 01555 endif 01556 endif 01557 01558 A = sum(s%hv(m1l:m2l,n1l)**1.5*s%dsv(m1l:m2l,n1l)) 01559 endif 01560 01561 #ifdef USEMPI 01562 call xmpi_allreduce(A,MPI_SUM) 01563 #endif 01564 01565 ! compute distribution of discharge over discharge cells according to: 01566 ! vv = CONST * hv^0.5 01567 ! qy = CONST * hv^1.5 01568 ! Q = CONST * hv^1.5 * dsv 01569 01570 if (indomain) then 01571 CONST = qnow/max(A,par%eps) 01572 s%vv(m1l:m2l,n1l) = CONST*sqrt(s%hv(m1l:m2l,n1l)) 01573 s%qy(m1l:m2l,n1l) = CONST*sqrt(s%hv(m1l:m2l,n1l))*s%hv(m1l:m2l,n1l) 01574 endif 01575 01576 elseif (m1.eq.m2) then 01577 01578 if (indomain) then 01579 isborder = (m1l.eq.1 .and. xmpi_istop) .or. (m1l.eq.s%nx .and. xmpi_isbot) 01580 n1l=max(n1l,jminl) 01581 n2l=min(n2l,jmaxl) 01582 01583 ! make sure discharge is an inflow at border and water levels are defined 01584 if (isborder) then 01585 if (m1l.gt.1) then 01586 qnow = -1*qnow 01587 s%hu(m1l,n1l:n2l) = s%hh(m1l,n1l:n2l) 01588 elseif (m1l.lt.s%nx) then 01589 s%hu(m1l,n1l:n2l) = s%hh(m1l+1,n1l:n2l) 01590 endif 01591 endif 01592 01593 A = sum(s%hu(m1l,n1l:n2l)**1.5*s%dnu(m1l,n1l:n2l)) 01594 endif 01595 01596 #ifdef USEMPI 01597 call xmpi_allreduce(A,MPI_SUM) 01598 #endif 01599 01600 ! compute distribution of discharge over discharge cells according to: 01601 ! uu = CONST * hu^0.5 01602 ! qx = CONST * hu^1.5 01603 ! Q = CONST * hu^1.5 * dnu 01604 01605 if (indomain) then 01606 CONST = qnow/max(A,par%eps) 01607 s%uu(m1l,n1l:n2l) = CONST*sqrt(s%hu(m1l,n1l:n2l)) 01608 s%qx(m1l,n1l:n2l) = CONST*sqrt(s%hu(m1l,n1l:n2l))*s%hu(m1l,n1l:n2l) 01609 endif 01610 endif 01611 endif 01612 enddo 01613 end subroutine discharge_boundary_h 01614 01615 subroutine discharge_boundary_v(s,par) 01616 use params, only: parameters 01617 use spaceparamsdef 01618 use interp, only: linear_interp 01619 #ifdef USEMPI 01620 use xmpi_module 01621 #endif 01622 01623 implicit none 01624 01625 type(spacepars),target :: s 01626 type(parameters) :: par 01627 01628 integer :: i,indx 01629 integer :: m1,n1 01630 integer :: m1l,n1l 01631 real*8 :: qnow 01632 integer :: ishift,jshift 01633 01634 01635 #ifdef USEMPI 01636 ishift = s%is(xmpi_rank+1)-1 01637 jshift = s%js(xmpi_rank+1)-1 01638 #else 01639 ishift = 0 01640 jshift = 0 01641 #endif 01642 01643 do i = 1,par%ndischarge 01644 if (s%pntdisch(i).eq.1) then 01645 call linear_interp(s%tdisch,s%qdisch(:,i),par%ntdischarge,par%t,qnow,indx) 01646 01647 m1 = s%pdisch(i,1) 01648 n1 = s%pdisch(i,2) 01649 01650 ! convert to local coordinates 01651 m1l = min(max(m1-ishift,1),s%nx) 01652 n1l = min(max(n1-jshift,1),s%ny) 01653 01654 ! add mass 01655 s%zs(m1l,n1l) = s%zs(m1l,n1l)+qnow*par%dt*s%dsdnzi(m1l,n1l) 01656 endif 01657 enddo 01658 end subroutine discharge_boundary_v 01659 01660 ! 01661 !============================================================================== 01662 subroutine velocity_Boundary(bcfile,ug,vg,zg,wg,dug,dvg,nyg,t,isSet_U,isSet_V,isSet_Z,isSet_W,isSet_dU, & 01663 isSet_dV,isSet_Q,force_init,bcst) 01664 !============================================================================== 01665 ! 01666 01667 !------------------------------------------------------------------------------- 01668 ! DECLARATIONS 01669 !------------------------------------------------------------------------------- 01670 01671 ! 01672 !-------------------------- PURPOSE ---------------------------- 01673 ! 01674 ! Reads boundary values for the velocity from a file. Used as forcing in the 01675 ! nonh module 01676 01677 !-------------------------- METHOD ---------------------------- 01678 01679 01680 !-------------------------- DEPENDENCIES ---------------------------- 01681 01682 use params 01683 use xmpi_module, only: Halt_Program, xmaster 01684 use spaceparamsdef, only: spacepars 01685 use logging_module, only: writelog, report_file_read_error 01686 use readkey_module, only: lowercase 01687 use mnemmodule, only: maxnamelen 01688 01689 implicit none 01690 01691 include 'nh_pars.inc' !Default precision etc. 01692 01693 !-------------------------- ARGUMENTS ---------------------------- 01694 ! 01695 character(len=*) ,intent(in) :: bcfile 01696 integer(kind=iKind) ,intent(in) :: nyg 01697 real(kind=rKind),dimension(nyg+1) :: ug,vg,zg,wg,dug,dvg 01698 real(kind=rKind) ,intent(in) :: t 01699 logical,intent(out) :: isSet_U,isSet_V,isSet_Z,isSet_W,isSet_dU,isSet_dV,isSet_Q 01700 logical,intent(in),optional :: force_init 01701 real(kind=rKind),intent(in),optional :: bcst 01702 ! 01703 01704 !-------------------------- LOCAL VARIABLES ---------------------------- 01705 character(len=iFileNameLen),save :: filename_U = '' 01706 integer ,save :: unit_U = iUnit_U 01707 01708 logical,save :: lVarU = .false. 01709 logical,save :: lIsEof = .false. 01710 logical,save :: lNH_boun_U = .false. 01711 logical,save :: initialize = .true. 01712 logical :: lExists,lOpened 01713 character(slen) :: string 01714 character(len=2),allocatable,dimension(:),save :: header 01715 integer(kind=ikind) :: iAllocErr,ier 01716 integer(kind=ikind),save :: nvar 01717 01718 integer(kind=ikind),save :: iZ = 0 01719 integer(kind=ikind),save :: iU = 0 01720 integer(kind=ikind),save :: iV = 0 01721 integer(kind=ikind),save :: idU = 0 01722 integer(kind=ikind),save :: idV = 0 01723 integer(kind=ikind),save :: iQ = 0 01724 integer(kind=ikind),save :: idQ = 0 01725 integer(kind=ikind),save :: iW = 0 01726 integer(kind=ikind),save :: iT = 0 01727 integer(kind=ikind) :: i 01728 01729 logical :: fi_local 01730 real(kind=rKind) :: bcstarttime 01731 01732 01733 real(kind=rKind),allocatable,dimension(:,:) :: tmp 01734 01735 real(kind=rKind),allocatable,dimension(:),save :: u0 ! 01736 real(kind=rKind),allocatable,dimension(:),save :: u1 ! 01737 01738 real(kind=rKind),allocatable,dimension(:),save :: v0 ! 01739 real(kind=rKind),allocatable,dimension(:),save :: v1 ! 01740 01741 real(kind=rKind),allocatable,dimension(:),save :: z0 ! 01742 real(kind=rKind),allocatable,dimension(:),save :: z1 ! 01743 01744 real(kind=rKind),allocatable,dimension(:),save :: w0 ! 01745 real(kind=rKind),allocatable,dimension(:),save :: w1 ! 01746 01747 real(kind=rKind),allocatable,dimension(:),save :: du0 ! 01748 real(kind=rKind),allocatable,dimension(:),save :: du1 ! 01749 01750 real(kind=rKind),allocatable,dimension(:),save :: dv0 ! 01751 real(kind=rKind),allocatable,dimension(:),save :: dv1 ! 01752 01753 real(kind=rKind),save :: t0 = 0. 01754 real(kind=rKind),save :: t1 = 0. 01755 01756 !------------------------------------------------------------------------------- 01757 ! IMPLEMENTATION 01758 !------------------------------------------------------------------------------- 01759 01760 if (xmaster) then 01761 if (present(bcst)) then 01762 bcstarttime = bcst 01763 else 01764 bcstarttime = 0.d0 01765 endif 01766 01767 if (present(force_init)) then 01768 fi_local = force_init 01769 else 01770 fi_local = .false. 01771 endif 01772 01773 if (bcfile/=filename_U .or. fi_local) then 01774 filename_U=trim(bcfile) 01775 initialize = .true. 01776 endif 01777 ! 01778 ! INITIALIZE NEW FILE 01779 ! 01780 if (initialize) then 01781 initialize = .false. 01782 01783 !Does the file exist and is it already opened (reuse files?) 01784 inquire(file=trim(filename_U),EXIST=lExists,OPENED=lOpened) 01785 01786 ! Close previous file if needed 01787 if(lOpened) then 01788 close(unit_U) 01789 endif 01790 01791 ! Deallocate previous arrays if needed 01792 if (allocated(header)) deallocate(header) 01793 if (allocated(tmp)) deallocate(tmp) 01794 if (allocated(u0)) deallocate(u0) 01795 if (allocated(u1)) deallocate(u1) 01796 if (allocated(v0)) deallocate(v0) 01797 if (allocated(v1)) deallocate(v1) 01798 if (allocated(z0)) deallocate(z0) 01799 if (allocated(z1)) deallocate(z1) 01800 if (allocated(w0)) deallocate(w0) 01801 if (allocated(w1)) deallocate(w1) 01802 if (allocated(du0)) deallocate(du0) 01803 if (allocated(dv0)) deallocate(dv0) 01804 if (allocated(du1)) deallocate(du1) 01805 if (allocated(dv1)) deallocate(dv1) 01806 01807 if (lExists) then 01808 !If exists, open 01809 open(unit_U,file=filename_U,access='SEQUENTIAL',action='READ',form='FORMATTED') 01810 lNH_boun_U = .true. 01811 else 01812 lNH_boun_U = .false. 01813 call writelog('lswe','','Error bc file does not exist: ',trim(filename_U)) 01814 call halt_program 01815 return 01816 endif 01817 01818 01819 !SCALAR OR VECTOR INPUT? 01820 read(unit_U,fmt=*,iostat=ier) string 01821 if (ier .ne. 0) then 01822 call report_file_read_error(filename_U) 01823 endif 01824 string = lowercase(string) 01825 if (string == 'scalar') then 01826 lVaru = .false. 01827 elseif (string == 'vector') then 01828 lVaru = .true. 01829 else 01830 01831 endif 01832 01833 !READ NUMBER OF VARIABLES 01834 read(unit_U,fmt=*,iostat=ier) nvar 01835 if (ier .ne. 0) then 01836 call report_file_read_error(filename_U) 01837 endif 01838 01839 !READ HEADER LINES 01840 allocate(header(nvar)) 01841 read(unit_U,fmt=*,iostat=ier) header 01842 if (ier .ne. 0) then 01843 call report_file_read_error(filename_U) 01844 endif 01845 01846 do i=1,nvar 01847 if (header(i) == 'z' .or. header(i) == 'Z' .or. & 01848 header(i) == 'zs' .or. header(i) == 'ZS' .or. & 01849 header(i) == 'Zs') then 01850 iZ = i; 01851 elseif (header(i) == 't' .or. header(i) == 'T') then 01852 iT = i; 01853 elseif (header(i) == 'u' .or. header(i) == 'U') then 01854 iU = i; 01855 elseif (header(i) == 'v' .or. header(i) == 'V') then 01856 iV = i; 01857 elseif (header(i) == 'w' .or. header(i) == 'W') then 01858 iW = i; 01859 elseif (header(i) == 'dU' .or. header(i) == 'DU') then 01860 idU = i; 01861 elseif (header(i) == 'dV' .or. header(i) == 'DV') then 01862 idV = i; 01863 elseif (header(i) == 'q' .or. header(i) == 'Q') then 01864 iQ = i; 01865 iU = i; 01866 elseif (header(i) == 'dq' .or. header(i) == 'dQ') then 01867 idQ = i; 01868 idU = i 01869 endif 01870 enddo 01871 01872 !Allocate arrays at old and new timelevel 01873 if (iU > 0) then 01874 allocate(u0(nyg+1),stat=iAllocErr); if (iAllocErr /= 0) call Halt_program 01875 allocate(u1(nyg+1),stat=iAllocErr); if (iAllocErr /= 0) call Halt_program 01876 u0 = 0.0d0 01877 u1 = 0.0d0 01878 endif 01879 01880 if (iV > 0) then 01881 allocate(v0(nyg+1),stat=iAllocErr); if (iAllocErr /= 0) call Halt_program 01882 allocate(v1(nyg+1),stat=iAllocErr); if (iAllocErr /= 0) call Halt_program 01883 v0 = 0.0d0 01884 v1 = 0.0d0 01885 endif 01886 01887 if (iZ > 0) then 01888 allocate(z0(nyg+1),stat=iAllocErr); if (iAllocErr /= 0) call Halt_program 01889 allocate(z1(nyg+1),stat=iAllocErr); if (iAllocErr /= 0) call Halt_program 01890 z0 = 0.0d0 01891 z1 = 0.0d0 01892 endif 01893 01894 if (iW > 0) then 01895 allocate(w0(nyg+1),stat=iAllocErr); if (iAllocErr /= 0) call Halt_program 01896 allocate(w1(nyg+1),stat=iAllocErr); if (iAllocErr /= 0) call Halt_program 01897 w0 = 0.0d0 01898 w1 = 0.0d0 01899 endif 01900 01901 if (idU > 0) then 01902 allocate(du0(nyg+1),stat=iAllocErr); if (iAllocErr /= 0) call Halt_program 01903 allocate(du1(nyg+1),stat=iAllocErr); if (iAllocErr /= 0) call Halt_program 01904 du0 = 0.0d0 01905 du1 = 0.0d0 01906 endif 01907 if (idV > 0) then 01908 allocate(dv0(nyg+1),stat=iAllocErr); if (iAllocErr /= 0) call Halt_program 01909 allocate(dv1(nyg+1),stat=iAllocErr); if (iAllocErr /= 0) call Halt_program 01910 dv0 = 0.0d0 01911 dv1 = 0.0d0 01912 endif 01913 01914 allocate(tmp(nyg+1,nvar-1)) 01915 !Read two first timelevels 01916 call velocity_Boundary_read(t0,tmp,Unit_U,lVaru,lIsEof,nvar) 01917 t0 = t0+bcstarttime 01918 if (iU >0) u0 = tmp(:,iU-1) 01919 if (iV >0) v0 = tmp(:,iV-1) 01920 if (iZ >0) z0 = tmp(:,iZ-1) 01921 if (iW >0) w0 = tmp(:,iW-1) 01922 if (idU >0) dU0 = tmp(:,idU-1) 01923 if (idV >0) dV0 = tmp(:,idV-1) 01924 ! 01925 if (lIsEof) then 01926 t1=t0 01927 if (iU >0) u1=u0 01928 if (iV >0) v1=v0 01929 if (iZ >0) z1=z0 01930 if (iW >0) w1=w0 01931 if (idU >0) du1=du0 01932 if (idV >0) dv1=dv0 01933 else 01934 call velocity_Boundary_read(t1,tmp,Unit_U,lVaru,lIsEof,nvar) 01935 t1 = t1+bcstarttime 01936 if (iU >0) u1 = tmp(:,iU-1) 01937 if (iV >0) v1 = tmp(:,iV-1) 01938 if (iZ >0) z1 = tmp(:,iZ-1) 01939 if (iW >0) w1 = tmp(:,iW-1) 01940 if (idu >0) du1 = tmp(:,idU-1) 01941 if (idv >0) dv1 = tmp(:,idV-1) 01942 endif 01943 deallocate(tmp) 01944 if (iU > 0) then 01945 ug = u0 + (u1-u0)*(t-t0)/(t1-t0) 01946 isSet_U = .true. 01947 else 01948 ug = 0 01949 isSet_U = .false. 01950 endif 01951 01952 if (iV > 0) then 01953 vg = v0 + (v1-v0)*(t-t0)/(t1-t0) 01954 isSet_V = .true. 01955 else 01956 vg = 0 01957 isSet_V = .false. 01958 endif 01959 01960 if (iZ > 0) then 01961 zg = z0 + (z1-z0)*(t-t0)/(t1-t0) 01962 isSet_Z = .true. 01963 else 01964 zg = 0 01965 isSet_Z = .false. 01966 endif 01967 if (iW > 0) then 01968 wg = w0 + (w1-w0)*(t-t0)/(t1-t0) 01969 isSet_W = .true. 01970 else 01971 wg = 0 01972 isSet_W = .false. 01973 endif 01974 01975 if (idu > 0) then 01976 dug = dU0 + (dU1-dU0)*(t-t0)/(t1-t0) 01977 isSet_dU = .true. 01978 else 01979 dug = 0 01980 isSet_dU = .false. 01981 endif 01982 01983 if (idv > 0) then 01984 dvg = dV0 + (dV1-dV0)*(t-t0)/(t1-t0) 01985 isSet_dV = .true. 01986 else 01987 dvg = 0 01988 isSet_dV = .false. 01989 endif 01990 01991 if (iQ > 0) then 01992 ! 01993 isSet_Q = .true. 01994 else 01995 isSet_Q = .false. 01996 ! 01997 endif 01998 01999 02000 return 02001 endif ! initialize 02002 ! 02003 ! REGULAR TIMESTEP READ BOUNDARY CONDITIONS 02004 ! 02005 if (lNH_boun_U) then 02006 if (.not. lIsEof) then 02007 allocate(tmp(nyg+1,nvar-1)) 02008 !If current time not located within interval read next line 02009 do while (.not. (t>=t0 .and. t<t1)) 02010 t0 = t1 02011 if (iU >0) u0 = u1 02012 if (iV >0) v0 = v1 02013 if (iZ >0) z0 = z1 02014 if (iW >0) w0 = w1 02015 if (idu >0) du0 = du1 02016 if (idv >0) dv0 = dv1 02017 call velocity_Boundary_read(t1,tmp,Unit_U,lVaru,lIsEof,nvar) 02018 t1 = t1+bcstarttime 02019 if (iU >0) u1 = tmp(:,iU-1) 02020 if (iV >0) v1 = tmp(:,iV-1) 02021 if (iZ >0) z1 = tmp(:,iZ-1) 02022 if (iW >0) w1 = tmp(:,iW-1) 02023 if (idu >0) du1 = tmp(:,idu-1) 02024 if (idv >0) dv1 = tmp(:,idv-1) 02025 if (lIsEof) exit !Exit on end of file condition 02026 end do 02027 deallocate(tmp) 02028 02029 if(lIsEof) then 02030 if (iU >0) ug = u1 02031 if (iV >0) vg = v1 02032 if (iZ >0) zg = z1 02033 if (iW >0) wg = w1 02034 if (idu >0) dug = du1 02035 if (idv >0) dvg = dv1 02036 if (iQ > 0) then 02037 isSet_Q = .true. 02038 else 02039 isSet_Q = .false. 02040 endif 02041 return 02042 else 02043 !Linear interpolation of u in time 02044 if (iU > 0) then 02045 ug = u0 + (u1-u0)*(t-t0)/(t1-t0) 02046 isSet_U = .true. 02047 else 02048 ug = 0 02049 isSet_U = .false. 02050 endif 02051 02052 if (iV > 0) then 02053 vg = v0 + (v1-v0)*(t-t0)/(t1-t0) 02054 isSet_V = .true. 02055 else 02056 vg = 0 02057 isSet_V = .false. 02058 endif 02059 02060 if (iZ > 0) then 02061 zg = z0 + (z1-z0)*(t-t0)/(t1-t0) 02062 isSet_Z = .true. 02063 else 02064 zg = 0 02065 isSet_Z = .false. 02066 endif 02067 if (iW > 0) then 02068 wg = w0 + (w1-w0)*(t-t0)/(t1-t0) 02069 isSet_W = .true. 02070 else 02071 wg = 0 02072 isSet_W = .false. 02073 endif 02074 if (idu > 0) then 02075 dug = dU0 + (dU1-dU0)*(t-t0)/(t1-t0) 02076 isSet_dU = .true. 02077 else 02078 dug = 0 02079 isSet_dU = .false. 02080 endif 02081 if (idv > 0) then 02082 dvg = dV0 + (dV1-dV0)*(t-t0)/(t1-t0) 02083 isSet_dV = .true. 02084 else 02085 dvg = 0 02086 isSet_dV = .false. 02087 endif 02088 ! 02089 if (iQ > 0) then 02090 ! 02091 isSet_Q = .true. 02092 ! 02093 else 02094 isSet_Q = .false. 02095 endif 02096 ! 02097 endif 02098 else 02099 !If end of file the last value which was available is used until the end of the computation 02100 if (iU > 0) ug = u1 02101 if (iV > 0) vg = v1 02102 if (iZ > 0) zg = z1 02103 if (idu > 0) dug = dU1 02104 if (idv > 0) dvg = dV1 02105 if (iW > 0) wg = w1 02106 if (iQ > 0) then 02107 isSet_Q = .true. 02108 else 02109 isSet_Q = .false. 02110 endif 02111 endif ! .not. lIsEof 02112 endif ! lNH_boun_U 02113 endif ! xmaster 02114 end subroutine velocity_Boundary 02115 02116 !============================================================================== 02117 subroutine velocity_Boundary_read(t,vector,iUnit,isvec,iseof,nvar) 02118 !============================================================================== 02119 02120 !-------------------------- PURPOSE ---------------------------- 02121 02122 ! Reads scalar/vector for timeseries 02123 02124 !------------------------------------------------------------------------------- 02125 ! DECLARATIONS 02126 !------------------------------------------------------------------------------- 02127 02128 02129 02130 ! - NONE - 02131 02132 !-------------------------- DEPENDENCIES ---------------------------- 02133 02134 use xmpi_module, only: Halt_Program 02135 use logging_module, only: report_file_read_error 02136 02137 implicit none 02138 02139 include 'nh_pars.inc' !Default precision etc. 02140 !-------------------------- ARGUMENTS ---------------------------- 02141 02142 real(kind=rKind),intent(inout) :: t 02143 real(kind=rKind),intent(inout),dimension(:,:) :: vector 02144 integer(kind=iKind),intent(in) :: iUnit 02145 integer(kind=iKind),intent(in) :: nvar 02146 logical ,intent(in) :: isvec 02147 logical ,intent(out) :: iseof 02148 02149 !-------------------------- LOCAL VARIABLES ---------------------------- 02150 02151 real(kind=rKind) :: scalar(nvar-1) 02152 integer(kind=iKind) :: i 02153 integer(kind=iKind) :: ioStat 02154 character(len=slen) :: ername 02155 02156 !------------------------------------------------------------------------------- 02157 ! IMPLEMENTATION 02158 !------------------------------------------------------------------------------- 02159 02160 02161 iseof = .false. 02162 02163 if (isvec) then 02164 read(iUnit,ioStat=ioStat,fmt=*,end=9000) t,vector 02165 else 02166 read(iUnit,ioStat=ioStat,fmt=*,end=9000) t,scalar 02167 do i=1,nvar-1 02168 vector(:,i) = scalar(i) 02169 enddo 02170 endif 02171 02172 if (iostat /= 0) then 02173 inquire(iUnit,name=ername) 02174 call report_file_read_error(ername) 02175 endif 02176 02177 return 02178 9000 iseof = .true. 02179 02180 end subroutine velocity_Boundary_read 02181 02182 end module boundaryconditions
1.7.6.1
