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