<HTML> <BODY BGCOLOR=#ccccdd LINK=#0000aa VLINK=#0000ff ALINK=#ff0000 ><BASE TARGET="bottom_target"><PRE>
<A NAME='DA_VERIF_TOOLS'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_VERIF_TOOLS' TARGET='top_target'><IMG SRC="../../gif/bar_purple.gif" border=0></A>
module da_verif_tools 1
implicit none
! Define some private constants
real, private, parameter :: pi = 3.1415926
real, private, parameter :: deg_per_rad = 180.0/pi
real, private, parameter :: rad_per_deg = pi / 180.0
logical, private, parameter :: print_detail_map = .false.
integer, private, parameter :: stdout = 6
! Mean Earth Radius in m. The value below is consistent
! with NCEP's routines and grids.
! real, public , parameter :: earth_radius_m = 6371200.0 ! from Brent
real, public , parameter :: earth_radius_m = 6370000.0
real, public , parameter :: radians_per_degree = pi / 180.0
! Define public parameters
! Projection codes for proj_info structure:
integer, public, parameter :: PROJ_LATLON = 0
integer, public, parameter :: PROJ_MERC = 1
integer, public, parameter :: PROJ_LC = 3
integer, public, parameter :: PROJ_PS = 5
! Define data structures to define various projections
type proj_info
integer :: code ! integer code for projection type
real :: lat1 ! SW latitude (1,1) in degrees (-90->90N)
real :: lon1 ! SW longitude (1,1) in degrees (-180->180E)
real :: dx ! Grid spacing in meters at truelats, used
! only for ps, lc, and merc projections
real :: dlat ! Lat increment for lat/lon grids
real :: dlon ! Lon increment for lat/lon grids
real :: stdlon ! Longitude parallel to y-axis (-180->180E)
real :: truelat1 ! First true latitude (all projections)
real :: truelat2 ! Second true lat (LC only)
real :: hemi ! 1 for NH, -1 for SH
real :: cone ! Cone factor for LC projections
real :: polei ! Computed i-location of pole point
real :: polej ! Computed j-location of pole point
real :: rsw ! Computed radius to SW corner
real :: rebydx ! Earth radius divided by dx
real :: knowni ! X-location of known lat/lon
real :: knownj ! Y-location of known lat/lon
real :: latinc ! Latitude increments in degrees
real :: loninc ! Longitude increments in degrees
logical :: init ! Flag to indicate if this struct is
! ready for use
end type proj_info
type(proj_info) :: map_info
character(len=10000),private :: message(50)
contains
<A NAME='DA_LLXY_LATLON'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_LLXY_LATLON' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_llxy_latlon(lat, lon, proj, x, y) 2,2
!-----------------------------------------------------------------------
! Purpose: Compute the x/y location of a lat/lon on a LATLON grid.
!-----------------------------------------------------------------------
implicit none
real, intent(in) :: lat
real, intent(in) :: lon
type(proj_info), intent(in) :: proj
real, intent(out) :: x
real, intent(out) :: y
real :: deltalat
real :: deltalon
real :: lon360
real :: latinc
real :: loninc
! Extract the latitude and longitude increments for this grid
! (e.g., 2.5 deg for NCEP reanalysis) from the proj structure, where
! loninc is saved in the stdlon tag and latinc is saved in truelat1
latinc = proj%truelat1
loninc = proj%stdlon
! Compute deltalat and deltalon as the difference between the input
! lat/lon and the origin lat/lon
deltalat = lat - proj%lat1
! To account for issues around the dateline, convert the incoming
! longitudes to be 0->360.0
if (lon < 0) then
lon360 = lon + 360.0
else
lon360 = lon
end if
deltalon = lon360 - proj%lon1
! Compute x/y
x = deltalon/loninc + 1.0
y = deltalat/latinc + 1.0
end subroutine da_llxy_latlon
<A NAME='DA_LLXY_LC'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_LLXY_LC' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_llxy_lc(lat, lon, proj, x, y) 2,2
!-----------------------------------------------------------------------
! Purpose: compute the geographical latitude and longitude values
! to the cartesian x/y on a Lambert Conformal projection.
!-----------------------------------------------------------------------
implicit none
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
real, intent(out) :: x ! Cartesian X coordinate
real, intent(out) :: y ! Cartesian Y coordinate
real :: arg
real :: deltalon
real :: tl1r
real :: rm
real :: ctl1r
! Compute deltalon between known longitude and standard lon and ensure
! it is not in the cut zone
deltalon = lon - proj%stdlon
if (deltalon > +180.0) deltalon = deltalon - 360.0
if (deltalon < -180.0) deltalon = deltalon + 360.0
! 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.0*proj%hemi-lat)*rad_per_deg/2.0) / &
TAN((90.0*proj%hemi-proj%truelat1)*rad_per_deg/2.0))**proj%cone
arg = proj%cone*(deltalon*rad_per_deg)
x = proj%polei + proj%hemi * rm * Sin(arg)
y = 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 == -1.0) then
x = 2.0 - x
y = 2.0 - y
end if
end subroutine da_llxy_lc
<A NAME='DA_LLXY_MERC'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_LLXY_MERC' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_llxy_merc(lat, lon, proj, x, y) 2,2
!-----------------------------------------------------------------------
! Purpose: Compute x,y coordinate from lat lon for mercator projection
!-----------------------------------------------------------------------
implicit none
real, intent(in) :: lat
real, intent(in) :: lon
type(proj_info),intent(in) :: proj
real,intent(out) :: x
real,intent(out) :: y
real :: deltalon
deltalon = lon - proj%lon1
if (deltalon < -180.0) deltalon = deltalon + 360.0
if (deltalon > 180.0) deltalon = deltalon - 360.0
x = 1.0 + (deltalon/(proj%dlon*deg_per_rad))
y = 1.0 + (ALOG(TAN(0.5*((lat + 90.0) * rad_per_deg)))) / &
proj%dlon - proj%rsw
end subroutine da_llxy_merc
<A NAME='DA_LLXY_PS'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_LLXY_PS' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_llxy_ps(lat,lon,proj,x,y) 2,2
!-----------------------------------------------------------------------
! Purpose: Given latitude (-90 to 90), longitude (-180 to 180), and the
! standard polar-stereographic projection information via the
! public proj structure, this routine returns the x/y indices which
! if within the domain range from 1->nx and 1->ny, respectively.
!-----------------------------------------------------------------------
implicit none
real, intent(in) :: lat
real, intent(in) :: lon
type(proj_info),intent(in) :: proj
real, intent(out) :: x !(x-index)
real, intent(out) :: y !(y-index)
real :: reflon
real :: scale_top
real :: ala
real :: alo
real :: rm
reflon = proj%stdlon + 90.0
! Compute numerator term of map scale factor
scale_top = 1.0 + proj%hemi * Sin(proj%truelat1 * rad_per_deg)
! Find radius to desired point
ala = lat * rad_per_deg
rm = proj%rebydx * COS(ala) * scale_top/(1.0 + proj%hemi *Sin(ala))
alo = (lon - reflon) * rad_per_deg
x = proj%polei + rm * COS(alo)
y = proj%polej + proj%hemi * rm * Sin(alo)
end subroutine da_llxy_ps
<A NAME='DA_LLXY_WRF'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_LLXY_WRF' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_llxy_wrf(proj, lat, lon, x, y) 5,12
!-----------------------------------------------------------------------
! Purpose: Converts input lat/lon values to the cartesian (x, y) value
! for the given projection.
!-----------------------------------------------------------------------
implicit none
type(proj_info), intent(in) :: proj
real, intent(in) :: lat
real, intent(in) :: lon
real, intent(out) :: x
real, intent(out) :: y
if (.NOT.proj%init) then
write(stdout,*) "da_llxy_wrf.inc"// &
"You have not called map_set for this projection!"
stop
end if
select case(proj%code)
case(PROJ_LATLON)
call da_llxy_latlon
(lat,lon,proj,x,y)
case(PROJ_MERC)
call da_llxy_merc(lat,lon,proj,x,y)
x = x + proj%knowni - 1.0
y = y + proj%knownj - 1.0
case(PROJ_PS)
call da_llxy_ps(lat,lon,proj,x,y)
case(PROJ_LC)
call da_llxy_lc(lat,lon,proj,x,y)
x = x + proj%knowni - 1.0
y = y + proj%knownj - 1.0
case default
write(unit=message(1),fmt='(A,I2)') &
'Unrecognized map projection code: ', proj%code
write(stdout,*) "da_llxy_wrf.inc"//message(1:1)
stop
end select
end subroutine da_llxy_wrf
<A NAME='DA_XYLL'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_XYLL' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_xyll(proj, xx, yy, lat, lon) 2,12
!-----------------------------------------------------------------------
! Purpose: Computes geographical latitude and longitude for a given (i,j)
! point in a grid with a projection of proj
!-----------------------------------------------------------------------
implicit none
type(proj_info), intent(in) :: proj
real, intent(in) :: xx
real, intent(in) :: yy
real, intent(out) :: lat
real, intent(out) :: lon
real :: x, y
if (.NOT.proj%init) then
write(stdout,*) "da_xyll.inc"// &
"You have not called map_set for this projection!"
stop
end if
x = xx
y = yy
select case (proj%code)
case (PROJ_LATLON)
call da_xyll_latlon(x, y, proj, lat, lon)
case (PROJ_MERC)
x = xx - proj%knowni + 1.0
y = yy - proj%knownj + 1.0
call da_xyll_merc(x, y, proj, lat, lon)
case (PROJ_PS)
call da_xyll_ps(x, y, proj, lat, lon)
case (PROJ_LC)
x = xx - proj%knowni + 1.0
y = yy - proj%knownj + 1.0
call da_xyll_lc(x, y, proj, lat, lon)
case default
write(unit=message(1),fmt='(A,I2)') &
"Unrecognized map projection code: ", proj%code
write(stdout,*) "da_xyll.inc"//message(1:1)
stop
end select
end subroutine da_xyll
<A NAME='DA_XYLL_LATLON'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_XYLL_LATLON' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_xyll_latlon(x, y, proj, lat, lon) 2,2
!-----------------------------------------------------------------------
! Purpose: Compute the lat/lon location of an x/y on a LATLON grid.
!-----------------------------------------------------------------------
implicit none
real, intent(in) :: x
real, intent(in) :: y
type(proj_info), intent(in) :: proj
real, intent(out) :: lat
real, intent(out) :: lon
real :: deltalat
real :: deltalon
real :: latinc
real :: loninc
! Extract the latitude and longitude increments for this grid
! (e.g., 2.5 deg for NCEP reanalysis) from the proj structure, where
! loninc is saved in the stdlon tag and latinc is saved in truelat1
latinc = proj%truelat1
loninc = proj%stdlon
! Compute deltalat and deltalon
deltalat = (x-1.0)*latinc
deltalon = (y-1.0)*loninc
lat = proj%lat1 + deltalat
lon = proj%lon1 + deltalon
if ((ABS(lat) > 90.0).OR.(ABS(deltalon) > 360.0)) then
! Off the earth for this grid
lat = -999.0
lon = -999.0
else
lon = lon + 360.0
lon = MOD(lon,360.0)
if (lon > 180.0) lon = lon -360.0
end if
end subroutine da_xyll_latlon
<A NAME='DA_XYLL_LC'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_XYLL_LC' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_xyll_lc(x, y, proj, lat, lon) 2,2
! subroutine da_to convert from the (x,y) cartesian coordinate to the
! geographical latitude and longitude for a Lambert Conformal projection.
implicit none
real, intent(in) :: x ! Cartesian X coordinate
real, intent(in) :: y ! Cartesian Y coordinate
type(proj_info),intent(in) :: proj ! Projection info structure
real, intent(out) :: lat ! Latitude (-90->90 deg N)
real, intent(out) :: lon ! Longitude (-180->180 E)
real :: inew
real :: jnew
real :: r
real :: chi,chi1,chi2
real :: r2
real :: xx
real :: yy
chi1 = (90.0 - proj%hemi*proj%truelat1)*rad_per_deg
chi2 = (90.0 - proj%hemi*proj%truelat2)*rad_per_deg
! See if we are in the southern hemispere and flip the indices
! if we are.
if (proj%hemi == -1.0) then
inew = -x + 2.0
jnew = -y + 2.0
else
inew = x
jnew = y
end if
! Compute radius**2 to i/j location
xx = inew - proj%polei
yy = proj%polej - jnew
r2 = (xx*xx + yy*yy)
r = sqrt(r2)/proj%rebydx
! Convert to lat/lon
if (r2 == 0.0) then
lat = proj%hemi * 90.0
lon = proj%stdlon
else
! Longitude
lon = proj%stdlon + deg_per_rad * ATAN2(proj%hemi*xx,yy)/proj%cone
lon = MOD(lon+360.0, 360.0)
! Latitude. Latitude determined by solving an equation adapted
! from:
! Maling, D.H., 1973: Coordinate Systems and Map Projections
! Equations #20 in Appendix I.
if (chi1 == chi2) then
chi = 2.0*ATAN((r/TAN(chi1))**(1.0/proj%cone) * TAN(chi1*0.5))
else
chi = 2.0*ATAN((r*proj%cone/Sin(chi1))**(1.0/proj%cone) * TAN(chi1*0.5))
end if
lat = (90.0-chi*deg_per_rad)*proj%hemi
end if
if (lon > +180.0) lon = lon - 360.0
if (lon < -180.0) lon = lon + 360.0
end subroutine da_xyll_lc
<A NAME='DA_XYLL_MERC'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_XYLL_MERC' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_xyll_merc(x, y, proj, lat, lon) 2,2
!-----------------------------------------------------------------------
! Compute the lat/lon from i/j for mercator projection
!-----------------------------------------------------------------------
implicit none
real,intent(in) :: x
real,intent(in) :: y
type(proj_info),intent(in) :: proj
real, intent(out) :: lat
real, intent(out) :: lon
lat = 2.0*ATAN(EXP(proj%dlon*(proj%rsw + y-1.0)))*deg_per_rad - 90.0
lon = (x-1.0)*proj%dlon*deg_per_rad + proj%lon1
if (lon > 180.0) lon = lon - 360.0
if (lon < -180.0) lon = lon + 360.0
end subroutine da_xyll_merc
<A NAME='DA_XYLL_PS'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_XYLL_PS' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_xyll_ps(x, y, proj, lat, lon) 2,2
! This is the inverse subroutine da_of llij_ps. It returns the
! latitude and longitude of an x/y point given the projection info
! structure.
implicit none
real, intent(in) :: x ! Column
real, intent(in) :: y ! Row
type (proj_info), intent(in) :: proj
real, intent(out) :: lat ! -90 -> 90 North
real, intent(out) :: lon ! -180 -> 180 East
real :: reflon
real :: scale_top
real :: xx,yy
real :: gi2, r2
real :: arccos
! Compute the reference longitude by rotating 90 degrees to the east
! to find the longitude line parallel to the positive x-axis.
reflon = proj%stdlon + 90.0
! Compute numerator term of map scale factor
scale_top = 1.0 + proj%hemi * Sin(proj%truelat1 * rad_per_deg)
! Compute radius to point of interest
xx = x - proj%polei
yy = (y - proj%polej) * proj%hemi
r2 = xx**2 + yy**2
! Now the magic code
if (r2 == 0.0) then
lat = proj%hemi * 90.0
lon = reflon
else
gi2 = (proj%rebydx * scale_top)**2.0
lat = deg_per_rad * proj%hemi * ASin((gi2-r2)/(gi2+r2))
arccos = ACOS(xx/sqrt(r2))
if (yy > 0) then
lon = reflon + deg_per_rad * arccos
else
lon = reflon - deg_per_rad * arccos
end if
end if
! Convert to a -180 -> 180 East convention
if (lon > 180.0) lon = lon - 360.0
if (lon < -180.0) lon = lon + 360.0
end subroutine da_xyll_ps
<A NAME='DA_SET_LC'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_SET_LC' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_set_lc(proj) 2,4
!-----------------------------------------------------------------------
! Purpose: Initialize the remaining items in the proj structure for a
! lambert conformal grid.
!-----------------------------------------------------------------------
implicit none
type(proj_info), intent(inout) :: proj
real :: arg
real :: deltalon1
real :: tl1r
real :: ctl1r
! Compute cone factor
call da_lc_cone(proj%truelat1, proj%truelat2, proj%cone)
if (print_detail_map) then
write(unit=stdout, fmt='(A,F8.6)') 'Computed cone factor: ', proj%cone
end if
! Compute longitude differences and ensure we stay out of the
! forbidden "cut zone"
deltalon1 = proj%lon1 - proj%stdlon
if (deltalon1 .gt. +180.0) deltalon1 = deltalon1 - 360.0
if (deltalon1 .lt. -180.0) deltalon1 = deltalon1 + 360.0
! Convert truelat1 to radian and compute COS for later use
tl1r = proj%truelat1 * rad_per_deg
ctl1r = COS(tl1r)
! Compute the radius to our known lower-left (SW) corner
proj%rsw = proj%rebydx * ctl1r/proj%cone * &
(TAN((90.0*proj%hemi-proj%lat1)*rad_per_deg/2.0) / &
TAN((90.0*proj%hemi-proj%truelat1)*rad_per_deg/2.0))**proj%cone
! Find pole point
arg = proj%cone*(deltalon1*rad_per_deg)
proj%polei = 1.0 - proj%hemi * proj%rsw * Sin(arg)
proj%polej = 1.0 + proj%rsw * COS(arg)
if (print_detail_map) then
write(unit=stdout,fmt='(A,2F10.3)') 'Computed pole i/j = ', proj%polei, proj%polej
end if
end subroutine da_set_lc
<A NAME='DA_SET_PS'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_SET_PS' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_set_ps(proj) 2,2
! Initializes a polar-stereographic map projection from the partially
! filled proj structure. This routine computes the radius to the
! southwest corner and computes the i/j location of the pole for use
! in llij_ps and ijll_ps.
implicit none
type(proj_info), intent(inout) :: proj
real :: ala1
real :: alo1
real :: reflon
real :: scale_top
! To define the cone factor for polar stereographic projection
proj%cone = 1.0
reflon = proj%stdlon + 90.0
! Compute numerator term of map scale factor
scale_top = 1.0 + proj%hemi * Sin(proj%truelat1 * rad_per_deg)
! Compute radius to lower-left (SW) corner
ala1 = proj%lat1 * rad_per_deg
proj%rsw = proj%rebydx*COS(ala1)*scale_top/(1.0+proj%hemi*Sin(ala1))
! Find the pole point
alo1 = (proj%lon1 - reflon) * rad_per_deg
proj%polei = proj%knowni - proj%rsw * COS(alo1)
proj%polej = proj%knownj - proj%hemi * proj%rsw * Sin(alo1)
if (print_detail_map) then
write(unit=stdout,fmt='(A,2F10.1)') 'Computed (I,J) of pole point: ',proj%polei,proj%polej
end if
end subroutine da_set_ps
<A NAME='DA_MAP_INIT'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_MAP_INIT' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_map_init(proj) 2,2
!-----------------------------------------------------------------------
! Purpose: Initializes the map projection structure to missing values
!-----------------------------------------------------------------------
implicit none
type(proj_info), intent(inout) :: proj
proj%lat1 = -999.9
proj%lon1 = -999.9
proj%dx = -999.9
proj%stdlon = -999.9
proj%truelat1 = -999.9
proj%truelat2 = -999.9
proj%hemi = 0.0
proj%cone = -999.9
proj%polei = -999.9
proj%polej = -999.9
proj%rsw = -999.9
proj%knowni = -999.9
proj%knownj = -999.9
proj%latinc = -999.9
proj%loninc = -999.9
proj%init = .false.
end subroutine da_map_init
<A NAME='DA_MAP_SET'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_MAP_SET' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_map_set(proj_code,lat1,lon1,knowni,knownj,dx,stdlon,truelat1,truelat2,latinc,loninc,proj) 3,16
! Given a partially filled proj_info structure, this routine computes
! polei, polej, rsw, and cone (if LC projection) to complete the
! structure. This allows us to eliminate redundant calculations when
! calling the coordinate conversion routines multiple times for the
! same map.
! This will generally be the first routine called when a user wants
! to be able to use the coordinate conversion routines, and it
! will call the appropriate subroutines based on the
! proj%code which indicates which projection type this is.
implicit none
integer, intent(in) :: proj_code
real, intent(in) :: lat1
real, intent(in) :: lon1
real, intent(in) :: dx
real, intent(in) :: stdlon
real, intent(in) :: truelat1
real, intent(in) :: truelat2
real, intent(in) :: knowni , knownj
real, intent(in) :: latinc, loninc
type(proj_info), intent(out) :: proj
! First, check for validity of mandatory variables in proj
if (ABS(lat1) > 90.0) then
write(stdout,*) "da_map_set.inc"// &
"Latitude of origin corner required as follows: -90N <= lat1 < = 90N"
stop
end if
if (ABS(lon1) > 180.0) then
write(stdout,*) "da_map_set.inc"// &
"Longitude of origin required as follows: -180E <= lon1 <= 180W"
stop
end if
if ((dx .LE. 0.0).AND.(proj_code .NE. PROJ_LATLON)) then
write(stdout,*) "da_map_set.inc"// &
"Require grid spacing (dx) in meters be positive!"
stop
end if
if ((ABS(stdlon) > 180.0).AND.(proj_code .NE. PROJ_MERC)) then
write(stdout,*) "da_map_set.inc"// &
"Need orientation longitude (stdlon) as: -180E <= lon1 <= 180W"
stop
end if
if (proj%code .NE. PROJ_LATLON .and. ABS(truelat1)>90.0) then
write(stdout,*) "da_map_set.inc"// &
"Set true latitude 1 for all projections!"
stop
end if
call da_map_init(proj)
proj%code = proj_code
proj%lat1 = lat1
proj%lon1 = lon1
proj%knowni = knowni
proj%knownj = knownj
proj%dx = dx
proj%stdlon = stdlon
proj%truelat1 = truelat1
proj%truelat2 = truelat2
proj%latinc = latinc
proj%loninc = loninc
if (proj%code .NE. PROJ_LATLON) then
proj%dx = dx
if (truelat1 < 0.0) then
proj%hemi = -1.0
else
proj%hemi = 1.0
end if
proj%rebydx = earth_radius_m / dx
end if
pick_proj: select case(proj%code)
case(PROJ_PS)
if (print_detail_map) then
write(unit=stdout,fmt='(A)') 'Setting up POLAR STEREOGRAPHIC map...'
end if
call da_set_ps(proj)
case(PROJ_LC)
if (print_detail_map) then
write(unit=stdout,fmt='(A)') 'Setting up LAMBERT CONFORMAL map...'
end if
if (ABS(proj%truelat2) > 90.0) then
if (print_detail_map) then
write(unit=stdout,fmt='(A)') 'Second true latitude not set, assuming a tangent'
write(unit=stdout,fmt='(A,F10.3)') 'projection at truelat1: ', proj%truelat1
end if
proj%truelat2=proj%truelat1
end if
call da_set_lc(proj)
case (PROJ_MERC)
if (print_detail_map) then
write(unit=stdout,fmt='(A)') 'Setting up MERCATOR map...'
end if
call da_set_merc(proj)
case (PROJ_LATLON)
if (print_detail_map) then
write(unit=stdout,fmt='(A)') 'Setting up CYLinDRICAL EQUIDISTANT LATLON map...'
end if
! Convert lon1 to 0->360 notation
if (proj%lon1 < 0.0) proj%lon1 = proj%lon1 + 360.0
proj%cone = 1.0
case default
write(unit=message(1),fmt='(A,I2)') 'Unknown projection code: ', proj%code
write(stdout,*) "da_map_set.inc"//message(1:1)
stop
end select pick_proj
proj%init = .true.
end subroutine da_map_set
<A NAME='DA_SET_MERC'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_SET_MERC' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_set_merc(proj) 2,2
!--------------------------------------------------------------------------
! Purpose: Sets up the remaining basic elements for the mercator projection
!--------------------------------------------------------------------------
implicit none
type(proj_info), intent(inout) :: proj
real :: clain
! Preliminary variables
clain = COS(rad_per_deg*proj%truelat1)
proj%dlon = proj%dx / (earth_radius_m * clain)
! Compute distance from equator to origin, and store in the
! proj%rsw tag.
proj%rsw = 0.0
if (proj%lat1 .NE. 0.0) then
proj%rsw = (alog(tan(0.5*((proj%lat1+90.)*rad_per_deg))))/proj%dlon
end if
end subroutine da_set_merc
<A NAME='DA_LC_CONE'><A href='../../html_code/verif_obs/da_verif_tools.f90.html#DA_LC_CONE' TARGET='top_target'><IMG SRC="../../gif/bar_red.gif" border=0></A>
subroutine da_lc_cone(truelat1, truelat2, cone) 2,2
!------------------------------------------------------------------------
! Purpose: compute the cone factor of a Lambert Conformal projection
!------------------------------------------------------------------------
implicit none
real, intent(in) :: truelat1 ! (-90 -> 90 degrees N)
real, intent(in) :: truelat2 ! " " " " "
real, intent(out) :: cone
! First, see if this is a secant or tangent projection. For tangent
! projections, truelat1 = truelat2 and the cone is tangent to the
! Earth's surface at this latitude. For secant projections, the cone
! intersects the Earth's surface at each of the distinctly different
! latitudes
if (abs(truelat1-truelat2) > 0.1) then
cone = alog10(cos(truelat1*rad_per_deg)) - &
alog10(cos(truelat2*rad_per_deg))
cone = cone /(alog10(tan((45.0 - abs(truelat1)/2.0) * rad_per_deg)) - &
alog10(tan((45.0 - abs(truelat2)/2.0) * rad_per_deg)))
else
cone = sin(abs(truelat1)*rad_per_deg)
end if
end subroutine da_lc_cone
end module da_verif_tools