module morphevolution_1 implicit none save #ifdef HAVE_CONFIG_H #include "config.h" #endif contains subroutine transus_1(s,par) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Copyright (C) 2007 UNESCO-IHE, WL|Delft Hydraulics and Delft University ! ! Dano Roelvink, Ap van Dongeren, Ad Reniers, Jamie Lescinski, ! ! Jaap van Thiel de Vries, Robert McCall ! ! ! ! d.roelvink@unesco-ihe.org ! ! UNESCO-IHE Institute for Water Education ! ! P.O. Box 3015 ! ! 2601 DA Delft ! ! The Netherlands ! ! ! ! This library is free software; you can redistribute it and/or ! ! modify it under the terms of the GNU Lesser General Public ! ! License as published by the Free Software Foundation; either ! ! version 2.1 of the License, or (at your option) any later version. ! ! ! ! This library is distributed in the hope that it will be useful, ! ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU ! ! Lesser General Public License for more details. ! ! ! ! You should have received a copy of the GNU Lesser General Public ! ! License along with this library; if not, write to the Free Software ! ! Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 ! ! USA ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! use params use spaceparams use xmpi_module use interp use paramsconst use morphevolution ! use vsmumod implicit none type(spacepars),target :: s type(parameters) :: par integer :: i,isig integer :: j,jg real*8 :: exp_ero real*8,dimension(:),allocatable,save :: chain,cumchain real*8,dimension(:,:),allocatable,save :: vmag2,uau,uav,um,vm,ueu_sed,uev_sed,veu_sed,vev_sed real*8,dimension(:,:),allocatable,save :: ccvt,dcdz,dsigt,aref real*8,dimension(:,:),allocatable,save :: cc,ccb,cu,cv,Sus,Svs real*8,dimension(:,:),allocatable,save :: cub,cvb,Sub,Svb,pbbedu,pbbedv real*8,dimension(:,:),allocatable,save :: suq3d,svq3d,eswmax,eswbed,sigs,deltas real*8,dimension(:,:,:),allocatable,save :: dsig,ccv,sdif,cuq3d,cvq3d,fac real*8,dimension(:,:),allocatable,save :: sinthm,costhm real*8 :: delta,delta_x,shields,ftheta,psi_x,Sbtot ! Lodewijk: for direction of sediment transport (bed slope effect) real*8 :: Ssmtot, dzbds, Sbmtot ! Lodewijk: for magnitude of sediment transport (bed slope effect) real*8 :: depthcheck real*8,dimension(:,:),allocatable,save :: trans_inf !Huub: concept model; !include 's.ind' !include 's.inp' if (.not. allocated(vmag2)) then allocate(vmag2 (s%nx+1,s%ny+1)) allocate(uau (s%nx+1,s%ny+1)) allocate(uav (s%nx+1,s%ny+1)) allocate(ueu_sed (s%nx+1,s%ny+1)) allocate(uev_sed (s%nx+1,s%ny+1)) allocate(veu_sed (s%nx+1,s%ny+1)) allocate(vev_sed (s%nx+1,s%ny+1)) allocate(cu (s%nx+1,s%ny+1)) allocate(cv (s%nx+1,s%ny+1)) allocate(cc (s%nx+1,s%ny+1)) allocate(ccb (s%nx+1,s%ny+1)) allocate(fac (s%nx+1,s%ny+1,par%ngd)) allocate(Sus (s%nx+1,s%ny+1)) allocate(Svs (s%nx+1,s%ny+1)) allocate(cub (s%nx+1,s%ny+1)) allocate(cvb (s%nx+1,s%ny+1)) allocate(Sub (s%nx+1,s%ny+1)) allocate(Svb (s%nx+1,s%ny+1)) allocate(pbbedu (s%nx+1,s%ny+1)) allocate(pbbedv (s%nx+1,s%ny+1)) allocate(ccvt (s%nx+1,s%ny+1)) allocate(dcdz (s%nx+1,s%ny+1)) allocate(dsigt (s%nx+1,s%ny+1)) allocate(dsig(s%nx+1,s%ny+1,par%kmax)) allocate(ccv(s%nx+1,s%ny+1,par%kmax)) allocate(sdif(s%nx+1,s%ny+1,par%kmax)) allocate(um (s%nx+1,s%ny+1)) allocate(vm (s%nx+1,s%ny+1)) allocate(deltas(s%nx+1,s%ny+1)) allocate(sigs(s%nx+1,s%ny+1)) allocate(eswmax(s%nx+1,s%ny+1)) allocate(eswbed(s%nx+1,s%ny+1)) allocate(suq3d(s%nx+1,s%ny+1)) allocate(svq3d(s%nx+1,s%ny+1)) allocate(cuq3d(s%nx+1,s%ny+1,par%kmax)) allocate(cvq3d(s%nx+1,s%ny+1,par%kmax)) allocate(aref(s%nx+1,s%ny+1)) allocate(chain(par%kmax)) allocate(cumchain(par%kmax)) allocate (sinthm(s%nx+1,s%ny+1)) allocate (costhm(s%nx+1,s%ny+1)) allocate(trans_inf(s%nx+1,s%ny+1)) !Huub: concept model; delta_x = 0.d0 ! Lodewijk shields = 0.d0 ! Lodewijk ftheta = 0.d0 ! Lodewijk psi_x = 0.d0 ! Lodewijk Sbtot = 0.d0 ! Lodewijk delta = 0.d0 ! Lodewijk uau = 0.d0 uav = 0.d0 um = 0.d0 vm = 0.d0 chain = 0.0d0 cumchain = 0.0d0 fac = 1.d0 exp_ero = 0.d0 ! generate sigma grid shape... do isig=2,par%kmax chain(isig) = par%sigfac**(isig-1) cumchain(isig) = cumchain(isig-1)+chain(isig) enddo endif ! use eulerian velocities vmag2 = s%ue**2+s%ve**2 cu = 0.0d0 cv = 0.0d0 cub = 0.0d0 cvb = 0.0d0 s%dcsdx = 0.0d0 s%dcsdy = 0.0d0 sinthm = sin(s%thetamean-s%alfaz) costhm = cos(s%thetamean-s%alfaz) ! short wave runup if (par%swrunup==1 .and. par%struct==1) then call hybrid(s,par) endif ! compute turbulence due to wave breaking if (par%lwt==1 .or. par%turb == TURB_BORE_AVERAGED .or. par%turb == TURB_WAVE_AVERAGED) then call waveturb(s,par) endif if (par%swave==1) then ! include wave skewness and assymetry in sediment advection velocity if (par%waveform==WAVEFORM_RUESSINK_VANRIJN)then call RvR(s,par) elseif (par%waveform==WAVEFORM_VANTHIEL) then call vT(s,par) endif endif call sedtransform(s,par) !Calculate suspended and bed load transport sand. call mccall_vanrijn_1(s,par) !Calculate bed load transport gravel if (par%subsus==1) then call rijper_vanrijn(s,par) !Calculate bed transport of sand under gravel endif ! ! Do loop to differentiate between gravel and sand formulations (Huub) Does not work for concept composite model. ! do jg = 1,par%ngd ! if (s%D50(jg)>0.002) then ! select case (par%formG) ! case (FORM_NIELSEN2006) ! call Nielsen2006(s,par) ! ! case (FORM_MCCALL_VANRIJN) ! call mccall_vanrijn(s,par) ! end select ! else ! endif ! end do ! ! calculate equilibrium concentration/sediment source ! select case (par%form) ! case (FORM_SOULSBY_VANRIJN,FORM_VANTHIEL_VANRIJN,FORM_VANRIJN1993) ! ! Soulsby van Rijn and Van Thiel de Vries & Reniers 2008 formulations ! call sedtransform(s,par) ! case (FORM_NIELSEN2006) ! call Nielsen2006(s,par) ! if (par%bulk==0) then ! return ! endif ! case (FORM_MCCALL_VANRIJN) ! call mccall_vanrijn(s,par) ! if (par%bulk==0) then ! return ! endif ! end select ! compute reduction factor for sediment sources due to presence of hard layers jg = 1 do j=1,s%ny+1 do i=1,s%nx+1 exp_ero = par%morfac*par%dt/(1.d0-par%por)*s%hh(i,j)*(s%ceqsg(i,j,jg)*s%pbbed(i,j,1,jg)/s%Tsg(i,j,jg) & + s%ceqbg(i,j,jg)*s%pbbed(i,j,1,jg)/par%dt) ! limit erosion to available sediment on top of hard layer wwvv changed tiny into epsilon fac(i,j,jg) = min(1.d0,s%structdepth(i,j)*s%pbbed(i,j,1,jg)/max(epsilon(0.d0),exp_ero) ) !if (fac(i,j,jg)*exp_ero > dzbed(i,j,1)*pbbed(i,j,1,jg)) then ! limit erosion to available sand in top layer ! fac(i,j,jg) = min(fac(i,j,jg),dzbed(i,j,1)*pbbed(i,j,1,jg)/max(tiny(0.d0),exp_ero) ) ! write(*,*)'WARNING: expected erosion from top layer is larger than available sand in top layer' !endif enddo enddo ! compute diffusion coefficient s%Dc = par%facDc*(par%nuh+par%nuhfac*s%hh*(s%DR/par%rho)**(1.d0/3.d0)) jg = 1 cc = s%ccg(:,:,jg) if (s%D50(jg)>0.002d0) then ! RJ: set ceqsg to zero for gravel. ! Dano: try without this fix cc = 0.d0 ! Can be used to test total transport mode endif ! ! X-direction ! ! Get ua in u points and split out in u and v direction uau(1:s%nx,:) = 0.5*(s%ua(1:s%nx,:)*costhm(1:s%nx,:)+s%ua(2:s%nx+1,:)*costhm(1:s%nx,:)) uav(1:s%nx,:) = 0.5*(s%ua(1:s%nx,:)*sinthm(1:s%nx,:)+s%ua(2:s%nx+1,:)*sinthm(1:s%nx,:)) if (par%nz>1) then ueu_sed(1:s%nx,:) = 0.5*(s%ue_sed(1:s%nx,:)+s%ue_sed(2:s%nx+1,:)) else ueu_sed=s%ueu endif veu_sed(1:s%nx,:) = 0.5*(s%ve_sed(1:s%nx,:)+s%ve_sed(2:s%nx+1,:)) if (xmpi_isbot) then veu_sed(s%nx+1,:) = veu_sed(s%nx,:) endif ! Compute vmagu including ua ! s%vmagu = sqrt((s%uu+uau)**2+(s%vu+uav)**2) s%vmagu = sqrt((ueu_sed+uau)**2+(veu_sed+uav)**2) ! sediment advection velocity for suspended load and bed load respectively ! REMARK: when vreps does not equal vv; no mass conservation ! s%ureps = s%ueu+uau ! s%urepb = s%ueu+uau ! RJ maybe reduce this velocity? s%ureps = ueu_sed+uau s%urepb = ueu_sed+uau ! do j=1,s%ny+1 do i=1,s%nx if(s%ureps(i,j)>0.d0) then ! test cu(i,j)=cc(i,j) cu(i,j)=par%thetanum*cc(i,j)+(1.d0-par%thetanum)*cc(min(i+1,s%nx),j) cub(i,j)=par%thetanum*s%ceqbg(i,j,jg)+(1.d0-par%thetanum)& *s%ceqbg(min(i+1,s%nx),j,jg) !cub(i,j)=par%thetanum*ccb(i,j)+(1.d0-par%thetanum)*ccb(min(i+1,nx),j) elseif(s%ureps(i,j)<0.d0) then cu(i,j)=par%thetanum*cc(i+1,j)+(1.d0-par%thetanum)*cc(max(i,2),j) cub(i,j)=par%thetanum*s%ceqbg(i+1,j,jg)+(1.d0-par%thetanum)& *s%ceqbg(max(i,2),j,jg) !cub(i,j)=par%thetanum*ccb(i+1,j)+(1.d0-par%thetanum)*ccb(max(i,2),j) else cu(i,j)=0.5d0*(cc(i,j)+cc(i+1,j)) cub(i,j)=0.5d0*(s%ceqbg(i,j,jg)+s%ceqbg(i+1,j,jg)) !cub(i,j)=0.5d0*(ccb(i,j)+ccb(i+1,j)) endif s%dcsdx(i,j)=(cc(i+1,j)-cc(i,j))/s%dsu(i,j) enddo enddo ! wwvv dcdx(nx:1,:) is still untouched, correct this ofr the parallel case cu(s%nx+1,:) = cc(s%nx+1,:) !Robert ! wwvv fix this in parallel case ! wwvv in parallel version, there will be a discrepancy between the values ! of dzbdx(nx+1,:). !wwvv so fix that ! Sus = 0.d0 Sub = 0.d0 ! ! suspended load, Lodewijk: no bed slope effect (yet) Sus=par%sus*(cu*s%ureps*s%hu-s%Dc*s%hu*s%dcsdx)*s%wetu ! bed load, Lodewijk: no bed slope effect (yet) Sub=par%bed*(cub*s%urepb*s%hu)*s%wetu ! ! Originally bed slope effect of XBeach : par%bdslpeffmag = 1 ! Original one, but only on bed load : par%bdslpeffmag = 2 if (par%bdslpeffmag == BDSLPEFFMAG_ROELV_TOTAL) then if (par%bermslope>0) then where (s%H/s%hu>1.0d0.or.(s%hu<1.d0.and. par%wavemodel == WAVEMODEL_STATIONARY)) Sus = Sus-par%sus*(10.d0*par%facsl*cu*s%vmagu*s%hu*(s%dzbdx-par%bermslope))*s%wetu elsewhere Sus = Sus-par%sus*(par%facsl*cu*s%vmagu*s%hu*s%dzbdx)*s%wetu end where else Sus = Sus-par%sus*(par%facsl*cu*s%vmagu*s%hu*s%dzbdx)*s%wetu endif endif if (par%bdslpeffmag == BDSLPEFFMAG_ROELV_TOTAL .or. par%bdslpeffmag == BDSLPEFFMAG_ROELV_BED) then if (par%bermslope>0) then where (s%H/s%hu>1.0d0 .or. (s%hu<1.d0.and. par%wavemodel == WAVEMODEL_STATIONARY)) Sub = Sub-par%bed*(10.0d0*par%facsl*cub*s%vmagu*s%hu*(s%dzbdx-par%bermslope))*s%wetu elsewhere Sub = Sub-par%bed*(par%facsl*cub*s%vmagu*s%hu*s%dzbdx)*s%wetu end where else Sub = Sub-par%bed*(par%facsl*cub*s%vmagu*s%hu*s%dzbdx)*s%wetu endif endif ! ! ! Y-direction ! ! Jaap: get ua in v points and split out in u and v direction if (s%ny>0) then uau(:,1:s%ny) = 0.5*(s%ua(:,1:s%ny)*costhm(:,1:s%ny)+s%ua(:,2:s%ny+1)*costhm(:,1:s%ny)) uav(:,1:s%ny) = 0.5*(s%ua(:,1:s%ny)*sinthm(:,1:s%ny)+s%ua(:,2:s%ny+1)*sinthm(:,1:s%ny)) uau(:,s%ny+1) = uau(:,s%ny) ! Jaap uav(:,s%ny+1) = uav(:,s%ny) ! Jaap if (par%nz>1) then vev_sed(:,1:s%ny) = 0.5*(s%ve_sed(:,1:s%ny)+s%ve_sed(:,2:s%ny+1)) else vev_sed=s%vev endif uev_sed(:,1:s%ny) = 0.5*(s%ue_sed(:,1:s%ny)+s%ue_sed(:,2:s%ny+1)) if (xmpi_isright) then uev_sed(:,1:s%ny+1) = uev_sed(:,1:s%ny) endif else uau=s%ua*costhm uav=s%ua*sinthm uev_sed=s%ue_sed vev_sed=s%ve_sed endif ! Jaap: compute vmagv including ua ! s%vmagv = sqrt((s%uv+uau)**2+(s%vv+uav)**2) s%vmagv = sqrt((uev_sed+uau)**2+(vev_sed+uav)**2) ! sediment advection velocity for suspended load and bed load respectively ! REMARK: when vreps does not equal vv; no mass conservation ! s%vreps = s%vev+uav ! s%vrepb = s%vev+uav ! RJ maybe reduce this velocity? Should be s%vv instead of s%vev? s%vreps = vev_sed+uav s%vrepb = vev_sed+uav ! RJ maybe reduce this velocity? Should be s%vv instead of s%vev? ! if (s%ny>0) then do j=1,s%ny do i=1,s%nx+1 if(s%vreps(i,j)>0) then cv(i,j)=par%thetanum*cc(i,j)+(1.d0-par%thetanum)*cc(i,min(j+1,s%ny)) cvb(i,j)=par%thetanum*s%pbbed(i,j,1,jg)*s%ceqbg(i,j,jg)+(1.d0-par%thetanum)& *s%pbbed(i,min(j+1,s%ny),1,jg)*s%ceqbg(i,min(j+1,s%ny),jg) !cvb(i,j)=par%thetanum*ccb(i,j)+(1.d0-par%thetanum)*ccb(i,min(j+1,ny)) elseif(s%vreps(i,j)<0) then cv(i,j)=par%thetanum*cc(i,j+1)+(1.d0-par%thetanum)*cc(i,max(j,2)) cvb(i,j)=par%thetanum*s%pbbed(i,j+1,1,jg)*s%ceqbg(i,j+1,jg)+(1.d0-par%thetanum)& *s%pbbed(i,max(j,2),1,jg)*s%ceqbg(i,max(j,2),jg) !cvb(i,j)=par%thetanum*ccb(i,j+1)+(1.d0-par%thetanum)*ccb(i,max(j,2)) else cv(i,j)=0.5d0*(cc(i,j)+cc(i,j+1)) !Jaap: cc instead of cv cvb(i,j)=0.5d0*(s%pbbed(i,j,1,jg)*s%ceqbg(i,j,jg)+s%pbbed(i,j+1,1,jg)*s%ceqbg(i,j+1,jg)) !cvb(i,j)=0.5d0*(ccb(i,j)+ccb(i,j+1)) end if s%dcsdy(i,j)=(cc(i,j+1)-cc(i,j))/s%dnv(i,j) !Jaap end do end do ! wwvv dcdy(:,ny+1) is not filled in, so in parallel case: cv(:,s%ny+1) = cc(:,s%ny+1) !Robert ! wwvv in parallel version, there will be a discrepancy between the values ! of dzbdy(:,ny+1). ! wwvv so fix that else cv = cc cvb = s%ceqbg(:,:,jg) endif ! s%ny>0 ! ! Compute sedimnent transport in v-direction ! Svs = 0.d0 Svb = 0.d0 ! Suspended load Svs=par%sus*(cv*s%vreps*s%hv-s%Dc*s%hv*s%dcsdy)*s%wetv ! Bed load Svb=par%bed*(cvb*s%vrepb*s%hv)*s%wetv ! ! Originally bed slope magnitude effect of XBeach : par%bdslpeffmag = 1 ! Original one, but only on bed load : par%bdslpeffmag = 2 if (par%bdslpeffmag == BDSLPEFFMAG_ROELV_TOTAL) then Svs = Svs-par%sus*(par%facsl*cv*s%vmagv*s%hv*s%dzbdy)*s%wetv endif if (par%bdslpeffmag == BDSLPEFFMAG_ROELV_TOTAL .or. par%bdslpeffmag == BDSLPEFFMAG_ROELV_BED ) then Svb = Svb-par%bed*(par%facsl*cvb*s%vmagv*s%hv*s%dzbdy)*s%wetv endif ! ! ! Bed slope magnitude effect (as Souslby intended) and change direction transport (see Van Rijn 1993 (section 7.2.6)) ! ! if (par%bdslpeffmag == BDSLPEFFMAG_SOULS_TOTAL .or. par%bdslpeffmag == BDSLPEFFMAG_SOULS_BED) then do j=1,s%ny+1 do i=1,s%nx+1 if ((dabs(Sub(i,j)) > 0.000001d0) .or. (dabs(Svb(i,j)) > 0.000001d0)) then Sbmtot = dsqrt( Sub(i,j)**2.d0 + Svb(i,j)**2.d0 ) dzbds = s%dzbdx(i,j)*Sub(i,j)/Sbmtot + s%dzbdy(i,j)*Svb(i,j)/Sbmtot ! dzbdn = s%dzbdx*Svb/Sbtot + s%dzbdy*Sub/Sbtot Sub(i,j) = Sub(i,j)*(1.d0 - par%facsl*dzbds) Svb(i,j) = Svb(i,j)*(1.d0 - par%facsl*dzbds) ! endif if (((dabs(Sus(i,j)) > 0.000001d0) .or. (dabs(Svs(i,j)) > 0.000001d0)) & .and. par%bdslpeffmag == BDSLPEFFMAG_SOULS_TOTAL) then Ssmtot = dsqrt( Sus(i,j)**2.d0 + Svs(i,j)**2.d0 ) dzbds = s%dzbdx(i,j)*Sus(i,j)/Ssmtot + s%dzbdy(i,j)*Svs(i,j)/Ssmtot ! dzbdn = s%dzbdx*Svb/Sbtot + s%dzbdy*Sub/Sbtot Sus(i,j) = Sus(i,j)*(1.d0 - par%facsl*dzbds) Svs(i,j) = Svs(i,j)*(1.d0 - par%facsl*dzbds) endif enddo enddo endif ! ! Lodewijk: modify the direction of the bed load transport based on the bed slope, see Van Rijn 1993 (section 7.2.6) if (par%bdslpeffdir == BDSLPEFFDIR_TALMON) then do j=1,s%ny+1 do i=1,s%nx+1 if (((dabs(s%urepb(i,j)) > 0.0001d0) .or. (dabs(s%vrepb(i,j)) > 0.0001d0)) & .and. ((dabs(s%taubx(i,j)) > 0.0001d0) .or. (dabs(s%tauby(i,j)) > 0.0001d0))) then if (s%urepb(i,j) < 0.d0) then delta_x = datan(s%vrepb(i,j)/s%urepb(i,j))+par%px ! Angle between fluid velocity vector and the s%x-axis else delta_x = datan(s%vrepb(i,j)/s%urepb(i,j)) ! Angle between fluid velocity vector and the s%x-axis endif delta = (par%rhos-par%rho)/par%rho shields = sqrt(s%taubx(i,j)**2 + s%tauby(i,j)**2)/(delta*par%rho*par%g*s%D50(jg)) ! shields = (urepb(i,j)**2.d0+vrepb(i,j)**2.d0)*s%cf(i,j)/(par%g*D50(jg)*delta) ftheta = 1.d0/(9.d0*(s%D50(jg)/s%hh(i,j))**0.3d0*sqrt(shields)) ! Talmon psi_x = datan( (dsin(delta_x)-ftheta*s%dzbdy(i,j)) / (dcos(delta_x)-ftheta*s%dzbdx(i,j)) ) psi_x = par%bdslpeffdirfac*(psi_x - delta_x)+delta_x Sbtot = dsqrt( Sub(i,j)**2.d0 + Svb(i,j)**2.d0 ) ! Magnitude of sediment transport without direction modifcation ! Decompose the sediment transport again, know with the knowledge of the direction of the sediment transport vector Sub(i,j) = Sbtot * dcos(psi_x) Svb(i,j) = Sbtot * dsin(psi_x) else Sub(i,j) = 0.d0 Svb(i,j) = 0.d0 endif enddo enddo endif ! ! ! do j=1,s%ny+1 do i=1,s%nx if(Sub(i,j)>0.d0) then pbbedu(i,j) = s%pbbed(i,j,1,jg) elseif(Sub(i,j)<0.d0) then pbbedu(i,j)= s%pbbed(i+1,j,1,jg) else pbbedu(i,j)=0.5d0*(s%pbbed(i,j,1,jg)+s%pbbed(i+1,j,1,jg)) endif enddo enddo ! !Sub = pbbedu*Sub ! do j=1,s%ny do i=1,s%nx+1 if(Svb(i,j)>0) then pbbedv(i,j)=s%pbbed(i,j,1,jg) else if(Svb(i,j)<0) then pbbedv(i,j)=s%pbbed(i,j+1,1,jg) else pbbedv(i,j)=0.5d0*(s%pbbed(i,j,1,jg)+s%pbbed(i,j+1,1,jg)) end if end do end do ! Svb = pbbedv*Svb ! ! BRJ: implicit concentration update (compute sources first, sink must be computed after updating actual sed.conc.) ! if (s%ny>0) then do j=2,s%ny do i=2,s%nx ! Changed to hh from hold by RJ (13072009) !**2/max(hh(i,j),par%hmin) s%ero(i,j,jg) = fac(i,j,jg)*s%hh(i,j)*s%ceqsg(i,j,jg)*s%pbbed(i,j,1,jg)/s%Tsg(i,j,jg) ! depo_ex(i,j,jg) = max(hold(i,j),0.01d0)*cc(i,j)/Tsg(i,j,jg) ! BRJ: the volume in the water column is updated and not the volume concentration. cc(i,j) = (par%dt*s%Tsg(i,j,jg))/(par%dt+s%Tsg(i,j,jg))* & (s%hold(i,j)*cc(i,j)/par%dt -((Sus(i,j)*s%dnu(i,j)-Sus(i-1,j)*s%dnu(i-1,j)+& Svs(i,j)*s%dsv(i,j)-Svs(i,j-1)*s%dsv(i,j-1))*s%dsdnzi(i,j)-& s%ero(i,j,jg))) cc(i,j)=max(cc(i,j),0.0d0) ! Jaap: negative cc's are possible... cc(i,j)=min(cc(i,j),par%cmax*s%hh(i,j)) s%depo_ex(i,j,jg) = cc(i,j)/s%Tsg(i,j,jg) enddo enddo else j=1 do i=2,s%nx ! depthcheck=s%z_g(i,j)-s%overlap(i,j) ! if(depthcheck>par%depth_thr) then !threshold ! s%trans_inf(i,j)=0 ! elseif (depthcheck>0 .and. depthcheck<=par%depth_thr) then ! s%trans_inf(i,j)=1-depthcheck/par%depth_thr ! else ! s%trans_inf(i,j)=1 ! end if ! ! s%depo_ex(i,j,jg) = cc(i,j)/s%Tsg(i,j,jg) ! s%ero(i,j,jg) = fac(i,j,jg)*s%hh(i,j)*s%ceqsg(i,j,jg)*s%trans_inf(i,j)/s%Tsg(i,j,jg) ! if(s%ero(i,j,jg)>s%depo_ex(i,j,jg)) then ! s%ero(i,j,jg)=min(s%depo_ex(i,j,jg),s%ero(i,j,jg)) ! endif !! depo_ex(i,j,jg) = max(hold(i,j),0.01d0)*cc(i,j)/Tsg(i,j,jg) !! BRJ: the volume in the water column is updated and not the volume concentration. !cc(i,j) = (par%dt*s%Tsg(i,j,jg))/(par%dt+s%Tsg(i,j,jg))* & !(s%hold(i,j)*cc(i,j)/par%dt -((Sus(i,j)*s%dnu(i,j)-Sus(i-1,j)*s%dnu(i-1,j))*s%dsdnzi(i,j)-& !s%ero(i,j,jg))) ! !cc(i,j)=max(cc(i,j),0.0d0) ! Jaap: negative cc's are possible... !cc(i,j)=min(cc(i,j),par%cmax*s%hh(i,j)) ! !s%depo_ex(i,j,jg) = cc(i,j)/s%Tsg(i,j,jg) ! Changed to hh from hold by RJ (13072009) !**2/max(hh(i,j),par%hmin) s%ero(i,j,jg) = fac(i,j,jg)*s%hh(i,j)*s%ceqsg(i,j,jg)*s%pbbed(i,j,1,jg)/s%Tsg(i,j,jg) ! depo_ex(i,j,jg) = max(hold(i,j),0.01d0)*cc(i,j)/Tsg(i,j,jg) ! BRJ: the volume in the water column is updated and not the volume concentration. cc(i,j) = (par%dt*s%Tsg(i,j,jg))/(par%dt+s%Tsg(i,j,jg))* & (s%hold(i,j)*cc(i,j)/par%dt -((Sus(i,j)*s%dnu(i,j)-Sus(i-1,j)*s%dnu(i-1,j))*s%dsdnzi(i,j)-& s%ero(i,j,jg))) cc(i,j)=max(cc(i,j),0.0d0) ! Jaap: negative cc's are possible... cc(i,j)=min(cc(i,j),par%cmax*s%hh(i,j)) s%depo_ex(i,j,jg) = cc(i,j)/s%Tsg(i,j,jg) ! s%depo_ex(i,j,jg) = cc(i,j)/s%Tsg(i,j,jg) ! s%ero(i,j,jg) = fac(i,j,jg)*s%hh(i,j)*s%ceqsg(i,j,jg)*s%pbbed(i,j,1,jg)/s%Tsg(i,j,jg) ! if (s%ero(i,j,jg)>depo ! s%ero(i,j,jg)=minval(dep,ero) enddo endif cc = cc/s%hh ! do lateral boundaries... if(xmpi_istop)then cc(1,:)=cc(2,:) s%ero(1,:,jg)=s%ero(2,:,jg) s%depo_ex(1,:,jg)=s%depo_ex(2,:,jg) endif if(xmpi_isleft .and. s%ny>0)then cc(:,1)=cc(:,2) s%ero(:,1,jg)=s%ero(:,2,jg) s%depo_ex(:,1,jg)=s%depo_ex(:,2,jg) endif if(xmpi_istop)then cc(s%nx+1,:)=cc(s%nx+1-1,:) s%ero(s%nx+1,:,jg)=s%ero(s%nx,:,jg) s%depo_ex(s%nx+1,:,jg)=s%depo_ex(s%nx,:,jg) endif if(xmpi_isright .and. s%ny>0)then cc(:,s%ny+1)=cc(:,s%ny+1-1) s%ero(:,s%ny+1,jg)=s%ero(:,s%ny,jg) s%depo_ex(:,s%ny+1,jg)=s%depo_ex(:,s%ny,jg) endif ! wwvv fix the first and last rows and columns of cc in parallel case #ifdef USEMPI call xmpi_shift_ee(cc) #endif ! ! wwvv border columns and rows of ccg Svg and Sug have to be communicated s%ccg(:,:,jg) = cc s%Svsg(:,:,jg) = Svs s%Susg(:,:,jg) = Sus s%Svbg(:,:,jg) = Svb s%Subg(:,:,jg) = Sub ! number of sediment fractions s%vmag=sqrt(max(vmag2,par%umin)) end subroutine transus_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine bed_update_conceptmodel(s,par) use params use spaceparams use xmpi_module use math_tools use paramsconst use morphevolution use interp use morphevolution implicit none type(spacepars),target :: s type(parameters) :: par integer :: i,j,j1,jg,ii,ie,id,je,jd,jdz,ndz, hinterland integer , dimension(:,:,:),allocatable,save :: indSus,indSub,indSvs,indSvb real*8 :: dzb,dzmax,dzt,dzleft,sdz,dzavt,fac,Savailable,dAfac real*8 , dimension(:,:),allocatable,save :: Sout,hav,tempexchange real*8 , dimension(par%ngd) :: edg,edg1,edg2,dzg real*8 , dimension(:),pointer :: dz real*8 , dimension(:,:),pointer :: pb integer :: n_aval real*8,save :: delta real*8 , dimension(par%ngd) :: flux_sand,flux_gravel, flux_subsand !Huub: concept model real*8 , dimension(:,:),allocatable,save :: dz_g !Huub: concept model real*8 , dimension(:,:),allocatable,save :: dz_s !Huub: concept model real*8 , dimension(:,:),allocatable,save :: q_s, q_g, q_ss !Huub: concept model real*8 , dimension(:,:,:),allocatable,save :: Suss_init, Susg_init !Huub: concept model real*8 :: depthcheck, z_gnew, dz_sextra, zb_init !Huub: concept model real*8 :: overlap_space, overlap_init ,counter_init, mass1, mass2, mass_check, sign_check !Huub: concept model !include 's.ind' !include 's.inp' if (.not. allocated(Sout)) then allocate(Sout(s%nx+1,s%ny+1)) allocate(hav(s%nx+1,s%ny+1)) allocate(indSus(s%nx+1,s%ny+1,par%ngd)) allocate(indSub(s%nx+1,s%ny+1,par%ngd)) allocate(indSvs(s%nx+1,s%ny+1,par%ngd)) allocate(indSvb(s%nx+1,s%ny+1,par%ngd)) allocate(q_g(s%nx+1,s%ny+1)) !Huub: concept model; allocate(q_s(s%nx+1,s%ny+1)) !Huub: concept model; allocate(q_ss(s%nx+1,s%ny+1)) !Huub: concept model; allocate(dz_g(s%nx+1,s%ny+1)) !Huub: concept model; allocate(dz_s(s%nx+1,s%ny+1)) !Huub: concept model; allocate(Suss_init(s%nx+1,s%ny+1,par%ngd)) !Huub: concept model; allocate(Susg_init(s%nx+1,s%ny+1,par%ngd)) !Huub: concept model; if (par%gwflow==1) then allocate(tempexchange(s%nx+1,s%ny+1)) endif delta = (par%rhos-par%rho)/par%rho endif ! Super fast 1D if (s%ny==0) then j1 = 1 else j1 = 2 endif s%dzbnow = 0.d0 dzb = 0.d0 if (par%t>=par%morstart .and. par%t < par%morstop .and. par%morfac > .999d0) then jg=1 do j=1,s%ny+1 do i=2,s%nx if(s%ureps(i,j)>0.d0) then depthcheck=s%z_g(i,j)-s%overlap(i,j) if(depthcheck>par%depth_thr) then s%trans_inf(i,j)=0 elseif (depthcheck>0 .and. depthcheck<=par%depth_thr) then s%trans_inf(i,j)=1-depthcheck/par%depth_thr else s%trans_inf(i,j)=1 endif s%Susg(i,j,jg)=s%Susg(i,j,jg)*s%trans_inf(i,j) s%Subg(i,j,jg)=s%Subg(i,j,jg)*s%trans_inf(i,j) s%Subg_g(i,j,2)=s%Subg_g(i,j,2)*(1-s%trans_inf(i,j)) elseif (s%ureps(i,j)<0.d0) then depthcheck=s%z_g(i+1,j)-s%overlap(i+1,j) if(depthcheck>par%depth_thr) then s%trans_inf(i,j)=0 elseif (depthcheck>0 .and. depthcheck<=par%depth_thr) then s%trans_inf(i,j)=1-depthcheck/par%depth_thr else s%trans_inf(i,j)=1 endif s%Susg(i,j,jg)=s%Susg(i,j,jg)*s%trans_inf(i,j) s%Subg(i,j,jg)=s%Subg(i,j,jg)*s%trans_inf(i,j) s%Subg_g(i,j,2)=s%Subg_g(i,j,2)*(1-s%trans_inf(i,j)) endif if (s%gwu(i,j)>0.d0) then depthcheck=s%z_g(i,j)-s%overlap(i,j) if(depthcheck>par%depth_thr) then s%trans_inf_gwf(i,j)=0 elseif (depthcheck>0 .and. depthcheck<=par%depth_thr) then s%trans_inf_gwf(i,j)=1-depthcheck/par%depth_thr else s%trans_inf_gwf(i,j)=1 endif s%Suss(i,j,jg)=s%Suss(i,j,jg)*(1-s%trans_inf_gwf(i,j)) elseif (s%gwu(i,j)<0.d0) then depthcheck=s%z_g(i+1,j)-s%overlap(i+1,j) if(depthcheck>par%depth_thr) then s%trans_inf_gwf(i,j)=0 elseif (depthcheck>0 .and. depthcheck<=par%depth_thr) then s%trans_inf_gwf(i,j)=1-depthcheck/par%depth_thr else s%trans_inf_gwf(i,j)=1 endif s%Suss(i,j,jg)=s%Suss(i,j,jg)*(1-s%trans_inf_gwf(i,j)) endif enddo enddo !Huub: make Suss more diffuse Suss_init=s%Suss do j=1,s%ny+1 do i=2,s%nx if(Suss_init(i,j,jg)>0.d0) then s%Suss(i,j,jg)=0.25d0*Suss_init(i-1,j,jg)+0.5d0*Suss_init(i,j,jg)+0.25d0*Suss_init(i+1,j,jg) else s%Suss(i,j,jg)=0.25d0*Suss_init(i-1,j,jg)+0.5d0*Suss_init(i,j,jg)+0.25d0*Suss_init(i+1,j,jg) endif enddo enddo !!Huub: make Susg more diffuse !Susg_init=s%Susg ! ! do j=1,s%ny+1 ! do i=2,s%nx ! if(Susg_init(i,j,jg)>0.d0) then ! s%Susg(i,j,jg)=0.25d0*Susg_init(i-1,j,jg)+0.5d0*Susg_init(i,j,jg)+0.25d0*Susg_init(i+1,j,jg) ! else ! s%Susg(i,j,jg)=0.25d0*Susg_init(i-1,j,jg)+0.5d0*Susg_init(i,j,jg)+0.25d0*Susg_init(i+1,j,jg) ! endif ! enddo ! enddo ! ! bed_predict ! ! reduce sediment transports when hard layer comes to surface ! this step is mainly necessary at the transition from hard layers to sand !if (par%struct == 1) then ! Huub: is zero for concept modelB ! jg = 1 indSus = 0 indSub = 0 indSvs = 0 indSvb = 0 Sout = 0.d0 do j=j1,s%ny+1 do i=2,s%nx+1 ! ! fluxes at i,j ! if (s%Subg(i,j,jg) > 0.d0) then ! bed load s%u-direction ! indSub(i,j,jg) = 1 ! Sout(i,j) = Sout(i,j) + s%Subg(i,j,jg)*s%dnu(i,j) ! endif ! ! fluxes at i-1,j ! if (s%Subg(i-1,j,jg) < 0.d0 ) then ! bed load s%u-direction ! Sout(i,j) = Sout(i,j) - s%Subg(i-1,j,jg)*s%dnu(i-1,j) ! endif ! fluxes at i,j if (s%Subg_g(i,j,2) > 0.d0) then ! bed load s%u-direction indSub(i,j,2) = 1 Sout(i,j) = Sout(i,j) + s%Subg_g(i,j,2)*s%dnu(i,j) endif ! fluxes at i-1,j if (s%Subg_g(i-1,j,2) < 0.d0 ) then ! bed load s%u-direction Sout(i,j) = Sout(i,j) - s%Subg_g(i-1,j,2)*s%dnu(i-1,j) endif ! if (par%sourcesink==0) then ! ! fluxes at i,j ! if (s%Susg(i,j,jg) > 0.d0 ) then ! suspended load s%u-direction ! indSus(i,j,jg) = 1 ! Sout(i,j) = Sout(i,j) + s%Susg(i,j,jg)*s%dnu(i,j) ! endif ! ! fluxes at i-1,j ! if (s%Susg(i-1,j,jg) < 0.d0 ) then ! suspended load s%u-direction ! Sout(i,j) = Sout(i,j) - s%Susg(i-1,j,jg)*s%dnu(i-1,j) ! endif ! endif enddo enddo ! if (s%ny>0) then ! do j=j1,s%ny+1 ! do i=2,s%nx+1 ! if (s%Svbg(i,j,jg) > 0.d0 ) then ! bed load s%v-direction ! indSvb(i,j,jg) = 1 ! Sout(i,j) = Sout(i,j) + s%Svbg(i,j,jg)*s%dsv(i,j) ! endif ! ! fluxes at i,j-1 ! if (s%Svbg(i,j-1,jg) < 0.d0 ) then ! bed load s%v-direction ! Sout(i,j) = Sout(i,j) - s%Svbg(i,j-1,jg)*s%dsv(i,j-1) ! endif ! if (par%sourcesink==0) then ! if (s%Svsg(i,j,jg ) > 0.d0 ) then ! suspended load s%v-direction ! indSvs(i,j,jg) = 1 ! Sout(i,j) = Sout(i,j) + s%Svsg(i,j,jg)*s%dsv(i,j) ! endif ! ! fluxes at i,j-1 ! if (s%Svsg(i,j-1,jg) < 0.d0 ) then ! suspended load s%v-direction ! Sout(i,j) = Sout(i,j) - s%Svsg(i,j-1,jg)*s%dsv(i,j-1) ! endif ! endif ! sourcesink = 0 ! enddo !s%nx+1 ! enddo !s%ny+1 ! endif !s%ny>0 ! ! do j=j1,s%ny+1 do i=2,s%nx+1 Savailable = s%z_g(i,j)/par%morfac/par%dt*(1.d0-par%por)/s%dsdnzi(i,j) ! reduction factor for cell outgoing sediment transports wwvv changed tiny into epsilon fac = min(1.d0,Savailable/max(Sout(i,j),epsilon(0.d0)) ) ! fix sediment transports for the presence of a hard layer; remind indSus etc are 1 in cases of cell outgoing transports ! updated S oell outgoing transports cell incoming transports if (fac < 1.d0)then s%Subg_g(i,j,2) = fac*indSub(i,j,2)*s%Subg_g(i,j,2) + (1-indSub(i,j,2))*s%Subg_g(i,j,2) s%Subg_g(i-1,j,2) = fac*(1-indSub(i-1,j,2))*s%Subg_g(i-1,j,2) + indSub(i-1,j,2)*s%Subg_g(i-1,j,2) ! s%Subg(i,j,jg) = fac*indSub(i,j,jg)*s%Subg(i,j,jg) + (1-indSub(i,j,jg))*s%Subg(i,j,jg) ! s%Subg(i-1,j,jg) = fac*(1-indSub(i-1,j,jg))*s%Subg(i-1,j,jg) + indSub(i-1,j,jg)*s%Subg(i-1,j,jg) ! if (s%ny>0) then ! s%Svbg(i,j,jg) = fac*indSvb(i,j,jg)*s%Svbg(i,j,jg) + (1-indSvb(i,j,jg))*s%Svbg(i,j,jg) ! s%Svbg(i,j-1,jg) = fac*(1-indSvb(i,j-1,jg))*s%Svbg(i,j-1,jg) + indSvb(i,j-1,jg)*s%Svbg(i,j-1,jg) ! endif ! if (par%sourcesink==0) then ! s%Susg(i,j,jg) = fac*indSus(i,j,jg)*s%Susg(i,j,jg) + (1-indSus(i,j,jg))*s%Susg(i,j,jg) ! s%Susg(i-1,j,jg) = fac*(1-indSus(i-1,j,jg))*s%Susg(i-1,j,jg) + indSus(i-1,j,jg)*s%Susg(i-1,j,jg) ! if (s%ny>0) then ! s%Svsg(i,j,jg) = fac*indSvs(i,j,jg)*s%Svsg(i,j,jg) + (1-indSvs(i,j,jg))*s%Svsg(i,j,jg) ! s%Svsg(i,j-1,jg) = fac*(1-indSvs(i,j-1,jg))*s%Svsg(i,j-1,jg) + indSvs(i,j-1,jg)*s%Svsg(i,j-1,jg) ! endif !s%ny = 0 ! endif ! sourcesink = 0 endif !fac<1.d0 enddo ! s%nx+1 enddo !s%ny + 1 !endif !struct == 1 ! j=1 do i=2, s%nx !Huub: loop to identify fluxes, which is used to see if fluxes are positive or negative. flux_sand=( s%Susg(i,j,:)*s%dnu(i,j)-s%Susg(i-1,j,:)*s%dnu(i-1,j) + s%Subg(i,j,:)*s%dnu(i,j) & -s%Subg(i-1,j,:)*s%dnu(i-1,j) + s%Suss(i,j,:)*s%dnu(i,j)-s%Suss(i-1,j,:)*s%dnu(i-1,j))*s%dsdnzi(i,j) !Huub: transport gradient of sand, see page 128 Kingsday release !flux_sand=(s%Suss(i,j,:)*s%dnu(i,j)-s%Suss(i-1,j,:)*s%dnu(i-1,j))*s%dsdnzi(i,j) !Huub: transport gradient of sand, see page 128 Kingsday release q_s(i,j)=flux_sand(1)!*s%trans_inf(i,j) !No transport for trans_inf==0 flux_gravel=(s%Subg_g(i,j,2)*s%dnu(i,j)-s%Subg_g(i-1,j,2)*s%dnu(i-1,j))*s%dsdnzi(i,j) !Huub: gravelflux same numerical implementation as for sand. q_g(i,j)=flux_gravel(1)!*(1-s%trans_inf(i,j)) !Only transport for trans_inf<1 !q_g(i,j)=0 flux_subsand=(s%Suss(i,j,:)*s%dnu(i,j)-s%Suss(i-1,j,:)*s%dnu(i-1,j))*s%dsdnzi(i,j) ! flux_subsand=(s%Sus_s(i,j,:)*s%dnu(i,j)-s%Sus_s(i-1,j,:)*s%dnu(i-1,j)) q_ss(i,j)=flux_subsand(1) !q_ss(i,j)=0 end do do i=2, s%nx !Huub: first calculate bed level change of gravel dz_g(i,j)=(par%morfac*par%dt/(1.d0-par%por_fr2))*(q_g(i,j)) !Normal exner with gravel porosity z_gnew=s%z_g(i,j) - dz_g(i,j) if(z_gnew<0)then dz_g(i,j)=s%z_g(i,j) endif z_gnew=s%z_g(i,j) - dz_g(i,j) s%z_g(i,j)=max(z_gnew,0.0d0) if (s%counter(i,j)>0) then !Calculations to put gravel underneath sand // z_g and overlap must change // z_s also changes counter_init=s%counter(i,j) !can go away? mass1=counter_init*(1-par%por)/par%morfac/par%dt mass2=-dz_g(i,j)*par%por_fr2*(1-par%por)/par%morfac/par%dt !mass of sand that could fit into gravel // -dz_g, because accretion gives negative dz_g's mass_check=mass1-mass2 if (mass_check>0 .and. dz_g(i,j)<0) then !HUUB: FIXFIXFIX (done 11/01 // 17/01 dz_g<0 because accretion is negative. ) dz_sextra=dz_g(i,j) + ((mass_check)/(1-par%por))*par%morfac*par%dt ! dz_g + . mass_check.. --> because accretion flux is reduced and dz_g is negative s%counter(i,j)=((mass_check)/(1.d0-par%por))*par%morfac*par%dt s%overlap(i,j)=s%z_g(i,j) elseif (mass_check<=0 .and. dz_g(i,j)<0) then !masscheck<0 or ==0 Huub: 17/01 dz_g<0 because accretion is negative. dz_sextra=-(mass1*par%morfac*par%dt/(1-par%por))/par%por_fr2 !HUUB: CHECK (done 11/01 done 17/01) s%counter(i,j)=0 s%overlap(i,j)=s%z_g(i,j)+par%morfac*par%dt*mass_check/((1-par%por)*par%por_fr2) !HUUB: CHECK (done 11/01 ()?) else dz_sextra=0 endif s%z_s(i,j)=s%z_s(i,j)-dz_sextra !Huub: CHECK!!! endif depthcheck=(s%z_g(i,j)-s%overlap(i,j)) overlap_space=depthcheck if (depthcheck>par%depth_thr) then !quite some gravel present if (q_s(i,j)<0) then !negative is accretion dz_s(i,j)=(par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j)))/par%por_fr2 s%overlap(i,j)=s%overlap(i,j)-dz_s(i,j) depthcheck=(s%z_g(i,j)-s%overlap(i,j)) if (depthcheck<0) then dz_s(i,j)=par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j)+overlap_space*par%por_fr2*(1-par%por)/par%morfac/par%dt & )-overlap_space !Zie afleiding Huub s%overlap(i,j)=s%z_g(i,j) s%counter(i,j)=dz_s(i,j)-overlap_space endif elseif (q_s(i,j)>0) then !positive is erosion (kan ook door randvoorwaarden) dz_s(i,j)=(par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j)))/par%por_fr2 s%overlap(i,j)=s%overlap(i,j)-dz_s(i,j) if (s%overlap(i,j)<0) then dz_s(i,j)=par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j)-(s%z_g(i,j)-overlap_space)*& par%por_fr2*(1-par%por)/par%morfac/par%dt) + (s%z_g(i,j)-overlap_space) s%overlap(i,j)=0 endif else !q_s==0 dz_s(i,j)=0 endif else !depthcheck =par%depth_thr of 0) then !to fill last bit of gravel sign_check=q_s(i,j) +depthcheck*par%por_fr2*(1-par%por)/par%morfac/par%dt if (sign_check>0) then !if flux whas not enough to fill gravel dz_s(i,j)=(par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j)))/par%por_fr2 s%overlap(i,j)=s%overlap(i,j)-dz_s(i,j) else dz_s(i,j)=par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j) +depthcheck*par%por_fr2*(1-par%por)/par%morfac/par%dt) - depthcheck s%overlap(i,j)=s%z_g(i,j) s%counter(i,j)=dz_s(i,j)-depthcheck endif else dz_s(i,j)=par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j)) if(s%z_g(i,j)==s%overlap(i,j) .and. s%z_g(i,j)>0) then !Huub; mogelijk toevoegen if counter > 0 s%counter(i,j)=s%counter(i,j)-dz_s(i,j) endif endif elseif (q_s(i,j)>0) then !erosion if(s%counter(i,j)==0 .and. s%overlap(i,j)==0) then dz_s(i,j)=par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j)) !elseif(s%counter(i,j)==0 .and. s%overlap(i,j)>0) then ! dz_s(i,j)=(par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j)))/par%por_fr2 ! s%overlap(i,j)=s%overlap(i,j)-dz_s(i,j) ! ! sign_check=s%z_g(i,j)-dz_s(i,j) ! ! if(sign_check<0) then ! dz_s(i,j)=par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j) - s%z_g(i,j)*par%por_fr2*(1-par%por)/par%morfac/par%dt) + s%z_g(i,j) ! s%overlap(i,j)=0 ! endif ! elseif (depthcheck>0) then ! tiny gravel present dz_s(i,j)=(par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j)))/par%por_fr2 if (dz_s(i,j)>s%overlap(i,j)) then dz_s(i,j)=(par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j)-s%overlap(i,j)*(1-par%por)*par%por_fr2/par%morfac/par%dt)) & !HUUB FIXEN + s%overlap(i,j) !HUUB FIXEN (18/01 gefixt.) s%overlap(i,j)=0 else s%overlap(i,j)=s%overlap(i,j)-dz_s(i,j) endif else !Huub; thus if counter>0 (and overlap =0(?) or overlap > 0) dz_s(i,j)=par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j)) if (dz_s(i,j)>s%counter(i,j)) then dz_s(i,j)=((par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j)-s%counter(i,j)*(1-par%por)/par%morfac/par%dt))/par%por_fr2) & !HUUB FIXEN + s%counter(i,j) !HUUB FIXEN (18/01 gefixt.) overlap_init=s%overlap(i,j) s%overlap(i,j)=s%overlap(i,j)-dz_s(i,j)-s%counter(i,j) if (s%overlap(i,j)<0) then dz_s(i,j)=par%morfac*par%dt/(1.d0-par%por)*(q_s(i,j)-s%counter(i,j)*(1-par%por)/par%morfac/par%dt-& overlap_init*par%por_fr2*(1-par%por)/par%morfac/par%dt) + overlap_init + s%counter(i,j) s%overlap(i,j)=0 endif s%counter(i,j)=0 else s%counter(i,j)=s%counter(i,j)-dz_s(i,j) endif endif else dz_s(i,j)=0 endif endif s%z_s(i,j)=s%z_s(i,j)-dz_s(i,j) zb_init=s%zb(i,j) s%zb(i,j)=s%z_s(i,j)+(s%z_g(i,j)-s%overlap(i,j)) s%dzbnow(i,j) = s%dzbnow(i,j)-(zb_init-s%zb(i,j)) s%dzbdt(i,j) = s%dzbnow(i,j)/par%dt ! s%sedero(i,j) = s%sedero(i,j)-sum(dzg) ! s%structdepth(i,j) = max(0.d0,s%structdepth(i,j)-sum(dzg)) enddo !if (s%ny>0) then ! do j=2,s%ny ! do i=2,s%nx ! ! ! bed level changes per fraction in this morphological time step in meters sand including pores ! ! positive in case of erosion ! if (par%sourcesink==0) then ! dzg=par%morfac*par%dt/(1.d0-par%por)*( & ! Dano, dz from sus transport gradients ! ( s%Susg(i,j,:)*s%dnu(i,j)-s%Susg(i-1,j,:)*s%dnu(i-1,j) +& ! s%Svsg(i,j,:)*s%dsv(i,j)-s%Svsg(i,j-1,:)*s%dsv(i,j-1) +& ! ! dz from bed load transport gradients ! s%Subg(i,j,:)*s%dnu(i,j)-s%Subg(i-1,j,:)*s%dnu(i-1,j)+& ! s%Svbg(i,j,:)*s%dsv(i,j)-s%Svbg(i,j-1,:)*s%dsv(i,j-1) )*s%dsdnzi(i,j) ) ! elseif (par%sourcesink==1) then ! dzg=par%morfac*par%dt/(1.d0-par%por)*( & ! s%ero(i,j,:)-s%depo_ex(i,j,:) +& ! ( s%Subg(i,j,:)*s%dnu(i,j)-s%Subg(i-1,j,:)*s%dnu(i-1,j)+& ! s%Svbg(i,j,:)*s%dsv(i,j)-s%Svbg(i,j-1,:)*s%dsv(i,j-1) )*s%dsdnzi(i,j) ) ! endif ! ! if (par%ngd==1) then ! Simple bed update in case one fraction ! ! s%zb(i,j) = s%zb(i,j)-sum(dzg) ! s%dzbnow(i,j) = s%dzbnow(i,j)-sum(dzg) ! naamgeveing? ! s%dzbdt(i,j) = s%dzbnow(i,j)/par%dt ! s%sedero(i,j) = s%sedero(i,j)-sum(dzg) ! s%structdepth(i,j) = max(0.d0,s%structdepth(i,j)-sum(dzg)) ! ! else ! ! erosion/deposition rate of sand mass (m/s) ! ! positive in case of erosion ! edg = dzg*(1.d0-par%por)/par%dt ! ! ! dz=>s%dzbed(i,j,:) ! pb=>s%pbbed(i,j,:,:) ! ! call update_fractions(par,s,i,j,dz,pb,edg,sum(dzg)) ! ! endif ! ! enddo ! s%nx+1 ! enddo ! s%ny+1 !else ! j=1 ! do i=2,s%nx ! ! bed level changes per fraction in this morphological time step in meters sand including pores ! ! positive in case of erosion ! if (par%sourcesink==0) then ! dzg=par%morfac*par%dt/(1.d0-par%por)*( & ! Dano, dz from sus transport gradients ! ( s%Susg(i,j,:)*s%dnu(i,j)-s%Susg(i-1,j,:)*s%dnu(i-1,j) +& ! s%Subg(i,j,:)*s%dnu(i,j)-s%Subg(i-1,j,:)*s%dnu(i-1,j) )*s%dsdnzi(i,j) +& ! (s%Svsg(i,j,:)+s%Svbg(i,j,:))*par%lsgrad) ! elseif (par%sourcesink==1) then ! dzg=par%morfac*par%dt/(1.d0-par%por)*( & ! s%ero(i,j,:)-s%depo_ex(i,j,:) +& ! ( s%Subg(i,j,:)*s%dnu(i,j)-s%Subg(i-1,j,:)*s%dnu(i-1,j) )*s%dsdnzi(i,j) +& ! (s%Svsg(i,j,:)+s%Svbg(i,j,:))*par%lsgrad) ! endif ! ! if (par%ngd==1) then ! Simple bed update in case one fraction ! ! s%zb(i,j) = s%zb(i,j)-sum(dzg) ! s%dzbnow(i,j) = s%dzbnow(i,j)-sum(dzg) ! s%dzbdt(i,j) = s%dzbnow(i,j)/par%dt ! s%sedero(i,j) = s%sedero(i,j)-sum(dzg) ! s%structdepth(i,j) = max(0.d0,s%structdepth(i,j)-sum(dzg)) ! ! else ! multiple fractions... ! ! erosion/deposition rate of sand mass (m/s) ! ! positive in case of erosion ! edg = dzg*(1.d0-par%por)/par%dt ! ! dz=>s%dzbed(i,j,:) ! pb=>s%pbbed(i,j,:,:) ! ! call update_fractions(par,s,i,j,dz,pb,edg,sum(dzg)) ! ! endif ! ! enddo ! s%nx+1 !endif !s%ny>0 #ifdef USEMPI call xmpi_shift_ee(s%zb) #endif ! In case of groundwater, we need to update groundwater levels and surface water if(par%gwflow==1) then where (s%wetz==1 .and. s%dzbnow>0.d0) ! accretion in wet areas s%zs = s%zs + s%dzbnow*(1.d0-par%por) ! zs = zs + dzbnow - dzbnow*par%por s%gwlevel = s%gwlevel+s%dzbnow s%infil = s%infil + s%dzbnow*par%por/par%dt elsewhere (s%wetz==1 .and. s%dzbnow<0.d0) ! erosion in wet areas ! maximum water leaving groundwater = gwlevel(not updated) - zb(updated) ! note that exfiltration is negative infil tempexchange = min((s%zb-s%gwlevel)*par%por,0.d0) s%zs = s%zs + s%dzbnow - tempexchange s%gwlevel = s%gwlevel + tempexchange/par%por s%infil = s%infil + tempexchange endwhere else ! only need to update water levels where (s%wetz==1) s%zs = s%zs+s%dzbnow endwhere endif !Huub: next part extracted from avalanching s%dzbdx=0.d0 s%dzbdy=0.d0 do j=1,s%ny+1 do i=1,s%nx s%dzbdx(i,j)=(s%zb(i+1,j)-s%zb(i,j))/s%dsu(i,j) enddo enddo do j=1,s%ny do i=1,s%nx+1 s%dzbdy(i,j)=(s%zb(i,j+1)-s%zb(i,j))/s%dnv(i,j) enddo enddo ! ! bed boundary conditions ! ! Dano: the following Neumann boundaries introduce a mass error equal to ! the bed level change due to the b.c. times the cell area. if(xmpi_isleft .and. s%ny>0) then s%zb(:,1) = s%zb(:,2) s%dzbdt(:,1) = s%dzbdt(:,2) s%dzbnow(:,1) = s%dzbnow(:,2) s%sedero(:,1) = s%sedero(:,2) s%structdepth(:,1) = s%structdepth(:,2) s%pbbed(:,1,:,:)=s%pbbed(:,2,:,:) s%z0bed(:,1)=s%z0bed(:,2) s%dzbed(:,1,:)=s%dzbed(:,2,:) s%z_s(:,1)= s%z_s(:,2) s%z_g(:,1)= s%z_g(:,2) s%overlap(:,1)= s%overlap(:,2) s%counter(:,1)= s%counter(:,2) endif if(xmpi_isright .and. s%ny>0) then s%zb(:,s%ny+1) = s%zb(:,s%ny) s%dzbdt(:,s%ny+1) = s%dzbdt(:,s%ny) s%dzbnow(:,s%ny+1) = s%dzbnow(:,s%ny) s%sedero(:,s%ny+1) = s%sedero(:,s%ny) s%structdepth(:,s%ny+1) = s%structdepth(:,s%ny) s%pbbed(:,s%ny+1,:,:)=s%pbbed(:,s%ny,:,:) s%z0bed(:,s%ny+1)=s%z0bed(:,s%ny) s%dzbed(:,s%ny+1,:)=s%dzbed(:,s%ny,:) s%z_s(:,s%ny+1)= s%z_s(:,s%ny) s%z_g(:,s%ny+1)= s%z_g(:,s%ny) s%overlap(:,s%ny+1)= s%overlap(:,s%ny) s%counter(:,s%ny+1)= s%counter(:,s%ny) endif ! Robert: in parallel version bed update must take place on internal boundaries: #ifdef USEMPI call xmpi_shift_ee(s%zb) #endif ! Update representative sed.diameter at the bed for flow friction and output if (par%ngd>1) then do j=j1,max(s%ny,1) do i=2,s%nx depthcheck=s%z_g(i,j)-s%overlap(i,j) if(depthcheck>par%depth_thr) then s%D50top = s%D50(2) s%D90top = s%D90(2) else s%D50top = s%D50(1) s%D90top = s%D90(1) endif enddo enddo endif endif ! if par%t>par%morstart end subroutine bed_update_conceptmodel !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine mccall_vanrijn_1(s,par) use params use spaceparams use xmpi_module use math_tools use paramsconst use morphevolution IMPLICIT NONE type(spacepars),target :: s type(parameters) :: par real*8 :: omegap real*8,save :: phirad,delta,khlim,reposerad,reposedzdx integer :: i,j,jg real*8,dimension(:,:),allocatable,save :: dudtsmooth real*8,dimension(:,:),allocatable,save :: fsed real*8,dimension(:,:),allocatable,save :: Arms,umeanupd,uvarupd real*8,dimension(:,:),allocatable,save :: umeanupdphi,uvarupdphi real*8,dimension(:,:),allocatable,save :: shields,qsedu,qsedutemp,dist real*8,dimension(:,:),allocatable,save :: blphi,facbl,facrw,facslp,facrwf real*8,dimension(:,:),allocatable,save :: ulocal,ulocalold,philocal real*8,dimension(:,:),allocatable,save :: dcfinl,dcfl,cffac real*8,dimension(:) ,allocatable,save :: shieldscrit,dstar real*8,dimension(:,:),allocatable,save :: signShields,thetacrlocal real*8,dimension(:,:),allocatable,save :: phishields,nEF real*8,dimension(:,:),allocatable,save :: dzbdxf real*8 :: Te,kvis,Sster,cc1,cc2,wster real*8 , dimension(:),allocatable,save :: w if (.not. allocated(dudtsmooth)) then allocate(dudtsmooth(s%nx+1,s%ny+1)) allocate(fsed(s%nx+1,s%ny+1)) allocate(shields(s%nx+1,s%ny+1)) allocate(qsedu(s%nx+1,s%ny+1)) allocate(qsedutemp(s%nx+1,s%ny+1)) allocate(dist(s%nx+1,s%nx+1)) allocate (blphi (s%nx+1,s%ny+1)) allocate (facbl (s%nx+1,s%ny+1)) allocate (facrw (s%nx+1,s%ny+1)) allocate (facrwf(s%nx+1,s%ny+1)) allocate (facslp (s%nx+1,s%ny+1)) allocate (ulocal (s%nx+1,s%ny+1)) allocate (ulocalold (s%nx+1,s%ny+1)) allocate(philocal(s%nx+1,s%ny+1)) allocate(umeanupdphi(s%nx+1,s%ny+1)) allocate(uvarupdphi(s%nx+1,s%ny+1)) allocate(dcfinl(s%nx+1,s%ny+1)) allocate(dcfl(s%nx+1,s%ny+1)) allocate(cffac(s%nx+1,s%ny+1)) cffac = 0.d0 allocate(Arms(s%nx+1,s%ny+1)) allocate(umeanupd(s%nx+1,s%ny+1)) allocate(uvarupd(s%nx+1,s%ny+1)) umeanupd = 0.d0 uvarupd = 0.d0 allocate(signShields(s%nx+1,s%ny+1)) signShields = 0.d0 allocate(thetacrlocal(s%nx+1,s%ny+1)) thetacrlocal = par%thetcr allocate(shieldscrit(par%ngd)) allocate(dstar(par%ngd)) allocate(dzbdxf(s%nx+1,s%ny+1)) if (par%formG==FORM_WILCOCK_CROW) then allocate(phishields(s%nx+1,s%ny+1)) elseif (par%formG==FORM_ENGELUND_FREDSOE) then allocate(nEF(s%nx+1,s%ny+1)) endif if(par%bulk==1) then Te = 20.d0 kvis = 4.d0/(20.d0+Te)*1d-5 ! Van rijn, 1993 allocate(w(par%ngd)) do jg=1,par%ngd Sster = s%D50(jg)/(4*kvis)*sqrt((par%rhos/par%rho-1)*par%g*s%D50(jg)) cc1 = 1.06d0*tanh(0.064d0*Sster*exp(-7.5d0/Sster**2)) cc2 = 0.22d0*tanh(2.34d0*Sster**(-1.18d0)*exp(-0.0064d0*Sster**2)) wster = cc1+cc2*Sster w(jg) = wster*sqrt((par%rhos/par%rho-1.d0)*par%g*s%D50(jg)) enddo endif fsed = 0.d0 phirad = par%phit/180*par%px delta = (par%rhos-par%rho)/par%rho ulocalold = s%uu ulocal = s%uu philocal = phirad umeanupdphi = 0.d0 uvarupdphi = 0.d0 ! ! find water depth for which kh(based on Trep) ~ 0.1 khlim = 0.01d0*par%Trep**2*par%g/4/par%px**2 ! ! Angle of repose reposerad = par%reposeangle*par%px/180 reposedzdx = tan(reposerad) do i=1,s%nx+1 dist(:,i) = abs(s%xu(i,1)-s%xu(:,1)) dist(:,i) = dist(:,i)/0.25d0 dist(:,i) = max(dist(:,i),1.d0) dist(:,i) = 1.d0-dist(:,i) enddo jg=2 ! Make jg=2, because all trnasports are done in j=2. (differentitation is done with different fluxes) dstar(jg) = (par%g*(par%rhos/par%rho-1)/1.d-6/1.d-6)**(1.d0/3)*par%D50(jg) !Huub use par%D50(2), because this is the second D50 (gravel) if(par%thetcr<0.d0) then shieldscrit(jg) = 0.3d0/(1+1.2d0*dstar(jg))+0.055d0*(1-exp(-0.02d0*dstar(jg))) else shieldscrit(jg) = par%thetcr endif endif !do j=1,s%ny+1 ! do i=1,nx ! Lloc = par%Trep*sqrt(par%g*max(0.1d0,hh(i,j)))/8 ! irange1 = max(1, i-nint(Lloc/dsu(i,j)))! staggered grid means lower should be at least zb(i) ! irange2 = min(s%nx+1,i+nint(Lloc/dsu(i,j))) ! irange2 = max(i+1,irange2) ! staggered grid means upper should be at least zb(i+1) ! dzbdxf(i,j) = (zb(irange2,j)-zb(irange1,j))/(sdist(irange2,j)-sdist(irange1,j)) ! enddo !enddo !dzbdxf(s%nx+1,:) = dzbdxf(nx,:) dzbdxf = s%dzbdx jg=2 ! Replace instead of do loop (2nd fraction is gravel) (Huub) ! compute shields parameter if (par%gwflow==1 .and. par%inclrelweight==1) then !facrw = 0.5d0*max(infil/Kz,-1.d0) facrw = 0.5d0*s%infil/s%Kzinf else facrw = 0.d0 endif where(s%wetu==1) ! try adding (groundwater) head gradient shields = (abs(s%taubx-par%incldzdx*s%dzsdx*par%rho*par%g*s%D50(jg))) /(par%rho*par%g*s%D50(jg)*max(delta+facrw,1e-6)) signShields = sign(1.d0,s%taubx-par%incldzdx*s%dzsdx*par%rho*par%g*s%D50(jg)) elsewhere shields = 0.d0 signShields = 0.d0 endwhere ! ! compute bed slope effect if (par%slopecorr==SLOPECORR_NIELSEN) then where(s%wetu==1) where (s%dzbdx>0.d0) shields = max(shields/cos(min(atan(dzbdxf),reposerad))-min(dzbdxf,reposedzdx)*signShields,0.d0) elsewhere(s%dzbdx<0.d0) shields = max(shields/cos(max(atan(dzbdxf),-reposerad))-max(dzbdxf,-reposedzdx)*signShields,0.d0) endwhere endwhere elseif (par%slopecorr==SLOPECORR_FREDSOE_DEIGAARD) then where(s%wetu==1) where (s%dzbdx>0.d0) shields = shields*cos(atan(min(dzbdxf,reposerad)))* & max(1.d0-signShields*min(dzbdxf/reposedzdx,1.d0),0.d0) elsewhere(s%dzbdx<0.d0) shields = shields*cos(atan(max(dzbdxf,-reposerad)))* & max(1.d0-signShields*max(dzbdxf/reposedzdx,-1.d0),0.d0) endwhere endwhere endif !! compute streaming effect !if(par%streaming==1) then ! ! Note assumes wave propagation is in direction of positive s. Modify to negative ! ! contribution in case wave direction is in negative s direction ! Arms = sqrt(2.d0)/omegap*sqrt(uvarupd) ! where(s%wetu==1) ! where(uvarupd>1d0*abs(umeanupd)) ! fe = exp(5.5d0*(min(170*sqrt(max(shields-0.05d0,0.d0))*s%D50(jg),s%hh)/Arms)**0.2-6.3d0) ! elsewhere(uvarupd<0.25d0*abs(umeanupd)) ! fe = 0.d0 ! elsewhere ! fe = ((uvarupd/abs(umeanupd))-0.25d0)/0.75d0 * & ! exp(5.5d0*(min(170*sqrt(max(shields-0.05d0,0.d0))*s%D50(jg),s%hh)/Arms)**0.2-6.3d0) ! endwhere ! ! shields = shields + (2d0/3/par%px*fe*Arms**3*omegap**3/(par%g*s%hh)) / & ! ((delta+facrw)*par%g*s%D50(jg)) * & ! signShields ! shields = max(shields,0.d0) ! endwhere !endif thetacrlocal = shieldscrit(jg) ! compute sediment transport select case (par%formG) case (FORM_MPM) par%Ctrans = 8.d0 where(s%wetu==1 .and. shields>thetacrlocal) qsedu = (par%Ctrans*(shields-thetacrlocal)**1.5d0)* & sqrt(delta*par%g*s%D50(jg)**3)*signShields elsewhere qsedu = 0.d0 endwhere case (FORM_WONG_PARKER) par%Ctrans = 3.97d0 where(s%wetu==1 .and. shields>thetacrlocal) qsedu = (par%Ctrans*(shields-thetacrlocal)**1.5d0)* & sqrt(delta*par%g*s%D50(jg)**3)*signShields elsewhere qsedu = 0.d0 endwhere case (FORM_FL_VB) ! Fernandez Luque, R., and van Beek, R.(1976). Erosion and Transport of Bed-Load ! Sediment. J. Hydraul, Res., 14(2), 127-144 par%Ctrans = 5.7d0 where(s%wetu==1 .and. shields>thetacrlocal) qsedu = (par%Ctrans*(shields-thetacrlocal)**1.5d0)* & sqrt(delta*par%g*s%D50(jg)**3)*signShields elsewhere qsedu = 0.d0 endwhere case (FORM_WILCOCK_CROW) where(s%wetu==1) phishields = shields/thetacrlocal where(phishields>=1.35d0) qsedu = 14*(1.d0-0.894d0/sqrt(phishields))**4.5d0*shields**1.5*sqrt(delta*par%g)*s%D50(jg)**1.5*signShields elsewhere qsedu = 0.002d0*phishields**7.5d0*shields**1.5*sqrt(delta*par%g)*s%D50(jg)**1.5*signShields endwhere elsewhere qsedu = 0.d0 endwhere case (FORM_ENGELUND_FREDSOE) where(s%wetu==1 .and. shields>thetacrlocal) ! note: pi/6 ~ 0.5236, dynamic friction ~ 0.5-0.65 (Fredsoe & Deigaard) nEF = (1.d0+(0.5236d0*0.5d0/(shields-thetacrlocal))**4)**-0.25d0 qsedu = 5*nEF*(sqrt(shields)-0.7d0*sqrt(thetacrlocal))*sqrt(delta*par%g*s%D50(jg)**3)*signShields elsewhere qsedu = 0.d0 endwhere case (FORM_FREDSOE_DEIGAARD) ! Equation 7.59 where(s%wetu==1 .and. shields>thetacrlocal) ! note: 30/pi = 9.5493, mu = 0.65 -> 30/pi/mu = 14.6912 qsedu = 14.6912d0*(shields-thetacrlocal)*(sqrt(shields)-0.7d0*sqrt(thetacrlocal)) * & sqrt(delta*par%g*s%D50(jg)**3)*signShields elsewhere qsedu = 0.d0 endwhere case(FORM_MCCALL_VANRIJN) ! Eq 10. note: devision by rhos to get m2/s transport rate par%Ctrans = 0.5d0 where(s%wetu==1 .and. shields>thetacrlocal) qsedu = par%Ctrans*s%D50(jg)*dstar(jg)**(-0.3d0)*sqrt(abs(s%taubx)/par%rho)* & ((shields-thetacrlocal)/thetacrlocal)*signShields elsewhere qsedu = 0.d0 endwhere end select where(qsedu>0.d0) qsedu = qsedu * par%uprushfac elsewhere qsedu = qsedu * par%backwashfac endwhere !qsedu = min(qsedu,ustar*100*D50(jg)*(1.d0-par%por)) ! minimum porosity ~0.60 for fluidized bed where(qsedu>0.d0) qsedu = min(qsedu,(1.d0-0.60d0)*(abs(s%qx))) elsewhere(qsedu<0.d0) qsedu = max(qsedu,(1.d0-0.60d0)*(-abs(s%qx))) endwhere if (par%bulk==0) then if(par%thetanum<1.d0) then do i=2,s%nx s%Subg_g(i,:,jg) = (1.d0-par%thetanum)/2*(qsedu(i-1,:)+qsedu(i+1,:)) + & !Huub: bed transport gravel par%thetanum * qsedu(i,:) enddo if(xmpi_istop) s%Subg(1,:,jg) = qsedu(i,:) if(xmpi_isbot) s%Subg(s%nx+1,:,jg) = qsedu(s%nx+1,:) else s%Subg_g(:,:,jg) = qsedu endif ! To Do: fix this properly for MPI do j=1,s%ny+1 do i=1,s%nx+1 if(s%Subg_g(i,j,jg)>0.d0) then s%Subg_g(i,j,jg) = s%Subg_g(i,j,jg)*s%sedcal(jg) !*s%pbbed(i,j,1,jg) else s%Subg_g(i,j,jg) = s%Subg_g(i,j,jg)*s%sedcal(jg) !*s%pbbed(min(i+1,s%nx+1),j,1,jg) endif enddo enddo else s%Tsg(:,:,jg) = max(par%tsfac*s%hu/w(jg),par%Tsmin) s%ceqsg(:,:,jg) = qsedu/s%hu/s%uu s%ceqbg(:,:,jg) = 0.d0 endif end subroutine mccall_vanrijn_1 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine rijper_vanrijn(s,par) use params use spaceparams use xmpi_module use math_tools use interp use paramsconst use morphevolution implicit none type(spacepars),target :: s type(parameters) :: par integer :: i,isig integer :: j,jg real*8,save :: delta, psi_cr0, phi_r real*8 , dimension(:,:),allocatable ,save :: vmg,Cd,Asb,dhdx,dhdy,Ts,hfac,hloc real*8 , dimension(:,:),allocatable ,save :: urms2,Ucr,Ucrc,Ucrw,term1,B2,vero,Ucrb,Ucrs real*8,dimension(:) ,allocatable,save :: shieldscrit,dstar real*8,dimension(:,:),allocatable,save :: signShields,thetacrlocal,gwu_sr,signgwu,signpsi real*8 , dimension(:,:),allocatable ,save :: shields,qsr,phi,psi_cr,psi if (.not. allocated(vmg)) then allocate (vmg (s%nx+1,s%ny+1)) allocate (term1 (s%nx+1,s%ny+1)) allocate (B2 (s%nx+1,s%ny+1)) allocate (Cd (s%nx+1,s%ny+1)) allocate (Asb (s%nx+1,s%ny+1)) allocate (Ucr (s%nx+1,s%ny+1)) allocate (Ucrc (s%nx+1,s%ny+1)) allocate (Ucrw (s%nx+1,s%ny+1)) allocate (urms2 (s%nx+1,s%ny+1)) allocate (hloc (s%nx+1,s%ny+1)) allocate (Ts (s%nx+1,s%ny+1)) allocate (shieldscrit (par%ngd)) allocate (dstar (par%ngd)) allocate (shields(s%nx+1,s%ny+1)) allocate (qsr (s%nx+1,s%ny+1)) allocate (gwu_sr (s%nx+1,s%ny+1)) allocate(signShields(s%nx+1,s%ny+1)) signShields = 0.d0 allocate(thetacrlocal(s%nx+1,s%ny+1)) thetacrlocal = par%thetcr allocate(signgwu(s%nx+1,s%ny+1)) signgwu = 0.d0 allocate(signpsi(s%nx+1,s%ny+1)) signpsi = 0.d0 allocate(psi(s%nx+1,s%ny+1)) allocate(psi_cr(s%nx+1,s%ny+1)) allocate(phi(s%nx+1,s%ny+1)) delta = (par%rhos-par%rho)/par%rho psi_cr0 = par%psi_cr0 phi_r = par%reposeangle*par%px/180 jg=1 ! Replace instead of do loop (2nd fraction is gravel) (Huub) dstar(jg) = (par%g*(par%rhos/par%rho-1)/1.d-6/1.d-6)**(1.d0/3)*par%D50(jg) if(par%thetcr<0.d0) then shieldscrit(jg) = 0.3d0/(1+1.2d0*dstar(jg))+0.055d0*(1-exp(-0.02d0*dstar(jg))) else shieldscrit(jg) = par%thetcr endif endif jg=1 ! Replace instead of do loop (2nd fraction is gravel) (Huub) !First limit the gwu maximum speed (can get very high) where (s%gwu>par%gwu_max) gwu_sr=par%gwu_max signgwu=sign(1.d0,gwu_sr) elsewhere (-s%gwu>par%gwu_max) gwu_sr=-par%gwu_max elsewhere gwu_sr=s%gwu signgwu=sign(1.d0,gwu_sr) endwhere where (s%z_g>0 .and. s%z_s/=s%gwlevel) ! where (s%wetu==1) Huub: correct? 12-01 --> changed to s%z_g>0 s%taubx_gwf=par%cf_ss*par%rho*(gwu_sr/par%por_fr2)**2*signgwu !100= s%cfu elsewhere s%taubx_gwf = 0.d0 endwhere !do j=1,s%ny+1 ! do i=1,s%nx ! cfu_gwf(i,j)= (par%g*((s%gwheight(i,j)+s%gwheight(i+1,j))/2)/(par%kz**2))*(s%gwheight(i,j)-s%gwheight(i+1))*s%dsu(i,j)) !enddo !enddo where(s%z_g>0) ! HUUB, misschien werkt dit niet! ! try adding (groundwater) head gradient shields = (abs(s%taubx_gwf)) /(par%rho*par%g*s%D50(jg)*max(delta,1e-6)) signShields = sign(1.d0,s%taubx_gwf) elsewhere shields = 0.d0 signShields = 0.d0 endwhere !! MPM !par%Ctrans = 8.d0 !where(s%wetu==1 .and. shields>thetacrlocal) ! qsedu = (par%Ctrans*(shields-thetacrlocal)**1.5d0)* & ! sqrt(delta*par%g*s%D50(jg)**3)*signShields !elsewhere ! qsedu = 0.d0 !endwhere thetacrlocal = shieldscrit(jg) ! McCall_Van Rijn par%Ctrans = 0.5d0 where(s%z_g>0 .and. shields>thetacrlocal) qsr = par%qsr_cal*par%Ctrans*s%D50(jg)*dstar(jg)**(-0.3d0)*sqrt(abs(s%taubx_gwf)/par%rho)* & ((shields-thetacrlocal)/thetacrlocal)*signShields elsewhere qsr = 0.d0 endwhere if (par%jacobsen==1) then do j=1,s%ny+1 do i=2,s%nx phi(i,j)= atan(((s%zb(i,j)-s%zb(i+1,1))/s%dsu(i,j))) !phi should be in be radians !phi(i,j)= atan(((s%zb(i,j)-s%zb(i+1,1))/s%dsu(i,j))*par%px/180) !phi should be in be radians psi_cr(i,j)=psi_cr0* (sin(phi_r-phi(i,j)))/(sin(phi_r)) enddo enddo where(s%z_g>0) ! HUUB, misschien werkt dit niet! ! try adding (groundwater) head gradient psi = ((gwu_sr)**2)/(par%g*s%D50(jg)*max(delta,1e-6)) signpsi = sign(1.d0,s%gwu) elsewhere psi = 0.d0 signpsi = 0.d0 endwhere !Jacbosen where(s%z_g>0 .and. psi>psi_cr) qsr = par%C_1 * (psi - psi_cr)**(1.5d0) * (gwu_sr) /(abs(gwu_sr)) elsewhere qsr = 0.d0 endwhere endif where(qsr>0.d0) qsr = min(qsr,(1.d0-0.60d0)*(abs(s%gwqx))) elsewhere(qsr<0.d0) qsr = max(qsr,(1.d0-0.60d0)*(-abs(s%gwqx))) endwhere if(par%thetanum<1.d0) then do i=2,s%nx s%Suss(i,:,jg) = (1.d0-par%thetanum)/2*(qsr(i-1,:)+qsr(i+1,:)) + & par%thetanum * qsr(i,:) enddo if(xmpi_istop) s%Suss(1,:,jg) = qsr(i,:) if(xmpi_isbot) s%Suss(s%nx+1,:,jg) = qsr(s%nx+1,:) else s%Suss(:,:,jg) = qsr endif ! To Do: fix this properly for MPI !do j=1,s%ny+1 ! do i=1,s%nx+1 ! if(s%Suss(i,j,jg)>0.d0) then ! s%Suss(i,j,jg) = s%Suss(i,j,jg)*s%sedcal(jg) ! else ! s%Suss(i,j,jg) = s%Suss(i,j,jg)*s%sedcal(jg) ! endif ! enddo !enddo end subroutine end module morphevolution_1