subroutine da_psichi_to_uv(psi, chi, coefx,coefy, u, v)

   !---------------------------------------------------------------------------
   !  Purpose: Calculate wind components u and v from psi and chi.
   !
   !  Method:  u = coefx * (-dpsi/dy + dchi/dx)
   !           v = coefy * ( dpsi/dx + dchi/dy)
   !
   !  Assumptions: Unstaggered grid.
   !               Lateral boundary conditions - dpsi/dn, dchi/dn = 0 (FCT)
   !               grid size may or may not be equal
   !
   !    Updated for Analysis on Arakawa-C grid
   !    Author: Syed RH Rizvi,  MMM/ESSL/NCAR,  Date: 10/22/2008
   !----------------------------------------------------------------------------

   implicit none
   
   real, intent(inout) :: psi(ims:ime,jms:jme,kms:kme) ! Stream function
   real, intent(inout) :: chi(ims:ime,jms:jme,kms:kme) ! Velocity potential
   real, intent(in)    :: coefx(ims:ime,jms:jme)       ! Multiplicative coeff.
   real, intent(in)    :: coefy(ims:ime,jms:jme)       ! Multiplicative coeff.
   real, intent(out)   :: u(ims:ime,jms:jme,kms:kme)   ! u wind comp.
   real, intent(out)   :: v(ims:ime,jms:jme,kms:kme)   ! v wind comp.


   integer            :: i, j, k                      ! Loop counters.
   integer            :: is, ie                       ! 1st dim. end points.
   integer            :: js, je                       ! 2nd dim. end points.
#ifdef A2C
   real               :: psi1, psi2, chi1, chi2
#endif

   if (trace_use) call da_trace_entry("da_psichi_to_uv")

   !---------------------------------------------------------------------------
   ! [1.0] For Global application, set Wast-Eest Periodic boundary
   !---------------------------------------------------------------------------

   ! Computation to check for edge of domain:
   is = its
   ie = ite
   js = jts
   je = jte
   if (jts == jds) js = jds+1
   if (jte == jde) je = jde-1

#ifdef A2C
   if (fg_format == fg_format_kma_global) then
#else
   if (global) then 
#endif
      call da_set_boundary_3d(psi)
      call da_set_boundary_3d(chi)
   else
      if (its == ids) is = ids+1
      if (ite == ide) ie = ide-1
   end if


   !$OMP PARALLEL DO &
   !$OMP PRIVATE (i, j, k)
   do k = kts, kte
      !------------------------------------------------------------------------
      !  [2.0] Compute u, v at interior points (2nd order central finite diffs):
      !------------------------------------------------------------------------
      do j = js, je
         do i = is, ie
#ifdef A2C
          if ((fg_format==fg_format_wrf_arw_regional) .or. &
              (fg_format==fg_format_wrf_arw_global  ) ) then
            psi1 = 0.5*(psi(i,j+1,k) + psi(i-1,j+1,k)) 
            psi2 = 0.5*(psi(i,j-1,k) + psi(i-1,j-1,k)) 
            chi1 = 2.0*(chi(i,j,k) - chi(i-1,j,k))
            u(i,j,k) = coefy(i,j)*(psi2 - psi1 ) + coefx(i,j)*chi1 

            psi1 = 0.5*(psi(i+1,j-1,k) + psi(i+1,j,k)) 
            psi2 = 0.5*(psi(i-1,j-1,k) + psi(i-1,j,k)) 
            chi1 = 2.0*(chi(i,j,k) - chi(i,j-1,k))
            v(i,j,k) = coefx(i,j)*(psi1 - psi2 ) + coefy(i,j)*chi1                              

          else if (fg_format == fg_format_kma_global) then
#endif
            u(i,j,k) = -coefy(i,j)*(psi(i  ,j+1,k) - psi(i  ,j-1,k)) + &
                        coefx(i,j)*(chi(i+1,j  ,k) - chi(i-1,j  ,k)) 

            v(i,j,k) = coefx(i,j)*(psi(i+1,j  ,k) - psi(i-1,j  ,k))  + &
                       coefy(i,j)*(chi(i  ,j+1,k) - chi(i  ,j-1,k)) 
#ifdef A2C
          else
             write(unit=message(1),fmt='(A,I5)') &
                   "Wrong choice for fg_format = ",fg_format
             call da_error(__FILE__,__LINE__,message(1:1))
          end if
#endif
         end do
      end do

#ifdef A2C
      if (fg_format == fg_format_kma_global) cycle
#else
      if (global) cycle
#endif
      !------------------------------------------------------------------------
      ! [3.0] Compute u, v at domain boundaries:
      !------------------------------------------------------------------------

      ! [3.1] Western boundaries:

      if (its == ids) then
         i = its
         do j = js, je
#ifdef A2C
          if ((fg_format==fg_format_wrf_arw_regional) .or. &
              (fg_format==fg_format_wrf_arw_global  ) ) then
          u(i,j,k) = 2.0*u(i+1,j,k) - u(i+2,j,k)
          v(i,j,k) = 2.0*v(i+1,j,k) - v(i+2,j,k)
          
          else if (fg_format == fg_format_kma_global) then
#endif
            u(i,j,k) = -coefy(i,j)*(psi(i  ,j+1,k) - psi(i,j-1,k)) + &
                        coefx(i,j)*(chi(i+2,j  ,k) - chi(i,j  ,k))  

            v(i,j,k) = coefx(i,j)*(psi(i+2,j  ,k) - psi(i,j  ,k))  + &
                       coefy(i,j)*(chi(i  ,j+1,k) - chi(i,j-1,k)) 
#ifdef A2C
          else
             write(unit=message(1),fmt='(A,I5)') &
                   "Wrong choice for fg_format = ",fg_format
             call da_error(__FILE__,__LINE__,message(1:1))
          end if
#endif
         end do
      end if

      ! [3.2] Eastern boundaries:

      if (ite == ide) then
         i = ite
         do j = js, je
#ifdef A2C
         if ((fg_format==fg_format_wrf_arw_regional) .or. &
              (fg_format==fg_format_wrf_arw_global  ) ) then
          u(i,j,k) = 2.0*u(i-1,j,k) - u(i-2,j,k)
          v(i,j,k) = 2.0*v(i-1,j,k) - v(i-2,j,k)
         else if (fg_format == fg_format_kma_global) then
#endif
            u(i,j,k) = -coefy(i,j)*(psi(i,j+1,k) - psi(i  ,j-1,k)) + &
                        coefx(i,j)*(chi(i,j  ,k) - chi(i-2,j  ,k)) 

            v(i,j,k) = coefx(i,j)*(psi(i,j  ,k) - psi(i-2,j  ,k))+ &
                       coefy(i,j)*(chi(i,j+1,k) - chi(i  ,j-1,k)) 
#ifdef A2C
         else
             write(unit=message(1),fmt='(A,I5)') &
                   "Wrong choice for fg_format = ",fg_format
             call da_error(__FILE__,__LINE__,message(1:1))
         end if
#endif
         end do
      end if

      ! [3.3] Southern boundaries:

      if (jts == jds) then
         j = jts
         do i = is, ie
#ifdef A2C
          if ((fg_format==fg_format_wrf_arw_regional) .or. &
              (fg_format==fg_format_wrf_arw_global  ) ) then
          u(i,j,k) = 2.0*u(i,j+1,k) - u(i,j+2,k)
          v(i,j,k) = 2.0*v(i,j+1,k) - v(i,j+2,k)
         else if (fg_format == fg_format_kma_global) then
#endif
            u(i,j,k) = -coefy(i,j)*(psi(i  ,j+2,k) - psi(i  ,j,k)) + &
                        coefx(i,j)*(chi(i+1,j  ,k) - chi(i-1,j,k))  

            v(i,j,k) = coefx(i,j)*(psi(i+1,j  ,k) - psi(i-1,j,k)) + &
                       coefy(i,j)*(chi(i  ,j+2,k) - chi(i  ,j,k)) 
#ifdef A2C
         else
             write(unit=message(1),fmt='(A,I5)') &
                   "Wrong choice for fg_format = ",fg_format
             call da_error(__FILE__,__LINE__,message(1:1))
         end if
#endif
         end do
      end if

      ! [3.4] Northern boundaries:

      if (jte == jde) then
         j = jte
         do i = is, ie
#ifdef A2C
         if ((fg_format==fg_format_wrf_arw_regional) .or. &
             (fg_format==fg_format_wrf_arw_global  ) ) then
          u(i,j,k) = 2.0*u(i,j-1,k) - u(i,j-2,k)
          v(i,j,k) = 2.0*v(i,j-1,k) - v(i,j-2,k)
         else if (fg_format == fg_format_kma_global) then
#endif
            u(i,j,k) = -coefy(i,j)*(psi(i  ,j,k) - psi(i  ,j-2,k)) + &
                        coefx(i,j)*(chi(i+1,j,k) - chi(i-1,j  ,k))  

            v(i,j,k) = coefx(i,j)*(psi(i+1,j,k) - psi(i-1,j  ,k))+ &
                       coefy(i,j)*(chi(i  ,j,k) - chi(i  ,j-2,k)) 
#ifdef A2C
         else
             write(unit=message(1),fmt='(A,I5)') &
                   "Wrong choice for fg_format = ",fg_format
             call da_error(__FILE__,__LINE__,message(1:1))
         end if
#endif
         end do
      end if

      !------------------------------------------------------------------------
      ! [4.0] Corner points (assume average of surrounding points - poor?):
      !------------------------------------------------------------------------

      ! [4.1] Bottom-left point:

      if (its == ids .AND. jts == jds) then
         u(its,jts,k) = 0.5 * (u(its+1,jts,k) + u(its,jts+1,k))
         v(its,jts,k) = 0.5 * (v(its+1,jts,k) + v(its,jts+1,k))
      end if

      ! [4.2] Top-left point:

#ifdef A2C
      if (its == ids .AND. jte == jde) then
         u(its,jte,k) = 0.5 * (u(its+1,jte,k) + u(its,jte-1,k))
         v(its,jte,k) = 0.5 * (v(its+1,jte,k) + v(its,jte-1,k))
      end if
#else
      if (ite == ide .AND. jts == jds) then
         u(ite,jts,k) = 0.5 * (u(ite-1,jts,k) + u(ite,jts+1,k))
         v(ite,jts,k) = 0.5 * (v(ite-1,jts,k) + v(ite,jts+1,k))
      endif
#endif

      ! [4.3] Bottom-right point:

#ifdef A2C
      if (ite == ide .AND. jts == jds) then
         u(ite,jts,k) = 0.5 * (u(ite-1,jts,k) + u(ite,jts+1,k))
         v(ite,jts,k) = 0.5 * (v(ite-1,jts,k) + v(ite,jts+1,k))
      end if
#else
      if (its == ids .AND. jte == jde) then
         u(its,jte,k) = 0.5 * (u(its+1,jte,k) + u(its,jte-1,k))
         v(its,jte,k) = 0.5 * (v(its+1,jte,k) + v(its,jte-1,k))
      end if
#endif

      ! [4.4] Top-right point:

      if (ite == ide .AND. jte == jde) then
         u(ite,jte,k) = 0.5 * (u(ite-1,jte,k) + u(ite,jte-1,k))
         v(ite,jte,k) = 0.5 * (v(ite-1,jte,k) + v(ite,jte-1,k))
      end if
   end do
   !$OMP END PARALLEL DO

   !---------------------------------------------------------------------------
   ! [5.0] For Global application, set Wast-Eest Periodic boundary
   !---------------------------------------------------------------------------

#ifdef A2C
   if (fg_format == fg_format_kma_global) then
#else
   if (global) then
#endif
      call da_set_boundary_3d(u)
      call da_set_boundary_3d(v)
   end if

   if (trace_use) call da_trace_exit("da_psichi_to_uv")

end subroutine da_psichi_to_uv