Apply radiation conditions to 3D u,v at open boundaries.
1523 type(ocean_grid_type),
intent(inout) :: G
1524 type(ocean_OBC_type),
pointer :: OBC
1525 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(inout) :: u_new
1528 real,
dimension(SZIB_(G),SZJ_(G),SZK_(G)),
intent(in) :: u_old
1529 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(inout) :: v_new
1532 real,
dimension(SZI_(G),SZJB_(G),SZK_(G)),
intent(in) :: v_old
1533 real,
intent(in) :: dt
1535 real :: dhdt, dhdx, dhdy, gamma_u, gamma_v, gamma_2
1536 real :: cff, Cx, Cy, tau
1537 real :: rx_max, ry_max
1538 real :: rx_new, rx_avg
1539 real :: ry_new, ry_avg
1540 real :: cff_new, cff_avg
1541 real,
pointer,
dimension(:,:,:) :: rx_tangential=>null()
1542 real,
pointer,
dimension(:,:,:) :: ry_tangential=>null()
1543 real,
pointer,
dimension(:,:,:) :: cff_tangential=>null()
1544 real,
parameter :: eps = 1.0e-20
1545 type(OBC_segment_type),
pointer :: segment => null()
1546 integer :: i, j, k, is, ie, js, je, nz, n
1547 integer :: is_obc, ie_obc, js_obc, je_obc
1549 is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
1551 if (.not.
associated(obc))
return
1553 if (.not.(obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)) &
1559 do n=1,obc%number_of_segments
1560 segment=>obc%segment(n)
1561 if (.not. segment%on_pe) cycle
1562 if (segment%is_E_or_W .and. segment%radiation)
then
1565 do j=segment%HI%jsd,segment%HI%jed
1566 segment%rx_normal(i,j,k) = obc%rx_normal(i,j,k)
1569 elseif (segment%is_N_or_S .and. segment%radiation)
then
1572 do i=segment%HI%isd,segment%HI%ied
1573 segment%ry_normal(i,j,k) = obc%ry_normal(i,j,k)
1577 if (segment%is_E_or_W .and. segment%oblique)
then
1580 do j=segment%HI%jsd,segment%HI%jed
1581 segment%rx_normal(i,j,k) = obc%rx_normal(i,j,k)
1582 segment%ry_normal(i,j,k) = obc%ry_normal(i,j,k)
1583 segment%cff_normal(i,j,k) = obc%cff_normal(i,j,k)
1586 elseif (segment%is_N_or_S .and. segment%oblique)
then
1589 do i=segment%HI%isd,segment%HI%ied
1590 segment%rx_normal(i,j,k) = obc%rx_normal(i,j,k)
1591 segment%ry_normal(i,j,k) = obc%ry_normal(i,j,k)
1592 segment%cff_normal(i,j,k) = obc%cff_normal(i,j,k)
1598 gamma_u = obc%gamma_uv ; gamma_v = obc%gamma_uv
1599 rx_max = obc%rx_max ; ry_max = obc%rx_max
1600 do n=1,obc%number_of_segments
1601 segment=>obc%segment(n)
1602 if (.not. segment%on_pe) cycle
1603 if (segment%oblique)
call gradient_at_q_points(g,segment,u_new,v_new)
1604 if (segment%direction == obc_direction_e)
then
1606 if (i<g%HI%IscB) cycle
1607 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed
1608 if (segment%radiation)
then
1609 dhdt = u_old(i-1,j,k)-u_new(i-1,j,k)
1610 dhdx = u_new(i-1,j,k)-u_new(i-2,j,k)
1612 if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max)
1613 rx_avg = (1.0-gamma_u)*segment%rx_normal(i,j,k) + gamma_u*rx_new
1614 segment%rx_normal(i,j,k) = rx_avg
1618 segment%normal_vel(i,j,k) = (u_new(i,j,k) + rx_avg*u_new(i-1,j,k)) / (1.0+rx_avg)
1621 obc%rx_normal(i,j,k) = segment%rx_normal(i,j,k)
1622 elseif (segment%oblique)
then
1623 dhdt = u_old(i-1,j,k)-u_new(i-1,j,k)
1624 dhdx = u_new(i-1,j,k)-u_new(i-2,j,k)
1625 if (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) > 0.0)
then
1626 dhdy = segment%grad_normal(j-1,1,k)
1627 elseif (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) == 0.0)
then
1630 dhdy = segment%grad_normal(j,1,k)
1632 if (dhdt*dhdx < 0.0) dhdt = 0.0
1634 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
1635 ry_new = min(cff,max(dhdt*dhdy,-cff))
1636 rx_avg = (1.0-gamma_u)*segment%rx_normal(i,j,k) + gamma_u*rx_new
1637 ry_avg = (1.0-gamma_u)*segment%ry_normal(i,j,k) + gamma_u*ry_new
1638 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
1639 segment%rx_normal(i,j,k) = rx_avg
1640 segment%ry_normal(i,j,k) = ry_avg
1641 segment%cff_normal(i,j,k) = cff_avg
1642 segment%normal_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + rx_avg*u_new(i-1,j,k)) - &
1643 (max(ry_avg,0.0)*segment%grad_normal(j-1,2,k) + min(ry_avg,0.0)*segment%grad_normal(j,2,k))) / &
1647 obc%rx_normal(i,j,k) = segment%rx_normal(i,j,k)
1648 obc%ry_normal(i,j,k) = segment%ry_normal(i,j,k)
1649 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
1650 elseif (segment%gradient)
then
1651 segment%normal_vel(i,j,k) = u_new(i-1,j,k)
1653 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then
1655 if (dhdt*dhdx <= 0.0)
then
1656 tau = segment%Velocity_nudging_timescale_in
1658 tau = segment%Velocity_nudging_timescale_out
1660 gamma_2 = dt / (tau + dt)
1661 segment%normal_vel(i,j,k) = (1 - gamma_2) * segment%normal_vel(i,j,k) + &
1662 gamma_2 * segment%nudged_normal_vel(i,j,k)
1665 if (segment%radiation_tan .or. segment%radiation_grad)
then
1667 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1669 rx_tangential(i,segment%HI%JsdB,k) = segment%rx_normal(i,segment%HI%jsd,k)
1670 rx_tangential(i,segment%HI%JedB,k) = segment%rx_normal(i,segment%HI%jed,k)
1671 do j=segment%HI%JsdB+1,segment%HI%JedB-1
1672 rx_tangential(i,j,k) = 0.5*(segment%rx_normal(i,j,k) + segment%rx_normal(i,j+1,k))
1675 if (segment%radiation_tan)
then
1676 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1677 rx_avg = rx_tangential(i,j,k)
1678 segment%tangential_vel(i,j,k) = (v_new(i,j,k) + rx_avg*v_new(i-1,j,k)) / (1.0+rx_avg)
1681 if (segment%nudged_tan)
then
1682 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1684 if (rx_tangential(i,j,k) <= 0.0)
then
1685 tau = segment%Velocity_nudging_timescale_in
1687 tau = segment%Velocity_nudging_timescale_out
1689 gamma_2 = dt / (tau + dt)
1690 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
1691 gamma_2 * segment%nudged_tangential_vel(i,j,k)
1694 if (segment%radiation_grad)
then
1695 js_obc = max(segment%HI%JsdB,g%jsd+1)
1696 je_obc = min(segment%HI%JedB,g%jed-1)
1697 do k=1,nz ;
do j=js_obc,je_obc
1698 rx_avg = rx_tangential(i,j,k)
1708 segment%tangential_grad(i,j,k) = ((v_new(i,j,k) - v_new(i-1,j,k))*g%IdxBu(i-1,j) + &
1709 rx_avg*(v_new(i-1,j,k) - v_new(i-2,j,k))*g%IdxBu(i-2,j)) / (1.0+rx_avg)
1712 if (segment%nudged_grad)
then
1713 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1715 if (rx_tangential(i,j,k) <= 0.0)
then
1716 tau = segment%Velocity_nudging_timescale_in
1718 tau = segment%Velocity_nudging_timescale_out
1720 gamma_2 = dt / (tau + dt)
1721 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
1722 gamma_2 * segment%nudged_tangential_grad(i,j,k)
1725 deallocate(rx_tangential)
1727 if (segment%oblique_tan .or. segment%oblique_grad)
then
1729 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1730 allocate(ry_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1731 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1733 rx_tangential(i,segment%HI%JsdB,k) = segment%rx_normal(i,segment%HI%jsd,k)
1734 rx_tangential(i,segment%HI%JedB,k) = segment%rx_normal(i,segment%HI%jed,k)
1735 ry_tangential(i,segment%HI%JsdB,k) = segment%ry_normal(i,segment%HI%jsd,k)
1736 ry_tangential(i,segment%HI%JedB,k) = segment%ry_normal(i,segment%HI%jed,k)
1737 cff_tangential(i,segment%HI%JsdB,k) = segment%cff_normal(i,segment%HI%jsd,k)
1738 cff_tangential(i,segment%HI%JedB,k) = segment%cff_normal(i,segment%HI%jed,k)
1739 do j=segment%HI%JsdB+1,segment%HI%JedB-1
1740 rx_tangential(i,j,k) = 0.5*(segment%rx_normal(i,j,k) + segment%rx_normal(i,j+1,k))
1741 ry_tangential(i,j,k) = 0.5*(segment%ry_normal(i,j,k) + segment%ry_normal(i,j+1,k))
1742 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i,j+1,k))
1745 if (segment%oblique_tan)
then
1746 do k=1,nz ;
do j=segment%HI%JsdB+1,segment%HI%JedB-1
1747 rx_avg = rx_tangential(i,j,k)
1748 ry_avg = ry_tangential(i,j,k)
1749 cff_avg = cff_tangential(i,j,k)
1750 segment%tangential_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + rx_avg*v_new(i-1,j,k)) - &
1751 (max(ry_avg,0.0)*segment%grad_tan(j,2,k) + min(ry_avg,0.0)*segment%grad_tan(j+1,2,k))) / &
1755 if (segment%nudged_tan)
then
1756 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1758 if (rx_tangential(i,j,k) <= 0.0)
then
1759 tau = segment%Velocity_nudging_timescale_in
1761 tau = segment%Velocity_nudging_timescale_out
1763 gamma_2 = dt / (tau + dt)
1764 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
1765 gamma_2 * segment%nudged_tangential_vel(i,j,k)
1768 if (segment%oblique_grad)
then
1769 js_obc = max(segment%HI%JsdB,g%jsd+1)
1770 je_obc = min(segment%HI%JedB,g%jed-1)
1771 do k=1,nz ;
do j=segment%HI%JsdB+1,segment%HI%JedB-1
1772 rx_avg = rx_tangential(i,j,k)
1773 ry_avg = ry_tangential(i,j,k)
1774 cff_avg = cff_tangential(i,j,k)
1775 segment%tangential_grad(i,j,k) = ((cff_avg*(v_new(i,j,k) - v_new(i-1,j,k))*g%IdxBu(i-1,j) &
1776 + rx_avg*(v_new(i-1,j,k) - v_new(i-2,j,k))*g%IdxBu(i-2,j)) - &
1777 (max(ry_avg,0.0)*segment%grad_gradient(j,2,k) + min(ry_avg,0.0)*segment%grad_gradient(j+1,2,k))) / &
1781 if (segment%nudged_grad)
then
1782 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1784 if (rx_tangential(i,j,k) <= 0.0)
then
1785 tau = segment%Velocity_nudging_timescale_in
1787 tau = segment%Velocity_nudging_timescale_out
1789 gamma_2 = dt / (tau + dt)
1790 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
1791 gamma_2 * segment%nudged_tangential_grad(i,j,k)
1794 deallocate(rx_tangential)
1795 deallocate(ry_tangential)
1796 deallocate(cff_tangential)
1800 if (segment%direction == obc_direction_w)
then
1802 if (i>g%HI%IecB) cycle
1803 do k=1,nz ;
do j=segment%HI%jsd,segment%HI%jed
1804 if (segment%radiation)
then
1805 dhdt = u_old(i+1,j,k)-u_new(i+1,j,k)
1806 dhdx = u_new(i+1,j,k)-u_new(i+2,j,k)
1808 if (dhdt*dhdx > 0.0) rx_new = min( (dhdt/dhdx), rx_max)
1809 rx_avg = (1.0-gamma_u)*segment%rx_normal(i,j,k) + gamma_u*rx_new
1810 segment%rx_normal(i,j,k) = rx_avg
1814 segment%normal_vel(i,j,k) = (u_new(i,j,k) + rx_avg*u_new(i+1,j,k)) / (1.0+rx_avg)
1817 obc%rx_normal(i,j,k) = segment%rx_normal(i,j,k)
1818 elseif (segment%oblique)
then
1819 dhdt = u_old(i+1,j,k)-u_new(i+1,j,k)
1820 dhdx = u_new(i+1,j,k)-u_new(i+2,j,k)
1821 if (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) > 0.0)
then
1822 dhdy = segment%grad_normal(j-1,1,k)
1823 elseif (dhdt*(segment%grad_normal(j,1,k) + segment%grad_normal(j-1,1,k)) == 0.0)
then
1826 dhdy = segment%grad_normal(j,1,k)
1828 if (dhdt*dhdx < 0.0) dhdt = 0.0
1830 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
1831 ry_new = min(cff,max(dhdt*dhdy,-cff))
1832 rx_avg = (1.0-gamma_u)*segment%rx_normal(i,j,k) + gamma_u*rx_new
1833 ry_avg = (1.0-gamma_u)*segment%ry_normal(i,j,k) + gamma_u*ry_new
1834 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
1835 segment%rx_normal(i,j,k) = rx_avg
1836 segment%ry_normal(i,j,k) = ry_avg
1837 segment%cff_normal(i,j,k) = cff_avg
1838 segment%normal_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + rx_avg*u_new(i+1,j,k)) - &
1839 (max(ry_avg,0.0)*segment%grad_normal(j-1,2,k) + min(ry_avg,0.0)*segment%grad_normal(j,2,k))) / &
1843 obc%rx_normal(i,j,k) = segment%rx_normal(i,j,k)
1844 obc%ry_normal(i,j,k) = segment%ry_normal(i,j,k)
1845 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
1846 elseif (segment%gradient)
then
1847 segment%normal_vel(i,j,k) = u_new(i+1,j,k)
1849 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then
1851 if (dhdt*dhdx <= 0.0)
then
1852 tau = segment%Velocity_nudging_timescale_in
1854 tau = segment%Velocity_nudging_timescale_out
1856 gamma_2 = dt / (tau + dt)
1857 segment%normal_vel(i,j,k) = (1 - gamma_2) * segment%normal_vel(i,j,k) + &
1858 gamma_2 * segment%nudged_normal_vel(i,j,k)
1861 if (segment%radiation_tan .or. segment%radiation_grad)
then
1863 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1865 rx_tangential(i,segment%HI%JsdB,k) = segment%rx_normal(i,segment%HI%jsd,k)
1866 rx_tangential(i,segment%HI%JedB,k) = segment%rx_normal(i,segment%HI%jed,k)
1867 do j=segment%HI%JsdB+1,segment%HI%JedB-1
1868 rx_tangential(i,j,k) = 0.5*(segment%rx_normal(i,j,k) + segment%rx_normal(i,j+1,k))
1871 if (segment%radiation_tan)
then
1872 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1873 rx_avg = rx_tangential(i,j,k)
1874 segment%tangential_vel(i,j,k) = (v_new(i+1,j,k) + rx_avg*v_new(i+2,j,k)) / (1.0+rx_avg)
1877 if (segment%nudged_tan)
then
1878 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1880 if (rx_tangential(i,j,k) <= 0.0)
then
1881 tau = segment%Velocity_nudging_timescale_in
1883 tau = segment%Velocity_nudging_timescale_out
1885 gamma_2 = dt / (tau + dt)
1886 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
1887 gamma_2 * segment%nudged_tangential_vel(i,j,k)
1890 if (segment%radiation_grad)
then
1891 js_obc = max(segment%HI%JsdB,g%jsd+1)
1892 je_obc = min(segment%HI%JedB,g%jed-1)
1893 do k=1,nz ;
do j=js_obc,je_obc
1894 rx_avg = rx_tangential(i,j,k)
1904 segment%tangential_grad(i,j,k) = ((v_new(i+2,j,k) - v_new(i+1,j,k))*g%IdxBu(i+1,j) + &
1905 rx_avg*(v_new(i+3,j,k) - v_new(i+2,j,k))*g%IdxBu(i+2,j)) / (1.0+rx_avg)
1908 if (segment%nudged_grad)
then
1909 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1911 if (rx_tangential(i,j,k) <= 0.0)
then
1912 tau = segment%Velocity_nudging_timescale_in
1914 tau = segment%Velocity_nudging_timescale_out
1916 gamma_2 = dt / (tau + dt)
1917 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
1918 gamma_2 * segment%nudged_tangential_grad(i,j,k)
1921 deallocate(rx_tangential)
1923 if (segment%oblique_tan .or. segment%oblique_grad)
then
1925 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1926 allocate(ry_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1927 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
1929 rx_tangential(i,segment%HI%JsdB,k) = segment%rx_normal(i,segment%HI%jsd,k)
1930 rx_tangential(i,segment%HI%JedB,k) = segment%rx_normal(i,segment%HI%jed,k)
1931 ry_tangential(i,segment%HI%JsdB,k) = segment%ry_normal(i,segment%HI%jsd,k)
1932 ry_tangential(i,segment%HI%JedB,k) = segment%ry_normal(i,segment%HI%jed,k)
1933 cff_tangential(i,segment%HI%JsdB,k) = segment%cff_normal(i,segment%HI%jsd,k)
1934 cff_tangential(i,segment%HI%JedB,k) = segment%cff_normal(i,segment%HI%jed,k)
1935 do j=segment%HI%JsdB+1,segment%HI%JedB-1
1936 rx_tangential(i,j,k) = 0.5*(segment%rx_normal(i,j,k) + segment%rx_normal(i,j+1,k))
1937 ry_tangential(i,j,k) = 0.5*(segment%ry_normal(i,j,k) + segment%ry_normal(i,j+1,k))
1938 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i,j+1,k))
1941 if (segment%oblique_tan)
then
1942 do k=1,nz ;
do j=segment%HI%JsdB+1,segment%HI%JedB-1
1943 rx_avg = rx_tangential(i,j,k)
1944 ry_avg = ry_tangential(i,j,k)
1945 cff_avg = cff_tangential(i,j,k)
1946 segment%tangential_vel(i,j,k) = ((cff_avg*v_new(i+1,j,k) + rx_avg*v_new(i+2,j,k)) - &
1947 (max(ry_avg,0.0)*segment%grad_tan(j,2,k) + min(ry_avg,0.0)*segment%grad_tan(j+1,2,k))) / &
1951 if (segment%nudged_tan)
then
1952 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1954 if (rx_tangential(i,j,k) <= 0.0)
then
1955 tau = segment%Velocity_nudging_timescale_in
1957 tau = segment%Velocity_nudging_timescale_out
1959 gamma_2 = dt / (tau + dt)
1960 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
1961 gamma_2 * segment%nudged_tangential_vel(i,j,k)
1964 if (segment%oblique_grad)
then
1965 js_obc = max(segment%HI%JsdB,g%jsd+1)
1966 je_obc = min(segment%HI%JedB,g%jed-1)
1967 do k=1,nz ;
do j=segment%HI%JsdB+1,segment%HI%JedB-1
1968 rx_avg = rx_tangential(i,j,k)
1969 ry_avg = ry_tangential(i,j,k)
1970 cff_avg = cff_tangential(i,j,k)
1971 segment%tangential_grad(i,j,k) = ((cff_avg*(v_new(i+2,j,k) - v_new(i+1,j,k))*g%IdxBu(i+1,j) &
1972 + rx_avg*(v_new(i+3,j,k) - v_new(i+2,j,k))*g%IdxBu(i+2,j)) - &
1973 (max(ry_avg,0.0)*segment%grad_gradient(j,2,k) + min(ry_avg,0.0)*segment%grad_gradient(j+1,2,k))) / &
1977 if (segment%nudged_grad)
then
1978 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
1980 if (rx_tangential(i,j,k) <= 0.0)
then
1981 tau = segment%Velocity_nudging_timescale_in
1983 tau = segment%Velocity_nudging_timescale_out
1985 gamma_2 = dt / (tau + dt)
1986 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
1987 gamma_2 * segment%nudged_tangential_grad(i,j,k)
1990 deallocate(rx_tangential)
1991 deallocate(ry_tangential)
1992 deallocate(cff_tangential)
1996 if (segment%direction == obc_direction_n)
then
1998 if (j<g%HI%JscB) cycle
1999 do k=1,nz ;
do i=segment%HI%isd,segment%HI%ied
2000 if (segment%radiation)
then
2001 dhdt = v_old(i,j-1,k)-v_new(i,j-1,k)
2002 dhdy = v_new(i,j-1,k)-v_new(i,j-2,k)
2004 if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max)
2005 ry_avg = (1.0-gamma_v)*segment%ry_normal(i,j,k) + gamma_v*ry_new
2006 segment%ry_normal(i,j,k) = ry_avg
2010 segment%normal_vel(i,j,k) = (v_new(i,j,k) + ry_avg*v_new(i,j-1,k)) / (1.0+ry_avg)
2013 obc%ry_normal(i,j,k) = segment%ry_normal(i,j,k)
2014 elseif (segment%oblique)
then
2015 dhdt = v_old(i,j-1,k)-v_new(i,j-1,k)
2016 dhdy = v_new(i,j-1,k)-v_new(i,j-2,k)
2017 segment%ry_normal(i,j,k) = ry_avg
2018 if (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) > 0.0)
then
2019 dhdx = segment%grad_normal(i-1,1,k)
2020 elseif (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) == 0.0)
then
2023 dhdx = segment%grad_normal(i,1,k)
2025 if (dhdt*dhdy < 0.0) dhdt = 0.0
2027 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2028 rx_new = min(cff,max(dhdt*dhdx,-cff))
2029 rx_avg = (1.0-gamma_u)*segment%rx_normal(i,j,k) + gamma_u*rx_new
2030 ry_avg = (1.0-gamma_u)*segment%ry_normal(i,j,k) + gamma_u*ry_new
2031 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
2032 segment%rx_normal(i,j,k) = rx_avg
2033 segment%ry_normal(i,j,k) = ry_avg
2034 segment%cff_normal(i,j,k) = cff_avg
2035 segment%normal_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + ry_avg*v_new(i,j-1,k)) - &
2036 (max(rx_avg,0.0)*segment%grad_normal(i-1,2,k) + min(rx_avg,0.0)*segment%grad_normal(i,2,k))) / &
2040 obc%rx_normal(i,j,k) = segment%rx_normal(i,j,k)
2041 obc%ry_normal(i,j,k) = segment%ry_normal(i,j,k)
2042 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
2043 elseif (segment%gradient)
then
2044 segment%normal_vel(i,j,k) = v_new(i,j-1,k)
2046 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then
2048 if (dhdt*dhdy <= 0.0)
then
2049 tau = segment%Velocity_nudging_timescale_in
2051 tau = segment%Velocity_nudging_timescale_out
2053 gamma_2 = dt / (tau + dt)
2054 segment%normal_vel(i,j,k) = (1 - gamma_2) * segment%normal_vel(i,j,k) + &
2055 gamma_2 * segment%nudged_normal_vel(i,j,k)
2058 if (segment%radiation_tan .or. segment%radiation_grad)
then
2060 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2062 rx_tangential(segment%HI%IsdB,j,k) = segment%ry_normal(segment%HI%isd,j,k)
2063 rx_tangential(segment%HI%IedB,j,k) = segment%ry_normal(segment%HI%ied,j,k)
2064 do i=segment%HI%IsdB+1,segment%HI%IedB-1
2065 rx_tangential(i,j,k) = 0.5*(segment%ry_normal(i,j,k) + segment%ry_normal(i+1,j,k))
2068 if (segment%radiation_tan)
then
2069 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2070 rx_avg = rx_tangential(i,j,k)
2071 segment%tangential_vel(i,j,k) = (u_new(i,j,k) + rx_avg*u_new(i,j-1,k)) / (1.0+rx_avg)
2074 if (segment%nudged_tan)
then
2075 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2077 if (rx_tangential(i,j,k) <= 0.0)
then
2078 tau = segment%Velocity_nudging_timescale_in
2080 tau = segment%Velocity_nudging_timescale_out
2082 gamma_2 = dt / (tau + dt)
2083 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
2084 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2087 if (segment%radiation_grad)
then
2088 is_obc = max(segment%HI%IsdB,g%isd+1)
2089 ie_obc = min(segment%HI%IedB,g%ied-1)
2090 do k=1,nz ;
do i=is_obc,ie_obc
2091 rx_avg = rx_tangential(i,j,k)
2101 segment%tangential_grad(i,j,k) = ((u_new(i,j,k) - u_new(i,j-1,k))*g%IdyBu(i,j-1) + &
2102 rx_avg*(u_new(i,j-1,k) - u_new(i,j-2,k))*g%IdyBu(i,j-2)) / (1.0+rx_avg)
2105 if (segment%nudged_grad)
then
2106 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2108 if (ry_tangential(i,j,k) <= 0.0)
then
2109 tau = segment%Velocity_nudging_timescale_in
2111 tau = segment%Velocity_nudging_timescale_out
2113 gamma_2 = dt / (tau + dt)
2114 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
2115 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2118 deallocate(rx_tangential)
2120 if (segment%oblique_tan .or. segment%oblique_grad)
then
2122 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2123 allocate(ry_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2124 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2126 rx_tangential(segment%HI%IsdB,j,k) = segment%ry_normal(segment%HI%isd,j,k)
2127 rx_tangential(segment%HI%IedB,j,k) = segment%ry_normal(segment%HI%ied,j,k)
2128 ry_tangential(segment%HI%IsdB,j,k) = segment%rx_normal(segment%HI%isd,j,k)
2129 ry_tangential(segment%HI%IedB,j,k) = segment%rx_normal(segment%HI%ied,j,k)
2130 cff_tangential(segment%HI%IsdB,j,k) = segment%cff_normal(segment%HI%isd,j,k)
2131 cff_tangential(segment%HI%IedB,j,k) = segment%cff_normal(segment%HI%ied,j,k)
2132 do i=segment%HI%IsdB+1,segment%HI%IedB-1
2133 rx_tangential(i,j,k) = 0.5*(segment%ry_normal(i,j,k) + segment%ry_normal(i+1,j,k))
2134 ry_tangential(i,j,k) = 0.5*(segment%rx_normal(i,j,k) + segment%rx_normal(i+1,j,k))
2135 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i+1,j,k))
2138 if (segment%oblique_tan)
then
2139 do k=1,nz ;
do i=segment%HI%IsdB+1,segment%HI%IedB-1
2140 rx_avg = rx_tangential(i,j,k)
2141 ry_avg = ry_tangential(i,j,k)
2142 cff_avg = cff_tangential(i,j,k)
2143 segment%tangential_vel(i,j,k) = ((cff_avg*u_new(i,j,k) + rx_avg*u_new(i,j-1,k)) - &
2144 (max(ry_avg,0.0)*segment%grad_tan(i,2,k) + min(ry_avg,0.0)*segment%grad_tan(i+1,2,k))) / &
2148 if (segment%nudged_tan)
then
2149 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2151 if (ry_tangential(i,j,k) <= 0.0)
then
2152 tau = segment%Velocity_nudging_timescale_in
2154 tau = segment%Velocity_nudging_timescale_out
2156 gamma_2 = dt / (tau + dt)
2157 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
2158 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2161 if (segment%oblique_grad)
then
2162 is_obc = max(segment%HI%IsdB,g%isd+1)
2163 ie_obc = min(segment%HI%IedB,g%ied-1)
2164 do k=1,nz ;
do i=segment%HI%IsdB+1,segment%HI%IedB-1
2165 rx_avg = rx_tangential(i,j,k)
2166 ry_avg = ry_tangential(i,j,k)
2167 cff_avg = cff_tangential(i,j,k)
2168 segment%tangential_grad(i,j,k) = ((cff_avg*(u_new(i,j,k) - u_new(i,j-1,k))*g%IdyBu(i,j-1) &
2169 + rx_avg*(u_new(i,j-1,k) - u_new(i,j-2,k))*g%IdyBu(i,j-2)) - &
2170 (max(ry_avg,0.0)*segment%grad_gradient(i,2,k) + min(ry_avg,0.0)*segment%grad_gradient(i+1,2,k))) / &
2174 if (segment%nudged_grad)
then
2175 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2177 if (ry_tangential(i,j,k) <= 0.0)
then
2178 tau = segment%Velocity_nudging_timescale_in
2180 tau = segment%Velocity_nudging_timescale_out
2182 gamma_2 = dt / (tau + dt)
2183 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
2184 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2187 deallocate(rx_tangential)
2188 deallocate(ry_tangential)
2189 deallocate(cff_tangential)
2193 if (segment%direction == obc_direction_s)
then
2195 if (j>g%HI%JecB) cycle
2196 do k=1,nz ;
do i=segment%HI%isd,segment%HI%ied
2197 if (segment%radiation)
then
2198 dhdt = v_old(i,j+1,k)-v_new(i,j+1,k)
2199 dhdy = v_new(i,j+1,k)-v_new(i,j+2,k)
2201 if (dhdt*dhdy > 0.0) ry_new = min( (dhdt/dhdy), ry_max)
2202 ry_avg = (1.0-gamma_v)*segment%ry_normal(i,j,k) + gamma_v*ry_new
2203 segment%ry_normal(i,j,k) = ry_avg
2207 segment%normal_vel(i,j,k) = (v_new(i,j,k) + ry_avg*v_new(i,j+1,k)) / (1.0+ry_avg)
2210 obc%ry_normal(i,j,k) = segment%ry_normal(i,j,k)
2211 elseif (segment%oblique)
then
2212 dhdt = v_old(i,j+1,k)-v_new(i,j+1,k)
2213 dhdy = v_new(i,j+1,k)-v_new(i,j+2,k)
2214 if (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) > 0.0)
then
2215 dhdx = segment%grad_normal(i-1,1,k)
2216 elseif (dhdt*(segment%grad_normal(i,1,k) + segment%grad_normal(i-1,1,k)) == 0.0)
then
2219 dhdx = segment%grad_normal(i,1,k)
2221 if (dhdt*dhdy < 0.0) dhdt = 0.0
2223 cff_new = max(dhdx*dhdx + dhdy*dhdy, eps)
2224 rx_new = min(cff,max(dhdt*dhdx,-cff))
2225 rx_avg = (1.0-gamma_u)*segment%rx_normal(i,j,k) + gamma_u*rx_new
2226 ry_avg = (1.0-gamma_u)*segment%ry_normal(i,j,k) + gamma_u*ry_new
2227 cff_avg = (1.0-gamma_u)*segment%cff_normal(i,j,k) + gamma_u*cff_new
2228 segment%rx_normal(i,j,k) = rx_avg
2229 segment%ry_normal(i,j,k) = ry_avg
2230 segment%cff_normal(i,j,k) = cff_avg
2231 segment%normal_vel(i,j,k) = ((cff_avg*v_new(i,j,k) + ry_avg*v_new(i,j+1,k)) - &
2232 (max(rx_avg,0.0)*segment%grad_normal(i-1,2,k) + min(rx_avg,0.0)*segment%grad_normal(i,2,k))) / &
2236 obc%rx_normal(i,j,k) = segment%rx_normal(i,j,k)
2237 obc%ry_normal(i,j,k) = segment%ry_normal(i,j,k)
2238 obc%cff_normal(i,j,k) = segment%cff_normal(i,j,k)
2239 elseif (segment%gradient)
then
2240 segment%normal_vel(i,j,k) = v_new(i,j+1,k)
2242 if ((segment%radiation .or. segment%oblique) .and. segment%nudged)
then
2244 if (dhdt*dhdy <= 0.0)
then
2245 tau = segment%Velocity_nudging_timescale_in
2247 tau = segment%Velocity_nudging_timescale_out
2249 gamma_2 = dt / (tau + dt)
2250 segment%normal_vel(i,j,k) = (1 - gamma_2) * segment%normal_vel(i,j,k) + &
2251 gamma_2 * segment%nudged_normal_vel(i,j,k)
2254 if (segment%radiation_tan .or. segment%radiation_grad)
then
2256 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2258 rx_tangential(segment%HI%IsdB,j,k) = segment%ry_normal(segment%HI%isd,j,k)
2259 rx_tangential(segment%HI%IedB,j,k) = segment%ry_normal(segment%HI%ied,j,k)
2260 do i=segment%HI%IsdB+1,segment%HI%IedB-1
2261 rx_tangential(i,j,k) = 0.5*(segment%ry_normal(i,j,k) + segment%ry_normal(i+1,j,k))
2264 if (segment%radiation_tan)
then
2265 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2266 rx_avg = rx_tangential(i,j,k)
2267 segment%tangential_vel(i,j,k) = (u_new(i,j+1,k) + rx_avg*u_new(i,j+2,k)) / (1.0+rx_avg)
2270 if (segment%nudged_tan)
then
2271 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2273 if (rx_tangential(i,j,k) <= 0.0)
then
2274 tau = segment%Velocity_nudging_timescale_in
2276 tau = segment%Velocity_nudging_timescale_out
2278 gamma_2 = dt / (tau + dt)
2279 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
2280 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2283 if (segment%radiation_grad)
then
2284 is_obc = max(segment%HI%IsdB,g%isd+1)
2285 ie_obc = min(segment%HI%IedB,g%ied-1)
2286 do k=1,nz ;
do i=is_obc,ie_obc
2287 rx_avg = rx_tangential(i,j,k)
2297 segment%tangential_grad(i,j,k) = ((u_new(i,j+2,k) - u_new(i,j+1,k))*g%IdyBu(i,j+1) + &
2298 rx_avg*(u_new(i,j+3,k) - u_new(i,j+2,k))*g%IdyBu(i,j+2)) / (1.0+rx_avg)
2301 if (segment%nudged_grad)
then
2302 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2304 if (ry_tangential(i,j,k) <= 0.0)
then
2305 tau = segment%Velocity_nudging_timescale_in
2307 tau = segment%Velocity_nudging_timescale_out
2309 gamma_2 = dt / (tau + dt)
2310 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
2311 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2314 deallocate(rx_tangential)
2316 if (segment%oblique_tan .or. segment%oblique_grad)
then
2318 allocate(rx_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2319 allocate(ry_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2320 allocate(cff_tangential(segment%HI%IsdB:segment%HI%IedB,segment%HI%JsdB:segment%HI%JedB,nz))
2322 rx_tangential(segment%HI%IsdB,j,k) = segment%ry_normal(segment%HI%isd,j,k)
2323 rx_tangential(segment%HI%IedB,j,k) = segment%ry_normal(segment%HI%ied,j,k)
2324 ry_tangential(segment%HI%IsdB,j,k) = segment%rx_normal(segment%HI%isd,j,k)
2325 ry_tangential(segment%HI%IedB,j,k) = segment%rx_normal(segment%HI%ied,j,k)
2326 cff_tangential(segment%HI%IsdB,j,k) = segment%cff_normal(segment%HI%isd,j,k)
2327 cff_tangential(segment%HI%IedB,j,k) = segment%cff_normal(segment%HI%ied,j,k)
2328 do i=segment%HI%IsdB+1,segment%HI%IedB-1
2329 rx_tangential(i,j,k) = 0.5*(segment%ry_normal(i,j,k) + segment%ry_normal(i+1,j,k))
2330 ry_tangential(i,j,k) = 0.5*(segment%rx_normal(i,j,k) + segment%rx_normal(i+1,j,k))
2331 cff_tangential(i,j,k) = 0.5*(segment%cff_normal(i,j,k) + segment%cff_normal(i+1,j,k))
2334 if (segment%oblique_tan)
then
2335 do k=1,nz ;
do i=segment%HI%IsdB+1,segment%HI%IedB-1
2336 rx_avg = rx_tangential(i,j,k)
2337 ry_avg = ry_tangential(i,j,k)
2338 cff_avg = cff_tangential(i,j,k)
2339 segment%tangential_vel(i,j,k) = ((cff_avg*u_new(i,j+1,k) + rx_avg*u_new(i,j+2,k)) - &
2340 (max(ry_avg,0.0)*segment%grad_tan(i,2,k) + min(ry_avg,0.0)*segment%grad_tan(i+1,2,k))) / &
2344 if (segment%nudged_tan)
then
2345 do k=1,nz ;
do i=segment%HI%IsdB,segment%HI%IedB
2347 if (ry_tangential(i,j,k) <= 0.0)
then
2348 tau = segment%Velocity_nudging_timescale_in
2350 tau = segment%Velocity_nudging_timescale_out
2352 gamma_2 = dt / (tau + dt)
2353 segment%tangential_vel(i,j,k) = (1 - gamma_2) * segment%tangential_vel(i,j,k) + &
2354 gamma_2 * segment%nudged_tangential_vel(i,j,k)
2357 if (segment%oblique_grad)
then
2358 is_obc = max(segment%HI%IsdB,g%isd+1)
2359 ie_obc = min(segment%HI%IedB,g%ied-1)
2360 do k=1,nz ;
do i=segment%HI%IsdB+1,segment%HI%IedB-1
2361 rx_avg = rx_tangential(i,j,k)
2362 ry_avg = ry_tangential(i,j,k)
2363 cff_avg = cff_tangential(i,j,k)
2364 segment%tangential_grad(i,j,k) = ((cff_avg*(u_new(i,j+2,k) - u_new(i,j+1,k))*g%IdyBu(i,j+1) &
2365 + rx_avg*(u_new(i,j+3,k) - u_new(i,j+2,k))*g%IdyBu(i,j+2)) - &
2366 (max(ry_avg,0.0)*segment%grad_gradient(i,2,k) + min(ry_avg,0.0)*segment%grad_gradient(i+1,2,k))) / &
2370 if (segment%nudged_grad)
then
2371 do k=1,nz ;
do j=segment%HI%JsdB,segment%HI%JedB
2373 if (ry_tangential(i,j,k) <= 0.0)
then
2374 tau = segment%Velocity_nudging_timescale_in
2376 tau = segment%Velocity_nudging_timescale_out
2378 gamma_2 = dt / (tau + dt)
2379 segment%tangential_grad(i,j,k) = (1 - gamma_2) * segment%tangential_grad(i,j,k) + &
2380 gamma_2 * segment%nudged_tangential_grad(i,j,k)
2383 deallocate(rx_tangential)
2384 deallocate(ry_tangential)
2385 deallocate(cff_tangential)
2391 call open_boundary_apply_normal_flow(obc, g, u_new, v_new)
2393 call pass_vector(u_new, v_new, g%Domain, clock=id_clock_pass)