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(in)    :: 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)
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(0,cv_size, be, ep, cv, iv, vp, vv, xbx, y, grid, &
61       config_flags)
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