direct_transform Subroutine

  subroutine direct_transform(stretch_factor, lon_p, lat_p, lon, lat)
    real(ESMF_KIND_R8),    intent(in):: stretch_factor !< Stretching factor
    real(ESMF_KIND_R8),    intent(in):: lon_p, lat_p   !< center location of the target face, radian
!  0 <= lon <= 2*pi ;    -pi/2 <= lat <= pi/2
    real(ESMF_KIND_R8), intent(inout) :: lon(:,:), lat(:,:)
!
    real(ESMF_KIND_R8)  :: lat_t, sin_p, cos_p, sin_lat, cos_lat, sin_o, p2, two_pi
    real(ESMF_KIND_R8)  :: c2p1, c2m1
    integer:: i, j

    p2 = 0.5d0*pi
    two_pi = 2.d0*pi

    c2p1 = 1.d0 + stretch_factor*stretch_factor
    c2m1 = 1.d0 - stretch_factor*stretch_factor

    sin_p = sin(lat_p)
    cos_p = cos(lat_p)

    do j=1,size(lon,2)
       do i=1,size(lon,1)
          if ( abs(c2m1) > 1.d-7 ) then
               sin_lat = sin(lat(i,j))
               lat_t = asin( (c2m1+c2p1*sin_lat)/(c2p1+c2m1*sin_lat) )
          else         ! no stretching
               lat_t = lat(i,j)
          endif
          sin_lat = sin(lat_t)
          cos_lat = cos(lat_t)
            sin_o = -(sin_p*sin_lat + cos_p*cos_lat*cos(lon(i,j)))
          if ( (1.-abs(sin_o)) < 1.d-7 ) then    ! poles
               lon(i,j) = 0.d0
               lat(i,j) = sign( p2, sin_o )
          else
               lat(i,j) = asin( sin_o )
               lon(i,j) = lon_p + atan2( -cos_lat*sin(lon(i,j)),   &
                          -sin_lat*cos_p+cos_lat*sin_p*cos(lon(i,j)))
               if ( lon(i,j) < 0.d0 ) then
                    lon(i,j) = lon(i,j) + two_pi
               elseif( lon(i,j) >= two_pi ) then
                    lon(i,j) = lon(i,j) - two_pi
               endif
          endif
       enddo
    enddo

  end subroutine direct_transform