angles.F90

Go to the documentation of this file.

module angles

  use gengrid_kinds, only : dbl_kind, int_kind
  use grdvars,       only : debug

  implicit none

contains

  subroutine find_angq(iind,jind,xangCt,anglet,angle)

    integer       , intent(in)  :: iind(2),jind(2)
    real(dbl_kind), intent(in)  :: xangCt(:)
    real(dbl_kind), intent(in)  :: anglet(:,:)
    real(dbl_kind), intent(out) :: angle(:,:)

    ! local variables
    integer :: i,j

    real(dbl_kind) :: angle_0, angle_w, angle_s, angle_sw
    real(dbl_kind) :: p25 = 0.25

    !---------------------------------------------------------------------
    ! find the angle on corners using the same relationship CICE uses
    ! internally to calculate angles on Ct using angles on Bu
    !
    !           w-----------------0 Ct(i+1,j+1)
    !           |                 |
    !        ----------Bu(i,j)---------- Bu lies on seam at j=nj
    !           |                 |
    !   Ct(i,j) sw----------------s
    !
    !---------------------------------------------------------------------

    angle = 0.0
    do j = jind(1)+1,jind(2)
       do i = iind(1),iind(2)-1
          if (j .lt. jind(2)) then
             angle_0  = anglet(i+1,j+1)
             angle_w  = anglet(i,  j+1)
             angle_s  = anglet(i+1,j  )
             angle_sw = anglet(i  ,j  )
          else
             angle_0  = xangct(i+1  )
             angle_w  = xangct(i    )
             angle_s  = anglet(i+1,j)
             angle_sw = anglet(i,  j)
          end if
          angle(i,j) = atan2(p25*(sin(angle_0) + sin(angle_w) + sin(angle_s) + sin(angle_sw)), &
                             p25*(cos(angle_0) + cos(angle_w) + cos(angle_s) + cos(angle_sw)))

          if (abs(angle(i,j)) .le. 1.0e-10)angle(i,j) = 0.0
       enddo
    enddo

  end subroutine find_angq

  subroutine find_angchk(iind,jind,angle,angchk)

    integer       , intent(in)  :: iind(2),jind(2)
    real(dbl_kind), intent(in)  :: angle(:,:)
    real(dbl_kind), intent(out) :: angchk(:,:)

    ! local variables
    integer :: i,j
    real(dbl_kind) :: angle_0, angle_w, angle_s, angle_sw
    real(dbl_kind) :: p25 = 0.25

    !---------------------------------------------------------------------
    ! check: calculate anglet from angle on corners as CICE does internally.
    ! since angle changes sign between CICE and MOM6, (-1)*angchk ~ anglet
    !
    !               w-----------------0 Bu(i,j)
    !               |                 |
    !               |     Ct(i,j)     |
    !               |                 |
    !   Bu(i-1,j-1) sw----------------s
    !
    !---------------------------------------------------------------------

    angchk = 0.0
    do j = jind(1)+1,jind(2)
       do i = iind(1)+1,iind(2)
          angle_0  = angle(i  ,j  )
          angle_w  = angle(i-1,j  )
          angle_s  = angle(i,  j-1)
          angle_sw = angle(i-1,j-1)
          angchk(i,j) = atan2(p25*(sin(angle_0) + sin(angle_w) + sin(angle_s) + sin(angle_sw)), &
                              p25*(cos(angle_0) + cos(angle_w) + cos(angle_s) + cos(angle_sw)))
       enddo
    enddo

  end subroutine find_angchk


  subroutine find_ang(iind,jind,lonBu,latBu,lonCt,anglet)

    integer       , intent(in)  :: iind(2),jind(2)
    real(dbl_kind), intent(in)  :: lonBu(:,:)
    real(dbl_kind), intent(in)  :: latBu(:,:)
    real(dbl_kind), intent(in)  :: lonCt(:,:)
    real(dbl_kind), intent(out) :: anglet(:,:)

    ! local variables
    integer :: i,j,m,n
    integer :: ii,jj

    ! from geolonB fix in MOM6
    real(dbl_kind) :: len_lon ! The periodic range of longitudes, usually 360 degrees.
    real(dbl_kind) :: pi_720deg ! One quarter the conversion factor from degrees to radians.
    real(dbl_kind) :: lonB(2,2)  ! The longitude of a point, shifted to have about the same value.
    real(dbl_kind) :: lon_scale = 0.0

    !---------------------------------------------------------------------
    ! rotation angle for "use_bugs" = false case from MOM6
    ! src/initialization/MOM_shared_initialization.F90 but allow for not
    ! having halo values
    ! note this does not reproduce sin_rot,cos_rot found in MOM6 output
    ! differences are ~O 10-6
    !---------------------------------------------------------------------

    anglet = 0.0
    pi_720deg = atan(1.0) / 180.0
    len_lon = 360.0
    do j=jind(1),jind(2); do i = iind(1),iind(2)
       do n=1,2 ; do m=1,2
          jj = j+n-2; ii = i+m-2
          if(jj .eq. 0)jj = 1
          if(ii .eq. 0)ii = iind(2)
          lonb(m,n) = modulo_around_point(lonbu(ii,jj), lonct(i,j), len_lon)
          !  lonB(m,n) = modulo_around_point(LonBu(I+m-2,J+n-2), LonCt(i,j), len_lon)
       enddo; enddo
       jj = j-1; ii = i-1
       if(jj .eq. 0)jj = 1
       if(ii .eq. 0)ii = iind(2)
       lon_scale = cos(pi_720deg*((latbu(ii,jj) + latbu(i,j)) + &
            (latbu(i,jj) + latbu(ii,j)) ) )
       anglet(i,j) = atan2(lon_scale*((lonb(1,2) - lonb(2,1)) + (lonb(2,2) - lonb(1,1))), &
            (latbu(ii,j) - latbu(i,jj)) + &
            (latbu(i,j) - latbu(ii,jj)) )

       !lon_scale = cos(pi_720deg*((LatBu(I-1,J-1) + LatBu(I,J)) + &
       !                           (LatBu(I,J-1) + LatBu(I-1,J)) ) )
       !anglet(i,j) = atan2(lon_scale*((lonB(1,2) - lonB(2,1)) + (lonB(2,2) - lonB(1,1))), &
       !               (LatBu(I-1,J) - LatBu(I,J-1)) + &
       !               (LatBu(I,J) - LatBu(I-1,J-1)) )
    enddo; enddo

  end subroutine find_ang
  ! -----------------------------------------------------------------------------
  function modulo_around_point(x, xc, Lx) result(x_mod)
    use gengrid_kinds, only : dbl_kind

    real(dbl_kind), intent(in) :: x
    real(dbl_kind), intent(in) :: xc
    real(dbl_kind), intent(in) :: lx
    real(dbl_kind) :: x_mod

    if (lx > 0.0) then
       x_mod = modulo(x - (xc - 0.5*lx), lx) + (xc - 0.5*lx)
    else
       x_mod = x
    endif
  end function modulo_around_point
end module angles