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