da_check_vptox_adjoint.inc
References to this file elsewhere.
1 subroutine da_check_vptox_adjoint(ne, xb, be, ep, xp, vp, xa)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Test Vp to X routine and adjoint for compatibility.
5 !
6 ! Method: Standard adjoint test: < x, x > = < v_adj, v >.
7 !---------------------------------------------------------------------------
8
9 implicit none
10
11 integer, intent(in) :: ne ! Ensemble size.
12 type (xb_type), intent(in) :: xb ! first guess.
13 type (be_type), intent(in) :: be ! background errors.
14 type (ep_type), intent(in) :: ep ! Ensemble perturbation type.
15 type (xpose_type), intent(inout) :: xp ! Dimensions and xpose buffers.
16 type (vp_type),intent(inout) :: vp ! grdipt/level cv (local).
17 type (x_type), intent(inout) :: xa ! grad_x(jo)
18
19 real :: adj_par_lhs ! < x, x >
20 real :: adj_par_rhs ! < v_adj, v >
21 real :: adj_sum_lhs ! < x, x >
22 real :: adj_sum_rhs ! < v_adj, v >
23 real :: vp2_v1(ims:ime,jms:jme,kms:kme)
24 real :: vp2_v2(ims:ime,jms:jme,kms:kme)
25 real :: vp2_v3(ims:ime,jms:jme,kms:kme)
26 real :: vp2_v4(ims:ime,jms:jme,kms:kme)
27 real :: vp2_v5(ims:ime,jms:jme,kms:kme)
28 real :: vp2_alpha(ims:ime,jms:jme,1:ne)
29
30 if (trace_use) call da_trace_entry("da_check_vptox_adjoint")
31
32 !-------------------------------------------------------------------------
33 ! [1.0] Initialise:
34 !-------------------------------------------------------------------------
35
36 write(unit=*, fmt='(/a/)') 'da_check_vptox_adjoint: Test Results:'
37
38 call da_zero_x(xa)
39
40 vp2_v1(:,:,:) = vp % v1(:,:,:)
41 vp2_v2(:,:,:) = vp % v2(:,:,:)
42
43 call da_psichi_to_uv(vp % v1, vp % v2, xb % coefx, &
44 xb % coefy , xa % u, xa % v)
45
46 adj_par_lhs = sum(xa % u(its:ite,jts:jte,:)**2) / typical_u_rms**2
47 adj_par_lhs = sum(xa % v(its:ite,jts:jte,:)**2) / typical_v_rms**2 + &
48 adj_par_lhs
49
50 xa % u(:,:,:) = xa % u(:,:,:) / typical_u_rms**2
51 xa % v(:,:,:) = xa % v(:,:,:) / typical_v_rms**2
52
53 vp%v1(:,:,:)=0.
54 vp%v2(:,:,:)=0.
55
56 call da_psichi_to_uv_adj(xa % u, xa % v, xb % coefx, &
57 xb % coefy, vp % v1, vp % v2)
58
59 adj_par_rhs = sum(vp % v1(its:ite,jts:jte,:) * vp2_v1(its:ite,jts:jte,:))
60 adj_par_rhs = sum(vp % v2(its:ite,jts:jte,:) * vp2_v2(its:ite,jts:jte,:)) + &
61 adj_par_rhs
62
63 adj_sum_lhs = wrf_dm_sum_real(adj_par_lhs)
64 adj_sum_rhs = wrf_dm_sum_real(adj_par_rhs)
65 write(unit=stdout,fmt='(A,2F10.2)') &
66 'TEST_COVERAGE_da_check_vptox_adjoint A: adj_sum_lhs,adj_sum_rhs = ', &
67 adj_sum_lhs,adj_sum_rhs
68
69 if (rootproc) then
70 write(unit=stdout, fmt='(/a/)') &
71 ' da_check_da_psichi_to_uv: Adjoint Test Results:'
72
73 write(unit=stdout, fmt='(/)')
74 write(unit=stdout, fmt='(a,1pe22.14)') &
75 'Whole Domain: < u_v, u_v > = ', adj_sum_lhs, &
76 'Whole Domain: < psi_chi, psi_chi_adj > = ', adj_sum_rhs
77 end if
78
79 vp%v1(:,:,:) = vp2_v1(:,:,:)
80 vp%v2(:,:,:) = vp2_v2(:,:,:)
81
82 call da_zero_x(xa)
83
84 vp2_v1(:,:,:) = vp % v1(:,:,:)
85 vp2_v2(:,:,:) = vp % v2(:,:,:)
86 vp2_v3(:,:,:) = vp % v3(:,:,:)
87 vp2_v4(:,:,:) = vp % v4(:,:,:)
88 vp2_v5(:,:,:) = vp % v5(:,:,:)
89 if (be % ne > 0) vp2_alpha(:,:,:) = vp % alpha(:,:,:)
90
91 !-------------------------------------------------------------------------
92 ! [2.0] Perform x = U vp transform:
93 !-------------------------------------------------------------------------
94
95 call da_transform_vptox(xb, vp, xp, xa, be, ep)
96
97 !-------------------------------------------------------------------------
98 ! [3.0] Calculate LHS of adjoint test equation:
99 !-------------------------------------------------------------------------
100
101 ! xa % u(:,:,:) = 0.0
102 ! xa % v(:,:,:) = 0.0
103 ! xa % t(:,:,:) = 0.0
104 ! xa % q(:,:,:) = 0.0
105 ! xa%psfc(:,:) = 0.0
106
107 ! xa % p(:,:,:) = 0.0
108 ! xa % rho(:,:,:) = 0.0
109 ! xa % w(:,:,:) = 0.0
110 ! xa % wh(:,:,:) = 0.0
111 ! xa % rh(:,:,:) = 0.0
112 ! xa % qt(:,:,:) = 0.0
113 ! xa % qcw(:,:,:) = 0.0
114 ! xa % qrn(:,:,:) = 0.0
115
116 adj_par_lhs = sum(xa%u(its:ite,jts:jte,:)**2)/typical_u_rms**2
117 adj_par_lhs = sum(xa%v(its:ite,jts:jte,:)**2)/typical_v_rms**2 + adj_par_lhs
118 adj_par_lhs = sum(xa%t(its:ite,jts:jte,:)**2)/typical_t_rms**2 + adj_par_lhs
119 adj_par_lhs = sum(xa%q(its:ite,jts:jte,:)**2)/typical_q_rms**2 + adj_par_lhs
120 adj_par_lhs = sum(xa%psfc(its:ite,jts:jte)**2)/typical_p_rms**2 + adj_par_lhs
121
122 adj_par_lhs = sum(xa%p(its:ite,jts:jte,:)**2)/typical_p_rms**2 + adj_par_lhs
123 adj_par_lhs = sum(xa%rho(its:ite,jts:jte,:)**2)/typical_rho_rms**2 + &
124 adj_par_lhs
125
126 if (Use_RadarObs) then
127 adj_par_lhs = adj_par_lhs &
128 + sum(xa % wh (its:ite, jts:jte, kts:kte)**2)/typical_w_rms**2
129 else
130 adj_par_lhs = adj_par_lhs &
131 + sum(xa % w (its:ite, jts:jte, kts:kte)**2)/typical_w_rms**2
132 end if
133
134 if (cv_options_hum == 2) then
135 adj_par_lhs = sum(xa % rh(its:ite,jts:jte,:)**2) / &
136 typical_rh_rms**2 + adj_par_lhs
137 end if
138
139 if (cv_options_hum == 3) then
140 adj_par_lhs = sum(xa % qcw(its:ite,jts:jte,:)**2) / &
141 typical_qcw_rms**2 + adj_par_lhs
142 adj_par_lhs = sum(xa % qrn(its:ite,jts:jte,:)**2) / &
143 typical_qrn_rms**2 + adj_par_lhs
144 adj_par_lhs = sum(xa % qt (its:ite,jts:jte,:)**2) / &
145 typical_q_rms**2 + adj_par_lhs
146 end if
147
148 !-------------------------------------------------------------------------
149 ! [4.0] Rescale input to adjoint routine:
150 !-------------------------------------------------------------------------
151
152 xa % u(:,:,:) = xa % u(:,:,:) / typical_u_rms**2
153 xa % v(:,:,:) = xa % v(:,:,:) / typical_v_rms**2
154 xa % t(:,:,:) = xa % t(:,:,:) / typical_t_rms**2
155 xa % q(:,:,:) = xa % q(:,:,:) / typical_q_rms**2
156 xa%psfc(:,:) = xa%psfc(:,:) / typical_p_rms**2
157
158 xa % p(:,:,:) = xa % p(:,:,:) / typical_p_rms**2
159 xa % rho(:,:,:) = xa % rho(:,:,:) / typical_rho_rms**2
160
161 if (Use_RadarObs) then
162 xa %wh(:,:,:) = xa %wh(:,:,:) / typical_w_rms**2
163 xa % w(:,:,:) = 0.0
164 else
165 xa %w (:,:,:) = xa %w (:,:,:) / typical_w_rms**2
166 end if
167
168 if (cv_options_hum == 2) then
169 xa % rh(:,:,:) = xa % rh(:,:,:) / typical_rh_rms**2
170 end if
171
172 if (cv_options_hum == 3) then
173 xa % qcw(:,:,:) = xa % qcw(:,:,:) / typical_qcw_rms**2
174 xa % qrn(:,:,:) = xa % qrn(:,:,:) / typical_qrn_rms**2
175 xa % qt (:,:,:) = xa % qt (:,:,:) / typical_q_rms**2
176 end if
177
178 !-------------------------------------------------------------------------
179 ! [5.0] Perform adjoint operation:
180 !-------------------------------------------------------------------------
181
182 call da_zero_vp_type (vp)
183
184 call da_transform_vptox_adj(xb, xa, vp, be, ep, xp)
185
186 !-------------------------------------------------------------------------
187 ! [6.0] Calculate RHS of adjoint test equation:
188 !-------------------------------------------------------------------------
189
190 adj_par_rhs = sum(vp % v1(its:ite,jts:jte,:) * vp2_v1(its:ite,jts:jte,:))
191 adj_par_rhs = sum(vp % v2(its:ite,jts:jte,:) * vp2_v2(its:ite,jts:jte,:)) + &
192 adj_par_rhs
193 adj_par_rhs = sum(vp % v3(its:ite,jts:jte,:) * vp2_v3(its:ite,jts:jte,:)) + &
194 adj_par_rhs
195 adj_par_rhs = sum(vp % v4(its:ite,jts:jte,:) * vp2_v4(its:ite,jts:jte,:)) + &
196 adj_par_rhs
197 adj_par_rhs = sum(vp % v5(its:ite,jts:jte,:) * vp2_v5(its:ite,jts:jte,:)) + &
198 adj_par_rhs
199
200 if (be % ne > 0) then
201 adj_par_rhs = sum(vp % alpha(its:ite,jts:jte,:) * &
202 vp2_alpha(its:ite,jts:jte,:)) + adj_par_rhs
203 end if
204
205 !-------------------------------------------------------------------------
206 ! [7.0] Print output:
207 !-------------------------------------------------------------------------
208
209 adj_sum_lhs = wrf_dm_sum_real(adj_par_lhs)
210 adj_sum_rhs = wrf_dm_sum_real(adj_par_rhs)
211 write (unit=stdout,fmt='(A,2F10.2)') &
212 'TEST_COVERAGE_da_check_vptox_adjoint B: adj_sum_lhs,adj_sum_rhs = ', &
213 adj_sum_lhs,adj_sum_rhs
214
215 if (rootproc) then
216 write(unit=stdout, fmt='(/)')
217 write(unit=stdout, fmt='(a,1pe22.14)') &
218 'Whole Domain: < x, x > = ', adj_sum_lhs, &
219 'Whole Domain: < vp_adj, vp > = ', adj_sum_rhs
220 end if
221
222 vp % v1(:,:,:) = vp2_v1(:,:,:)
223 vp % v2(:,:,:) = vp2_v2(:,:,:)
224 vp % v3(:,:,:) = vp2_v3(:,:,:)
225 vp % v4(:,:,:) = vp2_v4(:,:,:)
226 vp % v5(:,:,:) = vp2_v5(:,:,:)
227 if (be % ne > 0) vp % alpha(:,:,:) = vp2_alpha(:,:,:)
228
229 write(unit=stdout, fmt='(/a/)') &
230 'End of da_check_vptox_adjoint.'
231
232 if (trace_use) call da_trace_exit("da_check_vptox_adjoint")
233
234 end subroutine da_check_vptox_adjoint
235
236