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 602112 : subroutine c_sw(grid, u, v, pt, delp, &
64 602112 : u2, v2, &
65 602112 : 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 602112 : real(r8), pointer:: sc(:)
125 602112 : real(r8), pointer:: dc(:,:)
126 :
127 602112 : real(r8), pointer:: cosp(:)
128 602112 : real(r8), pointer:: acosp(:)
129 602112 : real(r8), pointer:: cose(:)
130 :
131 602112 : real(r8), pointer:: dxdt(:)
132 602112 : real(r8), pointer:: dxe(:)
133 602112 : real(r8), pointer:: rdxe(:)
134 602112 : real(r8), pointer:: dtdx2(:)
135 602112 : real(r8), pointer:: dtdx4(:)
136 602112 : real(r8), pointer:: dtxe5(:)
137 602112 : real(r8), pointer:: dycp(:)
138 602112 : real(r8), pointer:: cye(:)
139 :
140 602112 : real(r8), pointer:: fc(:)
141 :
142 602112 : real(r8), pointer:: sinlon(:)
143 602112 : real(r8), pointer:: coslon(:)
144 602112 : real(r8), pointer:: sinl5(:)
145 602112 : real(r8), pointer:: cosl5(:)
146 :
147 1204224 : real(r8) :: fx(grid%im,grid%jfirst:grid%jlast)
148 1204224 : real(r8) :: xfx(grid%im,grid%jfirst:grid%jlast)
149 1204224 : real(r8) :: tm2(grid%im,grid%jfirst:grid%jlast)
150 :
151 1204224 : real(r8) :: va(grid%im,grid%jfirst-1:grid%jlast)
152 :
153 1204224 : real(r8) :: wk4(grid%im+2,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d)
154 :
155 1204224 : real(r8) :: wk1(grid%im,grid%jfirst-1:grid%jlast+1)
156 :
157 1204224 : real(r8) :: cry(grid%im,grid%jfirst-1:grid%jlast+1)
158 1204224 : real(r8) :: fy(grid%im,grid%jfirst-1:grid%jlast+1)
159 :
160 1204224 : real(r8) :: ymass(grid%im,grid%jfirst: grid%jlast+1)
161 1204224 : real(r8) :: yfx(grid%im,grid%jfirst: grid%jlast+1)
162 :
163 1204224 : real(r8) :: crx(grid%im,grid%jfirst-grid%ng_c:grid%jlast+grid%ng_c)
164 1204224 : real(r8) :: vort_u(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
165 1204224 : real(r8) :: vort(grid%im,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_d)
166 :
167 1204224 : real(r8) :: fxjv(grid%im,grid%jfirst-1:grid%jn2g0)
168 1204224 : real(r8) :: p1dv(grid%im,grid%jfirst-1:grid%jn2g0)
169 1204224 : real(r8) :: cx1v(grid%im,grid%jfirst-1:grid%jn2g0)
170 :
171 1204224 : real(r8) :: qtmp(-grid%im/3:grid%im+grid%im/3)
172 1204224 : real(r8) :: qtmpv(-grid%im/3:grid%im+grid%im/3, grid%jfirst-1:grid%jn2g0)
173 1204224 : real(r8) :: slope(-grid%im/3:grid%im+grid%im/3)
174 1204224 : real(r8) :: al(-grid%im/3:grid%im+grid%im/3)
175 1204224 : real(r8) :: ar(-grid%im/3:grid%im+grid%im/3)
176 1204224 : real(r8) :: a6(-grid%im/3:grid%im+grid%im/3)
177 :
178 : real(r8) :: p1ke, p2ke
179 :
180 1204224 : logical :: ffsl(grid%jm)
181 1204224 : 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 602112 : im = grid%im
194 602112 : jm = grid%jm
195 602112 : jfirst = grid%jfirst
196 602112 : jlast = grid%jlast
197 :
198 602112 : jn2g0 = grid%jn2g0
199 :
200 602112 : ng_c = grid%ng_c
201 602112 : ng_d = grid%ng_d
202 602112 : ng_s = grid%ng_s
203 :
204 602112 : rcap = grid%rcap
205 :
206 602112 : zt_c = grid%zt_c
207 602112 : dydt = grid%dydt
208 602112 : dtdy5 = grid%dtdy5
209 :
210 602112 : sc => grid%sc
211 602112 : dc => grid%dc
212 :
213 602112 : cosp => grid%cosp
214 602112 : acosp => grid%acosp
215 602112 : cose => grid%cose
216 :
217 602112 : dxdt => grid%dxdt
218 602112 : dxe => grid%dxe
219 602112 : rdxe => grid%rdxe
220 602112 : dtdx2 => grid%dtdx2
221 602112 : dtdx4 => grid%dtdx4
222 602112 : dtxe5 => grid%dtxe5
223 602112 : dycp => grid%dycp
224 602112 : cye => grid%cye
225 602112 : fc => grid%fc
226 :
227 602112 : sinlon => grid%sinlon
228 602112 : coslon => grid%coslon
229 602112 : sinl5 => grid%sinl5
230 602112 : cosl5 => grid%cosl5
231 :
232 :
233 : ! Set loop limits
234 :
235 602112 : im2 = im/2
236 :
237 602112 : js2g0 = max(2,jfirst)
238 602112 : js2gc = max(2,jfirst-ng_c) ! NG lats on S (starting at 2)
239 602112 : jn1gc = min(jm,jlast+ng_c) ! ng_c lats on N (ending at jm)
240 602112 : js2g1 = max(2,jfirst-1)
241 602112 : jn1g1 = min(jm,jlast+1)
242 602112 : jn2gc = min(jm-1,jlast+ng_c) ! NG latitudes on N (ending at jm-1)
243 602112 : js2gc1 = max(2,jfirst-ng_c+1) ! NG-1 latitudes on S (starting at 2)
244 :
245 : ! KE at poles
246 602112 : if ( jfirst-ng_d <= 1 ) then
247 18816 : p1ke = D0_125*(u2(1, 1)**2 + v2(1, 1)**2)
248 : endif
249 :
250 602112 : if ( jlast+ng_d >= jm ) then
251 18816 : p2ke = D0_125*(u2(1,jm)**2 + v2(1,jm)**2)
252 : endif
253 :
254 602112 : if ( jfirst /= 1 ) then
255 171291456 : do i=1,im
256 171291456 : cry(i,jfirst-1) = dtdy5*vc(i,jfirst-1)
257 : enddo
258 :
259 : endif
260 :
261 2991744 : do j=js2g0,jn1g1 ! ymass needed on NS
262 691205760 : do i=1,im
263 688214016 : cry(i,j) = dtdy5*vc(i,j)
264 690603648 : ymass(i,j) = cry(i,j)*cose(j)
265 : enddo
266 : enddo
267 :
268 : ! New va definition
269 :
270 602112 : 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 2982336 : do j=js2g1,jn2g0 ! va needed on S (for YCC, iv==1)
280 688486848 : do i=1,im
281 687884736 : 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 4760448 : do j=js2gc,jn2gc ! ffsl needed on N*sg S*sg
289 1201759104 : do i=1,im
290 1201759104 : crx(i,j) = uc(i,j)*dtdx2(j)
291 : enddo
292 4158336 : ffsl(j) = .false.
293 4760448 : if( cosp(j) < zt_c ) then
294 171981934 : do i=1,im
295 171981934 : if( abs(crx(i,j)) > D1_0 ) then
296 8534 : ffsl(j) = .true.
297 : #if ( !defined UNICOSMP ) || ( !defined NEC_SX )
298 8534 : 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 1204224 : call tp2c( ptk, va(1,jfirst), delpf(1,jfirst-ng_c), &
309 602112 : crx(1,jfirst-ng_c), cry(1,jfirst), &
310 : im, jm, iord, jord, ng_c, xfx, &
311 : yfx, ffsl, rcap, acosp, &
312 602112 : crx(1,jfirst), ymass, cosp, &
313 2408448 : 0, jfirst, jlast)
314 :
315 2389632 : do j=js2g0,jn2g0 ! xfx not ghosted
316 2389632 : if( ffsl(j) ) then
317 1233163 : do i=1,im
318 1233163 : 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 1204224 : call tp2c(tm2 ,va(1,jfirst), pt(1,jfirst-ng_c), &
328 602112 : crx(1,jfirst-ng_c), cry(1,jfirst), &
329 : im, jm, iord, jord, ng_c, fx, &
330 : fy(1,jfirst), ffsl, rcap, acosp, &
331 1806336 : xfx, yfx, cosp, 1, jfirst, jlast)
332 :
333 : ! use wk4, crx as work arrays
334 602112 : call pft2d(ptk(1,js2g0), sc, &
335 : dc, im, jn2g0-js2g0+1, &
336 1204224 : wk4, crx )
337 : call pft2d(tm2(1,js2g0), sc, &
338 : dc, im, jn2g0-js2g0+1, &
339 602112 : wk4, crx )
340 :
341 2408448 : do j=jfirst,jlast
342 522633216 : do i=1,im
343 520224768 : ptk(i,j) = delp(i,j) + ptk(i,j)
344 522031104 : 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 1204224 : call ycc(im, jm, fy, vc(1,jfirst-2), va(1,jfirst-1), &
353 1806336 : va(1,jfirst-1), jord, 1, jfirst, jlast)
354 :
355 2982336 : do j=js2g1,jn2g0
356 :
357 687884736 : do i=1,im
358 687884736 : cx1v(i,j) = dtdx4(j)*u2(i,j)
359 : enddo
360 :
361 2380224 : sldv(j) = .false.
362 2380224 : if( cosp(j) < zt_c ) then
363 99587514 : do i=1,im
364 99587514 : if( abs(cx1v(i,j)) > D1_0 ) then
365 4267 : sldv(j) = .true.
366 : #if ( !defined UNICOSMP ) || ( !defined NEC_SX )
367 4267 : exit
368 : #endif
369 : endif
370 : enddo
371 : endif
372 :
373 2380224 : p1dv(im,j) = uc(1,j)
374 686106624 : do i=1,im-1
375 685504512 : 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 602112 : jfirst-1, jn2g0, jfirst-1, jn2g0)
386 :
387 2982336 : do j=js2g1,jn2g0
388 688486848 : do i=1,im
389 687884736 : wk1(i,j) = dxdt(j)*fxjv(i,j) + dydt*fy(i,j)
390 : enddo
391 : enddo
392 :
393 602112 : if ( jfirst == 1 ) then
394 2718912 : do i=1,im
395 2718912 : wk1(i,1) = p1ke
396 : enddo
397 : endif
398 :
399 602112 : if ( jlast == jm ) then
400 2718912 : do i=1,im
401 2718912 : wk1(i,jm) = p2ke
402 : enddo
403 : endif
404 :
405 : ! crx redefined
406 4177152 : do j=js2gc1,jn1gc
407 3575040 : crx(1,j) = dtxe5(j)*u(im,j)
408 1030213632 : do i=2,im
409 1029611520 : crx(i,j) = dtxe5(j)*u(i-1,j)
410 : enddo
411 : enddo
412 :
413 602112 : if ( jfirst /=1 ) then
414 171291456 : do i=1,im
415 171291456 : cry(i,jfirst-1) = dtdy5*v(i,jfirst-1)
416 : enddo
417 : endif
418 :
419 2408448 : do j=jfirst,jlast
420 522633216 : do i=1,im
421 520224768 : cry(i,j) = dtdy5*v(i,j)
422 522031104 : ymass(i,j) = cry(i,j)*cosp(j) ! ymass actually unghosted
423 : enddo
424 : enddo
425 :
426 2399040 : do j=js2g0,jlast
427 519914304 : do i=1,im
428 519312192 : 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 602112 : if ( jfirst-ng_d <= 1 ) then
435 5437824 : do i=1,im
436 5437824 : vort_u(i,1) = D0_0
437 : enddo
438 : endif
439 :
440 4760448 : do j=js2gc,jn2gc
441 1202361216 : do i=1,im
442 1201759104 : vort_u(i,j) = uc(i,j)*cosp(j)
443 : enddo
444 : enddo
445 :
446 602112 : if ( jlast+ng_d >= jm ) then
447 5437824 : do i=1,im
448 5437824 : vort_u(i,jm) = D0_0
449 : enddo
450 : endif
451 :
452 4177152 : do j=js2gc1,jn1gc
453 : ! The computed absolute vorticity on C-Grid is assigned to vort
454 7150080 : vort(1,j) = fc(j) + (vort_u(1,j-1)-vort_u(1,j))*cye(j) + &
455 10725120 : (vc(1,j) - vc(im,j))*rdxe(j)
456 :
457 1030213632 : do i=2,im
458 1026036480 : vort(i,j) = fc(j) + (vort_u(i,j-1)-vort_u(i,j))*cye(j) + &
459 2055648000 : (vc(i,j) - vc(i-1,j))*rdxe(j)
460 : enddo
461 : enddo
462 :
463 4177152 : do j=js2gc1,jn1gc ! ffsl needed on N*ng S*(ng-1)
464 3575040 : ffsl(j) = .false.
465 4177152 : if( cose(j) < zt_c ) then
466 155034110 : do i=1,im
467 155034110 : if( abs(crx(i,j)) > D1_0 ) then
468 10325 : ffsl(j) = .true.
469 : #if ( !defined UNICOSMP ) || ( !defined NEC_SX )
470 10325 : exit
471 : #endif
472 : endif
473 : enddo
474 : endif
475 : enddo
476 :
477 1204224 : call tpcc( tm2, ymass, vort(1,jfirst-ng_d), crx(1,jfirst-ng_c), &
478 602112 : cry(1,jfirst), im, jm, ng_c, ng_d, &
479 : iord, jord, fx, fy(1,jfirst), ffsl, cose, &
480 2408448 : jfirst, jlast, slope, qtmp, al, ar, a6 )
481 :
482 2389632 : do j=js2g0,jn2g0
483 1787520 : uc(1,j) = uc(1,j) + dtdx2(j)*(wk1(im,j)-wk1(1,j)) + dycp(j)*fy(1,j)
484 515407872 : do i=2,im
485 514805760 : 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 2399040 : do j=js2g0,jlast
489 517515264 : do i=1,im-1
490 517515264 : vc(i,j) = vc(i,j) + dtdy5*(wk1(i,j-1)-wk1(i,j))-dxe(j)*fx(i+1,j)
491 : enddo
492 2399040 : vc(im,j) = vc(im,j) + dtdy5*(wk1(im,j-1)-wk1(im,j))-dxe(j)*fx(1,j)
493 : enddo
494 : !EOC
495 602112 : 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 602112 : subroutine d_sw( grid, u, v, uc, vc, &
506 602112 : pt, delp, delpf, cx3, cy3, &
507 602112 : mfx, mfy, cdx, cdy, &
508 602112 : cdxde, cdxdp, cdyde, cdydp, & !ldel2 variables
509 602112 : cdxdiv, cdydiv, cdx4, cdy4, cdtau4, &
510 : ldiv2, ldiv4, ldel2, &
511 : iord, jord, tiny, am_correction, &
512 602112 : 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 602112 : real(r8), pointer:: sc(:)
603 602112 : real(r8), pointer:: dc(:,:)
604 602112 : real(r8), pointer:: se(:)
605 602112 : real(r8), pointer:: de(:,:)
606 :
607 602112 : real(r8), pointer :: cosp(:)
608 602112 : real(r8), pointer :: acosp(:)
609 602112 : real(r8), pointer :: cose(:)
610 :
611 602112 : real(r8), pointer :: sinlon(:)
612 602112 : real(r8), pointer :: coslon(:)
613 602112 : real(r8), pointer :: sinl5(:)
614 602112 : real(r8), pointer :: cosl5(:)
615 :
616 602112 : real(r8), pointer :: dtdx(:)
617 602112 : real(r8), pointer :: dtdxe(:)
618 602112 : real(r8), pointer :: dx(:)
619 602112 : real(r8), pointer :: rdx(:)
620 602112 : real(r8), pointer :: cy(:)
621 602112 : real(r8), pointer :: dyce(:)
622 602112 : real(r8), pointer :: dtxe5(:)
623 602112 : real(r8), pointer :: txe5(:)
624 :
625 602112 : real(r8), pointer :: f0(:)
626 :
627 1204224 : real(r8) fx(grid%im,grid%jfirst:grid%jlast)
628 1204224 : real(r8) xfx(grid%im,grid%jfirst:grid%jlast)
629 :
630 : !for del2 damping
631 1204224 : real(r8) :: wku(grid%im,grid%jfirst-1:grid%jlast+1)
632 1204224 : real(r8) :: wkv(grid%im,grid%jfirst-1:grid%jlast+1)
633 :
634 : !for div4 damping
635 1204224 : real(r8) :: wkdiv4(grid%im+2,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_s)
636 1204224 : real(r8) :: wk2div4(grid%im+1,grid%jfirst-grid%ng_s:grid%jlast+grid%ng_s)
637 :
638 1204224 : real(r8) wk1(grid%im,grid%jfirst-1:grid%jlast+1)
639 :
640 1204224 : real(r8) cry(grid%im,grid%jfirst-1:grid%jlast+1)
641 1204224 : real(r8) fy(grid%im,grid%jfirst-2:grid%jlast+2)!halo changed for div4
642 :
643 1204224 : real(r8) ymass(grid%im,grid%jfirst: grid%jlast+1)
644 1204224 : real(r8) yfx(grid%im,grid%jfirst: grid%jlast+1)
645 :
646 1204224 : real(r8) va(grid%im,grid%jfirst-1:grid%jlast)
647 1204224 : real(r8) ub(grid%im,grid%jfirst: grid%jlast+1)
648 :
649 : ! AM correction
650 1204224 : real(r8) :: dfx(grid%im,grid%jfirst:grid%jlast)
651 1204224 : real(r8) :: dfy(grid%im,grid%jfirst-2:grid%jlast+2)
652 1204224 : real(r8) :: dvdx(grid%im,grid%jfirst-grid%ng_d:grid%jlast+grid%ng_d)
653 :
654 1204224 : 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 1204224 : real(r8) fxjv(grid%im,grid%jfirst-1:grid%jn1g1)
661 1204224 : real(r8) qtmpv(-grid%im/3:grid%im+grid%im/3, grid%jfirst-1:grid%jn1g1)
662 1204224 : real(r8) slope(-grid%im/3:grid%im+grid%im/3)
663 1204224 : real(r8) al(-grid%im/3:grid%im+grid%im/3)
664 1204224 : real(r8) ar(-grid%im/3:grid%im+grid%im/3)
665 1204224 : real(r8) a6(-grid%im/3:grid%im+grid%im/3)
666 :
667 : real(r8) c1, c2
668 1204224 : real(r8) uanp(grid%im), uasp(grid%im), vanp(grid%im), vasp(grid%im)
669 : real(r8) un, vn, us, vs, r2im
670 :
671 1204224 : 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 1204224 : logical ffsl(grid%jm)
675 1204224 : 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 602112 : nq = grid%nq
690 :
691 602112 : im = grid%im
692 602112 : jm = grid%jm
693 602112 : jfirst = grid%jfirst
694 602112 : jlast = grid%jlast
695 602112 : ng_d = grid%ng_d
696 602112 : ng_s = grid%ng_s
697 602112 : js2g0 = grid%js2g0
698 602112 : jn1g1 = grid%jn1g1
699 :
700 602112 : rcap = grid%rcap
701 602112 : zt_d = grid%zt_d
702 602112 : tdy5 = grid%tdy5
703 602112 : rdy = grid%rdy
704 602112 : dtdy = grid%dtdy
705 602112 : dtdy5 = grid%dtdy5
706 :
707 602112 : sc => grid%sc
708 602112 : dc => grid%dc
709 602112 : se => grid%se
710 602112 : de => grid%de
711 :
712 602112 : cosp => grid%cosp
713 602112 : acosp => grid%acosp
714 602112 : cose => grid%cose
715 :
716 602112 : sinlon => grid%sinlon
717 602112 : coslon => grid%coslon
718 602112 : sinl5 => grid%sinl5
719 602112 : cosl5 => grid%cosl5
720 :
721 602112 : dtdx => grid%dtdx
722 602112 : dtdxe => grid%dtdxe
723 602112 : dx => grid%dx
724 602112 : rdx => grid%rdx
725 602112 : cy => grid%cy
726 602112 : dyce => grid%dyce
727 602112 : dtxe5 => grid%dtxe5
728 602112 : txe5 => grid%txe5
729 :
730 602112 : f0 => grid%f0
731 :
732 : ! Set loop limits
733 :
734 602112 : jn2g0 = min(jm-1,jlast)
735 602112 : jn2g1 = min(jm-1,jlast+1)
736 602112 : jn2g2 = min(jm-1,jlast+2)
737 :
738 602112 : js2gd = max(2,jfirst-ng_d) ! NG latitudes on S (starting at 1)
739 602112 : jn2gd = min(jm-1,jlast+ng_d) ! NG latitudes on S (ending at jm-1)
740 602112 : jn1gd = min(jm,jlast+ng_d) ! NG latitudes on N (ending at jm)
741 602112 : js2gs = max(2,jfirst-ng_s)
742 602112 : jn1gs = min(jm,jlast+ng_s) ! NG latitudes on N (ending at jm)
743 :
744 : ! Get C-grid U-wind at poles.
745 602112 : im2 = im/2
746 602112 : r2im = 0.5_r16/real(im,r16)
747 :
748 602112 : if ( jfirst <= 1 ) then
749 : !
750 : ! Treat SP
751 : !
752 2709504 : do i=1,im-1
753 2700096 : uasp(i) = uc(i,2) + uc(i+1,2)
754 2709504 : vasp(i) = vc(i,2) + vc(i,3)
755 : enddo
756 9408 : uasp(im) = uc(im,2) + uc(1,2)
757 9408 : vasp(im) = vc(im,2) + vc(im,3)
758 :
759 : ! Projection at SP
760 :
761 9408 : us = D0_0
762 9408 : vs = D0_0
763 1364160 : do i=1,im2
764 2709504 : us = us + (uasp(i+im2)-uasp(i))*sinlon(i) &
765 4064256 : + (vasp(i)-vasp(i+im2))*coslon(i)
766 : vs = vs + (uasp(i+im2)-uasp(i))*coslon(i) &
767 1364160 : + (vasp(i+im2)-vasp(i))*sinlon(i)
768 : enddo
769 9408 : us = us*r2im
770 9408 : vs = vs*r2im
771 :
772 : ! get U-wind at SP
773 :
774 1364160 : do i=1,im2
775 1354752 : uc(i, 1) = -us*sinl5(i) - vs*cosl5(i)
776 1364160 : uc(i+im2, 1) = -uc(i, 1)
777 : enddo
778 :
779 : endif
780 :
781 602112 : if ( jlast >= jm ) then
782 : !
783 : ! Treat NP
784 : !
785 2709504 : do i=1,im-1
786 2700096 : uanp(i) = uc(i,jm-1) + uc(i+1,jm-1)
787 2709504 : vanp(i) = vc(i,jm-1) + vc(i,jm)
788 : enddo
789 9408 : uanp(im) = uc(im,jm-1) + uc(1,jm-1)
790 9408 : vanp(im) = vc(im,jm-1) + vc(im,jm)
791 :
792 : ! Projection at NP
793 :
794 9408 : un = D0_0
795 9408 : vn = D0_0
796 1364160 : do i=1,im2
797 2709504 : un = un + (uanp(i+im2)-uanp(i))*sinlon(i) &
798 4064256 : + (vanp(i+im2)-vanp(i))*coslon(i)
799 : vn = vn + (uanp(i)-uanp(i+im2))*coslon(i) &
800 1364160 : + (vanp(i+im2)-vanp(i))*sinlon(i)
801 : enddo
802 9408 : un = un*r2im
803 9408 : vn = vn*r2im
804 :
805 : ! get U-wind at NP
806 :
807 1364160 : do i=1,im2
808 1354752 : uc(i,jm) = -un*sinl5(i) + vn*cosl5(i)
809 1364160 : uc(i+im2,jm) = -uc(i,jm)
810 : enddo
811 :
812 : endif
813 :
814 5927040 : do j=js2gd,jn2gd ! crx needed on N*ng S*ng
815 1539506304 : do i=1,im
816 1538904192 : crx(i,j) = dtdx(j)*uc(i,j)
817 : enddo
818 : enddo
819 :
820 5927040 : do j=js2gd,jn2gd ! ffsl needed on N*ng S*ng
821 5324928 : ffsl(j) = .false.
822 5927040 : if( cosp(j) < zt_d ) then
823 470136736 : do i=1,im
824 470136736 : if( abs(crx(i,j)) > D1_0 ) then
825 36763 : ffsl(j) = .true.
826 : #if ( !defined UNICOSMP ) || ( !defined NEC_SX )
827 36763 : exit
828 : #endif
829 : endif
830 : enddo
831 : endif
832 : enddo
833 :
834 2991744 : do j=js2g0,jn1g1 ! cry, ymass needed on N
835 691205760 : do i=1,im
836 688214016 : cry(i,j) = dtdy*vc(i,j)
837 690603648 : ymass(i,j) = cry(i,j)*cose(j)
838 : enddo
839 : enddo
840 :
841 2389632 : do j=js2g0,jn2g0 ! No ghosting
842 517195392 : do i=1,im
843 516593280 : if( cry(i,j)*cry(i,j+1) > D0_0 ) then
844 491226349 : if( cry(i,j) > D0_0 ) then
845 248110843 : va(i,j) = cry(i,j)
846 : else
847 243115506 : va(i,j) = cry(i,j+1) ! cry ghosted on N
848 : endif
849 : else
850 23579411 : va(i,j) = D0_0
851 : endif
852 : enddo
853 : enddo
854 :
855 : ! transport polar filtered delp
856 1204224 : 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 602112 : rcap, acosp,crx(1,jfirst), ymass, &
860 2408448 : 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 2408448 : do j = jfirst, jlast
877 522633216 : do i = 1, im
878 520224768 : ddp(i,j) = 0.0_r8
879 522031104 : duc(i,j) = 0.0_r8
880 : end do
881 : end do
882 :
883 602112 : 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 602112 : if( nq > 0 ) then
894 2389632 : do j=js2g0,jn2g0 ! No ghosting needed
895 517195392 : do i=1,im
896 514805760 : cx3(i,j) = cx3(i,j) + crx(i,j)
897 516593280 : mfx(i,j) = mfx(i,j) + xfx(i,j)
898 : enddo
899 : enddo
900 :
901 2399040 : do j=js2g0,jlast ! No ghosting needed
902 519914304 : do i=1,im
903 517515264 : cy3(i,j) = cy3(i,j) + cry(i,j)
904 519312192 : mfy(i,j) = mfy(i,j) + yfx(i,j)
905 : enddo
906 : enddo
907 : endif
908 2389632 : do j=js2g0,jn2g0 ! No ghosting needed
909 2389632 : if( ffsl(j) ) then
910 4499730 : do i=1,im
911 4499730 : 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 2408448 : do j=jfirst,jlast
918 522633216 : do i=1,im
919 : ! SAVE old delp: pressure thickness ~ "air density"
920 520224768 : wk1(i,j) = delp(i,j)
921 522031104 : delp(i,j) = wk1(i,j) + ub(i,j)
922 : enddo
923 : enddo
924 :
925 : ! pt Advection
926 1204224 : call tp2c(ub(1,jfirst),va(1,jfirst),pt(1,jfirst-ng_d), &
927 : crx(1,jfirst-ng_d),cry(1,jfirst), &
928 602112 : im,jm,iord,jord,ng_d,fx,fy(1,jfirst), &
929 : ffsl, rcap, acosp, &
930 1806336 : xfx, yfx(1,jfirst), cosp, 1, jfirst,jlast)
931 :
932 : ! Update pt.
933 2408448 : do j=jfirst,jlast
934 522633216 : do i=1,im
935 522031104 : 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 2991744 : do j=js2g0,jn1g1 ! ub needed on N
945 2389632 : ub(1,j) = dtdy5*(vc(1,j) + vc(im,j))
946 688816128 : do i=2,im
947 688214016 : ub(i,j) = dtdy5*(vc(i,j) + vc(i-1,j))
948 : enddo
949 : enddo
950 :
951 1204224 : call ytp(im, jm, fy(1,jfirst), v(1,jfirst-ng_d), ub(1,jfirst), &
952 1806336 : ub(1,jfirst), ng_d, jord, 1, jfirst, jlast)
953 : ! End using ub as v (CFL) on B-grid
954 :
955 2991744 : do j=js2g0,jn1g1 ! ub needed on N
956 691205760 : do i=1,im
957 690603648 : 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 2991744 : do j=js2g0,jn1g1 ! wk1 needed on N
963 2389632 : sldv(j) = .false.
964 2991744 : if( cose(j) < zt_d ) then
965 217004260 : do i=1,im
966 217004260 : if( abs(ub(i,j)) > D1_0 ) then ! ub ghosted on N
967 25463 : sldv(j) = .true.
968 : #if ( !defined UNICOSMP ) || ( !defined NEC_SX )
969 25463 : 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 602112 : jfirst, jlast+1, jfirst-1, jn1g1)
982 2991744 : do j=js2g0,jn1g1 ! wk1 needed on N
983 691205760 : do i=1,im
984 690603648 : 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 602112 : 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 602112 : if (ldiv4) then
1045 : !
1046 : ! filter velocity components for stability
1047 : !
1048 602112 : call pft2d(u(1,js2gd), grid%sediv4, grid%dediv4, im, jn1gs-js2gd+1, &
1049 1204224 : wkdiv4, wk2div4 )
1050 :
1051 602112 : call pft2d(v(1,js2gs), grid%scdiv4, grid%dcdiv4, im, jn2gd-js2gs+1, &
1052 1204224 : wkdiv4, wk2div4 )
1053 :
1054 : !**************************************************************************
1055 : !
1056 : ! div4 damping
1057 : !
1058 : !**************************************************************************
1059 :
1060 4760448 : do j=max(2,jfirst-2), min(jm-1,grid%jlast+2) ! fy need on NS (below)
1061 1202361216 : do i=1,im
1062 1201759104 : fy(i,j) = v(i,j)*cosp(j) ! v ghosted on NS at least
1063 : enddo
1064 : enddo
1065 :
1066 4177152 : do j=max(2,jfirst-1),min(jm,grid%jlast+2)
1067 : ! i=1
1068 3575040 : uc(1,j) = u(im,j) - u(1,j) ! u ghosted on N at least
1069 1030213632 : do i=2,im
1070 1029611520 : uc(i,j) = u(i-1,j) - u(i,j)
1071 : enddo
1072 : enddo
1073 : !
1074 : ! compute divergence
1075 : !
1076 602112 : if ( jfirst == 1 ) then
1077 : ! j=2
1078 2718912 : do i=1,im
1079 2718912 : div(i,2) = - cdydiv(2)*fy(i, 2) + cdxdiv(2)*uc(i,2)
1080 : enddo
1081 : endif
1082 :
1083 4158336 : do j=max(3,jfirst-1),min(jm-1,grid%jlast+2) ! wk1 needed on N (after TP2D)
1084 1028350848 : do i=1,im
1085 1027748736 : div(i,j) = cdydiv(j)*(fy(i,j-1) - fy(i,j)) + cdxdiv(j)*uc(i,j)
1086 : enddo
1087 : enddo
1088 :
1089 602112 : if ( jlast == jm ) then
1090 2718912 : do i=1,im
1091 2718912 : div(i,jm) = cdydiv(jm)*fy(i,jm-1) + cdxdiv(jm)*uc(i,jm)
1092 : enddo
1093 : endif
1094 :
1095 602112 : if ( jfirst == 1 ) then
1096 9408 : i=1
1097 9408 : j=2
1098 18816 : deldiv = cdx4(j) * (div(i+1,j )-D2_0*div(i,j)+div(im ,j ))+&
1099 18816 : cdy4(j) * (cosp(j)*(div(i,j+1)-div(i,j)))
1100 9408 : wk1(i,j) = wk1(i,j) +cdtau4(j)*deldiv
1101 2700096 : do i=2,im-1
1102 8072064 : deldiv = cdx4(j) * (div(i+1,j )-D2_0*div(i,j)+div(i-1,j ))+&
1103 10762752 : cdy4(j) * (cosp(j )*(div(i ,j+1)-div(i ,j)))
1104 2700096 : wk1(i,j) = wk1(i,j) + cdtau4(j)*deldiv
1105 : enddo
1106 9408 : i=im
1107 9408 : deldiv = cdx4(j) * (div(1 ,j )-D2_0*div(i,j)+div(i-1,j ))+&
1108 18816 : cdy4(j) * (cosp(j )*(div(i,j+1)-div(i,j)))
1109 9408 : wk1(i,j) = wk1(i,j) + cdtau4(j)*deldiv
1110 : endif
1111 :
1112 2972928 : do j=max(3,jfirst),jn2g1 ! wk1 needed on N (after TP2D)
1113 2370816 : i=1
1114 4741632 : deldiv = cdx4(j) * (div(i+1,j )-D2_0*div(i,j)+div(im ,j ))+&
1115 : cdy4(j) * ( &
1116 2370816 : cosp(j )*(div(i ,j+1)-div(i,j ))-&
1117 9483264 : cosp(j-1)*(div(i ,j )-div(i,j-1)))
1118 2370816 : wk1(i,j) = wk1(i,j) +cdtau4(j)*deldiv
1119 680424192 : do i=2,im-1
1120 2034160128 : 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 2712213504 : cosp(j-1)*(div(i ,j )-div(i ,j-1)))
1124 680424192 : wk1(i,j) = wk1(i,j) + cdtau4(j)*deldiv
1125 : enddo
1126 2370816 : i=im
1127 2370816 : 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 4741632 : cosp(j-1)*(div(i ,j )-div(i,j-1)))
1131 2972928 : wk1(i,j) = wk1(i,j) + cdtau4(j)*deldiv
1132 : enddo
1133 :
1134 602112 : if ( jlast == jm ) then
1135 9408 : i=1
1136 9408 : j = jm
1137 18816 : deldiv = cdx4(j) * (div(i+1,j )-D2_0*div(i,j)+div(im,j ))+&
1138 28224 : cdy4(j) * (-cosp(j-1)*(div(i,j)-div(i,j-1)))
1139 9408 : wk1(i,j) = wk1(i,j) +cdtau4(j)*deldiv
1140 2700096 : do i=2,im-1
1141 8072064 : deldiv = cdx4(j) * (div(i+1,j )-D2_0*div(i,j)+div(i-1,j ))+&
1142 10762752 : cdy4(j) * (-cosp(j-1)*(div(i,j)-div(i,j-1)))
1143 2700096 : wk1(i,j) = wk1(i,j) + cdtau4(j)*deldiv
1144 : enddo
1145 9408 : i=im
1146 9408 : j=jm
1147 9408 : deldiv = cdx4(j) * (div(1,j )-D2_0*div(i,j)+div(i-1,j ))+&
1148 18816 : cdy4(j) * (-cosp(j-1)*(div(i,j)-div(i,j-1)))
1149 9408 : wk1(i,j) = wk1(i,j) +cdtau4(j)*deldiv
1150 : endif
1151 : end if
1152 :
1153 870653952 : wku(:,:) = D0_0
1154 870653952 : wkv(:,:) = D0_0
1155 602112 : 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 5945856 : do j=js2gd,jn1gd
1244 1544944128 : do i=1,im
1245 1544342016 : delpf(i,j) = u(i,j)*cose(j) ! u ghosted on N*ng S*ng
1246 : enddo
1247 : enddo
1248 :
1249 602112 : if ( jfirst-ng_d <= 1 ) then
1250 18816 : c1 = D0_0
1251 5437824 : do i=1,im
1252 5437824 : c1 = c1 + delpf(i,2)
1253 : end do
1254 18816 : c1 = -c1*rdy*rcap
1255 5437824 : do i=1,im
1256 5437824 : uc(i,1) = c1
1257 : enddo
1258 : endif
1259 :
1260 602112 : if ( jlast+ng_d >= jm ) then
1261 18816 : c2 = D0_0
1262 5437824 : do i=1,im
1263 5437824 : c2 = c2 + delpf(i,jm)
1264 : end do
1265 18816 : c2 = c2*rdy*rcap
1266 5437824 : do i=1,im
1267 5437824 : uc(i,jm) = c2
1268 : enddo
1269 : else
1270 : ! This is an attempt to avoid ghosting u on N*(ng+1)
1271 168572544 : do i=1,im
1272 168572544 : uc(i,jn2gd) = D1E30
1273 : enddo
1274 : endif
1275 :
1276 5343744 : do j=js2gd, min(jm-1,jlast+ng_d-1)
1277 1365590016 : do i=1,im-1
1278 4082545152 : uc(i ,j) = ( delpf(i,j) - delpf(i ,j+1)) * cy(j) + &
1279 5448135168 : (v(i+1,j) - v(i ,j)) * rdx(j)
1280 : enddo
1281 14224896 : uc(im,j) = (delpf(im,j) - delpf(im,j+1)) * cy(j) + &
1282 19568640 : (v( 1 ,j) - v(im,j)) * rdx(j)
1283 : enddo
1284 :
1285 : ! uc is relative vorticity at this point
1286 :
1287 5964672 : do j=max(1,jfirst-ng_d), jn1gd
1288 1550381952 : do i=1,im
1289 1549779840 : uc(i,j) = uc(i,j) + f0(j)
1290 : ! uc is absolute vorticity
1291 : enddo
1292 : enddo
1293 :
1294 1204224 : call tp2d(va(1,jfirst), uc(1,jfirst-ng_d), crx(1,jfirst-ng_d), &
1295 602112 : cry(1,jfirst), im, jm, iord, jord, ng_d, fx, &
1296 602112 : fy (1,jfirst), ffsl, crx(1,jfirst), &
1297 1806336 : ymass, cosp, 0, jfirst, jlast)
1298 :
1299 4816896 : do j = jfirst-2, jlast+2
1300 1218674688 : do i = 1, im
1301 1218072576 : vf (i,j) = 0.0_r8
1302 : end do
1303 : end do
1304 :
1305 : ! AM correction
1306 602112 : 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 2399040 : do j = js2g0, jlast
1358 517515264 : do i=1,im-1
1359 517515264 : uc(i ,j) = dtdxe(j)*(wk1(i ,j)-wk1(i+1,j)) +dyce(j)*fy(i ,j) +wku(i ,j)
1360 : enddo
1361 2399040 : uc(im,j) = dtdxe(j)*(wk1(im,j)-wk1( 1,j)) +dyce(j)*fy(im,j) +wku(im,j)
1362 : enddo
1363 :
1364 2389632 : do j = js2g0, jn2g0
1365 517195392 : do i=1,im
1366 516593280 : 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 602112 : end subroutine d_sw
1371 : !-----------------------------------------------------------------------
1372 :
1373 : !-----------------------------------------------------------------------
1374 : !BOP
1375 : ! !IROUTINE: d2a2c_winds --- Interpolate winds
1376 : !
1377 : ! !INTERFACE:
1378 602112 : 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 1204224 : real(r8) uanp(grid%im), uasp(grid%im), vanp(grid%im), vasp(grid%im), r2im
1421 :
1422 602112 : real(r8), pointer:: sinlon(:)
1423 602112 : real(r8), pointer:: coslon(:)
1424 602112 : real(r8), pointer:: sinl5(:)
1425 602112 : real(r8), pointer:: cosl5(:)
1426 :
1427 : ! AM correction
1428 602112 : real(r8), pointer:: cosp(:)
1429 602112 : real(r8), pointer:: acosp(:)
1430 602112 : 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 602112 : im = grid%im
1440 602112 : jm = grid%jm
1441 602112 : jfirst = grid%jfirst
1442 602112 : jlast = grid%jlast
1443 :
1444 602112 : ng_c = grid%ng_c
1445 602112 : ng_d = grid%ng_d
1446 602112 : ng_s = grid%ng_s
1447 :
1448 602112 : im2 = im/2
1449 :
1450 602112 : js1gc = max(1,jfirst-ng_c) ! ng_c lats on S (starting at 1)
1451 602112 : jn1gc = min(jm,jlast+ng_c) ! ng_c latitudes on N (ending at jm)
1452 602112 : js2gc = max(2,jfirst-ng_c) ! ng_c lats on S (starting at 2)
1453 602112 : jn2gc = min(jm-1,jlast+ng_c) ! ng_c latitudes on N (ending at jm-1)
1454 602112 : js2gs = max(2,jfirst-ng_s) ! ng_s latitudes on S (starting at 2)
1455 602112 : jn2gsm1= min(jm-1,jlast+ng_s-1) ! ng_s-1 latitudes on N (ending at jm-1)
1456 602112 : js2gd = max(2,jfirst-ng_d) ! ng_d latitudes on S (starting at 2)
1457 602112 : jn2gd = min(jm-1,jlast+ng_d) ! ng_d latitudes on N (ending at jm-1)
1458 602112 : js2g2 = max(2,jfirst-2) ! 2 latitudes on S (starting at 2)
1459 602112 : jn1g2 = min(jm,jlast+2) ! 2 latitudes on N (ending at jm)
1460 :
1461 602112 : sinlon => grid%sinlon
1462 602112 : coslon => grid%coslon
1463 602112 : sinl5 => grid%sinl5
1464 602112 : cosl5 => grid%cosl5
1465 :
1466 : ! AM correction
1467 602112 : cosp => grid%cosp
1468 602112 : acosp => grid%acosp
1469 602112 : cose => grid%cose
1470 :
1471 : ! Get D-grid V-wind at the poles.
1472 :
1473 602112 : r2im = 0.5_r16/real(im,r16)
1474 :
1475 602112 : if ( jfirst-ng_d <= 1 ) then
1476 :
1477 : !
1478 : ! Treat SP
1479 : !
1480 5419008 : do i=1,im-1
1481 5400192 : uasp(i) = u(i,2) + u(i,3)
1482 5419008 : vasp(i) = v(i,2) + v(i+1,2)
1483 : enddo
1484 :
1485 18816 : uasp(im) = u(im,2) + u(im,3)
1486 18816 : vasp(im) = v(im,2) + v(1,2)
1487 :
1488 : ! Projection at SP
1489 :
1490 18816 : us = D0_0
1491 18816 : vs = D0_0
1492 :
1493 2728320 : do i=1,im2
1494 5419008 : us = us + (uasp(i+im2)-uasp(i))*sinlon(i) &
1495 8128512 : + (vasp(i)-vasp(i+im2))*coslon(i)
1496 : vs = vs + (uasp(i+im2)-uasp(i))*coslon(i) &
1497 2728320 : + (vasp(i+im2)-vasp(i))*sinlon(i)
1498 : enddo
1499 18816 : us = us*r2im
1500 18816 : vs = vs*r2im
1501 :
1502 : ! get V-wind at SP
1503 :
1504 2728320 : do i=1,im2
1505 2709504 : v(i, 1) = us*cosl5(i) - vs*sinl5(i)
1506 2728320 : v(i+im2,1) = -v(i, 1)
1507 : enddo
1508 :
1509 : endif
1510 :
1511 602112 : if ( jlast+ng_d >= jm ) then
1512 :
1513 : !
1514 : ! Treat NP
1515 : !
1516 5419008 : do i=1,im-1
1517 5400192 : uanp(i) = u(i,jm-1) + u(i,jm)
1518 5419008 : vanp(i) = v(i,jm-1) + v(i+1,jm-1)
1519 : enddo
1520 :
1521 18816 : uanp(im) = u(im,jm-1) + u(im,jm)
1522 18816 : vanp(im) = v(im,jm-1) + v(1,jm-1)
1523 :
1524 : ! Projection at NP
1525 :
1526 18816 : un = D0_0
1527 18816 : vn = D0_0
1528 2728320 : do i=1,im2
1529 5419008 : un = un + (uanp(i+im2)-uanp(i))*sinlon(i) &
1530 8128512 : + (vanp(i+im2)-vanp(i))*coslon(i)
1531 : vn = vn + (uanp(i)-uanp(i+im2))*coslon(i) &
1532 2728320 : + (vanp(i+im2)-vanp(i))*sinlon(i)
1533 : enddo
1534 18816 : un = un*r2im
1535 18816 : vn = vn*r2im
1536 :
1537 : ! get V-wind at NP
1538 :
1539 2728320 : do i=1,im2
1540 2709504 : v(i,jm) = -un*cosl5(i) - vn*sinl5(i)
1541 2728320 : v(i+im2,jm) = -v(i,jm)
1542 : enddo
1543 :
1544 : endif
1545 :
1546 1566695424 : ua(:,:) = D0_0
1547 1566695424 : va(:,:) = D0_0
1548 :
1549 5927040 : do j=js2gs,jn2gd
1550 1533579264 : do i=1,im-1
1551 1533579264 : va(i,j) = v(i,j) + v(i+1,j)
1552 : enddo
1553 5927040 : va(im,j) = v(im,j) + v(1,j)
1554 : enddo
1555 :
1556 602112 : 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 5343744 : do j = js2gd, jn2gsm1 ! WS: should be safe since u defined to jn2gs
1564 1370933760 : do i = 1, im
1565 1370331648 : 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 602112 : if ( reset_winds .and. met_rlx > D0_0 ) then
1574 391673856 : ua(:,:) = (D1_0-met_rlx)*ua(:,:) + met_rlx*( D2_0*u_cen(:,:) )
1575 391673856 : va(:,:) = (D1_0-met_rlx)*va(:,:) + met_rlx*( D2_0*v_cen(:,:) )
1576 : endif
1577 :
1578 602112 : if ( jfirst-ng_d <= 1 ) then
1579 : ! Projection at SP
1580 : us = D0_0
1581 : vs = D0_0
1582 :
1583 :
1584 2728320 : do i=1,im2
1585 8128512 : us = us + (ua(i+im2,2)-ua(i ,2))*sinlon(i) &
1586 8128512 : + (va(i ,2)-va(i+im2,2))*coslon(i)
1587 : vs = vs + (ua(i+im2,2)-ua(i ,2))*coslon(i) &
1588 2728320 : + (va(i+im2,2)-va(i ,2))*sinlon(i)
1589 : enddo
1590 :
1591 18816 : us = us/im
1592 18816 : vs = vs/im
1593 :
1594 : ! SP
1595 2728320 : do i=1,im2
1596 2709504 : ua(i,1) = -us*sinlon(i) - vs*coslon(i)
1597 2709504 : va(i,1) = us*coslon(i) - vs*sinlon(i)
1598 2709504 : ua(i+im2,1) = -ua(i,1)
1599 2728320 : va(i+im2,1) = -va(i,1)
1600 : enddo
1601 :
1602 : endif
1603 :
1604 602112 : if ( jlast+ng_d >= jm ) then
1605 :
1606 : ! Projection at NP
1607 : un = D0_0
1608 : vn = D0_0
1609 :
1610 2728320 : j = jm-1
1611 2728320 : do i=1,im2
1612 8128512 : un = un + (ua(i+im2,j)-ua(i ,j))*sinlon(i) &
1613 10838016 : + (va(i+im2,j)-va(i ,j))*coslon(i)
1614 : vn = vn + (ua(i ,j)-ua(i+im2,j))*coslon(i) &
1615 2728320 : + (va(i+im2,j)-va(i ,j))*sinlon(i)
1616 : enddo
1617 :
1618 18816 : un = un/im
1619 18816 : vn = vn/im
1620 :
1621 : ! NP
1622 2728320 : do i=1,im2
1623 2709504 : ua(i,jm) = -un*sinlon(i) + vn*coslon(i)
1624 2709504 : va(i,jm) = -un*coslon(i) - vn*sinlon(i)
1625 2709504 : ua(i+im2,jm) = -ua(i,jm)
1626 2728320 : va(i+im2,jm) = -va(i,jm)
1627 : enddo
1628 :
1629 : endif
1630 :
1631 :
1632 : ! A -> C
1633 4779264 : do j=js1gc,jn1gc ! uc needed N*ng_c S*ng_c, ua defined at poles
1634 4177152 : uc(1,j) = D0_25*(ua(1,j)+ua(im,j))
1635 1203621888 : do i=2,im
1636 1203019776 : uc(i,j) = D0_25*(ua(i,j)+ua(i-1,j))
1637 : enddo
1638 : enddo
1639 :
1640 602112 : 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 4769856 : do j = js2g2, jn1g2
1648 1205080128 : do i = 1, im
1649 1204478016 : vc(i,j) = D0_25*(va(i,j) + va(i,j-1))
1650 : end do
1651 : end do
1652 : end if
1653 :
1654 602112 : if ( jfirst==1 ) then
1655 2718912 : do i=1,im
1656 2718912 : vc(i,1) = D0_0 ! Needed to avoid undefined values in VC
1657 : enddo
1658 : endif
1659 : !EOC
1660 602112 : end subroutine d2a2c_winds
1661 : !-----------------------------------------------------------------------
1662 : end module sw_core
|
