wrfv2_fire/var/da/da_recursive_filter/da_transform_through_rf_adj.inc

   1 subroutine da_transform_through_rf_adj(grid, mz,rf_alpha, val, field)
   2
   3    !---------------------------------------------------------------------------
   4    ! Purpose: Control variable transform through the recursive filter.
   5    !
   6    ! Method: 1) Transform to nondimensional v_hat space.
   7    !         2) Perform recursive filter in non-dimensional space (ds = 1).
   8    !         3) Scale recursive filter output.
   9    !         4) Transform filtered field to dimensional space.
  10    !         5) Optionally scale by background error.
  11    !---------------------------------------------------------------------------
  12
  13    implicit none
  14
  15    type(domain),     intent(inout) :: grid
  16    integer,          intent(in)    :: mz                             ! Vertical truncation.
  17    real,             intent(in)    :: rf_alpha(mz)                   ! RF scale parameter.
  18    real,             intent(in)    :: val(jds:jde,mz)                ! Error standard deviation.
  19    real,             intent(inout) :: field(ims:ime,jms:jme,kms:kme) ! Field to be transformed.
  20
  21    integer :: rf_passes_over_two    ! rf_passes / 2
  22    integer :: i, j, m, n, pass, ij  ! Loop counters.
  23    real    :: p_x(ims:ime,jms:jme)  ! sqrt(Grid box area).
  24    real    :: val_j(grid%xp%jtsy:grid%xp%jtey)
  25    real    :: val_i(grid%xp%itsx:grid%xp%itex)
  26
  27    !-------------------------------------------------------------------------
  28    ! [1.0]: Initialise:
  29    !-------------------------------------------------------------------------
  30
  31    if (trace_use_dull) call da_trace_entry("da_transform_through_rf_adj")
  32
  33    rf_passes_over_two = rf_passes / 2
  34
  35    ! [1.1] Define inner product (square root of grid box area):
  36    !$OMP PARALLEL DO &
  37    !$OMP PRIVATE (ij, i, j)
  38    do ij = 1 , grid%num_tiles
  39       do j = grid%j_start(ij), grid%j_end(ij)
  40          do i = its, ite
  41             p_x(i,j) = sqrt(grid%xb%grid_box_area(i,j))
  42          end do
  43       end do
  44    end do
  45    !$OMP END PARALLEL DO
  46
  47    !-------------------------------------------------------------------------
  48    ! [4.0]: Perform 1D recursive filter in y-direction:
  49    !-------------------------------------------------------------------------
  50
  51    ! [4.3] Optionally scale by background error:
  52
  53    ! be_s % val = Gridpoint standard deviation - only required for
  54    ! vert_corr = vert_corr_1 as scaling is performed in vertical transform
  55    ! for vert_corr = vert_corr_2:
  56
  57    if (vert_corr == vert_corr_1) then
  58       do m = 1, mz
  59          do i = its, ite
  60             field(i,jts:jte,m) = field(i,jts:jte,m) * val(jts:jte,m)
  61          end do
  62       end do
  63    end if
  64
  65    ! [4.2] Transform filtered field to dimensional space:
  66
  67    !$OMP PARALLEL DO &
  68    !$OMP PRIVATE (ij ,m, j, i)
  69    do ij = 1 , grid%num_tiles
  70       do m = 1, mz
  71          do j = grid%j_start(ij), grid%j_end(ij)
  72             do i = its, ite
  73                grid%xp%v1z(i,j,m) = field(i,j,m) * p_x(i,j)
  74             end do
  75          end do
  76       end do
  77    end do
  78    !$OMP END PARALLEL DO
  79
  80    ! [4.1] Apply (i',j',k -> i',j,k') (grid%xp%v1z -> grid%xp%v1y)
  81    ! convert vertical column to y-stripe
  82
  83    call da_transpose_z2y (grid)
  84
  85    !-------------------------------------------------------------------------
  86    ! [3.0]: Perform 1D recursive filter in y-direction:
  87    !-------------------------------------------------------------------------
  88
  89    !  [3.2] Apply 1D filter in y direction:
  90
  91    n=grid%xp%jtey-grid%xp%jtsy+1
  92    !$OMP PARALLEL DO &
  93    !$OMP PRIVATE (m, i, val_j, pass, j)
  94    do m = grid%xp%ktsy, min(grid%xp%ktey, mz)
  95       do i = grid%xp%itsy, grid%xp%itey
  96          do j = grid%xp%jtsy, grid%xp%jtey
  97             val_j(j) = grid%xp%v1y(i,j,m)
  98          end do
  99          do pass = rf_passes_over_two, 1, -1
 100             call da_recursive_filter_1d_adj(pass, rf_alpha(m), val_j, n)
 101          end do
 102          do j = grid%xp%jtsy, grid%xp%jtey
 103             grid%xp%v1y(i,j,m) = val_j(j)
 104          end do
 105       end do
 106    end do
 107    !$OMP END PARALLEL DO
 108
 109    ! [3.1] Apply (i',j,k' -> i,j',k') (grid%xp%v1y -> grid%xp%v1x)
 110    ! convert from y-stripe to x-stripe
 111
 112    call da_transpose_y2x (grid)
 113
 114    !-------------------------------------------------------------------------
 115    ! [2.0]: Perform 1D recursive filter in x-direction:
 116    !-------------------------------------------------------------------------
 117
 118    ! [2.2] Apply 1D filter in x direction:
 119
 120    n = grid%xp%itex-grid%xp%itsx+1
 121
 122    !$OMP PARALLEL DO &
 123    !$OMP PRIVATE ( m, j, pass, i, val_i)
 124    do m = grid%xp%ktsx, min(grid%xp%ktex,mz)
 125       do j = grid%xp%jtsx, grid%xp%jtex
 126          do i = grid%xp%itsx, grid%xp%itex
 127             val_i(i) = grid%xp%v1x(i,j,m)
 128          end do
 129          do pass = rf_passes_over_two, 1, -1
 130             call da_recursive_filter_1d_adj(pass, rf_alpha(m), val_i, n)
 131          end do
 132          do i = grid%xp%itsx, grid%xp%itex
 133             grid%xp%v1x(i,j,m) = val_i(i)
 134          end do
 135       end do
 136    end do
 137    !$OMP END PARALLEL DO
 138
 139    ! [2.1] Apply (i,j',k' -> i',j',k) (grid%xp%v1x -> grid%xp%v1z)
 140    ! convert from x-stripe to vertical column
 141
 142    call da_transpose_x2z (grid)
 143
 144    !-------------------------------------------------------------------------
 145    ! [1.0]: Initialise:
 146    !-------------------------------------------------------------------------
 147
 148    ! [1.2] Transform to nondimensional v_hat space:
 149
 150    !$OMP PARALLEL DO &
 151    !$OMP PRIVATE (ij ,m, i, j)
 152    do ij = 1 , grid%num_tiles
 153       do m = 1, mz
 154          do j = grid%j_start(ij), grid%j_end(ij)
 155             do i = its, ite
 156                field(i,j,m) = grid%xp%v1z(i,j,m) / p_x(i,j)
 157             end do
 158          end do
 159       end do
 160    end do
 161    !$OMP END PARALLEL DO
 162
 163    if (trace_use_dull) call da_trace_exit("da_transform_through_rf_adj")
 164
 165 end subroutine da_transform_through_rf_adj
 166
 167