Apply radiation conditions to 3D u,v at open boundaries.
2119 type(ocean_grid_type),
intent(inout) :: G
2120 type(ocean_OBC_type),
pointer :: OBC
2121 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u_new
2124 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u_old
2125 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v_new
2128 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v_old
2129 type(unit_scale_type),
intent(in) :: US
2130 real,
intent(in) :: dt
2132 real :: dhdt, dhdx, dhdy
2133 real :: gamma_u, gamma_2
2135 real :: rx_max, ry_max
2136 real :: rx_new, rx_avg
2137 real :: ry_new, ry_avg
2138 real :: cff_new, cff_avg
2139 real,
allocatable,
dimension(:,:,:) :: &
2148 type(OBC_segment_type),
pointer :: segment => null()
2149 integer :: i, j, k, is, ie, js, je, m, nz, n
2150 integer :: is_obc, ie_obc, js_obc, je_obc
2152 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
2154 if (.not.
associated(obc))
return 2156 if (.not.(obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)) &
2159 eps = 1.0e-20*us%m_s_to_L_T**2
2164 if (obc%gamma_uv < 1.0)
then 2165 do n=1,obc%number_of_segments
2166 segment=>obc%segment(n)
2167 if (.not. segment%on_pe) cycle
2168 if (segment%is_E_or_W .and. segment%radiation)
then 2171 do j=segment%HI%jsd,segment%HI%jed
2172 segment%rx_norm_rad(i,j,k) = obc%rx_normal(i,j,k)
2175 elseif (segment%is_N_or_S .and. segment%radiation)
then 2178 do i=segment%HI%isd,segment%HI%ied
2179 segment%ry_norm_rad(i,j,k) = obc%ry_normal(i,j,k)
2183 if (segment%is_E_or_W .and. segment%oblique)
then 2186 do j=segment%HI%jsd,segment%HI%jed
2187 segment%rx_norm_obl(i,j,k) = obc%rx_oblique(i,j,k)
2188 segment%ry_norm_obl(i,j,k) = obc%ry_oblique(i,j,k)
2189 segment%cff_normal(i,j,k) = obc%cff_normal(i,j,k)
2192 elseif (segment%is_N_or_S .and. segment%oblique)
then 2195 do i=segment%HI%isd,segment%HI%ied
2196 segment%rx_norm_obl(i,j,k) = obc%rx_oblique(i,j,k)
2197 segment%ry_norm_obl(i,j,k) = obc%ry_oblique(i,j,k)
2198 segment%cff_normal(i,j,k) = obc%cff_normal(i,j,k)
2206 do n=1,obc%number_of_segments
2207 segment=>obc%segment(n)
2208 if (
associated(segment%tr_Reg))
then 2209 if (segment%is_E_or_W)
then 2212 if (
associated(segment%tr_Reg%Tr(m)%tres))
then 2214 do j=segment%HI%jsd,segment%HI%jed
2215 segment%tr_Reg%Tr(m)%tres(i,j,k) = obc%tres_x(i,j,k,m)
2223 if (
associated(segment%tr_Reg%Tr(m)%tres))
then 2225 do i=segment%HI%isd,segment%HI%ied
2226 segment%tr_Reg%Tr(m)%tres(i,j,k) = obc%tres_y(i,j,k,m)
2235 gamma_u = obc%gamma_uv
2236 rx_max = obc%rx_max ; ry_max = obc%rx_max
2237 do n=1,obc%number_of_segments
2238 segment=>obc%segment(n)
2239 if (.not. segment%on_pe) cycle
2240 if (segment%oblique)
call gradient_at_q_points(g, segment, u_new(:,:,:), v_new(:,:,:))
2241 if (segment%direction == obc_direction_e)
then 2243 if (i<g%HI%IscB) cycle
2244 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed
2245 if (segment%radiation)
then 2246 dhdt = (u_old(i-1,j,k) - u_new(i-1,j,k))
2247 dhdx = (u_new(i-1,j,k) - u_new(i-2,j,k))
2249 if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max)
2250 if (gamma_u < 1.0)
then 2251 rx_avg = (1.0-gamma_u)*segment%rx_norm_rad(i,j,k) + gamma_u*rx_new
2255 segment%rx_norm_rad(i,j,k) = rx_avg
2259 segment%normal_vel(i,j,k) = (u_new(i,j,k) + rx_avg*u_new(i-1,j,k)) / (1.0+rx_avg)
2262 if (gamma_u < 1.0)
then 2263 obc%rx_normal(i,j,k) = segment%rx_norm_rad(i,j,k)
2265 elseif (segment%oblique)
then 2266 dhdt = (u_old(i-1,j,k) - u_new(i-1,j,k))
2267 dhdx = (u_new(i-1,j,k) - u_new(i-2,j,k))
2268 if (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) > 0.0)
then 2269 dhdy = segment%grad_normal(j-1,1,k)
2270 elseif (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) == 0.0)
then 2273 dhdy = segment%grad_normal(j,1,k)
2275 if (dhdt*dhdx < 0.0) dhdt = 0.0
2276 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2277 rx_new = min(dhdt*dhdx, cff_new*rx_max)
2278 ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
2279 if (gamma_u < 1.0)
then 2280 rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
2281 ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
2282 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
2288 segment%rx_norm_obl(i,j,k) = rx_avg
2289 segment%ry_norm_obl(i,j,k) = ry_avg
2290 segment%cff_normal(i,j,k) = cff_avg
2291 segment%normal_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + rx_avg*u_new(i-1,j,k)) - &
2292 (max(ry_avg,0.0)*segment%grad_normal(j-1,2,k) + &
2293 min(ry_avg,0.0)*segment%grad_normal(j,2,k))) / &
2295 if (gamma_u < 1.0)
then 2298 obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
2299 obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
2300 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
2302 elseif (segment%gradient)
then 2303 segment%normal_vel(i,j,k) = u_new(i-1,j,k)
2305 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then 2307 if (dhdt*dhdx <= 0.0)
then 2308 tau = segment%Velocity_nudging_timescale_in
2310 tau = segment%Velocity_nudging_timescale_out
2312 gamma_2 = dt / (tau + dt)
2313 segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
2314 gamma_2 * segment%nudged_normal_vel(i,j,k)
2317 if (segment%radiation_tan .or. segment%radiation_grad)
then 2319 allocate(rx_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2321 if (gamma_u < 1.0)
then 2322 rx_tang_rad(i,segment%HI%JsdB,k) = segment%rx_norm_rad(i,segment%HI%jsd,k)
2323 rx_tang_rad(i,segment%HI%JedB,k) = segment%rx_norm_rad(i,segment%HI%jed,k)
2324 do j=segment%HI%JsdB+1,segment%HI%JedB-1
2325 rx_tang_rad(i,j,k) = 0.5*(segment%rx_norm_rad(i,j,k) + segment%rx_norm_rad(i,j+1,k))
2328 do j=segment%HI%JsdB,segment%HI%JedB
2329 dhdt = v_old(i,j,k)-v_new(i,j,k)
2330 dhdx = v_new(i,j,k)-v_new(i-1,j,k)
2331 rx_tang_rad(i,j,k) = 0.0
2332 if (dhdt*dhdx > 0.0) rx_tang_rad(i,j,k) = min( (dhdt/dhdx), rx_max)
2336 if (segment%radiation_tan)
then 2337 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2338 rx_avg = rx_tang_rad(i,j,k)
2339 segment%tangential_vel(i,j,k) = (v_new(i,j,k) + rx_avg*v_new(i-1,j,k)) / (1.0+rx_avg)
2342 if (segment%nudged_tan)
then 2343 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2345 if (rx_tang_rad(i,j,k) <= 0.0)
then 2346 tau = segment%Velocity_nudging_timescale_in
2348 tau = segment%Velocity_nudging_timescale_out
2350 gamma_2 = dt / (tau + dt)
2351 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2352 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2355 if (segment%radiation_grad)
then 2356 js_obc = max(segment%HI%JsdB,g%jsd+1)
2357 je_obc = min(segment%HI%JedB,g%jed-1)
2358 do k=1,nz ;
do j=js_obc,je_obc
2359 rx_avg = rx_tang_rad(i,j,k)
2369 segment%tangential_grad(i,j,k) = ((v_new(i,j,k) - v_new(i-1,j,k))*g%IdxBu(i-1,j) + &
2370 rx_avg*(v_new(i-1,j,k) - v_new(i-2,j,k))*g%IdxBu(i-2,j)) / (1.0+rx_avg)
2373 if (segment%nudged_grad)
then 2374 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2376 if (rx_tang_rad(i,j,k) <= 0.0)
then 2377 tau = segment%Velocity_nudging_timescale_in
2379 tau = segment%Velocity_nudging_timescale_out
2381 gamma_2 = dt / (tau + dt)
2382 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2383 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2386 deallocate(rx_tang_rad)
2388 if (segment%oblique_tan .or. segment%oblique_grad)
then 2390 allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2391 allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2392 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2394 if (gamma_u < 1.0)
then 2395 rx_tang_obl(i,segment%HI%JsdB,k) = segment%rx_norm_obl(i,segment%HI%jsd,k)
2396 rx_tang_obl(i,segment%HI%JedB,k) = segment%rx_norm_obl(i,segment%HI%jed,k)
2397 ry_tang_obl(i,segment%HI%JsdB,k) = segment%ry_norm_obl(i,segment%HI%jsd,k)
2398 ry_tang_obl(i,segment%HI%JedB,k) = segment%ry_norm_obl(i,segment%HI%jed,k)
2399 cff_tangential(i,segment%HI%JsdB,k) = segment%cff_normal(i,segment%HI%jsd,k)
2400 cff_tangential(i,segment%HI%JedB,k) = segment%cff_normal(i,segment%HI%jed,k)
2401 do j=segment%HI%JsdB+1,segment%HI%JedB-1
2402 rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i,j+1,k))
2403 ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i,j+1,k))
2404 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i,j+1,k))
2407 do j=segment%HI%JsdB,segment%HI%JedB
2408 dhdt = v_old(i,j,k)-v_new(i,j,k)
2409 dhdx = v_new(i,j,k)-v_new(i-1,j,k)
2410 if (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) > 0.0)
then 2411 dhdy = segment%grad_tan(j,1,k)
2412 elseif (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) == 0.0)
then 2415 dhdy = segment%grad_tan(j+1,1,k)
2417 if (dhdt*dhdx < 0.0) dhdt = 0.0
2418 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2419 rx_new = min(dhdt*dhdx, cff_new*rx_max)
2420 ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
2421 rx_tang_obl(i,j,k) = rx_new
2422 ry_tang_obl(i,j,k) = ry_new
2423 cff_tangential(i,j,k) = cff_new
2427 if (segment%oblique_tan)
then 2428 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2429 rx_avg = rx_tang_obl(i,j,k)
2430 ry_avg = ry_tang_obl(i,j,k)
2431 cff_avg = cff_tangential(i,j,k)
2432 segment%tangential_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + rx_avg*v_new(i-1,j,k)) - &
2433 (max(ry_avg,0.0)*segment%grad_tan(j,2,k) + &
2434 min(ry_avg,0.0)*segment%grad_tan(j+1,2,k))) / &
2438 if (segment%nudged_tan)
then 2439 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2441 if (rx_tang_obl(i,j,k) <= 0.0)
then 2442 tau = segment%Velocity_nudging_timescale_in
2444 tau = segment%Velocity_nudging_timescale_out
2446 gamma_2 = dt / (tau + dt)
2447 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2448 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2451 if (segment%oblique_grad)
then 2452 js_obc = max(segment%HI%JsdB,g%jsd+1)
2453 je_obc = min(segment%HI%JedB,g%jed-1)
2454 do k=1,nz ;
do j=segment%HI%JsdB+1,segment%HI%JedB-1
2455 rx_avg = rx_tang_obl(i,j,k)
2456 ry_avg = ry_tang_obl(i,j,k)
2457 cff_avg = cff_tangential(i,j,k)
2458 segment%tangential_grad(i,j,k) = &
2459 ((cff_avg*(v_new(i,j,k) - v_new(i-1,j,k))*g%IdxBu(i-1,j) + &
2460 rx_avg*(v_new(i-1,j,k) - v_new(i-2,j,k))*g%IdxBu(i-2,j)) - &
2461 (max(ry_avg,0.0)*segment%grad_gradient(j,2,k) + &
2462 min(ry_avg,0.0)*segment%grad_gradient(j+1,2,k)) ) / &
2466 if (segment%nudged_grad)
then 2467 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2469 if (rx_tang_obl(i,j,k) <= 0.0)
then 2470 tau = segment%Velocity_nudging_timescale_in
2472 tau = segment%Velocity_nudging_timescale_out
2474 gamma_2 = dt / (tau + dt)
2475 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2476 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2479 deallocate(rx_tang_obl)
2480 deallocate(ry_tang_obl)
2481 deallocate(cff_tangential)
2485 if (segment%direction == obc_direction_w)
then 2487 if (i>g%HI%IecB) cycle
2488 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed
2489 if (segment%radiation)
then 2490 dhdt = (u_old(i+1,j,k) - u_new(i+1,j,k))
2491 dhdx = (u_new(i+1,j,k) - u_new(i+2,j,k))
2493 if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max)
2494 if (gamma_u < 1.0)
then 2495 rx_avg = (1.0-gamma_u)*segment%rx_norm_rad(i,j,k) + gamma_u*rx_new
2499 segment%rx_norm_rad(i,j,k) = rx_avg
2503 segment%normal_vel(i,j,k) = (u_new(i,j,k) + rx_avg*u_new(i+1,j,k)) / (1.0+rx_avg)
2504 if (gamma_u < 1.0)
then 2507 obc%rx_normal(i,j,k) = segment%rx_norm_rad(i,j,k)
2509 elseif (segment%oblique)
then 2510 dhdt = (u_old(i+1,j,k) - u_new(i+1,j,k))
2511 dhdx = (u_new(i+1,j,k) - u_new(i+2,j,k))
2512 if (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) > 0.0)
then 2513 dhdy = segment%grad_normal(j-1,1,k)
2514 elseif (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) == 0.0)
then 2517 dhdy = segment%grad_normal(j,1,k)
2519 if (dhdt*dhdx < 0.0) dhdt = 0.0
2521 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2522 rx_new = min(dhdt*dhdx, cff_new*rx_max)
2523 ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
2524 if (gamma_u < 1.0)
then 2525 rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
2526 ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
2527 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
2533 segment%rx_norm_obl(i,j,k) = rx_avg
2534 segment%ry_norm_obl(i,j,k) = ry_avg
2535 segment%cff_normal(i,j,k) = cff_avg
2536 segment%normal_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + rx_avg*u_new(i+1,j,k)) - &
2537 (max(ry_avg,0.0)*segment%grad_normal(j-1,2,k) + &
2538 min(ry_avg,0.0)*segment%grad_normal(j,2,k))) / &
2540 if (gamma_u < 1.0)
then 2543 obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
2544 obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
2545 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
2547 elseif (segment%gradient)
then 2548 segment%normal_vel(i,j,k) = u_new(i+1,j,k)
2550 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then 2552 if (dhdt*dhdx <= 0.0)
then 2553 tau = segment%Velocity_nudging_timescale_in
2555 tau = segment%Velocity_nudging_timescale_out
2557 gamma_2 = dt / (tau + dt)
2558 segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
2559 gamma_2 * segment%nudged_normal_vel(i,j,k)
2562 if (segment%radiation_tan .or. segment%radiation_grad)
then 2564 allocate(rx_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2566 if (gamma_u < 1.0)
then 2567 rx_tang_rad(i,segment%HI%JsdB,k) = segment%rx_norm_rad(i,segment%HI%jsd,k)
2568 rx_tang_rad(i,segment%HI%JedB,k) = segment%rx_norm_rad(i,segment%HI%jed,k)
2569 do j=segment%HI%JsdB+1,segment%HI%JedB-1
2570 rx_tang_rad(i,j,k) = 0.5*(segment%rx_norm_rad(i,j,k) + segment%rx_norm_rad(i,j+1,k))
2573 do j=segment%HI%JsdB,segment%HI%JedB
2574 dhdt = v_old(i+1,j,k)-v_new(i+1,j,k)
2575 dhdx = v_new(i+1,j,k)-v_new(i+2,j,k)
2576 rx_tang_rad(i,j,k) = 0.0
2577 if (dhdt*dhdx > 0.0) rx_tang_rad(i,j,k) = min( (dhdt/dhdx), rx_max)
2581 if (segment%radiation_tan)
then 2582 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2583 rx_avg = rx_tang_rad(i,j,k)
2584 segment%tangential_vel(i,j,k) = (v_new(i+1,j,k) + rx_avg*v_new(i+2,j,k)) / (1.0+rx_avg)
2587 if (segment%nudged_tan)
then 2588 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2590 if (rx_tang_rad(i,j,k) <= 0.0)
then 2591 tau = segment%Velocity_nudging_timescale_in
2593 tau = segment%Velocity_nudging_timescale_out
2595 gamma_2 = dt / (tau + dt)
2596 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2597 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2600 if (segment%radiation_grad)
then 2601 js_obc = max(segment%HI%JsdB,g%jsd+1)
2602 je_obc = min(segment%HI%JedB,g%jed-1)
2603 do k=1,nz ;
do j=js_obc,je_obc
2604 rx_avg = rx_tang_rad(i,j,k)
2614 segment%tangential_grad(i,j,k) = ((v_new(i+2,j,k) - v_new(i+1,j,k))*g%IdxBu(i+1,j) + &
2615 rx_avg*(v_new(i+3,j,k) - v_new(i+2,j,k))*g%IdxBu(i+2,j)) / (1.0+rx_avg)
2618 if (segment%nudged_grad)
then 2619 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2621 if (rx_tang_rad(i,j,k) <= 0.0)
then 2622 tau = segment%Velocity_nudging_timescale_in
2624 tau = segment%Velocity_nudging_timescale_out
2626 gamma_2 = dt / (tau + dt)
2627 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2628 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2631 deallocate(rx_tang_rad)
2633 if (segment%oblique_tan .or. segment%oblique_grad)
then 2635 allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2636 allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2637 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2639 if (gamma_u < 1.0)
then 2640 rx_tang_obl(i,segment%HI%JsdB,k) = segment%rx_norm_obl(i,segment%HI%jsd,k)
2641 rx_tang_obl(i,segment%HI%JedB,k) = segment%rx_norm_obl(i,segment%HI%jed,k)
2642 ry_tang_obl(i,segment%HI%JsdB,k) = segment%ry_norm_obl(i,segment%HI%jsd,k)
2643 ry_tang_obl(i,segment%HI%JedB,k) = segment%ry_norm_obl(i,segment%HI%jed,k)
2644 cff_tangential(i,segment%HI%JsdB,k) = segment%cff_normal(i,segment%HI%jsd,k)
2645 cff_tangential(i,segment%HI%JedB,k) = segment%cff_normal(i,segment%HI%jed,k)
2646 do j=segment%HI%JsdB+1,segment%HI%JedB-1
2647 rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i,j+1,k))
2648 ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i,j+1,k))
2649 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i,j+1,k))
2652 do j=segment%HI%JsdB,segment%HI%JedB
2653 dhdt = v_old(i+1,j,k)-v_new(i+1,j,k)
2654 dhdx = v_new(i+1,j,k)-v_new(i+2,j,k)
2655 if (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) > 0.0)
then 2656 dhdy = segment%grad_tan(j,1,k)
2657 elseif (dhdt*(segment%grad_tan(j,1,k) + segment%grad_tan(j+1,1,k)) == 0.0)
then 2660 dhdy = segment%grad_tan(j+1,1,k)
2662 if (dhdt*dhdx < 0.0) dhdt = 0.0
2663 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2664 rx_new = min(dhdt*dhdx, cff_new*rx_max)
2665 ry_new = min(cff_new,max(dhdt*dhdy,-cff_new))
2666 rx_tang_obl(i,j,k) = rx_new
2667 ry_tang_obl(i,j,k) = ry_new
2668 cff_tangential(i,j,k) = cff_new
2672 if (segment%oblique_tan)
then 2673 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2674 rx_avg = rx_tang_obl(i,j,k)
2675 ry_avg = ry_tang_obl(i,j,k)
2676 cff_avg = cff_tangential(i,j,k)
2677 segment%tangential_vel(i,j,k) = ((cff_avg*v_new(i+1,j,k) + rx_avg*v_new(i+2,j,k)) - &
2678 (max(ry_avg,0.0)*segment%grad_tan(j,2,k) + &
2679 min(ry_avg,0.0)*segment%grad_tan(j+1,2,k))) / &
2683 if (segment%nudged_tan)
then 2684 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2686 if (rx_tang_obl(i,j,k) <= 0.0)
then 2687 tau = segment%Velocity_nudging_timescale_in
2689 tau = segment%Velocity_nudging_timescale_out
2691 gamma_2 = dt / (tau + dt)
2692 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2693 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2696 if (segment%oblique_grad)
then 2697 js_obc = max(segment%HI%JsdB,g%jsd+1)
2698 je_obc = min(segment%HI%JedB,g%jed-1)
2699 do k=1,nz ;
do j=segment%HI%JsdB+1,segment%HI%JedB-1
2700 rx_avg = rx_tang_obl(i,j,k)
2701 ry_avg = ry_tang_obl(i,j,k)
2702 cff_avg = cff_tangential(i,j,k)
2703 segment%tangential_grad(i,j,k) = &
2704 ((cff_avg*(v_new(i+2,j,k) - v_new(i+1,j,k))*g%IdxBu(i+1,j) + &
2705 rx_avg*(v_new(i+3,j,k) - v_new(i+2,j,k))*g%IdxBu(i+2,j)) - &
2706 (max(ry_avg,0.0)*segment%grad_gradient(j,2,k) + &
2707 min(ry_avg,0.0)*segment%grad_gradient(j+1,2,k))) / &
2711 if (segment%nudged_grad)
then 2712 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2714 if (rx_tang_obl(i,j,k) <= 0.0)
then 2715 tau = segment%Velocity_nudging_timescale_in
2717 tau = segment%Velocity_nudging_timescale_out
2719 gamma_2 = dt / (tau + dt)
2720 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2721 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2724 deallocate(rx_tang_obl)
2725 deallocate(ry_tang_obl)
2726 deallocate(cff_tangential)
2730 if (segment%direction == obc_direction_n)
then 2732 if (j<g%HI%JscB) cycle
2733 do k=1,nz ;
do i=segment%HI%isd,segment%HI%ied
2734 if (segment%radiation)
then 2735 dhdt = (v_old(i,j-1,k) - v_new(i,j-1,k))
2736 dhdy = (v_new(i,j-1,k) - v_new(i,j-2,k))
2738 if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max)
2739 if (gamma_u < 1.0)
then 2740 ry_avg = (1.0-gamma_u)*segment%ry_norm_rad(i,j,k) + gamma_u*ry_new
2744 segment%ry_norm_rad(i,j,k) = ry_avg
2748 segment%normal_vel(i,j,k) = (v_new(i,j,k) + ry_avg*v_new(i,j-1,k)) / (1.0+ry_avg)
2749 if (gamma_u < 1.0)
then 2752 obc%ry_normal(i,j,k) = segment%ry_norm_rad(i,j,k)
2754 elseif (segment%oblique)
then 2755 dhdt = (v_old(i,j-1,k) - v_new(i,j-1,k))
2756 dhdy = (v_new(i,j-1,k) - v_new(i,j-2,k))
2757 if (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) > 0.0)
then 2758 dhdx = segment%grad_normal(i-1,1,k)
2759 elseif (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) == 0.0)
then 2762 dhdx = segment%grad_normal(i,1,k)
2764 if (dhdt*dhdy < 0.0) dhdt = 0.0
2765 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2766 ry_new = min(dhdt*dhdy, cff_new*ry_max)
2767 rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
2768 if (gamma_u < 1.0)
then 2769 rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
2770 ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
2771 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
2777 segment%rx_norm_obl(i,j,k) = rx_avg
2778 segment%ry_norm_obl(i,j,k) = ry_avg
2779 segment%cff_normal(i,j,k) = cff_avg
2780 segment%normal_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + ry_avg*v_new(i,j-1,k)) - &
2781 (max(rx_avg,0.0)*segment%grad_normal(i-1,2,k) +&
2782 min(rx_avg,0.0)*segment%grad_normal(i,2,k))) / &
2784 if (gamma_u < 1.0)
then 2787 obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
2788 obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
2789 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
2791 elseif (segment%gradient)
then 2792 segment%normal_vel(i,j,k) = v_new(i,j-1,k)
2794 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then 2796 if (dhdt*dhdy <= 0.0)
then 2797 tau = segment%Velocity_nudging_timescale_in
2799 tau = segment%Velocity_nudging_timescale_out
2801 gamma_2 = dt / (tau + dt)
2802 segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
2803 gamma_2 * segment%nudged_normal_vel(i,j,k)
2806 if (segment%radiation_tan .or. segment%radiation_grad)
then 2808 allocate(ry_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2810 if (gamma_u < 1.0)
then 2811 ry_tang_rad(segment%HI%IsdB,j,k) = segment%ry_norm_rad(segment%HI%isd,j,k)
2812 ry_tang_rad(segment%HI%IedB,j,k) = segment%ry_norm_rad(segment%HI%ied,j,k)
2813 do i=segment%HI%IsdB+1,segment%HI%IedB-1
2814 ry_tang_rad(i,j,k) = 0.5*(segment%ry_norm_rad(i,j,k) + segment%ry_norm_rad(i+1,j,k))
2817 do i=segment%HI%IsdB,segment%HI%IedB
2818 dhdt = u_old(i,j-1,k)-u_new(i,j-1,k)
2819 dhdy = u_new(i,j-1,k)-u_new(i,j-2,k)
2820 ry_tang_rad(i,j,k) = 0.0
2821 if (dhdt*dhdy > 0.0) ry_tang_rad(i,j,k) = min( (dhdt/dhdy), rx_max)
2825 if (segment%radiation_tan)
then 2826 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2827 ry_avg = ry_tang_rad(i,j,k)
2828 segment%tangential_vel(i,j,k) = (u_new(i,j,k) + ry_avg*u_new(i,j-1,k)) / (1.0+ry_avg)
2831 if (segment%nudged_tan)
then 2832 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2834 if (ry_tang_rad(i,j,k) <= 0.0)
then 2835 tau = segment%Velocity_nudging_timescale_in
2837 tau = segment%Velocity_nudging_timescale_out
2839 gamma_2 = dt / (tau + dt)
2840 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2841 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2844 if (segment%radiation_grad)
then 2845 is_obc = max(segment%HI%IsdB,g%isd+1)
2846 ie_obc = min(segment%HI%IedB,g%ied-1)
2847 do k=1,nz ;
do i=is_obc,ie_obc
2848 ry_avg = ry_tang_rad(i,j,k)
2858 segment%tangential_grad(i,j,k) = ((u_new(i,j,k) - u_new(i,j-1,k))*g%IdyBu(i,j-1) + &
2859 ry_avg*(u_new(i,j-1,k) - u_new(i,j-2,k))*g%IdyBu(i,j-2)) / (1.0+ry_avg)
2862 if (segment%nudged_grad)
then 2863 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2865 if (ry_tang_rad(i,j,k) <= 0.0)
then 2866 tau = segment%Velocity_nudging_timescale_in
2868 tau = segment%Velocity_nudging_timescale_out
2870 gamma_2 = dt / (tau + dt)
2871 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2872 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2875 deallocate(ry_tang_rad)
2877 if (segment%oblique_tan .or. segment%oblique_grad)
then 2879 allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2880 allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2881 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2883 if (gamma_u < 1.0)
then 2884 rx_tang_obl(segment%HI%IsdB,j,k) = segment%rx_norm_obl(segment%HI%isd,j,k)
2885 rx_tang_obl(segment%HI%IedB,j,k) = segment%rx_norm_obl(segment%HI%ied,j,k)
2886 ry_tang_obl(segment%HI%IsdB,j,k) = segment%ry_norm_obl(segment%HI%isd,j,k)
2887 ry_tang_obl(segment%HI%IedB,j,k) = segment%ry_norm_obl(segment%HI%ied,j,k)
2888 cff_tangential(segment%HI%IsdB,j,k) = segment%cff_normal(segment%HI%isd,j,k)
2889 cff_tangential(segment%HI%IedB,j,k) = segment%cff_normal(segment%HI%ied,j,k)
2890 do i=segment%HI%IsdB+1,segment%HI%IedB-1
2891 rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i+1,j,k))
2892 ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i+1,j,k))
2893 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i+1,j,k))
2896 do i=segment%HI%IsdB,segment%HI%IedB
2897 dhdt = u_old(i,j,k)-u_new(i,j,k)
2898 dhdy = u_new(i,j,k)-u_new(i,j-1,k)
2899 if (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) > 0.0)
then 2900 dhdx = segment%grad_tan(i,1,k)
2901 elseif (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) == 0.0)
then 2904 dhdx = segment%grad_tan(i+1,1,k)
2906 if (dhdt*dhdy < 0.0) dhdt = 0.0
2907 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2908 ry_new = min(dhdt*dhdy, cff_new*ry_max)
2909 rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
2910 rx_tang_obl(i,j,k) = rx_new
2911 ry_tang_obl(i,j,k) = ry_new
2912 cff_tangential(i,j,k) = cff_new
2916 if (segment%oblique_tan)
then 2917 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2918 rx_avg = rx_tang_obl(i,j,k)
2919 ry_avg = ry_tang_obl(i,j,k)
2920 cff_avg = cff_tangential(i,j,k)
2921 segment%tangential_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + ry_avg*u_new(i,j-1,k)) - &
2922 (max(rx_avg,0.0)*segment%grad_tan(i,2,k) + &
2923 min(rx_avg,0.0)*segment%grad_tan(i+1,2,k))) / &
2927 if (segment%nudged_tan)
then 2928 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2930 if (ry_tang_obl(i,j,k) <= 0.0)
then 2931 tau = segment%Velocity_nudging_timescale_in
2933 tau = segment%Velocity_nudging_timescale_out
2935 gamma_2 = dt / (tau + dt)
2936 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
2937 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2940 if (segment%oblique_grad)
then 2941 is_obc = max(segment%HI%IsdB,g%isd+1)
2942 ie_obc = min(segment%HI%IedB,g%ied-1)
2943 do k=1,nz ;
do i=segment%HI%IsdB+1,segment%HI%IedB-1
2944 rx_avg = rx_tang_obl(i,j,k)
2945 ry_avg = ry_tang_obl(i,j,k)
2946 cff_avg = cff_tangential(i,j,k)
2947 segment%tangential_grad(i,j,k) = &
2948 ((cff_avg*(u_new(i,j,k) - u_new(i,j-1,k))*g%IdyBu(i,j-1) + &
2949 ry_avg*(u_new(i,j-1,k) - u_new(i,j-2,k))*g%IdyBu(i,j-2)) - &
2950 (max(rx_avg,0.0)*segment%grad_gradient(i,2,k) + &
2951 min(rx_avg,0.0)*segment%grad_gradient(i+1,2,k))) / &
2955 if (segment%nudged_grad)
then 2956 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2958 if (ry_tang_obl(i,j,k) <= 0.0)
then 2959 tau = segment%Velocity_nudging_timescale_in
2961 tau = segment%Velocity_nudging_timescale_out
2963 gamma_2 = dt / (tau + dt)
2964 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
2965 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2968 deallocate(rx_tang_obl)
2969 deallocate(ry_tang_obl)
2970 deallocate(cff_tangential)
2974 if (segment%direction == obc_direction_s)
then 2976 if (j>g%HI%JecB) cycle
2977 do k=1,nz ;
do i=segment%HI%isd,segment%HI%ied
2978 if (segment%radiation)
then 2979 dhdt = (v_old(i,j+1,k) - v_new(i,j+1,k))
2980 dhdy = (v_new(i,j+1,k) - v_new(i,j+2,k))
2982 if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max)
2983 if (gamma_u < 1.0)
then 2984 ry_avg = (1.0-gamma_u)*segment%ry_norm_rad(i,j,k) + gamma_u*ry_new
2988 segment%ry_norm_rad(i,j,k) = ry_avg
2992 segment%normal_vel(i,j,k) = (v_new(i,j,k) + ry_avg*v_new(i,j+1,k)) / (1.0+ry_avg)
2993 if (gamma_u < 1.0)
then 2996 obc%ry_normal(i,j,k) = segment%ry_norm_rad(i,j,k)
2998 elseif (segment%oblique)
then 2999 dhdt = (v_old(i,j+1,k) - v_new(i,j+1,k))
3000 dhdy = (v_new(i,j+1,k) - v_new(i,j+2,k))
3001 if (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) > 0.0)
then 3002 dhdx = segment%grad_normal(i-1,1,k)
3003 elseif (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) == 0.0)
then 3006 dhdx = segment%grad_normal(i,1,k)
3008 if (dhdt*dhdy < 0.0) dhdt = 0.0
3010 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
3011 ry_new = min(dhdt*dhdy, cff_new*ry_max)
3012 rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
3013 if (gamma_u < 1.0)
then 3014 rx_avg = (1.0-gamma_u)*segment%rx_norm_obl(i,j,k) + gamma_u*rx_new
3015 ry_avg = (1.0-gamma_u)*segment%ry_norm_obl(i,j,k) + gamma_u*ry_new
3016 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
3022 segment%rx_norm_obl(i,j,k) = rx_avg
3023 segment%ry_norm_obl(i,j,k) = ry_avg
3024 segment%cff_normal(i,j,k) = cff_avg
3025 segment%normal_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + ry_avg*v_new(i,j+1,k)) - &
3026 (max(rx_avg,0.0)*segment%grad_normal(i-1,2,k) + &
3027 min(rx_avg,0.0)*segment%grad_normal(i,2,k))) / &
3029 if (gamma_u < 1.0)
then 3032 obc%rx_oblique(i,j,k) = segment%rx_norm_obl(i,j,k)
3033 obc%ry_oblique(i,j,k) = segment%ry_norm_obl(i,j,k)
3034 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
3036 elseif (segment%gradient)
then 3037 segment%normal_vel(i,j,k) = v_new(i,j+1,k)
3039 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then 3041 if (dhdt*dhdy <= 0.0)
then 3042 tau = segment%Velocity_nudging_timescale_in
3044 tau = segment%Velocity_nudging_timescale_out
3046 gamma_2 = dt / (tau + dt)
3047 segment%normal_vel(i,j,k) = (1.0 - gamma_2) * segment%normal_vel(i,j,k) + &
3048 gamma_2 * segment%nudged_normal_vel(i,j,k)
3051 if (segment%radiation_tan .or. segment%radiation_grad)
then 3053 allocate(ry_tang_rad(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
3055 if (gamma_u < 1.0)
then 3056 ry_tang_rad(segment%HI%IsdB,j,k) = segment%ry_norm_rad(segment%HI%isd,j,k)
3057 ry_tang_rad(segment%HI%IedB,j,k) = segment%ry_norm_rad(segment%HI%ied,j,k)
3058 do i=segment%HI%IsdB+1,segment%HI%IedB-1
3059 ry_tang_rad(i,j,k) = 0.5*(segment%ry_norm_rad(i,j,k) + segment%ry_norm_rad(i+1,j,k))
3062 do i=segment%HI%IsdB,segment%HI%IedB
3063 dhdt = u_old(i,j+1,k)-u_new(i,j+1,k)
3064 dhdy = u_new(i,j+1,k)-u_new(i,j+2,k)
3065 ry_tang_rad(i,j,k) = 0.0
3066 if (dhdt*dhdy > 0.0) ry_tang_rad(i,j,k) = min( (dhdt/dhdy), rx_max)
3070 if (segment%radiation_tan)
then 3071 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
3072 ry_avg = ry_tang_rad(i,j,k)
3073 segment%tangential_vel(i,j,k) = (u_new(i,j+1,k) + ry_avg*u_new(i,j+2,k)) / (1.0+ry_avg)
3076 if (segment%nudged_tan)
then 3077 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
3079 if (ry_tang_rad(i,j,k) <= 0.0)
then 3080 tau = segment%Velocity_nudging_timescale_in
3082 tau = segment%Velocity_nudging_timescale_out
3084 gamma_2 = dt / (tau + dt)
3085 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
3086 gamma_2 * segment%nudged_tangential_vel(i,j,k)
3089 if (segment%radiation_grad)
then 3090 is_obc = max(segment%HI%IsdB,g%isd+1)
3091 ie_obc = min(segment%HI%IedB,g%ied-1)
3092 do k=1,nz ;
do i=is_obc,ie_obc
3093 ry_avg = ry_tang_rad(i,j,k)
3103 segment%tangential_grad(i,j,k) = ((u_new(i,j+2,k) - u_new(i,j+1,k))*g%IdyBu(i,j+1) + &
3104 ry_avg*(u_new(i,j+3,k) - u_new(i,j+2,k))*g%IdyBu(i,j+2)) / (1.0+ry_avg)
3107 if (segment%nudged_grad)
then 3108 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
3110 if (ry_tang_rad(i,j,k) <= 0.0)
then 3111 tau = segment%Velocity_nudging_timescale_in
3113 tau = segment%Velocity_nudging_timescale_out
3115 gamma_2 = dt / (tau + dt)
3116 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
3117 gamma_2 * segment%nudged_tangential_grad(i,j,k)
3120 deallocate(ry_tang_rad)
3122 if (segment%oblique_tan .or. segment%oblique_grad)
then 3124 allocate(rx_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
3125 allocate(ry_tang_obl(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
3126 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
3128 if (gamma_u < 1.0)
then 3129 rx_tang_obl(segment%HI%IsdB,j,k) = segment%rx_norm_obl(segment%HI%isd,j,k)
3130 rx_tang_obl(segment%HI%IedB,j,k) = segment%rx_norm_obl(segment%HI%ied,j,k)
3131 ry_tang_obl(segment%HI%IsdB,j,k) = segment%ry_norm_obl(segment%HI%isd,j,k)
3132 ry_tang_obl(segment%HI%IedB,j,k) = segment%ry_norm_obl(segment%HI%ied,j,k)
3133 cff_tangential(segment%HI%IsdB,j,k) = segment%cff_normal(segment%HI%isd,j,k)
3134 cff_tangential(segment%HI%IedB,j,k) = segment%cff_normal(segment%HI%ied,j,k)
3135 do i=segment%HI%IsdB+1,segment%HI%IedB-1
3136 rx_tang_obl(i,j,k) = 0.5*(segment%rx_norm_obl(i,j,k) + segment%rx_norm_obl(i+1,j,k))
3137 ry_tang_obl(i,j,k) = 0.5*(segment%ry_norm_obl(i,j,k) + segment%ry_norm_obl(i+1,j,k))
3138 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i+1,j,k))
3141 do i=segment%HI%IsdB,segment%HI%IedB
3142 dhdt = u_old(i,j+1,k)-u_new(i,j+1,k)
3143 dhdy = u_new(i,j+1,k)-u_new(i,j+2,k)
3144 if (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) > 0.0)
then 3145 dhdx = segment%grad_tan(i,1,k)
3146 elseif (dhdt*(segment%grad_tan(i,1,k) + segment%grad_tan(i+1,1,k)) == 0.0)
then 3149 dhdx = segment%grad_tan(i+1,1,k)
3151 if (dhdt*dhdy < 0.0) dhdt = 0.0
3152 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
3153 ry_new = min(dhdt*dhdy, cff_new*ry_max)
3154 rx_new = min(cff_new,max(dhdt*dhdx,-cff_new))
3155 rx_tang_obl(i,j,k) = rx_new
3156 ry_tang_obl(i,j,k) = ry_new
3157 cff_tangential(i,j,k) = cff_new
3161 if (segment%oblique_tan)
then 3162 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
3163 rx_avg = rx_tang_obl(i,j,k)
3164 ry_avg = ry_tang_obl(i,j,k)
3165 cff_avg = cff_tangential(i,j,k)
3166 segment%tangential_vel(i,j,k) = ((cff_avg*u_new(i,j+1,k) + ry_avg*u_new(i,j+2,k)) - &
3167 (max(rx_avg,0.0)*segment%grad_tan(i,2,k) + &
3168 min(rx_avg,0.0)*segment%grad_tan(i+1,2,k)) ) / &
3172 if (segment%nudged_tan)
then 3173 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
3175 if (ry_tang_obl(i,j,k) <= 0.0)
then 3176 tau = segment%Velocity_nudging_timescale_in
3178 tau = segment%Velocity_nudging_timescale_out
3180 gamma_2 = dt / (tau + dt)
3181 segment%tangential_vel(i,j,k) = (1.0 - gamma_2) * segment%tangential_vel(i,j,k) + &
3182 gamma_2 * segment%nudged_tangential_vel(i,j,k)
3185 if (segment%oblique_grad)
then 3186 is_obc = max(segment%HI%IsdB,g%isd+1)
3187 ie_obc = min(segment%HI%IedB,g%ied-1)
3188 do k=1,nz ;
do i=segment%HI%IsdB+1,segment%HI%IedB-1
3189 rx_avg = rx_tang_obl(i,j,k)
3190 ry_avg = ry_tang_obl(i,j,k)
3191 cff_avg = cff_tangential(i,j,k)
3192 segment%tangential_grad(i,j,k) = &
3193 ((cff_avg*(u_new(i,j+2,k) - u_new(i,j+1,k))*g%IdyBu(i,j+1) + &
3194 ry_avg*(u_new(i,j+3,k) - u_new(i,j+2,k))*g%IdyBu(i,j+2)) - &
3195 (max(rx_avg,0.0)*segment%grad_gradient(i,2,k) + &
3196 min(rx_avg,0.0)*segment%grad_gradient(i+1,2,k))) / &
3200 if (segment%nudged_grad)
then 3201 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
3203 if (ry_tang_obl(i,j,k) <= 0.0)
then 3204 tau = segment%Velocity_nudging_timescale_in
3206 tau = segment%Velocity_nudging_timescale_out
3208 gamma_2 = dt / (tau + dt)
3209 segment%tangential_grad(i,j,k) = (1.0 - gamma_2) * segment%tangential_grad(i,j,k) + &
3210 gamma_2 * segment%nudged_tangential_grad(i,j,k)
3213 deallocate(rx_tang_obl)
3214 deallocate(ry_tang_obl)
3215 deallocate(cff_tangential)
3221 call open_boundary_apply_normal_flow(obc, g, u_new, v_new)
3223 call pass_vector(u_new, v_new, g%Domain, clock=id_clock_pass)