var/da/da_test/da_check_vtoy_adjoint.inc

   1 subroutine da_check_vtoy_adjoint(cv_size,grid, config_flags, vp, vv, xbx, be, ep, 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    integer,                    intent(in)    :: cv_size
  12    type(grid_config_rec_type), intent(inout) :: config_flags
  13    type(domain),               intent(inout) :: grid
  14    type(vp_type),              intent(inout) :: vv    ! Grdipt/EOF CV.
  15    type(vp_type),              intent(inout) :: vp    ! Grdipt/level CV.
  16    type(xbx_type),             intent(inout) :: xbx   ! Header & non-gridded vars.
  17    type(be_type),              intent(in)    :: be    ! background error structure.
  18    type(ep_type),              intent(in)    :: ep     ! ensemble perturbation structure.
  19    type(iv_type),              intent(inout) :: iv    ! ob. increment vector.
  20    type(y_type),               intent(inout) :: y     ! y = h (xa)
  21
  22    real    :: cv(1:cv_size)          ! Test control variable.
  23    real    :: cv_2(1:cv_size)
  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_dull) call da_trace_entry("da_check_vtoy_adjoint")
  30
  31    write(unit=stdout, fmt='(/a/a)') &
  32         '       da_check_vtoy_adjoint:',&
  33         '---------------------------------------'
  34
  35    call random_number(cv(:))
  36    cv(:) = cv(:) - 0.5
  37
  38    cv_2(1:cv_size) = cv(1:cv_size)
  39
  40    call da_zero_x(grid%xa)
  41    call da_zero_vp_type(vp)
  42    call da_zero_vp_type(vv)
  43
  44    call da_transform_vtoy(cv_size, be, ep, cv, iv, vp, vv, xbx, y, grid, config_flags, vp, vv)
  45
  46    !-------------------------------------------------------------------------
  47    ! [3.0] Calculate LHS of adjoint test equation and
  48    !        Rescale input to adjoint routine :
  49    !-------------------------------------------------------------------------
  50
  51    call da_get_y_lhs_value(iv, y, partial_lhs, pertile_lhs, adj_sum_lhs)
  52
  53    cv = 0.0
  54
  55    ! WHY?
  56    ! call da_zero_vp_type(vp)
  57    ! call da_zero_vp_type(vv)
  58    ! call da_zero_x(grid%xa)
  59
  60    call da_transform_vtoy_adj(cv_size, be, ep, cv, iv, vp, vv, xbx, y, grid, &
  61       config_flags, .true., vp, vv)
  62
  63    adj_rhs = sum(cv(1:cv_size) * cv_2(1:cv_size))
  64
  65    !-------------------------------------------------------------------------
  66    !      Print output:
  67    !-------------------------------------------------------------------------
  68
  69 #ifdef DM_PARALLEL
  70    if (global) then
  71       adj_sum_rhs = adj_rhs
  72    else
  73       call mpi_allreduce(adj_rhs, adj_sum_rhs, 1, true_mpi_real, mpi_sum, &
  74                        comm, ierr)
  75    end if
  76 #else
  77    adj_sum_rhs = adj_rhs
  78    adj_sum_lhs = partial_lhs
  79 #endif
  80
  81 #ifdef DM_PARALLEL
  82    if (rootproc) then
  83       write(unit=stdout, fmt='(A,1pe22.14)') &
  84       'Whole Domain  < y, y     > = ', adj_sum_lhs
  85       write(unit=stdout, fmt='(A,1pe22.14)') &
  86          'Whole Domain  < v, v_adj > = ', adj_sum_rhs
  87    end if
  88 #else
  89    write(unit=stdout, fmt='(A,1pe22.14)') &
  90       'Whole Domain  < y, y     > = ', adj_sum_lhs
  91    write(unit=stdout, fmt='(A,1pe22.14)') &
  92       'Whole Domain  < v, v_adj > = ', adj_sum_rhs
  93 #endif
  94
  95    if (trace_use_dull) call da_trace_exit("da_check_vtoy_adjoint")
  96
  97 end subroutine da_check_vtoy_adjoint
  98
  99