da_transform_xtoy_rttov_adj.inc
References to this file elsewhere.
1 subroutine da_transform_xtoy_rttov_adj ( iv, jo_grad_y, jo_grad_x )
2
3 !---------------------------------------------------------------------------
4 ! Purpose: transform gradient from obs space to model grid space.
5 !
6 ! METHOD: jo_grad_x = H^T jo_grad_y = - H^T R^-1 ( d - H delta_x )
7 ! 1. input gradient in obs space and reference state of RTTOV
8 ! 2. call adjoint of RTM
9 ! 3. adjoint of interpolation from model grid to obs loc
10 !---------------------------------------------------------------------------
11
12 implicit none
13
14 !#include "rttov_v2q.interface"
15
16 type (x_type), intent(inout) :: jo_grad_x !
17 type (y_type), intent(in) :: jo_grad_y ! H' delta_x
18 type (iv_type), intent(in) :: iv ! O-B structure.
19
20 #ifdef RTTOV
21
22 integer :: i, j, k,c ! Index dimension.
23 integer :: nlevels ! Number of obs levels.
24 integer :: num_rad ! Number of radiance obs
25 #ifdef DM_PARALLEL
26 integer :: status(mpi_status_size) ! MPI status.
27 #endif
28
29 real, allocatable :: model_mr(:)
30
31
32 integer :: nchan, n
33
34 type(con_vars_type), allocatable :: con_vars(:), con_vars_ad(:)
35 type(aux_vars_type), allocatable :: aux_vars(:)
36
37 ! evenly distributed versions
38 type(con_vars_type), allocatable :: d_con_vars(:), d_con_vars_ad(:)
39 type(aux_vars_type), allocatable :: d_aux_vars(:)
40 real, allocatable :: d_tb(:,:)
41 integer :: d_num_rad, l_num_rad,e,s,len,nl
42 real, allocatable :: temp_t(:,:),temp_q(:,:), temp_ps(:), rtemp(:,:)
43 real, allocatable :: temp_tb(:,:)
44 integer, allocatable :: itemp(:,:)
45
46 if (iv%num_inst < 1) return
47
48 if (trace_use) call da_trace_entry("da_transform_xtoy_rttov_adj")
49
50 do i = 1, iv%num_inst ! loop for sensor
51
52 if (iv%instid(i)%num_rad_glo < 1) cycle
53
54 num_rad = iv%instid(i)%num_rad
55 nlevels = iv%instid(i)%nlevels
56 nchan = iv%instid(i)%nchan
57
58 if (num_rad > 0) then
59
60 allocate (model_mr(1:nlevels))
61
62 allocate (con_vars(num_rad))
63 allocate (con_vars_ad(num_rad))
64 allocate (aux_vars(num_rad))
65
66 !---------------------------------------------------------------
67 ! [1.0] assign tb = R^-1 Re :
68 !---------------------------------------------------------------
69
70 !---------------------------------------------
71 ! [2.0] get input of da_rttov_ad
72 !---------------------------------------------
73
74 do n = 1, num_rad
75 con_vars(n) % nlevels = nlevels
76
77 allocate (con_vars(n) % t(nlevels))
78 allocate (con_vars(n) % q(nlevels))
79 allocate (con_vars_ad(n) % t(nlevels))
80 allocate (con_vars_ad(n) % q(nlevels))
81
82 con_vars(n) % t(:) = 0.0
83 con_vars(n) % q(:) = 0.0
84 con_vars_ad(n) % t(:) = 0.0
85 con_vars_ad(n) % q(:) = 0.0
86
87 con_vars(n) % t(:) = iv%instid(i)%t(:,n)
88 con_vars(n) % q(:) = iv%instid(i)%mr(:,n)
89 con_vars(n) % ps = iv%instid(i)%ps(n)
90
91 aux_vars(n) % t2m = iv%instid(i)%t2m(n)
92 aux_vars(n) % q2m = iv%instid(i)%mr2m(n)
93 aux_vars(n) % u10 = iv%instid(i)%u10(n)
94 aux_vars(n) % v10 = iv%instid(i)%v10(n)
95
96 aux_vars(n) % surftype = iv%instid(i)%landsea_mask(n)
97 aux_vars(n) % surft = iv%instid(i)%ts(n)
98 ! aux_vars(n) % fastem(:) = 0.0
99
100 aux_vars(n) % satzen = iv%instid(i)%satzen(n)
101 aux_vars(n) % satazi = iv%instid(i)%satazi(n)
102 end do
103 end if
104
105 !-------------------------------------------------
106 ! [2.9] Distribute evenly across processors
107 !-------------------------------------------------
108
109 d_num_rad=num_tovs_after(i,myproc+1)
110
111 if (d_num_rad > 0) then
112 allocate (d_con_vars(d_num_rad))
113 allocate (d_con_vars_ad(d_num_rad))
114 allocate (d_aux_vars(d_num_rad))
115 allocate (d_tb(nchan,d_num_rad))
116 d_tb(:,:) = 0.0
117 d_con_vars(:) % nlevels = nlevels
118 do n = 1, d_num_rad
119 allocate (d_con_vars(n) % t(nlevels))
120 allocate (d_con_vars(n) % q(nlevels))
121 d_con_vars(n) % t(:) = 0.0
122 d_con_vars(n) % q(:) = 0.0
123 allocate (d_con_vars_ad(n) % t(nlevels))
124 allocate (d_con_vars_ad(n) % q(nlevels))
125 d_con_vars_ad(n) % t(:) = 0.0
126 d_con_vars_ad(n) % q(:) = 0.0
127 end do
128
129 ! Fill up with data that stays local
130
131 l_num_rad=min(num_rad,d_num_rad)
132
133 if (l_num_rad > 0) then
134 do n = 1, l_num_rad
135 d_con_vars(n) % t(:) = con_vars(n) % t(:)
136 d_con_vars(n) % q(:) = con_vars(n) % q(:)
137 ! d_aux_vars(n) % fastem(:) = 0.0
138 end do
139 d_con_vars(1:l_num_rad) % nlevels = con_vars(1:l_num_rad) % nlevels
140 d_con_vars(1:l_num_rad) % ps = con_vars(1:l_num_rad) % ps
141 d_aux_vars(1:l_num_rad) % t2m = aux_vars(1:l_num_rad) % t2m
142 d_aux_vars(1:l_num_rad) % q2m = aux_vars(1:l_num_rad) % q2m
143 d_aux_vars(1:l_num_rad) % u10 = aux_vars(1:l_num_rad) % u10
144 d_aux_vars(1:l_num_rad) % v10 = aux_vars(1:l_num_rad) % v10
145 d_aux_vars(1:l_num_rad) % surftype = aux_vars(1:l_num_rad) % surftype
146 d_aux_vars(1:l_num_rad) % surft = aux_vars(1:l_num_rad) % surft
147 d_aux_vars(1:l_num_rad) % satzen = aux_vars(1:l_num_rad) % satzen
148 d_aux_vars(1:l_num_rad) % satazi = aux_vars(1:l_num_rad) % satazi
149
150 d_tb(:,1:l_num_rad) = jo_grad_y%instid(i)%tb(:,1:l_num_rad)
151 end if
152
153 ! Get data from elsewhere
154 #ifdef DM_PARALLEL
155 do c=1,tovs_copy_count(i)
156 if (tovs_send_pe(i,c)==myproc) then
157 s=tovs_send_start(i,c)
158 len=tovs_send_count(i,c)
159 e=s+len-1
160 nl=len*nlevels
161 allocate(temp_t(nlevels,len))
162 do j=1,len
163 temp_t(:,j)=con_vars(s+j-1) % t(:)
164 end do
165 call mpi_send( temp_t,nl, true_mpi_real, tovs_recv_pe(i,c), &
166 c*8+1, comm, ierr)
167
168 allocate(temp_q(nlevels,len))
169 do j=1,len
170 temp_q(:,j)=con_vars(s+j-1) % q(:)
171 end do
172 call mpi_send (temp_q,nl, true_mpi_real, tovs_recv_pe(i,c), &
173 c*8+2, comm, ierr)
174
175 allocate (temp_tb(nchan,len))
176 temp_tb(:,:)=jo_grad_y%instid(i)%tb(:,s:e)
177 call mpi_send (temp_tb,len*nchan, &
178 true_mpi_real, tovs_recv_pe(i,c), c*8+3, comm, ierr)
179
180 allocate (rtemp(len,8))
181 rtemp(:,1)= con_vars(s:e) % ps
182 rtemp(:,2)= aux_vars(s:e) % t2m
183 rtemp(:,3)= aux_vars(s:e) % q2m
184 rtemp(:,4)= aux_vars(s:e) % u10
185 rtemp(:,5)= aux_vars(s:e) % v10
186 rtemp(:,6)= aux_vars(s:e) % surft
187 rtemp(:,7)= aux_vars(s:e) % satzen
188 rtemp(:,8)= aux_vars(s:e) % satazi
189 call mpi_send(rtemp,len*8, true_mpi_real, tovs_recv_pe(i,c), &
190 c*8+4, comm, ierr)
191
192 allocate (itemp(len,2))
193 itemp(:,1)= con_vars(s:e) % nlevels
194 itemp(:,2)= aux_vars(s:e) % surftype
195 call mpi_send(itemp,len*2, mpi_integer, tovs_recv_pe(i,c), &
196 c*8+5, comm, ierr)
197 end if
198
199 if (tovs_recv_pe(i,c)==myproc) then
200 s=tovs_recv_start(i,c)
201 len=tovs_send_count(i,c)
202 e=s+len-1
203 nl=len*nlevels
204 allocate(temp_t(nlevels,len))
205 call mpi_recv(temp_t,nl, true_mpi_real, tovs_send_pe(i,c), &
206 c*8+1, comm, status, ierr)
207 do j=1,len
208 d_con_vars(s+j-1) % t(:)=temp_t(:,j)
209 end do
210
211 allocate(temp_q(nlevels,len))
212 call mpi_recv(temp_q,nl, true_mpi_real, tovs_send_pe(i,c), &
213 c*8+2, comm, status, ierr)
214 do j=1,len
215 d_con_vars(s+j-1) % q(:)=temp_q(:,j)
216 end do
217
218 allocate (temp_tb(nchan,len))
219 call mpi_recv(temp_tb,len*nchan, true_mpi_real, &
220 tovs_send_pe(i,c), c*8+3, comm, status, ierr)
221 d_tb(:,s:e)=temp_tb(:,:)
222
223 allocate (rtemp(len,8))
224 call mpi_recv(rtemp,len*8, true_mpi_real, &
225 tovs_send_pe(i,c), c*8+4, comm, status, ierr)
226 d_con_vars(s:e) % ps = rtemp(:,1)
227 d_aux_vars(s:e) % t2m = rtemp(:,2)
228 d_aux_vars(s:e) % q2m = rtemp(:,3)
229 d_aux_vars(s:e) % u10 = rtemp(:,4)
230 d_aux_vars(s:e) % v10 = rtemp(:,5)
231 d_aux_vars(s:e) % surft = rtemp(:,6)
232 d_aux_vars(s:e) % satzen = rtemp(:,7)
233 d_aux_vars(s:e) % satazi = rtemp(:,8)
234
235 allocate (itemp(len,2))
236 call mpi_recv(itemp,len*2, mpi_integer, &
237 tovs_send_pe(i,c), c*8+5, comm, status, ierr)
238 d_con_vars(s:e) % nlevels = itemp(:,1)
239 d_aux_vars(s:e) % surftype = itemp(:,2)
240
241 end if
242
243 if (tovs_send_pe(i,c)==myproc) then
244 deallocate (temp_tb)
245 deallocate (temp_t)
246 deallocate (temp_q)
247 deallocate (itemp)
248 deallocate (rtemp)
249 end if
250 if (tovs_recv_pe(i,c)==myproc) then
251 deallocate (temp_tb)
252 deallocate (temp_t)
253 deallocate (temp_q)
254 deallocate (rtemp)
255 deallocate (itemp)
256 end if
257 end do
258 #endif
259
260 if (d_num_rad > 0) then
261 if (tovs_batch) then
262 call da_rttov_ad (i, nchan, d_num_rad, d_con_vars, &
263 d_aux_vars, d_con_vars_ad, d_tb)
264 else
265 do n=1,d_num_rad
266 call da_rttov_ad (i, nchan, 1, d_con_vars(n:n), &
267 d_aux_vars(n:n), d_con_vars_ad(n:n), d_tb(:,n:n))
268 end do
269 end if
270 end if
271
272 ! Transfer data back
273
274 if (l_num_rad > 0) then
275 do n = 1, l_num_rad
276 con_vars_ad(n) % t(:) = d_con_vars_ad(n) % t(:)
277 con_vars_ad(n) % q(:) = d_con_vars_ad(n) % q(:)
278 end do
279
280 con_vars_ad(1:l_num_rad) % ps = d_con_vars_ad(1:l_num_rad) % ps
281 end if
282 end if ! d_num_rad > 0
283
284 ! Return the data to other processors. Note the meaning of send_pe and
285 ! recv_pe is swapped here
286
287 nl=nlevels
288 #ifdef DM_PARALLEL
289 do c=1,tovs_copy_count(i)
290 if (tovs_recv_pe(i,c)==myproc) then
291 s=tovs_recv_start(i,c)
292 len=tovs_send_count(i,c)
293 e=s+len-1
294 nl=len*nlevels
295 allocate(temp_t(nlevels,len))
296 do j=1,len
297 temp_t(:,j)=d_con_vars_ad(s+j-1) % t(:)
298 end do
299 call mpi_send (temp_t,nl, true_mpi_real, tovs_send_pe(i,c), &
300 c*8+6, comm, ierr)
301
302 allocate(temp_q(nlevels,len))
303 do j=1,len
304 temp_q(:,j)=d_con_vars_ad(s+j-1) % q(:)
305 end do
306 call mpi_send (temp_q,nl, true_mpi_real, tovs_send_pe(i,c), &
307 c*8+7, comm, ierr)
308
309 allocate(rtemp(len,1))
310 rtemp(:,1) = d_con_vars_ad(s:e) % ps
311 call mpi_send (rtemp,len, true_mpi_real, tovs_send_pe(i,c), &
312 c*8+8, comm, ierr)
313 end if
314
315 if (tovs_send_pe(i,c)==myproc) then
316 s=tovs_send_start(i,c)
317 len=tovs_send_count(i,c)
318 e=s+len-1
319 nl=len*nlevels
320 allocate(temp_t(nlevels,len))
321 call mpi_recv (temp_t,nl, true_mpi_real, tovs_recv_pe(i,c), &
322 c*8+6, comm, status, ierr)
323 do j=1,len
324 con_vars_ad(s+j-1) % t(:)=temp_t(:,j)
325 end do
326
327 allocate(temp_q(nlevels,len))
328 call mpi_recv (temp_q,nl, true_mpi_real, tovs_recv_pe(i,c), &
329 c*8+7, comm, status, ierr)
330 do j=1,len
331 con_vars_ad(s+j-1) % q(:)=temp_q(:,j)
332 end do
333
334 allocate(rtemp(len,1))
335 call mpi_recv (rtemp,len, true_mpi_real, tovs_recv_pe(i,c), &
336 c*8+8, comm, status, ierr)
337 con_vars_ad(s:e) % ps=rtemp(:,1)
338 end if
339
340 if (tovs_recv_pe(i,c)==myproc) then
341 deallocate (temp_t)
342 deallocate (temp_q)
343 deallocate (rtemp)
344 end if
345 if (tovs_send_pe(i,c)==myproc) then
346 deallocate (temp_t)
347 deallocate (temp_q)
348 deallocate (rtemp)
349 end if
350 end do
351 #endif
352 if (d_num_rad > 0) then
353 do n=1,d_num_rad
354 deallocate (d_con_vars(n) % t)
355 deallocate (d_con_vars(n) % q)
356 deallocate (d_con_vars_ad(n) % t)
357 deallocate (d_con_vars_ad(n) % q)
358 end do
359
360 deallocate (d_tb)
361 deallocate (d_con_vars)
362 deallocate (d_aux_vars)
363 deallocate (d_con_vars_ad)
364 end if
365
366 ! adjoint of convert to hPa
367
368 if (num_rad > 0) then
369 con_vars_ad(:)% ps = con_vars_ad(:)%ps * 0.01
370 end if
371
372 do n=1,num_rad
373 ! 4.2 scale transform
374
375 do k=1, nlevels
376 model_mr(k) = con_vars_ad(n) % q(k)
377
378 if (iv%instid(i)%info%zk(k,n) <= 0.0) then
379 con_vars_ad(n)%t(k) = 0.0 !coefs(i) % ref_prfl_t(k,gas_id_watervapour)
380 con_vars_ad(n)%q(k) = 0.0 !coefs(i) % ref_prfl_mr(k,gas_id_watervapour)
381 else
382 ! adjoint of q(kg/kg) to ppmv
383
384 con_vars_ad(n)%q(k) = model_mr(k) * q2ppmv
385 ! call rttov_v2q (&
386 ! gas_unit_ppmv, &! in
387 ! model_mr(k), &! in
388 ! gas_id_watervapour, &! in
389 ! model_mr(k), &! in
390 ! con_vars_ad(n)%q(k) ) ! inout
391 end if
392 end do
393 end do
394
395 if (num_rad > 0) then
396 allocate(temp_t(nlevels,num_rad))
397 do n=1,num_rad
398 temp_t(:,n) = con_vars_ad(n)% t(:)
399 end do
400 call da_interp_lin_3d_adj (jo_grad_x%t, iv%instid(i)%info, temp_t)
401 deallocate(temp_t)
402
403 allocate(temp_q(nlevels,num_rad))
404 do n=1,num_rad
405 temp_q(:,n) = con_vars_ad(n)% q(:)
406 end do
407 call da_interp_lin_3d_adj (jo_grad_x%q, iv%instid(i)%info, temp_q)
408 deallocate(temp_q)
409
410 allocate(temp_ps(num_rad))
411 do n=1,num_rad
412 temp_ps(n) = con_vars_ad(n)% ps
413 end do
414 call da_interp_lin_2d_adj (jo_grad_x% psfc, iv%instid(i)%info, 1, temp_ps)
415 deallocate(temp_ps)
416 end if
417
418 if (num_rad > 0) then
419 do n = 1, num_rad
420 deallocate (con_vars(n) % t)
421 deallocate (con_vars(n) % q)
422 deallocate (con_vars_ad(n) % t)
423 deallocate (con_vars_ad(n) % q)
424 end do
425
426 deallocate (con_vars)
427 deallocate (aux_vars)
428 deallocate (con_vars_ad)
429 deallocate (model_mr)
430 end if
431 end do ! end loop for sensor
432
433 if (trace_use) call da_trace_exit("da_transform_xtoy_rttov_adj")
434 #else
435 call da_error(__FILE__,__LINE__, &
436 (/"Must compile with $RTTOV option for radiances"/))
437 #endif
438
439 end subroutine da_transform_xtoy_rttov_adj
440