module da_verif_tools
   
   
   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

subroutine da_llxy_latlon(lat, lon, proj, x, y)

   !----------------------------------------------------------------------- 
   ! 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


subroutine da_llxy_lc(lat, lon, proj, x, y)

   !-----------------------------------------------------------------------
   ! 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


subroutine da_llxy_merc(lat, lon, proj, x, y)

   !-----------------------------------------------------------------------
   ! 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


subroutine da_llxy_ps(lat,lon,proj,x,y)

   !-----------------------------------------------------------------------
   ! 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


subroutine da_llxy_wrf(proj, lat, lon, x, y)

   !-----------------------------------------------------------------------
   ! 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


subroutine da_xyll(proj, xx, yy, lat, lon)

   !-----------------------------------------------------------------------
   ! 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


subroutine da_xyll_latlon(x, y, proj, lat, lon)

   !-----------------------------------------------------------------------
   ! 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


subroutine da_xyll_lc(x, y, proj, lat, lon)

   ! 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


subroutine da_xyll_merc(x, y, proj, lat, lon)

   !-----------------------------------------------------------------------
   ! 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


subroutine da_xyll_ps(x, y, proj, lat, lon)

   ! 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


subroutine da_set_lc(proj)

   !-----------------------------------------------------------------------
   ! 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                             


subroutine da_set_ps(proj)

   ! 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


subroutine da_map_init(proj)

   !-----------------------------------------------------------------------
   ! 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


subroutine da_map_set(proj_code,lat1,lon1,knowni,knownj,dx,stdlon,truelat1,truelat2,latinc,loninc,proj)
   ! 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


subroutine da_set_merc(proj)
  
   !--------------------------------------------------------------------------
   ! 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


subroutine da_lc_cone(truelat1, truelat2, cone)

   !------------------------------------------------------------------------
   ! 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