var/da/da_tools/da_set_lc.inc

   1 subroutine da_set_lc(proj)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: Initialize the remaining items in the proj structure for a
   5    ! lambert conformal grid.
   6    !-----------------------------------------------------------------------
   7
   8    implicit none
   9
  10    type(proj_info), intent(inout) :: proj
  11
  12    real :: arg
  13    real :: deltalon1
  14    real :: tl1r
  15    real :: ctl1r
  16
  17    if (trace_use_dull) call da_trace_entry("da_set_lc")
  18
  19    ! Compute cone factor
  20    call da_lc_cone(proj%truelat1, proj%truelat2, proj%cone)
  21    if (print_detail_map) then
  22       write(unit=stdout, fmt='(A,F8.6)') 'Computed cone factor: ', proj%cone
  23    end if
  24    ! Compute longitude differences and ensure we stay out of the
  25    ! forbidden "cut zone"
  26    deltalon1 = proj%lon1 - proj%stdlon
  27    if (deltalon1 .gt. +180.0) deltalon1 = deltalon1 - 360.0
  28    if (deltalon1 .lt. -180.0) deltalon1 = deltalon1 + 360.0
  29
  30    ! Convert truelat1 to radian and compute COS for later use
  31    tl1r = proj%truelat1 * rad_per_deg
  32    ctl1r = COS(tl1r)
  33
  34    ! Compute the radius to our known lower-left (SW) corner
  35    proj%rsw = proj%rebydx * ctl1r/proj%cone * &
  36            (TAN((90.0*proj%hemi-proj%lat1)*rad_per_deg/2.0) / &
  37             TAN((90.0*proj%hemi-proj%truelat1)*rad_per_deg/2.0))**proj%cone
  38
  39    ! Find pole point
  40    arg = proj%cone*(deltalon1*rad_per_deg)
  41    proj%polei = 1.0 - proj%hemi * proj%rsw * Sin(arg)
  42    proj%polej = 1.0 + proj%rsw * COS(arg)
  43    if (print_detail_map) then
  44       write(unit=stdout,fmt='(A,2F10.3)') 'Computed pole i/j = ', proj%polei, proj%polej
  45    end if
  46
  47    if (trace_use_dull) call da_trace_exit("da_set_lc")
  48
  49 end subroutine da_set_lc
  50
  51