da_transfer_kmatoxb.inc
References to this file elsewhere.
1 subroutine da_transfer_kmatoxb(xbx, grid)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Transfers fields from KMA to first guess (xb) 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 integer :: is, ie, js, je, ks, ke
15 real :: tmpvar, earth_radius_sq,conv
16 real :: tmpp,tmp_p,tmpps,tmp_ps
17 ! real :: rgh_fac(grid%xp%ims:grid%xp%ime,grid%xp%jms:grid%xp%jme)
18 character(len=19) :: current_date
19
20 real :: TV(grid%xp%kms:grid%xp%kme)
21 real :: ALPHA(grid%xp%kms:grid%xp%kme)
22 real :: PHALF(grid%xp%kms:grid%xp%kme)
23 real :: PHALFL(grid%xp%kms:grid%xp%kme)
24
25 real :: pu, pd, coef, delp, hydro, rgasg, shgt
26
27 real, dimension(grid%xp%jds:grid%xp%jde) :: loc_latc_mean
28
29 real, allocatable :: arrayglobal(:,:)
30
31 if (trace_use) call da_trace_entry("da_transfer_kmatoxb")
32
33 conv = pi/180.0
34 earth_radius_sq = earth_radius*1000.0 * earth_radius*1000.0 * &
35 conv*(360.0/(coarse_ix-1))*(180.0/(coarse_jy-2))*conv
36 COEF=0.6080
37 RGASG = gas_constant/gravity
38
39 !---------------------------------------------------------------------------
40 ! Set xb array range indices for processor subdomain.
41 !---------------------------------------------------------------------------
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 % ds = grid%dx
51 grid%xb % kma_a = grid%kma_a
52 grid%xb % kma_b = grid%kma_b
53
54 if (print_detail_xb) then
55 write(unit=stdout, fmt='(/a/)') &
56 'lvl kma_a kma_b'
57
58 do k=ks,ke
59 write(unit=stdout, fmt='(i3,8e20.12)') k, grid%xb%kma_a(k), grid%xb%kma_b(k)
60 end do
61 end if
62
63 grid%xb % mix = grid%xp % ide - grid%xp % ids + 1
64 grid%xb % mjy = grid%xp % jde - grid%xp % jds + 1
65 grid%xb % mkz = grid%xp % kde - grid%xp % kds + 1
66
67 !---------------------------------------------------------------------------
68 ! KMA-specific fields:
69 !---------------------------------------------------------------------------
70
71 ! Fix ptop as 0.4 hPa.
72 ptop = 40.0
73 ps0 = 100000.0
74 tis0 = 0.0
75 tlp = 0.0
76
77 grid%xb % ptop = ptop
78 grid%xb % ps0 = ps0
79 grid%xb % tlp = tlp
80
81 !---------------------------------------------------------------------------
82 ! Convert KMA fields to xb:
83 !---------------------------------------------------------------------------
84
85 grid%xb%lat(is:ie,js:je) = grid%xlat(is:ie,js:je)
86 grid%xb%lon(is:ie,js:je) = grid%xlong(is:ie,js:je)
87
88 ! Fill halo region for u and v.
89 call wrf_dm_halo(grid%xp%domdesc,grid%xp%comms,grid%xp%halo_id3)
90 tmpvar = 1.0/real(grid%xp%ide-grid%xp%ids+1)
91
92 !---------------------------------------------------------------
93 ! transfer u,v,t & q(specific humidity in g/g)
94 !---------------------------------------------------------------
95
96 do k=ks,ke
97 do j=js,je
98 do i=is,ie
99 grid%xb%u(i,j,k) = grid%em_u_2(i,j,k)
100 grid%xb%v(i,j,k) = grid%em_v_2(i,j,k)
101 grid%xb%t(i,j,k) = grid%em_t_2(i,j,k)
102 grid%xb%q(i,j,k) = grid%moist(i,j,k,P_QV)
103 end do
104 end do
105 end do
106
107 !---------------------------------------------------------------------------
108 ! Fix wind at Poles
109 !---------------------------------------------------------------------------
110
111 call da_get_vpoles(grid%xb%u,grid%xb%v, &
112 grid%xp%ids,grid%xp%ide,grid%xp%jds,grid%xp%jde, &
113 grid%xp%ims,grid%xp%ime,grid%xp%jms,grid%xp%jme,grid%xp%kms,grid%xp%kme, &
114 grid%xp%its,grid%xp%ite,grid%xp%jts,grid%xp%jte,grid%xp%kts,grid%xp%kte )
115
116 if (print_detail_xb) then
117 write(unit=stdout,fmt='(A,2I5)') 'grid%xp%its, grid%xp%ite=', grid%xp%its, grid%xp%ite
118 write(unit=stdout,fmt='(A,2I5)') 'grid%xp%jts, grid%xp%jte=', grid%xp%jts, grid%xp%jte
119 write(unit=stdout,fmt='(A,2I5)') 'grid%xp%kts, grid%xp%kte=', grid%xp%kts, grid%xp%kte
120 write(unit=stdout,fmt='(A,2I5)') 'grid%xp%ims, grid%xp%ime=', grid%xp%ims, grid%xp%ime
121 write(unit=stdout,fmt='(A,2I5)') 'grid%xp%jms, grid%xp%jme=', grid%xp%jms, grid%xp%jme
122 write(unit=stdout,fmt='(A,2I5)') 'grid%xp%kms, grid%xp%kme=', grid%xp%kms, grid%xp%kme
123 write(unit=stdout,fmt='(A,2I5)') 'grid%xp%ids, grid%xp%ide=', grid%xp%ids, grid%xp%ide
124 write(unit=stdout,fmt='(A,2I5)') 'grid%xp%jds, grid%xp%jde=', grid%xp%jds, grid%xp%jde
125 write(unit=stdout,fmt='(A,2I5)') 'grid%xp%kds, grid%xp%kde=', grid%xp%kds, grid%xp%kde
126
127 write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%u, dim=1)=', size(grid%xb%u, dim=1)
128 write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%u, dim=2)=', size(grid%xb%u, dim=2)
129 write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%u, dim=3)=', size(grid%xb%u, dim=3)
130 write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%v, dim=1)=', size(grid%xb%v, dim=1)
131 write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%v, dim=2)=', size(grid%xb%v, dim=2)
132 write(unit=stdout,fmt='(A,I5)') 'size(grid%xb%v, dim=3)=', size(grid%xb%v, dim=3)
133
134 write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%u, dim=1)=', size(grid%xa%u, dim=1)
135 write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%u, dim=2)=', size(grid%xa%u, dim=2)
136 write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%u, dim=3)=', size(grid%xa%u, dim=3)
137 write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%v, dim=1)=', size(grid%xa%v, dim=1)
138 write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%v, dim=2)=', size(grid%xa%v, dim=2)
139 write(unit=stdout,fmt='(A,I5)') 'size(grid%xa%v, dim=3)=', size(grid%xa%v, dim=3)
140
141 write(unit=stdout,fmt='(A,I5)') 'size(grid%em_u_2, dim=1)=', size(grid%em_u_2, dim=1)
142 write(unit=stdout,fmt='(A,I5)') 'size(grid%em_u_2, dim=2)=', size(grid%em_u_2, dim=2)
143 write(unit=stdout,fmt='(A,I5)') 'size(grid%em_u_2, dim=3)=', size(grid%em_u_2, dim=3)
144 write(unit=stdout,fmt='(A,I5)') 'size(grid%em_v_2, dim=1)=', size(grid%em_v_2, dim=1)
145 write(unit=stdout,fmt='(A,I5)') 'size(grid%em_v_2, dim=2)=', size(grid%em_v_2, dim=2)
146 write(unit=stdout,fmt='(A,I5)') 'size(grid%em_v_2, dim=3)=', size(grid%em_v_2, dim=3)
147 end if
148
149 ALPHA(ke) = LOG(2.0)
150
151 do j=js,je
152 do i=is,ie
153 grid%xb%cori(i,j) = grid%f(i,j)
154 grid%xb%psfc(i,j) = grid%psfc(i,j)
155 grid%xb%terr(i,j) = grid%ht(i,j)
156
157 ! Zhiquan Liu add some RTTOV variables
158 !--------------------------------------
159 grid%xb%t2(i,j) = grid%t2(i,j)
160 grid%xb%q2(i,j) = grid%q2(i,j)
161 grid%xb%u10(i,j) = grid%u10(i,j)
162 grid%xb%v10(i,j) = grid%v10(i,j)
163 grid%xb%tsk(i,j) = grid%tsk(i,j)
164 ! currently KMA guess have no SST, replaced by TSK
165 ! grid%xb%tgrn(i,j) = grid%sst(i,j)
166 grid%xb%tgrn(i,j) = grid%tsk(i,j)
167 grid%xb%landmask(i,j) = grid%landmask(i,j)
168 grid%xb%xland(i,j) = grid%xland(i,j)
169
170 grid%xb%smois(i,j) = grid%smois(i,j,1)
171 grid%xb%tslb(i,j) = grid%tslb(i,j,1)
172 grid%xb%xice(i,j) = grid%xice(i,j)
173 grid%xb%ivgtyp(i,j) = grid%ivgtyp(i,j)
174 grid%xb%isltyp(i,j) = grid%isltyp(i,j)
175 grid%xb%vegfra(i,j) = grid%vegfra(i,j)
176 ! convert snow water equivalent kg/m2 to snow depth(mm)
177 grid%xb%snowh(i,j) = grid%snow(i,j)*10.
178
179 !---------------------------------------------------------------------
180 ! Syed RH Rizvi
181 !
182 ! The following code is to be activated after
183 ! getting sst, land use, land mask and snow information for KMA grid
184 ! This infor needs to be packed in "KMA2NETCDF" convertor
185 !
186 !---------------------------------------------------------------------
187
188 ! grid%xb%tgrn(i,j) = grid%sst(i,j)
189 ! if (grid%xb%tgrn(i,j) < 100.0) grid%xb%tgrn(i,j) = tmn(i,j)
190 ! grid%xb%lanu(i,j) = grid%lu_index(i,j)
191 ! grid%xb%landmask(i,j) = grid%landmask(i,j)
192 ! grid%xb%xland(i,j) = grid%xland(i,j)
193 ! grid%xb%snow(i,j) = grid%snowc(i,j)
194
195 ! Since data is on full levels get full level pr & ht.
196 ! Note p = A + B * Psfc formula needs Psfc to be in hPa
197
198 do K = ks,ke-1
199 PU = grid%xb%KMA_A(K+1) + grid%xb%KMA_B(K+1)*grid%xb%psfc(i,j)/100.0
200 PD = grid%xb%KMA_A(K ) + grid%xb%KMA_B(K )*grid%xb%psfc(i,j)/100.0
201 if (PU == PD)then
202 message(1)='PU, PD equal and so denominator will be zero'
203 write(unit=message(2),fmt='(A,3I5)') ' i, j, k = ',i,j,k
204 write(unit=message(3),fmt='(A,2F10.2)') &
205 ' KMA_A ',grid%KMA_A(k),grid%KMA_A(K+1)
206 write(unit=message(4),fmt='(A,2F10.2)') &
207 ' KMA_B ',grid%KMA_B(k),grid%KMA_B(K+1)
208 write(unit=message(5),fmt='(A,2F10.2)') &
209 ' psfc(Pa) = ',grid%xb%psfc(i,j)
210 call da_error(__FILE__,__LINE__,message(1:5))
211 end if
212 grid%xb%p(i,j,k)= 100.0 * exp ((PD*LOG(PD)-PU*LOG(PU))/(PD-PU) -1.0)
213 end do
214
215 grid%xb%p(i,j,ke) = 100.0*(grid%xb%KMA_A(ke)+grid%xb%KMA_B(ke)*grid%xb%psfc(i,j)/100.0)/2.0
216
217 do k=ks,ke
218 if (grid%xb%t(i,j,k) <= 0.0) then
219 write(unit=message(1),fmt='(A,3I5,F10.2)') &
220 'Problem in temp = ',i,j,k,grid%xb%t(i,j,k)
221 call da_error(__FILE__,__LINE__,message(1:1))
222 end if
223
224 grid%xb%rho(i,j,k)=grid%xb%p(i,j,k)/(gas_constant* &
225 (1.0+COEF*grid%xb%q(I,J,K))*grid%xb%t(I,J,K))
226 end do
227
228 ! compute full level height
229
230 do K = ks,ke
231 PHALF(K) = grid%xb%KMA_A(K) + grid%xb%KMA_B(K)*grid%xb%psfc(I,J)/100.0
232 TV (K) = (1.0+COEF*grid%xb%q(I,J,K))*grid%xb%t(I,J,K)*rgasg
233 end do
234
235 do K = ks,ke-1
236 DELP = PHALF(K) - PHALF(K+1)
237 PHALFL(K)= LOG(PHALF(K)/PHALF(K+1))
238 ALPHA(K) = 1.0-PHALF(K+1)*PHALFL(K)/DELP
239 end do
240
241 SHGT = grid%xb%terr(i,j)
242 do K = ks, ke
243 grid%xb%h(I,J,K) = SHGT + ALPHA(K)*TV(K)
244 end do
245
246 HYDRO = 0.0
247 do K = ks+1, ke
248 HYDRO = HYDRO + TV(K-1)*PHALFL(K-1)
249 grid%xb%h(I,J,K) = grid%xb%h(I,J,K) + HYDRO
250 end do
251 end do
252 end do
253
254 !---------------------------------------------------------------------------
255 ! Sigma values are needed for interpolating
256 ! background error statistics in vertical
257 ! Fix sigmah values based on global average surface pressure
258 ! and level wise pressure
259 !---------------------------------------------------------------------------
260
261 tmp_ps = sum(grid%xb%psfc(is:ie,js:je)) /real((grid%xp%ide-grid%xp%ids+1)*(grid%xp%jde-grid%xp%jds+1))
262
263 tmpps = wrf_dm_sum_real(tmp_ps)
264
265 do k=ks,ke
266 tmp_p = sum(grid%xb%p(is:ie,js:je,k)) /real((grid%xp%ide-grid%xp%ids+1)*(grid%xp%jde-grid%xp%jds+1))
267 tmpp = wrf_dm_sum_real(tmp_p)
268
269 ! 0.01 is added to see that sigmah should not become negative
270 grid%xb%sigmah(ke+ks-k) = (tmpp - ptop+0.01) / (tmpps- ptop+0.01)
271 if (grid%xb%sigmah(ke+ks-k) < 0.0) then
272 write(unit=message(1),fmt='(A,F10.2)') &
273 ' average surface pressure = ',tmpps
274 write(unit=message(2),fmt='(A,I5,A,F10.2)') &
275 ' average pressure at this level= ',k,' = ',tmpp
276 write(unit=message(3),fmt='(A,I5,A,F10.2)') &
277 ' negative sigmah(',ke+ks-k,') = ',grid%xb%sigmah(ke+ks-k)
278 call da_error(__FILE__,__LINE__,message(1:3))
279 end if
280 end do
281
282 ! Fix ztop
283
284 grid%xb%ztop = grid%xb%h(is,js,ke)
285
286 if (print_detail_xb) then
287 write(unit=stdout, fmt='(/5a/)') &
288 'lvl h p t'
289
290 do k=ks,ke
291 write(unit=stdout, fmt='(i3,8e20.12)') k, &
292 grid%xb%h(is,js,k), grid%xb%p(is,js,k), grid%xb%t(is,js,k)
293 end do
294
295 write (unit=stdout,fmt='(A)') ' '
296 write (unit=stdout,fmt='(A,I5)') 'grid%xb%u(is,je,ke)=', grid%xb%u(is,je,ke)
297 write (unit=stdout,fmt='(A,I5)') 'grid%xb%v(ie,js,ke)=', grid%xb%v(ie,js,ke)
298 write (unit=stdout,fmt='(A,I5)') 'grid%xb%w(is,js,ke)=', grid%xb%w(is,js,ke)
299 write (unit=stdout,fmt='(A,I5)') 'grid%xb%t(is,js,ke)=', grid%xb%t(is,js,ke)
300 write (unit=stdout,fmt='(A,I5)') 'grid%xb%p(is,js,ke)=', grid%xb%p(is,js,ke)
301 write (unit=stdout,fmt='(A,I5)') 'grid%xb%q(is,js,ke)=', grid%xb%q(is,js,ke)
302 write (unit=stdout,fmt='(A,I5)') 'grid%xb%hf(is,js,ke)=', grid%xb%hf(is,js,ke)
303 write (unit=stdout,fmt='(A,I5)') 'grid%xb%map_factor(is,js)=', grid%xb%map_factor(is,js)
304 write (unit=stdout,fmt='(A,I5)') 'grid%xb%cori(is,js)=', grid%xb%cori(is,js)
305 write (unit=stdout,fmt='(A,I5)') 'grid%xb%tgrn(is,js)=', grid%xb%tgrn(is,js)
306 write (unit=stdout,fmt='(A,I5)') 'grid%xb%lat(is,js)=', grid%xb%lat(is,js)
307 write (unit=stdout,fmt='(A,I5)') 'grid%xb%lon(is,js)=', grid%xb%lon(is,js)
308 write (unit=stdout,fmt='(A,I5)') 'grid%xb%terr(is,js)=', grid%xb%terr(is,js)
309 write (unit=stdout,fmt='(A,I5)') 'grid%xb%snow(is,js)=', grid%xb%snow(is,js)
310 write (unit=stdout,fmt='(A,I5)') 'grid%xb%lanu(is,js)=', grid%xb%lanu(is,js)
311 write (unit=stdout,fmt='(A,I5)') 'grid%xb%landmask(is,js)=', grid%xb%landmask(is,js)
312 write (unit=stdout,fmt='(A)') ' '
313 end if
314
315 write(current_date(1:19), fmt='(i4.4, 5(a1, i2.2))') &
316 grid%start_year, '-', &
317 grid%start_month, '-', &
318 grid%start_day, '_', &
319 grid%start_hour, ':', &
320 grid%start_minute,':', &
321 grid%start_second
322
323 write(unit=stdout,fmt='(2A)') 'Current date is ',current_date
324
325 !---------------------------------------------------------------------------
326 ! Syed RH Rizvi
327 !
328 ! Following code for calculating roughness length needs to be activated
329 ! after getting land use data for KMA grid
330 !---------------------------------------------------------------------------
331
332 !
333 ! Set proper value for "mminlu" if using KMA info
334
335 !xbx % mminlu = 'USGS'
336 !call da_roughness_from_lanu(19, xbx % mminlu, current_date, grid%xp, &
337 ! grid%xb % lanu(grid%xp%ims,grid%xp%jms), grid%xb % rough(grid%xp%ims,grid%xp%jms))
338
339 do j=js,je
340 do i=is,ie
341 if (grid%xb%ztop < grid%xb%hf(i,j,ke+1)) grid%xb%ztop = grid%xb%hf(i,j,ke+1)
342 ! Calculate grid_box area and vertical inner product for use in
343 ! vertical transform:
344 grid%xb % grid_box_area(i,j) = earth_radius_sq * cos(grid%xlat(i,j)*conv)
345
346 if (vertical_ip == 1) then
347 ! Vertical inner product is sqrt(Delta p):
348 do k=1,grid%xb%mkz
349 grid%xb % vertical_inner_product(i,j,k) = &
350 sqrt(grid%xb%p(i,j,k)-grid%xb%p(i,j,k+1))
351 end do
352 else if (vertical_ip == 2) then
353 ! Vertical inner product is Delta p:
354 do k=1,grid%xb%mkz
355 grid%xb % vertical_inner_product(i,j,k) = grid%xb%p(i,j,k)-grid%xb%p(i,j,k+1)
356 end do
357 end if
358
359 !------------------------------------------------------------------------------
360 ! Syed RH Rizvi
361 !
362 ! Following code for calculating roughness length needs to be activated
363 ! to calculate surface variables (10-m wind, and 2-m T, Q) ,
364 ! After testing KMA-WRFVAR for SFC_ASSIM_OPTIONS = 2
365 !
366 !------------------------------------------------------------------------------
367
368 ! call da_sfc_wtq(grid%xb%psfc(i,j), grid%xb%tgrn(i,j), &
369 ! grid%xb%p(i,j,ks), grid%xb%t(i,j,ks), grid%xb%q(i,j,ks), &
370 ! grid%xb%u(i,j,ks), grid%xb%v(i,j,ks), &
371 ! grid%xb%p(i,j,ks+1), grid%xb%t(i,j,ks+1), grid%xb%q(i,j,ks+1), &
372 ! grid%xb%h(i,j,ks), grid%xb%rough(i,j),grid%xb%landmask(i,j), &
373 ! grid%xb%u10(i,j), grid%xb%v10(i,j), grid%xb%t2(i,j), grid%xb%q2(i,j), &
374 ! grid%xb%regime(i,j))
375
376 ! write (unit=stdout,fmt='(7I5)') &
377 ! i,j,grid%xb%psfc(i,j),grid%xb%t2(i,j),grid%xb%q2(i,j), &
378 ! grid%xb%u10(i,j),grid%xb%v10(i,j)
379 ! ------------------------------------------------------------------------------
380
381 end do
382 end do
383
384 !---------------------------------------------------------------------------
385 ! Calculate saturation vapour pressure and relative humidity:
386 !---------------------------------------------------------------------------
387
388 do j=js,je
389 do k=ks,ke
390 do i=is,ie
391 call da_tpq_to_rh(grid%xb % t(i,j,k), grid%xb % p(i,j,k), &
392 grid%xb % q(i,j,k), &
393 grid%xb %es(i,j,k), grid%xb %qs(i,j,k), grid%xb %rh(i,j,k))
394 end do
395 end do
396 end do
397
398 !---------------------------------------------------------------------------
399 ! Calculate dew point temperature:
400 !---------------------------------------------------------------------------
401
402 call da_trh_to_td(grid%xb % rh, grid%xb % t, grid%xb % td, grid%xp)
403
404 if (print_detail_xb) then
405 i=is; j=js; k=ks
406
407 write (unit=stdout,fmt='(A,3I5)') 'i,j,k=', i,j,k
408 write (unit=stdout,fmt='(A,F10.2)') 'grid%xb % td(i,j,k)=', grid%xb % td(i,j,k)
409 write (unit=stdout,fmt='(A,F10.2)') 'grid%xb % es(i,j,k)=', grid%xb % es(i,j,k)
410 write (unit=stdout,fmt='(A,F10.2)') 'grid%xb % rh(i,j,k)=', grid%xb % rh(i,j,k)
411 write (unit=stdout,fmt='(A,F10.2)') 'grid%xb % qs(i,j,k)=', grid%xb % qs(i,j,k)
412 write (unit=stdout,fmt='(A)') ' '
413 end if
414
415 !---------------------------------------------------------------------------
416 ! Sea level pressure and total precipitable water
417 !---------------------------------------------------------------------------
418
419 if (print_detail_xb) then
420 write(unit=stdout, fmt='(3a, i8)') &
421 'file:', __FILE__, ', line:', __LINE__
422
423 write(unit=stdout, fmt='(a)') &
424 'Sea level pressure and total precipitable water'
425 end if
426
427 !---------------------------------------------------------------------------
428 ! Syed RH Rizvi
429 !
430 ! Following code for calculating roughness length needs to be activated
431 ! for calculating roughness length if sea level pressure is desired
432 !---------------------------------------------------------------------------
433
434 ! call da_wrf_tpq_2_slp (grid%xb)
435
436 call da_integrat_dz(grid%xb)
437 !---------------------------------------------------------------------------
438 ! Syed RH Rizvi
439 !
440 ! Following code for calculating roughness length needs to be activated
441 ! for surface wind speed for SFC_ASSIM_OPTIONS = 2
442 !
443 ! It is not working because KMA terrain height output from wave to grid
444 ! transform is negative at many grid points
445 !---------------------------------------------------------------------------
446
447 ! if (print_detail_xb) then
448 ! write(unit=stdout, fmt='(3a, i8)') &
449 ! 'file:', __FILE__, ', line:', __LINE__
450 ! write(unit=stdout, fmt='(a)') &
451 ! 'Surface Wind speed'
452 ! end if
453
454 ! tmpvar = log(10.0/0.0001)
455 ! do j=js,je
456 ! do i=is,ie
457 ! if (grid%xb%hf(i,j,ks) <= 0) then
458 ! write(unit=message(1), fmt=*) ' zero hf at i/j ',i,j,' ks= ',ks
459 ! call da_error(__FILE__,__LINE__,message(1:1))
460 ! end if
461 ! rgh_fac(i,j) = 1.0/log(grid%xb%hf(i,j,ks)/0.0001)
462
463 ! grid%xb%speed(i,j) = sqrt(grid%xb%u(i,j,ks)*grid%xb%u(i,j,ks) &
464 ! + grid%xb%v(i,j,ks)*grid%xb%v(i,j,ks) + 1.0e-6) &
465 ! *tmpvar*rgh_fac(i,j)
466 ! end do
467 ! end do
468
469 !---------------------------------------------------------------------------
470 ! Syed RH Rizvi
471 !
472 ! Following code for calculating roughness length needs to be activated
473 ! if SSMI brightness temperature are used
474 !---------------------------------------------------------------------------
475
476 ! if (print_detail_xb) then
477 ! write(unit=stdout, fmt='(3a, i8)') &
478 ! 'file:', __FILE__, ', line:', __LINE__
479
480 ! write(unit=stdout, fmt='(a)') &
481 ! 'Brightness temperature'
482 ! end if
483
484 ! call da_transform_xtotb(xa)
485
486 !---------------------------------------------------------------------------
487 ! Calculate means for later use in setting up background errors.
488 !---------------------------------------------------------------------------
489
490 allocate (xbx % latc_mean(grid%xp%jds:grid%xp%jde))
491
492 tmpvar = 1.0/real(grid%xp%ide-grid%xp%ids+1)
493 loc_latc_mean(:) = 0.0
494
495 ! Bitwise-exact reduction preserves operation order of serial code for
496 ! testing, at the cost of much-increased run-time. Turn it off when not
497 ! testing. This will always be .false. for a serial or 1-process MPI run.
498 if (testing_dm_exact) then
499 allocate(arrayglobal(grid%xp%ids:grid%xp%ide, grid%xp%jds:grid%xp%jde))
500 call da_wv_patch_to_global(grid%xb%lat, arrayglobal, grid%xp%domdesc, &
501 grid%xp%ids, grid%xp%ide, grid%xp%jds, grid%xp%jde, &
502 grid%xp%ims, grid%xp%ime, grid%xp%jms, grid%xp%jme, &
503 grid%xp%its, grid%xp%ite, grid%xp%jts, grid%xp%jte)
504 do j=grid%xp%jds,grid%xp%jde
505 loc_latc_mean(j) = tmpvar*sum(arrayglobal(grid%xp%ids:grid%xp%ide, j))
506 end do
507 deallocate(arrayglobal)
508 ! Broadcast result from monitor to other tasks.
509 call wrf_dm_bcast_real(loc_latc_mean, (grid%xp%jde-grid%xp%jds+1))
510 xbx % latc_mean = loc_latc_mean
511 else
512 do j=js,je
513 loc_latc_mean(j) = tmpvar*sum(grid%xb % lat(is:ie, j))
514 end do
515 call wrf_dm_sum_reals (loc_latc_mean, xbx % latc_mean)
516 end if
517
518 write(unit=stdout,fmt=*) 'js je is ie= ', js,je,is,ie
519
520 ! do j=js,je
521 ! do i=is,ie
522 ! write(unit=stdout,fmt='(2i4,3f12.2,e12.4,2f12.2)') &
523 ! j,i,grid%xb%psfc(i,j),grid%xb%tsk(i,j),grid%xb%t2(i,j), &
524 ! grid%xb%q2(i,j),grid%xb%u10(i,j),grid%xb%v10(i,j)
525 ! end do
526 ! end do
527
528 if (trace_use) call da_trace_exit("da_transfer_kmatoxb")
529
530 end subroutine da_transfer_kmatoxb
531
532