wrfv2_fire/var/da/da_recursive_filter/da_transform_through_rf.inc

   1 subroutine da_transform_through_rf(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")
  32
  33    rf_passes_over_two = rf_passes / 2
  34
  35    !$OMP PARALLEL DO &
  36    !$OMP PRIVATE ( ij, i, j )
  37    do ij = 1 , grid%num_tiles
  38       ! [1.1] Define inner product (square root of grid box area):
  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    !$OMP PARALLEL DO &
  48    !$OMP PRIVATE ( ij, m, i, j )
  49    do ij = 1 , grid%num_tiles
  50       do m = 1, mz
  51          do j = grid%j_start(ij), grid%j_end(ij)
  52             do i = its, ite
  53                grid%xp%v1z(i,j,m) = 0.0
  54             end do
  55          end do
  56       end do
  57    end do
  58    !$OMP END PARALLEL DO
  59
  60       ! [1.2] Transform to nondimensional v_hat space:
  61
  62    !$OMP PARALLEL DO &
  63    !$OMP PRIVATE ( ij, m, i, j )
  64    do ij = 1 , grid%num_tiles
  65       do m = 1, mz
  66          do j = grid%j_start(ij), grid%j_end(ij)
  67             do i = its, ite
  68                grid%xp%v1z(i,j,m) = field(i,j,m) / p_x(i,j)
  69             end do
  70          end do
  71       end do
  72    end do
  73    !$OMP END PARALLEL DO
  74
  75    !-------------------------------------------------------------------------
  76    ! [2.0]: Perform 1D recursive filter in x-direction:
  77    !-------------------------------------------------------------------------
  78
  79    ! [2.1] Apply (i',j',k -> i,j',k') (grid%xp%v1z -> grid%xp%v1x)
  80    ! convert from vertical column to x-stripe
  81
  82    call da_transpose_z2x (grid)
  83
  84    ! [2.2] Apply 1D filter in x direction:
  85
  86    n = grid%xp%itex - grid%xp%itsx + 1
  87    !$OMP PARALLEL DO &
  88    !$OMP PRIVATE ( m, j, pass, val_i, i )
  89    do m = grid%xp%ktsx, min(grid%xp%ktex,mz)
  90       do j = grid%xp%jtsx, grid%xp%jtex
  91          do i = grid%xp%itsx, grid%xp%itex
  92             val_i(i) = grid%xp%v1x(i,j,m)
  93          end do
  94          do pass = 1, rf_passes_over_two
  95             call da_recursive_filter_1d(pass, rf_alpha(m), val_i, n)
  96          end do
  97          do i = grid%xp%itsx, grid%xp%itex
  98             grid%xp%v1x(i,j,m) = val_i(i)
  99          end do
 100       end do
 101    end do
 102    !$OMP END PARALLEL DO
 103
 104    !-------------------------------------------------------------------------
 105    ! [3.0]: Perform 1D recursive filter in y-direction:
 106    !-------------------------------------------------------------------------
 107
 108    ! [3.1] Apply (i, j' ,k' -> i', j ,k') (grid%xp%v1x -> grid%xp%v1y)
 109    ! convert from vertical column to y-stripe
 110
 111    call da_transpose_x2y (grid)
 112
 113    ! [3.2] Apply 1D filter in y direction:
 114
 115    n = grid%xp%jtey - grid%xp%jtsy + 1
 116    !$OMP PARALLEL DO &
 117    !$OMP PRIVATE ( m, i, pass, val_j, j )
 118    do m = grid%xp%ktsy, min(grid%xp%ktey,mz)
 119       do i = grid%xp%itsy, grid%xp%itey
 120          do j = grid%xp%jtsy, grid%xp%jtey
 121             val_j(j) = grid%xp%v1y(i,j,m)
 122          end do
 123          do pass = 1, rf_passes_over_two
 124             call da_recursive_filter_1d(pass, rf_alpha(m), val_j, n)
 125          end do
 126          do j = grid%xp%jtsy, grid%xp%jtey
 127             grid%xp%v1y(i,j,m) = val_j(j)
 128          end do
 129       end do
 130    end do
 131    !$OMP END PARALLEL DO
 132
 133    !-------------------------------------------------------------------------
 134    ! [4.0]: Perform 1D recursive filter in y-direction:
 135    !-------------------------------------------------------------------------
 136
 137    ! [4.1] Apply (i',j,k' -> i',j',k) (grid%xp%v1y -> grid%xp%v1z)
 138    ! convert from y-stripe to vertical column.
 139
 140    call da_transpose_y2z (grid)
 141
 142    ! [4.2] Transform filtered field to dimensional space:
 143
 144    !$OMP PARALLEL DO &
 145    !$OMP PRIVATE ( ij, m, i, j )
 146    do ij = 1 , grid%num_tiles
 147       do m = 1, mz
 148          do j = grid%j_start(ij), grid%j_end(ij)
 149    !     do j = jts, jte
 150             do i = its, ite
 151                field(i,j,m) = grid%xp%v1z(i,j,m) * p_x(i,j)
 152             end do
 153          end do
 154       end do
 155    end do
 156    !$OMP END PARALLEL DO
 157
 158    ! [4.3] Optionally scale by background error:
 159
 160    ! be_s % val = Gridpoint standard deviation - only required for
 161    ! vert_corr = vert_corr_1 as scaling is performed in vertical transform
 162    ! for vert_corr = vert_corr_2:
 163
 164    if (vert_corr == vert_corr_1) then
 165       do m = 1, mz
 166          do i = its, ite
 167             field(i,jts:jte,m) = field(i,jts:jte,m) * val(jts:jte,m)
 168          end do
 169       end do
 170    end if
 171
 172    if (trace_use_dull) call da_trace_exit("da_transform_through_rf")
 173
 174 end subroutine da_transform_through_rf
 175
 176