subroutine da_psichi_to_uv(psi, chi, coefx,coefy, u, v) 3,11
!---------------------------------------------------------------------------
! 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