From 6c03672fb8bcb4f891995e6d404c8c3982142f68 Mon Sep 17 00:00:00 2001 From: Zhiyun Du Date: Wed, 18 Mar 2026 15:33:13 -0400 Subject: [PATCH 1/2] Add support for 3D sediment bed variable output --- src/Core/schism_glbl.F90 | 5 ++ src/Core/schism_io.F90 | 103 ++++++++++++++++++++++++- src/Core/scribe_io.F90 | 156 ++++++++++++++++++++++++++++++++++++-- src/Hydro/misc_subs.F90 | 51 +++++++++++++ src/Hydro/schism_init.F90 | 61 +++++++++++++-- src/Hydro/schism_step.F90 | 96 +++++++++++++++++++++-- src/Sediment/sed_init.F90 | 13 ++++ src/Sediment/sed_mod.F90 | 3 + src/Sediment/sediment.F90 | 28 +++++++ 9 files changed, 498 insertions(+), 18 deletions(-) diff --git a/src/Core/schism_glbl.F90 b/src/Core/schism_glbl.F90 index b36bf40e1..d26a9a606 100644 --- a/src/Core/schism_glbl.F90 +++ b/src/Core/schism_glbl.F90 @@ -153,11 +153,16 @@ module schism_glbl integer,save :: ihfskip,nrec,nspool,ifile,ifile_len, & &noutput,it_main,iths_main,id_out_var(2000),id_out_ww3(100),ncount_2dnode, & &ncount_2delem,ncount_2dside,ncount_3dnode,ncount_3delem,ncount_3dside,nsend_varout + + ! added by Zhiyun Du + integer,save :: ncount_3dbed + integer,save,allocatable :: srqst7(:) real(rkind),save :: time_stamp !simulation time in sec !Send var buffers real(4),save,allocatable :: varout_3dnode(:,:,:),varout_3delem(:,:,:),varout_3dside(:,:,:) real(4),save,allocatable :: varout_2dnode(:,:),varout_2delem(:,:),varout_2dside(:,:) + real(4),save,allocatable :: varout_3dbed(:,:,:) ! add by Zhiyun Du character(len=48),save,allocatable :: outfile_ns(:) character(len=48),save :: a_48 character(len=16),save :: a_16 diff --git a/src/Core/schism_io.F90 b/src/Core/schism_io.F90 index 01dea1443..b4ab252a9 100644 --- a/src/Core/schism_io.F90 +++ b/src/Core/schism_io.F90 @@ -36,7 +36,8 @@ module schism_io &one_dim,two_dim,time_dim,time_dims(1),itime_id,ele_dims(2),x_dims(1), & &y_dims(1),z_dims(1),var2d_dims(2),var3d_dims(3),var4d_dims(4),dummy_dim(1), & &data_start_1d(1),data_start_2d(2),data_start_3d(3),data_start_4d(4), & - &data_count_1d(1),data_count_2d(2),data_count_3d(3),data_count_4d(4) + &data_count_1d(1),data_count_2d(2),data_count_3d(3),data_count_4d(4), & + &nbed_dim ! added by Zhiyun Du integer,save,public :: ncid_schism_io public :: write_obe @@ -360,6 +361,106 @@ subroutine writeout_nc(varid,var_nm,i23d,idim1,idim2,outvar1,outvar2) end subroutine writeout_nc +!=============================================================================== +! add by Zhiyun Du, to output 3D bed variables +subroutine writeout_nc_bed(varid,var_nm,idim1,idim2,outvar1) +!------------------------------------------------------------------------------- +! Netcdf output for element-based bed-layer variables. +! Can be called from any routine, but make sure that +! the calling routine is called inside the main time loop +! exactly ONCE per step! +! +! This routine is intended for variables like: +! bed(Nbed,nea,MBEDP) +! bed(Nbed,nea,i_tracer) +! +! Inputs: +! var_nm: name of the output variable in nc file +! idim1: vertical bed dimension, should be Nbed +! idim2: horizontal dimension, should be nea +! outvar1(idim1,idim2): output array, e.g. bed(:,:,ithck) +! In/out: +! varid: first call will define the variable and return its ID; +! later calls reuse that ID +! +! Output variable dimensions in netCDF: +! (nbed, nele, time) +!------------------------------------------------------------------------------- + implicit none + + character(len=*), intent(in) :: var_nm + integer, intent(in) :: idim1,idim2 + real(rkind), intent(in) :: outvar1(idim1,idim2) + integer, intent(inout) :: varid + + character(len=1000) :: var_nm2 + integer :: i,iret,iret2,irec,len_var,idim2p + real*4 :: a1d(1) + +!----- Return if not output step + if(mod(it_main,nspool)/=0) return + + if(idim1<=0) call parallel_abort('writeout_nc_bed: should be Nbed for idim1') + + irec=(it_main-(ifile-1)*ihfskip)/nspool + if(irec<=0) call parallel_abort('writeout_nc_bed: irec<=0') + + var_nm2=adjustl(var_nm) + len_var=len_trim(var_nm2) + +!----- Dump time + a1d(1)=real(time_stamp) + data_start_1d(1)=irec + data_count_1d(1)=1 + iret=nf90_put_var(ncid_schism_io,itime_id,a1d,data_start_1d,data_count_1d) + +!----- Element-based output + idim2p=ne + if(idim20) then + allocate(var3dbed(Nbed,ne_max,nproc_compute), stat=istat) + if(istat/=0) call parallel_abort('scribe_init: var3dbed') + allocate(var3dbed_gb(Nbed,ne_global), stat=istat) + if(istat/=0) call parallel_abort('scribe_init: var3dbed_gb') + endif +#endif ! call mpi_barrier(comm_scribe,ierr) if(myrank_scribe==0) write(16,*)'finished scribe_init:',myrank_schism @@ -565,6 +587,7 @@ subroutine scribe_step(it) enddo !j itmp5=istart_sed_3dnode+ntrs(5) !index of iof_sed so far + ! total SSC if(iof_sed(itmp5+1)==1) call scribe_recv_write(it,1,1,itotal,icount_out_name) ! itotal=itotal+1 ! icount_out_name=icount_out_name+1 !index into out_name @@ -582,6 +605,15 @@ subroutine scribe_step(it) ! endif !myrank_schism ! endif !iof_ itmp5=itmp5+1 +! output bed thickness 3D, Zhiyun Du + if(iof_sed(itmp5+1)==1) call scribe_recv_write_bed(it,1,itotal,icount_out_name) + itmp5=itmp5+1 + +! bed fraction for each sediment class + do j=1,ntrs(5) + if(iof_sed(itmp5+j)==1) call scribe_recv_write_bed(it,1,itotal,icount_out_name) + enddo + itmp5=itmp5+ntrs(5) #endif /*USE_SED*/ #ifdef USE_ECO @@ -1020,6 +1052,87 @@ subroutine nc_writeout3D(imode,it,idim1,idim2,var3d_gb2,vname,i23d) end subroutine nc_writeout3D +!=============================================================================== +! add a new subroutine to write 3d vars in the bed, Zhiyun Du + subroutine nc_writeout3D_bed(it,idim1,idim2,var3d_bed_gb,vname,i23d) + implicit none + + ! idim1 = Nbed + ! idim2 = ne_global + integer, intent(in) :: it,idim1,idim2,i23d + real(4), intent(in) :: var3d_bed_gb(idim1,idim2) + character(len=*), intent(in) :: vname + + integer :: irec,iret,ivarid,chunks(3) + integer,save :: ivar_id_bed = -999 + character(len=140) :: fname + character(len=12) :: ifile_char + real(rkind) :: a1d(1) + + ! Check inputs + if(idim1/=Nbed) call parallel_abort('nc_writeout3D_bed: idim1/=Nbed') + if(idim2/=ne_global) call parallel_abort('nc_writeout3D_bed: idim2/=ne_global') + + chunks(1)=idim1 + chunks(2)=idim2 + chunks(3)=1 + if(4.d0*chunks(1)*chunks(2)*chunks(3)>3.d9) call parallel_abort('nc_writeout3D_bed: chunk size') + + if(mod(it-nspool,ihfskip)==0) then + ifile=(it-1)/ihfskip+1 + write(ifile_char,'(i12)') ifile + fname=trim(adjustl(out_dir))//'/'//trim(adjustl(vname))//'_'//trim(adjustl(ifile_char))//'.nc' + iret=nf90_create(trim(adjustl(fname)),OR(NF90_NETCDF4,NF90_CLOBBER),ncid_schism_3d) + + call fill_header_static(iof_ugrid,ncid_schism_3d,itime_id,node_dim, & + &nele_dim,nedge_dim,four_dim,nv_dim,one_dim,two_dim,time_dim) + + iret=nf90_redef(ncid_schism_3d) + + ! define bed-layer dimension if not already in header + iret=nf90_def_dim(ncid_schism_3d,'nSCHISM_bed_layers',Nbed,nbed_dim) + if(iret/=NF90_NOERR .and. iret/=NF90_EEXIST) then + call parallel_abort('nc_writeout3D_bed: def nbed_dim') + endif + + var3d_dims(1)=nbed_dim + var3d_dims(2)=nele_dim + var3d_dims(3)=time_dim + + iret=nf90_def_var(ncid_schism_3d,trim(adjustl(vname)),NF90_FLOAT,var3d_dims,ivar_id_bed) + if(iret.ne.NF90_NOERR) call parallel_abort('nc_writeout3D_bed: var_dims') + + iret=nf90_put_att(ncid_schism_3d,ivar_id_bed,'i23d',i23d) + iret=nf90_put_att(ncid_schism_3d,ivar_id_bed,'missing_value',NF90_FILL_FLOAT) + + iret=nf90_def_var_chunking(ncid_schism_3d,ivar_id_bed,NF90_CHUNKED,chunks) + iret=nf90_def_var_deflate(ncid_schism_3d,ivar_id_bed,0,1,4) + call add_mesh_attributes(ncid_schism_3d, ivar_id_bed) + iret=nf90_enddef(ncid_schism_3d) + endif + + irec=(it-(ifile-1)*ihfskip)/nspool + if(irec<=0) call parallel_abort('nc_writeout3D_bed: irec<=0') + + a1d(1)=dble(it*dt) + data_start_1d(1)=irec + data_count_1d(1)=1 + iret=nf90_put_var(ncid_schism_3d,itime_id,a1d,data_start_1d,data_count_1d) + if(iret/=NF90_NOERR) call parallel_abort('nc_writeout3D_bed: put time') + + data_start_3d=(/1,1,irec/) + data_count_3d=(/Nbed,idim2,1/) + if(ivar_id_bed <= 0) call parallel_abort('nc_writeout3D_bed: ivar_id_bed not set') + iret=nf90_put_var(ncid_schism_3d,ivar_id_bed,real(var3d_bed_gb),data_start_3d,data_count_3d) + if(iret/=NF90_NOERR) call parallel_abort('nc_writeout3D_bed: put 3D bed var') + + if(mod(it,ihfskip)==0) then + iret=nf90_close(ncid_schism_3d) + ivar_id_bed = -999 + endif + + end subroutine nc_writeout3D_bed + !=============================================================================== ! Declare static and time vars and fill in static vars ! Outputs: dimension IDs and var ID for time array (itime_id0) @@ -1358,6 +1471,39 @@ subroutine scribe_recv_write(it,imode,ivs,itotal,icount_out_name) end subroutine scribe_recv_write +! ============================================================================== +! added a new subroutine to rev and write 3d vars in sediment bed, Zhiyun Du + subroutine scribe_recv_write_bed(it,ivs,itotal,icount_out_name) + implicit none + integer, intent(in) :: it,ivs + integer, intent(inout) :: itotal,icount_out_name + character(len=20) :: varname3 + integer :: i,m,irank,ierr + + if(ivs/=1.and.ivs/=2) call parallel_abort('scribe_recv_write_bed: ivs') + + do m=1,ivs + itotal=itotal+1 + icount_out_name=icount_out_name+1 + irank=nproc_schism-itotal + + if(myrank_schism==irank) then + do i=1,nproc_compute + call mpi_irecv(var3dbed(:,:,i),ne(i)*Nbed,MPI_REAL4,i-1,200+itotal,comm_schism,rrqst2(i),ierr) + enddo + call mpi_waitall(nproc_compute,rrqst2,MPI_STATUSES_IGNORE,ierr) + + do i=1,nproc_compute + var3dbed_gb(1:Nbed,ielg(1:ne(i),i)) = var3dbed(1:Nbed,1:ne(i),i) + enddo + + varname3=out_name(icount_out_name) + call nc_writeout3D_bed(it,Nbed,ne_global,var3dbed_gb,varname3,iout_23d(icount_out_name)) + endif + enddo + + end subroutine scribe_recv_write_bed + !=============================================================================== subroutine add_mesh_attributes(ncid, varid) implicit none diff --git a/src/Hydro/misc_subs.F90 b/src/Hydro/misc_subs.F90 index 7d8dee09d..ff890cd02 100644 --- a/src/Hydro/misc_subs.F90 +++ b/src/Hydro/misc_subs.F90 @@ -6353,6 +6353,57 @@ subroutine savensend3D_scribe(icount,imode,ivs,nvrt0,npes,savevar1,savevar2) end subroutine savensend3D_scribe +! add by Zhiyun Du, to save 3d var in sedbed and send to scribe +! Save temp 3D sediment-bed vars and send to scribes + subroutine savensend3D_bed_scribe(icount,ivs,Nbed0,nees,savevar1,savevar2) + use schism_glbl, only : rkind,ne,nsend_varout, & + &varout_3dbed,ncount_3dbed,srqst7 + use schism_msgp, only : nscribes,nproc_schism,comm_schism,parallel_abort + + implicit none + include 'mpif.h' + + ! ivs: 1 scalar; 2 vector + ! nees: resident element number, should be ne + integer, intent(in) :: ivs,Nbed0,nees + ! icount: global counter + integer, intent(inout) :: icount + real(rkind), intent(in) :: savevar1(Nbed0,nees) + real(rkind), optional, intent(in) :: savevar2(Nbed0,nees) + + integer :: j,ierr + + ! Check + if(nees/=ne) call parallel_abort('savensend3D_bed_scribe: nees/=ne') + if(ivs/=1.and.ivs/=2) call parallel_abort('savensend3D_bed_scribe: ivs') + + if(ivs==2.and..not.present(savevar2)) then + call parallel_abort('savensend3D_bed_scribe: missing vector component') + endif + + do j=1,ivs ! scalar/vector + icount=icount+1 + nsend_varout=nsend_varout+1 + if(nsend_varout>nscribes.or.icount>ncount_3dbed) then + write(*,*) 'ERROR in savensend3D_bed_scribe before send' + write(*,*) 'icount=',icount,' ncount_3dbed=',ncount_3dbed + write(*,*) 'nsend_varout=',nsend_varout,' nscribes=',nscribes + call parallel_abort('savensend3D_bed_scribe: too many sends') + endif + + if(j==1) then + varout_3dbed(:,:,icount)=savevar1(:,1:nees) + else + varout_3dbed(:,:,icount)=savevar2(:,1:nees) + endif + + call mpi_isend(varout_3dbed(:,1:ne,icount),ne*Nbed0,MPI_REAL4, & + & nproc_schism-nsend_varout,200+nsend_varout,comm_schism, & + & srqst7(nsend_varout),ierr) + enddo !j + + end subroutine savensend3D_bed_scribe + !dir$ attributes forceinline :: signa2 function signa2(x1,x2,x3,y1,y2,y3) !------------------------------------------------------------------------------- diff --git a/src/Hydro/schism_init.F90 b/src/Hydro/schism_init.F90 index 816dbc1a6..b37995674 100644 --- a/src/Hydro/schism_init.F90 +++ b/src/Hydro/schism_init.F90 @@ -442,9 +442,10 @@ subroutine schism_init(iorder,indir,iths,ntime) !To add new outputs, simply make sure the iof_* has sufficient size, and !just add the output statements in _step and flags in param.nml (same !order). Flags for modules other than hydro are only used inside USE_* +! changed to iof_sed(3*sed_class+40) by Zhiyun Du, original is 20 if(iorder==0) then allocate(iof_hydro(40),iof_wwm(40),iof_gen(max(1,ntracer_gen)),iof_age(max(1,ntracer_age)),level_age(ntracer_age/2), & - &iof_sed(3*sed_class+20),iof_eco(max(1,eco_class)),iof_icm_core(17),iof_icm_silica(2),iof_icm_zb(2),iof_icm_ph(4), & + &iof_sed(3*sed_class+40),iof_eco(max(1,eco_class)),iof_icm_core(17),iof_icm_silica(2),iof_icm_zb(2),iof_icm_ph(4), & &iof_icm_srm(4),iof_cos(20),iof_fib(5),iof_sed2d(14),iof_ice(10),iof_mice(10),iof_ana(20),iof_marsh(2),iof_dvd(max(1,ntrs(12))), & !dim of srqst7 increased to account for 2D elem/side etc &srqst7(nscribes+10),veg_vert_z(nbins_veg_vert+1),veg_vert_scale_cd(nbins_veg_vert+1), & @@ -6313,7 +6314,7 @@ subroutine schism_init(iorder,indir,iths,ntime) #endif /*USE_WWM*/ #ifdef USE_SED - do i=7,13 + do i=7,16 ! modified by Zhiyun Du, original 7,13 if(iof_sed(i)/=0) then ncount_2dnode=ncount_2dnode+1 iout_23d(ncount_2dnode)=1 @@ -6332,18 +6333,24 @@ subroutine schism_init(iorder,indir,iths,ntime) out_name(ncount_2dnode)='sedErosionalFlux' case(13) out_name(ncount_2dnode)='sedDepositionalFlux' + case(14) + out_name(ncount_2dnode)='bedWaveStress' + case(15) + out_name(ncount_2dnode)='bedCurrentStress' + case(16) + out_name(ncount_2dnode)='bedWaveCurrentStress' end select endif enddo !i - if(iof_sed(14)/=0) then + if(iof_sed(17)/=0) then ! modified by Zhiyun Du, original 14 ncount_2dnode=ncount_2dnode+2 iout_23d(ncount_2dnode-1:ncount_2dnode)=1 out_name(ncount_2dnode-1)='sedBedloadTransportX' out_name(ncount_2dnode)='sedBedloadTransportY' endif !iof - ised_out_sofar=14 !set output flag index so far + ised_out_sofar=17 !set output flag index so far, modified by Zhiyun Du , original 14 do i=1,ntrs(5) if(iof_sed(i+ised_out_sofar)==1) then !vectors write(ifile_char,'(i12)')i @@ -6660,6 +6667,31 @@ subroutine schism_init(iorder,indir,iths,ntime) ised_out_sofar=ised_out_sofar+1 #endif /*USE_SED*/ +! added by Zhiyun Du, two 3d vars in sediment bed +ncount_3dbed = 0 +#ifdef USE_SED + istart_sed_3dnode = istart_sed_3dnode + if(iof_sed(ised_out_sofar+1)==1) then + ncount_3dbed=ncount_3dbed+1 + counter_out_name=counter_out_name+1 + iout_23d(counter_out_name)=2 + out_name(counter_out_name)='BedThickness_3d' + endif + ised_out_sofar=ised_out_sofar+1 + + do i=1,ntrs(5) + if(iof_sed(i+ised_out_sofar)==1) then + write(ifile_char,'(i12)')i + ifile_char=adjustl(ifile_char); itmp2=len_trim(ifile_char) + ncount_3dbed=ncount_3dbed+1 + counter_out_name=counter_out_name+1 + iout_23d(counter_out_name)=2 + out_name(counter_out_name)='BedFraction_3d_'//ifile_char(1:itmp2) + endif !iof + enddo !i + ised_out_sofar=ised_out_sofar+ntrs(5) !index for iof_sed so far +#endif /*USE_SED*/ + #ifdef USE_ECO do i=1,ntrs(6) if(iof_eco(i)==1) then @@ -6854,12 +6886,17 @@ subroutine schism_init(iorder,indir,iths,ntime) allocate(varout_3delem(nvrt,ne,ncount_3delem),stat=istat) if(istat/=0) call parallel_abort('INIT: 3delem') endif + + if(ncount_3dbed>0) then + allocate(varout_3dbed(Nbed,ne,ncount_3dbed),stat=istat) + if(istat/=0) call parallel_abort('INIT: 3delem_bed') + endif endif !iorder !... Send basic time info to scribes. Make sure the send vars are not altered ! during non-block sends/recv ! Min # of scribes required (all 2D (nodes/elem/side) vars share 1 scribe) - noutvars=ncount_3dnode+ncount_3delem+ncount_3dside+1 + noutvars=ncount_3dnode+ncount_3delem+ncount_3dside+1+ncount_3dbed ! add the last by Zhiyun Du if (noutvars > nscribes) then write(errmsg, '(A,I0,A,A,I0,A)') 'INIT: Too few scribes (', nscribes , '). ', & ' Please specify atleast equal to number of output variables (', noutvars, ')' @@ -6874,6 +6911,12 @@ subroutine schism_init(iorder,indir,iths,ntime) call mpi_send(ncount_2dnode,1,itype,nproc_schism-i,102,comm_schism,ierr) call mpi_send(nc_out,1,itype,nproc_schism-i,103,comm_schism,ierr) call mpi_send(nvrt,1,itype,nproc_schism-i,104,comm_schism,ierr) + + ! add send Nbed by Zhiyun Du, use tag 141 that has not been used +#ifdef USE_SED + call mpi_send(Nbed,1,itype,nproc_schism-i,141,comm_schism,ierr) +#endif + call mpi_send(np_global,1,itype,nproc_schism-i,105,comm_schism,ierr) call mpi_send(ne_global,1,itype,nproc_schism-i,106,comm_schism,ierr) call mpi_send(ns_global,1,itype,nproc_schism-i,107,comm_schism,ierr) @@ -6890,6 +6933,12 @@ subroutine schism_init(iorder,indir,iths,ntime) call mpi_send(ncount_3dnode,1,itype,nproc_schism-i,117,comm_schism,ierr) call mpi_send(ncount_3dside,1,itype,nproc_schism-i,118,comm_schism,ierr) call mpi_send(ncount_3delem,1,itype,nproc_schism-i,119,comm_schism,ierr) + + ! add send ncount_3dbed by Zhiyun Du, use tag 123 +#ifdef USE_SED + call mpi_send(ncount_3dbed,1,itype,nproc_schism-i,123,comm_schism,ierr) +#endif + call mpi_send(iout_23d,counter_out_name,itype,nproc_schism-i,120,comm_schism,ierr) call mpi_send(h0,1,rtype,nproc_schism-i,121,comm_schism,ierr) call mpi_send(ntrs,natrm,itype,nproc_schism-i,122,comm_schism,ierr) @@ -6902,7 +6951,7 @@ subroutine schism_init(iorder,indir,iths,ntime) call mpi_send(iof_gen,max(1,ntrs(3)),itype,nproc_schism-i,130,comm_schism,ierr) call mpi_send(iof_age,max(1,ntrs(4)),itype,nproc_schism-i,131,comm_schism,ierr) - call mpi_send(iof_sed,3*ntrs(5)+20,itype,nproc_schism-i,132,comm_schism,ierr) + call mpi_send(iof_sed,3*ntrs(5)+40,itype,nproc_schism-i,132,comm_schism,ierr) ! changed to +40 by Zhiyun Du call mpi_send(iof_eco,max(1,ntrs(6)),itype,nproc_schism-i,133,comm_schism,ierr) call mpi_send(iof_dvd,max(1,ntrs(12)),itype,nproc_schism-i,134,comm_schism,ierr) call mpi_send(istart_sed_3dnode,1,itype,nproc_schism-i,135,comm_schism,ierr) diff --git a/src/Hydro/schism_step.F90 b/src/Hydro/schism_step.F90 index cb0f0fec2..185536636 100644 --- a/src/Hydro/schism_step.F90 +++ b/src/Hydro/schism_step.F90 @@ -63,7 +63,8 @@ subroutine schism_step(it) USE sed_mod, only : Wsed,Srho,Nbed,MBEDP,bedldu,bedldv,bed,bottom, & &bed_frac,mcoefd,bed_fracn,bed_d50n,bed_taun,& &bedforms_rough,bed_rough,izcr,izsw,izwr,izbld, & - &bed,ithck,iaged,ntr_l,Sd50,eroflxn,depflxn,poron,Qaccun,Qaccvn + &bed,ithck,iaged,ntr_l,Sd50,eroflxn,depflxn,poron,Qaccun,Qaccvn, & + &tau_c_n,tau_w_n,tau_wc_n ! added by Zhiyun Du #endif #ifdef USE_SED2D @@ -129,6 +130,9 @@ subroutine schism_step(it) !ncid_schout,ncid_schout2,ncid_schout3,ncid_schout4,ncid_schout5, & !&ncid_schout6,ncid_schout7,ncid_schout_2,ncid_schout2_2,ncid_schout3_2,ncid_schout4_2,ncid_schout5_2, & !&ncid_schout6_2,ncid_schout7_2, + integer :: icount_bed ! added by Zhiyun Du + real(rkind), allocatable :: tmp_bed(:,:), tmp_bed2(:,:) ! added by Zhiyun Du + integer :: nbed_send_expected, nbed_send_done ! added by Zhiyun Du ! integer :: nstp,nnew !Tsinghua group !1120:close real(rkind) :: cwtmp,cwtmp2,cwtmp3,wtmp1,wtmp2,time,ramp,rampbc,rampwind,rampwafo,dzdx,dzdy, & &dudz,dvdz,dudx,dudx2,dvdx,dvdx2,dudy,dudy2,dvdy,dvdy2, & @@ -9089,11 +9093,22 @@ subroutine schism_step(it) &'SED_eroflx',1,1,npa,eroflxn) ![kg/m/m/s] if(iof_sed(13)==1) call writeout_nc(id_out_var(noutput+17), & &'SED_depflx',1,1,npa,depflxn) ![kg/m/m/s] + +! add below, Zhiyun Du, 03/12/2026 if(iof_sed(14)==1) call writeout_nc(id_out_var(noutput+18), & + &'wave_bed_stress',1,1,npa,tau_w_n*rho0) ![Pa] + if(iof_sed(15)==1) call writeout_nc(id_out_var(noutput+19), & + &'current_bed_stress',1,1,npa,tau_c_n*rho0) ![Pa] + if(iof_sed(16)==1) call writeout_nc(id_out_var(noutput+20), & + &'wave_current_bed_stress',1,1,npa,tau_wc_n*rho0) ![Pa] + + if(iof_sed(17)==1) call writeout_nc(id_out_var(noutput+21), & &'SED_qbdl_acc',1,1,npa,Qaccun,Qaccvn) ![[kg/m/s]] - noutput=noutput+14 - icount=14 !offset + + + noutput=noutput+17 ! original is noutput+14 + icount=17 !offset ! original value is 14 do i=1,ntrs(5) write(it_char,'(i72)')i @@ -9127,8 +9142,25 @@ subroutine schism_step(it) if(iof_sed(icount+1)==1) call writeout_nc(id_out_var(noutput+icount+5), & &'SED_TSC',2,nvrt,npa,total_sus_conc) + noutput=noutput+1 + icount=icount+1 +! added by Zhiyun Du, note "writeout_nc_bed", output bed thickness + if(iof_sed(icount+1)==1) call writeout_nc_bed(id_out_var(noutput+icount+5), & + &'SED_bed_thickness_layer',Nbed,nea,bed(:,:,ithck)) noutput=noutput+1 + icount=icount+1 + +! added by Zhiyun Du, output sed fraction + do i=1,ntrs(5) + write(it_char,'(i72)')i + it_char=adjustl(it_char); lit=len_trim(it_char) + itmp=irange_tr(1,5)+i-1 !global tracer # + if(iof_sed(icount+i)==1) call writeout_nc_bed(id_out_var(noutput+icount+i+4), & + &'SED_bedfrac_layer_'//it_char(1:lit),Nbed,nea,bed_frac(:,:,i)) + enddo !i + noutput=noutput+ntrs(5) + icount=icount+ntrs(5) #endif /*USE_SED*/ #ifdef USE_ECO @@ -9592,7 +9624,7 @@ subroutine schism_step(it) #endif /*USE_WWM*/ #ifdef USE_SED - do i=7,13 + do i=7,16 ! change by Zhiyun Du, original is 7,13 if(iof_sed(i)/=0) then icount=icount+1 if(icount>ncount_2dnode) call parallel_abort('STEP: icount>nscribes(2.0)') @@ -9611,18 +9643,24 @@ subroutine schism_step(it) varout_2dnode(icount,:)=eroflxn(1:np) case(13) varout_2dnode(icount,:)=depflxn(1:np) + case(14) ! add by Zhiyun Du + varout_2dnode(icount,:)=tau_w_n(1:np)*rho0 + case(15) ! add by Zhiyun Du + varout_2dnode(icount,:)=tau_c_n(1:np)*rho0 + case(16) ! add by Zhiyun Du + varout_2dnode(icount,:)=tau_wc_n(1:np)*rho0 end select endif enddo !i - if(iof_sed(14)/=0) then + if(iof_sed(17)/=0) then ! modified by Zhiyun Du, original is 14 icount=icount+2 if(icount>ncount_2dnode) call parallel_abort('STEP: icount>nscribes(2.11)') varout_2dnode(icount-1,:)=Qaccun(1:np) varout_2dnode(icount,:)=Qaccvn(1:np) endif - ised_out_sofar=14 !set output flag index so far + ised_out_sofar=17 !set output flag index so far, original is 14, modified by Zhiyun Du do i=1,ntrs(5) if(iof_sed(i+ised_out_sofar)==1) then !vectors icount=icount+2 @@ -9969,6 +10007,7 @@ subroutine schism_step(it) ! nsend_varout=nsend_varout+1 ! if(nsend_varout>nscribes.or.icount>ncount_3dnode) call parallel_abort('STEP: icount>nscribes(1.81)') itmp=irange_tr(1,5)+i-1 !tracer # +! write(*,*) 'icount_ssc_tracers = ', icount call savensend3D_scribe(icount,1,1,nvrt,np,tr_nd(itmp,:,1:np)) ! varout_3dnode(:,:,icount)=tr_nd(itmp,:,1:np) ! call mpi_isend(varout_3dnode(:,1:np,icount),np*nvrt,MPI_REAL4,nproc_schism-nsend_varout, & @@ -9986,8 +10025,53 @@ subroutine schism_step(it) ! call mpi_isend(varout_3dnode(:,1:np,icount),np*nvrt,MPI_REAL4,nproc_schism-nsend_varout, & ! &200+nsend_varout,comm_schism,srqst7(nsend_varout),ierr) endif + + ised_out_sofar=ised_out_sofar+1 ! added by Zhiyun Du #endif /*USE_SED*/ +! added by Zhiyun Du, to save 3d vars in sediment bed +icount_bed = 0 +#ifdef USE_SED + nbed_send_expected = ncount_3dbed + nbed_send_done = 0 + + ! Allocate temp buffers once if we will send anything + if ( (iof_sed(ised_out_sofar+1)==1) .or. any(iof_sed(ised_out_sofar+2:ised_out_sofar+1+ntrs(5))==1) ) then + if (.not. allocated(tmp_bed)) allocate(tmp_bed(Nbed, ne)) + if (.not. allocated(tmp_bed2)) allocate(tmp_bed2(Nbed, ne)) ! optional second component (future vector) + endif + + ! bed thickness all layers (3D bed) + if (iof_sed(ised_out_sofar+1)==1) then + tmp_bed(:,:) = bed(1:Nbed, 1:ne, ithck) + call savensend3D_bed_scribe(icount_bed, 1, Nbed, ne, tmp_bed) + nbed_send_done = nbed_send_done + 1 + endif + ised_out_sofar = ised_out_sofar + 1 + + ! bed fraction all layers by class (3D bed) + do i=1,ntrs(5) + if (iof_sed(i+ised_out_sofar)==1) then + tmp_bed(:,:) = bed_frac(1:Nbed, 1:ne, i) + call savensend3D_bed_scribe(icount_bed, 1, Nbed, ne, tmp_bed) + nbed_send_done = nbed_send_done + 1 + endif + enddo + ised_out_sofar = ised_out_sofar + ntrs(5) + + ! Sanity check, may not needed + if (nbed_send_expected /= nbed_send_done) then + write(*,*) 'ERROR bedout: expected bed sends=', nbed_send_expected, ' but did=', nbed_send_done + write(*,*) 'ERROR bedout: icount_bed(after)=', icount_bed + call flush(6) + call parallel_abort('bedout: send count mismatch (ncount_3dbed vs actual)') + endif + + ! (Optional) deallocate temps here, or keep them allocated for performance: + if (allocated(tmp_bed)) deallocate(tmp_bed) + if (allocated(tmp_bed2)) deallocate(tmp_bed2) +#endif + #ifdef USE_ECO do i=1,ntrs(6) if(iof_eco(i)==1) then diff --git a/src/Sediment/sed_init.F90 b/src/Sediment/sed_init.F90 index b21f94803..82dbda94c 100644 --- a/src/Sediment/sed_init.F90 +++ b/src/Sediment/sed_init.F90 @@ -219,6 +219,13 @@ SUBROUTINE sed_alloc() IF(i/=0) CALL parallel_abort('sed_alloc: lbc_sed allocation failure') ALLOCATE(bc_sed(npa),stat=i) IF(i/=0) CALL parallel_abort('sed_alloc: bc_sed allocation failure') + ! add below by Zhiyun Du + ALLOCATE(tau_c_n(npa),stat=i) + IF(i/=0) CALL parallel_abort('sed_alloc: tau_c_n allocation failure') + ALLOCATE(tau_w_n(npa),stat=i) + IF(i/=0) CALL parallel_abort('sed_alloc: tau_w_n allocation failure') + ALLOCATE(tau_wc_n(npa),stat=i) + IF(i/=0) CALL parallel_abort('sed_alloc: tau_wc_n allocation failure') ! BM ALLOCATE(poron(npa),stat=i) @@ -312,6 +319,9 @@ SUBROUTINE sed_alloc() bed_d50n(:) = IniVal bed_taun(:) = IniVal bed_rough(:) = IniVal + tau_c_n(:) = IniVal ! add by Zhiyun Du + tau_w_n(:) = IniVal ! add by Zhiyun Du + tau_wc_n(:) = IniVal ! add by Zhiyun Du ! BM poron(:) = IniVal @@ -781,6 +791,9 @@ SUBROUTINE sed_init bed_taun(:) = 0.0d0 bed_fracn(:,:) = 0.0d0 bdfc(:) = 0.0d0 + tau_c_n(:) = 0.0d0 ! add by Zhiyun Du + tau_w_n(:) = 0.0d0 ! add by Zhiyun Du + tau_wc_n(:) = 0.0d0 ! add by Zhiyun Du DO i=1,nea DO j=1,i34(i) DO ised=1,ntr_l diff --git a/src/Sediment/sed_mod.F90 b/src/Sediment/sed_mod.F90 index 744b6c72c..34f9f89d3 100644 --- a/src/Sediment/sed_mod.F90 +++ b/src/Sediment/sed_mod.F90 @@ -267,6 +267,9 @@ MODULE sed_mod REAL(rkind), ALLOCATABLE :: bed_taun(:) !Bottom shear stress (m^2/s/s) REAL(rkind), ALLOCATABLE :: bed_rough(:) !Apparent Roughness length (bedform prediction) REAL(rkind), ALLOCATABLE :: imnp(:) !BM: morphological ramp value (-); imnp(npa) + REAL(rkind), ALLOCATABLE :: tau_c_n(:) !Bottom current shear stress (m^2/s/s), add by Zhiyun Du + REAL(rkind), ALLOCATABLE :: tau_w_n(:) !Bottom wave shear stress (m^2/s/s),add by Zhiyun Du + REAL(rkind), ALLOCATABLE :: tau_wc_n(:) !Bottom current-wave shear stress (m^2/s/s),add by Zhiyun Du ! WWM variables defined at element centres REAL(rkind), ALLOCATABLE :: hs(:) !Significant wave height from WWM (m) ; (hs(nea) diff --git a/src/Sediment/sediment.F90 b/src/Sediment/sediment.F90 index d647d3431..5657d2405 100644 --- a/src/Sediment/sediment.F90 +++ b/src/Sediment/sediment.F90 @@ -1573,6 +1573,34 @@ subroutine sediment(it,moitn,mxitn,rtol,dave,tot_bedmass) ENDDO !i CALL exchange_p2d(bed_taun(:)) +! Added by Zhiyun Du, interpolate tau_c,tau_w,tau_wc to node, for output + tau_c_n = 0.0d0 + tau_w_n = 0.0d0 + tau_wc_n = 0.0d0 + DO i=1,np + IF(idry(i)==1) CYCLE + ta = 0.0d0 + DO j=1,nne(i) + ie=indel(j,i) + IF(idry_e(ie)==0)THEN + ta=ta+area(ie) + tau_c_n(i) = tau_c_n(i) + tau_c(ie) * area(ie) + tau_w_n(i) = tau_w_n(i) + tau_w(ie) * area(ie) + tau_wc_n(i) = tau_wc_n(i) + tau_wc(ie) * area(ie) + ENDIF + ENDDO !j + IF(ta==0.d0)THEN + CALL parallel_abort('SEDIMENT: elem2nod (2)') + ELSE + tau_c_n(i) = tau_c_n(i)/ta + tau_w_n(i) = tau_w_n(i)/ta + tau_wc_n(i) = tau_wc_n(i)/ta + ENDIF + ENDDO !i + CALL exchange_p2d(tau_c_n(:)) + CALL exchange_p2d(tau_w_n(:)) + CALL exchange_p2d(tau_wc_n(:)) + !--------------------------------------------------------------------- ! - BCG: If only one bed layer, re-initialize bed thickness to initial ! thickness and update bed_mass accordingly From cdf8d90532d5920860d52be6e5b578058d0e459f Mon Sep 17 00:00:00 2001 From: Zhiyun Du Date: Wed, 18 Mar 2026 16:22:09 -0400 Subject: [PATCH 2/2] add updated param.nml for 3D bed output --- sample_inputs/param.nml | 159 ++++++++++++++++++++++++---------------- 1 file changed, 94 insertions(+), 65 deletions(-) diff --git a/sample_inputs/param.nml b/sample_inputs/param.nml index aed51ed9f..4cfb99041 100644 --- a/sample_inputs/param.nml +++ b/sample_inputs/param.nml @@ -20,25 +20,25 @@ ibc = 0 !Baroclinic option ibtp = 1 - rnday = 30 !total run time in days - dt = 100. !Time step in sec + rnday = 1459 !total run time in days + dt = 150. !Time step in sec ! Grid for WWM (USE_WWM) msc2 = 24 !same as msc in .nml ... for consitency check between SCHISM and WWM - mdc2 = 30 !same as mdc in .nml + mdc2 = 24 !same as mdc in .nml ! Define # of tracers in tracer modules (if enabled) ntracer_gen = 2 !user defined module (USE_GEN) ntracer_age = 4 !age calculation (USE_AGE). Must be =2*N where N is # of age tracers - sed_class = 5 !SED3D (USE_SED) + sed_class = 7 !SED3D (USE_SED) eco_class = 27 !EcoSim (USE_ECO): must be between [25,60] ! # of vertical bins in vegetation module (used only if iveg/=0) nbins_veg_vert = 2 ! Global output controls - nspool = 36 !output step spool - ihfskip = 864 !stack spool; every ihfskip steps will be put into 1_*, 2_*, etc... + nspool = 24 !output step spool + ihfskip = 576 !stack spool; every ihfskip steps will be put into 1_*, 2_*, etc... / &OPT @@ -90,7 +90,7 @@ loadtide_coef = 0.1 !only used if iloadtide=2,3 (for '3', the default should be 0.12) ! Starting time - start_year = 2000 !int + start_year = 2018 !int start_month = 1 !int start_day = 1 !int start_hour = 0 !double @@ -102,13 +102,13 @@ ! Notes for lon/lat: make sure hgrid.ll and grid in sflux are consistent in ! longitude range! !----------------------------------------------------------------------- - ics = 1 !Coordinate option + ics = 2 !Coordinate option !----------------------------------------------------------------------- ! Hotstart option. 0: cold start; 1: hotstart with time reset to 0; 2: ! continue from the step in hotstart.nc !----------------------------------------------------------------------- - ihot = 0 + ihot = 1 !----------------------------------------------------------------------- ! Equation of State type used @@ -153,7 +153,7 @@ ! indvel=0, ishapiro=1 (shapiro0=0.5), ihorcon=inter_mom=0. ! For applications that include eddying regime, refer to the manual. !----------------------------------------------------------------------- - ihorcon = 0 + ihorcon = 2 hvis_coef0 = 0.025 !const. diffusion # if ihorcon/=0; <=0.025 for ihorcon=2, <=0.125 for ihorcon=1 ! cdh = 0.01 !needed only if ihorcon/=0; land friction coefficient - not active yet @@ -201,8 +201,8 @@ ! force expression. ! To get SCHISM-only result withno feedback from WWM, compile without WWM. !----------------------------------------------------------------------- - icou_elfe_wwm = 0 - nstep_wwm = 1 !call WWM every this many time steps + icou_elfe_wwm = 1 + nstep_wwm = 3 !call WWM every this many time steps iwbl = 0 !wave boundary layer formulation (used only if USE_WMM and !icou_elfe_wwm/=0 and nchi=1. If icou_elfe_wwm=0, set iwbl=0): !1-modified Grant-Madsen formulation; 2-Soulsby (1997) @@ -294,7 +294,7 @@ ! If USE_NWM_BMI is on with if_source/=0, some parts of the reading and some b.c. will be bypassed (but the inputs ! will still be needed). !----------------------------------------------------------------------- - if_source = 0 + if_source = -1 dramp_ss = 2 !needed if if_source/=0; ramp-up period in days for source/sinks (no ramp-up if <=0) !---------------------------------------------------------------------- @@ -353,7 +353,7 @@ ! the latitude specified in CPP projection ('sfea0') is used. If ncor=1 and ics=2, ! Coriolis is calculated from local latitude, and 'sfea0' is not used. !----------------------------------------------------------------------- - ncor = 0 !should usually be 1 if ics=2 + ncor = 1 !should usually be 1 if ics=2 rlatitude = 46 !if ncor=-1 coricoef = 0 !if ncor=0 @@ -424,7 +424,7 @@ ! Max. horizontal velocity magnitude, used mainly to prevent problem in ! bulk aerodynamic module !----------------------------------------------------------------------- - rmaxvel = 5. + rmaxvel = 10. !----------------------------------------------------------------------- ! Following parameters control backtracking @@ -480,7 +480,7 @@ ! in backtracking as above). dtb_max/min are the max/min steps allowed - ! actual step is calculated adaptively based on local gradient. !----------------------------------------------------------------------- - nadv = 1 + nadv = 2 dtb_max = 30. !in sec dtb_min = 10. @@ -510,7 +510,7 @@ eps2_tvd_imp = 1.e-14 !Optional hybridized ELM transport for efficiency - ielm_transport = 0 !1: turn on + ielm_transport = 1 !1: turn on max_subcyc = 10 !used only if ielm_transport/=0. Max # of subcycling per time step in transport allowed !if itr_met = 4, the following parameters are needed @@ -561,7 +561,7 @@ ! If iwind_form=-3, the stress is calculated according to the param. of Donelan et al. (1993) based on the wave age. ! In all cases, if USE_ICE the stress in ice-covered portion is calculated by ICE routine. !----------------------------------------------------------------------- - nws = 0 + nws = 2 wtiminc = 150. !time step for atmos. forcing. Default: same as dt drampwind = 1. !ramp-up period in days for wind (no ramp-up if <=0) iwindoff = 0 !needed only if nws/=0; '1': needs windfactor.gr3 @@ -584,8 +584,8 @@ ! hmin_airsea_ex, hmin_salt_ex: ! 0.2 m is recommended for both "1" and "2" !----------------------------------------------------------------------- - ihconsv = 0 !heat exchange option - isconsv = 0 !evaporation/precipitation model + ihconsv = 1 !heat exchange option + isconsv = 1 !evaporation/precipitation model i_hmin_airsea_ex = 2 ! no effect if ihconsv=0 hmin_airsea_ex = 0.2 ![m], no effect if ihconsv=0 i_hmin_salt_ex = 2 ! no effect if isconsv=0 @@ -632,8 +632,8 @@ ! Code will ignore junk values (<=-99) inside [MOD]_nu.nc, so 1 way to avoid ! nudging for a tracer is to set its nudged values to -9999 !----------------------------------------------------------------------- - inu_tr(1) = 0 !T - inu_tr(2) = 0 !S + inu_tr(1) = 2 !T + inu_tr(2) = 2 !S inu_tr(3) = 0 !GEN inu_tr(4) = 0 !Age inu_tr(5) = 0 !SED3D @@ -679,8 +679,8 @@ ! If error occurs like 'bktrk_subs: overflow' or 'MAIN: nbtrk > mxnbt' ! gradually increasing these will solve the problem !----------------------------------------------------------------------- - s1_mxnbt = 0.5 - s2_mxnbt = 3.5 + s1_mxnbt = 1.0 + s2_mxnbt = 4.0 !----------------------------------------------------------------------- ! Flag for harmonic analysis for elevation. If used , need to turn on USE_HA @@ -834,8 +834,8 @@ !----------------------------------------------------------------------- ! Option for hotstart outputs !----------------------------------------------------------------------- - nhot = 0 !1: output *_hotstart every 'nhot_write' steps - nhot_write = 8640 !must be a multiple of ihfskip if nhot=1 + nhot = 1 !1: output *_hotstart every 'nhot_write' steps + nhot_write = 576 !must be a multiple of ihfskip if nhot=1 !----------------------------------------------------------------------- ! Station output option. If iout_sta/=0, need output skip (nspool_sta) and @@ -862,15 +862,15 @@ iof_hydro(11) = 0 !evaporation rate [kg/m/m/s] {evaporationRate} 2D iof_hydro(12) = 0 !precipitation rate [kg/m/m/s] {precipitationRate} 2D iof_hydro(13) = 0 !Bottom stress vector [kg/m/s^2(Pa)] {bottomStressX,Y} 2D vector - iof_hydro(14) = 0 !wind velocity vector [m/s] {windSpeedX,Y} 2D vector + iof_hydro(14) = 1 !wind velocity vector [m/s] {windSpeedX,Y} 2D vector iof_hydro(15) = 0 !wind stress vector [m^2/s/s] {windStressX,Y} 2D vector - iof_hydro(16) = 0 !depth-averaged vel vector [m/s] {depthAverageVelX,Y} 2D vector + iof_hydro(16) = 1 !depth-averaged vel vector [m/s] {depthAverageVelX,Y} 2D vector iof_hydro(17) = 0 !vertical velocity [m/s] {verticalVelocity} 3D - iof_hydro(18) = 0 !water temperature [C] {temperature} 3D - iof_hydro(19) = 0 !water salinity [PSU] {salinity} 3D + iof_hydro(18) = 1 !water temperature [C] {temperature} 3D + iof_hydro(19) = 1 !water salinity [PSU] {salinity} 3D iof_hydro(20) = 0 !water density [kg/m^3] {waterDensity} 3D - iof_hydro(21) = 0 !vertical eddy diffusivity [m^2/s] {diffusivity} 3D - iof_hydro(22) = 0 !vertical eddy viscosity [m^2/s] {viscosity} 3D + iof_hydro(21) = 1 !vertical eddy diffusivity [m^2/s] {diffusivity} 3D + iof_hydro(22) = 1 !vertical eddy viscosity [m^2/s] {viscosity} 3D iof_hydro(23) = 0 !turbulent kinetic energy {turbulentKineticEner} 3D iof_hydro(24) = 0 !turbulent mixing length [m] {mixingLength} 3D ! iof_hydro(25) = 1 !z-coord {zCoordinates} 3D - this flag should be on for visIT etc @@ -879,7 +879,7 @@ iof_hydro(28) = 0 !vertical vel. @elem [m/s] {verticalVelAtElement} 3D iof_hydro(29) = 0 !T @prism centers [C] {temperatureAtElement} 3D iof_hydro(30) = 0 !S @prism centers [PSU] {salinityAtElement} 3D - iof_hydro(31) = 0 !Barotropic pressure gradient force vector (m.s-2) @side centers {pressure_gradient} 2D vector + iof_hydro(31) = 1 !Barotropic pressure gradient force vector (m.s-2) @side centers {pressure_gradient} 2D vector !----------------------------------------------------------------------- ! Outputs from optional modules. Only uncomment these if the USE_* is on @@ -892,15 +892,15 @@ !----------------------------------------------------------------------- ! Outputs from WWM (USE_WWM must be on in Makefile) !----------------------------------------------------------------------- -! iof_wwm(1) = 0 !sig. height (m) {sigWaveHeight} 2D -! iof_wwm(2) = 0 !Mean average period (sec) - TM01 {meanWavePeriod} 2D + iof_wwm(1) = 1 !sig. height (m) {sigWaveHeight} 2D + iof_wwm(2) = 1 !Mean average period (sec) - TM01 {meanWavePeriod} 2D ! iof_wwm(3) = 0 !Zero down crossing period for comparison with buoy (s) - TM02 {zeroDowncrossPeriod} 2D ! iof_wwm(4) = 0 !Average period of wave runup/overtopping - TM10 {TM10} 2D ! iof_wwm(5) = 0 !Mean wave number (1/m) {meanWaveNumber} 2D ! iof_wwm(6) = 0 !Mean wave length (m) {meanWaveLength} 2D -! iof_wwm(7) = 0 !Mean average energy transport direction (degr) - MWD in NDBC? {meanWaveDirection} 2D + iof_wwm(7) = 1 !Mean average energy transport direction (degr) - MWD in NDBC? {meanWaveDirection} 2D ! iof_wwm(8) = 0 !Mean directional spreading (degr) {meanDirSpreading} 2D -! iof_wwm(9) = 0 !Discrete peak period (sec) - Tp {peakPeriod} 2D + iof_wwm(9) = 1 !Discrete peak period (sec) - Tp {peakPeriod} 2D ! iof_wwm(10) = 0 !Continuous peak period based on higher order moments (sec) {continuousPeakPeriod} 2D ! iof_wwm(11) = 0 !Peak phase vel. (m/s) {peakPhaseVel} 2D ! iof_wwm(12) = 0 !Peak n-factor {peakNFactor} 2D @@ -911,7 +911,7 @@ ! iof_wwm(17) = 0 !Peak directional spreading {peakSpreading} 2D ! iof_wwm(18) = 0 !Discrete peak direction (radian?) {discretePeakDirectio} 2D ! iof_wwm(19) = 0 !Orbital vel. (m/s) {orbitalVelocity} 2D -! iof_wwm(20) = 0 !RMS Orbital vel. (m/s) {rmsOrbitalVelocity} 2D + iof_wwm(20) = 1 !RMS Orbital vel. (m/s) {rmsOrbitalVelocity} 2D ! iof_wwm(21) = 0 !Bottom excursion period (sec?) {bottomExcursionPerio} 2D ! iof_wwm(22) = 0 !Bottom wave period (sec) {bottomWavePeriod} 2D ! iof_wwm(23) = 0 !Uresell number based on peak period {UresellNumber} 2D @@ -921,7 +921,7 @@ ! iof_wwm(27) = 0 !Roller energy dissipation rate (W/m²) @nodes {Drol} 2D ! iof_wwm(28) = 0 !Total wave energy dissipation rate by depth-induced breaking (W/m²) @nodes {wave_sbrtot} 2D -! iof_wwm(29) = 0 !Total wave energy dissipation rate by bottom friction (W/m²) @nodes {wave_sbftot} 2D + iof_wwm(29) = 1 !Total wave energy dissipation rate by bottom friction (W/m²) @nodes {wave_sbftot} 2D ! iof_wwm(30) = 0 !Total wave energy dissipation rate by whitecapping (W/m²) @nodes {wave_sdstot} 2D ! iof_wwm(31) = 0 !Total wave energy dissipation rate by vegetation (W/m²) @nodes {wave_svegtot} 2D ! iof_wwm(32) = 0 !Total wave energy input rate from atmospheric forcing (W/m²) @nodes {wave_sintot} 2D @@ -951,32 +951,61 @@ ! Specific outputs in SED3D (USE_SED must be on in Makefile; ! otherwise these are not needed) !----------------------------------------------------------------------- -! iof_sed(1) = 0 ! total bed thickness @elem (m) {sedBedThickness} 2D -! iof_sed(2) = 0 ! total bed age over all layers @elem (sec) {sedBedAge} 2D -! iof_sed(3) = 0 ! Sediment transport roughness length @elem (m) (sedTransportRough) {z0st} 2D -! iof_sed(4) = 0 !current-ripples roughness length @elem (m) (sedRoughCurrentRippl) {z0cr} 2D -! iof_sed(5) = 0 !sand-waves roughness length (m) @elem (z0sw_elem) {sedRoughSandWave} 2D -! iof_sed(6) = 0 !wave-ripples roughness length @elem (m) (z0wr_elem) {sedRoughWaveRipple} 2D - -! iof_sed(7) = 0 !bottom depth _change_ from init. condition (m) {sedDepthChange} 2D -! iof_sed(8) = 0 ! Bed median grain size in the active layer (mm) {sedD50} 2D -! iof_sed(9) = 0 ! Bottom shear stress (Pa) {sedBedStress} 2D -! iof_sed(10) = 0 ! Bottom roughness lenghth (mm) {sedBedRoughness} 2D -! iof_sed(11) = 0 ! sediment porosity in the top layer (-) {sedPorocity} 2D -! iof_sed(12) = 0 ! erosion flux for suspended transport (kg/m/m/s) {sedErosionalFlux} 2D -! iof_sed(13) = 0 ! deposition flux for suspended transport (kg/m/m/s) {sedDepositionalFlux} 2D -! iof_sed(14) = 0 ! Bedload transport rate vector due to wave acceleration (kg/m/s) {sedBedloadTransportX,Y} 2D vector - -! Example of using 2 classes -! iof_sed(15) = 0 !Bedload transport rate _vector_ (kg.m-1.s-1) for 1st tracer (one output per class) {sedBedload[X,Y]_1} 2D vector -! iof_sed(16) = 0 !Bedload transport of 2nd class {sedBedFraction_[X,Y]_2} 2D vector -! iof_sed(17) = 0 !Bed fraction 1st tracer (one output per class) [-] {sedBedFraction_1} 2D -! iof_sed(18) = 0 !Bed fraction of 2nd class {sedBedFraction_2} 2D - -! iof_sed(19) = 0 !conc. of 1st class (one output need by each class) [g/L] {sedConcentration_1} 3D -! iof_sed(20) = 0 !conc. of 2nd class {sedConcentration_2} 3D - -! iof_sed(21) = 0 !total suspended concentration (g/L) {totalSuspendedLoad} 3D + iof_sed(1) = 1 ! total bed thickness @elem (m) {sedBedThickness} 2D + iof_sed(2) = 0 ! total bed age over all layers @elem (sec) {sedBedAge} 2D + iof_sed(3) = 0 ! Sediment transport roughness length @elem (m) (sedTransportRough) {z0st} 2D + iof_sed(4) = 0 !current-ripples roughness length @elem (m) (sedRoughCurrentRippl) {z0cr} 2D + iof_sed(5) = 0 !sand-waves roughness length (m) @elem (z0sw_elem) {sedRoughSandWave} 2D + iof_sed(6) = 0 !wave-ripples roughness length @elem (m) (z0wr_elem) {sedRoughWaveRipple} 2D + + iof_sed(7) = 0 !bottom depth _change_ from init. condition (m) {sedDepthChange} 2D + iof_sed(8) = 0 ! Bed median grain size in the active layer (mm) {sedD50} 2D + iof_sed(9) = 0 ! Bottom shear stress (Pa) {sedBedStress} 2D + iof_sed(10) = 0 ! Bottom roughness lenghth (mm) {sedBedRoughness} 2D + iof_sed(11) = 0 ! sediment porosity in the top layer (-) {sedPorocity} 2D + iof_sed(12) = 1 ! erosion flux for suspended transport (kg/m/m/s) {sedErosionalFlux} 2D + iof_sed(13) = 1 ! deposition flux for suspended transport (kg/m/m/s) {sedDepositionalFlux} 2D + iof_sed(14) = 1 ! Bottom wave shear stress (Pa) {sedBedStress} 2D + iof_sed(15) = 1 ! Bottom current shear stress (Pa) {sedBedStress} 2D + iof_sed(16) = 1 ! Bottom wave-current shear stress (Pa) {sedBedStress} 2D + iof_sed(17) = 0 ! Bedload transport rate vector due to wave acceleration (kg/m/s) {sedBedloadTransportX,Y} 2D vector + +! Example of using 7 classes + iof_sed(18) = 0 !Bedload transport rate _vector_ (kg.m-1.s-1) for 1st tracer (one output per class) {sedBedload[X,Y]_1} 2D vector + iof_sed(19) = 0 !Bedload transport of 2nd class {sedBedFraction_[X,Y]_2} 2D vector + iof_sed(20) = 0 !Bedload transport of 3rd class {sedBedFraction_[X,Y]_2} 2D vector + iof_sed(21) = 0 !Bedload transport of 4th class {sedBedFraction_[X,Y]_2} 2D vector + iof_sed(22) = 0 !Bedload transport of 5th class {sedBedFraction_[X,Y]_2} 2D vector + iof_sed(23) = 0 !Bedload transport of 6th class {sedBedFraction_[X,Y]_2} 2D vector + iof_sed(24) = 0 !Bedload transport of 7th class {sedBedFraction_[X,Y]_2} 2D vector + + iof_sed(25) = 0 !Bed fraction 1st tracer (one output per class) [-] {sedBedFraction_1} 2D + iof_sed(26) = 0 !Bed fraction of 2nd class {sedBedFraction_2} 2D + iof_sed(27) = 0 !Bed fraction of 3rd class {sedBedFraction_2} 2D + iof_sed(28) = 0 !Bed fraction of 4th class {sedBedFraction_2} 2D + iof_sed(29) = 0 !Bed fraction of 5th class {sedBedFraction_2} 2D + iof_sed(30) = 0 !Bed fraction of 6th class {sedBedFraction_2} 2D + iof_sed(31) = 0 !Bed fraction of 7th class {sedBedFraction_2} 2D + + iof_sed(32) = 1 !conc. of 1st class (one output need by each class) [g/L] {sedConcentration_1} 3D + iof_sed(33) = 1 !conc. of 2nd class {sedConcentration_2} 3D + iof_sed(34) = 1 !conc. of 3rd class {sedConcentration_2} 3D + iof_sed(35) = 1 !conc. of 4th class {sedConcentration_2} 3D + iof_sed(36) = 1 !conc. of 5th class {sedConcentration_2} 3D + iof_sed(37) = 1 !conc. of 6th class {sedConcentration_2} 3D + iof_sed(38) = 1 !conc. of 7th class {sedConcentration_2} 3D + + iof_sed(39) = 0 !total suspended concentration (g/L) {totalSuspendedLoad} 3D + + iof_sed(40) = 1 !bed thickness for all layers (m) 3D + + iof_sed(41) = 1 !bed fraction for all layers 1st tracer 3D + iof_sed(42) = 1 !bed fraction for all layers 2nd tracer 3D + iof_sed(43) = 1 !bed fraction for all layers 3rd tracer 3D + iof_sed(44) = 1 !bed fraction for all layers 4th tracer 3D + iof_sed(45) = 1 !bed fraction for all layers 5th tracer 3D + iof_sed(46) = 1 !bed fraction for all layers 6th tracer 3D + iof_sed(47) = 1 !bed fraction for all layers 7th tracer 3D !----------------------------------------------------------------------- ! EcoSim outputs