!/===========================================================================/ ! Copyright (c) 2007, The University of Massachusetts Dartmouth ! Produced at the School of Marine Science & Technology ! Marine Ecosystem Dynamics Modeling group ! All rights reserved. ! ! FVCOM has been developed by the joint UMASSD-WHOI research team. For ! details of authorship and attribution of credit please see the FVCOM ! technical manual or contact the MEDM group. ! ! ! This file is part of FVCOM. For details, see http://fvcom.smast.umassd.edu ! The full copyright notice is contained in the file COPYRIGHT located in the ! root directory of the FVCOM code. This original header must be maintained ! in all distributed versions. ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" ! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, ! THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR ! PURPOSE ARE DISCLAIMED. ! !/---------------------------------------------------------------------------/ ! CVS VERSION INFORMATION ! $Id$ ! $Name$ ! $Revision$ !/===========================================================================/ !======================================================================= ! FVCOM Sediment Module ! ! Copyright: 2005(c) ! ! THIS IS A DEMONSTRATION RELEASE. THE AUTHOR(S) MAKE NO REPRESENTATION ! ABOUT THE SUITABILITY OF THIS SOFTWARE FOR ANY OTHER PURPOSE. IT IS ! PROVIDED "AS IS" WITHOUT EXPRESSED OR IMPLIED WARRANTY. ! ! THIS ORIGINAL HEADER MUST BE MAINTAINED IN ALL DISTRIBUTED ! VERSIONS. ! ! Contact: G. Cowles ! School for Marine Science and Technology, Umass-Dartmouth ! ! Based on the Community Sediment Transport Model (CSTM) as implemented ! in ROMS by J. Warner (USGS) ! ! Comments: Sediment Dynamics Module ! ! Current FVCOM dependency ! ! init_sed.F: - user defined sediment model initial conditions ! mod_ncdio.F: - netcdf output includes concentration/bottom/bed fields ! fvcom.F: - main calls sediment setup ! internal.F: - calls sediment advance ! ! History ! Feb 7, 2008: added initialization of bottom(:,:) to 0 (w/ T. Hamada) ! : fixed loop bounds in hot start and archive for conc (w/ T. Hamada) ! : added comments describing theoretical bases of dynamics ! Feb 14,2008: added non-constant settling velocity for cohesive sediments (w/ T. Hamada) ! : updated vertical flux routine to handle non-constant vertical velocity (w/ T. Hamada) ! : added a user-defined routine to calculate settling velocity based on concentration (w/ T. Hamada) ! : added a user-defined routine to calculate erosion for a general case (w/ T. Hamada) ! ! PLEASE NOTE!!!!!!!!!!! ! Do NOT USE INTEL FORTRAN COMPILER VERSION 11.0 IT HAS KNOWN BUGS WHEN DEALING WITH TYPES WITH ALLOCATABLE ! COMPONENTS. YOU WILL SEE WEIRD BEHAVIOR. VERSION 11.1 IS OK. ! ! ! Later ! 1.) Modify vertical flux routines to work with general vertical coordinate ! 2.) Add divergence term for bedload transport calc ! 3.) Add ripple roughness calculation ! 4.) Add morphological change (bathymetry + vertical velocity condition) ! 5.) Eliminate excess divisions and recalcs ! !======================================================================= Module Mod_Sed_CSTMS #if defined (SEDIMENT) && (CSTMS_SED) Use Mod_Par Use Mod_Prec Use Mod_Types Use Mod_wd Use Control, only : seddis Use all_vars,only : CNSTNT,UNIFORM,SEDIMENT_PARAMETER_TYPE # if defined (WAVE_CURRENT_INTERACTION) USE MOD_WAVE_CURRENT_INTERACTION # endif # if defined (MULTIPROCESSOR) Use Mod_Par # endif # if defined (FLUID_MUD) USE MOD_FLUID_MUD # endif Use Mod_CSTMS_vars Use Mod_FlocMod implicit none contains !========================================================================== ! Allocate Data, Initialize Concentration Fields and Bed Parameters !========================================================================== Subroutine Setup_Sed(fname,Nin,Nin_gl,N_sed,N_sed_Max,Sed_Names) use control, only : msr,casename,input_dir use lims, only : m Use Mod_Utils, only: pstop implicit none integer, intent(in) :: Nin,Nin_gl character(len=80) :: fname integer, intent(out) :: N_sed integer, intent(in ) :: N_sed_max character(len=20) :: Sed_Names(N_sed_max) !-Local-------------------------- logical :: fexist integer :: ised integer :: i,k if(dbg_set(dbg_sbr)) write(ipt,*) "Start: setup_sed" Nobc = Nin Nobc_gl = Nin_gl !---------------------------------------------------------------- ! check arguments and ensure necessary files exist !---------------------------------------------------------------- !initialize sediment model iteration counter sed_its = 0 !ensure sediment parameter file exists if(fname == "none" .or. fname == "NONE" .or. fname == "None")then sedfile = "./"//trim(input_dir)//"/"//trim(casename)//'_sediment.inp' else sedfile = "./"//trim(input_dir)//"/"//trim(fname) endif inquire(file=trim(sedfile),exist=fexist) if(.not.fexist)then write(*,*)'sediment parameter file: ',trim(sedfile),' does not exist' write(*,*)'stopping' call pstop end if !---------------------------------------------------------------- ! read in sediment parameter file !---------------------------------------------------------------- call read_sed_params ! J. Ge for fluid mud layer # if defined (FLUID_MUD) !---------------------------------------------------------------- ! read in fluid mud layer parameter file !---------------------------------------------------------------- call initialize_fluid_mud # endif ! J. Ge for fluid mud layer !---------------------------------------------------------------- ! allocate data !---------------------------------------------------------------- call alloc_sed_vars call alloc_sedbed if(SED_FLOCS)then CALL allocate_sedflocs CALL initialize_sedflocs_param (SEDFLOCS % f_mass, & & SEDFLOCS % f_diam, & & SEDFLOCS % f_g1_sh, & & SEDFLOCS % f_g1_ds, & & SEDFLOCS % f_g3, & & SEDFLOCS % f_l1_sh, & & SEDFLOCS % f_l1_ds, & & SEDFLOCS % f_coll_prob_sh, & & SEDFLOCS % f_coll_prob_ds, & & SEDFLOCS % f_l3) end if !---------------------------------------------------------------- ! initialize fields ! these will be overwritten if this is a restart case !---------------------------------------------------------------- call init_sed !---------------------------------------------------------------- ! setup open boundary condition nudging values !---------------------------------------------------------------- if(sed_nudge) call setup_sed_obc !---------------------------------------------------------------- ! setup point source river forcing data !---------------------------------------------------------------- ! if(numqbc_gl .and. sed_source) call setup_sed_PTsource ! if(sed_source) call setup_sed_PTsource !-------------------------------------------------- !Calculate Bottom Thicknesses !-------------------------------------------------- do i=1,m do k=1,nbed bottom(i,nthck) = bottom(i,nthck) + bed(i,k,ithck) end do end do !---------------------------------------------------------------- ! report sediment parameters and statistics to screen !---------------------------------------------------------------- call report_sed_setup !---------------------------------------------------------------- ! update bottom properties !---------------------------------------------------------------- call update_bottom_properties !---------------------------------------------------------------- ! set arguments for sed names and N_sed for FVCOM to use ! in setting up river forcing !---------------------------------------------------------------- N_sed = Nst do i=1,Nst sed_names(i) = sed(i)%sname end do if(dbg_set(dbg_sbr)) write(ipt,*) "End: setup_sed" End Subroutine Setup_Sed !======================================================================= ! ! ! This routine it is the main driver for the sediment-transport ! ! model. Currently, it includes calls to the following routines: ! ! ! ! * Vertical settling of sediment in the water column. ! ! * Erosive and depositional flux interactions of sediment ! ! between water column and the bed. ! ! * Transport of multiple grain sizes. ! ! * Bed layer stratigraphy. ! ! * Bed morphology. ! ! * Bedload based on Meyer Peter Mueller. ! ! * Bedload based on Soulsby combined waves + currents ! ! (p166 Soulsby 1997) ! ! * Bedload slope term options: Nemeth et al, 2006, Coastal ! ! Engineering, v 53, p 265-275; Lesser et al, 2004, Coastal ! ! Engineering, v 51, p 883-915. ! ! ! ! * Seawater/sediment vertical level distribution: ! ! ! ! W-level RHO-level ! ! ! ! N _________ ! ! | | ! ! | N | ! ! N-1 |_________| S ! ! | | E ! ! | N-1 | A ! ! 2 |_________| W ! ! | | A ! ! | 2 | T ! ! 1 |_________| E ! ! | | R ! ! | 1 | ! ! 0 |_________|_____ bathymetry ! ! |/////////| ! ! | 1 | ! ! 1 |_________| S ! ! | | E ! ! | 2 | D ! ! 2 |_________| I ! ! | | M ! ! | Nbed-1 | E ! ! Nbed-1 |_________| N ! ! | | T ! ! | Nbed | ! ! Nbed |_________| ! ! ! ! References: ! ! ! ! Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G. ! ! Arango, 2008: Development of a three-dimensional, regional, ! ! coupled wave, current, and sediment-transport model, Computers ! ! & Geosciences, 34, 1284-1306. ! ! ! !======================================================================= ! !*********************************************************************** !======================================================================= ! Advance Sediment Model by time step DTSED !======================================================================= Subroutine Advance_Sed(DTin,Tin,taub_in) Use Input_Util Use Scalar Use Control, only : ireport,iint,msr,par Use Lims, only : m,mt,nt,kbm1,numqbc,kb,nprocs,myid Use All_Vars,only : BACKWARD_ADVECTION,BACKWARD_STEP ! Use Mod_OBCS,only : iobcn # if defined (MULTIPROCESSOR) Use Mod_Par # endif ! ! Imported variable declarations. ! Implicit None real(sp), intent(in ) :: DTin,Tin real(sp), intent(in ) :: taub_in(m) integer :: ierr if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Advance_Sed_CSTMS " !------------------------------------------------------ !Check for Initialization !------------------------------------------------------ if(iint < sed_start) return nstp = 1+MOD(iint-sed_start,2) nnew = 3-nstp !------------------------------------------------------ !Increment Sed Model Counter and Report Sed Model Init !------------------------------------------------------ sed_its = sed_its + 1 if(msr .and. sed_its == 1)then write(*,*)'========Sed Model Initiated=======' endif !------------------------------------------------------ !Set up Model Forcing -> time step and bottom shear !Convert taux/tauy to node-based values !Convert Model tau/rho to tau [N/m^2] !------------------------------------------------------ DTsed = DTin T_model = Tin # if defined (WAVE_CURRENT_INTERACTION) if(MB_BBL_USE.or.SG_BBL_USE.or.SSW_BBL_USE)then taub = taucwmax*rho_water else taub = taub_in*rho_water endif # else taub = taub_in*rho_water # endif !------------------------------------------------------ !set bed age to current model time if first sed call !------------------------------------------------------ if(sed_its ==1)then bed(:,:,iaged) = T_model endif !------------------------------------------------------ !Calculate Active Layer Thickness --> bottom(:,iactv) !------------------------------------------------------ call calc_active_layer !----------------------------------------------------------------------- ! Compute sediment bedload transport. !----------------------------------------------------------------------- ! if(BEDLOAD)then CALL calc_sed_bedload endif ! !----------------------------------------------------------------------- ! Compute sediment suspended transport. !----------------------------------------------------------------------- if(SUSLOAD)then ! !----------------------------------------------------------------------- ! Compute sediment flocculation !----------------------------------------------------------------------- ! if(SED_FLOCS)then CALL calc_sed_flocmod endif ! !----------------------------------------------------------------------- ! Compute sediment vertical settling. !----------------------------------------------------------------------- ! CALL calc_sed_settling ! !----------------------------------------------------------------------- ! Compute bed-water column exchanges: erosion and deposition. !----------------------------------------------------------------------- ! CALL calc_sed_fluxes ! !----------------------------------------------------------------------- ! Compute sediment due to biomass. !----------------------------------------------------------------------- ! if(SED_BIOMASS)then CALL calc_sed_biomass endif !----------------------------------------------------------------------- ! Compute fluid mud dynamics. !----------------------------------------------------------------------- ! # if defined (FLUID_MID) call calc_fluid_mud # endif endif ! end of suspended load ! !----------------------------------------------------------------------- ! Compute sediment bed stratigraphy. !----------------------------------------------------------------------- ! if (COHESIVE_BED .or. MIXED_BED)then CALL calc_sed_bed_cohesive elseif(NONCOHESIVE_BED2)then CALL calc_sed_bed2 else CALL calc_sed_bed endif ! !----------------------------------------------------------------------- ! Compute sediment bed biodiffusivity. !----------------------------------------------------------------------- ! if(SED_BIODIFF)then CALL calc_sed_biodiff endif ! !----------------------------------------------------------------------- ! Compute sediment surface layer properties. !----------------------------------------------------------------------- ! CALL Update_Bottom_Properties ! !----------------------------------------------------------------------- ! Horizontal Advection of Sediment Concentration !----------------------------------------------------------------------- ! CALL calc_sed_advection ! !------------------------------------------------------ !Vertical Diffusion of Sediment Concentrations !------------------------------------------------------ ! CALL calc_sed_diffusion !------------------------------------------------------ ! Set Open Boundary Conditions !------------------------------------------------------ CALL calc_sed_bcond !------------------------------------------------------ ! Update Depth Delta and Morphology Array ! this array is positive if material has accumulated, ! negative if eroded !------------------------------------------------------ call update_thickness_delta !------------------------------------------------------ ! Update Concentration Variables to next time level !------------------------------------------------------ ! CALL finalize_sed_step !------------------------------------------------------ ! Report the status of sediment !------------------------------------------------------ ! if(mod(sed_its,n_report)==0)call sed_report RETURN END SUBROUTINE Advance_Sed !======================================================================= !Read Sediment Parameters From Sediment Input File !======================================================================= Subroutine Read_Sed_Params Use Input_Util Implicit None integer linenum,i,k1,iscan real(sp) :: ftemp character(len=120) :: stemp real(sp),allocatable:: ncstmp(:) real(sp),allocatable:: nnstmp(:) real(sp),allocatable:: nsttmp(:) logical, allocatable:: nswitch(:) character(len=80),allocatable:: strtmp(:) ! ! Imported variable declarations ! ! ! Local variable declarations. ! integer :: iTrcStr, iTrcEnd integer :: ifield, igrid, itracer, itrc, nline !----------------------------------------------------------------------- ! Initialize. !----------------------------------------------------------------------- ! if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Read_Sed_Params " ! !----------------------------------------------------------------------- ! Read in cohesive and non-cohesive model parameters. !----------------------------------------------------------------------- ! !read in interation interval for reporting Call Get_Val(n_report,sedfile,'N_REPORT',line=linenum,echo=.true.) !read in Depositional bed layer thickness criteria to create a new layer Call Get_Val(newlayer_thick,sedfile,'NEWLAYER_THICK',line=linenum,echo=.true.) !read in number of cohesive sediment class Call Get_Val(ncs,sedfile,'NCS',line=linenum,echo=.true.) !read in number of non-cohesive sediment class Call Get_Val(nns,sedfile,'NNS',line=linenum,echo=.true.) !total number of classes for cohesive and non-cohesive sediment NST = NCS+NNS if(nst < 1)then write(*,*)'sediment input file must have at least one sediment class' write(*,*)'currently nst = ',nst stop endif NSED = NST ! to be consistent with original output code. allocate(sed(nst)) !J. Ge for tracer advection Allocate(sed0(nst)) Allocate(sed2(nst)) !J. Ge for tracer advection !read in switches for bedload calculation Call Get_Val(BEDLOAD,sedfile,'BEDLOAD',line=linenum,echo=.true.) !read in switches for suspended sediment calculation Call Get_Val(SUSLOAD,sedfile,'SUSLOAD',line=linenum,echo=.true.) !read in minumum sediment density Call Get_Val(min_Srho,sedfile,'MIN_SRHO',line=linenum,echo=.false.) !read in switches for 1-D sediment calculation Call Get_Val(oned_model,sedfile,'SED_ONED',line=linenum,echo=.true.) !read in switches for cohesive bed layer stratigraphy Call Get_Val(COHESIVE_BED,sedfile,'COHESIVE_BED',line=linenum,echo=.true.) if(COHESIVE_BED .and. NNS>0 )then write(*,*)'cohesive bed dynamics has been specified, however,the' write(*,*)'number of non-cohesive sediment classes are not-zero' write(*,*)'currently nns = ',nns stop end if !read in switches for non-cohesive bed layer stratigraphy Call Get_Val(NONCOHESIVE_BED2,sedfile,'NONCOHESIVE_BED2',line=linenum,echo=.true.) !read in switches for mixed-bed layer stratigraphy Call Get_Val(MIXED_BED,sedfile,'MIXED_BED',line=linenum,echo=.true.) !read in switches for morphology calculation Call Get_Val(SED_MORPH,sedfile,'SED_MORPH',line=linenum,echo=.true.) !read in switches for floc bed layer stratigraphy Call Get_Val(SED_FLOCS,sedfile,'SED_FLOCS',line=linenum,echo=.true.) !read in switches for sediment biomass due to vegation growth Call Get_Val(SED_BIOMASS,sedfile,'SED_BIOMASS',line=linenum,echo=.true.) !read in switches for sediment bed layer mixing (biodiffusion) Call Get_Val(SED_BIODIFF,sedfile,'SED_BIODIFF',line=linenum,echo=.true.) !read in switches for sediment defloccuation Call Get_Val(SED_DEFLOC,sedfile,'SED_DEFLOC',line=linenum,echo=.true.) if(SED_BIOMASS)then !read in switches for bottom seagrass influences Call Get_Val(SEAGRASS_BOTTOM ,sedfile,'SEAGRASS_BOTTOM ',line=linenum,echo=.true.) !read in switches for bottom seagrass sink Call Get_Val(SEAGRASS_SINK,sedfile,'SEAGRASS_SINK ',line=linenum,echo=.true.) endif if(MIXED_BED)then !read in cohesive transition Call Get_Val(transC,sedfile,'TRANSC',line=linenum,echo=.true.) !read in noncohesive transition Call Get_Val(transN,sedfile,'TRANSN',line=linenum,echo=.true.) endif if(BEDLOAD)then !read in switches for Meyer-Peter Mueller formulation Call Get_Val(BEDLOAD_MPM,sedfile,'BEDLOAD_MPM',line=linenum,echo=.true.) !read in switches for Soulsby formulation Call Get_Val(BEDLOAD_SOULSBY,sedfile,'BEDLOAD_SOULSBY',line=linenum,echo=.true.) !read in switches for memeth slope calculation Call Get_Val(SLOPE_NEMETH,sedfile,'SLOPE_NEMETH',line=linenum,echo=.true.) !read in switches for lesser slope calculation Call Get_Val(SLOPE_LESSER,sedfile,'SLOPE_LESSER',line=linenum,echo=.true.) endif !read in switch for constant tau cd Call Get_Val(sed_tau_cd_const,sedfile,'SED_TAU_CD_CONST',line=linenum,echo=.true.) !read in switch for linear tau cd Call Get_Val(sed_tau_cd_lin,sedfile,'SED_TAU_CD_LIN',line=linenum,echo=.true.) !read in Bed load transport rate coefficient Call Get_Val(bedload_coeff,sedfile,'BEDLOAD_COEFF',line=linenum,echo=.true.) !read in linear continuation for vertical settling Call Get_Val(linear_continuation,sedfile,'LINEAR_CONTINUATION',line=linenum,echo=.true.) !read in neumann continuation for vertical settling Call Get_Val(neumann,sedfile,'NEUMANN',line=linenum,echo=.true.) if(SED_BIODIFF)then !read in setting biodiffusivity profile Call Get_Val(DB_PROFILE,sedfile,'DB_PROFILE',line=linenum,echo=.true.) !read in Maximum biodiffusivity Call Get_Val(Dbmx,sedfile,'DBMAX',line=linenum,echo=.true.) !read in Minimum biodiffusivity Call Get_Val(Dbmm,sedfile,'DBMIN',line=linenum,echo=.true.) !read in Depth of maximum biodiffusivity Call Get_Val(Dbzs,sedfile,'DBZS',line=linenum,echo=.true.) !read in Depth of end of exponential biodiffusivity Call Get_Val(Dbzm,sedfile,'DBZM',line=linenum,echo=.true.) !read in Depth of minimum biodiffusivity Call Get_Val(Dbzp,sedfile,'DBZP',line=linenum,echo=.true.) endif if(ncs>0)then sed(1:ncs)%stype = 'cohesive' sed(ncs+1:nst)%stype = 'non-cohesive' allocate(ncstmp(ncs)) !read in sediment names for Cohesive sediment allocate(strtmp(ncs)) Call Get_Val_Array(strtmp,sedfile,'MUD_NAME',ncs,echo=.true.) sed(1:ncs)%sname=strtmp do i=1,ncs sed(i)%sname2=trim(strtmp(i))//'_bedload' end do deallocate(strtmp) !read in Median grain diameter (mm) for Cohesive sediment Call Get_Val_Array(ncstmp,sedfile,'MUD_SD50',ncs,echo=.true.) sed(1:ncs)%Sd50=ncstmp*0.001 !convert settling velocity to m/s from input mm/s !read in Sediment concentration (kg/m3). Call Get_Val_Array(ncstmp,sedfile,'MUD_CSED',ncs,echo=.true.) sed(1:ncs)%Csed_initial=ncstmp !read in Sediment grain density (kg/m3). Call Get_Val_Array(ncstmp,sedfile,'MUD_SRHO',ncs,echo=.true.) sed(1:ncs)%Srho=ncstmp !read in Particle settling velocity (mm/s). Call Get_Val_Array(ncstmp,sedfile,'MUD_WSED',ncs,echo=.true.) sed(1:ncs)%Wsed=ncstmp*0.001 !convert settling velocity to m/s from input mm/s !read in Surface erosion rate (kg/m2/s). Call Get_Val_Array(ncstmp,sedfile,'MUD_ERATE',ncs,echo=.true.) sed(1:ncs)%Erate=ncstmp !read in Critical shear for erosion and deposition (N/m2). Call Get_Val_Array(ncstmp,sedfile,'MUD_TAU_CE',ncs,echo=.true.) sed(1:ncs)%tau_ce=ncstmp Call Get_Val_Array(ncstmp,sedfile,'MUD_TAU_CD',ncs,echo=.true.) sed(1:ncs)%tau_cd=ncstmp !read in Porosity (nondimensional: 0.0-1.0): Vwater/(Vwater+Vsed). Call Get_Val_Array(ncstmp,sedfile,'MUD_POROS',ncs,echo=.true.) sed(1:ncs)%poros=ncstmp !read in Harmonic/biharmonic horizontal diffusion of tracer for nonlinear model Call Get_Val_Array(ncstmp,sedfile,'MUD_TNU2',ncs,echo=.true.) sed(1:ncs)%nl_tnu2=ncstmp Call Get_Val_Array(ncstmp,sedfile,'MUD_TNU4',ncs,echo=.true.) sed(1:ncs)%nl_tnu4=ncstmp ! read in Logical switches (TRUE/FALSE) to increase horizontal diffusivity ! of cohesive sediment trace in specific areas of the application ! domain (like sponge areas) allocate(nswitch(ncs)) Call Get_Val_Array(nswitch,sedfile,'MUD_Sponge',ncs,echo=.true.) sed(1:ncs)%LtracerSponge=nswitch !read in Vertical mixing coefficients for tracers in nonlinear model and ! basic state scale factor in adjoint-based algorithms. Call Get_Val_Array(ncstmp,sedfile,'MUD_AKT_BAK',ncs,echo=.true.) sed(1:ncs)%Akt_bak=ncstmp !read in Nudging/relaxation time scales, Call Get_Val_Array(ncstmp,sedfile,'MUD_TNUDG',ncs,echo=.true.) sed(1:ncs)%Tnudg=ncstmp !read in Morphological time scale factor (greater than or equal to 1.0) Call Get_Val_Array(ncstmp,sedfile,'MUD_MORPH_FAC',ncs,echo=.true.) sed(1:ncs)%morph_fac=ncstmp !read in Logical switches (TRUE/FALSE) to specify which variables to consider on ! tracers point Sources/Sinks (like river runoff). Call Get_Val_Array(nswitch,sedfile,'MUD_Ltracer',ncs,echo=.true.) sed(1:ncs)%LtracerSrc=nswitch !read in Logical switches (T/F) to process cohesive sediment tracer climatology. Call Get_Val_Array(nswitch,sedfile,'MUD_Ltclm',ncs,echo=.true.) sed(1:ncs)%LtracerCLM=nswitch !read in Logical switches (T/F) to activate nugding of cohesive !sediment tracer climatology. Call Get_Val_Array(nswitch,sedfile,'MUD_Tnudge',ncs,echo=.true.) sed(1:ncs)%LnudgeTCLM=nswitch deallocate(nswitch) endif if(COHESIVE_BED.or.MIXED_BED)then !read in minimum shear for erosion Call Get_Val(tcr_min,sedfile,'MUD_TAUCR_MIN',line=linenum,echo=.true.) !read in maximum shear for erosion Call Get_Val(tcr_max,sedfile,'MUD_TAUCR_MAX',line=linenum,echo=.true.) !read in Tau_crit profile slope Call Get_Val(tcr_slp,sedfile,'MUD_TAUCR_SLOPE',line=linenum,echo=.true.) !read in Tau_crit profile offset Call Get_Val(tcr_off,sedfile,'MUD_TAUCR_OFF',line=linenum,echo=.true.) !read in Tau_crit consolidation rate Call Get_Val(tcr_tim,sedfile,'MUD_TAUCR_TIME',line=linenum,echo=.true.) endif if(nns>0)then allocate(nswitch(nns)) allocate(nnstmp(nns)) !read in sediment names for non-Cohesive sediment allocate(strtmp(nns)) Call Get_Val_Array(strtmp,sedfile,'SAND_NAME',nns,echo=.true.) sed(ncs+1:ncs+nns)%sname=strtmp do i=1,nns sed(ncs+i)%sname2=trim(strtmp(i))//'_bedload' end do deallocate(strtmp) !read in Median grain diameter (mm) for non-cohesive sediment Call Get_Val_Array(nnstmp,sedfile,'SAND_SD50',nns,echo=.true.) sed(ncs+1:ncs+nns)%Sd50=nnstmp*.001 !convert mean grain diameter from mm to m !read in Sediment concentration (kg/m3) for non-cohesive sediment Call Get_Val_Array(nnstmp,sedfile,'SAND_CSED',nns,echo=.true.) sed(ncs+1:ncs+nns)%Csed_initial=nnstmp !read in grain density (kg/m3) for non-cohesive sediment Call Get_Val_Array(nnstmp,sedfile,'SAND_SRHO',nns,echo=.true.) sed(ncs+1:ncs+nns)%Srho=nnstmp !read in Particle settling velocity (mm/s) for non-cohesive sediment Call Get_Val_Array(nnstmp,sedfile,'SAND_WSED',nns,echo=.true.) sed(ncs+1:ncs+nns)%Wsed=nnstmp*.001 !convert settling velocity to m/s from input mm/s !read in Surface erosion rate (kg/m2/s) for non-cohesive sediment Call Get_Val_Array(nnstmp,sedfile,'SAND_ERATE',nns,echo=.true.) sed(ncs+1:ncs+nns)%Erate=nnstmp !read in Critical shear for erosion and deposition (N/m2). Call Get_Val_Array(nnstmp,sedfile,'SAND_TAU_CE',nns,echo=.true.) sed(ncs+1:ncs+nns)%tau_ce=nnstmp Call Get_Val_Array(nnstmp,sedfile,'SAND_TAU_CD',nns,echo=.true.) sed(ncs+1:ncs+nns)%tau_cd=nnstmp !read in Porosity (nondimensional: 0.0-1.0): Vwater/(Vwater+Vsed). Call Get_Val_Array(nnstmp,sedfile,'SAND_POROS',nns,echo=.true.) sed(ncs+1:ncs+nns)%poros=nnstmp !read in Harmonic/biharmonic horizontal diffusion of tracer Call Get_Val_Array(nnstmp,sedfile,'SAND_TNU2',nns,echo=.true.) sed(ncs+1:ncs+nns)%nl_tnu2=nnstmp Call Get_Val_Array(nnstmp,sedfile,'SAND_TNU4',nns,echo=.true.) sed(ncs+1:ncs+nns)%nl_tnu4=nnstmp !read in Logical switches (TRUE/FALSE) to increase horizontal diffusivity ! of noncohesive sediment trace Call Get_Val_Array(nswitch,sedfile,'SAND_Sponge',nns,echo=.true.) sed(ncs+1:ncs+nns)%LtracerSponge=nswitch !read in Vertical mixing coefficients for tracers in nonlinear model and ! basic state scale factor Call Get_Val_Array(nnstmp,sedfile,'SAND_AKT_BAK',nns,echo=.true.) sed(ncs+1:ncs+nns)%Akt_bak=nnstmp !read in Nudging/relaxation time scales, Call Get_Val_Array(nnstmp,sedfile,'SAND_TNUDG',nns,echo=.true.) sed(ncs+1:ncs+nns)%Tnudg=nnstmp !read in Morphological time scale factor (greater than or equal to 1.0) Call Get_Val_Array(nnstmp,sedfile,'SAND_MORPH_FAC',nns,echo=.true.) sed(ncs+1:ncs+nns)%morph_fac=nnstmp !read in Logical switches (TRUE/FALSE) to activate tracers point Sources/Sinks ! (like river runoff) and to specify which tracer variables to consider Call Get_Val_Array(nswitch,sedfile,'SAND_Ltsrc',nns,echo=.true.) sed(ncs+1:ncs+nns)%LtracerSrc=nswitch !read in Logical switches (TRUE/FALSE) to read and process ! noncohesive sediment tracers climatology Call Get_Val_Array(nswitch,sedfile,'SAND_Ltclm',nns,echo=.true.) sed(ncs+1:ncs+nns)%LtracerCLM=nswitch !read in Logical switches (TRUE/FALSE) to nudge the desired noncohesive sediment ! tracer climatology field. Call Get_Val_Array(nswitch,sedfile,'SAND_Tnudge',nns,echo=.true.) sed(ncs+1:ncs+nns)%LnudgeTCLM=nswitch endif !------------------------------------------------------------------------------ ! Flocculation Sediment Parameters, [1] values expected. !------------------------------------------------------------------------------ if(SED_FLOCS)then !read in l_ADS : boolean set to .true. if differential settling aggregation Call Get_Val(l_ADS,sedfile,'l_ADS',line=linenum,echo=.true.) !read in l_ASH : boolean set to .true. if shear aggregation Call Get_Val(l_ASH,sedfile,'l_ASH',line=linenum,echo=.true.) !read in l_COLLFRAG : boolean set to .true. if collision-induced fragmentation enable Call Get_Val(l_COLLFRAG,sedfile,'l_COLLFRAG',line=linenum,echo=.true.) !read in f_dp0 : Primary particle size (m), typically 4e-6 m Call Get_Val(f_dp0,sedfile,'f_dp0',line=linenum,echo=.true.) !read in f_nf : Floc fractal dimension, typically ranging from 1.6 to 2.6 Call Get_Val(f_nf,sedfile,'f_nf',line=linenum,echo=.true.) !read in f_dmax : (unused) Maximum diameter (m) Call Get_Val(f_dmax,sedfile,'f_dmax',line=linenum,echo=.true.) !read in f_nb_frag : number of fragments by shear erosion. If binary/ternary : 2. Call Get_Val(f_nb_frag,sedfile,'f_nb_frag',line=linenum,echo=.true.) !read in f_alpha : flocculation efficiency, ranging from 0 .to 1. Call Get_Val(f_alpha,sedfile,'f_alpha',line=linenum,echo=.true.) !read in f_beta : shear fragmentation rate Call Get_Val(f_beta,sedfile,'f_beta',line=linenum,echo=.true.) !read in f_ater : for ternary breakup, use 0.5, for binary : 0 Call Get_Val(f_ater,sedfile,'f_ater',line=linenum,echo=.true.) !read in f_ero_frac : fraction of the shear fragmentation term !transfered to shear erosion. Ranging from 0. (no erosion) to 1. !(all erosion) Call Get_Val(f_ero_frac,sedfile,'f_ero_frac',line=linenum,echo=.true.) !read in f_ero_nbfrag : Number of fragments induced by shear erosion. Call Get_Val(f_ero_nbfrag,sedfile,'f_ero_nbfrag',line=linenum,echo=.true.) !read in f_ero_iv : fragment size class (could be changed to a ! particle size or a particle distribution (INTEGER) Call Get_Val(f_ero_iv,sedfile,'f_ero_iv',line=linenum,echo=.true.) !read in f_collfragparam : Fragmentation rate for collision !-induced breakup Call Get_Val(f_collfragparam,sedfile,'f_collfragparam',line=linenum,echo=.true.) !read in f_clim : min concentration below which flocculation ! processes are not calculated Call Get_Val(f_clim,sedfile,'f_clim',line=linenum,echo=.true.) !read in l_testcase - if .TRUE. sets G(t) to values from Verney ! et al., 2011 lab experiment Call Get_Val(l_testcase,sedfile,'l_testcase',line=linenum,echo=.true.) !read in gls_p Call Get_Val(gls_p,sedfile,'GLS_P',line=linenum,echo=.true.) !read in gls_m Call Get_Val(gls_m,sedfile,'GLS_M',line=linenum,echo=.true.) !read in gls_n Call Get_Val(gls_n,sedfile,'GLS_N',line=linenum,echo=.true.) !read in gls_cmu0 Call Get_Val(gls_cmu0,sedfile,'GLS_CMU0',line=linenum,echo=.true.) !read in FLOC_TURB_DISS Call Get_Val(FLOC_TURB_DISS,sedfile,'FLOC_TURB_DISS',line=linenum,echo=.true.) !read in FLOC_BBL_DISS Call Get_Val(FLOC_BBL_DISS,sedfile,'FLOC_BBL_DISS',line=linenum,echo=.true.) endif if(SED_FLOCS.and.SED_DEFLOC)then allocate(mud_frac_eq(ncs)) !read in Equilibrium fractional class distribution (they should add ! up to 1). Call Get_Val_Array(mud_frac_eq,sedfile,'MUD_FRAC_EQ',ncs,echo=.true.) mud_frac_eq(1:ncs)=mud_frac_eq !read in Time scale of flocculation descomposition in bed Call Get_Val(t_dfloc,sedfile,'t_dfloc',line=linenum,echo=.true.) endif if(SED_BIOMASS)then !read in Number of hours for depth integration Call Get_Val(nTbiom,sedfile,'nTbiom',line=linenum,echo=.true.) !read in Call Get_Val(sgr_diam,sedfile,'SGR_DIAM',line=linenum,echo=.true.) !read in Call Get_Val(sgr_density,sedfile,'SGR_DENS',line=linenum,echo=.true.) !read in Call Get_Val(sgr_Hthres,sedfile,'SGR_HTHRES',line=linenum,echo=.true.) endif ! !----------------------------------------------------------------------- ! Report input parameters. !----------------------------------------------------------------------- ! !read in start interval for sed model Call Get_Val(sed_start,sedfile,'SED_START',line=linenum,echo=.true.) !hot start option Call Get_Val(sed_hot_start,sedfile,'SED_HOT_START',line=linenum,echo=.true.) !read in morphology parameters Call Get_Val(morpho_model,sedfile,'MORPHO_MODEL',line=linenum,echo=.true.) Call Get_Val(morpho_factor,sedfile,'MORPHO_FACTOR',line=linenum,echo=.true.) if(morpho_model)then Call Get_Val(morpho_incr,sedfile,'MORPHO_INCR',line=linenum,echo=.true.) Call Get_Val(morpho_strt,sedfile,'MORPHO_STRT',line=linenum,echo=.true.) else morpho_incr = 1 endif !read in number of bed layers Call Get_Val(nbed,sedfile,'NBED',line=linenum,echo=.true.) if(nbed < 1)then write(*,*)'sediment input file must have at least one bed layer' write(*,*)'currently nbed = ',nbed stop endif !read in logical for infinite sediment supply Call Get_Val(inf_bed,sedfile,'INF_BED',line=linenum,echo=.true.) !read in initial bed porosity allocate(init_bed_porosity(nbed)) if(nbed==1)then Call Get_Val(ftemp,sedfile,'INIT_BED_POROSITY',line=linenum,echo=.true.) init_bed_porosity(1) = ftemp else open(unit=999,file=trim(sedfile),form='formatted') iscan = scan_file(999,"INIT_BED_POROSITY",fvec = init_bed_porosity,nsze = nbed) if(iscan /= 0)then write(*,*)'problem reading INIT_BED_POROSITY from sediment param file' stop endif close(999) endif if(minval(init_bed_porosity) < 0. .or. maxval(init_bed_porosity) > 1)then write(*,*)'error in init_bed_porosity in sed param file' write(*,*)'you entered: ',init_bed_porosity write(*,*)'must be in range [0,1]' stop endif !read in initial bed thickness allocate(init_bed_thickness(nbed)) if(nbed==1)then Call Get_Val(ftemp,sedfile,'INIT_BED_THICKNESS',line=linenum,echo=.true.) init_bed_thickness(1) = ftemp else open(unit=999,file=trim(sedfile),form='formatted') iscan = scan_file(999,"INIT_BED_THICKNESS",fvec = init_bed_thickness,nsze = nbed) if(iscan /= 0)then write(*,*)'problem reading INIT_BED_THICKNESS from sediment param file' stop endif close(999) endif if(minval(init_bed_thickness) <= 0.)then write(*,*)'error in init_bed_thickness in sed param file' write(*,*)'you entered: ',init_bed_thickness write(*,*)'must be > 0' endif !read in initial bed age allocate(init_bed_age(nbed)) if(nbed==1)then Call Get_Val(ftemp,sedfile,'INIT_BED_AGE',line=linenum,echo=.true.) init_bed_age(1) = ftemp else open(unit=999,file=trim(sedfile),form='formatted') iscan = scan_file(999,"INIT_BED_AGE",fvec = init_bed_age,nsze = nbed) if(iscan /= 0)then write(*,*)'problem reading INIT_BED_AGE from sediment param file' stop endif close(999) endif if(minval(init_bed_age) <= 0.)then write(*,*)'error in init_bed_age in sed param file' write(*,*)'you entered: ',init_bed_age write(*,*)'must be > 0' endif !read in initial bed bio-diffusivity if(SED_BIODIFF)then allocate(init_bed_biodiff(nbed)) if(nbed==1)then Call Get_Val(ftemp,sedfile,'INIT_BED_BIODIFF',line=linenum,echo=.true.) init_bed_biodiff(1) = ftemp else open(unit=999,file=trim(sedfile),form='formatted') iscan = scan_file(999,"INIT_BED_BIODIFF",fvec = init_bed_biodiff,nsze = nbed) if(iscan /= 0)then write(*,*)'problem reading INIT_BED_BIODIFF from sediment param file' stop endif close(999) endif if(minval(init_bed_age) <= 0.)then write(*,*)'error in init_bed_bioddiff in sed param file' write(*,*)'you entered: ',init_bed_biodiff write(*,*)'must be > 0' endif end if !read Sediment critical stress for erosion allocate(init_bed_tau_crit(nbed)) if(nbed==1)then Call Get_Val(ftemp,sedfile,'INIT_BED_TAU_CRIT',line=linenum,echo=.true.) init_bed_tau_crit(1) = ftemp else open(unit=999,file=trim(sedfile),form='formatted') iscan = scan_file(999,"INIT_BED_TAU_CRIT",fvec = init_bed_tau_crit,nsze = nbed) if(iscan /= 0)then write(*,*)'problem reading INIT_BED_TAU_CRIT from sediment param file' stop endif close(999) endif if(minval(init_bed_tau_crit) <= 0.)then write(*,*)'error in init_bed_tau_crit in sed param file' write(*,*)'you entered: ',init_bed_tau_crit write(*,*)'must be > 0' endif !read in initial bed fraction allocate(init_bed_fraction(nst,nbed)) allocate(nsttmp(nst*nbed)) Call Get_Val_Array(nsttmp,sedfile,'INIT_BED_FRACTION',nst*nbed,echo=.true.) do i=1,nbed init_bed_fraction(1:nst,i)=nsttmp((i-1)*nst+1:(i-1)*nst+nst) if(sum(init_bed_fraction(1:nst,i))/=1.0)then write(*,*)'error in init_bed_fraction in sed param file' write(*,*)'also error in your chosen career' write(*,*)'in bed layer:',i write(*,*)'you entered: ',init_bed_fraction(1:nst,i) write(*,*)'must have the summary = 1' stop end if end do if(minval(init_bed_fraction) < 0. .or. maxval(init_bed_fraction) > 1.)then write(*,*)'error in init_bed_fraction in sed param file' write(*,*)'also error in your chosen career' write(*,*)'you entered: ',init_bed_fraction write(*,*)'must be >= 0 and <= 1' stop endif deallocate(nsttmp) !read in nudging switch Call Get_Val(sed_nudge,sedfile,'SED_NUDGE',line=linenum,echo=.true.) !read in nudging relaxation factor if(sed_nudge) then Call Get_Val(sed_alpha,sedfile,'SED_ALPHA',line=linenum,echo=.false.) Call Get_Val(sed_ramp ,sedfile,'SED_RAMP ',line=linenum,echo=.false.) endif !read in params controlling output Call Get_Val(sed_dumpbed,sedfile,'SED_DUMPBED',line=linenum,echo=.true.) Call Get_Val(sed_dumpbot,sedfile,'SED_DUMPBOT',line=linenum,echo=.true.) !hindered effect on vertical diffusion Call Get_Val(VERT_HINDERED,sedfile,'VERT_HINDERED',line=linenum,echo=.true.) !read in point source switch Call Get_Val(sed_source,sedfile,'SED_PTSOURCE',line=linenum,echo=.true.) !check values if(nst < 1)then write(*,*)'sediment input file must have at least one sediment class' write(*,*)'currently nst = ',nst stop endif if(nbed < 1)then write(*,*)'sediment input file must have at least one bed layer' write(*,*)'currently nbed = ',nbed stop endif if(morpho_factor < 0. )then write(*,*)'sediment morpho_factor must be positive' write(*,*)'currently is = ',morpho_factor stop endif if(morpho_incr < 0 )then write(*,*)'sediment morpho_incr must be positive' write(*,*)'currently is = ',morpho_incr stop endif if(morpho_strt < 0 )then write(*,*)'sediment morpho_strt must be non-negative' write(*,*)'currently is = ',morpho_strt stop endif if(dbg_set(dbg_sbr)) write(ipt,*) "End: Read_Sed_Params " End Subroutine Read_Sed_Params !======================================================================= ! Allocate Sediment Variables !======================================================================= Subroutine Alloc_Sed_Vars use lims, only: nt,mt,kb,numqbc,kbm1 ! use mod_obcs, only: iobcn implicit none integer i,ic,tmp1,tmp2,tmp3 if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Alloc_Sed_Vars " if(COHESIVE_BED.or.SED_BIODIFF.or. MIXED_BED)then MBEDP = 5 ! Bed properties else MBEDP = 4 ! Bed properties endif allocate(idSbed(MBEDP)) if(SEAGRASS_BOTTOM .and. MIXED_BED)then MBOTP = 34 ! Bottom properties elseif(SEAGRASS_BOTTOM .and. .not.(MIXED_BED))then MBOTP = 25 ! Bottom properties elseif(COHESIVE_BED.or.SED_BIODIFF)then MBOTP = 31 ! Bottom properties elseif(MIXED_BED)then MBOTP = 32 ! Bottom properties else MBOTP = 23 ! Bottom properties endif allocate(idBott(MBOTP)) if(SED_BIOMASS.and.SEAGRASS_BOTTOM.and.MIXED_BED)then isgrD = 33 ! seagrass shoot density isgrH = 34 ! seagrass height elseif(SED_BIOMASS.and.SEAGRASS_BOTTOM.and.(.not.(MIXED_BED)))then isgrD = 24 ! seagrass shoot density isgrH = 25 ! seagrass height endif n_bot_vars => MBOTP n_bed_chars => MBEDP ! !----------------------------------------------------------------------- ! Initialize tracer identification indices. !----------------------------------------------------------------------- ! ! Allocate various nested grid depended parameters ! if(COHESIVE_BED.or.MIXED_BED)then tcr_min = 0.0_sp tcr_max = 0.0_sp tcr_slp = 0.0_sp tcr_off = 0.0_sp tcr_tim = 0.0_sp endif if(MIXED_BED)then transC = 0.0_sp transN = 0.0_sp endif if(SED_BIOMASS)then sgr_diam = 0.0_sp sgr_density = 0.0_sp sgr_Hthres = 0.0_sp endif ! ! Allocate sediment tracers indices vectors. ! allocate ( idsed(MAX(1,NST)) ) allocate ( idmud(MAX(1,NCS)) ) allocate ( isand(MAX(1,NNS)) ) allocate ( idBmas(NST) ) allocate ( idfrac(NST) ) allocate ( idUbld(NST) ) allocate ( idVbld(NST) ) ! ! Set cohesive and noncohesive suspended sediment tracers ! identification indices. ! ic=0 DO i=1,NCS ic=ic+1 idmud(i)=ic idsed(i)=idmud(i) END DO DO i=1,NNS ic=ic+1 isand(i)=ic idsed(NCS+i)=isand(i) END DO if(dbg_set(dbg_sbr)) write(ipt,*) "End: Alloc_Sed_Vars " End Subroutine Alloc_Sed_Vars SUBROUTINE Alloc_sedbed Use Scalar Use Control, only : ireport,iint,msr,par Use Lims, only : m,mt,nt,kbm1,numqbc,kb,nprocs,myid # if defined (MULTIPROCESSOR) Use Mod_Par # endif implicit none integer i,tmp1,tmp2,tmp3 ! !----------------------------------------------------------------------- ! Allocate structure variables. !----------------------------------------------------------------------- ! if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Alloc_sedbed" ! !ensure arrays have nonzero dimension tmp1 = max(numqbc,1) tmp2 = max(Nobc+1,1) tmp3 = max(Nobc ,1) if(BEDLOAD.and.SED_DUMPBED)then allocate ( SEDBED% avgbedldu(0:mt,NST) ) ;SEDBED% avgbedldu = 0.0_sp allocate ( SEDBED% avgbedldv(0:mt,NST) ) ;SEDBED% avgbedldv = 0.0_sp endif allocate ( SEDBED % bed(0:mt,Nbed,MBEDP) ) ;SEDBED % bed = 0.0_sp allocate ( SEDBED % bed_frac(0:mt,Nbed,NST) ) ;SEDBED % bed_frac = 0.0_sp allocate ( SEDBED % bed_mass(0:mt,Nbed,2,NST) ) ;SEDBED % bed_mass = 0.0_sp if(SED_MORPH)then allocate ( SEDBED % bed_thick0(0:mt) ) ;SEDBED % bed_thick0 = 0.0_sp allocate ( SEDBED % bed_thick(0:mt,1:2) ) ;SEDBED % bed_thick = 0.0_sp endif if(BEDLOAD)then allocate ( SEDBED % bedldu(0:mt,NST) ) ;SEDBED % bedldu = 0.0_sp allocate ( SEDBED % bedldv(0:mt,NST) ) ;SEDBED % bedldv = 0.0_sp end if allocate ( SEDBED % bottom(0:mt,MBOTP) ) ;SEDBED % bottom = 0.0_sp bottom => SEDBED % bottom bed => SEDBED % bed do i=1,nst allocate(sed(i)%conc(0:mt,kb )) ; sed(i)%conc = 0.0_sp allocate(sed(i)%cnew(0:mt,kb )) ; sed(i)%cnew = 0.0_sp allocate(sed(i)%mass(0:mt,nbed)) ; sed(i)%mass = 0.0_sp allocate(sed(i)%frac(0:mt,nbed)) ; sed(i)%frac = 0.0_sp allocate(sed(i)%bflx(0:mt )) ; sed(i)%bflx = 0.0_sp allocate(sed(i)%eflx(0:mt )) ; sed(i)%eflx = 0.0_sp allocate(sed(i)%dflx(0:mt )) ; sed(i)%dflx = 0.0_sp allocate(sed(i)%cdis(500 )) ; sed(i)%cdis = 0.0_sp !JQI NOV2021 allocate(sed(i)%cflx(tmp3,kbm1)) ; sed(i)%cflx = 0.0_sp ! KURT GLAESEMANN - remove + 1 from dimension !JQI NOV2021 allocate(sed(i)%cobc(tmp3 )) ; sed(i)%cobc = 0.0_sp allocate(sed(i)%depm(0:mt )) ; sed(i)%depm = 0.0_sp allocate(sed(i)%wset0(0:mt,kb )) ; sed(i)%wset0 = 0.0_sp sed(i)%wset0 = sed(i)%Wsed allocate(sed(i)%t_cd(0:mt )) ; sed(i)%t_cd = 0.0_sp allocate(sed(i)%t_ce(0:mt )) ; sed(i)%t_ce = 0.0_sp allocate(sed(i)%rate(0:mt )) ; sed(i)%rate = 0.0_sp !J. Ge for tracer advection allocate(sed0(i)%conc(0:mt,kb )) ; sed0(i)%conc = 0.0 allocate(sed2(i)%conc(0:mt,kb )) ; sed2(i)%conc = 0.0 !J. Ge for tracer advection end do !allocate suspended sediment arrays allocate(csed(0:mt,kb)) ; csed = 0.0 !allocate bottom shear stress array allocate(taub(0:mt)) ; taub = 0.0_sp if(SUSLOAD)then allocate ( SEDBED % ero_flux(0:mt,NST) ) ;SEDBED % ero_flux = 0.0_sp allocate ( SEDBED % settling_flux(0:mt,NST) ) ;SEDBED % settling_flux = 0.0_sp endif if(SED_BIOMASS)then allocate ( SEDBED % Dstp_max(0:mt,24) ) ;SEDBED % Dstp_max = 0.1_sp if(SEAGRASS_SINK)then allocate (SgrN (0:mt)) ;SgrN = 0.0_sp end if endif if(dbg_set(dbg_sbr)) write(ipt,*) "End: Alloc_sedbed" RETURN END SUBROUTINE Alloc_sedbed !========================================================================== ! Update Bottom Properties due to changes in surface layer sed fractions ! Use geometric mean, good for log normal distribution !========================================================================== Subroutine Update_Bottom_Properties use lims, only: m implicit none integer :: i,ised real(sp) :: temp,sum1,sum2,sum3,sum4,eps if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Update_Bottom_Properties " ! initialize eps = epsilon(temp) ! update bottom properties using new bed fractions do i=1,m sum1 = 1.0 ; sum2 = 1.0 ; sum3 = 1.0 ; sum4 = 1.0 do ised=1,nst if(SEDIMENT_PARAMETER_TYPE/=UNIFORM)then sum1 = sum1*(sed(ised)%t_ce(i))**sedbed%bed_frac(i,1,ised) else sum1 = sum1*(sed(ised)%tau_ce)**sedbed%bed_frac(i,1,ised) end if sum2 = sum2*(sed(ised)%Sd50 )**sedbed%bed_frac(i,1,ised) sum3 = sum3*(sed(ised)%wsed )**sedbed%bed_frac(i,1,ised) sum4 = sum4*(sed(ised)%Srho )**sedbed%bed_frac(i,1,ised) end do bottom(i,itauc) = sum1 bottom(i,isd50) = sum2 bottom(i,iwsed) = sum3 bottom(i,idens) = max(sum4,min_Srho) end do if(dbg_set(dbg_sbr)) write(ipt,*) "End: Update_Bottom_Properties " End Subroutine Update_Bottom_Properties !========================================================================== ! Update Active Layer Thickness ! Use method of Harris and Wiberg (1997): ! -Approaches to quantifying long-term continental shelf sediment transport ! -with an example from the northern California STRESS mid-shelf site. ! [Continental Shelf REsearch, 17, 1389-1418] ! ! z_a = max[k_1*(tau_bot - tau_ce)*rho_0 , 0] + k_2 * D_50 ! where ! k_1 = .007 (empirical constant) ! k_2 = 6.0 (empirical constant) ! tau_bot = shear stress on the bed ! tau_ce = critical shear stress for erosion of the bed (averaged over ! sediment classes) ! D_50 = median grain diameter of the surface sediment ! note: our shear stresses are nondimensionalized by rho, so ! rho does not appear in actual calculation (below) !========================================================================== Subroutine Calc_Active_Layer use lims, only: m implicit none integer :: i real(sp) :: mfac mfac = dim(morpho_factor,1.0) !calculate active layer thickness do i=1,m bottom(i,iactv) = mfac*max(0.,.007*(taub(i)-bottom(i,itauc)))+6.0*bottom(i,isd50) end do End Subroutine Calc_Active_Layer !========================================================================== ! Report Sediment Setup To Screen !========================================================================== Subroutine Report_Sed_Setup use control, only: msr implicit none integer :: i if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Report_Sed_Setup " ! echo sediment model parameters to screen if(msr)then write(*,*) write(*,*)'------------------ Sediment Model Setup --------------------' write(*,*)'! nbed :',nbed write(*,*)'! nst :',nst write(*,*)'! n_report :',n_report if(sed_start > 0)then write(*,*)'! sed_start :',sed_start else write(*,*)'! sed_start : no delay' endif if(sed_hot_start)then write(*,*)'! sed_hot_start :',sed_hot_start endif write(*,*)'! min_Srho :',min_Srho if(sed_nudge)then write(*,*)'! open boundary nudging : active' write(*,*)'! nudging relax factor :',sed_alpha write(*,*)'! # nudging ramp its :',sed_ramp else write(*,*)'! open boundary nudging : not active' endif if(sed_source)then write(*,*)'! point source forcing : active' else write(*,*)'! point source forcing : not active' endif if(morpho_model)then write(*,*)'! morphodynamics : active' write(*,*)'! morpho_incr :',morpho_incr write(*,*)'! morpho_strt :',morpho_strt else write(*,*)'! morphodynamics : not active' endif write(*,*)'! morpho_factor :',morpho_factor if(bedload)then write(*,*)'! bedload dynamics : active' write(*,*)'! Critical Shields :',Shield_Cr_MPM write(*,*)'! Bedload Exponent :',Gamma_MPM write(*,*)'! Bedload Constant :',k_MPM write(*,*)'! Bedload Rate Coeff. :',bedload_rate if(bedload_smooth)then write(*,*)'! Bedload Smoother : active' endif else write(*,*)'! bedload dynamics : not active' endif if(susload)then write(*,*)'! susload dynamics : active' else write(*,*)'! susload dynamics : not active' endif if(inf_bed)then write(*,*)'! sediment supply : infinite' else write(*,*)'! sediment supply : finite' endif write(*,*)'! Settling CFL :',settle_cfl write(*,*)'!' write(*,*)'! class type Sd50 Wset tau_ce '//& 'tau_cd Erate Spor Srho' do i=1,nst write(*,'(1X,A1,1X,A12,1X,A6,6(F7.4,1X),F8.2)')'!', & trim(sed(i)%sname),trim(sed(i)%stype(1:6)), & sed(i)%Sd50,sed(i)%Wsed, & sed(i)%tau_ce,sed(i)%tau_cd,sed(i)%erate,sed(i)%poros, & sed(i)%Srho end do end if ! echo sediment model initial conditions to screen call sed_report write(*,*)'------------------------------------------------------------' write(*,*) if(dbg_set(dbg_sbr)) write(ipt,*) "End: Report_Sed_Setup " End Subroutine Report_Sed_Setup !========================================================================== ! Monitor Sediment Model: Compute and Write Statistics to Screen !========================================================================== Subroutine Sed_Report use control, only: msr ,par, mpi_fvcom_group use lims , only: m,mt,kbm1, msrid implicit none real(dp) :: tau_min,tau_max,tmp1(3),tmp2(3),dthck_max integer :: i,dim1,dim2,ierr # if defined (FLUID_MUD) real(dp) :: temp1(5),temp2(5),mud_max,ua_max,va_max,ua_min,va_min REAL(DP), DIMENSION(3) :: SBUF,RBUF1,RBUF2,RBUF3 real(dp) :: ESTOT,NSTOT # endif if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Sed_Report " !set up limits dim1 = mt dim2 = kbm1 !calculate max/min/rms concentrations do i=1,nst sed(i)%cmax = maxval(sed(i)%conc(1:dim1,1:dim2)) sed(i)%cmin = minval(sed(i)%conc(1:dim1,1:dim2)) sed(i)%crms = sum(dble(abs(sed(i)%conc(1:dim1,1:dim2)))) sed(i)%crms = sed(i)%crms/(dim1*dim2) end do dthck_max = maxval(abs(bottom(1:m,dthck))) tau_max = maxval(taub(1:m)) tau_min = minval(taub(1:m)) # if defined (FLUID_MUD) ESTOT = DBLE(NGL) NSTOT = DBLE(MGL) !calculate/max/min/rms mud layer mud_max = maxval(abs(D_FMLcm)) ua_max = maxval(ua_fml) ua_min = minval(ua_fml) va_max = maxval(va_fml) va_min = minval(va_fml) SBUF(1) = SUM(DBLE(UA_fml(1:n))) SBUF(2) = SUM(DBLE(VA_fml(1:n))) SBUF(3) = SUM(DBLE(D_FML(1:m))) RBUF1 = SBUF # endif # if defined (MULTIPROCESSOR) IF(PAR)THEN do i=1,nst tmp1(1:3) = (/sed(i)%cmax,sed(i)%cmin,sed(i)%crms/) CALL MPI_REDUCE(tmp1,tmp2,3,MPI_DP,MPI_MAX,MSRID-1,MPI_FVCOM_GROUP,IERR) if(msr)then sed(i)%cmax = tmp2(1) sed(i)%cmin = tmp2(2) sed(i)%crms = tmp2(3) endif end do tmp1(1:3) = (/dthck_max,tau_min,tau_max/) CALL MPI_REDUCE(tmp1(1),tmp2(1),1,MPI_DP,MPI_MAX,MSRID-1,MPI_FVCOM_GROUP,IERR) CALL MPI_REDUCE(tmp1(2),tmp2(2),1,MPI_DP,MPI_MIN,MSRID-1,MPI_FVCOM_GROUP,IERR) CALL MPI_REDUCE(tmp1(3),tmp2(3),1,MPI_DP,MPI_MAX,MSRID-1,MPI_FVCOM_GROUP,IERR) dthck_max = tmp2(1) tau_min = tmp2(2) tau_max = tmp2(3) # if defined (FLUID_MUD) temp1(1:5) = (/mud_max,ua_max,ua_min,va_max,va_min/) CALL MPI_REDUCE(temp1(1),temp2(1),1,MPI_DP,MPI_MAX,MSRID-1,MPI_FVCOM_GROUP,IERR) CALL MPI_REDUCE(temp1(2),temp2(2),1,MPI_DP,MPI_MAX,MSRID-1,MPI_FVCOM_GROUP,IERR) CALL MPI_REDUCE(temp1(3),temp2(3),1,MPI_DP,MPI_MIN,MSRID-1,MPI_FVCOM_GROUP,IERR) CALL MPI_REDUCE(temp1(4),temp2(4),1,MPI_DP,MPI_MAX,MSRID-1,MPI_FVCOM_GROUP,IERR) CALL MPI_REDUCE(temp1(5),temp2(5),1,MPI_DP,MPI_MIN,MSRID-1,MPI_FVCOM_GROUP,IERR) mud_max = temp2(1) ua_max = temp2(2) ua_min = temp2(3) va_max = temp2(4) va_min = temp2(5) CALL MPI_REDUCE(SBUF,RBUF1,3,MPI_DP,MPI_SUM,MSRID-1,MPI_FVCOM_GROUP,IERR) IF(ISNAN(RBUF1(1)).or.ISNAN(RBUF1(2)).or.ISNAN(RBUF1(3)))THEN CALL FATAL_ERROR("fluid mud caclulation fails with unstable iteration","Please check the model") END IF # endif ENDIF # endif ! write statistics to screen if(.not.msr)return write(*,* )'====================Sediment Model Stats======================' write(*,* )'! quantity : avg(g/L) max(g/L) ' do i=1,nst write(*,100)'!',trim(sed(i)%sname),' :',sed(i)%crms,sed(i)%cmax end do write(*,* )'! max sed thick change : ',dthck_max write(*,* )'! max bottom stress : ',tau_max write(*,* )'! min bottom stress : ',tau_min # if defined (FLUID_MUD) write(*,* )'! max mud thickness : ',mud_max write(*,* )'! max u-speed of mud : ',ua_max write(*,* )'! min u-speed of mud : ',ua_min write(*,* )'! max v-speed of mud : ',va_max write(*,* )'! min v-speed of mud : ',va_min write(*,* )'! avg u-speed of mud : ',RBUF1(1)/ESTOT write(*,* )'! avg v-speed of mud : ',RBUF1(2)/ESTOT write(*,* )'! avg thickness of mud : ',RBUF1(3)/NSTOT # endif 100 format(1x,a1,a20,a5,3f12.6) 101 format(1x,a26,f12.6) if(dbg_set(dbg_sbr)) write(ipt,*) "End: Sed_Report " End Subroutine Sed_Report !========================================================================== ! Setup open boundary nudging for sediment obc ! - check for existence of sed nudging file ! - read data (concentration of each sed type at each obc) !========================================================================== Subroutine Setup_Sed_OBC Use Control, only : casename,input_dir,msr,serial,par ! Use Mod_OBCs, only : iobcn,i_obc_gl,iobcn_gl Use Mod_Utils, only: pstop # if defined (MULTIPROCESSOR) Use Mod_Par, only : nlid # endif if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Setup_Sed_OBC " if(dbg_set(dbg_sbr)) write(ipt,*) "End: Setup_Sed_OBC " End Subroutine Setup_Sed_OBC !======================================================================= ! ! ! This routine computes sediment bedload transport using the Meyer- ! ! Peter and Muller (1948) formulation for unidirectional flow and ! ! Souksby and Damgaard (2005) algorithm that accounts for combined ! ! effect of currents and waves. ! ! ! ! References: ! ! ! ! Meyer-Peter, E. and R. Muller, 1948: Formulas for bedload transport ! ! In: Report on the 2nd Meeting International Association Hydraulic ! ! Research, Stockholm, Sweden, pp 39-64. ! ! ! ! Soulsby, R.L. and J.S. Damgaard, 2005: Bedload sediment transport ! ! in coastal waters, Coastal Engineering, 52 (8), 673-689. ! ! ! ! Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G. ! ! Arango, 2008: Development of a three-dimensional, regional, ! ! coupled wave, current, and sediment-transport model, Computers ! ! & Geosciences, 34, 1284-1306. ! ! ! !======================================================================= Subroutine calc_sed_bedload use control, only: msr ,par, mpi_fvcom_group use lims , only: m,mt,kbm1, msrid use all_vars # if defined (MULTIPROCESSOR) Use Mod_Par # endif implicit none ! ! Local variable declarations. ! integer :: i, ised, j, k,cnt,jj,cc,ibed integer :: n1,n2,n3 integer :: ia,ib,j1,j2 integer, allocatable ::IAB_LIMIT(:) real(sp), parameter :: eps = 1.0E-14_sp real(sp) :: cff1, cff2, cff3, cff4, cff5 real(sp) :: Dstp, bed_change real(sp), parameter :: one = 1.0_sp ! for BEDLOAD_MPM real(sp), dimension(0:mt) :: tau_w real(sp) :: sed_angle,max_slope real(sp) :: bedld, bedld_mass,poro real(sp) :: smgd, smgdr, osmgd, Umag real(sp) :: rhs_bed, Ra, phi, Clim,sed_slope real(sp) :: tau_x,tau_y,flowdir,bldx,bldy real(sp), allocatable :: cff(:) real(sp), dimension(0:mt) :: bld_dr real(sp), dimension(0:mt) :: dzdx real(sp), dimension(0:mt) :: dzdy real(sp), dimension(0:mt) :: a_slopex real(sp), dimension(0:mt) :: a_slopey real(sp), dimension(:),allocatable :: bflxe ! bedload scalar at cell center real(sp), dimension(:),allocatable :: bflxn ! bedload scalar at node point real(sp), dimension(0:nt) :: bustrcwmax_e real(sp), dimension(0:nt) :: bvstrcwmax_e real(sp), dimension(:), allocatable :: btx !x-component of bedload transport real(sp), dimension(:), allocatable :: bty !y-component of bedload transport real(sp), dimension(:), allocatable :: btxn !x-component of bedload transport real(sp), dimension(:), allocatable :: btyn !y-component of bedload transport real(sp), dimension(0:nt) :: dbdx_e !bed slope at cell face real(sp), dimension(0:nt) :: dbdy_e !bed slope at cell face real(sp), dimension(0:mt) :: dbdx_n !bed slope at node real(sp), dimension(0:mt) :: dbdy_n !bed slope at node ! for BEDLOAD_MPM real(sp), dimension(0:mt) :: angleu real(sp), dimension(0:mt) :: anglev ! for BEDLOAD_SOULSBY real(sp) :: theta_mean, theta_wav, w_asym real(sp) :: theta_max, theta_max1, theta_max2 real(sp) :: phi_x1, phi_x2, phi_x, phi_y real(sp) :: bedld_x, bedld_y, tau_cur, waven, wavec real(sp) :: dbdx,dbdy,beta,slopex,slopey real(sp), dimension(0:mt) :: phic real(sp), dimension(0:mt) :: phicw real(sp), dimension(0:mt) :: tau_wav real(sp), dimension(0:mt) :: tau_mean real(sp), parameter :: kdmax = 100.0_sp real(sp),allocatable,dimension(:,:)::temp real(sp),allocatable,dimension(:)::temp1 if(dbg_set(dbg_sbr)) write(ipt,*) "Start: calc_sed_bedload " allocate(btxn(0:mt)); allocate(btyn(0:mt)); allocate(btx(0:nt)); allocate(bty(0:nt)); allocate(bflxe(0:nt)); allocate(bflxn(0:mt)); ! !----------------------------------------------------------------------- ! Compute maximum bottom stress for MPM bedload or suspended load. !----------------------------------------------------------------------- ! if(BEDLOAD_MPM)then tau_w = taub /rho_water endif ! !----------------------------------------------------------------------- ! Compute bedload sediment transport. !----------------------------------------------------------------------- ! ! Compute some constant bed slope parameters. ! sed_angle=DTAN(33.0_sp*pi/180.0_sp) ! ! Compute angle between currents and waves (radians). ! do i=1,m # if defined (WAVE_CURRENT_INTERACTION) if(BEDLOAD_SOULSBY)then if(MB_BBL_USE.or.SG_BBL_USE.or.SSW_BBL_USE)then ! ! Compute angle between currents and waves, measure CCW from current ! direction toward wave vector. ! IF (bustrc(i).eq.0.0_sp) THEN phic(i)=0.5_sp*pi*SIGN(1.0_sp,bvstrc(i)) ELSE phic(i)=ATAN2(bvstrc(i),bustrc(i)) ENDIF phicw(i)=1.5_sp*pi-Dwave(i)-phic(i) !-angler(i) ! ! Compute stress components at rho points. ! tau_cur=SQRT(bustrc(i)**2+bvstrc(i)**2) tau_wav(i)=SQRT(bustrw(i)**2+bvstrw(j)**2) tau_mean(i)=tau_cur*(1.0_sp+1.2_sp*((tau_wav(i)/(tau_cur& &+tau_wav(i)+eps))**3.2_sp)) endif endif # endif ! if(BEDLOAD_MPM)then ! && !defined BSTRESS_UPWIND # if defined (WAVE_CURRENT_INTERACTION) if(MB_BBL_USE.or.SG_BBL_USE.or.SSW_BBL_USE)then cff1=bustrcwmax(i) cff2=bvstrcwmax(i) else cff1=wubot_n(i) cff2=wvbot_n(i) end if # else cff1=wubot_n(i) cff2=wvbot_n(i) # endif Umag=SQRT(cff1*cff1+cff2*cff2)+eps angleu(i)=cff1/Umag anglev(i)=cff2/Umag endif end do ! DO ised=NCS+1,NST btxn = 0.0_sp btyn = 0.0_sp btx = 0.0_sp bty = 0.0_sp bflxe = 0.0_sp bflxn = 0.0_sp smgd=(sed(ised)%Srho/rho0-1.0_sp)*grav*sed(ised)%Sd50 osmgd=1.0_sp/smgd smgdr=SQRT(smgd)*sed(ised)%Sd50*sed(ised)%Srho ! do i=1,m # if defined (WAVE_CURRENT_INTERACTION) if(BEDLOAD_SOULSBY)then if(MB_BBL_USE.or.SG_BBL_USE.or.SSW_BBL_USE)then ! ! Compute wave asymmetry factor, based on Fredosoe and Deigaard. ! Dstp=d(i) waven=2.0_sp*pi/(wlen(i)+eps) wavec=SQRT(grav/waven*tanh(waven*Dstp)) cff4=MIN(waven*Dstp,kdmax) cff1=-0.1875_sp*wavec*(waven*Dstp)**2/(SINH(cff4))**4 cff2=0.125_sp*grav*HSC1(i)**2/(wavec*Dstp+eps) cff3=pi*HSC1(i)/(Pwave_bot(i)*SINH(cff4)+eps) ! ! Compute wave asymmetry factor, based on the note of Soulsby ! cff1=MIN(0.375_sp*(HSC1(i)/Dstp)* & & ((waven*Dstp)/(SINH(cff4))**3),0.15_sp) w_asym=2.0_sp*cff1/(1.0_sp+cff1**2.0_sp) ! ! Compute nondimensional stresses. ! theta_wav=tau_wav(i)*osmgd+eps theta_mean=tau_mean(i)*osmgd ! cff1=theta_wav*(1.0_sp+w_asym) cff2=theta_wav*(1.0_sp-w_asym) theta_max1=SQRT((theta_mean+ & & cff1*COS(phicw(i)))**2+ & & (cff1*SIN(phicw(i)))**2) theta_max2=SQRT((theta_mean+ & & cff2*COS(phicw(i)+pi))**2+ & & (cff2*SIN(phicw(i)+pi))**2) theta_max=MAX(theta_max1,theta_max2) ! ! Motion initiation factor. ! cff3=0.5_sp*(1.0_sp+SIGN(1.0_sp, & & theta_max/sed(ised)%tau_ce-1.0_sp)) ! ! Calculate bed loads in direction of current and perpendicular ! direction. ! phi_x1=12.0_sp*SQRT(theta_mean)* & & MAX((theta_mean-sed(ised)%tau_ce),0.0_sp) phi_x2=12.0_sp*(0.9534_sp+0.1907*COS(2.0_sp*phicw(i)))* & & SQRT(theta_wav)*theta_mean+ & & 12.0_sp*(0.229_sp*w_asym*theta_wav**1.5_sp* & & COS(phicw(i))) ! phi_x=MAX(phi_x1,phi_x2) ! <- original IF (ABS(phi_x2).gt.phi_x1) THEN phi_x=phi_x2 ELSE phi_x=phi_x1 END IF bedld_x=phi_x*smgdr*cff3 ! cff5=theta_wav**1.5_sp+1.5_sp*(theta_mean**1.5_sp) phi_y=12.0_sp*0.1907_sp*theta_wav*theta_wav* & & (theta_mean*SIN(2.0_sp*phicw(i))+1.2_sp*w_asym* & & theta_wav*SIN(phicw(i)))/cff5*cff3 bedld_y=phi_y*smgdr ! ! Partition bedld into x- and y- directions. ! (btxn and btyn have dimensions of kg/m). ! btxn(i)=(bedld_x*COS(phic(i))-bedld_y*SIN(phic(i)))*DTsed btyn(i)=(bedld_x*SIN(phic(i))+bedld_y*COS(phic(i)))*DTsed end if end if # endif if(BEDLOAD_MPM)then ! ! Magnitude of bed load at rho points. Meyer-Peter Muller formulation. ! (BEDLD has dimensions of kg m-1 s-1). ! bedld=8.0_sp*(MAX((tau_w(i)*osmgd-0.047_sp), & & 0.0_sp)**1.5_sp)*smgdr ! ! Partition bedld into xi and eta directions, still at rho points. ! (btxn and btyn have dimensions of kg/m ). ! ! btxn(i)=angleu(i)*bedld*on_r(i)*DTsed ! original roms formular ! btyn(i)=anglev(i)*bedld*om_r(i)*DTsed ! original roms formular btxn(i) = angleu(i)*bedld*DTsed btyn(i) = anglev(i)*bedld*DTsed endif ! ! Correct for along-direction slope. Limit slope to 0.9*sed angle. ! ! cff1=0.5_sp*(1.0_sp+SIGN(1.0_sp,FX_r(i))) ! cff2=0.5_sp*(1.0_sp-SIGN(1.0_sp,FX_r(i))) ! cff3=0.5_sp*(1.0_sp+SIGN(1.0_sp,FE_r(i))) ! cff4=0.5_sp*(1.0_sp-SIGN(1.0_sp,FE_r(i))) ! if(SLOPE_NEMETH)then ! dzdx=(h(i+1,j)-h(i,j))/om_u(i+1,j)*cff1+ & ! & (h(i-1,j)-h(i,j))/om_u(i ,j)*cff2 ! dzdy=(h(i,j+1)-h(i,j))/on_v(i,j+1)*cff3+ & ! & (h(i,j-1)-h(i,j))/on_v(i ,j)*cff4 ! a_slopex=1.7_sp*dzdx ! a_slopey=1.7_sp*dzdy ! ! Add contriubiton of bed slope to bed load transport fluxes. ! ! FX_r(i,j)=FX_r(i,j)*(1.0_sp+a_slopex) ! FE_r(i,j)=FE_r(i,j)*(1.0_sp+a_slopey) ! ! elseif(SLOPE_LESSER)then ! dzdx=MIN(((h(i+1,j)-h(i ,j))/om_u(i+1,j)*cff1+ & ! & (h(i ,j)-h(i-1,j))/om_u(i ,j)*cff2),0.52_sp)* & ! & SIGN(1.0_sp,FX_r(i,j)) ! dzdy=MIN(((h(i,j+1)-h(i,j ))/on_v(i,j+1)*cff3+ & ! & (h(i,j )-h(i,j-1))/on_v(i ,j)*cff4),0.52_sp)* & ! & SIGN(1.0_sp,FE_r(i,j)) ! cff=DATAN(dzdx) ! a_slopex=sed_angle/(COS(cff)*(sed_angle-dzdx)) ! cff=DATAN(dzdy) ! a_slopey=sed_angle/(COS(cff)*(sed_angle-dzdy)) ! ! Add contriubiton of bed slope to bed load transport fluxes. ! ! FX_r(i,j)=FX_r(i,j)*a_slopex ! FE_r(i,j)=FE_r(i,j)*a_slopey ! endif end do ! calculate the bed slope at cell center do i=1,n n1 = nv(i,1) ; n2 = nv(i,2) ; n3 = nv(i,3) dbdx_e(i) = awx(i,1)*h(n1)+awx(i,2)*h(n2)+awx(i,3)*h(n3) dbdy_e(i) = awy(i,1)*h(n1)+awy(i,2)*h(n2)+awy(i,3)*h(n3) end do ! interpolate slope into node CALL E2N2D(dbdx_e,dbdx_n) CALL E2N2D(dbdy_e,dbdy_n) !----------CALCULATE DERIVATIVES OF DEPTH WITH X AND Y AT NODES----------------! !CALL DEPTH_GRADIENT if(SLOPE_NEMETH)then ! ! Correct for along-direction slope. Limit slope to 0.9*sed angle. ! bld_dr = sign(1.0_sp,btxn) dzdx=dbdx_n*bld_dr bld_dr = sign(1.0_sp,btyn) dzdy=dbdy_n*bld_dr a_slopex=1.7_sp*dzdx a_slopey=1.7_sp*dzdy ! Add contriubiton of bed slope to bed load transport fluxes. ! btxn=btxn*(1.0_sp+a_slopex) btyn=btyn*(1.0_sp+a_slopey) elseif(SLOPE_LESSER)then sed_slope = tan(sed_angle*3.14159/180.) max_slope = 0.8*sed_slope bld_dr = sign(1.0_sp,btxn) dzdx=MIN(dbdx_n*bld_dr,max_slope) bld_dr = sign(1.0_sp,btyn) dzdy=MIN(dbdy_n*bld_dr,max_slope) allocate(cff(0:mt)) cff=ATAN(dzdx) a_slopex=sed_angle/(COS(cff)*(sed_angle-dzdx)) cff=ATAN(dzdy) a_slopey=sed_angle/(COS(cff)*(sed_angle-dzdy)) deallocate(cff) ! ! Add contriubiton of bed slope to bed load transport fluxes. ! btxn=btxn*a_slopex btyn=btyn*a_slopey end if if(SED_MORPH)then ! ! Apply morphology factor. ! btxn=btxn*sed(ised)%morph_fac btyn=btyn*sed(ised)%morph_fac endif ! ! Apply bedload transport rate coefficient. Also limit ! bedload to the fraction of each sediment class. ! btxn=btxn*bedload_coeff*sedbed%bed_frac(:,1,ised) btyn=btyn*bedload_coeff*sedbed%bed_frac(:,1,ised) ! ! Limit bed load to available bed mass. ! ! ROMS does this, but doesn't make sense, mass is linked to divergence of ! transport, not horizontal fluxes , check with J. Warner ! bedld_mass=ABS(FX_r(i,j))+ABS(btyn(i,j))+eps ! FX_r=MIN(ABS(FX_r),sedbed%bed_mass(1:m,1,nstp,ised)*om_r*on_r*ABS(FX_r)/bedld_mass)*SIGN(1.0_sp,FX_r) ! btyn=MIN(ABS(btyn),sedbed%bed_mass(1:m,1,nstp,ised)*om_r*on_r*ABS(btyn)/bedld_mass)*SIGN(1.0_sp,btyn) !exchange btx/bty across interprocessor boundaries # if defined (MULTIPROCESSOR) if(par)then call aexchange(nc,myid,nprocs,btxn,btyn) endif # endif call n2e2d(btxn,btx) call n2e2d(btyn,bty) !add bed slope contribution to our horizontal fluxes !Eq 35, Warner etal, Comp & Geo, 2008 => qbl_slope = slopex,slopey !minimize the maximum downslope, but do not truncate upslope !upslope is limited by Lesser formula !range of slopefactor (slopex) = [.55,5.64] for .. !.. upslope(dhdx<0)=-infty -> downslope(dhdx>0) = infty !dbdx > 0 is a downslope (see bld_dr to enforce this) !sed_slope = tan(sed_angle*3.14159/180.) !max_slope = 0.8*sed_slope !do i=1,n ! n1 = nv(i,1) ; n2 = nv(i,2) ; n3 = nv(i,3) ! bld_dr = sign(1.0_sp,btx(i)) ! dbdx = min( (awx(i,1)*h(n1)+awx(i,2)*h(n2)+awx(i,3)*h(n3))*bld_dr ,max_slope) ! beta = atan(dbdx) ! slopex = sed_slope/( (sed_slope-dbdx)*cos(beta)) ! btx(i) = btx(i)*slopex ! bld_dr = sign(1.0_sp,bty(i)) ! dbdy = min( (awy(i,1)*h(n1)+awy(i,2)*h(n2)+awy(i,3)*h(n3))*bld_dr ,max_slope) ! beta = atan(dbdy) ! slopey = sed_slope/( (sed_slope-dbdy)*cos(beta)) ! bty(i) = bty(i)*slopey !end do # if defined (WAVE_CURRENT_INTERACTION) if(MB_BBL_USE.or.SG_BBL_USE.or.SSW_BBL_USE)then call n2e2d(bustrcwmax,bustrcwmax_e) call n2e2d(bvstrcwmax,bvstrcwmax_e) end if # endif !calculate the flux on the elements using first order upwind formulation do i=1,ne ia=iec(i,1) ib=iec(i,2) j1=ienode(i,1) j2=ienode(i,2) !bed stress in x and y averaged to an edge , used for upwind selection # if defined (WAVE_CURRENT_INTERACTION) if(MB_BBL_USE.or.SG_BBL_USE.or.SSW_BBL_USE)then tau_x = -0.5*(bustrcwmax_e(ia)+bustrcwmax_e(ib)) tau_y = -0.5*(bvstrcwmax_e(ia)+bvstrcwmax_e(ib)) else tau_x = 0.5*(wubot(ia)+wubot(ib)) tau_y = 0.5*(wvbot(ia)+wvbot(ib)) endif # else tau_x = 0.5*(wubot(ia)+wubot(ib)) tau_y = 0.5*(wvbot(ia)+wvbot(ib)) # endif !dot the bed stress vector with length-weighted edge normal !note bed stress opposes flow dir so we need neg sign !the vector (-dltyc,dltxc) points from ia to ib flowdir = -(-tau_x*dltyc(i) + tau_y*dltxc(i)) ! -(tau dot dL) !upwind select the x-dir and y-dir bedload transport bldx = ((one-sign(one,flowdir))*btx(ib)+(one+sign(one,flowdir))*btx(ia))*0.5_sp bldy = ((one-sign(one,flowdir))*bty(ib)+(one+sign(one,flowdir))*bty(ia))*0.5_sp !accum contributions to bedload transport divergence calculations at elements ia/ib !on either side of edge i ! => bedload in kg bflxe(ia) = bflxe(ia) - bldx*dltyc(i) + bldy*dltxc(i) bflxe(ib) = bflxe(ib) + bldx*dltyc(i) - bldy*dltxc(i) end do !limit fluxes to prevent bottom from breaking thru water surface. allocate(cff(0:nt)); cff = 0.0_sp do i=1,n ! Total thickness available. Dstp=MAX(d1(i),0.0_sp) ! Thickness change that wants to occur. n1 = nv(i,1) ; n2 = nv(i,2) ; n3 = nv(i,3) poro=(bed(n1,1,iporo)+bed(n2,1,iporo)+bed(n3,1,iporo))/3.0 bed_change= bflxe(i)/(sed(ised)%Srho*(1.0_sp-poro)) ! calculate the coefficients to make sure the change breaking ! the water surface cff(i)=MAX(bed_change-Dstp,0.0_sp) cff(i)=cff(i)/ABS(bed_change+eps) ! Limit that change to be less than available. bflxe(i) = bflxe(i)*(1-cff(i)) end do ! re-Limit the UV transport components that accumulated to bedload allocate(iab_limit(0:nt)); iab_limit = 0 ! flag to avoid multiple cutting do i=1,ne ia=iec(i,1) ib=iec(i,2) if(iab_limit(ia)==0)then btx(ia) = btx(ia)*(1-cff(ia)) bty(ia) = bty(ia)*(1-cff(ia)) iab_limit(ia)=1 endif if(iab_limit(ib)==0)then btx(ib) = btx(ib)*(1-cff(ib)) bty(ib) = bty(ib)*(1-cff(ib)) iab_limit(ib)=1 endif end do deallocate(iab_limit) deallocate(cff) !divide by element area to complete divergence => bedload in kg/m^2 !if there was a divergence of bedload transport, this should be a positive number do i=1,n bflxe(i) = bflxe(i)/art(i) end do !exchange bflxe across interprocessor boundaries # if defined (MULTIPROCESSOR) if(par)then call aexchange(ec,myid,nprocs,bflxe) call aexchange(ec,myid,nprocs,btx,bty) endif # endif ! interpolate element-based blfx to the nodes call e2n2d(bflxe,bflxn) call e2n2d(btx,btxn) call e2n2d(bty,btyn) # if defined (MULTIPROCESSOR) if(par)then call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,1,myid,nprocs,bflxn) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,1,myid,nprocs,btxn) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,1,myid,nprocs,btyn) call aexchange(nc,myid,nprocs,bflxn) call aexchange(nc,myid,nprocs,btxn,btyn) endif # endif ! evaluate change in bed properties. sedbed%bed_mass(1:m,1,nnew,ised)=MAX(sedbed%bed_mass(1:m,1,nstp,ised)& &-bflxn(1:m),0.0_sp) if(.not.SUSLOAD)then DO k=2,Nbed sedbed%bed_mass(1:m,k,nnew,ised)=sedbed%bed_mass(1:m,k,nstp,ised) END DO endif bed(1:m,1,ithck)=MAX(bed(1:m,1,ithck)- bflxn(1:m)/(sed(ised)%Srho*(1.0_sp& &-bed(1:m,1,iporo))),0.0_sp) ! !----------------------------------------------------------------------- ! Output bedload fluxes. !----------------------------------------------------------------------- ! sedbed%bedldu(1:m,ised)=btxn(1:m)/Dtsed sedbed%bedldv(1:m,ised)=btyn(1:m)/Dtsed END DO ! ! Need to update bed mass for the cohesive sediment types, becasue ! they did not partake in the bedload transport. ! DO ised=1,NCS sedbed%bed_mass(1:m,1,nnew,ised)=sedbed%bed_mass(1:m,1,nstp,ised) if(.not.SUSLOAD)then DO k=2,Nbed sedbed%bed_mass(1:m,k,nnew,ised)=sedbed%bed_mass(1:m,k,nstp,ised) END DO endif END DO ! ! Update mean surface properties. ! Sd50 must be positive definite, due to BBL routines. ! Srho must be >1000, due to (s-1) in BBL routines. ! DO i=1,m cff3=0.0_sp DO ised=1,NST cff3=cff3+sedbed%bed_mass(i,1,nnew,ised) END DO IF (cff3.eq.0.0_sp) THEN cff3=eps END IF DO ised=1,NST sedbed%bed_frac(i,1,ised)=sedbed%bed_mass(i,1,nnew,ised)/cff3 END DO End DO Call Update_Bottom_Properties deallocate(btx,bty,btxn,btyn,bflxe,bflxn) RETURN if(dbg_set(dbg_sbr)) write(ipt,*) "End: calc_sed_bedload " End Subroutine calc_sed_bedload !=================================================================== === ! ! ! This routine computes the vertical settling (sinking) of suspended ! ! sediment via a semi-Lagrangian advective flux algorithm. It uses a ! ! parabolic, vertical reconstructuion of the suspended sediment in ! ! the water column with PPT/WENO constraints to avoid oscillations. ! ! ! ! References: ! ! ! ! Colella, P. and P. Woodward, 1984: The piecewise parabolic method ! ! (PPM) for gas-dynamical simulations, J. Comp. Phys., 54, 174-201. ! ! ! ! Liu, X.D., S. Osher, and T. Chan, 1994: Weighted essentially ! ! nonoscillatory shemes, J. Comp. Phys., 115, 200-212. ! ! ! ! Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G. ! ! Arango, 2008: Development of a three-dimensional, regional, ! ! coupled wave, current, and sediment-transport model, Computers ! ! & Geosciences, 34, 1284-1306. ! ! ! !======================================================================= Subroutine calc_sed_settling implicit none ! ! Local variable declarations. ! integer :: i, indx, ised, j, k, ks real(sp) :: cff, cu, cffL, cffR, dltL, dltR integer, dimension(0:kb) :: ksource real(sp), dimension(0:kb) :: FC real(sp), dimension(0:kb) :: Hz_inv real(sp), dimension(0:kb) :: Hz_inv2 real(sp), dimension(0:kb) :: Hz_inv3 real(sp), dimension(0:kb) :: qc real(sp), dimension(0:kb) :: qR real(sp), dimension(0:kb) :: qL real(sp), dimension(0:kb) :: WR real(sp), dimension(0:kb) :: WL if(dbg_set(dbg_sbr)) write(ipt,*) "Start: calc_sed_settling " ! !----------------------------------------------------------------------- ! Add sediment vertical sinking (settling) term. !----------------------------------------------------------------------- ! ! Compute inverse thicknesses to avoid repeated divisions. ! ! DO ised=1,NST ! print*,'before settling:',ised,maxval(sed(ised)%conc),minval(sed(ised)%conc) ! end do do i=1,m DO k=1,kbm1 Hz_inv(k)=1.0_sp/(dt(i)*dz(i,k)) END DO DO k=2,kbm1 Hz_inv2(k)=1.0_sp/(dt(i)*dz(i,k)+dt(i)*dz(i,k-1)) END DO DO k=2,kbm1-1 Hz_inv3(k)=1.0_sp/(dt(i)*(dz(i,k-1)+dz(i,k)+dz(i,k+1))) END DO ! ! Copy concentration of suspended sediment into scratch array "qc" ! (q-central, restrict it to be positive) which is hereafter ! interpreted as a set of grid-box averaged values for sediment ! concentration. ! SED_LOOP: DO ised=1,NST DO k=1,kbm1 qc(k)=sed(ised)%conc(i,k) END DO ! !----------------------------------------------------------------------- ! Vertical sinking of suspended sediment. !----------------------------------------------------------------------- ! ! Reconstruct vertical profile of suspended sediment "qc" in terms ! of a set of parabolic segments within each grid box. Then, compute ! semi-Lagrangian flux due to sinking. ! DO k=2,kbm1 FC(k)=(qc(k-1)-qc(k))*Hz_inv2(k) END DO DO k=kbm1-1,2,-1 dltR=(dt(i)*dz(i,k))*FC(k) dltL=(dt(i)*dz(i,k))*FC(k+1) cff=dt(i)*(dz(i,k-1)+2.0_sp*dz(i,k)+dz(i,k+1)) cffR=cff*FC(k) cffL=cff*FC(k+1) ! ! Apply PPM monotonicity constraint to prevent oscillations within the ! grid box. ! IF ((dltR*dltL).le.0.0_sp) THEN dltR=0.0_sp dltL=0.0_sp ELSE IF (ABS(dltR).gt.ABS(cffL)) THEN dltR=cffL ELSE IF (ABS(dltL).gt.ABS(cffR)) THEN dltL=cffR END IF ! ! Compute right and left side values (qR,qL) of parabolic segments ! within grid box Hz(k); (WR,WL) are measures of quadratic variations. ! ! NOTE: Although each parabolic segment is monotonic within its grid ! box, monotonicity of the whole profile is not guaranteed, ! because qL(k+1)-qR(k) may still have different sign than ! qc(k+1)-qc(k). This possibility is excluded, after qL and qR ! are reconciled using WENO procedure. ! cff=(dltR-dltL)*Hz_inv3(k) dltR=dltR-cff*(dt(i)*dz(i,k-1)) dltL=dltL+cff*(dt(i)*dz(i,k+1)) qR(k)=qc(k)+dltR qL(k)=qc(k)-dltL WR(k)=(2.0_sp*dltR-dltL)**2 WL(k)=(dltR-2.0_sp*dltL)**2 END DO cff=1.0E-14_sp DO k=kbm1-1,3,-1 dltL=MAX(cff,WL(k )) dltR=MAX(cff,WR(k-1)) qR(k)=(dltR*qR(k)+dltL*qL(k-1))/(dltR+dltL) qL(k-1)=qR(k) END DO FC(1)=0.0_sp ! no-flux boundary condition if(LINEAR_CONTINUATION)then qL(1)=qR(2) qR(1)=2.0_sp*qc(1)-qL(1) elseif(NEUMANN)then qL(1)=qR(2) qR(1)=1.5_sp*qc(1)-0.5_sp*qL(1) else qR(1)=qc(1) ! default strictly monotonic qL(1)=qc(1) ! conditions qR(2)=qc(1) endif if(LINEAR_CONTINUATION)then qR(kbm1)=qL(kbm1-1) qL(kbm1)=2.0_sp*qc(kbm1)-qR(kbm1) elseif(NEUMANN)then qR(kbm1)=qL(kbm1-1) qL(kbm1)=1.5_sp*qc(kbm1)-0.5_sp*qR(kbm1) else qL(kbm1-1)=qc(kbm1) ! bottom grid boxes are qR(kbm1)=qc(kbm1) ! re-assumed to be qL(kbm1)=qc(kbm1) ! piecewise constant. endif ! ! Apply monotonicity constraint again, since the reconciled interfacial ! values may cause a non-monotonic behavior of the parabolic segments ! inside the grid box. ! DO k=kbm1,1,-1 dltR=qR(k)-qc(k) dltL=qc(k)-qL(k) cffR=2.0_sp*dltR cffL=2.0_sp*dltL IF ((dltR*dltL).lt.0.0_sp) THEN dltR=0.0_sp dltL=0.0_sp ELSE IF (ABS(dltR).gt.ABS(cffL)) THEN dltR=cffL ELSE IF (ABS(dltL).gt.ABS(cffR)) THEN dltL=cffR END IF qR(k)=qc(k)+dltR qL(k)=qc(k)-dltL END DO ! ! After this moment reconstruction is considered complete. The next ! stage is to compute vertical advective fluxes, FC. It is expected ! that sinking may occurs relatively fast, the algorithm is designed ! to be free of CFL criterion, which is achieved by allowing ! integration bounds for semi-Lagrangian advective flux to use as ! many grid boxes in upstream direction as necessary. ! ! In the two code segments below, WL is the z-coordinate of the ! departure point for grid box interface z_w with the same indices; ! FC is the finite volume flux; ksource(:,k) is index of vertical ! grid box which contains the departure point (restricted by kbm1). ! During the search: also add in content of whole grid boxes ! participating in FC. ! cff=DTsed*ABS(sed(ised)%Wsed) DO k=kbm1,1,-1 FC(k+1)=0.0_sp WL(k)=dt(i)*zz(i,k+1)+cff WR(k)=dt(i)*dz(i,k)*qc(k) ksource(k)=k END DO DO k=kbm1,1,-1 DO ks=k,2,-1 IF (WL(k).gt.dt(i)*zz(i,ks)) THEN ksource(k)=ks-1 FC(k+1)=FC(k+1)+WR(ks) END IF END DO END DO ! ! Finalize computation of flux: add fractional part. ! DO k=kbm1,1,-1 ks=ksource(k) cu=MIN(1.0_sp,(WL(k)-dt(i)*zz(i,ks+1))*Hz_inv(ks)) FC(k+1)=FC(k+1)+ & & dt(i)*dz(i,ks)*cu* & & (qL(ks)+ & & cu*(0.5_sp*(qR(ks)-qL(ks))- & & (1.5_sp-cu)* & & (qR(ks)+qL(ks)-2.0_sp*qc(ks)))) END DO DO k=kbm1,1,-1 sed(ised)%conc(i,k)=sed(ised)%conc(i,k)+(FC(k)-FC(k+1)) END DO sedbed%settling_flux(i,ised)=FC(kb) END DO SED_LOOP END DO ! DO ised=1,NST ! print*,'after settling:',ised,maxval(sed(ised)%conc),minval(sed(ised)%conc) ! end do if(dbg_set(dbg_sbr)) write(ipt,*) "End: calc_sed_settling " RETURN End Subroutine calc_sed_settling !======================================================================= ! ! ! This computes sediment bed and water column exchanges: deposition, ! ! resuspension, and erosion. ! ! ! ! References: ! ! ! ! Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G. ! ! Arango, 2008: Development of a three-dimensional, regional, ! ! coupled wave, current, and sediment-transport model, Computers ! ! & Geosciences, 34, 1284-1306. ! ! ! !======================================================================= Subroutine calc_sed_fluxes implicit none ! ! Local variable declarations. ! integer :: Ksed, i, indx, ised, j, k, ks integer :: bnew, cnt,jj,cc real(sp) :: cff, cff1, cff2, cff3, cff4, cff5 real(sp), dimension(0:mt) :: tau_w if(dbg_set(dbg_sbr)) write(ipt,*) "Start: calc_sed_fluxes" if(BEDLOAD)then bnew=nnew else bnew=nstp endif ! !----------------------------------------------------------------------- ! Compute sediment deposition, resuspension, and erosion. !----------------------------------------------------------------------- ! if(BEDLOAD_MPM.or.SUSLOAD)then tau_w = taub endif # if !defined (FLUID_MUD) ! !----------------------------------------------------------------------- ! Sediment deposition and resuspension near the bottom. !----------------------------------------------------------------------- ! ! The deposition and resuspension of sediment on the bottom "bed" ! is due to precepitation settling_flux, already computed, and the ! resuspension (erosion, hence called ero_flux). The resuspension is ! applied to the bottom-most grid box value qc(:,1) so the total mass ! is conserved. Restrict "ero_flux" so that "bed" cannot go negative ! after both fluxes are applied. ! NODE_LOOP : DO i=1,m # if defined (WET_DRY) if(iswetn(i)/=1)cycle ! no calculation when dry node # endif SED_LOOP: DO ised=1,NST ! ! Calculate critical shear stress in Pa ! if(COHESIVE_BED)then ! cff =1.0_sp/bed(i,1,ibtcr) cff=1.0_sp/sed(ised)%tau_ce elseif(MIXED_BED)then cff = MAX(bottom(i,idprp)*bed(i,1,ibtcr)+ & & (1.0_sp-bottom(i,idprp))*sed(ised)%tau_ce, & & sed(ised)%tau_ce) cff=1.0_sp/cff else cff=1.0_sp/sed(ised)%tau_ce endif ! ! Compute erosion, ero_flux (kg/m2). ! cff1=(1.0_sp-bed(i,1,iporo))*sedbed%bed_frac(i,1,ised) cff2=DTsed*sed(ised)%Erate*cff1 cff3=sed(ised)%Srho*cff1 cff4=sedbed%bed_mass(i,1,bnew,ised) cff5=sedbed%settling_flux(i,ised) !CRS if(SED_TAU_CD_CONST.or.SED_TAU_CD_LIN)then if (tau_w(i).GT.sed(ised)%tau_cd) then cff5=0.0_sp endif endif if(SED_TAU_CD_LIN)then if (tau_w(i).LE.sed(ised)%tau_cd.AND. & & sed(ised)%tau_cd.GT.0.0_sp) then cff5=(1.0_sp-tau_w(i)/sed(ised)%tau_cd)*cff5 endif endif sedbed%ero_flux(i,ised)= & & MIN(MAX(0.0_sp,cff2*(cff*tau_w(i)-1.0_sp)), & & MIN(cff3*bottom(i,iactv),cff4)+ & & cff5) ! ! Update global tracer variables (mT units) for erosive flux. ! sed(ised)%conc(i,kbm1)=sed(ised)%conc(i,kbm1)+sedbed%ero_flux(i,ised) if (SED_TAU_CD_CONST .or. SED_TAU_CD_LIN)then sed(ised)%conc(i,kbm1)=sed(ised)%conc(i,kbm1) + & & (sedbed%settling_flux(i,ised)-cff5) sedbed%settling_flux(i,ised)=cff5 endif END DO SED_LOOP END DO NODE_LOOP # endif ! J. Ge for fluid mud layer # if defined (FLUID_MUD) Settling = 0.0_SP do ised=1,nst do i=1,m # if defined (WET_DRY) if(iswetn(i)==1)then # endif Settling(i) = Settling(i) + sed(ised)%dflx(i) # if defined (WET_DRY) else Settling(i) = 0.0_sp end if # endif end do end do # endif if(dbg_set(dbg_sbr)) write(ipt,*) "End: calc_sed_fluxes" return End Subroutine calc_sed_fluxes !======================================================================= ! ! ! This computes sediment biomass due to vegation growth ! ! ! ! References: ! ! ! ! ! !======================================================================= Subroutine calc_sed_biomass implicit none ! ! Local variable declarations. ! integer :: i, j, k, ised integer :: sstp, nbio_steps real(sp) :: cff, Dstp real(sp) :: sgr_kgmmol if(dbg_set(dbg_sbr)) write(ipt,*) "Start: calc_sed_biomass " ! !----------------------------------------------------------------------- ! Compute !----------------------------------------------------------------------- ! ! Compute number of model steps for each hour. ! nbio_steps=MAX(1,INT(3600.0_sp/DTsed)) ! ! Compute number of hourly values we need to save. ! If we want a 1 day avg, then need 24 values. ! NODE_LOOP : DO i=1,m # if defined (WET_DRY) if(iswetn(i)/=1)cycle ! no calculation when dry node # endif ! ! Only update the max depth once per hour. ! IF (MOD(iint,nbio_steps).eq.0) THEN ! ! Determine the index for placement of new value. ! ! sstp=1+MOD(iint-ntstart,24) sstp=1+MOD(iint,24) ! ! Save instantaneous depth at this instance and recompute max daily depth. ! Dstp=d(i) sedbed%Dstp_max(i,sstp)=Dstp cff=0.0_sp DO k=1,nTbiom cff=MAX(cff,sedbed%Dstp_max(i,k)) END DO bottom(i,imaxD)=cff END IF ! ! Seagrass as a bottom property ! if(SEAGRASS_BOTTOM)then if(SEAGRASS_SINK)then ! sgr_diam=0.01_sp ! SgrN has units of millimole_nitrogen meter-3 ! change moles to kg (2.8e-5 kg/millimole N2) switch to right formula ! cylinder height is mass over (density*pi*r*r) ! shoot height is cylinder height / shoot density sgr_kgmmol = 2.8e-5_sp ! sgr_density = 500.0_sp ! sgr_Hthres = 1.25_sp cff = SgrN(i)*sgr_kgmmol/ & & (0.25_sp*sgr_diam*sgr_diam*pi*sgr_density) ! bottom(i,isgrH) = cff/bottom(i,isgrD) ! IF (bottom(i,isgrH).gt.sgr_Hthres) THEN bottom(i,isgrD)=bottom(i,isgrD)* & & bottom(i,isgrH)/sgr_Hthres bottom(i,isgrH) = sgr_Hthres END IF else bottom(i,isgrH)=1.25_sp bottom(i,isgrD)=400.0_sp endif endif ! ! Update settling flux for depositing bio mass. ! ! ! Require (for now) that the first sed class be the new biomass. ! ised=1 cff=0.0_sp ! remove this line ! cff= funct( bottom(i,imaxD), DTsed) ! need real eq. in kg/m^2 sedbed%settling_flux(i,ised)=sedbed%settling_flux(i,ised)+cff END DO NODE_LOOP if(dbg_set(dbg_sbr)) write(ipt,*) "End: calc_sed_biomass " return End Subroutine calc_sed_biomass !======================================================================= ! ! ! This routine computes sediment bed layer stratigraphy. ! ! ! ! Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G. ! ! Arango, 2008: Development of a three-dimensional, regional, ! ! coupled wave, current, and sediment-transport model, Computers ! ! & Geosciences, 34, 1284-1306. ! ! ! !======================================================================= Subroutine calc_sed_bed_cohesive implicit none ! ! Local variable declarations. ! integer :: Ksed, i, ised,ibed, j, k, ks integer :: bnew, nnn, cnt,cc,jj real(sp), parameter :: eps = 1.0E-14_sp real(sp) :: cff, cff1, cff2, cff3,cff4 real(sp) :: thck_avail, thck_to_add real(sp), dimension(NST) :: nlysm real(sp), dimension(NST) :: dep_mass real(sp), dimension(0:mt) :: tau_w real(sp) :: pcoh real(sp) :: alpha, bzactv, frt, tcb_top, tcb_bot real(sp), dimension(Nbed) :: tcr real(sp), dimension(Nbed) :: bmz real(sp), dimension(NCS) :: masseq real(sp),allocatable,dimension(:,:)::temp if(dbg_set(dbg_sbr)) write(ipt,*) "Start: calc_sed_bed_cohesive" if(BEDLOAD)then bnew=nnew else bnew=nstp endif ! KLUDGE alert ! minlayer_thick = 0.0005 minlayer_thick = newlayer_thick ! !----------------------------------------------------------------------- ! Compute sediment bed layer stratigraphy. !----------------------------------------------------------------------- ! if(BEDLOAD_MPM.or.SUSLOAD)then tau_w = taub endif ! !----------------------------------------------------------------------- ! Update bed properties according to ero_flux and dep_flux. !----------------------------------------------------------------------- ! if(SUSLOAD)then I_LOOP : DO i=1,m ! ! The deposition and resuspension of sediment on the bottom "bed" ! is due to precipitation flux FC(:,0), already computed, and the ! resuspension (erosion, hence called ero_flux). The resuspension is ! applied to the bottom-most grid box value qc(:,1) so the total mass ! is conserved. Restrict "ero_flux" so that "bed" cannot go negative ! after both fluxes are applied. ! SED_LOOP: DO ised=1,NST dep_mass(ised)=0.0_sp if(SED_MORPH)then ! Apply morphology factor. sedbed%ero_flux(i,ised)=sedbed%ero_flux(i,ised)*sed(ised)%morph_fac sedbed%settling_flux(i,ised)=sedbed%settling_flux(i,ised)* & & sed(ised)%morph_fac endif ! Update bed mass arrays. sedbed%bed_mass(i,1,nnew,ised)=MAX(sedbed%bed_mass(i,1,bnew,ised)- & & (sedbed%ero_flux(i,ised)- & & sedbed%settling_flux(i,ised)), & & 0.0_sp) DO k=2,Nbed sedbed%bed_mass(i,k,nnew,ised)=sedbed%bed_mass(i,k,nstp,ised) END DO END DO SED_LOOP cff3=0.0_sp DO ised=1,NST cff3=cff3+sedbed%bed_mass(i,1,nnew,ised) END DO IF (cff3.eq.0.0_sp) THEN cff3=eps END IF bed(i,1,ithck)=0.0_sp DO ised=1,NST sedbed%bed_frac(i,1,ised)=sedbed%bed_mass(i,1,nnew,ised)/cff3 bed(i,1,ithck)=MAX(bed(i,1,ithck)+ & & sedbed%bed_mass(i,1,nnew,ised)/ & & (sed(ised)%Srho* & & (1.0_sp-bed(i,1,iporo))),0.0_sp) END DO END DO I_LOOP endif !/* SUSPLOAD section */ ! !----------------------------------------------------------------------- ! At this point, all deposition or erosion is complete, and ! has been added/subtracted to top layer. Thickness has NOT been corrected. !----------------------------------------------------------------------- ! I_LOOP_CB : DO i=1,m ! Calculate active layer thickness, bottom(i,iactv). ! (trunk version allows this to be zero...this has minimum of 6*D50) bottom(i,iactv)=MAX(0.0_sp, & & 0.007_sp* & & (tau_w(i)-bottom(i,itauc)))+ & & 6.0_sp*bottom(i,isd50) if(COHESIVE_BED)then bottom(i,iactv)=MAX(0.0_sp, & & 0.007_sp* & & (tau_w(i)-bed(i,1,ibtcr)))+ & & 6.0_sp*bottom(i,isd50) endif if(MIXED_BED)then cff1= MAX(0.0_sp, & & 0.007_sp* & & (tau_w(i)-bed(i,1,ibtcr)))+ & & 6.0_sp*bottom(i,isd50) cff2= MAX(0.0_sp, & & 0.007_sp* & & (tau_w(i)-bottom(i,itauc)))+ & & 6.0_sp*bottom(i,isd50) bottom(i,iactv)=MAX(cff1,cff2) endif if(SED_MORPH)then ! Apply morphology factor. bottom(i,iactv)=MAX(bottom(i,iactv)*sed(ised)%morph_fac, & & bottom(i,iactv)) endif if (COHESIVE_BED .or. MIXED_BED)then ! Find first layer with tc > tb ! Remember, the critical stresses apply to the TOP of each layer Ksed = 1 bzactv = 0.0_sp frt = 0.0_sp ! CRS cff1 = tau_w(i) tcb_top = bed(i,1,ibtcr) tcb_bot = bed(i,2,ibtcr) if (MIXED_BED)then ! Calc. tau crit for bottom of next layer ! Update cohesive property of seds in top layer cff3 = 0.0_sp cff4 = 1.0_sp DO ised=1,NCS cff3=cff3+sedbed%bed_frac(i,1,ised) cff4=cff4*sed(ised)%tau_ce**sedbed%bed_frac(i,1,ised) END DO DO ised=NCS+1,NST cff4=cff4*sed(ised)%tau_ce**sedbed%bed_frac(i,1,ised) ENDDO pcoh=min(max((cff3-transN)/(transC-transN) & & ,0.0_sp),1.0_sp) tcb_top = pcoh*bed(i,1,ibtcr)+(pcoh-1.0_sp)*cff4 !BF endif IF(cff1 .GT. tcb_top)THEN ! Calculate tcb_temp for next layer tcb_bot = bed(i,Ksed+1,ibtcr) if (MIXED_BED)then ! Recalculate cohesive fraction and mean grain tau crit ! Note that Ksed is used for grain props, and Ksed+1 for ! bed tau crit at bottom of layer Ksed cff3 = 0.0_sp cff4 = 1.0_sp DO ised=1,NCS cff3=cff3+sedbed%bed_frac(i,Ksed,ised) cff4=cff4*sed(ised)%tau_ce**sedbed%bed_frac(i,Ksed,ised) END DO DO ised=NCS+1,NST cff4=cff4*sed(ised)%tau_ce**sedbed%bed_frac(i,Ksed,ised) ENDDO ! Calculate cohesive behavior and blended tau crit pcoh=min(max((cff3-transN)/(transC-transN), & & 0.0_sp),1.0_sp) tcb_bot =pcoh*bed(i,Ksed+1,ibtcr)+(pcoh-1.0_sp)*cff4 !BF endif DO WHILE ( (Ksed.LE.(Nbed-1)) .AND. & & (cff1 .GT. tcb_bot)) !ALA Instead of adding entire 2nd layer, add just what is needed !ALA! Add entire layer !ALA bzactv = bzactv + bed(i,Ksed,ithck) ! Add thickness equivalent to difference between cff1 and tcb_bot IF (cff1.GT.bed(i,Ksed+1,ibtcr)) THEN bzactv = bzactv + bed(i,Ksed,ithck) tcb_top = tcb_bot tcb_bot = bed(i,Ksed+1,ibtcr) ELSE bzactv = bzactv + MIN(bed(i,Ksed,ithck), & & (MAX(0.0_sp,0.007_sp*(cff1-tcb_bot)))) tcb_top = cff1 tcb_bot = bed(i,Ksed+1,ibtcr) ENDIF if(MIXED_BED)then ! Recalculate cohesive fraction and mean grain tau crit cff3 = 0.0_sp cff4 = 1.0_sp DO ised=1,NCS cff3=cff3+sedbed%bed_frac(i,Ksed,ised) cff4=cff4*sed(ised)%tau_ce**sedbed%bed_frac(i,Ksed,ised) END DO DO ised=NCS+1,NST cff4=cff4*sed(ised)%tau_ce**sedbed%bed_frac(i,Ksed,ised) ENDDO ! Calculate cohesive behavior and blended tau crit pcoh=min(max((cff3-transN)/(transC-transN), & & 0.0_sp),1.0_sp) tcb_bot = pcoh*bed(i,Ksed+1,ibtcr)+(pcoh-1.0_sp)*cff4 !BF endif Ksed = Ksed+1 ENDDO frt = MAX(0.0_sp,(cff1-tcb_top) ) / & & MAX( eps, & & ( tcb_bot-tcb_top )) ENDIF bzactv = bzactv+frt*bed(i,Ksed,ithck) bzactv = MAX( bzactv, 6.0_sp*bottom(i,isd50) ) !CRS bottom(i,iactv)=min(bottom(i,iactv),bzactv) !?CRS ! bottom(i,iactv)=max(bottom(i,iactv),bzactv) !?CRS endif !/* defined COHESIVE_BED || defined MIXED_BED */ END DO I_LOOP_CB J_LOOP2 : DO i=1,m ! ! Calculate net deposition and erosion cff=0.0_sp cff2=0.0_sp DO ised=1,NST cff=cff+sedbed%settling_flux(i,ised) cff2=cff2+sedbed%ero_flux(i,ised) dep_mass(ised)=0.0_sp IF ((sedbed%ero_flux(i,ised)-sedbed%settling_flux(i,ised)).lt. & & 0.0_sp) THEN dep_mass(ised)=sedbed%settling_flux(i,ised)- & & sedbed%ero_flux(i,ised) END IF END DO bottom(i,idnet)=cff-cff2 IF ( cff-cff2.GT.0.0_sp) THEN ! NET deposition ! Deposition. Determine if we need to create a new bed layer if( COHESIVE_BED .or. MIXED_BED)then ! ! Calculate tau_crit of deposited bed ! ! Calculate new mass in top layer bmz(1) = 0.0_sp DO ised=1,NST bmz(1) = bmz(1)+sedbed%bed_mass(i,1,nnew,ised) END DO IF (Nbed.GT.1) THEN ! Average of cff1 and cff2, where ! cff1 = linear extension of previous tcr slope to new surface ! cff2 = minimum deposition cff1 = bed(i,1,ibtcr) - & & bottom(i,idnet) * & & (bed(i,2,ibtcr)-bed(i,1,ibtcr)) / & & (bmz(1)-bottom(i,idnet)) cff2 = MAX( tau_w(i) , tcr_min ) bed(i,1,ibtcr) = MIN(bed(i,1,ibtcr), & & MAX( 0.5_sp*(cff1+cff2), cff2 )) ! cff1 = bottom(i,idnet)/MAX(bmz(1),eps) ! cff2 = bed(i,1,ibtcr)+bed(i,2,ibtcr)-2*tcr_min ! bed(i,1,ibtcr) =bed(i,1,ibtcr) - cff1*cff2 ELSE ! TODO: Not sure what mud tau_crit should be for single-layer bed. ! Try weighted average of dep and old value cff1 = (bed(i,1,ibtcr)*(bmz(1)-bottom(i,idnet)) & & + cff2*bottom(i,idnet)) / bmz(1) bed(i,1,ibtcr) = MAX( cff1, cff2 ) END IF endif ! (no test for age here) bed(i,1,iaged)=T_model IF(bed(i,1,ithck).gt. & & MAX(bottom(i,iactv),newlayer_thick)) THEN ! Top layer is too thick IF (Nbed.gt.2) THEN IF(bed(i,2,ithck).lt.minlayer_thick) THEN ! Layer 2 is smaller than minimum size ! Instead of pushing down all layers, just combine top 2 layers cff=0.0_sp cff1=0.0_sp cff2=0.0_sp DO ised=1,NST cff =cff +dep_mass(ised) cff1=cff1+sedbed%bed_mass(i,1,nnew,ised) cff2=cff2+sedbed%bed_mass(i,2,nnew,ised) END DO if(COHESIVE_BED .or. MIXED_BED)then ! ! Assign new tau_crit at 2nd layer ! bed(i,2,ibtcr)= bed(i,2,ibtcr) - & & (cff1-cff)/(cff1+cff2)*(bed(i,3,ibtcr)-bed(i,1,ibtcr)) ! bed(i,2,ibtcr) = 0.5_sp*( bed(i,1,ibtcr)+ & ! & bed(i,2,ibtcr) ) endif ! ! Update bed mass DO ised=1,NST sedbed%bed_mass(i,2,nnew,ised)= & & MAX(sedbed%bed_mass(i,2,nnew,ised)+ & & sedbed%bed_mass(i,1,nnew,ised)- & & dep_mass(ised),0.0_sp) sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised) END DO ! ALA - average time and porosity ! ALA CHECK WITH CRS cff1 or cff1-cff for first layer bed(i,2,iaged)=(bed(i,1,iaged)*cff1+ & & bed(i,2,iaged)*cff2)/(cff1+cff2) bed(i,1,iaged)=T_model bed(i,2,iporo)=(bed(i,1,iporo)*cff1+ & & bed(i,2,iporo)*cff2)/(cff1+cff2) ! ALA CHECK WITH CRS POROSITY OF 1ST LAYER bed(i,1,iporo)=bed(i,1,iporo) ELSE ! Layer 2 is > minlayer thick, need another layer ! Combine bottom layers. cff1=0.0_sp cff2=0.0_sp DO ised=1,NST cff1=cff1+sedbed%bed_mass(i,Nbed-1,nnew,ised) cff2=cff2+sedbed%bed_mass(i,Nbed,nnew,ised) END DO bed(i,Nbed,iporo)= (bed(i,Nbed-1,iporo)*cff1+ & & bed(i,Nbed,iporo)*cff2)/(cff1+cff2) bed(i,Nbed,iaged)= (bed(i,Nbed-1,iaged)*cff1+ & & bed(i,Nbed,iaged)*cff2)/(cff1+cff2) if(COHESIVE_BED .or. MIXED_BED)then ! ! Assign tcrit at top of new bottom bed tcrit for Nbed-1 bed(i,Nbed,ibtcr)= bed(i,Nbed-1,ibtcr) ! bed(i,Nbed,ibtcr)=bed(i,Nbed-1,ibtcr) + & ! & cff2*(bed(i,Nbed,ibtcr)-bed(i,Nbed-1,ibtcr))/cff1 endif DO ised=1,NST sedbed%bed_mass(i,Nbed,nnew,ised)= & & sedbed%bed_mass(i,Nbed-1,nnew,ised)+ & & sedbed%bed_mass(i,Nbed ,nnew,ised) END DO ! ! Push layers down. DO k=Nbed-1,2,-1 bed(i,k,iporo)=bed(i,k-1,iporo) bed(i,k,iaged)=bed(i,k-1,iaged) DO ised =1,NST sedbed%bed_mass(i,k,nnew,ised)= & & sedbed%bed_mass(i,k-1,nnew,ised) END DO if(COHESIVE_BED .or. MIXED_BED)then bed(i,k,ibtcr)=bed(i,k-1,ibtcr) endif END DO ! Set new top parameters for top 2 layers if(COHESIVE_BED .or. MIXED_BED)then ! Tau_crit at top already determined ! Set tau_crit at 2nd layer cff=0.0_sp cff1=0.0_sp DO ised=1,NST cff =cff +dep_mass(ised) cff1=cff1+sedbed%bed_mass(i,1,nnew,ised) END DO cff2 = (bed(i,2,ibtcr)-bed(i,1,ibtcr))/cff1 bed(i,2,ibtcr) = bed(i,1,ibtcr)+cff*cff2 endif DO ised=1,NST sedbed%bed_mass(i,2,nnew,ised)= & & MAX(sedbed%bed_mass(i,2,nnew,ised)- & & dep_mass(ised),0.0_sp) sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised) END DO END IF ELSEIF (Nbed.eq.2) THEN ! NBED=2 ! if(COHESIVE_BED .or. MIXED_BED)then ! Tau_crit at top already determined ! Set tau_crit at 2nd layer cff=0.0_sp cff1=0.0_sp DO ised=1,NST cff =cff +dep_mass(ised) cff1=cff1+sedbed%bed_mass(i,1,nnew,ised) END DO cff2 = (bed(i,2,ibtcr)-bed(i,1,ibtcr))/cff1 bed(i,2,ibtcr) = bed(i,1,ibtcr)+cff*cff2 endif cff1=0.0_sp cff2=0.0_sp DO ised=1,NST cff1=cff1+sedbed%bed_mass(i,1,nnew,ised) cff2=cff2+sedbed%bed_mass(i,2,nnew,ised) END DO DO ised=1,NST sedbed%bed_mass(i,2,nnew,ised)= & & MAX(sedbed%bed_mass(i,2,nnew,ised)+ & & sedbed%bed_mass(i,1,nnew,ised)- & & dep_mass(ised),0.0_sp) sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised) END DO ! ALA - average time and porosity bed(i,2,iaged)=(bed(i,1,iaged)*cff1+ & & bed(i,2,iaged)*cff2)/(cff1+cff2) bed(i,1,iaged)=T_model bed(i,2,iporo)=(bed(i,1,iporo)*cff1+ & & bed(i,2,iporo)*cff2)/(cff1+cff2) ! ALA CHECK WITH CRS POROSITY OF 1ST LAYER bed(i,1,iporo)=bed(i,1,iporo) ELSE ! NBED=1 END IF ELSE ! Net deposition has occurred, but no new bed layer was created END IF ELSE ! Net erosion occurred if(COHESIVE_BED .or. MIXED_BED)then bmz(1) = 0.0_sp DO ised=1,NST bmz(1) = bmz(1)+sedbed%bed_mass(i,1,nnew,ised) END DO ! recalc tc for top of bed, based on linear ! interpolation and mass removed / orig. mass in top layer bed(i,1,ibtcr)=bed(i,1,ibtcr)+ & & MIN(1.0_sp,-bottom(i,idnet)/MAX(eps,bmz(1)))* & & MAX(0.0_sp,(bed(i,2,ibtcr)-bed(i,1,ibtcr))) endif bed(i,1,iaged)=T_model IF (Nbed.eq.2) THEN ! NBED=2 DO ised=1,NST sedbed%bed_mass(i,2,nnew,ised)= & & MAX(sedbed%bed_mass(i,2,nnew,ised)+ & & sedbed%bed_mass(i,1,nnew,ised)- & & dep_mass(ised),0.0_sp) sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised) END DO ELSEIF (Nbed.eq.1) THEN ! ALF NO NEED TO DO ANYTHING ELSE END IF END IF ! Recalculate thickness and fractions for all layers. DO k=1,Nbed cff3=0.0_sp DO ised=1,NST cff3=cff3+sedbed%bed_mass(i,k,nnew,ised) END DO IF (cff3.eq.0.0_sp) THEN cff3=eps END IF bed(i,k,ithck)=0.0_sp DO ised=1,NST sedbed%bed_frac(i,k,ised)=sedbed%bed_mass(i,k,nnew,ised)/cff3 bed(i,k,ithck)=MAX(bed(i,k,ithck)+ & & sedbed%bed_mass(i,k,nnew,ised)/ & & (sed(ised)%Srho* & & (1.0_sp-bed(i,k,iporo))),0.0_sp) END DO END DO END DO J_LOOP2 J_LOOP3 : DO i=1,m IF (bottom(i,iactv).gt.bed(i,1,ithck)) THEN IF (Nbed.eq.1) THEN bottom(i,iactv)=bed(i,1,ithck) ELSE thck_to_add=bottom(i,iactv)-bed(i,1,ithck) thck_avail=0.0_sp Ksed=1 ! initialize DO k=2,Nbed IF (thck_avail.lt.thck_to_add) THEN thck_avail=thck_avail+bed(i,k,ithck) Ksed=k END IF END DO ! ! Catch here if there was not enough bed material. ! IF (thck_avail.lt.thck_to_add) THEN bottom(i,iactv)=bed(i,1,ithck)+thck_avail thck_to_add=thck_avail END IF ! ! Update bed mass of top layer and fractional layer. ! cff2=MAX(thck_avail-thck_to_add,0.0_sp)/ & & MAX(bed(i,Ksed,ithck),eps) DO ised=1,NST cff1=0.0_sp DO k=1,Ksed cff1=cff1+sedbed%bed_mass(i,k,nnew,ised) END DO cff3=cff2*sedbed%bed_mass(i,Ksed,nnew,ised) sedbed%bed_mass(i,1 ,nnew,ised)=cff1-cff3 sedbed%bed_mass(i,Ksed,nnew,ised)=cff3 END DO ! ! Update thickness of fractional layer ksource_sed. ! bed(i,Ksed,ithck)=MAX(thck_avail-thck_to_add,0.0_sp) if(COHESIVE_BED .or. MIXED_BED)then ! ! Update tau_cr of fractional layer if(Nbed>2)then bed(i,Ksed,ibtcr) = bed(i,Ksed+1,ibtcr)- & & cff2*(bed(i,Ksed+1,ibtcr)-bed(i,Ksed,ibtcr)) end if endif ! ! Update bed fraction of top layer. ! cff3=0.0_sp DO ised=1,NST cff3=cff3+sedbed%bed_mass(i,1,nnew,ised) END DO IF (cff3.eq.0.0_sp) THEN cff3=eps END IF DO ised=1,NST sedbed%bed_frac(i,1,ised)=sedbed%bed_mass(i,1,nnew,ised)/cff3 END DO ! ! Upate bed thickness of top layer. ! bed(i,1,ithck)=bottom(i,iactv) ! ! Pull all layers closer to the surface. ! DO k=Ksed,Nbed ks=Ksed-2 bed(i,k-ks,ithck)=bed(i,k,ithck) bed(i,k-ks,iporo)=bed(i,k,iporo) bed(i,k-ks,iaged)=bed(i,k,iaged) if(COHESIVE_BED .or. MIXED_BED)then bed(i,k-ks,ibtcr)=bed(i,k,ibtcr) endif DO ised=1,NST sedbed%bed_frac(i,k-ks,ised)=sedbed%bed_frac(i,k,ised) sedbed%bed_mass(i,k-ks,nnew,ised)=sedbed%bed_mass(i,k,nnew,ised) END DO END DO ! ! Add new layers onto the bottom. Split what was in the bottom layer to ! fill these new empty cells. ("ks" is the number of new layers). ! ks=Ksed-2 cff=1.0_sp/REAL(ks+1,sp) if (COHESIVE_BED .or. MIXED_BED)then !alf cff1 = cff*(tcr_max-bed(i,Nbed-ks,ibtcr)) !alf Use the value from the bottom layer as max. cff1 = cff*(bed(i,Nbed,ibtcr)-bed(i,Nbed-ks,ibtcr)) ! Interpolate bottom layer to tau_crit_max ! TODO: should be the reference profile and not tau_crit_max DO k=Nbed-1,Nbed-ks,-1 bed(i,k,ibtcr)=bed(i,Nbed-ks,ibtcr)+ & & REAL(k-Nbed+ks,sp) * cff1 ENDDO endif ! ALA CHECK WITH CRS about bed_frac nnn=0 DO ised=1,NST nlysm(ised)=newlayer_thick*REAL(ks+1,sp)* & & (sed(ised)%Srho* & & (1.0_sp-bed(i,Nbed-ks,iporo)))* & & sedbed%bed_frac(i,Nbed-ks,ised) IF (ks.gt.0) THEN IF (sedbed%bed_mass(i,Nbed-ks,nnew,ised).gt. & & nlysm(ised)) THEN nnn=nnn+1 nlysm(ised)= newlayer_thick*REAL(ks,sp)* & & (sed(ised)%Srho* & & (1.0_sp-bed(i,Nbed-ks,iporo)))* & & sedbed%bed_frac(i,Nbed-ks,ised) END IF END IF END DO IF (nnn.eq.NST) THEN bed(i,Nbed,ithck)=bed(i,Nbed-ks,ithck)- & & newlayer_thick*REAL(ks,sp) DO ised=1,NST sedbed%bed_mass(i,Nbed,nnew,ised)= & & sedbed%bed_mass(i,Nbed-ks,nnew,ised)-nlysm(ised) END DO DO k=Nbed-1,Nbed-ks,-1 bed(i,k,ithck)=newlayer_thick bed(i,k,iaged)=bed(i,Nbed-ks,iaged) DO ised=1,NST sedbed%bed_frac(i,k,ised)=sedbed%bed_frac(i,Nbed-ks,ised) sedbed%bed_mass(i,k,nnew,ised)= & & nlysm(ised)/REAL(ks,sp) END DO END DO ELSE cff=1.0_sp/REAL(ks+1,sp) DO k=Nbed,Nbed-ks,-1 bed(i,k,ithck)=bed(i,Nbed-ks,ithck)*cff bed(i,k,iaged)=bed(i,Nbed-ks,iaged) DO ised=1,NST sedbed%bed_frac(i,k,ised)=sedbed%bed_frac(i,Nbed-ks,ised) sedbed%bed_mass(i,k,nnew,ised)= & & sedbed%bed_mass(i,Nbed-ks,nnew,ised)*cff END DO END DO END IF END IF ! Nbed > 1 END IF ! increase top bed layer if(MIXED_BED)then ! ! Update cohesive property of seds in top layer cff1 = 0.0_sp DO ised=1,NCS cff1=cff1+sedbed%bed_frac(i,1,ised) END DO bottom(i,idprp)=min(max((cff1-transN)/ & & (transC-transN),0.0_sp),1.0_sp) endif if(SED_FLOCS .and. SED_DEFLOC)then alpha = MIN(DTsed/t_dfloc,1.0_sp) DO k=2,Nbed cff3=0.0_sp DO ised=1,NST cff3=cff3+sedbed%bed_mass(i,k,nnew,ised) END DO cff1 = 0.0_sp DO ised=1,NCS cff1 = cff1+sedbed%bed_mass(i,k,nnew,ised) END DO DO ised=1,NCS masseq(ised)=mud_frac_eq(ised)*cff1 sedbed%bed_mass(i,k,nnew,ised)=alpha*masseq(ised)+ & & sedbed%bed_mass(i,k,nnew,ised)*(1.0_sp-alpha) sedbed%bed_frac(i,k,ised)=sedbed%bed_mass(i,k,nnew,ised)/cff3 END DO END DO endif if(COHESIVE_BED .or. MIXED_BED)then ! ! Key to this algorithm: "mass available depth" (mad) [kg/m2] ! present tau crit profile (tcp) ! representative tau crit profile (tcr) ! Compute depth and mass depth of sediment column ! Compute cumulative depth ! Compute mass depth bmz(1) = 0.0_sp DO ised=1,NST bmz(1) = bmz(1)+sedbed%bed_mass(i,1,nnew,ised) ENDDO DO k=2,Nbed bmz(k) = bmz(k-1) DO ised=1,NST bmz(k)=bmz(k)+sedbed%bed_mass(i,k,nnew,ised) ENDDO ENDDO ! Calculate representative critical shear stress profile ! Note that the values are for the TOP of the layer...we ! assume the bottom of the bottom layer has tcr = tcr_max tcr(1) = tcr_min DO k=2,Nbed tcr(k) = tcr_min IF (bmz(k-1).GT.eps) THEN ! tcr(k) = exp((log(bmz(k-1))- & ! & bottom(i,idoff))/ & ! & bottom(i,idslp)) tcr(k) = exp((log(bmz(k-1))- & & tcr_off)/ & & tcr_slp) ENDIF tcr(k) = MIN( MAX( tcr(k), tcr_min), tcr_max ) ENDDO ! ks=Ksed-2 ! DO k=Nbed,Nbed-ks,-1 ! bed(i,k,ibtcr)=tcr(k) ! ENDDO ! Relax tau crit profile bottom(i,k,ibtcr) ! towards representative profile tcr...100 x slower for "swelling" ! than consolidation ! TODO - make the factor an input parameter IF( bed(i,1,ibtcr).LE.tcr(1)) THEN ! alpha = MIN(DTsed/bottom(i,idtim),1.0_sp) alpha = MIN(DTsed/tcr_tim,1.0_sp) ELSE ! alpha = MIN(DTsed/(100.0_sp*bottom(i,idtim)),1.0_sp) alpha = MIN(DTsed/(100.0_sp*tcr_tim),1.0_sp) ENDIF DO k=1,Nbed bed(i,k,ibtcr)=bed(i,k,ibtcr)+ & & alpha*(tcr(k)-bed(i,k,ibtcr)) bed(i,k,ibtcr)= & & MIN( MAX( bed(i,k,ibtcr), tcr_min), tcr_max ) ENDDO !ALA Enforce last layer = tcr_max !ALA bed(i,Nbed,ibtcr)= tcr_max endif !/* cohesive bed */ END DO J_LOOP3 ! !----------------------------------------------------------------------- ! Store old bed thickness. !----------------------------------------------------------------------- ! if (SED_MORPH)then do i=1,m sedbed%bed_thick(i,nnew)=0.0_sp DO k=1,Nbed sedbed%bed_thick(i,nnew)=sedbed%bed_thick(i,nnew)+ & & bed(i,k,ithck) END DO end do # if defined (MULTIPROCESSOR) if(par)then call aexchange(nc,myid,nprocs,sedbed%bed_thick) endif # endif bottom(i,morph) = sedbed%bed_thick(i,nnew) endif # if defined (MULTIPROCESSOR) if(par)then allocate(temp(0:mt,1:nbed)) DO ised=1,NST temp = 0.0_sp temp = sedbed%bed_frac(:,:,ised) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) sedbed%bed_frac(:,:,ised) = temp temp = 0.0_sp temp = sedbed%bed_mass(:,:,nnew,ised) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) sedbed%bed_mass(:,:,nnew,ised) = temp end do do ibed=1,MBEDP temp = 0.0_sp temp = bed(:,:,ibed) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) bed(:,:,ibed)=temp end do deallocate(temp) endif # endif if(dbg_set(dbg_sbr)) write(ipt,*) "End: calc_sed_bed_cohesive" return end Subroutine calc_sed_bed_cohesive !======================================================================= ! ! ! This routine computes sediment bed layer stratigraphy. ! ! ! ! Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G. ! ! Arango, 2008: Development of a three-dimensional, regional, ! ! coupled wave, current, and sediment-transport model, Computers ! ! & Geosciences, 34, 1284-1306. ! ! ! !======================================================================= Subroutine calc_sed_bed_cohesive2 implicit none ! ! Local variable declarations. ! integer :: Kbed, i, ised,ibed, j, k, ks integer :: bnew, nnn, cnt,cc,jj real(sp), parameter :: eps = 1.0E-14_sp real(sp) :: cff, cff1, cff2, cff3,cff4 real(sp) :: thck_avail, thck_to_add real(sp), dimension(NST) :: nlysm real(sp), dimension(NST) :: dep_mass real(sp), dimension(0:mt) :: tau_w real(sp) :: pcoh real(sp) :: alpha, bzactv, frt, tcb_top, tcb_bot real(sp), dimension(Nbed) :: tcr real(sp), dimension(Nbed) :: bmz real(sp), dimension(NCS) :: masseq real(sp),allocatable,dimension(:,:)::temp if(dbg_set(dbg_sbr)) write(ipt,*) "Start: calc_sed_bed_cohesive" if(BEDLOAD)then bnew=nnew else bnew=nstp endif ! KLUDGE alert ! minlayer_thick = 0.0005 minlayer_thick = newlayer_thick ! !----------------------------------------------------------------------- ! Compute sediment bed layer stratigraphy. !----------------------------------------------------------------------- ! if(BEDLOAD_MPM.or.SUSLOAD)then tau_w = taub endif ! !----------------------------------------------------------------------- ! Update bed properties according to ero_flux and dep_flux. !----------------------------------------------------------------------- ! if(SUSLOAD)then NODE_LOOP : DO i=1,m ! ! The deposition and resuspension of sediment on the bottom "bed" ! is due to precipitation flux FC(:,0), already computed, and the ! resuspension (erosion, hence called ero_flux). The resuspension is ! applied to the bottom-most grid box value qc(:,1) so the total mass ! is conserved. Restrict "ero_flux" so that "bed" cannot go negative ! after both fluxes are applied. ! SED_LOOP: DO ised=1,NST dep_mass(ised)=0.0_sp if (SED_MORPH)then ! Apply morphology factor. sedbed%ero_flux(i,ised)=sedbed%ero_flux(i,ised)*sed(ised)%morph_fac sedbed%settling_flux(i,ised)=sedbed%settling_flux(i,ised)* & & sed(ised)%morph_fac endif ! Update bed mass arrays. sedbed%bed_mass(i,1,nnew,ised)=MAX(sedbed%bed_mass(i,1,bnew,ised)- & & (sedbed%ero_flux(i,ised)- & & sedbed%settling_flux(i,ised)), & & 0.0_sp) DO k=2,Nbed sedbed%bed_mass(i,k,nnew,ised)=sedbed%bed_mass(i,k,nstp,ised) END DO END DO SED_LOOP cff3=0.0_sp DO ised=1,NST cff3=cff3+sedbed%bed_mass(i,1,nnew,ised) END DO IF (cff3.eq.0.0_sp) THEN cff3=eps END IF bed(i,1,ithck)=0.0_sp DO ised=1,NST sedbed%bed_frac(i,1,ised)=sedbed%bed_mass(i,1,nnew,ised)/cff3 bed(i,1,ithck)=MAX(bed(i,1,ithck)+ & & sedbed%bed_mass(i,1,nnew,ised)/ & & (sed(ised)%Srho* & & (1.0_sp-bed(i,1,iporo))),0.0_sp) END DO END DO NODE_LOOP endif !/* SUSPLOAD section */ ! !----------------------------------------------------------------------- ! At this point, all deposition or erosion is complete, and ! has been added/subtracted to top layer. Thickness has NOT been corrected. !----------------------------------------------------------------------- ! NODE_LOOP_CB : DO i=1,m ! Calculate active layer thickness, bottom(i,iactv). ! (trunk version allows this to be zero...this has minimum of 6*D50) bottom(i,iactv)=MAX(0.0_sp, & & 0.007_sp* & & (tau_w(i)-bottom(i,itauc)))+ & & 6.0_sp*bottom(i,isd50) if (COHESIVE_BED)then bottom(i,iactv)=MAX(0.0_sp, & & 0.007_sp* & & (tau_w(i)-bed(i,1,ibtcr)))+ & & 6.0_sp*bottom(i,isd50) endif if (MIXED_BED)then cff1= MAX(0.0_sp, & & 0.007_sp* & & (tau_w(i)-bed(i,1,ibtcr)))+ & & 6.0_sp*bottom(i,isd50) cff2= MAX(0.0_sp, & & 0.007_sp* & & (tau_w(i)-bottom(i,itauc)))+ & & 6.0_sp*bottom(i,isd50) bottom(i,iactv)=MAX(cff1,cff2) endif if (SED_MORPH)then ! Apply morphology factor. bottom(i,iactv)=MAX(bottom(i,iactv)*sed(1)%morph_fac, & & bottom(i,iactv)) endif if (COHESIVE_BED .or. MIXED_BED)then ! Find first layer with tc > tb ! Remember, the critical stresses apply to the TOP of each layer Kbed = 1 bzactv = 0.0_sp frt = 0.0_sp ! CRS cff1 = tau_w(i) tcb_top = bed(i,1,ibtcr) tcb_bot = bed(i,2,ibtcr) if (MIXED_BED)then ! Calc. tau crit for bottom of next layer ! Update cohesive property of seds in top layer cff3 = 0.0_sp cff4 = 1.0_sp DO ised=1,NCS cff3=cff3+sedbed%bed_frac(i,1,ised) cff4=cff4*sed(ised)%tau_ce**sedbed%bed_frac(i,1,ised) END DO DO ised=NCS+1,NST cff4=cff4*sed(ised)%tau_ce**sedbed%bed_frac(i,1,ised) ENDDO pcoh=min(max((cff3-transN)/(transC-transN) & & ,0.0_sp),1.0_sp) tcb_top = pcoh*bed(i,1,ibtcr)+(pcoh-1.0_sp)*cff4 !BF endif IF(cff1 .GT. tcb_top)THEN ! Calculate tcb_temp for next layer tcb_bot = bed(i,Kbed+1,ibtcr) if (MIXED_BED)then ! Recalculate cohesive fraction and mean grain tau crit ! Note that Kbed is used for grain props, and Kbed+1 for ! bed tau crit at bottom of layer Kbed cff3 = 0.0_sp cff4 = 1.0_sp DO ised=1,NCS cff3=cff3+sedbed%bed_frac(i,Kbed,ised) cff4=cff4*sed(ised)%tau_ce**sedbed%bed_frac(i,Kbed,ised) END DO DO ised=NCS+1,NST cff4=cff4*sed(ised)%tau_ce**sedbed%bed_frac(i,Kbed,ised) ENDDO ! Calculate cohesive behavior and blended tau crit pcoh=min(max((cff3-transN)/(transC-transN), & & 0.0_sp),1.0_sp) tcb_bot =pcoh*bed(i,Kbed+1,ibtcr)+(pcoh-1.0_sp)*cff4 !BF endif DO WHILE ( (Kbed.LE.(Nbed-1)) .AND. & & (cff1 .GT. tcb_bot)) !ALA Instead of adding entire 2nd layer, add just what is needed !ALA! Add entire layer !ALA bzactv = bzactv + bed(i,Kbed,ithck) ! Add thickness equivalent to difference between cff1 and tcb_bot IF (cff1.GT.bed(i,Kbed+1,ibtcr)) THEN bzactv = bzactv + bed(i,Kbed,ithck) tcb_top = tcb_bot tcb_bot = bed(i,Kbed+1,ibtcr) ELSE bzactv = bzactv + MIN(bed(i,Kbed,ithck), & & (MAX(0.0_sp,0.007_sp*(cff1-tcb_bot)))) tcb_top = cff1 tcb_bot = bed(i,Kbed+1,ibtcr) ENDIF if (MIXED_BED)then ! Recalculate cohesive fraction and mean grain tau crit cff3 = 0.0_sp cff4 = 1.0_sp DO ised=1,NCS cff3=cff3+sedbed%bed_frac(i,Kbed,ised) cff4=cff4*sed(ised)%tau_ce**sedbed%bed_frac(i,Kbed,ised) END DO DO ised=NCS+1,NST cff4=cff4*sed(ised)%tau_ce**sedbed%bed_frac(i,Kbed,ised) ENDDO ! Calculate cohesive behavior and blended tau crit pcoh=min(max((cff3-transN)/(transC-transN), & & 0.0_sp),1.0_sp) tcb_bot = pcoh*bed(i,Kbed+1,ibtcr)+(pcoh-1.0_sp)*cff4 endif Kbed = Kbed+1 ENDDO frt = MAX(0.0_sp,(cff1-tcb_top) ) / & & MAX( eps, & & ( tcb_bot-tcb_top )) ENDIF bzactv = bzactv+frt*bed(i,Kbed,ithck) bzactv = MAX( bzactv, 6.0_sp*bottom(i,isd50) ) !CRS bottom(i,iactv)=min(bottom(i,iactv),bzactv) !?CRS ! bottom(i,iactv)=max(bottom(i,iactv),bzactv) !?CRS endif !/* defined COHESIVE_BED || defined MIXED_BED */ END DO NODE_LOOP_CB NODE_LOOP2 : DO i=1,m ! ! Calculate net deposition and erosion cff=0.0_sp cff2=0.0_sp DO ised=1,NST cff=cff+sedbed%settling_flux(i,ised) cff2=cff2+sedbed%ero_flux(i,ised) dep_mass(ised)=0.0_sp IF ((sedbed%ero_flux(i,ised)-sedbed%settling_flux(i,ised)).lt. & & 0.0_sp) THEN dep_mass(ised)=sedbed%settling_flux(i,ised)- & & sedbed%ero_flux(i,ised) END IF END DO bottom(i,idnet)=cff-cff2 IF ( cff-cff2.GT.0.0_sp) THEN ! NET deposition ! Deposition. Determine if we need to create a new bed layer if (COHESIVE_BED .or. MIXED_BED)then ! ! Calculate tau_crit of deposited bed ! ! Calculate new mass in top layer bmz(1) = 0.0_sp DO ised=1,NST bmz(1) = bmz(1)+sedbed%bed_mass(i,1,nnew,ised) END DO IF (Nbed.GT.1) THEN ! Average of cff1 and cff2, where ! cff1 = linear extension of previous tcr slope to new surface ! cff2 = minimum deposition cff1 = bed(i,1,ibtcr) - & & bottom(i,idnet) * & & (bed(i,2,ibtcr)-bed(i,1,ibtcr)) / & & (bmz(1)-bottom(i,idnet)) cff2 = MAX( tau_w(i) , tcr_min ) bed(i,1,ibtcr) = MIN(bed(i,1,ibtcr), & & MAX( 0.5_sp*(cff1+cff2), cff2 )) ! cff1 = bottom(i,idnet)/MAX(bmz(1),eps) ! cff2 = bed(i,1,ibtcr)+bed(i,2,ibtcr)-2*tcr_min ! bed(i,1,ibtcr) =bed(i,1,ibtcr) - cff1*cff2 ELSE ! TODO: Not sure what mud tau_crit should be for single-layer bed. ! Try weighted average of dep and old value cff1 = (bed(i,1,ibtcr)*(bmz(1)-bottom(i,idnet)) & & + cff2*bottom(i,idnet)) / bmz(1) cff2 = MAX( tau_w(i) , tcr_min ) bed(i,1,ibtcr) = MAX( cff1, cff2 ) END IF bed(i,1,ibtcr) = 0.1 !neeed be removed!! endif ! (no test for age here) bed(i,1,iaged)=T_model IF(bed(i,1,ithck).gt. & & MAX(bottom(i,iactv),newlayer_thick)) THEN ! Top layer is too thick IF (Nbed.gt.2) THEN IF(bed(i,2,ithck).lt.minlayer_thick) THEN ! Layer 2 is smaller than minimum size ! Instead of pushing down all layers, just combine top 2 layers cff=0.0_sp cff1=0.0_sp cff2=0.0_sp DO ised=1,NST cff =cff +dep_mass(ised) cff1=cff1+sedbed%bed_mass(i,1,nnew,ised) cff2=cff2+sedbed%bed_mass(i,2,nnew,ised) END DO if (COHESIVE_BED .or. MIXED_BED)then ! ! Assign new tau_crit at 2nd layer ! bed(i,2,ibtcr)= bed(i,2,ibtcr) - & & (cff1-cff)/(cff1+cff2)*(bed(i,3,ibtcr)-bed(i,1,ibtcr)) ! bed(i,2,ibtcr) = 0.5_sp*( bed(i,1,ibtcr)+ & ! & bed(i,2,ibtcr) ) endif ! ! Update bed mass DO ised=1,NST sedbed%bed_mass(i,2,nnew,ised)= & & MAX(sedbed%bed_mass(i,2,nnew,ised)+ & & sedbed%bed_mass(i,1,nnew,ised)- & & dep_mass(ised),0.0_sp) sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised) END DO ! ALA - average time and porosity ! ALA CHECK WITH CRS cff1 or cff1-cff for first layer bed(i,2,iaged)=(bed(i,1,iaged)*cff1+ & & bed(i,2,iaged)*cff2)/(cff1+cff2) bed(i,1,iaged)=T_model bed(i,2,iporo)=(bed(i,1,iporo)*cff1+ & & bed(i,2,iporo)*cff2)/(cff1+cff2) ! ALA CHECK WITH CRS POROSITY OF 1ST LAYER bed(i,1,iporo)=bed(i,1,iporo) ELSE ! Layer 2 is > minlayer thick, need another layer ! Combine bottom layers. cff1=0.0_sp cff2=0.0_sp DO ised=1,NST cff1=cff1+sedbed%bed_mass(i,Nbed-1,nnew,ised) cff2=cff2+sedbed%bed_mass(i,Nbed,nnew,ised) END DO bed(i,Nbed,iporo)= & & (bed(i,Nbed-1,iporo)*cff1+ & & bed(i,Nbed,iporo)*cff2)/(cff1+cff2) bed(i,Nbed,iaged)= & & (bed(i,Nbed-1,iaged)*cff1+ & & bed(i,Nbed,iaged)*cff2)/(cff1+cff2) if (COHESIVE_BED.or.MIXED_BED)then ! ! Assign tcrit at top of new bottom bed tcrit for Nbed-1 bed(i,Nbed,ibtcr)= bed(i,Nbed-1,ibtcr) ! bed(i,Nbed,ibtcr)=bed(i,Nbed-1,ibtcr) + & ! & cff2*(bed(i,Nbed,ibtcr)-bed(i,Nbed-1,ibtcr))/cff1 endif DO ised=1,NST sedbed%bed_mass(i,Nbed,nnew,ised)= & & sedbed%bed_mass(i,Nbed-1,nnew,ised)+ & & sedbed%bed_mass(i,Nbed ,nnew,ised) END DO ! ! Push layers down. DO k=Nbed-1,2,-1 bed(i,k,iporo)=bed(i,k-1,iporo) bed(i,k,iaged)=bed(i,k-1,iaged) DO ised =1,NST sedbed%bed_mass(i,k,nnew,ised)= & & sedbed%bed_mass(i,k-1,nnew,ised) END DO if (COHESIVE_BED.or.MIXED_BED)then bed(i,k,ibtcr)=bed(i,k-1,ibtcr) endif END DO ! Set new top parameters for top 2 layers if (COHESIVE_BED.or.MIXED_BED)then ! Tau_crit at top already determined ! Set tau_crit at 2nd layer cff=0.0_sp cff1=0.0_sp DO ised=1,NST cff =cff +dep_mass(ised) cff1=cff1+sedbed%bed_mass(i,1,nnew,ised) END DO cff2 = (bed(i,2,ibtcr)-bed(i,1,ibtcr))/cff1 bed(i,2,ibtcr) = bed(i,1,ibtcr)+cff*cff2 endif DO ised=1,NST sedbed%bed_mass(i,2,nnew,ised)= & & MAX(sedbed%bed_mass(i,2,nnew,ised)- & & dep_mass(ised),0.0_sp) sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised) END DO END IF ELSEIF (Nbed.eq.2) THEN ! ! NBED=2 ! if (COHESIVE_BED .and. MIXED_BED)then ! Tau_crit at top already determined ! Set tau_crit at 2nd layer cff=0.0_sp cff1=0.0_sp DO ised=1,NST cff =cff +dep_mass(ised) cff1=cff1+sedbed%bed_mass(i,1,nnew,ised) END DO cff2 = (bed(i,2,ibtcr)-bed(i,1,ibtcr))/cff1 bed(i,2,ibtcr) = bed(i,1,ibtcr)+cff*cff2 endif cff1=0.0_sp cff2=0.0_sp DO ised=1,NST cff1=cff1+sedbed%bed_mass(i,1,nnew,ised) cff2=cff2+sedbed%bed_mass(i,2,nnew,ised) END DO DO ised=1,NST sedbed%bed_mass(i,2,nnew,ised)= & & MAX(sedbed%bed_mass(i,2,nnew,ised)+ & & sedbed%bed_mass(i,1,nnew,ised)- & & dep_mass(ised),0.0_sp) sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised) END DO ! ALA - average time and porosity bed(i,2,iaged)=(bed(i,1,iaged)*cff1+ & & bed(i,2,iaged)*cff2)/(cff1+cff2) bed(i,1,iaged)=T_model bed(i,2,iporo)=(bed(i,1,iporo)*cff1+ & & bed(i,2,iporo)*cff2)/(cff1+cff2) ! ALA CHECK WITH CRS POROSITY OF 1ST LAYER bed(i,1,iporo)=bed(i,1,iporo) ELSE ! NBED=1 END IF ELSE ! Net deposition has occurred, but no new bed layer was created END IF ELSE ! Net erosion occurred if (COHESIVE_BED .or. MIXED_BED)then IF (Nbed.ge.2) THEN bmz(1) = 0.0_sp DO ised=1,NST bmz(1) = bmz(1)+sedbed%bed_mass(i,1,nnew,ised) END DO ! recalc tc for top of bed, based on linear ! interpolation and mass removed / orig. mass in top layer bed(i,1,ibtcr)=bed(i,1,ibtcr)+ & & MIN(1.0_sp,-bottom(i,idnet)/MAX(eps,bmz(1)))* & & MAX(0.0_sp,(bed(i,2,ibtcr)-bed(i,1,ibtcr))) endif endif bed(i,1,iaged)=T_model IF (Nbed.eq.2) THEN ! NBED=2 DO ised=1,NST sedbed%bed_mass(i,2,nnew,ised)= & & MAX(sedbed%bed_mass(i,2,nnew,ised)+ & & sedbed%bed_mass(i,1,nnew,ised)- & & dep_mass(ised),0.0_sp) sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised) END DO ELSEIF (Nbed.eq.1) THEN ! ALF NO NEED TO DO ANYTHING ELSE END IF END IF ! Recalculate thickness and fractions for all layers. DO k=1,Nbed cff3=0.0_sp DO ised=1,NST cff3=cff3+sedbed%bed_mass(i,k,nnew,ised) END DO IF (cff3.eq.0.0_sp) THEN cff3=eps END IF bed(i,k,ithck)=0.0_sp DO ised=1,NST sedbed%bed_frac(i,k,ised)=sedbed%bed_mass(i,k,nnew,ised)/cff3 bed(i,k,ithck)=MAX(bed(i,k,ithck)+ & & sedbed%bed_mass(i,k,nnew,ised)/ & & (sed(ised)%Srho* & & (1.0_sp-bed(i,k,iporo))),0.0_sp) END DO END DO END DO NODE_LOOP2 NODE_LOOP3 : DO i=1,m IF (bottom(i,iactv).gt.bed(i,1,ithck)) THEN IF (Nbed.eq.1) THEN bottom(i,iactv)=bed(i,1,ithck) ELSE thck_to_add=bottom(i,iactv)-bed(i,1,ithck) thck_avail=0.0_sp Kbed=1 ! initialize DO k=2,Nbed IF (thck_avail.lt.thck_to_add) THEN thck_avail=thck_avail+bed(i,k,ithck) Kbed=k END IF END DO ! ! Catch here if there was not enough bed material. ! IF (thck_avail.lt.thck_to_add) THEN bottom(i,iactv)=bed(i,1,ithck)+thck_avail thck_to_add=thck_avail END IF ! ! Update bed mass of top layer and fractional layer. ! cff2=MAX(thck_avail-thck_to_add,0.0_sp)/ & & MAX(bed(i,Kbed,ithck),eps) DO ised=1,NST cff1=0.0_sp DO k=1,Kbed cff1=cff1+sedbed%bed_mass(i,k,nnew,ised) END DO cff3=cff2*sedbed%bed_mass(i,Kbed,nnew,ised) sedbed%bed_mass(i,1 ,nnew,ised)=cff1-cff3 sedbed%bed_mass(i,Kbed,nnew,ised)=cff3 END DO ! ! Update thickness of fractional layer ksource_sed. ! bed(i,Kbed,ithck)=MAX(thck_avail-thck_to_add,0.0_sp) if (COHESIVE_BED .or. MIXED_BED)then ! ! Update tau_cr of fractional layer if(Kbed 1 END IF ! increase top bed layer if (MIXED_BED)then ! ! Update cohesive property of seds in top layer cff1 = 0.0_sp DO ised=1,NCS cff1=cff1+sedbed%bed_frac(i,1,ised) END DO bottom(i,idprp)=min(max((cff1-transN)/ & & (transC-transN),0.0_sp),1.0_sp) endif if (SED_FLOCS .and. SED_DEFLOC)then alpha = MIN(DTsed/t_dfloc,1.0_sp) DO k=2,Nbed cff3=0.0_sp DO ised=1,NST cff3=cff3+sedbed%bed_mass(i,k,nnew,ised) END DO cff1 = 0.0_sp DO ised=1,NCS cff1 = cff1+sedbed%bed_mass(i,k,nnew,ised) END DO DO ised=1,NCS masseq(ised)=mud_frac_eq(ised)*cff1 sedbed%bed_mass(i,k,nnew,ised)=alpha*masseq(ised)+ & & sedbed%bed_mass(i,k,nnew,ised)*(1.0_sp-alpha) sedbed%bed_frac(i,k,ised)=sedbed%bed_mass(i,k,nnew,ised)/cff3 END DO END DO endif if (COHESIVE_BED.or. MIXED_BED)then ! ! Key to this algorithm: "mass available depth" (mad) [kg/m2] ! present tau crit profile (tcp) ! representative tau crit profile (tcr) ! Compute depth and mass depth of sediment column ! Compute cumulative depth ! Compute mass depth bmz(1) = 0.0_sp DO ised=1,NST bmz(1) = bmz(1)+sedbed%bed_mass(i,1,nnew,ised) ENDDO DO k=2,Nbed bmz(k) = bmz(k-1) DO ised=1,NST bmz(k)=bmz(k)+sedbed%bed_mass(i,k,nnew,ised) ENDDO ENDDO ! Calculate representative critical shear stress profile ! Note that the values are for the TOP of the layer...we ! assume the bottom of the bottom layer has tcr = tcr_max tcr(1) = tcr_min DO k=2,Nbed tcr(k) = tcr_min IF (bmz(k-1).GT.eps) THEN ! tcr(k) = exp((log(bmz(k-1))- & ! & bottom(i,idoff))/ & ! & bottom(i,idslp)) tcr(k) = exp((log(bmz(k-1))- & & tcr_off)/ & & tcr_slp) ENDIF tcr(k) = MIN( MAX( tcr(k), tcr_min), tcr_max ) ENDDO ! ks=Kbed-2 ! DO k=Nbed,Nbed-ks,-1 ! bed(i,k,ibtcr)=tcr(k) ! ENDDO ! Relax tau crit profile bottom(i,k,ibtcr) ! towards representative profile tcr...100 x slower for "swelling" ! than consolidation ! TODO - make the factor an input parameter IF( bed(i,1,ibtcr).LE.tcr(1)) THEN ! alpha = MIN(DTsed/bottom(i,idtim),1.0_sp) alpha = MIN(DTsed/tcr_tim,1.0_sp) ELSE ! alpha = MIN(DTsed/(100.0_sp*bottom(i,idtim)),1.0_sp) alpha = MIN(DTsed/(100.0_sp*tcr_tim),1.0_sp) ENDIF DO k=1,Nbed bed(i,k,ibtcr)=bed(i,k,ibtcr)+ & & alpha*(tcr(k)-bed(i,k,ibtcr)) bed(i,k,ibtcr)= & & MIN( MAX( bed(i,k,ibtcr), tcr_min), tcr_max ) ENDDO !ALA Enforce last layer = tcr_max !ALA bed(i,Nbed,ibtcr)= tcr_max endif !/* cohesive bed */ END DO NODE_LOOP3 ! !----------------------------------------------------------------------- ! Store old bed thickness. !----------------------------------------------------------------------- ! if (SED_MORPH)then do i=1,m sedbed%bed_thick(i,nnew)=0.0_sp DO k=1,Nbed sedbed%bed_thick(i,nnew)=sedbed%bed_thick(i,nnew)+ & & bed(i,k,ithck) END DO end do # if defined (MULTIPROCESSOR) if(par)then call aexchange(nc,myid,nprocs,sedbed%bed_thick) endif # endif bottom(i,morph) = sedbed%bed_thick(i,nnew) endif # if defined (MULTIPROCESSOR) if(par)then allocate(temp(0:mt,1:nbed)) DO ised=1,NST temp = 0.0_sp temp = sedbed%bed_frac(:,:,ised) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) sedbed%bed_frac(:,:,ised) = temp temp = 0.0_sp temp = sedbed%bed_mass(:,:,nnew,ised) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) sedbed%bed_mass(:,:,nnew,ised) = temp end do do ibed=1,MBEDP temp = 0.0_sp temp = bed(:,:,ibed) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) bed(:,:,ibed)=temp end do deallocate(temp) endif # endif if(dbg_set(dbg_sbr)) write(ipt,*) "End: calc_sed_bed_cohesive" return End Subroutine calc_sed_bed_cohesive2 !======================================================================= ! ! ! This routine computes sediment bed layer stratigraphy. ! ! ! ! Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G. ! ! Arango, 2008: Development of a three-dimensional, regional, ! ! coupled wave, current, and sediment-transport model, Computers ! ! & Geosciences, 34, 1284-1306. ! ! ! !======================================================================= Subroutine calc_sed_bed2 implicit none ! ! Local variable declarations. ! integer :: Ksed, i, ised,ibed, j, k, ks integer :: bnew, nnn, cnt,cc,jj real(sp), parameter :: eps = 1.0E-14_sp real(sp) :: cff, cff1, cff2, cff3 real(sp) :: thck_avail, thck_to_add real(sp), dimension(NST) :: nlysm real(sp), dimension(NST) :: dep_mass real(sp), dimension(0:mt) :: tau_w real(sp),allocatable,dimension(:,:)::temp if(dbg_set(dbg_sbr)) write(ipt,*) "Start: calc_sed_bed2" if(BEDLOAD)then bnew=nnew else bnew=nstp endif ! KLUDGE alert ! minlayer_thick = 0.0005 minlayer_thick = newlayer_thick ! !----------------------------------------------------------------------- ! Compute sediment bed layer stratigraphy. !----------------------------------------------------------------------- ! if(BEDLOAD_MPM.or.SUSLOAD)then tau_w = taub endif ! !----------------------------------------------------------------------- ! Update bed properties according to ero_flux and dep_flux. !----------------------------------------------------------------------- ! if (SUSLOAD)then NODE_LOOP : DO i=1,m ! ! The deposition and resuspension of sediment on the bottom "bed" ! is due to precipitation flux FC(:,0), already computed, and the ! resuspension (erosion, hence called ero_flux). The resuspension is ! applied to the bottom-most grid box value qc(:,1) so the total mass ! is conserved. Restrict "ero_flux" so that "bed" cannot go negative ! after both fluxes are applied. ! SED_LOOP: DO ised=1,NST dep_mass(ised)=0.0_sp if (SED_MORPH)then ! Apply morphology factor. sedbed%ero_flux(i,ised)=sedbed%ero_flux(i,ised)*sed(ised)%morph_fac sedbed%settling_flux(i,ised)=sedbed%settling_flux(i,ised)* & & sed(ised)%morph_fac endif ! Update bed mass arrays. sedbed%bed_mass(i,1,nnew,ised)=MAX(sedbed%bed_mass(i,1,bnew,ised)- & & (sedbed%ero_flux(i,ised)- & & sedbed%settling_flux(i,ised)), & & 0.0_sp) DO k=2,Nbed sedbed%bed_mass(i,k,nnew,ised)=sedbed%bed_mass(i,k,nstp,ised) END DO END DO SED_LOOP cff3=0.0_sp DO ised=1,NST cff3=cff3+sedbed%bed_mass(i,1,nnew,ised) END DO IF (cff3.eq.0.0_sp) THEN cff3=eps END IF bed(i,1,ithck)=0.0_sp DO ised=1,NST sedbed%bed_frac(i,1,ised)=sedbed%bed_mass(i,1,nnew,ised)/cff3 bed(i,1,ithck)=MAX(bed(i,1,ithck)+ & & sedbed%bed_mass(i,1,nnew,ised)/ & & (sed(ised)%Srho* & & (1.0_sp-bed(i,1,iporo))),0.0_sp) END DO END DO NODE_LOOP endif !/* SUSPLOAD section */ ! !----------------------------------------------------------------------- ! At this point, all deposition or erosion is complete, and ! has been added/subtracted to top layer. Thickness has NOT been corrected. !----------------------------------------------------------------------- ! NODE_LOOP2 : DO i=1,m ! Calculate active layer thickness, bottom(i,iactv). ! (trunk version allows this to be zero...this has minimum of 6*D50) bottom(i,iactv)=MAX(0.0_sp, & & 0.007_sp* & & (tau_w(i)-bottom(i,itauc)))+ & & 6.0_sp*bottom(i,isd50) if(SED_MORPH)then ! Apply morphology factor. bottom(i,iactv)=MAX(bottom(i,iactv)*sed(1)%morph_fac, & & bottom(i,iactv)) endif ! ! Calculate net deposition and erosion cff=0.0_sp cff2=0.0_sp DO ised=1,NST cff=cff+sedbed%settling_flux(i,ised) cff2=cff2+sedbed%ero_flux(i,ised) dep_mass(ised)=0.0_sp IF ((sedbed%ero_flux(i,ised)-sedbed%settling_flux(i,ised)).lt. & & 0.0_sp) THEN dep_mass(ised)=sedbed%settling_flux(i,ised)- & & sedbed%ero_flux(i,ised) END IF END DO IF ( cff-cff2.GT.0.0_sp) THEN ! NET depostion ! Deposition. Determine if we need to create a new bed layer ! (no test for age here) bed(i,1,iaged)=T_model IF(bed(i,1,ithck).gt. & & MAX(bottom(i,iactv),newlayer_thick)) THEN ! Top layer is too thick IF (Nbed.gt.2) THEN IF(bed(i,2,ithck).lt.minlayer_thick) THEN ! Layer 2 is smaller than minimum size ! Instead of pushing down all layers, just combine top 2 layers cff=0.0_sp cff1=0.0_sp cff2=0.0_sp DO ised=1,NST cff =cff +dep_mass(ised) cff1=cff1+sedbed%bed_mass(i,1,nnew,ised) cff2=cff2+sedbed%bed_mass(i,2,nnew,ised) END DO ! Update bed mass DO ised=1,NST sedbed%bed_mass(i,2,nnew,ised)= & & MAX(sedbed%bed_mass(i,2,nnew,ised)+ & & sedbed%bed_mass(i,1,nnew,ised)- & & dep_mass(ised),0.0_sp) sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised) END DO ! ALA - average T_model and porosity ! ALA CHECK WITH CRS cff1 or cff1-cff for first layer bed(i,2,iaged)=(bed(i,1,iaged)*cff1+ & & bed(i,2,iaged)*cff2)/(cff1+cff2) bed(i,1,iaged)=T_model bed(i,2,iporo)=(bed(i,1,iporo)*cff1+ & & bed(i,2,iporo)*cff2)/(cff1+cff2) ! ALA CHECK WITH CRS POROSITY OF 1ST LAYER bed(i,1,iporo)=bed(i,1,iporo) ELSE ! Layer 2 is > minlayer thick, need another layer ! Combine bottom layers. cff1=0.0_sp cff2=0.0_sp DO ised=1,NST cff1=cff1+sedbed%bed_mass(i,Nbed-1,nnew,ised) cff2=cff2+sedbed%bed_mass(i,Nbed,nnew,ised) END DO bed(i,Nbed,iporo)= & & (bed(i,Nbed-1,iporo)*cff1+ & & bed(i,Nbed,iporo)*cff2)/(cff1+cff2) bed(i,Nbed,iaged)= & & (bed(i,Nbed-1,iaged)*cff1+ & & bed(i,Nbed,iaged)*cff2)/(cff1+cff2) DO ised=1,NST sedbed%bed_mass(i,Nbed,nnew,ised)= & & sedbed%bed_mass(i,Nbed-1,nnew,ised)+ & & sedbed%bed_mass(i,Nbed ,nnew,ised) END DO ! ! Push layers down. DO k=Nbed-1,2,-1 bed(i,k,iporo)=bed(i,k-1,iporo) bed(i,k,iaged)=bed(i,k-1,iaged) DO ised =1,NST sedbed%bed_mass(i,k,nnew,ised)= & & sedbed%bed_mass(i,k-1,nnew,ised) END DO END DO ! Set new top parameters for top 2 layers DO ised=1,NST sedbed%bed_mass(i,2,nnew,ised)= & & MAX(sedbed%bed_mass(i,2,nnew,ised)- & & dep_mass(ised),0.0_sp) sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised) END DO END IF ELSEIF (Nbed.eq.2) THEN ! NBED=2 cff1=0.0_sp cff2=0.0_sp DO ised=1,NST cff1=cff1+sedbed%bed_mass(i,1,nnew,ised) cff2=cff2+sedbed%bed_mass(i,2,nnew,ised) END DO DO ised=1,NST sedbed%bed_mass(i,2,nnew,ised)= & & MAX(sedbed%bed_mass(i,2,nnew,ised)+ & & sedbed%bed_mass(i,1,nnew,ised)- & & dep_mass(ised),0.0_sp) sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised) END DO ! ALA - average T_model and porosity bed(i,2,iaged)=(bed(i,1,iaged)*cff1+ & & bed(i,2,iaged)*cff2)/(cff1+cff2) bed(i,1,iaged)=T_model bed(i,2,iporo)=(bed(i,1,iporo)*cff1+ & & bed(i,2,iporo)*cff2)/(cff1+cff2) ! ALA CHECK WITH CRS POROSITY OF 1ST LAYER bed(i,1,iporo)=bed(i,1,iporo) ELSE ! NBED=1 END IF ELSE ! Net deposition has occured, but no new bed layer was created END IF ELSE ! Net erosion occurred bed(i,1,iaged)=T_model IF (Nbed.eq.2) THEN ! NBED=2 DO ised=1,NST sedbed%bed_mass(i,2,nnew,ised)= & & MAX(sedbed%bed_mass(i,2,nnew,ised)+ & & sedbed%bed_mass(i,1,nnew,ised)- & & dep_mass(ised),0.0_sp) sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised) END DO ELSEIF (Nbed.eq.1) THEN ! ALF NO NEED TO DO ANYTHING ELSE END IF END IF ! Recalculate thickness and fractions for all layers. DO k=1,Nbed cff3=0.0_sp DO ised=1,NST cff3=cff3+sedbed%bed_mass(i,k,nnew,ised) END DO IF (cff3.eq.0.0_sp) THEN cff3=eps END IF bed(i,k,ithck)=0.0_sp DO ised=1,NST sedbed%bed_frac(i,k,ised)=sedbed%bed_mass(i,k,nnew,ised)/cff3 bed(i,k,ithck)=MAX(bed(i,k,ithck)+ & & sedbed%bed_mass(i,k,nnew,ised)/ & & (sed(ised)%Srho* & & (1.0_sp-bed(i,k,iporo))),0.0_sp) END DO END DO END DO NODE_LOOP2 NODE_LOOP3 : DO i=1,m IF (bottom(i,iactv).gt.bed(i,1,ithck)) THEN IF (Nbed.eq.1) THEN bottom(i,iactv)=bed(i,1,ithck) ELSE thck_to_add=bottom(i,iactv)-bed(i,1,ithck) thck_avail=0.0_sp Ksed=1 ! initialize DO k=2,Nbed IF (thck_avail.lt.thck_to_add) THEN thck_avail=thck_avail+bed(i,k,ithck) Ksed=k END IF END DO ! ! Catch here if there was not enough bed material. ! IF (thck_avail.lt.thck_to_add) THEN bottom(i,iactv)=bed(i,1,ithck)+thck_avail thck_to_add=thck_avail END IF ! ! Update bed mass of top layer and fractional layer. ! cff2=MAX(thck_avail-thck_to_add,0.0_sp)/ & & MAX(bed(i,Ksed,ithck),eps) DO ised=1,NST cff1=0.0_sp DO k=1,Ksed cff1=cff1+sedbed%bed_mass(i,k,nnew,ised) END DO cff3=cff2*sedbed%bed_mass(i,Ksed,nnew,ised) sedbed%bed_mass(i,1 ,nnew,ised)=cff1-cff3 sedbed%bed_mass(i,Ksed,nnew,ised)=cff3 END DO ! ! Update thickness of fractional layer ksource_sed. ! bed(i,Ksed,ithck)=MAX(thck_avail-thck_to_add,0.0_sp) ! ! Update bed fraction of top layer. ! cff3=0.0_sp DO ised=1,NST cff3=cff3+sedbed%bed_mass(i,1,nnew,ised) END DO IF (cff3.eq.0.0_sp) THEN cff3=eps END IF DO ised=1,NST sedbed%bed_frac(i,1,ised)=sedbed%bed_mass(i,1,nnew,ised)/cff3 END DO ! ! Upate bed thickness of top layer. ! bed(i,1,ithck)=bottom(i,iactv) ! ! Pull all layers closer to the surface. ! DO k=Ksed,Nbed ks=Ksed-2 bed(i,k-ks,ithck)=bed(i,k,ithck) bed(i,k-ks,iporo)=bed(i,k,iporo) bed(i,k-ks,iaged)=bed(i,k,iaged) DO ised=1,NST sedbed%bed_frac(i,k-ks,ised)=sedbed%bed_frac(i,k,ised) sedbed%bed_mass(i,k-ks,nnew,ised)=sedbed%bed_mass(i,k,nnew,ised) END DO END DO ! ! Add new layers onto the bottom. Split what was in the bottom layer to ! fill these new empty cells. ("ks" is the number of new layers). ! ks=Ksed-2 ! ALA CHECK WITH CRS about bed_frac nnn=0 DO ised=1,NST nlysm(ised)=newlayer_thick*REAL(ks+1,sp)* & & (sed(ised)%Srho* & & (1.0_sp-bed(i,Nbed-ks,iporo)))* & & sedbed%bed_frac(i,Nbed-ks,ised) IF (ks.gt.0) THEN IF (sedbed%bed_mass(i,Nbed-ks,nnew,ised).gt. & & nlysm(ised)) THEN nnn=nnn+1 nlysm(ised)= & & newlayer_thick*REAL(ks,sp)* & & (sed(ised)%Srho* & & (1.0_sp-bed(i,Nbed-ks,iporo)))* & & sedbed%bed_frac(i,Nbed-ks,ised) END IF END IF END DO IF (nnn.eq.NST) THEN bed(i,Nbed,ithck)=bed(i,Nbed-ks,ithck)- & & newlayer_thick*REAL(ks,sp) DO ised=1,NST sedbed%bed_mass(i,Nbed,nnew,ised)= & & sedbed%bed_mass(i,Nbed-ks,nnew,ised)-nlysm(ised) END DO DO k=Nbed-1,Nbed-ks,-1 bed(i,k,ithck)=newlayer_thick bed(i,k,iaged)=bed(i,Nbed-ks,iaged) DO ised=1,NST sedbed%bed_frac(i,k,ised)=sedbed%bed_frac(i,Nbed-ks,ised) sedbed%bed_mass(i,k,nnew,ised)= & & nlysm(ised)/REAL(ks,sp) END DO END DO ELSE cff=1.0_sp/REAL(ks+1,sp) DO k=Nbed,Nbed-ks,-1 bed(i,k,ithck)=bed(i,Nbed-ks,ithck)*cff bed(i,k,iaged)=bed(i,Nbed-ks,iaged) DO ised=1,NST sedbed%bed_frac(i,k,ised)=sedbed%bed_frac(i,Nbed-ks,ised) sedbed%bed_mass(i,k,nnew,ised)= & & sedbed%bed_mass(i,Nbed-ks,nnew,ised)*cff END DO END DO END IF END IF ! Nbed > 1 END IF ! increase top bed layer END DO NODE_LOOP3 ! !----------------------------------------------------------------------- ! Store old bed thickness. !----------------------------------------------------------------------- ! if (SED_MORPH)then do i=1,m sedbed%bed_thick(i,nnew)=0.0_sp DO k=1,Nbed sedbed%bed_thick(i,nnew)=sedbed%bed_thick(i,nnew)+ & & bed(i,k,ithck) END DO end do # if defined (MULTIPROCESSOR) if(par)then call aexchange(nc,myid,nprocs,sedbed%bed_thick) endif # endif bottom(i,morph) = sedbed%bed_thick(i,nnew) endif ! !----------------------------------------------------------------------- ! Apply periodic or gradient boundary conditions to property arrays. !----------------------------------------------------------------------- ! # if defined (MULTIPROCESSOR) if(par)then allocate(temp(0:mt,1:nbed)) DO ised=1,NST temp = 0.0_sp temp = sedbed%bed_frac(:,:,ised) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) sedbed%bed_frac(:,:,ised) = temp temp = 0.0_sp temp = sedbed%bed_mass(:,:,nnew,ised) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) sedbed%bed_mass(:,:,nnew,ised) = temp end do do ibed=1,MBEDP temp = 0.0_sp temp = bed(:,:,ibed) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) bed(:,:,ibed)=temp end do deallocate(temp) endif # endif if(dbg_set(dbg_sbr)) write(ipt,*) "End: calc_sed_bed2" return End Subroutine calc_sed_bed2 !======================================================================= ! ! ! This routine computes sediment bed layer stratigraphy. ! ! ! ! Warner, J.C., C.R. Sherwood, R.P. Signell, C.K. Harris, and H.G. ! ! Arango, 2008: Development of a three-dimensional, regional, ! ! coupled wave, current, and sediment-transport model, Computers ! ! & Geosciences, 34, 1284-1306. ! ! ! !======================================================================= Subroutine calc_sed_bed implicit none ! ! Local variable declarations. ! integer :: Ksed, i, ised,ibed, j, k, ks integer :: bnew, nnn, cnt,cc,jj real(sp), parameter :: eps = 1.0E-14_sp real(sp) :: cff, cff1, cff2, cff3 real(sp) :: thck_avail, thck_to_add real(sp), dimension(NST) :: dep_mass real(sp), dimension(0:mt) :: tau_w real(sp),allocatable,dimension(:,:)::temp if(dbg_set(dbg_sbr)) write(ipt,*) "Start: calc_sed_bed" if(BEDLOAD)then bnew=nnew else bnew=nstp endif ! !----------------------------------------------------------------------- ! Compute sediment bed layer stratigraphy. !----------------------------------------------------------------------- ! if(BEDLOAD_MPM.or.SUSLOAD)then tau_w = taub endif ! !----------------------------------------------------------------------- ! Update bed properties according to ero_flux and dep_flux. !----------------------------------------------------------------------- ! if(SUSLOAD)then NODE_LOOP : DO i=1,m SED_LOOP: DO ised=1,NST ! ! The deposition and resuspension of sediment on the bottom "bed" ! is due to precepitation flux FC(:,0), already computed, and the ! resuspension (erosion, hence called ero_flux). The resuspension is ! applied to the bottom-most grid box value qc(:,1) so the total mass ! is conserved. Restrict "ero_flux" so that "bed" cannot go negative ! after both fluxes are applied. ! dep_mass(ised)=0.0_sp if( SED_MORPH)then ! ! Apply morphology factor. ! sedbed%ero_flux(i,ised)=sedbed%ero_flux(i,ised)*sed(ised)%morph_fac sedbed%settling_flux(i,ised)=sedbed%settling_flux(i,ised)* & & sed(ised)%morph_fac ! endif IF ((sedbed%ero_flux(i,ised)-sedbed%settling_flux(i,ised)).lt. & & 0.0_sp) THEN ! ! If first time step of deposit, then store deposit material in ! temporary array, dep_mass. ! IF ((T_model.gt.(bed(i,1,iaged)+1.1_sp*DTsed)).and. & & (bed(i,1,ithck).gt.newlayer_thick)) THEN dep_mass(ised)=sedbed%settling_flux(i,ised)- & & sedbed%ero_flux(i,ised) END IF bed(i,1,iaged)=T_model END IF ! ! Update bed mass arrays. ! sedbed%bed_mass(i,1,nnew,ised)=MAX(sedbed%bed_mass(i,1,bnew,ised)- & & (sedbed%ero_flux(i,ised)- & & sedbed%settling_flux(i,ised)), & & 0.0_sp) DO k=2,Nbed sedbed%bed_mass(i,k,nnew,ised)=sedbed%bed_mass(i,k,nstp,ised) END DO END DO SED_LOOP ! ! If first time step of deposit, create new layer and combine bottom ! two bed layers. ! cff=0.0_sp ! ! Determine if deposition ocurred here. ! IF (Nbed.gt.1) THEN DO ised=1,NST cff=cff+dep_mass(ised) END DO IF (cff.gt.0.0_sp) THEN ! ! Combine bottom layers. ! bed(i,Nbed,iporo)=0.5_sp*(bed(i,Nbed-1,iporo)+ & & bed(i,Nbed,iporo)) bed(i,Nbed,iaged)=0.5_sp*(bed(i,Nbed-1,iaged)+ & & bed(i,Nbed,iaged)) DO ised=1,NST sedbed%bed_mass(i,Nbed,nnew,ised)= & & sedbed%bed_mass(i,Nbed-1,nnew,ised)+ & & sedbed%bed_mass(i,Nbed ,nnew,ised) END DO ! ! Push layers down. ! DO k=Nbed-1,2,-1 bed(i,k,iporo)=bed(i,k-1,iporo) bed(i,k,iaged)=bed(i,k-1,iaged) DO ised =1,NST sedbed%bed_mass(i,k,nnew,ised)=sedbed%bed_mass(i,k-1,nnew,ised) END DO END DO ! ! Set new top layer parameters. ! DO ised=1,NST sedbed%bed_mass(i,2,nnew,ised)=MAX(sedbed%bed_mass(i,2,nnew,ised)-& & dep_mass(ised),0.0_sp) sedbed%bed_mass(i,1,nnew,ised)=dep_mass(ised) END DO END IF END IF !NBED=1 ! ! Recalculate thickness and fractions for all layers. ! DO k=1,Nbed cff3=0.0_sp DO ised=1,NST cff3=cff3+sedbed%bed_mass(i,k,nnew,ised) END DO IF (cff3.eq.0.0_sp) THEN cff3=eps END IF bed(i,k,ithck)=0.0_sp DO ised=1,NST sedbed%bed_frac(i,k,ised)=sedbed%bed_mass(i,k,nnew,ised)/cff3 bed(i,k,ithck)=MAX(bed(i,k,ithck)+ & & sedbed%bed_mass(i,k,nnew,ised)/ & & (sed(ised)%Srho* & & (1.0_sp-bed(i,k,iporo))),0.0_sp) END DO END DO END DO NODE_LOOP ! ! End of Suspended Sediment only section. ! endif ! ! Ensure top bed layer thickness is greater or equal than active layer ! thickness. If need to add sed to top layer, then entrain from lower ! levels. Create new layers at bottom to maintain Nbed. ! NODE_LOOP2 : DO i=1,m ! ! Calculate active layer thickness, bottom(i,iactv). ! bottom(i,iactv)=MAX(0.0_sp, & & 0.007_sp* & & (tau_w(i)-bottom(i,itauc)))+ & & 6.0_sp*bottom(i,isd50) ! if (SED_MORPH)then ! ! Apply morphology factor. ! bottom(i,iactv)=MAX(bottom(i,iactv)*sed(1)%morph_fac, & & bottom(i,iactv)) endif ! IF (bottom(i,iactv).gt.bed(i,1,ithck)) THEN IF (Nbed.eq.1) THEN bottom(i,iactv)=bed(i,1,ithck) ELSE thck_to_add=bottom(i,iactv)-bed(i,1,ithck) thck_avail=0.0_sp Ksed=1 ! initialize DO k=2,Nbed IF (thck_avail.lt.thck_to_add) THEN thck_avail=thck_avail+bed(i,k,ithck) Ksed=k END IF END DO ! ! Catch here if there was not enough bed material. ! IF (thck_avail.lt.thck_to_add) THEN bottom(i,iactv)=bed(i,1,ithck)+thck_avail thck_to_add=thck_avail END IF ! ! Update bed mass of top layer and fractional layer. ! cff2=MAX(thck_avail-thck_to_add,0.0_sp)/ & & MAX(bed(i,Ksed,ithck),eps) DO ised=1,NST cff1=0.0_sp DO k=1,Ksed cff1=cff1+sedbed%bed_mass(i,k,nnew,ised) END DO cff3=cff2*sedbed%bed_mass(i,Ksed,nnew,ised) sedbed%bed_mass(i,1 ,nnew,ised)=cff1-cff3 sedbed%bed_mass(i,Ksed,nnew,ised)=cff3 END DO ! ! Update thickness of fractional layer ksource_sed. ! bed(i,Ksed,ithck)=MAX(thck_avail-thck_to_add,0.0_sp) ! ! Upate bed fraction of top layer. ! cff3=0.0_sp DO ised=1,NST cff3=cff3+sedbed%bed_mass(i,1,nnew,ised) END DO IF (cff3.eq.0.0_sp) THEN cff3=eps END IF DO ised=1,NST sedbed%bed_frac(i,1,ised)=sedbed%bed_mass(i,1,nnew,ised)/cff3 END DO ! ! Upate bed thickness of top layer. ! bed(i,1,ithck)=bottom(i,iactv) ! ! Pull all layers closer to the surface. ! DO k=Ksed,Nbed ks=Ksed-2 bed(i,k-ks,ithck)=bed(i,k,ithck) bed(i,k-ks,iporo)=bed(i,k,iporo) bed(i,k-ks,iaged)=bed(i,k,iaged) DO ised=1,NST sedbed%bed_frac(i,k-ks,ised)=sedbed%bed_frac(i,k,ised) sedbed%bed_mass(i,k-ks,nnew,ised)=sedbed%bed_mass(i,k,nnew,ised) END DO END DO ! ! Add new layers onto the bottom. Split what was in the bottom layer to ! fill these new empty cells. ("ks" is the number of new layers). ! ks=Ksed-2 cff=1.0_sp/REAL(ks+1,sp) DO k=Nbed,Nbed-ks,-1 bed(i,k,ithck)=bed(i,Nbed-ks,ithck)*cff bed(i,k,iaged)=bed(i,Nbed-ks,iaged) DO ised=1,NST sedbed%bed_frac(i,k,ised)=sedbed%bed_frac(i,Nbed-ks,ised) sedbed%bed_mass(i,k,nnew,ised)= & & sedbed%bed_mass(i,Nbed-ks,nnew,ised)*cff END DO END DO END IF ! Nbed > 1 END IF ! increase top bed layer END DO NODE_LOOP2 ! !----------------------------------------------------------------------- ! Store old bed thickness. !----------------------------------------------------------------------- ! if (SED_MORPH)then do i=1,m sedbed%bed_thick(i,nnew)=0.0_sp DO k=1,Nbed sedbed%bed_thick(i,nnew)=sedbed%bed_thick(i,nnew)+ & & bed(i,k,ithck) END DO end do # if defined (MULTIPROCESSOR) if(par)then call aexchange(nc,myid,nprocs,sedbed%bed_thick) endif # endif bottom(i,morph) = sedbed%bed_thick(i,nnew) endif # if defined (MULTIPROCESSOR) if(par)then allocate(temp(0:mt,1:nbed)) DO ised=1,NST temp = 0.0_sp temp = sedbed%bed_frac(:,:,ised) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) sedbed%bed_frac(:,:,ised) = temp temp = 0.0_sp temp = sedbed%bed_mass(:,:,nnew,ised) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) sedbed%bed_mass(:,:,nnew,ised) = temp end do do ibed=1,MBEDP temp = 0.0_sp temp = bed(:,:,ibed) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) bed(:,:,ibed)=temp end do deallocate(temp) endif # endif if(dbg_set(dbg_sbr)) write(ipt,*) "End: calc_sed_bed" return End Subroutine calc_sed_bed !======================================================================= ! ! ! This routine computes sediment bed layer mixing (biodiffusion) ! ! ! ! ! !======================================================================= Subroutine calc_sed_biodiff implicit none ! ! Local variable declarations. ! integer :: i, j, k, ised,ibed real(sp), parameter :: eps = 1.0E-14_sp real(sp), parameter :: cmy2ms = 3.1688765E-12_sp ! multiply cm2/yr by this to get m2/s real(sp) :: cff, cff1, cff2, cff3 real(sp), dimension(0:kb) :: dep_mass integer :: iu,il,lp,ii ! real(sp) :: rtemp, Dbmx, Dbmm, zs, zm, zp real(sp) :: rtemp real(sp), dimension(Nbed) :: zb real(sp), dimension(Nbed) :: zc real(sp), dimension(Nbed) :: dzui real(sp), dimension(Nbed) :: dzli real(sp), dimension(Nbed) :: dzmi real(sp), dimension(Nbed) :: Db real(sp), dimension(Nbed) :: Dc real(sp), dimension(Nbed) :: a real(sp), dimension(Nbed) :: d real(sp), dimension(Nbed) :: b real(sp), dimension(Nbed) :: cc real(sp), dimension(Nbed) :: dd real(sp),allocatable,dimension(:,:)::temp if(dbg_set(dbg_sbr)) write(ipt,*) "Start: calc_sed_biodiff " ! !----------------------------------------------------------------------- ! Compute sediment bed layer stratigraphy. !----------------------------------------------------------------------- ! ! print *, 'Mixing...' NODE_LOOP : DO i=1,m ! ! Set mixing coefficient profile ! (hardwire uniform mixing) ! ! DO k=1,Nbed ! Db(k)=1.0E-8_sp ! ENDDO IF (Nbed.GT.2) THEN ! Compute cumulative depth zb(1)=bed(i,1,ithck) DO k=2,Nbed zb(k)=zb(k-1)+bed(i,k,ithck) END DO ! Compute depths of bed centers zc(1)=0.5_sp*(bed(i,1,ithck)) DO k=2,Nbed zc(k)=zb(k-1)+0.5_sp*bed(i,k,ithck) END DO ! ! Set biodiffusivity profile if (DB_PROFILE)then ! Depth-varying biodiffusivity profile ! !ALF Dbmx = bottom(i,idbmx) !ALF Dbmm = bottom(i,idbmm) !ALF zs = bottom(i,idbzs) !ALF zm = bottom(i,idbzm) !ALF zp = bottom(i,idbzp) DO k=1,Nbed Db(k)= Dbmx ! IF( zb(k).GT.Dbzp)THEN ! should be .GE. ? IF( zb(k).GE.Dbzp)THEN Db(k)=0.0_sp ELSE IF((zb(k) .GT. Dbzs).AND. & & (zb(k) .LE. Dbzm)) THEN rtemp= LOG(Dbmm/Dbmx)/(-Dbzm-Dbzs) Db(k)=Dbmx*exp(rtemp*(-zb(k)-Dbzs)) ELSEIF((zb(k).GT.Dbzm).AND. & & (zb(k).LT.Dbzp) ) THEN Db(k)=(Dbmm-(Dbmm/(Dbzp-Dbzm)))* & & (zb(k)-Dbzm) ENDIF ENDIF END DO else ! Uniform biodiffusivity profile at max value ! DO k=1,Nbed Db(k)=Dbmx ! Db(k)=bottom(i,idbmx) ENDDO endif !/* defined DB_PROFILE */ ! write(*,*) 'Db',Db ! Calculate finite differences dzui(1)=1.0_sp/(zc(2)-zc(1)) dzli(1)=1.E35_sp ! should not be needed dzmi(1)=1.0_sp/bed(i,1,ithck) DO k=2,Nbed-1 dzui(k)=1.0_sp/(zc(k+1)-zc(k)) dzli(k)=1.0_sp/(zc(k)-zc(k-1)) !dzmi(k)=1.0_sp/(zb(k+1)-zb(k)) ! equivalent: dzmi(k)=1.0/bed(i,k,ithck) ! Tridiagonal terms b(k)= -dtsed*dzmi(k)*Db(k-1)*dzli(k) d(k)=1.0_sp+dtsed*dzmi(k)* & & ( Db(k-1)*dzli(k)+Db(k)*dzui(k) ) a(k)= -dtsed*dzmi(k)*Db(k)*dzui(k) ENDDO dzui(Nbed)=1.0E35_sp ! should not be needed dzli(Nbed)=1.0_sp/(zc(Nbed)-zc(Nbed-1)) dzmi(Nbed)=1.0_sp/bed(i,Nbed,ithck) ! No-flux boundary conditions b(1) = 999.9_sp ! should not be needed d(1)= 1.0_sp +dtsed*dzmi(1)*Db(1)*dzui(1) a(1)= -dtsed*dzmi(1)*Db(1)*dzui(1) b(Nbed)= -dtsed*dzmi(Nbed)*Db(Nbed-1)*dzli(Nbed) d(Nbed)=1.0_sp+dtsed*dzmi(Nbed)*Db(Nbed-1)*dzli(Nbed) a(Nbed)=999.9_sp ! should not be needed ! ! Calculate mixing for each size fraction DO ised=1,NST ! ...make working copies DO k=1,Nbed cc(k) = sedbed%bed_frac(i,k,ised) dd(k)= d(k) ENDDO ! Solve a tridiagonal system of equations using Thomas' algorithm ! Anderson, Tannehill, and Pletcher (1984) pp. 549-550 ! ...establish upper triangular matrix il = 1 iu = Nbed lp = il+1 DO k = lp,iu rtemp = b(k)/dd(k-1) dd(k)= dd(k)-rtemp*a(k-1); cc(k)= cc(k)-rtemp*cc(k-1); ENDDO ! ...back substitution cc(iu) = cc(iu)/dd(iu) DO k = lp,iu ii = iu-k+il; cc(ii) = (cc(ii)-a(ii)*cc(ii+1))/dd(ii); ENDDO ! ...solution stored in cc; copy out DO k = 1,Nbed sedbed%bed_frac(i,k,ised)=cc(k) ENDDO ENDDO ! TODO - Mix porosity or assign it as f(depth)? ! TODO - Mix age? ! Recompute bed masses DO k=1,Nbed ! write(*,*) i,k,(bed_frac(i,k,ised),ised=1,NST) ! debugging: ensure fracs add up to 1 cff3 = 0.0_sp DO ised=1,NST cff3 = cff3+sedbed%bed_frac(i,k,ised) ENDDO if( abs(1.0_sp-cff3).GT.1e-6 ) & & write(*,*) 'error: sum_frac: ',cff3 cff3=0.0_sp DO ised=1,NST cff3=cff3+sedbed%bed_mass(i,k,nnew,ised) ENDDO DO ised=1,NST sedbed%bed_mass(i,k,nnew,ised)=sedbed%bed_frac(i,k,ised)*cff3 ENDDO ENDDO END IF !NBED.GT.2 END DO NODE_LOOP # if defined (MULTIPROCESSOR) if(par)then allocate(temp(0:mt,1:nbed)) DO ised=1,NST temp = 0.0_sp temp = sedbed%bed_frac(:,:,ised) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) sedbed%bed_frac(:,:,ised) = temp temp = 0.0_sp temp = sedbed%bed_mass(:,:,nnew,ised) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) sedbed%bed_mass(:,:,nnew,ised) = temp end do do ibed=1,MBEDP temp = 0.0_sp temp = bed(:,:,ibed) call node_match(1,nbn,bn_mlt,bn_loc,bnc,mt,nbed,myid,nprocs,temp) call aexchange(nc,myid,nprocs,temp) bed(:,:,ibed)=temp end do deallocate(temp) endif # endif if(dbg_set(dbg_sbr)) write(ipt,*) "End: calc_sed_biodiff" return End Subroutine calc_sed_biodiff ! !----------------------------------------------------------------------- ! Horizontal Advection of Sediment Concentration !----------------------------------------------------------------------- ! Subroutine calc_fluid_mud implicit none integer :: i if(dbg_set(dbg_sbr)) write(ipt,*) "Start: calc_fluid_mud " # if defined (FLUID_MUD) call fluid_mud_source_sink(DTsed,sed(1)%Srho,bottom(:,bflag)& &,bed(:,1,ithck),bed(:,1,iporo),sedbed%frac(:,1,1),taub,sed(1)%t_ce,sed(1)%rate) call internal_step_fluid_mud do ised=1,nst do i=1,m # if defined (WET_DRY) if(iswetn(i)==1)then # endif sed(ised)%conc(i,kbm1) = sed(ised)%conc(i,kbm1) + Entrainment(i)/(d(i)*dz(i,kbm1)) sed(ised)%dflx(i) = dewatering(i) sed(ised)%eflx(i) = BedErosion(i) # if defined (WET_DRY) else sed(ised)%conc(i,kbm1) = 0.0_sp sed(ised)%dflx(i) = 0.0_sp sed(ised)%eflx(i) = 0.0_sp end if # endif end do end do # endif if(dbg_set(dbg_sbr)) write(ipt,*) "End: calc_fluid_mud " Return End Subroutine calc_fluid_mud ! !----------------------------------------------------------------------- ! Horizontal Advection of Sediment Concentration !----------------------------------------------------------------------- ! Subroutine Calc_Sed_Advection Use Scalar implicit none integer :: i,k,ised,d_cdis,d_cflx if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Calc_Sed_Advection " if(.not. oned_model)then d_cdis = 1 if(numqbc > 0) d_cdis = numqbc d_cflx = 1 if(Nobc > 0) d_cflx = Nobc do i=1,nst !set discharge sediment concentrations if(numqbc > 0 .and. sed_source)then sed(i)%cdis(1:numqbc) = seddis(1:numqbc,i) endif !3D advection ! call adv_scal(sed(i)%conc, & ! sed(i)%cnew, & ! d_cdis, & ! sed(i)%cdis(1:numqbc), & ! d_cflx, & ! sed(i)%cflx, & ! DTsed ,sed_source) !J. Ge for tracer advection call adv_scal(sed(i)%conc,sed0(i)%conc,sed2(i)%conc, & sed(i)%cnew, & d_cdis, & sed(i)%cdis(1:numqbc), & d_cflx, & sed(i)%cflx, & DTsed ,sed_source) !J. Ge for tracer advection call fct_sed(sed(i)%conc,sed(i)%cnew) end do else do i=1,nst sed(i)%cnew = sed(i)%conc end do endif if(dbg_set(dbg_sbr)) write(ipt,*) "End: Calc_Sed_Advection " return End Subroutine Calc_Sed_Advection ! !------------------------------------------------------ !Vertical Diffusion of Sediment Concentrations !------------------------------------------------------ ! Subroutine Calc_Sed_Diffusion Use Scalar implicit none integer :: i if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Calc_Sed_Diffusion " !Jianzhong GE 03/05/2013 if(vert_hindered)call hindered_diffusion !Jianzhong GE 03/05/2013 do i=1,nst call vdif_scal(sed(i)%cnew,DTsed) end do if(dbg_set(dbg_sbr)) write(ipt,*) "End: Calc_Sed_Diffusion " return End Subroutine Calc_Sed_Diffusion !Jianzhong Ge 03/05/2013 Subroutine Hindered_Diffusion use lims, only: m use all_vars,only: kh,kbm1 implicit none real(sp) :: sed_conc(1:m,1:kbm1) !mass concentration real(sp) :: vol_conc(1:m,1:kbm1) !volume concentration integer :: i,ii,k real,parameter :: cs0=0.65 real :: phi if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Hindered_Diffusion" sed_conc = 0.0 do ii=1,nst sed_conc=sed_conc+sed(ii)%cnew end do do i=1,m do k=1,kbm1 if(sed_conc(i,k)>2.5)then ! vol_conc(i,k)=sed_conc(i,k)/2650.0 ! phi=1+(vol_conc(i,k)/Cs0)**0.8-2.0*(vol_conc(i,k)/Cs0)**0.4 ! kh(i,k)=kh(i,k)*phi**4.0 kh(i,k)=kh(i,k)*0.02 ! print*,kh(i,k),phi**4.0,vol_conc(i,k),sed_conc(i,k) end if end do end do if(dbg_set(dbg_sbr)) write(ipt,*) "End: Hindered_Diffusion" return End Subroutine Hindered_Diffusion !Jianzhong Ge 03/05/2013 !------------------------------------------------------ ! Set Boundary Conditions !------------------------------------------------------ Subroutine Calc_Sed_Bcond Use Scalar implicit none integer :: i real(sp) :: temp if(dbg_set(dbg_sbr)) write(ipt,*) "Start: Calc_Sed_Bcond" !------------------------------------------------------ !Exchange Concentration on Processor Boundaries !------------------------------------------------------ # if defined (MULTIPROCESSOR) if(par)then do i=1,nst call aexchange(nc,myid,nprocs,sed(i)%cnew) end do endif # endif !------------------------------------------------------ ! Set Open Boundary Conditions (nudging to background) !------------------------------------------------------ if(.not. oned_model)then if(sed_ramp == 0) sed_ramp = 1 temp = min(float(sed_its/sed_ramp), 1.0)*sed_alpha if(.not.sed_nudge) temp = 0.0 do i=1,nst call bcond_scal_OBC(sed(i)%conc, & sed(i)%cnew, & sed(i)%cflx, & sed(i)%cobc, & DTsed, & temp) end do endif !------------------------------------------------------ ! Point Source Boundary Conditions !------------------------------------------------------ if(sed_source .and. .not. oned_model)then do i=1,nst call bcond_scal_PTsource(sed(i)%conc, & sed(i)%cnew, & sed(i)%cdis) end do endif if(dbg_set(dbg_sbr)) write(ipt,*) "End: Calc_Sed_Bcond" return End Subroutine Calc_Sed_Bcond !========================================================================== ! Update Total Sediment Thickness for Reporting !========================================================================== Subroutine Update_Thickness_Delta use lims, only: m implicit none integer :: i,k,ised !shift last thickness to lthck bottom(:,lthck) = bottom(:,nthck) bottom(:,nthck) = 0.0 !calculate new thickness do i=1,m do k=1,nbed bottom(i,nthck) = bottom(i,nthck) + bed(i,k,ithck) end do end do !calculate new delta do i=1,m bottom(i,dthck) = bottom(i,dthck) + bottom(i,nthck) - bottom(i,lthck) end do !calculate new total mass do i=1,m bottom(i,tmass) = 0.0 do ised=1,nsed do k=1,nbed bottom(i,tmass) = bottom(i,tmass) + sedbed%bed_mass(i,k,nnew,ised) end do end do end do End Subroutine Update_Thickness_Delta Subroutine finalize_sed_step implicit none integer :: i,ised if(dbg_set(dbg_sbr)) write(ipt,*) "Start: finalize_sed_step" !------------------------------------------------------ ! Update Concentration Variables to next time level !------------------------------------------------------ do ised=1,nst sed(ised)%conc = sed(ised)%cnew end do !------------------------------------------------------ ! store the fraction at bed surface !------------------------------------------------------ do ised=1,nst sed(ised)%frac(1:m,1:nbed)=sedbed%bed_frac(1:m,1:nbed,ised) end do !------------------------------------------------------ ! Limit to Non-Negative (modify to Venkat Limiter) !------------------------------------------------------ do i=1,nst sed(i)%conc = max(sed(i)%cnew,0.0) end do !------------------------------------------------------ ! Calculate total sediment concentration (all classes) !------------------------------------------------------ csed = 0.0 do i=1,nst csed = csed + sed(i)%conc end do !------------------------------------------------------ !Exchange Concentration on Processor Boundaries !------------------------------------------------------ # if defined (MULTIPROCESSOR) if(par)then do i=1,nst call aexchange(nc,myid,nprocs,sed(i)%cnew) end do call aexchange(nc,myid,nprocs,taub) endif # endif !J. Ge for tracer advection IF(BACKWARD_ADVECTION==.TRUE.)THEN IF(BACKWARD_STEP==1)THEN do i=1,nst SED0(i)%conc = sed(i)%cnew end do ELSEIF(BACKWARD_STEP==2)THEN do i=1,nst SED2(i)%conc = SED0(i)%conc SED0(i)%conc = sed(i)%cnew end do ENDIF ENDIF !J. Ge for tracer advection if(dbg_set(dbg_sbr)) write(ipt,*) "End: finalize_sed_step" Return End Subroutine finalize_sed_step SUBROUTINE UPDATE_DEPTH !------------------------------------------------------------------------------| USE ALL_VARS USE MOD_OBCS IMPLICIT NONE REAL(SP) :: TEMP, DX, DY INTEGER :: I,I1,K,J1,J2,CC LOGICAL :: ISONB2 !------------------------------------------------------------------------------| IF(DBG_SET(DBG_SBR)) write(ipt,*) "START: UPDATE_DEPTH" IF(DBG_SET(DBG_SBR)) write(ipt,*) "END: UPDATE_DEPTH" END SUBROUTINE UPDATE_DEPTH # endif End Module Mod_Sed_CSTMS