da_uvprho_to_w_adj.inc
References to this file elsewhere.
1 subroutine da_uvprho_to_w_adj(grid)
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 (domain), intent(inout) :: grid
13
14 integer :: is, ie ! 1st dim. end points.
15 integer :: js, je ! 2nd dim. end points.
16
17 integer :: I,J,K
18
19 real, dimension(ims:ime,jms:jme,kms:kme) :: URHO, VRHO
20 real, dimension(ims:ime,jms:jme,kms:kme) :: DIV, WZ
21 real :: TERM3
22
23 if (trace_use) call da_trace_entry("da_uvprho_to_w_adj")
24
25 ! initialize to zero for some variables because of the adjoint requirements.
26
27 WZ(:,:,:) = 0.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 grid%xa%w(I,J,K+1)=grid%xa%w(I,J,K+1)+grid%xa%w(I,J,K)
39 WZ(I,J,K)=grid%xa%w(I,J,K)*(grid%xb%hf(I,J,K)-grid%xb%hf(I,J,K+1))
40 grid%xa%w(I,J,K)=0.0
41 end do
42 grid%xa%w(I,J,kte+1)=0.0
43 end do
44 end do
45
46 call da_w_adjustment_adj(grid%xb,WZ)
47
48 ! Divide by constant
49
50 WZ(its:ite,jts:jte,kts:kte)=WZ(its:ite,jts:jte,kts:kte)/(GAMMA*grid%xb%p(its:ite,jts:jte,kts:kte))
51
52 ! Term 4: Derivative of basic vertical velocity with respect to z.
53
54 do J=jts,jte
55 do I=its,ite
56 do K=kte,kts,-1
57 grid%xa%p(I,J,K)=grid%xa%p(I,J,K)-WZ(I,J,K)*GAMMA* &
58 (grid%xb%w(I,J,K+1)-grid%xb%w(I,J,K))/ &
59 (grid%xb%hf(I,J,K+1)-grid%xb%hf(I,J,K))
60 end do
61 end do
62 end do
63
64 ! Term 3.2: Vertical integration of the basic mass divergence
65
66 do J=jts,jte
67 do I=its,ite
68 do K=kts,kte-1
69 TERM3=TERM3+WZ(I,J,K)
70 DIV(I,J,K+1)=DIV(I,J,K+1)+ &
71 TERM3*GRAVITY*0.5*(grid%xb%h(I,J,K+1)-grid%xb%h(I,J,K))
72 DIV(I,J,K) =DIV(I,J,K)+ &
73 TERM3*GRAVITY*0.5*(grid%xb%h(I,J,K+1)-grid%xb%h(I,J,K))
74 end do
75 TERM3=0.0
76 end do
77 end do
78
79 call da_uv_to_divergence_adj(grid, URHO,VRHO, DIV)
80
81 ! Computation to check for edge of domain:
82 if (test_transforms) then
83 is = its-1; ie = ite+1; js = jts-1; je = jte+1
84 if (its == ids) is = ids; if (ite == ide) ie = ide
85 if (jts == jds) js = jds; if (jte == jde) je = jde
86 else
87 is = its
88 ie = ite
89 js = jts
90 je = jte
91 end if
92
93 grid%xa%rho(is:ie,js:je,kts:kte)=grid%xa%rho(is:ie,js:je,kts:kte)+VRHO(is:ie,js:je,kts:kte)*grid%xb%v(is:ie,js:je,kts:kte)
94 grid%xa%rho(is:ie,js:je,kts:kte)=grid%xa%rho(is:ie,js:je,kts:kte)+URHO(is:ie,js:je,kts:kte)*grid%xb%u(is:ie,js:je,kts:kte)
95
96 URHO(:,:,:)=0.0
97 VRHO(:,:,:)=0.0
98
99 ! Term 3.1: Vertical integration of the perturbed mass divergence
100
101 do J=jts,jte
102 do I=its,ite
103 do K=kts,kte-1
104 TERM3=TERM3+WZ(I,J,K)
105 DIV(I,J,K+1)=DIV(I,J,K+1)+ &
106 TERM3*GRAVITY*0.5*(grid%xb%h(I,J,K+1)-grid%xb%h(I,J,K))
107 DIV(I,J,K) =DIV(I,J,K)+ &
108 TERM3*GRAVITY*0.5*(grid%xb%h(I,J,K+1)-grid%xb%h(I,J,K))
109 end do
110 TERM3=0.0
111 end do
112 end do
113
114 call da_uv_to_divergence_adj(grid, URHO,VRHO, DIV)
115
116 grid%xa%v(is:ie,js:je,kts:kte)=grid%xa%v(is:ie,js:je,kts:kte)+VRHO(is:ie,js:je,kts:kte)*grid%xb%rho(is:ie,js:je,kts:kte)
117 grid%xa%u(is:ie,js:je,kts:kte)=grid%xa%u(is:ie,js:je,kts:kte)+URHO(is:ie,js:je,kts:kte)*grid%xb%rho(is:ie,js:je,kts:kte)
118
119 URHO(:,:,:)=0.0
120 VRHO(:,:,:)=0.0
121
122 ! Term 2.2: Divergence term from basic wind
123
124 call da_uv_to_divergence(grid%xb, grid%xb%u,grid%xb%v, DIV)
125
126 grid%xa%p(its:ite,jts:jte,kts:kte)=grid%xa%p(its:ite,jts:jte,kts:kte)-WZ(its:ite,jts:jte,kts:kte)*GAMMA*DIV(its:ite,jts:jte,kts:kte)
127
128 ! Term 2.1: Divergence term from perturbed wind
129
130 DIV(its:ite,jts:jte,kts:kte)=-WZ(its:ite,jts:jte,kts:kte)*GAMMA*grid%xb%p(its:ite,jts:jte,kts:kte) ! DIV redefined
131
132 call da_uv_to_divergence_adj(grid, grid%xa%u,grid%xa%v, DIV)
133
134 ! Computation to check for edge of domain:
135 is = its
136 ie = ite
137 js = jts
138 je = jte
139 if (its == ids) is = ids+1
140 if (ite == ide) ie = ide-1
141 if (jts == jds) js = jds+1
142 if (jte == jde) je = jde-1
143
144 ! Term 1.2: Basic pressure advection along the perturbed wind
145
146 if (jte == jde) then
147 j = jte
148 do K=kts,kte
149 do I=its, ite
150 WZ(I,J-1,K)=WZ(I,J-1,K)+WZ(I,J,K)
151 end do
152 end do
153 end if
154
155 if (jts == jds) then
156 j = jts
157 do K=kts,kte
158 do I=its, ite
159 WZ(I,J+1,K)=WZ(I,J+1,K)+WZ(I,J,K)
160 end do
161 end do
162 end if
163
164 if (ite == ide) then
165 i = ite
166 do K=kts,kte
167 do J=js,je
168 WZ(I-1,J,K)=WZ(I-1,J,K)+WZ(I,J,K)
169 end do
170 end do
171 end if
172
173 if (its == ids) then
174 i = its
175 do K=kts,kte
176 do J=js,je
177 WZ(I+1,J,K)=WZ(I+1,J,K)+WZ(I,J,K)
178 end do
179 end do
180 end if
181
182 do K=kts,kte
183 do J=js,je
184 do I=is,ie
185 grid%xa%v(I,J,K)=grid%xa%v(I,J,K)-WZ(I,J,K)* &
186 (grid%xb%p(I,J+1,K)-grid%xb%p(I,J-1,K))*grid%xb%coefy(I,J)
187 grid%xa%u(I,J,K)=grid%xa%u(I,J,K)-WZ(I,J,K)* &
188 (grid%xb%p(I+1,J,K)-grid%xb%p(I-1,J,K))*grid%xb%coefx(I,J)
189 end do
190 end do
191 end do
192
193 !-------------------------------------------------------------------------
194 ! Computation to check for edge of domain:
195 ! This is only for adjoint, as we have to cross the processor boundary
196 ! to get the contribution.
197
198 is = its - 1
199 ie = ite + 1
200 js = jts - 1
201 je = jte + 1
202
203 grid%xp%v1z(its:ite, jts:jte, kts:kte) = wz(its:ite, jts:jte, kts:kte)
204
205 #ifdef DM_PARALLEL
206 #include "HALO_BAL_EQN_ADJ.inc"
207 #endif
208
209 if (its == ids) then
210 is = ids+1
211 else
212 wz(is, js:je, kts:kte) = grid%xp%v1z(is, js:je, kts:kte)
213 end if
214
215 if (ite == ide) then
216 ie = ide-1
217 else
218 wz(ie, js:je, kts:kte) = grid%xp%v1z(ie, js:je, kts:kte)
219 end if
220
221 if (jts == jds) then
222 js = jds+1
223 else
224 wz(is:ie, js, kts:kte) = grid%xp%v1z(is:ie, js, kts:kte)
225 end if
226
227 if (jte == jde) then
228 je = jde-1
229 else
230 wz(is:ie, je, kts:kte) = grid%xp%v1z(is:ie, je, kts:kte)
231 end if
232
233 ! Term 1.1: Perturbed pressure advection along the basic wind
234
235 do K=kts,kte
236 do J=js,je
237 do I=is,ie
238 grid%xa%p(I,J+1,K)=grid%xa%p(I,J+1,K)-WZ(I,J,K)*grid%xb%v(I,J,K)*grid%xb%coefy(I,J)
239 grid%xa%p(I,J-1,K)=grid%xa%p(I,J-1,K)+WZ(I,J,K)*grid%xb%v(I,J,K)*grid%xb%coefy(I,J)
240 grid%xa%p(I+1,J,K)=grid%xa%p(I+1,J,K)-WZ(I,J,K)*grid%xb%u(I,J,K)*grid%xb%coefx(I,J)
241 grid%xa%p(I-1,J,K)=grid%xa%p(I-1,J,K)+WZ(I,J,K)*grid%xb%u(I,J,K)*grid%xb%coefx(I,J)
242 end do
243 end do
244 end do
245
246 WZ(:,:,:) = 0.0
247
248 if (trace_use) call da_trace_exit("da_uvprho_to_w_adj")
249
250 end subroutine da_uvprho_to_w_adj
251
252