|
XBeach
|
00001 !============================================================================== 00002 ! FLOW_SECONDORDER_MODULE 00003 !============================================================================== 00004 ! 00005 ! DATE AUTHOR CHANGES 00006 ! 00007 ! october 2009 Pieter Bart Smit New module 00008 ! 00009 !****************************************************************************** 00010 ! INTERFACE 00011 !****************************************************************************** 00012 00013 module flow_secondorder_module 00014 00015 implicit none 00016 save 00017 private 00018 00019 !----------------------------- PARAMETERS ----------------------------------- 00020 00021 include 'nh_pars.inc' !Default precision etc. 00022 00023 !----------------------------- VARIABLES ----------------------------------- 00024 00025 00026 real(kind=rKind),dimension(:,:),allocatable :: wrk1 !Temporary storage for: 00027 ! (i ) du-velocity at i ,j+1/2 00028 ! (ii) dv-velocity at i ,j+1/2 00029 ! (iii) dqx at i+1/2,j 00030 real(kind=rKind),dimension(:,:),allocatable :: wrk2 !Temporary storage for: 00031 ! (i ) du-velocity at i+1/2,j 00032 ! (ii) dv-velocity at i+1/2,j 00033 ! (iii) dqx at i ,j+1/2 00034 00035 logical :: initialized = .false. 00036 00037 public flow_secondorder_advUV 00038 public flow_secondorder_advW 00039 public flow_secondorder_con 00040 public minmod 00041 public flow_secondorder_huhv 00042 ! 00043 !****************************************************************************** 00044 ! SUBROUTINES/FUNCTIONS 00045 !****************************************************************************** 00046 00047 contains 00048 00049 subroutine flow_secondorder_init(s) 00050 ! 00051 00052 !------------------------------------------------------------------------------- 00053 ! DECLARATIONS 00054 !------------------------------------------------------------------------------- 00055 00056 ! 00057 !-------------------------- PURPOSE ---------------------------- 00058 ! 00059 ! Initializes the resources needed for the second order mcCormack scheme 00060 ! 00061 ! 00062 !-------------------------- DEPENDENCIES ---------------------------- 00063 ! 00064 ! 00065 ! -- MODULES -- 00066 use spaceparams 00067 00068 ! 00069 !-------------------------- ARGUMENTS ---------------------------- 00070 ! 00071 type(spacepars) ,intent(inout) :: s 00072 00073 !Allocate resources 00074 allocate ( wrk1(s%nx+1,s%ny+1)) 00075 allocate ( wrk2(s%nx+1,s%ny+1)) 00076 00077 wrk1 = 0.0_rKind 00078 wrk2 = 0.0_rKind 00079 00080 initialized = .true. 00081 end subroutine flow_secondorder_init 00082 00083 ! 00084 !============================================================================== 00085 subroutine flow_secondorder_advUV(s,par,uu_old,vv_old) 00086 !============================================================================== 00087 ! 00088 00089 ! DATE AUTHOR CHANGES 00090 ! 00091 ! November 2010 Pieter Bart Smit New Subroutine 00092 ! November 2014 Pieter Bart Smit Updated for curvilinear and mpi 00093 00094 !------------------------------------------------------------------------------- 00095 ! DECLARATIONS 00096 !------------------------------------------------------------------------------- 00097 00098 ! 00099 !-------------------------- PURPOSE ---------------------------- 00100 ! 00101 ! Calculates second order correction to the advection terms for U and V 00102 ! 00103 00104 ! 00105 !-------------------------- DEPENDENCIES ---------------------------- 00106 ! 00107 use spaceparams 00108 use params, only: parameters 00109 00110 ! 00111 !-------------------------- METHOD ---------------------------- 00112 ! 00113 ! First a prediction for the velocity on the new timelevel is done in flow_timestep 00114 ! using an explicit first order euler step. We use this prediction (s%uu,s%vv), and the velocity 00115 ! at the beginning of the step (uu_old, vv_old), to do a slope limited correction. 00116 ! 00117 00118 ! 00119 !-------------------------- ARGUMENTS ---------------------------- 00120 ! 00121 type(spacepars) ,intent(inout) :: s 00122 type(parameters) ,intent(in) :: par 00123 00124 real(kind=rkind),dimension(s%nx+1,s%ny+1),intent(in) :: vv_old !"Old" v-velocity 00125 real(kind=rkind),dimension(s%nx+1,s%ny+1),intent(in) :: uu_old !"Old" u-velocity 00126 ! 00127 00128 !-------------------------- LOCAL VARIABLES ---------------------------- 00129 ! 00130 integer(kind=iKind) :: iww !i-2 :Two points 'west' of current point i 00131 integer(kind=iKind) :: iw !i-1 :One ... 00132 integer(kind=iKind) :: i !Current point 00133 integer(kind=iKind) :: ie !i+1 00134 integer(kind=iKind) :: iee !i+2 00135 integer(kind=iKind) :: jnn !j-2 00136 integer(kind=iKind) :: jn !j-1 00137 integer(kind=iKind) :: j !Current point 00138 integer(kind=iKind) :: js !j+1 00139 integer(kind=iKind) :: jss !j+2 00140 integer(kind=iKind) :: jmin1d,jmax1d,jmin,jmax,imin,imax ! index for superfast1D 00141 00142 real(kind=rKind) :: mindepth ! Near the dry/wet interface the higher order interpolations 00143 ! can cause unwanted effects. To avoid this any extrapolation/interpolation 00144 ! is disabled when the lowest surface elevation within the molecule is lower then 00145 ! the highers bottom elevation. 00146 00147 real(kind=rKind) :: delta1 ! "Central" difference 00148 real(kind=rKind) :: delta2 ! "Upwind" difference 00149 00150 real(kind=rKind) :: qx,qy !discharge south 00151 real(kind=rKind) :: fac !A factor 00152 00153 !------------------------------------------------------------------------------- 00154 ! IMPLEMENTATION 00155 !------------------------------------------------------------------------------- 00156 00157 !Initialize/allocate arrays on first entry 00158 if (.not. initialized) then 00159 call flow_secondorder_init(s) 00160 endif 00161 00162 ! 00163 imin = 2 00164 if (xmpi_istop) then 00165 ! 00166 imin = 3 00167 ! 00168 endif 00169 00170 imax = s%nx 00171 if (xmpi_isbot) then 00172 ! 00173 imax = s%nx-1 00174 ! 00175 endif 00176 00177 jmin = 2 00178 if (xmpi_isleft) then 00179 ! 00180 jmin = 3 00181 ! 00182 endif 00183 00184 jmax = s%ny 00185 if (xmpi_isright) then 00186 ! 00187 jmax = s%ny-1 00188 ! 00189 endif 00190 00191 ! 00192 if (s%ny>0) then 00193 ! 00194 jmin1d = 2 00195 jmax1d = s%ny 00196 ! 00197 else 00198 ! 00199 jmin1d = 1 00200 jmax1d = 1 00201 ! 00202 endif 00203 00204 ! 00205 !-- calculate qx * du in z-points -- 00206 ! 00207 do j=jmin1d,jmax1d 00208 ! 00209 do i=2,s%nx 00210 ! 00211 wrk1(i,j) = 0.0_rKind 00212 ie = min(i+1,imax+1) 00213 iee = min(i+2,imax+1) 00214 iw = max(i-1,1) 00215 iww = max(i-2,1) 00216 mindepth = minval(s%zs(iww:iee,j))-maxval(s%zb(iww:iee,j)) 00217 00218 qx = .5 * ( s%qx(i,j) + s%qx(iw,j) ) * cos( s%alfau( i ,j ) - s%alfau( iw,j) ) 00219 ! 00220 if (mindepth > par%eps) then 00221 ! 00222 if ( qx > 0.d0 .and. i>imin ) then 00223 ! 00224 delta1 = (s%uu(i,j) -uu_old(iw ,j)) / s%dsz(i,j) 00225 delta2 = (s%uu(iw,j)-uu_old(iww,j)) / s%dsz(iw,j) 00226 wrk1(i,j) = 0.5d0*s%dsu(iw,j) * minmod(delta1,delta2) * qx 00227 ! 00228 elseif ( qx < 0.d0 .and. i < imax + 1 ) then 00229 ! 00230 delta1 = (uu_old(i,j) -s%uu(iw ,j)) / s%dsz(i,j) 00231 delta2 = (uu_old(ie,j)-s%uu(i,j)) / s%dsz(ie,j) 00232 wrk1(i,j) = -0.5d0*s%dsu(i,j) * minmod(delta1,delta2) *qx 00233 ! 00234 endif 00235 ! 00236 endif 00237 ! 00238 enddo 00239 ! 00240 enddo 00241 00242 ! 00243 !-- calculate qy * du in c-points -- 00244 ! 00245 if ( s%ny > 0 ) then 00246 ! 00247 do j=2,s%ny 00248 ! 00249 do i=2,s%nx-1 00250 ! 00251 wrk2(i,j) = 0.0_rkind 00252 js = min(j+1,jmax+1) 00253 jss = min(j+2,jmax+1) 00254 jn = max(j-1,1) 00255 mindepth = minval(s%zs(i:i+1,jn:jss))-maxval(s%zb(i:i+1,jn:jss)) 00256 ! 00257 if ((mindepth > par%eps)) then 00258 ! 00259 qy = ( s%qy(i+1,j) + s%qy(i,j) ) * .5d0 * cos( s%alfau( i , j ) - s%alfau( i , j + 1 ) ) 00260 ! 00261 if ( qy > 0.d0 .and. j > jmin - 1 ) then 00262 ! 00263 delta1 = (s%uu(i,js) -uu_old(i ,j )) / s%dnv(i,j) 00264 delta2 = (s%uu(i,j) -uu_old(i ,jn)) / s%dnv(i,jn) 00265 wrk2(i,j) = 0.5d0*s%dnv(i,j)*minmod(delta1,delta2) * qy 00266 ! 00267 elseif ( qy < 0.d0 .and. j < jmax ) then 00268 ! 00269 delta1 = (uu_old(i,js) -s%uu(i ,j)) / s%dnv(i,j) 00270 delta2 = (uu_old(i,jss)-s%uu(i,js)) / s%dnv(i,js) 00271 wrk2(i,j) = -0.5d0*s%dnv(i,j)*minmod(delta1,delta2) * qy 00272 ! 00273 endif 00274 endif 00275 ! 00276 enddo 00277 ! 00278 enddo 00279 ! 00280 wrk2(:,1 ) = 0.0_rKind 00281 ! 00282 endif 00283 00284 !CORRECTION TO U 00285 if (s%ny>0) then 00286 ! 00287 do j=jmin_uu,jmax_uu 00288 ! 00289 do i=imin_uu,imax_uu 00290 ! 00291 fac = par%dt / s%hum(i,j) * s%dsdnui(i,j) 00292 s%uu(i,j) = s%uu(i,j) - fac * ( s%dnz(i+1,j) * wrk1(i+1,j) - s%dnz(i, j) * wrk1(i ,j ) & 00293 + s%dsc(i,j ) * wrk2(i ,j) - s%dsc(i,j-1) * wrk2(i ,j-1) ) 00294 ! 00295 enddo 00296 ! 00297 enddo 00298 ! 00299 else 00300 ! 00301 do i=imin_uu,imax_uu 00302 ! 00303 s%uu(i,1) = s%uu(i,1)-par%dt/s%hum(i,1)*( ( wrk1(i+1,1) - wrk1(i ,1 ) ) / s%dsu(i,1) ) 00304 ! 00305 enddo 00306 ! 00307 endif 00308 00309 !== SECOND ORDER EXPLICIT CORRECTION TO V == 00310 if (s%ny>2) then 00311 ! 00312 do j=2,s%ny 00313 ! 00314 do i=2,s%nx 00315 ! 00316 wrk1(i,j) = 0. 00317 js = min(j+1,jmax+1) 00318 jss = min(j+2,jmax+1) 00319 jn = max(j-1,1) 00320 jnn = max(j-2,1) 00321 mindepth = minval(s%zs(i,jnn:jss))-maxval(s%zb(i,jnn:jss)) 00322 ! 00323 qy = 0.5d0 * ( s%qy(i,j)+s%qy(i,jn) ) * cos( s%alfav(i,j) - s%alfav(i,jn) ) * s%dsz(i,j) 00324 ! 00325 if (mindepth > par%eps) then 00326 ! 00327 if (( qy > 0.0) .and. (j>jmin)) then 00328 ! 00329 delta1 = (s%vv(i,j ) - vv_old(i ,jn )) / s%dnz(i,j) 00330 delta2 = (s%vv(i,jn) - vv_old(i ,jnn)) / s%dnz(i,jn) 00331 wrk1(i,j) = 0.5d0*s%dnv(i,jn)*minmod(delta1,delta2)*qy 00332 ! 00333 elseif ( qy < -0.0d0 .and. j<jmax+1) then 00334 ! 00335 delta1 = (vv_old(i,j ) - s%vv(i,jn)) / s%dnz(i,j) 00336 delta2 = (vv_old(i,js) - s%vv(i,j )) / s%dnz(i,js) 00337 wrk1(i,j) = -0.5d0*s%dnv(i,j)*minmod(delta1,delta2)*qy 00338 ! 00339 endif 00340 ! 00341 endif 00342 enddo 00343 enddo 00344 00345 do j=2,s%ny-1 00346 ! 00347 do i=2,s%nx 00348 ! 00349 wrk2(i,j) = 0.0d0 00350 ie = min(i+1,s%nx) 00351 iee = min(i+2,s%nx) 00352 iw = max(i-1,1) 00353 mindepth = minval(s%zs(iw:iee,j:j+1))-maxval(s%zb(iw:iee,j:j+1)) 00354 ! 00355 if (mindepth > par%eps) then 00356 ! 00357 qx = .5d0 * ( s%qx(i,j+1) + s%qx(i,j) ) * cos( s%alfav( i , j ) - s%alfav( i + 1 , j ) ) * s%dsc(i,j) 00358 ! 00359 if ( qx > par%Umin .and. i > imin - 1) then 00360 delta1 = (s%vv(ie,j) - vv_old(i ,j )) / s%dsu(i,1) 00361 delta2 = (s%vv(i ,j) - vv_old(iw,j)) / s%dsu(iw,1) 00362 wrk2(i,j) = 0.5d0*s%dsu(i,1)*minmod(delta1,delta2)*qx 00363 elseif ( qx < -par%Umin .and. i < imax ) then 00364 delta1 = (vv_old(ie,j) -s%vv(i ,j)) / s%dsu(ie,1) 00365 delta2 = (vv_old(iee,j)-s%vv(ie,j)) / s%dsu(ie,1) 00366 wrk2(i,j) = -0.5d0*s%dsu(i,1)*minmod(delta1,delta2)*qx 00367 endif 00368 endif 00369 enddo 00370 enddo 00371 wrk2(1,:) = 0.0_rKind 00372 00373 !CORRECTION TO V 00374 do j=jmin_vv,jmax_vv 00375 ! 00376 do i=imin_vv,imax_vv 00377 ! 00378 fac = par%dt / s%hvm(i,j) * s%dsdnvi(i,j) 00379 s%vv(i,j) = s%vv(i,j) - fac*( wrk2(i,j) - wrk2(i-1,j) + wrk1(i,j+1) - wrk1(i ,j) ) 00380 ! 00381 enddo 00382 ! 00383 enddo 00384 endif 00385 ! 00386 end subroutine flow_secondorder_advUV 00387 00388 ! 00389 !============================================================================== 00390 subroutine flow_secondorder_advW(s,par,w,w_old,mskZ) 00391 !============================================================================== 00392 ! 00393 00394 ! DATE AUTHOR CHANGES 00395 ! 00396 ! November 2010 Pieter Bart Smit New Subroutine 00397 ! November 2014 Pieter Bart Smit Modified for nonhq3d 00398 00399 !------------------------------------------------------------------------------- 00400 ! DECLARATIONS 00401 !------------------------------------------------------------------------------- 00402 00403 ! 00404 !-------------------------- PURPOSE ---------------------------- 00405 ! 00406 ! Calculates second order correction to the advection terms for W. Only used 00407 ! in combination WITH the non-hydrostatic module 00408 00409 !-------------------------- DEPENDENCIES ---------------------------- 00410 ! 00411 00412 ! -- MODULES -- 00413 use spaceparams 00414 use params, only: parameters 00415 use xmpi_module 00416 00417 00418 !-------------------------- ARGUMENTS ---------------------------- 00419 ! 00420 00421 type(spacepars) ,intent(inout) :: s 00422 type(parameters) ,intent(in) :: par 00423 00424 real(kind=rKind),dimension(s%nx+1,s%ny+1),intent(inout) :: w !The velocity at the star level 00425 real(kind=rKind),dimension(s%nx+1,s%ny+1),intent(in) :: w_old !The vertical velocity at the n-level 00426 integer(kind=iKind),dimension(s%nx+1,s%ny+1),intent(in) :: mskZ !A mask that determines whether or not the w-equation is active 00427 00428 ! 00429 !-------------------------- LOCAL VARIABLES ---------------------------- 00430 ! 00431 integer(kind=iKind) :: iw !i-1 :One ... 00432 integer(kind=iKind) :: i !Current point 00433 integer(kind=iKind) :: ie !i+1 00434 integer(kind=iKind) :: iee !i+2 00435 integer(kind=iKind) :: jn !j+1 00436 integer(kind=iKind) :: jnn !j+2 00437 integer(kind=iKind) :: j !Current point 00438 integer(kind=iKind) :: js !j-1 00439 integer(kind=iKind) :: jmin,jmax,imin,imax 00440 00441 real(kind=rKind) :: rfac 00442 real(kind=rKind) :: mindepth 00443 real(kind=rKind) :: delta1 00444 real(kind=rKind) :: delta2,qx,qy 00445 00446 !------------------------------------------------------------------------------- 00447 ! IMPLEMENTATION 00448 !------------------------------------------------------------------------------- 00449 00450 00451 ! 00452 imin = 1 00453 if (xmpi_istop) then 00454 ! 00455 imin = 2 00456 ! 00457 endif 00458 00459 imax = s%nx 00460 if (xmpi_isbot) then 00461 ! 00462 imax = s%nx-1 00463 ! 00464 endif 00465 00466 jmin = 1 00467 if (xmpi_isleft) then 00468 ! 00469 jmin = 2 00470 ! 00471 endif 00472 00473 jmax = s%ny 00474 if (xmpi_isright) then 00475 ! 00476 jmax = s%ny-1 00477 ! 00478 endif 00479 ! 00480 00481 !Initialize/allocate arrays on first entry 00482 if (.not. initialized) then 00483 ! 00484 call flow_secondorder_init(s) 00485 ! 00486 endif 00487 00488 if (s%ny > 0) then 00489 ! 00490 ! 2D Codepath 00491 ! 00492 00493 ! 00494 ! Calculate the updated fluxes of w-mom. in the s-dir 00495 ! 00496 wrk1 = 0.0d0 00497 do j=2,s%ny 00498 ! 00499 do i=1,s%nx 00500 ! 00501 if ( mskZ(i,j) == 0 ) cycle 00502 00503 wrk1(i,j) = 0.0_rKind 00504 ! 00505 ie = min(i+1,imax+1) 00506 iee = min(i+2,imax+1) 00507 iw = max(i-1,1) 00508 00509 rfac = real( mskZ( ie , j ) * mskZ( iee , j ) * mskZ( iw , j ), kind=rKind ) 00510 00511 mindepth = minval(s%zs(iw:iee,j))-maxval(s%zb(iw:iee,j)) 00512 qx = s%qx(i,j) * s%dnu( i,j ) * rfac 00513 ! 00514 if (mindepth > par%eps) then 00515 ! 00516 if ( qx > 0.0_rKind .and. i>imin ) then 00517 ! 00518 delta1 = (w(ie,j ) - w_old(i ,j )) / s%dsu(i,j) 00519 delta2 = (w(i,j ) - w_old(iw ,j )) / s%dsu(iw,j) 00520 wrk1(i,j) = 0.5d0*s%dsu(i,j)*minmod(delta1,delta2) * qx 00521 ! 00522 elseif (qx < 0.0_rKind .and. i<imax) then 00523 ! 00524 delta1 = (w_old(ie ,j) - w(i ,j )) / s%dsu(i,j) 00525 delta2 = (w_old(iee,j ) - w(ie,j )) / s%dsu(ie,j) 00526 wrk1(i,j) = -0.5d0*s%dsu(i,j)*minmod(delta1,delta2)*qx 00527 ! 00528 endif 00529 ! 00530 endif 00531 ! 00532 enddo 00533 ! 00534 enddo 00535 00536 ! 00537 ! Calculate the updated fluxes of w-mom. in the n-dir 00538 ! 00539 wrk2 = 0. 00540 do j=1,s%ny 00541 ! 00542 do i=2,s%nx 00543 ! 00544 if ( mskZ(i,j) == 0 ) cycle 00545 00546 wrk2(i,j) = 0.0_rKind 00547 ! 00548 jn = min(j+1,jmax+1) 00549 jnn = min(j+2,jmax+1) 00550 js = max(j-1,1) 00551 00552 rfac = real( mskZ( i , jn ) * mskZ( i , jnn ) * mskZ( i , js ), kind=rKind ) 00553 00554 mindepth = minval(s%zs(i,js:jnn))-maxval(s%zb(i,js:jnn)) 00555 qy = s%qy(i,j) * s%dsv( i,j ) * rfac 00556 ! 00557 if (mindepth > par%eps) then 00558 ! 00559 if ( qy > 0.0_rKind .and. j>jmin ) then 00560 ! 00561 delta1 = (w(i,jn ) - w_old(i ,j )) / s%dnv(i,j) 00562 delta2 = (w(i,j ) - w_old(i ,js )) / s%dnv(i,js) 00563 wrk2(i,j) = 0.5d0*s%dnv(i,j)*minmod(delta1,delta2) * qy 00564 ! 00565 elseif (qy < 0.0_rKind .and. j<jmax) then 00566 ! 00567 delta1 = (w_old(i ,jn) - w(i ,j )) / s%dnv(i,j) 00568 delta2 = (w_old(i,jnn ) - w(i,jn )) / s%dnv(i,jn) 00569 wrk2(i,j) = -0.5d0*s%dnv(i,j)*minmod(delta1,delta2)*qy 00570 ! 00571 endif 00572 ! 00573 endif 00574 ! 00575 enddo 00576 ! 00577 enddo 00578 ! 00579 00580 ! 00581 if (xmpi_istop) then 00582 ! 00583 wrk1(1,:) = 0.0d0 00584 ! 00585 endif 00586 if (xmpi_isbot) then 00587 ! 00588 wrk1(s%nx,:) = 0.0d0 00589 ! 00590 endif 00591 if (xmpi_isleft) then 00592 ! 00593 wrk2(:,1) = 0.0d0 00594 ! 00595 endif 00596 if (xmpi_isright) then 00597 ! 00598 wrk2(:,s%ny) =0. 00599 ! 00600 endif 00601 ! 00602 00603 ! 00604 ! Update w-mom 00605 ! 00606 do j=jmin_zs,jmax_zs 00607 ! 00608 do i=imin_zs,imax_zs 00609 ! 00610 if ( mskZ(i,j) == 0 ) cycle 00611 ! 00612 w(i,j) = w(i,j) - par%dt * s%dsdnzi(i,j) *( wrk1(i,j) + wrk2(i,j) - wrk1(i-1,j) - wrk2(i,j-1) ) /s%hh(i,j) 00613 ! 00614 enddo 00615 ! 00616 enddo 00617 ! 00618 else 00619 ! 00620 ! 1D Codepath 00621 ! 00622 j=1 00623 do i=2,s%nx 00624 ! 00625 if ( mskZ(i,j) == 0 ) cycle 00626 00627 wrk1(i,j) = 0.0_rKind 00628 ! 00629 ie = min(i+1,s%nx) 00630 iee = min(i+2,s%nx) 00631 iw = max(i-1,1) 00632 00633 rfac = real( mskZ( ie , j ) * mskZ( iee , j ) * mskZ( iw , j ), kind=rKind ) 00634 00635 mindepth = minval(s%zs(iw:iee,j))-maxval(s%zb(iw:iee,j)) 00636 qx = s%qx(i,j) * rfac 00637 ! 00638 if (mindepth > par%eps) then 00639 ! 00640 if ( qx > 0.0_rKind .and. i>2 ) then 00641 ! 00642 delta1 = (w(ie,j ) - w_old(i ,j )) / s%dsu(i,j) 00643 delta2 = (w(i,j ) - w_old(iw ,j )) / s%dsu(iw,j) 00644 wrk1(i,j) = 0.5d0*s%dsu(i,j)*minmod(delta1,delta2) * qx 00645 ! 00646 elseif (qx < 0.0_rKind .and. i<s%nx-1) then 00647 ! 00648 delta1 = (w_old(ie ,j) - w(i ,j )) / s%dsu(i,j) 00649 delta2 = (w_old(iee,j ) - w(ie,j )) / s%dsu(ie,j) 00650 wrk1(i,j) = -0.5d0*s%dsu(i,j)*minmod(delta1,delta2)*qx 00651 ! 00652 endif 00653 ! 00654 endif 00655 ! 00656 enddo 00657 00658 ! 00659 if (xmpi_istop) then 00660 ! 00661 wrk1(1,:) = 0.0d0 00662 ! 00663 endif 00664 if (xmpi_isbot) then 00665 ! 00666 wrk1(s%nx,:) = 0.0d0 00667 ! 00668 endif 00669 ! 00670 ! Update w-mom to ** level 00671 ! 00672 do i=imin_zs,imax_zs 00673 ! 00674 if ( mskZ(i,1) == 0 ) cycle 00675 ! 00676 w(i,1) = w(i,1)-par%dt/s%hh(i,1)* ( wrk1(i,1)- wrk1(i-1,1) )/ s%dsz(i,1) 00677 ! 00678 enddo 00679 ! 00680 endif 00681 ! 00682 end subroutine flow_secondorder_advW 00683 00684 ! 00685 !============================================================================== 00686 subroutine flow_secondorder_con(s,par,zs_old) 00687 !============================================================================== 00688 ! 00689 ! DATE AUTHOR CHANGES 00690 ! 00691 ! November 2010 Pieter Bart Smit New Subroutine 00692 00693 !------------------------------------------------------------------------------- 00694 ! DECLARATIONS 00695 !------------------------------------------------------------------------------- 00696 00697 ! 00698 !-------------------------- PURPOSE ---------------------------- 00699 ! 00700 ! Calculates second order correction to the continuity terms 00701 ! 00702 00703 !-------------------------- DEPENDENCIES ---------------------------- 00704 ! 00705 use spaceparams 00706 use params, only: parameters 00707 00708 !-------------------------- ARGUMENTS ---------------------------- 00709 ! 00710 type(spacepars) ,intent(inout) :: s 00711 type(parameters) ,intent(in) :: par 00712 00713 real(kind=rKind),dimension(s%nx+1,s%ny+1),intent(in) :: zs_old 00714 ! 00715 00716 !-------------------------- LOCAL VARIABLES ---------------------------- 00717 ! 00718 00719 integer(kind=iKind) :: iw !i-1 :One ... 00720 integer(kind=iKind) :: i !Current point 00721 integer(kind=iKind) :: ie !i+1 00722 integer(kind=iKind) :: iee !i+2 00723 integer(kind=iKind) :: jn !j-1 00724 integer(kind=iKind) :: j !Current point 00725 integer(kind=iKind) :: js !j+1 00726 integer(kind=iKind) :: jss !j+2 00727 integer(kind=iKind) :: jmin,jmax !index for superfast1D 00728 00729 real(kind=rKind) :: mindepth 00730 real(kind=rKind) :: delta1 00731 real(kind=rKind) :: delta2 00732 !------------------------------------------------------------------------------- 00733 ! IMPLEMENTATION 00734 !------------------------------------------------------------------------------- 00735 00736 if (s%ny>0) then 00737 jmin = 2 00738 jmax = s%ny 00739 else 00740 jmin = 1 00741 jmax = 1 00742 endif 00743 00744 !== SECOND ORDER EXPLICIT CORRECTION TO CONTINUITY == 00745 !Initialize/allocate arrays on first entry 00746 if (.not. initialized) then 00747 call flow_secondorder_init(s) 00748 endif 00749 !return 00750 !correction to mass flux qx 00751 do j=jmin,jmax 00752 do i=2,s%nx-1 00753 wrk1(i,j) = 0.0_rkind 00754 ie = min(i+1,s%nx) 00755 iee = min(i+2,s%nx) 00756 iw = max(i-1,2) 00757 mindepth = minval(s%zs(iw:iee,j))-maxval(s%zb(iw:iee,j)) 00758 if (mindepth>par%eps) then 00759 if (s%uu(i,j) > par%umin .and. i>2 ) then 00760 delta1 = (s%zs(ie,j)- zs_old(i,j ))/ s%dsu(i,1) 00761 delta2 = (s%zs(i,j) - zs_old(iw,j))/ s%dsu(i-1,1) 00762 wrk1(i,j) = s%uu(i,j)*0.5d0*s%dsu(i,1)*minmod(delta1,delta2) 00763 elseif (s%uu(i,j) < -par%umin .and. i<s%nx-1) then 00764 delta1 = (zs_old(iee,j) - s%zs(ie,j)) / s%dsu(i+1,1) 00765 delta2 = (zs_old(ie ,j) - s%zs(i ,j)) / s%dsu(i,1) 00766 wrk1(i,j) = - s%uu(i,j)*0.5d0*s%dsu(i,1)*minmod(delta1,delta2) 00767 endif 00768 endif 00769 enddo 00770 enddo 00771 wrk1(1 ,:) = 0.0_rkind 00772 wrk1(s%nx,:) = 0.0_rkind 00773 00774 !Correction to mass flux qy 00775 if (s%ny > 2) then 00776 do j=2,s%ny-1 00777 do i=2,s%nx 00778 wrk2(i,j) = 0.0_rKind 00779 js = min(j+1,s%ny) 00780 jss = min(j+2,s%ny) 00781 jn = max(j-1,2) 00782 mindepth = minval(s%zs(i,jn:jss))-maxval(s%zb(i,jn:jss)) 00783 if (mindepth> par%eps) then 00784 if (s%vv(i,j) > par%Umin .and. j>2) then 00785 delta1 = (s%zs(i,js) - zs_old(i,j ))/ s%dnv(1,j) 00786 delta2 = (s%zs(i,j) - zs_old(i,jn))/ s%dnv(1,j-1) 00787 wrk2(i,j) = s%vv(i,j)*0.5d0*s%dnv(1,j)*minmod(delta1,delta2) 00788 elseif (s%vv(i,j) < -par%Umin .and. j<s%ny-1) then 00789 delta1 = (zs_old(i,jss) - s%zs(i,js)) / s%dnv(1,j+1) 00790 delta2 = (zs_old(i ,js) - s%zs(i ,j)) / s%dnv(1,j) 00791 wrk2(i,j) = s%vv(i,j)*0.5d0*s%dnv(1,j)*minmod(delta1,delta2) 00792 endif 00793 endif 00794 enddo 00795 enddo 00796 wrk2(:,1 ) = 0.0_rKind 00797 wrk2(:,s%ny) = 0.0_rKind 00798 else 00799 wrk2 = 0.0_rKind 00800 endif 00801 00802 !Update waterlevels 00803 if (s%ny>0) then 00804 do j=jmin_zs,jmax_zs 00805 do i=imin_zs,imax_zs 00806 s%zs(i,j) = s%zs(i,j)-par%dt*( (wrk1(i,j)-wrk1(i-1,j))/ s%dsz(i,1) & 00807 + (wrk2(i,j)-wrk2(i,j-1))/ s%dnz(1,j) ) 00808 enddo 00809 enddo 00810 else 00811 do i=imin_zs,imax_zs 00812 s%zs(i,1) = s%zs(i,1)-par%dt*( (wrk1(i,1)-wrk1(i-1,1))/ s%dsz(i,1) ) 00813 enddo 00814 endif 00815 00816 !Update fluxes 00817 if (s%ny>0) then 00818 s%qx(imin_uu:imax_uu,jmin_uu:jmax_uu) = s%qx(imin_uu:imax_uu,jmin_uu:jmax_uu) +wrk1(imin_uu:imax_uu,jmin_uu:jmax_uu) 00819 s%qy(imin_vv:imax_vv,jmin_vv:jmax_vv) = s%qy(imin_vv:imax_vv,jmin_vv:jmax_vv) +wrk2(imin_vv:imax_vv,jmin_vv:jmax_vv) 00820 else 00821 s%qx(imin_uu:imax_uu,jmin_uu:jmax_uu) = s%qx(imin_uu:imax_uu,jmin_uu:jmax_uu) +wrk1(imin_uu:imax_uu,jmin_uu:jmax_uu) 00822 endif 00823 00824 end subroutine flow_secondorder_con 00825 00826 ! 00827 !============================================================================== 00828 subroutine flow_secondorder_huhv(s,par) 00829 !============================================================================== 00830 ! 00831 ! DATE AUTHOR CHANGES 00832 ! 00833 ! November 2014 Pieter Bart Smit New Subroutine 00834 00835 !------------------------------------------------------------------------------- 00836 ! DECLARATIONS 00837 !------------------------------------------------------------------------------- 00838 00839 ! 00840 !-------------------------- PURPOSE ---------------------------- 00841 ! 00842 ! Calculates second order correction to the waterlevels 00843 ! 00844 00845 !-------------------------- DEPENDENCIES ---------------------------- 00846 ! 00847 use spaceparams 00848 use params, only: parameters 00849 use xmpi_module 00850 00851 !-------------------------- ARGUMENTS ---------------------------- 00852 ! 00853 type(spacepars) ,intent(inout) :: s 00854 type(parameters) ,intent(in) :: par 00855 00856 ! real(kind=rKind),dimension(s%nx+1,s%ny+1),intent(in) :: zs_old 00857 ! 00858 00859 !-------------------------- LOCAL VARIABLES ---------------------------- 00860 ! 00861 00862 integer(kind=iKind) :: iw !i-1 :One ... 00863 integer(kind=iKind) :: i !Current point 00864 integer(kind=iKind) :: ie !i+1 00865 integer(kind=iKind) :: iee !i+2 00866 integer(kind=iKind) :: j !Current point 00867 integer(kind=iKind) :: jmin,jmax,imin,imax 00868 real(kind=rKind) :: mindepth 00869 real(kind=rKind) :: delta1 00870 real(kind=rKind) :: delta2 00871 !------------------------------------------------------------------------------- 00872 ! IMPLEMENTATION 00873 !------------------------------------------------------------------------------- 00874 00875 ! 00876 imin = 1 00877 if (xmpi_istop) then 00878 ! 00879 imin = 2 00880 ! 00881 endif 00882 00883 imax = s%nx 00884 if (xmpi_isbot) then 00885 ! 00886 imax = s%nx-1 00887 ! 00888 endif 00889 00890 jmin = 1 00891 if (xmpi_isleft) then 00892 ! 00893 jmin = 2 00894 ! 00895 endif 00896 00897 jmax = s%ny 00898 if (xmpi_isright) then 00899 ! 00900 jmax = s%ny-1 00901 ! 00902 endif 00903 00904 if (s%ny > 0) then 00905 ! 00906 ! 2DH Code 00907 ! 00908 do j= jmin_uu , jmax_uu 00909 ! 00910 do i = imin_uu, imax_uu ! Robert: was 1 , s%nx 00911 ! 00912 ie = min(i+1,imax+1) 00913 iee = min(i+2,imax+1) 00914 iw = max(i-1,1) 00915 ! 00916 mindepth = minval(s%zs(iw:iee,j))-maxval(s%zb(iw:iee,j)) 00917 ! 00918 if (mindepth>par%eps) then 00919 ! 00920 if (s%uu(i,j) > par%umin .and. i > imin ) then 00921 ! 00922 delta1 = (s%zs(ie,j)- s%zs(i,j ))/ s%dsu(i,1) 00923 delta2 = (s%zs(i,j) - s%zs(iw,j))/ s%dsu(i-1,1) 00924 s%hu(i,j) = s%hu(i,j) + 0.5d0*s%dsu(i,1)*minmod(delta1,delta2) 00925 ! 00926 elseif (s%uu(i,j) < -par%umin .and. i < imax ) then 00927 ! 00928 delta1 = (s%zs(iee,j) - s%zs(ie,j)) / s%dsu(i+1,1) 00929 delta2 = (s%zs(ie ,j) - s%zs(i ,j)) / s%dsu(i,1) 00930 s%hu(i,j) = s%hu(i,j) - 0.5d0*s%dsu(i,1)*minmod(delta1,delta2) 00931 ! 00932 endif 00933 ! 00934 endif 00935 ! 00936 enddo 00937 ! 00938 enddo 00939 00940 do j= jmin_vv,jmax_vv ! Robert: was 1,s%ny 00941 ! 00942 ie = min(j+1,jmax+1) 00943 iee = min(j+2,jmax+1) 00944 iw = max(j-1,1) 00945 ! 00946 do i=imin_vv,imax_vv 00947 ! 00948 mindepth = minval(s%zs(i,iw:iee))-maxval(s%zb(i,iw:iee)) 00949 ! 00950 if (mindepth>par%eps) then 00951 ! 00952 if (s%vv(i,j) > par%umin .and. j>jmin ) then 00953 ! 00954 delta1 = (s%zs(i,ie)- s%zs(i,j ))/ s%dnv(i,j) 00955 delta2 = (s%zs(i,j) - s%zs(i,iw))/ s%dnv(i,j-1) 00956 s%hv(i,j) = s%hv(i,j) + 0.5d0*s%dnv(i,j)*minmod(delta1,delta2) 00957 ! 00958 elseif (s%vv(i,j) < -par%umin .and. j<jmax) then 00959 ! 00960 delta1 = (s%zs(i,iee) - s%zs(i,ie)) / s%dnv(i,j+1) 00961 delta2 = (s%zs(i ,ie) - s%zs(i ,j)) / s%dnv(i,j) 00962 s%hv(i,j) = s%hv(i,j) - 0.5d0*s%dnv(i,j)*minmod(delta1,delta2) 00963 ! 00964 endif 00965 ! 00966 endif 00967 ! 00968 enddo 00969 enddo 00970 else 00971 ! 00972 ! 1DH Code 00973 ! 00974 j= 1 00975 do i=imin_uu, imax_uu ! Robert: was 1,s%nx-1 00976 ! 00977 ie = min(i+1,imax+1) 00978 iee = min(i+2,imax+1) 00979 iw = max(i-1,1) 00980 mindepth = minval(s%zs(iw:iee,j))-maxval(s%zb(iw:iee,j)) 00981 ! 00982 if (mindepth>par%eps) then 00983 ! 00984 if (s%uu(i,j) > par%umin .and. i>imin ) then 00985 ! 00986 delta1 = (s%zs(ie,j)- s%zs(i,j ))/ s%dsu(i,1) 00987 delta2 = (s%zs(i,j) - s%zs(iw,j))/ s%dsu(i-1,1) 00988 s%hu(i,j) = s%hu(i,j) + 0.5d0*s%dsu(i,1)*minmod(delta1,delta2) 00989 ! 00990 elseif (s%uu(i,j) < -par%umin .and. i<imax) then 00991 ! 00992 delta1 = (s%zs(iee,j) - s%zs(ie,j)) / s%dsu(i+1,1) 00993 delta2 = (s%zs(ie ,j) - s%zs(i ,j)) / s%dsu(i,1) 00994 s%hu(i,j) = s%hu(i,j) - 0.5d0*s%dsu(i,1)*minmod(delta1,delta2) 00995 endif 00996 ! 00997 endif 00998 ! 00999 enddo 01000 ! 01001 endif 01002 01003 end subroutine flow_secondorder_huhv 01004 ! 01005 !============================================================================== 01006 real(kind=rKind) pure function minmod(delta1,delta2) 01007 !============================================================================== 01008 ! 01009 ! DATE AUTHOR CHANGES 01010 ! 01011 ! November 2010 Pieter Bart Smit New Subroutine 01012 ! 01013 !------------------------------------------------------------------------------- 01014 ! DECLARATIONS 01015 !------------------------------------------------------------------------------- 01016 01017 ! 01018 !-------------------------- PURPOSE ---------------------------- 01019 ! 01020 ! Limits succesive gradients using the TVD minmod limiter. 01021 ! 01022 !-------------------------- METHOD ---------------------------- 01023 ! 01024 ! - MINMOD LIMITER - 01025 ! 01026 ! The minmod slope limiter essentiallycorrects the first order upwind 01027 ! estimate for hu/u/v with a second order upwind or central approximation. 01028 ! If the succesive gradients are opposite in sign the correction is set to zero to 01029 ! avoid unwanted oscillations. 01030 ! 01031 ! Note: Declared as pure to explicitly denote that the function has no side effects. 01032 ! If this is not supported in the compiler the pure attribute can be safely removed. 01033 ! 01034 ! ( see for instance Hirch [2001] for a better introduction to limiters. ) 01035 01036 ! 01037 !-------------------------- DEPENDENCIES ---------------------------- 01038 ! 01039 ! - NONE - 01040 ! 01041 !-------------------------- ARGUMENTS ---------------------------- 01042 ! 01043 real(kind=rKind),intent(in) :: delta1 ! first gradient 01044 real(kind=rKind),intent(in) :: delta2 ! secpond gradient 01045 ! 01046 !-------------------------- LOCAL VARIABLES ---------------------------- 01047 ! 01048 if (delta1*delta2 <= 0.0_rKind) then 01049 minmod = 0.0_rKind 01050 return 01051 endif 01052 01053 if (delta1 > 0.0_rKind) then 01054 minmod = min(delta1,delta2) 01055 elseif (delta1 < 0.0_rKind) then 01056 minmod = max(delta1,delta2) 01057 else 01058 minmod = 0.0_rKind 01059 endif 01060 end function minmod 01061 01062 01063 end module flow_secondorder_module
1.7.6.1
