<HTML> <BODY BGCOLOR=#ccccdd LINK=#0000aa VLINK=#0000ff ALINK=#ff0000 ><BASE TARGET="bottom_target"><PRE>
<A NAME='DA_BALANCE_CYCLOTERM'><A href='../../html_code/dynamics/da_balance_cycloterm.inc.html#DA_BALANCE_CYCLOTERM' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_balance_cycloterm (xb, k, term_x, term_y),2
!---------------------------------------------------------------------------
! Purpose: Calculates cyclostrophic term in balance equation.
!
! Method: Term is rho (u.grad) u on a single level.
!---------------------------------------------------------------------------
implicit none
type(xb_type), intent(in) :: xb ! First guess structure.
integer, intent(in) :: k ! Model level.
real, intent(inout) :: term_x(:,:) ! x component of term.
real, intent(inout) :: term_y(:,:) ! y component of term.
integer :: i, j ! Loop counters.
integer :: is, ie ! 1st dim. end points.
integer :: js, je ! 2nd dim. end points.
real :: data(ims:ime,jms:jme) ! Temporary storage.
real :: varb
if (trace_use) call da_trace_entry
("da_balance_cycloterm")
!---------------------------------------------------------------------------
! [1.0] Initialise:
!---------------------------------------------------------------------------
! Computation to check for edge of domain:
is = its; ie = ite; js = jts; je = jte
if (.not. global .and. its==ids) is = ids+1
if (.not. global .and. ite==ide) ie = ide-1
if (jts==jds) js = jds+1
if (jte==jde) je = jde-1
!---------------------------------------------------------------------------
! [2.0] Calculate term_x = rho M (u.du/dx + v.du/dy):
!---------------------------------------------------------------------------
! [2.1] Interior points:
do j = js, je
do i = is, ie
data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)*(xb%u(i+1,j,k) - &
xb%u(i-1,j,k)) + xb%v(i,j,k) * xb%coefy(i,j)*(xb%u(i,j+1,k) - &
xb%u(i,j-1,k))
end do
end do
if (.NOT. global) then ! For global only interior points
! [2.2] Bottom boundaries:
if (its==ids) then
i = its
do j = js, je
varb = -3.0*xb%u(i,j,k)+4.0*xb%u(i+1,j,k)-xb%u(i+2,j,k)
data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * varb + &
xb%v(i,j,k) * xb%coefy(i,j) * (xb%u(i,j+1,k) - xb%u(i,j-1,k))
end do
end if
! [2.3] Top boundaries:
if (ite==ide) then
i = ite
do j = js, je
varb = 3.0*xb%u(i,j,k)-4.0*xb%u(i-1,j,k)+xb%u(i-2,j,k)
data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * varb + &
xb%v(i,j,k) * xb%coefy(i,j) * (xb%u(i,j+1,k) - xb%u(i,j-1,k))
end do
end if
! [2.4] Left boundaries:
if (jts==jds) then
j = jts
do i = is, ie
varb = -3.0*xb%u(i,j,k)+4.0*xb%u(i,j+1,k)-xb%u(i,j+2,k)
data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * (xb%u(i+1,j,k) - &
xb%u(i-1,j,k)) + xb%v(i,j,k) * xb%coefy(i,j) * varb
end do
end if
! [2.5] Right boundaries:
if (jte==jde) then
j = jte
do i = is, ie
varb = 3.0*xb%u(i,j,k)-4.0*xb%u(i,j-1,k)+xb%u(i,j-2,k)
data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * (xb%u(i+1,j,k) - &
xb%u(i-1,j,k)) + xb%v(i,j,k) * xb%coefy(i,j) * varb
end do
end if
! [2.6] Corner points:
if (its==ids .AND. jts==jds) then
data(its,jts) = 0.5 * (data(its,jts+1) + data(its+1,jts))
end if
if (ite==ide .AND. jts==jds) then
data(ite,jts) = 0.5 * (data(ite-1,jts) + data(ite,jts+1))
end if
if (its==ids .AND. jte==jde) then
data(its,jde) = 0.5 * (data(its,jde-1) + data(its+1,jde))
end if
if (ite==ide .AND. jte==jde) then
data(ite,jte) = 0.5 * (data(ite-1,jte) + data(ite,jte-1))
end if
end if ! not global
! [2.7] Multiply by rho and add to term_x:
term_x(its:ite,jts:jte) = xb%rho(its:ite,jts:jte,k)*data(its:ite,jts:jte) + term_x(its:ite,jts:jte)
!---------------------------------------------------------------------------
! [3.0] Calculate term_y = rho M (u.dv/dx + v.dv/dy):
!---------------------------------------------------------------------------
! [3.1] Interior points:
do j = js, je
do i = is, ie
data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)*(xb%v(i+1,j,k) - xb%v(i-1,j,k)) + &
xb%v(i,j,k) * xb%coefy(i,j)*(xb%v(i,j+1,k) - xb%v(i,j-1,k))
end do
end do
if (.NOT. global) then ! For global only interior points
! [3.2] Bottom boundaries:
if (its==ids) then
i = its
do j = js, je
varb = -3.0*xb%v(i,j,k)+4.0*xb%v(i+1,j,k)-xb%v(i+2,j,k)
data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)* varb + &
xb%v(i,j,k) * xb%coefy(i,j)* (xb%v(i,j+1,k) - xb%v(i,j-1,k))
end do
end if
! [3.3] Top boundaries:
if (ite==ide) then
i = ite
do j = js, je
varb = 3.0*xb%v(i,j,k)-4.0*xb%v(i-1,j,k)+xb%v(i-2,j,k)
data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)* varb + &
xb%v(i,j,k) * xb%coefy(i,j)* (xb%v(i,j+1,k) - xb%v(i,j-1,k))
end do
end if
! [3.4] Left boundaries:
if (jts==jds) then
j = jts
do i = is, ie
varb = -3.0*xb%v(i,j,k)+4.0*xb%v(i,j+1,k)-xb%v(i,j+2,k)
data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)* (xb%v(i+1,j,k) - xb%v(i-1,j,k)) + &
xb%v(i,j,k) * xb%coefy(i,j)* varb
end do
end if
! [3.5] Right boundaries:
if (jte==jde) then
j = jte
do i = is, ie
varb = 3.0*xb%v(i,j,k)-4.0*xb%v(i,j+1,k)+xb%v(i,j+2,k)
data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)* (xb%v(i+1,j,k) - xb%v(i-1,j,k)) + &
xb%v(i,j,k) * xb%coefy(i,j)* varb
end do
end if
! [3.6] Corner points:
if (its==ids .AND. jts==jds) then
data(its,jts) = 0.5 * (data(its,jts+1) + data(its+1,jts))
end if
if (ite==ide .AND. jts==jds) then
data(ite,jts) = 0.5 * (data(ite-1,jts) + data(ite,jts+1))
end if
if (its==ids .AND. jte==jde) then
data(its,jde) = 0.5 * (data(its,jde-1) + data(its+1,jde))
end if
if (ite==ide .AND. jte==jde) then
data(ite,jte) = 0.5 * (data(ite-1,jte) + data(ite,jte-1))
end if
end if ! not global
! [3.7] Multiply by rho and add to term_y
term_y(its:ite,jts:jte)=xb%rho(its:ite,jts:jte,k)* data(its:ite,jts:jte) + &
term_y(its:ite,jts:jte)
if (trace_use) call da_trace_exit("da_balance_cycloterm")
end subroutine da_balance_cycloterm