1 subroutine da_transform_xtoy_rttov (grid, iv, y )
3 !---------------------------------------------------------------------------
4 ! Purpose: transform from analysis increment to
5 ! pertubation radiance.
7 ! METHOD: delta_y = H delta_x
8 ! 1. input reference state of RTTOV_TL
9 ! 2. interpolate analysis increment to obs location
11 !---------------------------------------------------------------------------
15 type (domain), intent(in) :: grid
16 type (y_type), intent(inout) :: y ! H' delta_x
17 type (iv_type), intent(in) :: iv ! O-B structure.
21 integer :: i,j,k,c ! Index dimension.
22 integer :: nlevels ! Number of obs levels.
24 real, allocatable :: model_t(:,:) ! Model value t at ob location.
25 real, allocatable :: model_q(:,:) ! Model value q(kg/kg) at ob location.
26 real, allocatable :: model_mr(:,:) ! Model value mr(ppmv) at ob location.
27 real, allocatable :: model_psfc(:)
29 integer :: num_rad, nchan, n , l
31 type(con_vars_type), allocatable :: con_vars(:), con_vars_tl(:)
32 type(aux_vars_type), allocatable :: aux_vars(:), aux_vars_tl(:)
34 ! evenly distributed versions
35 type(con_vars_type), allocatable :: d_con_vars(:), d_con_vars_tl(:)
36 type(aux_vars_type), allocatable :: d_aux_vars(:), d_aux_vars_tl(:)
37 real, allocatable :: d_tb(:,:)
39 integer :: d_num_rad, l_num_rad,e,s,len,nl
40 real, allocatable :: temp_t(:,:),temp_q(:,:)
41 real, allocatable :: temp_t_tl(:,:),temp_q_tl(:,:),rtemp(:,:)
42 real, allocatable :: temp_tb(:,:)
43 real, allocatable :: temp_t_jac(:,:,:), temp_q_jac(:,:,:), temp_ps_jac(:,:)
44 integer, allocatable :: itemp(:,:)
48 integer :: status(mpi_status_size) ! MPI status.
51 if (iv%num_inst < 1) return
53 if (trace_use) call da_trace_entry("da_transform_xtoy_rttov")
55 do i= 1, iv%num_inst ! loop for sensor
57 if (iv%instid(i)%num_rad_glo < 1) cycle
59 num_rad = iv%instid(i)%num_rad
60 nlevels = iv%instid(i)%nlevels
61 nchan = iv%instid(i)%nchan
64 allocate (con_vars(num_rad))
65 allocate (con_vars_tl(num_rad))
66 allocate (aux_vars(num_rad))
67 allocate (aux_vars_tl(num_rad))
69 allocate(model_t(nlevels,num_rad))
70 allocate(model_q(nlevels,num_rad))
71 allocate(model_mr(nlevels,num_rad))
72 allocate(model_psfc(num_rad))
79 call da_interp_lin_3d (grid%xa%t, iv%instid(i)%info, model_t)
80 call da_interp_lin_3d (grid%xa%q, iv%instid(i)%info, model_q)
82 call da_interp_lin_2d (grid%xa%psfc, iv%instid(i)%info, 1, model_psfc)
84 model_psfc(:) = 0.01*model_psfc(:) ! convert to hPa
85 con_vars(:) % nlevels = nlevels
86 con_vars_tl(:) % nlevels = nlevels
87 aux_vars_tl(:) % t2m = 0.0
88 aux_vars_tl(:) % q2m = 0.0
89 aux_vars_tl(:) % u10 = 0.0
90 aux_vars_tl(:) % v10 = 0.0
91 ! aux_vars(:) % fastem(:) = 0.0
92 aux_vars_tl(:) % surft = 0.0
93 ! aux_vars_tl(:) % fastem(:) = 0.0
97 if ( iv%instid(i)%info%zk(k,n) <= 0.0 ) then
98 model_t(k,n) = 0.0 !coefs(i) % ref_prfl_t(k,gas_id_watervapour)
99 model_mr(k,n) = 0.0 !coefs(i) % ref_prfl_mr(k,gas_id_watervapour)
101 ! model_mr(k,n) = model_q(k,n) * 1.60771704e+6 ! convert q(kg/kg) to ppmv
102 model_mr(k,n) = model_q(k,n) * q2ppmv
106 allocate (con_vars(n) % t(nlevels))
107 allocate (con_vars(n) % q(nlevels))
108 con_vars(n) % t(:) = iv%instid(i)%t (:,n)
109 con_vars(n) % q(:) = iv%instid(i)%mr(:,n)
110 con_vars(n) % ps = iv%instid(i)%ps(n)
112 if ( use_rttov_kmatrix ) then
113 allocate (con_vars(n) % t_jac(nchan,nlevels))
114 allocate (con_vars(n) % q_jac(nchan,nlevels))
115 allocate (con_vars(n) % ps_jac(nchan))
117 con_vars(n) % t_jac(:,k) = iv%instid(i)%t_jacobian(:,k,n)
118 con_vars(n) % q_jac(:,k) = iv%instid(i)%q_jacobian(:,k,n)
120 con_vars(n) % ps_jac(:)= iv%instid(i)%ps_jacobian(:,n)
123 allocate (con_vars_tl(n) % t(nlevels))
124 allocate (con_vars_tl(n) % q(nlevels))
126 con_vars_tl(n) % t(:) = model_t (:,n)
127 con_vars_tl(n) % q(:) = model_mr (:,n)
128 con_vars_tl(n) % ps = model_psfc(n)
130 aux_vars(n) % t2m = iv%instid(i)%t2m(n)
131 aux_vars(n) % q2m = iv%instid(i)%mr2m(n)
132 aux_vars(n) % u10 = iv%instid(i)%u10(n)
133 aux_vars(n) % v10 = iv%instid(i)%v10(n)
134 aux_vars(n) % surftype = iv%instid(i)%surftype(n)
135 aux_vars(n) % surft = iv%instid(i)%ts(n)
136 aux_vars(n) % satzen = iv%instid(i)%satzen(n)
137 aux_vars(n) % satazi = iv%instid(i)%satazi(n)
138 aux_vars(n) % solzen = iv%instid(i)%solzen(n)
139 aux_vars(n) % solazi = iv%instid(i)%solazi(n)
140 aux_vars(n) % elevation = iv%instid(i)%elevation(n)
141 aux_vars(n) % rlat = iv%instid(i)%info%lat(1,n)
148 deallocate(model_psfc)
151 ! [1.3] Call RTM TL model
153 d_num_rad=num_tovs_after(i,myproc+1)
155 if (d_num_rad > 0) then
156 allocate (d_con_vars(d_num_rad))
157 allocate (d_con_vars_tl(d_num_rad))
158 allocate (d_aux_vars(d_num_rad))
159 allocate (d_aux_vars_tl(d_num_rad))
160 allocate (d_tb(nchan,d_num_rad))
161 d_con_vars(:) % nlevels = nlevels
164 allocate (d_con_vars(n) % t(nlevels))
165 allocate (d_con_vars(n) % q(nlevels))
167 if ( use_rttov_kmatrix ) then
168 allocate (d_con_vars(n) % t_jac(nchan,nlevels))
169 allocate (d_con_vars(n) % q_jac(nchan,nlevels))
170 allocate (d_con_vars(n) % ps_jac(nchan))
173 allocate (d_con_vars_tl(n) % t(nlevels))
174 allocate (d_con_vars_tl(n) % q(nlevels))
178 ! Fill up with data that stays local
179 l_num_rad=Min(num_rad,d_num_rad)
181 if (l_num_rad > 0) then
182 d_con_vars(1:l_num_rad) % nlevels = con_vars(1:l_num_rad) % nlevels
183 d_con_vars_tl(1:l_num_rad) % nlevels = con_vars_tl(1:l_num_rad) % nlevels
186 d_con_vars(n) % t(:) = con_vars(n) % t(:)
187 d_con_vars(n) % q(:) = con_vars(n) % q(:)
188 if ( use_rttov_kmatrix ) then
189 d_con_vars(n) % t_jac(:,:) = con_vars(n) % t_jac(:,:)
190 d_con_vars(n) % q_jac(:,:) = con_vars(n) % q_jac(:,:)
191 d_con_vars(n) % ps_jac(:) = con_vars(n) % ps_jac(:)
193 d_con_vars_tl(n) % t(:) = con_vars_tl(n) % t(:)
194 d_con_vars_tl(n) % q(:) = con_vars_tl(n) % q(:)
195 ! d_aux_vars(n) % fastem(:) = 0.0
198 d_con_vars(1:l_num_rad) % ps = con_vars(1:l_num_rad) % ps
199 d_con_vars_tl(1:l_num_rad) % ps = con_vars_tl(1:l_num_rad) % ps
201 d_aux_vars(1:l_num_rad) % t2m = aux_vars(1:l_num_rad) % t2m
202 d_aux_vars(1:l_num_rad) % q2m = aux_vars(1:l_num_rad) % q2m
203 d_aux_vars(1:l_num_rad) % u10 = aux_vars(1:l_num_rad) % u10
204 d_aux_vars(1:l_num_rad) % v10 = aux_vars(1:l_num_rad) % v10
205 d_aux_vars(1:l_num_rad) % surftype = aux_vars(1:l_num_rad) % surftype
206 d_aux_vars(1:l_num_rad) % surft = aux_vars(1:l_num_rad) % surft
207 ! d_aux_vars(1:l_num_rad) % fastem(:) = aux_vars(1:l_num_rad) % fastem(:)
208 d_aux_vars(1:l_num_rad) % satzen = aux_vars(1:l_num_rad) % satzen
209 d_aux_vars(1:l_num_rad) % satazi = aux_vars(1:l_num_rad) % satazi
210 d_aux_vars(1:l_num_rad) % solzen = aux_vars(1:l_num_rad) % solzen
211 d_aux_vars(1:l_num_rad) % solazi = aux_vars(1:l_num_rad) % solazi
212 d_aux_vars(1:l_num_rad) % elevation = aux_vars(1:l_num_rad) % elevation
213 d_aux_vars(1:l_num_rad) % rlat = aux_vars(1:l_num_rad) % rlat
215 d_aux_vars_tl(1:l_num_rad) % t2m = aux_vars_tl(1:l_num_rad) % t2m
216 d_aux_vars_tl(1:l_num_rad) % q2m = aux_vars_tl(1:l_num_rad) % q2m
217 d_aux_vars_tl(1:l_num_rad) % u10 = aux_vars_tl(1:l_num_rad) % u10
218 d_aux_vars_tl(1:l_num_rad) % v10 = aux_vars_tl(1:l_num_rad) % v10
219 d_aux_vars_tl(1:l_num_rad) % surft = aux_vars_tl(1:l_num_rad) % surft
220 ! d_aux_vars_tl(1:l_numrad) % fastem(:) = aux_vars_tl(1:l_num_rad) % fastem(:)
224 ! Get data from elsewhere
226 do c=1,tovs_copy_count(i)
227 if (tovs_send_pe(i,c)==myproc) then
228 s=tovs_send_start(i,c)
229 len=tovs_send_count(i,c)
233 allocate(temp_t(nlevels,len))
235 temp_t(:,j)=con_vars(s+j-1) % t(:)
237 call mpi_send( temp_t,nl, true_mpi_real, tovs_recv_pe(i,c), &
240 allocate(temp_t_tl(nlevels,len))
242 temp_t_tl(:,j)=con_vars_tl(s+j-1) % t(:)
244 call mpi_send( temp_t_tl,nl, true_mpi_real, tovs_recv_pe(i,c), &
247 allocate(temp_q(nlevels,len))
249 temp_q(:,j)=con_vars(s+j-1) % q(:)
251 call mpi_send (temp_q,nl, true_mpi_real, tovs_recv_pe(i,c), &
254 allocate(temp_q_tl(nlevels,len))
256 temp_q_tl(:,j)=con_vars_tl(s+j-1) % q(:)
258 call mpi_send (temp_q_tl,nl, true_mpi_real, tovs_recv_pe(i,c), &
261 allocate (rtemp(len,18))
262 rtemp(:,1)= con_vars(s:e) % ps
263 rtemp(:,2)= con_vars_tl(s:e) % ps
264 rtemp(:,3)= aux_vars(s:e) % t2m
265 rtemp(:,4)= aux_vars(s:e) % q2m
266 rtemp(:,5)= aux_vars(s:e) % u10
267 rtemp(:,6)= aux_vars(s:e) % v10
268 rtemp(:,7)= aux_vars(s:e) % surft
269 rtemp(:,8)= aux_vars(s:e) % satzen
270 rtemp(:,9)= aux_vars(s:e) % satazi
271 rtemp(:,10)= aux_vars(s:e) % solzen
272 rtemp(:,11)= aux_vars(s:e) % solazi
273 rtemp(:,12)= aux_vars(s:e) % elevation
274 rtemp(:,13)= aux_vars(s:e) % rlat
275 rtemp(:,14)= aux_vars_tl(s:e) % t2m
276 rtemp(:,15)= aux_vars_tl(s:e) % q2m
277 rtemp(:,16)= aux_vars_tl(s:e) % u10
278 rtemp(:,17)= aux_vars_tl(s:e) % v10
279 rtemp(:,18)= aux_vars_tl(s:e) % surft
280 call mpi_send (rtemp,len*18, true_mpi_real, tovs_recv_pe(i,c), &
283 allocate (itemp(len,2))
284 itemp(:,1)= con_vars(s:e) % nlevels ! aux_vars_tl identical
285 itemp(:,2)= aux_vars(s:e) % surftype ! aux_vars_tl identical
286 call mpi_send (itemp,len*2, mpi_integer, tovs_recv_pe(i,c), &
291 deallocate (temp_t_tl)
292 deallocate (temp_q_tl)
296 if ( use_rttov_kmatrix ) then
297 allocate(temp_t_jac(nchan,nlevels,len))
300 temp_t_jac(:,k,j)=con_vars(s+j-1) % t_jac(:,k)
303 call mpi_send( temp_t_jac, nl*nchan, true_mpi_real, tovs_recv_pe(i,c), &
306 allocate(temp_q_jac(nchan,nlevels,len))
309 temp_q_jac(:,k,j)=con_vars(s+j-1) % q_jac(:,k)
312 call mpi_send( temp_q_jac, nl*nchan, true_mpi_real, tovs_recv_pe(i,c), &
315 allocate(temp_ps_jac(nchan,len))
317 temp_ps_jac(:,j)=con_vars(s+j-1) % ps_jac(:)
319 call mpi_send( temp_ps_jac, len*nchan, true_mpi_real, tovs_recv_pe(i,c), &
322 deallocate (temp_t_jac)
323 deallocate (temp_q_jac)
324 deallocate (temp_ps_jac)
327 if (tovs_recv_pe(i,c)==myproc) then
328 s=tovs_recv_start(i,c)
329 len=tovs_send_count(i,c)
333 allocate(temp_t(nlevels,len))
334 call mpi_recv (temp_t,nl, true_mpi_real, tovs_send_pe(i,c), &
335 c*10+1, comm, status, ierr)
337 d_con_vars(s+j-1) % t(:)=temp_t(:,j)
340 allocate(temp_t_tl(nlevels,len))
341 call mpi_recv (temp_t_tl,nl, true_mpi_real, tovs_send_pe(i,c), &
342 c*10+2, comm, status, ierr)
344 d_con_vars_tl(s+j-1) % t(:)=temp_t_tl(:,j)
347 allocate(temp_q(nlevels,len))
348 call mpi_recv (temp_q,nl, true_mpi_real, tovs_send_pe(i,c), &
349 c*10+3, comm, status, ierr)
351 d_con_vars(s+j-1) % q(:)=temp_q(:,j)
354 allocate(temp_q_tl(nlevels,len))
355 call mpi_recv (temp_q_tl,nl, true_mpi_real, tovs_send_pe(i,c), &
356 c*10+4, comm, status, ierr)
358 d_con_vars_tl(s+j-1) % q(:)=temp_q_tl(:,j)
361 allocate (rtemp(len,18))
362 call mpi_recv (rtemp,len*18, true_mpi_real, tovs_send_pe(i,c), &
363 c*10+5, comm, status, ierr)
364 d_con_vars(s:e) % ps = rtemp(:,1)
365 d_con_vars_tl(s:e) % ps = rtemp(:,2)
366 d_aux_vars(s:e) % t2m = rtemp(:,3)
367 d_aux_vars(s:e) % q2m = rtemp(:,4)
368 d_aux_vars(s:e) % u10 = rtemp(:,5)
369 d_aux_vars(s:e) % v10 = rtemp(:,6)
370 d_aux_vars(s:e) % surft = rtemp(:,7)
371 d_aux_vars(s:e) % satzen = rtemp(:,8)
372 d_aux_vars(s:e) % satazi = rtemp(:,9)
373 d_aux_vars(s:e) % solzen = rtemp(:,10)
374 d_aux_vars(s:e) % solazi = rtemp(:,11)
375 d_aux_vars(s:e) % elevation = rtemp(:,12)
376 d_aux_vars(s:e) % rlat = rtemp(:,13)
377 d_aux_vars_tl(s:e) % t2m = rtemp(:,14)
378 d_aux_vars_tl(s:e) % q2m = rtemp(:,15)
379 d_aux_vars_tl(s:e) % u10 = rtemp(:,16)
380 d_aux_vars_tl(s:e) % v10 = rtemp(:,17)
381 d_aux_vars_tl(s:e) % surft = rtemp(:,18)
383 allocate (itemp(len,2))
384 call mpi_recv (itemp,len*2, mpi_integer, tovs_send_pe(i,c), &
385 c*10+6, comm, status, ierr)
387 d_con_vars(s:e) % nlevels = itemp(:,1)
388 d_aux_vars(s:e) % surftype = itemp(:,2)
390 d_con_vars_tl(s:e) % nlevels = d_con_vars(s:e) % nlevels
391 d_aux_vars_tl(s:e) % surftype = d_aux_vars(s:e) % surftype
395 deallocate (temp_t_tl)
396 deallocate (temp_q_tl)
400 if ( use_rttov_kmatrix ) then
401 allocate(temp_t_jac(nchan,nlevels,len))
402 call mpi_recv (temp_t_jac,nl*nchan, true_mpi_real, tovs_send_pe(i,c), &
403 c*10+7, comm, status, ierr)
406 d_con_vars(s+j-1) % t_jac(:,k) = temp_t_jac(:,k,j)
410 allocate(temp_q_jac(nchan,nlevels,len))
411 call mpi_recv (temp_q_jac,nl*nchan, true_mpi_real, tovs_send_pe(i,c), &
412 c*10+8, comm, status, ierr)
415 d_con_vars(s+j-1) % q_jac(:,k) = temp_q_jac(:,k,j)
419 allocate(temp_ps_jac(nchan,len))
420 call mpi_recv (temp_ps_jac,len*nchan, true_mpi_real, tovs_send_pe(i,c), &
421 c*10+9, comm, status, ierr)
423 d_con_vars(s+j-1) % ps_jac(:) = temp_ps_jac(:,j)
426 deallocate (temp_t_jac)
427 deallocate (temp_q_jac)
428 deallocate (temp_ps_jac)
433 if (d_num_rad > 0) then
436 call da_rttov_tl (i, nchan, d_num_rad, d_con_vars, &
437 d_aux_vars, d_con_vars_tl, d_aux_vars_tl, d_tb)
440 if (.not. use_rttov_kmatrix) then
444 call da_rttov_tl (i, nchan, 1, d_con_vars(n:n), &
445 d_aux_vars(n:n), d_con_vars_tl(n:n), d_aux_vars_tl(n:n), &
448 !$OMP END PARALLEL DO
454 d_tb(l,n) = d_tb(l,n) + &
455 d_con_vars(n)%t_jac(l,k) * d_con_vars_tl(n)%t(k) + &
456 d_con_vars(n)%q_jac(l,k) * d_con_vars_tl(n)%q(k)
458 d_tb(l,n) = d_tb(l,n) + &
459 d_con_vars(n)%ps_jac(l) * d_con_vars_tl(n)%ps
468 ! Return the local data
469 if (l_num_rad > 0) then
470 y%instid(i)%tb(:,1:l_num_rad) = d_tb(:,1:l_num_rad)
473 ! Return the data to other processors. Note the meaning of send_pe and recv_pe is
477 do c=1,tovs_copy_count(i)
478 if (tovs_recv_pe(i,c)==myproc) then
479 s=tovs_recv_start(i,c)
480 len=tovs_send_count(i,c)
482 allocate(temp_tb(nchan,len))
483 temp_tb(:,:) = d_tb(:,s:e)
484 call mpi_send (temp_tb,len*nchan, true_mpi_real, tovs_send_pe(i,c), c*10+10, comm, ierr)
487 if (tovs_send_pe(i,c)==myproc) then
488 s=tovs_send_start(i,c)
489 len=tovs_send_count(i,c)
491 allocate(temp_tb(nchan,len))
492 call mpi_recv (temp_tb,len*nchan, true_mpi_real, &
493 tovs_recv_pe(i,c), c*10+10, comm, status, ierr)
494 y%instid(i)%tb(:,s:e)=temp_tb(:,:)
500 if (d_num_rad > 0) then
502 deallocate (d_con_vars(n) % t)
503 deallocate (d_con_vars(n) % q)
504 deallocate (d_con_vars_tl(n) % t)
505 deallocate (d_con_vars_tl(n) % q)
506 if ( use_rttov_kmatrix ) then
507 deallocate (d_con_vars(n) % t_jac)
508 deallocate (d_con_vars(n) % q_jac)
509 deallocate (d_con_vars(n) % ps_jac)
514 deallocate (d_con_vars)
515 deallocate (d_aux_vars)
516 deallocate (d_aux_vars_tl)
517 deallocate (d_con_vars_tl)
520 !-------------------------------------------------------------------
522 !-------------------------------------------------------------------
524 if (num_rad > 0) then
526 deallocate (con_vars(n) % t)
527 deallocate (con_vars(n) % q)
528 deallocate (con_vars_tl(n) % t)
529 deallocate (con_vars_tl(n) % q)
530 if ( use_rttov_kmatrix ) then
531 deallocate (con_vars(n) % t_jac)
532 deallocate (con_vars(n) % q_jac)
533 deallocate (con_vars(n) % ps_jac)
536 deallocate (con_vars)
537 deallocate (aux_vars)
538 deallocate (con_vars_tl)
539 deallocate (aux_vars_tl)
541 end do ! end loop for sensor
543 if (trace_use) call da_trace_exit("da_transform_xtoy_rttov")
545 call da_error(__FILE__,__LINE__, &
546 (/"Must compile with $RTTOV option for radiances"/))
549 end subroutine da_transform_xtoy_rttov