var/da/da_tools/da_llxy_merc_new.inc

   1 subroutine da_llxy_merc_new(proj, info)
   2
   3    !-----------------------------------------------------------------------
   4    ! Purpose: Compute x/y coordinate from lat lon for mercator projection
   5    !-----------------------------------------------------------------------
   6
   7    implicit none
   8
   9    type(proj_info), intent(in)    :: proj
  10    type(infa_type), intent(inout) :: info
  11
  12    real :: deltalon
  13    integer :: n
  14
  15 if (trace_use) call da_trace_entry("da_llxy_merc_new")
  16
  17 ! FAST
  18
  19 !   where (lon(:,:) - proj%lon1 < -180.0)
  20 !      x(:,:) = 1.0 + (lon(:,:) - proj%lon1 + 360.0)/(proj%dlon*deg_per_rad))
  21 !   elsewhere (lon(:,:) - proj%lon1 > 180.0)
  22 !      x(:,:) = 1.0 + (lon(:,:) - proj%lon1 - 360.0)/(proj%dlon*deg_per_rad))
  23 !   else
  24 !      x(:,:) = 1.0 + (lon(:,:) - proj%lon1)/(proj%dlon*deg_per_rad))
  25 !   end where
  26
  27 !   y(:,:) = 1.0 + (ALOG(TAN(0.5*((lat(:,: + 90.0) * rad_per_deg)))) / proj%dlon - proj%rsw
  28
  29 ! SLOW
  30
  31    do n=lbound(info%lat,2),ubound(info%lat,2)
  32       deltalon = info%lon(1,n) - proj%lon1
  33       if (deltalon < -180.0) deltalon = deltalon + 360.0
  34       if (deltalon > 180.0) deltalon = deltalon - 360.0
  35       info%x(:,n) = 1.0 + (deltalon/(proj%dlon*deg_per_rad))
  36       info%y(:,n) = 1.0 + (ALOG(TAN(0.5*((info%lat(1,n) + 90.0) * rad_per_deg)))) / proj%dlon - proj%rsw
  37    end do
  38
  39 if (trace_use) call da_trace_exit("da_llxy_merc_new")
  40
  41
  42 end subroutine da_llxy_merc_new
  43
  44