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