da_transfer_xatowrf.inc
References to this file elsewhere.
1 subroutine da_transfer_xatowrf(grid)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Convert analysis increments into WRF increments
5 !---------------------------------------------------------------------------
6
7 implicit none
8
9 type(domain), intent(inout) :: grid
10
11 integer :: i, j, k
12
13 integer :: is, ie, js, je, ks, ke
14
15 real :: sdmd, s1md
16
17 ! arrays to hold wrf increments on the c-grid
18
19 real, dimension(grid%xp%ims:grid%xp%ime,grid%xp%jms:grid%xp%jme, &
20 grid%xp%kms:grid%xp%kme) :: &
21 u_cgrid, v_cgrid, q_cgrid, ph_cgrid
22
23 real, dimension(grid%xp%ims:grid%xp%ime,grid%xp%jms:grid%xp%jme) :: mu_cgrid
24
25 real :: t_full , p_full, rho_full, q_full, ph_full , ph_xb_hd, &
26 qvf1, qvf2, qvf1_b, qvf2_b
27
28 if (trace_use) call da_trace_entry("da_transfer_xatowrf")
29
30 !------------------------------------------------------------------------
31 ! set array range indices for processor subdomain.et array range indices for processor subdomain.
32 !------------------------------------------------------------------------
33
34 is = grid%xp % its
35 ie = grid%xp % ite
36 js = grid%xp % jts
37 je = grid%xp % jte
38 ks = grid%xp % kts
39 ke = grid%xp % kte
40
41 ! To keep the background PH perturbation:
42
43 do j=js,je
44 do i=is,ie
45 do k=ks, ke+1
46 ph_cgrid(i,j,k) = grid%em_ph_2(i,j,k)
47 end do
48 end do
49 end do
50
51 !---------------------------------------------------------------------------
52 ! [1.0] Get the mixing ratio of moisture first, as it is easy.
53 !---------------------------------------------------------------------------
54
55 do k=ks,ke
56 do j=js,je
57 do i=is,ie
58 if ((grid%xb%q(i,j,k)+grid%xa%q(i,j,k)) < 0.0) then
59 q_cgrid(i,j,k) =-grid%xb%q(i,j,k)/(1.0 - grid%xb%q(i,j,k))**2
60 else
61 q_cgrid(i,j,k) = grid%xa%q(i,j,k)/(1.0 - grid%xb%q(i,j,k))**2
62 end if
63 end do
64 end do
65 end do
66
67 !---------------------------------------------------------------------------
68 ! [2.0] compute increments of dry-column air mass per unit area
69 !---------------------------------------------------------------------------
70
71 do j=js,je
72 do i=is,ie
73 sdmd=0.0
74 s1md=0.0
75 do k=ks,ke
76 sdmd=sdmd+q_cgrid(i,j,k)*grid%em_dnw(k)
77 s1md=s1md+(1.0+grid%moist(i,j,k,P_QV))*grid%em_dnw(k)
78 end do
79
80 mu_cgrid(i,j)=-(grid%xa%psfc(i,j)+grid%xb%psac(i,j)*sdmd)/s1md
81 end do
82 end do
83
84 !---------------------------------------------------------------------------
85 ! [3.0] compute pressure increments
86 !---------------------------------------------------------------------------
87
88 ! Tangent linear code for grid%xa%p (based on WRF "real.init.code")
89 ! developed by Y.-R. Guo 05/13/2004:
90
91 do j=js,je
92 do i=is,ie
93 k = ke
94 qvf1 = 0.5*(q_cgrid(i,j,k)+q_cgrid(i,j,k))
95 qvf1_b = 0.5*(grid%moist(i,j,k,P_QV)+grid%moist(i,j,k,P_QV))
96 qvf2 = - qvf1 / ((1.+qvf1_b)*(1.+qvf1_b))
97 qvf2_b = 1./(1.+qvf1)
98 qvf1 = qvf1*qvf2_b + qvf1_b*qvf2
99 qvf1_b = qvf1_b*qvf2_b
100 grid%xa%p(i,j,k) = (-0.5/grid%em_rdnw(k)) * &
101 ((mu_cgrid(i,j)+qvf1*grid%em_mub(i,j)) / qvf2_b &
102 -(grid%em_mu_2(i,j)+qvf1_b*grid%em_mub(i,j))*qvf2/(qvf2_b*qvf2_b))
103
104 do k = ke-1,1,-1
105 qvf1 = 0.5*(q_cgrid(i,j,k)+q_cgrid(i,j,k+1))
106 qvf1_b = 0.5*(grid%moist(i,j,k,P_QV)+grid%moist(i,j,k+1,P_QV))
107 qvf2 = - qvf1 / ((1.+qvf1_b)*(1.+qvf1_b))
108 qvf2_b = 1./(1.+qvf1_b)
109 qvf1 = qvf1*qvf2_b + qvf1_b*qvf2
110 qvf1_b = qvf1_b*qvf2_b
111 grid%xa%p(i,j,k) = grid%xa%p(i,j,k+1) &
112 - (1./grid%em_rdn(k+1)) * &
113 ((mu_cgrid(i,j)+qvf1*grid%em_mub(i,j)) / qvf2_b &
114 -(grid%em_mu_2(i,j)+qvf1_b*grid%em_mub(i,j))*qvf2/(qvf2_b*qvf2_b))
115 end do
116 end do
117 end do
118
119 ! Adjust grid%xa%q to make grid%xb%q + grid%xa%q > 0.0
120
121 if (check_rh == 2) then
122 ! Shu-Hua~s TPW conservation:
123 call da_check_rh(grid%xb, grid%xa, grid%xp)
124 else if (check_rh == 1) then
125 ! Simple resetting to max/min values:
126 call da_check_rh_simple(grid%xb, grid%xa, grid%xp)
127 end if
128
129 do k=ks,ke
130 do j=js,je
131 do i=is,ie
132 q_cgrid(i,j,k) = grid%xa%q(i,j,k)/(1.0 - grid%xb%q(i,j,k))**2
133 end do
134 end do
135 end do
136
137 !---------------------------------------------------------------------------
138 ! [4.0] Convert temperature increments into theta increments
139 ! Evaluate also the increments of (1/rho) and geopotential
140 !---------------------------------------------------------------------------
141
142 if (print_detail_xa) then
143 write(unit=stdout, fmt='(a, e24.12)') &
144 'sum(abs(grid%xa%t(is:ie,js:je,ks:ke)))=', &
145 sum(abs(grid%xa%t(is:ie,js:je,ks:ke))), &
146 'sum(abs(grid%xa%p(is:ie,js:je,ks:ke)))=', &
147 sum(abs(grid%xa%p(is:ie,js:je,ks:ke))), &
148 'sum(abs(grid%xb%t(is:ie,js:je,ks:ke)))=', &
149 sum(abs(grid%xb%t(is:ie,js:je,ks:ke))), &
150 'sum(abs(grid%xb%p(is:ie,js:je,ks:ke)))=', &
151 sum(abs(grid%xb%p(is:ie,js:je,ks:ke))), &
152 'sum(abs(grid%em_t_2 (is:ie,js:je,ks:ke)))=', &
153 sum(abs(grid%em_t_2 (is:ie,js:je,ks:ke)))
154
155 write(unit=stdout, fmt='(2(2x, a, e20.12))') &
156 'maxval(grid%xa%u(is:ie,js:je,ks:ke))=', &
157 maxval(grid%xa%u(is:ie,js:je,ks:ke)), &
158 'minval(grid%xa%u(is:ie,js:je,ks:ke))=', &
159 minval(grid%xa%u(is:ie,js:je,ks:ke)), &
160 'maxval(grid%xa%v(is:ie,js:je,ks:ke))=', &
161 maxval(grid%xa%v(is:ie,js:je,ks:ke)), &
162 'minval(grid%xa%v(is:ie,js:je,ks:ke))=', &
163 minval(grid%xa%v(is:ie,js:je,ks:ke)), &
164 'maxval(grid%xa%t(is:ie,js:je,ks:ke))=', &
165 maxval(grid%xa%t(is:ie,js:je,ks:ke)), &
166 'minval(grid%xa%t(is:ie,js:je,ks:ke))=', &
167 minval(grid%xa%t(is:ie,js:je,ks:ke)), &
168 'maxval(grid%xa%q(is:ie,js:je,ks:ke))=', &
169 maxval(grid%xa%q(is:ie,js:je,ks:ke)), &
170 'minval(grid%xa%q(is:ie,js:je,ks:ke))=', &
171 minval(grid%xa%q(is:ie,js:je,ks:ke)), &
172 'maxval(grid%xa%p(is:ie,js:je,ks:ke))=', &
173 maxval(grid%xa%p(is:ie,js:je,ks:ke)), &
174 'minval(grid%xa%p(is:ie,js:je,ks:ke))=', &
175 minval(grid%xa%p(is:ie,js:je,ks:ke)), &
176 'maxval(grid%xa%psfc(is:ie,js:je)) =', &
177 maxval(grid%xa%psfc(is:ie,js:je)), &
178 'minval(grid%xa%psfc(is:ie,js:je)) =', &
179 minval(grid%xa%psfc(is:ie,js:je))
180 end if
181
182 do j=js,je
183 do i=is,ie
184
185 ph_full = grid%ht(i,j) * gravity
186 ph_xb_hd = grid%ht(i,j) * gravity
187 do k = ks, ke
188 ! To obtain all of the full fields: t, p, q(mixing ratio), rho
189 t_full = grid%xa%t(i,j,k) + grid%xb%t(i,j,k)
190 p_full = grid%xa%p(i,j,k) + grid%xb%p(i,j,k)
191 q_full = grid%moist(i,j,k,P_QV) + q_cgrid(i,j,k)
192
193 ! Note: According to WRF, this is the dry air density used to
194 ! compute the geopotential height:
195 rho_full = p_full / (gas_constant*t_full*(1.0+q_full/rd_over_rv))
196
197 ! To compute the theta increment with the full fields:
198 grid%em_t_2(i,j,k) = t_full*(ps0/p_full)**kappa - ts0
199
200 ! The full field of analysis ph:
201 ph_full = ph_full &
202 - grid%xb%dnw(k) * (grid%xb%psac(i,j)+mu_cgrid(i,j)) / rho_full
203
204 ! background hydrostatic phi:
205 ph_xb_hd = ph_xb_hd &
206 - grid%xb%dnw(k) * grid%xb%psac(i,j) / grid%xb%rho(i,j,k)
207
208 ! The analysis perturbation = Hydro_ph - base_ph + nonhydro_xb_ph:
209 grid%em_ph_2(i,j,k+1) = ph_full - grid%em_phb(i,j,k+1) &
210 + (grid%xb%hf(i,j,k+1)*gravity - ph_xb_hd)
211 end do
212 end do
213 end do
214
215 ! To compute the geopotential height increment:
216
217 do k=ks, ke+1
218 do j=js,je
219 do i=is,ie
220 ph_cgrid(i,j,k) = grid%em_ph_2(i,j,k) - ph_cgrid(i,j,k)
221 end do
222 end do
223 end do
224
225 ! ========================
226 ! Write out the increment:
227 ! ========================
228
229 if (write_increments) then
230 write(unit=stdout,fmt='(/"Write out increment for plotting......")')
231 call da_write_increments (grid%xp, grid%xb, grid%xa, q_cgrid, mu_cgrid, ph_cgrid)
232 end if
233
234 ! CONVERT FROM A-GRID TO C-GRID
235
236 !TBH: NOTE that grid%xp%halo_id3 = HALO_PSICHI_UV_ADJ which is currently defined
237 !TBH: in the Registry as "dyn_em 24:grid%xa%u,grid%xa%v,grid%xa%psfc". Clearly it is not
238 !TBH: necessary to update halos in grid%xa%psfc here! Also, 24-point stencil is
239 !TBH: too thick, 9-point should suffice. Apparently, grid%xp%halo_id3 is used
240 !TBH: in many places! This needs to be fixed.
241
242 ! Fill the halo region for u and v.
243
244 call wrf_dm_halo(grid%xp%domdesc,grid%xp%comms,grid%xp%halo_id3)
245
246 ! Fill the boundary
247
248 ! The southern boundary
249 if (js == grid%xp%jds) then
250 grid%xa%v(is:ie,js-1,ks:ke)=2.0*grid%xa%v(is:ie,js ,ks:ke) &
251 - grid%xa%v(is:ie,js+1,ks:ke)
252 end if
253
254 ! The northern boundary
255 if (je == grid%xp%jde) then
256 grid%xa%v(is:ie,je+1,ks:ke)=2.0*grid%xa%v(is:ie,je ,ks:ke) &
257 - grid%xa%v(is:ie,je-1,ks:ke)
258 end if
259
260 ! The western boundary
261 if (is == grid%xp%ids) then
262 grid%xa%u(is-1,js:je,ks:ke)=2.0*grid%xa%u(is ,js:je,ks:ke) &
263 - grid%xa%u(is+1,js:je,ks:ke)
264 end if
265
266 ! The eastern boundary
267 if (ie == grid%xp%ide) then
268 grid%xa%u(ie+1,js:je,ks:ke)=2.0*grid%xa%u(ie ,js:je,ks:ke) &
269 - grid%xa%u(ie-1,js:je,ks:ke)
270 end if
271
272 do k=ks,ke
273 do j=js,je+1
274 do i=is,ie+1
275 u_cgrid(i,j,k)=0.5*(grid%xa%u(i-1,j ,k)+grid%xa%u(i,j,k))
276 v_cgrid(i,j,k)=0.5*(grid%xa%v(i ,j-1,k)+grid%xa%v(i,j,k))
277 end do
278 end do
279 end do
280
281 !------------------------------------------------------------------------
282 ! For later plot and comparation Purpose only, zero out the unused var.
283 !------------------------------------------------------------------------
284
285 ! The northern boundary
286 if (je == grid%xp%jde) then
287 u_cgrid(is:ie+1,je+1,ks:ke)=0.0
288 end if
289
290 ! The eastern boundary
291 if (ie == grid%xp%ide) then
292 v_cgrid(ie+1,js:je+1,ks:ke)=0.0
293 end if
294
295 !---------------------------------------------------------------------------
296 ! [5.0] add increment to the original guess and update "grid"
297 !---------------------------------------------------------------------------
298
299 do j=js,je
300 do i=is,ie
301 grid%em_mu_2(i,j) = grid%em_mu_2(i,j) + mu_cgrid(i,j)
302 grid%em_mu0(i,j) = grid%em_mub(i,j) + grid%em_mu_2(i,j)
303 grid%em_w_2(i,j,ke+1)= grid%em_w_2(i,j,ke+1) + grid%xa%w(i,j,ke+1)
304 end do
305
306 do k=ks,ke
307 do i=is,ie
308 grid%em_u_2(i,j,k) = grid%em_u_2(i,j,k) + u_cgrid(i,j,k)
309 grid%em_v_2(i,j,k) = grid%em_v_2(i,j,k) + v_cgrid(i,j,k)
310 grid%em_w_2(i,j,k) = grid%em_w_2(i,j,k) + grid%xa%w(i,j,k)
311 grid%moist(i,j,k,P_QV) = grid%moist(i,j,k,P_QV) + q_cgrid(i,j,k)
312 ! make sure qv is positive.
313 if (grid%moist(i,j,k,P_QV) < 1.0e-6) grid%moist(i,j,k,P_QV) = 1.0e-6
314
315 if (size(grid%moist,dim=4) >= 4) then
316 grid%moist(i,j,k,p_qc) = grid%moist(i,j,k,p_qc) + grid%xa%qcw(i,j,k)
317 grid%moist(i,j,k,p_qr) = grid%moist(i,j,k,p_qr) + grid%xa%qrn(i,j,k)
318 if (grid%moist(i,j,k,p_qc) < 0.0) grid%moist(i,j,k,p_qc) = 0.0
319 if (grid%moist(i,j,k,p_qr) < 0.0) grid%moist(i,j,k,p_qr) = 0.0
320 end if
321
322 if (size(grid%moist,dim=4) >= 6) then
323 grid%moist(i,j,k,p_qi) = grid%moist(i,j,k,p_qi) + grid%xa%qci(i,j,k)
324 grid%moist(i,j,k,p_qs) = grid%moist(i,j,k,p_qs) + grid%xa%qsn(i,j,k)
325 if (grid%moist(i,j,k,p_qi) < 0.0) grid%moist(i,j,k,p_qi) = 0.0
326 if (grid%moist(i,j,k,p_qs) < 0.0) grid%moist(i,j,k,p_qs) = 0.0
327 end if
328
329 if (size(grid%moist,dim=4) >= 7) then
330 grid%moist(i,j,k,p_qg) = grid%moist(i,j,k,p_qg) + grid%xa%qgr(i,j,k)
331 if (grid%moist(i,j,k,p_qg) < 0.0) grid%moist(i,j,k,p_qg) = 0.0
332 end if
333 end do
334 end do
335 end do
336
337 ! The northern boundary
338 if (je == grid%xp%jde) then
339 j=je+1
340 do k=ks,ke
341 do i=is,ie
342 grid%em_v_2(i,j,k) = grid%em_v_2(i,j,k) + v_cgrid(i,j,k)
343 end do
344 end do
345 end if
346
347 ! The eastern boundary
348 if (ie == grid%xp%ide) then
349 i=ie+1
350 do k=ks,ke
351 do j=js,je
352 grid%em_u_2(i,j,k) = grid%em_u_2(i,j,k) + u_cgrid(i,j,k)
353 end do
354 end do
355 end if
356
357 if (print_detail_xa) then
358 write(unit=stdout, fmt=*) 'simple variables:'
359
360 if (ie == grid%xp%ide) then
361 write (unit=stdout,fmt=*) ' '
362
363 do k=ks+5,ke,10
364 do j=js,je,10
365 write(unit=stdout, fmt=*) &
366 ' grid%em_u_2(', grid%xp%ide+1, ',', j, ',', k, ')=', &
367 grid%em_u_2(grid%xp%ide+1,j,k)
368 end do
369 write(unit=stdout, fmt=*) ' '
370 end do
371 end if
372
373 if (je == grid%xp%jde) then
374 write(unit=stdout, fmt=*) ' '
375
376 do k=ks+5,ke,10
377 do i=is,ie,10
378 write(unit=stdout, fmt=*) &
379 ' grid%em_v_2(', i, ',', grid%xp%jde+1, ',', k, ')=', &
380 grid%em_v_2(i, grid%xp%jde+1,k)
381 end do
382 write(unit=stdout, fmt=*) ' '
383 end do
384 end if
385
386 write(unit=stdout, fmt='(2(2x, a, e20.12))') &
387 'maxval(mu_cgrid(is:ie,js:je)) =', &
388 maxval(mu_cgrid(is:ie,js:je)), &
389 'minval(mu_cgrid(is:ie,js:je)) =', &
390 minval(mu_cgrid(is:ie,js:je)), &
391 'maxval(u_cgrid(is:ie,js:je,ks:ke)) =', &
392 maxval(u_cgrid(is:ie,js:je,ks:ke)), &
393 'minval(u_cgrid(is:ie,js:je,ks:ke)) =', &
394 minval(u_cgrid(is:ie,js:je,ks:ke)), &
395 'maxval(v_cgrid(is:ie,js:je,ks:ke)) =', &
396 maxval(v_cgrid(is:ie,js:je,ks:ke)), &
397 'minval(v_cgrid(is:ie,js:je,ks:ke)) =', &
398 minval(v_cgrid(is:ie,js:je,ks:ke)), &
399 'maxval(q_cgrid(is:ie,js:je,ks:ke)) =', &
400 maxval(q_cgrid(is:ie,js:je,ks:ke)), &
401 'minval(q_cgrid(is:ie,js:je,ks:ke)) =', &
402 minval(q_cgrid(is:ie,js:je,ks:ke))
403
404 do k=ks,ke
405 write(unit=stdout, fmt='(a, i3)') 'k=', k
406
407 write(unit=stdout, fmt='(2(2x, a, e20.12))') &
408 'maxval(u_cgrid(is:ie,js:je,k)) =', maxval(u_cgrid(is:ie,js:je,k)), &
409 'minval(u_cgrid(is:ie,js:je,k)) =', minval(u_cgrid(is:ie,js:je,k)), &
410 'maxval(v_cgrid(is:ie,js:je,k)) =', maxval(v_cgrid(is:ie,js:je,k)), &
411 'minval(v_cgrid(is:ie,js:je,k)) =', minval(v_cgrid(is:ie,js:je,k)), &
412 'maxval(q_cgrid(is:ie,js:je,k)) =', maxval(q_cgrid(is:ie,js:je,k)), &
413 'minval(q_cgrid(is:ie,js:je,k)) =', minval(q_cgrid(is:ie,js:je,k))
414 end do
415 end if
416
417 if (trace_use) call da_trace_exit("da_transfer_xatowrf")
418
419 end subroutine da_transfer_xatowrf
420
421