var/da/da_tools/da_llxy_lc_new.inc

   1 subroutine da_llxy_lc_new(proj, info)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: compute the geographical latitude and longitude values
   5    ! to the cartesian x/y on a Lambert Conformal projection.
   6    !-----------------------------------------------------------------------
   7
   8    implicit none
   9
  10    type(proj_info), intent(in)    :: proj     ! Projection info structure
  11    type(infa_type), intent(inout) :: info
  12
  13
  14    real    :: tl1r
  15 !   real    :: temp1
  16 !   real    :: temp2
  17    real    :: ctl1r
  18    real    :: deltalon, rm, arg
  19    integer :: n
  20
  21    if (trace_use) call da_trace_entry("da_llxy_lc_new")
  22
  23 ! FAST
  24
  25    ! Convert truelat1 to radian and compute COS for later use
  26 !   tl1r = proj%truelat1 * rad_per_deg
  27 !   ctl1r = COS(tl1r)
  28 !   temp1 = TAN((90.0*proj%hemi-proj%truelat1)*rad_per_deg/2.0)
  29 !   temp2 = proj%rebydx * ctl1r/proj%cone
  30
  31    ! Compute deltalon between known longitude and standard lon and ensure
  32    ! it is not in the cut zone
  33
  34 !   where (lon - proj%stdlon > +180.0)
  35 !      x = proj%polei + proj%hemi * (temp2 * (TAN((90.0*proj%hemi-lat)*rad_per_deg/2.0) / temp1)**proj%cone &
  36 !         * SIN(proj%cone*((lon - proj%stdlon-360.0)*rad_per_deg))
  37 !      y = proj%polej - (temp2 * (TAN((90.0*proj%hemi-lat)*rad_per_deg/2.0) / temp1)**proj%cone &
  38 !         * COS(proj%cone*((lon - proj%stdlon-360.0)*rad_per_deg))
  39 !   elsewhere  (lon - proj%stdlon - -180.0)
  40 !      x = proj%polei + proj%hemi * (temp2 * (TAN((90.0*proj%hemi-lat)*rad_per_deg/2.0) / temp1) &
  41 !         * SIN(proj%cone*((lon - proj%stdlon+360.0)*rad_per_deg))
  42 !      y = proj%polej - (temp2 * (TAN((90.0*proj%hemi-lat)*rad_per_deg/2.0) / temp1)**proj%cone &
  43 !         * COS(proj%cone*((lon - proj%stdlon+360.0)*rad_per_deg))
  44 !   else
  45 !      x = proj%polei + proj%hemi * (temp2 * (TAN((90.0*proj%hemi-lat)*rad_per_deg/2.0) / temp1) &
  46 !         * SIN(proj%cone*((lon - proj%stdlon)*rad_per_deg))
  47 !      y = proj%polej - (temp2 * (TAN((90.0*proj%hemi-lat)*rad_per_deg/2.0) / temp1)**proj%cone &
  48 !         * COS(proj%cone*((lon - proj%stdlon)*rad_per_deg))
  49 !   end where
  50
  51    ! Finally, if we are in the southern hemisphere, flip the i/j
  52    ! values to a coordinate system where (1,1) is the SW corner
  53    ! (what we assume) which is different than the original NCEP
  54    ! algorithms which used the NE corner as the origin in the
  55    ! southern hemisphere (left-hand vs. right-hand coordinate?)
  56 !   if (proj%hemi == -1.0) then
  57 !      x(:,:) = 2.0 - i(:,:)
  58 !      y(:,:) = 2.0 - j(:,:)
  59 !   end if
  60
  61 ! SLOW
  62
  63    do n=lbound(info%lat,2), ubound(info%lat,2)
  64       ! Compute deltalon between known longitude and standard lon and ensure
  65       ! it is not in the cut zone
  66       deltalon = info%lon(1,n) - proj%stdlon
  67       if (deltalon > +180.0) deltalon = deltalon - 360.0
  68       if (deltalon < -180.0) deltalon = deltalon + 360.0
  69
  70       ! Convert truelat1 to radian and compute COS for later use
  71       tl1r = proj%truelat1 * rad_per_deg
  72       ctl1r = COS(tl1r)
  73
  74       ! Radius to desired point
  75       rm = proj%rebydx * ctl1r/proj%cone * &
  76           (TAN((90.0*proj%hemi-info%lat(1,n))*rad_per_deg/2.0) / &
  77            TAN((90.0*proj%hemi-proj%truelat1)*rad_per_deg/2.0))**proj%cone
  78
  79       arg = proj%cone*(deltalon*rad_per_deg)
  80       info%x(:,n) = proj%polei + proj%hemi * rm * Sin(arg)
  81       info%y(:,n) = proj%polej - rm * COS(arg)
  82
  83       ! Finally, if we are in the southern hemisphere, flip the i/j
  84       ! values to a coordinate system where (1,1) is the SW corner
  85       ! (what we assume) which is different than the original NCEP
  86       ! algorithms which used the NE corner as the origin in the
  87       ! southern hemisphere (left-hand vs. right-hand coordinate?)
  88       if (proj%hemi == -1.0) then
  89          info%x(:,n) = 2.0 - info%x(:,n)
  90          info%y(:,n) = 2.0 - info%y(:,n)
  91       end if
  92    end do
  93
  94    if (trace_use) call da_trace_exit("da_llxy_lc_new")
  95
  96 end subroutine da_llxy_lc_new
  97
  98