da_uvprho_to_w_lin.inc
References to this file elsewhere.
1 subroutine da_uvprho_to_w_lin(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
26 ! Computation to check for edge of domain:
27 is = its
28 ie = ite
29 js = jts
30 je = jte
31 if (its == ids) is = ids+1
32 if (ite == ide) ie = ide-1
33 if (jts == jds) js = jds+1
34 if (jte == jde) je = jde-1
35
36 WZ(:,:,:) = 0.0
37 ! Term 1.1: perturbed pressure advection along the basic wind
38 do K=kts,kte
39 do J=js,je
40 do I=is,ie
41 WZ(I,J,K)=WZ(I,J,K)-xb%u(I,J,K)*(xa%p(I+1,J,K)-xa%p(I-1,J,K))* &
42 xb%coefx(I,J)
43 WZ(I,J,K)=WZ(I,J,K)-xb%v(I,J,K)*(xa%p(I,J+1,K)-xa%p(I,J-1,K))* &
44 xb%coefy(I,J)
45 end do
46 end do
47 end do
48
49 ! Term 1.2: Basic pressure advection along the perturbed wind
50
51 do K=kts,kte
52 do J=js,je
53 do I=is,ie
54 WZ(I,J,K)=WZ(I,J,K)-xa%u(I,J,K)*(xb%p(I+1,J,K)-xb%p(I-1,J,K))* &
55 xb%coefx(I,J)
56 WZ(I,J,K)=WZ(I,J,K)-xa%v(I,J,K)*(xb%p(I,J+1,K)-xb%p(I,J-1,K))* &
57 xb%coefy(I,J)
58 end do
59 end do
60 end do
61
62 ! Dealing the laterial boundary because of the advection.
63 ! boundary too simple? (It is the same as fill in interpf, fill can be used)
64
65 if (its == ids) then
66 i = its
67 do K=kts,kte
68 do J=js,je
69 WZ(I,J,K)=WZ(I+1,J,K)
70 end do
71 end do
72 end if
73
74 if (ite == ide) then
75 i = ite
76 do K=kts,kte
77 do J=js,je
78 WZ(I,J,K)=WZ(I-1,J,K)
79 end do
80 end do
81 end if
82
83 if (jts == jds) then
84 j = jts
85 do K=kts,kte
86 do I=its, ite
87 WZ(I,J,K)=WZ(I,J+1,K)
88 end do
89 end do
90 end if
91
92 if (jte == jde) then
93 j = jte
94 do K=kts,kte
95 do I=its, ite
96 WZ(I,J,K)=WZ(I,J-1,K)
97 end do
98 end do
99 end if
100
101 ! Term 2.1: Divergence term from perturbed wind
102
103 call da_uv_to_divergence(xb, xa%u, xa%v, DIV)
104
105 do K=kts,kte
106 do J=jts,jte
107 do I=its,ite
108 WZ(I,J,K)=WZ(I,J,K)-GAMMA*xb%p(I,J,K)*DIV(I,J,K)
109 end do
110 end do
111 end do
112
113 ! Term 2.2: Divergence term from basic wind
114
115 call da_uv_to_divergence(xb, xb%u, xb%v, DIV)
116
117 do K=kts,kte
118 do J=jts,jte
119 do I=its,ite
120 WZ(I,J,K)=WZ(I,J,K)-GAMMA*xa%p(I,J,K)*DIV(I,J,K)
121 end do
122 end do
123 end do
124
125 ! Computation to check for edge of domain:
126 is = its-1; ie = ite+1; js = jts-1; je = jte+1
127 if (its == ids) is = ids; if (ite == ide) ie = ide
128 if (jts == jds) js = jds; if (jte == jde) je = jde
129
130 ! Term 3.1: Vertical integration of the perturbed mass divergence
131
132 do K=kts,kte
133 do J=js,je
134 do I=is,ie
135 URHO(I,J,K)=xb%rho(I,J,K)*xa%u(I,J,K)
136 VRHO(I,J,K)=xb%rho(I,J,K)*xa%v(I,J,K)
137 end do
138 end do
139 end do
140
141 call da_uv_to_divergence(xb, URHO, VRHO, DIV)
142
143 do J=jts,jte
144 do I=its,ite
145 TERM3=0.0
146
147 do K=kte-1,kts,-1
148 TERM3=TERM3+GRAVITY*(DIV(I,J,K+1)+DIV(I,J,K))*0.5 &
149 *(xb%h(I,J,K+1)-xb%h(I,J,K))
150 WZ(I,J,K)=WZ(I,J,K)+TERM3
151 end do
152 end do
153 end do
154
155 ! Term 3.2: Vertical integration of the basic mass divergence
156
157 do K=kts,kte
158 do J=js,je
159 do I=is,ie
160 URHO(I,J,K)=xa%rho(I,J,K)*xb%u(I,J,K)
161 VRHO(I,J,K)=xa%rho(I,J,K)*xb%v(I,J,K)
162 end do
163 end do
164 end do
165
166 call da_uv_to_divergence(xb, URHO, VRHO, DIV)
167
168 do J=jts,jte
169 do I=its,ite
170 TERM3=0.0
171
172 do K=kte-1,kts,-1
173 TERM3=TERM3+GRAVITY*(DIV(I,J,K+1)+DIV(I,J,K))*0.5 &
174 *(xb%h(I,J,K+1)-xb%h(I,J,K))
175 WZ(I,J,K)=WZ(I,J,K)+TERM3
176 end do
177 end do
178 end do
179
180 ! Term 4: Derivative of basic vertical velocity with respect to z.
181
182 do J=jts,jte
183 do I=its,ite
184 do K=kts,kte
185 WZ(I,J,K)=WZ(I,J,K)-GAMMA*xa%p(I,J,K)*(xb%w(I,J,K+1)-xb%w(I,J,K))/ &
186 (xb%hf(I,J,K+1)-xb%hf(I,J,K))
187 end do
188 end do
189 end do
190
191 ! Divided by constant
192
193 do K=kts,kte
194 do J=jts,jte
195 do I=its,ite
196 WZ(I,J,K)=WZ(I,J,K)/(GAMMA*xb%p(I,J,K))
197 end do
198 end do
199 end do
200
201 ! integration to calculate the vertical velocity
202
203 call da_w_adjustment_lin(xb,xa%w,WZ)
204
205 do J=jts,jte
206 do I=its,ite
207 xa%w(I,J,kte+1)=0.0
208 do K=kte,kts,-1
209 xa%w(I,J,K)=xa%w(I,J,K+1)+ &
210 WZ(I,J,K)*(xb%hf(I,J,K)-xb%hf(I,J,K+1))
211 end do
212 end do
213 end do
214
215 end subroutine da_uvprho_to_w_lin
216
217