a suggestion about wrf_latlon_to_ij

recently, I used the function wrf_latlon_to_ij.
when the data points become very large, for example 300*500, this function
needs long time to process.
I read the source code ./ni/src/lib/nfpfort/wrf_user.f, and found this
function use the least distance method.
I think it should use projection information to transform a latitude
longitude pair into a i j pair on grid,
because the wrf grid is a regular projection.
I give a example fortran program as following:

SUBROUTINE llij_lc( lat, lon, proj, i, j)

  ! Subroutine to compute the geographical latitude and longitude values
  ! to the cartesian x/y on a Lambert Conformal projection.

    IMPLICIT NONE

    ! Input Args
    REAL, INTENT(IN) :: lat ! Latitude (-90->90 deg N)
    REAL, INTENT(IN) :: lon ! Longitude (-180->180 E)
    TYPE(proj_info),INTENT(IN) :: proj ! Projection info structure

    ! Output Args
    REAL, INTENT(OUT) :: i ! Cartesian X coordinate
    REAL, INTENT(OUT) :: j ! Cartesian Y coordinate

    ! Locals
    REAL :: arg
    REAL :: deltalon
    REAL :: tl1r
    REAL :: rm
    REAL :: ctl1r

    ! BEGIN CODE

    ! Compute deltalon between known longitude and standard lon and ensure
    ! it is not in the cut zone
    deltalon = lon - proj%stdlon
    IF (deltalon .GT. +180.) deltalon = deltalon - 360.
    IF (deltalon .LT. -180.) deltalon = deltalon + 360.

    ! Convert truelat1 to radian and compute COS for later use
    tl1r = proj%truelat1 * rad_per_deg
    ctl1r = COS(tl1r)

    ! Radius to desired point
    rm = proj%rebydx * ctl1r/proj%cone * &
         (TAN((90.*proj%hemi-lat)*rad_per_deg/2.) / &
          TAN((90.*proj%hemi-proj%truelat1)*rad_per_deg/2.))**proj%cone

    arg = proj%cone*(deltalon*rad_per_deg)
    i = proj%polei + proj%hemi * rm * SIN(arg)
    j = proj%polej - rm * COS(arg)

    ! Finally, if we are in the southern hemisphere, flip the i/j
    ! values to a coordinate system where (1,1) is the SW corner
    ! (what we assume) which is different than the original NCEP
    ! algorithms which used the NE corner as the origin in the
    ! southern hemisphere (left-hand vs. right-hand coordinate?)
    IF (proj%hemi .EQ. -1.) THEN
      i = 2. - i
      j = 2. - j
    ENDIF
    RETURN
  END SUBROUTINE llij_lc

_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Jun 04 2008 - 04:00:39 MDT