da_check_vvtovp_adjoint.inc
References to this file elsewhere.
1 subroutine da_check_vvtovp_adjoint(ne, xb, be, vv, vp)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Test Vv to Vp routine and adjoint for compatibility.
5 !
6 ! Method: Standard adjoint test: < Vp, Vp > = < Vv_adj, Vv >.
7 !---------------------------------------------------------------------------
8
9 implicit none
10
11 integer, intent(in) :: ne ! Ensemble size.
12 type (xb_type), intent(in) :: xb ! first guess (local).
13 type (be_type), intent(in) :: be ! background error structure.
14 type (vp_type), intent(inout) :: vv ! CV(i,j,m).
15 type (vp_type), intent(inout) :: vp ! CV(i,j,k)
16
17 real :: adj_par_lhs ! < x, x >
18 real :: adj_par_rhs ! < v_adj, v >
19 real :: adj_sum_lhs ! < x, x >
20 real :: adj_sum_rhs ! < v_adj, v >
21
22 real :: vv2_v1(ims:ime,jms:jme,kms:kme)
23 real :: vv2_v2(ims:ime,jms:jme,kms:kme)
24 real :: vv2_v3(ims:ime,jms:jme,kms:kme)
25 real :: vv2_v4(ims:ime,jms:jme,kms:kme)
26 real :: vv2_v5(ims:ime,jms:jme,kms:kme)
27 real :: vv2_alpha(ims:ime,jms:jme,1:ne)
28
29 if (trace_use) call da_trace_entry("da_check_vvtovp_adjoint")
30
31 !----------------------------------------------------------------------
32 ! [1.0] Initialise:
33 !----------------------------------------------------------------------
34
35 write(unit=stdout, fmt='(/a/)') 'da_check_vvtovp_adjoint: Test Results:'
36
37 !----------------------------------------------------------------------
38 ! [2.0] Perform Vp = U_v Vv transform:
39 !----------------------------------------------------------------------
40
41 call da_vertical_transform('u', be, &
42 xb % vertical_inner_product, &
43 vv, vp)
44
45 !----------------------------------------------------------------------
46 ! [3.0] Calculate LHS of adjoint test equation:
47 !----------------------------------------------------------------------
48
49 adj_par_lhs = sum(vp % v1(its:ite,jts:jte,kts:kte)**2) * inv_typ_vp1_sumsq &
50 + sum(vp % v2(its:ite,jts:jte,kts:kte)**2) * inv_typ_vp2_sumsq &
51 + sum(vp % v3(its:ite,jts:jte,kts:kte)**2) * inv_typ_vp3_sumsq &
52 + sum(vp % v4(its:ite,jts:jte,kts:kte)**2) * inv_typ_vp4_sumsq &
53 + sum(vp % v5(its:ite,jts:jte,kts:kte)**2) * inv_typ_vp5_sumsq
54
55 if (be % ne > 0) then
56 adj_par_lhs = adj_par_lhs + &
57 sum(vp % alpha(its:ite,jts:jte,1:be%ne)**2) * inv_typ_vpalpha_sumsq
58 end if
59
60 !----------------------------------------------------------------------
61 ! [4.0] Rescale input to adjoint routine:
62 !----------------------------------------------------------------------
63
64 vp % v1(its:ite,jts:jte,kts:kte) = vp % v1(its:ite,jts:jte,kts:kte) * &
65 inv_typ_vp1_sumsq
66 vp % v2(its:ite,jts:jte,kts:kte) = vp % v2(its:ite,jts:jte,kts:kte) * &
67 inv_typ_vp2_sumsq
68 vp % v3(its:ite,jts:jte,kts:kte) = vp % v3(its:ite,jts:jte,kts:kte) * &
69 inv_typ_vp3_sumsq
70 vp % v4(its:ite,jts:jte,kts:kte) = vp % v4(its:ite,jts:jte,kts:kte) * &
71 inv_typ_vp4_sumsq
72 vp % v5(its:ite,jts:jte,kts:kte) = vp % v5(its:ite,jts:jte,kts:kte) * &
73 inv_typ_vp5_sumsq
74
75 if (be % ne > 0) then
76 vp % alpha(its:ite,jts:jte,1:be%ne) = &
77 vp % alpha(its:ite,jts:jte,1:be%ne) * inv_typ_vpalpha_sumsq
78 end if
79
80 !----------------------------------------------------------------------
81 ! [5.0] Perform adjoint operation:
82 !----------------------------------------------------------------------
83
84 vv2_v1(its:ite,jts:jte,1:be%v1%mz) = vv % v1(its:ite,jts:jte,1:be%v1%mz)
85 vv2_v2(its:ite,jts:jte,1:be%v2%mz) = vv % v2(its:ite,jts:jte,1:be%v2%mz)
86 vv2_v3(its:ite,jts:jte,1:be%v3%mz) = vv % v3(its:ite,jts:jte,1:be%v3%mz)
87 vv2_v4(its:ite,jts:jte,1:be%v4%mz) = vv % v4(its:ite,jts:jte,1:be%v4%mz)
88 vv2_v5(its:ite,jts:jte,1:be%v5%mz) = vv % v5(its:ite,jts:jte,1:be%v5%mz)
89
90 if (be % ne > 0) then
91 vv2_alpha(its:ite,jts:jte,1:be%ne) = vv % alpha(its:ite,jts:jte,1:be%ne)
92 end if
93
94 call da_vertical_transform('u_adj', be, &
95 xb % vertical_inner_product, &
96 vv, vp)
97
98 !----------------------------------------------------------------------
99 ! [6.0] Calculate RHS of adjoint test equation:
100 !----------------------------------------------------------------------
101
102 adj_par_rhs = 0.0
103 if (be % v1 % mz > 0) &
104 adj_par_rhs = sum(vv % v1(its:ite,jts:jte,1:be%v1%mz) * &
105 vv2_v1(its:ite,jts:jte,1:be%v1%mz)) + adj_par_rhs
106 if (be % v2 % mz > 0) &
107 adj_par_rhs = sum(vv % v2(its:ite,jts:jte,1:be%v2%mz) * &
108 vv2_v2(its:ite,jts:jte,1:be%v2%mz)) + adj_par_rhs
109 if (be % v3 % mz > 0) &
110 adj_par_rhs = sum(vv % v3(its:ite,jts:jte,1:be%v3%mz) * &
111 vv2_v3(its:ite,jts:jte,1:be%v3%mz)) + adj_par_rhs
112 if (be % v4 % mz > 0) &
113 adj_par_rhs = sum(vv % v4(its:ite,jts:jte,1:be%v4%mz) * &
114 vv2_v4(its:ite,jts:jte,1:be%v4%mz)) + adj_par_rhs
115 if (be % v5 % mz == 1) &
116 adj_par_rhs = sum(vv % v5(its:ite,jts:jte,1:be%v5%mz) * &
117 vv2_v5(its:ite,jts:jte,1:be%v5%mz)) + adj_par_rhs
118 if (be % ne > 0) &
119 adj_par_rhs = sum(vv % alpha(its:ite,jts:jte,1:be%ne) * &
120 vv2_alpha(its:ite,jts:jte,1:be%ne)) + adj_par_rhs
121
122 !----------------------------------------------------------------------
123 ! [7.0] Print output:
124 !----------------------------------------------------------------------
125
126 write(unit=stdout, fmt='(a,1pe22.14)') &
127 'Single domain < vp, vp > = ', adj_par_lhs, &
128 'Single domain < Vv_adj, Vv > = ', adj_par_rhs
129
130 adj_sum_lhs = wrf_dm_sum_real(adj_par_lhs)
131 adj_sum_rhs = wrf_dm_sum_real(adj_par_rhs)
132
133 write(unit=stdout,fmt='(A,2F15.2)') &
134 'TEST_COVERAGE_da_check_vvtovp_adjoint: adj_sum_lhs,adj_sum_rhs = ', &
135 adj_sum_lhs,adj_sum_rhs
136
137 if (rootproc) then
138 write(unit=stdout, fmt='(/)')
139 write(unit=stdout, fmt='(a,1pe22.14)') &
140 'Whole Domain: < Vp, Vp > = ', adj_sum_lhs, &
141 'Whole Domain: < Vv_adj, Vv > = ', adj_sum_rhs
142 end if
143
144 vv % v1(its:ite,jts:jte,1:be%v1%mz) = vv2_v1(its:ite,jts:jte,1:be%v1%mz)
145 vv % v2(its:ite,jts:jte,1:be%v2%mz) = vv2_v2(its:ite,jts:jte,1:be%v2%mz)
146 vv % v3(its:ite,jts:jte,1:be%v3%mz) = vv2_v3(its:ite,jts:jte,1:be%v3%mz)
147 vv % v4(its:ite,jts:jte,1:be%v4%mz) = vv2_v4(its:ite,jts:jte,1:be%v4%mz)
148 vv % v5(its:ite,jts:jte,1:be%v5%mz) = vv2_v5(its:ite,jts:jte,1:be%v5%mz)
149
150 if (be % ne > 0) then
151 vv % alpha(its:ite,jts:jte,1:be%ne) = vv2_alpha(its:ite,jts:jte,1:be%ne)
152 end if
153
154 write(unit=stdout, fmt='(/a/)') 'da_check_vvtovp_adjoint: Test Finished.'
155
156 if (trace_use) call da_trace_exit("da_check_vvtovp_adjoint")
157
158 end subroutine da_check_vvtovp_adjoint
159
160