da_uvprho_to_w_adj.inc
References to this file elsewhere.
1 subroutine da_uvprho_to_w_adj(xb, xa, xp)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Calculates vertical velocity increments from Richardson's Eq.
5 !
6 ! Method: Richardson's Eq., which
7 ! combines continuity Eq., thermodynamic Eq. and hrdrostatic Eq.
8 !---------------------------------------------------------------------------
9
10 implicit none
11
12 type (xb_type), intent(in) :: xb ! First guess structure.
13 type (x_type), intent(inout) :: xa ! increment structure.
14 type (xpose_type), intent(inout) :: xp ! domain decomposition vars.
15
16 integer :: is, ie ! 1st dim. end points.
17 integer :: js, je ! 2nd dim. end points.
18
19 integer :: I,J,K
20
21 real, dimension(ims:ime,jms:jme,kms:kme) :: URHO, VRHO
22 real, dimension(ims:ime,jms:jme,kms:kme) :: DIV, WZ
23 real :: TERM3
24
25 if (trace_use) call da_trace_entry("da_uvprho_to_w_adj")
26
27 !initialize to zero for some variables because of the adjoint requirements.
28
29 WZ(:,:,:)=0.
30 URHO(:,:,:)=0.0
31 VRHO(:,:,:)=0.0
32 DIV(:,:,:)=0.0
33 TERM3=0.0
34
35 ! integration to calculate the vertical velocity
36
37 do J=jts,jte
38 do I=its,ite
39 do K=kts,kte
40 xa%w(I,J,K+1)=xa%w(I,J,K+1)+xa%w(I,J,K)
41 WZ(I,J,K)=xa%w(I,J,K)*(xb%hf(I,J,K)-xb%hf(I,J,K+1))
42 xa%w(I,J,K)=0.0
43 end do
44 xa%w(I,J,kte+1)=0.0
45 end do
46 end do
47
48 call da_w_adjustment_adj(xb,WZ)
49
50 ! Divide by constant
51
52 do K=kts,kte
53 do J=jts,jte
54 do I=its,ite
55 WZ(I,J,K)=WZ(I,J,K)/(GAMMA*xb%p(I,J,K))
56 end do
57 end do
58 end do
59
60 ! Term 4: Derivative of basic vertical velocity with respect to z.
61
62 do J=jts,jte
63 do I=its,ite
64 do K=kte,kts,-1
65 xa%p(I,J,K)=xa%p(I,J,K)-WZ(I,J,K)*GAMMA* &
66 (xb%w(I,J,K+1)-xb%w(I,J,K))/ &
67 (xb%hf(I,J,K+1)-xb%hf(I,J,K))
68 end do
69 end do
70 end do
71
72 ! Term 3.2: Vertical integration of the basic mass divergence
73
74 do J=jts,jte
75 do I=its,ite
76 do K=kts,kte-1
77 TERM3=TERM3+WZ(I,J,K)
78 DIV(I,J,K+1)=DIV(I,J,K+1)+ &
79 TERM3*GRAVITY*0.5*(xb%h(I,J,K+1)-xb%h(I,J,K))
80 DIV(I,J,K) =DIV(I,J,K)+ &
81 TERM3*GRAVITY*0.5*(xb%h(I,J,K+1)-xb%h(I,J,K))
82 end do
83 TERM3=0.0
84 end do
85 end do
86
87 call da_uv_to_divergence_adj(xp, xb, URHO,VRHO, DIV)
88
89 ! Computation to check for edge of domain:
90 if (Testing_WRFVAR) then
91 is = its-1; ie = ite+1; js = jts-1; je = jte+1
92 if (its == ids) is = ids; if (ite == ide) ie = ide
93 if (jts == jds) js = jds; if (jte == jde) je = jde
94 else
95 is = its
96 ie = ite
97 js = jts
98 je = jte
99 end if
100
101 do K=kts,kte
102 do J=js,je
103 do I=is,ie
104 xa%rho(I,J,K)=xa%rho(I,J,K)+VRHO(I,J,K)*xb%v(I,J,K)
105 xa%rho(I,J,K)=xa%rho(I,J,K)+URHO(I,J,K)*xb%u(I,J,K)
106 end do
107 end do
108 end do
109
110 URHO(:,:,:)=0.0
111 VRHO(:,:,:)=0.0
112
113 ! Term 3.1: Vertical integration of the perturbed mass divergence
114
115 do J=jts,jte
116 do I=its,ite
117 do K=kts,kte-1
118 TERM3=TERM3+WZ(I,J,K)
119 DIV(I,J,K+1)=DIV(I,J,K+1)+ &
120 TERM3*GRAVITY*0.5*(xb%h(I,J,K+1)-xb%h(I,J,K))
121 DIV(I,J,K) =DIV(I,J,K)+ &
122 TERM3*GRAVITY*0.5*(xb%h(I,J,K+1)-xb%h(I,J,K))
123 end do
124 TERM3=0.0
125 end do
126 end do
127
128 call da_uv_to_divergence_adj(xp, xb, URHO,VRHO, DIV)
129
130 do K=kts,kte
131 do J=js,je
132 do I=is,ie
133 xa%v(I,J,K)=xa%v(I,J,K)+VRHO(I,J,K)*xb%rho(I,J,K)
134 xa%u(I,J,K)=xa%u(I,J,K)+URHO(I,J,K)*xb%rho(I,J,K)
135 end do
136 end do
137 end do
138
139 URHO(:,:,:)=0.0
140 VRHO(:,:,:)=0.0
141
142 ! Term 2.2: Divergence term from basic wind
143
144 call da_uv_to_divergence(xb, xb%u,xb%v, DIV)
145
146 do K=kts,kte
147 do J=jts,jte
148 do I=its,ite
149 xa%p(I,J,K)=xa%p(I,J,K)-WZ(I,J,K)*GAMMA*DIV(I,J,K)
150 end do
151 end do
152 end do
153
154 ! Term 2.1: Divergence term from perturbed wind
155
156 do K=kts,kte
157 do J=jts,jte
158 do I=its,ite
159 DIV(I,J,K)=-WZ(I,J,K)*GAMMA*xb%p(I,J,K) ! DIV redefined
160 end do
161 end do
162 end do
163
164 call da_uv_to_divergence_adj(xp, xb, xa%u,xa%v, DIV)
165
166 ! Computation to check for edge of domain:
167 is = its
168 ie = ite
169 js = jts
170 je = jte
171 if (its == ids) is = ids+1
172 if (ite == ide) ie = ide-1
173 if (jts == jds) js = jds+1
174 if (jte == jde) je = jde-1
175
176 ! Term 1.2: Basic pressure advection along the perturbed wind
177
178 if (jte == jde) then
179 j = jte
180 do K=kts,kte
181 do I=its, ite
182 WZ(I,J-1,K)=WZ(I,J-1,K)+WZ(I,J,K)
183 end do
184 end do
185 end if
186
187 if (jts == jds) then
188 j = jts
189 do K=kts,kte
190 do I=its, ite
191 WZ(I,J+1,K)=WZ(I,J+1,K)+WZ(I,J,K)
192 end do
193 end do
194 end if
195
196 if (ite == ide) then
197 i = ite
198 do K=kts,kte
199 do J=js,je
200 WZ(I-1,J,K)=WZ(I-1,J,K)+WZ(I,J,K)
201 end do
202 end do
203 end if
204
205 if (its == ids) then
206 i = its
207 do K=kts,kte
208 do J=js,je
209 WZ(I+1,J,K)=WZ(I+1,J,K)+WZ(I,J,K)
210 end do
211 end do
212 end if
213
214 do K=kts,kte
215 do J=js,je
216 do I=is,ie
217 xa%v(I,J,K)=xa%v(I,J,K)-WZ(I,J,K)* &
218 (xb%p(I,J+1,K)-xb%p(I,J-1,K))*xb%coefy(I,J)
219 xa%u(I,J,K)=xa%u(I,J,K)-WZ(I,J,K)* &
220 (xb%p(I+1,J,K)-xb%p(I-1,J,K))*xb%coefx(I,J)
221 end do
222 end do
223 end do
224
225 !-------------------------------------------------------------------------
226 ! Computation to check for edge of domain:
227 ! This is only for adjoint, as we have to cross the processor boundary
228 ! to get the contribution.
229
230 is = its - 1
231 ie = ite + 1
232 js = jts - 1
233 je = jte + 1
234
235 xp%v1z(its:ite, jts:jte, kts:kte) = wz(its:ite, jts:jte, kts:kte)
236
237 call wrf_dm_halo(xp%domdesc,xp%comms,xp%halo_id2)
238
239 if (its == ids) then
240 is = ids+1
241 else
242 wz(is, js:je, kts:kte) = xp%v1z(is, js:je, kts:kte)
243 end if
244
245 if (ite == ide) then
246 ie = ide-1
247 else
248 wz(ie, js:je, kts:kte) = xp%v1z(ie, js:je, kts:kte)
249 end if
250
251 if (jts == jds) then
252 js = jds+1
253 else
254 wz(is:ie, js, kts:kte) = xp%v1z(is:ie, js, kts:kte)
255 end if
256
257 if (jte == jde) then
258 je = jde-1
259 else
260 wz(is:ie, je, kts:kte) = xp%v1z(is:ie, je, kts:kte)
261 end if
262
263 ! Term 1.1: Perturbed pressure advection along the basic wind
264
265 do K=kts,kte
266 do J=js,je
267 do I=is,ie
268 xa%p(I,J+1,K)=xa%p(I,J+1,K)-WZ(I,J,K)*xb%v(I,J,K)*xb%coefy(I,J)
269 xa%p(I,J-1,K)=xa%p(I,J-1,K)+WZ(I,J,K)*xb%v(I,J,K)*xb%coefy(I,J)
270 xa%p(I+1,J,K)=xa%p(I+1,J,K)-WZ(I,J,K)*xb%u(I,J,K)*xb%coefx(I,J)
271 xa%p(I-1,J,K)=xa%p(I-1,J,K)+WZ(I,J,K)*xb%u(I,J,K)*xb%coefx(I,J)
272 end do
273 end do
274 end do
275
276 WZ(:,:,:) = 0.0
277
278 if (trace_use) call da_trace_exit("da_uvprho_to_w_adj")
279
280 end subroutine da_uvprho_to_w_adj
281
282