da_check_xtoy_adjoint.inc
References to this file elsewhere.
1 subroutine da_check_xtoy_adjoint(grid, config_flags, iv, 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 (domain), intent(inout) :: grid
12 type(grid_config_rec_type), intent(inout) :: config_flags
13 type (iv_type), intent(in) :: iv ! ob. increment vector.
14 type (y_type), intent(inout) :: y ! y = h (grid%xa)
15
16 real :: adj_ttl_lhs ! < y, y >
17 real :: adj_ttl_rhs ! < x, x_adj >
18
19 real :: partial_lhs ! < y, y >
20 real :: partial_rhs ! < x, x_adj >
21
22 real :: pertile_lhs ! < y, y >
23 real :: pertile_rhs ! < x, x_adj >
24
25 real, dimension(ims:ime, jms:jme, kms:kme) :: xa2_u, xa2_v, xa2_t, &
26 xa2_p, xa2_q, xa2_rh
27 real, dimension(ims:ime, jms:jme, kms:kme) :: xa2_w
28 real, dimension(ims:ime, jms:jme) :: xa2_psfc
29 real, dimension(ims:ime, jms:jme, kms:kme) :: xa2_qcw, xa2_qrn
30 integer :: nobwin, ndynopt
31 character(len=4) :: filnam
32 integer :: wrf_done_unit
33
34 if (trace_use) call da_trace_entry("da_check_xtoy_adjoint")
35
36 write (unit=stdout, fmt='(/a)') ' da_check_xtoy_adjoint: Adjoint Test Results:'
37
38 !----------------------------------------------------------------------
39 ! [1.0] Initialise:
40 !----------------------------------------------------------------------
41
42 partial_lhs = 0.0
43 pertile_lhs = 0.0
44
45 #ifdef DM_PARALLEL
46 #include "HALO_XA.inc"
47 #endif
48
49 ! grid%xa%u(ims:ime, jms:jme, kms:kme) = 0.0
50 ! grid%xa%v(ims:ime, jms:jme, kms:kme) = 0.0
51 ! grid%xa%w(ims:ime, jms:jme, kms:kme) = 0.0
52 ! grid%xa%t(ims:ime, jms:jme, kms:kme) = 0.0
53 ! grid%xa%p(ims:ime, jms:jme, kms:kme) = 0.0
54 ! grid%xa%q(ims:ime, jms:jme, kms:kme) = 0.0
55 ! grid%xa%rh(ims:ime, jms:jme, kms:kme) = 0.0
56 ! grid%xa%psfc(ims:ime, jms:jme) = 0.0
57
58 xa2_u(ims:ime, jms:jme, kms:kme) = grid%xa%u(ims:ime, jms:jme, kms:kme)
59 xa2_v(ims:ime, jms:jme, kms:kme) = grid%xa%v(ims:ime, jms:jme, kms:kme)
60 xa2_t(ims:ime, jms:jme, kms:kme) = grid%xa%t(ims:ime, jms:jme, kms:kme)
61 xa2_p(ims:ime, jms:jme, kms:kme) = grid%xa%p(ims:ime, jms:jme, kms:kme)
62 xa2_q(ims:ime, jms:jme, kms:kme) = grid%xa%q(ims:ime, jms:jme, kms:kme)
63 xa2_w(ims:ime, jms:jme, kms:kme) = grid%xa%w(ims:ime, jms:jme, kms:kme)
64 xa2_rh(ims:ime, jms:jme, kms:kme)= grid%xa%rh(ims:ime, jms:jme, kms:kme)
65 xa2_psfc(ims:ime, jms:jme) = grid%xa%psfc(ims:ime, jms:jme)
66
67 xa2_qcw(ims:ime, jms:jme, kms:kme) = grid%xa%qcw(ims:ime, jms:jme, kms:kme)
68 xa2_qrn(ims:ime, jms:jme, kms:kme) = grid%xa%qrn(ims:ime, jms:jme, kms:kme)
69
70 if (var4d) then
71 call da_transfer_xatowrftl(grid, config_flags, 'tl01')
72
73 #ifdef DM_PARALLEL
74 call da_system("da_run_wrfplus_tl.ksh pre")
75 if (rootproc) then
76 call da_system("rm -rf wrf_done")
77 call da_system("touch wrf_go_ahead")
78 call da_get_unit(wrf_done_unit)
79 do while (.true.)
80 open(wrf_done_unit,file="wrf_done",status="old",err=303)
81 close(wrf_done_unit)
82 exit
83 303 continue
84 call da_system("sleep 1")
85 end do
86 call da_free_unit(wrf_done_unit)
87 end if
88 call wrf_get_dm_communicator ( comm )
89 call mpi_barrier(comm, ierr)
90 call da_system("da_run_wrfplus_tl.ksh post")
91 #else
92 call da_system("da_run_wrfplus_tl.ksh")
93 #endif
94 end if
95
96 do nobwin=1, num_fgat_time
97
98 if (var4d) then
99 write(filnam,'(a2,i2.2)') 'tl',nobwin
100 call da_transfer_wrftltoxa(grid, config_flags, filnam)
101 end if
102
103 call da_pt_to_rho_lin(grid)
104
105 if (sfc_assi_options == 2) then
106 call da_transform_xtowtq (grid)
107 #ifdef DM_PARALLEL
108 #include "HALO_SFC_XA.inc"
109 #endif
110 end if
111
112 if (use_ssmt1obs .or. use_ssmt2obs .or. use_gpspwobs .or. &
113 use_gpsrefobs .or. &
114 use_ssmitbobs .or. use_ssmiretrievalobs) then
115
116 ! Now do something for PW
117 call da_transform_xtotpw(grid)
118
119 ! GPS Refractivity:
120 if (use_gpsrefobs) then
121 call da_transform_xtogpsref_lin(grid)
122 end if
123
124 if (use_ssmt1obs .or. use_ssmt2obs .or. &
125 use_ssmitbobs .or. use_ssmiretrievalobs) then
126 if (global) then
127 call da_error(__FILE__,__LINE__, &
128 (/"grid%xb%speed is not available, see da_transfer_kmatoxb.inc"/))
129 end if
130 call da_transform_xtoseasfcwind_lin(grid)
131 end if
132
133 if (use_ssmitbobs) call da_transform_xtotb_lin (grid)
134
135 #ifdef DM_PARALLEL
136 #include "HALO_SSMI_XA.inc"
137 #endif
138 end if
139
140 !----------------------------------------------------------------------
141 ! [2.0] Perform y = Hx transform:
142 !----------------------------------------------------------------------
143
144 call da_transform_xtoy (grid, iv, y)
145
146 !----------------------------------------------------------------------
147 ! [3.0] Calculate LHS of adjoint test equation and
148 ! Rescale input to adjoint routine :
149 !----------------------------------------------------------------------
150
151 if (iv%info(sound)%nlocal > 0) then
152 call da_check_xtoy_adjoint_sound (iv, y, partial_lhs, pertile_lhs)
153 call da_check_xtoy_adjoint_sonde_sfc (iv, y, partial_lhs, pertile_lhs)
154 end if
155
156 if (iv%info(synop)%nlocal > 0) call da_check_xtoy_adjoint_synop (iv, y, partial_lhs, pertile_lhs)
157 if (iv%info(geoamv)%nlocal > 0) call da_check_xtoy_adjoint_geoamv (iv, y, partial_lhs, pertile_lhs)
158 if (iv%info(polaramv)%nlocal > 0) call da_check_xtoy_adjoint_polaramv (iv, y, partial_lhs, pertile_lhs)
159 if (iv%info(airep)%nlocal > 0) call da_check_xtoy_adjoint_airep (iv, y, partial_lhs, pertile_lhs)
160 if (iv%info(pilot)%nlocal > 0) call da_check_xtoy_adjoint_pilot (iv, y, partial_lhs, pertile_lhs)
161 if (iv%info(radar)%nlocal > 0) call da_check_xtoy_adjoint_radar (iv, y, partial_lhs, pertile_lhs)
162 if (iv%info(satem)%nlocal > 0) call da_check_xtoy_adjoint_satem (iv, y, partial_lhs, pertile_lhs)
163 if (iv%info(metar)%nlocal > 0) call da_check_xtoy_adjoint_metar (iv, y, partial_lhs, pertile_lhs)
164 if (iv%info(ships)%nlocal > 0) call da_check_xtoy_adjoint_ships (iv, y, partial_lhs, pertile_lhs)
165 if (iv%info(gpspw)%nlocal > 0) call da_check_xtoy_adjoint_gpspw (iv, y, partial_lhs, pertile_lhs)
166 if (iv%info(gpsref)%nlocal > 0) call da_check_xtoy_adjoint_gpsref (iv, y, partial_lhs, pertile_lhs)
167 if (iv%info(ssmi_tb)%nlocal > 0) call da_check_xtoy_adjoint_ssmi_tb (iv, y, partial_lhs, pertile_lhs)
168 if (iv%info(ssmi_rv)%nlocal > 0) call da_check_xtoy_adjoint_ssmi_rv (iv, y, partial_lhs, pertile_lhs)
169 if (iv%info(ssmt2)%nlocal > 0) call da_check_xtoy_adjoint_ssmt1 (iv, y, partial_lhs, pertile_lhs)
170 if (iv%info(ssmt2)%nlocal > 0) call da_check_xtoy_adjoint_ssmt2 (iv, y, partial_lhs, pertile_lhs)
171 if (iv%info(qscat)%nlocal > 0) call da_check_xtoy_adjoint_qscat (iv, y, partial_lhs, pertile_lhs)
172 if (iv%info(profiler)%nlocal > 0) call da_check_xtoy_adjoint_profiler (iv, y, partial_lhs, pertile_lhs)
173 if (iv%info(buoy)%nlocal > 0) call da_check_xtoy_adjoint_buoy (iv, y, partial_lhs, pertile_lhs)
174 if (iv%info(bogus)%nlocal > 0) call da_check_xtoy_adjoint_bogus (iv, y, partial_lhs, pertile_lhs)
175 if (iv%num_inst > 0) call da_check_xtoy_adjoint_rad (iv, y, partial_lhs, pertile_lhs)
176
177 !----------------------------------------------------------------------
178 ! [5.0] Perform adjoint operation:
179 !----------------------------------------------------------------------
180
181 call da_zero_x (grid%xa)
182
183 call da_transform_xtoy_adj (grid, iv, y, grid%xa)
184
185
186 if (use_ssmt1obs .or. use_ssmt2obs .or. use_gpspwobs .or. &
187 use_gpsrefobs .or. &
188 use_ssmitbobs .or. use_ssmiretrievalobs) then
189
190 if (use_ssmitbobs) call da_transform_xtotb_adj (grid)
191
192 ! for PW
193 call da_transform_xtotpw_adj (grid)
194
195 ! GPS Refractivity:
196 if (use_gpsrefobs) call da_transform_xtogpsref_adj (grid)
197
198 if (use_ssmt1obs .or. use_ssmt2obs .or. &
199 use_ssmitbobs .or. use_ssmiretrievalobs) then
200 if (global) then
201 call da_error(__FILE__,__LINE__, &
202 (/"grid%xb%speed is not available, see da_transfer_kmatoxb.inc"/))
203 end if
204 call da_transform_xtoseasfcwind_adj (grid)
205 end if
206 end if
207
208 ! Now do something for surface variables
209 if (sfc_assi_options == 2) then
210 call da_transform_xtowtq_adj (grid)
211 end if
212
213 call da_pt_to_rho_adj (grid)
214
215 if (var4d) then
216
217 grid%g_u_2 = 0.0
218 grid%g_v_2 = 0.0
219 grid%g_w_2 = 0.0
220 grid%g_t_2 = 0.0
221 grid%g_ph_2 = 0.0
222 grid%g_mu_2 = 0.0
223 grid%g_moist = 0.0
224
225 write(unit=filnam,fmt='(a2,i2.2)') 'af',nobwin
226
227 call da_transfer_wrftltoxa_adj(grid, config_flags, filnam)
228
229 !!!!!! exclude TL&AD integration!!!!!!!
230 #ifdef MYOK
231 ndynopt = grid%dyn_opt
232 grid%dyn_opt = DYN_EM_TL
233 call nl_set_dyn_opt (1 , DYN_EM_TL)
234 call da_med_initialdata_input(grid , config_flags, filnam)
235 grid%a_u_2 = grid%g_u_2
236 grid%a_v_2 = grid%g_v_2
237 grid%a_w_2 = grid%g_w_2
238 grid%a_t_2 = grid%g_t_2
239 grid%a_ph_2 = grid%g_ph_2
240 grid%a_mu_2 = grid%g_mu_2
241 grid%a_moist = grid%g_moist
242
243 grid%dyn_opt = DYN_EM_AD
244 call nl_set_dyn_opt (1 , DYN_EM_AD)
245 call da_med_initialdata_output(grid , config_flags, 'gr01')
246
247 grid%dyn_opt = ndynopt
248 call nl_set_dyn_opt (1 , DYN_EM)
249 grid%g_u_2 = 0.0
250 grid%g_v_2 = 0.0
251 grid%g_w_2 = 0.0
252 grid%g_t_2 = 0.0
253 grid%g_ph_2 = 0.0
254 grid%g_mu_2 = 0.0
255 grid%g_moist = 0.0
256 grid%a_u_2 = 0.0
257 grid%a_v_2 = 0.0
258 grid%a_w_2 = 0.0
259 grid%a_t_2 = 0.0
260 grid%a_ph_2 = 0.0
261 grid%a_mu_2 = 0.0
262 grid%a_moist = 0.0
263 #endif
264 !!!!!!
265 end if
266 end do
267
268 if (var4d) then
269
270 #ifndef TLDF
271 ndynopt = grid%dyn_opt
272 grid%dyn_opt = DYN_EM_TL
273 call nl_set_dyn_opt (1 , DYN_EM_TL)
274
275 call da_med_initialdata_input(grid , config_flags, 'tldf')
276 grid%g_u_2 = 0.0
277 grid%g_v_2 = 0.0
278 grid%g_w_2 = 0.0
279 grid%g_t_2 = 0.0
280 grid%g_ph_2 = 0.0
281 grid%g_mu_2 = 0.0
282 grid%g_moist = 0.0
283
284 call med_hist_out(grid , 3 , config_flags)
285
286 grid%dyn_opt = ndynopt
287 call nl_set_dyn_opt (1 , DYN_EM)
288 #endif
289
290 #ifdef DM_PARALLEL
291 call da_system("da_run_wrfplus_ad.ksh pre")
292 if (rootproc) then
293 call da_system("rm -rf wrf_done")
294 call da_system("touch wrf_go_ahead")
295 call da_get_unit(wrf_done_unit)
296 do while (.true.)
297 open(wrf_done_unit,file="wrf_done",status="old",err=304)
298 close(wrf_done_unit)
299 exit
300 304 continue
301 call da_system("sleep 1")
302 end do
303 call da_free_unit(wrf_done_unit)
304 end if
305 call wrf_get_dm_communicator ( comm )
306 call mpi_barrier(comm, ierr)
307 call da_system("da_run_wrfplus_ad.ksh post")
308 #else
309 call da_system("da_run_wrfplus_ad.ksh")
310 #endif
311
312 grid%xa%u(its:ite, jts:jte, kts:kte) = 0.0
313 grid%xa%v(its:ite, jts:jte, kts:kte) = 0.0
314 grid%xa%w(its:ite, jts:jte, kts:kte+1) = 0.0
315 grid%xa%q(its:ite, jts:jte, kts:kte) = 0.0
316 grid%xa%t(its:ite, jts:jte, kts:kte) = 0.0
317 grid%xa%psfc(its:ite, jts:jte) = 0.0
318
319 call da_transfer_xatowrftl_adj(grid, config_flags, 'gr01')
320 grid%a_u_2 = grid%xa%u
321 grid%a_v_2 = grid%xa%v
322 grid%a_w_2 = grid%xa%w
323 grid%a_t_2 = grid%xa%t
324 grid%a_ph_2 = grid%xa%p
325 grid%a_mu_2 = grid%xa%psfc
326 grid%a_moist(:,:,:,p_a_qv) = grid%xa%q(:,:,:)
327
328 grid%dyn_opt = DYN_EM_AD
329 call nl_set_dyn_opt (1 , DYN_EM_AD)
330 call da_med_initialdata_output(grid , config_flags, 'gr01')
331
332 grid%dyn_opt = ndynopt
333 call nl_set_dyn_opt (1 , DYN_EM)
334
335 end if
336
337 pertile_rhs = sum (grid%xa%u(ims:ime, jms:jme, kms:kme) * xa2_u(ims:ime, jms:jme, kms:kme)) &
338 + sum (grid%xa%v(ims:ime, jms:jme, kms:kme) * xa2_v(ims:ime, jms:jme, kms:kme)) &
339 + sum (grid%xa%w(ims:ime, jms:jme, kms:kme) * xa2_w(ims:ime, jms:jme, kms:kme)) &
340 + sum (grid%xa%t(ims:ime, jms:jme, kms:kme) * xa2_t(ims:ime, jms:jme, kms:kme)) &
341 + sum (grid%xa%p(ims:ime, jms:jme, kms:kme) * xa2_p(ims:ime, jms:jme, kms:kme)) &
342 + sum (grid%xa%q(ims:ime, jms:jme, kms:kme) * xa2_q(ims:ime, jms:jme, kms:kme)) &
343 + sum (grid%xa%rh(ims:ime, jms:jme, kms:kme)* xa2_rh(ims:ime, jms:jme, kms:kme)) &
344 + sum (grid%xa%psfc(ims:ime, jms:jme) * xa2_psfc(ims:ime, jms:jme))
345 pertile_rhs = pertile_rhs &
346 + sum (grid%xa%qcw(ims:ime, jms:jme, kms:kme) * xa2_qcw(ims:ime, jms:jme, kms:kme)) &
347 + sum (grid%xa%qrn(ims:ime, jms:jme, kms:kme) * xa2_qrn(ims:ime, jms:jme, kms:kme))
348
349 !----------------------------------------------------------------------
350 ! [6.0] Calculate RHS of adjoint test equation:
351 !----------------------------------------------------------------------
352
353 partial_rhs = sum (grid%xa%u(its:ite, jts:jte, kts:kte) * xa2_u(its:ite, jts:jte, kts:kte)) &
354 + sum (grid%xa%v(its:ite, jts:jte, kts:kte) * xa2_v(its:ite, jts:jte, kts:kte)) &
355 + sum (grid%xa%w(its:ite, jts:jte, kts:kte+1) * xa2_w(its:ite, jts:jte, kts:kte+1)) &
356 + sum (grid%xa%t(its:ite, jts:jte, kts:kte) * xa2_t(its:ite, jts:jte, kts:kte)) &
357 + sum (grid%xa%p(its:ite, jts:jte, kts:kte) * xa2_p(its:ite, jts:jte, kts:kte)) &
358 + sum (grid%xa%q(its:ite, jts:jte, kts:kte) * xa2_q(its:ite, jts:jte, kts:kte)) &
359 + sum (grid%xa%rh(its:ite, jts:jte, kts:kte)* xa2_rh(its:ite, jts:jte, kts:kte)) &
360 + sum (grid%xa%psfc(its:ite, jts:jte) * xa2_psfc(its:ite, jts:jte))
361
362 partial_rhs = partial_rhs &
363 + sum (grid%xa%qcw(its:ite, jts:jte, kts:kte) * xa2_qcw(its:ite, jts:jte, kts:kte)) &
364 + sum (grid%xa%qrn(its:ite, jts:jte, kts:kte) * xa2_qrn(its:ite, jts:jte, kts:kte))
365
366 !----------------------------------------------------------------------
367 ! [7.0] Print output:
368 !----------------------------------------------------------------------
369
370 write (unit=stdout, fmt='(A,1pe22.14)') &
371 ' Tile < y, y > = ', pertile_lhs, &
372 ' Tile < x, x_adj > = ', pertile_rhs
373
374 adj_ttl_lhs = wrf_dm_sum_real (partial_lhs)
375 adj_ttl_rhs = wrf_dm_sum_real (partial_rhs)
376
377 if (rootproc) then
378 write (unit=stdout, fmt='(A,1pe22.14)') ' Whole Domain < y, y > = ', adj_ttl_lhs
379 write (unit=stdout, fmt='(A,1pe22.14)') ' Whole Domain < x, x_adj > = ', adj_ttl_rhs
380 end if
381
382 if (trace_use) call da_trace_exit("da_check_xtoy_adjoint")
383
384 end subroutine da_check_xtoy_adjoint
385
386