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 !-------------------------------------------------------------------------
31 ! [1.0] Initialise:
32 !-------------------------------------------------------------------------
33
34 write(unit=*, fmt='(/a/)') 'da_check_vptox_adjoint: Test Results:'
35
36 call da_zero_x(xa)
37
38 vp2_v1(:,:,:) = vp % v1(:,:,:)
39 vp2_v2(:,:,:) = vp % v2(:,:,:)
40
41 call da_psichi_to_uv(vp % v1, vp % v2, xb % coefx, &
42 xb % coefy , xa % u, xa % v)
43
44 adj_par_lhs = sum(xa % u(its:ite,jts:jte,:)**2) / typical_u_rms**2
45 adj_par_lhs = sum(xa % v(its:ite,jts:jte,:)**2) / typical_v_rms**2 + &
46 adj_par_lhs
47
48 xa % u(:,:,:) = xa % u(:,:,:) / typical_u_rms**2
49 xa % v(:,:,:) = xa % v(:,:,:) / typical_v_rms**2
50
51 vp%v1(:,:,:)=0.
52 vp%v2(:,:,:)=0.
53
54 call da_psichi_to_uv_adj(xa % u, xa % v, xb % coefx, &
55 xb % coefy, vp % v1, vp % v2)
56
57 adj_par_rhs = sum(vp % v1(its:ite,jts:jte,:) * vp2_v1(its:ite,jts:jte,:))
58 adj_par_rhs = sum(vp % v2(its:ite,jts:jte,:) * vp2_v2(its:ite,jts:jte,:)) + &
59 adj_par_rhs
60
61 adj_sum_lhs = wrf_dm_sum_real(adj_par_lhs)
62 adj_sum_rhs = wrf_dm_sum_real(adj_par_rhs)
63 write(unit=stdout,fmt='(A,2F10.2)') &
64 'TEST_COVERAGE_da_check_vptox_adjoint A: adj_sum_lhs,adj_sum_rhs = ', &
65 adj_sum_lhs,adj_sum_rhs
66
67 if (rootproc) then
68 write(unit=stdout, fmt='(/a/)') &
69 ' da_check_da_psichi_to_uv: Adjoint Test Results:'
70
71 write(unit=stdout, fmt='(/)')
72 write(unit=stdout, fmt='(a,1pe22.14)') &
73 'Whole Domain: < u_v, u_v > = ', adj_sum_lhs, &
74 'Whole Domain: < psi_chi, psi_chi_adj > = ', adj_sum_rhs
75 end if
76
77 vp%v1(:,:,:) = vp2_v1(:,:,:)
78 vp%v2(:,:,:) = vp2_v2(:,:,:)
79
80 call da_zero_x(xa)
81
82 vp2_v1(:,:,:) = vp % v1(:,:,:)
83 vp2_v2(:,:,:) = vp % v2(:,:,:)
84 vp2_v3(:,:,:) = vp % v3(:,:,:)
85 vp2_v4(:,:,:) = vp % v4(:,:,:)
86 vp2_v5(:,:,:) = vp % v5(:,:,:)
87 if (be % ne > 0) vp2_alpha(:,:,:) = vp % alpha(:,:,:)
88
89 !-------------------------------------------------------------------------
90 ! [2.0] Perform x = U vp transform:
91 !-------------------------------------------------------------------------
92
93 call da_transform_vptox(xb, vp, xp, xa, be, ep)
94
95 !-------------------------------------------------------------------------
96 ! [3.0] Calculate LHS of adjoint test equation:
97 !-------------------------------------------------------------------------
98
99 ! xa % u(:,:,:) = 0.0
100 ! xa % v(:,:,:) = 0.0
101 ! xa % t(:,:,:) = 0.0
102 ! xa % q(:,:,:) = 0.0
103 ! xa%psfc(:,:) = 0.0
104
105 ! xa % p(:,:,:) = 0.0
106 ! xa % rho(:,:,:) = 0.0
107 ! xa % w(:,:,:) = 0.0
108 ! xa % wh(:,:,:) = 0.0
109 ! xa % rh(:,:,:) = 0.0
110 ! xa % qt(:,:,:) = 0.0
111 ! xa % qcw(:,:,:) = 0.0
112 ! xa % qrn(:,:,:) = 0.0
113
114 adj_par_lhs = sum(xa%u(its:ite,jts:jte,:)**2)/typical_u_rms**2
115 adj_par_lhs = sum(xa%v(its:ite,jts:jte,:)**2)/typical_v_rms**2 + adj_par_lhs
116 adj_par_lhs = sum(xa%t(its:ite,jts:jte,:)**2)/typical_t_rms**2 + adj_par_lhs
117 adj_par_lhs = sum(xa%q(its:ite,jts:jte,:)**2)/typical_q_rms**2 + adj_par_lhs
118 adj_par_lhs = sum(xa%psfc(its:ite,jts:jte)**2)/typical_p_rms**2 + adj_par_lhs
119
120 adj_par_lhs = sum(xa%p(its:ite,jts:jte,:)**2)/typical_p_rms**2 + adj_par_lhs
121 adj_par_lhs = sum(xa%rho(its:ite,jts:jte,:)**2)/typical_rho_rms**2 + &
122 adj_par_lhs
123
124 if (Use_RadarObs) then
125 adj_par_lhs = adj_par_lhs &
126 + sum(xa % wh (its:ite, jts:jte, kts:kte)**2)/typical_w_rms**2
127 else
128 adj_par_lhs = adj_par_lhs &
129 + sum(xa % w (its:ite, jts:jte, kts:kte)**2)/typical_w_rms**2
130 end if
131
132 if (cv_options_hum == 2) then
133 adj_par_lhs = sum(xa % rh(its:ite,jts:jte,:)**2) / &
134 typical_rh_rms**2 + adj_par_lhs
135 end if
136
137 if (cv_options_hum == 3) then
138 adj_par_lhs = sum(xa % qcw(its:ite,jts:jte,:)**2) / &
139 typical_qcw_rms**2 + adj_par_lhs
140 adj_par_lhs = sum(xa % qrn(its:ite,jts:jte,:)**2) / &
141 typical_qrn_rms**2 + adj_par_lhs
142 adj_par_lhs = sum(xa % qt (its:ite,jts:jte,:)**2) / &
143 typical_q_rms**2 + adj_par_lhs
144 end if
145
146 !-------------------------------------------------------------------------
147 ! [4.0] Rescale input to adjoint routine:
148 !-------------------------------------------------------------------------
149
150 xa % u(:,:,:) = xa % u(:,:,:) / typical_u_rms**2
151 xa % v(:,:,:) = xa % v(:,:,:) / typical_v_rms**2
152 xa % t(:,:,:) = xa % t(:,:,:) / typical_t_rms**2
153 xa % q(:,:,:) = xa % q(:,:,:) / typical_q_rms**2
154 xa%psfc(:,:) = xa%psfc(:,:) / typical_p_rms**2
155
156 xa % p(:,:,:) = xa % p(:,:,:) / typical_p_rms**2
157 xa % rho(:,:,:) = xa % rho(:,:,:) / typical_rho_rms**2
158
159 if (Use_RadarObs) then
160 xa %wh(:,:,:) = xa %wh(:,:,:) / typical_w_rms**2
161 xa % w(:,:,:) = 0.0
162 else
163 xa %w (:,:,:) = xa %w (:,:,:) / typical_w_rms**2
164 end if
165
166 if (cv_options_hum == 2) then
167 xa % rh(:,:,:) = xa % rh(:,:,:) / typical_rh_rms**2
168 end if
169
170 if (cv_options_hum == 3) then
171 xa % qcw(:,:,:) = xa % qcw(:,:,:) / typical_qcw_rms**2
172 xa % qrn(:,:,:) = xa % qrn(:,:,:) / typical_qrn_rms**2
173 xa % qt (:,:,:) = xa % qt (:,:,:) / typical_q_rms**2
174 end if
175
176 !-------------------------------------------------------------------------
177 ! [5.0] Perform adjoint operation:
178 !-------------------------------------------------------------------------
179
180 call da_zero_vp_type (vp)
181
182 call da_transform_vptox_adj(xb, xa, vp, be, ep, xp)
183
184 !-------------------------------------------------------------------------
185 ! [6.0] Calculate RHS of adjoint test equation:
186 !-------------------------------------------------------------------------
187
188 adj_par_rhs = sum(vp % v1(its:ite,jts:jte,:) * vp2_v1(its:ite,jts:jte,:))
189 adj_par_rhs = sum(vp % v2(its:ite,jts:jte,:) * vp2_v2(its:ite,jts:jte,:)) + &
190 adj_par_rhs
191 adj_par_rhs = sum(vp % v3(its:ite,jts:jte,:) * vp2_v3(its:ite,jts:jte,:)) + &
192 adj_par_rhs
193 adj_par_rhs = sum(vp % v4(its:ite,jts:jte,:) * vp2_v4(its:ite,jts:jte,:)) + &
194 adj_par_rhs
195 adj_par_rhs = sum(vp % v5(its:ite,jts:jte,:) * vp2_v5(its:ite,jts:jte,:)) + &
196 adj_par_rhs
197
198 if (be % ne > 0) then
199 adj_par_rhs = sum(vp % alpha(its:ite,jts:jte,:) * &
200 vp2_alpha(its:ite,jts:jte,:)) + adj_par_rhs
201 end if
202
203 !-------------------------------------------------------------------------
204 ! [7.0] Print output:
205 !-------------------------------------------------------------------------
206
207 adj_sum_lhs = wrf_dm_sum_real(adj_par_lhs)
208 adj_sum_rhs = wrf_dm_sum_real(adj_par_rhs)
209 write (unit=stdout,fmt='(A,2F10.2)') &
210 'TEST_COVERAGE_da_check_vptox_adjoint B: adj_sum_lhs,adj_sum_rhs = ', &
211 adj_sum_lhs,adj_sum_rhs
212
213 if (rootproc) then
214 write(unit=stdout, fmt='(/)')
215 write(unit=stdout, fmt='(a,1pe22.14)') &
216 'Whole Domain: < x, x > = ', adj_sum_lhs, &
217 'Whole Domain: < vp_adj, vp > = ', adj_sum_rhs
218 end if
219
220 vp % v1(:,:,:) = vp2_v1(:,:,:)
221 vp % v2(:,:,:) = vp2_v2(:,:,:)
222 vp % v3(:,:,:) = vp2_v3(:,:,:)
223 vp % v4(:,:,:) = vp2_v4(:,:,:)
224 vp % v5(:,:,:) = vp2_v5(:,:,:)
225 if (be % ne > 0) vp % alpha(:,:,:) = vp2_alpha(:,:,:)
226
227 write(unit=stdout, fmt='(/a/)') &
228 'End of da_check_vptox_adjoint.'
229
230 end subroutine da_check_vptox_adjoint
231
232