subroutine da_balance_cycloterm (xb, k, term_x, term_y) !--------------------------------------------------------------------------- ! 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