var/da/da_dynamics/da_psichi_to_uv.inc

   1 subroutine da_psichi_to_uv(psi, chi, coefx,coefy, u, v)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Calculate wind components u and v from psi and chi.
   5    !
   6    !  Method:  u = coefx * (-dpsi/dy + dchi/dx)
   7    !           v = coefy * ( dpsi/dx + dchi/dy)
   8    !
   9    !  Assumptions: Unstaggered grid.
  10    !               Lateral boundary conditions - dpsi/dn, dchi/dn = 0 (FCT)
  11    !               grid size may or may not be equal
  12    !
  13    !    Updated for Analysis on Arakawa-C grid
  14    !    Author: Syed RH Rizvi,  MMM/ESSL/NCAR,  Date: 10/22/2008
  15    !----------------------------------------------------------------------------
  16
  17    implicit none
  18
  19    real, intent(inout) :: psi(ims:ime,jms:jme,kms:kme) ! Stream function
  20    real, intent(inout) :: chi(ims:ime,jms:jme,kms:kme) ! Velocity potential
  21    real, intent(in)    :: coefx(ims:ime,jms:jme)       ! Multiplicative coeff.
  22    real, intent(in)    :: coefy(ims:ime,jms:jme)       ! Multiplicative coeff.
  23    real, intent(out)   :: u(ims:ime,jms:jme,kms:kme)   ! u wind comp.
  24    real, intent(out)   :: v(ims:ime,jms:jme,kms:kme)   ! v wind comp.
  25
  26
  27    integer            :: i, j, k                      ! Loop counters.
  28    integer            :: is, ie                       ! 1st dim. end points.
  29    integer            :: js, je                       ! 2nd dim. end points.
  30 #ifdef A2C
  31    real               :: psi1, psi2, chi1, chi2
  32 #endif
  33
  34    if (trace_use) call da_trace_entry("da_psichi_to_uv")
  35
  36    !---------------------------------------------------------------------------
  37    ! [1.0] For Global application, set Wast-Eest Periodic boundary
  38    !---------------------------------------------------------------------------
  39
  40    ! Computation to check for edge of domain:
  41    is = its
  42    ie = ite
  43    js = jts
  44    je = jte
  45    if (jts == jds) js = jds+1
  46    if (jte == jde) je = jde-1
  47
  48 #ifdef A2C
  49    if (fg_format == fg_format_kma_global) then
  50 #else
  51    if (global) then
  52 #endif
  53       call da_set_boundary_3d(psi)
  54       call da_set_boundary_3d(chi)
  55    else
  56       if (its == ids) is = ids+1
  57       if (ite == ide) ie = ide-1
  58    end if
  59
  60
  61    !$OMP PARALLEL DO &
  62    !$OMP PRIVATE (i, j, k)
  63    do k = kts, kte
  64       !------------------------------------------------------------------------
  65       !  [2.0] Compute u, v at interior points (2nd order central finite diffs):
  66       !------------------------------------------------------------------------
  67       do j = js, je
  68          do i = is, ie
  69 #ifdef A2C
  70           if ((fg_format==fg_format_wrf_arw_regional) .or. &
  71               (fg_format==fg_format_wrf_arw_global  ) ) then
  72             psi1 = 0.5*(psi(i,j+1,k) + psi(i-1,j+1,k))
  73             psi2 = 0.5*(psi(i,j-1,k) + psi(i-1,j-1,k))
  74             chi1 = 2.0*(chi(i,j,k) - chi(i-1,j,k))
  75             u(i,j,k) = coefy(i,j)*(psi2 - psi1 ) + coefx(i,j)*chi1
  76
  77             psi1 = 0.5*(psi(i+1,j-1,k) + psi(i+1,j,k))
  78             psi2 = 0.5*(psi(i-1,j-1,k) + psi(i-1,j,k))
  79             chi1 = 2.0*(chi(i,j,k) - chi(i,j-1,k))
  80             v(i,j,k) = coefx(i,j)*(psi1 - psi2 ) + coefy(i,j)*chi1
  81
  82           else if (fg_format == fg_format_kma_global) then
  83 #endif
  84             u(i,j,k) = -coefy(i,j)*(psi(i  ,j+1,k) - psi(i  ,j-1,k)) + &
  85                         coefx(i,j)*(chi(i+1,j  ,k) - chi(i-1,j  ,k))
  86
  87             v(i,j,k) = coefx(i,j)*(psi(i+1,j  ,k) - psi(i-1,j  ,k))  + &
  88                        coefy(i,j)*(chi(i  ,j+1,k) - chi(i  ,j-1,k))
  89 #ifdef A2C
  90           else
  91              write(unit=message(1),fmt='(A,I5)') &
  92                    "Wrong choice for fg_format = ",fg_format
  93              call da_error(__FILE__,__LINE__,message(1:1))
  94           end if
  95 #endif
  96          end do
  97       end do
  98
  99 #ifdef A2C
 100       if (fg_format == fg_format_kma_global) cycle
 101 #else
 102       if (global) cycle
 103 #endif
 104       !------------------------------------------------------------------------
 105       ! [3.0] Compute u, v at domain boundaries:
 106       !------------------------------------------------------------------------
 107
 108       ! [3.1] Western boundaries:
 109
 110       if (its == ids) then
 111          i = its
 112          do j = js, je
 113 #ifdef A2C
 114           if ((fg_format==fg_format_wrf_arw_regional) .or. &
 115               (fg_format==fg_format_wrf_arw_global  ) ) then
 116           u(i,j,k) = 2.0*u(i+1,j,k) - u(i+2,j,k)
 117           v(i,j,k) = 2.0*v(i+1,j,k) - v(i+2,j,k)
 118
 119           else if (fg_format == fg_format_kma_global) then
 120 #endif
 121             u(i,j,k) = -coefy(i,j)*(psi(i  ,j+1,k) - psi(i,j-1,k)) + &
 122                         coefx(i,j)*(chi(i+2,j  ,k) - chi(i,j  ,k))
 123
 124             v(i,j,k) = coefx(i,j)*(psi(i+2,j  ,k) - psi(i,j  ,k))  + &
 125                        coefy(i,j)*(chi(i  ,j+1,k) - chi(i,j-1,k))
 126 #ifdef A2C
 127           else
 128              write(unit=message(1),fmt='(A,I5)') &
 129                    "Wrong choice for fg_format = ",fg_format
 130              call da_error(__FILE__,__LINE__,message(1:1))
 131           end if
 132 #endif
 133          end do
 134       end if
 135
 136       ! [3.2] Eastern boundaries:
 137
 138       if (ite == ide) then
 139          i = ite
 140          do j = js, je
 141 #ifdef A2C
 142          if ((fg_format==fg_format_wrf_arw_regional) .or. &
 143               (fg_format==fg_format_wrf_arw_global  ) ) then
 144           u(i,j,k) = 2.0*u(i-1,j,k) - u(i-2,j,k)
 145           v(i,j,k) = 2.0*v(i-1,j,k) - v(i-2,j,k)
 146          else if (fg_format == fg_format_kma_global) then
 147 #endif
 148             u(i,j,k) = -coefy(i,j)*(psi(i,j+1,k) - psi(i  ,j-1,k)) + &
 149                         coefx(i,j)*(chi(i,j  ,k) - chi(i-2,j  ,k))
 150
 151             v(i,j,k) = coefx(i,j)*(psi(i,j  ,k) - psi(i-2,j  ,k))+ &
 152                        coefy(i,j)*(chi(i,j+1,k) - chi(i  ,j-1,k))
 153 #ifdef A2C
 154          else
 155              write(unit=message(1),fmt='(A,I5)') &
 156                    "Wrong choice for fg_format = ",fg_format
 157              call da_error(__FILE__,__LINE__,message(1:1))
 158          end if
 159 #endif
 160          end do
 161       end if
 162
 163       ! [3.3] Southern boundaries:
 164
 165       if (jts == jds) then
 166          j = jts
 167          do i = is, ie
 168 #ifdef A2C
 169           if ((fg_format==fg_format_wrf_arw_regional) .or. &
 170               (fg_format==fg_format_wrf_arw_global  ) ) then
 171           u(i,j,k) = 2.0*u(i,j+1,k) - u(i,j+2,k)
 172           v(i,j,k) = 2.0*v(i,j+1,k) - v(i,j+2,k)
 173          else if (fg_format == fg_format_kma_global) then
 174 #endif
 175             u(i,j,k) = -coefy(i,j)*(psi(i  ,j+2,k) - psi(i  ,j,k)) + &
 176                         coefx(i,j)*(chi(i+1,j  ,k) - chi(i-1,j,k))
 177
 178             v(i,j,k) = coefx(i,j)*(psi(i+1,j  ,k) - psi(i-1,j,k)) + &
 179                        coefy(i,j)*(chi(i  ,j+2,k) - chi(i  ,j,k))
 180 #ifdef A2C
 181          else
 182              write(unit=message(1),fmt='(A,I5)') &
 183                    "Wrong choice for fg_format = ",fg_format
 184              call da_error(__FILE__,__LINE__,message(1:1))
 185          end if
 186 #endif
 187          end do
 188       end if
 189
 190       ! [3.4] Northern boundaries:
 191
 192       if (jte == jde) then
 193          j = jte
 194          do i = is, ie
 195 #ifdef A2C
 196          if ((fg_format==fg_format_wrf_arw_regional) .or. &
 197              (fg_format==fg_format_wrf_arw_global  ) ) then
 198           u(i,j,k) = 2.0*u(i,j-1,k) - u(i,j-2,k)
 199           v(i,j,k) = 2.0*v(i,j-1,k) - v(i,j-2,k)
 200          else if (fg_format == fg_format_kma_global) then
 201 #endif
 202             u(i,j,k) = -coefy(i,j)*(psi(i  ,j,k) - psi(i  ,j-2,k)) + &
 203                         coefx(i,j)*(chi(i+1,j,k) - chi(i-1,j  ,k))
 204
 205             v(i,j,k) = coefx(i,j)*(psi(i+1,j,k) - psi(i-1,j  ,k))+ &
 206                        coefy(i,j)*(chi(i  ,j,k) - chi(i  ,j-2,k))
 207 #ifdef A2C
 208          else
 209              write(unit=message(1),fmt='(A,I5)') &
 210                    "Wrong choice for fg_format = ",fg_format
 211              call da_error(__FILE__,__LINE__,message(1:1))
 212          end if
 213 #endif
 214          end do
 215       end if
 216
 217       !------------------------------------------------------------------------
 218       ! [4.0] Corner points (assume average of surrounding points - poor?):
 219       !------------------------------------------------------------------------
 220
 221       ! [4.1] Bottom-left point:
 222
 223       if (its == ids .AND. jts == jds) then
 224          u(its,jts,k) = 0.5 * (u(its+1,jts,k) + u(its,jts+1,k))
 225          v(its,jts,k) = 0.5 * (v(its+1,jts,k) + v(its,jts+1,k))
 226       end if
 227
 228       ! [4.2] Top-left point:
 229
 230 #ifdef A2C
 231       if (its == ids .AND. jte == jde) then
 232          u(its,jte,k) = 0.5 * (u(its+1,jte,k) + u(its,jte-1,k))
 233          v(its,jte,k) = 0.5 * (v(its+1,jte,k) + v(its,jte-1,k))
 234       end if
 235 #else
 236       if (ite == ide .AND. jts == jds) then
 237          u(ite,jts,k) = 0.5 * (u(ite-1,jts,k) + u(ite,jts+1,k))
 238          v(ite,jts,k) = 0.5 * (v(ite-1,jts,k) + v(ite,jts+1,k))
 239       endif
 240 #endif
 241
 242       ! [4.3] Bottom-right point:
 243
 244 #ifdef A2C
 245       if (ite == ide .AND. jts == jds) then
 246          u(ite,jts,k) = 0.5 * (u(ite-1,jts,k) + u(ite,jts+1,k))
 247          v(ite,jts,k) = 0.5 * (v(ite-1,jts,k) + v(ite,jts+1,k))
 248       end if
 249 #else
 250       if (its == ids .AND. jte == jde) then
 251          u(its,jte,k) = 0.5 * (u(its+1,jte,k) + u(its,jte-1,k))
 252          v(its,jte,k) = 0.5 * (v(its+1,jte,k) + v(its,jte-1,k))
 253       end if
 254 #endif
 255
 256       ! [4.4] Top-right point:
 257
 258       if (ite == ide .AND. jte == jde) then
 259          u(ite,jte,k) = 0.5 * (u(ite-1,jte,k) + u(ite,jte-1,k))
 260          v(ite,jte,k) = 0.5 * (v(ite-1,jte,k) + v(ite,jte-1,k))
 261       end if
 262    end do
 263    !$OMP END PARALLEL DO
 264
 265    !---------------------------------------------------------------------------
 266    ! [5.0] For Global application, set Wast-Eest Periodic boundary
 267    !---------------------------------------------------------------------------
 268
 269 #ifdef A2C
 270    if (fg_format == fg_format_kma_global) then
 271 #else
 272    if (global) then
 273 #endif
 274       call da_set_boundary_3d(u)
 275       call da_set_boundary_3d(v)
 276    end if
 277
 278    if (trace_use) call da_trace_exit("da_psichi_to_uv")
 279
 280 end subroutine da_psichi_to_uv
 281
 282