da_transfer_wrftoxb.inc
References to this file elsewhere.
1 subroutine da_transfer_wrftoxb(xbx, grid)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Transfers fields from WRF to first guess structure.
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 type (xbx_type), intent(inout) :: xbx ! Header & non-gridded vars.
10
11 type(domain), intent(inout) :: grid
12
13 integer :: i, j, k
14
15 real :: theta, tmpvar
16
17 real, dimension(ims:ime,jms:jme) :: rgh_fac
18
19 character(len=19) :: current_date
20
21 real :: loc_psac_mean
22
23 real, dimension(jds:jde) :: loc_latc_mean
24
25 integer :: size2d
26
27 real, dimension(kms:kme) :: DDT
28
29 real :: qvf1, cvpm, cpovcv, ppb, ttb, albn, aln, height
30 real, allocatable :: arrayglobal(:,:)
31
32 !---------------------------------------------------------------------------
33 ! Set xb array range indices for processor subdomain.
34 !---------------------------------------------------------------------------
35
36 if (trace_use) call da_trace_entry("da_transfer_wrftoxb")
37
38 grid%xb % map = grid%map_proj
39 grid%xb % ds = grid%dx
40
41 grid%xb % mix = grid%xp % ide - grid%xp % ids + 1
42 grid%xb % mjy = grid%xp % jde - grid%xp % jds + 1
43 grid%xb % mkz = grid%xp % kde - grid%xp % kds + 1
44
45 ! WHY?
46 !rizvi xbx%big_header%bhi(5,5) = grid%start_year
47 !rizvi xbx%big_header%bhi(6,5) = grid%start_month
48 !rizvi xbx%big_header%bhi(7,5) = grid%start_day
49 !rizvi xbx%big_header%bhi(8,5) = grid%start_hour
50
51 !---------------------------------------------------------------------------
52 ! WRF-specific fitelds:
53 !---------------------------------------------------------------------------
54
55 ptop = grid%p_top
56
57 grid%xb%sigmaf(kte+1) = grid%em_znw(kte+1)
58
59 grid%xb%znw(kte+1) = grid%em_znw(kte+1)
60 grid%xb%znu(kte+1) = 0.0
61
62 do k=kts,kte
63 grid%xb%sigmah(k) = grid%em_znu(k)
64 grid%xb%sigmaf(k) = grid%em_znw(k)
65
66 grid%xb%znu(k) = grid%em_znu(k)
67 grid%xb%znw(k) = grid%em_znw(k)
68 grid%xb%dn(k) = grid%em_dn(k)
69 grid%xb%dnw(k) = grid%em_dnw(k)
70 end do
71
72 grid%xb % ptop = ptop
73
74 !---------------------------------------------------------------------------
75 ! Convert WRF fitelds to xb:
76 !---------------------------------------------------------------------------
77
78 if (print_detail_xb) then
79 write(unit=stdout, fmt='(3a, i8)') &
80 'file:', __FILE__, ', line:', __LINE__
81
82 write(unit=stdout, fmt=*) 'its,ite=', its,ite
83 write(unit=stdout, fmt=*) 'jts,jte=', jts,jte
84 write(unit=stdout, fmt=*) 'kts,kte=', kts,kte
85
86 write(unit=stdout, fmt='(/5a/)') &
87 'lvl dnw dn rdnw rdn'
88
89 do k=kts,kte+1
90 write(unit=stdout, fmt='(i3,8f16.8)') k, &
91 grid%em_dnw(k), grid%em_dn(k), grid%em_rdnw(k), grid%em_rdn(k)
92 end do
93
94 write(unit=stdout, fmt='(/5a/)') &
95 'lvl znu znw rdnw rdn'
96
97 do k=kts,kte+1
98 write(unit=stdout, fmt='(i3,8f16.8)') k, &
99 grid%xb%sigmah(k), grid%xb%sigmaf(k), grid%em_rdnw(k), &
100 grid%em_rdn(k)
101 end do
102
103 write(unit=stdout, fmt='(/5a/)') &
104 'lvl phb ph_2'
105
106 do k=kts,kte
107 write(unit=stdout, fmt='(i3,8e20.12)') k, &
108 grid%em_phb(its,jts,k), grid%em_ph_2(its,jts,k)
109 end do
110
111 write(unit=stdout, fmt=*) 'simple variables:'
112
113 if (jte == jde) then
114 write(unit=stdout, fmt=*) ' '
115
116 do k=kts+5,kte,10
117 do i=its,ite,10
118 write(unit=stdout, fmt=*) &
119 ' grid%em_v_2(', i, ',', jde+1, ',', k, ')=', &
120 grid%em_v_2(i, jde+1,k)
121 end do
122 write(unit=stdout, fmt=*) ' '
123 end do
124 end if
125
126 if (ite == ide) then
127 write(unit=stdout, fmt=*) ' '
128
129 do k=kts+5,kte,10
130 do j=jts,jte,10
131 write(unit=stdout, fmt=*) &
132 ' grid%em_u_2(', ide+1, ',', j, ',', k, ')=', &
133 grid%em_u_2(ide+1,j,k)
134 end do
135 write(unit=stdout, fmt=*) ' '
136 end do
137 end if
138
139 write(unit=stdout, fmt=*) 'simple variables:'
140
141 write(unit=stdout,fmt=*) &
142 ' grid%em_u_1(its,jts,kts)=', grid%em_u_1(its,jts,kts)
143 write(unit=stdout,fmt=*) &
144 ' grid%em_v_1(its,jts,kts)=', grid%em_v_1(its,jts,kts)
145 write(unit=stdout,fmt=*) &
146 ' grid%em_w_1(its,jts,kts)=', grid%em_w_1(its,jts,kts)
147 write(unit=stdout,fmt=*) &
148 ' grid%em_t_1(its,jts,kts)=', grid%em_t_1(its,jts,kts)
149 write(unit=stdout,fmt=*) &
150 ' grid%em_ph_1(its,jts,kts)=',grid%em_ph_1(its,jts,kts)
151
152
153 write(unit=stdout,fmt=*) &
154 ' grid%em_u_2(its,jte,kts)=', grid%em_u_2(its,jte,kts)
155 write(unit=stdout,fmt=*) &
156 ' grid%em_v_2(ite,jts,kts)=', grid%em_v_2(ite,jts,kts)
157 write(unit=stdout,fmt=*) &
158 ' grid%em_w_2(its,jts,kts)=', grid%em_w_2(its,jts,kts)
159 write(unit=stdout,fmt=*) &
160 ' grid%em_t_2(its,jts,kts)=', grid%em_t_2(its,jts,kts)
161 write(unit=stdout,fmt=*) &
162 ' grid%em_ph_2(its,jts,kts)=',grid%em_ph_2(its,jts,kts)
163 write(unit=stdout,fmt=*) &
164 ' grid%em_phb(its,jts,kts)=', grid%em_phb(its,jts,kts)
165
166 write(unit=stdout, fmt=*) &
167 ' grid%sm31,grid%em31,grid%sm32,grid%em32, grid%sm33,grid%em33=', &
168 grid%sm31,grid%em31,grid%sm32,grid%em32,grid%sm33,grid%em33
169
170 write(unit=stdout, fmt=*) ' grid%p_top=', grid%p_top
171 write(unit=stdout, fmt=*) ' grid%em_znu(kts)=', grid%em_znu(kts)
172 write(unit=stdout, fmt=*) ' grid%em_mu0(its,jts)=', grid%em_mu0(its,jts)
173 write(unit=stdout, fmt=*) ' grid%em_mub(its,jts)=', grid%em_mub(its,jts)
174 write(unit=stdout, fmt=*) ' grid%em_mu_2(its,jts)=', &
175 grid%em_mu_2(its,jts)
176
177 write(unit=stdout, fmt=*) ' hbot(its,jts)=', grid%hbot(its,jts)
178 write(unit=stdout, fmt=*) ' htop(its,jts)=', grid%htop(its,jts)
179
180 write(unit=stdout, fmt=*) ' grid%p_top=', grid%p_top
181 write(unit=stdout, fmt=*) ' num_moist=', num_moist
182 write(unit=stdout, fmt=*) ' P_QV=', P_QV
183
184 write(unit=stdout, fmt=*) ' moist(its,jts,kts,2)=', &
185 grid%moist(its,jts,kts,2)
186 write(unit=stdout, fmt=*) ' '
187 end if
188
189 !---------------------------------------------------------------
190 ! Need this to exchange values in the halo region.
191 ! grid%xa%u and grid%xa%v are used as temporary arrays and so
192 ! it is easy to use the existing exchange scheme.
193 !
194 ! Note, this is needed as u_2 and v_2 has no guarantee
195 ! the most east column, and the most north row are
196 ! properly initailized for each tile.
197 !---------------------------------------------------------------
198
199 do j=jts,jte
200 do k=kts,kte
201 do i=its,ite+1
202 grid%xa%u(i,j,k) = grid%em_u_2(i,j,k)
203 end do
204 end do
205 end do
206
207 do j=jts,jte+1
208 do k=kts,kte
209 do i=its,ite
210 grid%xa%v(i,j,k) = grid%em_v_2(i,j,k)
211 end do
212 end do
213 end do
214
215 ! Fill the halo region for u and v.
216
217 #ifdef DM_PARALLEL
218 #include "HALO_PSICHI_UV_ADJ.inc"
219 #endif
220
221
222 if (print_detail_xb) then
223 write(unit=stdout, fmt=*) &
224 ' ids,ide,jds,jde,kds,kde=', ids,ide,jds,jde,kds,kde
225 write(unit=stdout, fmt=*) &
226 ' its,ite,jts,jte,kts,kte=', its,ite,jts,jte,kts,kte
227 write(unit=stdout, fmt=*) &
228 ' ims,ime,jms,jme,kms,kme=', ims,ime,jms,jme,kms,kme
229
230 write(unit=stdout, fmt=*) &
231 ' lbound(grid%xb%u)=', lbound(grid%xb%u)
232 write(unit=stdout, fmt=*) &
233 ' lbound(grid%xb%v)=', lbound(grid%xb%v)
234 write(unit=stdout, fmt=*) &
235 ' lbound(grid%em_u_2)=', lbound(grid%em_u_2)
236 write(unit=stdout, fmt=*) &
237 ' lbound(grid%em_v_2)=', lbound(grid%em_v_2)
238 write(unit=stdout, fmt=*) &
239 ' ubound(grid%xb%u)=', ubound(grid%xb%u)
240 write(unit=stdout, fmt=*) &
241 ' ubound(grid%xb%v)=', ubound(grid%xb%v)
242 write(unit=stdout, fmt=*) &
243 ' ubound(grid%em_u_2)=', ubound(grid%em_u_2)
244 write(unit=stdout, fmt=*) &
245 ' ubound(grid%em_v_2)=', ubound(grid%em_v_2)
246 end if
247
248 do j=jts,jte
249 k = kte+1
250
251 do i=its,ite
252 grid%em_p(i,j,k) = 0.0
253 grid%xb%map_factor(i,j) = grid%msft(i,j)
254 grid%xb%cori(i,j) = grid%f(i,j)
255 grid%xb%tgrn(i,j) = grid%sst(i,j)
256 if (grid%xb%tgrn(i,j) < 100.0) &
257 grid%xb%tgrn(i,j) = grid%tmn(i,j)
258 grid%xb%lat(i,j) = grid%xlat(i,j)
259 grid%xb%lon(i,j) = grid%xlong(i,j)
260 grid%xb%terr(i,j) = grid%ht(i,j)
261 grid%xb%snow(i,j) = grid%snowc(i,j)
262 grid%xb%lanu(i,j) = grid%lu_index(i,j)
263 grid%xb%landmask(i,j) = grid%landmask(i,j)
264 grid%xb%xland(i,j) = grid%xland(i,j)
265 ! Z. Liu below are variables used by RTTOV
266 grid%xb%tsk(i,j) = grid%tsk(i,j)
267 grid%xb%smois(i,j) = grid%smois(i,j,1)
268 grid%xb%tslb(i,j) = grid%tslb(i,j,1)
269 grid%xb%xice(i,j) = grid%xice(i,j)
270 grid%xb%ivgtyp(i,j) = grid%ivgtyp(i,j)
271 grid%xb%isltyp(i,j) = grid%isltyp(i,j)
272 grid%xb%vegfra(i,j) = grid%vegfra(i,j)
273 grid%xb%snowh(i,j) = grid%snowh(i,j)*1000.0 ! meter to mm
274 end do
275
276 ! WHY?
277 ! Adapted the code from "real.init.code" by Y.-R. Guo 05/13/2004:
278
279 ! do i=its,ite
280 ! k = kte
281 ! qvf1 = 0.5*(grid%moist(i,j,k,P_QV)+grid%moist(i,j,k,P_QV))
282 ! qvf2 = 1.0/(1.0+qvf1)
283 ! qvf1 = qvf1*qvf2
284 ! grid%xb%p(i,j,k) = -0.5*(grid%mu_2(i,j)+qvf1* &
285 ! grid%mub(i,j))/grid%em_rdnw(k)/qvf2
286
287 ! do k = kte-1,1,-1
288 ! qvf1 = 0.5*(grid%moist(i,j,k,P_QV)+grid%moist(i,j,k+1,P_QV))
289 ! qvf2 = 1.0/(1.0+qvf1)
290 ! qvf1 = qvf1*qvf2
291 ! grid%p(i,j,k) = grid%p(i,j,k+1) - &
292 ! (grid%mu_2(i,j)+qvf1*grid%mub(i,j))/qvf2/rdn(k+1)
293 ! end do
294 ! end do
295
296 ! Adapted the code from WRF module_big_step_utilitites_em.F ----
297 ! subroutine calc_p_rho_phi Y.-R. Guo (10/20/2004)
298
299 cvpm = - (1.0 - gas_constant/cp)
300 cpovcv = cp / (cp - gas_constant)
301
302 do k=kts,kte
303 do i=its,ite
304 ! The base specific volume (from real.init.code)
305 ppb = grid%em_znu(k) * grid%em_mub(i,j) + ptop
306 ttb = (base_temp + base_lapse*log(ppb/base_pres)) * &
307 (base_pres/ppb)**kappa
308 albn = (gas_constant/base_pres) * ttb * (ppb/base_pres)**cvpm
309
310 qvf1 = 1.0 + grid%moist(i,j,k,P_QV) / rd_over_rv
311 aln = -1.0 / (grid%em_mub(i,j)+grid%em_mu_2(i,j)) * &
312 (albn*grid%em_mu_2(i,j) + grid%em_rdnw(k) * &
313 (grid%em_ph_2(i,j,k+1) - grid%em_ph_2(i,j,k)))
314 ! total pressure:
315 grid%xb%p(i,j,k) = base_pres * &
316 ((gas_constant*(t0+grid%em_t_2(i,j,k))*qvf1) / &
317 (base_pres*(aln+albn)))**cpovcv
318 ! total density
319 grid%xb%rho(i,j,k)= 1.0 / (albn+aln)
320 ! pressure purtubation:
321 grid%em_p(i,j,k) = grid%xb%p(i,j,k) - ppb
322 end do
323 end do
324
325 do k=kts,kte+1
326 do i=its,ite
327 grid%xb%hf(i,j,k) = (grid%em_phb(i,j,k)+grid%em_ph_2(i,j,k))/gravity
328 grid%xa%w (i,j,k) = grid%em_w_2(i,j,k)
329 grid%xb%w (i,j,k) = grid%em_w_2(i,j,k)
330 end do
331 end do
332
333 do k=kts,kte
334 do i=its,ite
335 grid%xb%u(i,j,k) = 0.5*(grid%xa%u(i,j,k)+grid%xa%u(i+1,j,k))
336 grid%xb%v(i,j,k) = 0.5*(grid%xa%v(i,j,k)+grid%xa%v(i,j+1,k))
337 grid%xb%wh(i,j,k)= 0.5*(grid%xb%w(i,j,k)+grid%xb%w(i,j,k+1))
338 grid%xb%h(i,j,k) = 0.5*(grid%xb%hf(i,j,k)+grid%xb%hf(i,j,k+1))
339
340 grid%xb%q(i,j,k) = grid%moist(i,j,k,P_QV)
341
342 theta = t0 + grid%em_t_2(i,j,k)
343 grid%xb%t(i,j,k) = theta*(grid%xb%p(i,j,k)/base_pres)**kappa
344
345 ! Convert to specific humidity from mixing ratio of water vapor:
346 grid%xb%q(i,j,k)=grid%xb%q(i,j,k)/(1.0+grid%xb%q(i,j,k))
347
348 if (grid%xb%q(i,j,k) < 1.0e-6) &
349 grid%xb%q(i,j,k) = 1.0e-6
350
351 ! Background qrn needed for radar radial velocity assmiilation:
352
353 if (size(grid%moist,dim=4) >= 4) then
354 grid%xb%qcw(i,j,k) = grid%moist(i,j,k,p_qc)
355 grid%xb%qrn(i,j,k) = grid%moist(i,j,k,p_qr)
356 grid%xb%qt (i,j,k) = grid%xb%q(i,j,k) + grid%xb%qcw(i,j,k) + &
357 grid%xb%qrn(i,j,k)
358 end if
359
360 if (size(grid%moist,dim=4) >= 6) then
361 grid%xb%qci(i,j,k) = grid%moist(i,j,k,p_qi)
362 grid%xb%qsn(i,j,k) = grid%moist(i,j,k,p_qs)
363 end if
364
365 if (size(grid%moist,dim=4) >= 7) then
366 grid%xb%qgr(i,j,k) = grid%moist(i,j,k,p_qg)
367 end if
368 end do
369 end do
370
371 do i=its,ite
372 grid%xb%psac(i,j) = grid%em_mub(i,j)+grid%em_mu_2(i,j)
373 grid%xb%psfc(i,j) = grid%em_mub(i,j)+grid%em_p(i,j,kts)+grid%p_top
374
375 if (grid%xb%tgrn(i,j) < 100.0) &
376 grid%xb%tgrn(i,j) = grid%xb%t(i,j,kts)+ &
377 0.0065*(grid%xb%h(i,j,kts)-grid%xb%hf(i,j,kts))
378 end do
379 end do
380
381 grid%xb%ztop = grid%xb%hf(its,jts,kte+1)
382
383 if (print_detail_xb) then
384 write(unit=stdout, fmt=*) ' '
385 if (print_detail_xb) then
386 write(unit=stdout, fmt='(/5a/)') &
387 'lvl h p t'
388
389 do k=kts,kte
390 write(unit=stdout, fmt='(i3,8e20.12)') k, &
391 grid%xb%h(its,jts,k), grid%xb%p(its,jts,k), grid%xb%t(its,jts,k)
392 end do
393 end if
394
395 write(unit=stdout,fmt=*) ' '
396 write(unit=stdout,fmt=*) 'grid%xb%u(its,jte,kte)=', grid%xb%u(its,jte,kte)
397 write(unit=stdout,fmt=*) 'grid%xb%v(ite,jts,kte)=', grid%xb%v(ite,jts,kte)
398 write(unit=stdout,fmt=*) 'grid%xb%w(its,jts,kte)=', grid%xb%w(its,jts,kte)
399 write(unit=stdout,fmt=*) 'grid%xb%t(its,jts,kte)=', grid%xb%t(its,jts,kte)
400 write(unit=stdout,fmt=*) 'grid%xb%p(its,jts,kte)=', grid%xb%p(its,jts,kte)
401 write(unit=stdout,fmt=*) 'grid%xb%q(its,jts,kte)=', grid%xb%q(its,jts,kte)
402 write(unit=stdout,fmt=*) 'grid%xb%h(its,jts,kte)=', grid%xb%h(its,jts,kte)
403 write(unit=stdout,fmt=*) &
404 'grid%xb%hf(its,jts,kte)=', grid%xb%hf(its,jts,kte)
405 write(unit=stdout,fmt=*) &
406 'grid%xb%map_factor(its,jts)=', grid%xb%map_factor(its,jts)
407 write(unit=stdout,fmt=*) 'grid%xb%cori(its,jts)=', grid%xb%cori(its,jts)
408 write(unit=stdout,fmt=*) 'grid%xb%tgrn(its,jts)=', grid%xb%tgrn(its,jts)
409 write(unit=stdout,fmt=*) 'grid%xb%lat(its,jts)=', grid%xb%lat(its,jts)
410 write(unit=stdout,fmt=*) 'grid%xb%lon(its,jts)=', grid%xb%lon(its,jts)
411 write(unit=stdout,fmt=*) 'grid%xb%terr(its,jts)=', grid%xb%terr(its,jts)
412 write(unit=stdout,fmt=*) 'grid%xb%snow(its,jts)=', grid%xb%snow(its,jts)
413 write(unit=stdout,fmt=*) 'grid%xb%lanu(its,jts)=', grid%xb%lanu(its,jts)
414 write(unit=stdout,fmt=*) &
415 'grid%xb%landmask(its,jts)=', grid%xb%landmask(its,jts)
416 write(unit=stdout,fmt=*) '(ite,jte)=', ite,jte
417 write(unit=stdout,fmt=*) 'grid%xb%lat(ite,jte)=', grid%xb%lat(ite,jte)
418 write(unit=stdout,fmt=*) 'grid%xb%lon(ite,jte)=', grid%xb%lon(ite,jte)
419 write(unit=stdout,fmt=*) ' '
420 end if
421
422 !---------------------------------------------------------------------------
423 ! [3.0] Calculate vertical inner product for use in vertical transform:
424 !---------------------------------------------------------------------------
425
426 if (vertical_ip == vertical_ip_sqrt_delta_p) then
427 ! Vertical inner product is sqrt(Delta p):
428 do k=kts,kte
429 grid%xb % vertical_inner_product(its:ite,jts:jte,k) = &
430 sqrt(grid%xb % psac(its:ite,jts:jte) * grid%xb%sigmah(k))
431 end do
432 else if (vertical_ip == vertical_ip_delta_p) then
433
434 ! Vertical inner product is Delta p:
435 do k=1,grid%xb%mkz
436 grid % xb % vertical_inner_product(its:ite,jts:jte,k) = &
437 grid % xb % psac(its:ite,jts:jte) * grid%xb%sigmah(k)
438 end do
439 end if
440
441 !---------------------------------------------------------------------------
442 ! Roughness
443 !---------------------------------------------------------------------------
444
445 current_date = 'yyyy-mm-dd_hh:mm:ss'
446
447 write(current_date(1:19), fmt='(i4.4, 5(a1, i2.2))') &
448 grid%start_year, '-', &
449 grid%start_month, '-', &
450 grid%start_day, '_', &
451 grid%start_hour, ':', &
452 grid%start_minute, ':', &
453 grid%start_second
454
455 xbx % mminlu = 'USGS'
456
457 call da_roughness_from_lanu(19, xbx % mminlu, current_date, &
458 grid%xb % lanu, grid%xb % rough)
459
460 !---------------------------------------------------------------------------
461 ! Calculate 1/grid box areas:
462 !---------------------------------------------------------------------------
463
464 if (print_detail_xb) then
465 write(unit=stdout, fmt='(/a, e24.16)') &
466 'grid%xb % ds=', grid%xb % ds
467
468 write(unit=stdout, fmt='(a, e24.16/)') &
469 'grid%xb % map_factor(its,jts)=', grid%xb % map_factor(its,jts)
470 end if
471
472 do j=jts,jte
473 do i=its,ite
474 if (grid%xb%ztop < grid%xb%hf(i,j,kte+1)) &
475 grid%xb%ztop = grid%xb%hf(i,j,kte+1)
476
477 tmpvar = grid%xb%ds / grid%xb%map_factor(i,j)
478
479 grid%xb % grid_box_area(i,j) = tmpvar*tmpvar
480
481 ! Calculate surface variable(wind, moisture, temperature)
482 ! sfc variables: 10-m wind, and 2-m T, Q, at cross points
483
484 height = grid%xb%h(i,j,kts) - grid%xb%terr(i,j)
485 if (height <= 0.0) then
486 message(1) = "Negative height found"
487 write (unit=message(2),FMT='(2I6,A,F10.2,A,F10.2)') &
488 i,j,' ht = ',grid%xb%h(i,j,kts) ,' terr = ',grid%xb%terr(i,j)
489 call da_error(__FILE__,__LINE__, message(1:2))
490 end if
491
492 call da_sfc_wtq(grid%xb%psfc(i,j), grid%xb%tgrn(i,j), &
493 grid%xb%p(i,j,kts), grid%xb%t(i,j,kts), grid%xb%q(i,j,kts), &
494 grid%xb%u(i,j,kts), grid%xb%v(i,j,kts), &
495 grid%xb%p(i,j,kts+1), grid%xb%t(i,j,kts+1), grid%xb%q(i,j,kts+1), &
496 height, grid%xb%rough(i,j),grid%xb%xland(i,j), &
497 grid%xb%u10(i,j), grid%xb%v10(i,j), grid%xb%t2(i,j), &
498 grid%xb%q2(i,j), grid%xb%regime(i,j))
499 end do
500 end do
501
502 !---------------------------------------------------------------------------
503 ! Calculate saturation vapour pressure and relative humidity:
504 !---------------------------------------------------------------------------
505
506 do j=jts,jte
507 do k=kts,kte
508 do i=its,ite
509 call da_tpq_to_rh(grid%xb % t(i,j,k), grid%xb % p(i,j,k), &
510 grid%xb % q(i,j,k), grid%xb %es(i,j,k), grid%xb %qs(i,j,k), &
511 grid%xb %rh(i,j,k))
512 end do
513 end do
514 end do
515
516 !---------------------------------------------------------------------------
517 ! Calculate dew point temperature:
518 !---------------------------------------------------------------------------
519
520 call da_trh_to_td(grid%xb % rh, grid%xb % t, grid%xb % td)
521
522 if (print_detail_xb) then
523 i=its; j=jts; k=kts
524
525 write(unit=stdout, fmt=*) 'i,j,k=', i,j,k
526 write(unit=stdout, fmt=*) 'grid%xb % td(i,j,k)=', grid%xb % td(i,j,k)
527 write(unit=stdout, fmt=*) 'grid%xb % es(i,j,k)=', grid%xb % es(i,j,k)
528 write(unit=stdout, fmt=*) 'grid%xb % rh(i,j,k)=', grid%xb % rh(i,j,k)
529 write(unit=stdout, fmt=*) 'grid%xb % qs(i,j,k)=', grid%xb % qs(i,j,k)
530 write(unit=stdout, fmt=*) ' '
531 end if
532
533 !---------------------------------------------------------------------------
534 ! Sea level pressure and total precipitable water
535 !---------------------------------------------------------------------------
536
537 call da_wrf_tpq_2_slp (grid)
538
539 ! WHY?
540 ! do j = jts,jte
541 ! do i = its,ite
542 ! call da_tpq_to_slp(grid%xb%t(i,j,:), grid%xb%q(i,j,:), &
543 ! grid%xb%p(i,j,:), grid%xb%terr(i,j), &
544 ! grid%xb%psfc(i,j), grid%xb%slp(i,j), grid%xp)
545 ! end do
546 ! end do
547
548 call da_integrat_dz(grid)
549
550 !---------------------------------------------------------------------------
551 ! Surface wind speed
552 !---------------------------------------------------------------------------
553
554 tmpvar = log(10.0/0.0001)
555
556 do j=jts,jte
557 do i=its,ite
558 height = grid%xb%h(i,j,kts) - grid%xb%terr(i,j)
559 rgh_fac(i,j) = 1.0/log(height/0.0001)
560
561 grid%xb%speed(i,j) = sqrt(grid%xb%u(i,j,kts)*grid%xb%u(i,j,kts) &
562 + grid%xb%v(i,j,kts)*grid%xb%v(i,j,kts) + 1.0e-6) &
563 *tmpvar*rgh_fac(i,j)
564 end do
565 end do
566
567 !---------------------------------------------------------------------------
568 ! Brightness temperature SH Chen
569 !---------------------------------------------------------------------------
570
571 if (use_ssmitbobs) &
572 call da_transform_xtotb(grid)
573
574 !---------------------------------------------------------------------------
575 ! GPS Refractivity linked by Y.-R. Guo 05/28/2004
576 !---------------------------------------------------------------------------
577
578 call da_transform_xtogpsref(grid)
579
580 !---------------------------------------------------------------------------
581 ! Ground-based GPS ZTD must follow the GPS Refractivity calculation.
582 !---------------------------------------------------------------------------
583
584 ! WHY?
585 ! if (use_gpsztdobs) call da_transform_xtoztd(grid%xb)
586
587 !---------------------------------------------------------------------------
588 ! Calculate means for later use in setting up background errors.
589 !---------------------------------------------------------------------------
590
591 ! WHY?
592 ! if (.not. associated(xbx % latc_mean)) then
593 allocate (xbx % latc_mean(jds:jde))
594 if (trace_use) call da_trace("da_transfer_wrftoxb",&
595 message="allocated xbx%latc_mean")
596 ! end if
597
598 size2d = (ide-ids+1)*(jde-jds+1)
599
600 tmpvar = 1.0/real(size2d)
601
602 ! Bitwitse-exact reduction preserves operation order of serial code for
603 ! testing, at the cost of much-increased run-time. Turn it off when not
604 ! testing. Thits will always be .false. for a serial or 1-process MPI run.
605 if (test_dm_exact) then
606 allocate(arrayglobal(ids:ide, jds:jde))
607 call da_patch_to_global(grid,grid%xb%psac, arrayglobal)
608 loc_psac_mean = tmpvar*sum(arrayglobal(ids:ide,jds:jde))
609 deallocate(arrayglobal)
610 else
611 loc_psac_mean = tmpvar*sum(grid%xb % psac(its:ite,jts:jte))
612 end if
613
614 tmpvar = 1.0/real(ide-ids+1)
615
616 if (test_dm_exact) then
617 allocate(arrayglobal(ids:ide, jds:jde))
618 call da_patch_to_global(grid,grid%xb%lat, arrayglobal)
619 do j=jds,jde
620 loc_latc_mean(j) = tmpvar*sum(arrayglobal(ids:ide, j))
621 end do
622 deallocate(arrayglobal)
623 else
624 do j=jts,jte
625 loc_latc_mean(j) = tmpvar*sum(grid%xb % lat(its:ite, j))
626 end do
627 end if
628
629 if (test_dm_exact) then
630 ! Broadcast result from monitor to other tasks.
631 call wrf_dm_bcast_real(loc_psac_mean, 1)
632 xbx % psac_mean = loc_psac_mean
633 ! Broadcast result from monitor to other tasks.
634 call wrf_dm_bcast_real(loc_latc_mean, (jde-jds+1))
635 xbx % latc_mean = loc_latc_mean
636 else
637 xbx % psac_mean = wrf_dm_sum_real(loc_psac_mean)
638 call wrf_dm_sum_reals(loc_latc_mean, xbx % latc_mean)
639 end if
640
641 if (print_detail_xb) then
642 ! write(unit=stdout, fmt=*) 'loc_psac_mean =', loc_psac_mean
643 write(unit=stdout, fmt=*) 'xbx % psac_mean=', xbx % psac_mean
644
645 ! write(unit=stdout, fmt=*) 'loc_latc_mean =', loc_latc_mean(jts)
646 write(unit=stdout, fmt=*) 'xbx % latc_mean=', xbx % latc_mean(jts)
647 end if
648
649 ! Calculate time step from one dimensional cloud model parameterization
650
651 if (dt_cloud_model) then
652 do j = jts, jte
653 do i = its, ite
654 call da_cloud_model (grid%xb%t(I,J,:), grid%xb%p(I,J,:), &
655 grid%xb%q(I,J,:), grid%xb%qcw(I,J,:), grid%xb%qrn(I,J,:), &
656 grid%xb%h(I,J,:), grid%xb%hf(I,J,:), ddt, kts, kte)
657
658 do k = kts, kte
659 grid%xb%delt(i,j,k) = DDT(k)
660 end do
661 end do
662 end do
663 end if
664
665 if (trace_use) call da_trace_exit("da_transfer_wrftoxb")
666
667 end subroutine da_transfer_wrftoxb
668