|
XBeach
|
00001 module initialize_module 00002 use typesandkinds 00003 implicit none 00004 private 00005 public setbathy_init, grid_bathy, drifter_init, wave_init, sed_init, flow_init, discharge_init, hotstart_init_1, hotstart_init_2 00006 save 00007 integer imin_ee,imax_ee,jmin_ee,jmax_ee 00008 integer imin_uu,imax_uu,jmin_uu,jmax_uu 00009 integer imin_vv,imax_vv,jmin_vv,jmax_vv 00010 integer imin_zs,imax_zs,jmin_zs,jmax_zs 00011 00012 contains 00013 subroutine grid_bathy(s,par) 00014 00015 use params, only: parameters 00016 use spaceparamsdef 00017 use spaceparams, only: gridprops 00018 use xmpi_module 00019 use general_mpi_module 00020 use logging_module, only: writelog, report_file_read_error 00021 use typesandkinds, only: slen 00022 use paramsconst 00023 #ifdef USEMPI 00024 use mpi 00025 #endif 00026 00027 implicit none 00028 00029 type(spacepars) :: s 00030 type(parameters) :: par 00031 00032 integer :: i,ier,idum,n,m,ier2,ier3,dum 00033 integer :: j 00034 integer :: itheta 00035 real*8 :: degrad 00036 character(slen) :: line 00037 00038 s%nx = par%nx 00039 s%ny = par%ny 00040 s%dx = par%dx 00041 s%dy = par%dy 00042 s%xori = par%xori 00043 s%yori = par%yori 00044 s%alfa = par%alfa 00045 s%posdwn = par%posdwn 00046 s%vardx = par%vardx 00047 00048 if (s%alfa.lt.0) then 00049 s%alfa = 360.d0+s%alfa 00050 endif 00051 00052 s%alfa = s%alfa*atan(1.0d0)/45.d0 00053 ! Robert: huh? 00054 s%posdwn = s%posdwn*sign(s%posdwn,1.d0) 00055 ! end huh? 00056 00057 if(.not. (xmaster .or. xomaster)) return 00058 if(.true.) then 00059 allocate(s%x(1:s%nx+1,1:s%ny+1)) 00060 allocate(s%y(1:s%nx+1,1:s%ny+1)) 00061 allocate(s%xz(1:s%nx+1,1:s%ny+1)) 00062 allocate(s%yz(1:s%nx+1,1:s%ny+1)) 00063 allocate(s%xu(1:s%nx+1,1:s%ny+1)) 00064 allocate(s%yu(1:s%nx+1,1:s%ny+1)) 00065 allocate(s%xv(1:s%nx+1,1:s%ny+1)) 00066 allocate(s%yv(1:s%nx+1,1:s%ny+1)) 00067 allocate(s%zb(1:s%nx+1,1:s%ny+1)) 00068 allocate(s%zb0(1:s%nx+1,1:s%ny+1)) 00069 allocate(s%dzbdx(1:s%nx+1,1:s%ny+1)) 00070 allocate(s%dzbdy(1:s%nx+1,1:s%ny+1)) 00071 allocate(s%dsu(1:s%nx+1,1:s%ny+1)) 00072 allocate(s%dsv(1:s%nx+1,1:s%ny+1)) 00073 allocate(s%dsz(1:s%nx+1,1:s%ny+1)) 00074 allocate(s%dsc(1:s%nx+1,1:s%ny+1)) 00075 allocate(s%dnu(1:s%nx+1,1:s%ny+1)) 00076 allocate(s%dnv(1:s%nx+1,1:s%ny+1)) 00077 allocate(s%dnz(1:s%nx+1,1:s%ny+1)) 00078 allocate(s%dnc(1:s%nx+1,1:s%ny+1)) 00079 allocate(s%dsdnui(1:s%nx+1,1:s%ny+1)) 00080 allocate(s%dsdnvi(1:s%nx+1,1:s%ny+1)) 00081 allocate(s%dsdnzi(1:s%nx+1,1:s%ny+1)) 00082 allocate(s%alfau(1:s%nx+1,1:s%ny+1)) 00083 allocate(s%alfav(1:s%nx+1,1:s%ny+1)) 00084 allocate(s%alfaz(1:s%nx+1,1:s%ny+1)) 00085 allocate(s%sdist(1:s%nx+1,1:s%ny+1)) 00086 allocate(s%ndist(1:s%nx+1,1:s%ny+1)) 00087 endif 00088 00089 call writelog('l','' ,'------------------------------------') 00090 call writelog('ls','','Building Grid and Bathymetry') 00091 call writelog('l','', '------------------------------------') 00092 00093 ! 00094 ! Create grid and bathymetry 00095 ! 00096 ! Jaap make switch here to read XBeach or Delft3D format respectively 00097 ! 00098 ! wv in the following select case construct, s%zb, s%x and s%y are determined 00099 ! and they are read from file 00100 ! we let xmaster read these entities and send them to xomaster. 00101 ! because xomaster has also a need for these entities 00102 ! 00103 00104 if (xmaster) then 00105 select case(par%gridform) 00106 case(GRIDFORM_XBEACH) 00107 select case(s%vardx) 00108 case(0) 00109 if (par%setbathy .ne. 1 .and. par%hotstart .ne. 1) then 00110 open(31,file=par%depfile) 00111 do j=1,s%ny+1 00112 read(31,*,iostat=ier)(s%zb(i,j),i=1,s%nx+1) 00113 if (ier .ne. 0) then 00114 call report_file_read_error(par%depfile) 00115 endif 00116 end do 00117 close(31) 00118 endif 00119 do j=1,s%ny+1 00120 do i=1,s%nx+1 00121 s%x(i,j)=(i-1)*s%dx 00122 s%y(i,j)=(j-1)*s%dy 00123 end do 00124 end do 00125 case (1) ! Robert keep vardx == 1 for backwards compatibility?? 00126 if (par%setbathy .ne. 1 .and. par%hotstart .ne. 1) then 00127 open (31,file=par%depfile) 00128 read (31,*,iostat=ier)((s%zb(i,j),i=1,s%nx+1),j=1,s%ny+1) 00129 if (ier .ne. 0) then 00130 call report_file_read_error(par%depfile) 00131 endif 00132 close(31) 00133 endif 00134 00135 open (32,file=par%xfile) 00136 read (32,*,iostat=ier)((s%x(i,j),i=1,s%nx+1),j=1,s%ny+1) 00137 if (ier .ne. 0) then 00138 call report_file_read_error(par%xfile) 00139 endif 00140 close(32) 00141 00142 if (s%ny>0 .and. par%yfile/=' ') then 00143 open (33,file=par%yfile) 00144 read (33,*,iostat=ier)((s%y(i,j),i=1,s%nx+1),j=1,s%ny+1) 00145 if (ier .ne. 0) then 00146 call report_file_read_error(par%yfile) 00147 endif 00148 close(33) 00149 elseif (s%ny==0 .and. par%yfile/=' ') then 00150 open (33,file=par%yfile) 00151 read (33,*,iostat=ier)((s%y(i,j),i=1,s%nx+1),j=1,s%ny+1) 00152 if (ier .ne. 0) then 00153 call report_file_read_error(par%yfile) 00154 endif 00155 close(33) 00156 else 00157 s%y = 0.d0 00158 end if 00159 !endif 00160 case default 00161 call writelog('esl','','Invalid value for vardx: ',par%vardx) 00162 call halt_program 00163 end select 00164 00165 case (GRIDFORM_DELFT3D) 00166 ! 00167 ! Gridfile 00168 ! 00169 open(31,file=par%xyfile,status='old',iostat=ier) 00170 if (ier .ne. 0) then 00171 call report_file_read_error(par%xyfile) 00172 endif 00173 ! http://oss.deltares.nl/documents/183920/185723/Delft3D-FLOW_User_Manual.pdf section A.2.2 00174 ! skip comment text in file... 00175 do 00176 read(31,'(a)',iostat=ier)line 00177 if (line(1:1)/='*') then 00178 exit 00179 endif 00180 enddo 00181 read(31,*,iostat=ier2) idum,idum 00182 ! new grid format now specifies missing value, so catch this error 00183 if (ier2 .ne. 0) then 00184 read(31,*,iostat=ier2) idum,idum 00185 endif 00186 read(31,*,iostat=ier3) idum,idum,idum 00187 ! if any iostat is still /= 0 then there is an error reading the file 00188 if (ier+ier2+ier3 .ne. 0) then 00189 call report_file_read_error(par%xyfile) 00190 endif 00191 00192 read(31,*,iostat=ier) & 00193 (line,dum,(s%x(m,n),m=1,s%nx+1),n=1,s%ny+1), & 00194 (line,dum,(s%y(m,n),m=1,s%nx+1),n=1,s%ny+1) 00195 if (ier .ne. 0) then 00196 call report_file_read_error(par%xyfile) 00197 endif 00198 00199 close(31) 00200 ! 00201 ! Depfile 00202 ! 00203 if (par%setbathy .ne. 1) then 00204 open(33,file=par%depfile,status='old') 00205 do n=1,s%ny+1 00206 read(33,*,iostat=ier)(s%zb(m,n),m=1,s%nx+1) 00207 if (ier .ne. 0) then 00208 call report_file_read_error(par%depfile) 00209 endif 00210 enddo 00211 close(33) 00212 endif 00213 end select 00214 endif 00215 ! xmaster 00216 00217 degrad=par%px/180.d0 00218 ! send zb, x, y to xomaster 00219 00220 #ifdef USEMPI 00221 ! wwvv todo use xmpi_send 00222 if(xmaster) then 00223 ! MPI_SEND(BUF, COUNT, DATATYPE, DEST, TAG, COMM, IERROR) 00224 ! <type> BUF(*) 00225 ! INTEGER COUNT, DATATYPE, DEST, TAG, COMM, IERROR 00226 call MPI_Send(s%zb, size(s%zb), MPI_DOUBLE_PRECISION, xmpi_omaster, 1, xmpi_ocomm, ier) 00227 call MPI_Send(s%x, size(s%x), MPI_DOUBLE_PRECISION, xmpi_omaster, 2, xmpi_ocomm, ier) 00228 call MPI_Send(s%y, size(s%y), MPI_DOUBLE_PRECISION, xmpi_omaster, 3, xmpi_ocomm, ier) 00229 else 00230 ! MPI_RECV(BUF, COUNT, DATATYPE, SOURCE, TAG, COMM, STATUS, IERROR) 00231 ! <type> BUF(*) 00232 ! INTEGER COUNT, DATATYPE, SOURCE, TAG, COMM 00233 ! INTEGER STATUS(MPI_STATUS_SIZE), IERROR 00234 call MPI_Recv(s%zb, size(s%zb), MPI_DOUBLE_PRECISION, xmpi_imaster, 1, xmpi_ocomm, MPI_STATUS_IGNORE, ier) 00235 call MPI_Recv(s%x, size(s%x), MPI_DOUBLE_PRECISION, xmpi_imaster, 2, xmpi_ocomm, MPI_STATUS_IGNORE, ier) 00236 call MPI_Recv(s%y, size(s%y), MPI_DOUBLE_PRECISION, xmpi_imaster, 3, xmpi_ocomm, MPI_STATUS_IGNORE, ier) 00237 endif 00238 #endif 00239 00240 if(.true.) then 00241 00242 s%zb=-s%zb*s%posdwn 00243 ! Make sure that at the lateral boundaries the bathymetry is alongshore uniform 00244 if (s%ny>0) then 00245 s%zb(:,1) = s%zb(:,2) 00246 s%zb(:,s%ny+1) = s%zb(:,s%ny) 00247 endif 00248 s%zb(1,:)=s%zb(2,:) 00249 s%zb(s%nx+1,:)=s%zb(s%nx,:) 00250 00251 call gridprops (s) 00252 00253 s%zb0 = s%zb 00254 if(par%swave==1) then 00255 ! 00256 ! Specify theta-grid 00257 ! 00258 ! 00259 ! from Nautical wave directions in degrees to Cartesian wave directions in radian !!! 00260 ! 00261 ! s%theta0=(1.5d0*par%px-s%alfa)-par%dir0*atan(1.d0)/45.d0 ! Updated in waveparams.f90 for instat 4,5,6,7 00262 s%theta0=(1.5d0*par%px)-par%dir0*atan(1.d0)/45.d0 ! Updated in waveparams.f90 for instat 4,5,6,7 00263 do while(s%theta0<-par%px) 00264 s%theta0=s%theta0+2.d0*par%px 00265 enddo 00266 do while(s%theta0>par%px) 00267 s%theta0=s%theta0-2.d0*par%px 00268 enddo 00269 00270 !degrad=par%px/180.d0 00271 ! 00272 if (par%thetanaut==1) then 00273 s%thetamin=(270-par%thetamax)*degrad 00274 s%thetamax=(270-par%thetamin)*degrad 00275 else 00276 ! rotate theta grid to world coordinates for backwards compatibility 00277 s%thetamin=par%thetamin+s%alfa/degrad 00278 s%thetamax=par%thetamax+s%alfa/degrad 00279 00280 s%thetamin=s%thetamin*degrad 00281 s%thetamax=s%thetamax*degrad 00282 endif 00283 00284 ! try and fix whatever strange angles given in params.txt 00285 s%thetamin = mod(s%thetamin,2.d0*par%px) 00286 s%thetamax = mod(s%thetamax,2.d0*par%px) 00287 ! thetamin should always be smaller than thetamax 00288 if(s%thetamin>=s%thetamax) then 00289 if (s%thetamax>=0.d0) then 00290 do while(s%thetamin>=s%thetamax) 00291 s%thetamin = s%thetamin-2.d0*par%px 00292 enddo 00293 else 00294 do while(s%thetamin>s%thetamax) 00295 s%thetamax = s%thetamax+2.d0*par%px 00296 enddo 00297 endif 00298 elseif(s%thetamax>s%thetamin+2.d0*par%px) then 00299 do while(s%thetamax>s%thetamin+2.d0*par%px) ! note, most should already be captured by mod statements above, 00300 ! but can still occur under strange conditions 00301 s%thetamin = s%thetamin+2.d0*par%px 00302 enddo 00303 endif 00304 00305 !if (s%thetamax>par%px) then 00306 ! s%thetamax=s%thetamax-2*par%px 00307 ! s%thetamin=s%thetamin-2*par%px 00308 !endif 00309 !if (s%thetamin<-par%px) then 00310 ! s%thetamax=s%thetamax+2*par%px 00311 ! s%thetamin=s%thetamin+2*par%px 00312 !endif 00313 00314 00315 if(par%single_dir==0) then 00316 s%dtheta=par%dtheta*degrad 00317 s%ntheta=nint((s%thetamax-s%thetamin)/s%dtheta) 00318 else 00319 s%dtheta=s%thetamax-s%thetamin 00320 s%ntheta=1 00321 endif 00322 else 00323 s%dtheta=2*par%px 00324 s%ntheta = 1 00325 ! below needed sometimes (like combination wavemodel = nonh and wavebctype = params) 00326 s%theta0=(1.5d0*par%px)-par%dir0*atan(1.d0)/45.d0 00327 do while(s%theta0<-par%px) 00328 s%theta0=s%theta0+2.d0*par%px 00329 enddo 00330 do while(s%theta0>par%px) 00331 s%theta0=s%theta0-2.d0*par%px 00332 enddo 00333 endif 00334 00335 ! Always allocate room incase of output request and memory sharing 00336 allocate(s%theta(1:s%ntheta)) 00337 allocate(s%thet(1:s%nx+1,1:s%ny+1,1:s%ntheta)) 00338 allocate(s%costh(1:s%nx+1,1:s%ny+1,1:s%ntheta)) 00339 allocate(s%sinth(1:s%nx+1,1:s%ny+1,1:s%ntheta)) 00340 00341 if (par%single_dir==1) then 00342 s%dtheta_s=par%dtheta_s*degrad 00343 s%ntheta_s=nint((s%thetamax-s%thetamin)/s%dtheta_s) 00344 allocate(s%theta_s(1:s%ntheta_s)) 00345 allocate(s%thet_s(1:s%nx+1,1:s%ny+1,1:s%ntheta_s)) 00346 allocate(s%costh_s(1:s%nx+1,1:s%ny+1,1:s%ntheta_s)) 00347 allocate(s%sinth_s(1:s%nx+1,1:s%ny+1,1:s%ntheta_s)) 00348 else 00349 s%dtheta_s=2*par%px 00350 s%ntheta_s=0 00351 allocate(s%theta_s(1:s%ntheta)) 00352 allocate(s%thet_s(1:s%nx+1,1:s%ny+1,1:s%ntheta)) 00353 allocate(s%costh_s(1:s%nx+1,1:s%ny+1,1:s%ntheta)) 00354 allocate(s%sinth_s(1:s%nx+1,1:s%ny+1,1:s%ntheta)) 00355 endif 00356 00357 ! Always allocate room incase of output request and memory sharing 00358 00359 00360 if (par%swave==1) then 00361 do itheta=1,s%ntheta 00362 s%theta(itheta)=mod(s%thetamin+s%dtheta/2+s%dtheta*(itheta-1),2*par%px) 00363 end do 00364 00365 do itheta=1,s%ntheta 00366 do j=1,s%ny+1 00367 do i=1,s%nx+1 00368 s%thet(i,j,itheta) = s%theta(itheta) 00369 s%costh(i,j,itheta)=cos(mod(s%theta(itheta)-s%alfaz(i,j),2*par%px)) 00370 s%sinth(i,j,itheta)=sin(mod(s%theta(itheta)-s%alfaz(i,j),2*par%px)) 00371 enddo 00372 enddo 00373 enddo 00374 00375 if (par%single_dir==1) then 00376 do itheta=1,s%ntheta_s 00377 s%theta_s(itheta)=mod(s%thetamin+s%dtheta_s/2+s%dtheta_s*(itheta-1),2*par%px) 00378 end do 00379 00380 do itheta=1,s%ntheta_s 00381 do j=1,s%ny+1 00382 do i=1,s%nx+1 00383 s%thet_s(i,j,itheta) = s%theta_s(itheta) 00384 s%costh_s(i,j,itheta)=cos(mod(s%theta_s(itheta)-s%alfaz(i,j),2*par%px)) 00385 s%sinth_s(i,j,itheta)=sin(mod(s%theta_s(itheta)-s%alfaz(i,j),2*par%px)) 00386 enddo 00387 enddo 00388 enddo 00389 endif 00390 endif 00391 endif 00392 00393 !if (xmaster) then 00394 if(.true.) then 00395 ! Initialize dzbdx, dzbdy 00396 do j=1,s%ny+1 00397 do i=1,s%nx 00398 s%dzbdx(i,j)=(s%zb(i+1,j)-s%zb(i,j))/s%dsu(i,j) 00399 enddo 00400 enddo 00401 ! dummy, needed to keep compiler happy 00402 s%dzbdx(s%nx+1,:)=s%dzbdx(s%nx,:) 00403 00404 do j=1,s%ny 00405 do i=1,s%nx+1 00406 s%dzbdy(i,j)=(s%zb(i,j+1)-s%zb(i,j))/s%dnv(i,j) 00407 enddo 00408 enddo 00409 if (s%ny>0) then 00410 s%dzbdy(:,s%ny+1)=s%dzbdy(:,s%ny) 00411 endif 00412 endif 00413 00414 end subroutine grid_bathy 00415 00416 subroutine setbathy_init(s,par) 00417 use params, only: parameters 00418 use spaceparamsdef, only: spacepars 00419 use filefunctions, only: create_new_fid 00420 use logging_module, only: report_file_read_error 00421 use interp, only: linear_interp 00422 use xmpi_module, only: xmaster 00423 00424 type(spacepars) :: s 00425 type(parameters) :: par 00426 00427 integer :: i,j,it 00428 integer :: ier,fid,dummy 00429 00430 if(.not. xmaster) return 00431 00432 if (par%setbathy==1) then 00433 allocate(s%setbathy(s%nx+1,s%ny+1,par%nsetbathy)) 00434 allocate(s%tsetbathy(par%nsetbathy)) 00435 ! start file read 00436 fid = create_new_fid() 00437 open (fid,file=par%setbathyfile) 00438 do it=1,par%nsetbathy 00439 read(fid,*,iostat=ier)s%tsetbathy(it) 00440 if (ier .ne. 0) then 00441 call report_file_read_error(par%setbathyfile) 00442 endif 00443 do j=1,s%ny+1 00444 read(fid,*,iostat=ier)(s%setbathy(i,j,it),i=1,s%nx+1) 00445 if (ier .ne. 0) then 00446 call report_file_read_error(par%setbathyfile) 00447 endif 00448 enddo 00449 enddo 00450 close(fid) 00451 ! Interpolate initial bathymetry 00452 do j=1,s%ny+1 00453 do i=1,s%nx+1 00454 call LINEAR_INTERP(s%tsetbathy,s%setbathy(i,j,:),par%nsetbathy, & 00455 0.d0,s%zb(i,j),dummy) 00456 enddo 00457 enddo 00458 s%zb0 = s%zb 00459 else 00460 ! give MPI bcast a memory address 00461 allocate(s%setbathy(0,0,0)) 00462 allocate(s%tsetbathy(0)) 00463 endif 00464 end subroutine setbathy_init 00465 00466 subroutine wave_init (s,par) 00467 use params, only: parameters 00468 use spaceparamsdef, only: spacepars 00469 use logging_module, only: report_file_read_error 00470 use xmpi_module, only: xmaster 00471 use wave_functions_module, only: dispersion 00472 00473 use paramsconst 00474 00475 IMPLICIT NONE 00476 00477 type(spacepars),target :: s 00478 type(parameters) :: par 00479 00480 integer :: itheta, i, j, ier 00481 logical :: exists 00482 00483 if(.not. xmaster) return 00484 00485 allocate(s%thetamean(1:s%nx+1,1:s%ny+1)) 00486 allocate(s%Fx(1:s%nx+1,1:s%ny+1)) 00487 allocate(s%Fy(1:s%nx+1,1:s%ny+1)) 00488 allocate(s%Sxx(1:s%nx+1,1:s%ny+1)) 00489 allocate(s%Sxy(1:s%nx+1,1:s%ny+1)) 00490 allocate(s%Syy(1:s%nx+1,1:s%ny+1)) 00491 allocate(s%n(1:s%nx+1,1:s%ny+1)) 00492 allocate(s%H(1:s%nx+1,1:s%ny+1)) 00493 allocate(s%cgx(1:s%nx+1,1:s%ny+1,1:s%ntheta)) 00494 allocate(s%cgy(1:s%nx+1,1:s%ny+1,1:s%ntheta)) 00495 allocate(s%cx(1:s%nx+1,1:s%ny+1,1:s%ntheta)) 00496 allocate(s%cy(1:s%nx+1,1:s%ny+1,1:s%ntheta)) 00497 allocate(s%ctheta(1:s%nx+1,1:s%ny+1,1:s%ntheta)) 00498 allocate(s%sigt(1:s%nx+1,1:s%ny+1,1:s%ntheta)) 00499 allocate(s%ee(1:s%nx+1,1:s%ny+1,1:s%ntheta)) 00500 allocate(s%rr(1:s%nx+1,1:s%ny+1,1:s%ntheta)) 00501 if (par%single_dir==1) then ! Robert: always allocate these variables 00502 allocate(s%cgx_s(1:s%nx+1,1:s%ny+1,1:s%ntheta_s)) 00503 allocate(s%cgy_s(1:s%nx+1,1:s%ny+1,1:s%ntheta_s)) 00504 allocate(s%ctheta_s(1:s%nx+1,1:s%ny+1,1:s%ntheta_s)) 00505 allocate(s%ee_s(1:s%nx+1,1:s%ny+1,1:s%ntheta_s)) 00506 else 00507 allocate(s%cgx_s(0,0,0)) 00508 allocate(s%cgy_s(0,0,0)) 00509 allocate(s%ctheta_s(0,0,0)) 00510 allocate(s%ee_s(0,0,0)) 00511 endif 00512 allocate(s%sigm(1:s%nx+1,1:s%ny+1)) 00513 allocate(s%c(1:s%nx+1,1:s%ny+1)) 00514 allocate(s%cg(1:s%nx+1,1:s%ny+1)) 00515 allocate(s%k(1:s%nx+1,1:s%ny+1)) 00516 allocate(s%ui(1:s%nx+1,1:s%ny+1)) 00517 allocate(s%vi(1:s%nx+1,1:s%ny+1)) 00518 allocate(s%E(1:s%nx+1,1:s%ny+1)) 00519 allocate(s%R(1:s%nx+1,1:s%ny+1)) 00520 allocate(s%urms(1:s%nx+1,1:s%ny+1)) 00521 allocate(s%D(1:s%nx+1,1:s%ny+1)) 00522 allocate(s%Df(1:s%nx+1,1:s%ny+1)) 00523 allocate(s%fw(1:s%nx+1,1:s%ny+1)) 00524 allocate(s%Dp(1:s%nx+1,1:s%ny+1)) 00525 allocate(s%Qb(1:s%nx+1,1:s%ny+1)) 00526 allocate(s%ust(1:s%nx+1,1:s%ny+1)) 00527 allocate(s%uwf(1:s%nx+1,1:s%ny+1)) 00528 allocate(s%vwf(1:s%nx+1,1:s%ny+1)) 00529 allocate(s%ustr(1:s%nx+1,1:s%ny+1)) 00530 allocate(s%usd(1:s%nx+1,1:s%ny+1)) 00531 allocate(s%bi(1:s%ny+1)) 00532 allocate(s%DR(1:s%nx+1,1:s%ny+1)) 00533 allocate(s%umwci (1:s%nx+1,1:s%ny+1)) 00534 allocate(s%vmwci (1:s%nx+1,1:s%ny+1)) 00535 allocate(s%zswci (1:s%nx+1,1:s%ny+1)) 00536 allocate(s%BR(1:s%nx+1,1:s%ny+1)) 00537 allocate(s%hhw(1:s%nx+1,1:s%ny+1)) 00538 allocate(s%hhws(1:s%nx+1,1:s%ny+1)) 00539 allocate(s%uws(1:s%nx+1,1:s%ny+1)) 00540 allocate(s%vws(1:s%nx+1,1:s%ny+1)) 00541 allocate(s%hhwcins(1:s%nx+1,1:s%ny+1)) 00542 allocate(s%uwcins(1:s%nx+1,1:s%ny+1)) 00543 allocate(s%vwcins(1:s%nx+1,1:s%ny+1)) 00544 !allocate(s%tm(1:s%nx+1,1:s%ny+1)) 00545 ! 00546 ! Initial condition 00547 ! 00548 s%ee = 0.d0 00549 s%thetamean = 0.d0 00550 s%Fx = 0.d0 00551 s%Fy = 0.d0 00552 s%Sxx = 0.d0 00553 s%Sxy = 0.d0 00554 s%Syy = 0.d0 00555 s%n = 0.d0 00556 s%H = 0.d0 00557 s%cgx = 0.d0 00558 s%cgy = 0.d0 00559 s%cx = 0.d0 00560 s%cy = 0.d0 00561 s%ctheta = 0.d0 00562 if (par%single_dir==1) then 00563 s%ee_s = 0.d0 00564 s%cgx_s = 0.d0 00565 s%cgy_s = 0.d0 00566 s%ctheta_s = 0.d0 00567 endif 00568 s%sigt = 0.d0 00569 s%rr = 0.d0 00570 s%sigm = 0.d0 00571 s%c = 0.d0 00572 s%cg = 0.d0 00573 s%k = 0.d0 00574 s%ui = 0.d0 00575 s%vi = 0.d0 00576 s%E = 0.d0 00577 s%R = 0.d0 00578 s%urms = 0.d0 00579 s%D = 0.d0 00580 s%Qb = 0.d0 00581 s%ust = 0.d0 00582 s%uwf = 0.d0 00583 s%vwf = 0.d0 00584 s%ustr = 0.d0 00585 s%usd = 0.d0 00586 s%bi = 0.d0 00587 s%DR = 0.d0 00588 s%Df = 0.d0 00589 s%BR = par%Beta 00590 s%hhw = s%hh 00591 s%hhws = s%hh 00592 s%uws = s%u 00593 s%vws = s%v 00594 s%hhwcins = s%hh 00595 s%uwcins = s%u 00596 s%vwcins = s%v 00597 s%isSet_Vbc = 0 00598 !s%tm = 0.d0 00599 00600 00601 ! introduce intrinsic frequencies for wave action 00602 if ( par%wbctype==WBCTYPE_PARAMETRIC .or. & 00603 par%wbctype==WBCTYPE_JONS_TABLE .or. & 00604 par%wbctype==WBCTYPE_SWAN .or. & 00605 par%wbctype==WBCTYPE_VARDENS .or. & 00606 par%wbctype==WBCTYPE_REUSE .or. & 00607 par%wbctype==WBCTYPE_TS_NONH & 00608 ) par%Trep=10.d0 00609 !Robert 00610 ! incorrect values are computed below for instat = 4/5/6/7 00611 ! in this case right values are computed in wave params.f90 00612 if ( par%wbctype==WBCTYPE_PARAMETRIC .or. & 00613 par%wbctype==WBCTYPE_JONS_TABLE .or. & 00614 par%wbctype==WBCTYPE_SWAN .or. & 00615 par%wbctype==WBCTYPE_VARDENS) then 00616 if(xmaster) call spectral_wave_init (s,par) ! only used by xmaster 00617 endif 00618 00619 do itheta=1,s%ntheta 00620 s%sigt(:,:,itheta) = 2*par%px/par%Trep 00621 end do 00622 s%sigm = sum(s%sigt,3)/s%ntheta 00623 call dispersion(par,s,s%hh) 00624 ! 00625 inquire(file=par%wavfricfile,exist=exists) 00626 if ((exists)) then 00627 open(723,file=par%wavfricfile) 00628 do j=1,s%ny+1 00629 read(723,*,iostat=ier)(s%fw(i,j),i=1,s%nx+1) 00630 if (ier .ne. 0) then 00631 call report_file_read_error(par%wavfricfile) 00632 endif 00633 enddo 00634 close(723) 00635 else 00636 s%fw=par%wavfriccoef 00637 endif 00638 00639 00640 end subroutine wave_init 00641 00642 00643 subroutine spectral_wave_init(s,par) 00644 use params, only: parameters 00645 use spaceparamsdef, only: spacepars 00646 use filefunctions, only: create_new_fid 00647 use logging_module, only: writelog 00648 use spectral_wave_bc_module, only: lastwaveelevation, n_index_loc, bcfiles, bccount, nspectra, reuseall, spectrumendtime 00649 use xmpi_module, only: halt_program 00650 use typesandkinds, only: slen 00651 00652 type(spacepars) :: s 00653 type(parameters) :: par 00654 00655 integer :: fid,err 00656 integer :: i 00657 integer,dimension(1) :: minlocation 00658 character(slen) :: testline 00659 real*8,dimension(:),allocatable :: xspec,yspec,mindist 00660 real*8 :: mindistr 00661 00662 call writelog('l','','--------------------------------') 00663 call writelog('l','','Initializing spectral wave boundary conditions ') 00664 ! Initialize that wave boundary conditions need to be calculated (first time at least) 00665 ! Stored and defined in spectral_wave_bc_module 00666 reuseall = .false. 00667 ! Initialize the number of times wave boundary conditions have been generated. 00668 ! Stored and defined in spectral_wave_bc_module 00669 bccount = 0 00670 ! Initialize bcendtime to zero. 00671 ! Stored and defined in spectral_wave_bc_module 00672 spectrumendtime = 0.d0 00673 ! Initialise lastwaveheight to zero 00674 ! Stored and defined in spectral_wave_bc_module 00675 allocate(lastwaveelevation(s%ny+1,s%ntheta)) 00676 lastwaveelevation = 0.d0 00677 00678 if (par%nspectrumloc<1) then 00679 call writelog('ewls','','number of boundary spectra (''nspectrumloc'') may not be less than 1') 00680 call halt_program 00681 endif 00682 00683 ! open location list file 00684 fid = create_new_fid() 00685 open(fid,file=par%bcfile,status='old',form='formatted') 00686 ! check for LOCLIST 00687 read(fid,*)testline 00688 if (trim(testline)=='LOCLIST') then 00689 allocate(n_index_loc(par%nspectrumloc)) ! stored and defined in spectral_wave_bc_module 00690 allocate(bcfiles(par%nspectrumloc)) ! stored and defined in spectral_wave_bc_module 00691 allocate(xspec(par%nspectrumloc)) 00692 allocate(yspec(par%nspectrumloc)) 00693 allocate(mindist(s%ny+1)) 00694 do i=1,par%nspectrumloc 00695 ! read x,y and file name per location 00696 read(fid,*,IOSTAT=err)xspec(i),yspec(i),bcfiles(i)%fname 00697 bcfiles(i)%listline = 0 00698 if (err /= 0) then 00699 ! something has gone wrong during the read of this file 00700 call writelog('ewls','a,i0,a,a)','error reading line ',i+1,' of file ',par%bcfile) 00701 call writelog('ewls','','check file for format errors and ensure the number of ',& 00702 'lines is equal to nspectrumloc') 00703 call halt_program 00704 endif 00705 enddo 00706 ! convert x,y locations of the spectra to the nearest grid points on s=1 boundary 00707 do i=1,par%nspectrumloc 00708 mindist=sqrt((xspec(i)-s%xz(1,:))**2+(yspec(i)-s%yz(1,:))**2) 00709 ! look up the location of the found minimum 00710 minlocation=minloc(mindist) 00711 ! minimum distance 00712 mindistr = mindist(minlocation(1)) 00713 ! store location 00714 n_index_loc(i) = minlocation(1) 00715 call writelog('ls','(a,i0,a,i0)','Spectrum ',i,' placed at n = ',minlocation(1)) 00716 call writelog('ls','(a,f0.3)','Distance spectrum to grid: ',mindistr) 00717 enddo 00718 nspectra = par%nspectrumloc ! stored and defined in spectral_wave_bc_module 00719 deallocate(xspec,yspec,mindist) 00720 else 00721 if (par%nspectrumloc==1) then 00722 allocate(n_index_loc(1)) ! stored and defined in spectral_wave_bc_module 00723 allocate(bcfiles(1)) ! stored and defined in spectral_wave_bc_module 00724 n_index_loc = 1 00725 bcfiles(1)%fname = par%bcfile 00726 bcfiles(1)%listline = 0 ! for files that have multiple lines, set listline to 0 00727 nspectra = 1 ! stored and defined in spectral_wave_bc_module 00728 else 00729 call writelog('ewls','','if nspectrumloc>1 then bcfile should contain spectra locations with LOCLIST header') 00730 close(fid) 00731 call halt_program 00732 endif 00733 endif 00734 close(fid) 00735 00736 call writelog('l','','--------------------------------') 00737 00738 end subroutine spectral_wave_init 00739 00740 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00741 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00742 00743 subroutine flow_init (s,par) 00744 use params, only: parameters 00745 use spaceparamsdef 00746 use logging_module, only: writelog, report_file_read_error 00747 use xmpi_module 00748 use compute_tide_module, only: tide_init 00749 use paramsconst 00750 00751 00752 IMPLICIT NONE 00753 00754 type(spacepars),target :: s 00755 type(parameters), intent(in) :: par 00756 00757 integer*4 :: i,j,ig,indt 00758 logical :: exists 00759 logical :: offshoreregime 00760 integer :: indoff,indbay,ier 00761 00762 real*8,dimension(:),allocatable :: xzs0,yzs0,szs0 00763 real*8,dimension(:,:),allocatable :: vmagvold,vmaguold 00764 real*8 :: flowerr 00765 00766 00767 if(.not. xmaster) return 00768 00769 if (par%hotstart .ne. 1) then 00770 allocate(s%zs(1:s%nx+1,1:s%ny+1)) 00771 allocate(s%uu(1:s%nx+1,1:s%ny+1)) 00772 allocate(s%vv(1:s%nx+1,1:s%ny+1)) 00773 endif 00774 00775 allocate(s%dzsdt(1:s%nx+1,1:s%ny+1)) 00776 allocate(s%dzsdx(1:s%nx+1,1:s%ny+1)) 00777 allocate(s%dzsdy(1:s%nx+1,1:s%ny+1)) 00778 allocate(s%dzbdt(1:s%nx+1,1:s%ny+1)) 00779 allocate(s%dzbnow(1:s%nx+1,1:s%ny+1)) 00780 allocate(s%qx(1:s%nx+1,1:s%ny+1)) 00781 allocate(s%qy(1:s%nx+1,1:s%ny+1)) 00782 allocate(s%sedero(1:s%nx+1,1:s%ny+1)) 00783 allocate(s%ueu(1:s%nx+1,1:s%ny+1)) 00784 allocate(s%vev(1:s%nx+1,1:s%ny+1)) 00785 allocate(s%vmagu(1:s%nx+1,1:s%ny+1)) 00786 allocate(s%vmagv(1:s%nx+1,1:s%ny+1)) 00787 allocate(s%vmageu(1:s%nx+1,1:s%ny+1)) 00788 allocate(s%vmagev(1:s%nx+1,1:s%ny+1)) 00789 allocate(s%u(1:s%nx+1,1:s%ny+1)) 00790 allocate(s%v(1:s%nx+1,1:s%ny+1)) 00791 allocate(s%ue(1:s%nx+1,1:s%ny+1)) 00792 allocate(s%ve(1:s%nx+1,1:s%ny+1)) 00793 allocate(s%ue_sed(1:s%nx+1,1:s%ny+1)) 00794 allocate(s%ve_sed(1:s%nx+1,1:s%ny+1)) 00795 allocate(s%hold(1:s%nx+1,1:s%ny+1)) 00796 allocate(s%wetu(1:s%nx+1,1:s%ny+1)) 00797 allocate(s%wetv(1:s%nx+1,1:s%ny+1)) 00798 allocate(s%wetz(1:s%nx+1,1:s%ny+1)) 00799 allocate(s%wete(1:s%nx+1,1:s%ny+1)) 00800 allocate(s%hh(1:s%nx+1,1:s%ny+1)) 00801 allocate(s%hu(1:s%nx+1,1:s%ny+1)) 00802 allocate(s%hv(1:s%nx+1,1:s%ny+1)) 00803 allocate(s%hum(1:s%nx+1,1:s%ny+1)) 00804 allocate(s%hvm(1:s%nx+1,1:s%ny+1)) 00805 allocate(s%vu(1:s%nx+1,1:s%ny+1)) 00806 allocate(s%uv(1:s%nx+1,1:s%ny+1)) 00807 allocate(s%maxzs(1:s%nx+1,1:s%ny+1)) 00808 allocate(s%minzs(1:s%nx+1,1:s%ny+1)) 00809 allocate(s%taubx(1:s%nx+1,1:s%ny+1)) 00810 allocate(s%tauby(1:s%nx+1,1:s%ny+1)) 00811 allocate(s%ws(1:s%nx+1,1:s%ny+1)) 00812 allocate(s%wscrit(1:s%nx+1,1:s%ny+1)) 00813 allocate(s%wb(1:s%nx+1,1:s%ny+1)) 00814 allocate(s%nuh(1:s%nx+1,1:s%ny+1)) 00815 allocate(s%pres(1:s%nx+1,1:s%ny+1)) 00816 allocate(s%ph(1:s%nx+1,1:s%ny+1)) 00817 allocate(s%wi(2,1:s%ny+1)) 00818 allocate(s%dUi(2,1:s%ny+1)) 00819 allocate(s%dVi(2,1:s%ny+1)) 00820 allocate(s%zi(2,1:s%ny+1)) 00821 allocate(s%bedfriccoef(1:s%nx+1,1:s%ny+1)) 00822 allocate(s%cf(1:s%nx+1,1:s%ny+1)) 00823 allocate(s%cfu(1:s%nx+1,1:s%ny+1)) 00824 allocate(s%cfv(1:s%nx+1,1:s%ny+1)) 00825 allocate(s%zs0(1:s%nx+1,1:s%ny+1)) 00826 allocate(s%zs1(1:s%nx+1,1:s%ny+1)) ! water level minus tide 00827 allocate(s%dzs0dn(1:s%nx+1,1:s%ny+1)) 00828 allocate(s%zs0fac(1:s%nx+1,1:s%ny+1,2)) 00829 allocate(s%wm(1:s%nx+1,1:s%ny+1)) 00830 allocate(s%umean(1:s%nx+1,1:s%ny+1)) 00831 allocate(s%vmean(1:s%nx+1,1:s%ny+1)) 00832 allocate(s%ur(1:s%nx+1,1:s%ny+1)) 00833 allocate(s%xyzs01(2)) 00834 allocate(s%xyzs02(2)) 00835 allocate(s%xyzs03(2)) 00836 allocate(s%xyzs04(2)) 00837 00838 allocate(s%dU(1:s%nx+1,1:s%ny+1)) 00839 allocate(s%dV(1:s%nx+1,1:s%ny+1)) 00840 00841 allocate(s%sigz(1:par%nz)) 00842 allocate(s%uz(1:s%nx+1,1:s%ny+1,1:par%nz)) 00843 allocate(s%vz(1:s%nx+1,1:s%ny+1,1:par%nz)) 00844 allocate(s%ustz(1:s%nx+1,1:s%ny+1,1:par%nz)) 00845 allocate(s%nutz(1:s%nx+1,1:s%ny+1,1:par%nz)) 00846 00847 allocate(szs0(1:2)) 00848 allocate(xzs0(1:2)) 00849 allocate(yzs0(1:2)) 00850 allocate(vmagvold(1:s%nx+1,1:s%ny+1)) 00851 allocate(vmaguold(1:s%nx+1,1:s%ny+1)) 00852 00853 allocate(s%ududx(1:s%nx+1,1:s%ny+1)) 00854 allocate(s%vdvdy(1:s%nx+1,1:s%ny+1)) 00855 allocate(s%udvdx(1:s%nx+1,1:s%ny+1)) 00856 allocate(s%vdudy(1:s%nx+1,1:s%ny+1)) 00857 allocate(s%viscu(1:s%nx+1,1:s%ny+1)) 00858 allocate(s%viscv(1:s%nx+1,1:s%ny+1)) 00859 ! nonh hydrostatic flow initialisation 00860 allocate(s%breaking(1:s%nx+1,1:s%ny+1)) 00861 00862 ! vegetation initialization 00863 !allocate(s%vegtype(1:s%nx+1,1:s%ny+1)) 00864 ! allocate(s%nsecveg(1:s%nx+1,1:s%ny+1)) 00865 !allocate(s%Cdveg(1:s%nx+1,1:s%ny+1,1:3)) 00866 !allocate(s%ahveg(1:s%nx+1,1:s%ny+1,1:3)) 00867 !allocate(s%bveg(1:s%nx+1,1:s%ny+1,1:3)) 00868 !allocate(s%Nveg(1:s%nx+1,1:s%ny+1,1:3)) 00869 00870 ! added bed roughness due to second order effects 00871 allocate(s%taubx_add(1:s%nx+1,1:s%ny+1)) 00872 allocate(s%tauby_add(1:s%nx+1,1:s%ny+1)) 00873 00874 allocate(s%hstokes(1:s%nx+1,1:s%ny+1)) 00875 00876 ! Just to be sure! 00877 if(par%hotstart .ne. 1) then 00878 s%zs = 0.d0 00879 s%uu = 0.d0 00880 s%vv = 0.d0 00881 endif 00882 s%dzsdt = 0.0d0 00883 s%dzsdx = 0.0d0 00884 s%dzsdy = 0.0d0 00885 s%dzbdt = 0.0d0 00886 s%dzbnow = 0.d0 00887 s%qx = 0.0d0 00888 s%qy = 0.0d0 00889 s%sedero = 0.0d0 00890 s%ueu = 0.0d0 00891 s%vev = 0.0d0 00892 s%vmagu = 0.0d0 00893 s%vmagv = 0.0d0 00894 s%vmageu = 0.0d0 00895 s%vmagev = 0.0d0 00896 s%u = 0.0d0 00897 s%v = 0.0d0 00898 s%ue = 0.0d0 00899 s%ve = 0.0d0 00900 s%ue_sed = 0.0d0 00901 s%ve_sed = 0.0d0 00902 s%hold = 0.0d0 00903 s%wetu = 0 00904 s%wetv = 0 00905 s%wetz = 0 00906 s%wete = 0 00907 s%hh = 0.0d0 00908 s%hu = 0.0d0 00909 s%hv = 0.0d0 00910 s%hum = 0.0d0 00911 s%hvm = 0.0d0 00912 s%vu = 0.0d0 00913 s%uv = 0.0d0 00914 s%dU = 0.0d0 00915 s%dV = 0.0d0 00916 s%maxzs = 0.0d0 00917 s%minzs = 0.0d0 00918 s%taubx = 0.0d0 00919 s%tauby = 0.0d0 00920 s%ws = 0.0d0 00921 s%wb = 0.0d0 00922 s%nuh = 0.0d0 00923 s%pres = 0.0d0 00924 s%ph = 0.0d0 00925 s%wi = 0.0d0 00926 s%zi = 0.0d0 00927 s%cf = 0.0d0 00928 s%zs0 = 0.0d0 00929 s%zs0fac = 0.0d0 00930 s%wm = 0.0d0 00931 s%umean = 0.0d0 00932 s%vmean = 0.0d0 00933 s%ur = 0.0d0 00934 s%xyzs01 = 0.0d0 00935 s%xyzs02 = 0.0d0 00936 s%xyzs03 = 0.0d0 00937 s%xyzs04 = 0.0d0 00938 00939 szs0 = 0.0d0 00940 xzs0 = 0.0d0 00941 yzs0 = 0.0d0 00942 00943 ! TODO: do this properly.... 00944 ! All variables above, should be initialized below (for all cells) 00945 s%zs0 = 0.0d0 00946 s%dzs0dn = 0.d0 00947 s%ue = 0.0d0 00948 s%ve = 0.0d0 00949 s%ws = 0.0d0 00950 s%wscrit = 0.d0 00951 s%wb = 0.0d0 00952 s%pres = 0.0d0 00953 s%zi = 0.0d0 00954 s%wi = 0.0d0 00955 s%nuh = 0.0d0 00956 !s%cf = par%cf 00957 s%wm =0.d0 00958 s%zs0fac = 0.0d0 00959 00960 s%ududx = 0.d0 00961 s%vdvdy = 0.d0 00962 s%udvdx = 0.d0 00963 s%vdudy = 0.d0 00964 s%viscu = 0.d0 00965 s%viscv = 0.d0 00966 00967 s%taubx_add = 0.0d0 00968 s%tauby_add = 0.0d0 00969 00970 s%hstokes = 0.d0 00971 00972 if (par%wavemodel==WAVEMODEL_NONH) then 00973 s%breaking = 0 00974 endif 00975 00976 !if (par%vegetation==1) then 00977 ! s%vegtype = 0 00978 ! s%Cdveg = 0.d0 00979 ! s%ahveg = 0.d0 00980 ! s%bveg = 0.d0 00981 ! s%Nveg = 0.d0 00982 !endif 00983 00984 ! 00985 ! set-up tide and surge waterlevels 00986 call tide_init(s,par) 00987 00988 inquire(file=par%bedfricfile,exist=exists) 00989 if ((exists)) then 00990 open(723,file=par%bedfricfile) 00991 do j=1,s%ny+1 00992 read(723,*,iostat=ier)(s%bedfriccoef(i,j),i=1,s%nx+1) 00993 if (ier .ne. 0) then 00994 call report_file_read_error(par%bedfricfile) 00995 endif 00996 enddo 00997 close(723) 00998 else 00999 s%bedfriccoef=par%bedfriccoef 01000 endif 01001 01002 ! 01003 ! set zs, hh, wetu, wetv, wetz 01004 ! 01005 01006 s%zs0 = max(s%zs0,s%zb) 01007 ! cjaap: replaced par%hmin by par%eps 01008 if (par%hotstart .ne. 1) s%zs=max(s%zb+par%eps,s%zs0) 01009 s%hh=max(s%zs-s%zb,par%eps) 01010 01011 01012 !Initialize hu correctly to prevent spurious initial flow (Pieter) 01013 do j=1,s%ny+1 01014 do i=1,s%nx 01015 s%hu(i,j) = max(s%zs(i,j),s%zs(i+1,j))-max(s%zb(i,j),s%zb(i+1,j)) 01016 enddo 01017 enddo 01018 s%hu(s%nx+1,:)=s%hu(s%nx,:) 01019 01020 if (s%ny>0) then 01021 do j=1,s%ny 01022 do i=1,s%nx+1 01023 s%hv(i,j) = max(s%zs(i,j),s%zs(i,j+1))-max(s%zb(i,j),s%zb(i,j+1)) 01024 enddo 01025 enddo 01026 s%hv(:,s%ny+1)=s%hv(:,s%ny) 01027 else 01028 s%hv(:,1)=s%zs(:,1)-s%zb(:,1) 01029 endif 01030 01031 s%hum(1:s%nx,:) = 0.5d0*(s%hh(1:s%nx,:)+s%hh(2:s%nx+1,:)) 01032 s%hum(s%nx+1,:)=s%hh(s%nx+1,:) 01033 01034 ! R+L: Why are these variable reinitialise? 01035 s%dzsdt=0.d0 01036 s%dzsdx=0.d0 01037 s%dzsdy=0.d0 01038 s%dzbdt=0.d0 01039 s%uu=0.d0 01040 s%u=0.d0 01041 s%vv=0.d0 01042 s%v=0.d0 01043 s%vu=0.d0 01044 s%uv=0.d0 01045 s%ueu=0.d0 01046 s%vev=0.d0 01047 s%qx=0.d0 01048 s%qy=0.d0 01049 s%sedero=0.d0 01050 s%vmagu=0.d0 01051 s%vmagv=0.d0 01052 s%vmageu=0.d0 01053 s%vmagev=0.d0 01054 s%taubx=0.d0 01055 s%tauby=0.d0 01056 s%maxzs=-999.d0 01057 s%minzs=999.d0 01058 where(s%zs>s%zb+par%eps) 01059 s%wetz=1 01060 s%wetu=1 01061 s%wetv=1 01062 s%wete=1 01063 elsewhere 01064 s%wetz=0 01065 s%wetu=0 01066 s%wetv=0 01067 s%wete=0 01068 endwhere 01069 ! 01070 ! Start with initial velocities based on balance between water level gradient 01071 ! and bed friction term 01072 if (par%hotstartflow==1) then 01073 call writelog('ls','','Calculating stationary flow field for initial condition') 01074 ! water level gradients 01075 do j=1,s%ny+1 01076 do i=2,s%nx 01077 s%dzsdx(i,j)=(s%zs0(i+1,j)-s%zs0(i,j))/s%dsu(i,j) 01078 end do 01079 end do 01080 do j=1,s%ny ! Dano need to get correct slope on boundary y=0 01081 do i=1,s%nx+1 01082 s%dzsdy(i,j)=(s%zs0(i,j+1)-s%zs0(i,j))/s%dnv(i,j) 01083 end do 01084 end do 01085 ! Water depth in u,v-points mean for wind forcing 01086 do j=1,s%ny+1 01087 do i=1,s%nx+1 !Ap 01088 s%hum(i,j)=max(.5d0*(s%hh(i,j)+s%hh(min(s%nx,i)+1,j)),par%eps) 01089 end do 01090 end do 01091 do j=1,s%ny+1 01092 do i=1,s%nx+1 01093 s%hvm(i,j)=max(.5d0*(s%hh(i,j)+s%hh(i,min(s%ny,j)+1)),par%eps) 01094 end do 01095 end do 01096 ! residual error 01097 flowerr = huge(0.d0) 01098 vmagvold = 0.d0 01099 vmaguold = 0.d0 01100 do while (flowerr > 0.00001d0) 01101 ! 01102 ! Balance in longshore 01103 vmagvold=0.5d0*(vmagvold+sqrt(s%uv**2+s%vv**2)) ! mean needed for convergence 01104 ! solve v-balance of pressure gradient, wind forcing and bed friction 01105 where (s%wetv==1) 01106 s%vv = s%hv/s%cf/max(vmagvold,0.000001d0) & 01107 *(-par%g*s%dzsdy+par%rhoa*par%Cd*s%windnv*sqrt(s%windsu**2+s%windnv**2)/(par%rho*s%hvm)) ! Kees: wind correction 01108 elsewhere 01109 s%vv = 0.d0 01110 endwhere 01111 ! update vmagev 01112 ! u velocity in v points 01113 if (s%ny>0) then 01114 s%uv(2:s%nx,1:s%ny)= .25d0*(s%uu(1:s%nx-1,1:s%ny)+s%uu(2:s%nx,1:s%ny)+ & 01115 s%uu(1:s%nx-1,2:s%ny+1)+s%uu(2:s%nx,2:s%ny+1)) 01116 ! boundaries? 01117 ! wwvv and what about uv(:,1) ? 01118 if(xmpi_isright) then 01119 s%uv(:,s%ny+1) = s%uv(:,s%ny) 01120 endif 01121 else 01122 s%uv(2:s%nx,1)= .5d0*(s%uu(1:s%nx-1,1)+s%uu(2:s%nx,1)) 01123 endif !s%ny>0 01124 s%vmagev = sqrt(s%uv**2+s%vv**2) 01125 ! 01126 ! Balance in cross shore 01127 vmaguold= 0.5d0*(vmaguold+sqrt(s%uu**2+s%vu**2)) ! mean needed for convergence 01128 ! Solve balance of forces 01129 where (s%wetu==1) 01130 s%uu = s%hu/s%cf/max(vmaguold,0.000001d0) & 01131 *(-par%g*s%dzsdx+par%rhoa*par%Cd*s%windsu*sqrt(s%windsu**2+s%windnv**2)/(par%rho*s%hum)) ! Kees: wind correction 01132 elsewhere 01133 s%uu = 0.d0 01134 endwhere 01135 ! update vmageu 01136 if (s%ny>0) then 01137 s%vu(1:s%nx,2:s%ny)= 0.25d0*(s%vv(1:s%nx,1:s%ny-1)+s%vv(1:s%nx,2:s%ny)+ & 01138 s%vv(2:s%nx+1,1:s%ny-1)+s%vv(2:s%nx+1,2:s%ny)) 01139 if(xmpi_isleft) then 01140 s%vu(:,1) = s%vu(:,2) 01141 endif 01142 if(xmpi_isright) then 01143 s%vu(:,s%ny+1) = s%vu(:,s%ny) 01144 endif 01145 else 01146 s%vu(1:s%nx,1)= 0.5d0*(s%vv(1:s%nx,1)+s%vv(2:s%nx+1,1)) 01147 endif !s%ny>0 01148 s%vmageu = sqrt(s%uu**2+s%vu**2) 01149 ! 01150 ! Check residual error 01151 flowerr = max(maxval(abs(s%vmagev-vmagvold)),maxval(abs(s%vmageu-vmaguold))) 01152 enddo 01153 ! calculate all derivatives 01154 s%ueu=s%uu 01155 s%vev=s%vv 01156 s%qx=s%uu*s%hu 01157 s%qy=s%vv*s%hv 01158 s%vmagu=s%vmageu 01159 s%vmagv=s%vmagev 01160 s%u(2:s%nx,:)=0.5d0*(s%uu(1:s%nx-1,:)+s%uu(2:s%nx,:)) 01161 if(xmpi_istop) then 01162 s%u(1,:)=s%uu(1,:) 01163 endif 01164 if(xmpi_isbot) then 01165 s%u(s%nx+1,:)=s%u(s%nx,:) 01166 endif 01167 if (s%ny>0) then 01168 s%v(:,2:s%ny)=0.5d0*(s%vv(:,1:s%ny-1)+s%vv(:,2:s%ny)) 01169 if(xmpi_isleft) then 01170 s%v(:,1)=s%vv(:,1) 01171 endif 01172 if(xmpi_isright) then 01173 s%v(:,s%ny+1)=s%v(:,s%ny) 01174 endif 01175 s%v(s%nx+1,:)=s%v(s%nx,:) 01176 else ! Dano 01177 s%v=s%vv 01178 endif !s%ny>0 01179 endif 01180 ! 01181 ! Initialize for tide instant boundary condition 01182 ! 01183 if (par%tidetype==TIDETYPE_INSTANT) then 01184 ! RJ: 22-09-2010 01185 ! Check for whole domain whether a grid cell should be associated with 01186 ! 1) offshore tide and surge 01187 ! 2) bay tide and surge 01188 ! 3) no tide and surge 01189 ! 4) weighted tide and surge (for completely wet arrays) 01190 ! relative weight of offshore boundary and bay boundary for each grid point is stored in zs0fac 01191 ! 01192 s%zs0fac = 0.d0 01193 do j = 1,s%ny+1 01194 offshoreregime = .true. 01195 indoff = s%nx+1 ! ind of last point (starting at offshore boundary) that should be associated with offshore boundary 01196 indbay = 1 ! ind of first point (starting at offshore boundary) that should be associated with bay boundary 01197 do i = 1,s%nx+1 01198 if (offshoreregime .and. s%wetz(i,j)==0) then 01199 indoff = max(i-1,1) 01200 offshoreregime = .false. 01201 endif 01202 if (s%wetz(i,j)==0 .and. s%wetz(min(i+1,s%nx+1),j)==1) then 01203 indbay = min(i+1,s%nx+1) 01204 endif 01205 enddo 01206 01207 if (indbay==1 .and. indoff==s%nx+1) then ! in case of completely wet arrays linear interpolation for s%zs0fac 01208 ! Dano: don't know how to fix this for curvilinear 01209 ! zs0fac(:,j,2) = (xz-xz(1))/(xz(s%nx+1)-xz(1)) 01210 ! zs0fac(:,j,1) = 1-zs0fac(:,j,2) 01211 else ! in all other cases we assume three regims offshore, dry and bay 01212 s%zs0fac(1:indoff,j,1) = 1.d0 01213 s%zs0fac(1:indoff,j,2) = 0.d0 01214 if (indbay > 1) then 01215 s%zs0fac(indoff+1:indbay-1,j,1) = 0.d0 01216 s%zs0fac(indbay:s%nx+1,j,1) = 0.d0 01217 s%zs0fac(indoff+1:indbay-1,j,2) = 0.d0 01218 s%zs0fac(indbay:s%nx+1,j,2) = 1.d0 01219 endif 01220 endif 01221 enddo 01222 01223 endif ! tidetype = instant water level boundary 01224 01225 end subroutine flow_init 01226 01227 01228 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01229 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01230 01231 01232 subroutine sed_init (s,par) 01233 use params, only: parameters 01234 use spaceparamsdef, only: spacepars 01235 use xmpi_module, only: xmaster 01236 use logging_module, only: writelog, report_file_read_error 01237 use typesandkinds, only: slen 01238 01239 IMPLICIT NONE 01240 01241 type(spacepars),target :: s 01242 type(parameters) :: par 01243 01244 integer :: i,j,m,jg,start,ier 01245 character(slen) :: fnameg 01246 character(len=4) :: tempc 01247 real*8 :: tempr 01248 01249 01250 if(.not. xmaster) return 01251 01252 allocate(s%ccg(1:s%nx+1,1:s%ny+1,par%ngd)) 01253 allocate(s%dcbdy(1:s%nx+1,1:s%ny+1)) 01254 allocate(s%dcbdx(1:s%nx+1,1:s%ny+1)) 01255 allocate(s%dcsdy(1:s%nx+1,1:s%ny+1)) 01256 allocate(s%dcsdx(1:s%nx+1,1:s%ny+1)) 01257 allocate(s%Tsg(1:s%nx+1,1:s%ny+1,par%ngd)) 01258 allocate(s%Susg(1:s%nx+1,1:s%ny+1,par%ngd)) 01259 allocate(s%Svsg(1:s%nx+1,1:s%ny+1,par%ngd)) 01260 allocate(s%Subg(1:s%nx+1,1:s%ny+1,par%ngd)) 01261 allocate(s%Svbg(1:s%nx+1,1:s%ny+1,par%ngd)) 01262 allocate(s%vmag(1:s%nx+1,1:s%ny+1)) 01263 allocate(s%ceqsg(1:s%nx+1,1:s%ny+1,par%ngd)) 01264 allocate(s%ceqbg(1:s%nx+1,1:s%ny+1,par%ngd)) 01265 !if (par%dilatancy==1) allocate(s%D15(1:par%ngd)) ! Lodewijk 01266 !if (par%dilatancy==1) then 01267 ! allocate(s%D15(1:par%ngd)) ! Lodewijk 01268 !else 01269 ! allocate(s%D15(1)) ! wwvv: s%D15 will be distributed, so it must exist 01270 !endif 01271 allocate(s%D15(1:par%ngd)) ! Robert: this is distributed according to stencil in spaceparams.tmpl 01272 allocate(s%D50(1:par%ngd)) 01273 allocate(s%D90(1:par%ngd)) 01274 allocate(s%D50top(1:s%nx+1,1:s%ny+1)) 01275 allocate(s%D90top(1:s%nx+1,1:s%ny+1)) 01276 allocate(s%sedcal(1:par%ngd)) 01277 allocate(s%ucrcal(1:par%ngd)) 01278 allocate(s%nd(1:s%nx+1,1:s%ny+1)) 01279 allocate(s%dzbed(1:s%nx+1,1:s%ny+1,1:par%nd)) 01280 allocate(s%pbbed(1:s%nx+1,1:s%ny+1,1:par%nd,1:par%ngd)) 01281 allocate(s%z0bed(1:s%nx+1,1:s%ny+1)) 01282 allocate(s%ureps(1:s%nx+1,1:s%ny+1)) 01283 allocate(s%urepb(1:s%nx+1,1:s%ny+1)) 01284 allocate(s%vreps(1:s%nx+1,1:s%ny+1)) 01285 allocate(s%vrepb(1:s%nx+1,1:s%ny+1)) 01286 allocate(s%ero(1:s%nx+1,1:s%ny+1,1:par%ngd)) 01287 allocate(s%depo_ex(1:s%nx+1,1:s%ny+1,1:par%ngd)) 01288 allocate(s%depo_im(1:s%nx+1,1:s%ny+1,1:par%ngd)) 01289 allocate(s%kb(1:s%nx+1,1:s%ny+1)) 01290 allocate(s%Tbore(1:s%nx+1,1:s%ny+1)) 01291 allocate(s%ua(1:s%nx+1,1:s%ny+1)) 01292 allocate(s%dzav(1:s%nx+1,1:s%ny+1)) 01293 allocate(s%Sk(1:s%nx+1,1:s%ny+1)) 01294 allocate(s%As(1:s%nx+1,1:s%ny+1)) 01295 allocate(s%kturb(1:s%nx+1,1:s%ny+1)) 01296 allocate(s%rolthick(1:s%nx+1,1:s%ny+1)) 01297 allocate(s%Sutot(1:s%nx+1,1:s%ny+1)) ! Only really for easy output 01298 allocate(s%Svtot(1:s%nx+1,1:s%ny+1)) ! Only really for easy output 01299 allocate(s%cctot(1:s%nx+1,1:s%ny+1)) ! Only really for easy output 01300 allocate(s%runup(1:s%ny+1)) 01301 allocate(s%Hrunup(1:s%ny+1)) 01302 allocate(s%xHrunup(1:s%ny+1)) 01303 allocate(s%istruct(1:s%ny+1)) 01304 allocate(s%iwl(1:s%ny+1)) 01305 allocate(s%strucslope(1:s%ny+1)) 01306 allocate(s%Dc(1:s%nx+1,1:s%ny+1)) 01307 allocate(s%ccz(1:s%nx+1,1:s%ny+1,1:par%kmax)) ! Q3D: concentration profile 01308 allocate(s%ca(1:s%nx+1,1:s%ny+1)) ! Q3D: reference concentration 01309 allocate(s%refA(1:s%nx+1,1:s%ny+1)) ! Q3D: reference level 01310 01311 ! Initialize so structures can be implemented more easily 01312 s%pbbed = 0.d0 01313 ! 01314 ! Set grain size(s) 01315 ! 01316 if (par%dilatancy==1) s%D15 = par%D15(1:par%ngd) 01317 s%D50 = par%D50(1:par%ngd) 01318 s%D90 = par%D90(1:par%ngd) 01319 s%sedcal = par%sedcal(1:par%ngd) 01320 s%ucrcal = par%ucrcal(1:par%ngd) 01321 01322 01323 if (par%ngd==1) then 01324 01325 ! No multi sediment, but we do need some data to keep the script running 01326 01327 s%pbbed(:,:,:,1)=1.d0 ! set sand fraction everywhere, not structure fraction (if exist) which is still 0.d0 01328 par%nd_var=2 01329 01330 s%dzbed(:,:,1:par%nd_var-1) = max(par%dzg1,10.d0) 01331 s%dzbed(:,:,par%nd_var) = max(par%dzg2,10.d0) 01332 s%dzbed(:,:,par%nd_var+1:par%nd) = max(par%dzg3,10.d0) 01333 01334 else 01335 01336 ! Fill pbed en dzbed 01337 ! 01338 do jg=1,par%ngd 01339 write(tempc,'(i4)')jg 01340 start=4-floor(log10(real(jg))) 01341 write(fnameg,'(a,a,a)')'gdist',tempc(start:4),'.inp' 01342 open(31,file=fnameg) 01343 do m=1,par%nd 01344 do j=1,s%ny+1 01345 read(31,*,iostat=ier)(s%pbbed(i,j,m,jg),i=1,s%nx+1) 01346 if (ier .ne. 0) then 01347 call report_file_read_error(fnameg) 01348 endif 01349 enddo 01350 enddo 01351 close(31) 01352 enddo 01353 ! Rework pbbed so that sum fractions = 1 01354 do m=1,par%nd 01355 do j=1,s%ny+1 !Jaap instead of 2:s%ny 01356 do i=1,s%nx+1 !Jaap instead of 2:s%nx 01357 01358 tempr=sum(s%pbbed(i,j,m,1:par%ngd)) 01359 if (abs(1.d0-tempr)>0.000001d0) then 01360 ! Maybe fix this warning if in combination with structures 01361 call writelog('lws','(a,i0,a,i0,a,i0,a)',' Warning: Resetting sum of sediment fractions in point (',& 01362 i,',',j,') layer ,',m,& 01363 ' to equal unity.') 01364 if (tempr<=tiny(0.d0)) then ! In case cell has zero sediment (i.e. only hard structure) 01365 s%pbbed(i,j,m,:)=1.d0/dble(par%ngd) 01366 else 01367 s%pbbed(i,j,m,:)=s%pbbed(i,j,m,:)/tempr 01368 endif 01369 endif 01370 enddo 01371 enddo 01372 enddo 01373 ! boundary neumann --> Jaap not necessary already done in loop above 01374 01375 ! sediment thickness 01376 s%dzbed(:,:,1:par%nd_var-1) = par%dzg1 01377 s%dzbed(:,:,par%nd_var) = par%dzg2 01378 s%dzbed(:,:,par%nd_var+1:par%nd) = par%dzg3 01379 endif 01380 01381 ! Initialize representative sed.diameter at the bed for flow friction and output 01382 do j=1,s%ny+1 01383 do i=1,s%nx+1 01384 s%D50top(i,j) = sum(s%pbbed(i,j,1,:)*s%D50) 01385 s%D90top(i,j) = sum(s%pbbed(i,j,1,:)*s%D90) 01386 enddo 01387 enddo 01388 ! 01389 ! Set non-erodable layer 01390 ! 01391 allocate(s%structdepth(s%nx+1,s%ny+1)) 01392 01393 s%structdepth = 100.d0 01394 01395 if (par%struct==1 .and. par%hotstart==0) then 01396 !call readkey('params.txt','ne_layer',fnameh) 01397 !open(31,file=fnameh) 01398 open(31,file=par%ne_layer) 01399 01400 do j=1,s%ny+1 01401 read(31,*,iostat=ier)(s%structdepth(i,j),i=1,s%nx+1) 01402 if (ier .ne. 0) then 01403 call report_file_read_error(par%ne_layer) 01404 endif 01405 end do 01406 01407 close(31) 01408 01409 endif 01410 01411 ! bottom of sediment model 01412 s%z0bed = s%zb - sum(s%dzbed,DIM=3) 01413 01414 s%nd = par%nd 01415 01416 s%ureps = 0.d0 01417 s%vreps = 0.d0 01418 s%ccg = 0.d0 01419 s%ceqbg = 0.d0 01420 s%ceqsg = 0.d0 01421 s%Susg = 0.d0 01422 s%Svsg = 0.d0 01423 s%Subg = 0.d0 01424 s%Svbg = 0.d0 01425 s%dcsdx = 0.d0 01426 s%dcsdy = 0.d0 01427 s%dcbdx = 0.d0 01428 s%dcbdy = 0.d0 01429 s%ero = 0.d0 01430 s%depo_im = 0.d0 01431 s%depo_ex = 0.d0 01432 s%kb = 0.d0 01433 s%Tbore = 0.d0 01434 s%ua = 0.d0 01435 s%dzav = 0.d0 01436 s%Sk = 0.d0 01437 s%As = 0.d0 01438 s%kturb = 0.d0 01439 s%rolthick = 0.d0 01440 s%Sutot = 0.d0 01441 s%Svtot = 0.d0 01442 s%cctot = 0.d0 01443 s%runup = 0.d0 01444 s%Hrunup = 0.d0 01445 s%xHrunup = 0.d0 01446 s%istruct = s%nx+1 01447 s%strucslope = 0.d0 01448 s%Dc = 0.d0 01449 01450 ! Initialize dzbdx, dzbdy 01451 do j=1,s%ny+1 01452 do i=1,s%nx 01453 s%dzbdx(i,j)=(s%zb(i+1,j)-s%zb(i,j))/s%dsu(i,j) 01454 enddo 01455 enddo 01456 ! dummy, needed to keep compiler happy 01457 s%dzbdx(s%nx+1,:)=s%dzbdx(s%nx,:) 01458 01459 if (s%ny>0) then 01460 do j=1,s%ny 01461 do i=1,s%nx+1 01462 s%dzbdy(i,j)=(s%zb(i,j+1)-s%zb(i,j))/s%dnv(i,j) 01463 enddo 01464 enddo 01465 s%dzbdy(:,s%ny+1)=s%dzbdy(:,s%ny) 01466 else 01467 s%dzbdy=0.d0 01468 endif 01469 01470 end subroutine sed_init 01471 01472 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01473 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01474 01475 subroutine discharge_init(s, par) 01476 01477 use params, only: parameters 01478 use spaceparamsdef, only: spacepars 01479 use xmpi_module, only: xmaster 01480 use logging_module, only: report_file_read_error 01481 01482 01483 implicit none 01484 01485 type(spacepars),target :: s 01486 type(parameters) :: par 01487 01488 integer :: i,j 01489 integer :: io 01490 integer :: m1,m2,n1,n2 01491 real*8, dimension(:),allocatable :: xdb,ydb,xde,yde 01492 integer,dimension(2) :: mnb,mne 01493 01494 01495 if(.not. xmaster) return 01496 01497 io = 0 01498 01499 allocate(xdb (par%ndischarge) ) 01500 allocate(ydb (par%ndischarge) ) 01501 allocate(xde (par%ndischarge) ) 01502 allocate(yde (par%ndischarge) ) 01503 01504 allocate(s%pntdisch (1:par%ndischarge) ) 01505 allocate(s%pdisch (1:par%ndischarge , 1:4) ) 01506 allocate(s%tdisch (1:par%ntdischarge) ) 01507 allocate(s%qdisch (1:par%ntdischarge , 1:par%ndischarge) ) 01508 01509 s%pntdisch = 0 01510 s%tdisch = 0.d0 01511 s%pdisch = 0 01512 s%qdisch = 0.d0 01513 01514 if (par%ndischarge>0) then 01515 01516 ! read discharge locations 01517 open(10,file=par%disch_loc_file) 01518 do i=1,par%ndischarge 01519 read(10,*,IOSTAT=io) xdb(i),ydb(i),xde(i),yde(i) 01520 if (io .ne. 0) then 01521 call report_file_read_error(par%disch_loc_file) 01522 endif 01523 ! distinguish between horizontal and vertical discharge 01524 if (xdb(i).eq.xde(i) .and. ydb(i).eq.yde(i)) then 01525 s%pntdisch(i) = 1 01526 else 01527 s%pntdisch(i) = 0 01528 endif 01529 01530 enddo 01531 close(10) 01532 01533 if (par%ntdischarge>0) then 01534 01535 ! read time series 01536 open(10,file=par%disch_timeseries_file) 01537 do i=1,par%ntdischarge 01538 read(10,*,IOSTAT=io) s%tdisch(i),(s%qdisch(i,j),j=1,par%ndischarge) 01539 if (io .ne. 0) then 01540 call report_file_read_error(par%disch_timeseries_file) 01541 endif 01542 enddo 01543 close(10) 01544 endif 01545 endif 01546 01547 ! initialise each discharge location 01548 do i=1,par%ndischarge 01549 01550 ! dxd = abs(xde(i)-xdb(i)) 01551 ! dyd = abs(yde(i)-ydb(i)) 01552 mnb = minloc(sqrt((s%xz-xdb(i))**2+(s%yz-ydb(i))**2)) 01553 mne = minloc(sqrt((s%xz-xde(i))**2+(s%yz-yde(i))**2)) 01554 01555 ! convert discharge location to cell indices depending on type of discharge: 01556 ! point discharge, in v-direction or in u-direction 01557 if (s%pntdisch(i).eq.1) then 01558 01559 ! point discharge (no orientation, no added momentum, just mass) 01560 01561 ! mnb = minloc(sqrt((s%xz-xdb(i))**2+(s%yz-ydb(i))**2)) 01562 ! mne = minloc(sqrt((s%xz-xde(i))**2+(s%yz-yde(i))**2)) 01563 01564 s%pdisch(i,:) = (/mnb(1),mnb(2),0,0/) 01565 ! elseif (dxd.gt.dyd) then 01566 elseif (mnb(1).ne.mne(1)) then 01567 01568 ! discharge through v-points 01569 01570 ! mnb = minloc(sqrt((s%xv-xdb(i))**2+(s%yv-ydb(i))**2)) 01571 ! mne = minloc(sqrt((s%xv-xde(i))**2+(s%yv-yde(i))**2)) 01572 01573 m1 = minval((/mnb(1),mne(1)/)) 01574 m2 = maxval((/mnb(1),mne(1)/)) 01575 n1 = nint(0.5*(mnb(2)+mne(2))) 01576 01577 if (n1.lt.1) n1 = 1 01578 if (n1.gt.s%ny) n1 = s%ny 01579 01580 s%pdisch(i,:) = (/m1,n1,m2,n1/) 01581 else 01582 01583 ! discharge through u-points 01584 01585 ! mnb = minloc(sqrt((s%xu-xdb(i))**2+(s%yu-ydb(i))**2)) 01586 ! mne = minloc(sqrt((s%xu-xde(i))**2+(s%yu-yde(i))**2)) 01587 01588 m1 = nint(0.5*(mnb(1)+mne(1))) 01589 n1 = minval((/mnb(2),mne(2)/)) 01590 n2 = maxval((/mnb(2),mne(2)/)) 01591 01592 if (m1.lt.1) m1 = 1 01593 if (m1.gt.s%nx) m1 = s%nx 01594 01595 s%pdisch(i,:) = (/m1,n1,m1,n2/) 01596 endif 01597 enddo 01598 01599 ! incorporate morfac 01600 if (par%morfacopt == 1) then 01601 s%tdisch = s%tdisch/max(par%morfac,1.d0) 01602 endif 01603 end subroutine discharge_init 01604 01605 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01606 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 01607 01608 subroutine drifter_init(s, par) 01609 01610 use params, only: parameters 01611 use spaceparamsdef, only: spacepars 01612 use logging_module, only: report_file_read_error 01613 use xmpi_module 01614 01615 implicit none 01616 01617 type(spacepars),target :: s 01618 type(parameters) :: par 01619 01620 integer :: i,ier 01621 real*8 :: xdrift,ydrift 01622 real*8 :: ds,dn 01623 integer,dimension(2) :: mn 01624 01625 01626 if(.not. (xmaster .or. xomaster)) return 01627 01628 allocate(s%idrift (par%ndrifter)) 01629 allocate(s%jdrift (par%ndrifter)) 01630 allocate(s%tdriftb (par%ndrifter)) 01631 allocate(s%tdrifte (par%ndrifter)) 01632 01633 if (par%ndrifter>0) then 01634 01635 ! read drifter file 01636 if(xmaster) then 01637 open(10,file=par%drifterfile) 01638 do i=1,par%ndrifter 01639 read(10,*,iostat=ier)xdrift,ydrift,s%tdriftb(i),s%tdrifte(i) 01640 if (ier .ne. 0) then 01641 call report_file_read_error(par%drifterfile) 01642 endif 01643 01644 mn = minloc(sqrt((s%xz-xdrift)**2+(s%yz-ydrift)**2)) 01645 01646 ds = (xdrift - s%xz(mn(1),mn(2)))*cos(s%alfaz(mn(1),mn(2))) & 01647 +(ydrift - s%yz(mn(1),mn(2)))*sin(s%alfaz(mn(1),mn(2))) 01648 dn = -(xdrift - s%xz(mn(1),mn(2)))*sin(s%alfaz(mn(1),mn(2))) & 01649 +(ydrift - s%yz(mn(1),mn(2)))*cos(s%alfaz(mn(1),mn(2))) 01650 01651 s%idrift(i) = mn(1) + ds/s%dsu(mn(1),mn(2)) 01652 s%jdrift(i) = mn(2) + dn/s%dnv(mn(1),mn(2)) 01653 enddo 01654 close(10) 01655 01656 ! incorporate morfac 01657 if (par%morfacopt == 1) then 01658 s%tdriftb = s%tdriftb/max(par%morfac,1.d0) 01659 s%tdrifte = s%tdrifte/max(par%morfac,1.d0) 01660 endif 01661 endif 01662 01663 #ifdef USEMPI 01664 call xmpi_send(xmpi_imaster, xmpi_omaster,s%idrift) 01665 call xmpi_send(xmpi_imaster, xmpi_omaster,s%jdrift) 01666 call xmpi_send(xmpi_imaster, xmpi_omaster,s%tdriftb) 01667 call xmpi_send(xmpi_imaster, xmpi_omaster,s%tdrifte) 01668 #endif 01669 01670 endif 01671 end subroutine drifter_init 01672 01673 subroutine hotstart_init_1(s,par) 01674 use params 01675 use spaceparams 01676 use filefunctions 01677 use logging_module 01678 use paramsconst 01679 01680 implicit none 01681 01682 type(spacepars),target :: s 01683 type(parameters) :: par 01684 01685 integer :: i,j,ier 01686 01687 allocate(s%zs(1:s%nx+1,1:s%ny+1)) 01688 allocate(s%uu(1:s%nx+1,1:s%ny+1)) 01689 allocate(s%vv(1:s%nx+1,1:s%ny+1)) 01690 ! Robert: not zb, that is allocated in grid_bathy called earlier 01691 01692 ! all models 01693 call lowlevel_read_hotstart_2d(s%zs,s%nx,s%ny,'zs',par%hotstartfileno) 01694 call lowlevel_read_hotstart_2d(s%zb,s%nx,s%ny,'zb',par%hotstartfileno) 01695 call lowlevel_read_hotstart_2d(s%uu,s%nx,s%ny,'uu',par%hotstartfileno) 01696 call lowlevel_read_hotstart_2d(s%vv,s%nx,s%ny,'vv',par%hotstartfileno) 01697 s%zb0 = s%zb 01698 end subroutine hotstart_init_1 01699 01700 subroutine hotstart_init_2(s,par) 01701 use params 01702 use spaceparams 01703 use filefunctions 01704 use logging_module 01705 use paramsconst 01706 01707 implicit none 01708 01709 type(spacepars),target :: s 01710 type(parameters) :: par 01711 01712 integer :: i,j,ier 01713 ! groundwater parameters 01714 if(par%gwflow==1) then 01715 call lowlevel_read_hotstart_2d(s%gwlevel,s%nx,s%ny,'gwlevel',par%hotstartfileno) 01716 call lowlevel_read_hotstart_2d(s%gwhead,s%nx,s%ny,'gwhead',par%hotstartfileno) 01717 if(par%gwnonh==1) call lowlevel_read_hotstart_2d(s%gwcurv,s%nx,s%ny,'gwcurv',par%hotstartfileno) 01718 endif 01719 01720 ! wave energy 01721 if (par%swave==1) then 01722 call lowlevel_read_hotstart_3d(s%ee,s%nx,s%ny,s%ntheta,'ee',par%hotstartfileno) 01723 call lowlevel_read_hotstart_3d(s%rr,s%nx,s%ny,s%ntheta,'rr',par%hotstartfileno) 01724 endif 01725 01726 ! single-dir 01727 if (par%single_dir==1) then 01728 call lowlevel_read_hotstart_3d(s%ee_s,s%nx,s%ny,s%ntheta_s,'ee_s',par%hotstartfileno) 01729 endif 01730 01731 ! nonh parameters 01732 if (par%wavemodel == WAVEMODEL_NONH) then 01733 call lowlevel_read_hotstart_2d_int(s%breaking,s%nx,s%ny,'breaking',par%hotstartfileno) 01734 call lowlevel_read_hotstart_2d(s%wb,s%nx,s%ny,'wb',par%hotstartfileno) 01735 call lowlevel_read_hotstart_2d(s%ws,s%nx,s%ny,'ws',par%hotstartfileno) 01736 endif 01737 01738 ! sediment transport 01739 if (par%sedtrans==1 .and. (par%sus==1 .or. par%bulk==1)) then 01740 call lowlevel_read_hotstart_3d(s%ccg,s%nx,s%ny,par%ngd,'ccg',par%hotstartfileno) 01741 endif 01742 01743 ! turbulence 01744 if (par%turb==1) then 01745 call lowlevel_read_hotstart_2d(s%kturb,s%nx,s%ny,'kturb',par%hotstartfileno) 01746 endif 01747 01748 ! structures 01749 if(par%struct==1) then 01750 call lowlevel_read_hotstart_2d(s%structdepth,s%nx,s%ny,'structdepth',par%hotstartfileno) 01751 endif 01752 01753 ! nhplus 01754 if (par%nonhq3d==1) then 01755 call lowlevel_read_hotstart_2d(s%dU,s%nx,s%ny,'dU',par%hotstartfileno) 01756 call lowlevel_read_hotstart_2d(s%dV,s%nx,s%ny,'dV',par%hotstartfileno) 01757 endif 01758 01759 ! umean and vmean 01760 if ((par%flow==1).or.(par%wavemodel==WAVEMODEL_NONH)) then 01761 call lowlevel_read_hotstart_2d(s%umean,s%nx,s%ny,'umean',par%hotstartfileno) 01762 call lowlevel_read_hotstart_2d(s%vmean,s%nx,s%ny,'vmean',par%hotstartfileno) 01763 endif 01764 01765 end subroutine hotstart_init_2 01766 01767 subroutine lowlevel_read_hotstart_2d(var,nx,ny,varname,ith) 01768 use filefunctions 01769 01770 integer,intent(in) :: nx,ny,ith 01771 real*8,dimension(nx+1,ny+1) :: var 01772 character(*), intent(in) :: varname 01773 01774 integer :: unit,i,j,ierr 01775 character(128) :: fname,rowfmt 01776 01777 write(fname,'(a,a,I6.6,a)') 'hotstart_',trim(varname),ith,'.dat' 01778 call check_file_exist(fname) 01779 call check_file_length(fname,nx+1,ny+1) 01780 01781 unit = create_new_fid() 01782 write(rowfmt,'(a,I12,a)') '(',ny+1,'(1X,D))' 01783 open(unit,file=fname,form='FORMATTED', iostat=ierr,status='OLD',action='READ') 01784 do j=1,ny+1 01785 read(unit,FMT=rowfmt,iostat=ierr) (var(i,j), i=1,nx+1) 01786 enddo 01787 close(unit) 01788 01789 end subroutine lowlevel_read_hotstart_2d 01790 01791 subroutine lowlevel_read_hotstart_2d_int(var,nx,ny,varname,ith) 01792 use filefunctions 01793 01794 integer,intent(in) :: nx,ny,ith 01795 integer,dimension(nx+1,ny+1) :: var 01796 character(*), intent(in) :: varname 01797 01798 integer :: unit,i,j,ierr 01799 character(128) :: fname,rowfmt 01800 01801 write(fname,'(a,a,I6.6,a)') 'hotstart_',trim(varname),ith,'.dat' 01802 call check_file_exist(fname) 01803 call check_file_length(fname,nx+1,ny+1) 01804 01805 unit = create_new_fid() 01806 write(rowfmt,'(a,I12,a)') '(',ny+1,'(1X,I4))' 01807 open(unit,file=fname,form='FORMATTED', iostat=ierr) 01808 do j=1,ny+1 01809 read(unit,FMT=rowfmt,iostat=ierr) (var(i,j), i=1,nx+1) 01810 enddo 01811 close(unit) 01812 01813 end subroutine lowlevel_read_hotstart_2d_int 01814 01815 subroutine lowlevel_read_hotstart_3d(var,nx,ny,nd,varname,ith) 01816 use filefunctions 01817 01818 integer,intent(in) :: nx,ny,nd,ith 01819 real*8,dimension(nx+1,ny+1,nd) :: var 01820 character(*), intent(in) :: varname 01821 01822 integer :: unit,i,j,k,ierr 01823 character(128) :: fname,rowfmt 01824 01825 write(fname,'(a,a,I6.6,a)') 'hotstart_',trim(varname),ith,'.dat' 01826 call check_file_exist(fname) 01827 call check_file_length(fname,nx+1,ny+1,nd) 01828 01829 unit = create_new_fid() 01830 write(rowfmt,'(a,I12,a)') '(',ny+1,'(1X,D))' 01831 open(unit,file=fname,form='FORMATTED', iostat=ierr) 01832 do k=1,nd 01833 do j=1,ny+1 01834 read(unit,FMT=rowfmt,iostat=ierr) (var(i,j,k), i=1,nx+1) 01835 enddo 01836 enddo 01837 close(unit) 01838 01839 end subroutine lowlevel_read_hotstart_3d 01840 01841 end module initialize_module
1.7.6.1
