da_balance_cycloterm.inc
References to this file elsewhere.
1 subroutine da_balance_cycloterm (xb, k, term_x, term_y)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Calculates cyclostrophic term in balance equation.
5 !
6 ! Method: Term is rho (u.grad) u on a single level.
7 !---------------------------------------------------------------------------
8
9 implicit none
10
11 type(xb_type), intent(in) :: xb ! First guess structure.
12 integer, intent(in) :: k ! Model level.
13 real, intent(inout) :: term_x(:,:) ! x component of term.
14 real, intent(inout) :: term_y(:,:) ! y component of term.
15
16 integer :: i, j ! Loop counters.
17 integer :: is, ie ! 1st dim. end points.
18 integer :: js, je ! 2nd dim. end points.
19
20 real :: data(ims:ime,jms:jme) ! Temporary storage.
21
22 real :: varb
23
24 if (trace_use) call da_trace_entry("da_balance_cycloterm")
25
26 !---------------------------------------------------------------------------
27 ! [1.0] Initialise:
28 !---------------------------------------------------------------------------
29
30 ! Computation to check for edge of domain:
31 is = its; ie = ite; js = jts; je = jte
32 if (.not. global .and. its==ids) is = ids+1
33 if (.not. global .and. ite==ide) ie = ide-1
34 if (jts==jds) js = jds+1
35 if (jte==jde) je = jde-1
36
37 !---------------------------------------------------------------------------
38 ! [2.0] Calculate term_x = rho M (u.du/dx + v.du/dy):
39 !---------------------------------------------------------------------------
40
41 ! [2.1] Interior points:
42
43 do j = js, je
44 do i = is, ie
45 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)*(xb%u(i+1,j,k) - &
46 xb%u(i-1,j,k)) + xb%v(i,j,k) * xb%coefy(i,j)*(xb%u(i,j+1,k) - &
47 xb%u(i,j-1,k))
48 end do
49 end do
50
51 if (.NOT. global) then ! For global only interior points
52
53 ! [2.2] Bottom boundaries:
54
55 if (its==ids) then
56 i = its
57
58 do j = js, je
59 varb = -3.0*xb%u(i,j,k)+4.0*xb%u(i+1,j,k)-xb%u(i+2,j,k)
60
61 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * varb + &
62 xb%v(i,j,k) * xb%coefy(i,j) * (xb%u(i,j+1,k) - xb%u(i,j-1,k))
63 end do
64 end if
65
66 ! [2.3] Top boundaries:
67
68 if (ite==ide) then
69 i = ite
70
71 do j = js, je
72 varb = 3.0*xb%u(i,j,k)-4.0*xb%u(i-1,j,k)+xb%u(i-2,j,k)
73
74 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * varb + &
75 xb%v(i,j,k) * xb%coefy(i,j) * (xb%u(i,j+1,k) - xb%u(i,j-1,k))
76 end do
77 end if
78
79 ! [2.4] Left boundaries:
80
81 if (jts==jds) then
82 j = jts
83
84 do i = is, ie
85 varb = -3.0*xb%u(i,j,k)+4.0*xb%u(i,j+1,k)-xb%u(i,j+2,k)
86
87 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * (xb%u(i+1,j,k) - &
88 xb%u(i-1,j,k)) + xb%v(i,j,k) * xb%coefy(i,j) * varb
89 end do
90 end if
91
92 ! [2.5] Right boundaries:
93
94 if (jte==jde) then
95 j = jte
96
97 do i = is, ie
98 varb = 3.0*xb%u(i,j,k)-4.0*xb%u(i,j-1,k)+xb%u(i,j-2,k)
99
100 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * (xb%u(i+1,j,k) - &
101 xb%u(i-1,j,k)) + xb%v(i,j,k) * xb%coefy(i,j) * varb
102 end do
103 end if
104
105 ! [2.6] Corner points:
106
107 if (its==ids .AND. jts==jds) then
108 data(its,jts) = 0.5 * (data(its,jts+1) + data(its+1,jts))
109 end if
110
111 if (ite==ide .AND. jts==jds) then
112 data(ite,jts) = 0.5 * (data(ite-1,jts) + data(ite,jts+1))
113 end if
114
115 if (its==ids .AND. jte==jde) then
116 data(its,jde) = 0.5 * (data(its,jde-1) + data(its+1,jde))
117 end if
118
119 if (ite==ide .AND. jte==jde) then
120 data(ite,jte) = 0.5 * (data(ite-1,jte) + data(ite,jte-1))
121 end if
122 end if ! not global
123
124 ! [2.7] Multiply by rho and add to term_x:
125
126 term_x(its:ite,jts:jte) = xb%rho(its:ite,jts:jte,k)*data(its:ite,jts:jte) + term_x(its:ite,jts:jte)
127
128 !---------------------------------------------------------------------------
129 ! [3.0] Calculate term_y = rho M (u.dv/dx + v.dv/dy):
130 !---------------------------------------------------------------------------
131
132 ! [3.1] Interior points:
133
134 do j = js, je
135 do i = is, ie
136 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)*(xb%v(i+1,j,k) - xb%v(i-1,j,k)) + &
137 xb%v(i,j,k) * xb%coefy(i,j)*(xb%v(i,j+1,k) - xb%v(i,j-1,k))
138 end do
139 end do
140
141
142 if (.NOT. global) then ! For global only interior points
143
144 ! [3.2] Bottom boundaries:
145
146 if (its==ids) then
147 i = its
148
149 do j = js, je
150 varb = -3.0*xb%v(i,j,k)+4.0*xb%v(i+1,j,k)-xb%v(i+2,j,k)
151
152 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)* varb + &
153 xb%v(i,j,k) * xb%coefy(i,j)* (xb%v(i,j+1,k) - xb%v(i,j-1,k))
154 end do
155 end if
156
157 ! [3.3] Top boundaries:
158
159 if (ite==ide) then
160 i = ite
161
162 do j = js, je
163 varb = 3.0*xb%v(i,j,k)-4.0*xb%v(i-1,j,k)+xb%v(i-2,j,k)
164
165 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)* varb + &
166 xb%v(i,j,k) * xb%coefy(i,j)* (xb%v(i,j+1,k) - xb%v(i,j-1,k))
167 end do
168 end if
169
170 ! [3.4] Left boundaries:
171
172 if (jts==jds) then
173 j = jts
174
175 do i = is, ie
176 varb = -3.0*xb%v(i,j,k)+4.0*xb%v(i,j+1,k)-xb%v(i,j+2,k)
177
178 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)* (xb%v(i+1,j,k) - xb%v(i-1,j,k)) + &
179 xb%v(i,j,k) * xb%coefy(i,j)* varb
180 end do
181 end if
182
183 ! [3.5] Right boundaries:
184
185 if (jte==jde) then
186 j = jte
187
188 do i = is, ie
189 varb = 3.0*xb%v(i,j,k)-4.0*xb%v(i,j+1,k)+xb%v(i,j+2,k)
190
191 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)* (xb%v(i+1,j,k) - xb%v(i-1,j,k)) + &
192 xb%v(i,j,k) * xb%coefy(i,j)* varb
193 end do
194 end if
195
196 ! [3.6] Corner points:
197
198 if (its==ids .AND. jts==jds) then
199 data(its,jts) = 0.5 * (data(its,jts+1) + data(its+1,jts))
200 end if
201
202 if (ite==ide .AND. jts==jds) then
203 data(ite,jts) = 0.5 * (data(ite-1,jts) + data(ite,jts+1))
204 end if
205
206 if (its==ids .AND. jte==jde) then
207 data(its,jde) = 0.5 * (data(its,jde-1) + data(its+1,jde))
208 end if
209
210 if (ite==ide .AND. jte==jde) then
211 data(ite,jte) = 0.5 * (data(ite-1,jte) + data(ite,jte-1))
212 end if
213 end if ! not global
214
215 ! [3.7] Multiply by rho and add to term_y
216
217 term_y(its:ite,jts:jte)=xb%rho(its:ite,jts:jte,k)* data(its:ite,jts:jte) + &
218 term_y(its:ite,jts:jte)
219
220 if (trace_use) call da_trace_exit("da_balance_cycloterm")
221
222 end subroutine da_balance_cycloterm
223
224