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