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