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