da_transfer_xatowrftl.inc
References to this file elsewhere.
1 subroutine da_transfer_xatowrftl(grid, config_flags, filnam)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Convert analysis increments into WRFTL increments
5 ! (following xatowrf, but only keep the increments)
6 !---------------------------------------------------------------------------
7
8 implicit none
9
10 type(domain), intent(inout) :: grid
11 type(grid_config_rec_type), intent(inout) :: config_flags
12
13 character*4, intent(in) :: filnam
14
15 integer :: i, j, k
16 integer :: is, ie, js, je, ks, ke
17 real :: sdmd, s1md
18 real :: rho_cgrid
19
20 real :: g_press(grid%xp%ims:grid%xp%ime,grid%xp%jms:grid%xp%jme, &
21 grid%xp%kms:grid%xp%kme)
22
23 integer ndynopt
24
25 if (trace_use) call da_trace_entry("da_transfer_xatowrftl")
26
27 is=grid%xp%its
28 ie=grid%xp%ite
29 js=grid%xp%jts
30 je=grid%xp%jte
31 ks=grid%xp%kts
32 ke=grid%xp%kte
33
34 !---------------------------------------------------------------------------
35 ! [1.0] Get the mixing ratio of moisture first, as it is easy.
36 !---------------------------------------------------------------------------
37
38 do k=ks,ke
39 do j=js,je
40 do i=is,ie
41 grid%g_moist(i,j,k,P_G_QV)=grid%xa%q(i,j,k)/(1.0-grid%xb%q(i,j,k))**2
42 end do
43 end do
44 end do
45
46 !---------------------------------------------------------------------------
47 ! [2.0] COMPUTE increments of dry-column air mass per unit area
48 !---------------------------------------------------------------------------
49
50 do j=js,je
51 do i=is,ie
52 sdmd=0.0
53 s1md=0.0
54 do k=ks,ke
55 sdmd=sdmd+grid%g_moist(i,j,k,P_G_QV)*grid%em_dnw(k)
56 s1md=s1md+(1.0+grid%moist(i,j,k,P_QV))*grid%em_dnw(k)
57 end do
58 grid%em_g_mu_2(i,j)=-(grid%xa%psfc(i,j)+grid%xb%psac(i,j)*sdmd)/s1md
59 end do
60 end do
61
62 !---------------------------------------------------------------------------
63 ! [3.0] compute pressure increments (for computing theta increments)
64 !---------------------------------------------------------------------------
65
66 do j=js,je
67 do i=is,ie
68 g_press(i,j,ke+1)=0.0
69 do k=ke,ks,-1
70 g_press(i,j,k)=g_press(i,j,k+1) &
71 -(grid%em_g_mu_2(i,j)*(1.0+grid%moist(i,j,k,P_QV)) &
72 +(grid%em_mu_2(i,j)+grid%em_mub(i,j))*grid%g_moist(i,j,k,P_G_QV))* &
73 grid%em_dn(k)
74 grid%xa%p(i,j,k)=0.5*(g_press(i,j,k)+g_press(i,j,k+1))
75 end do
76 end do
77 end do
78
79 !---------------------------------------------------------------------------
80 ! [4.0] convert temperature increments into theta increments
81 ! evaluate also the increments of (1/rho) and geopotential
82 !---------------------------------------------------------------------------
83
84 do k=ks,ke
85 do j=js,je
86 do i=is,ie
87 grid%em_g_t_2(i,j,k)=(300.0+grid%em_t_2(i,j,k))* &
88 (grid%xa%t(i,j,k)/grid%xb%t(i,j,k) &
89 -kappa*grid%xa%p(i,j,k)/grid%xb%p(i,j,k))
90 end do
91 end do
92 end do
93 do j=js,je
94 do i=is,ie
95 grid%em_g_ph_2(i,j,ks)=0.
96 do k=ks,ke
97 rho_cgrid=grid%xb%rho(i,j,k) &
98 *(grid%xa%p(i,j,k)/grid%xb%p(i,j,k) &
99 -grid%xa%t(i,j,k)/grid%xb%t(i,j,k) &
100 -0.61*grid%xa%q(i,j,k)/(1.+0.61*grid%xb%q(i,j,k)))
101 grid%em_g_ph_2(i,j,k+1)=grid%em_g_ph_2(i,j,k) &
102 -(g_press(i,j,k+1)-g_press(i,j,k) &
103 +(grid%em_ph_2(i,j,k+1)-grid%em_ph_2(i,j,k))*rho_cgrid) &
104 /grid%xb%rho(i,j,k)
105 end do
106 end do
107 end do
108
109 !---------------------------------------------------------------------------
110 ! [5.0] convert from a-grid to c-grid
111 !---------------------------------------------------------------------------
112
113 #ifdef DM_PARALLEL
114 ! Fill the halo region for u and v.
115
116 call wrf_dm_halo(grid%xp%domdesc,grid%xp%comms,grid%xp%halo_id3)
117
118 ! The southern boundary (fill A-GRID boundaries)
119 ! To keep the gradient, A(0) = 2A(1)-A(2)
120 if (js == grid%xp%jds) &
121 grid%xa%v(is:ie,js-1,ks:ke)=2.0*grid%xa%v(is:ie,js,ks:ke)-grid%xa%v(is:ie,js+1,ks:ke)
122
123 ! The western boundary (fill A-GRID boundaries)
124 ! To keep the gradient, A(0) = 2A(1)-A(2)
125 if (is == grid%xp%ids) &
126 grid%xa%u(is-1,js:je,ks:ke)=2.0*grid%xa%u(is,js:je,ks:ke)-grid%xa%u(is+1,js:je,ks:ke)
127
128 do k=ks,ke
129 do j=js,je
130 do i=is,ie
131 grid%em_g_u_2(i,j,k)=0.5*(grid%xa%u(i-1,j ,k)+grid%xa%u(i,j,k))
132 grid%em_g_v_2(i,j,k)=0.5*(grid%xa%v(i ,j-1,k)+grid%xa%v(i,j,k))
133 end do
134 end do
135 end do
136
137 ! The northern boundary (compute C-GRID boundaries)
138 ! To keep the gradient, A(N+1) = 2A(N)-A(N-1)
139 ! and on C-Grid, this will lead to C(N+1)=(A(N)+A(N+1))/2=(3A(N)-A(N-1))/2
140 if (je == grid%xp%jde) &
141 grid%em_g_v_2(is:ie,je+1,ks:ke)=&
142 (3.0*grid%xa%v(is:ie,je,ks:ke)-grid%xa%v(is:ie,je-1,ks:ke))/2.
143
144 ! The eastern boundary (compute C-GRID boundaries)
145 ! To keep the gradient, A(N+1) = 2A(N)-A(N-1)
146 ! and on C-Grid, this will lead to C(N+1)=(A(N)+A(N+1))/2=(3A(N)-A(N-1))/2
147 if (ie == grid%xp%ide) &
148 grid%em_g_u_2(ie+1,js:je,ks:ke)=&
149 (3.0*grid%xa%u(ie,js:je,ks:ke)- grid%xa%u(ie-1,js:je,ks:ke))/2.
150
151 #else
152
153 do k=ks,ke
154 do j=js,je
155 do i=is+1,ie
156 grid%em_g_u_2(i,j,k)=0.5*(grid%xa%u(i-1,j,k)+grid%xa%u(i,j,k))
157 end do
158 end do
159 do j=js+1,je
160 do i=is,ie
161 grid%em_g_v_2(i,j,k)=0.5*(grid%xa%v(i,j-1,k)+grid%xa%v(i,j,k))
162 end do
163 end do
164 end do
165
166 ! To keep the gradient, A(N+1) = 2A(N)-A(N-1)
167 ! and on C-Grid, this will lead to C(N+1)=(A(N)+A(N+1))/2=(3A(N)-A(N-1))/2
168
169 ! The eastern boundary
170 grid%em_g_u_2(ie+1,js:je,ks:ke)=(3.0*grid%xa%u(ie,js:je,ks:ke)-grid%xa%u(ie-1,js:je,ks:ke))/2.
171
172 ! The northern boundary
173 grid%em_g_v_2(is:ie,je+1,ks:ke)=(3.0*grid%xa%v(is:ie,je,ks:ke)-grid%xa%v(is:ie,je-1,ks:ke))/2.
174
175 ! To keep the gradient, A(0) = 2A(1)-A(2)
176 ! and on C-Grid, this will lead to C(1)=(A(0)+A(1))/2=(3A(1)-A(2))/2
177
178 ! The western boundary
179 grid%em_g_u_2(is,js:je,ks:ke)=(3.0*grid%xa%u(is,js:je,ks:ke)-grid%xa%u(is+1,js:je,ks:ke))/2.
180
181 ! The southern boundary
182 grid%em_g_v_2(is:ie,js,ks:ke)=(3.0*grid%xa%v(is:ie,js,ks:ke)-grid%xa%v(is:ie,js+1,ks:ke))/2.
183
184 #endif
185
186 !---------------------------------------------------------------------------
187 ! [6.0] save OTHERinCREMENT
188 !---------------------------------------------------------------------------
189
190 do j=js,je
191 do k=ks,ke+1
192 do i=is,ie
193 grid%em_g_w_2(i,j,k)=grid%xa%w(i,j,k)
194 end do
195 end do
196 end do
197
198 !---------------------------------------------------------------------------
199 ! [5.0] save inCREMENT
200 !---------------------------------------------------------------------------
201
202 #ifdef VAR4D_MICROPHYSICS
203 ! New code needed once we introduce the microphysics to 4dvar in 2008
204 if (size(moist,dim=4) >= 4) then
205 do k=ks,ke
206 do j=js,je
207 do i=is,ie
208 g_moist(i,j,k,p_qcw) = grid%xa%qcw(i,j,k)
209 g_moist(i,j,k,p_qrn) = grid%xa%qrn(i,j,k)
210 end do
211 end do
212 end do
213 end if
214
215 if (size(moist,dim=4) >= 6) then
216 do k=ks,ke
217 do j=js,je
218 do i=is,ie
219 g_moist(i,j,k,p_qci) = grid%xa%qci(i,j,k)
220 g_moist(i,j,k,p_qsn) = grid%xa%qsn(i,j,k)
221 end do
222 end do
223 end do
224 end if
225
226 if (size(moist,dim=4) >= 7) then
227 do k=ks,ke
228 do j=js,je
229 do i=is,ie
230 g_moist(i,j,k,p_qgr) = grid%xa%qgr(i,j,k)
231 end do
232 end do
233 end do
234 end if
235
236 #endif
237
238 !---------------------------------------------------------------------------
239 ! [7.0] output
240 !---------------------------------------------------------------------------
241
242 ndynopt = grid%dyn_opt
243 grid%dyn_opt = DYN_EM_TL
244 call nl_set_dyn_opt (1 , DYN_EM_TL)
245
246 call da_med_initialdata_output(grid , config_flags, filnam)
247
248 grid%dyn_opt = ndynopt
249 call nl_set_dyn_opt (1 , DYN_EM)
250
251 if (trace_use) call da_trace_exit("da_transfer_xatowrftl")
252
253 end subroutine da_transfer_xatowrftl
254
255