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