Line data    Source code

       1             : module mass_matrix_mod
       2             :   use shr_kind_mod,   only: r8=>shr_kind_r8
       3             :   use dimensions_mod, only: np, nelemd
       4             :   use quadrature_mod, only: quadrature_t, gauss ,gausslobatto
       5             :   use element_mod,    only: element_t
       6             :   use parallel_mod,   only: parallel_t
       7             :   use edge_mod,       only: edgevpack, edgevunpack,                           &
       8             :        freeedgebuffer,initedgebuffer
       9             :   use edgetype_mod,   only: edgebuffer_t
      10             :   use bndry_mod,      only: bndry_exchange
      11             : 
      12             : implicit none
      13             : private
      14             : 
      15             :   public :: mass_matrix
      16             : 
      17             : contains
      18             : 
      19             : ! ===========================================
      20             : ! mass_matrix:
      21             : !
      22             : ! Compute the mass matrix for each element...
      23             : ! ===========================================
      24             : 
      25        3072 :   subroutine mass_matrix(par,elem)
      26             : 
      27             :     type (parallel_t),intent(in) :: par
      28             :     type (element_t) :: elem(:)
      29             : 
      30        3072 :     type (EdgeBuffer_t)    :: edge
      31             : 
      32             :     real(kind=r8)  da                     ! area element
      33             : 
      34             :     type (quadrature_t) :: gp
      35             : 
      36             :     integer ii
      37             :     integer i,j
      38             :     integer kptr
      39             :     integer iptr
      40             : 
      41             :     ! ===================
      42             :     ! begin code
      43             :     ! ===================
      44             : 
      45        3072 :     call initEdgeBuffer(par,edge,elem,1,nthreads=1)
      46             : 
      47             :     ! =================================================
      48             :     ! mass matrix on the velocity grid
      49             :     ! =================================================    
      50             : 
      51        3072 :     gp=gausslobatto(np)
      52             :  
      53       24672 :     do ii=1,nelemd
      54      108000 :        do j=1,np
      55      453600 :           do i=1,np
      56             :               ! MNL: metric term for map to reference element is now in metdet!
      57      345600 :              elem(ii)%mp(i,j)=gp%weights(i)*gp%weights(j)
      58      432000 :              elem(ii)%rmp(i,j)=elem(ii)%mp(i,j)
      59             :           end do
      60             :        end do
      61             : 
      62       21600 :        kptr=0
      63       24672 :        call edgeVpack(edge,elem(ii)%rmp,1,kptr,ii)
      64             : 
      65             :     end do
      66             : 
      67             :     ! ==============================
      68             :     ! Insert boundary exchange here
      69             :     ! ==============================
      70             : 
      71        3072 :     call bndry_exchange(par,edge,location='mass_matrix #1')
      72             : 
      73       24672 :     do ii=1,nelemd
      74             : 
      75       21600 :        kptr=0
      76       21600 :        call edgeVunpack(edge,elem(ii)%rmp,1,kptr,ii)
      77             : 
      78      111072 :        do j=1,np
      79      453600 :           do i=1,np
      80      432000 :              elem(ii)%rmp(i,j)=1.0_r8/elem(ii)%rmp(i,j)
      81             :           end do
      82             :        end do
      83             : 
      84             :     end do
      85             : !$OMP BARRIER
      86             : 
      87        3072 :     deallocate(gp%points)
      88        3072 :     deallocate(gp%weights)
      89             : 
      90             :     ! =============================================
      91             :     ! compute spherical element mass matrix
      92             :     ! =============================================
      93       24672 :     do ii=1,nelemd
      94      108000 :        do j=1,np
      95      453600 :           do i=1,np
      96      345600 :              elem(ii)%spheremp(i,j)=elem(ii)%mp(i,j)*elem(ii)%metdet(i,j)
      97      432000 :              elem(ii)%rspheremp(i,j)=elem(ii)%spheremp(i,j)
      98             :           end do
      99             :        end do
     100       21600 :        kptr=0
     101       24672 :        call edgeVpack(edge,elem(ii)%rspheremp,1,kptr,ii)
     102             :     end do
     103        3072 :     call bndry_exchange(par,edge,location='mass_matrix #2')
     104       24672 :     do ii=1,nelemd
     105       21600 :        kptr=0
     106       21600 :        call edgeVunpack(edge,elem(ii)%rspheremp,1,kptr,ii)
     107      111072 :        do j=1,np
     108      453600 :           do i=1,np
     109      432000 :              elem(ii)%rspheremp(i,j)=1.0_r8/elem(ii)%rspheremp(i,j)
     110             :           end do
     111             :        end do
     112             :     end do
     113             : !$OMP BARRIER
     114             : 
     115        3072 :     call FreeEdgeBuffer(edge)
     116             :        
     117        3072 :   end subroutine mass_matrix
     118             : 
     119             : end module mass_matrix_mod
     120             :