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