wrfv2_fire/var/da/da_dynamics/da_psichi_to_uv_adj.inc

   1 subroutine da_psichi_to_uv_adj(u, v, coefx, coefy, psi, chi)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Adjoint code of da_psichi_to_uv
   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) :: u(ims:ime,jms:jme,kms:kme)   ! u wind comp.
  20    real, intent(inout) :: v(ims:ime,jms:jme,kms:kme)   ! v wind comp.
  21    real, intent(inout) :: psi(ims:ime,jms:jme,kms:kme) ! Stream function
  22    real, intent(inout) :: chi(ims:ime,jms:jme,kms:kme) ! Velocity potential
  23    real, intent(in)    :: coefx(ims:ime,jms:jme)       ! Multiplicative coeff.
  24    real, intent(in)    :: coefy(ims:ime,jms:jme)       ! Multiplicative coeff.
  25
  26    integer :: i, j, k                      ! Loop counters.
  27    integer :: is, ie                       ! 1st dim. end points.
  28    integer :: js, je                       ! 2nd dim. end points.
  29 #ifdef A2C
  30    real               :: psi1, psi2, chi1, chi2
  31 #endif
  32
  33    if (trace_use) call da_trace_entry("da_psichi_to_uv_adj")
  34
  35    !---------------------------------------------------------------------------
  36    ! [5.0] For Global application, set Wast-Eest Periodic boundary
  37    !---------------------------------------------------------------------------
  38
  39 #ifdef A2C
  40    if (fg_format == fg_format_kma_global) then
  41 #else
  42    if (global) then
  43 #endif
  44       call da_set_boundary_3d(u)
  45       call da_set_boundary_3d(v)
  46    end if
  47
  48    !---------------------------------------------------------------------------
  49    ! Computation to check for edge of domain:
  50    !---------------------------------------------------------------------------
  51
  52    is = its-1
  53    ie = ite+1
  54    js = jts-1
  55    je = jte+1
  56    if (jts == jds) js = jds+1
  57    if (jte == jde) je = jde-1
  58
  59 #ifdef A2C
  60    if (.not. (fg_format == fg_format_kma_global) ) then
  61 #else
  62    if (.not. global) then
  63 #endif
  64       if (its == ids) is = ids+1
  65       if (ite == ide) ie = ide-1
  66
  67       !$OMP PARALLEL DO &
  68       !$OMP PRIVATE ( k )
  69       do k = kts, kte
  70
  71          !---------------------------------------------------------------------
  72          ! [4.0] Corner points (assume average of surrounding points - poor?):
  73          !---------------------------------------------------------------------
  74
  75          ! [4.1] Bottom-left point:
  76
  77          if (its == ids .AND. jts == jds) then
  78             u(its+1,jts,k) = u(its+1,jts,k) + 0.5 * u(its,jts,k)
  79             u(its,jts+1,k) = u(its,jts+1,k) + 0.5 * u(its,jts,k)
  80             v(its+1,jts,k) = v(its+1,jts,k) + 0.5 * v(its,jts,k)
  81             v(its,jts+1,k) = v(its,jts+1,k) + 0.5 * v(its,jts,k)
  82          end if
  83
  84         ! [4.2] Top-left point:
  85
  86 #ifdef A2C
  87          if (its == ids .AND. jte == jde) then
  88             u(its+1,jte,k) = u(its+1,jte,k) + 0.5 * u(its,jte,k)
  89             u(its,jte-1,k) = u(its,jte-1,k) + 0.5 * u(its,jte,k)
  90             v(its+1,jte,k) = v(its+1,jte,k) + 0.5 * v(its,jte,k)
  91             v(its,jte-1,k) = v(its,jte-1,k) + 0.5 * v(its,jte,k)
  92          end if
  93 #else
  94          if (ite == ide .AND. jts == jds) then
  95             u(ite-1,jts,k) = u(ite-1,jts,k) + 0.5 * u(ite,jts,k)
  96             u(ite,jts+1,k) = u(ite,jts+1,k) + 0.5 * u(ite,jts,k)
  97             v(ite-1,jts,k) = v(ite-1,jts,k) + 0.5 * v(ite,jts,k)
  98             v(ite,jts+1,k) = v(ite,jts+1,k) + 0.5 * v(ite,jts,k)
  99          end if
 100 #endif
 101
 102          ! [4.3] Bottom-right point:
 103
 104 #ifdef A2C
 105          if (ite == ide .AND. jts == jds) then
 106             u(ite-1,jts,k) = u(ite-1,jts,k) + 0.5 * u(ite,jts,k)
 107             u(ite,jts+1,k) = u(ite,jts+1,k) + 0.5 * u(ite,jts,k)
 108             v(ite-1,jts,k) = v(ite-1,jts,k) + 0.5 * v(ite,jts,k)
 109             v(ite,jts+1,k) = v(ite,jts+1,k) + 0.5 * v(ite,jts,k)
 110          end if
 111
 112 #else
 113          if (its == ids .AND. jte == jde) then
 114             u(its+1,jte,k) = u(its+1,jte,k) + 0.5 * u(its,jte,k)
 115             u(its,jte-1,k) = u(its,jte-1,k) + 0.5 * u(its,jte,k)
 116             v(its+1,jte,k) = v(its+1,jte,k) + 0.5 * v(its,jte,k)
 117             v(its,jte-1,k) = v(its,jte-1,k) + 0.5 * v(its,jte,k)
 118          end if
 119 #endif
 120          ! [4.4] Top-right point:
 121
 122          if (ite == ide .AND. jte == jde) then
 123             u(ite-1,jte,k) = u(ite-1,jte,k) + 0.5 * u(ite,jte,k)
 124             u(ite,jte-1,k) = u(ite,jte-1,k) + 0.5 * u(ite,jte,k)
 125             v(ite-1,jte,k) = v(ite-1,jte,k) + 0.5 * v(ite,jte,k)
 126             v(ite,jte-1,k) = v(ite,jte-1,k) + 0.5 * v(ite,jte,k)
 127          end if
 128       end do
 129       !$OMP END PARALLEL DO
 130    end if
 131
 132
 133    ! [3.0] Compute u, v at domain boundaries:
 134
 135    !$OMP PARALLEL DO &
 136    !$OMP PRIVATE ( i, j, k, message )
 137    do k = kts, kte
 138 #ifdef A2C
 139    if (.not. (fg_format == fg_format_kma_global) ) then
 140 #else
 141    if (.not. global) then
 142 #endif
 143          ! [3.4] Northern boundaries:
 144
 145          if (jte == jde) then
 146             j = jte
 147
 148             do i = ie, is, -1
 149 #ifdef A2C
 150               if ((fg_format==fg_format_wrf_arw_regional) .or. &
 151                   (fg_format==fg_format_wrf_arw_global  ) ) then
 152              u(i,j-1,k) = u(i,j-1,k) + 2.0*u(i,j,k)
 153              v(i,j-1,k) = v(i,j-1,k) + 2.0*v(i,j,k)
 154              u(i,j-2,k) = u(i,j-2,k) -     u(i,j,k)
 155              v(i,j-2,k) = v(i,j-2,k) -     v(i,j,k)
 156
 157              else if (fg_format == fg_format_kma_global) then
 158 #endif
 159                psi(i+1,j,k) = psi(i+1,j,k) + coefx(i,j) * v(i,j,k)
 160                psi(i-1,j,k) = psi(i-1,j,k) - coefx(i,j) * v(i,j,k)
 161                chi(i,j  ,k) = chi(i,j  ,k) + coefy(i,j) * v(i,j,k)
 162                chi(i,j-2,k) = chi(i,j-2,k) - coefy(i,j) * v(i,j,k)
 163
 164                psi(i,j  ,k) = psi(i,j  ,k) - coefy(i,j) * u(i,j,k)
 165                psi(i,j-2,k) = psi(i,j-2,k) + coefy(i,j) * u(i,j,k)
 166                chi(i+1,j,k) = chi(i+1,j,k) + coefx(i,j) * u(i,j,k)
 167                chi(i-1,j,k) = chi(i-1,j,k) - coefx(i,j) * u(i,j,k)
 168 #ifdef A2C
 169              else
 170                  write(unit=message(1),fmt='(A,I5)') &
 171                        "Wrong choice for fg_format = ",fg_format
 172                  call da_error(__FILE__,__LINE__,message(1:1))
 173              end if
 174 #endif
 175             end do
 176          end if
 177          ! [3.3] Southern boundaries:
 178
 179          if (jts == jds) then
 180             j = jts
 181
 182             do i = ie, is, -1
 183
 184 #ifdef A2C
 185              if ((fg_format==fg_format_wrf_arw_regional) .or. &
 186                  (fg_format==fg_format_wrf_arw_global  ) ) then
 187              u(i,j+1,k) = u(i,j+1,k) + 2.0*u(i,j,k)
 188              v(i,j+1,k) = v(i,j+1,k) + 2.0*v(i,j,k)
 189              u(i,j+2,k) = u(i,j+2,k) -     u(i,j,k)
 190              v(i,j+2,k) = v(i,j+2,k) -     v(i,j,k)
 191
 192              else if (fg_format == fg_format_kma_global) then
 193 #endif
 194
 195                psi(i+1,j,k) = psi(i+1,j,k) + coefx(i,j) * v(i,j,k)
 196                psi(i-1,j,k) = psi(i-1,j,k) - coefx(i,j) * v(i,j,k)
 197                chi(i,j+2,k) = chi(i,j+2,k) + coefy(i,j) * v(i,j,k)
 198                chi(i,j  ,k) = chi(i,j  ,k) - coefy(i,j) * v(i,j,k)
 199
 200                psi(i,j+2,k) = psi(i,j+2,k) - coefy(i,j) * u(i,j,k)
 201                psi(i,j  ,k) = psi(i,j  ,k) + coefy(i,j) * u(i,j,k)
 202                chi(i+1,j,k) = chi(i+1,j,k) + coefx(i,j) * u(i,j,k)
 203                chi(i-1,j,k) = chi(i-1,j,k) - coefx(i,j) * u(i,j,k)
 204 #ifdef A2C
 205              else
 206                write(unit=message(1),fmt='(A,I5)') &
 207                      "Wrong choice for fg_format = ",fg_format
 208                call da_error(__FILE__,__LINE__,message(1:1))
 209              end if
 210 #endif
 211             end do
 212          end if
 213
 214          ! [3.2] Eastern boundaries:
 215
 216          if (ite == ide) then
 217             i = ite
 218             do j = je, js, -1
 219 #ifdef A2C
 220              if ((fg_format==fg_format_wrf_arw_regional) .or. &
 221                  (fg_format==fg_format_wrf_arw_global  ) ) then
 222              u(i-1,j,k) = u(i-1,j,k) + 2.0*u(i,j,k)
 223              v(i-1,j,k) = v(i-1,j,k) + 2.0*v(i,j,k)
 224              u(i-2,j,k) = u(i-2,j,k) -     u(i,j,k)
 225              v(i-2,j,k) = v(i-2,j,k) -     v(i,j,k)
 226              else if (fg_format == fg_format_kma_global) then
 227 #endif
 228                psi(i  ,j,k) = psi(i  ,j,k) + coefx(i,j) * v(i,j,k)
 229                psi(i-2,j,k) = psi(i-2,j,k) - coefx(i,j) * v(i,j,k)
 230                chi(i,j+1,k) = chi(i,j+1,k) + coefy(i,j) * v(i,j,k)
 231                chi(i,j-1,k) = chi(i,j-1,k) - coefy(i,j) * v(i,j,k)
 232
 233                psi(i,j+1,k) = psi(i,j+1,k) - coefy(i,j) * u(i,j,k)
 234                psi(i,j-1,k) = psi(i,j-1,k) + coefy(i,j) * u(i,j,k)
 235                chi(i  ,j,k) = chi(i  ,j,k) + coefx(i,j) * u(i,j,k)
 236                chi(i-2,j,k) = chi(i-2,j,k) - coefx(i,j) * u(i,j,k)
 237 #ifdef A2C
 238              else
 239                write(unit=message(1),fmt='(A,I5)') &
 240                      "Wrong choice for fg_format = ",fg_format
 241                call da_error(__FILE__,__LINE__,message(1:1))
 242              end if
 243 #endif
 244             end do
 245          end if
 246
 247          ! [3.1] Western boundaries:
 248          if (its == ids) then
 249             i = its
 250
 251             do j = je, js, -1
 252 #ifdef A2C
 253              if ((fg_format==fg_format_wrf_arw_regional) .or. &
 254                  (fg_format==fg_format_wrf_arw_global  ) ) then
 255              u(i+1,j,k) = u(i+1,j,k) + 2.0*u(i,j,k)
 256              v(i+1,j,k) = v(i+1,j,k) + 2.0*v(i,j,k)
 257              u(i+2,j,k) = u(i+2,j,k) -     u(i,j,k)
 258              v(i+2,j,k) = v(i+2,j,k) -     v(i,j,k)
 259              else if (fg_format == fg_format_kma_global) then
 260 #endif
 261                psi(i+2,j,k) = psi(i+2,j,k) + coefx(i,j) * v(i,j,k)
 262                psi(i  ,j,k) = psi(i  ,j,k) - coefx(i,j) * v(i,j,k)
 263                chi(i,j+1,k) = chi(i,j+1,k) + coefy(i,j) * v(i,j,k)
 264                chi(i,j-1,k) = chi(i,j-1,k) - coefy(i,j) * v(i,j,k)
 265
 266                psi(i,j+1,k) = psi(i,j+1,k) - coefy(i,j) * u(i,j,k)
 267                psi(i,j-1,k) = psi(i,j-1,k) + coefy(i,j) * u(i,j,k)
 268                chi(i+2,j,k) = chi(i+2,j,k) + coefx(i,j) * u(i,j,k)
 269                chi(i,  j,k) = chi(i,  j,k) - coefx(i,j) * u(i,j,k)
 270 #ifdef A2C
 271              else
 272                  write(unit=message(1),fmt='(A,I5)') &
 273                        "Wrong choice for fg_format = ",fg_format
 274                  call da_error(__FILE__,__LINE__,message(1:1))
 275              end if
 276 #endif
 277             end do
 278          end if
 279    end if
 280
 281       !------------------------------------------------------------------------
 282       ! [2.0] Compute u, v at interior points (2nd order central finite diffs):
 283       !------------------------------------------------------------------------
 284       do j = je, js, -1
 285          do i = ie, is, -1
 286 #ifdef A2C
 287             if ((fg_format==fg_format_wrf_arw_regional) .or. &
 288                 (fg_format==fg_format_wrf_arw_global  ) ) then
 289
 290             psi1 = coefx(i,j)*v(i,j,k)
 291             psi2 =-coefx(i,j)*v(i,j,k)
 292             chi1 = coefy(i,j)*v(i,j,k)
 293             chi(i,j,k)   = chi(i,j,k) + 2.0*chi1
 294             chi(i,j-1,k) = chi(i,j-1,k) - 2.0*chi1
 295             psi(i-1,j-1,k) = psi(i-1,j-1,k) +0.5*psi2
 296             psi(i-1,j  ,k) = psi(i-1,j  ,k) +0.5*psi2
 297             psi(i+1,j-1,k) = psi(i+1,j-1,k) + 0.5*psi1
 298             psi(i+1,j  ,k) = psi(i+1,j  ,k) + 0.5*psi1
 299
 300             psi1 = -coefy(i,j)*u(i,j,k)
 301             psi2 =  coefy(i,j)*u(i,j,k)
 302             chi1 =  coefx(i,j)*u(i,j,k)
 303             chi(i,j,k)    = chi(i,j,k)    + 2.0*chi1
 304             chi(i-1,j,k)  = chi(i-1,j,k)  - 2.0*chi1
 305             psi(i,j-1,k)  = psi(i,j-1,k)  + 0.5*psi2
 306             psi(i-1,j-1,k)= psi(i-1,j-1,k)+ 0.5*psi2
 307             psi(i,j+1,k)  = psi(i,j+1,k)  + 0.5*psi1
 308             psi(i-1,j+1,k)= psi(i-1,j+1,k)+ 0.5*psi1
 309           else if (fg_format == fg_format_kma_global) then
 310 #endif
 311             psi(i+1,j,k) = psi(i+1,j,k) + coefx(i,j) * v(i,j,k)
 312             psi(i-1,j,k) = psi(i-1,j,k) - coefx(i,j) * v(i,j,k)
 313             chi(i,j+1,k) = chi(i,j+1,k) + coefy(i,j) * v(i,j,k)
 314             chi(i,j-1,k) = chi(i,j-1,k) - coefy(i,j) * v(i,j,k)
 315
 316             psi(i,j+1,k) = psi(i,j+1,k) - coefy(i,j) * u(i,j,k)
 317             psi(i,j-1,k) = psi(i,j-1,k) + coefy(i,j) * u(i,j,k)
 318             chi(i+1,j,k) = chi(i+1,j,k) + coefx(i,j) * u(i,j,k)
 319             chi(i-1,j,k) = chi(i-1,j,k) - coefx(i,j) * u(i,j,k)
 320 #ifdef A2C
 321           else
 322               write(unit=message(1),fmt='(A,I5)') &
 323                     "Wrong choice for fg_format = ",fg_format
 324               call da_error(__FILE__,__LINE__,message(1:1))
 325           end if
 326 #endif
 327          end do
 328       end do
 329    end do    !  loop over levels
 330    !$OMP END PARALLEL DO
 331
 332    !---------------------------------------------------------------------------
 333    ! [5.0] For Global application, set Wast-Eest Periodic boundary
 334    !---------------------------------------------------------------------------
 335
 336 #ifdef A2C
 337    if( fg_format == fg_format_kma_global ) then
 338 #else
 339    if (global) then
 340 #endif
 341        call da_set_boundary_3d(psi)
 342        call da_set_boundary_3d(chi)
 343    end if
 344
 345    if (trace_use) call da_trace_exit("da_psichi_to_uv_adj")
 346
 347 end subroutine da_psichi_to_uv_adj