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