Re: a suggestion about wrf_latlon_to_ij

Hi Yunfei,

Thanks for your suggestion. I'll pass this onto the WRF developers,
since they are the ones who gave us the Fortran code for the NCL
version of this function.

--Mary

On Wed, 4 Jun 2008, yunfei zhang wrote:

> 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 - 10:28:17 MDT