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