! Copyright 2014 College of William and Mary ! ! Licensed under the Apache License, Version 2.0 (the "License"); ! you may not use this file except in compliance with the License. ! You may obtain a copy of the License at ! ! http://www.apache.org/licenses/LICENSE-2.0 ! ! Unless required by applicable law or agreed to in writing, software ! distributed under the License is distributed on an "AS IS" BASIS, ! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. ! See the License for the specific language governing permissions and ! limitations under the License. !-------------------------------------------------------------------- SUBROUTINE sed2d_main(it) !-------------------------------------------------------------------- ! Main sed2d routine ! ! Authors: Guillaume Dodet (guillaume.dodet@univ-brest.fr) ! Thomas Guerin (thomas.guerin@univ-lr.fr) ! ! History: ! 03/2013 - G.Dodet: - Corrected bugs in node centered FE-method; ! - Modified filtering section; ! - Added space-variable d50 input and Cdsed ! output; ! 04/2013 - G.Dodet: - Added nskip parameter to skip first iterations ! before updating bed level; ! - Split in 2 routines for each numerical method; ! - Added diffusive filter on total transport; ! 06/2013 - G.Dodet: - Added filter on velocity before computing ! sediment fluxes to prevent development of ! oscillations along the coastline induced by ! wave-force/barotropic pressure gradient ! inconsistency (temporary solution); ! - Added diffusive filter on sediment fluxes ! - Added bedform predictors and new roughness; ! - Added timers; ! 07/2013 - G.Dodet: - Converted output transport rate in kg/m/s; ! - G.Dodet: Added iterative shapiro filter for currents ! - T.Guerin: Added Camenen and Larson (2011) ! 04/2014 - T.Guerin: - Removed part for getting offshore wave ! parameters (needed for Camemen and Larson ! formula and wave asymmetry calculation) and ! replaced it by local parameters computation ! - Added initialization of Cd_e, z0_e ! - Removed morphological time step (dtsed2d) to ! be consistent with multi-class mode (will be ! adapted in the future) ! 10/2016 - T.Guerin: Merging single-class and multi-class routines ! (one single routine is also kept for the ! different numerical schemes) ! 11/2016 - T.Guerin: Adding WENO scheme to solve the Exner equation ! 04/2017 - T.Guerin: Added morphological ramp factor ! 05/2017 - T.Guerin: - Added Wu an Lin (2014) formula ! - Added TR2004 formula (van Rijn et al., 2004) !-------------------------------------------------------------------- USE schism_glbl, ONLY : Cdp,dldxy,dp,dt,eta2,idry,idry_e,idry_s,ielg,iegl,& &ipgl,indel,iplg,isdel,isidenode,isidenei2, & &elside,nea,nne,i34,elnode,np,npa,ns,nsa,pi, & &out_wwm,su2,sv2,timer_ns,rhosed,rkind,kbs,zs,nvrt,& &ihydraulics,isblock_sd,errmsg,time_stamp,area, & &in_dir,out_dir,len_in_dir,len_out_dir USE schism_msgp, ONLY : comm,exchange_e2d,exchange_p2d,exchange_s2d, & &ierr,myrank,parallel_abort,rtype,parallel_finalize USE hydraulic_structures, ONLY : nhtblocks USE sed2d_mod, ONLY : Cdsed,cflsed,d_class,diffac,dpdxy,dpdxy_e, & &dtsed2d,h0_sed,idrag,idsed2d,ifilt,imeth, & &imorpho,irough,iskip,islope,itrans,nskip, & &poro,qav,qb,qb_e,qdt_e,qfilter,qramp,qs,qs_e, & &qtot,qtot_e,ufilter,vc_area, & &transfac,z0_e,z0cr_e,z0sw_e,z0wr_e,u2_tmp,v2_tmp,& &bed_delta,d50,d90,F_class,h_inf, & &h_lim_max,h_lim_min,h_top,h2,nb_class, & &nb_layer,morfac,ised_dump,imnp USE sed2d_friction USE sed2d_transport USE sed2d_filter USE sed2d_morpho IMPLICIT NONE INCLUDE 'mpif.h' !- Arguments -------------------------------------------------------- INTEGER, INTENT(IN) :: it !- Local variables -------------------------------------------------- INTEGER :: i,id,iel,inode,iside,j,k,k2,n,neigh,rank_tmp REAL(rkind) :: beta,D_star,D_star_class,d50_e,d90_e,hs,htot, & &kpeak,rampfac,sum_F_class,suru,surv,tau_cr, & &tau_cr_class,theta_cr,theta_cr_class,tmp,tp,Ucr_c, & &Ucr_c_class,Ucr_w,Ucr_w_class,u_star,udir,ue,uorb, & &uorbp,utmp,vtmp,ve,wdir,wdirc,wdirs,wlpeak,wtmp,z0 REAL(rkind), DIMENSION(npa) :: dp_filt,dp_tmp,dp1,neigh2,sum_dh REAL(rkind), DIMENSION(nea) :: Cd_e,tmp_e REAL(rkind), DIMENSION(nsa) :: su2_tmp,su2_tmp1,su2_tmp2,sv2_tmp, & &sv2_tmp1,sv2_tmp2,tmp_s REAL(rkind), DIMENSION(npa,2) :: qdt REAL(rkind), DIMENSION(nb_class) :: F_class_e REAL(rkind), DIMENSION(npa,nb_class) :: dh_class REAL(rkind), DIMENSION(nea,2,nb_class) :: qb_class_e,qs_class_e, & &qtot_class_e !- Dumping-related variables --------------------------------------- INTEGER,save :: ne_dump INTEGER, allocatable :: ie_dump(:) REAL(rkind),save :: t_dump !time in dumping option REAL(rkind), allocatable :: vol_dump(:) logical, save :: first_call=.true. !-------------------------------------------------------------------- !- Compute ramp factor IF(it0) THEN rampfac = TANH(2.d0*(it-iskip)*dt/qramp) ELSE rampfac = 1.d0 ENDIF ENDIF #ifdef INCLUDE_TIMING wtmp = mpi_wtime() !start of timer #endif !- Side velocities filter before computing transport ---------------- !- Final vel used is s[u,v]2_tmp su2_tmp=0 sv2_tmp=0 do i=1,nsa if(idry_s(i)==1) cycle do k=kbs(i),nvrt-1 su2_tmp(i)=su2_tmp(i)+(su2(k,i)+su2(k+1,i))/2*(zs(k+1,i)-zs(k,i)) sv2_tmp(i)=sv2_tmp(i)+(sv2(k,i)+sv2(k+1,i))/2*(zs(k+1,i)-zs(k,i)) enddo !k htot=zs(nvrt,i)-zs(kbs(i),i) if(htot<=0) CALL parallel_abort('SED2D: htot<=0') su2_tmp(i)=su2_tmp(i)/htot sv2_tmp(i)=sv2_tmp(i)/htot enddo !i IF(ufilter /= 0) THEN su2_tmp1 = su2_tmp !su2(1,:) sv2_tmp1 = sv2_tmp !sv2(1,:) su2_tmp2 = su2_tmp !su2(1,:) sv2_tmp2 = sv2_tmp !sv2(1,:) !- Skip if the element has a side at the wet-dry interface DO i = 1,nsa IF(isdel(1,i)>0 .AND. isdel(2,i)>0) THEN !then resident and internal IF(idry_e(isdel(1,i))+idry_e(isdel(2,i))==1) THEN !wet/dry interface su2_tmp1(i) = 0.d0 sv2_tmp1(i) = 0.d0 ELSE su2_tmp1(i) = su2_tmp(i) !su2(1,i) sv2_tmp1(i) = sv2_tmp(i) !sv2(1,i) ENDIF ENDIF ENDDO !nsa CALL exchange_s2d(su2_tmp1) CALL exchange_s2d(sv2_tmp1) !- Iterative shapiro filter applied on the current field !JZ: lon/lat; DO k = 1,ufilter DO i = 1,ns if(isdel(2,i)==0.or.idry_s(i)==1) CYCLE if(ihydraulics/=0.and.nhtblocks>0) then if(isblock_sd(1,i)/=0) cycle endif suru = 0.d0 surv = 0.d0 DO j = 1,4 id = isidenei2(j,i) suru = suru+su2_tmp1(id) surv = surv+sv2_tmp1(id) ENDDO su2_tmp2(i) = su2_tmp1(i)+0.25d0*(suru-4.d0*su2_tmp1(i)) sv2_tmp2(i) = sv2_tmp1(i)+0.25d0*(surv-4.d0*sv2_tmp1(i)) ENDDO !ns DO i = 1,ns if(isdel(2,i)==0.or.idry_s(i)==1) CYCLE if(ihydraulics/=0.and.nhtblocks>0) then if(isblock_sd(1,i)/=0) cycle endif su2_tmp(i) = su2_tmp2(i) sv2_tmp(i) = sv2_tmp2(i) ENDDO CALL exchange_s2d(su2_tmp) CALL exchange_s2d(sv2_tmp) su2_tmp1 = su2_tmp sv2_tmp1 = sv2_tmp ENDDO !k=1,ufilter IF(myrank==0) WRITE(16,*)'done filtering velocities (sed2d)' ENDIF !ufilter/=0 #ifdef INCLUDE_TIMING timer_ns(11) = timer_ns(11)+mpi_wtime()-wtmp !end of timer #endif !- Dumping/dredging if(ised_dump/=0) then !For 1st call (including hot), init. read if(first_call) then !Time stamps in this file must be one of the time steps open(18,file=in_dir(1:len_in_dir)//'sed_dump.in',status='old') read(18,*) do read(18,*,iostat=k)t_dump,ne_dump !time in sec if(k/=0) then ised_dump=0 !reset exit endif if(t_dump>=time_stamp) exit read(18,*) !vol enddo endif !first_call if(ised_dump/=0.and.abs(t_dump-time_stamp)<1.e-4) then !in case end of file allocate(ie_dump(ne_dump),vol_dump(ne_dump),stat=k) if(k/=0) call parallel_abort('SED2D: alloc (9)') read(18,*)(ie_dump(k),vol_dump(k),k=1,ne_dump) if(myrank==0) write(16,*)'SED2D,start dumping at time:',real(t_dump),ne_dump !Modify depth but not bottom() or bed_frac do k=1,ne_dump iel=ie_dump(k) !global index if(iegl(iel)%rank==myrank) then i=iegl(iel)%id !local index tmp=vol_dump(k)/area(i) !m do j=1,i34(i) dp(elnode(j,i))=dp(elnode(j,i))-tmp enddo !j endif !iegl enddo !k deallocate(ie_dump,vol_dump) !Prep for next record read(18,*,iostat=k)t_dump,ne_dump if(k/=0) ised_dump=0 !reset endif !ised_dump/=0 endif !ised_dump/=0 !-------------------------------------------------------------------- !- Compute sediment transport at element center for each grain size !- class in the top layer of the sediment column !-------------------------------------------------------------------- !- Initialization Cd_e(:) = 0.d0 !Drag coefficient (-) dpdxy_e(:,:) = 0.d0 !Bed slope (m/m) z0_e(:) = 0.d0 !Roughness length (m) qb_class_e(:,:,:) = 0.d0 !Bed load transport for each grain size class (m3/s/m) qs_class_e(:,:,:) = 0.d0 !Suspended load transport for each grain size class (m3/s/m) qtot_class_e(:,:,:) = 0.d0 !Total transport for each grain size class (m3/s/m) qb_e(:,:) = 0.d0 !Total bed load transport (sum over grain size classes) (m3/s/m) qs_e(:,:) = 0.d0 !Total suspended load transport (sum over grain size classes) (m3/s/m) qtot_e(:,:) = 0.d0 !Total total transport (sum over grain size classes) (m3/s/m) !- Start loop over elements DO i = 1,nea #ifdef INCLUDE_TIMING wtmp = mpi_wtime() !start of timer #endif !- Initialization d50_e = 0.d0 !d50 (m) d90_e = 0.d0 !d90 (m) htot = 0.d0 !Total depth (m) ue = 0.d0 !X-velocity (m/s) ve = 0.d0 !Y-velocity (m/s) uorb = 0.d0 !Init orbital velocity computed by wwm (m/s) hs = 0.d0 !Significant wave height computed by wwm(m) tp = 0.d0 !Peak period computed by wwm (s) wlpeak = 0.d0 !Wave length associated to peak period by wwm(m) uorbp = 0.d0 !Peak orbital velocity (m/s) wdir = 0.d0 !Mean average energy transport direction (rad) wdirc = 0.d0 !cosine value needed for calculating mean wave direction in each element wdirs = 0.d0 !sine value needed for calculating mean wave direction in each element dp1 = dp !Store depth before bed update for output purpose DO j = 1,i34(i) inode = elnode(j,i) htot = htot+(eta2(inode)+dp(inode))/i34(i) !Total water depth ENDDO IF(idry_e(i)==1.OR.htot<=h0_sed) CYCLE DO j = 1,i34(i) inode = elnode(j,i) iside = elside(j,i) !- Velocity ue = ue+su2_tmp(iside)/i34(i) ve = ve+sv2_tmp(iside)/i34(i) !- d50 and d90 corresponding to the top layer d50_e = d50_e+d50(inode,1)/i34(i) d90_e = d90_e+d90(inode,1)/i34(i) !- bed slope dpdxy_e(i,1) = dpdxy_e(i,1)+dp(inode)*dldxy(j,1,i) dpdxy_e(i,2) = dpdxy_e(i,2)+dp(inode)*dldxy(j,2,i) !- Wave parameters #ifdef USE_WWM uorb = uorb+out_wwm(inode,22)/i34(i) hs = hs+out_wwm(inode,1)/i34(i) tp = tp+out_wwm(inode,12)/i34(i) wlpeak = wlpeak+out_wwm(inode,17)/i34(i) wdirc = wdirc+(COS((270.d0-out_wwm(inode,9))*pi/180.d0))/i34(i) wdirs = wdirs+(SIN((270.d0-out_wwm(inode,9))*pi/180.d0))/i34(i) #endif !- Drag coefficient IF(idrag==1)THEN Cd_e(i) = Cd_e(i)+Cdp(inode)/i34(i) ENDIF ENDDO! j = 1,i34 #ifdef USE_WWM wdir = ATAN2(wdirs,wdirc) #endif #ifdef INCLUDE_TIMING timer_ns(4) = timer_ns(4)+mpi_wtime()-wtmp !end of timer wtmp = mpi_wtime() !start of timer #endif !- Compute bed slope in streamwise direction used for slope effect on ! total transport (Soulsby, 1997). Positive if flows runs uphill. IF(islope == 1) THEN udir = ATAN2(ve,ue) beta = -(dpdxy_e(i,1)*COS(udir)+dpdxy_e(i,2)*SIN(udir)) beta = MIN(beta,0.6d0) !May lead to negative transport otherwise ELSE beta = 0.d0 ENDIF !- Compute threshold parameters, roughness and drag coefficients ! related to the global sediment characteristics (i.e. d50 and d90) CALL compute_thresh(d50_e,d90_e,htot,tp,D_star,Ucr_c,Ucr_w, & &theta_cr,tau_cr) SELECT CASE(irough) CASE(1) !Skin friction only z0_e(i) = d50_e/12.d0 CASE(2) !Bedform associated roughness (Soulsby, 1997) CALL bedform_predictor_s97(d50_e,htot,ue,ve,uorb,tp, & &tau_cr,theta_cr,z0cr_e(i),z0wr_e(i),z0sw_e(i),z0_e(i)) CASE(3) !Bedform associated roughness (Van-Rijn, 2007) CALL bedform_predictor_vr07(d50_e,htot,ue,ve,uorb, & &z0cr_e(i),z0wr_e(i),z0sw_e(i),z0_e(i)) END SELECT IF(IABS(idrag).GT.1) CALL compute_drag(d50_e,htot,z0_e(i),ue,ve,idrag,Cd_e(i)) #ifdef INCLUDE_TIMING timer_ns(5) = timer_ns(5)+mpi_wtime()-wtmp !end of timer wtmp = mpi_wtime() !start of timer #endif #ifdef USE_DEBUG !- Test if parameters are physically correct if(htot<0.d0 .or. SQRT(ue*ue+ve*ve)>7.5d0 .or. d50_e>0.05d0.or. & &d90_e>0.10d0 .or. SQRT(dpdxy_e(i,1)**2.d0+dpdxy_e(i,2)**2.d0)& &>5.d0 .or. uorb>10.d0 .or. hs>15.d0 .or. tp>30.d0 .or. & &wlpeak>1500.d0 .or. Cd_e(i)>10.d0 .or. z0>0.5d0) then write(12,*)'Warning (1) for it:',it,' and el: ', ielg(i) write(12,*)'Check following values: ' write(12,*)'total depth = ',htot write(12,*)'current velocity = ',SQRT(ue*ue+ve*ve) write(12,*)'d50 = ',d50_e write(12,*)'d90 = ',d90_e write(12,*)'slope = ',SQRT(dpdxy_e(i,1)**2.d0+dpdxy_e(i,2)**2.d0) write(12,*)'orbital velocity = ',uorb write(12,*)'hs = ',hs write(12,*)'tp = ',tp write(12,*)'wlpeak = ',wlpeak write(12,*)'Cd = ',Cd_e(i) write(12,*)'z0 = ',z0_e(i) ! CALL parallel_abort('Sed2d: chech value in outputs/nonfatal') endif #endif !- Start loop over grain size classes DO k=1,nb_class !- Compute sediment class fractions at element center if (k==1) then F_class_e(:) = 0.d0 do j=1,i34(i) inode = elnode(j,i) do k2 = 1,nb_class F_class_e(k2) = F_class_e(k2)+F_class(inode,k2,1)/i34(i) enddo enddo !j=1,i34(i) endif !- Compute threshold parameters related to grain size class k CALL compute_thresh(d_class(k),d90_e,htot,tp,D_star_class, & &Ucr_c_class,Ucr_w_class,theta_cr_class,tau_cr_class) !- Compute sediment transport for grain size class k ! NB: we always keep Cd computed from the d50 SELECT CASE(itrans) CASE(1) !Engelund-Hansen (1967) CALL eh67(Cd_e(i),d_class(k),ue,ve,beta,qtot_class_e(i,:,k)) CASE(2) !Ackers and White (1973) CALL aw73(Cd_e(i),d_class(k),ue,ve,htot,beta,D_star_class, & &qtot_class_e(i,:,k)) CASE(3) !Soulsby - Van Rijn (Soulsby, 1997) CALL svr97_bedl(Cd_e(i),d_class(k),ue,ve,htot,dpdxy_e(i,:),& &Ucr_c_class,tau_cr_class,uorb,qb_class_e(i,:,k)) CALL svr97_susp(Cd_e(i),d_class(k),ue,ve,D_star_class, & &Ucr_c_class,uorb,qs_class_e(i,:,k)) CASE(4) !Van-Rijn (2007) CALL vr07_bedl(Cd_e(i),d_class(k),ue,ve,htot,dpdxy_e(i,:),uorb,& &Ucr_c_class,Ucr_w_class,tau_cr_class,qb_class_e(i,:,k)) CALL vr07_susp(d_class(k),ue,ve,uorb,D_star_class,Ucr_c_class, & &Ucr_w_class,qs_class_e(i,:,k)) CASE(5) !Camenen and Larson (2011) CALL cl11_tot(d_class(k),D_star_class,dpdxy_e(i,:),hs,htot,& &theta_cr_class,tp,ue,ve,wdir,wlpeak, & &qb_class_e(i,:,k),qs_class_e(i,:,k)) CASE(6) !Wu and Lin (2014) CALL wl14_tot(d_class(k),d_class(:),d50_e,d90_e,D_star_class, & &F_class_e(:),htot,dpdxy_e(i,:),ue,ve,hs,tp,wdir,wlpeak,& &qb_class_e(i,:,k),qs_class_e(i,:,k)) CASE(7) !TRANSPOR2004 (van Rijn et al., 2004) CALL tr04_tot(d_class(k),d50_e,d90_e,htot,dpdxy_e(i,:),ue,ve,hs,& &tp,wdir,wlpeak,qb_class_e(i,:,k),qs_class_e(i,:,k)) END SELECT !- Apply slope effect and compute total transport for grain size class k if (itrans==3.or.itrans==4.or.itrans==5.or.itrans==6) then qb_class_e(i,:,k) = qb_class_e(i,:,k)*(1.d0-1.6d0*beta) qs_class_e(i,:,k) = qs_class_e(i,:,k)*(1.d0-1.6d0*beta) qtot_class_e(i,:,k) = qb_class_e(i,:,k)+qs_class_e(i,:,k) elseif (itrans==7) then !NB: no slope effect with tr04 formula because already done in tr04_tot qtot_class_e(i,:,k) = qb_class_e(i,:,k)+qs_class_e(i,:,k) endif !- Multiply the transport by the corresponding grain size fraction qb_class_e(i,:,k) = F_class_e(k)*qb_class_e(i,:,k) qs_class_e(i,:,k) = F_class_e(k)*qs_class_e(i,:,k) qtot_class_e(i,:,k) = F_class_e(k)*qtot_class_e(i,:,k) !- Apply ramp, correction, and diffusion factors if required qtot_class_e(i,1,k) = rampfac*transfac*(qtot_class_e(i,1,k)+(1.d0-poro)*diffac* & &abs(qtot_class_e(i,1,k))*dpdxy_e(i,1)) qtot_class_e(i,2,k) = rampfac*transfac*(qtot_class_e(i,2,k)+(1.d0-poro)*diffac* & &abs(qtot_class_e(i,2,k))*dpdxy_e(i,2)) !- Sum over grain size classes qb_e(i,:) = qb_e(i,:) + qb_class_e(i,:,k) qs_e(i,:) = qs_e(i,:) + qs_class_e(i,:,k) qtot_e(i,:) = qtot_e(i,:) + qtot_class_e(i,:,k) #ifdef INCLUDE_TIMING timer_ns(6) = timer_ns(6)+mpi_wtime()-wtmp !end of timer wtmp = mpi_wtime() !start of timer #endif ENDDO !k=1,nb_class ENDDO !i=1,nea #ifdef INCLUDE_TIMING timer_ns(7) = timer_ns(7)+mpi_wtime()-wtmp !end of timer wtmp = mpi_wtime() !start of timer #endif !-------------------------------------------------------------------- !- Compute bed change at nodes for each grain size class !-------------------------------------------------------------------- DO k=1,nb_class !- Apply diffusive filter on total transport if required IF(qfilter == 1) THEN CALL sed2d_filter_diffu(qtot_class_e(:,1,k),tmp_e,nea) qtot_class_e(:,1,k) = tmp_e CALL sed2d_filter_diffu(qtot_class_e(:,2,k),tmp_e,nea) qtot_class_e(:,2,k) = tmp_e CALL exchange_e2d(qtot_class_e(:,1,k)) CALL exchange_e2d(qtot_class_e(:,2,k)) IF(myrank==0) WRITE(16,*)'done filtering transport rate (sed2d) for grain size class:',k !Recompute total transport (sum over grain size classes) if(k==1) qtot_e(:,:) = 0.d0 qtot_e(:,1) = qtot_e(:,1) + qtot_class_e(:,1,k) qtot_e(:,2) = qtot_e(:,2) + qtot_class_e(:,2,k) ENDIF !- Time-integrated sediment transport (m^2) for grain size class k ! NB: qdt_e is related to the transport of grain size class k ! but not total transport over all grain size classes qdt_e(:,:) = 0.d0 !initialization qdt_e(:,1) = qtot_class_e(:,1,k)*dt qdt_e(:,2) = qtot_class_e(:,2,k)*dt !- Apply morphological factor qdt_e(:,1) = qdt_e(:,1)*morfac qdt_e(:,2) = qdt_e(:,2)*morfac IF (imorpho==1 .AND. it>iskip) THEN !Compute bottom evolution for grain size class k if (imeth==1) then CALL base_scheme(it) elseif (imeth==2) then CALL weno_scheme(it) endif !Save bottom evolution for grain size class k dh_class(:,k) = bed_delta(:) ENDIF ENDDO !k=1,nb_class !-------------------------------------------------------------------- !- Update bathymetry - !-------------------------------------------------------------------- IF(imorpho == 1 .AND. it>iskip) THEN sum_dh(:) = 0.d0 DO i=1,npa DO k=1,nb_class sum_dh(i) = sum_dh(i) + dh_class(i,k) ENDDO sum_dh(i) = sum_dh(i)*imnp(i) !apply morphological ramp factor dp(i) = dp(i) + sum_dh(i) !update bathy ENDDO ENDIF !-------------------------------------------------------------------- !- Update grain size fractions - !-------------------------------------------------------------------- IF (nb_class>1.and.imorpho==1.and.it>iskip) THEN DO i=1,npa IF (sum(dh_class(i,:))==0.d0) CYCLE !- Compute new sediment class fractions DO k=1,nb_class IF (sum_dh(i)>0.d0 .or. & &(sum_dh(i)==0.d0.and.dh_class(i,k)>0.d0)) THEN !EROSION case !YJZ: this check may be too restrictive? IF(sum_dh(i)>=h_lim_min) THEN write(errmsg,*)'negative thickness for layer #2: & &need to increase h_lim_min, or decrease time step, or decrease dtsed2d:',sum_dh(i),h_lim_min CALL parallel_abort(errmsg) ENDIF !New fractions for surface layer; draw from layer #2 so top layer !thickness stays same F_class(i,k,1)=(h_top*F_class(i,k,1)+sum_dh(i)*F_class(i,k,2)-dh_class(i,k))/h_top IF(F_class(i,k,1)<0.d0) F_class(i,k,1)=0.d0 ELSEIF (sum_dh(i)<0.d0 .or. & &(sum_dh(i)==0.d0.and.dh_class(i,k)<0.d0)) THEN !DEPOSITION case !YJZ: this check may be too restrictive? IF(ABS(sum_dh(i))>=h_lim_max) THEN write(errmsg,*)'layer #2 thickness too large: & &need to increase h_lim_max, or decrease time step, or decrease dtsed2d:',sum_dh(i),h_lim_max CALL parallel_abort(errmsg) ENDIF !New fractions for surface layer !Deposit into layer #1 (-dh), and swap out same amount to layer #2 F_class(i,k,1)=((h_top+sum_dh(i))*F_class(i,k,1)-dh_class(i,k))/h_top IF(F_class(i,k,1)<0.d0) F_class(i,k,1)=0.d0 !New sub-surface layer fractions F_class(i,k,2)=(h2(i)*F_class(i,k,2)-sum_dh(i)*F_class(i,k,1))/(h2(i)-sum_dh(i)) !>=0 ENDIF !sum_dh(i) sign ENDDO !k=1,nb_class !- Update layer #2 IF (sum_dh(i)>0.d0) THEN !EROSION case !New sub-surface layer thickness h2(i) = h2(i) - sum_dh(i) !should >0 b/cos sum_dh(i)<=h_lim_min IF (h2(i)0 h2(i) = h2(i) + h_inf !Adding new layer => updating fractions - implies infinite supply !of sediment DO n=3,nb_layer-1 F_class(i,:,n) = F_class(i,:,n+1) ENDDO ENDIF !h2(i) ELSEIF (sum_dh(i)<0.d0) THEN !DEPOSITION !New sub-surface layer thickness h2(i) = h2(i) - sum_dh(i) !sum_dh(i)<0 IF (h2(i)>h_lim_max) THEN !Splitting sub-surface layer into 2 layers !New sub-surface layer thickness h2(i) = h2(i) - h_inf if(h2(i) updating fractions DO n=1,nb_layer-2 F_class(i,:,nb_layer-n+1) = F_class(i,:,nb_layer-n) ENDDO ENDIF !h2(i)>h_lim_max ENDIF !sum_dh(i) !- Compute sum of fractions and correct fractions if sum /= 1 DO n=1,nb_layer sum_F_class = 0.d0 DO k=1,nb_class sum_F_class = sum_F_class + F_class(i,k,n) ENDDO IF (sum_F_class/=1.d0) THEN if(sum_F_class==0) then if (n==1 .and. imnp(i)==0) then !non-erodable surface layer case F_class(i,:,n) = 0.d0 F_class(i,nb_class,n) = 1.d0 else write(errmsg,*)'SED2D: all eroded; ',n,iplg(i) call parallel_abort(errmsg) endif else F_class(i,:,n) = F_class(i,:,n)/sum_F_class endif ENDIF ENDDO !nb_layer !- Compute new d50 (weighted geometric mean) DO n=1,nb_layer d50(i,n) = d_class(1)**F_class(i,1,n) DO k=2,nb_class d50(i,n) = d50(i,n) * d_class(k)**F_class(i,k,n) ENDDO ENDDO !nb_layer !- Compute new d90 (to be improved) d90(i,:) = 2.5d0*d50(i,:) ENDDO !i=1,npa ENDIF !imorpho=1 .and. it>iskip #ifdef INCLUDE_TIMING timer_ns(9) = timer_ns(9)+mpi_wtime()-wtmp !end of timer wtmp = mpi_wtime() !start of timer #endif !-------------------------------------------------------------------- !- Apply user-defined filter - !-------------------------------------------------------------------- IF (imorpho==1.AND.it>iskip.AND.ifilt>0) THEN SELECT CASE(ifilt) CASE(1) CALL sed2d_filter_extrema(dp,dp_filt) CASE(2) CALL sed2d_filter_slope(dp,dp_filt) CASE(3) CALL sed2d_filter_slope(dp,dp_tmp) CALL sed2d_filter_extrema(dp_tmp,dp_filt) END SELECT ! prevent filters from changing depths where imnp=0 do i = 1, npa if (imnp(i) .eq. 1) dp(i) = dp_filt(i) end do CALL exchange_p2d(dp) IF(myrank==0)WRITE(16,*)'done filtering new bathy. (sed2d)' ENDIF !-------------------------------------------------------------------- !- Check if new bathymetry /= NaN - !-------------------------------------------------------------------- DO i = 1,npa IF(dp(i)/=dp(i)) THEN WRITE(errmsg,*)'Sed2d: NaN in depth at step:',it,' and node: ', iplg(i),';dp = ',dp(i) CALL parallel_abort(errmsg) ENDIF ENDDO #ifdef INCLUDE_TIMING timer_ns(10) = timer_ns(10)+mpi_wtime()-wtmp !end of timer wtmp = mpi_wtime() !start of timer #endif !-------------------------------------------------------------------- !- Interpolate variables at node for output purpose and also for ! passing back Cdsed to schism_step !-------------------------------------------------------------------- Cdsed = 0.d0 cflsed = 0.d0 dpdxy = 0.d0 qav = 0.d0 qdt = 0.d0 qb = 0.d0 qs = 0.d0 qtot = 0.d0 DO i = 1,np IF(idry(i) == 1) CYCLE neigh = 0 DO j = 1,nne(i) iel = indel(j,i) Cdsed(i) = Cdsed(i)+Cd_e(iel) dpdxy(i,:) = dpdxy(i,:)+dpdxy_e(iel,:) qav(i,:) = qav(i,:)+qtot_e(iel,:)*dt qdt(i,:) = qdt(i,:)+qtot_e(iel,:)*dt qb(i,:) = qb(i,:)+qb_e(iel,:) qs(i,:) = qs(i,:)+qs_e(iel,:) qtot(i,:) = qtot(i,:)+qtot_e(iel,:) neigh = neigh+1 ENDDO !j if(neigh==0.or.eta2(i)+dp1(i)==0) call parallel_abort('SED2D: div. by 0 (12)') Cdsed(i) = Cdsed(i)/neigh dpdxy(i,1) = dpdxy(i,1)/neigh dpdxy(i,2) = dpdxy(i,2)/neigh qav(i,:) = rhosed*qav(i,:)/(neigh*dt*dtsed2d) !(kg/m/s) qdt(i,:) = qdt(i,:)/neigh !(m3/m) qb(i,:) = rhosed*qb(i,:)/neigh !(kg/m/s) qs(i,:) = rhosed*qs(i,:)/neigh !(kg/m/s) qtot(i,:) = rhosed*qtot(i,:)/neigh !(kg/m/s) cflsed(i) = 3.d0*SQRT(qdt(i,1)*qdt(i,1)+qdt(i,2)*qdt(i,2))/ & &((eta2(i)+dp1(i))*SQRT(vc_area(i))) ENDDO !np CALL exchange_p2d(Cdsed(:)) DO i=1,2 CALL exchange_p2d(dpdxy(:,i)) CALL exchange_p2d(qav(:,i)) CALL exchange_p2d(qdt(:,i)) CALL exchange_p2d(qb(:,i)) CALL exchange_p2d(qs(:,i)) CALL exchange_p2d(qtot(:,i)) ENDDO CALL exchange_p2d(cflsed(:)) u2_tmp = 0.d0 v2_tmp = 0.d0 neigh2 = 0.d0 DO i = 1,nsa DO j = 1,2 inode = isidenode(j,i) IF(idry(inode) == 1) CYCLE u2_tmp(inode) = u2_tmp(inode)+su2_tmp(i) v2_tmp(inode) = v2_tmp(inode)+sv2_tmp(i) neigh2(inode) = neigh2(inode)+1 ENDDO !j=1,2 ENDDO !nsa CALL exchange_p2d(u2_tmp(:)) !check if necessary CALL exchange_p2d(v2_tmp(:)) !check if necessary CALL exchange_p2d(neigh2(:)) !check if necessary DO i = 1,npa IF(idry(i) == 1) CYCLE if(neigh2(i)==0) call parallel_abort('sed2d_main_node: impossible (1)') u2_tmp(i) = u2_tmp(i)/neigh2(i) v2_tmp(i) = v2_tmp(i)/neigh2(i) ENDDO !npa tmp=sum(dp) WRITE(12,*)'SED2D, it=',it,' ; sum of depth=',tmp IF(myrank==0) WRITE(16,*)'done computing total sediment load (sed2d) and & &exiting sed2d_main; check nonfatal* for sum of depths' if(tmp/=tmp) call parallel_abort('sed2d_main_node: NaN in sum') first_call=.false. #ifdef INCLUDE_TIMING timer_ns(8) = timer_ns(8)+mpi_wtime()-wtmp !end of timer #endif END SUBROUTINE sed2d_main