diff options
Diffstat (limited to 'src/include/upwind_second2.h')
| -rw-r--r-- | src/include/upwind_second2.h | 102 |
1 files changed, 51 insertions, 51 deletions
diff --git a/src/include/upwind_second2.h b/src/include/upwind_second2.h index d2e3c17..2f7fdb0 100644 --- a/src/include/upwind_second2.h +++ b/src/include/upwind_second2.h @@ -4,90 +4,90 @@ ! $Header$ ! If this is an active cell... -if ( eh_mask(i,j,k) .ge. 0 ) then +if ( eh_mask(i,j,k,l) .ge. 0 ) then ! If this is not a boundary cell... - if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. & - ( .not. btest(eh_mask(i,j,k),1) ) ) then + if ( ( .not. btest(eh_mask(i,j,k,l),0) ) .and. & + ( .not. btest(eh_mask(i,j,k,l),1) ) ) then ! If the cell to the left is not a boundary cell... - if ( ( .not. btest(eh_mask(i-1,j,k),0) ) .or. & - ( eh_mask(i-1,j,k) .eq. -2 ) ) then + if ( ( .not. btest(eh_mask(i-1,j,k,l),0) ) .or. & + ( eh_mask(i-1,j,k,l) .eq. -2 ) ) then ! Calculate left handed one sided derivative - al = idx * ( three * f(i,j,k) - four * f(i-1,j,k) + f(i-2,j,k) ) + al = idx * ( three * f(i,j,k,l) - four * f(i-1,j,k,l) + f(i-2,j,k,l) ) else al = zero end if ! If the cell to the right is not a boundary cell... - if ( ( .not. btest(eh_mask(i+1,j,k),1) ) .or. & - ( eh_mask(i+1,j,k) .eq. -2 ) ) then + if ( ( .not. btest(eh_mask(i+1,j,k,l),1) ) .or. & + ( eh_mask(i+1,j,k,l) .eq. -2 ) ) then ! Calculate right handed one sided derivative - ar = idx * ( - three * f(i,j,k) + four * f(i+1,j,k) - f(i+2,j,k) ) + ar = idx * ( - three * f(i,j,k,l) + four * f(i+1,j,k,l) - f(i+2,j,k,l) ) else ar = zero end if ! Else if the cell is a left boundary cell. - else if ( btest(eh_mask(i,j,k),0) ) then + else if ( btest(eh_mask(i,j,k,l),0) ) then ! calculate only right handed one sided difference. al = zero - ar = idx * ( - three * f(i,j,k) + four * f(i+1,j,k) - f(i+2,j,k) ) + ar = idx * ( - three * f(i,j,k,l) + four * f(i+1,j,k,l) - f(i+2,j,k,l) ) ! Else it must be a right boundary cell. else ! So calculate only the left handed one sided difference. - al = idx * ( three * f(i,j,k) - four * f(i-1,j,k) + f(i-2,j,k) ) + al = idx * ( three * f(i,j,k,l) - four * f(i-1,j,k,l) + f(i-2,j,k,l) ) ar = zero end if ! Do the same for the y-derivative. - if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. & - ( .not. btest(eh_mask(i,j,k),3) ) ) then - if ( ( .not. btest(eh_mask(i,j-1,k),2) ) .or. & - ( eh_mask(i,j-1,k) .eq. -2 ) ) then - bl = idy * ( three * f(i,j,k) - four * f(i,j-1,k) + f(i,j-2,k) ) + if ( ( .not. btest(eh_mask(i,j,k,l),2) ) .and. & + ( .not. btest(eh_mask(i,j,k,l),3) ) ) then + if ( ( .not. btest(eh_mask(i,j-1,k,l),2) ) .or. & + ( eh_mask(i,j-1,k,l) .eq. -2 ) ) then + bl = idy * ( three * f(i,j,k,l) - four * f(i,j-1,k,l) + f(i,j-2,k,l) ) else bl = zero end if - if ( ( .not. btest(eh_mask(i,j+1,k),3) ) .or. & - ( eh_mask(i,j+1,k) .eq. -2 ) ) then - br = idy * ( - three * f(i,j,k) + four * f(i,j+1,k) - f(i,j+2,k) ) + if ( ( .not. btest(eh_mask(i,j+1,k,l),3) ) .or. & + ( eh_mask(i,j+1,k,l) .eq. -2 ) ) then + br = idy * ( - three * f(i,j,k,l) + four * f(i,j+1,k,l) - f(i,j+2,k,l) ) else br = zero end if - else if ( btest(eh_mask(i,j,k),2) ) then + else if ( btest(eh_mask(i,j,k,l),2) ) then bl = zero - br = idy * ( - three * f(i,j,k) + four * f(i,j+1,k) - f(i,j+2,k) ) + br = idy * ( - three * f(i,j,k,l) + four * f(i,j+1,k,l) - f(i,j+2,k,l) ) else - bl = idy * ( three * f(i,j,k) - four * f(i,j-1,k) + f(i,j-2,k) ) + bl = idy * ( three * f(i,j,k,l) - four * f(i,j-1,k,l) + f(i,j-2,k,l) ) br = zero end if ! And the z-derivative - if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. & - ( .not. btest(eh_mask(i,j,k),5) ) ) then - if ( ( .not. btest(eh_mask(i,j,k-1),4) ) .or. & - ( eh_mask(i,j,k-1) .eq. -2 ) ) then - cl = idz * ( three * f(i,j,k) - four * f(i,j,k-1) + f(i,j,k-2) ) + if ( ( .not. btest(eh_mask(i,j,k,l),4) ) .and. & + ( .not. btest(eh_mask(i,j,k,l),5) ) ) then + if ( ( .not. btest(eh_mask(i,j,k-1,l),4) ) .or. & + ( eh_mask(i,j,k-1,l) .eq. -2 ) ) then + cl = idz * ( three * f(i,j,k,l) - four * f(i,j,k-1,l) + f(i,j,k-2,l) ) else cl = zero end if - if ( ( .not. btest(eh_mask(i,j,k+1),5) ) .or. & - ( eh_mask(i,j,k+1) .eq. -2 ) ) then - cr = idz * ( - three * f(i,j,k) + four * f(i,j,k+1) - f(i,j,k+2) ) + if ( ( .not. btest(eh_mask(i,j,k+1,l),5) ) .or. & + ( eh_mask(i,j,k+1,l) .eq. -2 ) ) then + cr = idz * ( - three * f(i,j,k,l) + four * f(i,j,k+1,l) - f(i,j,k+2,l) ) else cr = zero end if - else if ( btest(eh_mask(i,j,k),4) ) then + else if ( btest(eh_mask(i,j,k,l),4) ) then cl = zero - cr = idz * ( - three * f(i,j,k) + four * f(i,j,k+1) - f(i,j,k+2) ) + cr = idz * ( - three * f(i,j,k,l) + four * f(i,j,k+1,l) - f(i,j,k+2,l) ) else - cl = idz * ( three * f(i,j,k) - four * f(i,j,k-1) + f(i,j,k-2) ) + cl = idz * ( three * f(i,j,k,l) - four * f(i,j,k-1,l) + f(i,j,k-2,l) ) cr = zero end if @@ -112,46 +112,46 @@ if ( eh_mask(i,j,k) .ge. 0 ) then ! If f>0 choose the correct one sided derivative depending on the ! magnitude of the negative and positive parts - if ( f(i,j,k) .gt. 0 ) then + if ( f(i,j,k,l) .gt. 0 ) then if ( abs(alplus) .gt. abs(arminus) ) then - dfx(i,j,k) = alplus + dfx(i,j,k,l) = alplus else - dfx(i,j,k) = arminus + dfx(i,j,k,l) = arminus end if if ( abs(blplus) .gt. abs(brminus) ) then - dfy(i,j,k) = blplus + dfy(i,j,k,l) = blplus else - dfy(i,j,k) = brminus + dfy(i,j,k,l) = brminus end if if ( abs(clplus) .gt. abs(crminus) ) then - dfz(i,j,k) = clplus + dfz(i,j,k,l) = clplus else - dfz(i,j,k) = crminus + dfz(i,j,k,l) = crminus end if endif ! Ditto if f<0. - if ( f(i,j,k) .lt. 0 ) then + if ( f(i,j,k,l) .lt. 0 ) then if ( abs(alminus) .gt. abs(arplus) ) then - dfx(i,j,k) = alminus + dfx(i,j,k,l) = alminus else - dfx(i,j,k) = arplus + dfx(i,j,k,l) = arplus end if if ( abs(blminus) .gt. abs(brplus) ) then - dfy(i,j,k) = blminus + dfy(i,j,k,l) = blminus else - dfy(i,j,k) = brplus + dfy(i,j,k,l) = brplus end if if ( abs(clminus) .gt. abs(crplus) ) then - dfz(i,j,k) = clminus + dfz(i,j,k,l) = clminus else - dfz(i,j,k) = crplus + dfz(i,j,k,l) = crplus end if end if ! If the cell is not active set the derivatives to zero. else - dfx(i,j,k) = zero - dfy(i,j,k) = zero - dfz(i,j,k) = zero + dfx(i,j,k,l) = zero + dfy(i,j,k,l) = zero + dfz(i,j,k,l) = zero end if |