Line data Source code
1 : module derivative_mod
2 : use shr_kind_mod, only: r8=>shr_kind_r8
3 : use cam_abortutils, only: endrun
4 : use dimensions_mod, only : np, nc, npdg, nelemd, nlev
5 : use quadrature_mod, only : quadrature_t, gauss, gausslobatto,legendre, jacobi
6 : ! needed for spherical differential operators:
7 : use physconst, only: ra
8 : use element_mod, only : element_t
9 : use control_mod, only : hypervis_scaling, hypervis_power
10 : use perf_mod, only : t_startf, t_stopf
11 :
12 : implicit none
13 : private
14 :
15 : type, public :: derivative_t
16 : real (kind=r8) :: Dvv(np,np)
17 : real (kind=r8) :: Dvv_diag(np,np)
18 : real (kind=r8) :: Dvv_twt(np,np)
19 : real (kind=r8) :: Mvv_twt(np,np) ! diagonal matrix of GLL weights
20 : real (kind=r8) :: Mfvm(np,nc+1)
21 : real (kind=r8) :: Cfvm(np,nc)
22 : real (kind=r8) :: legdg(np,np)
23 : end type derivative_t
24 :
25 : type, public :: derivative_stag_t
26 : real (kind=r8) :: D(np,np)
27 : real (kind=r8) :: M(np,np)
28 : real (kind=r8) :: Dpv(np,np)
29 : real (kind=r8) :: D_twt(np,np)
30 : real (kind=r8) :: M_twt(np,np)
31 : real (kind=r8) :: M_t(np,np)
32 : end type derivative_stag_t
33 :
34 : real (kind=r8), allocatable :: integration_matrix(:,:)
35 : real (kind=r8), allocatable :: integration_matrix_physgrid(:,:)
36 : real (kind=r8), allocatable :: boundary_interp_matrix(:,:,:)
37 :
38 : ! ======================================
39 : ! Public Interfaces
40 : ! ======================================
41 :
42 : public :: subcell_integration
43 : public :: subcell_dss_fluxes
44 : public :: subcell_div_fluxes
45 : public :: subcell_Laplace_fluxes
46 : public :: allocate_subcell_integration_matrix_cslam !for consistent se-cslam algorithm
47 : public :: allocate_subcell_integration_matrix_physgrid !for integration se basis functions over physgrid control volumes
48 :
49 : public :: derivinit
50 :
51 : public :: gradient
52 : public :: gradient_wk
53 : public :: vorticity
54 : public :: divergence
55 :
56 : public :: interpolate_gll2fvm_corners
57 : public :: interpolate_gll2fvm_points
58 : public :: remap_phys2gll
59 :
60 :
61 : interface divergence
62 : module procedure divergence_nonstag
63 : module procedure divergence_stag
64 : end interface
65 :
66 : interface gradient
67 : module procedure gradient_str_nonstag
68 : module procedure gradient_str_stag
69 : end interface
70 :
71 : interface gradient_wk
72 : module procedure gradient_wk_nonstag
73 : module procedure gradient_wk_stag
74 : end interface
75 :
76 : public :: v2pinit
77 :
78 : private :: dmatinit
79 : private :: dvvinit
80 : private :: dpvinit
81 :
82 : ! these routines compute spherical differential operators as opposed to
83 : ! the gnomonic coordinate operators above. Vectors (input or output)
84 : ! are always expressed in lat-lon coordinates
85 : !
86 : ! note that weak derivatives (integrated by parts form) can be defined using
87 : ! contra or co-variant test functions, so
88 : !
89 : public :: gradient_sphere
90 : public :: gradient_sphere_wk_testcov
91 : public :: gradient_sphere_wk_testcontra ! only used for debugging
92 : public :: ugradv_sphere
93 : public :: vorticity_sphere
94 : public :: vorticity_sphere_diag
95 : public :: divergence_sphere
96 : public :: curl_sphere
97 : public :: curl_sphere_wk_testcov
98 : public :: divergence_sphere_wk
99 : public :: laplace_sphere_wk
100 : public :: vlaplace_sphere_wk
101 : public :: vlaplace_sphere_wk_mol
102 : public :: element_boundary_integral
103 : public :: edge_flux_u_cg
104 : public :: gll_to_dgmodal
105 : public :: dgmodal_to_gll
106 :
107 : public :: limiter_optim_iter_full
108 :
109 : contains
110 :
111 : ! ==========================================
112 : ! derivinit:
113 : !
114 : ! Initialize the matrices for taking
115 : ! derivatives and interpolating
116 : ! ==========================================
117 :
118 1536 : subroutine derivinit(deriv,fvm_corners, fvm_points)
119 : type (derivative_t) :: deriv
120 : real (kind=r8),optional :: fvm_corners(nc+1)
121 : real (kind=r8),optional :: fvm_points(nc)
122 :
123 : ! Local variables
124 : type (quadrature_t) :: gp ! Quadrature points and weights on pressure grid
125 :
126 : real (kind=r8) :: dmat(np,np)
127 : real (kind=r8) :: dpv(np,np)
128 : real (kind=r8) :: v2p(np,np)
129 : real (kind=r8) :: p2v(np,np)
130 : real (kind=r8) :: dvv(np,np)
131 : real (kind=r8) :: dvv_diag(np,np)
132 : real (kind=r8) :: v2v(np,np)
133 : real (kind=r8) :: xnorm
134 : integer i,j
135 :
136 : ! ============================================
137 : ! initialize matrices in r8 precision
138 : ! and transfer results into r8
139 : ! floating point precision
140 : ! ============================================
141 :
142 1536 : gp=gausslobatto(np)
143 :
144 : ! Legendre polynomials of degree npdg-1, on the np GLL grid:
145 1536 : if (npdg>np) call endrun( 'FATAL ERROR: npdg>np')
146 1536 : if (npdg>0 .and. npdg<np) then
147 : ! in this case, we use a DG basis of Legendre polynomials
148 : ! stored at the GLL points.
149 0 : do i=1,np
150 0 : deriv%legdg(1:npdg,i) = legendre(gp%points(i),npdg-1)
151 : end do
152 : ! normalize
153 0 : do j=1,npdg
154 0 : xnorm=sqrt(sum(deriv%legdg(j,:)*deriv%legdg(j,:)*gp%weights(:)))
155 0 : deriv%legdg(j,:)=deriv%legdg(j,:)/xnorm
156 : enddo
157 : endif
158 1536 : call dvvinit(dvv,gp)
159 32256 : deriv%Dvv(:,:) = dvv(:,:)
160 :
161 7680 : do i=1,np
162 32256 : do j=1,np
163 30720 : if (i.eq.j) then
164 6144 : deriv%dvv_diag(i,j) = dvv(i,j)
165 : else
166 18432 : deriv%dvv_diag(i,j) = 0.0_r8
167 : endif
168 : end do
169 : end do
170 :
171 1536 : v2v = 0.0_r8
172 7680 : do i=1,np
173 7680 : v2v(i,i) = gp%weights(i)
174 : end do
175 :
176 7680 : do i=1,np
177 32256 : do j=1,np
178 30720 : dvv(j,i) = dvv(j,i)*gp%weights(i)
179 : end do
180 : end do
181 :
182 32256 : deriv%Dvv_twt = TRANSPOSE(dvv)
183 32256 : deriv%Mvv_twt = v2v
184 1536 : if (present(fvm_corners)) &
185 1536 : call v2pinit(deriv%Mfvm,gp%points,fvm_corners,np,nc+1)
186 :
187 1536 : if (present(fvm_points)) &
188 1536 : call v2pinit(deriv%Cfvm,gp%points,fvm_points,np,nc)
189 :
190 : ! notice we deallocate this memory here even though it was allocated
191 : ! by the call to gausslobatto.
192 1536 : deallocate(gp%points)
193 1536 : deallocate(gp%weights)
194 :
195 1536 : end subroutine derivinit
196 :
197 : ! =======================================
198 : ! dmatinit:
199 : !
200 : ! Compute rectangular v->p
201 : ! derivative matrix (dmat)
202 : ! =======================================
203 :
204 : subroutine dmatinit(dmat)
205 :
206 : real (kind=r8) :: dmat(np,np)
207 :
208 : ! Local variables
209 :
210 : type (quadrature_t) :: gll
211 : type (quadrature_t) :: gs
212 :
213 : integer i,j
214 : real(kind=r8) fact,f1,f2
215 : real(kind=r8) func0,func1
216 : real(kind=r8) dis,c0,c1
217 :
218 : real(kind=r8) :: leg(np,np)
219 : real(kind=r8) :: jac(0:np-1)
220 : real(kind=r8) :: djac(0:np-1)
221 :
222 : c0 = 0.0_r8
223 : c1 = 1.0_r8
224 :
225 : gll= gausslobatto(np)
226 : gs = gauss(np)
227 :
228 : ! =============================================================
229 : ! Compute Legendre polynomials on Gauss-Lobatto grid (velocity)
230 : ! =============================================================
231 :
232 : do i=1,np
233 : leg(:,i) = legendre(gll%points(i),np-1)
234 : end do
235 :
236 : ! ================================================================
237 : ! Derivatives of velocity cardinal functions on pressure grid
238 : ! d(i,j) = D(j,i) = D' (D-transpose) since D(i,j) = dh_j(x_i)/dx
239 : ! ================================================================
240 :
241 : fact = np*(np-1)
242 :
243 : do j=1,np
244 : call jacobi(np-1,gs%points(j),c0,c0,jac(0:np-1),djac(0:np-1))
245 : func0 = jac(np-1)
246 : func1 = djac(np-1)
247 : f1 = fact*func0
248 : f2 = (c1 - gs%points(j))*(c1 + gs%points(j)) * func1
249 : do i = 1, np
250 : if ( gs%points(j) /= gll%points(i) ) then
251 : dis = gs%points(j) - gll%points(i)
252 : dmat(i,j) = func0 / ( leg(np,i)*dis ) + f2 / (fact*leg(np,i)*dis*dis)
253 : else
254 : dmat(i,j) = c0
255 : endif
256 : end do
257 : end do
258 :
259 : deallocate(gll%points)
260 : deallocate(gll%weights)
261 :
262 : deallocate(gs%points)
263 : deallocate(gs%weights)
264 :
265 : end subroutine dmatinit
266 :
267 : ! =======================================
268 : ! dpvinit:
269 : !
270 : ! Compute rectangular p->v
271 : ! derivative matrix (dmat)
272 : ! for strong gradients
273 : ! =======================================
274 :
275 : subroutine dpvinit(dmat)
276 :
277 : real (kind=r8) :: dmat(np,np)
278 :
279 : ! Local variables
280 :
281 : type (quadrature_t) :: gll
282 : type (quadrature_t) :: gs
283 :
284 : integer i,j
285 : real(kind=r8) dis,c0,c1
286 :
287 : real(kind=r8) :: legv(0:np,np)
288 : real(kind=r8) :: dlegv(0:np,np)
289 :
290 : real(kind=r8) :: leg(0:np)
291 : real(kind=r8) :: dleg(0:np)
292 :
293 : c0 = 0.0_r8
294 : c1 = 1.0_r8
295 :
296 : gll= gausslobatto(np)
297 : gs = gauss(np)
298 :
299 : ! =============================================================
300 : ! Compute Legendre polynomials on Gauss-Lobatto grid (velocity)
301 : ! =============================================================
302 :
303 : do i=1,np
304 : call jacobi(np,gll%points(i),c0,c0,legv(0:np,i),dlegv(0:np,i))
305 : end do
306 :
307 : ! ================================================================
308 : ! Derivatives of velocity cardinal functions on pressure grid
309 : ! d(i,j) = D(j,i) = D' (D-transpose) since D(i,j) = dh_j(x_i)/dx
310 : ! ================================================================
311 :
312 : do j=1,np
313 : call jacobi(np,gs%points(j),c0,c0,leg(0:np),dleg(0:np))
314 : do i = 1, np
315 : if ( gs%points(j) /= gll%points(i) ) then
316 : dis = gll%points(i) - gs%points(j)
317 : dmat(j,i) = dlegv(np,i)/( dleg(np)*dis ) - legv(np,i)/ (dleg(np)*dis*dis)
318 : else
319 : dmat(j,i) = c0
320 : endif
321 : end do
322 : end do
323 :
324 : deallocate(gll%points)
325 : deallocate(gll%weights)
326 :
327 : deallocate(gs%points)
328 : deallocate(gs%weights)
329 :
330 : end subroutine dpvinit
331 :
332 : ! =======================================
333 : ! v2pinit:
334 : ! Compute interpolation matrix from gll(1:n1) -> gs(1:n2)
335 : ! =======================================
336 3072 : subroutine v2pinit(v2p,gll,gs,n1,n2)
337 : integer :: n1,n2
338 : real(kind=r8) :: v2p(n1,n2)
339 6144 : real(kind=r8) :: v2p_new(n1,n2)
340 : real(kind=r8) :: gll(n1),gs(n2)
341 : ! Local variables
342 :
343 : integer i,j,k,m,l
344 : real(kind=r8) fact,f1, sum
345 : real(kind=r8) func0,func1
346 :
347 6144 : real(kind=r8) :: leg(n1,n1)
348 : real(kind=r8) :: jac(0:n1-1)
349 : real(kind=r8) :: djac(0:n1-1)
350 : real(kind=r8) :: c0,c1
351 :
352 : type(quadrature_t) :: gll_pts
353 6144 : real(kind=r8) :: leg_out(n1,n2)
354 6144 : real(kind=r8) :: gamma(n1)
355 :
356 3072 : c0 = 0.0_r8
357 3072 : c1 = 1.0_r8
358 :
359 : ! ==============================================================
360 : ! Compute Legendre polynomials on Gauss-Lobatto grid (velocity)
361 : ! ==============================================================
362 :
363 3072 : fact = -n1*(n1-1)
364 15360 : do i=1,n1
365 12288 : leg(:,i) = legendre(gll(i),n1-1)
366 15360 : leg(n1,i) = fact * leg(n1,i)
367 : end do
368 :
369 : ! ===================================================
370 : ! Velocity cardinal functions on pressure grid
371 : ! ===================================================
372 : ! NEW VERSION, with no division by (gs(j)-gll(i)):
373 :
374 : ! compute legendre polynomials at output points:
375 3072 : gll_pts = gausslobatto(n1)
376 :
377 3072 : fact = -n1*(n1-1)
378 13824 : do i=1,n2
379 10752 : leg_out(:,i) = legendre(gs(i),n1-1)
380 13824 : leg_out(n1,i) = fact * leg_out(n1,i)
381 : end do
382 :
383 :
384 : ! compute gamma: (normalization factor for inv(leg)
385 15360 : do m=1,n1
386 12288 : gamma(m)=0
387 61440 : do i=1,n1
388 61440 : gamma(m)=gamma(m)+leg(m,i)*leg(m,i)*gll_pts%weights(i)
389 : enddo
390 15360 : gamma(m)=1/gamma(m)
391 : enddo
392 :
393 : ! compute product of leg_out * inv(leg):
394 13824 : do j=1,n2 ! this should be fvm points
395 56832 : do l=1,n1 ! this should be GLL points
396 43008 : sum=0
397 215040 : do k=1,n1 ! number of polynomials = number of GLL points
398 215040 : sum=sum + leg_out(k,j)*gamma(k)*leg(k,l)
399 : enddo
400 53760 : v2p_new(l,j) = gll_pts%weights(l)*sum
401 : enddo
402 : enddo
403 3072 : deallocate(gll_pts%points)
404 3072 : deallocate(gll_pts%weights)
405 :
406 56832 : v2p=v2p_new
407 3072 : end subroutine v2pinit
408 :
409 :
410 :
411 : ! =======================================
412 : ! dvvinit:
413 : !
414 : ! Compute rectangular v->v
415 : ! derivative matrix (dvv)
416 : ! =======================================
417 :
418 1536 : subroutine dvvinit(dvv,gll)
419 :
420 : real(kind=r8) :: dvv(np,np)
421 : type (quadrature_t) :: gll
422 :
423 : ! Local variables
424 :
425 : real(kind=r8) :: leg(np,np)
426 : real(kind=r8) :: c0,c1,c4
427 :
428 : integer i,j
429 :
430 1536 : c0 = 0.0_r8
431 1536 : c1 = 1.0_r8
432 1536 : c4 = 4.0_r8
433 :
434 7680 : do i=1,np
435 7680 : leg(:,i) = legendre(gll%points(i),np-1)
436 : end do
437 :
438 32256 : dvv(:,:) = c0
439 7680 : do j=1,np
440 15360 : do i=1,j-1
441 15360 : dvv(j,i) = (c1/(gll%points(i)-gll%points(j)))*leg(np,i)/leg(np,j)
442 : end do
443 6144 : dvv(j,j) = c0
444 16896 : do i=j+1,np
445 15360 : dvv(j,i) = (c1/(gll%points(i)-gll%points(j)))*leg(np,i)/leg(np,j)
446 : end do
447 : end do
448 :
449 :
450 1536 : dvv(np,np) = + np*(np-1)/c4
451 1536 : dvv(1,1) = - np*(np-1)/c4
452 :
453 1536 : end subroutine dvvinit
454 :
455 : ! ================================================
456 : ! divergence_stag:
457 : !
458 : ! Compute divergence (maps v grid -> p grid)
459 : ! ================================================
460 :
461 0 : subroutine divergence_stag(v,deriv,div)
462 :
463 : real(kind=r8), intent(in) :: v(np,np,2)
464 : type (derivative_stag_t), intent(in) :: deriv
465 : real(kind=r8), intent(out) :: div(np,np)
466 :
467 : ! Local
468 :
469 : integer i
470 : integer j
471 : integer l
472 :
473 : real(kind=r8) sumx00
474 : real(kind=r8) sumy00
475 :
476 : real(kind=r8) :: vtemp(np,np,2)
477 :
478 :
479 0 : do j=1,np
480 0 : do l=1,np
481 :
482 : sumx00=0.0d0
483 : sumy00=0.0d0
484 : !DIR$ UNROLL(NP)
485 0 : do i=1,np
486 0 : sumx00 = sumx00 + deriv%D(i,l )*v(i,j ,1)
487 0 : sumy00 = sumy00 + deriv%M(i,l )*v(i,j ,2)
488 : enddo
489 0 : vtemp(j ,l ,1) = sumx00
490 0 : vtemp(j ,l ,2) = sumy00
491 : enddo
492 : enddo
493 0 : do j=1,np
494 0 : do i=1,np
495 : sumx00=0.0d0
496 : sumy00=0.0d0
497 : !DIR$ UNROLL(NP)
498 0 : do l=1,np
499 0 : sumx00 = sumx00 + deriv%M(l,j )*vtemp(l,i ,1)
500 0 : sumy00 = sumy00 + deriv%D(l,j )*vtemp(l,i ,2)
501 : enddo
502 0 : div(i ,j ) = sumx00 + sumy00
503 :
504 : enddo
505 : enddo
506 :
507 0 : end subroutine divergence_stag
508 :
509 : ! ================================================
510 : ! divergence_nonstag:
511 : !
512 : ! Compute divergence (maps v->v)
513 : ! ================================================
514 :
515 0 : subroutine divergence_nonstag(v,deriv,div)
516 :
517 : real(kind=r8), intent(in) :: v(np,np,2)
518 : type (derivative_t), intent(in) :: deriv
519 :
520 : real(kind=r8), intent(out) :: div(np,np)
521 :
522 : ! Local
523 :
524 : integer i
525 : integer j
526 : integer l
527 :
528 : real(kind=r8) :: dudx00
529 : real(kind=r8) :: dvdy00
530 :
531 : real(kind=r8) :: vvtemp(np,np)
532 :
533 0 : do j=1,np
534 0 : do l=1,np
535 : dudx00=0.0d0
536 : dvdy00=0.0d0
537 : !DIR$ UNROLL(NP)
538 0 : do i=1,np
539 0 : dudx00 = dudx00 + deriv%Dvv(i,l )*v(i,j ,1)
540 0 : dvdy00 = dvdy00 + deriv%Dvv(i,l )*v(j ,i,2)
541 : end do
542 :
543 0 : div(l ,j ) = dudx00
544 0 : vvtemp(j ,l ) = dvdy00
545 : end do
546 : end do
547 0 : do j=1,np
548 0 : do i=1,np
549 0 : div(i,j)=div(i,j)+vvtemp(i,j)
550 : end do
551 : end do
552 :
553 0 : end subroutine divergence_nonstag
554 :
555 : ! ================================================
556 : ! gradient_wk_stag:
557 : !
558 : ! Compute the weak form gradient:
559 : ! maps scalar field on the pressure grid to the
560 : ! velocity grid
561 : ! ================================================
562 :
563 0 : function gradient_wk_stag(p,deriv) result(dp)
564 :
565 : type (derivative_stag_t), intent(in) :: deriv
566 : real(kind=r8), intent(in) :: p(np,np)
567 :
568 : real(kind=r8) :: dp(np,np,2)
569 :
570 : ! Local
571 :
572 : integer i
573 : integer j
574 : integer l
575 :
576 : real(kind=r8) sumx00,sumx01
577 : real(kind=r8) sumy00,sumy01
578 :
579 : real(kind=r8) :: vtempt(np,np,2)
580 :
581 0 : do j=1,np
582 0 : do l=1,np
583 : sumx00=0.0d0
584 : sumy00=0.0d0
585 : !DIR$ UNROLL(NP)
586 0 : do i=1,np
587 0 : sumx00 = sumx00 + deriv%D_twt(i,l )*p(i,j )
588 0 : sumy00 = sumy00 + deriv%M_twt(i,l )*p(i,j )
589 : enddo
590 0 : vtempt(j ,l ,1) = sumx00
591 0 : vtempt(j ,l ,2) = sumy00
592 : enddo
593 : enddo
594 0 : do j=1,np
595 0 : do i=1,np
596 : sumx00=0.0d0
597 : sumy00=0.0d0
598 : !DIR$ UNROLL(NP)
599 0 : do l=1,np
600 0 : sumx00 = sumx00 + deriv%M_twt(l,j )*vtempt(l,i ,1)
601 0 : sumy00 = sumy00 + deriv%D_twt(l,j )*vtempt(l,i ,2)
602 : enddo
603 0 : dp(i ,j ,1) = sumx00
604 0 : dp(i ,j ,2) = sumy00
605 : enddo
606 : enddo
607 :
608 :
609 0 : end function gradient_wk_stag
610 :
611 : ! ================================================
612 : ! gradient_wk_nonstag:
613 : !
614 : ! Compute the weak form gradient:
615 : ! maps scalar field on the Gauss-Lobatto grid to the
616 : ! weak gradient on the Gauss-Lobbatto grid
617 : ! ================================================
618 :
619 0 : function gradient_wk_nonstag(p,deriv) result(dp)
620 :
621 : type (derivative_t), intent(in) :: deriv
622 : real(kind=r8), intent(in) :: p(np,np)
623 :
624 : real(kind=r8) :: dp(np,np,2)
625 :
626 : ! Local
627 :
628 : integer i
629 : integer j
630 : integer l
631 :
632 : real(kind=r8) sumx00
633 : real(kind=r8) sumy00
634 :
635 : real(kind=r8) :: vvtempt(np,np,2)
636 :
637 0 : do j=1,np
638 0 : do l=1,np
639 : sumx00=0.0d0
640 : sumy00=0.0d0
641 : !DIR$ UNROLL(NP)
642 0 : do i=1,np
643 0 : sumx00 = sumx00 + deriv%Dvv_twt(i,l )*p(i,j )
644 0 : sumy00 = sumy00 + deriv%Mvv_twt(i,l )*p(i,j )
645 : end do
646 0 : vvtempt(j ,l ,1) = sumx00
647 0 : vvtempt(j ,l ,2) = sumy00
648 : end do
649 : end do
650 :
651 0 : do j=1,np
652 0 : do i=1,np
653 : sumx00=0.0d0
654 : sumy00=0.0d0
655 : !DIR$ UNROLL(NP)
656 0 : do l=1,np
657 0 : sumx00 = sumx00 + deriv%Mvv_twt(l,j )*vvtempt(l,i ,1)
658 0 : sumy00 = sumy00 + deriv%Dvv_twt(l,j )*vvtempt(l,i ,2)
659 : end do
660 0 : dp(i ,j ,1) = sumx00
661 0 : dp(i ,j ,2) = sumy00
662 : end do
663 : end do
664 0 : end function gradient_wk_nonstag
665 :
666 : ! ================================================
667 : ! gradient_str_stag:
668 : !
669 : ! Compute the *strong* form gradient:
670 : ! maps scalar field on the pressure grid to the
671 : ! velocity grid
672 : ! ================================================
673 :
674 0 : subroutine gradient_str_stag(p,deriv,dp)
675 :
676 : type (derivative_stag_t), intent(in) :: deriv
677 : real(kind=r8), intent(in) :: p(np,np)
678 :
679 : real(kind=r8), intent(out) :: dp(np,np,2)
680 :
681 : ! Local
682 :
683 : integer i
684 : integer j
685 : integer l
686 :
687 : real(kind=r8) sumx00
688 : real(kind=r8) sumy00
689 :
690 : real(kind=r8) :: vtempt(np,np,2)
691 0 : do j=1,np
692 0 : do l=1,np
693 : sumx00=0.0d0
694 : sumy00=0.0d0
695 : !DIR$ UNROLL(NP)
696 0 : do i=1,np
697 0 : sumx00 = sumx00 + deriv%Dpv(i,l )*p(i,j )
698 0 : sumy00 = sumy00 + deriv%M_t(i,l )*p(i,j )
699 : enddo
700 0 : vtempt(j ,l ,1) = sumx00
701 0 : vtempt(j ,l ,2) = sumy00
702 : enddo
703 : enddo
704 0 : do j=1,np
705 0 : do i=1,np
706 : sumx00=0.0d0
707 : sumy00=0.0d0
708 : !DIR$ UNROLL(NP)
709 0 : do l=1,np
710 0 : sumx00 = sumx00 + deriv%M_t(l,j )*vtempt(l,i ,1)
711 0 : sumy00 = sumy00 + deriv%Dpv(l,j )*vtempt(l,i ,2)
712 : enddo
713 0 : dp(i ,j ,1) = sumx00
714 0 : dp(i ,j ,2) = sumy00
715 : enddo
716 : enddo
717 :
718 0 : end subroutine gradient_str_stag
719 :
720 : ! ================================================
721 : ! gradient_str_nonstag:
722 : !
723 : ! Compute the *strong* gradient on the velocity grid
724 : ! of a scalar field on the velocity grid
725 : ! ================================================
726 :
727 0 : subroutine gradient_str_nonstag(s,deriv,ds)
728 :
729 : type (derivative_t), intent(in) :: deriv
730 : real(kind=r8), intent(in) :: s(np,np)
731 : real(kind=r8), intent(out) :: ds(np,np,2)
732 :
733 : integer i
734 : integer j
735 : integer l
736 : real(kind=r8) :: dsdx00,dsdx01
737 : real(kind=r8) :: dsdy00,dsdy01
738 0 : do j=1,np
739 0 : do l=1,np
740 : dsdx00=0.0d0
741 : dsdy00=0.0d0
742 : !DIR$ UNROLL(NP)
743 0 : do i=1,np
744 0 : dsdx00 = dsdx00 + deriv%Dvv(i,l )*s(i,j )
745 0 : dsdy00 = dsdy00 + deriv%Dvv(i,l )*s(j ,i)
746 : end do
747 0 : ds(l ,j ,1) = dsdx00
748 0 : ds(j ,l ,2) = dsdy00
749 : end do
750 : end do
751 0 : end subroutine gradient_str_nonstag
752 :
753 : ! ================================================
754 : ! vorticity:
755 : !
756 : ! Compute the vorticity of the velocity field on the
757 : ! velocity grid
758 : ! ================================================
759 :
760 0 : subroutine vorticity(v,deriv,vort)
761 :
762 : type (derivative_t), intent(in) :: deriv
763 : real(kind=r8), intent(in) :: v(np,np,2)
764 :
765 : real(kind=r8), intent(out) :: vort(np,np)
766 :
767 : integer i
768 : integer j
769 : integer l
770 :
771 : real(kind=r8) :: dvdx00,dvdx01
772 : real(kind=r8) :: dudy00,dudy01
773 :
774 : real(kind=r8) :: vvtemp(np,np)
775 0 : do j=1,np
776 0 : do l=1,np
777 : dudy00=0.0d0
778 : dvdx00=0.0d0
779 : !DIR$ UNROLL(NP)
780 0 : do i=1,np
781 0 : dvdx00 = dvdx00 + deriv%Dvv(i,l )*v(i,j ,2)
782 0 : dudy00 = dudy00 + deriv%Dvv(i,l )*v(j ,i,1)
783 : enddo
784 0 : vort(l ,j ) = dvdx00
785 0 : vvtemp(j ,l ) = dudy00
786 : enddo
787 : enddo
788 0 : do j=1,np
789 0 : do i=1,np
790 0 : vort(i,j)=vort(i,j)-vvtemp(i,j)
791 : end do
792 : end do
793 :
794 0 : end subroutine vorticity
795 :
796 : ! ================================================
797 : ! interpolate_gll2fvm_points:
798 : !
799 : ! shape funtion interpolation from data on GLL grid to cellcenters on physics grid
800 : ! Author: Christoph Erath
801 : ! ================================================
802 0 : function interpolate_gll2fvm_points(v,deriv) result(p)
803 :
804 : real(kind=r8), intent(in) :: v(np,np)
805 : type (derivative_t) :: deriv
806 : real(kind=r8) :: p(nc,nc)
807 :
808 : ! Local
809 : integer i
810 : integer j
811 : integer l
812 :
813 : real(kind=r8) sumx00,sumx01
814 : real(kind=r8) sumx10,sumx11
815 : real(kind=r8) vtemp(np,nc)
816 :
817 0 : do j=1,np
818 0 : do l=1,nc
819 : sumx00=0.0d0
820 : !DIR$ UNROLL(NP)
821 0 : do i=1,np
822 0 : sumx00 = sumx00 + deriv%Cfvm(i,l )*v(i,j )
823 : enddo
824 0 : vtemp(j ,l) = sumx00
825 : enddo
826 : enddo
827 0 : do j=1,nc
828 0 : do i=1,nc
829 : sumx00=0.0d0
830 : !DIR$ UNROLL(NP)
831 0 : do l=1,np
832 0 : sumx00 = sumx00 + deriv%Cfvm(l,j )*vtemp(l,i)
833 : enddo
834 0 : p(i ,j ) = sumx00
835 : enddo
836 : enddo
837 0 : end function interpolate_gll2fvm_points
838 : ! ================================================
839 : ! interpolate_gll2fvm_corners:
840 : !
841 : ! shape funtion interpolation from data on GLL grid to physics grid
842 : !
843 : ! ================================================
844 0 : function interpolate_gll2fvm_corners(v,deriv) result(p)
845 :
846 : real(kind=r8), intent(in) :: v(np,np)
847 : type (derivative_t), intent(in) :: deriv
848 : real(kind=r8) :: p(nc+1,nc+1)
849 :
850 : ! Local
851 : integer i
852 : integer j
853 : integer l
854 :
855 : real(kind=r8) sumx00,sumx01
856 : real(kind=r8) sumx10,sumx11
857 : real(kind=r8) vtemp(np,nc+1)
858 :
859 0 : do j=1,np
860 0 : do l=1,nc+1
861 : sumx00=0.0d0
862 : !DIR$ UNROLL(NP)
863 0 : do i=1,np
864 0 : sumx00 = sumx00 + deriv%Mfvm(i,l )*v(i,j )
865 : enddo
866 0 : vtemp(j ,l) = sumx00
867 : enddo
868 : enddo
869 0 : do j=1,nc+1
870 0 : do i=1,nc+1
871 : sumx00=0.0d0
872 : !DIR$ UNROLL(NP)
873 0 : do l=1,np
874 0 : sumx00 = sumx00 + deriv%Mfvm(l,j )*vtemp(l,i)
875 : enddo
876 0 : p(i ,j ) = sumx00
877 : enddo
878 : enddo
879 0 : end function interpolate_gll2fvm_corners
880 :
881 :
882 : ! ================================================
883 : ! remap_phys2gll:
884 : !
885 : ! interpolate to an equally spaced (in reference element coordinate system)
886 : ! "physics" grid to the GLL grid
887 : !
888 : ! 1st order, monotone, conservative
889 : ! MT initial version 2013
890 : ! ================================================
891 0 : function remap_phys2gll(pin,nphys) result(pout)
892 : integer :: nphys
893 : real(kind=r8), intent(in) :: pin(nphys*nphys)
894 : real(kind=r8) :: pout(np,np)
895 :
896 : ! Local
897 : integer, save :: nphys_init=0
898 : integer, save :: nintersect
899 : real(kind=r8),save,pointer :: acell(:) ! arrivial cell index of i'th intersection
900 : real(kind=r8),save,pointer :: dcell(:) ! departure cell index of i'th intersection
901 : real(kind=r8),save,pointer :: delta(:) ! length of i'th intersection
902 : real(kind=r8),save,pointer :: delta_a(:) ! length of arrival cells
903 : integer in_i,in_j,ia,ja,id,jd,count,i,j
904 : logical :: found
905 :
906 : real(kind=r8) :: tol = 1.0e-13_r8
907 : real(kind=r8) :: weight,x1,x2,dx
908 0 : real(kind=r8) :: gll_edges(np+1),phys_edges(nphys+1)
909 : type(quadrature_t) :: gll_pts
910 0 : if (nphys_init/=nphys) then
911 : ! setup (most be done on masterthread only) since all data is static
912 : ! MT: move barrier inside if loop - we dont want a barrier every regular call
913 : !OMP BARRIER
914 : !OMP MASTER
915 0 : nphys_init=nphys
916 : ! find number of intersections
917 0 : nintersect = np+nphys-1 ! max number of possible intersections
918 0 : allocate(acell(nintersect))
919 0 : allocate(dcell(nintersect))
920 0 : allocate(delta(nintersect))
921 0 : allocate(delta_a(np))
922 :
923 : ! compute phys grid cell edges on [-1,1]
924 0 : do i=1,nphys+1
925 0 : dx = 2d0/nphys
926 0 : phys_edges(i)=-1 + (i-1)*dx
927 : enddo
928 :
929 : ! compute GLL cell edges on [-1,1]
930 0 : gll_pts = gausslobatto(np)
931 0 : gll_edges(1)=-1
932 0 : do i=2,np
933 0 : gll_edges(i) = gll_edges(i-1) + gll_pts%weights(i-1)
934 : enddo
935 0 : gll_edges(np+1)=1
936 0 : delta_a=gll_pts%weights
937 0 : deallocate(gll_pts%points)
938 0 : deallocate(gll_pts%weights)
939 :
940 0 : count=0
941 0 : x1=-1
942 0 : do while ( abs(x1-1) > tol )
943 : ! find point x2 closet to x1 and x2>x1:
944 0 : x2 = 1.1_r8
945 0 : do ia=2,np+1
946 0 : if (gll_edges(ia)>x1) then
947 0 : if ( ( gll_edges(ia)-x1) < (x2-x1) ) then
948 0 : x2=gll_edges(ia)
949 : endif
950 : endif
951 : enddo
952 0 : do id=2,nphys+1
953 0 : if (phys_edges(id)>x1) then
954 0 : if ( ( phys_edges(id)-x1) < (x2-x1) ) then
955 0 : x2=phys_edges(id)
956 : endif
957 : endif
958 : enddo
959 0 : print *,'x2=',x2
960 0 : if (x2>1+tol) call endrun('ERROR: did not find next intersection point')
961 0 : if (x2<=x1) call endrun('ERROR: next intersection point did not advance')
962 0 : count=count+1
963 0 : if (count>nintersect) call endrun('ERROR: search failuer: nintersect was too small')
964 0 : delta(count)=x2-x1
965 :
966 0 : found=.false.
967 0 : do ia=1,np
968 0 : if (gll_edges(ia) <= x1+tol .and. x2-tol <= gll_edges(ia+1)) then
969 0 : found=.true.
970 0 : acell(count)=ia
971 : endif
972 : enddo
973 0 : if (.not. found) call endrun('ERROR: interval search problem')
974 :
975 0 : found=.false.
976 0 : do id=1,nphys
977 0 : if (phys_edges(id) <= x1+tol .and. x2-tol <= phys_edges(id+1)) then
978 0 : found=.true.
979 0 : dcell(count)=id
980 : endif
981 : enddo
982 0 : if (.not. found) call endrun('ERROR: interval search problem')
983 0 : x1=x2
984 : enddo
985 : ! reset to actual number of intersections
986 0 : nintersect=count
987 : !OMP END MASTER
988 : !OMP BARRIER
989 : endif
990 :
991 0 : pout=0
992 0 : do in_i = 1,nintersect
993 0 : do in_j = 1,nintersect
994 0 : ia = acell(in_i)
995 0 : ja = acell(in_j)
996 0 : id = dcell(in_i)
997 0 : jd = dcell(in_j)
998 : ! mass in intersection region: value*area_intersect
999 : ! value_arrival = value*area_intersect/area_arrival
1000 0 : weight = ( delta(in_i)*delta(in_j) ) / ( delta_a(ia)*delta_a(ja))
1001 : ! accumulate contribution from each intersection region:
1002 0 : pout(ia,ja) = pout(ia,ja) + weight*pin(id+(jd-1)*nphys)
1003 : enddo
1004 : enddo
1005 :
1006 0 : end function remap_phys2gll
1007 :
1008 : !----------------------------------------------------------------
1009 :
1010 : !DIR$ ATTRIBUTES FORCEINLINE :: gradient_sphere
1011 16277274000 : subroutine gradient_sphere(s,deriv,Dinv,ds)
1012 : !
1013 : ! input s: scalar
1014 : ! output ds: spherical gradient of s, lat-lon coordinates
1015 : !
1016 :
1017 : type (derivative_t), intent(in) :: deriv
1018 : real(kind=r8), intent(in), dimension(np,np,2,2) :: Dinv
1019 : real(kind=r8), intent(in) :: s(np,np)
1020 : real(kind=r8), intent(out) :: ds(np,np,2)
1021 :
1022 : integer i
1023 : integer j
1024 : integer l
1025 :
1026 : real(kind=r8) :: dsdx00, dsdy00
1027 : real(kind=r8) :: v1(np,np),v2(np,np)
1028 :
1029 81386370000 : do j=1,np
1030 >34182*10^7 : do l=1,np
1031 : dsdx00=0.0d0
1032 : dsdy00=0.0d0
1033 : !DIR$ UNROLL(NP)
1034 >13021*10^8 : do i=1,np
1035 >10417*10^8 : dsdx00 = dsdx00 + deriv%Dvv(i,l )*s(i,j )
1036 >13021*10^8 : dsdy00 = dsdy00 + deriv%Dvv(i,l )*s(j ,i)
1037 : end do
1038 >26043*10^7 : v1(l ,j ) = dsdx00*ra
1039 >32554*10^7 : v2(j ,l ) = dsdy00*ra
1040 : end do
1041 : end do
1042 : ! convert covarient to latlon
1043 81386370000 : do j=1,np
1044 >34182*10^7 : do i=1,np
1045 >26043*10^7 : ds(i,j,1)=Dinv(i,j,1,1)*v1(i,j) + Dinv(i,j,2,1)*v2(i,j)
1046 >32554*10^7 : ds(i,j,2)=Dinv(i,j,1,2)*v1(i,j) + Dinv(i,j,2,2)*v2(i,j)
1047 : enddo
1048 : enddo
1049 :
1050 16277274000 : end subroutine gradient_sphere
1051 :
1052 :
1053 2477919600 : function curl_sphere_wk_testcov(s,deriv,elem) result(ds)
1054 : !
1055 : ! integrated-by-parts gradient, w.r.t. COVARIANT test functions
1056 : ! input s: scalar (assumed to be s*khat)
1057 : ! output ds: weak curl, lat/lon coordinates
1058 : !
1059 : ! starting with:
1060 : ! PHIcov1 = (PHI,0) covariant vector
1061 : ! PHIcov2 = (0,PHI) covariant vector
1062 : !
1063 : ! ds1 = integral[ PHIcov1 dot curl(s*khat) ]
1064 : ! ds2 = integral[ PHIcov2 dot curl(s*khat) ]
1065 : ! integrate by parts:
1066 : ! ds1 = integral[ vor(PHIcov1) * s ]
1067 : ! ds2 = integral[ vor(PHIcov1) * s ]
1068 : !
1069 : ! PHIcov1 = (PHI^mn,0)
1070 : ! PHIcov2 = (0,PHI^mn)
1071 : ! vorticity() acts on covariant vectors:
1072 : ! ds1 = sum wij g s_ij 1/g ( (PHIcov1_2)_x - (PHIcov1_1)_y )
1073 : ! = -sum wij s_ij d/dy (PHI^mn )
1074 : ! for d/dy component, only sum over i=m
1075 : ! = -sum w_mj s_mj d( PHI^n)(j)
1076 : ! j
1077 : !
1078 : ! ds2 = sum wij g s_ij 1/g ( (PHIcov2_2)_x - (PHIcov2_1)_y )
1079 : ! = +sum wij s_ij d/dx (PHI^mn )
1080 : ! for d/dx component, only sum over j=n
1081 : ! = +sum w_in s_in d( PHI^m)(i)
1082 : ! i
1083 : !
1084 : type (derivative_t), intent(in) :: deriv
1085 : type (element_t), intent(in) :: elem
1086 : real(kind=r8), intent(in) :: s(np,np)
1087 :
1088 : real(kind=r8) :: ds(np,np,2)
1089 :
1090 : integer i,j,l,m,n
1091 : real(kind=r8) :: dscontra(np,np,2)
1092 :
1093 2477919600 : dscontra=0
1094 12389598000 : do n=1,np
1095 52036311600 : do m=1,np
1096 : !DIR$ UNROLL(NP)
1097 >20814*10^7 : do j=1,np
1098 : ! phi(n)_y sum over second index, 1st index fixed at m
1099 >15858*10^7 : dscontra(m,n,1)=dscontra(m,n,1)-(elem%mp(m,j)*s(m,j)*deriv%Dvv(n,j) )*ra
1100 : ! phi(m)_x sum over first index, second index fixed at n
1101 >19823*10^7 : dscontra(m,n,2)=dscontra(m,n,2)+(elem%mp(j,n)*s(j,n)*deriv%Dvv(m,j) )*ra
1102 : enddo
1103 : enddo
1104 : enddo
1105 :
1106 : ! convert contra -> latlon
1107 12389598000 : do j=1,np
1108 52036311600 : do i=1,np
1109 39646713600 : ds(i,j,1)=(elem%D(i,j,1,1)*dscontra(i,j,1) + elem%D(i,j,1,2)*dscontra(i,j,2))
1110 49558392000 : ds(i,j,2)=(elem%D(i,j,2,1)*dscontra(i,j,1) + elem%D(i,j,2,2)*dscontra(i,j,2))
1111 : enddo
1112 : enddo
1113 2477919600 : end function curl_sphere_wk_testcov
1114 :
1115 :
1116 2477919600 : function gradient_sphere_wk_testcov(s,deriv,elem) result(ds)
1117 : !
1118 : ! integrated-by-parts gradient, w.r.t. COVARIANT test functions
1119 : ! input s: scalar
1120 : ! output ds: weak gradient, lat/lon coordinates
1121 : ! ds = - integral[ div(PHIcov) s ]
1122 : !
1123 : ! PHIcov1 = (PHI^mn,0)
1124 : ! PHIcov2 = (0,PHI^mn)
1125 : ! div() acts on contra components, so convert test function to contra:
1126 : ! PHIcontra1 = metinv PHIcov1 = (a^mn,b^mn)*PHI^mn
1127 : ! a = metinv(1,1) b=metinv(2,1)
1128 : !
1129 : ! ds1 = sum wij g s_ij 1/g ( g a PHI^mn)_x + ( g b PHI^mn)_y )
1130 : ! = sum wij s_ij ag(m,n) d/dx( PHI^mn ) + bg(m,n) d/dy( PHI^mn)
1131 : ! i,j
1132 : ! for d/dx component, only sum over j=n
1133 : ! = sum w_in s_in ag(m,n) d( PHI^m)(i)
1134 : ! i
1135 : ! for d/dy component, only sum over i=m
1136 : ! = sum w_mj s_mj bg(m,n) d( PHI^n)(j)
1137 : ! j
1138 : !
1139 : !
1140 : ! This formula is identical to gradient_sphere_wk_testcontra, except that
1141 : ! g(m,n) is replaced by a(m,n)*g(m,n)
1142 : ! and we have two terms for each componet of ds
1143 : !
1144 : !
1145 : type (derivative_t), intent(in) :: deriv
1146 : type (element_t), intent(in) :: elem
1147 : real(kind=r8), intent(in) :: s(np,np)
1148 :
1149 : real(kind=r8) :: ds(np,np,2)
1150 :
1151 : integer i,j,l,m,n
1152 : real(kind=r8) :: dscontra(np,np,2)
1153 :
1154 :
1155 2477919600 : dscontra=0
1156 12389598000 : do n=1,np
1157 52036311600 : do m=1,np
1158 : !DIR$ UNROLL(NP)
1159 >20814*10^7 : do j=1,np
1160 >31717*10^7 : dscontra(m,n,1)=dscontra(m,n,1)-(&
1161 >15858*10^7 : (elem%mp(j,n)*elem%metinv(m,n,1,1)*elem%metdet(m,n)*s(j,n)*deriv%Dvv(m,j) ) +&
1162 : (elem%mp(m,j)*elem%metinv(m,n,2,1)*elem%metdet(m,n)*s(m,j)*deriv%Dvv(n,j) ) &
1163 >47576*10^7 : ) *ra
1164 :
1165 : dscontra(m,n,2)=dscontra(m,n,2)-(&
1166 : (elem%mp(j,n)*elem%metinv(m,n,1,2)*elem%metdet(m,n)*s(j,n)*deriv%Dvv(m,j) ) +&
1167 : (elem%mp(m,j)*elem%metinv(m,n,2,2)*elem%metdet(m,n)*s(m,j)*deriv%Dvv(n,j) ) &
1168 >19823*10^7 : ) *ra
1169 : enddo
1170 : enddo
1171 : enddo
1172 : ! convert contra -> latlon
1173 12389598000 : do j=1,np
1174 52036311600 : do i=1,np
1175 39646713600 : ds(i,j,1)=(elem%D(i,j,1,1)*dscontra(i,j,1) + elem%D(i,j,1,2)*dscontra(i,j,2))
1176 49558392000 : ds(i,j,2)=(elem%D(i,j,2,1)*dscontra(i,j,1) + elem%D(i,j,2,2)*dscontra(i,j,2))
1177 : enddo
1178 : enddo
1179 :
1180 2477919600 : end function gradient_sphere_wk_testcov
1181 :
1182 :
1183 0 : function gradient_sphere_wk_testcontra(s,deriv,elem) result(ds)
1184 : !
1185 : ! integrated-by-parts gradient, w.r.t. CONTRA test functions
1186 : ! input s: scalar
1187 : ! output ds: weak gradient, lat/lon coordinates
1188 : !
1189 : ! integral[ div(phivec) s ] = sum spheremp()* divergence_sphere(phivec) *s
1190 : ! ds1 = above formual with phivec=(PHI,0) in CONTRA coordinates
1191 : ! ds2 = above formual with phivec=(0,PHI) in CONTRA coordinates
1192 : !
1193 : ! PHI = (phi,0)
1194 : ! s1 = sum w_ij s_ij g_ij 1/g_ij ( g_ij PHI^mn )x
1195 : ! = sum w_ij s_ij g_mn dx(PHI^mn)_ij
1196 : ! ij
1197 : ! because x derivative is zero for j<>n, only have to sum over j=n
1198 : ! s1(m,n) = sum w_i,n g_mn dx(PHI^m)_i,n s_i,n
1199 : ! i
1200 : !
1201 : type (derivative_t), intent(in) :: deriv
1202 : type (element_t), intent(in) :: elem
1203 : real(kind=r8), intent(in) :: s(np,np)
1204 :
1205 : real(kind=r8) :: ds(np,np,2)
1206 :
1207 : integer i,j,l,m,n
1208 : real(kind=r8) :: dscov(np,np,2)
1209 :
1210 0 : dscov=0
1211 0 : do n=1,np
1212 0 : do m=1,np
1213 : !DIR$ UNROLL(NP)
1214 0 : do j=1,np
1215 : ! phi(m)_x sum over first index, second index fixed at n
1216 0 : dscov(m,n,1)=dscov(m,n,1)-(elem%mp(j,n)*elem%metdet(m,n)*s(j,n)*deriv%Dvv(m,j) )*ra
1217 : ! phi(n)_y sum over second index, 1st index fixed at m
1218 0 : dscov(m,n,2)=dscov(m,n,2)-(elem%mp(m,j)*elem%metdet(m,n)*s(m,j)*deriv%Dvv(n,j) )*ra
1219 : enddo
1220 : enddo
1221 : enddo
1222 :
1223 : ! convert covariant -> latlon
1224 0 : ds(:,:,1)=elem%Dinv(:,:,1,1)*dscov(:,:,1) + elem%Dinv(:,:,2,1)*dscov(:,:,2)
1225 0 : ds(:,:,2)=elem%Dinv(:,:,1,2)*dscov(:,:,1) + elem%Dinv(:,:,2,2)*dscov(:,:,2)
1226 :
1227 0 : end function gradient_sphere_wk_testcontra
1228 :
1229 0 : function ugradv_sphere(u,v,deriv,elem) result(ugradv)
1230 : !
1231 : ! input: vectors u and v (latlon coordinates)
1232 : ! output: vector [ u dot grad ] v (latlon coordinates)
1233 : !
1234 : type (derivative_t), intent(in) :: deriv
1235 : type (element_t), intent(in) :: elem
1236 : real(kind=r8), intent(in) :: u(np,np,2)
1237 : real(kind=r8), intent(in) :: v(np,np,2)
1238 :
1239 : real(kind=r8) :: ugradv(np,np,2)
1240 : real(kind=r8) :: dum_cart(np,np,3)
1241 : real(kind=r8) :: temp(np,np,2)
1242 :
1243 : integer :: component
1244 :
1245 : ! latlon -> cartesian
1246 0 : do component=1,3
1247 : ! Summing along the third dimension is a sum over components for each point.
1248 : ! (This is just a faster way of doing a dot product for each grid point,
1249 : ! since reindexing the inputs to use the intrinsic effectively would be
1250 : ! just asking for trouble.)
1251 0 : dum_cart(:,:,component)=sum( elem%vec_sphere2cart(:,:,component,:)*v(:,:,:) ,3)
1252 : end do
1253 :
1254 : ! Do ugradv on the cartesian components.
1255 0 : do component=1,3
1256 : ! Dot u with the gradient of each component
1257 0 : call gradient_sphere(dum_cart(:,:,component),deriv,elem%Dinv,temp)
1258 0 : dum_cart(:,:,component) = sum( u(:,:,:) * temp,3)
1259 : enddo
1260 :
1261 : ! cartesian -> latlon
1262 0 : do component=1,2
1263 : ! vec_sphere2cart is its own pseudoinverse.
1264 0 : ugradv(:,:,component) = sum(dum_cart(:,:,:)*elem%vec_sphere2cart(:,:,:,component), 3)
1265 : end do
1266 :
1267 0 : end function ugradv_sphere
1268 :
1269 :
1270 :
1271 0 : function curl_sphere(s,deriv,elem) result(ds)
1272 : !
1273 : ! input s: scalar (assumed to be s khat)
1274 : ! output curl(s khat) vector in lat-lon coordinates
1275 : !
1276 : ! This subroutine can be used to compute divergence free velocity fields,
1277 : ! since div(ds)=0
1278 : !
1279 : ! first compute:
1280 : ! curl(s khat) = (1/jacobian) ( ds/dy, -ds/dx ) in contra-variant coordinates
1281 : ! then map to lat-lon
1282 : !
1283 : type (derivative_t), intent(in) :: deriv
1284 : type (element_t), intent(in) :: elem
1285 : real(kind=r8), intent(in) :: s(np,np)
1286 :
1287 : real(kind=r8) :: ds(np,np,2)
1288 :
1289 : integer i
1290 : integer j
1291 : integer l
1292 :
1293 : real(kind=r8) :: dsdx00
1294 : real(kind=r8) :: dsdy00
1295 : real(kind=r8) :: v1(np,np),v2(np,np)
1296 :
1297 0 : do j=1,np
1298 0 : do l=1,np
1299 : dsdx00=0.0d0
1300 : dsdy00=0.0d0
1301 : !DIR$ UNROLL(NP)
1302 0 : do i=1,np
1303 0 : dsdx00 = dsdx00 + deriv%Dvv(i,l )*s(i,j )
1304 0 : dsdy00 = dsdy00 + deriv%Dvv(i,l )*s(j ,i)
1305 : end do
1306 0 : v2(l ,j ) = -dsdx00*ra
1307 0 : v1(j ,l ) = dsdy00*ra
1308 : end do
1309 : end do
1310 : ! convert contra -> latlon *and* divide by jacobian
1311 0 : do j=1,np
1312 0 : do i=1,np
1313 0 : ds(i,j,1)=(elem%D(i,j,1,1)*v1(i,j) + elem%D(i,j,1,2)*v2(i,j))/elem%metdet(i,j)
1314 0 : ds(i,j,2)= (elem%D(i,j,2,1)*v1(i,j) + elem%D(i,j,2,2)*v2(i,j))/elem%metdet(i,j)
1315 : enddo
1316 : enddo
1317 :
1318 0 : end function curl_sphere
1319 :
1320 :
1321 : !--------------------------------------------------------------------------
1322 :
1323 :
1324 :
1325 5362718400 : subroutine divergence_sphere_wk(v,deriv,elem,div)
1326 : !
1327 : ! input: v = velocity in lat-lon coordinates
1328 : ! ouput: div(v) spherical divergence of v, integrated by parts
1329 : !
1330 : ! Computes -< grad(psi) dot v >
1331 : ! (the integrated by parts version of < psi div(v) > )
1332 : !
1333 : ! note: after DSS, divergence_sphere() and divergence_sphere_wk()
1334 : ! are identical to roundoff, as theory predicts.
1335 : !
1336 : real(kind=r8), intent(in) :: v(np,np,2) ! in lat-lon coordinates
1337 : type (derivative_t), intent(in) :: deriv
1338 : type (element_t), intent(in) :: elem
1339 : real(kind=r8),intent(out) :: div(np,np)
1340 :
1341 : ! Local
1342 :
1343 : integer i,j,m,n
1344 :
1345 : real(kind=r8) :: vtemp(np,np,2)
1346 : real(kind=r8) :: ggtemp(np,np,2)
1347 : real(kind=r8) :: gtemp(np,np,2)
1348 : real(kind=r8) :: psi(np,np)
1349 : real(kind=r8) :: xtmp
1350 :
1351 : ! latlon- > contra
1352 26813592000 : do j=1,np
1353 >11261*10^7 : do i=1,np
1354 85803494400 : vtemp(i,j,1)=(elem%Dinv(i,j,1,1)*v(i,j,1) + elem%Dinv(i,j,1,2)*v(i,j,2))
1355 >10725*10^7 : vtemp(i,j,2)=(elem%Dinv(i,j,2,1)*v(i,j,1) + elem%Dinv(i,j,2,2)*v(i,j,2))
1356 : enddo
1357 : enddo
1358 :
1359 26813592000 : do n=1,np
1360 >11261*10^7 : do m=1,np
1361 :
1362 85803494400 : div(m,n)=0
1363 : !DIR$ UNROLL(NP)
1364 >45046*10^7 : do j=1,np
1365 >34321*10^7 : div(m,n)=div(m,n)-(elem%spheremp(j,n)*vtemp(j,n,1)*deriv%Dvv(m,j) &
1366 : +elem%spheremp(m,j)*vtemp(m,j,2)*deriv%Dvv(n,j)) &
1367 >77223*10^7 : * ra
1368 : enddo
1369 :
1370 : end do
1371 : end do
1372 :
1373 5362718400 : end subroutine divergence_sphere_wk
1374 :
1375 :
1376 :
1377 0 : function element_boundary_integral(v,deriv,elem) result(result)
1378 : !
1379 : ! input: v = velocity in lat-lon coordinates
1380 : ! ouput: result(i,j) = contour integral of PHI_ij * v dot normal
1381 : ! where PHI_ij = cardinal function at i,j GLL point
1382 : !
1383 : ! this routine is used just to check spectral element integration by parts identities
1384 : !
1385 : real(kind=r8), intent(in) :: v(np,np,2) ! in lat-lon coordinates
1386 : type (derivative_t), intent(in) :: deriv
1387 : type (element_t), intent(in) :: elem
1388 : real(kind=r8) :: result(np,np)
1389 :
1390 : ! Local
1391 : real(kind=r8) :: ucontra(np,np,2) ! in lat-lon coordinates
1392 : integer i,j
1393 :
1394 : ! latlon->contra
1395 0 : do j=1,np
1396 0 : do i=1,np
1397 0 : ucontra(i,j,1)=(elem%Dinv(i,j,1,1)*v(i,j,1) + elem%Dinv(i,j,1,2)*v(i,j,2))
1398 0 : ucontra(i,j,2)=(elem%Dinv(i,j,2,1)*v(i,j,1) + elem%Dinv(i,j,2,2)*v(i,j,2))
1399 : enddo
1400 : enddo
1401 :
1402 : ! note: GLL weights weight(i) = Mvv_twt(i,i)
1403 0 : result=0
1404 0 : j=1
1405 0 : do i=1,np
1406 0 : result(i,j)=result(i,j)-deriv%Mvv_twt(i,i)*elem%metdet(i,j)*ucontra(i,j,2)*ra
1407 : enddo
1408 :
1409 0 : j=np
1410 0 : do i=1,np
1411 0 : result(i,j)=result(i,j)+deriv%Mvv_twt(i,i)*elem%metdet(i,j)*ucontra(i,j,2)*ra
1412 : enddo
1413 :
1414 0 : i=1
1415 0 : do j=1,np
1416 0 : result(i,j)=result(i,j)-deriv%Mvv_twt(j,j)*elem%metdet(i,j)*ucontra(i,j,1)*ra
1417 : enddo
1418 :
1419 0 : i=np
1420 0 : do j=1,np
1421 0 : result(i,j)=result(i,j)+deriv%Mvv_twt(j,j)*elem%metdet(i,j)*ucontra(i,j,1)*ra
1422 : enddo
1423 0 : end function element_boundary_integral
1424 :
1425 :
1426 :
1427 0 : function edge_flux_u_cg( v,p,pedges, deriv, elem, u_is_contra) result(result)
1428 : !
1429 : !
1430 : ! input: v = velocity in contra or lat-lon coordinates (CONTINUIOUS)
1431 : ! p = scalar on this element
1432 : ! pedges = scalar edge data from neighbor elements
1433 : !
1434 : ! ouput: result(i,j) = contour integral of PHI_ij * pstar * v dot normal
1435 : ! where PHI_ij = cardinal function at i,j GLL point
1436 : ! pstar = centered or other flux
1437 : !
1438 : real(kind=r8), intent(in) :: v(np,np,2)
1439 : real(kind=r8), intent(in) :: p(np,np)
1440 : real(kind=r8), intent(in) :: pedges(0:np+1,0:np+1)
1441 : type (derivative_t), intent(in) :: deriv
1442 : type (element_t), intent(in) :: elem
1443 : real(kind=r8) :: result(np,np)
1444 : logical :: u_is_contra
1445 :
1446 : ! Local
1447 : real(kind=r8) :: ucontra(np,np,2) ! in lat-lon coordinates
1448 : real(kind=r8) :: flux,pstar
1449 : integer i,j
1450 :
1451 :
1452 0 : result=0
1453 :
1454 :
1455 0 : if (u_is_contra) then
1456 0 : ucontra=v
1457 : else
1458 : ! latlon->contra
1459 0 : do j=1,np
1460 0 : do i=1,np
1461 0 : ucontra(i,j,1)=(elem%Dinv(i,j,1,1)*v(i,j,1) + elem%Dinv(i,j,1,2)*v(i,j,2))
1462 0 : ucontra(i,j,2)=(elem%Dinv(i,j,2,1)*v(i,j,1) + elem%Dinv(i,j,2,2)*v(i,j,2))
1463 : enddo
1464 : enddo
1465 : endif
1466 : ! upwind
1467 0 : do i=1,np
1468 0 : j=1
1469 0 : pstar=p(i,j)
1470 0 : if (ucontra(i,j,2)>0) pstar=pedges(i,0)
1471 0 : flux = -pstar*ucontra(i,j,2)*( deriv%Mvv_twt(i,i)*elem%metdet(i,j)*ra)
1472 0 : result(i,j)=result(i,j)+flux
1473 :
1474 0 : j=np
1475 0 : pstar=p(i,j)
1476 0 : if (ucontra(i,j,2)<0) pstar=pedges(i,np+1)
1477 0 : flux = pstar*ucontra(i,j,2)* ( deriv%Mvv_twt(i,i)*elem%metdet(i,j)*ra)
1478 0 : result(i,j)=result(i,j)+flux
1479 : enddo
1480 :
1481 0 : do j=1,np
1482 0 : i=1
1483 0 : pstar=p(i,j)
1484 0 : if (ucontra(i,j,1)>0) pstar=pedges(0,j)
1485 0 : flux = -pstar*ucontra(i,j,1)* ( deriv%Mvv_twt(j,j)*elem%metdet(i,j)*ra)
1486 0 : result(i,j)=result(i,j)+flux
1487 :
1488 0 : i=np
1489 0 : pstar=p(i,j)
1490 0 : if (ucontra(i,j,1)<0) pstar=pedges(np+1,j)
1491 0 : flux = pstar*ucontra(i,j,1)* ( deriv%Mvv_twt(j,j)*elem%metdet(i,j)*ra)
1492 0 : result(i,j)=result(i,j)+flux
1493 : end do
1494 :
1495 0 : end function edge_flux_u_cg
1496 :
1497 :
1498 : !DIR$ ATTRIBUTES FORCEINLINE :: vorticity_sphere
1499 4503891600 : subroutine vorticity_sphere(v,deriv,elem,vort)
1500 : !
1501 : ! input: v = velocity in lat-lon coordinates
1502 : ! ouput: spherical vorticity of v
1503 : !
1504 :
1505 : type (derivative_t), intent(in) :: deriv
1506 : type (element_t), intent(in) :: elem
1507 : real(kind=r8), intent(in) :: v(np,np,2)
1508 :
1509 : real(kind=r8), intent(out) :: vort(np,np)
1510 :
1511 : integer i
1512 : integer j
1513 : integer l
1514 :
1515 : real(kind=r8) :: dvdx00,dudy00
1516 : real(kind=r8) :: vco(np,np,2)
1517 : real(kind=r8) :: vtemp(np,np)
1518 :
1519 : ! convert to covariant form
1520 22519458000 : do j=1,np
1521 94581723600 : do i=1,np
1522 72062265600 : vco(i,j,1)=(elem%D(i,j,1,1)*v(i,j,1) + elem%D(i,j,2,1)*v(i,j,2))
1523 90077832000 : vco(i,j,2)=(elem%D(i,j,1,2)*v(i,j,1) + elem%D(i,j,2,2)*v(i,j,2))
1524 : enddo
1525 : enddo
1526 :
1527 22519458000 : do j=1,np
1528 94581723600 : do l=1,np
1529 :
1530 : dudy00=0.0d0
1531 : dvdx00=0.0d0
1532 :
1533 : !DIR$ UNROLL(NP)
1534 >36031*10^7 : do i=1,np
1535 >28824*10^7 : dvdx00 = dvdx00 + deriv%Dvv(i,l )*vco(i,j ,2)
1536 >36031*10^7 : dudy00 = dudy00 + deriv%Dvv(i,l )*vco(j ,i,1)
1537 : enddo
1538 :
1539 72062265600 : vort(l ,j ) = dvdx00
1540 90077832000 : vtemp(j ,l ) = dudy00
1541 : enddo
1542 : enddo
1543 :
1544 22519458000 : do j=1,np
1545 94581723600 : do i=1,np
1546 90077832000 : vort(i,j)=(vort(i,j)-vtemp(i,j))*(elem%rmetdet(i,j)*ra)
1547 : end do
1548 : end do
1549 :
1550 4503891600 : end subroutine vorticity_sphere
1551 :
1552 0 : function vorticity_sphere_diag(v,deriv,elem) result(vort)
1553 : !
1554 : ! input: v = velocity in lat-lon coordinates
1555 : ! ouput: diagonal component of spherical vorticity of v
1556 : !
1557 :
1558 : type (derivative_t), intent(in) :: deriv
1559 : type (element_t), intent(in) :: elem
1560 : real(kind=r8), intent(in) :: v(np,np,2)
1561 :
1562 : real(kind=r8) :: vort(np,np)
1563 :
1564 : integer i
1565 : integer j
1566 : integer l
1567 :
1568 : real(kind=r8) :: dvdx00,dudy00
1569 : real(kind=r8) :: vco(np,np,2)
1570 : real(kind=r8) :: vtemp(np,np)
1571 : real(kind=r8) :: rdx
1572 : real(kind=r8) :: rdy
1573 :
1574 : ! convert to covariant form
1575 :
1576 0 : do j=1,np
1577 0 : do i=1,np
1578 0 : vco(i,j,1)=(elem%D(i,j,1,1)*v(i,j,1) + elem%D(i,j,2,1)*v(i,j,2))
1579 0 : vco(i,j,2)=(elem%D(i,j,1,2)*v(i,j,1) + elem%D(i,j,2,2)*v(i,j,2))
1580 : enddo
1581 : enddo
1582 :
1583 :
1584 0 : do j=1,np
1585 0 : do l=1,np
1586 : dudy00=0.0d0
1587 : dvdx00=0.0d0
1588 : !DIR$ UNROLL(NP)
1589 0 : do i=1,np
1590 0 : dvdx00 = dvdx00 + deriv%Dvv_diag(i,l)*vco(i,j ,2)
1591 0 : dudy00 = dudy00 + deriv%Dvv_diag(i,l)*vco(j ,i,1)
1592 : enddo
1593 0 : vort(l ,j) = dvdx00
1594 0 : vtemp(j ,l) = dudy00
1595 : enddo
1596 : enddo
1597 :
1598 0 : do j=1,np
1599 0 : do i=1,np
1600 0 : vort(i,j)=(vort(i,j)-vtemp(i,j))*(elem%rmetdet(i,j)*ra)
1601 : end do
1602 : end do
1603 :
1604 0 : end function vorticity_sphere_diag
1605 :
1606 : !DIR$ ATTRIBUTES FORCEINLINE :: divergence_sphere
1607 6597676800 : subroutine divergence_sphere(v,deriv,elem,div)
1608 : !
1609 : ! input: v = velocity in lat-lon coordinates
1610 : ! ouput: div(v) spherical divergence of v
1611 : !
1612 :
1613 :
1614 : real(kind=r8), intent(in) :: v(np,np,2) ! in lat-lon coordinates
1615 : type (derivative_t), intent(in) :: deriv
1616 : type (element_t), intent(in) :: elem
1617 : real(kind=r8), intent(out) :: div(np,np)
1618 :
1619 : ! Local
1620 :
1621 : integer i
1622 : integer j
1623 : integer l
1624 :
1625 : real(kind=r8) :: dudx00
1626 : real(kind=r8) :: dvdy00
1627 : real(kind=r8) :: gv(np,np,2),vvtemp(np,np)
1628 :
1629 : ! convert to contra variant form and multiply by g
1630 32988384000 : do j=1,np
1631 >13855*10^7 : do i=1,np
1632 >10556*10^7 : gv(i,j,1)=elem%metdet(i,j)*(elem%Dinv(i,j,1,1)*v(i,j,1) + elem%Dinv(i,j,1,2)*v(i,j,2))
1633 >13195*10^7 : gv(i,j,2)=elem%metdet(i,j)*(elem%Dinv(i,j,2,1)*v(i,j,1) + elem%Dinv(i,j,2,2)*v(i,j,2))
1634 : enddo
1635 : enddo
1636 :
1637 : ! compute d/dx and d/dy
1638 32988384000 : do j=1,np
1639 >13855*10^7 : do l=1,np
1640 : dudx00=0.0d0
1641 : dvdy00=0.0d0
1642 : !DIR$ UNROLL(NP)
1643 >52781*10^7 : do i=1,np
1644 >42225*10^7 : dudx00 = dudx00 + deriv%Dvv(i,l )*gv(i,j ,1)
1645 >52781*10^7 : dvdy00 = dvdy00 + deriv%Dvv(i,l )*gv(j ,i,2)
1646 : end do
1647 >10556*10^7 : div(l ,j ) = dudx00
1648 >13195*10^7 : vvtemp(j ,l ) = dvdy00
1649 : end do
1650 : end do
1651 :
1652 32988384000 : do j=1,np
1653 >13855*10^7 : do i=1,np
1654 >13195*10^7 : div(i,j)=(div(i,j)+vvtemp(i,j))*(elem%rmetdet(i,j)*ra)
1655 : end do
1656 : end do
1657 :
1658 6597676800 : end subroutine divergence_sphere
1659 :
1660 :
1661 : !DIR$ ATTRIBUTES FORCEINLINE :: laplace_sphere_wk
1662 5362718400 : subroutine laplace_sphere_wk(s,deriv,elem,laplace,var_coef,mol_nu)
1663 : !
1664 : ! input: s = scalar
1665 : ! ouput: -< grad(PHI), grad(s) > = weak divergence of grad(s)
1666 : ! note: for this form of the operator, grad(s) does not need to be made C0
1667 : !
1668 : real(kind=r8), intent(in) :: s(np,np)
1669 : type (derivative_t), intent(in) :: deriv
1670 : type (element_t), intent(in) :: elem
1671 : real(kind=r8) :: laplace(np,np)
1672 : logical, intent(in) :: var_coef
1673 : real(kind=r8), intent(in), optional :: mol_nu(np,np) !variable nu (e.g. molecular diffusion)
1674 :
1675 : real(kind=r8) :: laplace2(np,np)
1676 : integer i,j
1677 :
1678 : ! Local
1679 : real(kind=r8) :: grads(np,np,2), oldgrads(np,np,2)
1680 :
1681 5362718400 : call gradient_sphere(s,deriv,elem%Dinv,grads)
1682 : !
1683 : ! molecular diffusion coefficient
1684 : !
1685 5362718400 : if (present(mol_nu)) then
1686 0 : if (var_coef) &
1687 0 : call endrun('ERROR: this option is only for non-grid dependent variable viscosity')
1688 0 : grads(:,:,1) = grads(:,:,1)*mol_nu(:,:)
1689 0 : grads(:,:,2) = grads(:,:,2)*mol_nu(:,:)
1690 : end if
1691 :
1692 5362718400 : if (var_coef) then
1693 5269212000 : if (hypervis_power/=0 ) then
1694 : ! scalar viscosity with variable coefficient
1695 0 : grads(:,:,1) = grads(:,:,1)*elem%variable_hyperviscosity(:,:)
1696 0 : grads(:,:,2) = grads(:,:,2)*elem%variable_hyperviscosity(:,:)
1697 5269212000 : else if (hypervis_scaling /=0 ) then
1698 : ! tensor hv, (3)
1699 0 : oldgrads=grads
1700 0 : do j=1,np
1701 0 : do i=1,np
1702 0 : grads(i,j,1) = oldgrads(i,j,1)*elem%tensorVisc(i,j,1,1) + &
1703 0 : oldgrads(i,j,2)*elem%tensorVisc(i,j,1,2)
1704 : grads(i,j,2) = oldgrads(i,j,1)*elem%tensorVisc(i,j,2,1) + &
1705 0 : oldgrads(i,j,2)*elem%tensorVisc(i,j,2,2)
1706 : end do
1707 : end do
1708 : else
1709 : ! do nothing: constant coefficient viscosity
1710 : endif
1711 : endif
1712 :
1713 : ! note: divergnece_sphere and divergence_sphere_wk are identical *after* bndry_exchange
1714 : ! if input is C_0. Here input is not C_0, so we should use divergence_sphere_wk().
1715 : ! laplace=divergence_sphere_wk(grads,deriv,elem)
1716 5362718400 : call divergence_sphere_wk(grads,deriv,elem,laplace)
1717 :
1718 5362718400 : end subroutine laplace_sphere_wk
1719 :
1720 : !DIR$ ATTRIBUTES FORCEINLINE :: vlaplace_sphere_wk
1721 2477919600 : subroutine vlaplace_sphere_wk(v,deriv,elem,undamprrcart,laplace,var_coef,nu_ratio)
1722 : !
1723 : ! input: v = vector in lat-lon coordinates
1724 : ! ouput: weak laplacian of v, in lat-lon coordinates
1725 : !
1726 : ! logic:
1727 : ! tensorHV: requires cartesian
1728 : ! nu_div/=nu: requires contra formulatino
1729 : !
1730 : ! One combination NOT supported: tensorHV and nu_div/=nu then abort
1731 : !
1732 : real(kind=r8), intent(in) :: v(np,np,2)
1733 : type (derivative_t), intent(in) :: deriv
1734 : type (element_t), intent(in) :: elem
1735 : logical, intent(in) :: undamprrcart
1736 : real(kind=r8), intent(out) :: laplace(np,np,2)
1737 : logical, optional, intent(in) :: var_coef
1738 : real(kind=r8), optional, intent(in) :: nu_ratio
1739 :
1740 :
1741 2477919600 : if (hypervis_scaling/=0 .and. var_coef) then
1742 : ! tensorHV is turned on - requires cartesian formulation
1743 0 : if (present(nu_ratio)) then
1744 0 : if (nu_ratio /= 1._r8) then
1745 0 : call endrun('ERROR: tensorHV can not be used with nu_div/=nu')
1746 : endif
1747 : endif
1748 0 : laplace=vlaplace_sphere_wk_cartesian(v,deriv,elem,var_coef,undamprrcart)
1749 : else
1750 : ! all other cases, use contra formulation:
1751 2477919600 : laplace=vlaplace_sphere_wk_contra(v,deriv,elem,var_coef,undamprrcart,nu_ratio)
1752 : endif
1753 :
1754 2477919600 : end subroutine vlaplace_sphere_wk
1755 : !
1756 : ! version of vlaplace_sphere_wk for molecular diffusion
1757 : !
1758 0 : subroutine vlaplace_sphere_wk_mol(v,deriv,elem,undamprrcart,mol_nu,laplace)
1759 : !
1760 : ! input: v = vector in lat-lon coordinates
1761 : ! ouput: weak laplacian of v, in lat-lon coordinates
1762 : !
1763 : real(kind=r8), intent(in) :: v(np,np,2)
1764 : type (derivative_t), intent(in):: deriv
1765 : type (element_t), intent(in) :: elem
1766 : logical, intent(in) :: undamprrcart
1767 : real(kind=r8), intent(in) :: mol_nu(np,np)
1768 : real(kind=r8), intent(out) :: laplace(np,np,2)
1769 :
1770 : real(kind=r8) :: vor(np,np),div(np,np)
1771 :
1772 : integer :: n,m
1773 :
1774 0 : call divergence_sphere(v,deriv,elem,div)
1775 0 : call vorticity_sphere(v,deriv,elem,vor)
1776 :
1777 0 : div = div*mol_nu(:,:)
1778 0 : vor = vor*mol_nu(:,:)
1779 :
1780 : laplace = gradient_sphere_wk_testcov(div,deriv,elem) - &
1781 0 : curl_sphere_wk_testcov(vor,deriv,elem)
1782 :
1783 0 : if (undamprrcart) then
1784 0 : do n=1,np
1785 0 : do m=1,np
1786 : ! add in correction so we dont damp rigid rotation
1787 0 : laplace(m,n,1)=laplace(m,n,1) + 2*elem%spheremp(m,n)*v(m,n,1)*(ra**2)
1788 0 : laplace(m,n,2)=laplace(m,n,2) + 2*elem%spheremp(m,n)*v(m,n,2)*(ra**2)
1789 : enddo
1790 : enddo
1791 : end if
1792 :
1793 :
1794 0 : end subroutine vlaplace_sphere_wk_mol
1795 :
1796 :
1797 :
1798 0 : function vlaplace_sphere_wk_cartesian(v,deriv,elem,var_coef,undamprrcart) result(laplace)
1799 : !
1800 : ! input: v = vector in lat-lon coordinates
1801 : ! ouput: weak laplacian of v, in lat-lon coordinates
1802 :
1803 : real(kind=r8), intent(in) :: v(np,np,2)
1804 : type (derivative_t), intent(in) :: deriv
1805 : type (element_t), intent(in) :: elem
1806 : logical, intent(in) :: var_coef
1807 : logical, intent(in) :: undamprrcart
1808 : real(kind=r8) :: laplace(np,np,2)
1809 : ! Local
1810 :
1811 : integer component
1812 : real(kind=r8) :: dum_cart(np,np,3)
1813 : real(kind=r8) :: dum_cart2(np,np)
1814 :
1815 :
1816 : ! latlon -> cartesian
1817 0 : do component=1,3
1818 0 : dum_cart2(:,:) = elem%vec_sphere2cart(:,:,component,1)*v(:,:,1) + &
1819 0 : elem%vec_sphere2cart(:,:,component,2)*v(:,:,2)
1820 : ! Do laplace on cartesian comps
1821 0 : call laplace_sphere_wk(dum_cart2,deriv,elem,dum_cart(:,:,component),var_coef)
1822 : enddo
1823 :
1824 : ! cartesian -> latlon
1825 0 : do component=1,2
1826 : ! vec_sphere2cart is its own pseudoinverse.
1827 0 : laplace(:,:,component) = dum_cart(:,:,1)*elem%vec_sphere2cart(:,:,1,component) + &
1828 : dum_cart(:,:,2)*elem%vec_sphere2cart(:,:,2,component) + &
1829 0 : dum_cart(:,:,3)*elem%vec_sphere2cart(:,:,3,component)
1830 : end do
1831 :
1832 0 : if (undamprrcart) then
1833 : ! add in correction so we dont damp rigid rotation
1834 0 : laplace(:,:,1)=laplace(:,:,1) + 2*elem%spheremp(:,:)*v(:,:,1)*(ra**2)
1835 0 : laplace(:,:,2)=laplace(:,:,2) + 2*elem%spheremp(:,:)*v(:,:,2)*(ra**2)
1836 : end if
1837 :
1838 0 : end function vlaplace_sphere_wk_cartesian
1839 :
1840 :
1841 :
1842 2477919600 : function vlaplace_sphere_wk_contra(v,deriv,elem,var_coef,undamprrcart,nu_ratio) result(laplace)
1843 : !
1844 : ! input: v = vector in lat-lon coordinates
1845 : ! ouput: weak laplacian of v, in lat-lon coordinates
1846 : !
1847 : ! du/dt = laplace(u) = grad(div) - curl(vor)
1848 : ! < PHI du/dt > = < PHI laplace(u) > PHI = covariant, u = contravariant
1849 : ! = < PHI grad(div) > - < PHI curl(vor) >
1850 : ! = grad_wk(div) - curl_wk(vor)
1851 : !
1852 : real(kind=r8), intent(in) :: v(np,np,2)
1853 : logical, intent(in) :: var_coef
1854 : type (derivative_t), intent(in) :: deriv
1855 : type (element_t), intent(in) :: elem
1856 : logical, intent(in) :: undamprrcart
1857 : real(kind=r8), optional, intent(in) :: nu_ratio
1858 :
1859 : real(kind=r8) :: laplace(np,np,2)
1860 :
1861 : ! Local
1862 :
1863 : integer i,j,l,m,n
1864 : real(kind=r8) :: vor(np,np),div(np,np)
1865 : real(kind=r8) :: v1,v2,div1,div2,vor1,vor2,phi_x,phi_y
1866 :
1867 2477919600 : call divergence_sphere(v,deriv,elem,div)
1868 2477919600 : call vorticity_sphere(v,deriv,elem,vor)
1869 :
1870 2477919600 : if (var_coef .and. hypervis_power/=0 ) then
1871 : ! scalar viscosity with variable coefficient
1872 0 : div = div*elem%variable_hyperviscosity(:,:)
1873 0 : vor = vor*elem%variable_hyperviscosity(:,:)
1874 : endif
1875 :
1876 52036311600 : if (present(nu_ratio)) div = nu_ratio*div
1877 :
1878 : laplace = gradient_sphere_wk_testcov(div,deriv,elem) - &
1879 >10655*10^7 : curl_sphere_wk_testcov(vor,deriv,elem)
1880 :
1881 2477919600 : if (undamprrcart) then
1882 12389598000 : do n=1,np
1883 52036311600 : do m=1,np
1884 : ! add in correction so we dont damp rigid rotation
1885 39646713600 : laplace(m,n,1)=laplace(m,n,1) + 2*elem%spheremp(m,n)*v(m,n,1)*(ra**2)
1886 49558392000 : laplace(m,n,2)=laplace(m,n,2) + 2*elem%spheremp(m,n)*v(m,n,2)*(ra**2)
1887 : enddo
1888 : enddo
1889 : end if
1890 2477919600 : end function vlaplace_sphere_wk_contra
1891 :
1892 0 : function gll_to_dgmodal(p,deriv) result(phat)
1893 : !
1894 : ! input: v = velocity in lat-lon coordinates
1895 : ! ouput: phat = Legendre coefficients
1896 : !
1897 : ! Computes < g dot p > = SUM g(i,j) p(i,j) w(i) w(j)
1898 : ! (the quadrature approximation on the *reference element* of the integral of p against
1899 : ! all Legendre polynomials up to degree npdg
1900 : !
1901 : ! for npdg < np, this routine gives the (exact) modal expansion of p/spheremp()
1902 : !
1903 : real(kind=r8), intent(in) :: p(np,np)
1904 : type (derivative_t), intent(in) :: deriv
1905 : real(kind=r8) :: phat(npdg,npdg)
1906 :
1907 : ! Local
1908 : integer i,j,m,n
1909 0 : real(kind=r8) :: A(np,npdg)
1910 0 : A=0
1911 0 : phat=0
1912 :
1913 : ! N^3 tensor product formulation:
1914 0 : do m=1,npdg
1915 0 : do j=1,np
1916 : !DIR$ UNROLL(NP)
1917 0 : do i=1,np
1918 0 : A(j,m)=A(j,m)+( p(i,j)*deriv%Mvv_twt(i,i)*deriv%Mvv_twt(j,j) )*deriv%legdg(m,i)
1919 : enddo
1920 : enddo
1921 : enddo
1922 :
1923 0 : do n=1,npdg
1924 0 : do m=1,npdg
1925 : !DIR$ UNROLL(NP)
1926 0 : do j=1,np
1927 0 : phat(m,n)=phat(m,n)+A(j,m)*deriv%legdg(n,j)
1928 : enddo
1929 : enddo
1930 : enddo
1931 :
1932 0 : end function
1933 :
1934 0 : function dgmodal_to_gll(phat,deriv) result(p)
1935 : !
1936 : ! input: phat = coefficients of Legendre expansion
1937 : ! ouput: p = sum expansion to evaluate phat at GLL points
1938 : !
1939 : real(kind=r8) :: p(np,np)
1940 : type (derivative_t), intent(in) :: deriv
1941 : real(kind=r8) :: phat(npdg,npdg)
1942 : ! Local
1943 : integer i,j,m,n
1944 0 : real(kind=r8) :: A(npdg,np)
1945 :
1946 0 : p(:,:)=0
1947 : ! tensor product version
1948 0 : A=0
1949 0 : do i=1,np
1950 0 : do n=1,npdg
1951 0 : do m=1,npdg
1952 0 : A(n,i)=A(n,i)+phat(m,n)*deriv%legdg(m,i)
1953 : enddo
1954 : enddo
1955 : enddo
1956 0 : do j=1,np
1957 0 : do i=1,np
1958 0 : do n=1,npdg
1959 0 : p(i,j) = p(i,j)+A(n,i)*deriv%legdg(n,j)
1960 : enddo
1961 : enddo
1962 : enddo
1963 :
1964 0 : end function dgmodal_to_gll
1965 :
1966 2072725200 : subroutine subcell_dss_fluxes(dss, p, n, metdet, C, fluxes)
1967 :
1968 : integer , intent(in) :: p
1969 : integer , intent(in) :: n
1970 : real (kind=r8), intent(in) :: dss (p,p)
1971 : real (kind=r8), intent(in) :: metdet (p,p)
1972 : real (kind=r8), intent(in) :: C (2,2,2)
1973 :
1974 : real (kind=r8) :: fluxes (n,n,4)
1975 :
1976 4145450400 : real (kind=r8) :: Bp(p,p)
1977 4145450400 : real (kind=r8) :: Tp(p,p)
1978 4145450400 : real (kind=r8) :: Lp(p,p)
1979 4145450400 : real (kind=r8) :: Rp(p,p)
1980 :
1981 4145450400 : real (kind=r8) :: B(n,n)
1982 4145450400 : real (kind=r8) :: T(n,n)
1983 4145450400 : real (kind=r8) :: L(n,n)
1984 2072725200 : real (kind=r8) :: R(n,n)
1985 :
1986 : integer :: i,j
1987 :
1988 >10985*10^7 : fluxes = 0
1989 :
1990 43527229200 : Bp = 0
1991 43527229200 : Tp = 0
1992 43527229200 : Rp = 0
1993 43527229200 : Lp = 0
1994 :
1995 12436351200 : Bp(:,1) = dss(:,1) ! bottom
1996 10363626000 : Tp(:,p) = dss(:,p) ! top
1997 10363626000 : Rp(p,:) = dss(p,:) ! right
1998 10363626000 : Lp(1,:) = dss(1,:) ! left
1999 :
2000 2072725200 : Bp(1,1) = C(1,1,2)
2001 2072725200 : Lp(1,1) = C(1,1,1)
2002 2072725200 : Bp(p,1) = C(2,1,2)
2003 2072725200 : Rp(p,1) = C(2,1,1)
2004 :
2005 2072725200 : Tp(1,p) = C(1,2,2)
2006 2072725200 : Lp(1,p) = C(1,2,1)
2007 2072725200 : Tp(p,p) = C(2,2,2)
2008 2072725200 : Rp(p,p) = C(2,2,1)
2009 :
2010 :
2011 2072725200 : call subcell_integration(Bp, p, n, metdet,B)
2012 2072725200 : call subcell_integration(Tp, p, n, metdet,T)
2013 2072725200 : call subcell_integration(Lp, p, n, metdet,L)
2014 2072725200 : call subcell_integration(Rp, p, n, metdet,R)
2015 :
2016 8290900800 : do i = 1,n
2017 26945427600 : do j = 1,n
2018 18654526800 : if (1<j) T(i,j) = T(i,j) + T(i,j-1)
2019 24872702400 : if (1<i) R(i,j) = R(i,j) + R(i-1,j)
2020 : end do
2021 : end do
2022 :
2023 8290900800 : do i = n,1,-1
2024 26945427600 : do j = n,1,-1
2025 18654526800 : if (j<n) B(i,j) = B(i,j) + B(i,j+1)
2026 24872702400 : if (i<n) L(i,j) = L(i,j) + L(i+1,j)
2027 : end do
2028 : end do
2029 :
2030 8290900800 : do i = 1,n
2031 26945427600 : do j = 1,n
2032 18654526800 : if (1==j) fluxes(i,j,1) = B(i,j)
2033 18654526800 : if (n==i) fluxes(i,j,2) = R(i,j)
2034 18654526800 : if (j==n) fluxes(i,j,3) = T(i,j)
2035 18654526800 : if (1==i) fluxes(i,j,4) = L(i,j)
2036 :
2037 18654526800 : if (1< j) fluxes(i,j,1) = B(i,j) - T(i,j-1)
2038 18654526800 : if (i< n) fluxes(i,j,2) = R(i,j) - L(i+1,j)
2039 18654526800 : if (j< n) fluxes(i,j,3) = T(i,j) - B(i,j+1)
2040 24872702400 : if (1< i) fluxes(i,j,4) = L(i,j) - R(i-1,j)
2041 : end do
2042 : end do
2043 :
2044 2072725200 : end subroutine subcell_dss_fluxes
2045 :
2046 810388800 : subroutine subcell_div_fluxes(u, p, n, metdet,fluxes)
2047 :
2048 : implicit none
2049 :
2050 : integer , intent(in) :: p
2051 : integer , intent(in) :: n
2052 : real (kind=r8), intent(in) :: u(p,p,2)
2053 : real (kind=r8), intent(in) :: metdet(p,p)
2054 :
2055 1620777600 : real (kind=r8) :: v(p,p,2)
2056 : real (kind=r8) :: fluxes(n,n,4)
2057 1620777600 : real (kind=r8) :: tb(n,p)
2058 1620777600 : real (kind=r8) :: lr(p,n)
2059 1620777600 : real (kind=r8) :: flux_l(n,n)
2060 1620777600 : real (kind=r8) :: flux_r(n,n)
2061 1620777600 : real (kind=r8) :: flux_b(n,n)
2062 1620777600 : real (kind=r8) :: flux_t(n,n)
2063 :
2064 : integer i,j
2065 :
2066 : if (.not.ALLOCATED(integration_matrix) .or. &
2067 810388800 : SIZE(integration_matrix,1).ne.n .or. &
2068 : SIZE(integration_matrix,2).ne.p) then
2069 0 : call endrun( 'FATAL ERROR: allocate_subcell_integration_matrix not called')
2070 : end if
2071 :
2072 17018164800 : v(:,:,1) = u(:,:,1)*metdet(:,:)
2073 17018164800 : v(:,:,2) = u(:,:,2)*metdet(:,:)
2074 :
2075 810388800 : tb = MATMUL(integration_matrix, v(:,:,2))
2076 810388800 : flux_b(:,:) = MATMUL(tb,TRANSPOSE(boundary_interp_matrix(:,1,:)))
2077 810388800 : flux_t(:,:) = MATMUL(tb,TRANSPOSE(boundary_interp_matrix(:,2,:)))
2078 :
2079 810388800 : lr = MATMUL(v(:,:,1),TRANSPOSE(integration_matrix))
2080 810388800 : flux_l(:,:) = MATMUL(boundary_interp_matrix(:,1,:),lr)
2081 810388800 : flux_r(:,:) = MATMUL(boundary_interp_matrix(:,2,:),lr)
2082 :
2083 10535054400 : fluxes(:,:,1) = -flux_b(:,:)*ra
2084 10535054400 : fluxes(:,:,2) = flux_r(:,:)*ra
2085 10535054400 : fluxes(:,:,3) = flux_t(:,:)*ra
2086 10535054400 : fluxes(:,:,4) = -flux_l(:,:)*ra
2087 :
2088 810388800 : end subroutine subcell_div_fluxes
2089 :
2090 1262336400 : subroutine subcell_Laplace_fluxes(u, deriv, elem, p, n, fluxes)
2091 :
2092 : implicit none
2093 :
2094 : integer , intent(in) :: p
2095 : integer , intent(in) :: n
2096 : type (derivative_t) , intent(in) :: deriv
2097 : type (element_t) , intent(in) :: elem
2098 : real (kind=r8), intent(in) :: u(p,p)
2099 :
2100 2524672800 : real (kind=r8) :: g(p,p,2)
2101 2524672800 : real (kind=r8) :: v(p,p,2)
2102 2524672800 : real (kind=r8) :: div(p,p,2)
2103 2524672800 : real (kind=r8) :: sub_int(n,n,2)
2104 : real (kind=r8) :: fluxes(n,n,4)
2105 :
2106 : integer i,j
2107 :
2108 1262336400 : call gradient_sphere(u,deriv,elem%Dinv,g)
2109 :
2110 26509064400 : v(:,:,1) = elem%Dinv(:,:,1,1)*g(:,:,1) + elem%Dinv(:,:,1,2)*g(:,:,2)
2111 26509064400 : v(:,:,2) = elem%Dinv(:,:,2,1)*g(:,:,1) + elem%Dinv(:,:,2,2)*g(:,:,2)
2112 :
2113 6311682000 : do j=1,p
2114 26509064400 : do i=1,p
2115 >10098*10^7 : div(i,j,1) = -SUM(elem%spheremp(:,j)*v(:,j,1)*deriv%Dvv(i,:))
2116 >10603*10^7 : div(i,j,2) = -SUM(elem%spheremp(i,:)*v(i,:,2)*deriv%Dvv(j,:))
2117 : end do
2118 : end do
2119 54280465200 : div = div * ra
2120 :
2121 26509064400 : div(:,:,1) = div(:,:,1) / elem%spheremp(:,:)
2122 26509064400 : div(:,:,2) = div(:,:,2) / elem%spheremp(:,:)
2123 :
2124 1262336400 : call subcell_integration(div(:,:,1), p, n, elem%metdet, sub_int(:,:,1))
2125 1262336400 : call subcell_integration(div(:,:,2), p, n, elem%metdet, sub_int(:,:,2))
2126 :
2127 5049345600 : do i=1,n
2128 12623364000 : do j=2,n
2129 7574018400 : sub_int(j,i,1) = sub_int(j,i,1) + sub_int(j-1,i,1)
2130 11361027600 : sub_int(i,j,2) = sub_int(i,j,2) + sub_int(i,j-1,2)
2131 : end do
2132 : end do
2133 :
2134 66903829200 : fluxes = 0
2135 5049345600 : do i=1,n
2136 16410373200 : do j=1,n
2137 11361027600 : if (1.lt.j) fluxes(i,j,1) = -sub_int(i, j-1,2)
2138 11361027600 : if (i.lt.n) fluxes(i,j,2) = sub_int(i, j, 1)
2139 11361027600 : if (j.lt.n) fluxes(i,j,3) = sub_int(i, j, 2)
2140 15148036800 : if (1.lt.i) fluxes(i,j,4) = -sub_int(i-1,j, 1)
2141 : end do
2142 : end do
2143 :
2144 1262336400 : end subroutine subcell_Laplace_fluxes
2145 :
2146 :
2147 : ! Given a field defined on the unit element, [-1,1]x[-1,1]
2148 : ! sample values, sampled_val, premultiplied by integration weights,
2149 : ! and a number, np, of Gauss-Lobatto-Legendre points. Divide
2150 : ! the square up into intervals by intervals sub-squares so that
2151 : ! there are now intervals**2 sub-cells. Integrate the
2152 : ! function defined by sampled_val over each of these
2153 : ! sub-cells and return the integrated values as an
2154 : ! intervals by intervals matrix.
2155 : !
2156 : ! Efficiency is obtained by computing and caching the appropriate
2157 : ! integration matrix the first time the function is called.
2158 11024812800 : subroutine subcell_integration(sampled_val, np, intervals, metdet,values,separate_physics_grid)
2159 :
2160 : implicit none
2161 :
2162 : integer , intent(in) :: np
2163 : integer , intent(in) :: intervals
2164 : real (kind=r8), intent(in) :: sampled_val(np,np)
2165 : real (kind=r8), intent(in) :: metdet(np,np)
2166 22049625600 : real (kind=r8) :: val(np,np)
2167 : real (kind=r8) :: values(intervals,intervals)
2168 :
2169 : logical, optional :: separate_physics_grid
2170 :
2171 : logical :: ltmp
2172 : integer i,j
2173 :
2174 11024812800 : if (present(separate_physics_grid)) then
2175 208656000 : ltmp = separate_physics_grid
2176 : else
2177 : ltmp = .false.
2178 : end if
2179 :
2180 :
2181 208656000 : if (ltmp) then
2182 : if (.not.ALLOCATED(integration_matrix_physgrid) .or. &
2183 0 : SIZE(integration_matrix_physgrid,1).ne.intervals .or. &
2184 : SIZE(integration_matrix_physgrid,2).ne.np) then
2185 0 : call endrun( 'FATAL ERROR: allocate_subcell_integration_matrix_physgrid not called')
2186 : end if
2187 0 : val = sampled_val * metdet
2188 : values = MATMUL(integration_matrix_physgrid, &
2189 0 : MATMUL(val,TRANSPOSE(integration_matrix_physgrid)))
2190 : else
2191 : if (.not.ALLOCATED(integration_matrix) .or. &
2192 11024812800 : SIZE(integration_matrix,1).ne.intervals .or. &
2193 : SIZE(integration_matrix,2).ne.np) then
2194 0 : call endrun( 'FATAL ERROR: allocate_subcell_integration_matrix not called')
2195 : end if
2196 : ! Multiply sampled values by spectral element weights
2197 >23152*10^7 : val = sampled_val * metdet
2198 :
2199 : ! Multiply the sampled values by the weighted jacobians.
2200 : ! Symmetry allows us to write this as J^t V J
2201 : ! where J is a vector.
2202 :
2203 : values = MATMUL(integration_matrix, &
2204 11024812800 : MATMUL(val,TRANSPOSE(integration_matrix)))
2205 :
2206 : end if
2207 11024812800 : end subroutine subcell_integration
2208 :
2209 :
2210 : ! Helper subroutine that will fill in a matrix needed to
2211 : ! integrate a function defined on the GLL points of a unit
2212 : ! square on sub-cells. So np is the number of integration
2213 : ! GLL points defined on the unit square (actually [-1,1]x[-1,1])
2214 : ! and intervals is the number to cut it up into, say a 3 by 3
2215 : ! set of uniform sub-cells. This function will fill the
2216 : ! subcell_integration matrix with the correct coefficients
2217 : ! to integrate over each subcell.
2218 1536 : subroutine allocate_subcell_integration_matrix_cslam(np, intervals)
2219 : !-----------------
2220 : !-----------------
2221 : use quadrature_mod, only : gausslobatto, quadrature_t
2222 :
2223 : implicit none
2224 :
2225 : integer , intent(in) :: np
2226 : integer , intent(in) :: intervals
2227 : real (kind=r8) :: values(intervals,intervals)
2228 :
2229 :
2230 : real(kind=r8), parameter :: zero = 0.0_r8, one=1.0_r8, two=2.0_r8
2231 :
2232 :
2233 3072 : real (kind=r8) :: sub_gll (intervals,np)
2234 :
2235 3072 : real (kind=r8) :: Lagrange_interp(intervals,np,np)
2236 : type (quadrature_t) :: gll
2237 :
2238 3072 : real (kind=r8) :: legrange_div(np)
2239 : real (kind=r8) :: a,b,x,y, x_j, x_i
2240 : real (kind=r8) :: r(1)
2241 : integer i,j,n,m
2242 :
2243 1536 : if (ALLOCATED(integration_matrix)) deallocate(integration_matrix)
2244 6144 : allocate(integration_matrix(intervals,np))
2245 1536 : if (ALLOCATED(boundary_interp_matrix)) deallocate(boundary_interp_matrix)
2246 6144 : allocate(boundary_interp_matrix(intervals,2,np))
2247 :
2248 1536 : gll = gausslobatto(np)
2249 :
2250 : ! The GLL (Gauss-Lobatto-Legendre) points are from [-1,1],
2251 : ! we have a bunch of sub-intervals defined by intervals that
2252 : ! go from [a,b] so we need to linearly map [-1,1] -> [a,b]
2253 : ! all the GLL points by y = (a/2)(1-x) + (b/2)(1+x)
2254 6144 : do i=1,intervals
2255 4608 : a = -one + (i-one)*two/intervals
2256 4608 : b = -one + i *two/intervals
2257 24576 : sub_gll(i,:) = (a+b)/two + gll%points(:)/intervals
2258 : end do
2259 :
2260 : ! Now to interpolate from the values at the input GLL
2261 : ! points to the sub-GLL points. Do this by Lagrange
2262 : ! interpolation. The jth Lagrange interpolating polynomial
2263 : ! for points x_i is
2264 : ! \prod_{i\ne j} (x-x_i)/(x_j-x_i)
2265 : ! These are then multiplied by the sampled values y_i
2266 : ! and summed.
2267 :
2268 : ! Save some time by pre-computing the denominitor. I think
2269 : ! this is OK since all the points are of order 1 so should
2270 : ! be well behaved.
2271 7680 : do n = 1,np
2272 6144 : x_j = gll%points(n)
2273 6144 : x = one
2274 30720 : do m = 1,np
2275 30720 : if (m.ne.n) then
2276 18432 : x_i = gll%points(m)
2277 18432 : x = x * (x_j-x_i)
2278 : endif
2279 : end do
2280 7680 : legrange_div(n)= x
2281 : end do
2282 6144 : do i=1,intervals
2283 24576 : do n=1,np
2284 18432 : x = sub_gll(i,n)
2285 96768 : do j = 1,np
2286 : y = one
2287 368640 : do m = 1,np
2288 368640 : if (m.ne.j) then
2289 221184 : x_i = gll%points(m)
2290 221184 : y = y * (x-x_i)
2291 : end if
2292 : end do
2293 92160 : Lagrange_interp(i,n,j) = y/legrange_div(j)
2294 : end do
2295 : end do
2296 : end do
2297 :
2298 : ! Integration is the GLL weights times Jacobians times
2299 : ! the interpolated values:
2300 : ! w^t I Y I^t w
2301 : ! where
2302 : ! w is GLL weights and Jacobians,
2303 : ! I is the Lagrange_interp matrix, and
2304 : ! Y is the coefficient matrix, sampled_val.
2305 : ! This can be written J Y J^t where
2306 : ! J = w^t I
2307 : ! J is integration_matrix
2308 6144 : do i=1,intervals
2309 24576 : integration_matrix(i,:) = MATMUL(gll%weights(:),Lagrange_interp(i,:,:))
2310 : end do
2311 :
2312 : ! There is still the Jacobian to consider. We are
2313 : ! integrating over [a,b] x [c,d] where
2314 : ! |b-a| = |d-c| = 2/Intervals
2315 : ! Multiply the weights appropriately given that
2316 : ! they are defined for a 2x2 square
2317 26112 : integration_matrix = integration_matrix/intervals
2318 :
2319 26112 : boundary_interp_matrix(:,1,:) = Lagrange_interp(:,1,:)
2320 26112 : boundary_interp_matrix(:,2,:) = Lagrange_interp(:,np,:)
2321 1536 : end subroutine allocate_subcell_integration_matrix_cslam
2322 :
2323 1536 : subroutine allocate_subcell_integration_matrix_physgrid(np, intervals)
2324 : !-----------------
2325 : !-----------------
2326 : use quadrature_mod, only : gausslobatto, quadrature_t
2327 :
2328 : implicit none
2329 :
2330 : integer , intent(in) :: np
2331 : integer , intent(in) :: intervals
2332 : real (kind=r8) :: values(intervals,intervals)
2333 :
2334 :
2335 : real(kind=r8), parameter :: zero = 0.0_r8, one=1.0_r8, two=2.0_r8
2336 :
2337 :
2338 3072 : real (kind=r8) :: sub_gll (intervals,np)
2339 :
2340 3072 : real (kind=r8) :: Lagrange_interp(intervals,np,np)
2341 : type (quadrature_t) :: gll
2342 :
2343 3072 : real (kind=r8) :: legrange_div(np)
2344 : real (kind=r8) :: a,b,x,y, x_j, x_i
2345 : real (kind=r8) :: r(1)
2346 : integer i,j,n,m
2347 :
2348 1536 : if (ALLOCATED(integration_matrix_physgrid)) deallocate(integration_matrix_physgrid)
2349 6144 : allocate(integration_matrix_physgrid(intervals,np))
2350 :
2351 1536 : gll = gausslobatto(np)
2352 :
2353 : ! The GLL (Gauss-Lobatto-Legendre) points are from [-1,1],
2354 : ! we have a bunch of sub-intervals defined by intervals that
2355 : ! go from [a,b] so we need to linearly map [-1,1] -> [a,b]
2356 : ! all the GLL points by y = (a/2)(1-x) + (b/2)(1+x)
2357 6144 : do i=1,intervals
2358 4608 : a = -one + (i-one)*two/intervals
2359 4608 : b = -one + i *two/intervals
2360 24576 : sub_gll(i,:) = (a+b)/two + gll%points(:)/intervals
2361 : end do
2362 :
2363 : ! Now to interpolate from the values at the input GLL
2364 : ! points to the sub-GLL points. Do this by Lagrange
2365 : ! interpolation. The jth Lagrange interpolating polynomial
2366 : ! for points x_i is
2367 : ! \prod_{i\ne j} (x-x_i)/(x_j-x_i)
2368 : ! These are then multiplied by the sampled values y_i
2369 : ! and summed.
2370 :
2371 : ! Save some time by pre-computing the denominitor. I think
2372 : ! this is OK since all the points are of order 1 so should
2373 : ! be well behaved.
2374 7680 : do n = 1,np
2375 6144 : x_j = gll%points(n)
2376 6144 : x = one
2377 30720 : do m = 1,np
2378 30720 : if (m.ne.n) then
2379 18432 : x_i = gll%points(m)
2380 18432 : x = x * (x_j-x_i)
2381 : endif
2382 : end do
2383 7680 : legrange_div(n)= x
2384 : end do
2385 6144 : do i=1,intervals
2386 24576 : do n=1,np
2387 18432 : x = sub_gll(i,n)
2388 96768 : do j = 1,np
2389 : y = one
2390 368640 : do m = 1,np
2391 368640 : if (m.ne.j) then
2392 221184 : x_i = gll%points(m)
2393 221184 : y = y * (x-x_i)
2394 : end if
2395 : end do
2396 92160 : Lagrange_interp(i,n,j) = y/legrange_div(j)
2397 : end do
2398 : end do
2399 : end do
2400 6144 : do i=1,intervals
2401 24576 : integration_matrix_physgrid(i,:) = MATMUL(gll%weights(:),Lagrange_interp(i,:,:))
2402 : end do
2403 26112 : integration_matrix_physgrid = integration_matrix_physgrid/intervals
2404 1536 : end subroutine allocate_subcell_integration_matrix_physgrid
2405 :
2406 :
2407 :
2408 0 : subroutine limiter_optim_iter_full(ptens,sphweights,minp,maxp,dpmass,kbeg,kend)
2409 : !
2410 : !The idea here is the following: We need to find a grid field which is closest
2411 : !to the initial field (in terms of weighted sum), but satisfies the min/max constraints.
2412 : !So, first we find values which do not satisfy constraints and bring these values
2413 : !to a closest constraint. This way we introduce some mass change (addmass),
2414 : !so, we redistribute addmass in the way that l2 error is smallest.
2415 : !This redistribution might violate constraints thus, we do a few iterations.
2416 : !
2417 : ! O. Guba ~2012 Documented in Guba, Taylor & St-Cyr, JCP 2014
2418 : ! I. Demeshko & M. Taylor 7/2015: Removed indirect addressing.
2419 : ! N. Lopez & M. Taylor 8/2015: Mass redistributon tweak which is better at
2420 : ! linear coorelation preservation
2421 : !
2422 : use dimensions_mod, only : np, np, nlev
2423 :
2424 : real (kind=r8), dimension(nlev), intent(inout) :: minp, maxp
2425 : real (kind=r8), dimension(np*np,nlev), intent(inout) :: ptens
2426 : real (kind=r8), dimension(np*np,nlev), intent(in), optional :: dpmass
2427 : real (kind=r8), dimension(np*np), intent(in) :: sphweights
2428 : integer, intent(in) :: kbeg, kend
2429 :
2430 : real (kind=r8), dimension(np,np) :: ptens_mass
2431 : integer k1, k, i, j, iter, weightsnum
2432 : real (kind=r8) :: addmass, weightssum, mass, sumc
2433 : real (kind=r8) :: x(np*np),c(np*np)
2434 : integer :: maxiter = np*np-1
2435 : real (kind=r8) :: tol_limiter = 5.0e-14_r8
2436 :
2437 0 : do k = kbeg, kend
2438 :
2439 0 : do k1=1,np*np
2440 0 : c(k1)=sphweights(k1)*dpmass(k1,k)
2441 0 : x(k1)=ptens(k1,k)/dpmass(k1,k)
2442 : enddo
2443 :
2444 0 : sumc=sum(c)
2445 0 : if (sumc <= 0 ) CYCLE ! this should never happen, but if it does, dont limit
2446 0 : mass=sum(c*x)
2447 :
2448 :
2449 :
2450 : ! relax constraints to ensure limiter has a solution:
2451 : ! This is only needed if runnign with the SSP CFL>1 or
2452 : ! due to roundoff errors
2453 0 : if( mass < minp(k)*sumc ) then
2454 0 : minp(k) = mass / sumc
2455 : endif
2456 0 : if( mass > maxp(k)*sumc ) then
2457 0 : maxp(k) = mass / sumc
2458 : endif
2459 :
2460 :
2461 :
2462 0 : do iter=1,maxiter
2463 :
2464 : addmass=0.0d0
2465 :
2466 0 : do k1=1,np*np
2467 0 : if((x(k1)>maxp(k))) then
2468 0 : addmass=addmass+(x(k1)-maxp(k))*c(k1)
2469 0 : x(k1)=maxp(k)
2470 : endif
2471 0 : if((x(k1)<minp(k))) then
2472 0 : addmass=addmass-(minp(k)-x(k1))*c(k1)
2473 0 : x(k1)=minp(k)
2474 : endif
2475 : enddo !k1
2476 :
2477 0 : if(abs(addmass)<=tol_limiter*abs(mass)) exit
2478 :
2479 0 : weightssum=0.0d0
2480 0 : if(addmass>0)then
2481 0 : do k1=1,np*np
2482 0 : if(x(k1)<maxp(k))then
2483 0 : weightssum=weightssum+c(k1)
2484 : endif
2485 : enddo !k1
2486 0 : do k1=1,np*np
2487 0 : if(x(k1)<maxp(k))then
2488 0 : x(k1)=x(k1)+addmass/weightssum
2489 : endif
2490 : enddo
2491 : else
2492 0 : do k1=1,np*np
2493 0 : if(x(k1)>minp(k))then
2494 0 : weightssum=weightssum+c(k1)
2495 : endif
2496 : enddo
2497 0 : do k1=1,np*np
2498 0 : if(x(k1)>minp(k))then
2499 0 : x(k1)=x(k1)+addmass/weightssum
2500 : endif
2501 : enddo
2502 : endif
2503 :
2504 :
2505 : enddo!end of iteration
2506 :
2507 0 : do k1=1,np*np
2508 0 : ptens(k1,k)=x(k1)
2509 : enddo
2510 :
2511 : enddo
2512 :
2513 0 : do k = kbeg, kend
2514 0 : do k1=1,np*np
2515 0 : ptens(k1,k)=ptens(k1,k)*dpmass(k1,k)
2516 : enddo
2517 : enddo
2518 :
2519 0 : end subroutine limiter_optim_iter_full
2520 :
2521 :
2522 :
2523 :
2524 :
2525 0 : end module derivative_mod
|
