angles.F90

Go to the documentation of this file.

module angles

  use gengrid_kinds, only : dbl_kind, int_kind
  use grdvars,       only : ni,nj,nx,ny
  use grdvars,       only : x,y,xsgp1,ysgp1,sg_maxlat
  use grdvars,       only : latbu,lonbu,lonct
  use grdvars,       only : angq,anglet
  use grdvars,       only : debug

  implicit none

contains
  subroutine find_angq

    ! local variables
    integer :: i,j,i1,i2,m,n

    ! pole locations on SG
    integer(int_kind) :: ipolesg(2)

    ! 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

    !---------------------------------------------------------------------
    ! to find angleq on seam, replicate supergrid values across seam
    !---------------------------------------------------------------------

    angq = 0.0
    xsgp1 = 0.0; ysgp1 = 0.0
    !pole on supergrid
    ipolesg = -1
    j = ny
    do i = 1,nx/2
       if(y(i,j) .eq. sg_maxlat)ipolesg(1) = i
    enddo
    do i = nx/2+1,nx
       if(y(i,j) .eq. sg_maxlat)ipolesg(2) = i
    enddo
    if(debug)print *,'poles found at ',ipolesg

    xsgp1(:,0:ny) = x(:,0:ny)
    ysgp1(:,0:ny) = y(:,0:ny)

    !check
    do i = ipolesg(1)-5,ipolesg(1)+5
       i2 = ipolesg(2)+(ipolesg(1)-i)+1
       if(debug)print *,i,i2
    enddo
    print *
    do i = ipolesg(2)-5,ipolesg(2)+5
       i2 = ipolesg(2)+(ipolesg(1)-i)+1
       if(debug)print *,i,i2
    enddo

    !replicate supergrid across pole
    do i = 1,nx
       i2 = ipolesg(2)+(ipolesg(1)-i)
       xsgp1(i,ny+1) = xsgp1(i2,ny)
       ysgp1(i,ny+1) = ysgp1(i2,ny)
    enddo

    !check
    j = ny+1
    i1 = ipolesg(1); i2 = ipolesg(2)-(ipolesg(1)-i1)
    print *,'replicate X across seam on SG'
    print *,xsgp1(i1-2,j),xsgp1(i2+2,j)
    print *,xsgp1(i1-1,j),xsgp1(i2+1,j)
    print *,xsgp1(i1,  j),xsgp1(i2,  j)
    print *,xsgp1(i1+1,j),xsgp1(i2-1,j)
    print *,xsgp1(i1+2,j),xsgp1(i2-2,j)

    print *,'replicate Y across seam on SG'
    print *,ysgp1(i1-2,j),ysgp1(i2+2,j)
    print *,ysgp1(i1-1,j),ysgp1(i2+1,j)
    print *,ysgp1(i1,  j),ysgp1(i2,  j)
    print *,ysgp1(i1+1,j),ysgp1(i2-1,j)
    print *,ysgp1(i1+2,j),ysgp1(i2-2,j)

    !---------------------------------------------------------------------
    ! rotation angle on supergrid vertices
    ! lonB: x(i-1,j-1) has same relationship to x(i,j) on SG as
    !       geolonT(i,j) has to geolonBu(i,j) on the reduced grid
    !---------------------------------------------------------------------

    ! constants as defined in MOM
    pi_720deg = atan(1.0) / 180.0
    len_lon = 360.0
    do j=1,ny ; do i=1,nx-1
       do n=1,2 ; do m=1,2
          lonb(m,n) = modulo_around_point(xsgp1(i+m-2,j+n-2), xsgp1(i-1,j-1), len_lon)
       enddo; enddo
       lon_scale    = cos(pi_720deg*(ysgp1(i-1,j-1) + ysgp1(i+1,j-1) + &
            ysgp1(i-1,j+1) + ysgp1(i+1,j+1)) )
       angq(i,j)    = atan2(lon_scale*((lonb(1,2) - lonb(2,1)) + (lonb(2,2) - lonb(1,1))), &
            ysgp1(i-1,j+1) + ysgp1(i+1,j+1) - &
            ysgp1(i-1,j-1) - ysgp1(i+1,j-1) )
    enddo; enddo

    !check
    if(debug) then
       j = ny
       i1 = ipolesg(1); i2 = ipolesg(2)-(ipolesg(1)-i1)
       print *,'angq along seam on SG'
       print *,angq(i1-2,j),angq(i2+2,j)
       print *,angq(i1-1,j),angq(i2+1,j)
       print *,angq(i1,  j),angq(i2,  j)
       print *,angq(i1+1,j),angq(i2-1,j)
       print *,angq(i1+2,j),angq(i2-2,j)
    end if

  end subroutine find_angq

  subroutine find_ang

    ! 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=1,nj; do i = 1,ni
       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 = ni
          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 = ni
       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