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