subroutine sud(dischy ,nst ,icreep ,betac ,mmax , & & nmax ,j ,nmmaxj ,nmmax ,kmax , & & lstsci ,nsrc ,lsecfl ,norow ,icx , & & icy ,dismmt ,irocol ,mnksrc , & & kfu ,kfv ,kfs ,kcs ,kspu , & & kadu ,kadv ,kcu ,kfumin ,kfumax , & & porosu ,s0 ,s1 ,u0 ,u1 , & & v1 ,w1 ,r0 ,qxk ,qyk , & & qzk ,guu ,gvv ,gvu ,gsqs , & & gud ,gvd ,gvz ,gsqiu ,dteu , & & circ2d ,circ3d ,disch , & & umdis ,umean ,hu ,hv ,dpu ,dzu1 , & & dpdksi ,thick ,sig ,dps ,taubpu , & & taubsu ,rho ,sumrho ,wsu ,fxw , & & wsbodyu ,idry ,crbc ,vicuv ,hu0 , & & vnu2d ,vicww ,rxx ,rxy ,dfu , & & deltau ,tp ,rlabda ,cfurou ,cfvrou , & & rttfu ,diapl ,rnpl , & & windsu ,patm ,fcorio ,evap ,ubrlsu , & & uwtypu ,hkru ,pship ,tgfsep ,a , & & b ,c ,d ,aa ,bb , & & cc ,dd ,tetau ,aak ,bbk , & & cck ,ddk ,d0 ,d0k ,bbka , & & bbkc ,ua ,ub ,soumud ,dis_nf , & & precip ,ustokes ,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: SUD evaluates/solves the implicitly coupled ! momentum and continuity equation at each ! half time step. ! Special approximation for pressure term, ! based on limiter to avoid artificial flow. ! Switch which makes it possible to use ! upwind-approach for wet cross section in shallow ! areas or if the model area contains structures. ! Method used: A.D.I.-scheme is used. ! Upwind-approach for wet cross section in shallow ! areas or if the model area contains structures. ! At 2D Weir points: ! - depth value DPU is not corrected like in general ! (3D) weir case (see CALDPU) ! - crest height is explicitly taken into account ! in drying check ! - 2D Turbulence model at depth points ! ! UA and UB are workarrays WRKB15 and WRKB16; Used only in CUCNP ! !!--pseudo code and references-------------------------------------------------- ! NONE !!--declarations---------------------------------------------------------------- use precision use flow2d3d_timers use globaldata use dfparall ! implicit none ! type(globdat),target :: gdp ! ! The following list of pointer parameters is used to point inside the gdp structure ! include 'flow_steps_f.inc' real(fp) , pointer :: eps integer , pointer :: maseva integer , pointer :: lundia integer , pointer :: ntstep real(fp) , pointer :: dco real(fp) , pointer :: drycrt real(fp) , pointer :: dryflc real(fp) , pointer :: hdt integer , pointer :: iter1 character(6) , pointer :: momsol character(8) , pointer :: dpuopt real(fp) , pointer :: rhow real(fp) , pointer :: ag integer , pointer :: iro logical , pointer :: wind logical , pointer :: culvert logical , pointer :: mudlay logical , pointer :: nfl logical , pointer :: zmodel logical , pointer :: wavcmp ! ! Global variables ! integer :: icreep ! Description and declaration in tricom.igs integer , intent(in) :: 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 , intent(in) :: icy !! Increment in the Y-dir. (see ICX) integer :: idry 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 :: lsecfl ! Description and declaration in dimens.igs integer :: lstsci ! Description and declaration in esm_alloc_int.f90 integer , intent(in) :: mmax ! Description and declaration in esm_alloc_int.f90 integer , intent(in) :: nmax ! 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 :: norow ! Description and declaration in esm_alloc_int.f90 integer :: nsrc ! Description and declaration in esm_alloc_int.f90 integer , intent(in) :: nst !! Time step number integer , dimension(5, norow) :: irocol ! Description and declaration in esm_alloc_int.f90 integer , dimension(7, nsrc) :: mnksrc ! Description and declaration in esm_alloc_int.f90 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) :: kcu ! 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) :: kfu ! Description and declaration in esm_alloc_int.f90 integer , dimension(gdp%d%nmlb:gdp%d%nmub) :: kfumax ! Description and declaration in esm_alloc_int.f90 integer , dimension(gdp%d%nmlb:gdp%d%nmub) :: kfumin ! Description and declaration in esm_alloc_int.f90 integer , dimension(gdp%d%nmlb:gdp%d%nmub) :: kfv ! Description and declaration in esm_alloc_int.f90 integer , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax) :: kspu ! Description and declaration in esm_alloc_int.f90 integer , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: kadu ! Description and declaration in esm_alloc_int.f90 integer , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: kadv ! Description and declaration in esm_alloc_int.f90 real(fp) :: betac ! Description and declaration in tricom.igs real(fp) , dimension(12, norow) :: crbc ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(4, norow) :: circ2d ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: a !! Internal work array, tridiagonal matrix water levels lower diagonal real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: aa !! Internal work array, coupling mean velocity with water level point in (N,M,K) left (down) real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: b !! Internal work array, tridiagonal matrix water levels main diagonal real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: bb !! Internal work array, coefficient mean velocity real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: c !! Internal work array, tridiagonal matrix water levels upper diagonal real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: cc !! Internal work array, coupling mean velocity with water level point right (upper) real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: d !! Internal work array, Right Hand side of the Continuity equation real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: d0 !! Internal work array, Explicit part of the Right Hand side real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: dd !! Internal work array, Right hand side of the momentum eq. at (N,M) real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: deltau ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: dfu ! Description and declaration in esm_alloc_real.f90 real(prec) , dimension(gdp%d%nmlb:gdp%d%nmub) :: dps ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: dpu ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: dteu ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: evap ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: fcorio ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: fxw ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: gsqiu ! 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) :: gud ! 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) :: gvd ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: gvu ! 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) :: gvz ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: hkru ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: hu0 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: hu ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: hv ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: patm ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: precip ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: pship ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: rlabda ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: s0 ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: s1 ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) , intent(in) :: soumud ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: taubpu ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: taubsu ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: tetau !! Factor for upwind approach S0 can be 0.0, 0.5 or 1.0 depending on value of HU, DCO, KSPU and UMEAN real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: tgfsep !! Water elev. induced by tide gen.force real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: tp ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: umean ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: uwtypu ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: vnu2d ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: windsu ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: wsu ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub) :: wsbodyu ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax) , intent(out) :: qzk ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax) :: vicww ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, 0:kmax) :: w1 ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, 3) :: cfurou ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, 3) :: cfvrou ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: aak !! Internal work array (in CUCNP & UZD) real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: bbk !! Internal work array (in CUCNP & UZD) real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: bbka !! Internal work array real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: bbkc !! Internal work array real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: cck !! Internal work array (in CUCNP & UZD) real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: d0k !! Internal work array real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: ddk !! Internal work array, diagonal space at (N,M,K) real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: diapl ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: dpdksi ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: dzu1 ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: porosu ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: qxk ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) , intent(in) :: qyk ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: rho ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: rnpl ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: rttfu ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: rxx ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: rxy ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: sumrho ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: u0 ! 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) :: ua real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: ub real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: ubrlsu ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: ustokes ! Description and declaration in trisol.igs 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+2) :: vicuv ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax) :: dis_nf ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(gdp%d%nmlb:gdp%d%nmub, kmax, lstsci) :: r0 ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(kmax) :: sig ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(kmax) :: thick ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(kmax, 2, norow) :: circ3d ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(nsrc) :: disch ! Description and declaration in esm_alloc_real.f90 real(fp) , dimension(nsrc) :: umdis ! Description and declaration in esm_alloc_real.f90 character(1) , dimension(nsrc) :: dismmt ! Description and declaration in esm_alloc_char.f90 character(8) :: dischy ! Description and declaration in tricom.igs ! ! Local variables ! integer :: ddb integer :: i integer :: icxy integer :: icol integer :: ierror integer :: intdir integer :: iter integer :: itr integer :: k integer :: kenm integer :: kk integer :: m integer :: mmaxddb integer :: n integer :: nhystp integer :: nm integer :: nmaxddb integer :: nmd integer :: nmf integer :: nmlu integer :: nmu logical :: error ! Flag for detection of closure error in mass-balance real(hp) :: bi real(hp) :: fac real(fp) :: epsomb real(fp) :: hdti real(fp) :: hnm real(fp) :: hucres real(fp) :: humean ! Mean value for H in U-points real(fp) :: drytrsh real(fp) :: pr real(fp) :: dxiu real(fp) :: dxid character(80) :: errtxt integer :: nm_pos ! indicating the array to be exchanged has nm index at the 2nd place, e.g., dbodsd(lsedtot,nm) ! !! executable statements ------------------------------------------------------- ! eps => gdp%gdconst%eps maseva => gdp%gdheat%maseva lundia => gdp%gdinout%lundia ntstep => gdp%gdinttim%ntstep dco => gdp%gdnumeco%dco drycrt => gdp%gdnumeco%drycrt dryflc => gdp%gdnumeco%dryflc hdt => gdp%gdnumeco%hdt iter1 => gdp%gdnumeco%iter1 momsol => gdp%gdnumeco%momsol dpuopt => gdp%gdnumeco%dpuopt rhow => gdp%gdphysco%rhow ag => gdp%gdphysco%ag iro => gdp%gdphysco%iro wind => gdp%gdprocs%wind culvert => gdp%gdprocs%culvert mudlay => gdp%gdprocs%mudlay nfl => gdp%gdprocs%nfl zmodel => gdp%gdprocs%zmodel wavcmp => gdp%gdprocs%wavcmp ! call timer_start(timer_sud_rest, gdp) ! ddb = gdp%d%ddbound nmaxddb = nmax + 2*gdp%d%ddbound mmaxddb = mmax + 2*gdp%d%ddbound hdti = 1.0_fp / hdt icxy = max(icx, icy) drytrsh = drycrt ! nm_pos = 1 ! if (idry == 1) then ! ! This is necessary because SUD can be repeated in case of drying ! in DRYCHK (DRYFLP <> NO) ! do nm = 1, nmmax umean(nm) = 0.0 if (kfu(nm)==1) then do k = 1, kmax umean(nm) = umean(nm) + thick(k)*u0(nm, k) enddo endif enddo call upwhu (j ,nmmaxj ,nmmax ,kmax ,icx , & & zmodel ,kcs ,kcu ,kspu ,dps , & & s0 ,dpu ,umean ,hu ,gdp ) endif do nm = 1, nmmax hu0(nm) = hu(nm) enddo ! ! AFTER CALCULATION OF HU FOR CUCNP THE VALUE SHOULD BE CORRECTED ! AND TETAU SHOULD BE CALCULATED FOR UPWIND APPROACH ! nmu = +icx do nm = 1, nmmax nmu = nmu + 1 hu(nm) = max(hu(nm), 0.01_fp) tetau(nm) = 0.5 if (kfu(nm)==1) then humean = 0.5*(s0(nm) + s0(nmu)) + dpu(nm) if (humean0 .or. dpuopt=='UPW' .or. momsol == 'flood ') then if (umean(nm)>=0.001) then tetau(nm) = 1.0 elseif (umean(nm)<= - 0.001) then tetau(nm) = 0.0 else tetau(nm) = 1.0 if (s0(nmu)>s0(nm)) tetau(nm) = 0.0 endif endif endif enddo call timer_stop(timer_sud_rest, gdp) ! call timer_start(timer_sud_cucnp, gdp) call cucnp(dischy ,icreep ,dpdksi ,s0 ,u0 , & & v1 ,w1 ,hu ,hv ,dps ,dpu , & & umean ,guu ,gvv ,gvu ,gsqs , & & gvd ,gud ,gvz ,gsqiu ,qxk , & & qyk ,disch ,umdis ,mnksrc ,dismmt ,j , & & nmmaxj ,nmmax ,kmax ,icx ,icy , & & nsrc ,lsecfl ,lstsci ,betac ,aak , & & bbk ,cck ,ddk ,bbka ,bbkc , & & thick ,sig ,rho ,sumrho ,vicuv , & & vnu2d ,vicww ,wsu ,fxw ,wsbodyu , & & rxx ,rxy ,kcs ,kcu ,kfu ,kfv , & & kfs ,kspu ,kadu ,kadv ,dfu ,deltau , & & tp ,rlabda ,cfurou ,cfvrou ,rttfu , & & r0 ,diapl ,rnpl ,taubpu ,taubsu , & & windsu ,patm ,fcorio ,ubrlsu ,uwtypu , & & hkru ,pship ,tgfsep ,dteu ,ua , & & ub ,ustokes ,.false. ,u1 ,s1 ,gdp ) call timer_stop(timer_sud_cucnp, gdp) ! ! INITIALISATION OF ITERATION OVER CONTINUITY EQUATION ! call timer_start(timer_sud_rest, gdp) do k = 1, kmax do nm = 1, nmmax if (kcs(nm) > 0) then d0k(nm, k) = gsqs(nm)*thick(k)*s0(nm)*hdti - qyk(nm, k) + qyk(nm - icy, k) endif enddo enddo ! ! IN LAYER 1 DUE TO PRECIPITATION/EVAPORATION ! FOR TIME DEPENDENT INPUT OR HEAT MODEL WITH SPECIAL REQUEST ! if (maseva>0) then do nm = 1, nmmax if (kcs(nm)==1) then d0k(nm, 1) = d0k(nm, 1) + precip(nm)*gsqs(nm) if (kfs(nm)==1) then d0k(nm, 1) = d0k(nm, 1) - (evap(nm)/rhow)*gsqs(nm) endif endif enddo endif ! ! ADDITION OF DISCHARGES (suction only permitted if the point isn't dry) ! do i = 1, nsrc nm = (mnksrc(5, i) + ddb) + ((mnksrc(4, i) - 1) + ddb)*icxy k = mnksrc(6, i) if (k .eq. -1) cycle kenm = min(1, kfu(nm) + kfu(nm - icx) + kfv(nm) + kfv(nm - icy)) if (kenm/=0 .or. disch(i)>=0.0) then if (k/=0) then d0k(nm, k) = d0k(nm, k) + disch(i) else do kk = 1, kmax d0k(nm, kk) = d0k(nm, kk) + disch(i)*thick(kk) enddo endif else write (errtxt, '(i0,i3)') nst, i call prterr(lundia ,'S208' ,trim(errtxt)) endif ! ! in case of an intake for an intake/outfall combination: ! if (mnksrc(7, i)>=2) then nm = (mnksrc(2, i) + ddb) + ((mnksrc(1, i) - 1) + ddb)*icxy k = mnksrc(3, i) kenm = min(1, kfu(nm) + kfu(nm - icx) + kfv(nm) + kfv(nm - icy)) if (kenm/=0 .or. -disch(i)>=0.0) then if (k/=0) then d0k(nm, k) = d0k(nm, k) - disch(i) else do kk = 1, kmax d0k(nm, kk) = d0k(nm, kk) - disch(i)*thick(kk) enddo endif ! ! in case of a culvert no warning generated ! elseif (mnksrc(7, i)/=3) then write (errtxt, '(i0,i3)') nst, i call prterr(lundia ,'S208' ,trim(errtxt)) else endif endif enddo ! ! ADDITION OF DISCHARGES from near field model ! if (nfl) then do nm = 1, nmmax do k = 1, kmax d0k(nm,k) = d0k(nm,k) + dis_nf(nm,k) enddo enddo endif ! ! add sources/sinks mud layer if mudlay == .true. ! if (mudlay) then do nm = 1, nmmax ! ! kfs mask array removed, since then the layer cannot increase ! any more after it has gone dry ! if (kcs(nm)==1) then d0k(nm, 1) = d0k(nm, 1) + gsqs(nm)*soumud(nm) endif enddo endif ! ! Initialise arrays a - dd for all (nm) ! a = 0.0 b = 1.0 c = 0.0 d0 = 0.0 aa = 0.0 bb = 1.0 cc = 0.0 dd = 0.0 ! do nm = 1, nmmax d(nm) = s0(nm) enddo ! do k = 1, kmax do nm = 1, nmmax if (kfu(nm)==1) then fac = real(porosu(nm,k),hp) * real(thick(k),hp) / real(bbk(nm,k),hp) aa(nm) = real(aa(nm),hp) + real(aak(nm,k),hp)*fac cc(nm) = real(cc(nm),hp) + real(cck(nm,k),hp)*fac dd(nm) = real(dd(nm),hp) + real(ddk(nm,k),hp)*fac endif enddo enddo do k = 1, kmax do nm = 1, nmmax if (kcs(nm)==1) d0(nm) = d0(nm) + d0k(nm, k) enddo enddo call timer_stop(timer_sud_rest, gdp) ! ! ITERATIVE LOOP OVER CURRENT ROW. CALCULATION OF CONTINUITY EQ. ! 9999 continue itr = 0 ! do iter = 1, iter1 ! ! BOUNDARY CONDITIONS ! call timer_start(timer_sud_cucbp, gdp) call cucbp(kmax ,norow ,icx , & & icy ,zmodel ,irocol ,kcs ,kfu , & & kfumin ,kfumax ,s0 ,u0 ,dpu , & & hu ,umean ,tetau ,guu ,gvu , & & dzu1 ,thick ,circ2d ,circ3d ,a , & & b ,c ,d ,aa ,bb , & & cc ,dd ,aak ,bbk ,cck , & & ddk ,crbc ,wavcmp ,gdp ) call timer_stop(timer_sud_cucbp, gdp) ! ! SET UP SYSTEM OF EQUATIONS FOR INTERIOR POINTS ! call timer_start(timer_sud_rest, gdp) if (momsol=='flood ') then nmd = -icx do nm = 1, nmmax nmd = nmd + 1 if (kcs(nm)==1) then dxid = hu0(nmd)*guu(nmd) dxiu = hu0(nm) *guu(nm) a(nm) = dxid*aa(nmd) b(nm) = hdti*gsqs(nm) + dxid*cc(nmd) - dxiu*aa(nm) c(nm) = -dxiu*cc(nm) d(nm) = d0(nm) - dxiu*dd(nm) + dxid*dd(nmd) endif enddo else nmd = -icx do nm = 1, nmmax nmd = nmd + 1 if (kcs(nm)==1) then a(nm) = guu(nmd)*(hu(nmd)*aa(nmd) - tetau(nmd)*dd(nmd)) b(nm) = hdti*gsqs(nm) & & + guu(nmd)*(hu(nmd)*cc(nmd) - (1.0 - tetau(nmd))*dd(nmd)) & & - guu(nm) *(hu(nm) *aa(nm) - tetau(nm) *dd(nm) ) c(nm) = -guu(nm)*(hu(nm)*cc(nm) - (1.0 - tetau(nm))*dd(nm)) d(nm) = d0(nm) - (guu(nm)*dpu(nm)*dd(nm) - guu(nmd)*dpu(nmd)*dd(nmd)) endif enddo endif call timer_stop(timer_sud_rest, gdp) ! ! Domain decomposition: ! Give Mapper chance to build the coupling equations ! Note that this is a two-stage process. First, coupling equations ! are built for the coupling points start+1;end-1 ! Secondly, coupling equations are built for the `end' coupling points, ! start and end. ! ! nhystp = nxtstp(d3dflow_build_adi_zeta, gdp) ! ! End of Domain decomposition addition ! ! !***SCALE ROWS OF MATRIX/RIGHT HAND SIDE VECTOR ! call timer_start(timer_sud_rowsc, gdp) do nm = 1, nmmax bi = 1.0_hp / real(b(nm),hp) a(nm) = real(a(nm),hp) * bi b(nm) = 1.0_fp c(nm) = real(c(nm),hp) * bi d(nm) = real(d(nm),hp) * bi enddo call timer_stop(timer_sud_rowsc, gdp) ! ! SOLUTION TRIDIAGONAL SYSTEM FOR THE WATERLEVELS ! if (nhystp==noneighbors) then ! ! Single domain case without domain decomposition ! The next piece of code in this IF-statement works for both serial and parallel runs ! In case of parallel runs twisted factorization technique is employed which is ! perfectly parallizable for two processors only. In case of more than 2 processors, ! this technique is combined with the block Jacobi approach at coupling points between ! pairs of "twisted" processors. Improvement in convergence is achieved by means of ! alternating the pairs of twisted processors at each iteration. ! if ( nproc > 2 ) then icol = mod(iter,2) else icol = 1 endif ! call timer_start(timer_sud_solve, gdp) if ( mod(inode,2) == icol ) then ! ! FORWARD SWEEP (elimination) ! ! Division by the pivot for nmf is not needed anymore ! because of row scaling ! do m = 2, mmaxddb nm = m*icx - icxy do n = 1, nmaxddb nm = nm + icy if (kcs(nm) > 0) then bi = 1.0_hp / (real(b(nm),hp) - real(a(nm),hp)*real(c(nm-icx),hp)) c(nm) = real(c(nm),hp) * bi d(nm) = (real(d(nm),hp) - real(a(nm),hp)*real(d(nm-icx),hp)) * bi endif enddo enddo else ! ! BACKWARD SWEEP (elimination) ! ! Division by the pivot for nmlu is not needed anymore ! because of row scaling ! do m = mmaxddb-1, 1, -1 nm = m*icx - icxy do n = 1, nmaxddb nm = nm + icy if (kcs(nm) > 0) then bi = 1.0_hp / (real(b(nm),hp) - real(c(nm),hp)*real(a(nm+icx),hp)) a(nm) = real(a(nm),hp) * bi d(nm) = (real(d(nm),hp) - real(c(nm),hp)*real(d(nm+icx),hp)) * bi endif enddo enddo endif ! ! exchange coefficients a, b, c and d with neighbours for parallel runs ! call dfexchg ( a, 1, 1, dfloat, nm_pos, gdp ) call dfexchg ( b, 1, 1, dfloat, nm_pos, gdp ) call dfexchg ( c, 1, 1, dfloat, nm_pos, gdp ) call dfexchg ( d, 1, 1, dfloat, nm_pos, gdp ) call dfsync(gdp) ! if ( mod(inode,2) == icol ) then ! ! FORWARD SWEEP in coupling points (elimination) ! do m = 1, mmaxddb nm = m*icx - icxy do n = 1, nmaxddb nm = nm + icy if (kcs(nm) == -1) then bi = 1.0_fp / (real(b(nm),hp) - real(a(nm),hp)*real(c(nm-icx),hp)) c (nm) = real(c(nm),hp) * bi d (nm) = (real(d(nm),hp) - real(a(nm),hp)*real(d(nm-icx),hp)) * bi s1(nm) = d(nm) endif enddo enddo ! ! BACKWARD SWEEP (substitution) ! nmlu = mmaxddb*icx - icxy do n = 1, nmaxddb nmlu = nmlu + icy if (kcs(nmlu) > 0) s1(nmlu) = d(nmlu) enddo do m = mmaxddb - 1, 1, -1 nm = m*icx - icxy do n = 1, nmaxddb nm = nm + icy if (kcs(nm) > 0) then d(nm) = d(nm) - c(nm)*d(nm + icx) s1(nm) = d(nm) endif enddo enddo else ! ! BACKWARD SWEEP in coupling points (elimination) ! do m = mmaxddb, 1, -1 nm = m*icx - icxy do n = 1, nmaxddb nm = nm + icy if (kcs(nm) == -1) then bi = 1.0_hp / (real(b(nm),hp) - real(c(nm),hp)*real(a(nm+icx),hp)) a (nm) = real(a(nm),hp) * bi d (nm) = (real(d(nm),hp) - real(c(nm),hp)*real(d(nm+icx),hp)) * bi s1(nm) = d(nm) endif enddo enddo ! ! FORWARD SWEEP (substitution) ! nmf = icx - icxy do n = 1, nmaxddb nmf = nmf + icy if (kcs(nmf) > 0) s1(nmf) = d(nmf) enddo do m = 2, mmaxddb nm = m*icx - icxy do n = 1, nmaxddb nm = nm + icy if (kcs(nm) > 0) then d (nm) = d(nm) - a(nm)*d(nm - icx) s1(nm) = d(nm) endif enddo enddo endif ! ! exchange s1 with neighbours for parallel runs ! call dfexchg ( s1, 1, 1, dfloat, nm_pos, gdp ) ! ! insert block Jacobi equation in coupling points ! do nm = 1, nmmax if ( kcs(nm) == -1 ) then a(nm) = 0.0 b(nm) = 1.0 c(nm) = 0.0 d(nm) = s1(nm) endif enddo call timer_stop(timer_sud_solve, gdp) else ! ! Domain decomposition: ! ! Wang solver for subdomains ! ! METHOD OF WANG ! ! First part of Wang's algoritm: pre-elimination: ! call timer_start(timer_sud_wangpre, gdp) call wangp1(s1 ,kcs ,irocol ,norow ,icx , & & icy ,j ,nmmaxj ,a ,b , & & c ,d ,gdp ) call timer_stop(timer_sud_wangpre, gdp) ! ! Now second part of Wang's algoritm: elimination of reduced ! system. This is carried out by a global mapper process. ! At end of global mapper process: ! at coupling point at the left side and ! at last computational point at the right side we have that ! a=c=0, b=1 and d=zeta=new water elevation ! ! Now Mapper builds and solves the reduced system of equation ! if (icx==1) then intdir = 0 else ! ! intdir = 1 corresponds to left_to_Right direction ! (see GAWS routines) ! intdir = 1 endif call timer_start(timer_sud_gwsslv, gdp) call gwsslv(intdir ) call timer_stop(timer_sud_gwsslv, gdp) ! ! Now third part of Wang's algoritm: back substitution ! call timer_start(timer_sud_wangback, gdp) call wangp3(s1 ,kcs ,irocol ,norow ,icx , & & icy ,j ,nmmaxj ,a ,b , & & c ,d ,gdp ) call timer_stop(timer_sud_wangback, gdp) ! ! in case of Hydra (DD method, Wang approach) array d does ! not contain the solution, but the right-hand side evaluation. ! Since array d (and also s1) is used in the remainder of SUD ! for the new water elevation, we have to copy s1 to d. ! do nm = 1, nmmax d(nm) = s1(nm) enddo endif ! ! TOTAL WATERDEPTH DRYING IN SUD ! call timer_start(timer_sud_rest, gdp) nmu = +icx do nm = 1, nmmax nmu = nmu + 1 if (kfu(nm)==1) then ! ! Special approach for 2D weir points: ! - depth value DPU is not corrected like in general (3D) weir case ! - crest height is explicitly taken into account in drying check ! hucres = 1E9 if (abs(kspu(nm, 0))==9) then if (umean(nm)>=0.001) then hucres = d(nm) + hkru(nm) elseif (umean(nm)<= - 0.001) then hucres = d(nm + icx) + hkru(nm) else hucres = max(d(nm + icx), d(nm)) + hkru(nm) endif endif hnm = tetau(nm)*d(nm) + (1. - tetau(nm))*d(nmu) + dpu(nm) ! ! CHECK FOR DRYING ! if (min(hnm, hucres)<=drytrsh) then aa(nm) = 0.0 bb(nm) = 1.0 cc(nm) = 0.0 dd(nm) = 0.0 kfu(nm) = 0 itr = 1 else if (momsol=='flood ') then bb(nm) = 1.0 else bb(nm) = hu(nm)/hnm endif hu(nm) = hnm endif endif enddo call timer_stop(timer_sud_rest, gdp) ! ! determine global maximum of 'itr' over all nodes ! Note: this enables to synchronize the repeating computation ! call dfreduce( itr, 1, dfint, dfmax, gdp ) ! ! REPEAT COMPUTATION IF POINT IS SET DRY ! FIRST RESET HU ! ! Domain decomposition: ! Synchronize on drying before finishing solve zeta ! ! Note that if iter 2 ) call dfmassc (s1 ,u1 ,qxk ,hu ,d0 , & & dpu ,porosu ,gsqs ,guu ,tetau , & & kcs ,kcu ,kfu ,thick ,nmmax , & & kmax ,icx ,gdp ) ! ! compute depth-averaged velocity ! do nm = 1, nmmax umean(nm) = 0.0 enddo do k = 1, kmax do nm = 1, nmmax umean(nm) = umean(nm) + thick(k)*u1(nm, k) enddo enddo ! ! Domain decomposition: ! nhystp = nxtstp(d3dflow_finish_wang, gdp) ! ! End of Domain decomposition addition ! if (kmax>1) then ! ! COMPUTATION VERTICAL VELOCITIES AND DISCHARGES ! ! Initialise arrays qzk and w1 for all (nm,k) ! qzk = 0.0 w1 = 0.0 ! do k = 1, kmax do nm = 1, nmmax if (kcs(nm)==1) then w1(nm, k) = w1(nm, k - 1) + thick(k)*s1(nm)*hdti & & + (qxk(nm, k) - qxk(nm - icx, k) & & - d0k(nm, k) ) / gsqs(nm) qzk(nm, k) = w1(nm, k)*gsqs(nm) endif enddo enddo ! ! exchange w1 with neighbours for parallel runs ! call dfexchg ( w1, 0, kmax, dfloat, nm_pos, gdp ) ! ! compute vertical discharge ! ! epsomb = max(eps, eps*hdt) ! error = .false. do nm = 1, nmmax if (abs(w1(nm, kmax))>epsomb) then error = .true. w1(nm, kmax) = 0.0 endif enddo ierror = 0 if (error) ierror = 1 call dfreduce( ierror, 1, dfint, dfmax, gdp ) error = ierror==1 if (error) then write (errtxt, '(a,e12.3,a,i0,a)') 'Mass closure error exceeds ', & & epsomb, ' after ', ntstep, ' timesteps.' call prterr(lundia, 'U190', trim(errtxt)) endif endif call timer_stop(timer_sud_veldisch, gdp) ! ! Optionally compute individual momentum terms for output ! if (gdp%gdflwpar%flwoutput%momentum) then ! pass hu or hu0 ? call timer_start(timer_sud_cucnp, gdp) call cucnp(dischy ,icreep ,dpdksi ,s0 ,u0 , & & v1 ,w1 ,hu ,hv ,dps ,dpu , & & umean ,guu ,gvv ,gvu ,gsqs , & & gvd ,gud ,gvz ,gsqiu ,qxk , & & qyk ,disch ,umdis ,mnksrc ,dismmt ,j , & & nmmaxj ,nmmax ,kmax ,icx ,icy , & & nsrc ,lsecfl ,lstsci ,betac ,aak , & & bbk ,cck ,ddk ,bbka ,bbkc , & & thick ,sig ,rho ,sumrho ,vicuv , & & vnu2d ,vicww ,wsu ,fxw ,wsbodyu , & & rxx ,rxy ,kcs ,kcu ,kfu ,kfv , & & kfs ,kspu ,kadu ,kadv ,dfu ,deltau , & & tp ,rlabda ,cfurou ,cfvrou ,rttfu , & & r0 ,diapl ,rnpl ,taubpu ,taubsu , & & windsu ,patm ,fcorio ,ubrlsu ,uwtypu , & & hkru ,pship ,tgfsep ,dteu ,ua , & & ub ,ustokes ,.true. ,u1 ,s1 ,gdp ) call timer_stop(timer_sud_cucnp, gdp) endif end subroutine sud