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 real, allocatable :: temp_tb(:,:)
44 integer, allocatable :: itemp(:,:)
45
46 #ifdef DM_PARALLEL
47 integer :: status(mpi_status_size) ! MPI status.
48 #endif
49
50 #if !defined(RTTOV)
51 call da_error(__FILE__,__LINE__, &
52 (/"Must compile with $RTTOV option for radiances"/))
53 #else
54
55 if (iv%num_inst < 1) return
56
57 if (trace_use) call da_trace_entry("da_transform_xtoy_rttov")
58
59 do i= 1, iv%num_inst ! loop for sensor
60
61 if (iv%instid(i)%num_rad_glo < 1) cycle
62
63 num_rad = iv%instid(i)%num_rad
64 nlevels = iv%instid(i)%nlevels
65 nchan = iv%instid(i)%nchan
66
67 if (num_rad > 0) then
68 allocate (con_vars(num_rad))
69 allocate (con_vars_tl(num_rad))
70 allocate (aux_vars(num_rad))
71 allocate (aux_vars_tl(num_rad))
72
73 allocate(model_t(nlevels,num_rad))
74 allocate(model_q(nlevels,num_rad))
75 allocate(model_mr(nlevels,num_rad))
76 allocate(model_psfc(num_rad))
77 model_t(:,:)=0.0
78 model_q(:,:)=0.0
79 model_mr(:,:)=0.0
80 model_psfc(:)=0.0
81
82 call da_interp_lin_3d_new( xa%t, xp, &
83 iv%instid(i)%loc_i, &
84 iv%instid(i)%loc_j, &
85 iv%instid(i)%loc_k, &
86 iv%instid(i)%loc_dx, &
87 iv%instid(i)%loc_dy, &
88 iv%instid(i)%loc_dz, &
89 iv%instid(i)%loc_dxm, &
90 iv%instid(i)%loc_dym, &
91 iv%instid(i)%loc_dzm, &
92 model_t, nlevels, num_rad)
93
94 call da_interp_lin_3d_new (xa%q, xp, &
95 iv%instid(i)%loc_i, &
96 iv%instid(i)%loc_j, &
97 iv%instid(i)%loc_k, &
98 iv%instid(i)%loc_dx, &
99 iv%instid(i)%loc_dy, &
100 iv%instid(i)%loc_dz, &
101 iv%instid(i)%loc_dxm, &
102 iv%instid(i)%loc_dym, &
103 iv%instid(i)%loc_dzm, &
104 model_q, nlevels, num_rad)
105
106 call da_interp_lin_2d_new (xa % psfc, xp%ims, xp%ime, xp%jms, xp%jme, &
107 iv%instid(i)%loc_i, iv%instid(i)%loc_j, &
108 iv%instid(i)%loc_dx, iv%instid(i)%loc_dy, &
109 iv%instid(i)%loc_dxm, iv%instid(i)%loc_dym, &
110 model_psfc,num_rad)
111
112 model_psfc(:) = 0.01*model_psfc(:) ! convert to hPa
113 con_vars(:) % nlevels = nlevels
114 con_vars_tl(:) % nlevels = nlevels
115 aux_vars_tl(:) % t2m = 0.0
116 aux_vars_tl(:) % q2m = 0.0
117 aux_vars_tl(:) % u10 = 0.0
118 aux_vars_tl(:) % v10 = 0.0
119 ! aux_vars(:) % fastem(:) = 0.
120 aux_vars_tl(:) % surft = 0.0
121 ! aux_vars_tl(:) % fastem(:) = 0.0
122
123 do n=1,num_rad
124 do k=1, nlevels
125 if ( iv%instid(i)%zk(k,n) <= 0.0 ) then
126 model_t(k,n) = 0.0 !coefs(i) % ref_prfl_t(k,gas_id_watervapour)
127 model_mr(k,n) = 0.0 !coefs(i) % ref_prfl_mr(k,gas_id_watervapour)
128 else
129 ! model_mr(k,n) = model_q(k,n) * 1.60771704e+6 ! convert q(kg/kg) to ppmv
130 model_mr(k,n) = model_q(k,n) * q2ppmv
131 ! call rttov_q2v ( &
132 ! gas_unit_specconc, &! in
133 ! model_q(k,n), &! in
134 ! gas_id_watervapour, &! in
135 ! model_q(k,n), &! in
136 ! model_mr(k,n) ) ! inout
137 end if
138 end do
139
140 allocate (con_vars(n) % t(nlevels))
141 allocate (con_vars(n) % q(nlevels))
142 con_vars(n) % t(:) = iv%instid(i)%t (:,n)
143 con_vars(n) % q(:) = iv%instid(i)%mr(:,n)
144 con_vars(n) % ps = iv%instid(i)%ps(n)
145
146 allocate (con_vars_tl(n) % t(nlevels))
147 allocate (con_vars_tl(n) % q(nlevels))
148
149 con_vars_tl(n) % t(:) = model_t (:,n)
150 con_vars_tl(n) % q(:) = model_mr (:,n)
151 con_vars_tl(n) % ps = model_psfc(n)
152
153 aux_vars(n) % t2m = iv%instid(i)%t2m(n)
154 aux_vars(n) % q2m = iv%instid(i)%mr2m(n)
155 aux_vars(n) % u10 = iv%instid(i)%u10(n)
156 aux_vars(n) % v10 = iv%instid(i)%v10(n)
157 aux_vars(n) % surftype = iv%instid(i)%landsea_mask(n)
158 aux_vars(n) % surft = iv%instid(i)%ts(n)
159 aux_vars(n) % satzen = iv%instid(i)%satzen(n)
160 aux_vars(n) % satazi = iv%instid(i)%satazi(n)
161
162 aux_vars_tl(n) % surftype = iv%instid(i)%landsea_mask(n)
163 aux_vars_tl(n) % satzen = iv%instid(i)%satzen(n)
164 aux_vars_tl(n) % satazi = iv%instid(i)%satazi(n)
165 end do
166
167 deallocate(model_t)
168 deallocate(model_q)
169 deallocate(model_mr)
170 deallocate(model_psfc)
171 end if
172
173 ! [1.3] Call RTM TL model
174
175 d_num_rad=num_tovs_after(i,myproc+1)
176
177 if (d_num_rad > 0) then
178 allocate (d_con_vars(d_num_rad))
179 allocate (d_con_vars_tl(d_num_rad))
180 allocate (d_aux_vars(d_num_rad))
181 allocate (d_aux_vars_tl(d_num_rad))
182 allocate (d_tb(nchan,d_num_rad))
183
184 d_con_vars(:) % nlevels = nlevels
185
186 do n = 1, d_num_rad
187 allocate (d_con_vars(n) % t(nlevels))
188 allocate (d_con_vars(n) % q(nlevels))
189
190 allocate (d_con_vars_tl(n) % t(nlevels))
191 allocate (d_con_vars_tl(n) % q(nlevels))
192 end do
193 end if
194
195 ! Fill up with data that stays local
196 l_num_rad=Min(num_rad,d_num_rad)
197
198 if (l_num_rad > 0) then
199 d_con_vars(1:l_num_rad) % nlevels = con_vars(1:l_num_rad) % nlevels
200 d_con_vars_tl(1:l_num_rad) % nlevels = con_vars_tl(1:l_num_rad) % nlevels
201
202 do n = 1, l_num_rad
203 d_con_vars(n) % t(:) = con_vars(n) % t(:)
204 d_con_vars(n) % q(:) = con_vars(n) % q(:)
205 d_con_vars_tl(n) % t(:) = con_vars_tl(n) % t(:)
206 d_con_vars_tl(n) % q(:) = con_vars_tl(n) % q(:)
207 ! d_aux_vars(n) % fastem(:) = 0.
208 end do
209
210 d_con_vars(1:l_num_rad) % ps = con_vars(1:l_num_rad) % ps
211 d_con_vars_tl(1:l_num_rad) % ps = con_vars_tl(1:l_num_rad) % ps
212
213 d_aux_vars(1:l_num_rad) % t2m = aux_vars(1:l_num_rad) % t2m
214 d_aux_vars(1:l_num_rad) % q2m = aux_vars(1:l_num_rad) % q2m
215 d_aux_vars(1:l_num_rad) % u10 = aux_vars(1:l_num_rad) % u10
216 d_aux_vars(1:l_num_rad) % v10 = aux_vars(1:l_num_rad) % v10
217 d_aux_vars(1:l_num_rad) % surftype = aux_vars(1:l_num_rad) % surftype
218 d_aux_vars(1:l_num_rad) % surft = aux_vars(1:l_num_rad) % surft
219 ! d_aux_vars(1:l_num_rad) % fastem(:) = aux_vars(1:l_num_rad) % fastem(:)
220 d_aux_vars(1:l_num_rad) % satzen = aux_vars(1:l_num_rad) % satzen
221 d_aux_vars(1:l_num_rad) % satazi = aux_vars(1:l_num_rad) % satazi
222
223 d_aux_vars_tl(1:l_num_rad) % t2m = aux_vars_tl(1:l_num_rad) % t2m
224 d_aux_vars_tl(1:l_num_rad) % q2m = aux_vars_tl(1:l_num_rad) % q2m
225 d_aux_vars_tl(1:l_num_rad) % u10 = aux_vars_tl(1:l_num_rad) % u10
226 d_aux_vars_tl(1:l_num_rad) % v10 = aux_vars_tl(1:l_num_rad) % v10
227 d_aux_vars_tl(1:l_num_rad) % surftype = aux_vars_tl(1:l_num_rad) % surftype
228 d_aux_vars_tl(1:l_num_rad) % surft = aux_vars_tl(1:l_num_rad) % surft
229 ! d_aux_vars_tl(1:l_numrad) % fastem(:) = aux_vars_tl(1:l_num_rad) % fastem(:)
230 d_aux_vars_tl(1:l_num_rad) % satzen = aux_vars_tl(1:l_num_rad) % satzen
231 d_aux_vars_tl(1:l_num_rad) % satazi = aux_vars_tl(1:l_num_rad) % satazi
232 end if
233
234 ! Get data from elsewhere
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 if (d_num_rad > 0) then
377 if (tovs_batch) then
378 call da_rttov_tl (i, nchan, d_num_rad, d_con_vars, &
379 d_aux_vars, d_con_vars_tl, d_aux_vars_tl, d_tb)
380 else
381 do n=1,d_num_rad
382 call da_rttov_tl (i, nchan, 1, d_con_vars(n:n), &
383 d_aux_vars(n:n), d_con_vars_tl(n:n), d_aux_vars_tl(n:n), &
384 d_tb(:,n:n))
385 end do
386 end if
387 end if
388
389 ! Transfer data back
390
391 ! Return the local data
392 if (l_num_rad > 0) then
393 y%instid(i)%tb(:,1:l_num_rad) = d_tb(:,1:l_num_rad)
394 end if
395
396 ! Return the data to other processors. Note the meaning of send_pe and recv_pe is
397 ! swapped here
398
399 nl=nlevels
400 #ifdef DM_PARALLEL
401 do c=1,tovs_copy_count(i)
402 if (tovs_recv_pe(i,c)==myproc) then
403 s=tovs_recv_start(i,c)
404 len=tovs_send_count(i,c)
405 e=s+len-1
406 nl=len*nlevels
407 allocate(temp_tb(nchan,len))
408 temp_tb(:,:) = d_tb(:,s:e)
409 call mpi_send (temp_tb,len*nchan, true_mpi_real, tovs_send_pe(i,c), c*7+7, comm, ierr)
410 deallocate(temp_tb)
411 end if
412 if (tovs_send_pe(i,c)==myproc) then
413 s=tovs_send_start(i,c)
414 len=tovs_send_count(i,c)
415 e=s+len-1
416 nl=len*nlevels
417 allocate(temp_tb(nchan,len))
418 call mpi_recv (temp_tb,len*nchan, true_mpi_real, &
419 tovs_recv_pe(i,c), c*7+7, comm, status, ierr)
420 y%instid(i)%tb(:,s:e)=temp_tb(:,:)
421 deallocate(temp_tb)
422 end if
423 end do
424 #endif
425
426 if (d_num_rad > 0) then
427 do n=1,d_num_rad
428 deallocate (d_con_vars(n) % t)
429 deallocate (d_con_vars(n) % q)
430 deallocate (d_con_vars_tl(n) % t)
431 deallocate (d_con_vars_tl(n) % q)
432 end do
433
434 deallocate (d_tb)
435 deallocate (d_con_vars)
436 deallocate (d_aux_vars)
437 deallocate (d_aux_vars_tl)
438 deallocate (d_con_vars_tl)
439 end if
440
441 !-------------------------------------------------------------------
442 ! [2.0] assign Hdx :
443 !-------------------------------------------------------------------
444
445 if (num_rad > 0) then
446 do n=1,num_rad
447 deallocate (con_vars(n) % t)
448 deallocate (con_vars(n) % q)
449 deallocate (con_vars_tl(n) % t)
450 deallocate (con_vars_tl(n) % q)
451 end do
452 deallocate (con_vars)
453 deallocate (aux_vars)
454 deallocate (con_vars_tl)
455 deallocate (aux_vars_tl)
456 end if
457 end do ! end loop for sensor
458
459 if (trace_use) call da_trace_exit("da_transform_xtoy_rttov")
460 #endif
461
462 end subroutine da_transform_xtoy_rttov
463
464