Line data Source code
1 : module sw_core
2 : !BOP
3 : !
4 : ! !MODULE: sw_core --- Utilities for solving the shallow-water equation
5 : !
6 : ! !USES:
7 : use dynamics_vars, only: T_FVDYCORE_GRID
8 : use shr_kind_mod, only : r8 => shr_kind_r8
9 :
10 : #ifdef NO_R16
11 : integer,parameter :: r16= selected_real_kind(12) ! 8 byte real
12 : #else
13 : integer,parameter :: r16= selected_real_kind(24) ! 16 byte real
14 : #endif
15 :
16 : !
17 : ! !PUBLIC MEMBER FUNCTIONS:
18 : public d2a2c_winds, c_sw, d_sw
19 : !
20 : ! !DESCRIPTION:
21 : !
22 : ! This module contains vertical independent part of the Lagrangian
23 : ! dynamics; in simpler terms, it solves the 2D shallow water equation
24 : ! (SWE).
25 : !
26 : ! \begin{tabular}{|l|l|} \hline \hline
27 : ! c_sw & \\ \hline
28 : ! d_sw &
29 : ! \end{tabular}
30 : !
31 : ! !REVISION HISTORY:
32 : ! 01.01.15 Lin Routines coalesced into this module
33 : ! 03.11.19 Sawyer Merged in CAM changes by Mirin
34 : ! 04.10.07 Sawyer ompinner now from dynamics_vars
35 : ! 05.03.25 Todling shr_kind_r8 can only be referenced once (MIPSpro-7.4.2)
36 : ! 05.05.25 Sawyer Merged CAM and GEOS5 versions (mostly CAM)
37 : ! 05.07.26 Worley Changes for Cray X1
38 : ! 05.07.05 Sawyer Interfaces of c_sw and d_sw simplified with grid
39 : ! 05.10.12 Worley More changes for Cray X1(E), avoiding array segment copying
40 : ! 06.01.18 Putman Allowed Y-dir courant number and mass flux to accumulate
41 : ! at jlast+1
42 : ! 06.09.06 Sawyer Isolated magic numbers as F90 parameters
43 : !
44 : !EOP
45 :
46 : ! Magic numbers used in this module
47 :
48 : real(r8), parameter :: D0_0 = 0.0_r8
49 : real(r8), parameter :: D0_125 = 0.125_r8
50 : real(r8), parameter :: D0_25 = 0.25_r8
51 : real(r8), parameter :: D0_5 = 0.5_r8
52 : real(r8), parameter :: D1_0 = 1.0_r8
53 : real(r8), parameter :: D2_0 = 2.0_r8
54 : real(r8), parameter :: D1E30 = 1.0e30_r8
55 :
56 : contains
57 :
58 : !-----------------------------------------------------------------------
59 : !BOP
60 : ! !IROUTINE: c_sw --- Solve the SWE on a C grid
61 : !
62 : ! !INTERFACE:
63 1505280 : subroutine c_sw(grid, u, v, pt, delp, &
64 1505280 : u2, v2, &
65 1505280 : uc, vc, ptc, delpf, ptk, &
66 : tiny, iord, jord, am_geom_crrct)
67 :
68 : ! Routine for shallow water dynamics on the C-grid
69 :
70 : ! !USES:
71 :
72 : use tp_core
73 : use pft_module, only : pft2d
74 :
75 : implicit none
76 :
77 : ! !INPUT PARAMETERS:
78 : type (T_FVDYCORE_GRID), intent(in) :: grid
79 : integer, intent(in):: iord
80 : integer, intent(in):: jord
81 : logical, intent(in):: am_geom_crrct
82 :
83 : real(r8), intent(in):: u2(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
84 : real(r8), intent(in):: v2(grid%im,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d)
85 :
86 : ! Prognostic variables:
87 : real(r8), intent(in):: u(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_s)
88 : real(r8), intent(in):: v(grid%im,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d)
89 : real(r8), intent(in):: pt(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
90 : real(r8), intent(in):: delp(grid%im,grid%jfirst:grid%jlast)
91 : real(r8), intent(in):: delpf(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
92 :
93 : real(r8), intent(in):: tiny
94 :
95 : ! !INPUT/OUTPUT PARAMETERS:
96 : real(r8), intent(inout):: uc(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
97 : real(r8), intent(inout):: vc(grid%im,grid%jfirst-2:grid%jlast+2 )
98 :
99 : ! !OUTPUT PARAMETERS:
100 : real(r8), intent(out):: ptc(grid%im,grid%jfirst:grid%jlast)
101 : real(r8), intent(out):: ptk(grid%im,grid%jfirst:grid%jlast)
102 :
103 : ! !DESCRIPTION:
104 : !
105 : ! Routine for shallow water dynamics on the C-grid
106 : !
107 : ! !REVISION HISTORY:
108 : ! WS 2003.11.19 Merged in CAM changes by Mirin
109 : ! WS 2004.10.07 Added ProTeX documentation
110 : ! WS 2005.07.01 Simplified interface by passing grid
111 : !
112 : !EOP
113 : !-----------------------------------------------------------------------
114 : !BOC
115 :
116 :
117 : !--------------------------------------------------------------
118 : ! Local
119 : real(r8) :: zt_c
120 : real(r8) :: dydt
121 : real(r8) :: dtdy5
122 : real(r8) :: rcap
123 :
124 1505280 : real(r8), pointer:: sc(:)
125 1505280 : real(r8), pointer:: dc(:,:)
126 :
127 1505280 : real(r8), pointer:: cosp(:)
128 1505280 : real(r8), pointer:: acosp(:)
129 1505280 : real(r8), pointer:: cose(:)
130 :
131 1505280 : real(r8), pointer:: dxdt(:)
132 1505280 : real(r8), pointer:: dxe(:)
133 1505280 : real(r8), pointer:: rdxe(:)
134 1505280 : real(r8), pointer:: dtdx2(:)
135 1505280 : real(r8), pointer:: dtdx4(:)
136 1505280 : real(r8), pointer:: dtxe5(:)
137 1505280 : real(r8), pointer:: dycp(:)
138 1505280 : real(r8), pointer:: cye(:)
139 :
140 1505280 : real(r8), pointer:: fc(:)
141 :
142 1505280 : real(r8), pointer:: sinlon(:)
143 1505280 : real(r8), pointer:: coslon(:)
144 1505280 : real(r8), pointer:: sinl5(:)
145 1505280 : real(r8), pointer:: cosl5(:)
146 :
147 3010560 : real(r8) :: fx(grid%im,grid%jfirst:grid%jlast)
148 3010560 : real(r8) :: xfx(grid%im,grid%jfirst:grid%jlast)
149 3010560 : real(r8) :: tm2(grid%im,grid%jfirst:grid%jlast)
150 :
151 3010560 : real(r8) :: va(grid%im,grid%jfirst-1:grid%jlast)
152 :
153 3010560 : real(r8) :: wk4(grid%im+2,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d)
154 :
155 3010560 : real(r8) :: wk1(grid%im,grid%jfirst-1:grid%jlast+1)
156 :
157 3010560 : real(r8) :: cry(grid%im,grid%jfirst-1:grid%jlast+1)
158 3010560 : real(r8) :: fy(grid%im,grid%jfirst-1:grid%jlast+1)
159 :
160 3010560 : real(r8) :: ymass(grid%im,grid%jfirst: grid%jlast+1)
161 3010560 : real(r8) :: yfx(grid%im,grid%jfirst: grid%jlast+1)
162 :
163 3010560 : real(r8) :: crx(grid%im,grid%jfirst-grid%ng_c:grid%jlast+grid%ng_c)
164 3010560 : real(r8) :: vort_u(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
165 3010560 : real(r8) :: vort(grid%im,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d)
166 :
167 3010560 : real(r8) :: fxjv(grid%im,grid%jfirst-1:grid%jn2g0)
168 3010560 : real(r8) :: p1dv(grid%im,grid%jfirst-1:grid%jn2g0)
169 3010560 : real(r8) :: cx1v(grid%im,grid%jfirst-1:grid%jn2g0)
170 :
171 3010560 : real(r8) :: qtmp(-grid%im/3:grid%im+grid%im/3)
172 3010560 : real(r8) :: qtmpv(-grid%im/3:grid%im+grid%im/3, grid%jfirst-1:grid%jn2g0)
173 3010560 : real(r8) :: slope(-grid%im/3:grid%im+grid%im/3)
174 3010560 : real(r8) :: al(-grid%im/3:grid%im+grid%im/3)
175 3010560 : real(r8) :: ar(-grid%im/3:grid%im+grid%im/3)
176 3010560 : real(r8) :: a6(-grid%im/3:grid%im+grid%im/3)
177 :
178 : real(r8) :: p1ke, p2ke
179 :
180 3010560 : logical :: ffsl(grid%jm)
181 3010560 : logical :: sldv(grid%jfirst-1:grid%jn2g0)
182 :
183 : integer :: i, j, im2
184 : integer :: js2g1, js2gc1, jn2gc, jn1g1, js2g0, js2gc, jn1gc
185 : integer :: im, jm, jfirst, jlast, jn2g0, ng_s, ng_c, ng_d
186 :
187 :
188 :
189 : !
190 : ! For convenience
191 : !
192 :
193 1505280 : im = grid%im
194 1505280 : jm = grid%jm
195 1505280 : jfirst = grid%jfirst
196 1505280 : jlast = grid%jlast
197 :
198 1505280 : jn2g0 = grid%jn2g0
199 :
200 1505280 : ng_c = grid%ng_c
201 1505280 : ng_d = grid%ng_d
202 1505280 : ng_s = grid%ng_s
203 :
204 1505280 : rcap = grid%rcap
205 :
206 1505280 : zt_c = grid%zt_c
207 1505280 : dydt = grid%dydt
208 1505280 : dtdy5 = grid%dtdy5
209 :
210 1505280 : sc => grid%sc
211 1505280 : dc => grid%dc
212 :
213 1505280 : cosp => grid%cosp
214 1505280 : acosp => grid%acosp
215 1505280 : cose => grid%cose
216 :
217 1505280 : dxdt => grid%dxdt
218 1505280 : dxe => grid%dxe
219 1505280 : rdxe => grid%rdxe
220 1505280 : dtdx2 => grid%dtdx2
221 1505280 : dtdx4 => grid%dtdx4
222 1505280 : dtxe5 => grid%dtxe5
223 1505280 : dycp => grid%dycp
224 1505280 : cye => grid%cye
225 1505280 : fc => grid%fc
226 :
227 1505280 : sinlon => grid%sinlon
228 1505280 : coslon => grid%coslon
229 1505280 : sinl5 => grid%sinl5
230 1505280 : cosl5 => grid%cosl5
231 :
232 :
233 : ! Set loop limits
234 :
235 1505280 : im2 = im/2
236 :
237 1505280 : js2g0 = max(2,jfirst)
238 1505280 : js2gc = max(2,jfirst-ng_c) ! NG lats on S (starting at 2)
239 1505280 : jn1gc = min(jm,jlast+ng_c) ! ng_c lats on N (ending at jm)
240 1505280 : js2g1 = max(2,jfirst-1)
241 1505280 : jn1g1 = min(jm,jlast+1)
242 1505280 : jn2gc = min(jm-1,jlast+ng_c) ! NG latitudes on N (ending at jm-1)
243 1505280 : js2gc1 = max(2,jfirst-ng_c+1) ! NG-1 latitudes on S (starting at 2)
244 :
245 : ! KE at poles
246 1505280 : if ( jfirst-ng_d <= 1 ) then
247 47040 : p1ke = D0_125*(u2(1, 1)**2 + v2(1, 1)**2)
248 : endif
249 :
250 1505280 : if ( jlast+ng_d >= jm ) then
251 47040 : p2ke = D0_125*(u2(1,jm)**2 + v2(1,jm)**2)
252 : endif
253 :
254 1505280 : if ( jfirst /= 1 ) then
255 428228640 : do i=1,im
256 428228640 : cry(i,jfirst-1) = dtdy5*vc(i,jfirst-1)
257 : enddo
258 :
259 : endif
260 :
261 7479360 : do j=js2g0,jn1g1 ! ymass needed on NS
262 1728014400 : do i=1,im
263 1720535040 : cry(i,j) = dtdy5*vc(i,j)
264 1726509120 : ymass(i,j) = cry(i,j)*cose(j)
265 : enddo
266 : enddo
267 :
268 : ! New va definition
269 :
270 1505280 : if (am_geom_crrct) then
271 0 : do j=js2g1,jn2g0 ! va needed on S (for YCC, iv==1)
272 0 : do i=1,im
273 : ! weight by cos
274 0 : va(i,j) = (cry(i,j)*cose(j)+cry(i,j+1)*cose(j+1))/(cose(j)+cose(j+1))
275 : end do
276 : end do
277 :
278 : else
279 7455840 : do j=js2g1,jn2g0 ! va needed on S (for YCC, iv==1)
280 1721217120 : do i=1,im
281 1719711840 : va(i,j) = 0.5_r8*(cry(i,j)+cry(i,j+1))
282 : end do
283 : end do
284 :
285 : end if
286 :
287 : ! SJL: Check if FFSL integer fluxes need to be computed
288 11901120 : do j=js2gc,jn2gc ! ffsl needed on N*sg S*sg
289 3004397760 : do i=1,im
290 3004397760 : crx(i,j) = uc(i,j)*dtdx2(j)
291 : enddo
292 10395840 : ffsl(j) = .false.
293 11901120 : if( cosp(j) < zt_c ) then
294 2644763281 : do i=1,im
295 2644763281 : if( abs(crx(i,j)) > D1_0 ) then
296 26147 : ffsl(j) = .true.
297 : #if ( !defined UNICOSMP ) || ( !defined NEC_SX )
298 26147 : exit
299 : #endif
300 : endif
301 : enddo
302 : endif
303 : enddo
304 :
305 : ! 2D transport of polar filtered delp (for computing fluxes!)
306 : ! Update is done on the unfiltered delp
307 :
308 3010560 : call tp2c( ptk, va(1,jfirst), delpf(1,jfirst-ng_c), &
309 1505280 : crx(1,jfirst-ng_c), cry(1,jfirst), &
310 : im, jm, iord, jord, ng_c, xfx, &
311 : yfx, ffsl, rcap, acosp, &
312 1505280 : crx(1,jfirst), ymass, cosp, &
313 6021120 : 0, jfirst, jlast)
314 :
315 5974080 : do j=js2g0,jn2g0 ! xfx not ghosted
316 5974080 : if( ffsl(j) ) then
317 3554700 : do i=1,im
318 3554700 : xfx(i,j) = xfx(i,j)/sign(max(abs(crx(i,j)),tiny),crx(i,j))
319 : enddo
320 : endif
321 : enddo
322 :
323 : ! pt-advection using pre-computed mass fluxes
324 : ! use tm2 below as the storage for pt increment
325 : ! WS 99.09.20 : pt, crx need on N*ng S*ng, yfx on N
326 :
327 3010560 : call tp2c(tm2 ,va(1,jfirst), pt(1,jfirst-ng_c), &
328 1505280 : crx(1,jfirst-ng_c), cry(1,jfirst), &
329 : im, jm, iord, jord, ng_c, fx, &
330 : fy(1,jfirst), ffsl, rcap, acosp, &
331 4515840 : xfx, yfx, cosp, 1, jfirst, jlast)
332 :
333 : ! use wk4, crx as work arrays
334 1505280 : call pft2d(ptk(1,js2g0), sc, &
335 : dc, im, jn2g0-js2g0+1, &
336 3010560 : wk4, crx )
337 : call pft2d(tm2(1,js2g0), sc, &
338 : dc, im, jn2g0-js2g0+1, &
339 1505280 : wk4, crx )
340 :
341 6021120 : do j=jfirst,jlast
342 1306583040 : do i=1,im
343 1300561920 : ptk(i,j) = delp(i,j) + ptk(i,j)
344 1305077760 : ptc(i,j) = (pt(i,j)*delp(i,j) + tm2(i,j))/ptk(i,j)
345 : enddo
346 : enddo
347 :
348 : !------------------
349 : ! Momentum equation
350 : !------------------
351 :
352 3010560 : call ycc(im, jm, fy, vc(1,jfirst-2), va(1,jfirst-1), &
353 4515840 : va(1,jfirst-1), jord, 1, jfirst, jlast)
354 :
355 7455840 : do j=js2g1,jn2g0
356 :
357 1719711840 : do i=1,im
358 1719711840 : cx1v(i,j) = dtdx4(j)*u2(i,j)
359 : enddo
360 :
361 5950560 : sldv(j) = .false.
362 5950560 : if( cosp(j) < zt_c ) then
363 1518981917 : do i=1,im
364 1518981917 : if( abs(cx1v(i,j)) > D1_0 ) then
365 15313 : sldv(j) = .true.
366 : #if ( !defined UNICOSMP ) || ( !defined NEC_SX )
367 15313 : exit
368 : #endif
369 : endif
370 : enddo
371 : endif
372 :
373 5950560 : p1dv(im,j) = uc(1,j)
374 1715266560 : do i=1,im-1
375 1713761280 : p1dv(i,j) = uc(i+1,j)
376 : enddo
377 :
378 : enddo
379 :
380 : call xtpv(im, sldv, fxjv, p1dv, cx1v, iord, cx1v, &
381 : cosp, 0, slope, qtmpv, al, ar, a6, &
382 : jfirst, jlast, js2g1, jn2g0, jm, &
383 : jfirst-1, jn2g0, jfirst-1, jn2g0, &
384 : jfirst-1, jn2g0, jfirst-1, jn2g0, &
385 1505280 : jfirst-1, jn2g0, jfirst-1, jn2g0)
386 :
387 7455840 : do j=js2g1,jn2g0
388 1721217120 : do i=1,im
389 1719711840 : wk1(i,j) = dxdt(j)*fxjv(i,j) + dydt*fy(i,j)
390 : enddo
391 : enddo
392 :
393 1505280 : if ( jfirst == 1 ) then
394 6797280 : do i=1,im
395 6797280 : wk1(i,1) = p1ke
396 : enddo
397 : endif
398 :
399 1505280 : if ( jlast == jm ) then
400 6797280 : do i=1,im
401 6797280 : wk1(i,jm) = p2ke
402 : enddo
403 : endif
404 :
405 : ! crx redefined
406 10442880 : do j=js2gc1,jn1gc
407 8937600 : crx(1,j) = dtxe5(j)*u(im,j)
408 2575534080 : do i=2,im
409 2574028800 : crx(i,j) = dtxe5(j)*u(i-1,j)
410 : enddo
411 : enddo
412 :
413 1505280 : if ( jfirst /=1 ) then
414 428228640 : do i=1,im
415 428228640 : cry(i,jfirst-1) = dtdy5*v(i,jfirst-1)
416 : enddo
417 : endif
418 :
419 6021120 : do j=jfirst,jlast
420 1306583040 : do i=1,im
421 1300561920 : cry(i,j) = dtdy5*v(i,j)
422 1305077760 : ymass(i,j) = cry(i,j)*cosp(j) ! ymass actually unghosted
423 : enddo
424 : enddo
425 :
426 5997600 : do j=js2g0,jlast
427 1299785760 : do i=1,im
428 1298280480 : tm2(i,j) = D0_5*(cry(i,j)+cry(i,j-1)) ! cry ghosted on S
429 : enddo
430 : enddo
431 :
432 : ! Compute absolute vorticity on the C-grid.
433 :
434 1505280 : if ( jfirst-ng_d <= 1 ) then
435 13594560 : do i=1,im
436 13594560 : vort_u(i,1) = D0_0
437 : enddo
438 : endif
439 :
440 11901120 : do j=js2gc,jn2gc
441 3005903040 : do i=1,im
442 3004397760 : vort_u(i,j) = uc(i,j)*cosp(j)
443 : enddo
444 : enddo
445 :
446 1505280 : if ( jlast+ng_d >= jm ) then
447 13594560 : do i=1,im
448 13594560 : vort_u(i,jm) = D0_0
449 : enddo
450 : endif
451 :
452 10442880 : do j=js2gc1,jn1gc
453 : ! The computed absolute vorticity on C-Grid is assigned to vort
454 17875200 : vort(1,j) = fc(j) + (vort_u(1,j-1)-vort_u(1,j))*cye(j) + &
455 26812800 : (vc(1,j) - vc(im,j))*rdxe(j)
456 :
457 2575534080 : do i=2,im
458 2565091200 : vort(i,j) = fc(j) + (vort_u(i,j-1)-vort_u(i,j))*cye(j) + &
459 5139120000 : (vc(i,j) - vc(i-1,j))*rdxe(j)
460 : enddo
461 : enddo
462 :
463 10442880 : do j=js2gc1,jn1gc ! ffsl needed on N*ng S*(ng-1)
464 8937600 : ffsl(j) = .false.
465 10442880 : if( cose(j) < zt_c ) then
466 2290390493 : do i=1,im
467 2290390493 : if( abs(crx(i,j)) > D1_0 ) then
468 29607 : ffsl(j) = .true.
469 : #if ( !defined UNICOSMP ) || ( !defined NEC_SX )
470 29607 : exit
471 : #endif
472 : endif
473 : enddo
474 : endif
475 : enddo
476 :
477 3010560 : call tpcc( tm2, ymass, vort(1,jfirst-ng_d), crx(1,jfirst-ng_c), &
478 1505280 : cry(1,jfirst), im, jm, ng_c, ng_d, &
479 : iord, jord, fx, fy(1,jfirst), ffsl, cose, &
480 6021120 : jfirst, jlast, slope, qtmp, al, ar, a6 )
481 :
482 5974080 : do j=js2g0,jn2g0
483 4468800 : uc(1,j) = uc(1,j) + dtdx2(j)*(wk1(im,j)-wk1(1,j)) + dycp(j)*fy(1,j)
484 1288519680 : do i=2,im
485 1287014400 : uc(i,j) = uc(i,j) + dtdx2(j)*(wk1(i-1,j)-wk1(i,j)) + dycp(j)*fy(i,j)
486 : enddo
487 : enddo
488 5997600 : do j=js2g0,jlast
489 1293788160 : do i=1,im-1
490 1293788160 : vc(i,j) = vc(i,j) + dtdy5*(wk1(i,j-1)-wk1(i,j))-dxe(j)*fx(i+1,j)
491 : enddo
492 5997600 : vc(im,j) = vc(im,j) + dtdy5*(wk1(im,j-1)-wk1(im,j))-dxe(j)*fx(1,j)
493 : enddo
494 : !EOC
495 1505280 : end subroutine c_sw
496 : !--------------------------------------------------------------------------
497 :
498 :
499 :
500 : !-----------------------------------------------------------------------
501 : !BOP
502 : ! !IROUTINE: d_sw --- Solve the SWE on a D grid
503 : !
504 : ! !INTERFACE:
505 1505280 : subroutine d_sw( grid, u, v, uc, vc, &
506 1505280 : pt, delp, delpf, cx3, cy3, &
507 1505280 : mfx, mfy, cdx, cdy, &
508 1505280 : cdxde, cdxdp, cdyde, cdydp, & !ldel2 variables
509 1505280 : cdxdiv, cdydiv, cdx4, cdy4, cdtau4, &
510 : ldiv2, ldiv4, ldel2, &
511 : iord, jord, tiny, am_correction, &
512 1505280 : ddp, duc, vf)
513 : !--------------------------------------------------------------------------
514 : ! Routine for shallow water dynamics on the D-grid
515 :
516 : ! !USES:
517 :
518 : use tp_core
519 : use pft_module, only : pft2d
520 :
521 : implicit none
522 :
523 : ! !INPUT PARAMETERS:
524 : type (T_FVDYCORE_GRID), intent(in) :: grid
525 : integer, intent(in):: iord
526 : integer, intent(in):: jord
527 : logical, intent(in) :: ldiv2,ldiv4,ldel2 !damping options
528 :
529 : ! Prognostic variables:
530 : real(r8), intent(inout):: u(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_s)
531 : real(r8), intent(inout):: v(grid%im,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d)
532 : ! Delta pressure
533 : real(r8), intent(inout):: delp(grid%im,grid%jfirst:grid%jlast)
534 : ! Potential temperature
535 : real(r8), intent(inout):: pt(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
536 :
537 : real(r8), intent(inout):: delpf(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
538 :
539 : real(r8), intent(in):: cdx (grid%js2g0:grid%jn1g1)
540 : real(r8), intent(in):: cdy (grid%js2g0:grid%jn1g1)
541 : !
542 : ! variables for div4 and del2 damping
543 : !
544 : real(r8), intent(in):: cdx4 (grid%js2g0:grid%jn1g1)
545 : real(r8), intent(in):: cdy4 (grid%js2g0:grid%jn1g1)
546 : real(r8), intent(in):: cdtau4(grid%js2g0:grid%jn1g1)
547 : real(r8), intent(in):: cdxde (grid%js2g0:grid%jn1g1)
548 : real(r8), intent(in):: cdxdp (grid%js2g0:grid%jn1g1)
549 : real(r8), intent(in):: cdydp (grid%js2g0:grid%jn1g1)
550 : real(r8), intent(in):: cdyde (grid%js2g0:grid%jn1g1)
551 : real(r8), intent(in):: cdxdiv(grid%jm)
552 : real(r8), intent(in):: cdydiv(grid%jm)
553 :
554 : real(r8), intent(in):: tiny
555 :
556 : ! !INPUT/OUTPUT PARAMETERS:
557 : real(r8), intent(inout):: uc(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
558 : real(r8), intent(inout):: vc(grid%im,grid%jfirst-2 :grid%jlast+2 )
559 : real(r8), intent(inout):: cx3(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)! Accumulated Courant number in X
560 : real(r8), intent(inout):: cy3(grid%im,grid%jfirst:grid%jlast+1) ! Accumulated Courant number in Y
561 : real(r8), intent(inout):: mfx(grid%im,grid%jfirst:grid%jlast) ! Mass flux in X (unghosted)
562 : real(r8), intent(inout):: mfy(grid%im,grid%jfirst:grid%jlast+1) ! Mass flux in Y
563 :
564 : ! AM correction
565 : logical, intent(in) :: am_correction ! logical switch for correction (generate out args)
566 : real(r8), intent(out) :: ddp(grid%im,grid%jfirst:grid%jlast)
567 : real(r8), intent(out) :: duc(grid%im,grid%jfirst:grid%jlast)
568 : real(r8), intent(out) :: vf(grid%im,grid%jfirst-2:grid%jlast+2 )
569 :
570 : ! !DESCRIPTION:
571 : !
572 : ! Routine for shallow water dynamics on the D-grid
573 : !
574 : ! !REVISION HISTORY:
575 : ! WS 2003.11.19 Merged in CAM changes by Mirin
576 : ! WS 2004.10.07 Added ProTeX documentation
577 : ! WS 2005.07.05 Simplified interface using grid
578 : !
579 : !EOP
580 : !-----------------------------------------------------------------------
581 : !BOC
582 :
583 :
584 : ! Local
585 : integer :: im
586 : integer :: jm
587 : integer :: jfirst
588 : integer :: jlast
589 : integer :: js2g0
590 : integer :: jn1g1
591 : integer :: ng_d
592 : integer :: ng_s
593 : integer :: nq
594 :
595 : real(r8) :: zt_d
596 : real(r8) :: tdy5
597 : real(r8) :: rdy
598 : real(r8) :: dtdy
599 : real(r8) :: dtdy5
600 : real(r8) :: rcap
601 :
602 1505280 : real(r8), pointer:: sc(:)
603 1505280 : real(r8), pointer:: dc(:,:)
604 1505280 : real(r8), pointer:: se(:)
605 1505280 : real(r8), pointer:: de(:,:)
606 :
607 1505280 : real(r8), pointer :: cosp(:)
608 1505280 : real(r8), pointer :: acosp(:)
609 1505280 : real(r8), pointer :: cose(:)
610 :
611 1505280 : real(r8), pointer :: sinlon(:)
612 1505280 : real(r8), pointer :: coslon(:)
613 1505280 : real(r8), pointer :: sinl5(:)
614 1505280 : real(r8), pointer :: cosl5(:)
615 :
616 1505280 : real(r8), pointer :: dtdx(:)
617 1505280 : real(r8), pointer :: dtdxe(:)
618 1505280 : real(r8), pointer :: dx(:)
619 1505280 : real(r8), pointer :: rdx(:)
620 1505280 : real(r8), pointer :: cy(:)
621 1505280 : real(r8), pointer :: dyce(:)
622 1505280 : real(r8), pointer :: dtxe5(:)
623 1505280 : real(r8), pointer :: txe5(:)
624 :
625 1505280 : real(r8), pointer :: f0(:)
626 :
627 3010560 : real(r8) fx(grid%im,grid%jfirst:grid%jlast)
628 3010560 : real(r8) xfx(grid%im,grid%jfirst:grid%jlast)
629 :
630 : !for del2 damping
631 3010560 : real(r8) :: wku(grid%im,grid%jfirst-1:grid%jlast+1)
632 3010560 : real(r8) :: wkv(grid%im,grid%jfirst-1:grid%jlast+1)
633 :
634 : !for div4 damping
635 3010560 : real(r8) :: wkdiv4(grid%im+2,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_s)
636 3010560 : real(r8) :: wk2div4(grid%im+1,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_s)
637 :
638 3010560 : real(r8) wk1(grid%im,grid%jfirst-1:grid%jlast+1)
639 :
640 3010560 : real(r8) cry(grid%im,grid%jfirst-1:grid%jlast+1)
641 3010560 : real(r8) fy(grid%im,grid%jfirst-2:grid%jlast+2)!halo changed for div4
642 :
643 3010560 : real(r8) ymass(grid%im,grid%jfirst: grid%jlast+1)
644 3010560 : real(r8) yfx(grid%im,grid%jfirst: grid%jlast+1)
645 :
646 3010560 : real(r8) va(grid%im,grid%jfirst-1:grid%jlast)
647 3010560 : real(r8) ub(grid%im,grid%jfirst: grid%jlast+1)
648 :
649 : ! AM correction
650 3010560 : real(r8) :: dfx(grid%im,grid%jfirst:grid%jlast)
651 3010560 : real(r8) :: dfy(grid%im,grid%jfirst-2:grid%jlast+2)
652 3010560 : real(r8) :: dvdx(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
653 :
654 3010560 : real(r8) crx(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
655 : #if defined(FILTER_MASS_FLUXES)
656 : real(r8) u2(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
657 : real(r8) v2(grid%im+2,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d)
658 : #endif
659 :
660 3010560 : real(r8) fxjv(grid%im,grid%jfirst-1:grid%jn1g1)
661 3010560 : real(r8) qtmpv(-grid%im/3:grid%im+grid%im/3, grid%jfirst-1:grid%jn1g1)
662 3010560 : real(r8) slope(-grid%im/3:grid%im+grid%im/3)
663 3010560 : real(r8) al(-grid%im/3:grid%im+grid%im/3)
664 3010560 : real(r8) ar(-grid%im/3:grid%im+grid%im/3)
665 3010560 : real(r8) a6(-grid%im/3:grid%im+grid%im/3)
666 :
667 : real(r8) c1, c2
668 3010560 : real(r8) uanp(grid%im), uasp(grid%im), vanp(grid%im), vasp(grid%im)
669 : real(r8) un, vn, us, vs, r2im
670 :
671 3010560 : real(r8) div (grid%im,grid%jfirst-1:grid%jlast+2) !for div4 damping
672 : real(r8) div4(grid%im,grid%jfirst-1:grid%jlast+1) !for div4 damping
673 :
674 3010560 : logical ffsl(grid%jm)
675 3010560 : logical sldv(grid%jfirst-1:grid%jn1g1)
676 :
677 : real(r8):: deldiv !for div4
678 :
679 :
680 : integer i, j
681 : integer js2gd, jn2g0, jn2g1, jn2gd, jn1gd
682 : integer jn2g2 !for extra halo for div4
683 : integer js2gs, jn1gs
684 : integer im2
685 :
686 : !
687 : ! For convenience
688 : !
689 1505280 : nq = grid%nq
690 :
691 1505280 : im = grid%im
692 1505280 : jm = grid%jm
693 1505280 : jfirst = grid%jfirst
694 1505280 : jlast = grid%jlast
695 1505280 : ng_d = grid%ng_d
696 1505280 : ng_s = grid%ng_s
697 1505280 : js2g0 = grid%js2g0
698 1505280 : jn1g1 = grid%jn1g1
699 :
700 1505280 : rcap = grid%rcap
701 1505280 : zt_d = grid%zt_d
702 1505280 : tdy5 = grid%tdy5
703 1505280 : rdy = grid%rdy
704 1505280 : dtdy = grid%dtdy
705 1505280 : dtdy5 = grid%dtdy5
706 :
707 1505280 : sc => grid%sc
708 1505280 : dc => grid%dc
709 1505280 : se => grid%se
710 1505280 : de => grid%de
711 :
712 1505280 : cosp => grid%cosp
713 1505280 : acosp => grid%acosp
714 1505280 : cose => grid%cose
715 :
716 1505280 : sinlon => grid%sinlon
717 1505280 : coslon => grid%coslon
718 1505280 : sinl5 => grid%sinl5
719 1505280 : cosl5 => grid%cosl5
720 :
721 1505280 : dtdx => grid%dtdx
722 1505280 : dtdxe => grid%dtdxe
723 1505280 : dx => grid%dx
724 1505280 : rdx => grid%rdx
725 1505280 : cy => grid%cy
726 1505280 : dyce => grid%dyce
727 1505280 : dtxe5 => grid%dtxe5
728 1505280 : txe5 => grid%txe5
729 :
730 1505280 : f0 => grid%f0
731 :
732 : ! Set loop limits
733 :
734 1505280 : jn2g0 = min(jm-1,jlast)
735 1505280 : jn2g1 = min(jm-1,jlast+1)
736 1505280 : jn2g2 = min(jm-1,jlast+2)
737 :
738 1505280 : js2gd = max(2,jfirst-ng_d) ! NG latitudes on S (starting at 1)
739 1505280 : jn2gd = min(jm-1,jlast+ng_d) ! NG latitudes on S (ending at jm-1)
740 1505280 : jn1gd = min(jm,jlast+ng_d) ! NG latitudes on N (ending at jm)
741 1505280 : js2gs = max(2,jfirst-ng_s)
742 1505280 : jn1gs = min(jm,jlast+ng_s) ! NG latitudes on N (ending at jm)
743 :
744 : ! Get C-grid U-wind at poles.
745 1505280 : im2 = im/2
746 1505280 : r2im = 0.5_r16/real(im,r16)
747 :
748 1505280 : if ( jfirst <= 1 ) then
749 : !
750 : ! Treat SP
751 : !
752 6773760 : do i=1,im-1
753 6750240 : uasp(i) = uc(i,2) + uc(i+1,2)
754 6773760 : vasp(i) = vc(i,2) + vc(i,3)
755 : enddo
756 23520 : uasp(im) = uc(im,2) + uc(1,2)
757 23520 : vasp(im) = vc(im,2) + vc(im,3)
758 :
759 : ! Projection at SP
760 :
761 23520 : us = D0_0
762 23520 : vs = D0_0
763 3410400 : do i=1,im2
764 6773760 : us = us + (uasp(i+im2)-uasp(i))*sinlon(i) &
765 10160640 : + (vasp(i)-vasp(i+im2))*coslon(i)
766 : vs = vs + (uasp(i+im2)-uasp(i))*coslon(i) &
767 3410400 : + (vasp(i+im2)-vasp(i))*sinlon(i)
768 : enddo
769 23520 : us = us*r2im
770 23520 : vs = vs*r2im
771 :
772 : ! get U-wind at SP
773 :
774 3410400 : do i=1,im2
775 3386880 : uc(i, 1) = -us*sinl5(i) - vs*cosl5(i)
776 3410400 : uc(i+im2, 1) = -uc(i, 1)
777 : enddo
778 :
779 : endif
780 :
781 1505280 : if ( jlast >= jm ) then
782 : !
783 : ! Treat NP
784 : !
785 6773760 : do i=1,im-1
786 6750240 : uanp(i) = uc(i,jm-1) + uc(i+1,jm-1)
787 6773760 : vanp(i) = vc(i,jm-1) + vc(i,jm)
788 : enddo
789 23520 : uanp(im) = uc(im,jm-1) + uc(1,jm-1)
790 23520 : vanp(im) = vc(im,jm-1) + vc(im,jm)
791 :
792 : ! Projection at NP
793 :
794 23520 : un = D0_0
795 23520 : vn = D0_0
796 3410400 : do i=1,im2
797 6773760 : un = un + (uanp(i+im2)-uanp(i))*sinlon(i) &
798 10160640 : + (vanp(i+im2)-vanp(i))*coslon(i)
799 : vn = vn + (uanp(i)-uanp(i+im2))*coslon(i) &
800 3410400 : + (vanp(i+im2)-vanp(i))*sinlon(i)
801 : enddo
802 23520 : un = un*r2im
803 23520 : vn = vn*r2im
804 :
805 : ! get U-wind at NP
806 :
807 3410400 : do i=1,im2
808 3386880 : uc(i,jm) = -un*sinl5(i) + vn*cosl5(i)
809 3410400 : uc(i+im2,jm) = -uc(i,jm)
810 : enddo
811 :
812 : endif
813 :
814 14817600 : do j=js2gd,jn2gd ! crx needed on N*ng S*ng
815 3848765760 : do i=1,im
816 3847260480 : crx(i,j) = dtdx(j)*uc(i,j)
817 : enddo
818 : enddo
819 :
820 14817600 : do j=js2gd,jn2gd ! ffsl needed on N*ng S*ng
821 13312320 : ffsl(j) = .false.
822 14817600 : if( cosp(j) < zt_d ) then
823 3379385079 : do i=1,im
824 3379385079 : if( abs(crx(i,j)) > D1_0 ) then
825 79787 : ffsl(j) = .true.
826 : #if ( !defined UNICOSMP ) || ( !defined NEC_SX )
827 79787 : exit
828 : #endif
829 : endif
830 : enddo
831 : endif
832 : enddo
833 :
834 7479360 : do j=js2g0,jn1g1 ! cry, ymass needed on N
835 1728014400 : do i=1,im
836 1720535040 : cry(i,j) = dtdy*vc(i,j)
837 1726509120 : ymass(i,j) = cry(i,j)*cose(j)
838 : enddo
839 : enddo
840 :
841 5974080 : do j=js2g0,jn2g0 ! No ghosting
842 1292988480 : do i=1,im
843 1291483200 : if( cry(i,j)*cry(i,j+1) > D0_0 ) then
844 1231265332 : if( cry(i,j) > D0_0 ) then
845 629235010 : va(i,j) = cry(i,j)
846 : else
847 602030322 : va(i,j) = cry(i,j+1) ! cry ghosted on N
848 : endif
849 : else
850 55749068 : va(i,j) = D0_0
851 : endif
852 : enddo
853 : enddo
854 :
855 : ! transport polar filtered delp
856 3010560 : call tp2c(ub(1,jfirst), va(1,jfirst), delpf(1,jfirst-ng_d), &
857 : crx(1,jfirst-ng_d),cry(1,jfirst),im,jm,iord,jord, &
858 : ng_d, xfx, yfx, ffsl, &
859 1505280 : rcap, acosp,crx(1,jfirst), ymass, &
860 6021120 : cosp, 0, jfirst, jlast)
861 :
862 : #if defined(FILTER_MASS_FLUXES)
863 : call pft2d( xfx(1,js2g0), sc, dc, im, jn2g0-js2g0+1, &
864 : v2, u2 )
865 : call pft2d( yfx(1,js2g0), se, de, im, jn1g1-js2g0+1, &
866 : v2, u2 )
867 : do j=js2g0,jn2g0
868 : do i=1,im-1
869 : ub(i,j) = xfx(i,j) - xfx(i+1,j) + (yfx(i,j)-yfx(i,j+1))*acosp(j)
870 : enddo
871 : ub(im,j) = xfx(im,j) - xfx(1,j) + (yfx(im,j)-yfx(im,j+1))*acosp(j)
872 : enddo
873 : #endif
874 :
875 : ! AM correction
876 6021120 : do j = jfirst, jlast
877 1306583040 : do i = 1, im
878 1300561920 : ddp(i,j) = 0.0_r8
879 1305077760 : duc(i,j) = 0.0_r8
880 : end do
881 : end do
882 :
883 1505280 : if (am_correction) then
884 0 : do j = js2g0, jn2g0
885 0 : do i = 1, im-1
886 0 : ddp( i,j) = (xfx(i+1,j) - xfx(i ,j))
887 : end do
888 0 : ddp(im,j) = (xfx( 1,j) - xfx(im,j))
889 : end do
890 : end if ! am_correction
891 :
892 : ! <<< Save necessary data for large time step tracer transport >>>
893 1505280 : if( nq > 0 ) then
894 5974080 : do j=js2g0,jn2g0 ! No ghosting needed
895 1292988480 : do i=1,im
896 1287014400 : cx3(i,j) = cx3(i,j) + crx(i,j)
897 1291483200 : mfx(i,j) = mfx(i,j) + xfx(i,j)
898 : enddo
899 : enddo
900 :
901 5997600 : do j=js2g0,jlast ! No ghosting needed
902 1299785760 : do i=1,im
903 1293788160 : cy3(i,j) = cy3(i,j) + cry(i,j)
904 1298280480 : mfy(i,j) = mfy(i,j) + yfx(i,j)
905 : enddo
906 : enddo
907 : endif
908 5974080 : do j=js2g0,jn2g0 ! No ghosting needed
909 5974080 : if( ffsl(j) ) then
910 9335567 : do i=1,im
911 9335567 : xfx(i,j) = xfx(i,j)/sign(max(abs(crx(i,j)),tiny),crx(i,j))
912 : enddo
913 : endif
914 : enddo
915 :
916 : ! Update delp
917 6021120 : do j=jfirst,jlast
918 1306583040 : do i=1,im
919 : ! SAVE old delp: pressure thickness ~ "air density"
920 1300561920 : wk1(i,j) = delp(i,j)
921 1305077760 : delp(i,j) = wk1(i,j) + ub(i,j)
922 : enddo
923 : enddo
924 :
925 : ! pt Advection
926 3010560 : call tp2c(ub(1,jfirst),va(1,jfirst),pt(1,jfirst-ng_d), &
927 : crx(1,jfirst-ng_d),cry(1,jfirst), &
928 1505280 : im,jm,iord,jord,ng_d,fx,fy(1,jfirst), &
929 : ffsl, rcap, acosp, &
930 4515840 : xfx, yfx(1,jfirst), cosp, 1, jfirst,jlast)
931 :
932 : ! Update pt.
933 6021120 : do j=jfirst,jlast
934 1306583040 : do i=1,im
935 1305077760 : pt(i,j) = (pt(i,j)*wk1(i,j)+ub(i,j)) / delp(i,j)
936 : enddo
937 : enddo
938 :
939 : ! Compute upwind biased kinetic energy at the four cell corners
940 :
941 : ! Start using ub as v (CFL) on B-grid (cell corners)
942 : ! ub here is average over updated C-grid points (involving
943 : ! 6 D-grid points) instead of 2 non-updated D-grid points
944 7479360 : do j=js2g0,jn1g1 ! ub needed on N
945 5974080 : ub(1,j) = dtdy5*(vc(1,j) + vc(im,j))
946 1722040320 : do i=2,im
947 1720535040 : ub(i,j) = dtdy5*(vc(i,j) + vc(i-1,j))
948 : enddo
949 : enddo
950 :
951 3010560 : call ytp(im, jm, fy(1,jfirst), v(1,jfirst-ng_d), ub(1,jfirst), &
952 4515840 : ub(1,jfirst), ng_d, jord, 1, jfirst, jlast)
953 : ! End using ub as v (CFL) on B-grid
954 :
955 7479360 : do j=js2g0,jn1g1 ! ub needed on N
956 1728014400 : do i=1,im
957 1726509120 : ub(i,j) = dtxe5(j)*(uc(i,j) + uc(i,j-1))
958 : ! uc will be used as work array after this point
959 : enddo
960 : enddo
961 :
962 7479360 : do j=js2g0,jn1g1 ! wk1 needed on N
963 5974080 : sldv(j) = .false.
964 7479360 : if( cose(j) < zt_d ) then
965 1522603628 : do i=1,im
966 1522603628 : if( abs(ub(i,j)) > D1_0 ) then ! ub ghosted on N
967 54999 : sldv(j) = .true.
968 : #if ( !defined UNICOSMP ) || ( !defined NEC_SX )
969 54999 : exit
970 : #endif
971 : endif
972 : enddo
973 : endif
974 : enddo
975 :
976 : call xtpv(im, sldv, fxjv, u, ub, iord, ub, cose, &
977 : 0, slope, qtmpv, al, ar, a6, &
978 : jfirst, jlast, js2g0, jn1g1, jm, &
979 : jfirst-1, jn1g1, jfirst-1, jn1g1, &
980 : jfirst-ng_d, jlast+ng_s, jfirst, jlast+1, &
981 1505280 : jfirst, jlast+1, jfirst-1, jn1g1)
982 7479360 : do j=js2g0,jn1g1 ! wk1 needed on N
983 1728014400 : do i=1,im
984 1726509120 : wk1(i,j) = txe5(j)*fxjv(i,j) + tdy5*fy(i,j) ! fy ghosted on N
985 : enddo
986 : enddo
987 :
988 : ! Add divergence damping to vector invariant form of the momentum eqn
989 : ! (absolute vorticity is damped by ffsl scheme, therefore divergence damping
990 : ! provides more consistent dissipation to divergent part of the flow)
991 :
992 : !--------------------------
993 : ! Perform divergence damping
994 : !--------------------------
995 :
996 1505280 : if (ldiv2) then
997 : !
998 : ! standard div2 damping
999 : !
1000 0 : do j=max(2,jfirst-1), jn2g1 ! fy need on NS (below)
1001 0 : do i=1,im
1002 : !
1003 : ! cosp is cosine(theta) at cell center discretized from the identity
1004 : !
1005 : ! cos(theta) = d(sin(theta))/d(theta)
1006 : !
1007 : ! as
1008 : !
1009 : ! cosp = (sine(j+1)-sine(j))/dp where dp = pi/(jm-1)
1010 : !
1011 0 : fy(i,j) = v(i,j)*cosp(j) ! v ghosted on NS at least
1012 : enddo
1013 : enddo
1014 :
1015 0 : do j=js2g0,jn1g1
1016 : ! i=1
1017 0 : uc(1,j) = u(im,j) - u(1,j) ! u ghosted on N at least
1018 0 : do i=2,im
1019 0 : uc(i,j) = u(i-1,j) - u(i,j)
1020 : enddo
1021 : enddo
1022 :
1023 0 : if ( jfirst == 1 ) then
1024 : ! j=2
1025 0 : do i=1,im
1026 0 : wk1(i,2) = wk1(i,2) - cdy(2)*fy(i, 2) + cdx(2)*uc(i,2)
1027 : enddo
1028 : endif
1029 :
1030 0 : do j=max(3,jfirst),jn2g1 ! wk1 needed on N (after TP2D)
1031 0 : do i=1,im
1032 0 : wk1(i,j) = wk1(i,j) + cdy(j)*(fy(i,j-1) - fy(i,j)) &
1033 0 : + cdx(j)*uc(i,j)
1034 : enddo
1035 : enddo
1036 :
1037 0 : if ( jlast == jm ) then
1038 0 : do i=1,im
1039 0 : wk1(i,jm) = wk1(i,jm) + cdy(jm)*fy(i,jm-1) + cdx(jm)*uc(i,jm)
1040 : enddo
1041 : endif
1042 : end if
1043 :
1044 1505280 : if (ldiv4) then
1045 : !
1046 : ! filter velocity components for stability
1047 : !
1048 1505280 : call pft2d(u(1,js2gd), grid%sediv4, grid%dediv4, im, jn1gs-js2gd+1, &
1049 3010560 : wkdiv4, wk2div4 )
1050 :
1051 1505280 : call pft2d(v(1,js2gs), grid%scdiv4, grid%dcdiv4, im, jn2gd-js2gs+1, &
1052 3010560 : wkdiv4, wk2div4 )
1053 :
1054 : !**************************************************************************
1055 : !
1056 : ! div4 damping
1057 : !
1058 : !**************************************************************************
1059 :
1060 11901120 : do j=max(2,jfirst-2), min(jm-1,grid%jlast+2) ! fy need on NS (below)
1061 3005903040 : do i=1,im
1062 3004397760 : fy(i,j) = v(i,j)*cosp(j) ! v ghosted on NS at least
1063 : enddo
1064 : enddo
1065 :
1066 10442880 : do j=max(2,jfirst-1),min(jm,grid%jlast+2)
1067 : ! i=1
1068 8937600 : uc(1,j) = u(im,j) - u(1,j) ! u ghosted on N at least
1069 2575534080 : do i=2,im
1070 2574028800 : uc(i,j) = u(i-1,j) - u(i,j)
1071 : enddo
1072 : enddo
1073 : !
1074 : ! compute divergence
1075 : !
1076 1505280 : if ( jfirst == 1 ) then
1077 : ! j=2
1078 6797280 : do i=1,im
1079 6797280 : div(i,2) = - cdydiv(2)*fy(i, 2) + cdxdiv(2)*uc(i,2)
1080 : enddo
1081 : endif
1082 :
1083 10395840 : do j=max(3,jfirst-1),min(jm-1,grid%jlast+2) ! wk1 needed on N (after TP2D)
1084 2570877120 : do i=1,im
1085 2569371840 : div(i,j) = cdydiv(j)*(fy(i,j-1) - fy(i,j)) + cdxdiv(j)*uc(i,j)
1086 : enddo
1087 : enddo
1088 :
1089 1505280 : if ( jlast == jm ) then
1090 6797280 : do i=1,im
1091 6797280 : div(i,jm) = cdydiv(jm)*fy(i,jm-1) + cdxdiv(jm)*uc(i,jm)
1092 : enddo
1093 : endif
1094 :
1095 1505280 : if ( jfirst == 1 ) then
1096 23520 : i=1
1097 23520 : j=2
1098 47040 : deldiv = cdx4(j) * (div(i+1,j )-D2_0*div(i,j)+div(im ,j ))+&
1099 47040 : cdy4(j) * (cosp(j)*(div(i,j+1)-div(i,j)))
1100 23520 : wk1(i,j) = wk1(i,j) +cdtau4(j)*deldiv
1101 6750240 : do i=2,im-1
1102 20180160 : deldiv = cdx4(j) * (div(i+1,j )-D2_0*div(i,j)+div(i-1,j ))+&
1103 26906880 : cdy4(j) * (cosp(j )*(div(i ,j+1)-div(i ,j)))
1104 6750240 : wk1(i,j) = wk1(i,j) + cdtau4(j)*deldiv
1105 : enddo
1106 23520 : i=im
1107 23520 : deldiv = cdx4(j) * (div(1 ,j )-D2_0*div(i,j)+div(i-1,j ))+&
1108 47040 : cdy4(j) * (cosp(j )*(div(i,j+1)-div(i,j)))
1109 23520 : wk1(i,j) = wk1(i,j) + cdtau4(j)*deldiv
1110 : endif
1111 :
1112 7432320 : do j=max(3,jfirst),jn2g1 ! wk1 needed on N (after TP2D)
1113 5927040 : i=1
1114 11854080 : deldiv = cdx4(j) * (div(i+1,j )-D2_0*div(i,j)+div(im ,j ))+&
1115 : cdy4(j) * ( &
1116 5927040 : cosp(j )*(div(i ,j+1)-div(i,j ))-&
1117 23708160 : cosp(j-1)*(div(i ,j )-div(i,j-1)))
1118 5927040 : wk1(i,j) = wk1(i,j) +cdtau4(j)*deldiv
1119 1701060480 : do i=2,im-1
1120 5085400320 : deldiv = cdx4(j) * (div(i+1,j )-D2_0*div(i,j)+div(i-1,j ))+&
1121 : cdy4(j) * ( &
1122 0 : cosp(j )*(div(i ,j+1)-div(i ,j ))-&
1123 6780533760 : cosp(j-1)*(div(i ,j )-div(i ,j-1)))
1124 1701060480 : wk1(i,j) = wk1(i,j) + cdtau4(j)*deldiv
1125 : enddo
1126 5927040 : i=im
1127 5927040 : deldiv = cdx4(j) * (div(1 ,j )-D2_0*div(i,j)+div(i-1,j ))+&
1128 : cdy4(j) * ( &
1129 0 : cosp(j )*(div(i ,j+1)-div(i,j ))-&
1130 11854080 : cosp(j-1)*(div(i ,j )-div(i,j-1)))
1131 7432320 : wk1(i,j) = wk1(i,j) + cdtau4(j)*deldiv
1132 : enddo
1133 :
1134 1505280 : if ( jlast == jm ) then
1135 23520 : i=1
1136 23520 : j = jm
1137 47040 : deldiv = cdx4(j) * (div(i+1,j )-D2_0*div(i,j)+div(im,j ))+&
1138 70560 : cdy4(j) * (-cosp(j-1)*(div(i,j)-div(i,j-1)))
1139 23520 : wk1(i,j) = wk1(i,j) +cdtau4(j)*deldiv
1140 6750240 : do i=2,im-1
1141 20180160 : deldiv = cdx4(j) * (div(i+1,j )-D2_0*div(i,j)+div(i-1,j ))+&
1142 26906880 : cdy4(j) * (-cosp(j-1)*(div(i,j)-div(i,j-1)))
1143 6750240 : wk1(i,j) = wk1(i,j) + cdtau4(j)*deldiv
1144 : enddo
1145 23520 : i=im
1146 23520 : j=jm
1147 23520 : deldiv = cdx4(j) * (div(1,j )-D2_0*div(i,j)+div(i-1,j ))+&
1148 47040 : cdy4(j) * (-cosp(j-1)*(div(i,j)-div(i,j-1)))
1149 23520 : wk1(i,j) = wk1(i,j) +cdtau4(j)*deldiv
1150 : endif
1151 : end if
1152 :
1153 2176634880 : wku(:,:) = D0_0
1154 2176634880 : wkv(:,:) = D0_0
1155 1505280 : if (ldel2) then
1156 : !**************************************************************************
1157 : !
1158 : ! regular del2 (Laplacian) damping
1159 : !
1160 : !**************************************************************************
1161 0 : if ( jfirst == 1 ) then
1162 0 : i=1
1163 0 : j=2
1164 0 : wku(i,j) = cdxde(j)* (u(i+1,j )-D2_0*u(i,j)+u(im ,j ))+&
1165 0 : cdyde(j)* (cosp(j )*(u(i ,j+1)-u(i ,j)))
1166 0 : wkv(i,j) = cdxdp(j) * (v(i+1,j )-D2_0*v(i,j)+v(im,j ))+&
1167 : cdydp(j) * ( &
1168 0 : cose(j+1)*(v(i ,j+1)-v(i,j ))-&
1169 0 : cose(j )*(v(i ,j ) ))
1170 : !line above: there is no v(i,j-1) since it is on the pole
1171 0 : do i=2,im-1
1172 0 : wku(i,j) = cdxde(j) * (u(i+1,j )-D2_0*u(i,j)+u(i-1,j ))+&
1173 0 : cdyde(j) * (cosp(j )*(u(i ,j+1)-u(i ,j)))
1174 : wkv(i,j) = cdxdp(j) * (v(i+1,j )-D2_0*v(i,j)+v(i-1,j ))+&
1175 : cdydp(j) * ( &
1176 0 : cose(j+1)*(v(i ,j+1)-v(i,j ))-&
1177 0 : cose(j )*(v(i ,j ) ))
1178 : enddo
1179 0 : i=im
1180 0 : wku(i,j) = cdxde(j) * (u(1 ,j )-D2_0*u(i,j)+u(i-1,j ))+&
1181 0 : cdyde(j) * (cosp(j )*(u(i ,j+1)-u(i ,j)))
1182 : wkv(i,j) = cdxdp(j) * (v(1,j )-D2_0*v(i,j)+v(i-1 ,j ))+&
1183 : cdydp(j) * ( &
1184 0 : cose(j+1)*(v(i ,j+1)-v(i,j ))-&
1185 0 : cose(j )*(v(i ,j ) ))
1186 : endif
1187 :
1188 0 : do j=max(3,jfirst),jn2g1 ! wk1 needed on N (after TP2D)
1189 0 : i=1
1190 0 : wku(i,j) = cdxde(j) * (u(i+1,j )-D2_0*u(i,j)+u(im ,j ))+&
1191 : cdyde(j) * ( &
1192 0 : cosp(j )*(u(i ,j+1)-u(i,j ))-&
1193 0 : cosp(j-1)*(u(i ,j )-u(i,j-1)))
1194 0 : wkv(i,j) = cdxdp(j) * (v(i+1,j )-D2_0*v(i,j)+v(im ,j ))+&
1195 : cdydp(j) * ( &
1196 0 : cose(j+1)*(v(i ,j+1)-v(i,j ))-&
1197 0 : cose(j )*(v(i ,j )-v(i,j-1)))
1198 0 : do i=2,im-1
1199 0 : wku(i,j) = cdxde(j) * (u(i+1,j )-D2_0*u(i,j)+u(i-1,j ))+&
1200 : cdyde(j) * ( &
1201 0 : cosp(j )*(u(i ,j+1)-u(i ,j ))-&
1202 0 : cosp(j-1)*(u(i ,j )-u(i ,j-1)))
1203 : wkv(i,j) = cdxdp(j) * (v(i+1,j )-D2_0*v(i,j)+v(i-1,j ))+&
1204 : cdydp(j) * ( &
1205 0 : cose(j+1)*(v(i ,j+1)-v(i ,j ))-&
1206 0 : cose(j )*(v(i ,j )-v(i ,j-1)))
1207 : enddo
1208 0 : i=im
1209 0 : wku(i,j) = cdxde(j) * (u(1 ,j )-D2_0*u(i,j)+u(i-1,j ))+&
1210 : cdyde(j) * ( &
1211 0 : cosp(j )*(u(i ,j+1)-u(i,j ))-&
1212 0 : cosp(j-1)*(u(i ,j )-u(i,j-1)))
1213 : wkv(i,j) = cdxdp(j) * (v(1 ,j )-D2_0*v(i,j)+v(i-1,j ))+&
1214 : cdydp(j) * ( &
1215 0 : cose(j+1)*(v(i ,j+1)-v(i,j ))-&
1216 0 : cose(j )*(v(i ,j )-v(i,j-1)))
1217 : enddo
1218 :
1219 0 : if ( jlast == jm ) then
1220 0 : i=1
1221 0 : j = jm
1222 0 : wku(i,jm) = cdxde(j) * (u(i+1,j )-D2_0*u(i,j)+u(im,j ))+&
1223 0 : cdyde(j) * (-cosp(j-1)*(u(i,j)-u(i,j-1)))
1224 0 : do i=2,im-1
1225 0 : wku(i,jm) = cdxde(j) * (u(i+1,j)-D2_0*u(i,j)+u(i-1,j))+&
1226 0 : cdyde(j) * (-cosp(j-1)*(u(i,j)-u(i,j-1)))
1227 : enddo
1228 0 : i=im
1229 0 : j=jm
1230 0 : wku(i,jm) = cdxde(j) * (u(1,j)-D2_0*u(i,j)+u(i-1,j))+&
1231 0 : cdyde(j) * (-cosp(j-1)*(u(i,j)-u(i,j-1)))
1232 : endif
1233 : end if
1234 :
1235 : !------------------------------------
1236 : ! End divergence damping computation
1237 : !------------------------------------
1238 :
1239 :
1240 : ! Compute Vorticity on the D grid
1241 : ! delpf used as work array
1242 :
1243 14864640 : do j=js2gd,jn1gd
1244 3862360320 : do i=1,im
1245 3860855040 : delpf(i,j) = u(i,j)*cose(j) ! u ghosted on N*ng S*ng
1246 : enddo
1247 : enddo
1248 :
1249 1505280 : if ( jfirst-ng_d <= 1 ) then
1250 47040 : c1 = D0_0
1251 13594560 : do i=1,im
1252 13594560 : c1 = c1 + delpf(i,2)
1253 : end do
1254 47040 : c1 = -c1*rdy*rcap
1255 13594560 : do i=1,im
1256 13594560 : uc(i,1) = c1
1257 : enddo
1258 : endif
1259 :
1260 1505280 : if ( jlast+ng_d >= jm ) then
1261 47040 : c2 = D0_0
1262 13594560 : do i=1,im
1263 13594560 : c2 = c2 + delpf(i,jm)
1264 : end do
1265 47040 : c2 = c2*rdy*rcap
1266 13594560 : do i=1,im
1267 13594560 : uc(i,jm) = c2
1268 : enddo
1269 : else
1270 : ! This is an attempt to avoid ghosting u on N*(ng+1)
1271 421431360 : do i=1,im
1272 421431360 : uc(i,jn2gd) = D1E30
1273 : enddo
1274 : endif
1275 :
1276 13359360 : do j=js2gd, min(jm-1,jlast+ng_d-1)
1277 3413975040 : do i=1,im-1
1278 10206362880 : uc(i ,j) = ( delpf(i,j) - delpf(i ,j+1)) * cy(j) + &
1279 13620337920 : (v(i+1,j) - v(i ,j)) * rdx(j)
1280 : enddo
1281 35562240 : uc(im,j) = (delpf(im,j) - delpf(im,j+1)) * cy(j) + &
1282 48921600 : (v( 1 ,j) - v(im,j)) * rdx(j)
1283 : enddo
1284 :
1285 : ! uc is relative vorticity at this point
1286 :
1287 14911680 : do j=max(1,jfirst-ng_d), jn1gd
1288 3875954880 : do i=1,im
1289 3874449600 : uc(i,j) = uc(i,j) + f0(j)
1290 : ! uc is absolute vorticity
1291 : enddo
1292 : enddo
1293 :
1294 3010560 : call tp2d(va(1,jfirst), uc(1,jfirst-ng_d), crx(1,jfirst-ng_d), &
1295 1505280 : cry(1,jfirst), im, jm, iord, jord, ng_d, fx, &
1296 1505280 : fy (1,jfirst), ffsl, crx(1,jfirst), &
1297 4515840 : ymass, cosp, 0, jfirst, jlast)
1298 :
1299 12042240 : do j = jfirst-2, jlast+2
1300 3046686720 : do i = 1, im
1301 3045181440 : vf (i,j) = 0.0_r8
1302 : end do
1303 : end do
1304 :
1305 : ! AM correction
1306 1505280 : if (am_correction) then
1307 :
1308 0 : if (jlast+ng_d >= jm) then
1309 0 : do i = 1, im
1310 0 : dvdx(i,jm) = 0.0_r8
1311 : end do
1312 : else
1313 0 : do i = 1, im
1314 0 : dvdx(i,jn2gd) = 0.0_r8
1315 : end do
1316 : end if
1317 :
1318 0 : if (jfirst-ng_d <= 1) then
1319 0 : do i = 1, im
1320 0 : dvdx(i,1) = 0.0_r8
1321 : end do
1322 : end if
1323 :
1324 0 : do j = js2gd, min(jm-1, jlast+ng_d-1)
1325 0 : do i = 1, im-1
1326 0 : dvdx( i,j) = (v(i+1,j) - v(i ,j)) * rdx(j)
1327 : end do
1328 0 : dvdx(im,j) = (v( 1,j) - v(im,j)) * rdx(j)
1329 : end do
1330 :
1331 0 : call tp2d(va(1,jfirst),dvdx(1,jfirst-ng_d), crx(1,jfirst-ng_d), &
1332 0 : cry(1,jfirst), im, jm, iord, jord, ng_d,dfx, &
1333 0 : dfy(1,jfirst), ffsl, crx(1,jfirst), &
1334 0 : ymass, cosp, 0, jfirst, jlast)
1335 :
1336 0 : do j = js2g0, jlast
1337 0 : do i = 1, im
1338 0 : duc(i,j) = dyce(j)*dfy(i,j)
1339 : end do
1340 : end do
1341 :
1342 0 : do j = js2g0, jlast
1343 0 : do i=1,im
1344 0 : vf(i,j) = dyce(j)*(fy(i,j)-dfy(i,j))
1345 : enddo
1346 : enddo
1347 :
1348 0 : do j = js2g0, jlast
1349 0 : do i = 1, im-1
1350 0 : duc( i,j) = dtdxe(j)*(wk1(i+1,j)-wk1(i ,j)) -duc( i,j) -wku( i,j)
1351 : end do
1352 0 : duc(im,j) = dtdxe(j)*(wk1( 1,j)-wk1(im,j)) -duc(im,j) -wku(im,j)
1353 : end do
1354 :
1355 : end if ! am_correction
1356 :
1357 5997600 : do j = js2g0, jlast
1358 1293788160 : do i=1,im-1
1359 1293788160 : uc(i ,j) = dtdxe(j)*(wk1(i ,j)-wk1(i+1,j)) +dyce(j)*fy(i ,j) +wku(i ,j)
1360 : enddo
1361 5997600 : uc(im,j) = dtdxe(j)*(wk1(im,j)-wk1( 1,j)) +dyce(j)*fy(im,j) +wku(im,j)
1362 : enddo
1363 :
1364 5974080 : do j = js2g0, jn2g0
1365 1292988480 : do i=1,im
1366 1291483200 : vc(i,j) = dtdy*(wk1(i,j)-wk1(i,j+1)) -dx(j)*fx(i,j) +wkv(i,j)
1367 : enddo
1368 : enddo
1369 :
1370 1505280 : end subroutine d_sw
1371 : !-----------------------------------------------------------------------
1372 :
1373 : !-----------------------------------------------------------------------
1374 : !BOP
1375 : ! !IROUTINE: d2a2c_winds --- Interpolate winds
1376 : !
1377 : ! !INTERFACE:
1378 1505280 : subroutine d2a2c_winds(grid, u, v, ua, va, uc, vc, u_cen, v_cen, &
1379 : reset_winds, met_rlx, am_geom_crrct)
1380 :
1381 : implicit none
1382 :
1383 : ! !PARAMETERS:
1384 : type (T_FVDYCORE_GRID), intent(in) :: grid
1385 :
1386 : real(r8), intent(in ):: u(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_s)
1387 : real(r8), intent(inout):: v(grid%im,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d)
1388 : real(r8), intent( out):: ua(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
1389 : real(r8), intent( out):: va(grid%im,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d)
1390 : real(r8), intent( out):: uc(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
1391 : real(r8), intent( out):: vc(grid%im,grid%jfirst-2:grid%jlast+2 )
1392 :
1393 : ! cell centered winds
1394 : logical , intent(in):: reset_winds
1395 : real(r8), intent(in):: u_cen(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
1396 : real(r8), intent(in):: v_cen(grid%im,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d)
1397 : real(r8), intent(in):: met_rlx
1398 : logical, intent(in):: am_geom_crrct
1399 :
1400 : ! !DESCRIPTION:
1401 : !
1402 : ! Calculate the cell-centered (A-grid) winds and the cell-wall perpendicular
1403 : ! (C-grid) winds from the cell-wall parallel (D-grid) winds.
1404 : !
1405 : ! This routine assumes that U and V have complete haloes! As a result,
1406 : ! the A-grid and C-grid results should have complete haloes from +/- ng_c
1407 : ! (which is generally smaller than ng_d). This feature has not been
1408 : ! thoroughly tested.
1409 : !
1410 : ! !REVISION HISTORY:
1411 : ! WP 2007.06.01 Creation
1412 : ! WS 2007.07.03 Added ProTeX documentation, removed unused vars.
1413 : ! WS 2009.04.15 Fixed numerous problems in indexing bounds
1414 : !
1415 : !EOP
1416 : !-----------------------------------------------------------------------
1417 : !BOC
1418 :
1419 : real(r8) us, vs, un, vn
1420 3010560 : real(r8) uanp(grid%im), uasp(grid%im), vanp(grid%im), vasp(grid%im), r2im
1421 :
1422 1505280 : real(r8), pointer:: sinlon(:)
1423 1505280 : real(r8), pointer:: coslon(:)
1424 1505280 : real(r8), pointer:: sinl5(:)
1425 1505280 : real(r8), pointer:: cosl5(:)
1426 :
1427 : ! AM correction
1428 1505280 : real(r8), pointer:: cosp(:)
1429 1505280 : real(r8), pointer:: acosp(:)
1430 1505280 : real(r8), pointer:: cose(:)
1431 :
1432 : integer :: i, j, im2
1433 : integer :: im, jm, jfirst, jlast, ng_s, ng_c, ng_d
1434 : integer :: jn1gc, js1gc, jn2gc, js2gc ! ng_c ghosted bounds
1435 : integer :: js2gd, jn2gd ! ng_d ghosted bounds
1436 : integer :: js2gs, jn2gsm1 ! ng_s ghosted bounds
1437 : integer :: js2g2, jn1g2 ! 2-lat ghosted bounds
1438 :
1439 1505280 : im = grid%im
1440 1505280 : jm = grid%jm
1441 1505280 : jfirst = grid%jfirst
1442 1505280 : jlast = grid%jlast
1443 :
1444 1505280 : ng_c = grid%ng_c
1445 1505280 : ng_d = grid%ng_d
1446 1505280 : ng_s = grid%ng_s
1447 :
1448 1505280 : im2 = im/2
1449 :
1450 1505280 : js1gc = max(1,jfirst-ng_c) ! ng_c lats on S (starting at 1)
1451 1505280 : jn1gc = min(jm,jlast+ng_c) ! ng_c latitudes on N (ending at jm)
1452 1505280 : js2gc = max(2,jfirst-ng_c) ! ng_c lats on S (starting at 2)
1453 1505280 : jn2gc = min(jm-1,jlast+ng_c) ! ng_c latitudes on N (ending at jm-1)
1454 1505280 : js2gs = max(2,jfirst-ng_s) ! ng_s latitudes on S (starting at 2)
1455 1505280 : jn2gsm1= min(jm-1,jlast+ng_s-1) ! ng_s-1 latitudes on N (ending at jm-1)
1456 1505280 : js2gd = max(2,jfirst-ng_d) ! ng_d latitudes on S (starting at 2)
1457 1505280 : jn2gd = min(jm-1,jlast+ng_d) ! ng_d latitudes on N (ending at jm-1)
1458 1505280 : js2g2 = max(2,jfirst-2) ! 2 latitudes on S (starting at 2)
1459 1505280 : jn1g2 = min(jm,jlast+2) ! 2 latitudes on N (ending at jm)
1460 :
1461 1505280 : sinlon => grid%sinlon
1462 1505280 : coslon => grid%coslon
1463 1505280 : sinl5 => grid%sinl5
1464 1505280 : cosl5 => grid%cosl5
1465 :
1466 : ! AM correction
1467 1505280 : cosp => grid%cosp
1468 1505280 : acosp => grid%acosp
1469 1505280 : cose => grid%cose
1470 :
1471 : ! Get D-grid V-wind at the poles.
1472 :
1473 1505280 : r2im = 0.5_r16/real(im,r16)
1474 :
1475 1505280 : if ( jfirst-ng_d <= 1 ) then
1476 :
1477 : !
1478 : ! Treat SP
1479 : !
1480 13547520 : do i=1,im-1
1481 13500480 : uasp(i) = u(i,2) + u(i,3)
1482 13547520 : vasp(i) = v(i,2) + v(i+1,2)
1483 : enddo
1484 :
1485 47040 : uasp(im) = u(im,2) + u(im,3)
1486 47040 : vasp(im) = v(im,2) + v(1,2)
1487 :
1488 : ! Projection at SP
1489 :
1490 47040 : us = D0_0
1491 47040 : vs = D0_0
1492 :
1493 6820800 : do i=1,im2
1494 13547520 : us = us + (uasp(i+im2)-uasp(i))*sinlon(i) &
1495 20321280 : + (vasp(i)-vasp(i+im2))*coslon(i)
1496 : vs = vs + (uasp(i+im2)-uasp(i))*coslon(i) &
1497 6820800 : + (vasp(i+im2)-vasp(i))*sinlon(i)
1498 : enddo
1499 47040 : us = us*r2im
1500 47040 : vs = vs*r2im
1501 :
1502 : ! get V-wind at SP
1503 :
1504 6820800 : do i=1,im2
1505 6773760 : v(i, 1) = us*cosl5(i) - vs*sinl5(i)
1506 6820800 : v(i+im2,1) = -v(i, 1)
1507 : enddo
1508 :
1509 : endif
1510 :
1511 1505280 : if ( jlast+ng_d >= jm ) then
1512 :
1513 : !
1514 : ! Treat NP
1515 : !
1516 13547520 : do i=1,im-1
1517 13500480 : uanp(i) = u(i,jm-1) + u(i,jm)
1518 13547520 : vanp(i) = v(i,jm-1) + v(i+1,jm-1)
1519 : enddo
1520 :
1521 47040 : uanp(im) = u(im,jm-1) + u(im,jm)
1522 47040 : vanp(im) = v(im,jm-1) + v(1,jm-1)
1523 :
1524 : ! Projection at NP
1525 :
1526 47040 : un = D0_0
1527 47040 : vn = D0_0
1528 6820800 : do i=1,im2
1529 13547520 : un = un + (uanp(i+im2)-uanp(i))*sinlon(i) &
1530 20321280 : + (vanp(i+im2)-vanp(i))*coslon(i)
1531 : vn = vn + (uanp(i)-uanp(i+im2))*coslon(i) &
1532 6820800 : + (vanp(i+im2)-vanp(i))*sinlon(i)
1533 : enddo
1534 47040 : un = un*r2im
1535 47040 : vn = vn*r2im
1536 :
1537 : ! get V-wind at NP
1538 :
1539 6820800 : do i=1,im2
1540 6773760 : v(i,jm) = -un*cosl5(i) - vn*sinl5(i)
1541 6820800 : v(i+im2,jm) = -v(i,jm)
1542 : enddo
1543 :
1544 : endif
1545 :
1546 3916738560 : ua(:,:) = D0_0
1547 3916738560 : va(:,:) = D0_0
1548 :
1549 14817600 : do j=js2gs,jn2gd
1550 3833948160 : do i=1,im-1
1551 3833948160 : va(i,j) = v(i,j) + v(i+1,j)
1552 : enddo
1553 14817600 : va(im,j) = v(im,j) + v(1,j)
1554 : enddo
1555 :
1556 1505280 : if (am_geom_crrct) then
1557 0 : do j = js2gd, jn2gsm1
1558 0 : do i = 1, im
1559 0 : ua(i,j) =(u(i,j)*cose(j) + u(i,j+1)*cose(j+1))/cosp(j) ! curl free -> curl free
1560 : end do
1561 : end do
1562 : else
1563 13359360 : do j = js2gd, jn2gsm1 ! WS: should be safe since u defined to jn2gs
1564 3427334400 : do i = 1, im
1565 3425829120 : ua(i,j) = u(i,j) + u(i,j+1) ! solid body -> solid body
1566 : end do
1567 : end do
1568 : end if
1569 : !
1570 : ! reset cell center winds to the offline meteorlogy data
1571 : !
1572 :
1573 1505280 : if ( reset_winds .and. met_rlx > D0_0 ) then
1574 0 : ua(:,:) = (D1_0-met_rlx)*ua(:,:) + met_rlx*( D2_0*u_cen(:,:) )
1575 0 : va(:,:) = (D1_0-met_rlx)*va(:,:) + met_rlx*( D2_0*v_cen(:,:) )
1576 : endif
1577 :
1578 1505280 : if ( jfirst-ng_d <= 1 ) then
1579 : ! Projection at SP
1580 : us = D0_0
1581 : vs = D0_0
1582 :
1583 :
1584 6820800 : do i=1,im2
1585 20321280 : us = us + (ua(i+im2,2)-ua(i ,2))*sinlon(i) &
1586 20321280 : + (va(i ,2)-va(i+im2,2))*coslon(i)
1587 : vs = vs + (ua(i+im2,2)-ua(i ,2))*coslon(i) &
1588 6820800 : + (va(i+im2,2)-va(i ,2))*sinlon(i)
1589 : enddo
1590 :
1591 47040 : us = us/im
1592 47040 : vs = vs/im
1593 :
1594 : ! SP
1595 6820800 : do i=1,im2
1596 6773760 : ua(i,1) = -us*sinlon(i) - vs*coslon(i)
1597 6773760 : va(i,1) = us*coslon(i) - vs*sinlon(i)
1598 6773760 : ua(i+im2,1) = -ua(i,1)
1599 6820800 : va(i+im2,1) = -va(i,1)
1600 : enddo
1601 :
1602 : endif
1603 :
1604 1505280 : if ( jlast+ng_d >= jm ) then
1605 :
1606 : ! Projection at NP
1607 : un = D0_0
1608 : vn = D0_0
1609 :
1610 6820800 : j = jm-1
1611 6820800 : do i=1,im2
1612 20321280 : un = un + (ua(i+im2,j)-ua(i ,j))*sinlon(i) &
1613 27095040 : + (va(i+im2,j)-va(i ,j))*coslon(i)
1614 : vn = vn + (ua(i ,j)-ua(i+im2,j))*coslon(i) &
1615 6820800 : + (va(i+im2,j)-va(i ,j))*sinlon(i)
1616 : enddo
1617 :
1618 47040 : un = un/im
1619 47040 : vn = vn/im
1620 :
1621 : ! NP
1622 6820800 : do i=1,im2
1623 6773760 : ua(i,jm) = -un*sinlon(i) + vn*coslon(i)
1624 6773760 : va(i,jm) = -un*coslon(i) - vn*sinlon(i)
1625 6773760 : ua(i+im2,jm) = -ua(i,jm)
1626 6820800 : va(i+im2,jm) = -va(i,jm)
1627 : enddo
1628 :
1629 : endif
1630 :
1631 :
1632 : ! A -> C
1633 11948160 : do j=js1gc,jn1gc ! uc needed N*ng_c S*ng_c, ua defined at poles
1634 10442880 : uc(1,j) = D0_25*(ua(1,j)+ua(im,j))
1635 3009054720 : do i=2,im
1636 3007549440 : uc(i,j) = D0_25*(ua(i,j)+ua(i-1,j))
1637 : enddo
1638 : enddo
1639 :
1640 1505280 : if (am_geom_crrct) then
1641 0 : do j = js2g2, jn1g2 ! vc needed N*2, S*2 (for ycc), va defined at poles
1642 0 : do i = 1, im
1643 0 : vc(i,j) = D0_25*(va(i,j)*cosp(j) + va(i,j-1)*cosp(j-1))/cose(j) ! div free -> div free
1644 : end do
1645 : end do
1646 : else
1647 11924640 : do j = js2g2, jn1g2
1648 3012700320 : do i = 1, im
1649 3011195040 : vc(i,j) = D0_25*(va(i,j) + va(i,j-1))
1650 : end do
1651 : end do
1652 : end if
1653 :
1654 1505280 : if ( jfirst==1 ) then
1655 6797280 : do i=1,im
1656 6797280 : vc(i,1) = D0_0 ! Needed to avoid undefined values in VC
1657 : enddo
1658 : endif
1659 : !EOC
1660 1505280 : end subroutine d2a2c_winds
1661 : !-----------------------------------------------------------------------
1662 : end module sw_core
|
