var/da/da_test/da_test_vtoy_transform.inc

   1 subroutine da_test_vtoy_transform(grid, config_flags, vp, vv, xbx, be, iv, y)
   2
   3    !---------------------------------------------------------------------------
   4    !  Purpose: Perform V To Y Adjoint transform test
   5    !
   6    !  Method:  Perform adjoint test on complete transform: <y,y> = <v_adj,v>.
   7    !---------------------------------------------------------------------------
   8
   9    implicit none
  10
  11    type(grid_config_rec_type), intent(inout) :: config_flags
  12    type(domain),               intent(inout) :: grid
  13    type(vp_type),              intent(inout) :: vv    ! Grdipt/EOF CV.
  14    type(vp_type),              intent(inout) :: vp    ! Grdipt/level CV.
  15    type(xbx_type),             intent(in)    :: xbx   ! Header & non-gridded vars.
  16    type(be_type),              intent(in)    :: be    ! background error structure.
  17    type(iv_type),              intent(inout) :: iv    ! ob. increment vector.
  18    type(y_type),               intent(inout) :: y     ! y = h (xa)
  19
  20
  21    real    :: cv(1:cv_size)          ! Test control variable.
  22    real    :: cv_2(1:cv_size)
  23    real    :: field(ims:ime, jms:jme)
  24    real    :: adj_sum_lhs               ! < y, y >
  25    real    :: adj_rhs,adj_sum_rhs       ! < v, v_adj >
  26    real    :: partial_lhs   ! < y, y >
  27    real    :: pertile_lhs   ! < y, y >
  28
  29    if (trace_use) call da_trace_entry("da_test_vtoy_transform")
  30
  31    write(unit=stdout, fmt='(a)') '       da_test_vtoy_transform:'
  32    write(unit=stdout, fmt='(a)') '---------------------------------------'
  33
  34    call random_number(cv(:))
  35    cv(:) = cv(:) - 0.5
  36
  37    cv_2(1:cv_size) = cv(1:cv_size)
  38
  39    call da_zero_x(grid%xa)
  40    call da_zero_vp_type(vp)
  41    call da_zero_vp_type(vv)
  42
  43    call da_transform_vtoy(be, cv, iv, vp, vv, xbx, y, grid, config_flags )
  44
  45    !-------------------------------------------------------------------------
  46    ! [3.0] Calculate LHS of adjoint test equation and
  47    !       Rescale input to adjoint routine :
  48    !-------------------------------------------------------------------------
  49
  50    call da_get_y_lhs_value( iv, y, partial_lhs, pertile_lhs, adj_sum_lhs )
  51
  52    cv = 0.0
  53
  54    ! call da_zero_vp_type(vp)
  55    ! call da_zero_vp_type(vv)
  56    ! call da_zero_x(grid%xa)
  57
  58    call da_transform_vtoy_adj(be, cv, iv, vp, vv, xbx, y, grid, config_flags, .true. )
  59
  60    adj_rhs = sum( cv(1:cv_size) * cv_2(1:cv_size) )
  61    !-------------------------------------------------------------------------
  62    ! Print output:
  63    !-------------------------------------------------------------------------
  64    ! FIX? with wrf_dm_sum_real
  65
  66 #ifdef DM_PARALLEL
  67   if( global ) then
  68      adj_sum_rhs = adj_rhs
  69   else
  70      call mpi_allreduce( adj_rhs, adj_sum_rhs, 1, true_mpi_real, mpi_sum, comm, ierr)
  71   end if
  72 #else
  73    adj_sum_rhs = adj_rhs
  74    adj_sum_lhs = partial_lhs
  75 #endif
  76
  77 #ifdef DM_PARALLEL
  78    if ( rootproc ) then
  79       write(unit=stdout, fmt='(A,1pe22.14)') &
  80          'Whole Domain  < y, y     > = ', adj_sum_lhs
  81       write(unit=stdout, fmt='(A,1pe22.14)') &
  82          'Whole Domain  < v, v_adj > = ', adj_sum_rhs
  83    end if
  84 #else
  85    write(unit=stdout, fmt='(A,1pe22.14)') &
  86       'Whole Domain  < y, y     > = ', adj_sum_lhs
  87    write(unit=stdout, fmt='(A,1pe22.14)') &
  88       'Whole Domain  < v, v_adj > = ', adj_sum_rhs
  89 #endif
  90
  91    if (trace_use) call da_trace_exit("da_test_vtoy_transform")
  92
  93 end subroutine da_test_vtoy_transform
  94
  95