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