scale_atm_dyn_dgm_globalsw.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scaleFElib.h"
module scale_atm_dyn_dgm_globalsw
  !-----------------------------------------------------------------------------
  !
  !++ Used modules
  !
  use scale_precision
  use scale_io
  use scale_prc
  use scale_prof
 
  use scale_const, only: &
    grav => const_grav
  
  use scale_sparsemat
  use scale_element_base, only: &
    elementbase2d
  use scale_element_quadrilateral, only: quadrilateralelement
  use scale_localmesh_2d, only: localmesh2d  
  use scale_mesh_base2d, only: meshbase2d
  use scale_localmeshfield_base, only: localmeshfield2d
  use scale_meshfield_base, only: meshfield2d
 
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedures
  !
  public :: atm_dyn_dgm_globalsw_init
  public :: atm_dyn_dgm_globalsw_final
  public :: atm_dyn_dgm_globalsw_cal_tend
 
  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  
  !-----------------------------------------------------------------------------
  !
  !++ Private procedures & variables
  !
  !-------------------
 
  integer, private, parameter :: VARS_H_ID      = 1
  integer, private, parameter :: VARS_U1_ID     = 2
  integer, private, parameter :: VARS_u2_ID     = 3
  integer, private, parameter :: PROG_VARS_NUM  = 3
  
  private :: cal_del_flux_dyn
 
contains
  subroutine atm_dyn_dgm_globalsw_init( mesh )
    implicit none
    class(meshbase2d), intent(in) :: mesh
    !--------------------------------------------
    return
  end subroutine atm_dyn_dgm_globalsw_init
 
 
  subroutine atm_dyn_dgm_globalsw_final()
    implicit none
    !--------------------------------------------
    return
  end subroutine atm_dyn_dgm_globalsw_final  
 
  !-------------------------------
 
!OCL SERIAL
  subroutine atm_dyn_dgm_globalsw_cal_tend( &
    h_dt, U_dt, V_dt,                       & ! (out)
    h_, u_, v_, hs_, u1_, u2_, coriolis,    & ! (in)
    dx, dy, sx, sy, lift, lmesh, elem  )
 
    implicit none
 
    class(localmesh2d), intent(in) :: lmesh
    class(elementbase2d), intent(in) :: elem
    type(sparsemat), intent(in) :: dx, dy, sx, sy, lift
    real(rp), intent(out) :: h_dt(elem%np,lmesh%nea)
    real(rp), intent(out) :: u_dt(elem%np,lmesh%nea)
    real(rp), intent(out) :: v_dt(elem%np,lmesh%nea)
    real(rp), intent(in)  :: h_(elem%np,lmesh%nea)
    real(rp), intent(in)  :: u_(elem%np,lmesh%nea)
    real(rp), intent(in)  :: v_(elem%np,lmesh%nea)
    real(rp), intent(in)  :: hs_(elem%np,lmesh%nea)
    real(rp), intent(in)  :: u1_(elem%np,lmesh%nea)
    real(rp), intent(in)  :: u2_(elem%np,lmesh%nea)
    real(rp), intent(in)  :: coriolis(elem%np,lmesh%nea)
 
    real(rp) :: fx(elem%np), fy(elem%np), liftdelflx(elem%np)
    real(rp) :: del_flux(elem%nfptot,lmesh%ne,prog_vars_num)
    real(rp) :: del_flux_aux(elem%nfptot,lmesh%ne,1)
    real(rp) :: u1_dt(elem%np), u2_dt(elem%np)
    real(rp) :: vor(elem%np), e(elem%np)
    integer :: ke       
    !------------------------------------------------------------------------
 
    call prof_rapstart('cal_dyn_tend_bndflux', 3)
    call cal_del_flux_dyn( del_flux, del_flux_aux,         & ! (out)
      h_, u_, v_, hs_, u1_, u2_,                           & ! (in)
      lmesh%Gsqrt, lmesh%GIJ(:,:,1,1), lmesh%GIJ(:,:,2,2), & ! (in)
      lmesh%normal_fn(:,:,1), lmesh%normal_fn(:,:,2),      & ! (in)
      lmesh%vmapM, lmesh%vmapP,                            & ! (in)
      lmesh, elem )                                          ! (in)
    call prof_rapend('cal_dyn_tend_bndflux', 3)
 
    !-----
    call prof_rapstart('cal_dyn_tend_interior', 3)
 
    !$omp parallel do private( &
    !$omp E, VOR,                          &
    !$omp u1_dt, u2_dt, Fx, Fy, LiftDelFlx )
    do ke = lmesh%NeS, lmesh%NeE
      !--
      call sparsemat_matmul(dx, u2_(:,ke), fx)
      call sparsemat_matmul(dy, u1_(:,ke), fy)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux_aux(:,ke,1), liftdelflx)
 
      vor(:) = ( &
          lmesh%Escale(:,ke,1,1) * fx(:) &
        - lmesh%Escale(:,ke,2,2) * fy(:) &
        + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
      
      e(:) = grav * ( hs_(:,ke) + h_(:,ke) )                        &
           + 0.5_rp * ( u1_(:,ke) * u_(:,ke) + u2_(:,ke) * v_(:,ke) )
      
      !-- h
      call sparsemat_matmul(dx, lmesh%Gsqrt(:,ke) * h_(:,ke) * u_(:,ke), fx)
      call sparsemat_matmul(dy, lmesh%Gsqrt(:,ke) * h_(:,ke) * v_(:,ke), fy)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,vars_h_id), liftdelflx)
 
      h_dt(:,ke) = - ( &
              lmesh%Escale(:,ke,1,1) * fx(:) &
            + lmesh%Escale(:,ke,2,2) * fy(:) &
            + liftdelflx(:) ) / lmesh%Gsqrt(:,ke)
      
      !-- u1
      call sparsemat_matmul(dx, e(:), fx)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,vars_u1_id), liftdelflx)
 
      u1_dt(:) = - (  &
            lmesh%Escale(:,ke,1,1) * fx(:)                                &
          + liftdelflx(:)                                                 &
          - lmesh%Gsqrt(:,ke) * ( coriolis(:,ke) + vor(:) ) * v_(:,ke)    )
      
      !-- u2
      call sparsemat_matmul(dy, e(:), fy)
      call sparsemat_matmul(lift, lmesh%Fscale(:,ke) * del_flux(:,ke,vars_u2_id), liftdelflx)
 
      u2_dt(:) = - (  &
            lmesh%Escale(:,ke,2,2) * fy(:)                                &
          + liftdelflx(:)                                                 &
          + lmesh%Gsqrt(:,ke) * ( coriolis(:,ke) + vor(:) ) * u_(:,ke)    )
    
      !--
      u_dt(:,ke) = lmesh%GIJ(:,ke,1,1) * u1_dt(:) + lmesh%GIJ(:,ke,1,2) * u2_dt(:) 
      v_dt(:,ke) = lmesh%GIJ(:,ke,2,1) * u1_dt(:) + lmesh%GIJ(:,ke,2,2) * u2_dt(:) 
    end do
 
    call prof_rapend('cal_dyn_tend_interior', 3)
 
    return
  end subroutine atm_dyn_dgm_globalsw_cal_tend
 
  !------
 
!OCL SERIAL
  subroutine cal_del_flux_dyn( del_flux, del_flux_aux,     &
    h_, U_, V_, hs_, u1_, u2_,                             &
    Gsqrt_, G11_, G22_, nx, ny, vmapM, vmapP, lmesh, elem  )
 
    implicit none
 
    class(localmesh2d), intent(in) :: lmesh
    class(elementbase2d), intent(in) :: elem  
    real(rp), intent(out) ::  del_flux(elem%nfptot,lmesh%ne,prog_vars_num)
    real(rp), intent(out) ::  del_flux_aux(elem%nfptot,lmesh%ne,1)
    real(rp), intent(in) ::  h_(elem%np*lmesh%nea)
    real(rp), intent(in) ::  u_(elem%np*lmesh%nea)  
    real(rp), intent(in) ::  v_(elem%np*lmesh%nea)  
    real(rp), intent(in) ::  hs_(elem%np*lmesh%nea)  
    real(rp), intent(in) ::  u1_(elem%np*lmesh%nea)  
    real(rp), intent(in) ::  u2_(elem%np*lmesh%nea)
    real(rp), intent(in) ::  gsqrt_(elem%np*lmesh%nea) 
    real(rp), intent(in) ::  g11_(elem%np*lmesh%ne) 
    real(rp), intent(in) ::  g22_(elem%np*lmesh%ne) 
    real(rp), intent(in) :: nx(elem%nfptot,lmesh%ne)
    real(rp), intent(in) :: ny(elem%nfptot,lmesh%ne)
    integer, intent(in) :: vmapm(elem%nfptot,lmesh%ne)
    integer, intent(in) :: vmapp(elem%nfptot,lmesh%ne)
    
    integer :: ke, ip(elem%nfptot), im(elem%nfptot)
    real(rp) :: velp(elem%nfptot), velm(elem%nfptot), alpha(elem%nfptot)
    real(rp) :: h_p(elem%nfptot), h_m(elem%nfptot)
    real(rp) :: e_p(elem%nfptot), e_m(elem%nfptot)
    real(rp) :: gii(elem%nfptot)
    !------------------------------------------------------------------------
 
 
    !$omp parallel do private( iM, iP, &
    !$omp alpha, VelM, VelP,           &
    !$omp h_P, h_M, E_P, E_M, Gii )
    do ke=lmesh%NeS, lmesh%NeE
      im(:) = vmapm(:,ke); ip(:) = vmapp(:,ke)
 
      h_m(:) = h_(im)
      h_p(:) = h_(ip)
      velm(:) = u_(im) * nx(:,ke) + v_(im) * ny(:,ke)
      velp(:) = u_(ip) * nx(:,ke) + v_(ip) * ny(:,ke)
      gii(:) = g11_(im) * abs(nx(:,ke)) + g22_(im) * abs(ny(:,ke))
 
      e_m(:) = 0.5_rp * ( u_(im) * u1_(im) + v_(im) * u2_(im) ) &
             + grav * ( h_m(:) + hs_(im) )
      e_p(:) = 0.5_rp * ( u_(ip) * u1_(ip) + v_(ip) * u2_(ip) ) &
             + grav * ( h_p(:) + hs_(ip) )
      
      alpha(:) = max( sqrt( gii(:) * grav * (h_m(:) + hs_(im)) ) + abs(velm(:)), &
                      sqrt( gii(:) * grav * (h_p(:) + hs_(ip)) ) + abs(velp(:))  )
      
      del_flux(:,ke,vars_h_id) = 0.5_rp * gsqrt_(im) * (     &
                      (h_p(:) * velp(:) - h_m(:) * velm(:) ) &
                    - alpha(:) * ( h_p(:) - h_m(:) )         )
 
      del_flux(:,ke,vars_u1_id) = 0.5_rp * (               &
                    ( e_p(:)  - e_m(:) ) * nx(:,ke)        &
                    - alpha(:) * ( u1_(ip) - u1_(im) )     )
 
      del_flux(:,ke,vars_u2_id) = 0.5_rp * (               &
                    ( e_p(:)  - e_m(:) ) * ny(:,ke)        &
                    - alpha(:) * ( u2_(ip) - u2_(im) )     )
      
      del_flux_aux(:,ke,1) = 0.5_rp * ( &
                     ( u2_(ip) - u2_(im) ) * nx(:,ke) &
                   - ( u1_(ip) - u1_(im) ) * ny(:,ke) )
    end do
 
    return
  end subroutine cal_del_flux_dyn
 
end module scale_atm_dyn_dgm_globalsw