module wave_impact_module contains !==========================================================================! ! ! Subroutine to determine wave erosion due to the impact of waves on a steep slope ! based on Mariotti and Fagherazzi (2010) or Roelvink and Reniers (2012) or new approach ! subroutine wave_impact(s,par,sed_input,local_wetz) ! use params use spaceparams use reflection_module ! implicit none type(parameters) :: par type(spacepars) :: s real*8, dimension(par%nx+1,par%ny+1,par%ngd), intent(out) :: sed_input ! [m3/m2/s] flux of material to fluid phase coming from dry cells real*8, dimension(par%ngd) :: temp_ero_dry ! [m3/m2/s] total eroded material from dry cells per number of fractions real*8, dimension(par%nx+1,par%ny+1) :: hloc real*8 :: Pcr ! [W/m] critical wave power threshold for erosion real*8 :: Pw ! [W/m] wave enery flux due to rollers and waves real*8 :: betaMF10 ! [kg/J] "erodibility" parameter from Mariotti and Fagherazzi (2010) real*8 :: betaMal11 ! [m3/J] "erodibility" parameter from Marani et al. (2011) !real*8 :: calwimp ! [-] calibration factor for the erodibility parameter real*8 :: beta ! [m3/J] or [kg/J] (depending on vol_mass) "erodibility" parameter real*8 :: ee ! [m3/m/s] erosion rate expressed in meter^2 of sediment (without pores) per second (m^2 are referred as a dx*dz quantity) integer, dimension(par%nx+1,par%ny+1), intent(in) :: local_wetz ! [-] mask to keep track of udpated wet/dry cells integer, dimension(par%nx+1) :: cellscarp ! [-] number of cells composing a scarp integer, dimension(par%nx+1) :: toe ! [-] vector containing the TOE position for each scarp identified integer, dimension(par%nx+1) :: top ! [-] vector containing the TOP position for each scarp identified integer :: counter ! [-] counter of the number of scarps integer :: subm ! [-] flag to check whether the domain is completely submerged (1) or not (0) integer :: vol_mass ! [-] flag to decide if Mal11 (1) or MF10 (0) approach is used integer :: i,j,jg,kk ! ! hloc = max(s%hh,par%eps) Pcr = 0.0d0 betaMF10 = 1.4E-04 betaMal11 = 1.2E-09 beta = 0.0d0 !calwimp = 20 temp_ero_dry = 0.d0 top = 0 toe = 0 cellscarp = 0 counter = 0 vol_mass = 1 s%ero_scarp = 0.0d0 ! ! ! Check if the domain is completely submerged (or completely dry) ! call submergence(s,par,local_wetz,subm) ! ! If the domain is not completely submerged or dry ! if (subm==0) then ! ! Transformation from [m3/m2/s] to [kg/m2/s] if necessary ! if (vol_mass==0) then do jg=1,par%ngd s%ero(:,:,jg) = s%ero(:,:,jg)*par%rhos end do beta = betaMF10 else beta = betaMal11 end if do j=1,par%ny+1 do i=1,par%nx-1 ! ! Identify the cell at the waterline (i is the last wet cell at waterline) ! TODO: ACCOUNT FOR THE CASE THERE IS A TIDAL POOL "SEPARATING" FLUID PHASE ! if (local_wetz(i,j)==1 .and. local_wetz(i+1,j)==0) then ! ! Calculate wave power at the waterline ! call wave_energy_flux(s,par,i,j,Pw) ! ! Determine erosion rate [m3/m/s] due to wave impact ! ee = max(0.d0,Pw-Pcr)*beta*(1-par%por)*par%calwimp ! ! If no erosion occurs, skip the following commands ! if (ee>0.d0) then ! ! Determine the retreat of the scarp ! if (par%retreat==1) then ! ! Retreat approach in which two cells at waterline are eroded ! proportionally to the water depth at the last wet cell ! call retreat_scheme_1(s,par,hloc,ee,temp_ero_dry,i,j) ! ! else ! ! Retreat approach in which erosion is proportional to ! the domain slope in order to represent the horizontal ! translation of the bank scarp. ! ! Identify the cells being part of the scarp ! call identify_scarp(s,par,cellscarp,toe,top,counter) ! ! if (sum(cellscarp)<=1) then ! ! In case no cell has slope>=scarpsl no scarp erosion occurs ! go to 10 ! else ! do kk=1,counter if (i>=toe(kk) .and. i<=top(kk)) then ! ! Retreat scheme is only applied to the scarp intersecting the waterline ! if (par%retreat==2) then ! ! Linear retreat scheme ! call retreat_scheme_2(s,par,ee,temp_ero_dry,toe(kk),top(kk),cellscarp(kk),i,j) ! else if (par%retreat==3) then ! ! Nonlinear retreat scheme ! call retreat_scheme_3(s,par,ee,temp_ero_dry,toe(kk),top(kk),cellscarp(kk),i,j) ! end if ! end if end do end if ! if statement for cellscarp = 0 end if ! if statement for retreat 1 or 2,3 ! ! Calculate flux of material eroded from dry cells as input for adv. diff. eq.; ! material is not added at a single cell but distributed over cells ! call distr_ero_dry(s,par,temp_ero_dry,sed_input,i,j) ! ! end if ! if statement for occurrence of erosion end if ! if to identify waterline end do ! do in x direction end do ! do in y direction end if ! if statement for submerged/dry domain ! ! Transformation from [kg/m2/s] to [m3/m2/s] in case MF10 approach is used ! 10 if (vol_mass==0) then s%ero = s%ero/par%rhos end if ! end subroutine wave_impact !==========================================================================! subroutine retreat_scheme_1(s,par,hloc,ee,temp_ero_dry,i,j) ! ! Subroutine to describe the lateral retreat assuming erosion interests the last wet cell ! and the first dry cell. Eosion rate is split among these two cells based on water depth ! at last wet cell ! use params use spaceparams ! implicit none type(parameters) :: par type(spacepars) :: s real*8, dimension(par%nx+1,par%ny+1), intent(in) :: hloc real*8, intent(in) :: ee ! [m3/m/s] erosion rate expressed in meter^2 of sediment (without pores) per second (m^2 are referred as a dx*dz quantity) integer, intent(in) :: i ! [-] position of the last wet cell integer, intent(in) :: j real*8, dimension(par%ngd), intent(out) :: temp_ero_dry ! [m3/m2/s] total eroded material from dry cells per number of fractions real*8, dimension(2) :: weight ! [-] weight function for erosion based on water depth real*8, dimension(par%ngd) :: cell_ero_dry ! [m3/m2/s] erosion rate at dry cell real*8, dimension(par%ngd) :: cell_ero_wet ! [m3/m2/s] erosion rate at wet cell real*8 :: global_mass ! [m3/m2/s] variable to check for mass (rate) conservation real*8 :: mass_diff ! [m3/m2/s] variable to check for mass (rate) conservation real*8 :: wave_L ! [m] wave length real*8 :: wave_L_temp ! [m] "cumulative" wave length real*8 :: x_temp real*8 :: local_slope real*8, dimension(par%nx+1,par%ngd) :: distributed_ero real*8, dimension(par%ngd) :: mass_diff2 integer :: jg,kk integer :: ncells integer :: counter ! ! ! Initialize variables ! global_mass = 0.0d0 mass_diff = 0.0d0 cell_ero_dry = 0.0d0 cell_ero_wet = 0.0d0 local_slope = 0.0d0 distributed_ero = 0.0d0 mass_diff2 = 0.d0 ! ! Split erosion rate among last wet cell and first dry cell ! weight(2) = min( hloc(i,j)/max(s%zb(i+1,j)-s%zb(i,j),par%eps), 1.0d0 ) weight(1) = abs(1.0d0-weight(2)) ! if ((weight(1)+weight(2))>1.1 .or. (weight(1)+weight(2))<0.9) then write(*,*) '!!...unconsistent impact erosion weights...' weight(1) = 0.5 weight(2) = 0.5 end if ! ! Update erosion for the last wet cell ! cell_ero_wet = ee/s%dsz(i,j)*s%pbbed(i,j,1,:)*weight(1) s%ero_scarp(i,j,:) = cell_ero_wet s%ero(i,j,:) = s%ero(i,j,:)+ee/s%dsz(i,j)*s%pbbed(i,j,1,:)*weight(1) ! ! Update erosion at first dry cell (and store eroded material for adv. diff. eq. input) ! cell_ero_dry = ee/s%dsz(i+1,j)*s%pbbed(i+1,j,1,:)*weight(2) s%ero_scarp(i+1,j,:) = cell_ero_dry s%ero(i+1,j,:) = s%ero(i+1,j,:)+cell_ero_dry temp_ero_dry = cell_ero_dry ! ! Check for mass conservation ! global_mass = sum(cell_ero_wet*s%dsz(i,j)+cell_ero_dry*s%dsz(i,j)) mass_diff = ee-global_mass ! ! In case cell at waterline is lower than previous cell, than material eroded at wet cell ! is split linearly decreasing over ! if (s%zb(i-1,j)>s%zb(i,j)) then local_slope = (s%zb(i,j)-s%zb(i-1,j))/s%dsu(i-1,j) counter = i-1 do while (local_slope<0.d0) local_slope = (s%zb(counter,j)-s%zb(counter-1,j))/s%dsu(counter-1,j) counter = counter-1 if (counter==0) then exit end if end do ncells = i-counter+1 s%ero(i,j,:) = s%ero(i,j,:)-cell_ero_wet do kk=1,ncells s%ero(i+1-kk,j,:) = s%ero(i+1-kk,j,:) + cell_ero_wet*s%dsz(i+1-kk,j)/sum(s%dsz(i-ncells+1:i,j)) distributed_ero(kk,:) = cell_ero_wet*s%dsz(i+1-kk,j)/sum(s%dsz(i-ncells+1:i,j)) end do !wave_L = 2.d0*par%px/s%k(i,j) !wave_L_temp = 0.d0 !counter = i !do while (wave_L>=wave_L_temp) ! wave_L_temp = wave_L_temp+s%dsz(counter,j) ! counter = counter-1 ! if (counter==0) then ! exit ! end if !end do !ncells = i-counter !s%ero(i,j,:) = s%ero(i,j,:)-cell_ero_wet !x_temp = s%dsz(i,j)/2 !do kk=1,ncells ! distributed_ero(kk,:) = 2.0d0*cell_ero_wet/wave_L_temp*(1-x_temp/wave_L_temp)*s%dsz(i+1-kk,j) ! s%ero(i+1-kk,j,:) = s%ero(i+1-kk,j,:)+2.0d0*cell_ero_wet/wave_L_temp*(1-x_temp/wave_L_temp)*s%dsz(i+1-kk,j) ! ! ! ! Increase x_temp in order to be always at cell center ! ! ! x_temp = x_temp+0.5d0*(s%dsz(i+1-kk,j)+s%dsz(i-kk,j)) ! ! !end do end if ! ! Check for mass erosion rate conservation after redistribution of eroded material ! do jg=1,par%ngd mass_diff2(jg) = sum(distributed_ero(:,jg))-cell_ero_wet(jg) end do ! ! end subroutine retreat_scheme_1 !==========================================================================! subroutine retreat_scheme_2(s,par,ee,temp_ero_dry,S_toe,S_top,S_cellscarp,i,j) ! ! Subroutine to describe the lateral retreat considering the bank is moving with a ! constant celerity proportional to the wave impact erosion rate ! use params use spaceparams ! implicit none type(parameters) :: par type(spacepars) :: s real*8, intent(in) :: ee ! [m3/m/s] erosion rate expressed in meter^2 of sediment (without pores) per second (m^2 are referred as a dx*dz quantity) integer, intent(in) :: S_toe ! [-] cell number referred to the TOE of the scarp intersecting the waterline integer, intent(in) :: S_top ! [-] cell number referred to the TOP of the scarp intersecting the waterline integer, intent(in) :: i ! [-] cell number at waterline integer, intent(in) :: j integer, intent(in) :: S_cellscarp ! [-] number of cells beaing part of the scarp real*8, dimension(par%ngd), intent(out) :: temp_ero_dry ! [m3/m2/s] total eroded material from dry cells per fractions real*8 :: scarp_celerity ! [m3/m2/s] average erosion rate over the scarp height real*8 :: cell_ero ! [m3/m2/s] erosion rate per cell real*8 :: mindt real*8 :: global_mass real*8 :: mass_diff !real*8 :: mydeltat ! [-] deltat in case we want a CFL=par%CFL integer :: jg,ii ! ! ! Determine scarp_celerity which represents the ! celerity of the scarp migrating horizontally ! scarp_celerity = 0.0d0 do ii=1,S_cellscarp scarp_celerity = scarp_celerity+(s%zb(S_toe+ii-1,j)-s%zb(S_toe+ii-2,j))/s%dsu(S_toe+ii-2,j)*s%dsz(S_toe+ii-1,j) end do scarp_celerity = ee/scarp_celerity ! ! Initialize some variables ! temp_ero_dry = 0.0d0 global_mass = 0.0d0 mass_diff = 0.0d0 mindt = 0.0d0 ! ! TODO: FIX PROBLEM WITH MASS CONSERVATION, WE HAVE DISCONTINUITIES AT THE BOUNDARIES OF THE SCARP ! WE NEED TO USE GODUNOV'S SCHEME??.. WRITE IT DOWN ! ! ! Control the value of "local" dt (it is to be formally sure; actually ! the opposite problem arises: par%dt is imposed by hydrodynamics, resulting ! in too low values of CFL for the bed update) ! mindt = par%CFL/(par%morfac*scarp_celerity)*minval(s%dsz(S_toe-1:S_top+1,j)) if (par%dt>mindt) then write(*,*) '!!...violated CFL in lateral retreat 2' end if ! do jg=1,par%ngd do ii=1,S_cellscarp ! ! Update erosion rate ! cell_ero = scarp_celerity*s%pbbed(S_toe+ii-1,j,1,jg)*(s%zb(S_toe+ii-1,j)-s%zb(S_toe+ii-2,j))/s%dsu(S_toe+ii-2,j) s%ero_scarp(S_toe+ii-1,j,jg) = cell_ero s%ero(S_toe+ii-1,j,jg) = s%ero(S_toe+ii-1,j,jg)+cell_ero ! ! Account for the material eroded from the dry cells ! if (S_toe+ii-1>i .and. S_toe+ii-1<=S_top) then temp_ero_dry(jg) = temp_ero_dry(jg)+cell_ero end if ! ! Check where the mass goes ! global_mass = global_mass+cell_ero*s%dsz(S_toe+ii-1,j) ! ! Hypothetical dt to have CFL=par%CFL ! !mydeltat = par%CFL*s%dsu(i,j)/scarp_celerity ! end do end do ! ! Check if mass is conserved during the retreat of the scarp, ! if mass_diff = 0 (in percentage) then erosion rate is correctly spread over the scarp ! mass_diff = 2.0d0*(global_mass-ee)/(global_mass+ee) if (abs(mass_diff)>=1.0E-6) then write(*,*) 'someone is stealing mass in retreat 2...' end if ! ! end subroutine retreat_scheme_2 !==========================================================================! subroutine retreat_scheme_3(s,par,ee,temp_ero_dry,S_toe,S_top,S_cellscarp,i,j) ! ! Subroutine to describe the lateral retreat considering the bank is moving with a ! nonuniform celerity that is function of the water depth in front of the scarp. ! Bed update update is described by a nonlinear advection equation and it is ! solved by Godunov's scheme (for details see Toro, 2009) ! use params use spaceparams ! implicit none type(parameters) :: par type(spacepars) :: s real*8, intent(in) :: ee ! [m3/m/s] erosion rate expressed in meter^2 of sediment (without pores) per second (m^2 are referred as a dx*dz quantity) ! integer, intent(in) :: S_toe ! [-] cell number referred to the TOE of the scarp at the waterline integer, intent(in) :: S_top ! [-] cell number referred to the TOP of the scarp at the waterline integer, intent(in) :: i integer, intent(in) :: j integer, intent(in) :: S_cellscarp ! [-] number of cells baing part of the scarp real*8, dimension(par%ngd), intent(out) :: temp_ero_dry ! [m3/m2/s] total eroded material from dry cells per fractions real*8 :: csc ! [m/s] constant used to describe non-uniform scarp celerity real*8 :: cell_ero ! [m3/m2/s] erosion rate per cell real*8 :: global_mass ! [m3/m2/s] variable to check if mass (per unit time) is conserved real*8 :: mass_diff ! (same as above) real*8, dimension(par%nx+1) :: z ! [m] local variable to apply the numerical scheme real*8, dimension(par%nx+1) :: F,FL,FR ! [m2/s] vectors of fluxes real*8, dimension(par%nx+1) :: lambdaL,lambdaR ! [m/s] vectors of characteristic curves speed real*8, dimension(par%nx+1) :: sRH ! [m/s] speed of the discontinuity (when the bank reaches a "vertical" slope) real*8 :: h,wl ! [m] local variables to apply the numerical scheme real*8 :: zp,ze ! [m] bank top and toe elevation in z coordinate framework real*8, dimension(2) :: maxlambda ! [m/s] maximum characteristic speed real*8 :: mindt real*8 :: checkFlux ! [m/s] variable to check if FL(i)=FR(i-1) integer :: jg,ii ! ! ! Initialize variables ! csc = 0.0d0 F = 0.0d0 FL = 0.0d0 FR = 0.0d0 lambdaL = 0.0d0 lambdaR = 0.0d0 maxlambda = 0.0d0 z = 0.0d0 h = 0.0d0 sRH = 0.0d0 temp_ero_dry = 0.0d0 global_mass = 0.0d0 mass_diff = 0.0d0 ! ! Determine a local coordinate system (perhaps it is not strictly necessary) ! !if (par%t>=1908.01360850388) then ! write(*,*), 'findwhy' !end if z = s%zb(:,j)-minval(s%zb(:,j)) h = (s%zs(i,j)-s%zb(S_toe-1,j)) - min(par%eps,abs(s%zs(i,j)-s%zb(i+1,j))) ze = z(S_toe-1) zp = z(S_top) wl = h+ze if (wl>=zp) then ! ! Since wl is mesured using s%zs at the bank toe cell, it is possible ! that the waterline elevation is lower than it. It means the cell at the top is ! dry but the wl elevation at the bank top is higher. In such a case the following ! value for csc is used. ! csc = ee/(zp-ze) else csc = ee*2.0d0 / (zp+h-ze) end if z(S_toe-2) = z(S_toe-1) z(S_top+1) = z(S_top) do ii=-1,S_cellscarp-1 ! ! Loop on cell interfaces ! checkFlux = 0.0d0 ! TODO: check what happens in case we are at the first cell of the scarp if (z(S_toe+ii+1)>=wl .and. z(S_toe+ii)>=wl) then ! ! Both cells are above the waterline ! ! Determine characteristics speed ! lambdaL(S_toe+ii) = lambda_emerged(csc,zp,ze,h,z(S_toe+ii)) lambdaR(S_toe+ii) = lambda_emerged(csc,zp,ze,h,z(S_toe+ii+1)) ! ! Fluxes ! FL(S_toe+ii) = flux_emerged(csc,zp,ze,h,z(S_toe+ii)) FR(S_toe+ii) = flux_emerged(csc,zp,ze,h,z(S_toe+ii+1)) ! ! Rankine-Hugoniot condition ! sRH(S_toe+ii) = rankine_hugoniot(FL(S_toe+ii),FR(S_toe+ii),z(S_toe+ii),z(S_toe+ii+1)) ! else if (z(S_toe+ii+1)>=wl .and. z(S_toe+ii)<=wl) then ! ! Mean water level is at the interface between two cells ! ! Determine characteristics speed ! lambdaL(S_toe+ii) = lambda_submerged(csc,ze,z(S_toe+ii)) lambdaR(S_toe+ii) = lambda_emerged(csc,zp,ze,h,z(S_toe+ii+1)) ! ! Fluxes ! FL(S_toe+ii) = flux_submerged(csc,ze,z(S_toe+ii)) FR(S_toe+ii) = flux_emerged(csc,zp,ze,h,z(S_toe+ii+1)) ! ! Rankine-Hugoniot condition ! sRH(S_toe+ii) = rankine_hugoniot(FL(S_toe+ii),FR(S_toe+ii),z(S_toe+ii),z(S_toe+ii+1)) ! else ! ! Both cells are below mean water level ! ! Determine characteristics speed ! lambdaL(S_toe+ii) = lambda_submerged(csc,ze,z(S_toe+ii)) !perhaps the call of these functions is redundant... lambdaR(S_toe+ii) = lambda_submerged(csc,ze,z(S_toe+ii)) ! ! Fluxes ! FL(S_toe+ii) = flux_submerged(csc,ze,z(S_toe+ii)) FR(S_toe+ii) = flux_submerged(csc,ze,z(S_toe+ii+1)) ! ! Rankine-Hugoniot condition ! sRH(S_toe+ii) = rankine_hugoniot(FL(S_toe+ii),FR(S_toe+ii),z(S_toe+ii),z(S_toe+ii+1)) ! end if ! ! Check for the accuracy of fluxes ! checkFlux = FL(S_toe+ii)-FR(S_toe+ii-1) if (abs(checkFlux)>1.0E-09) then write(*,*) 'wrong fluxes in retreat 3..?', checkFlux end if ! ! Behavior of characteristics curves and determination of the vector containing ! fluxes used to update the bed ! F(S_toe+ii) = getFlux(FL(S_toe+ii),FR(S_toe+ii),lambdaL(S_toe+ii),lambdaR(S_toe+ii),sRH(S_toe+ii)) ! end do ! ! Check the dt (it is to be formally sure; actually the opposite problem arises: par%dt ! is imposed by hydrodynamics, resulting a in too low value of CFL for the bed update) ! TODO: to reduce numerical diffusion update the bed not at every time step...?? ! maxlambda = max(abs(lambdaR),abs(lambdaL)) ! if it used max (min) instead of maxval (min) compiler error (???) if (maxval(maxlambda)>0.0d0) then mindt = minval(par%CFL*s%dsu(:,j)/(par%morfac*maxval(maxlambda))) if (par%dt>mindt) then write(*,*) '!!...violated CFL in lateral retreat 3' end if end if ! ! Boundary conditions ! F(S_toe-2) = F(S_toe-1) F(S_top+1) = F(S_top) ! ! Starto to update the bed (actually s%ero is updated, not s%zb) ! do jg=1,par%ngd do ii=0,S_cellscarp-1 ! ! Erosion rate at single cell ( it represents [z(n+1)-z(n)]/dt that is a negative quantity in ! case of erosion, but in XBeach [z(n+1)-z(n)]/dt = -ero ) ! cell_ero = -(F(S_toe+ii)-F(S_toe+ii-1))/s%dsz(S_toe+ii,j) * s%pbbed(S_toe+ii,j,1,jg) cell_ero = -cell_ero s%ero_scarp(S_toe+ii,j,jg) = cell_ero s%ero(S_toe+ii,j,jg) = s%ero(S_toe+ii,j,jg)+cell_ero ! ! Account for the material eroded from the dry cells ! if (S_toe+ii>i .and. S_toe+ii<=S_top) then temp_ero_dry(jg) = temp_ero_dry(jg)+cell_ero end if ! ! Check where the mass goes ! global_mass = global_mass+cell_ero*s%dsz(S_toe+ii,j) ! end do end do ! ! Check if mass is conserved during the retreat of the scarp. ! It would be sufficient to check if the output flux is equal to ee ! (and it is the case), but some errors can arise from teh subsequent ! procedures. mass_diff is expressed in percentage ! mass_diff = 2.0d0*(global_mass-ee)/(global_mass+ee) if (abs(mass_diff)>=1.0E-6) then write(*,*) 'someone is stealing mass in retreat 3...' end if ! ! end subroutine retreat_scheme_3 !==========================================================================! function lambda_submerged(csc,ze,zi) ! ! Characteristic speed submerged cell ! implicit none real*8, intent(in) :: csc real*8, intent(in) :: ze real*8, intent(in) :: zi real*8 :: lambda_submerged ! ! if (ze<=zi) then lambda_submerged = csc else lambda_submerged = 0.0d0 end if ! end function lambda_submerged !==========================================================================! function lambda_emerged(csc,zp,ze,h,zi) ! ! Characteristic speed emerged cell ! implicit none real*8, intent(in) :: csc real*8, intent(in) :: zp real*8, intent(in) :: ze real*8, intent(in) :: h real*8, intent(in) :: zi real*8 :: lambda_emerged ! ! if (zp>=zi) then if (abs(zp-ze-h)<1.0E-12) then lambda_emerged = 0.d0 else lambda_emerged = csc/(zp-ze-h)*(zp-zi) end if else lambda_emerged = 0.0d0 end if ! end function lambda_emerged !==========================================================================! function flux_submerged(csc,ze,zi) ! ! Flux submerged cell ! implicit none real*8, intent(in) :: csc real*8, intent(in) :: ze real*8, intent(in) :: zi real*8 :: flux_submerged ! ! if (zi>=ze) then flux_submerged = csc*(zi-ze) else flux_submerged = 0.0d0 end if ! end function flux_submerged !==========================================================================! function flux_emerged(csc,zp,ze,h,zi) ! ! Flux emerged cell ! implicit none real*8, intent(in) :: csc real*8, intent(in) :: zp real*8, intent(in) :: ze real*8, intent(in) :: h real*8, intent(in) :: zi real*8 :: Htot real*8 :: flux_emerged ! ! Htot = ze+h if (abs(zp-Htot)<1.0E-12) then flux_emerged = csc*h else if (zi<=zp) then flux_emerged = csc/(zp-Htot)*zi*(zp-zi/2.0d0) + csc/(zp-Htot)*Htot*(Htot/2.0d0-zp) + csc*h else flux_emerged = 0.5d0*csc/(zp-Htot)*zp**2.0d0 + csc/(zp-Htot)*Htot*(Htot/2.0d0-zp) + csc*h end if end if ! end function flux_emerged !==========================================================================! function rankine_hugoniot(fL,fR,zi,zii) ! ! Determine the speed of the discontinuity ! implicit none real*8, intent(in) :: fL real*8, intent(in) :: fR real*8, intent(in) :: zi real*8, intent(in) :: zii real*8 :: rankine_hugoniot ! ! if (zii==zi) then rankine_hugoniot = 0.0d0 else rankine_hugoniot = (fR-fL)/(zii-zi) end if ! end function rankine_hugoniot !==========================================================================! function getFlux(FL,FR,lambdaL,lambdaR,sRH) ! ! Check the characteristic curves behavior and get ! the right value for the flux at cell interface (i+1/2) ! implicit none real*8, intent(in) :: FL,FR real*8, intent(in) :: lambdaL,lambdaR real*8, intent(in) :: sRH real*8 :: getFlux ! ! if (lambdaL>lambdaR) then if (sRH>=0.0d0) then getFlux = FL else getFlux = FR end if else if (lambdaL>=0.0d0) then getFlux = FL elseif (lambdaR<=0.0d0) then getFlux = FR else getFlux = 0.0d0 end if end if ! end function getFlux !==========================================================================! subroutine identify_scarp(s,par,cellscarp,toe,top,counter) ! use params use spaceparams ! implicit none type(parameters) :: par type(spacepars) :: s real*8, dimension(par%nx+1,par%ny+1) :: dzbdxup !integer, dimension(par%nx+1,par%ny+1), intent(out) :: sl ! [-] flag to identify cells being part of the scarp integer, dimension(par%nx+1), intent(out) :: cellscarp ! [-] number of cells composing the scarp integer, dimension(par%nx+1), intent(out) :: toe ! [-] First cell being part of the scarp integer, dimension(par%nx+1), intent(out) :: top ! [-] Last cell being part of the scarp integer, intent(out) :: counter ! [-] identifier of the number of scarps integer :: ii,jj ! ! ! TODO: CONSIDER TO SET A THRESHOLD FOR SCARP HEIGHT ! Determine local slope of the domain ! dzbdxup = 0.0d0 top = 0 toe = 0 cellscarp = 0 counter = 0 ! do jj=1,par%ny+1 do ii=2,par%nx+1 ! ! The slope of the scarp is determined in an upwind manner ! dzbdxup(ii,jj)=(s%zb(ii,jj)-s%zb(ii-1,jj))/s%dsu(ii-1,jj) ! end do do ii=2,par%nx ! ! Identify cells being part of the scarp (slope >= scarpsl) ! if (dzbdxup(ii,jj)>=par%scarpsl .and. dzbdxup(ii-1,jj)=par%scarpsl .and. dzbdxup(ii+1,jj)=par%scarpsl .and. dzbdxup(ii-1,jj)>=par%scarpsl) then ! ! This is a cell within the scarp ! !sl(ii,jj) = 1 cellscarp(counter) = cellscarp(counter)+1 ! if (dzbdxup(ii,jj)>=par%scarpsl .and. dzbdxup(ii+1,jj)=wave_L_temp) wave_L_temp = wave_L_temp+s%dsz(counter,j) counter = counter-1 if (counter==0) then exit end if end do ncells = i-counter x_temp = s%dsz(i,j)/2 do kk=1,ncells sed_input(i+1-kk,j,jg) = 2.0d0*temp_ero_dry(jg)/wave_L_temp*(1-x_temp/wave_L_temp)*s%dsz(i+1-kk,j) ! ! Increase x_temp in order to be always at cell center ! x_temp = x_temp+0.5d0*(s%dsz(i+1-kk,j)+s%dsz(i-kk,j)) ! end do end do ! ! Check for mass erosion rate conservation ! do jg=1,par%ngd mass_diff(jg) = temp_ero_dry(jg)-sum(sed_input(:,:,jg)) end do if (sum(abs(mass_diff))>=1.0E-09) then write(*,*) 'error in redistributing eroded mass from dry cells...' end if ! ! end subroutine distr_ero_dry !==========================================================================! ! ! Subroutine to determine the wave energy flux approaching the bank and responsible for ! scarp erosion. Wave energy flux is averaged over a wave length from the waterline ! subroutine wave_energy_flux(s,par,i,j,Pw) ! use params use spaceparams ! implicit none type(parameters) :: par type(spacepars) :: s integer, intent(in) :: i integer, intent(in) :: j real*8, intent(out) :: Pw ! [W/m] wave energy flux comprises wave and roller contribution real*8 :: wave_L ! [m] wave length real*8 :: wave_L_temp ! [m] "cumulative" wave length real*8 :: cg_interface ! [m/s] group velocity at cell interface real*8 :: c_interface ! [m/s] phase (and roller) velocity at cell interface real*8 :: W_flux real*8 :: R_flux integer :: ncells ! [-] number of cells covering a wave length integer :: counter integer :: jg,kk ! ! ! Initialize wave energy flux ! Pw = 0.0d0 ! ! Determine wave length ! wave_L = 2.d0*par%px/s%k(i,j) wave_L_temp = 0.d0 counter = i do while (wave_L>=wave_L_temp) ! ! Increase wave_L_temp until wave_L is reached ! wave_L_temp = wave_L_temp+s%dsz(counter,j) counter = counter-1 if (counter==0) then exit end if end do ncells = i-counter wave_L_temp = 0.0d0 do kk=1,ncells ! ! Determine celerity at cell interface ! cg_interface = 0.5*(s%cg(i-kk+1,j)+s%cg(i-kk,j)) c_interface = 0.5*(s%c(i-kk+1,j)+s%c(i-kk,j)) ! ! Determine "weighted" energy fluxes at cell interface ! W_flux = s%E(i-kk,j)*cg_interface*s%dsu(i-kk,j) R_flux = s%R(i-kk,j)*c_interface*s%dsu(i-kk,j) ! ! Update overall wave energy flux ! Pw = Pw+W_flux+R_flux ! ! Update the new wave_L_temp ! wave_L_temp = wave_L_temp+s%dsu(i-kk,j) ! end do ! Pw = Pw/wave_L_temp ! ! end subroutine wave_energy_flux !==========================================================================! end module wave_impact_module