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