subroutine secrhs(s0 ,s1 ,dps ,u1 ,v1 , & & guu ,gvv ,gsqs ,j ,nmmaxj , & & nmmax ,kmax ,lstsci ,lsecfl ,icx , & & icy ,kfu ,kfv ,kfs ,kcs , & & xcor ,ycor ,sour ,sink ,cfurou , & & cfvrou ,fcorio ,curstr ,x3 ,x2y , & & xy2 ,y3 ,gdp ) !----- GPL --------------------------------------------------------------------- ! ! Copyright (C) Stichting Deltares, 2011-2014. ! ! 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$ ! $HeadURL$ !!--description----------------------------------------------------------------- ! ! Function: Computes righthandside terms for the transport ! equation of spiral motion intensity (secondary ! flow). ! 3-1-94 CF formula Ta aangepast. ! Method used: ! !!--pseudo code and references-------------------------------------------------- ! NONE !!--declarations---------------------------------------------------------------- use precision ! use globaldata ! implicit none ! type(globdat),target :: gdp ! ! The following list of pointer parameters is used to point inside the gdp structure ! real(fp) , pointer :: ag real(fp) , pointer :: vonkar integer , pointer :: lundia real(fp) , pointer :: chzmin real(fp) , pointer :: dryflc real(fp) , pointer :: rmincf ! ! Global variables ! integer :: icx !! Increment in the x-dir., if icx= nmax then computation proceeds in the x-dir. !! if icx=1 then computation proceeds in the y-dir. integer :: icy !! Increment in the y-dir. (see icx) integer :: j !! Begin pointer for arrays which have been transformed into 1d arrays. !! due to the shift in the 2nd (m-)index, j = -2*nmax + 1 integer :: kmax ! Description and declaration in esm_alloc_int.f90 integer , intent(in) :: lsecfl ! Description and declaration in dimens.igs integer , intent(in) :: lstsci ! Description and declaration in esm_alloc_int.f90 integer :: nmmax ! Description and declaration in dimens.igs integer :: nmmaxj ! Description and declaration in dimens.igs integer , dimension(gdp%d%nmlb:gdp%d%nmub) :: kcs ! Description and declaration in esm_alloc_int.f90 integer , dimension(gdp%d%nmlb:gdp%d%nmub) :: kfs ! Description and declaration in esm_alloc_int.f90 integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfu ! Description and declaration in esm_alloc_int.f90 integer , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: kfv ! Description and declaration in esm_alloc_int.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: curstr !! Internal work array real(prec), dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: dps ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: fcorio ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: gsqs ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: guu ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: gvv ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: s0 ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: s1 ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: x2y ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: x3 ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: xcor ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: xy2 ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: y3 ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: ycor ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, 3) , intent(in) :: cfurou ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, 3) , intent(in) :: cfvrou ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: u1 ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: v1 ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci), intent(out) :: sink ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci), intent(out) :: sour ! Description and declaration in esm_alloc_real.f90 ! ! ! Local variables ! integer :: ierr ! Error counter integer :: iexit ! Exit code when ALPHA < 0.5 integer :: k integer :: ndm integer :: nm integer :: nmd real(fp) :: alpha real(fp) :: bendeq ! Bend equilibrium for spiral motion real(fp) :: chezy ! Chezy coefficient at zeta point real(fp) :: corieq ! Coriolis equilibrium for spiral m. real(fp) :: geta2 real(fp) :: gksi2 real(fp) :: h0new ! Total waterdepth at new time real(fp) :: h0old ! Total waterdepth at old time real(fp) :: htrsh real(fp) :: riv real(fp) :: rminiv real(fp) :: tai real(fp) :: tanew real(fp) :: taold real(fp) :: umod ! Sqrt(uuu*uuu+vvv*vvv) real(fp) :: uuu ! U-velocity at zeta point real(fp) :: vvv ! V-velocity at zeta point character(15) :: errmsg ! Text string for error message ! ! !! executable statements ------------------------------------------------------- ! ! ! chzmin => gdp%gdnumeco%chzmin dryflc => gdp%gdnumeco%dryflc rmincf => gdp%gdnumeco%rmincf lundia => gdp%gdinout%lundia ag => gdp%gdphysco%ag vonkar => gdp%gdphysco%vonkar ! htrsh = 0.5*dryflc k = 1 ! ! compute curvature of streakline ! call curvat(u1 ,v1 ,gsqs ,guu ,gvv , & & j ,nmmaxj ,nmmax ,kmax ,icx , & & icy ,kcs ,kfs ,curstr ,x3 , & & x2y ,xy2 ,y3 ,gdp ) ! ! compute source and sink terms ! For SOUR the old time and for SINK the new time ! ierr = 0 nmd = -icx ndm = -icy ! do nm = 1, nmmax nmd = nmd + 1 ndm = ndm + 1 ! !-------Actual point should be active ! if ( (kfs(nm)==1) .and. (kcs(nm)==1) ) then h0old = max(htrsh, s0(nm) + real(dps(nm),fp)) h0new = max(htrsh, s1(nm) + real(dps(nm),fp)) uuu = 0.5*(u1(nm, k) + u1(nmd, k)) vvv = 0.5*(v1(nm, k) + v1(ndm, k)) umod = sqrt(uuu*uuu + vvv*vvv) gksi2 = gvv(nm) + gvv(ndm) geta2 = guu(nm) + guu(nmd) rminiv = 1.0/(rmincf*max(gksi2, geta2)) if (abs(curstr(nm))= .5 else SOURC & SINK become < 0. ! because TAH0 = 0.5*(1.-2.*ALPHA)*H0 should be > 0 ! if (chezy>chzmin) then alpha = sqrt(ag)/(vonkar*chezy) else alpha = sqrt(ag)/(vonkar*chzmin) endif if (alpha>=0.5) then ierr = ierr + 1 endif taold = 0.5*(1. - 2.*alpha)*h0old tanew = 0.5*(1. - 2.*alpha)*h0new tai = ((vonkar**2*alpha)*umod)/taold ! !---------Subtract Coriolis contribution ! sour(nm, k, lsecfl) = (bendeq - corieq)*tai tai = ((vonkar**2*alpha)*umod)/tanew sink(nm, k, lsecfl) = 1.*tai endif enddo ! !-----Number of errors ! if (ierr>0) then write (errmsg, '(a,i3,a)') 'in ', ierr, ' point(s)' call prterr(lundia ,'S230' ,errmsg ) ! ! !-------stop routine for DELFT3D ! iexit = 4 call d3stop(iexit ,gdp ) ! endif end subroutine secrhs