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 !---------------------------------------------------------------------------
25 ! [1.0] Initialise:
26 !---------------------------------------------------------------------------
27
28 ! Computation to check for edge of domain:
29 is = its; ie = ite; js = jts; je = jte
30 if (.not. global .and. its == ids ) is = ids+1
31 if (.not. global .and. ite == ide ) ie = ide-1
32 if (jts == jds ) js = jds+1
33 if (jte == jde ) je = jde-1
34
35 !---------------------------------------------------------------------------
36 ! [2.0] Calculate term_x = rho M ( u.du/dx + v.du/dy ):
37 !---------------------------------------------------------------------------
38
39 ! [2.1] Interior points:
40
41 do j = js, je
42 do i = is, ie
43 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j)*( xb%u(i+1,j,k) - &
44 xb%u(i-1,j,k) ) + xb%v(i,j,k) * xb%coefy(i,j)*( xb%u(i,j+1,k) - &
45 xb%u(i,j-1,k) )
46 end do
47 end do
48
49 if (.NOT. global) then ! For global only interior points
50
51 ! [2.2] Bottom boundaries:
52
53 if (its == ids ) then
54 i = its
55
56 do j = js, je
57 varb = -3.0*xb%u(i,j,k)+4.0*xb%u(i+1,j,k)-xb%u(i+2,j,k)
58
59 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * varb + &
60 xb%v(i,j,k) * xb%coefy(i,j) * ( xb%u(i,j+1,k) - xb%u(i,j-1,k) )
61 end do
62 end if
63
64 ! [2.3] Top boundaries:
65
66 if (ite == ide ) then
67 i = ite
68
69 do j = js, je
70 varb = 3.0*xb%u(i,j,k)-4.0*xb%u(i-1,j,k)+xb%u(i-2,j,k)
71
72 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * varb + &
73 xb%v(i,j,k) * xb%coefy(i,j) * ( xb%u(i,j+1,k) - xb%u(i,j-1,k) )
74 end do
75 end if
76
77 ! [2.4] Left boundaries:
78
79 if (jts == jds ) then
80 j = jts
81
82 do i = is, ie
83 varb = -3.0*xb%u(i,j,k)+4.0*xb%u(i,j+1,k)-xb%u(i,j+2,k)
84
85 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * ( xb%u(i+1,j,k) - &
86 xb%u(i-1,j,k) ) + xb%v(i,j,k) * xb%coefy(i,j) * varb
87 end do
88 end if
89
90 ! [2.5] Right boundaries:
91
92 if (jte == jde ) then
93 j = jte
94
95 do i = is, ie
96 varb = 3.0*xb%u(i,j,k)-4.0*xb%u(i,j-1,k)+xb%u(i,j-2,k)
97
98 data(i,j) = xb%u(i,j,k) * xb%coefx(i,j) * ( xb%u(i+1,j,k) - &
99 xb%u(i-1,j,k) ) + xb%v(i,j,k) * xb%coefy(i,j) * varb
100 end do
101 end if
102
103 ! [2.6] Corner points:
104
105 if (its == ids .AND. jts == jds ) then
106 data(its,jts) = 0.5 * ( data(its,jts+1) + data(its+1,jts) )
107 end if
108
109 if (ite == ide .AND. jts == jds ) then
110 data(ite,jts) = 0.5 * ( data(ite-1,jts) + data(ite,jts+1) )
111 end if
112
113 if (its == ids .AND. jte == jde ) then
114 data(its,jde) = 0.5 * ( data(its,jde-1) + data(its+1,jde) )
115 end if
116
117 if (ite == ide .AND. jte == jde ) then
118 data(ite,jte) = 0.5 * ( data(ite-1,jte) + data(ite,jte-1) )
119 end if
120
121 end if ! not global
122
123 ! [2.7] Multiply by rho and add to term_x:
124
125 term_x(its:ite,jts:jte) = xb%rho(its:ite,jts:jte,k)*data(its:ite,jts:jte) + &
126 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 end subroutine da_balance_cycloterm
221
222