subroutine da_balance_equation_lin(grid, xbx, u, v, phi_b)

   !---------------------------------------------------------------------------
   !  Purpose: Calculates balanced mass increment phi_b from wind increments.
   !
   !  Method:  Solve grad**2 Phi_b = - div[ k x rho f u + rho (u.grad ) u ] by
   !           1) Calculate RHS of balance equation in gridpt space.
   !           2) Solve Del**2 phi_b = RHS for phi_b using double (FCT).
   !
   !---------------------------------------------------------------------------

   implicit none

   type(domain),   intent(inout) :: grid
   type(xbx_type), intent(in)    :: xbx                            ! Header & non-gridded vars. 
   real,           intent(in)    :: u(ims:ime,jms:jme,kms:kme)     ! u wind comp.
   real,           intent(in)    :: v(ims:ime,jms:jme,kms:kme)     ! v wind comp.
   real,           intent(out)   :: phi_b(ims:ime,jms:jme,kms:kme) ! Balanced mass increment.

   integer :: i, j, k                 ! Loop counters.
   integer :: is, ie                  ! 1st dim. end points.
   integer :: js, je                  ! 2nd dim. end points.

   real    :: coefx(ims:ime,jms:jme)  ! Multiplicative coefficient.
   real    :: coefy(ims:ime,jms:jme)  ! Multiplicative coefficient.
   real    :: term_x(ims:ime,jms:jme) ! Balance eqn x term
   real    :: term_y(ims:ime,jms:jme) ! Balance eqn y term
   
   real    :: del2phi_b(ims:ime,jms:jme,kms:kme)  ! Del**2 phi_b/M**2

   if (trace_use) call da_trace_entry("da_balance_equation_lin")

   !---------------------------------------------------------------------------
   ! [1.0] Initialise iand set Multipilactive constants
   !---------------------------------------------------------------------------

   ! 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

   if (fg_format == fg_format_kma_global) then
      coefx = grid%xb%coefx
      coefy = grid%xb%coefy
   else if( fg_format == fg_format_wrf_arw_regional) then
      coefx = grid%xb%coefz
      coefy = coefx   
   else if (fg_format == fg_format_wrf_arw_global) then
      write (unit=message(1),fmt='(A,I3)') ' needs work for fg_format_wrf_arw_global  = ',fg_format
      call da_error(__FILE__,__LINE__,message(1:1))
   else if (fg_format == fg_format_wrf_nmm_regional) then
      write (unit=message(1),fmt='(A,I3)') ' needs work for fg_format_wrf_nmm_regional = ',fg_format
      call da_error(__FILE__,__LINE__,message(1:1))
   else
      write(unit=message(1),fmt='(A,I3)') 'Wrong FG_FORMAT = ',fg_format
      call da_error(__FILE__,__LINE__,message(1:1))
   end if

   ! [1.1] Multiplicative coefficent for conversion RHS->Del**2 phi_b/M**2:

   do k = kts, kte
      term_x(ims:ime,jms:jme) = 0.0
      term_y(ims:ime,jms:jme) = 0.0

      !---------------------------------------------------------------------------
      ! [2.0] Calculate RHS of balance equation in gridpt space:
      !---------------------------------------------------------------------------

      ! [2.1] Include geostrophic terms in balance eqn if requested:

      if (balance_type == balance_geo .OR. balance_type == balance_geocyc ) then
         call da_balance_geoterm_lin (grid%xb % cori, grid%xb % rho(:,:,k), u(:,:,k), v(:,:,k), &
            term_x, term_y)
      end if
      
      ! [2.2] Include cyclostrophic terms in balance eqn if requested:

      if (balance_type == balance_cyc .OR. balance_type == balance_geocyc ) then
         call da_balance_cycloterm_lin (grid%xb % rho(:,:,k), grid%xb % u(:,:,k),   &
            grid%xb % v(:,:,k), u(:,:,k), v(:,:,k), grid%xb % coefx(:,:), grid%xb % coefy(:,:), &
            term_x(:,:), term_y(:,:))
      end if
      
      ! [2.3] Take divergence to get Del**2 phi_b/M**2:

      do j = js, je
         do i = is, ie
          del2phi_b(i,j,k) = -coefx(i,j)*( term_x(i+1,j) - term_x(i-1,j)) &
                             -coefy(i,j)*( term_y(i,j+1) - term_y(i,j-1)) 
         end do
      end do

      ! [2.4] Del**2 Phi_b  boundary conditions:

      if (.not. global .and. its == ids ) del2phi_b(its,js:je,k) = 0.0
      if (.not. global .and. ite == ide ) del2phi_b(ite,js:je,k) = 0.0
      if (jts == jds ) del2phi_b(is:ie,jts,k) = 0.0
      if (jte == jde ) del2phi_b(is:ie,jte,k) = 0.0

      if (.not. global .and. (its == ids .and. jts == jds) ) del2phi_b(its,jts,k) = 0.0
      if (.not. global .and. (its == ids .and. jte == jde) ) del2phi_b(its,jte,k) = 0.0
      if (.not. global .and. (ite == ide .and. jts == jds) ) del2phi_b(ite,jts,k) = 0.0
      if (.not. global .and. (ite == ide .and. jte == jde) ) del2phi_b(ite,jte,k) = 0.0     
   end do

   !------------------------------------------------------------------------------
   !  [3.0] Solve Del**2 phi_b = RHS for phi_b:
   !------------------------------------------------------------------------------

   call da_solve_poissoneqn_fst(grid, xbx, del2phi_b, phi_b) 

   if (trace_use) call da_trace_exit("da_balance_equation_lin")

end subroutine da_balance_equation_lin