module spaceparams type spacepars real*8,dimension(:,:),allocatable :: x ! [m] x-coord. comp. grid (positive shoreward, perp. to coastline) real*8,dimension(:,:),allocatable :: y ! [m] y-coord. comp. grid real*8,dimension(:),allocatable :: xz ! [m] x-coord. comp. grid (positive shoreward, perp. to coastline) real*8,dimension(:),allocatable :: yz ! [m] y-coord. comp. grid real*8,dimension(:),allocatable :: xu ! [m] x-coord. in u points real*8,dimension(:),allocatable :: yv ! [m] y-coord. in v points real*8,dimension(:,:),allocatable :: xw ! [m] world x-coordinates real*8,dimension(:,:),allocatable :: yw ! [m] world y-coordinates real*8 :: dx !`[m] grid size x-direction real*8 :: dy ! [m] grid size y-direction real*8 :: xori ! [m] x-origin of grid in world coordinates real*8 :: yori ! [m] y-origin of grid in world coordinates real*8 :: alfa ! [rad] (deg on input) angle of grid w.r.t. East real*8 :: posdwn ! [-] depths defined positive downwards (1) or upwards(-1) integer :: nx ! [-] number of grid cells x-direction integer :: ny ! [-] number of grid cells y-direction real*8,dimension(:,:),allocatable :: zb ! [m] bed level real*8,dimension(:,:),allocatable :: zb0 ! [m] initial bed level real*8,dimension(:),allocatable :: theta ! [rad] wave angles directional distribution ! w.r.t. comp. x-axis integer :: ntheta ! [-] number of wave direction bins real*8 :: dtheta ! [rad] wave direction bin size real*8 :: theta0 ! [rad] mean incident wave angle real*8,dimension(:),allocatable :: cxsth ! [-] cos(theta) real*8,dimension(:),allocatable :: sxnth ! [-] sin(theta) real*8,dimension(:,:),allocatable :: thetamean ! [rad] mean wave angle real*8,dimension(:,:),allocatable :: Fx ! [N/m2] wave force x-direction real*8,dimension(:,:),allocatable :: Fy ! [N/m2] wave force y-direction real*8,dimension(:,:),allocatable :: Sxy ! [N/m] radiation stress real*8,dimension(:,:),allocatable :: Syy ! [N/m] radiation stress real*8,dimension(:,:),allocatable :: Sxx ! [N/m] radiation stress real*8,dimension(:,:),allocatable :: n ! [-] ratio group velocity/wave celerity real*8,dimension(:,:),allocatable :: H ! [m] wave height real*8,dimension(:,:,:),allocatable :: cgx ! [m/s] group velocity x-direction real*8,dimension(:,:,:),allocatable :: cgy ! [m/s] group velocity y-direction real*8,dimension(:,:,:),allocatable :: cx ! [m/s] wave celerity x-direction real*8,dimension(:,:,:),allocatable :: cy ! [m/s] wave celerity y-direction real*8,dimension(:,:,:),allocatable :: ctheta ! [rad/s] wave celerity theta-direction (refraction) real*8,dimension(:,:,:),allocatable :: ee ! [J/m2/rad] directionally distributed wave energy real*8,dimension(:,:,:),allocatable :: thet ! [rad] wave angles real*8,dimension(:,:,:),allocatable :: costhet ! [-] cos of wave angles real*8,dimension(:,:,:),allocatable :: sinthet ! [-] sin of wave angles real*8,dimension(:,:,:),allocatable :: sigt ! [rad/s] relative frequency real*8,dimension(:,:,:),allocatable :: rr ! [J/m2/rad] directionally distributed roller energy real*8,dimension(:,:),allocatable :: k ! [rad/m] wave number real*8,dimension(:,:),allocatable :: c ! [m/s] wave celerity real*8,dimension(:,:),allocatable :: cg ! [m/s] group velocity real*8,dimension(:,:),allocatable :: sigm ! [rad/s] mean frequency real*8,dimension(:,:),allocatable :: hh ! [m] water depth real*8,dimension(:,:),allocatable :: zs ! [m] water level real*8,dimension(:,:),allocatable :: zs0 ! [m] water level due to tide alone real*8,dimension(:),allocatable :: tideinpt ! [s] input time of input tidal signal real*8,dimension(:,:),allocatable :: tideinpz ! [m] input tidal signal real*8,dimension(:,:),allocatable :: dzsdt ! [m/s] rate of change water level real*8,dimension(:,:),allocatable :: dzbdt ! [m/s] rate of change bed level real*8,dimension(:,:),allocatable :: uu ! [m/s] (GLM) x-velocity in u-points real*8,dimension(:,:),allocatable :: vv ! [m/s] (GLM) y-velocity in v-points real*8,dimension(:,:),allocatable :: qx ! [m2/s] x-discharge in u-points real*8,dimension(:,:),allocatable :: qy ! [m2/s] y-discharge in u-points real*8,dimension(:,:),allocatable :: sedero ! [m] cum. sedimentation/erosion real*8,dimension(:,:),allocatable :: dcdx ! [kg/m3/m] concentration gradient x-dir. real*8,dimension(:,:),allocatable :: dcdy ! [kg/m3/m] concentration gradient y-dir. real*8,dimension(:,:),allocatable :: ui ! [m/s] incident bound wave velocity real*8,dimension(:,:),allocatable :: E ! [Nm/m2] wave energy real*8,dimension(:,:),allocatable :: R ! [Nm/m2] roller energy real*8,dimension(:,:),allocatable :: urms ! [m/s] orbital velocity real*8,dimension(:,:),allocatable :: D ! [W/m2] dissipation real*8,dimension(:,:),allocatable :: ust ! [m/s] Stokes drift real*8,dimension(:,:),allocatable :: tm ! [rad] mean wave direction real*8,dimension(:,:),allocatable :: ueu ! [m/s] Eulerian mean velocity x-dir. real*8,dimension(:,:),allocatable :: vev ! [m/s] Eulerian mean velocity y-dir. real*8,dimension(:,:),allocatable :: vmagu ! [m/s] (GLM) velocity magnitude u-points real*8,dimension(:,:),allocatable :: vmageu ! [m/s] (GLM) velocity magnitude u-points real*8,dimension(:,:),allocatable :: vmagv ! [m/s] (GLM) velocity magnitude v-points real*8,dimension(:,:),allocatable :: vmagev ! [m/s] (GLM) velocity magnitude v-points real*8,dimension(:,:),allocatable :: u ! [m/s] (GLM) x-velocity cell centre (for output) real*8,dimension(:,:),allocatable :: v ! [m/s] (GLM) y-velocity cell centre (for output) real*8,dimension(:,:),allocatable :: ue ! [m/s] Eulerian mean x-velocity cell centre (for output) real*8,dimension(:,:),allocatable :: ve ! [m/s] Eulerian mean y-velocity cell centre (for output) real*8,dimension(:,:),allocatable :: hold ! [m] water depth previous time step integer,dimension(:,:),allocatable :: wetu ! [-] mask wet/dry u-points integer,dimension(:,:),allocatable :: wetv ! [-] mask wet/dry v-points integer,dimension(:,:),allocatable :: wetz ! [-] mask wet/dry eta-points real*8,dimension(:,:),allocatable :: hu ! [m] water depth in u-points real*8,dimension(:,:),allocatable :: hv ! [m] water depth in v-points real*8,dimension(:,:),allocatable :: hum ! [m] water depth in u-points real*8,dimension(:,:),allocatable :: hvm ! [m] water depth in v-points !real*8,dimension(:,:),allocatable :: ceq ! [m3/m3] depth-averaged equilibrium concentration real*8,dimension(:,:),allocatable :: vmag ! [m/s] velocity magnitude in cell centre !real*8,dimension(:,:),allocatable :: Su ! [m2/s] sediment transport x-dir. (excluding pores) !real*8,dimension(:,:),allocatable :: Sv ! [m2/s] sediment transport y-dir. (excluding pores) !real*8,dimension(:,:),allocatable :: Ts ! [s] adaptation time scale !real*8,dimension(:,:),allocatable :: cc ! [m3/m3] depth-averaged concentration real*8,dimension(:,:,:),allocatable :: ccg ! [m3/m3] depth-averaged concentration for each sediment fraction real*8,dimension(:,:),allocatable :: uwf ! [m/s] x-comp. Stokes drift real*8,dimension(:,:),allocatable :: vwf ! [m/s] y-comp. Stokes drift real*8,dimension(:,:),allocatable :: ustr ! [m/s] return flow due to roller real*8,dimension(:,:),allocatable :: usd ! [m/s] return flow due to roller after breaker delay real*8,dimension(:),allocatable :: bi ! [m] incoming bound long wave real*8,dimension(:,:),allocatable :: DR ! [W/m2] roller energy dissipation real*8,dimension(:,:),allocatable :: umean ! [m/s] longterm mean velocity at bnds integer :: vardx ! [-] 0 = uniform grid size, 1 = variable grid size real*8,dimension(:,:),allocatable :: vu ! [m/s] y velocity in u points real*8,dimension(:,:),allocatable :: uv ! [m/s] x velocity in v points real*8,dimension(:,:,:,:),allocatable :: graindistr ! [-] fractional graindistribution for sediment classes real*8,dimension(:),allocatable :: D50 ! [m] D50 grain diameters for all sediment classses real*8,dimension(:),allocatable :: D90 ! [m] D90 grain diameters for all sediment classses real*8,dimension(:),allocatable :: sedcal ! [-] equilibrium sediment concentartion factor for each sediment class real*8,dimension(:,:,:),allocatable :: Tsg ! [s] sediment response time for each sediment class real*8,dimension(:,:,:),allocatable :: Sug ! [m2/s] sediment transport x-dir. for each sediment class (excluding pores) real*8,dimension(:,:,:),allocatable :: Svg ! [m2/s] sediment transport y-dir. for each sediment class (excluding pores) real*8,dimension(:,:,:),allocatable :: ceqg ! [m3/m3] depth-averaged equilibrium concentration for each sediment class real*8,dimension(:,:),allocatable :: ua ! [m/s] time averaged flow velocity due to wave assymetry real*8,dimension(:,:),allocatable :: BR ! [-] maximum wave surface slope used in roller dissipation formulation real*8,dimension(:,:),allocatable :: kb ! [m^2/s^2] near bed turbulence intensity due to depth induces breaking real*8,dimension(:,:),allocatable :: Tbore ! [s] wave period interval associated with breaking induced turbulence real*8,dimension(:,:),allocatable :: uon ! [m/s] onshore directed peak orbital velocity real*8,dimension(:,:),allocatable :: uoff ! [m/s] offshore directed peak orbital velocity real*8,dimension(:,:),allocatable :: dzav ! [m] total bed level change due to avalanching real*8,dimension(:,:),allocatable :: maxzs ! [m] maximum elevation in simulation real*8,dimension(:,:),allocatable :: minzs ! [m] minimum elevation in simulation end type contains subroutine grid_bathy(s,par) use params IMPLICIT NONE type(spacepars) :: s type(parameters) :: par character*80 :: fnameh,fnamex,fnamey integer :: i,iy integer :: j integer, external :: readkey_int integer :: itheta real*8, external :: readkey_dbl real*8 :: degrad,E,thetamin,thetamax ! Input file Keyword Default Minimum Maximum s%nx = readkey_int ('params.txt','nx', 50, 2, 10000) s%ny = readkey_int ('params.txt','ny', 2, 2, 10000) s%dx = readkey_dbl ('params.txt','dx', 0.d0, -1d9, 1d9) s%dy = readkey_dbl ('params.txt','dy', 0.d0, -1d9, 1d9) s%xori = readkey_dbl ('params.txt','xori', 0.d0, -1d9, 1d9) s%yori = readkey_dbl ('params.txt','yori', 0.d0, -1d9, 1d9) s%alfa = readkey_dbl ('params.txt','alfa', 0.d0, -360.d0, 360.d0) s%alfa=s%alfa*atan(1.)/45. s%posdwn= readkey_dbl ('params.txt','posdwn', 1.d0, -1.d0, 1.d0) s%vardx= readkey_int ('params.txt','vardx', 0, 0, 1) !Jaap allocate(s%x(1:s%nx+1,1:s%ny+1)) allocate(s%y(1:s%nx+1,1:s%ny+1)) allocate(s%xz(1:s%nx+1)) allocate(s%yz(1:s%ny+1)) allocate(s%xu(1:s%nx+1)) allocate(s%yv(1:s%ny+1)) allocate(s%xw(1:s%nx+1,1:s%ny+1)) allocate(s%yw(1:s%nx+1,1:s%ny+1)) allocate(s%zb(1:s%nx+1,1:s%ny+1)) allocate(s%zb0(1:s%nx+1,1:s%ny+1)) ! ! Create grid ! if(s%vardx==0)then call readkey('params.txt','depfile',fnameh) open(31,file=fnameh) do j=1,s%ny+1 read(31,*)(s%zb(i,j),i=1,s%nx+1) end do close(31) s%zb=-s%zb*s%posdwn do j=1,s%ny+1 do i=1,s%nx+1 s%x(i,j)=(i-1)*s%dx s%y(i,j)=(j-1)*s%dy end do end do elseif(s%vardx==1)then call readkey('params.txt','depfile',fnameh) call readkey('params.txt','xfile',fnamex) call readkey('params.txt','yfile',fnamey) open(31,file=fnameh) open(32,file=fnamex) open(33,file=fnamey) !do j = 1,s%ny+1 ! read(31,*)(s%zb(i,j),i=1,s%nx+1) ! read(32,*)(s%x(i,j),i=1,s%nx+1) ! read(33,*)(s%y(i,j),i=1,s%nx+1) !enddo read(31,*)((s%zb(i,j),i=1,s%nx+1),j=1,s%ny+1) read(32,*)((s%x(i,j),i=1,s%nx+1),j=1,s%ny+1) read(33,*)((s%y(i,j),i=1,s%nx+1),j=1,s%ny+1) close(31) close(32) close(33) if (abs(s%x(1,2)-s%x(1,1))>par%eps) then ! Apparently input grid is at an angle, therefore defined in world coordinates ! Find out grid orientation s%alfa=atan2(s%y(2,1)-s%y(1,1),s%x(2,1)-s%x(1,1)) s%xori=s%x(1,1) s%yori=s%y(1,1) s%xw=s%x s%yw=s%y s%x= cos(s%alfa)*(s%xw-s%xori)+sin(s%alfa)*(s%yw-s%yori) s%y=-sin(s%alfa)*(s%xw-s%xori)+cos(s%alfa)*(s%yw-s%yori) endif s%zb=-s%zb*s%posdwn endif s%xz = s%x(:,1); s%yz = s%y(1,:); s%xu(1:s%nx) = 0.5*(s%xz(1:s%nx)+s%xz(2:s%nx+1)) s%xu(s%nx+1) = s%xz(s%nx+1)+0.5*(s%xz(s%nx+1)-s%xz(s%nx)) s%yv(1:s%ny) = 0.5*(s%yz(1:s%ny)+s%yz(2:s%ny+1)) s%yv(s%ny+1) = s%yz(s%ny+1)+0.5*(s%yz(s%ny+1)-s%yz(s%ny)) s%xw=s%xori+s%x*cos(s%alfa)-s%y*sin(s%alfa) s%yw=s%yori+s%x*sin(s%alfa)+s%y*cos(s%alfa) s%zb0 = s%zb; ! ! Specify theta-grid ! ! ! from Nautical wave directions in degrees to Cartesian wave directions in radian !!! ! s%theta0=(1.5*par%px-s%alfa)-par%dir0*atan(1.)/45 ! Updated in waveparams.f90 for instat 4,5,6,7 degrad=par%px/180.; if (par%thetanaut==1) then thetamin=(270-par%thetamax)*degrad-s%alfa thetamax=(270-par%thetamin)*degrad-s%alfa if (thetamax>par%px) then thetamax=thetamax-2.d0*par%px thetamin=thetamin-2.d0*par%px endif if (thetamin<-par%px) then thetamax=thetamax+2.d0*par%px thetamin=thetamin+2.d0*par%px endif else thetamin=par%thetamin*degrad; thetamax=par%thetamax*degrad; endif s%dtheta=par%dtheta*degrad; s%ntheta=(thetamax-thetamin)/s%dtheta allocate(s%theta(1:s%ntheta)) allocate(s%cxsth(1:s%ntheta)) allocate(s%sxnth(1:s%ntheta)) do itheta=1,s%ntheta s%theta(itheta)=thetamin+s%dtheta/2+s%dtheta*(itheta-1) end do s%cxsth=cos(s%theta); s%sxnth=sin(s%theta); end subroutine end module