da_check_xtoy_adjoint.inc
References to this file elsewhere.
1 subroutine da_check_xtoy_adjoint(xb, xa, iv, xp, y)
2
3 !--------------------------------------------------------------------------
4 ! Purpose: Test observation operator transform and adjoint for compatibility.
5 !
6 ! Method: Standard adjoint test: < y, y > = < x, x_adj >.
7 !---------------------------------------------------------------------------
8
9 implicit none
10
11 type (xb_type), intent(in) :: xb ! first guess (local).
12 type (x_type), intent(inout) :: xa ! analysis increments (local).
13 type (ob_type), intent(in) :: iv ! ob. increment vector.
14 type (xpose_type), intent(in) :: xp ! Dimensions and xpose buffers(be).
15 type (y_type), intent(inout) :: y ! y = h (xa)
16
17 real :: adj_ttl_lhs ! < y, y >
18 real :: adj_ttl_rhs ! < x, x_adj >
19
20 real :: partial_lhs ! < y, y >
21 real :: partial_rhs ! < x, x_adj >
22
23 real :: pertile_lhs ! < y, y >
24 real :: pertile_rhs ! < x, x_adj >
25
26 real, dimension(ims:ime, jms:jme, kms:kme) :: xa2_u, xa2_v, xa2_t, &
27 xa2_p, xa2_q, xa2_rh
28 real, dimension(ims:ime, jms:jme, kms:kme) :: xa2_w
29 real, dimension(ims:ime, jms:jme) :: xa2_psfc
30 real, dimension(ims:ime, jms:jme, kms:kme) :: xa2_qcw, xa2_qrn
31
32 if (trace_use) call da_trace_entry("da_check_xtoy_adjoint")
33
34 write (unit=stdout, fmt='(/a)') ' da_check_xtoy_adjoint: Adjoint Test Results:'
35
36 !----------------------------------------------------------------------
37 ! [1.0] Initialise:
38 !----------------------------------------------------------------------
39
40 ! Exchange XA (u, v, t, p, q, psfc) halo region.
41 call wrf_dm_halo (xp%domdesc,xp%comms,xp%halo_id4)
42
43 ! xa%u(ims:ime, jms:jme, kms:kme) = 0.0
44 ! xa%v(ims:ime, jms:jme, kms:kme) = 0.0
45 ! xa%w(ims:ime, jms:jme, kms:kme) = 0.0
46 ! xa%t(ims:ime, jms:jme, kms:kme) = 0.0
47 ! xa%p(ims:ime, jms:jme, kms:kme) = 0.0
48 ! xa%q(ims:ime, jms:jme, kms:kme) = 0.0
49 ! xa%rh(ims:ime, jms:jme, kms:kme) = 0.0
50 ! xa%psfc(ims:ime, jms:jme) = 0.0
51
52 xa2_u(ims:ime, jms:jme, kms:kme) = xa%u(ims:ime, jms:jme, kms:kme)
53 xa2_v(ims:ime, jms:jme, kms:kme) = xa%v(ims:ime, jms:jme, kms:kme)
54 xa2_t(ims:ime, jms:jme, kms:kme) = xa%t(ims:ime, jms:jme, kms:kme)
55 xa2_p(ims:ime, jms:jme, kms:kme) = xa%p(ims:ime, jms:jme, kms:kme)
56 xa2_q(ims:ime, jms:jme, kms:kme) = xa%q(ims:ime, jms:jme, kms:kme)
57 xa2_w(ims:ime, jms:jme, kms:kme) = xa%w(ims:ime, jms:jme, kms:kme)
58 xa2_rh(ims:ime, jms:jme, kms:kme)= xa%rh(ims:ime, jms:jme, kms:kme)
59 xa2_psfc(ims:ime, jms:jme) = xa%psfc(ims:ime, jms:jme)
60
61 xa2_qcw(ims:ime, jms:jme, kms:kme) = xa%qcw(ims:ime, jms:jme, kms:kme)
62 xa2_qrn(ims:ime, jms:jme, kms:kme) = xa%qrn(ims:ime, jms:jme, kms:kme)
63
64 call da_pt_to_rho_lin(xb, xa, xp)
65
66 if (sfc_assi_options == 2) then
67 call da_transform_xtowtq (xp, xb, xa)
68 ! Exchange XA (surface variable: u10, v10, t2, q2) halo region.
69 call wrf_dm_halo(xp%domdesc,xp%comms,xp%halo_id6)
70 end if
71
72 if (use_ssmt1obs .or. use_ssmt2obs .or. use_GpspwObs .or. &
73 use_GpsrefObs .or. &
74 use_SsmiTbObs .or. use_SsmiRetrievalObs) then
75
76 ! Now do something for PW
77 call da_transform_xtotpw(xa, xb)
78
79 ! GPS Refractivity:
80 if (use_gpsrefobs) then
81 call da_transform_xtogpsref_lin(xa, xb, xp)
82 end if
83
84 if (use_ssmt1obs .or. use_ssmt2obs .or. &
85 use_SsmiTbObs .or. Use_SsmiRetrievalObs) then
86 if (global) then
87 call da_error(__FILE__,__LINE__, &
88 (/"xb%speed is not available, see da_transfer_kmatoxb.inc"/))
89 end if
90 call da_transform_xtoseasfcwind_lin(xa, xb)
91 end if
92
93 if (use_SsmiTbObs) call da_transform_xtotb_lin (xa, xb)
94
95 ! Exchange XA halo region.
96 call wrf_dm_halo(xp%domdesc,xp%comms,xp%halo_id8)
97 end if
98
99 !----------------------------------------------------------------------
100 ! [2.0] Perform y = Hx transform:
101 !----------------------------------------------------------------------
102
103 call da_transform_xtoy(xb, iv, xa, xp, y)
104
105 !----------------------------------------------------------------------
106 ! [3.0] Calculate LHS of adjoint test equation and
107 ! Rescale input to adjoint routine :
108 !----------------------------------------------------------------------
109
110 partial_lhs = 0.0
111 pertile_lhs = 0.0
112
113 if (use_SoundObs) then
114 call da_check_xtoy_adjoint_sound(iv, y, partial_lhs, pertile_lhs)
115 call da_check_xtoy_adjoint_sonde_sfc(iv, y, partial_lhs, pertile_lhs)
116 end if
117
118 if (use_SynopObs) &
119 call da_check_xtoy_adjoint_synop(iv, y, partial_lhs, pertile_lhs)
120
121 if (use_GeoAMVObs) &
122 call da_check_xtoy_adjoint_geoamv(iv, y, partial_lhs, pertile_lhs)
123
124 if (use_PolarAMVObs) &
125 call da_check_xtoy_adjoint_polaramv(iv, y, partial_lhs, pertile_lhs)
126
127 if (use_AirepObs) &
128 call da_check_xtoy_adjoint_airep(iv, y, partial_lhs, pertile_lhs)
129
130 if (use_PilotObs) &
131 call da_check_xtoy_adjoint_pilot(iv, y, partial_lhs, pertile_lhs)
132
133 if (use_RadarObs) &
134 call da_check_xtoy_adjoint_radar(iv, y, partial_lhs, pertile_lhs)
135
136 if (use_SatemObs) &
137 call da_check_xtoy_adjoint_satem(iv, y, partial_lhs, pertile_lhs)
138
139 if (use_MetarObs) &
140 call da_check_xtoy_adjoint_metar(iv, y, partial_lhs, pertile_lhs)
141
142 if (use_ShipsObs) &
143 call da_check_xtoy_adjoint_ships(iv, y, partial_lhs, pertile_lhs)
144
145 if (use_GpspwObs) &
146 call da_check_xtoy_adjoint_gpspw(iv, y, partial_lhs, pertile_lhs)
147
148 if (use_GpsrefObs) &
149 call da_check_xtoy_adjoint_gpsref(iv, y, partial_lhs, pertile_lhs)
150
151 if (use_SsmiTbObs .or. Use_SsmiRetrievalObs) &
152 call da_check_xtoy_adjoint_ssmi (iv, y, partial_lhs, pertile_lhs)
153
154 if (use_ssmt1obs) &
155 call da_check_xtoy_adjoint_ssmt1(iv, y, partial_lhs, pertile_lhs)
156
157 if (use_ssmt2obs) &
158 call da_check_xtoy_adjoint_ssmt2(iv, y, partial_lhs, pertile_lhs)
159
160 if (use_qscatobs) &
161 call da_check_xtoy_adjoint_qscat(iv, y, partial_lhs, pertile_lhs)
162
163 if (use_ProfilerObs) &
164 call da_check_xtoy_adjoint_profiler(iv, y, partial_lhs, pertile_lhs)
165
166 if (use_BuoyObs) &
167 call da_check_xtoy_adjoint_buoy(iv, y, partial_lhs, pertile_lhs)
168
169 if (use_BogusObs) &
170 call da_check_xtoy_adjoint_bogus(iv, y, partial_lhs, pertile_lhs)
171
172 if (use_rad) &
173 call da_check_xtoy_adjoint_rad(iv, y, partial_lhs, pertile_lhs)
174
175 !----------------------------------------------------------------------
176 ! [5.0] Perform adjoint operation:
177 !----------------------------------------------------------------------
178
179 call da_zero_x (xa)
180
181 call da_transform_xtoy_adj (xb, iv, xp, y, xa)
182
183
184 if (use_ssmt1obs .or. use_ssmt2obs .or. use_GpspwObs .or. &
185 use_GpsrefObs .or. &
186 use_SsmiTbObs .or. use_SsmiRetrievalObs) then
187
188 if (Use_SsmiTbObs) call da_transform_xtotb_adj (xa, xb)
189
190 ! for PW
191 call da_transform_xtotpw_adj (xa, xb)
192
193 ! GPS Refractivity:
194 if (use_GpsrefObs) call da_transform_xtogpsref_adj (xa, xb, xp)
195
196 if (use_ssmt1obs .or. use_ssmt2obs .or. &
197 use_SsmiTbObs .or. use_SsmiRetrievalObs) then
198 if (global) then
199 call da_error(__FILE__,__LINE__, &
200 (/"xb%speed is not available, see da_transfer_kmatoxb.inc"/))
201 end if
202 call da_transform_xtoseasfcwind_adj (xa, xb)
203 end if
204 end if
205
206 ! Now do something for surface variables
207 if (sfc_assi_options == 2) then
208 call da_transform_xtowtq_adj (xp, xb, xa)
209 end if
210
211 call da_pt_to_rho_adj (xb, xa)
212
213 pertile_rhs = sum (xa%u(ims:ime, jms:jme, kms:kme) * xa2_u(ims:ime, jms:jme, kms:kme)) &
214 + sum (xa%v(ims:ime, jms:jme, kms:kme) * xa2_v(ims:ime, jms:jme, kms:kme)) &
215 + sum (xa%w(ims:ime, jms:jme, kms:kme) * xa2_w(ims:ime, jms:jme, kms:kme)) &
216 + sum (xa%t(ims:ime, jms:jme, kms:kme) * xa2_t(ims:ime, jms:jme, kms:kme)) &
217 + sum (xa%p(ims:ime, jms:jme, kms:kme) * xa2_p(ims:ime, jms:jme, kms:kme)) &
218 + sum (xa%q(ims:ime, jms:jme, kms:kme) * xa2_q(ims:ime, jms:jme, kms:kme)) &
219 + sum (xa%rh(ims:ime, jms:jme, kms:kme)* xa2_rh(ims:ime, jms:jme, kms:kme)) &
220 + sum (xa%psfc(ims:ime, jms:jme) * xa2_psfc(ims:ime, jms:jme))
221 pertile_rhs = pertile_rhs &
222 + sum (xa%qcw(ims:ime, jms:jme, kms:kme) * xa2_qcw(ims:ime, jms:jme, kms:kme)) &
223 + sum (xa%qrn(ims:ime, jms:jme, kms:kme) * xa2_qrn(ims:ime, jms:jme, kms:kme))
224
225 !----------------------------------------------------------------------
226 ! [6.0] Calculate RHS of adjoint test equation:
227 !----------------------------------------------------------------------
228
229 partial_rhs = sum (xa%u(its:ite, jts:jte, kts:kte) * xa2_u(its:ite, jts:jte, kts:kte)) &
230 + sum (xa%v(its:ite, jts:jte, kts:kte) * xa2_v(its:ite, jts:jte, kts:kte)) &
231 + sum (xa%w(its:ite, jts:jte, kts:kte+1) * xa2_w(its:ite, jts:jte, kts:kte+1)) &
232 + sum (xa%t(its:ite, jts:jte, kts:kte) * xa2_t(its:ite, jts:jte, kts:kte)) &
233 + sum (xa%p(its:ite, jts:jte, kts:kte) * xa2_p(its:ite, jts:jte, kts:kte)) &
234 + sum (xa%q(its:ite, jts:jte, kts:kte) * xa2_q(its:ite, jts:jte, kts:kte)) &
235 + sum (xa%rh(its:ite, jts:jte, kts:kte)* xa2_rh(its:ite, jts:jte, kts:kte)) &
236 + sum (xa%psfc(its:ite, jts:jte) * xa2_psfc(its:ite, jts:jte))
237
238 partial_rhs = partial_rhs &
239 + sum (xa%qcw(its:ite, jts:jte, kts:kte) * xa2_qcw(its:ite, jts:jte, kts:kte)) &
240 + sum (xa%qrn(its:ite, jts:jte, kts:kte) * xa2_qrn(its:ite, jts:jte, kts:kte))
241
242 !----------------------------------------------------------------------
243 ! [7.0] Print output:
244 !----------------------------------------------------------------------
245
246 write (unit=stdout, fmt='(A,1pe22.14)') &
247 ' Tile < y, y > = ', pertile_lhs, &
248 ' Tile < x, x_adj > = ', pertile_rhs
249
250 adj_ttl_lhs = wrf_dm_sum_real (partial_lhs)
251 adj_ttl_rhs = wrf_dm_sum_real (partial_rhs)
252
253 if (rootproc) then
254 write (unit=stdout, fmt='(A,1pe22.14)') ' Whole Domain < y, y > = ', adj_ttl_lhs
255 write (unit=stdout, fmt='(A,1pe22.14)') ' Whole Domain < x, x_adj > = ', adj_ttl_rhs
256 end if
257
258 if (trace_use) call da_trace_exit("da_check_xtoy_adjoint")
259
260 end subroutine da_check_xtoy_adjoint
261
262