|
XBeach
|
00001 !============================================================================== 00002 ! MODULE NH_MAIN 00003 !============================================================================== 00004 00005 ! DATE AUTHOR CHANGES 00006 ! 00007 ! october 2009 Pieter Bart Smit New module 00008 ! october 2014 Pieter Bart Smit Included an updated pressure correction scheme/general maintenance 00009 00010 ! -- REMARKS -- 00011 ! There are "two" mutually exclusive implementations of the non-hydrostatic model 00012 ! in this file. The first is the original non-hydrostatic code (see Draft report) 00013 ! which has been left untouched (except for fixing some consistency errors). The 00014 ! second implementation is the reduced two-layer formulation. The new implementation 00015 ! uses two-layers in the vertical (controlled by a thickness parameter) and should 00016 ! give better results for shorter wave components. Moreover, it contains the previous 00017 ! code as a special case (thickness of the lower layer is set to 0). 00018 ! 00019 ! INCLUDED SUBROUTINES 00020 ! + nonh_cor 00021 ! This is the entry routine for the nonhydrostatic code when called from 00022 ! flowtimestep. Depending on the parameter nh_red it calls the old nh-routines 00023 ! or the new reduced routines. 00024 ! 00025 ! + nonh_init 00026 ! Initializes the various resources required. 00027 ! 00028 ! + nonh_1lay_cor 00029 ! The "old" non-hydrostatic correction routine, calculates the new dynamic 00030 ! pressure. 00031 ! 00032 ! + nonh_1lay_pred 00033 ! 00034 ! 00035 module nonh_module 00036 00037 implicit none 00038 save 00039 private 00040 00041 00042 ! If mpi is defined, the non-hydrostatic module is NOT included in the compilation 00043 ! to avoid unwanted side effects. 00044 00045 !****************************************************************************** 00046 ! INTERFACE 00047 !****************************************************************************** 00048 00049 00050 !----------------------------- PARAMETERS ----------------------------------- 00051 00052 include 'nh_pars.inc' !Default precision etc. 00053 00054 !----------------------------- VARIABLES ----------------------------------- 00055 00056 logical :: initialized = .false. !Indicates whether or not the nh_module has been initialized 00057 00058 ! 00059 ! -- Numerical parameter -- 00060 ! 00061 real (kind=rKind), allocatable, dimension(:,:) :: wcoef !The relative thickness of the lower layer in the reduced 00062 !two-layer model at z-points (optimisation parameter in "old model") 00063 ! 00064 ! -- Additional physical variables -- 00065 ! 00066 real (kind=rKind), allocatable, dimension(:,:) :: dp !Change in pressure between the old and 00067 !new timestep 00068 real (kind=rKind), allocatable, dimension(:,:) :: Wm !The mean vertical velocity in the top 00069 !layer (nonhq3d) 00070 real (kind=rKind), allocatable, dimension(:,:) :: Wm_old !The mean vertical velocity in the top 00071 !layer on the previous timestep (nonhq3d) 00072 real (kind=rKind), allocatable, dimension(:,:) :: Wm0 !The mean vertical velocity in the bottom 00073 !layer (nonhq3d) 00074 real (kind=rKind), allocatable, dimension(:,:) :: dU0 !The velocity (s-dir) difference between the top 00075 !and bottom layers at the previous timestep 00076 real (kind=rKind), allocatable, dimension(:,:) :: dV0 !The velocity (n-dir) difference between the top 00077 !and bottom layers at the previous timestep 00078 real (kind=rKind), allocatable, dimension(:,:) :: omega !The velocity vertical transport velocity 00079 00080 ! 00081 ! -- Arrays used for building the Poisson equation -- 00082 ! 00083 real (kind=rKind), allocatable, dimension(:,:,:) :: au !Pressure coeficients for u(i,j), au(0,i,j) refers to 00084 !s%nhpres(i,j) and au(1,i,j) to s%nhpres(i+1,j) 00085 real (kind=rKind), allocatable, dimension(:,:,:) :: av !Pressure coeficients for v(i,j), av(0,i,j) refers to 00086 !s%nhpres(i,j) and av(1,i,j) to s%nhpres(i,j+1) 00087 real (kind=rKind), allocatable, dimension(:,:,:) :: adu !Pressure coeficients for du(i,j), adu(0,i,j) refers to 00088 !s%nhpres(i,j) and au(1,i,j) to s%nhpres(i+1,j) 00089 real (kind=rKind), allocatable, dimension(:,:,:) :: adv !Pressure coeficients for dv(i,j), adv(0,i,j) refers to 00090 !s%nhpres(i,j) and adv(1,i,j) to s%nhpres(i,j+1) 00091 real (kind=rKind), allocatable, dimension(:,:,:) :: mat !5-diagonal Pressure matrix. 00092 ! Diagonal Pressure Point 00093 ! 1 i ,j 00094 ! 2 i-1,j 00095 ! 3 i+1,j 00096 ! 4 i ,j-1 00097 ! 5 i ,j+1 00098 real (kind=rKind), allocatable, dimension(:,:) :: rhs !RHS of the pressure matrix 00099 00100 real (kind=rKind), allocatable, dimension(:,:,:) :: aws !Pressure coeficients for ws(k,i,j) !**OBSOLETE IN nonhq3d** 00101 real (kind=rKind), allocatable, dimension(:,:,:) :: awb !Pressure coeficients for wb(k,i,j) !**OBSOLETE IN nonhq3d** 00102 real (kind=rKind), allocatable, dimension(:,:) :: aur !RHS for du(i,j) !**OBSOLETE IN nonhq3d** 00103 real (kind=rKind), allocatable, dimension(:,:) :: avr !RHS for dv(i,j) !**OBSOLETE IN nonhq3d** 00104 real (kind=rKind), allocatable, dimension(:,:) :: awbr !Pressure coeficients for wb(k,i,j) !**OBSOLETE IN nonhq3d** 00105 real (kind=rKind), allocatable, dimension(:,:) :: awsr !Pressure coeficients for ws(k,i,j) !**OBSOLETE IN nonhq3d** 00106 00107 ! 00108 ! -- Additional grid parameters -- 00109 ! 00110 real (kind=rKind), allocatable, dimension(:,:) :: zsu !free surface-level in u-velocity points 00111 real (kind=rKind), allocatable, dimension(:,:) :: zsv !free surface-level in v-velocity points 00112 real (kind=rKind), allocatable, dimension(:,:) :: zbu !free bed-level in u-velocity points 00113 real (kind=rKind), allocatable, dimension(:,:) :: zbv !free bed-level in v-velocity points 00114 00115 real (kind=rKind), allocatable, dimension(:) :: dxz !**OBSOLETE IN nonhq3d** 00116 real (kind=rKind), allocatable, dimension(:) :: dyz !**OBSOLETE IN nonhq3d** 00117 real (kind=rKind), allocatable, dimension(:) :: dxu !**OBSOLETE IN nonhq3d** 00118 real (kind=rKind), allocatable, dimension(:) :: dyv !**OBSOLETE IN nonhq3d** 00119 real (kind=rKind), allocatable, dimension(:) :: ddxz !**OBSOLETE IN nonhq3d** 00120 real (kind=rKind), allocatable, dimension(:) :: ddyz !**OBSOLETE IN nonhq3d** 00121 real (kind=rKind), allocatable, dimension(:) :: ddxu !**OBSOLETE IN nonhq3d** 00122 real (kind=rKind), allocatable, dimension(:) :: ddyv !**OBSOLETE IN nonhq3d** 00123 00124 integer(kind=iKind),allocatable,dimension(:,:) :: nonhU 00125 integer(kind=iKind),allocatable,dimension(:,:) :: nonhV 00126 integer(kind=iKind),allocatable,dimension(:,:) :: nonhZ 00127 real (kind=rKind), allocatable, dimension(:,:) :: lbreakcond 00128 00129 !--- PUBLIC SUBROUTINES --- 00130 public nonh_init 00131 public nonh_cor 00132 public nonh_init_wcoef 00133 00134 !--- PRIVATE SUBROUTINES 00135 00136 ! NONE 00137 contains 00138 ! 00139 !****************************************************************************** 00140 ! SUBROUTINES/FUNCTIONS 00141 !****************************************************************************** 00142 ! 00143 00144 ! 00145 !============================================================================== 00146 subroutine nonh_cor(s,par,ipredcor,uu0,vv0) 00147 !============================================================================== 00148 ! 00149 ! 00150 ! DATE AUTHOR CHANGES 00151 ! 00152 ! October 2014 Pieter Bart Smit new subroutine 00153 ! 00154 !------------------------------------------------------------------------------- 00155 ! DECLARATIONS 00156 !------------------------------------------------------------------------------- 00157 00158 !-------------------------- PURPOSE ---------------------------- 00159 00160 ! This is the new entry routine for the non-hydrostatic module. It is called twice during each timestep by the 00161 ! flow module, (I) (withe path=0) before the second-order corrections are applied (if applicable) to include the known pressure 00162 ! explicitly in the momentum equations, and (II) (path=1) just after the second order advection correction to calculate 00163 ! the new pressures implicitly. The codepath taken is specified with the path parameter (path=0 for step I, and path=1 for 00164 ! step II). 00165 ! 00166 ! Moreover, depending on the parameter nonhq3d, which activates the new reduced 2 layer model, again different codepaths 00167 ! are taken: 00168 ! 00169 ! IF nonhq3d=0 00170 ! call the "old" versions of the non-hydrostatic routines 00171 ! (path=0) nonh_1lay_pred (renamed from nonh_explicit) 00172 ! (path=1) nonh_1lay_cor (renamed from nonh_cor) 00173 ! IF nonhq3d=1 00174 ! call the "old" versions of the non-hydrostatic routine 00175 ! IF 1D 00176 ! (path=0) nonh_2lay_pred_2DV (renamed from nonh_explicit) 00177 ! (path=1) nonh_2lay_cor_2DV (renamed from nonh_cor) 00178 ! ELSE 00179 ! (path=0) nonh_2lay_pred_3D (renamed from nonh_explicit) 00180 ! (path=1) nonh_2lay_cor_3D (renamed from nonh_cor) 00181 ! 00182 !-------------------------- DEPENDENCIES ---------------------------- 00183 00184 use spaceparams 00185 use params, only: parameters 00186 use xmpi_module 00187 00188 !-------------------------- ARGUMENTS ---------------------------- 00189 00190 type(spacepars) ,intent(inout) :: s 00191 type(parameters),intent(in) :: par 00192 integer ,intent(in) :: ipredcor !Select if we want to apply the explicit predicition / or the implicit correction 00193 00194 real(kind=rkind),dimension(s%nx+1,s%ny+1),intent(in) :: uu0 00195 real(kind=rkind),dimension(s%nx+1,s%ny+1),intent(in) :: vv0 00196 00197 !------------------------------------------------------------------------------- 00198 ! IMPLEMENTATION 00199 !------------------------------------------------------------------------------- 00200 00201 select case ( par%nonhq3d ) 00202 ! 00203 case (1) 00204 ! 00205 ! The "new" reduced two-layer formulations are used 00206 ! 00207 if ( ipredcor == 0 ) then 00208 ! -- Predictor 00209 if ( s%ny == 0 ) then 00210 ! -- 2dv code path 00211 call nonh_2lay_pred_2dV(s,par,uu0) 00212 ! 00213 else 00214 ! -- "3d" code path 00215 call nonh_2lay_pred_3d(s,par,uu0,vv0) 00216 ! 00217 endif 00218 ! 00219 else 00220 ! 00221 if ( s%ny == 0 ) then 00222 ! 00223 call nonh_2lay_cor_2dV(s,par) 00224 ! 00225 else 00226 ! 00227 call nonh_2lay_cor_3d(s,par) 00228 ! 00229 endif 00230 ! 00231 endif 00232 ! 00233 case default 00234 ! 00235 ! The "old" code is used 00236 ! 00237 if ( ipredcor == 0 ) then 00238 ! -- Predictor 00239 call nonh_1lay_pred(s,par) 00240 ! 00241 else 00242 ! -- Corrector 00243 call nonh_1lay_cor(s,par) 00244 ! 00245 endif 00246 ! 00247 end select 00248 ! 00249 end subroutine nonh_cor 00250 00251 00252 ! 00253 !============================================================================== 00254 subroutine nonh_init(s,par) 00255 !============================================================================== 00256 ! 00257 00258 ! DATE AUTHOR CHANGES 00259 ! 00260 ! October 2010 Pieter Bart Smit New Subroutine 00261 00262 !------------------------------------------------------------------------------- 00263 ! DECLARATIONS 00264 !------------------------------------------------------------------------------- 00265 00266 ! 00267 !-------------------------- PURPOSE ---------------------------- 00268 ! 00269 ! Initializes nh subroutines 00270 ! 00271 !-------------------------- DEPENDENCIES ---------------------------- 00272 00273 use spaceparams 00274 use params, only: parameters 00275 00276 !-------------------------- ARGUMENTS ---------------------------- 00277 ! 00278 type(spacepars) ,intent(inout) :: s 00279 type(parameters),intent(in) :: par 00280 ! 00281 00282 !------------------------------------------------------------------------------- 00283 ! IMPLEMENTATION 00284 !------------------------------------------------------------------------------- 00285 00286 ! open(unit=PrintFileUnit,file=trim(PrintFileName)) 00287 00288 if(.not. xcompute) return 00289 00290 allocate(au (0:1,s%nx+1,s%ny+1)); au = 0.0_rKind 00291 allocate(av (0:1,s%nx+1,s%ny+1)); av = 0.0_rKind 00292 allocate(mat(5,s%nx+1,s%ny+1)) ; mat = 0.0_rKind 00293 allocate(rhs( s%nx+1,s%ny+1)) ; rhs = 0.0_rKind 00294 allocate(zbu( s%nx+1,s%ny+1)) ; zbu = 0.0_rKind 00295 allocate(zbv( s%nx+1,s%ny+1)) ; zbv = 0.0_rKind 00296 allocate(zsu( s%nx+1,s%ny+1)) ; zsu = 0.0_rKind 00297 allocate(zsv( s%nx+1,s%ny+1)) ; zsv = 0.0_rKind 00298 allocate(nonhU( s%nx+1,s%ny+1)); nonhU = 1 00299 allocate(nonhV( s%nx+1,s%ny+1)); nonhV = 1 00300 allocate(nonhZ( s%nx+1,s%ny+1)); nonhZ = 1 00301 allocate(lbreakcond(s%nx+1,s%ny+1)); lbreakcond = par%maxbrsteep 00302 allocate(dp(s%nx+1,s%ny+1)) ; dp = 0.0_rKind 00303 allocate(Wm(s%nx+1,s%ny+1)) ;Wm = 0.0_rKind 00304 allocate(Wm_old(s%nx+1,s%ny+1)) ;Wm_old = 0.0_rKind 00305 allocate(Wcoef(s%nx+1,s%ny+1)) ;Wcoef = 1.0_rKind 00306 allocate(aws (5,s%nx+1,s%ny+1)) ; aws = 0.0_rKind 00307 ! 00308 ! The reduced two-layer model 00309 ! 00310 if ( par%nonhq3d == 1 ) then 00311 ! 00312 allocate(omega(s%nx+1,s%ny+1)); omega = 0.0_rKind 00313 allocate(adu (0:1,s%nx+1,s%ny+1)); adu = 0.0_rKind 00314 allocate(adv (0:1,s%nx+1,s%ny+1)); adv = 0.0_rKind 00315 allocate(Wm0(s%nx+1,s%ny+1)) ;Wm0 = 0.0_rKind 00316 allocate(dU0(s%nx+1,s%ny+1)) ;dU0 = 0.0_rKind 00317 allocate(dV0(s%nx+1,s%ny+1)) ;dV0 = 0.0_rKind 00318 ! 00319 else 00320 ! 00321 ! These grid variables are obsolete in the new reduced model 00322 ! 00323 allocate(dxz ( s%nx+1)) ; dxz = 0.0_rKind 00324 allocate(dyz ( s%ny+1)) ; dyz = 0.0_rKind 00325 allocate(dxu ( s%nx+1)) ; dxu = 0.0_rKind 00326 allocate(dyv ( s%ny+1)) ; dyv = 0.0_rKind 00327 00328 allocate(ddxz ( s%nx+1)) ; ddxz = 0.0_rKind 00329 allocate(ddyz ( s%ny+1)) ; ddyz = 0.0_rKind 00330 allocate(ddxu ( s%nx+1)) ; ddxu = 0.0_rKind 00331 allocate(ddyv ( s%ny+1)) ; ddyv = 0.0_rKind 00332 allocate(awb (5,s%nx+1,s%ny+1)) ; awb = 0.0_rKind 00333 allocate(aur( s%nx+1,s%ny+1)) ; aur = 0.0_rKind 00334 allocate(avr( s%nx+1,s%ny+1)) ; avr = 0.0_rKind 00335 allocate(awbr(s%nx+1,s%ny+1)) ; awbr = 0.0_rKind 00336 allocate(awsr(s%nx+1,s%ny+1)) ; awsr = 0.0_rKind 00337 00338 !Initialize grid variables 00339 dyz = s%dnz(1,:) 00340 dyv = s%dnv(1,:) 00341 ddyz = 1.0_rKind/dyz 00342 ddyv = 1.0_rKind/dyv 00343 00344 dxu = s%dsu(:,1) 00345 dxz = s%dsz(:,1) 00346 ddxu = 1.0_rKind/dxu 00347 ddxz = 1.0_rKind/dxz 00348 ! 00349 endif 00350 00351 mat(1,1,:) = 1.0_rKind 00352 mat(1,s%nx+1,:) = 1.0_rKind 00353 mat(1,:,1) = 1.0_rKind 00354 mat(1,:,s%ny+1) = 1.0_rKind 00355 00356 initialized = .true. 00357 00358 if ( par%nonhq3d == 1 ) then 00359 ! 00360 wcoef = par%nhlay 00361 ! 00362 else 00363 ! 00364 call nonh_init_wcoef(s,par) 00365 ! 00366 endif 00367 ! 00368 00369 !Initialize levels at u/v points (zu,zbu etc.) 00370 call zuzv(s) 00371 00372 end subroutine nonh_init 00373 00374 00375 00376 00377 !============================================================================== 00378 ! "ORIGINAL" NON-HYDROSTATIC ROUTINES 00379 !============================================================================== 00380 ! 00381 ! 00382 ! 00383 !============================================================================== 00384 subroutine nonh_1lay_cor(s,par) 00385 !============================================================================== 00386 ! 00387 00388 ! DATE AUTHOR CHANGES 00389 ! 00390 ! October 2010 Pieter Bart Smit New Subroutine 00391 ! November 2010 Pieter Bart Smit Added explicit prediction for 2nd order scheme 00392 00393 !------------------------------------------------------------------------------- 00394 ! DECLARATIONS 00395 !------------------------------------------------------------------------------- 00396 00397 !-------------------------- PURPOSE ---------------------------- 00398 00399 ! Releases resources 00400 00401 !-------------------------- DEPENDENCIES ---------------------------- 00402 use spaceparams 00403 use params, only: parameters 00404 use solver_module, only: solver_solvemat 00405 use xmpi_module 00406 !-------------------------- ARGUMENTS ---------------------------- 00407 00408 type(spacepars) ,intent(inout) :: s 00409 type(parameters),intent(in) :: par 00410 00411 !-------------------------- LOCAL VARIABLES ---------------------------- 00412 00413 !Indices 00414 integer(kind=iKind) :: i !Index variables 00415 integer(kind=iKind) :: j !Index variables 00416 integer(kind=iKind) :: jmin,jmax !Index variables for superfast1D 00417 00418 ! logical 00419 logical :: sf1d !Switch for superfast1D 00420 00421 !Work 00422 00423 real(kind=rKind) :: dzs_e,dzs_w 00424 real(kind=rKind) :: dzs_s,dzs_n 00425 real(kind=rKind) :: dzb_e,dzb_w 00426 real(kind=rKind) :: dzb_n,dzb_s 00427 00428 00429 ! 00430 !------------------------------------------------------------------------------- 00431 ! IMPLEMENTATION 00432 !------------------------------------------------------------------------------- 00433 ! 00434 dzs_e = 0 00435 dzs_w = 0 00436 dzs_s = 0 00437 dzs_n = 0 00438 dzb_e = 0 00439 dzb_w = 0 00440 dzb_n = 0 00441 dzb_s = 0 00442 if (s%ny>0) then 00443 sf1d = .false. 00444 jmin = 2 00445 jmax = s%ny 00446 else 00447 sf1d = .true. 00448 jmin = 1 00449 jmax = 1 00450 endif 00451 00452 !call timer_start(timer_flow_nonh) 00453 if (.not. initialized) then 00454 call nonh_init(s,par) 00455 endif 00456 00457 call zuzv(s) 00458 00459 !Built pressure coefficients U 00460 aur = s%uu 00461 avr = s%vv 00462 00463 00464 aur(1,:) = s%uu(1,:) 00465 aur(s%nx,:) = s%uu(s%nx,:) 00466 00467 avr(:,1) = s%vv(:,1) 00468 if (s%ny>0) avr(:,s%ny) = s%vv(:,s%ny) 00469 00470 !Built pressure coefficients for W 00471 00472 !AW Bottom 00473 if (sf1d) then 00474 do i=2,s%nx 00475 if (nonhZ(i,1) == 1) then 00476 if (nonhU(i,1)*nonhU(i-1,1) == 1) then 00477 dzb_e = .5_rKind*(zbu(i,1)-zbu(i-1,1))*ddxz(i) 00478 dzb_w = dzb_e 00479 elseif (nonhU(i ,1) == 0) then 00480 dzb_e = .0_rKind 00481 dzb_w = .5_rKind*(s%zb(i,1)-zbu(i-1,1))*ddxz(i) 00482 elseif (nonhU(i-1,1) == 0) then 00483 dzb_e = .5_rKind*(zbu(i,1)-s%zb(i,1)) *ddxz(i) 00484 dzb_w = .0_rKind 00485 endif 00486 00487 if (nonhV(i,1) == 1) then 00488 dzb_s = .0_rKind 00489 dzb_n = dzb_s 00490 elseif (nonhV(i ,1 ) == 0) then 00491 dzb_s = .0_rKind 00492 dzb_n = .5_rKind*(s%zb(i,1)-zbv(i,1))*ddyz(1) 00493 endif 00494 00495 awb(1,i,1) = dzb_e * au(0,i ,1 )+dzb_w*au(1,i-1,1 ) & !main diagonal 00496 + dzb_s * av(0,i ,1 )+dzb_n*av(1,i ,1) 00497 awb(2,i,1) = dzb_w * au(0,i-1,1 ) !west 00498 awb(3,i,1) = dzb_e * au(1,i ,1 ) !east 00499 awb(4,i,1) = dzb_n * av(0,i ,1 ) !south 00500 awb(5,i,1) = dzb_s * av(1,i ,1 ) !north 00501 00502 awbr(i,1) = dzb_e*s%uu(i,1)+dzb_w*s%uu(i-1,1 ) & 00503 + dzb_s*s%vv(i,1)+dzb_n*s%vv(i ,1 ) 00504 else 00505 awb(:,i,1) = 0.0_rKind 00506 awbr(i,1) = 0.0_rKind 00507 endif 00508 enddo 00509 else 00510 do j=2,s%ny 00511 do i=2,s%nx 00512 if (nonhZ(i,j) == 1) then 00513 if (nonhU(i,j)*nonhU(i-1,j) == 1) then 00514 dzb_e = .5_rKind*(zbu(i,j)-zbu(i-1,j))*ddxz(i) 00515 ! FB: this was dzs_e, which was not defined. Assumed it was a typo. TODO: check. 00516 dzb_w = dzb_e 00517 elseif (nonhU(i ,j) == 0) then 00518 dzb_e = .0_rKind 00519 dzb_w = .5_rKind*(s%zb(i,j)-zbu(i-1,j))*ddxz(i) 00520 elseif (nonhU(i-1,j) == 0) then 00521 dzb_e = .5_rKind*(zbu(i,j)-s%zb(i,j)) *ddxz(i) 00522 dzb_w = .0_rKind 00523 endif 00524 00525 if (nonhV(i,j)*nonhV(i,j-1) == 1) then 00526 dzb_s = .5_rKind*(zbv(i,j)-zbv(i,j-1))*ddyz(j) 00527 dzb_n = dzb_s 00528 elseif (nonhV(i ,j ) == 0) then 00529 dzb_s = .0_rKind 00530 dzb_n = .5_rKind*(s%zb(i,j)-zbv(i,j-1))*ddyz(j) 00531 elseif (nonhV(i ,j-1) == 0) then 00532 dzb_s = .5_rKind*(zbv(i,j)-s%zb(i,j)) *ddyz(j) 00533 dzb_n = .0_rKind 00534 endif 00535 00536 awb(1,i,j) = dzb_e * au(0,i ,j )+dzb_w*au(1,i-1,j ) & !main diagonal 00537 + dzb_s * av(0,i ,j )+dzb_n*av(1,i ,j-1) 00538 awb(2,i,j) = dzb_w * au(0,i-1,j ) !west 00539 awb(3,i,j) = dzb_e * au(1,i ,j ) !east 00540 awb(4,i,j) = dzb_n * av(0,i ,j-1) !south 00541 awb(5,i,j) = dzb_s * av(1,i ,j ) !north 00542 00543 awbr(i,j) = dzb_e*s%uu(i,j)+dzb_w*s%uu(i-1,j ) & 00544 + dzb_s*s%vv(i,j)+dzb_n*s%vv(i ,j-1) 00545 else 00546 awb(:,i,j) = 0.0_rKind 00547 awbr(i,j) = 0.0_rKind 00548 endif 00549 enddo 00550 enddo 00551 endif 00552 00553 do j=jmin,jmax 00554 do i=2,s%nx 00555 if (nonhZ(i,j) == 1) then 00556 aws(1,i,j) = wcoef(i,j)*par%dt*2.0_rKind/s%hh(i,j)-awb(1,i,j) 00557 aws(2:5,i,j) = -awb(2:5,i,j) 00558 awsr(i,j) = 2.0_rKind * Wm(i,j)-awbr(i,j) 00559 else 00560 aws(:,i,j) = 0.0_rKind 00561 awsr(i,j) = 0.0_rKind 00562 endif 00563 enddo 00564 enddo 00565 00566 !Substitute in the continuity equation 00567 if (sf1d) then 00568 do i=2,s%nx 00569 if (nonhZ(i,1)==1) then 00570 if (nonhU(i,1)*nonhU(i-1,1) == 1) then 00571 dzs_e = .5_rKind*(zsu(i,1)-zsu(i-1,1))*dyz(1) 00572 dzs_w = dzs_e 00573 elseif (nonhU(i ,1) == 0) then 00574 dzs_e = .0_rKind 00575 dzs_w = .5_rKind*(s%zs(i,1)-zsu(i-1,1))*dyz(1) 00576 elseif (nonhU(i-1,1) == 0) then 00577 dzs_e = .5_rKind*(zsu(i,1)-s%zs(i,1)) *dyz(1) 00578 dzs_w = .0_rKind 00579 endif 00580 00581 ! wwvv problem the following construction make it impossible 00582 ! wwvv to determine if dzs_s gets a value 00583 ! wwvv (What other values than 0 and 1 can nonhV have) 00584 ! wwvv I put this kind of varibles to zero at the start 00585 ! wwvv of the subroutine 00586 if (nonhV(i,1) == 1) then 00587 dzs_s = .0_rKind 00588 dzs_n = dzs_s 00589 elseif (nonhV(i ,1 ) == 0) then 00590 dzs_s = .0_rKind 00591 dzs_n = .5_rKind*(s%zs(i,1)-zsv(i,1))*dxz(i) 00592 endif 00593 00594 00595 mat(1,i,1) = dyz(1)*( s%hu(i,1)*au(0,i,1) - s%hu(i-1,1 )*au(1,i-1, 1) ) & !subs U left/right face 00596 + dxz(i)*( s%hv(i,1)*av(0,i,1) - s%hv(i ,1 )*av(1,i , 1) ) & !subs V left/rigth face 00597 - dzs_e*au(0,i,1) - dzs_w*au(1,i-1,1) & !kin. boun. top 00598 - dzs_s*av(0,i,1) - dzs_n*av(1,i-1,1) & !kin. boun. top 00599 + dxz(i)*dyz(1)*aws(1,i,1) 00600 00601 mat(2,i,1) = - dyz(1)*s%hu(i-1,1)*au(0,i-1, 1) & 00602 - dzs_w*au(0,i-1,1) + dxz(i)*dyz(1)*aws(2,i,1) 00603 00604 mat(3,i,1) = dyz(1)*s%hu(i,1) * au(1,i, 1) & 00605 - dzs_e*au(1,i ,1) + dxz(i)*dyz(1)*aws(3,i,1) 00606 00607 mat(4,i,1) = - dxz(i)*s%hv(i,1) *av(0,i,1) & 00608 - dzs_n*av(0,i,1) + dxz(i)*dyz(1)*aws(4,i,1) 00609 00610 mat(5,i,1) = dxz(i)*s%hv(i,1) * av(1,i,1) & 00611 - dzs_s*av(1,i,1) + dxz(i)*dyz(1)*aws(5,i,1) 00612 00613 rhs(i,1) = - dyz(1)*( s%hu(i,1)*aur(i,1) - s%hu(i-1,1 )*aur(i-1, 1) ) & !subs U left/right face 00614 - dxz(i)*( s%hv(i,1)*avr(i,1) - s%hv(i ,1 )*avr(i , 1) ) & !subs V left/rigth face 00615 + dzs_e*aur(i,1) + dzs_w*aur(i-1,1) & !kin. boun. top 00616 + dzs_s*avr(i,1) + dzs_n*avr(i, 1) & !kin. boun. top 00617 - dxz(i)*dyz(1)*awsr(i,1) 00618 else 00619 mat(1 ,i,1) = 1.0_rKind 00620 mat(2:5,i,1) = 0.0_rKind 00621 rhs(i,1) = 0.0_rKind 00622 endif 00623 enddo 00624 else 00625 do j=2,s%ny 00626 do i=2,s%nx 00627 if (nonhZ(i,j)==1) then 00628 if (nonhU(i,j)*nonhU(i-1,j) == 1) then 00629 dzs_e = .5_rKind*(zsu(i,j)-zsu(i-1,j))*dyz(j) 00630 dzs_w = dzs_e 00631 elseif (nonhU(i ,j) == 0) then 00632 dzs_e = .0_rKind 00633 dzs_w = .5_rKind*(s%zs(i,j)-zsu(i-1,j))*dyz(j) 00634 elseif (nonhU(i-1,j) == 0) then 00635 dzs_e = .5_rKind*(zsu(i,j)-s%zs(i,j)) *dyz(j) 00636 dzs_w = .0_rKind 00637 endif 00638 00639 if (nonhV(i,j)*nonhV(i,j-1) == 1) then 00640 dzs_s = .5_rKind*(zsv(i,j)-zsv(i,j-1))*dxz(i) 00641 dzs_n = dzs_s 00642 elseif (nonhV(i ,j ) == 0) then 00643 dzs_s = .0_rKind 00644 dzs_n = .5_rKind*(s%zs(i,j)-zsv(i,j-1))*dxz(i) 00645 elseif (nonhV(i ,j-1) == 0) then 00646 dzs_s = .5_rKind*(zsv(i,j)-s%zs(i,j)) *dxz(i) 00647 dzs_n = .0_rKind 00648 endif 00649 00650 00651 mat(1,i,j) = dyz(j)*( s%hu(i,j)*au(0,i,j) - s%hu(i-1,j )*au(1,i-1, j) ) & !subs U left/right face 00652 + dxz(i)*( s%hv(i,j)*av(0,i,j) - s%hv(i ,j-1)*av(1,i ,j-1) ) & !subs V left/rigth face 00653 - dzs_e*au(0,i,j) - dzs_w*au(1,i-1,j) & !kin. boun. top 00654 - dzs_s*av(0,i,j) - dzs_n*av(1,i,j-1) & !kin. boun. top 00655 + dxz(i)*dyz(j)*aws(1,i,j) 00656 00657 mat(2,i,j) = - dyz(j)*s%hu(i-1,j)*au(0,i-1, j) & 00658 - dzs_w*au(0,i-1,j) + dxz(i)*dyz(j)*aws(2,i,j) 00659 00660 mat(3,i,j) = dyz(j)*s%hu(i,j) * au(1,i, j) & 00661 - dzs_e*au(1,i ,j) + dxz(i)*dyz(j)*aws(3,i,j) 00662 00663 mat(4,i,j) = - dxz(i)*s%hv(i,j-1) *av(0,i,j-1) & 00664 - dzs_n*av(0,i,j-1) + dxz(i)*dyz(j)*aws(4,i,j) 00665 00666 mat(5,i,j) = dxz(i)*s%hv(i,j) * av(1,i,j) & 00667 - dzs_s*av(1,i,j) + dxz(i)*dyz(j)*aws(5,i,j) 00668 00669 rhs(i,j) = - dyz(j)*( s%hu(i,j)*aur(i,j) - s%hu(i-1,j )*aur(i-1, j) ) & !subs U left/right face 00670 - dxz(i)*( s%hv(i,j)*avr(i,j) - s%hv(i ,j-1)*avr(i ,j-1) ) & !subs V left/rigth face 00671 + dzs_e*aur(i,j) + dzs_w*aur(i-1,j) & !kin. boun. top 00672 + dzs_s*avr(i,j) + dzs_n*avr(i,j-1) & !kin. boun. top 00673 - dxz(i)*dyz(j)*awsr(i,j) 00674 else 00675 mat(1 ,i,j) = 1.0_rKind 00676 mat(2:5,i,j) = 0.0_rKind 00677 rhs(i,j) = 0.0_rKind 00678 endif 00679 enddo 00680 enddo 00681 endif 00682 00683 !Solve matrix 00684 !call timer_start(timer_flow_nonh_solv) 00685 dp = 0.0_rKind 00686 call solver_solvemat( mat , rhs , dp , s%nx, s%ny,par) 00687 !call timer_stop(timer_flow_nonh_solv) 00688 00689 00690 s%pres = s%pres + dp 00691 00692 !Correct u/v/w 00693 00694 !U 00695 do j=jmin_uu,jmax_uu 00696 do i=imin_uu,imax_uu 00697 s%uu(i,j) = aur(i,j) + au(1,i,j)*dp(i+1,j)+au(0,i,j)*dp(i,j) 00698 enddo 00699 enddo 00700 00701 !v 00702 if (sf1d) then 00703 do i=imin_vv,imax_vv 00704 s%vv(i,1) = avr(i,1) + av(1,i,1)*dp(i,1)+av(0,i,1)*dp(i,1) 00705 enddo 00706 else 00707 do j=jmin_vv,jmax_vv 00708 do i=imin_vv,imax_vv 00709 s%vv(i,j) = avr(i,j) + av(1,i,j)*dp(i,j+1)+av(0,i,j)*dp(i,j) 00710 enddo 00711 enddo 00712 endif 00713 00714 !W 00715 if (sf1d) then 00716 do i=imin_zs,imax_zs 00717 if (nonhZ(i,1) == 1) then 00718 s%ws(i,1) = awsr(i,1) + dp(i , 1) * aws(1,i,1) & 00719 + dp(i-1,1) * aws(2,i,1) + dp(i+1,1 ) * aws(3,i,1) & 00720 + dp(i ,1) * aws(4,i,1) + dp(i , 1) * aws(5,i,1) 00721 00722 else 00723 !If the point is excluded in nh then get ws from the kinematic boundary condition 00724 if (s%wetU(i,1)*s%wetU(i-1,1) == 1) then 00725 dzs_e = .5_rKind*(zsu(i,1)-zsu(i-1,1))*ddxz(i) 00726 dzs_w = dzs_e 00727 elseif (s%wetU(i ,1) == 0) then 00728 dzs_e = .0_rKind 00729 dzs_w = .5_rKind*(s%zs(i,1)-zsu(i-1,1))*ddxz(i) 00730 elseif (s%wetU(i-1,1) == 0) then 00731 dzs_e = .5_rKind*(zsu(i,1)-s%zs(i,1)) *ddxz(i) 00732 dzs_w = .0_rKind 00733 else 00734 dzs_e = .0_rKind 00735 dzs_w = .0_rKind 00736 endif 00737 00738 if (s%wetV(i,1) == 1) then 00739 dzs_s = .0_rKind 00740 dzs_n = dzb_s 00741 elseif (s%wetV(i ,1 ) == 0) then 00742 dzs_s = .0_rKind 00743 dzs_n = .5_rKind*(s%zs(i,1)-zsv(i,1))*ddyz(1) 00744 endif 00745 00746 s%ws(i,1) = - (s%hu(i,1)*s%uu(i,1)-s%hu(i-1,1 )*s%uu(i-1,1 ))*ddxz(i) & 00747 - (s%hv(i,1)*s%vv(i,1)-s%hv(i ,1 )*s%vv(i ,1 ))*ddyz(1) & 00748 + dzs_e*s%uu(i,1)+dzs_w*s%uu(i-1,1) & 00749 + dzs_s*s%vv(i,1)+dzs_n*s%vv(i, 1) 00750 endif 00751 enddo 00752 do i=imin_zs,imax_zs 00753 if (nonhZ(i,1) == 1) then 00754 s%wb(i,1) = awbr(i,1) + dp(i , 1) * awb(1,i,1) & 00755 + dp(i-1,1) * awb(2,i,1) + dp(i+1,1 ) * awb(3,i,1) & 00756 + dp(i ,1) * awb(4,i,1) + dp(i ,1 ) * awb(5,i,1) 00757 else 00758 if (s%wetU(i,1)*s%wetU(i-1,1) == 1) then 00759 dzb_e = .5_rKind*(zbu(i,1)-zbu(i-1,1))*ddxz(i) 00760 dzb_w = dzs_e 00761 elseif (s%wetU(i ,1) == 0) then 00762 dzb_e = .0_rKind 00763 dzb_w = .5_rKind*(s%zb(i,1)-zbu(i-1,1))*ddxz(i) 00764 elseif (s%wetU(i-1,1) == 0) then 00765 dzb_e = .5_rKind*(zbu(i,1)-s%zb(i,1)) *ddxz(i) 00766 dzb_w = .0_rKind 00767 else 00768 dzb_e = .0_rKind 00769 dzb_w = .0_rKind 00770 endif 00771 00772 if (s%wetV(i,1) == 1) then 00773 dzb_s = .0_rKind 00774 dzb_n = dzb_s 00775 elseif (s%wetV(i ,1 ) == 0) then 00776 dzb_s = .0_rKind 00777 dzb_n = .5_rKind*(s%zb(i,1)-zbv(i,1))*ddyz(1) 00778 endif 00779 00780 s%wb(i,1) = + dzb_e*s%uu(i,1)+dzb_w*s%uu(i-1,1) & 00781 + dzb_s*s%vv(i,1)+dzb_n*s%vv(i,1) 00782 endif 00783 enddo 00784 !Assign boundaries 00785 #ifdef USEMPI 00786 call xmpi_shift_ee(s%wb) 00787 call xmpi_shift_ee(s%ws) 00788 if (xmpi_istop) then 00789 s%wb(1,:) = s%wb(2,:) 00790 endif 00791 if (xmpi_isbot) then 00792 s%ws(s%nx+1,:) = s%ws(s%nx,:) 00793 s%wb(s%nx+1,:) = s%wb(s%nx,:) 00794 endif 00795 #else 00796 s%ws(s%nx+1,:) = s%ws(s%nx,:) 00797 s%wb(1,:) = s%wb(2,:) 00798 s%wb(s%nx+1,:) = s%wb(s%nx,:) 00799 #endif 00800 else 00801 do j=jmin_zs,jmax_zs 00802 do i=imin_zs,imax_zs 00803 if (nonhZ(i,j) == 1) then 00804 s%ws(i,j) = awsr(i,j) + dp(i , j) * aws(1,i,j) & 00805 + dp(i-1,j) * aws(2,i,j) + dp(i+1,j ) * aws(3,i,j) & 00806 + dp(i ,j-1)* aws(4,i,j) + dp(i ,j+1) * aws(5,i,j) 00807 00808 else 00809 !If the point is excluded in nh then get ws from the kinematic boundary condition 00810 if (s%wetU(i,j)*s%wetU(i-1,j) == 1) then 00811 dzs_e = .5_rKind*(zsu(i,j)-zsu(i-1,j))*ddxz(i) 00812 dzs_w = dzs_e 00813 elseif (s%wetU(i ,j) == 0) then 00814 dzs_e = .0_rKind 00815 dzs_w = .5_rKind*(s%zs(i,j)-zsu(i-1,j))*ddxz(i) 00816 elseif (s%wetU(i-1,j) == 0) then 00817 dzs_e = .5_rKind*(zsu(i,j)-s%zs(i,j)) *ddxz(i) 00818 dzs_w = .0_rKind 00819 else 00820 dzs_e = .0_rKind 00821 dzs_w = .0_rKind 00822 endif 00823 00824 if (s%wetV(i,j)*s%wetV(i,j-1) == 1) then 00825 dzs_s = .5_rKind*(zsv(i,j)-zsv(i,j-1))*ddyz(j) 00826 dzs_n = dzb_s 00827 elseif (s%wetV(i ,j ) == 0) then 00828 dzs_s = .0_rKind 00829 dzs_n = .5_rKind*(s%zs(i,j)-zsv(i,j-1))*ddyz(j) 00830 elseif (s%wetV(i ,j-1) == 0) then 00831 dzs_s = .5_rKind*(zsv(i,j)-s%zs(i,j)) *ddyz(j) 00832 dzs_n = .0_rKind 00833 else 00834 dzs_e = .0_rKind 00835 dzs_w = .0_rKind 00836 endif 00837 00838 s%ws(i,j) = - (s%hu(i,j)*s%uu(i,j)-s%hu(i-1,j )*s%uu(i-1,j ))*ddxz(i) & 00839 - (s%hv(i,j)*s%vv(i,j)-s%hv(i ,j-1)*s%vv(i ,j-1))*ddyz(j) & 00840 + dzs_e*s%uu(i,j)+dzs_w*s%uu(i-1,j) & 00841 + dzs_s*s%vv(i,j)+dzs_n*s%vv(i,j-1) 00842 endif 00843 enddo 00844 enddo 00845 00846 do j=jmin_zs,jmax_zs 00847 do i=imin_zs,imax_zs 00848 if (nonhZ(i,j) == 1) then 00849 s%wb(i,j) = awbr(i,j) + dp(i , j) * awb(1,i,j) & 00850 + dp(i-1,j) * awb(2,i,j) + dp(i+1,j ) * awb(3,i,j) & 00851 + dp(i ,j-1)* awb(4,i,j) + dp(i ,j+1) * awb(5,i,j) 00852 else 00853 if (s%wetU(i,j)*s%wetU(i-1,j) == 1) then 00854 dzb_e = .5_rKind*(zbu(i,j)-zbu(i-1,j))*ddxz(i) 00855 dzb_w = dzs_e 00856 elseif (s%wetU(i ,j) == 0) then 00857 dzb_e = .0_rKind 00858 dzb_w = .5_rKind*(s%zb(i,j)-zbu(i-1,j))*ddxz(i) 00859 elseif (s%wetU(i-1,j) == 0) then 00860 dzb_e = .5_rKind*(zbu(i,j)-s%zb(i,j)) *ddxz(i) 00861 dzb_w = .0_rKind 00862 else 00863 dzb_e = .0_rKind 00864 dzb_w = .0_rKind 00865 endif 00866 00867 if (s%wetV(i,j)*s%wetV(i,j-1) == 1) then 00868 dzb_s = .5_rKind*(zbv(i,j)-zbv(i,j-1))*ddyz(j) 00869 dzb_n = dzb_s 00870 elseif (s%wetV(i ,j ) == 0) then 00871 dzb_s = .0_rKind 00872 dzb_n = .5_rKind*(s%zb(i,j)-zbv(i,j-1))*ddyz(j) 00873 elseif (s%wetV(i ,j-1) == 0) then 00874 dzb_s = .5_rKind*(zbv(i,j)-s%zb(i,j)) *ddyz(j) 00875 dzb_n = .0_rKind 00876 else 00877 dzb_e = .0_rKind 00878 dzb_w = .0_rKind 00879 endif 00880 00881 s%wb(i,j) = + dzb_e*s%uu(i,j)+dzb_w*s%uu(i-1,j) & 00882 + dzb_s*s%vv(i,j)+dzb_n*s%vv(i,j-1) 00883 endif 00884 enddo 00885 enddo 00886 #ifdef USEMPI 00887 call xmpi_shift_ee(s%wb) 00888 call xmpi_shift_ee(s%ws) 00889 if (xmpi_istop) then 00890 s%wb(1,:) = s%wb(2,:) 00891 endif 00892 if (xmpi_isbot) then 00893 s%ws(s%nx+1,:) = s%ws(s%nx,:) 00894 s%wb(s%nx+1,:) = s%wb(s%nx,:) 00895 endif 00896 if (xmpi_isleft) then 00897 s%ws(:,1) = s%ws(:,2) 00898 s%wb(:,1) = s%wb(:,2) 00899 endif 00900 if (xmpi_isright) then 00901 s%wb(:,s%ny+1) = s%wb(:,s%ny) 00902 s%ws(:,s%ny+1) = s%ws(:,s%ny) 00903 endif 00904 #else 00905 !Assign boundaries 00906 s%ws(:,1) = s%ws(:,2) 00907 s%ws(:,s%ny+1) = s%ws(:,s%ny) 00908 !s%ws(1,:) = s%ws(2,:) 00909 s%ws(s%nx+1,:) = s%ws(s%nx,:) 00910 00911 s%wb(:,1) = s%wb(:,2) 00912 s%wb(:,s%ny+1) = s%wb(:,s%ny) 00913 s%wb(1,:) = s%wb(2,:) 00914 s%wb(s%nx+1,:) = s%wb(s%nx,:) 00915 #endif 00916 endif 00917 00918 Wm_old = .5_rKind*(s%ws+s%wb) 00919 end subroutine nonh_1lay_cor 00920 00921 subroutine nonh_1lay_pred(s,par) 00922 !============================================================================== 00923 ! 00924 00925 ! DATE AUTHOR CHANGES 00926 ! 00927 ! November 2010 Pieter Bart Smit New Subroutine 00928 00929 00930 !------------------------------------------------------------------------------ 00931 ! DECLARATIONS 00932 !------------------------------------------------------------------------------ 00933 00934 !-------------------------- PURPOSE ---------------------------- 00935 ! 00936 ! Include the pressure explicitly in the predictor step 00937 ! 00938 !-------------------------- DEPENDENCIES ---------------------------- 00939 use spaceparams 00940 use params, only: parameters 00941 use flow_secondorder_module, only: flow_secondorder_advw 00942 !-------------------------- ARGUMENTS ---------------------------- 00943 00944 type(spacepars) ,intent(inout) :: s 00945 type(parameters),intent(in) :: par 00946 00947 !-------------------------- LOCAL VARIABLES ---------------------------- 00948 00949 !Indices 00950 integer(kind=iKind) :: i,ie,iee,iw !Index variables 00951 integer(kind=iKind) :: j,js 00952 integer(kind=iKind) :: jmin,jmax 00953 00954 real(kind=rKind) :: dwdx1 !Gradient of vertical velocity in x-dir at i+1/2,j 00955 real(kind=rKind) :: dwdx2 !Gradient of vertical velocity in x-dir at i-1/2,j 00956 real(kind=rKind) :: dwdy1 !Gradient of vertical velocity in x-dir at i ,j+1/2 00957 real(kind=rKind) :: dwdy2 !Gradient of vertical velocity in x-dir at i ,j-1/2 00958 real(kind=rKind) :: Vol 00959 real(kind=rKind) :: wmax,reformfac 00960 00961 if (.not. initialized) then 00962 call nonh_init(s,par) 00963 endif 00964 00965 if (s%ny>0) then 00966 jmin = 2 00967 jmax = s%ny 00968 else 00969 jmin = 1 00970 jmax = 1 00971 endif 00972 00973 ! 00974 ! Determine whether or not to use the nh correction method (dry/wet/breaking etc.) 00975 ! 00976 if(par%useXBeachGSettings==0 .or. par%nhbreaker<=1) then 00977 call nonh_break( s , par ) 00978 else 00979 ! (1) The point is dry 00980 ! (2) The relative wave length kd of the smallest possible wave (L=2dx) is smaller than kdmin 00981 ! (3) The interpolated waterlevel zs is below the bottom (steep cliffs with overwash situations) 00982 ! (4) Where Miche breaker criterium applies -> bores are hydrostatic 00983 ! max steepness = H/L = maxbrsteep 00984 ! dz/dx = maxbrsteep 00985 ! dz/dx = dz/dt/c = w/c = w/sqrt(gh) 00986 ! wmax = maxbrsteep*sqrt(gh) 00987 do j=1,s%ny+1 00988 do i=1,s%nx+1 00989 iw = max(i,i-1) 00990 ie = min(s%nx,i+1) 00991 iee = min(s%nx,i+2) 00992 00993 if ( (s%wetU(i,j)==1 ) & 00994 .and. (0.5_rKind*(s%zs(i,j) + s%zs(ie,j)) > zbu(i,j) ) & 00995 .and. ( s%dsu(i,1)*par%kdmin/par%px < s%hum(i,j) ) ) then 00996 nonhU(i,j) = 1 00997 else 00998 nonhU(i,j) = 0 00999 endif 01000 enddo 01001 enddo 01002 01003 if (s%ny>2) then 01004 do j=1,s%ny+1 01005 js = min(s%ny,j+1) 01006 do i=1,s%nx+1 01007 if ( (s%wetV(i,j)==1 ) & 01008 .and. (0.5_rKind*(s%zs(i,j) + s%zs(i,js)) > zbv(i,j) ) & 01009 .and. ( s%dnv(1,j)*par%kdmin/par%px < s%hvm(i,j) ) ) then 01010 nonhV(i,j) = 1 01011 else 01012 nonhV(i,j) = 0 01013 endif 01014 enddo 01015 enddo 01016 else 01017 nonhV = 0 01018 endif 01019 ! 01020 ! Determine if a velocity point will be included in the nonh pressure matrix, include if 01021 ! any of the surrounding velocity points is included. 01022 ! 01023 if (s%ny>0) then 01024 do j=2,s%ny 01025 do i=2,s%nx 01026 if (max(nonhV(i,j),nonhV(i,j-1),nonhU(i,j),nonhU(i-1,j)) > 0) then 01027 nonhZ(i,j) = 1 01028 else 01029 nonhZ(i,j) = 0 01030 endif 01031 enddo 01032 enddo 01033 else 01034 do i=2,s%nx 01035 if (max(nonhU(i,1),nonhU(i-1,1)) > 0) then 01036 nonhZ(i,1) = 1 01037 else 01038 nonhZ(i,1) = 0 01039 endif 01040 enddo 01041 endif 01042 01043 if (par%nhbreaker == 2) then 01044 ! First determine local breaker criterion 01045 lbreakcond = par%maxbrsteep 01046 if (s%ny==0) then 01047 do i=2,s%nx 01048 if (s%breaking(i,1)==1) then 01049 lbreakcond(i-1:i+1,1) = par%secbrsteep 01050 endif 01051 enddo 01052 else 01053 do j=jmin,jmax 01054 do i=2,s%nx 01055 if (s%breaking(i,j)==1) then 01056 lbreakcond(i-1:i+1,j-1:j+1) = par%secbrsteep 01057 endif 01058 enddo 01059 enddo 01060 endif 01061 #ifdef USEMPI 01062 call xmpi_shift_ee(lbreakcond) 01063 #endif 01064 ! Now find areas where main breaking criterion is exceeded 01065 do j=jmin,jmax 01066 do i=2,s%nx 01067 if (s%breaking(i,j)==1) then 01068 if(s%ws(i,j)<=0.d0) then 01069 s%breaking(i,j) = 0 01070 endif 01071 else 01072 wmax = lbreakcond(i,j)*sqrt(par%g*s%hh(i,j)) ! add current term in here too 01073 if (s%ws(i,j)>=wmax) then 01074 s%breaking(i,j) = 1 01075 endif 01076 endif 01077 enddo 01078 enddo 01079 ! turn off non-hydrostatic pressure correction in areas with breaking and increase viscosity 01080 do j=jmin,jmax 01081 do i=2,s%nx 01082 if (s%breaking(i,j)==1) then 01083 nonhZ(i,j) = 0 01084 s%pres(i,j) = 0.d0 01085 endif 01086 enddo 01087 enddo 01088 elseif (par%nhbreaker == 3) then 01089 ! First determine local breaker criterion 01090 lbreakcond = par%maxbrsteep 01091 if (s%ny==0) then 01092 do i=2,s%nx 01093 if (s%breaking(i,1)==1) then 01094 lbreakcond(i-1:i+1,1) = par%secbrsteep 01095 endif 01096 enddo 01097 else 01098 do j=jmin,jmax 01099 do i=2,s%nx 01100 if (s%breaking(i,j)==1) then 01101 lbreakcond(i-1:i+1,j-1:j+1) = par%secbrsteep 01102 endif 01103 enddo 01104 enddo 01105 endif 01106 ! Now find areas where main breaking criterion is exceeded 01107 do j=jmin,jmax 01108 do i=2,s%nx 01109 s%wscrit(i,j) = lbreakcond(i,j)*sqrt(par%g*s%hh(i,j)) ! add current term in here too 01110 if (s%breaking(i,j)==1) then 01111 if(s%ws(i,j)<=s%wscrit(i,j)) then 01112 s%breaking(i,j) = 0 01113 endif 01114 else 01115 if (s%ws(i,j)>=s%wscrit(i,j)) then 01116 s%breaking(i,j) = 1 01117 endif 01118 endif 01119 enddo 01120 enddo 01121 ! turn off non-hydrostatic pressure correction in areas with breaking and increase viscosity 01122 do j=jmin,jmax 01123 do i=2,s%nx 01124 if (s%breaking(i,j)==1) then 01125 nonhZ(i,j) = 0 01126 s%pres(i,j) = 0.d0 01127 endif 01128 enddo 01129 enddo 01130 endif 01131 endif 01132 01133 !Calculate explicit part average vertical momentum (advection) 01134 if (s%ny>0) then 01135 do j=2,s%ny 01136 do i=2,s%nx 01137 if (nonhZ(i,j) == 1) then 01138 Wm(i,j) = Wm_old(i,j) - par%dt*( ddxu(i-1)*max(s%qx(i-1,j ),0.0_rKind)*(Wm_old(i ,j ) & 01139 - Wm_old(i-1,j ))/s%hh(i,j) & 01140 + ddxu(i) *min(s%qx(i ,j ),0.0_rKind)*(Wm_old(i+1,j )-Wm_old(i ,j ))/s%hh(i,j) & 01141 + ddyv(j-1)*max(s%qy(i ,j-1),0.0_rKind)*(Wm_old(i ,j )-Wm_old(i ,j-1))/s%hh(i,j) & 01142 + ddyv(j )*min(s%qy(i ,j ),0.0_rKind)*(Wm_old(i ,j+1)-Wm_old(i ,j ))/s%hh(i,j) ) 01143 else 01144 Wm(i,j) = 0.0_rKind 01145 Wm_old(i,j) = 0.0_rKind 01146 s%ws(i,j) = 0.0_rKind 01147 s%wb(i,j) = 0.0_rKind 01148 s%pres(i,j) = 0.0_rKind 01149 endif 01150 enddo 01151 enddo 01152 else 01153 do i=2,s%nx 01154 if (nonhZ(i,1) == 1) then 01155 Wm(i,1) = Wm_old(i,1) - par%dt*( ddxu(i-1)*max(s%qx(i-1,1 ),0.0_rKind)*(Wm_old(i ,1 ) & 01156 -Wm_old(i-1,1 ))/s%hh(i,1) & 01157 + ddxu(i) *min(s%qx(i ,1 ),0.0_rKind)*(Wm_old(i+1,1 )-Wm_old(i ,1 ))/s%hh(i,1) ) 01158 else 01159 Wm(i,1) = 0.0_rKind 01160 Wm_old(i,1) = 0.0_rKind 01161 s%ws(i,1) = 0.0_rKind 01162 s%wb(i,1) = 0.0_rKind 01163 s%pres(i,1) = 0.0_rKind 01164 endif 01165 enddo 01166 endif 01167 01168 !Calculate explicit part vertical viscocity 01169 if (s%ny>0) then 01170 do j=2,s%ny 01171 do i=2,s%nx 01172 dwdx1 = .5d0*(s%nuh(i-1,j )+s%nuh(i,j ))*s%hu(i-1,j )*(Wm_old(i ,j )-Wm_old(i-1,j ))*ddxu(i-1) 01173 dwdx2 = .5d0*(s%nuh(i+1,j )+s%nuh(i,j ))*s%hu(i ,j )*(Wm_old(i+1,j )-Wm_old(i ,j ))*ddxu(i) 01174 dwdy1 = s%nuh(i ,j-1)*s%hu(i ,j-1)*(Wm_old(i ,j )-Wm_old(i ,j-1))*ddyv(j-1) 01175 dwdy2 = s%nuh(i ,j )*s%hu(i ,j )*(Wm_old(i ,j+1)-Wm_old(i ,j ))*ddyv(j) 01176 Wm(i,j) = Wm(i,j) + (1.0_rKind/s%hh(i,j))*par%dt*(dwdx2-dwdx1)*ddxz(i)*real(s%wetu(i,j)*s%wetu(i-1,j),rKind) 01177 + (1.0_rKind/s%hh(i,j))*par%dt*(dwdy2-dwdy1)*ddyz(j)*real(s%wetv(i,j)*s%wetv(i,j-1),rKind) 01178 enddo 01179 enddo 01180 else 01181 do i=2,s%nx 01182 dwdx1 = .5d0*(s%nuh(i-1,1 )+s%nuh(i,1 ))*s%hu(i-1,1 )*(Wm_old(i ,1 )-Wm_old(i-1,1 ))*ddxu(i-1) 01183 dwdx2 = .5d0*(s%nuh(i+1,1 )+s%nuh(i,1 ))*s%hu(i ,1 )*(Wm_old(i+1,1 )-Wm_old(i ,1 ))*ddxu(i) 01184 dwdy1 = 0.d0 01185 dwdy2 = 0.d0 01186 Wm(i,1) = Wm(i,1) + (1.0_rKind/s%hh(i,1))*par%dt*(dwdx2-dwdx1)*ddxz(i)*real(s%wetu(i,1)*s%wetu(i-1,1),rKind) 01187 enddo 01188 endif 01189 01190 do j=jmin,jmax 01191 do i=2,s%nx 01192 if (nonhU(i,j)==1) then 01193 vol = 0.5_rKind*par%dt/(s%hum(i,j)*dxu(i)) 01194 au(1,i,j) = - (s%zs(i+1,j) - s%zb(i ,j))*vol 01195 au(0,i,j) = + (s%zs(i ,j) - s%zb(i+1,j))*vol 01196 else 01197 au(1,i,j) = 0.0_rKind 01198 au(0,i,j) = 0.0_rKind 01199 endif 01200 enddo 01201 enddo 01202 if (xmpi_istop) then !wwvv: added test on istop 01203 au(:,1,:) = 0.0_rKind 01204 endif 01205 if (xmpi_isbot) then !wwvv: added test on isbot 01206 au(:,s%nx+1,:) = 0.0_rKind 01207 endif 01208 01209 01210 !Built pressure coefficients V 01211 !call timer_start(timer_flow_nonh_av) 01212 if (s%ny>2) then 01213 do j=2,s%ny 01214 do i=2,s%nx 01215 if (nonhV(i,j)==1)then 01216 vol = 0.5_rKind*par%dt/(s%hvm(i,j)*dyv(j)) 01217 av(1,i,j) = -(s%zs(i ,j+1) - s%zb(i ,j ))*vol 01218 av(0,i,j) = +(s%zs(i ,j ) - s%zb(i ,j+1))*vol 01219 avr(i,j) = s%vv(i,j) 01220 else 01221 av(1,i,j) = 0.0_rKind 01222 av(0,i,j) = 0.0_rKind 01223 endif 01224 enddo 01225 enddo 01226 if (xmpi_isleft) then !wwvv: added test on isleft 01227 av(:,:,1) = 0.0_rKind 01228 endif 01229 if (xmpi_isright) then !wwvv: added test on isright 01230 av(:,:,s%ny+1) = 0.0_rKind 01231 endif 01232 endif 01233 01234 !Include explicit approximation for pressure in s%uu and s%vv and Wm 01235 if (par%secorder == 1) then 01236 do j=jmin,jmax 01237 do i=imin_uu,imax_uu 01238 s%uu(i,j) = s%uu(i,j) + au(1,i,j) * s%pres(i+1,j) + au(0,i,j) * s%pres(i ,j) 01239 enddo 01240 enddo 01241 01242 if (s%ny>0) then 01243 do j=jmin_vv,jmax_vv 01244 do i=imin_vv,imax_vv 01245 s%vv(i,j) = s%vv(i,j) + av(1,i,j) * s%pres(i,j+1) + av(0,i,j) * s%pres(i,j ) 01246 enddo 01247 enddo 01248 else 01249 do i=imin_vv,imax_vv 01250 s%vv(i,1) = s%vv(i,1) + av(1,i,1) * s%pres(i,1) + av(0,i,1) * s%pres(i,1 ) 01251 enddo 01252 endif 01253 01254 do j=jmin_zs,jmax_zs 01255 do i=imin_zs,imax_zs 01256 if (nonhZ(i,j) == 1) then 01257 Wm(i,j) = Wm(i,j) + Wcoef(i,j)*par%dt * s%pres(i,j)/s%hh(i,j) 01258 endif 01259 enddo 01260 enddo 01261 call flow_secondorder_advW(s,par,Wm,Wm_old, nonhZ) 01262 endif 01263 01264 end subroutine nonh_1lay_pred 01265 01266 01267 01268 01269 !============================================================================== 01270 ! REDUCED NON-HYDROSTATIC ROUTINES 2DV 01271 !============================================================================== 01272 ! 01273 ! 01274 !============================================================================== 01275 subroutine nonh_2lay_cor_2dV(s,par) 01276 !============================================================================== 01277 ! 01278 01279 ! DATE AUTHOR CHANGES 01280 ! 01281 ! October 2014 Pieter Bart Smit New subroutine 01282 ! 01283 !------------------------------------------------------------------------------- 01284 ! DECLARATIONS 01285 !------------------------------------------------------------------------------- 01286 01287 !-------------------------- PURPOSE ---------------------------- 01288 01289 ! In this subroutine we calculate the nonh-hydrostatic pressures by solving 01290 ! a discrete poisson equation, and subsequently correct the velocity field 01291 ! so that the velocities at the new timelevel are divergence free. 01292 ! 01293 ! This subroutine uses the pressure coeficients calculated previously in 01294 ! the routine nonh_2lay_pred_2dV, and for each of the velocity components all 01295 ! processes have been accounted for (advection, bed-friction etc), including 01296 ! an explict estimation of the nonh press known from the previous timestep 01297 ! (unless we are calculating with first order accuracy). 01298 ! 01299 ! Subsequently we will take the following steps: 01300 ! 01301 ! 1) Calculate the water and bed levels at u-points (subroutine call to zuzv). 01302 ! 2) Built the poisson matrix by substituting the pressure dependencies of the 01303 ! momentum equations into the continuity equation. 01304 ! 3) Solve the resulting system (call to solver_solvemat). 01305 ! 4) Update the non-hydrostatic pressure. If second-order is active, dp is the 01306 ! difference between the old and the new pressure, so in that case we have 01307 ! 01308 ! s%press = s%press + dp 01309 ! 01310 ! else, it represents the new pressure directly, s%press = dp. 01311 ! 5) Correct the velocity components with the nonhydrostatic pressures. The 01312 ! new velocity field will be divergence free. 01313 ! 01314 ! After these steps, control is returned to flow_timestep. 01315 ! 01316 !-------------------------- DEPENDENCIES ---------------------------- 01317 ! 01318 use spaceparams 01319 use params, only: parameters 01320 use solver_module, only: solver_solvemat 01321 use xmpi_module 01322 ! 01323 !-------------------------- ARGUMENTS ---------------------------- 01324 ! 01325 type(spacepars) ,intent(inout) :: s 01326 type(parameters),intent(in) :: par 01327 ! 01328 !-------------------------- LOCAL PARAMETERS ---------------------------- 01329 ! 01330 integer(kind=iKind),parameter :: j = 1 !Index variables 01331 ! 01332 !-------------------------- LOCAL VARIABLES ---------------------------- 01333 ! 01334 01335 !Indices 01336 integer(kind=iKind) :: i,k !Index variables 01337 01338 !Work 01339 real(kind=rKind) :: dzsdx ! Free surface gradient (s-dir) 01340 real(kind=rKind) :: dzbdx ! Bed level gradient (s-dir) 01341 real(kind=rKind) :: alpha ! Parameter controlling how the velocity at the 01342 ! layer interface is determined ( =-wcoef for energy 01343 ! conservative behaviour). 01344 real(kind=rKind) :: alpha_u(0:1) ! The layer distribution parameter interpolated at 01345 ! the east (=1) and west (=0) faces. 01346 real(kind=rKind) :: hfU(0:1) ! Coefficient of the west U(i-1,j) (=0) or east U(i,j) (=1) velocity 01347 ! in the reduced continuity equation. 01348 real(kind=rKind) :: hfdU(0:1) ! Coefficient of the west dU(i-1,j) (=0) or east dU(i,j) (=1) velocity 01349 ! 01350 !------------------------------------------------------------------------------- 01351 ! IMPLEMENTATION 01352 !------------------------------------------------------------------------------- 01353 ! 01354 01355 ! 01356 ! ------ Substitute in the continuity equation ------ 01357 ! 01358 do i=2,s%nx 01359 ! 01360 if (nonhZ(i,j)==1) then 01361 ! 01362 do k = 0,1 01363 ! 01364 alpha_u(k) = .5d0 * (wcoef(i+k,j) + wcoef(i-1+k,j)) 01365 01366 hfU(k) = s%hu(i-1+k ,j) * ( 1.0d0 + alpha_u(k) ) 01367 hfdU(k) = s%hu(i-1+k ,j) * alpha_u(k) * ( 1.0d0 - alpha_u(k) ) 01368 ! 01369 enddo 01370 ! 01371 dzsdx = .5*( zsu(i,j) - zsu(i-1,j) ) 01372 dzbdx = .5*( zbu(i,j) + alpha_u(1) * s%hu(i,j) - zbu(i-1,j) - alpha_u(0) * s%hu(i-1,j) ) 01373 ! 01374 01375 01376 ! 01377 alpha = -wcoef(i,j) 01378 ! 01379 mat(1,i,j) = 2.0d0* s%dsz(i,j) * aws(1,i,j) & 01380 + hfU( 1)*au( 0,i,j) - hfU( 0)*au(1, i-1,j) & 01381 + hfdU(1)*adU(0,i,j) - hfdU(0)*adU(1,i-1,j) & 01382 - dzsdx * ( au(0,i,j) + au(1,i-1,j) - ( adU(0,i,j) + adU(1,i-1,j) )* wcoef(i,j) ) & 01383 - dzbdx * ( au(0,i,j) + au(1,i-1,j) + ( adU(0,i,j) + adU(1,i-1,j) )* alpha ) 01384 ! 01385 mat(2,i,j) = - hfU( 0)*au (0, i-1,j) - hfdU( 0)*adU(0, i-1,j) & 01386 - dzsdx * ( au(0,i-1,j) - adU(0,i-1,j) * wcoef(i,j) ) & 01387 - dzbdx * ( au(0,i-1,j) + adU(0,i-1,j) * alpha ) 01388 ! 01389 mat(3,i,j) = hfU(1)*au (1,i ,j) + hfdU(1)*adU(1,i,j) & 01390 - dzsdx * ( au(1,i,j) - adU(1,i,j) * wcoef(i,j) ) & 01391 - dzbdx * ( au(1,i,j) + adU(1,i,j) * alpha ) 01392 ! 01393 rhs(i,j) = - 2.0d0 * s%dsz(i,j) * Wm(i,j) & 01394 - hfU( 1)*s%uu( i ,j) + hfU( 0)*s%uu( i-1 ,j) & 01395 - hfdU(1)*s%dU( i ,j) + hfdU(0)*s%dU( i-1 ,j) & 01396 + dzsdx * ( s%uu(i,j) + s%uu(i-1,j) - ( s%dU(i,j) + s%dU(i-1,j) )* wcoef(i,j) ) & 01397 + dzbdx * ( s%uu(i,j) + s%uu(i-1,j) + ( s%dU(i,j) + s%dU(i-1,j) )* alpha ) 01398 else 01399 ! 01400 mat(1 ,i,j) = 1.0_rKind 01401 mat(2:5,i,j) = 0.0_rKind 01402 rhs(i,j) = 0.0_rKind 01403 ! 01404 endif 01405 ! 01406 enddo 01407 01408 ! 01409 !---------------- Solve Linear System ----------------- 01410 ! 01411 dp = 0.0_rKind 01412 call solver_solvemat( mat , rhs , dp , s%nx, s%ny,par) 01413 ! 01414 if (par%secorder == 1) then 01415 ! 01416 s%pres = s%pres + dp 01417 ! 01418 else 01419 ! 01420 s%pres = dp 01421 ! 01422 endif 01423 01424 ! 01425 !------------------ Correct u/v/w ------------------- 01426 ! 01427 !U 01428 do i=2,s%nx-1 01429 ! 01430 s%uu(i,j) = s%uu(i,j) + au(1,i,j)*dp(i+1,j)+au(0,i,j)*dp(i,j) 01431 ! 01432 enddo 01433 01434 !dU 01435 if (par%nhlay > 0. ) then 01436 ! 01437 do i=2,s%nx-1 01438 ! 01439 s%dU(i,j) = s%dU(i,j) + adU(1,i,j)*dp(i+1,j)+adU(0,i,j)*dp(i,j) 01440 ! 01441 enddo 01442 ! 01443 endif 01444 ! 01445 01446 #ifdef USEMPI 01447 ! Communicate velocities 01448 call xmpi_shift_ee(s%uu) 01449 call xmpi_shift_ee(s%dU) 01450 #endif 01451 01452 ! 01453 ! Calculate the surface vertical velocity from the kinematic condition 01454 ! 01455 ! 01456 do i=2,s%nx 01457 ! 01458 if ( s%wetz(i,j)==1 ) then 01459 ! 01460 dzsdx = .5*( zsu(i,j) - zsu(i-1,j ) ) / s%dsz(i,j) 01461 01462 s%ws(i,j) = - (s%hu(i,j)*s%uu(i,j) -s%hu(i-1,j )*s%uu(i-1,j )) / s%dsz(i,j) & 01463 + dzsdx * ( s%uu(i,j) + s%uu(i-1,j ) - ( s%dU(i,j) + s%dU(i-1,j ) ) * wcoef(i,j) ) 01464 ! 01465 else 01466 ! 01467 s%ws(i,j) = 0.0d0 01468 ! 01469 endif 01470 ! 01471 enddo 01472 ! 01473 01474 ! 01475 ! Calculate the bottom velocity from the kinematic condition 01476 ! 01477 ! 01478 do i=2,s%nx 01479 ! 01480 if ( s%wetz(i,j)==1 ) then 01481 ! 01482 dzbdx = .5*( zbu(i,j) - zbu(i-1,j ) ) / s%dsz(i,j) 01483 ! 01484 s%wb(i,j) = + dzbdx * ( s%uu(i,j) + s%uu(i-1,j ) + ( s%dU(i,j) + s%dU(i-1,j ) ) * (1-wcoef(i,j)) ) 01485 ! 01486 else 01487 ! 01488 s%wb(i,j) = 0.0d0 01489 ! 01490 endif 01491 ! 01492 enddo 01493 ! 01494 01495 ! 01496 do i=2,s%nx 01497 ! 01498 ! 01499 if ( nonhZ(i,j) == 1) then 01500 ! 01501 Wm(i,j) = Wm(i,j) + aws(1,i,j)*dp(i,j) 01502 ! 01503 else 01504 ! 01505 Wm(i,j) = 0.0d0 01506 ! 01507 endif 01508 ! 01509 enddo 01510 ! 01511 01512 if (par%nhlay > 0. ) then 01513 ! 01514 01515 do i=2,s%nx 01516 ! 01517 ! 01518 if ( nonhZ(i,j) == 1 ) then 01519 ! 01520 Wm0(i,j) = Wm(i,j) + 0.5d0* ( s%wb(i,j) - s%ws(i,j) ) 01521 ! 01522 else 01523 ! 01524 Wm0(i,j) = 0. 01525 ! 01526 endif 01527 ! 01528 enddo 01529 ! 01530 endif 01531 01532 01533 01534 #ifdef USEMPI 01535 call xmpi_shift_ee(Wm) 01536 call xmpi_shift_ee(s%dU) 01537 01538 if (xmpi_istop) then 01539 Wm(1,:) = Wm(2,:) 01540 endif 01541 if (xmpi_isbot) then 01542 Wm(s%nx+1,:) = Wm(s%nx,:) 01543 endif 01544 #else 01545 !Assign boundaries 01546 Wm(1,:) = Wm(2,:) 01547 Wm(s%nx+1,:) = Wm(s%nx,:) 01548 01549 #endif 01550 01551 ! 01552 ! Calculate the relative vertical advection 01553 ! 01554 do i=2,s%nx 01555 ! 01556 if (s%wetz(i,j) == 1) then 01557 ! 01558 omega(i,j) = - ( s%hu(i ,j) * dU0(i ,j) - s%hu(i-1,j) * dU0(i-1,j) ) / s%dsz( i , j ) 01559 ! 01560 else 01561 ! 01562 omega(i,j) = 0.0d0 01563 ! 01564 endif 01565 ! 01566 enddo 01567 01568 ! 01569 Wm_old = Wm 01570 !s%ws = s%dU 01571 dU0 = s%dU 01572 !1111111 01573 end subroutine nonh_2lay_cor_2dv 01574 01575 ! 01576 subroutine nonh_2lay_pred_2dV(s,par,uu0) 01577 !============================================================================== 01578 ! 01579 01580 ! DATE AUTHOR CHANGES 01581 ! 01582 ! October 2014 Pieter Bart Smit New Subroutine 01583 01584 01585 !------------------------------------------------------------------------------ 01586 ! DECLARATIONS 01587 !------------------------------------------------------------------------------ 01588 01589 !-------------------------- PURPOSE ---------------------------- 01590 ! 01591 ! Include the pressure explicitly in the predictor step. In this subroutine 01592 ! the following steps are taken: 01593 ! 01594 ! 1) Calculate whether or not breaking is active (handled in subroutine 01595 ! nonh_break ). 01596 ! 2) Calculate the coeficients of the nonh. pres. as they occur in the momentum 01597 ! equations. (these are use to calculate the explicit pressure contributions 01598 ! here, and the implicit contributions in the correction routine later). 01599 ! 3) Calculate the first order additional contributions to the advection due to 01600 ! the reduced two-layer formulation (if applicable). 01601 ! 4) Include an explicit estimate for the nonh press. in the momentum equations 01602 ! using the earlier calculated coef. (if second-order approximations for the 01603 ! advection are used, else this step is not necessary). 01604 ! 5) Include the contributions due to internal stress terms and bottom friction 01605 ! in the equation for the velocity difference (implicitly for stability). 01606 ! 01607 ! After these steps, control flows back to flow_timestep, and - if active - 01608 ! second-order approximations are calculated. Consequently, flow_timestep calls 01609 ! nonh_cor again and the nonhydrostatic computations are continued in 01610 ! nonh_2lay_cor. 01611 ! 01612 !-------------------------- DEPENDENCIES ---------------------------- 01613 use spaceparams 01614 use params, only: parameters 01615 use flow_secondorder_module, only: flow_secondorder_advw 01616 !-------------------------- ARGUMENTS ---------------------------- 01617 01618 type(spacepars) ,intent(inout) :: s 01619 type(parameters),intent(in) :: par 01620 real(kind=rkind),dimension(s%nx+1,s%ny+1),intent(in) :: uu0 !The velocity at the previous timestep 01621 01622 !-------------------------- LOCAL PARAMETERS ---------------------------- 01623 01624 integer(kind=iKind), parameter :: j = 1 !Y-index, constant in 1d 01625 01626 !-------------------------- LOCAL VARIABLES ---------------------------- 01627 !General 01628 real(kind=rKind) :: fac !A factor 01629 01630 !Indices 01631 integer(kind=iKind) :: i,k !Loop indices 01632 01633 !Used in calculation of pressure gradients 01634 real(kind=rKind) :: Vol !The volume of a momentum cell 01635 01636 !Used in calculation of advection 01637 real(kind=rKind) :: qx(0:1) !Discharge through the faces of a momentum cell (0: left face 01638 !1: right face) 01639 real(kind=rKind) :: faceval(0:1) !Value of dU at the faces of a momentum cell (0: left face 01640 !1: right face) 01641 01642 !Used in calculation of vertical exchange of momentum 01643 real(kind=rKind) :: fric !Bottom friction coeficient (= dt * cf * |U| / h) 01644 real(kind=rKind) :: stress !Coeficient used in calculation of stress 01645 real(kind=rKind) :: ub !The (Eulerian) velocity at the bottom 01646 real(kind=rKind) :: umag !The (Eulerian) magnitude of the velocity vector at the bottom 01647 real(kind=rKind) :: ve !The eddy viscosity used for vertical mixing 01648 real(kind=rKind) :: advec !Coeficient used in calculation of vertical advection 01649 real(kind=rKind) :: Amat,rhs !The 1*1 "matrix" and RHS of the vertical "system". 01650 01651 !------------------------------------------------------------------------------ 01652 ! IMPLEMENTATION 01653 !------------------------------------------------------------------------------ 01654 01655 ! 01656 ! Determine whether or not to use the nh correction method (dry/wet/breaking etc.) 01657 ! 01658 call nonh_break( s , par ) 01659 call zuzv(s) 01660 ! 01661 ! -- Calculate the pressure dependencies for U -- 01662 ! 01663 do i=2,s%nx 01664 ! 01665 if (nonhU(i,j)==1) then 01666 ! 01667 vol = 0.5d0 * par%dt/(s%hum(i,j)*s%dsu(i,j)) 01668 fac = (s%zb(i+1,j) - s%zb(i ,j)) * vol 01669 ! 01670 au(1,i,j) = - vol * (1.d0 + wcoef(i+1,j) ) * s%hh(i+1,j) - fac 01671 au(0,i,j) = + vol * (1.d0 + wcoef(i ,j) ) * s%hh(i,j) - fac 01672 ! 01673 else 01674 ! 01675 au(1,i,j) = 0.0_rKind 01676 au(0,i,j) = 0.0_rKind 01677 ! 01678 endif 01679 ! 01680 enddo 01681 01682 ! 01683 ! -- Calculate the pressure dependencies for dU -- 01684 ! 01685 01686 if ( par%nhlay > 0. ) then 01687 ! 01688 do i=2,s%nx 01689 ! 01690 if (nonhU(i,j)*(1-s%breaking(i,j))*(1-s%breaking(i-1,j))==1 ) then 01691 ! 01692 vol = 0.5_rKind*par%dt/(s%hum(i,j)*s%dsu(i,j)) 01693 ! 01694 fac = ( s%zs( i+1,j ) - s%zs(i,j) ) * vol / ( 1.0d0 - .5d0*wcoef(i,j) - .5d0*wcoef(i+1,j) ) & 01695 - ( wcoef( i+1,j ) - wcoef(i,j) ) * s%hum(i,j) * vol / (2.0d0 - wcoef(i,j) - wcoef(i+1,j) ) 01696 ! 01697 adU(1,i,j) = - s%hh(i+1,j)*vol + fac 01698 adU(0,i,j) = + s%hh(i,j)*vol + fac 01699 ! 01700 else 01701 ! 01702 adU(1,i,j) = 0.0_rKind 01703 adU(0,i,j) = 0.0_rKind 01704 ! 01705 endif 01706 ! 01707 enddo 01708 ! 01709 endif 01710 01711 ! 01712 ! -- Calculate the pressure dependencies for Wm -- 01713 ! 01714 01715 do i=2,s%nx 01716 ! 01717 01718 if (nonhZ(i,j) == 1) then 01719 ! 01720 aws(1,i,j) = par%dt/( s%hh(i,j) * ( 1.0d0 - wcoef( i , j ) ) ) 01721 ! 01722 else 01723 ! 01724 aws(1,i,j) = 0.0_rKind 01725 ! 01726 endif 01727 enddo 01728 ! 01729 ! 01730 !----------------- Advection Wm (FOU) ------------------- 01731 ! 01732 do i=2,s%nx 01733 ! 01734 if ( nonhZ(i,j) == 1) then 01735 ! 01736 01737 do k = 0,1 01738 ! 01739 qx(k) = ( s%qx(i - 1 + k , j ) ) 01740 01741 if (qx(k) > 0.) then 01742 ! 01743 faceval(k) = Wm_old( i - 1 + k ,j ) * real( nonhZ(i - 1 + k,j) , kind=rKind ) 01744 ! 01745 else 01746 ! 01747 faceval(k) = Wm_old( i + k ,j ) * real( nonhZ(i + k,j) , kind=rKind ) 01748 ! 01749 endif 01750 ! 01751 enddo 01752 01753 fac = par%dt / s%hh(i,j) / s%dsz(i,j) 01754 01755 Wm(i,j) = Wm_old(i,j) + fac * Wm_old(i,j)*( s%qx(i,j) - s%qx(i-1,j) ) - fac*( qx(1) * faceval(1) - qx(0) * faceval(0) ) 01756 ! 01757 else 01758 ! 01759 Wm(i,j) = 0.0_rKind 01760 Wm_old(i,j) = 0.0_rKind 01761 s%pres(i,j) = 0.0_rKind 01762 ! 01763 endif 01764 ! 01765 enddo 01766 ! 01767 !----------------- Advection dU (FOU) ------------------- 01768 ! 01769 if ( par%nhlay > 0. ) then 01770 ! 01771 do i=2,s%nx-1 01772 ! 01773 if ( s%wetu( i,j) == 1) then 01774 ! 01775 do k = 0,1 01776 ! 01777 qx(k) = .5 * ( s%qx(i - 1 + k , j ) + s%qx(i + k , j ) ) 01778 01779 if (qx(k) > 0.) then 01780 ! 01781 faceval(k) = dU0( i - 1 + k , j ) 01782 ! 01783 else 01784 ! 01785 faceval(k) = dU0( i+ k , j ) 01786 ! 01787 endif 01788 ! 01789 enddo 01790 01791 fac = par%dt / s%hum(i,j) / s%dsu(i,j) 01792 01793 s%dU(i,j) = dU0(i,j) + .5 * fac *dU0(i,j) * ( s%qx(i+1,j) - s%qx(i-1,j) ) & 01794 - fac * ( qx(1)*faceval(1) - qx(0)*faceval(0) ) 01795 ! 01796 else 01797 ! 01798 s%dU(i,j) = 0.0 01799 ! 01800 endif 01801 ! 01802 enddo 01803 ! 01804 endif 01805 ! 01806 01807 ! 01808 if (xmpi_istop) then 01809 ! 01810 au(:,1,:) = 0.0_rKind 01811 adU(:,1,:) = 0.0_rKind 01812 ! 01813 endif 01814 ! 01815 if (xmpi_isbot) then 01816 ! 01817 au(:,s%nx,:) = 0.0_rKind 01818 adU(:,s%nx,:) = 0.0_rKind 01819 ! 01820 endif 01821 ! 01822 if (par%nhlay > 0.) then 01823 ! 01824 do i=2,s%nx 01825 ! 01826 if ( nonhZ(i,j) == 1 ) then 01827 ! 01828 ! Vertical advection (First order upwind) 01829 ! 01830 01831 if ( omega(i,j) > 0.) then 01832 ! 01833 s%wm(i,j) = s%wm(i,j) + par%dt * omega(i,j) * wm0(i,j)/s%hh(i,j) * wcoef(i,j) 01834 ! 01835 else 01836 ! 01837 s%wm(i,j) = s%wm(i,j) + par%dt * omega(i,j) * wm_old(i,j)/s%hh(i,j) * wcoef(i,j) 01838 endif 01839 endif 01840 01841 enddo 01842 endif 01843 01844 !Include explicit approximation for pressure in s%uu du and Wm 01845 if (par%secorder == 1) then 01846 ! 01847 do i=2,s%nx-1 01848 ! 01849 s%uu(i,j) = s%uu(i,j) + au(1,i,j) * s%pres(i+1,j) + au(0,i,j) * s%pres(i ,j) 01850 ! 01851 enddo 01852 ! 01853 01854 if ( par%nhlay > 0. ) then 01855 ! 01856 do i=2,s%nx-1 01857 ! 01858 s%dU(i,j) = s%dU(i,j) + adU(1,i,j) * s%pres(i+1,j) + adU(0,i,j) * s%pres(i ,j) 01859 ! 01860 enddo 01861 ! 01862 endif 01863 01864 ! 01865 do i=2,s%nx 01866 ! 01867 Wm(i,j) = Wm(i,j) + aws(1,i,j) * s%pres(i,j) 01868 ! 01869 enddo 01870 ! 01871 01872 #ifdef USEMPI 01873 ! 01874 call xmpi_shift_ee(Wm) 01875 #endif 01876 ! 01877 call flow_secondorder_advW(s,par,Wm,Wm_old, nonhZ) 01878 ! 01879 #ifdef USEMPI 01880 ! 01881 call xmpi_shift_ee(Wm) 01882 ! 01883 #endif 01884 endif 01885 ! 01886 01887 if (par%nhlay > 0.) then 01888 ! 01889 ! 01890 ! 01891 ! 01892 !----------------- Bottom stress, internal stresses and vertical advection ------------------- 01893 ! 01894 ! Integration of these terms occurs 01895 ! implicity to ensure that they do not introduce any additional stability 01896 ! concerns. 01897 ! 01898 ! Moreover, the vertical exchange of momentum due to advection between the 01899 ! layers also influences dU (but not U), so we have 01900 ! 01901 ! the equation reads dU - dU* 01902 ! -------- = dU ( fric + stress ) + s%uu * fric - omega * U_face 01903 ! dt 01904 ! 01905 do i=2,s%nx-1 01906 ! 01907 if ( s%wetu(i,j) == 1 ) then 01908 ! 01909 ! 01910 ! Bottom friction 01911 ! 01912 ub = uu0(i,j) + (1.-par%nhlay) * dU0(i,j) 01913 umag = abs(ub) 01914 fric = par%dt * s%cfu(i,j) * umag / s%hu(i,j) 01915 01916 ! 01917 ! Internal stresses 01918 ! 01919 ve = 1.0e-4 !Vertical eddy viscocity 01920 stress = 2.* par%dt * ve / ( s%hum(i,j)**2 * (par%nhlay - par%nhlay**2) ) 01921 01922 ! 01923 ! Vertical advection 01924 ! 01925 advec = 0.5d0* ( omega(i,j) + omega(i+1,j) )* par%dt / s%hum(i,j) 01926 01927 ! 01928 ! Built rhs/diagonal 01929 ! 01930 fac = 0.5d0*( wcoef(i,j) + wcoef(i+1,j) ) 01931 rhs = s%dU(i,j) - uu0(i,j) * ( fric + advec ) 01932 Amat = 1.0d0 + (1.0d0-fac) * fric + stress & 01933 + (1.0d0 - fac) * max( advec , 0. ) & 01934 - fac * min( advec , 0. ) 01935 01936 ! 01937 ! "Solve" system 01938 ! 01939 s%dU(i,j) = rhs / Amat 01940 ! 01941 endif 01942 ! 01943 enddo 01944 ! 01945 endif 01946 01947 #ifdef USEMPI 01948 ! 01949 call xmpi_shift_ee(s%dU) 01950 call xmpi_shift_ee(s%dV) 01951 ! 01952 #endif 01953 ! 01954 end subroutine nonh_2lay_pred_2dv 01955 01956 01957 !============================================================================== 01958 ! REDUCED NON-HYDROSTATIC ROUTINES 3D 01959 !============================================================================== 01960 ! 01961 01962 01963 ! 01964 !============================================================================== 01965 subroutine nonh_2lay_cor_3d(s,par) 01966 !============================================================================== 01967 ! 01968 01969 ! DATE AUTHOR CHANGES 01970 ! 01971 ! October 2010 Pieter Bart Smit New Subroutine 01972 ! November 2010 Pieter Bart Smit Added explicit prediction for 2nd order scheme 01973 01974 !------------------------------------------------------------------------------- 01975 ! DECLARATIONS 01976 !------------------------------------------------------------------------------- 01977 01978 !-------------------------- PURPOSE ---------------------------- 01979 01980 ! In this subroutine we calculate the nonh-hydrostatic pressures by solving 01981 ! a discrete poisson equation, and subsequently correct the velocity field 01982 ! so that the velocities at the new timelevel are divergence free. 01983 ! 01984 ! This subroutine uses the pressure coeficients calculated previously in 01985 ! the routine nonh_2lay_pred_3d, and for each of the velocity components all 01986 ! processes have been accounted for (advection, bed-friction etc), including 01987 ! an explict estimation of the nonh press known from the previous timestep 01988 ! (unless we are calculating with first order accuracy). 01989 ! 01990 ! Subsequently we will take the following steps: 01991 ! 01992 ! 1) Calculate the water and bed levels at u-points (subroutine call to zuzv). 01993 ! 2) Built the poisson matrix by substituting the pressure dependencies of the 01994 ! momentum equations into the continuity equation. 01995 ! 3) Solve the resulting system (call to solver_solvemat). 01996 ! 4) Update the non-hydrostatic pressure. If second-order is active, dp is the 01997 ! difference between the old and the new pressure, so in that case we have 01998 ! 01999 ! s%press = s%press + dp 02000 ! 02001 ! else, it represents the new pressure directly, s%press = dp. 02002 ! 5) Correct the velocity components with the nonhydrostatic pressures. The 02003 ! new velocity field will be divergence free. 02004 ! 02005 ! After these steps, control is returned to flow_timestep. 02006 ! 02007 !-------------------------- DEPENDENCIES ---------------------------- 02008 use spaceparams 02009 use params, only: parameters 02010 use solver_module, only: solver_solvemat 02011 use xmpi_module 02012 !-------------------------- ARGUMENTS ---------------------------- 02013 02014 type(spacepars) ,intent(inout) :: s 02015 type(parameters),intent(in) :: par 02016 02017 !-------------------------- LOCAL VARIABLES ---------------------------- 02018 02019 !Indices 02020 integer(kind=iKind) :: i,j,k !Index variables 02021 02022 !Work 02023 real(kind=rKind) :: dzsdx ! Free surface gradient (s-dir) 02024 real(kind=rKind) :: dzbdx ! Bed level gradient (s-dir) 02025 real(kind=rKind) :: dzsdy ! Free surface gradient (n-dir) 02026 real(kind=rKind) :: dzbdy ! Bed level gradient (n-dir) 02027 real(kind=rKind) :: alpha ! Parameter controlling how the velocity at the 02028 ! layer interface is determined ( =-wcoef for energy 02029 ! conservative behaviour). 02030 real(kind=rKind) :: alpha_u(0:1) ! The layer distribution parameter interpolated at 02031 ! the east (=1) and west (=0) faces. 02032 real(kind=rKind) :: alpha_v(0:1) ! The layer distribution parameter interpolated at 02033 ! the north (=1) and south (=0) faces. 02034 real(kind=rKind) :: hfU(0:1) ! Coefficient of the west U(i-1,j) (=0) or east U(i,j) (=1) velocity 02035 ! in the reduced continuity equation. 02036 real(kind=rKind) :: hfdU(0:1) ! Coefficient of the west dU(i-1,j) (=0) or east dU(i,j) (=1) velocity 02037 ! difference in the reduced continuity equation. 02038 real(kind=rKind) :: hfV(0:1) ! Coefficient of the south V(i,j-1) (=0) or North V(i,j) (=1) velocity 02039 ! in the reduced continuity equation. 02040 real(kind=rKind) :: hfdV(0:1) ! Coefficient of the south dV(i,j-1) (=0) or North dV(i,j) (=1) velocity 02041 ! difference in the reduced continuity equation. 02042 02043 ! 02044 !------------------------------------------------------------------------------- 02045 ! IMPLEMENTATION 02046 !------------------------------------------------------------------------------- 02047 ! 02048 02049 ! 02050 !Substitute in the continuity equation 02051 ! 02052 do j=2,s%ny 02053 ! 02054 do i=2,s%nx 02055 ! 02056 if (nonhZ(i,j)==1) then 02057 ! 02058 do k = 0,1 02059 ! 02060 alpha_u(k) = .5d0 * (wcoef(i+k,j) + wcoef(i-1+k,j)) 02061 02062 hfU(k) = s%dsdnzi(i,j) * s%dnu(i-1+k,j) * s%hu(i-1+k ,j) * ( 1.0d0 + alpha_u(k) ) 02063 hfdU(k) = s%dsdnzi(i,j) * s%dnu(i-1+k,j) * s%hu(i-1+k ,j) * alpha_u(k) * ( 1.0d0 - alpha_u(k) ) 02064 ! 02065 enddo 02066 02067 do k = 0,1 02068 ! 02069 alpha_v(k) = .5d0 * (wcoef(i,j+k) + wcoef(i,j-1+k)) 02070 02071 hfV(k) = s%dsdnzi(i,j) * s%dsv(i,j-1+k) * s%hv(i,j-1+k ) * ( 1.0d0 + alpha_v(k) ) 02072 hfdV(k) = s%dsdnzi(i,j) * s%dsv(i,j-1+k) * s%hv(i,j-1+k ) * alpha_v(k) * ( 1.0d0 - alpha_v(k) ) 02073 ! 02074 enddo 02075 ! 02076 02077 if (nonhu(i,j)*nonhu(i-1,j) == 1 ) then 02078 ! 02079 dzsdx = .5*( zsu(i,j) - zsu(i-1,j ) ) / s%dsz(i,j) 02080 dzbdx = .5*( zbu(i,j) + alpha_u(1) * s%hu(i,j) - zbu(i-1,j ) - alpha_u(0) * s%hu(i-1,j ) ) / s%dsz(i,j) 02081 ! 02082 else 02083 dzsdx = 0.0d0 02084 dzbdx = 0.0 02085 endif 02086 02087 if (nonhv(i,j)*nonhv(i,j-1) == 1 ) then 02088 ! 02089 dzsdy = .5*( zsv(i,j) - zsv(i ,j-1) ) / s%dnz(i,j) 02090 dzbdy = .5*( zbv(i,j) + alpha_v(1) * s%hv(i,j) - zbv(i ,j-1) - alpha_v(0) * s%hv(i ,j-1) ) / s%dnz(i,j) 02091 else 02092 ! 02093 dzsdy = 0.0d0 02094 dzbdy = 0.0d0 02095 ! 02096 endif 02097 02098 alpha = -wcoef(i,j) 02099 ! 02100 ! Main diagonal ( i,j ) 02101 ! 02102 mat(1,i,j) = 2.0d0 * aws(1,i,j) & 02103 + hfU(1) * au(0, i ,j) - hfU(0) * au(1, i-1,j ) & 02104 + hfV(1) * av(0, i ,j) - hfV(0) * av(1, i ,j-1) & 02105 + hfdU(1) * adU(0,i ,j) - hfdU(0) * adU(1,i-1 ,j ) & 02106 + hfdV(1) * adV(0,i ,j) - hfdV(0) * adV(1,i ,j-1) & 02107 - dzsdx * ( au( 0,i,j) + au( 1,i-1,j ) - ( adU(0,i,j) + adU(1,i-1,j ) ) * wcoef(i,j) ) & 02108 - dzsdy * ( av( 0,i,j) + av( 1,i ,j-1) - ( adV(0,i,j) + adV(1,i ,j-1) ) * wcoef(i,j) ) & 02109 - dzbdx * ( au( 0,i,j) + au( 1,i-1,j ) + ( adU(0,i,j) + adU(1,i-1,j ) ) * alpha ) & 02110 - dzbdy * ( av( 0,i,j) + av( 1,i ,j-1) + ( adV(0,i,j) + adV(1,i ,j-1) ) * alpha ) 02111 ! 02112 ! West diagonal ( i-1,j ) 02113 ! 02114 mat(2,i,j) = - hfU(0) * au (0, i-1,j) & 02115 - hfdU(0) * adU(0, i-1,j) & 02116 - dzsdx * ( au(0,i-1,j) - adU(0,i-1,j) * wcoef(i,j) ) & 02117 - dzbdx * ( au(0,i-1,j) + adU(0,i-1,j) * alpha ) 02118 ! 02119 ! South diagonal (i,j-1 ) 02120 ! 02121 mat(4,i,j) = - hfV(0) * av (0, i,j-1) & 02122 - hfdV(0) * adV(0, i,j-1) & 02123 - dzsdy * ( av(0,i,j-1) - adV(0,i,j-1) * wcoef(i,j) ) & 02124 - dzbdy * ( av(0,i,j-1) + adV(0,i,j-1) * alpha ) 02125 ! 02126 ! East diagonal ( i+1,j ) 02127 ! 02128 mat(3,i,j) = hfU(1) * au (1,i ,j) & 02129 + hfdU(1) * adU(1,i ,j) & 02130 - dzsdx * ( au(1,i,j) - adU(1,i,j) * wcoef(i,j) ) & 02131 - dzbdx * ( au(1,i,j) + adU(1,i,j) * alpha ) 02132 ! 02133 ! North diagonal ( i,j+1 ) 02134 ! 02135 mat(5,i,j) = hfV(1) * av (1,i ,j) & 02136 + hfdV(1) * adV(1,i ,j) & 02137 - dzsdy * ( av(1,i,j) - adV(1,i,j) * wcoef(i,j) ) & 02138 - dzbdy * ( av(1,i,j) + adV(1,i,j) * alpha ) 02139 ! 02140 ! RHS 02141 ! 02142 rhs(i,j) = - 2.0d0 * Wm(i,j) & 02143 - hfU( 1 ) * s%uu( i ,j) + hfU( 0 ) * s%uu( i-1 ,j ) & 02144 - hfdU( 1 ) * s%dU( i ,j) + hfdU( 0 ) * s%dU( i-1 ,j ) & 02145 - hfV( 1 ) * s%vv( i ,j) + hfV( 0 ) * s%vv( i ,j-1) & 02146 - hfdV( 1 ) * s%dV( i ,j) + hfdV( 0 ) * s%dV( i ,j-1) & 02147 + dzsdx * ( s%uu(i,j) + s%uu(i-1,j ) - ( s%dU(i,j) + s%dU(i-1,j ) ) * wcoef(i,j) ) & 02148 + dzsdy * ( s%vv(i,j) + s%vv(i ,j-1) - ( s%dV(i,j) + s%dV(i ,j-1) ) * wcoef(i,j) ) & 02149 + dzbdx * ( s%uu(i,j) + s%uu(i-1,j ) + ( s%dU(i,j) + s%dU(i-1,j ) ) * alpha ) & 02150 + dzbdy * ( s%vv(i,j) + s%vv(i ,j-1) + ( s%dV(i,j) + s%dV(i ,j-1) ) * alpha ) 02151 ! 02152 else 02153 ! 02154 mat(1 ,i,j) = 1.0_rKind 02155 mat(2:5,i,j) = 0.0_rKind 02156 rhs(i,j) = 0.0_rKind 02157 ! 02158 endif 02159 enddo 02160 enddo 02161 02162 !------------------ Solve matrix ----------------- 02163 dp = 0.0_rKind 02164 call solver_solvemat( mat , rhs , dp , s%nx, s%ny,par) 02165 02166 if ( par%secorder == 1 ) then 02167 ! 02168 ! If second-order is active, we are calculating with the pressure difference 02169 ! 02170 s%pres = s%pres + dp 02171 ! 02172 else 02173 ! 02174 s%pres = dp 02175 ! 02176 endif 02177 02178 02179 !Correct u/v/w 02180 02181 !U 02182 do j=jmin_uu,jmax_uu 02183 ! 02184 do i=imin_uu,imax_uu 02185 ! 02186 s%uu(i,j) = s%uu(i,j) + au(1,i,j)*dp(i+1,j)+au(0,i,j)*dp(i,j) 02187 ! 02188 enddo 02189 ! 02190 enddo 02191 02192 !V 02193 do j=jmin_vv,jmax_vv 02194 ! 02195 do i=imin_vv,imax_vv 02196 ! 02197 s%vv(i,j) = s%vv(i,j) + av(1,i,j)*dp(i,j+1)+av(0,i,j)*dp(i,j) 02198 ! 02199 enddo 02200 ! 02201 enddo 02202 ! 02203 02204 ! 02205 if (par%nhlay > 0.) then 02206 ! 02207 !dU 02208 ! 02209 do j=jmin_uu,jmax_uu 02210 ! 02211 do i=imin_uu,imax_uu 02212 ! 02213 s%dU(i,j) = s%dU(i,j) + adU(1,i,j)*dp(i+1,j)+adU(0,i,j)*dp(i,j) 02214 ! 02215 enddo 02216 ! 02217 enddo 02218 ! 02219 !dV 02220 ! 02221 do j=jmin_vv,jmax_vv 02222 ! 02223 do i=imin_vv,imax_vv 02224 ! 02225 s%dV(i,j) = s%dV(i,j) + adV(1,i,j)*dp(i,j+1)+adV(0,i,j)*dp(i,j) 02226 ! 02227 enddo 02228 ! 02229 enddo 02230 ! 02231 endif 02232 ! 02233 02234 ! 02235 ! Communicate results... 02236 ! 02237 #ifdef USEMPI 02238 call xmpi_shift_ee(s%uu) 02239 call xmpi_shift_ee(s%vv) 02240 call xmpi_shift_ee(s%dU) 02241 call xmpi_shift_ee(s%dV) 02242 #endif 02243 02244 ! 02245 ! Calculate the surface vertical velocity from the kinematic condition 02246 ! 02247 do j=jmin_zs,jmax_zs 02248 ! 02249 do i=imin_zs,imax_zs 02250 ! 02251 if ( s%wetz(i,j)==1 ) then 02252 ! 02253 dzsdx = .5*( zsu(i,j) - zsu(i-1,j ) ) / s%dsz(i,j) 02254 dzsdy = .5*( zsv(i,j) - zsv(i ,j-1) ) / s%dnz(i,j) 02255 ! 02256 s%ws(i,j) = - (s%dnu(i,j)*s%hu(i,j)*s%uu(i,j) -s%dnu(i-1,j)*s%hu(i-1,j )*s%uu(i-1,j )) * s%dsdnzi(i,j) & 02257 - (s%dsv(i,j)*s%hv(i,j)*s%vv(i,j) -s%dsv(i,j-1)*s%hv( i,j-1)*s%vv(i ,j-1)) * s%dsdnzi(i,j) & 02258 + dzsdx * ( s%uu(i,j) + s%uu(i-1,j ) - ( s%dU(i,j) + s%dU(i-1,j ) ) * wcoef(i,j) ) & 02259 + dzsdy * ( s%vv(i,j) + s%vv(i ,j-1) - ( s%dV(i,j) + s%dV(i ,j-1) ) * wcoef(i,j) ) 02260 ! 02261 else 02262 ! 02263 s%ws(i,j) = 0.0d0 02264 ! 02265 endif 02266 ! 02267 enddo 02268 ! 02269 enddo 02270 02271 ! 02272 ! Calculate the bottom velocity from the kinematic condition 02273 ! 02274 do j=jmin_zs,jmax_zs 02275 ! 02276 do i=imin_zs,imax_zs 02277 ! 02278 if ( s%wetz(i,j)==1 ) then 02279 ! 02280 dzbdx = .5*( zbu(i,j) - zbu(i-1,j ) ) 02281 dzbdy = .5*( zbv(i,j) - zbv(i ,j-1) ) 02282 02283 s%wb(i,j) = + dzbdx * ( s%uu(i,j) + s%uu(i-1,j ) + ( s%dU(i,j) + s%dU(i-1,j ) ) * (1- wcoef(i,j)) ) / s%dsz(i,j) & 02284 + dzbdy * ( s%vv(i,j) + s%vv(i ,j-1) + ( s%dV(i,j) + s%dV(i ,j-1) ) * (1- wcoef(i,j)) ) / s%dnz(i,j) 02285 ! 02286 else 02287 ! 02288 s%wb(i,j) = 0.0d0 02289 ! 02290 endif 02291 ! 02292 enddo 02293 ! 02294 enddo 02295 02296 ! 02297 ! Calculate the mean vertical velocity in the top layer 02298 ! 02299 do j=jmin_zs,jmax_zs 02300 ! 02301 do i=imin_zs,imax_zs 02302 ! 02303 if ( nonhZ(i,j) == 1) then 02304 ! 02305 Wm(i,j) = Wm(i,j) + aws(1,i,j)*dp(i,j) 02306 ! 02307 else 02308 ! 02309 Wm(i,j) = 0.0d0 02310 ! 02311 endif 02312 ! 02313 enddo 02314 ! 02315 enddo 02316 02317 ! 02318 ! Calculate the relative vertical advection 02319 ! 02320 do j=jmin_zs,jmax_zs 02321 ! 02322 do i=imin_zs,imax_zs 02323 ! 02324 if (s%wetz(i,j)==1) then 02325 ! 02326 omega(i,j) = - ( s%dnu(i,j) * s%hu(i ,j) * dU0(i ,j) - s%dnu(i-1,j) * s%hu(i-1,j) * dU0(i-1,j) & 02327 + s%dsv(i,j) * s%hv(i ,j) * dV0(i ,j) - s%dsv(i,j-1) * s%hv(i,j-1) * dV0(i,j-1) ) * s%dsdnzi( i , j ) 02328 ! 02329 else 02330 ! 02331 omega(i,j) = 0.0d0 02332 ! 02333 endif 02334 ! 02335 enddo 02336 ! 02337 enddo 02338 02339 ! 02340 ! Calculate the mean vertical velocity in the lower layer 02341 ! 02342 do j=jmin_zs,jmax_zs 02343 ! 02344 do i=imin_zs,imax_zs 02345 ! 02346 ! 02347 if ( s%wetz(i,j)==1 ) then 02348 ! 02349 Wm0(i,j) = Wm(i,j) + .5d0 * ( s%wb(i,j) - s%ws(i,j) ) 02350 ! 02351 else 02352 ! 02353 Wm0(i,j) = 0. 02354 ! 02355 endif 02356 ! 02357 enddo 02358 ! 02359 enddo 02360 02361 02362 #ifdef USEMPI 02363 if (xmpi_istop) then 02364 Wm(1,:) = Wm(2,:) 02365 endif 02366 if (xmpi_isbot) then 02367 Wm(s%nx+1,:) = Wm(s%nx,:) 02368 endif 02369 if (xmpi_isleft) then 02370 Wm(1,:) = Wm(2,:) 02371 endif 02372 if (xmpi_isright) then 02373 Wm(:,s%ny+1) = Wm(:,s%ny) 02374 endif 02375 ! 02376 call xmpi_shift_ee(Wm) 02377 ! 02378 #else 02379 !Assign boundaries 02380 Wm(1,:) = Wm(2,:) 02381 Wm(s%nx+1,:) = Wm(s%nx,:) 02382 Wm(:,1) = Wm(:,2) 02383 Wm(:,s%ny+1) = Wm(:,s%ny) 02384 02385 #endif 02386 Wm_old = Wm 02387 dU0 = s%dU 02388 dV0 = s%dV 02389 ! 02390 end subroutine nonh_2lay_cor_3d 02391 02392 02393 subroutine nonh_2lay_pred_3d(s,par,uu0,vv0) 02394 !============================================================================== 02395 ! 02396 02397 ! DATE AUTHOR CHANGES 02398 ! 02399 ! October 2014 Pieter Bart Smit New Subroutine 02400 02401 02402 !------------------------------------------------------------------------------ 02403 ! DECLARATIONS 02404 !------------------------------------------------------------------------------ 02405 02406 !-------------------------- PURPOSE ---------------------------- 02407 ! 02408 ! Include the pressure explicitly in the predictor step. In this subroutine 02409 ! the following steps are taken: 02410 ! 02411 ! 1) Calculate whether or not breaking is active (handled in subroutine 02412 ! nonh_break ). 02413 ! 2) Calculate the coeficients of the nonh. pres. as they occur in the momentum 02414 ! equations. (these are use to calculate the explicit pressure contributions 02415 ! here, and the implicit contributions in the correction routine later). 02416 ! 3) Calculate the first order additional contributions to the advection due to 02417 ! the reduced two-layer formulation (if applicable). 02418 ! 4) Include an explicit estimate for the nonh press. in the momentum equations 02419 ! using the earlier calculated coef. (if second-order approximations for the 02420 ! advection are used, else this step is not necessary). 02421 ! 5) Include the contributions due to internal stress terms and bottom friction 02422 ! in the equation for the velocity difference (implicitly for stability). 02423 ! 02424 ! After these steps, control flows back to flow_timestep, and - if active - 02425 ! second-order approximations are calculated. Consequently, flow_timestep calls 02426 ! nonh_cor again and the nonhydrostatic computations are continued in 02427 ! nonh_2lay_cor. 02428 ! 02429 !-------------------------- DEPENDENCIES ---------------------------- 02430 use spaceparams 02431 use params, only: parameters 02432 use flow_secondorder_module, only: flow_secondorder_advw 02433 !-------------------------- ARGUMENTS ---------------------------- 02434 02435 type(spacepars) ,intent(inout) :: s 02436 type(parameters),intent(in) :: par 02437 real(kind=rkind),dimension(s%nx+1,s%ny+1),intent(in) :: uu0 02438 real(kind=rkind),dimension(s%nx+1,s%ny+1),intent(in) :: vv0 02439 02440 !-------------------------- VARIABLES ---------------------------- 02441 02442 !General 02443 real(kind=rKind) :: Fac !A factor 02444 02445 !Indices 02446 integer(kind=iKind) :: i,j,k !Loop indices 02447 02448 !Used in calculation of pressure gradients 02449 real(kind=rKind) :: Vol !The volume of a momentum cell 02450 02451 !Used in calculation of advection dU/dV/Wm 02452 real(kind=rKind) :: qx(0:1) !Discharge through the faces of a momentum cell (0: left face 02453 !1: right face) 02454 real(kind=rKind) :: qy(0:1) !Discharge through the faces of a momentum cell (0: bottom face 02455 !1: top face) 02456 real(kind=rKind) :: faceval_x(0:1) !Value of dU at the faces of a momentum cell (0: left face 02457 !1: right face) 02458 real(kind=rKind) :: faceval_y(0:1) !Value of dU at the faces of a momentum cell (0: left face 02459 !1: right face) 02460 real(kind=rKind) :: DU_face 02461 real(kind=rKind) :: DV_face 02462 real(kind=rKind) :: cosd,sind !cosine or sine of the difference angle 02463 real(kind=rKind) :: sig(0:1) !sign function 02464 02465 !Used in calculation of vertical exchange of momentum 02466 real(kind=rKind) :: fric !Bottom friction coeficient (= dt * cf * |U| / h) 02467 real(kind=rKind) :: stress !Coeficient used in calculation of stress 02468 real(kind=rKind) :: ub !The (Eulerian) u-velocity at the bottom 02469 real(kind=rKind) :: vb !The (Eulerian) v-velocity at the bottom 02470 real(kind=rKind) :: umag !The (Eulerian) magnitude of the velocity vector at the bottom 02471 real(kind=rKind) :: ve !The eddy viscosity used for vertical mixing 02472 real(kind=rKind) :: advec !Coeficient used in calculation of vertical advection 02473 real(kind=rKind) :: Amat,rhs !The 1*1 "matrix" and RHS of the vertical "system". 02474 02475 !------------------------------------------------------------------------------ 02476 ! IMPLEMENTATION 02477 !------------------------------------------------------------------------------ 02478 02479 ! 02480 ! Determine whether or not to use the nh correction method (dry/wet/breaking etc.) 02481 ! 02482 call nonh_break( s , par ) 02483 call zuzv(s) 02484 02485 ! 02486 ! Calculate the pressure dependencies for U 02487 ! 02488 do j=2,s%ny 02489 ! 02490 do i=1,s%nx 02491 ! 02492 if (nonhU(i,j)==1) then 02493 ! 02494 vol = 0.5d0 * par%dt/(s%hum(i,j)*s%dsu(i,j)) 02495 Fac = (s%zb(i+1,j) - s%zb(i ,j)) * vol 02496 au(1,i,j) = - vol * (1.d0 + wcoef(i,j) ) * s%hh(i+1,j) - Fac 02497 au(0,i,j) = + vol * (1.d0 + wcoef(i,j) ) * s%hh(i,j) - Fac 02498 ! 02499 else 02500 ! 02501 au(1,i,j) = 0.0_rKind 02502 au(0,i,j) = 0.0_rKind 02503 ! 02504 endif 02505 ! 02506 enddo 02507 ! 02508 enddo 02509 ! 02510 ! Calculate the pressure dependencies for V 02511 ! 02512 do j=1,s%ny 02513 ! 02514 do i=2,s%nx 02515 ! 02516 if (nonhV(i,j)==1) then 02517 ! 02518 vol = 0.5d0 * par%dt/(s%hvm(i,j)*s%dnv(i,j)) 02519 Fac = (s%zb(i,j+1) - s%zb(i ,j)) * vol 02520 av(1,i,j) = - vol * (1.d0 + wcoef(i,j+1) ) * s%hh(i,j+1) - Fac 02521 av(0,i,j) = + vol * (1.d0 + wcoef(i,j ) ) * s%hh(i,j) - Fac 02522 ! 02523 else 02524 ! 02525 av(1,i,j) = 0.0_rKind 02526 av(0,i,j) = 0.0_rKind 02527 ! 02528 endif 02529 ! 02530 enddo 02531 ! 02532 enddo 02533 02534 if ( par%nhlay > 0.0d0 ) then 02535 ! 02536 ! Calculate the pressure dependencies for dU 02537 ! 02538 do j=2,s%ny 02539 ! 02540 do i=1,s%nx 02541 ! 02542 if (nonhU(i,j)*(1-s%breaking(i,j))*(1-s%breaking(max(i-1,1),j))==1 ) then 02543 ! 02544 vol = 0.5_rKind*par%dt/(s%hum(i,j)*s%dsu(i,j)) 02545 02546 fac = ( s%zs( i+1,j ) - s%zs(i,j) ) * vol / ( 1.0d0 - .5d0*wcoef(i,j) - .5d0*wcoef(i+1,j) ) & 02547 - ( wcoef( i+1,j ) - wcoef(i,j) ) * s%hum(i,j) * vol / (2.0d0 - wcoef(i,j) - wcoef(i+1,j) ) 02548 02549 adU(1,i,j) = - s%hh(i+1,j)*vol + Fac 02550 adU(0,i,j) = + s%hh(i,j)*vol + Fac 02551 ! 02552 else 02553 ! 02554 adU(1,i,j) = 0.0_rKind 02555 adU(0,i,j) = 0.0_rKind 02556 ! 02557 endif 02558 ! 02559 enddo 02560 ! 02561 enddo 02562 ! 02563 ! Calculate the pressure dependencies for dV 02564 ! 02565 do j=1,s%ny 02566 ! 02567 do i=2,s%nx 02568 ! 02569 if (nonhV(i,j)*(1-s%breaking(i,j))*(1-s%breaking(i,max(j-1,1)))==1 ) then 02570 ! 02571 vol = 0.5_rKind*par%dt/(s%hvm(i,j)*s%dnv(i,j)) 02572 fac = ( s%zs( i,j ) - s%zs(i,j+1) ) * vol / ( 1.0d0 - .5d0*wcoef(i,j) - .5d0*wcoef(i,j+1) ) & 02573 - ( wcoef( i,j ) - wcoef(i,j+1) ) * s%hvm(i,j) * vol / (2.0d0 - wcoef(i,j) - wcoef(i,j+1) ) 02574 adV(1,i,j) = - s%hh(i,j+1)*vol + Fac 02575 adV(0,i,j) = + s%hh(i,j)*vol + Fac 02576 ! 02577 else 02578 ! 02579 adV(1,i,j) = 0.0_rKind 02580 adV(0,i,j) = 0.0_rKind 02581 ! 02582 endif 02583 ! 02584 enddo 02585 ! 02586 enddo 02587 ! 02588 endif 02589 02590 ! 02591 !Built pressure coefficients W 02592 ! 02593 do j=2,s%ny 02594 ! 02595 do i=2,s%nx 02596 ! 02597 if (nonhZ(i,j) == 1) then 02598 ! 02599 aws(1,i,j) = par%dt/( s%hh(i,j) * ( 1 - wcoef(i,j) ) ) 02600 ! 02601 else 02602 ! 02603 aws(1,i,j) = 0.0_rKind 02604 ! 02605 endif 02606 ! 02607 enddo 02608 ! 02609 enddo 02610 ! 02611 02612 ! 02613 if (xmpi_istop) then 02614 ! 02615 au(:,1,:) = 0.0_rKind 02616 adU(:,1,:) = 0.0_rKind 02617 ! 02618 endif 02619 ! 02620 if (xmpi_isbot) then 02621 ! 02622 au(:,s%nx,:) = 0.0_rKind 02623 adU(:,s%nx,:) = 0.0_rKind 02624 ! 02625 endif 02626 ! 02627 02628 ! 02629 if (xmpi_isleft) then 02630 ! 02631 av(:,:,1) = 0.0_rKind 02632 adV(:,:,1) = 0.0_rKind 02633 ! 02634 endif 02635 ! 02636 02637 ! 02638 if (xmpi_isright) then 02639 ! 02640 av(:,:,s%ny) = 0.0_rKind 02641 adV(:,:,s%ny) = 0.0_rKind 02642 ! 02643 endif 02644 02645 ! 02646 !----------------- Advection Wm (FOU) ------------------- 02647 ! 02648 ! d ( h U Wm ) d ( h V Wm ) 02649 ! ------------ + ------------ 02650 ! ds dn 02651 ! 02652 ! Advection according to Stelling & Duinmeijer 2003, momentum conservative 02653 ! 02654 do j=jmin_zs,jmax_zs 02655 ! 02656 do i=imin_zs,imax_zs 02657 ! 02658 if ( nonhZ(i,j)==1 ) then 02659 ! 02660 do k = 0,1 02661 ! 02662 qx(k) = s%qx(i - 1 + k , j ) * s%dnu(i - 1 + k , j ) 02663 qy(k) = s%qy(i , j - 1 + k) * s%dsv(i , j - 1 + k ) 02664 02665 if (qx(k) > 0.) then 02666 ! 02667 faceval_x(k) = Wm_old( i - 1 + k ,j ) * real( nonhZ(i - 1 + k,j) , kind=rKind ) 02668 ! 02669 else 02670 ! 02671 faceval_x(k) = Wm_old( i + k ,j ) * real( nonhZ(i + k,j) , kind=rKind ) 02672 ! 02673 endif 02674 02675 if (qy(k) > 0.) then 02676 ! 02677 faceval_y(k) = Wm_old( i ,j - 1 + k ) * real( nonhZ(i,j - 1 + k) , kind=rKind ) 02678 ! 02679 else 02680 ! 02681 faceval_y(k) = Wm_old( i,j + k ) * real( nonhZ(i,j + k) , kind=rKind ) 02682 ! 02683 endif 02684 ! 02685 enddo 02686 02687 fac = par%dt / s%hh(i,j) * s%dsdnzi(i,j) 02688 02689 Wm(i,j) = Wm_old(i,j) + fac * Wm_old(i,j) * ( qx(1) + qy(1) - qx(0) - qy(0) ) & 02690 - fac * ( qx(1) * faceval_x(1) + qy(1) * faceval_y(1) & 02691 - qx(0) * faceval_x(0) - qy(0) * faceval_y(0) ) 02692 ! 02693 else 02694 ! 02695 Wm(i,j) = 0.0_rKind 02696 Wm_old(i,j) = 0.0_rKind 02697 s%pres(i,j) = 0.0_rKind 02698 ! 02699 endif 02700 ! 02701 enddo 02702 ! 02703 enddo 02704 02705 if ( par%nhlay > 0.0d0 ) then 02706 ! 02707 ! 02708 !----------------- Advection dU (FOU) ------------------- 02709 ! 02710 ! d ( h U dU ) d ( h V dU ) 02711 ! ---------- + ---------- 02712 ! dx dy 02713 ! 02714 ! Advection according to Stelling & Duinmeijer 2003, momentum conservative (second term in their equation 24) 02715 ! 02716 sig(0) = -1 02717 sig(1) = 1 02718 ! 02719 do j=jmin_uu,jmax_uu 02720 ! 02721 do i=imin_uu,imax_uu 02722 ! 02723 if ( s%wetu( i,j) == 1) then 02724 ! 02725 ! Calculate transport velocities / face values 02726 ! 02727 do k = 0,1 02728 ! 02729 qx(k) = .5 * ( s%qx(i - 1 + k , j ) + s%qx(i + k , j ) ) * s%dnz(i+k,j) 02730 qy(k) = .5 * ( s%qy(i+1, j - 1 + k ) + s%qy(i, j - 1 + k ) ) * s%dsc(i,j+k) 02731 ! 02732 02733 dv_face = 0.5d0* ( dV0( i + k , j - 1 ) + dV0( i + k , j ) ) 02734 02735 cosd = cos( s%alfau(i+k,j) - s%alfau(i+k-1,j) ) 02736 sind = sig(k)*sin( s%alfau(i+k,j) - s%alfaz(i+k-1,j) ) 02737 ! 02738 if (qx(k) > 0.) then 02739 ! 02740 faceval_x(k) = dU0( i - 1 + k , j ) * cosd + dV_face * sind 02741 ! 02742 else 02743 ! 02744 faceval_x(k) = dU0( i+ k , j ) * cosd + dV_face * sind 02745 ! 02746 endif 02747 ! 02748 02749 ! 02750 cosd = cos( s%alfau(i,j+k-1) - s%alfau(i,j+k) ) 02751 sind = sig(k)*sin( s%alfau(i,j+k-1) - s%alfau(i,j+k) ) 02752 02753 dV_face = 0.5d0* ( dV0( i , j + k - 1) + dV0( i + 1 , j + k - 1 ) ) 02754 if (qy(k) > 0.) then 02755 ! 02756 faceval_y(k) = dU0( i , j -1 + k ) * cosd + dV_face*sind 02757 ! 02758 else 02759 ! 02760 faceval_y(k) = dU0( i , j + k) * cosd + dV_face*sind 02761 ! 02762 endif 02763 ! 02764 enddo 02765 02766 fac = par%dt / s%hum(i,j) * s%dsdnui(i,j) 02767 02768 s%dU(i,j) = dU0(i,j) + fac *dU0(i,j) * ( qx(1) + qy(1) - qx(0) - qy(0) ) & 02769 - fac *( qx(1)*faceval_x(1) + qy(1)*faceval_y(1) & 02770 - qx(0)*faceval_x(0) - qy(0)*faceval_y(0) ) 02771 ! 02772 else 02773 ! 02774 s%dU(i,j) = 0.0 02775 ! 02776 endif 02777 ! 02778 enddo 02779 ! 02780 enddo 02781 02782 ! 02783 !----------------- Advection dV (FOU) ------------------- 02784 ! 02785 ! d ( h V dV ) d ( h U dV ) 02786 ! ---------- + ---------- 02787 ! dy dx 02788 ! 02789 ! 02790 ! Advection according to Stelling & Duinmeijer 2003, momentum conservative (second term in their equation 25) 02791 ! 02792 do j=jmin_vv,jmax_vv 02793 ! 02794 do i=imin_vv,imax_vv 02795 ! 02796 if ( s%wetv( i,j) == 1) then 02797 ! 02798 ! Calculate transport velocities / face values 02799 ! 02800 do k = 0,1 02801 ! 02802 qy(k) = .5 * ( s%qy(i , j- 1 + k ) + s%qy(i , j + k ) ) * s%dsz(i,j+k) 02803 qx(k) = .5 * ( s%qx(i- 1 + k, j +1 ) + s%qx(i - 1 + k , j) ) * s%dnc(i+k,j) 02804 ! 02805 02806 cosd = cos( s%alfav(i+k-1,j) - s%alfav(i+k,j) ) 02807 sind = sig(k)*sin( s%alfav(i+k-1,j) - s%alfav(i+k,j) ) 02808 02809 dU_face = 0.5d0* ( dU0( i - 1 + k , j ) + dU0( i - 1 + k , j +1 ) ) 02810 if (qx(k) > 0.) then 02811 ! 02812 faceval_x(k) = dV0( i- 1 + k, j ) * cosd + dU_face * sind 02813 ! 02814 else 02815 ! 02816 faceval_x(k) = dV0( i + k, j ) *cosd + dU_face * sind 02817 ! 02818 endif 02819 ! 02820 02821 cosd = cos( s%alfav(i,j+k) - s%alfav(i,j+k-1) ) 02822 sind = sig(k)*sin( s%alfav(i,j+k) - s%alfav(i,j+k-1) ) 02823 ! 02824 dU_face = 0.5d0* ( dU0( i - 1 , j + k ) + dU0( i , j + k ) ) 02825 if (qy(k) > 0.) then 02826 ! 02827 faceval_y(k) = dV0( i , j -1 + k ) * cosd+ dU_face * sind 02828 ! 02829 else 02830 ! 02831 faceval_y(k) = dV0( i , j + k) * cosd + dU_face * sind 02832 ! 02833 endif 02834 ! 02835 enddo 02836 02837 fac = par%dt / s%hvm(i,j) * s%dsdnvi(i,j) 02838 02839 s%dV(i,j) = dV0(i,j) + fac *dV0(i,j) * ( qx(1) + qy(1) - qx(0) - qy(0) ) & 02840 - fac *( qx(1)*faceval_x(1) + qy(1)*faceval_y(1) & 02841 - qx(0)*faceval_x(0) - qy(0)*faceval_y(0) ) 02842 ! 02843 else 02844 ! 02845 s%dV(i,j) = 0.0 02846 ! 02847 endif 02848 ! 02849 enddo 02850 ! 02851 enddo 02852 ! 02853 endif 02854 02855 ! 02856 ! ------------- Include explicit approximations -------------------- 02857 ! 02858 if (par%secorder == 1) then 02859 ! 02860 ! Use explicit estimate for the nonhydrostatic pressure in U-momentum 02861 ! 02862 do j=jmin_uu,jmax_uu 02863 ! 02864 do i=imin_uu,imax_uu 02865 ! 02866 s%uu(i,j) = s%uu(i,j) + au(1,i,j) * s%pres(i+1,j) + au(0,i,j) * s%pres(i ,j) 02867 ! 02868 enddo 02869 ! 02870 enddo 02871 ! 02872 ! Use explicit estimate for the nonhydrostatic pressure in V-momentum 02873 ! 02874 do j=jmin_vv,jmax_vv 02875 ! 02876 do i=imin_vv,imax_vv 02877 ! 02878 s%vv(i,j) = s%vv(i,j) + av(1,i,j) * s%pres(i,j+1) + av(0,i,j) * s%pres(i ,j) 02879 ! 02880 enddo 02881 ! 02882 enddo 02883 02884 ! 02885 ! Use explicit estimate for the nonhydrostatic pressure in dU-momentum 02886 ! 02887 if ( par%nhlay > 0.0d0 ) then 02888 ! 02889 do j=jmin_uu,jmax_uu 02890 ! 02891 do i=imin_uu,imax_uu 02892 ! 02893 s%dU(i,j) = s%dU(i,j) + adU(1,i,j) * s%pres(i+1,j) + adU(0,i,j) * s%pres(i ,j) 02894 ! 02895 enddo 02896 ! 02897 enddo 02898 ! 02899 ! Use explicit estimate for the nonhydrostatic pressure in dV-momentum 02900 ! 02901 do j=jmin_vv,jmax_vv 02902 ! 02903 do i=imin_vv,imax_vv 02904 ! 02905 s%dV(i,j) = s%dV(i,j) + adV(1,i,j) * s%pres(i,j+1) + adV(0,i,j) * s%pres(i ,j) 02906 ! 02907 enddo 02908 ! 02909 enddo 02910 ! 02911 endif 02912 ! 02913 ! 02914 ! Use explicit estimate for the nonhydrostatic pressure in Wm-momentum 02915 ! 02916 do j=jmin_zs,jmax_zs 02917 ! 02918 do i=imin_zs,imax_zs 02919 ! 02920 Wm(i,j) = Wm(i,j) + aws(1,i,j) * s%pres(i,j) 02921 ! 02922 enddo 02923 ! 02924 enddo 02925 ! 02926 02927 #ifdef USEMPI 02928 ! 02929 call xmpi_shift_ee(Wm) 02930 #endif 02931 ! 02932 call flow_secondorder_advW(s,par,Wm,Wm_old, nonhZ) 02933 02934 #ifdef USEMPI 02935 call xmpi_shift_ee(Wm) 02936 #endif 02937 ! 02938 endif 02939 02940 ! 02941 if (par%nhlay > 0.) then 02942 ! 02943 ! 02944 ! 02945 ! 02946 !----------------- Bottom stress, internal stresses and vertical advection ------------------- 02947 ! 02948 ! Integration of bottom stress and internal stress terms occurs implicity to ensure that they 02949 ! do not introduce any additional stability concerns. 02950 ! 02951 ! Moreover, the vertical exchange of momentum due to advection between the 02952 ! layers also influences dU (but not U), so we have 02953 ! 02954 ! the equation reads dU - dU* 02955 ! -------- = dU ( fric + stress ) + s%uu * fric - omega * U_face 02956 ! dt 02957 ! To ensure that the system is solvable, we approximate U_face by it's upwind value. 02958 ! 02959 do j=jmin_uu,jmax_uu ! Robert: was do j=2,s%ny-1 02960 ! 02961 do i=imin_uu,imax_uu 02962 ! 02963 if ( s%wetu(i,j) == 1 .and. wcoef(i,j) > 0.0d0 ) then 02964 ! 02965 ! 02966 ! Bottom friction 02967 ! 02968 ub = uu0(i,j) + (1.- wcoef(i,j)) * s%dU(i,j) 02969 vb = s%vu(i,j) + .25 * (1.- wcoef(i,j)) * ( dV0(i,j) + dV0(i+1,j) + dV0(i,j-1) + dV0(i+1,j-1) ) 02970 umag = sqrt( ub**2 + vb**2 ) 02971 fric = par%dt * s%cfu(i,j) * umag / s%hu(i,j) 02972 02973 ! 02974 ! Internal stresses 02975 ! 02976 ve = 0.0001 !Vertical eddy viscocity 02977 stress = 2.* par%dt * ve / ( s%hum(i,j)**2 * (par%nhlay - par%nhlay**2) ) 02978 02979 ! 02980 ! Vertical advection 02981 ! 02982 advec = .5d0 * ( omega(i,j) + omega(i+1,j) ) * par%dt / s%hum(i,j) 02983 02984 ! 02985 ! Built rhs/diagonal 02986 ! 02987 fac = 0.5d0*( wcoef(i,j) + wcoef(i+1,j) ) 02988 rhs = s%dU(i,j) - uu0(i,j) * ( fric + advec ) 02989 02990 Amat = 1.0d0 + (1.0d0-fac) * fric + stress & 02991 + (1.0d0 - fac) * max( advec , 0. ) & 02992 - fac * min( advec , 0. ) 02993 02994 ! 02995 ! "Solve" system 02996 ! 02997 s%dU(i,j) = rhs / Amat 02998 ! 02999 endif 03000 ! 03001 enddo 03002 ! 03003 enddo 03004 ! 03005 03006 ! 03007 do j=jmin_vv,jmax_vv 03008 ! 03009 do i=imin_vv,imax_vv 03010 ! 03011 ! Velocity in the lower layer at u-point 03012 ! 03013 if (s%wetv(i,j) == 1 .and. wcoef(i,j) > 0.0d0 ) then 03014 ! 03015 ! Bottom friction 03016 ! 03017 vb = vv0(i,j) + (1.- wcoef(i,j)) * s%dV(i,j) 03018 ub = s%uv(i,j) + .25 * (1.- wcoef(i,j)) * ( dU0(i,j) + dU0(i,j+1) + dU0(i-1,j) + dU0(i-1,j+1) ) 03019 umag = sqrt( ub**2 + vb**2 ) 03020 fric = par%dt * s%cfu(i,j) * umag / s%hv(i,j) 03021 03022 ! 03023 ! Internal stresses 03024 ! 03025 ve = 0.0001 !Vertical eddy viscocity 03026 stress = 2.* par%dt * ve / ( s%hvm(i,j)**2 * (par%nhlay - par%nhlay**2) ) 03027 03028 ! 03029 ! Vertical advection 03030 ! 03031 advec = .5d0 * ( omega(i,j) + omega(i,j+1) ) * par%dt / s%hvm(i,j) 03032 03033 ! 03034 ! Built rhs/diagonal 03035 ! 03036 fac = 0.5d0*( wcoef(i,j) + wcoef(i+1,j) ) 03037 rhs = s%dV(i,j) - vv0(i,j) * ( fric + advec ) 03038 03039 Amat = 1.0d0 + (1.0d0-fac) * fric + stress & 03040 + (1.0d0 - fac) * max( advec , 0. ) & 03041 - fac * min( advec , 0. ) 03042 03043 ! 03044 ! "Solve" system 03045 ! 03046 s%dV(i,j) = rhs / Amat 03047 ! 03048 endif 03049 ! 03050 enddo 03051 ! 03052 enddo 03053 ! 03054 endif 03055 03056 #ifdef USEMPI 03057 !ARRAYS DU and DV need to be communicated 03058 call xmpi_shift_ee(s%dU) 03059 call xmpi_shift_ee(s%dV) 03060 #endif 03061 ! 03062 end subroutine nonh_2lay_pred_3d 03063 03064 subroutine nonh_break( s , par ) 03065 !-------------------------- PURPOSE ---------------------------- 03066 ! 03067 ! Include the pressure explicitly in the predictor step 03068 ! 03069 !-------------------------- DEPENDENCIES ---------------------------- 03070 use spaceparams 03071 use params 03072 use flow_secondorder_module 03073 use xmpi_module 03074 !-------------------------- ARGUMENTS ---------------------------- 03075 03076 type(spacepars) ,intent(inout) :: s 03077 type(parameters),intent(in) :: par 03078 03079 !-------------------------- LOCAL VARIABLES ---------------------------- 03080 03081 !Indices 03082 integer(kind=iKind) :: i,ie 03083 integer(kind=iKind) :: j,js 03084 integer(kind=iKind) :: jmin,jmax,ineighbour 03085 real(kind=rKind) :: wmax,dzdt 03086 03087 03088 if (s%ny>0) then 03089 jmin = 2 03090 jmax = s%ny 03091 else 03092 jmin = 1 03093 jmax = 1 03094 endif 03095 03096 ! 03097 ! Determine if a velocity point will be included in the nonh pressure matrix, do not include if: 03098 ! 03099 ! (1) The point is dry 03100 ! (2) The relative wave length kd of the smallest possible wave (L=2dx) is smaller than kdmin 03101 ! (3) The interpolated waterlevel zs is below the bottom (steep cliffs with overwash situations) 03102 ! (4) Where Miche breaker criterium applies -> bores are hydrostatic 03103 ! max steepness = H/L = maxbrsteep 03104 ! dz/dx = maxbrsteep 03105 ! dz/dx = dz/dt/c = w/c = w/sqrt(gh) 03106 ! wmax = maxbrsteep*sqrt(gh) 03107 03108 ! 03109 do j=1,s%ny+1 03110 ! 03111 do i=1,s%nx+1 03112 ! 03113 ie = min(s%nx,i+1) 03114 ! 03115 if ( (s%wetU(i,j)==1 ) .and. ( 0.5_rKind*( s%zs(i,j) + s%zs(ie,j) ) > zbu(i,j) ) ) then 03116 ! 03117 nonhU(i,j) = 1 03118 ! 03119 else 03120 ! 03121 nonhU(i,j) = 0 03122 ! 03123 endif 03124 ! 03125 enddo 03126 ! 03127 enddo 03128 03129 ! 03130 if (s%ny>2) then 03131 ! 03132 do j=1,s%ny+1 03133 ! 03134 js = min(s%ny,j+1) 03135 ! 03136 do i=1,s%nx+1 03137 ! 03138 if ( (s%wetV(i,j)==1 ) .and. (0.5_rKind*(s%zs(i,j) + s%zs(i,js)) > zbv(i,j) ) ) then 03139 ! 03140 nonhV(i,j) = 1 03141 ! 03142 else 03143 ! 03144 nonhV(i,j) = 0 03145 ! 03146 endif 03147 ! 03148 enddo 03149 ! 03150 enddo 03151 ! 03152 else 03153 ! 03154 nonhV = 0 03155 ! 03156 endif 03157 ! 03158 ! Determine if a velocity point will be included in the nonh pressure matrix, include if 03159 ! any of the surrounding velocity points is included. 03160 ! 03161 03162 ! 03163 do j=jmin_zs,jmax_zs 03164 ! 03165 js = max( j-1 , 1 ) 03166 ! 03167 do i=imin_zs,imax_zs 03168 ! 03169 if (max( nonhV(i,j) , nonhV(i,js) ,nonhU(i,j), nonhU(i-1,j)) > 0) then 03170 ! 03171 nonhZ(i,j) = 1 03172 ! 03173 else 03174 ! 03175 nonhZ(i,j) = 0 03176 ! 03177 endif 03178 ! 03179 enddo 03180 ! 03181 enddo 03182 ! 03183 03184 if (par%nhbreaker == 1) then 03185 ! 03186 ! Breaking is active 03187 ! 03188 do j = jmin_zs,jmax_zs 03189 ! 03190 js = max( j-1 , 1 ) 03191 ! 03192 do i =imin_zs,imax_zs 03193 ! 03194 dzdt = - ( s%uu(i,j) * s%hu(i,j) - s%uu(i-1,j) * s%hu(i-1,j) ) / s%dsz(i,j) 03195 ! 03196 ineighbour = 0 03197 if ( s%breaking(i-1,j)==1 .or. s%breaking(i+1,j) == 1 ) then 03198 ! 03199 ineighbour = 1 03200 ! 03201 endif 03202 03203 if ( s%ny > 1 ) then 03204 ! 03205 dzdt = dzdt - ( s%vv(i,j) * s%hv(i,j) - s%vv(i,j-1) * s%hv(i,j-1) ) / s%dnz(i,j) 03206 03207 if ( s%breaking(i,j-1) == 1 .or. s%breaking(i,j+1) == 1 ) then 03208 ! 03209 ineighbour = 1 03210 ! 03211 endif 03212 ! 03213 endif 03214 ! 03215 03216 03217 03218 wmax = sqrt(par%g*s%hh(i,j)) 03219 03220 if ( s%breaking(i,j) == 0 ) then 03221 ! 03222 if ( dzdt > par%maxbrsteep*wmax ) then 03223 ! 03224 s%breaking(i,j) = 2 03225 ! 03226 elseif (dzdt > par%reformsteep*wmax .and. ineighbour == 1) then 03227 ! 03228 s%breaking(i,j) = 2 03229 ! 03230 endif 03231 ! 03232 elseif ( s%breaking(i,j) == 1 ) then 03233 ! 03234 if ( dzdt < 0.0d0 ) then 03235 ! 03236 s%breaking(i,j) = -1 03237 ! 03238 endif 03239 ! 03240 endif 03241 ! 03242 enddo 03243 ! 03244 enddo 03245 03246 where (s%breaking==2) 03247 s%breaking = 1 03248 endwhere 03249 03250 where (s%breaking==-1) 03251 s%breaking = 0 03252 endwhere 03253 #ifdef USEMPI 03254 !ARRAYS communication 03255 call xmpi_shift_ee(s%breaking) 03256 #endif 03257 ! turn off non-hydrostatic pressure correction in areas with breaking and increase viscosity 03258 where (s%breaking/=0) 03259 nonhZ = 0 03260 s%pres = 0 03261 endwhere 03262 ! 03263 endif 03264 ! 03265 end subroutine nonh_break 03266 03267 ! 03268 !============================================================================== 03269 subroutine nonh_init_wcoef(s,par) 03270 !============================================================================== 03271 ! 03272 03273 ! DATE AUTHOR CHANGES 03274 ! 03275 ! March 2013 Robert McCall Move computation of wcoef to own function so callable separately 03276 03277 !------------------------------------------------------------------------------- 03278 ! DECLARATIONS 03279 !------------------------------------------------------------------------------- 03280 03281 !-------------------------- PURPOSE ---------------------------- 03282 03283 ! Recompute wcoef optimisation 03284 03285 !-------------------------- DEPENDENCIES ---------------------------- 03286 use spaceparams 03287 use params, only: parameters 03288 use wave_functions_module, only: dispersion 03289 03290 !-------------------------- ARGUMENTS ---------------------------- 03291 03292 type(spacepars) ,intent(in) :: s 03293 type(parameters),intent(in) :: par 03294 03295 !-------------------------- LOCAL VARIABLES ---------------------------- 03296 03297 integer(kind=iKind) :: i,j 03298 real (kind=rKind) :: d,sigma,k 03299 03300 03301 ! 03302 !------------------------------------------------------------------------------- 03303 ! IMPLEMENTATION 03304 !------------------------------------------------------------------------------- 03305 ! 03306 03307 if (allocated(wcoef)) then 03308 if ( par%nonhq3d == 1 ) then 03309 wcoef = par%nhlay 03310 else 03311 sigma = 2.0_rkind*par%px/par%Topt 03312 if (par%dispc <= 0.0_rKind) then 03313 !Calculate the optimum value of the dispersion coeficiant. 03314 do j=1,s%ny+1 03315 do i=1,s%nx+1 03316 d = max(s%zs0(i,j)-s%zb(i,j),par%eps) 03317 k = disper(sigma,d,par%g,2.0_rKind*par%px,0.0001_rKind) 03318 if (d>0.0_rKind) then 03319 wcoef(i,j) = 1.0_rKind / ( 4.0_rKind*par%g / (d*sigma**2) - 4.0_rKind / ((k*d)**2) ) 03320 else 03321 wcoef(i,j) = 1.d0 03322 endif 03323 enddo 03324 enddo 03325 else 03326 wcoef = par%dispc 03327 endif 03328 endif 03329 endif 03330 03331 03332 end subroutine nonh_init_wcoef 03333 ! 03334 03335 03336 03337 03338 subroutine zuzv(s) 03339 !============================================================================== 03340 ! 03341 03342 ! DATE AUTHOR CHANGES 03343 ! 03344 ! November 2010 Pieter Bart Smit New Subroutine 03345 03346 03347 !------------------------------------------------------------------------------ 03348 ! DECLARATIONS 03349 !------------------------------------------------------------------------------ 03350 03351 !-------------------------- PURPOSE ---------------------------- 03352 ! 03353 ! Interpolate bottom and free surface location to u/v points 03354 ! 03355 !-------------------------- DEPENDENCIES ---------------------------- 03356 use spaceparams 03357 use xmpi_module 03358 !-------------------------- ARGUMENTS ---------------------------- 03359 03360 type(spacepars) ,intent(inout) :: s 03361 !-------------------------- LOCAL VARIABLES ---------------------------- 03362 03363 integer(kind=iKind) :: i 03364 integer(kind=iKind) :: j 03365 03366 !------------------------------------------------------------------------------- 03367 ! IMPLEMENTATION 03368 !------------------------------------------------------------------------------- 03369 !Bottom location in U point 03370 do j=1,s%ny+1 03371 do i=1,s%nx 03372 zbu(i,j) = max(s%zb(i,j),s%zb(min(s%nx,i)+1,j)) 03373 enddo 03374 enddo 03375 #ifdef USEMPI 03376 if (xmpi_isbot) then 03377 zbu(s%nx+1,:) = s%zb(s%nx+1,:) 03378 endif 03379 #else 03380 zbu(s%nx+1,:) = s%zb(s%nx+1,:) 03381 #endif 03382 !Bottom location in V point 03383 if (s%ny>2) then 03384 do j=1,s%ny 03385 do i=1,s%nx+1 03386 zbv(i,j) = max(s%zb(i,j),s%zb(i,min(s%ny,j)+1)) 03387 enddo 03388 enddo 03389 zbv(:,s%ny+1) = s%zb(:,s%ny+1) 03390 else 03391 zbv = s%zb 03392 endif 03393 03394 !Free surface location in u-point 03395 do j=1,s%ny+1 03396 do i=1,s%nx 03397 ! 03398 zsu(i,j) = s%hu(i,j) + zbu(i,j) 03399 ! 03400 enddo 03401 enddo 03402 03403 !Free surface location in v-point 03404 if (s%ny>2) then 03405 do j=1,s%ny 03406 do i=1,s%nx+1 03407 zsv(i,j) = s%hv(i,j) + zbv(i,j) 03408 enddo 03409 enddo 03410 else 03411 do j=1,s%ny+1 03412 zsv(:,j) = s%zs(:,min(2,s%ny+1)) 03413 enddo 03414 endif 03415 ! 03416 end subroutine zuzv 03417 03418 ! 03419 !============================================================================== 03420 real(kind=rKind) function disper(w,d,g,pi2,accuracy) 03421 !============================================================================== 03422 ! 03423 03424 03425 !------------------------------------------------------------------------------- 03426 ! DECLARATIONS 03427 !------------------------------------------------------------------------------- 03428 03429 !-------------------------- PURPOSE ---------------------------- 03430 03431 ! Calculate k for a given intrinsic frequency w and depth d. First use the fenton 03432 ! approximation and then iterate for better accuracy (when necessary) 03433 03434 03435 !-------------------------- ARGUMENTS ---------------------------- 03436 03437 real(kind=rKind),intent(in) :: w 03438 real(kind=rKind),intent(in) :: d 03439 real(kind=rKind),intent(in) :: g 03440 real(kind=rKind),intent(in) :: pi2 03441 real(kind=rKind),intent(in) :: accuracy 03442 03443 !-------------------------- LOCAL VARIABLES ---------------------------- 03444 03445 real(kind=rKind) :: k 03446 real(kind=rKind) :: alpha 03447 real(kind=rKind) :: L,Ldeep,pi2d,w2 03448 real(kind=rKind) :: error 03449 03450 !------------------------------------------------------------------------------- 03451 ! IMPLEMENTATION 03452 !------------------------------------------------------------------------------- 03453 w2 = w**2 03454 alpha = d*w2/g 03455 k = alpha*(1.0_rKind/sqrt(tanh(alpha)))/d 03456 L = pi2/k 03457 Ldeep = g*pi2/(w2) 03458 pi2d = pi2*d 03459 03460 error = abs(g*k*tanh(k*d)-w2)/w2 03461 do while (error > accuracy) 03462 L = Ldeep*tanh(pi2d/L) 03463 k = pi2/L 03464 error = abs(g*k*tanh(k*d)-w2)/w2 03465 enddo 03466 03467 disper = k 03468 end function disper 03469 03470 03471 end module nonh_module
1.7.6.1
