shelf_triangular_festuff Module Reference

Data Types
type	ice_shelf_cs

Functions/Subroutines
subroutine	matrix_diagonal_triangle (CS, u_diagonal, v_diagonal)

Function/Subroutine Documentation

matrix_diagonal_triangle()

subroutine shelf_triangular_festuff::matrix_diagonal_triangle	(	type(ice_shelf_cs), pointer	CS,
		real, dimension (:,:), intent(inout)	u_diagonal,
		real, dimension (:,:), intent(inout)	v_diagonal
	)

Definition at line 165 of file shelf_triangular_FEstuff.F90.

  type(ice_shelf_cs),    pointer       :: cs
  real, dimension (:,:), intent(inout) :: u_diagonal, v_diagonal

! returns the diagonal entries of the matrix for a Jacobi preconditioning

  real, pointer, dimension (:,:)       :: umask, vmask, &
                          nu_lower, nu_upper, beta_lower, beta_upper, hmask
  type(ocean_grid_type), pointer :: g
  integer ::  i, j, is, js, cnt, isc, jsc, iec, jec
  real :: a, n, ux, uy, vx, vy, eps_min, domain_width, dxh, dyh, dxdyh

  g => cs%grid

!   if (G%symmetric) then
!     isym=1
!   else
!     isym=0
!   endif


  isc = g%isc ; jsc = g%jsc ; iec = g%iec ; jec = g%jec

  umask => cs%umask ; vmask => cs%vmask ; hmask => cs%hmask
  nu_lower => cs%ice_visc_lower_tri ; nu_upper => cs%ice_visc_upper_tri
  beta_lower => cs%taub_beta_eff_lower_tri ; beta_upper => cs%taub_beta_eff_upper_tri

  do i=isc-1,iec+1  ; do j=jsc-1,jec+1 ; if (hmask(i,j) .eq. 1) then
    dxh = g%dxT(i,j)
    dyh = g%dyT(i,j)
    dxdyh = g%areaT(i,j)

    if (umask(i,j-1) .eq. 1) then ! this (bot right) is a degree of freedom node

      ux = 1./dxh ; uy = 0./dyh
      vx = 0. ; vy = 0.

      u_diagonal(i,j-1) = u_diagonal(i,j-1) + &
          .5 * dxdyh * nu_lower(i,j) * ((4*ux+2*vy) * (1./dxh) + (uy+vy) * (0./dyh))

      u_diagonal(i,j-1) = u_diagonal(i,j-1) + &
          beta_lower(i,j) * dxdyh * 1./24

      ux = 0. ; uy = 0.
      vx = 1./dxh ; vy = 0./dyh

      v_diagonal(i,j-1) = v_diagonal(i,j-1) + &
          .5 * dxdyh * nu_lower(i,j) * ((uy+vx) * (1./dxh) + (4*vy+2*ux) * (0./dyh))

      v_diagonal(i,j-1) = v_diagonal(i,j-1) + &
          beta_lower(i,j) * dxdyh * 1./24

      ux = 0./dxh ; uy = -1./dyh
      vx = 0. ; vy = 0.

      u_diagonal(i,j-1) = u_diagonal(i,j-1) + &
          .5 * dxdyh * nu_upper(i,j) * ((4*ux+2*vy) * (0./dxh) + (uy+vy) * (-1./dyh))

      u_diagonal(i,j-1) = u_diagonal(i,j-1) + &
          beta_upper(i,j) * dxdyh * 1./24

      vx = 0./dxh ; vy = -1./dyh
      ux = 0. ; uy = 0.

      v_diagonal(i,j-1) = v_diagonal(i,j-1) + &
          .5 * dxdyh * nu_upper(i,j) * ((uy+vx) * (0./dxh) + (4*vy+2*ux) * (-1./dyh))

      v_diagonal(i,j-1) = v_diagonal(i,j-1) + &
          beta_upper(i,j) * dxdyh * 1./24

    endif

    if (umask(i-1,j) .eq. 1) then ! this (top left) is a degree of freedom node

      ux = 0./dxh ; uy = 1./dyh
      vx = 0. ; vy = 0.

      u_diagonal(i-1,j) = u_diagonal(i-1,j) + &
          .5 * dxdyh * nu_lower(i,j) * ((4*ux+2*vy) * (0./dxh) + (uy+vy) * (1./dyh))

      u_diagonal(i,j-1) = u_diagonal(i,j-1) + &
          beta_lower(i,j) * dxdyh * 1./24

      ux = 0. ; uy = 0.
      vx = 0./dxh ; vy = 1./dyh

      v_diagonal(i-1,j) = v_diagonal(i-1,j) + &
          .5 * dxdyh * nu_lower(i,j) * ((uy+vx) * (0./dxh) + (4*vy+2*ux) * (1./dyh))

      v_diagonal(i,j-1) = v_diagonal(i,j-1) + &
          beta_lower(i,j) * dxdyh * 1./24

      ux = -1./dxh ; uy = 0./dyh
      vx = 0. ; vy = 0.

      u_diagonal(i-1,j) = u_diagonal(i-1,j) + &
          .5 * dxdyh * nu_upper(i,j) * ((4*ux+2*vy) * (-1./dxh) + (uy+vy) * (0./dyh))

      u_diagonal(i,j-1) = u_diagonal(i,j-1) + &
          beta_upper(i,j) * dxdyh * 1./24

      vx = -1./dxh ; vy = 0./dyh
      ux = 0. ; uy = 0.

      v_diagonal(i-1,j) = v_diagonal(i-1,j) + &
          .5 * dxdyh * nu_upper(i,j) * ((uy+vx) * (-1./dxh) + (4*vy+2*ux) * (0./dyh))

      v_diagonal(i,j-1) = v_diagonal(i,j-1) + &
          beta_upper(i,j) * dxdyh * 1./24

    endif

    if (umask(i-1,j-1) .eq. 1) then ! this (bot left) is a degree of freedom node

      ux = -1./dxh ; uy = -1./dyh
      vx = 0. ; vy = 0.

      u_diagonal(i-1,j-1) = u_diagonal(i-1,j-1) + &
          .5 * dxdyh * nu_upper(i,j) * ((4*ux+2*vy) * (-1./dxh) + (uy+vy) * (-1./dyh))

      u_diagonal(i-1,j-1) = u_diagonal(i-1,j-1) + &
          beta_lower(i,j) * dxdyh * 1./24

      vx = -1./dxh ; vy = -1./dyh
      ux = 0. ; uy = 0.

      v_diagonal(i-1,j-1) = v_diagonal(i-1,j-1) + &
          .5 * dxdyh * nu_upper(i,j) * ((uy+vx) * (-1./dxh) + (4*vy+2*ux) * (-1./dyh))

      v_diagonal(i-1,j-1) = v_diagonal(i-1,j-1) + &
          beta_lower(i,j) * dxdyh * 1./24
    endif

    if (umask(i,j) .eq. 1) then ! this (top right) is a degree of freedom node

      ux = 1./ dxh ; uy = 1./dyh
      vx = 0. ; vy = 0.

      u_diagonal(i,j) = u_diagonal(i,j) + &
          .5 * dxdyh * nu_upper(i,j) * ((4*ux+2*vy) * (1./dxh) + (uy+vy) * (1./dyh))

      u_diagonal(i,j) = u_diagonal(i,j) + &
          beta_upper(i,j) * dxdyh * 1./24

      vx = 1./ dxh ; vy = 1./dyh
      ux = 0. ; uy = 0.

      v_diagonal(i,j) = v_diagonal(i,j) + &
          .5 * dxdyh * nu_upper(i,j) * ((uy+vx) * (1./dxh) + (4*vy+2*ux) * (1./dyh))

      v_diagonal(i,j) = v_diagonal(i,j) + &
          beta_upper(i,j) * dxdyh * 1./24

    endif
  endif ; enddo ; enddo