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