subroutine eqtran(sig ,thick ,kmax ,ws ,ltur , & & frac ,sigmol ,dicww ,lundia ,taucr0 , & & rksrs ,i2d3d ,lsecfl ,spirint ,suspfrac , & & tetacr ,concin , & & dzduu ,dzdvv ,ubot ,tauadd ,sus , & & bed ,susw ,bedw ,espir ,wave , & & scour ,ubot_from_com ,camax ,eps , & & iform ,par ,numintpar ,numrealpar,numstrpar , & & dllfunc ,dllhandle ,intpar ,realpar ,strpar , & !output: & aks ,caks ,taurat ,seddif ,rsedeq , & & kmaxsd ,conc2d ,sbcu ,sbcv ,sbwu , & & sbwv ,sswu ,sswv ,dss ,caks_ss3d , & & aks_ss3d ,ust2 ,t_relax ,error ) !----- GPL --------------------------------------------------------------------- ! ! Copyright (C) Stichting Deltares, 2011-2017. ! ! This program is free software: you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by ! the Free Software Foundation version 3. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License ! along with this program. If not, see . ! ! contact: delft3d.support@deltares.nl ! Stichting Deltares ! P.O. Box 177 ! 2600 MH Delft, The Netherlands ! ! All indications and logos of, and references to, "Delft3D" and "Deltares" ! are registered trademarks of Stichting Deltares, and remain the property of ! Stichting Deltares. All rights reserved. ! !------------------------------------------------------------------------------- ! $Id: eqtran.f90 7398 2017-07-19 14:14:21Z jagers $ ! $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/7545/src/utils_gpl/morphology/packages/morphology_kernel/src/eqtran.f90 $ !!--description----------------------------------------------------------------- ! ! !!--pseudo code and references-------------------------------------------------- ! NONE !!--declarations---------------------------------------------------------------- use precision use message_module, only: write_error use mathconsts, only: pi, ee use iso_c_binding, only: c_char use morphology_data_module ! implicit none ! ! Call variables ! integer(pntrsize) , intent(in) :: dllhandle integer , intent(in) :: i2d3d integer , intent(in) :: iform integer , intent(in) :: kmax ! Description and declaration in esm_alloc_int.f90 integer , intent(in) :: lsecfl ! Description and declaration in esm_alloc_int.f90 integer , intent(in) :: ltur ! Description and declaration in esm_alloc_int.f90 integer , intent(in) :: lundia ! Description and declaration in inout.igs integer , intent(in) :: numintpar integer , intent(in) :: numrealpar integer , intent(in) :: numstrpar integer , dimension(numintpar) , intent(inout) :: intpar real(fp) , intent(in) :: bed real(fp) , intent(in) :: bedw real(fp) , intent(in) :: camax real(fp) , dimension(kmax) , intent(inout) :: concin real(fp) , dimension(0:kmax) , intent(in) :: dicww ! Description and declaration in esm_alloc_real.f90 real(fp) , intent(in) :: dzduu ! Description and declaration in esm_alloc_real.f90 real(fp) , intent(in) :: dzdvv ! Description and declaration in esm_alloc_real.f90 real(fp) , intent(in) :: eps real(fp) , intent(in) :: espir real(fp) , intent(in) :: frac ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(30) , intent(inout) :: par real(fp) , intent(in) :: rksrs ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(kmax) , intent(in) :: sig ! Description and declaration in esm_alloc_real.f90 real(fp) , intent(in) :: sigmol ! Description and declaration in esm_alloc_real.f90 real(fp) , intent(in) :: spirint ! Spiral flow intensity real(fp) , intent(in) :: sus real(fp) , intent(in) :: susw real(fp) , intent(in) :: tauadd real(fp) , intent(in) :: taucr0 real(fp) , intent(in) :: tetacr real(fp) , dimension(kmax) , intent(in) :: thick ! Description and declaration in esm_alloc_real.f90 real(fp) , intent(in) :: ubot ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(0:kmax) , intent(in) :: ws ! Description and declaration in esm_alloc_real.f90 real(hp) , dimension(numrealpar), intent(inout) :: realpar logical , intent(in) :: scour logical , intent(in) :: suspfrac ! suspended sediment fraction logical , intent(in) :: ubot_from_com logical , intent(in) :: wave character(256) , intent(in) :: dllfunc character(256), dimension(numstrpar), intent(inout) :: strpar ! logical , intent(out) :: error integer , intent(out) :: kmaxsd real(fp) , intent(inout):: aks ! out parameter for Van Rijn, in parameter for others real(fp) , intent(out) :: aks_ss3d real(fp) , intent(out) :: caks real(fp) , intent(out) :: caks_ss3d real(fp) , intent(out) :: conc2d real(fp) , intent(out) :: dss real(fp), dimension(kmax) , intent(out) :: rsedeq real(fp) , intent(out) :: sbcu real(fp) , intent(out) :: sbcv real(fp) , intent(out) :: sbwu real(fp) , intent(out) :: sbwv real(fp), dimension(0:kmax) , intent(out) :: seddif real(fp) , intent(out) :: sswu real(fp) , intent(out) :: sswv real(fp) , intent(out) :: t_relax real(fp) , intent(out) :: taurat real(fp) , intent(out) :: ust2 ! ! Local variables ! integer :: i integer(pntrsize) :: ierror_ptr integer :: k integer(pntrsize), external :: perf_function_eqtran real(fp) :: ag real(fp) :: alphaspir real(fp) :: avgu real(fp) :: bakdif real(fp) :: cesus real(fp) :: chezy real(fp) :: cosa real(fp) :: d10 real(fp) :: d15 real(fp) :: d90 real(fp) :: dg real(fp) :: dgsd real(fp) :: di50 real(fp) :: dstar real(fp) :: fac3d2d real(fp) :: h1 real(fp) :: hidexp real(fp) :: hrms real(fp) :: mudfrac real(fp) :: rhosol real(fp) :: rhowat real(fp) :: rlabda real(fp) :: sag real(fp) :: salinity real(fp) :: sandfrac real(fp) :: sbot real(fp) :: sina real(fp) :: ssus real(fp) :: sscu real(fp) :: sscv real(fp) :: taub real(fp) :: teta real(fp) :: tp real(fp) :: txg real(fp) :: tyg real(fp) :: u real(fp) :: umod real(fp) :: uorb real(fp) :: ustarc real(fp) :: utot real(fp) :: uuu real(fp) :: v real(fp) :: vicmol real(fp) :: vonkar real(fp) :: vvv real(fp) :: z0cur real(fp) :: z0rou real(fp) :: zumod real(fp) :: kwtur real(fp) :: dzbdt real(fp) :: dzdx real(fp) :: dzdy real(fp) :: poros real(fp) :: ua real(fp) :: va real(fp) :: uamg ! ! Interface to dll is in High precision! ! real(hp) :: cesus_dll real(hp) :: sbc_dll real(hp) :: sbcu_dll real(hp) :: sbcv_dll real(hp) :: sbwu_dll real(hp) :: sbwv_dll real(hp) :: ssus_dll real(hp) :: sswu_dll real(hp) :: sswv_dll real(hp) :: t_relax_dll character(1024) :: errmsg character(256) :: message ! Contains message from internal or external transport formula character(kind=c_char) :: message_c(257) ! C- version of "message", including C_NULL_CHAR ! Calling perf_function_eqtran with "message" caused problems ! Solved by using "message_c" logical :: equi_conc ! equilibrium concentration given (instead of susp. transport rate) logical :: sbc_total ! total bed load given (instead of m,n components) logical :: sus_total ! total suspended load given (instead of m,n components) ! !! executable statements ------------------------------------------------------- ! ierror_ptr = 0 error = .false. ! utot = real(realpar(RP_EFVLM),fp) u = real(realpar(RP_EFUMN),fp) v = real(realpar(RP_EFVMN),fp) uuu = real(realpar(RP_UCHAR),fp) vvv = real(realpar(RP_VCHAR),fp) umod = real(realpar(RP_VELCH),fp) zumod = real(realpar(RP_ZVLCH),fp) h1 = real(realpar(RP_DEPTH),fp) chezy = real(realpar(RP_CHEZY),fp) hrms = real(realpar(RP_HRMS) ,fp) tp = real(realpar(RP_TPEAK),fp) teta = real(realpar(RP_TETA) ,fp) rlabda = real(realpar(RP_RLAMB),fp) uorb = real(realpar(RP_UORB) ,fp) kwtur = real(realpar(RP_KWTUR),fp) dzbdt = real(realpar(RP_BLCHG),fp) dzdx = real(realpar(RP_DZDX) ,fp) dzdy = real(realpar(RP_DZDY) ,fp) di50 = real(realpar(RP_D50) ,fp) dss = real(realpar(RP_DSS) ,fp) dstar = real(realpar(RP_DSTAR),fp) d10 = real(realpar(RP_D10MX),fp) d15 = real(realpar(RP_D15MX),fp) d90 = real(realpar(RP_D90MX),fp) mudfrac = real(realpar(RP_MUDFR),fp) hidexp = real(realpar(RP_HIDEX),fp) !ws = real(realpar(RP_SETVL),fp) rhosol = real(realpar(RP_RHOSL),fp) rhowat = real(realpar(RP_RHOWT),fp) salinity = real(realpar(RP_SALIN),fp) !temp = real(realpar(RP_TEMP) ,fp) ag = real(realpar(RP_GRAV) ,fp) vicmol = real(realpar(RP_VICML),fp) taub = real(realpar(RP_TAUB) ,fp) !ubed = real(realpar(RP_UBED ),fp) !vbed = real(realpar(RP_VBED ),fp) !velb = real(realpar(RP_VELBD),fp) !zvelb = real(realpar(RP_ZVLBD),fp) vonkar = real(realpar(RP_VNKAR),fp) z0cur = real(realpar(RP_Z0CUR),fp) z0rou = real(realpar(RP_Z0ROU),fp) ustarc = real(realpar(RP_USTAR),fp) dg = real(realpar(RP_DG) ,fp) dgsd = real(realpar(RP_DGSD) ,fp) sandfrac = real(realpar(RP_SNDFR),fp) ! cesus = 0.0_fp sbot = 0.0_fp sbcu = 0.0_fp sbcv = 0.0_fp ssus = 0.0_fp sscu = 0.0_fp sscv = 0.0_fp sbwu = 0.0_fp sbwv = 0.0_fp sswu = 0.0_fp sswv = 0.0_fp ua = 0.0_fp va = 0.0_fp uamg = 0.0_fp sag = sqrt(ag) bakdif = vicmol / sigmol ! if (iform > 0) then if (suspfrac) then ! ! By default entrainment and deposition (re)move sediment into/from ! the bottom-most layer. ! kmaxsd = kmax ! ! Use diffusivity of turbulence model as vertical sediment diffusion ! coefficient. In the future, we may OPTIONALLY enable Van Rijn's ! analytical 1984, 1993 or 2004 formulations here. ! do k = 0, kmax seddif(k) = dicww(k) enddo ! seddif(kmax) = vonkar*z0rou*ustarc endif endif ! equi_conc = .false. sbc_total = .false. sus_total = .false. if (iform == -1) then ! ! Van Rijn 1993 ! call tram1(numrealpar,realpar ,wave ,par , & & kmax ,bed , & & tauadd ,taucr0 ,aks ,eps ,camax , & & frac ,sig ,thick ,ws , & & dicww ,ltur , & & kmaxsd ,taurat ,caks , & & seddif ,sigmol ,rsedeq ,scour ,bedw , & & susw ,sbcu ,sbcv ,sbwu ,sbwv , & & sswu ,sswv ,conc2d ,error , & & message ) ! if (error) then call write_error(message, unit=lundia) return endif ! caks_ss3d = caks aks_ss3d = aks ! sbc_total = .false. sus_total = .false. elseif (iform == -2) then ! ! Van Rijn 2004 ! call tram2(numrealpar,realpar ,wave ,i2d3d ,par , & & kmax ,bed ,dzduu ,dzdvv ,rksrs , & & tauadd ,taucr0 ,aks ,eps ,camax , & & frac ,sig ,thick ,ws , & & dicww ,ltur ,aks_ss3d , & & kmaxsd ,taurat ,caks ,caks_ss3d ,concin , & & seddif ,sigmol ,rsedeq ,scour ,bedw , & & susw ,sbcu ,sbcv ,sbwu ,sbwv , & & sswu ,sswv ,tetacr ,conc2d ,error , & & message ) ! if (error) then call write_error(message, unit=lundia) return endif ! sbc_total = .false. sus_total = .false. elseif (iform == 1) then ! ! Engelund-Hansen ! call tranb1(utot ,di50 ,chezy ,h1 ,par , & & sbot ,ssus ) ! sbc_total = .true. sus_total = .true. elseif (iform == 2) then ! ! Meyer-Peter-Muller ! call tranb2(utot ,di50 ,d90 ,chezy ,h1 , & & par ,hidexp ,sbot ,ssus ) ! sbc_total = .true. sus_total = .true. elseif (iform == 3) then ! ! Ackers-White ! call tranb3(utot ,d90 ,chezy ,h1 ,par , & & sbot ,ssus ) ! sbc_total = .true. sus_total = .true. elseif (iform == 4) then ! ! general relation for bed load ! call tranb4(utot ,di50 ,chezy ,par ,hidexp , & & sbot ,ssus ) ! sbc_total = .true. sus_total = .true. elseif (iform == 5) then ! ! Bijker ! call tranb5(u ,v ,di50 ,d90 ,chezy , & & h1 ,hrms ,tp ,teta ,par , & & dzduu ,dzdvv ,sbcu ,sbcv ,sscu , & & sscv ,cesus ,vonkar ) ! sbc_total = .false. sus_total = .false. elseif (iform == 6) then ! ! Bailard ! errmsg = 'Bailard method is disabled' call write_error(errmsg, unit=lundia) error = .true. return ! !call tranb6(utot ,u ,v ,chezy ,h1 , & ! & hrms ,tp ,teta ,diss ,dzduu , & ! & dzdvv ,par ,sbcu ,sbcv ,sscu , & ! & sscv ) ! sbc_total = .false. sus_total = .false. elseif (iform == 7) then ! ! Van Rijn (1984, modified) ! call tranb7(utot ,di50 ,d90 ,h1 ,par , & & sbot ,ssus ,vonkar ,mudfrac) ! sbc_total = .true. sus_total = .true. elseif (iform == 8) then ! ! Van Rijn / Ribberink (1994) ! errmsg = 'Van Rijn/Ribberink (1994) method is disabled' call write_error(errmsg, unit=lundia) error = .true. return ! !call tranb8(u ,v ,hrms ,h1 ,teta , & ! & tp ,di50 ,d90 ,diss ,dzduu , & ! & dzdvv ,par ,sbcu ,sbcv ,sscu , & ! & sscv ) ! sbc_total = .false. sus_total = .false. elseif (iform == 9) then ! ! Silt module ! errmsg = 'Original Delft3D-MOR Silt module is disabled' call write_error(errmsg, unit=lundia) error = .true. return ! !call tranb9(utot ,h1 ,alfs ,sbot ,ssus ) ! sbc_total = .true. sus_total = .true. elseif (iform == 10) then ! ! Ashida and Michiue ! errmsg = 'Ashida and Michiue method is disabled' call write_error(errmsg, unit=lundia) error = .true. return ! !call trab10(utot ,di50 ,chezy ,h1 ,cosa , & ! & sina ,dzduu ,dzdvv ,par ,sbot , & ! & ssus ) ! sbc_total = .true. sus_total = .true. elseif (iform == 11) then ! ! Soulsby and Van Rijn ! call trab11(u ,v ,hrms ,h1 ,tp , & & di50 ,par ,sbcu ,sbcv ,sscu , & & sscv ,ubot ,vonkar ,ubot_from_com ) ! sbc_total = .false. sus_total = .false. elseif (iform == 12) then ! ! Soulsby ! call trab12(u ,v ,hrms ,h1 ,tp , & & teta ,di50 ,par ,sbcu ,sbcv , & & sscu ,sscv ,ubot ,vonkar ,ubot_from_com) ! sbc_total = .false. sus_total = .false. elseif (iform == 13) then ! ! test transport (Wang) Fredsoe ! call tran9t(utot ,di50 ,d90 ,chezy ,h1 , & & ustarc ,par ,sbot ,ssus ) ! sbc_total = .true. sus_total = .true. elseif (iform == 14) then ! ! generalized Ashida and Michiue ! call trab14(utot ,di50 ,chezy ,par ,hidexp , & & sbot ,ssus ) ! sbc_total = .true. sus_total = .true. elseif (iform == 16) then ! ! Wilcock & Crowe ! call trabwc(utot ,di50 ,taub ,par ,sbot , & & ssus ,dg ,sandfrac ,chezy ) ! sbc_total = .true. sus_total = .true. elseif (iform == 17) then ! ! Modified Wilcock & Crowe ! call trabwc2(utot ,di50 ,taub ,par ,sbot , & & ssus ,dg ,dgsd ,chezy ) ! sbc_total = .true. sus_total = .true. elseif (iform == 18) then ! ! Gaeuman et al. (development of Wilcock & Crowe) ! call trabg(utot ,di50 ,taub ,par ,sbot , & & ssus ,dg ,dgsd ,chezy ) ! sbc_total = .true. sus_total = .true. elseif (iform == 19) then ! ! van Thiel / Van Rijn (2008) ! call trab19(u ,v ,hrms ,rlabda ,teta ,h1 ,tp , & & di50 ,d15 ,d90 ,par ,dzbdt ,vicmol ,poros , & & chezy ,dzdx ,dzdy ,sbcu ,sbcv ,sscu ,sscv , & & ua ,va ,ubot ,kwtur ,vonkar ,ubot_from_com ) ! uamg = sqrt(ua*ua+va*va) realpar(RP_UAU) = real(ua ,hp) ! needed for suspended transport realpar(RP_VAU) = real(va ,hp) ! sbc_total = .false. sus_total = .false. elseif (iform == 20) then ! ! Soulsby / Van Rijn with XBeach adaptations ! call trab20(u ,v ,hrms ,rlabda ,teta ,h1 ,tp , & & di50 ,d15 ,d90 ,par ,dzbdt ,vicmol ,poros , & & chezy ,dzdx ,dzdy ,sbcu ,sbcv ,sscu ,sscv , & & ua ,va ,ubot ,kwtur ,vonkar ,ubot_from_com ) ! uamg = sqrt(ua*ua+va*va) realpar(RP_UAU) = real(ua ,hp) ! needed for suspended transport realpar(RP_VAU) = real(va ,hp) ! sbc_total = .false. sus_total = .false. elseif (iform == 21) then ! ! Nino and Garcia (1998) ! call tranb21(utot, di50, d90, h1, par, sbot, ssus) ! sbc_total = .true. sus_total = .true. elseif (iform == 15) then ! ! User defined formula in DLL ! Input parameters are passed via realpar/intpar/strpar-arrays that have ! been mostly filled in calling routine ! realpar(RP_UMEAN) = real(u ,hp) realpar(RP_VMEAN) = real(v ,hp) realpar(RP_VELMN) = real(utot ,hp) ! Initialisation of output variables of user defined transport formulae ! sbc_total = .false. ! may be changed by user defined formulae sus_total = .true. ! always true for user defined formulae sbc_dll = 0.0_hp sbcu_dll = 0.0_hp sbcv_dll = 0.0_hp sbwu_dll = 0.0_hp sbwv_dll = 0.0_hp ! equi_conc = .false. cesus_dll = 0.0_hp ssus_dll = 0.0_hp sswu_dll = 0.0_hp sswv_dll = 0.0_hp t_relax_dll = 0.0_hp message = ' ' do i=1,256 message_c(i) = message(i:i) enddo message_c(257) = C_NULL_CHAR ! ! psem/vsem is used to be sure this works fine in DD calculations ! call psemlun ierror_ptr = perf_function_eqtran(dllhandle , dllfunc , & intpar , numintpar , & realpar , numrealpar , & strpar , numstrpar , & sbc_total, sbc_dll , sbcu_dll , & sbcv_dll , sbwu_dll , sbwv_dll , & equi_conc, cesus_dll, ssus_dll , & sswu_dll , sswv_dll , t_relax_dll , & message_c) message = transfer(message_c(1:256), message) call vsemlun if (ierror_ptr /= 0) then errmsg = 'Cannot find function "'//trim(dllfunc)//'" in dynamic library.' call write_error(errmsg, unit=lundia) error = .true. return endif if (message /= ' ') then errmsg = 'Message from user defined transport formulae '//trim(dllfunc)//' :' call write_error(errmsg, unit=lundia) write (lundia,'(a,a )') ' ', trim(message) write (lundia,'(a )') ' ' error = .true. return endif ! ! Output parameters ! sbot = real(sbc_dll ,fp) sbcu = real(sbcu_dll,fp) sbcv = real(sbcv_dll,fp) sbwu = real(sbwu_dll,fp) sbwv = real(sbwv_dll,fp) ! cesus = real(cesus_dll,fp) ssus = real(ssus_dll ,fp) sswu = real(sswu_dll ,fp) sswv = real(sswv_dll ,fp) ! t_relax = real(t_relax_dll,fp) else errmsg = 'Transport formula not recognized' call write_error(errmsg, unit=lundia) error = .true. return endif ! if (iform > 0) then ! ! Change from volume to mass concentrations/transport rates. ! ! Van Rijn 1993 or Van Rijn 2004 include rhosol in the transport ! rates, the other formulae don't. So, multiply the transport ! rates and concentrations by rhosol now. ! if (sbc_total) then sbot = sbot * rhosol else sbcu = sbcu * rhosol sbcv = sbcv * rhosol endif ! cesus = cesus * rhosol if (sus_total) then ssus = ssus * rhosol else sscu = sscu * rhosol sscv = sscv * rhosol ssus = sqrt(sscu**2 + sscv**2) endif ! sbwu = sbwu * rhosol sbwv = sbwv * rhosol sswu = sswu * rhosol sswv = sswv * rhosol endif ! ! If only bed load transport magnitude is given, then the bed load ! should be oriented based on the near bed velocity. The near bed ! velocity is given by either the velocity (uuu,vvv) in 3D or ! the depth averaged velocity (uuu,vvv) corrected for spiral flow in 2D. ! ust2 = ustarc**2 if (sbc_total) then if (umod > 0.0_fp) then ! ! Correct bed load transport direction for spiral flow intensity ! if (lsecfl == 0) then alphaspir = 0.0_fp else alphaspir = sqrt(ag) / 0.4_fp / chezy alphaspir = 12.5_fp * espir * (1.0_fp-alphaspir) ! 12.5 = 2.0_fp/(von Karman)^2 alphaspir = alphaspir * spirint / umod endif txg = ust2 * (uuu + alphaspir*vvv) / umod tyg = ust2 * (vvv - alphaspir*uuu) / umod ust2 = sqrt(txg**2 + tyg**2) if (ust2 > eps) then cosa = txg / ust2 sina = tyg / ust2 else cosa = 0.0_fp sina = 0.0_fp endif sbcu = sbot * cosa sbcv = sbot * sina else sbcu = 0.0_fp sbcv = 0.0_fp endif endif ! ! Adjust for calibration factors ! sbcu = bed * sbcu sbcv = bed * sbcv sbwu = bedw * sbwu sbwv = bedw * sbwv sswu = susw * sswu sswv = susw * sswv ! if (iform>0) then cesus = sus * cesus ssus = sus * ssus sscu = sus * sscu sscv = sus * sscv else caks = sus * caks caks_ss3d= sus * caks_ss3d do k = 1, kmax rsedeq(k) = sus * rsedeq(k) enddo conc2d = sus * conc2d endif ! if (suspfrac) then ! ! For this fraction a 2D or 3D advection-diffusion will be solved. ! In 2D case we need to provide: conc2d ! In 3D case we need to provide: aks and caks ! if (iform <= 0) then ! ! Van Rijn 1993 or 2004 formula ! NOTE: conc2d, aks, and caks have already been computed in Van Rijn ! specific manner. ! else ! ! If we are not using Van Rijn 1993 or Van Rijn 2004 ! then we still need to provide values for conc2d, aks, and caks. ! if (equi_conc) then ! ! Concentration given by transport formula ! else ! ! Suspended transport rate given by transport formula, ! derive concentration. Add non-zero uamg for van Thiel ! and XBeach like Soulsby / Van Rijn (iform=19, 20) ! cesus = ssus / (utot+uamg+eps) / h1 endif ! ! Concentration needs to be multiplied by frac to match Van Rijn ! formulae. ! conc2d = cesus * frac ! ! Convert depth averaged concentration to reference concentration ! at distance aks from bed. ! kmaxsd = kmax call factor3d2d(kmax ,aks ,kmaxsd ,sig , & & thick ,seddif ,ws ,bakdif , & & z0rou ,h1 ,fac3d2d ) caks = conc2d / (fac3d2d+eps) / rhosol rsedeq = 0.0_fp ! aks_ss3d = aks caks_ss3d = caks endif else ! .not.suspfrac ! ! Note: in case of bedload sediment type, the suspended load is added to the ! bed load to compute the total load. Bed slope effect will be applied to ! both the bed- and the suspended part. ! ! The effect of frac will be included in the bed load at a later stage. ! if (iform <= 0) then ! ! Van Rijn 1993 or 2004 formula ! NOTE: conc2d is based on caks which has been multiplied by ! frac above. Since sbcu/v will be multiplied by frac later (again) ! there is the risk of double correction. So, undo the multiplication. ! if (frac>0.0_fp) then sbcu = sbcu + conc2d * h1 * u / frac sbcv = sbcv + conc2d * h1 * v / frac endif elseif (equi_conc) then ! ! Concentration given by transport formula ! sbcu = sbcu + cesus * h1 * u sbcv = sbcv + cesus * h1 * v elseif (sus_total) then ! ! Total suspended transport rate given by transport formula, ! assume that it is transported in depth-averaged direction ! sbcu = sbcu + ssus * u / (utot+eps) sbcv = sbcv + ssus * v / (utot+eps) else ! ! Suspended transport rate components given by transport ! formula, add them individually to the bed load ! sbcu = sbcu + sscu sbcv = sbcv + sscv endif endif end subroutine eqtran