var/da/da_tools/da_map_set.inc

   1 subroutine da_map_set(proj_code,lat1,lon1,knowni,knownj,dx,stdlon,truelat1,truelat2,latinc,loninc,proj)
   2    ! Given a partially filled proj_info structure, this routine computes
   3    ! polei, polej, rsw, and cone (if LC projection) to complete the
   4    ! structure.  This allows us to eliminate redundant calculations when
   5    ! calling the coordinate conversion routines multiple times for the
   6    ! same map.
   7    ! This will generally be the first routine called when a user wants
   8    ! to be able to use the coordinate conversion routines, and it
   9    ! will call the appropriate subroutines based on the
  10    ! proj%code which indicates which projection type  this is.
  11
  12    implicit none
  13
  14    integer,         intent(in)  :: proj_code
  15    real,            intent(in)  :: lat1
  16    real,            intent(in)  :: lon1
  17    real,            intent(in)  :: dx
  18    real,            intent(in)  :: stdlon
  19    real,            intent(in)  :: truelat1
  20    real,            intent(in)  :: truelat2
  21    real,            intent(in)  :: knowni , knownj
  22    real,            intent(in)  :: latinc, loninc
  23    type(proj_info), intent(out) :: proj
  24
  25    if (trace_use_dull) call da_trace_entry("da_map_set")
  26
  27    ! First, check for validity of mandatory variables in proj
  28    if (ABS(lat1) > 90.0) then
  29       call da_error(__FILE__,__LINE__, &
  30          (/"Latitude of origin corner required as follows: -90N <= lat1 < = 90N"/))
  31    end if
  32    if (ABS(lon1) > 180.0) then
  33       call da_error(__FILE__,__LINE__, &
  34          (/"Longitude of origin required as follows: -180E <= lon1 <= 180W"/))
  35    end if
  36    if ((dx .LE. 0.0).AND.(proj_code .NE. PROJ_LATLON)) then
  37       call da_error(__FILE__,__LINE__, &
  38          (/"Require grid spacing (dx) in meters be positive!"/))
  39    end if
  40    if ((ABS(stdlon) > 180.0).AND.(proj_code .NE. PROJ_MERC)) then
  41       call da_error(__FILE__,__LINE__, &
  42          (/"Need orientation longitude (stdlon) as: -180E <= lon1 <= 180W"/))
  43    end if
  44    if (proj%code .NE. PROJ_LATLON .and. ABS(truelat1)>90.0) then
  45       call da_error(__FILE__,__LINE__, &
  46          (/"Set true latitude 1 for all projections!"/))
  47    end if
  48
  49    call da_map_init(proj)
  50    proj%code  = proj_code
  51    proj%lat1 = lat1
  52    proj%lon1 = lon1
  53    proj%knowni = knowni
  54    proj%knownj = knownj
  55    proj%dx    = dx
  56    proj%stdlon = stdlon
  57    proj%truelat1 = truelat1
  58    proj%truelat2 = truelat2
  59    proj%latinc   = latinc
  60    proj%loninc   = loninc
  61    if (proj%code .NE. PROJ_LATLON) then
  62       proj%dx = dx
  63       if (truelat1 < 0.0) then
  64          proj%hemi = -1.0
  65       else
  66          proj%hemi = 1.0
  67       end if
  68       proj%rebydx = earth_radius_m / dx
  69    end if
  70    pick_proj: select case(proj%code)
  71
  72       case(PROJ_PS)
  73          if (print_detail_map) then
  74             write(unit=stdout,fmt='(A)') 'Setting up POLAR STEREOGRAPHIC map...'
  75          end if
  76          call da_set_ps(proj)
  77
  78       case(PROJ_LC)
  79          if (print_detail_map) then
  80             write(unit=stdout,fmt='(A)') 'Setting up LAMBERT CONFORMAL map...'
  81          end if
  82          if (ABS(proj%truelat2) > 90.0) then
  83             if (print_detail_map) then
  84                write(unit=stdout,fmt='(A)') 'Second true latitude not set, assuming a tangent'
  85                write(unit=stdout,fmt='(A,F10.3)') 'projection at truelat1: ', proj%truelat1
  86             end if
  87             proj%truelat2=proj%truelat1
  88          end if
  89          call da_set_lc(proj)
  90
  91       case (PROJ_MERC)
  92          if (print_detail_map) then
  93             write(unit=stdout,fmt='(A)') 'Setting up MERCATOR map...'
  94          end if
  95          call da_set_merc(proj)
  96
  97       case (PROJ_LATLON)
  98          if (print_detail_map) then
  99             write(unit=stdout,fmt='(A)') 'Setting up CYLINDRICAL EQUIDISTANT LATLON map...'
 100          end if
 101          ! Convert lon1 to 0->360 notation
 102          if (proj%lon1 < 0.0) proj%lon1 = proj%lon1 + 360.0
 103
 104          proj%cone = 1.0
 105       case default
 106          write(unit=message(1),fmt='(A,I2)') 'Unknown projection code: ', proj%code
 107          call da_error(__FILE__,__LINE__,message(1:1))
 108    end select pick_proj
 109    proj%init = .true.
 110
 111    if (trace_use_dull) call da_trace_exit("da_map_set")
 112
 113 end subroutine da_map_set
 114
 115