subroutine weightsa (proj,inter,stlon,stlat,dlon,dlat,end_index,end_index_n,lat,lon,wiss,index,lcross,wagdot,lmore)

  implicit none

  integer, dimension(4) :: end_index, end_index_n

  real, dimension(end_index(1),end_index(2)) :: lat, lon
  real, dimension(end_index_n(1),end_index_n(2),4) :: wiss
  real, dimension(end_index_n(1),end_index_n(2),1) :: wagdot
  real :: dlon, dlat, stlat, stlon
  integer :: inter
  integer, dimension(end_index_n(1),end_index_n(2),4,2) :: index

  integer :: i, j, k, l, m
  real :: swiss
  double precision :: arg, arglat, arglon
!  double precision :: arga, lata, latb
  double precision :: rearth, pi, fr, rf, xlon, xlat
  real :: ill,jll
  real :: iall,jall,ibll,jbll
!  double precision :: argab

  logical :: lcross
  logical :: lmore

  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
     LOGICAL          :: init     ! Flag to indicate if this struct is
                                  ! ready for use
  END TYPE proj_info

  TYPE(proj_info) :: proj

  rearth=6370.
  pi=3.14159265358
  fr=pi/180.
  rf=180./pi

  if (lmore) write(*,'(/,8x, "Subroutine weightsa: Interpolation method: ",I7,/)') inter

  select case (inter)

  case (0)

  if (lcross) then
  do i=1,end_index_n(1)
    xlon=dreal(stlon+(i-1)*dlon)
  do j=1,end_index_n(2)
    xlat=dreal(stlat+(j-1)*dlat)

            call latlon_to_ij(proj, sngl(xlat), sngl(xlon), jll, ill)
!            print*, 'Titi_0', ill, jll

            l=aint(ill)
            k=aint(jll)

            if (l.lt.1) then
               write(*,'(/,8x, "Error: start_lat must be greater",/)')
               stop
            endif
            if (l.ge.end_index(1)) then
               write(*,'(/,8x, "Error: end_lat must be lesser",/)')
               stop
            endif
            if (k.lt.1) then
               write(*,'(/,8x, "Error: start_lon must be greater",/)')
               stop
            endif
            if (k.ge.end_index(2)) then
               write(*,'(/,8x, "Error: end_lon must be lesser",/)')
               stop
            endif


            index(i,j,1,1)=l
            index(i,j,1,2)=k
            wiss(i,j,1)=ill
            wiss(i,j,2)=jll

!            if (lmore) print* , 'Hola case (0)', index(i,j,1,1), index(i,j,1,2), wiss(i,j,1), wiss(i,j,2)
!            stop 0 
  enddo
  enddo
!  stop
  else
  do i=1,end_index_n(1)
   xlon=dreal(stlon+(i-1)*dlon)
  do j=1,end_index_n(2)
    xlat=dreal(stlat+(j-1)*dlat)

            call latlon_to_ij(proj, sngl(xlat), sngl(xlon), jll, ill)
            l=aint(ill)
            k=aint(jll)
!            print*, 'Titi_2', ill, jll
!            stop 122
!            print*, xlat, xlon, l,'<',end_index(1), k, '<', end_index(2)
            if (l.lt.1) then
                write(*,'(/,8x, "Error: start_lat must be greater",/)')
                stop
            endif
            if (l.ge.end_index(1)) then
                write(*,'(/,8x, "Error: end_lat must be lesser",/)')
                stop
            endif
            if (k.lt.1) then
                write(*,'(/,8x, "Error: start_lon must be greater",/)')
                stop
            endif
            if (k.ge.end_index(2)) then
                write(*,'(/,8x, "Error: end_lon must be lesser",/)')
                stop
            endif

            index(i,j,1,1)=l
            index(i,j,1,2)=k
            wiss(i,j,1)=ill
            wiss(i,j,2)=jll

!            if (lmore) print* , 'Hola case (0)', index(i,j,1,1), index(i,j,1,2), wiss(i,j,1), wiss(i,j,2)
!            stop 0 

! It computes the angles between the grids for the point of lon/lat
            call latlon_to_ij(proj,sngl(xlat),sngl(xlon),jall,iall)
            call latlon_to_ij(proj,sngl(xlat+0.01),sngl(xlon),jbll,ibll)
            wagdot(i,j,1)=sngl(rf*datan(dreal((jbll-jall)/(ibll-iall))))
!            print*, 'Arga', i, j, jbll,'-',jall,ibll,'-',iall,wagdot(i,j,1)

  enddo
  enddo
  endif

!    stop

  case (3)

  if (lcross) then
  do i=1,end_index_n(1)
    xlon=dreal(stlon+(i-1)*dlon)
  do j=1,end_index_n(2)
    xlat=dreal(stlat+(j-1)*dlat)

            call latlon_to_ij(proj, sngl(xlat), sngl(xlon), jll, ill)
!            print*, 'Titi_0', ill, jll

!            do m=1,4
!            call latlon_to_ij(proj, lat(index(i,j,m,1),index(i,j,m,2)),&
!                 lon(index(i,j,m,1),index(i,j,m,2)), jll, ill)
!            print*, 'Titi_1', ill, jll
!            xlat=41.38
!            xlon=2.12
!            call latlon_to_ij(proj, sngl(xlat), sngl(xlon), jll, ill)
!            print*, 'Titi_1', ill, jll
!            print*, 'Titi_1', aint(ill), aint(jll)
!            enddo
!            stop 121


            l=aint(ill)
            k=aint(jll)

            if (l.lt.1) then
               write(*,'(/,8x, "Error: start_lat must be greater",/)')
               stop
            endif
            if (l.ge.end_index(1)) then
               write(*,'(/,8x, "Error: end_lat must be lesser",/)')
               stop
            endif
            if (k.lt.1) then
               write(*,'(/,8x, "Error: start_lon must be greater",/)')
               stop
            endif
            if (k.ge.end_index(2)) then
               write(*,'(/,8x, "Error: end_lon must be lesser",/)')
               stop
            endif

            index(i,j,1,1)=l
            index(i,j,1,2)=k

            index(i,j,2,1)=l
            index(i,j,2,2)=k+1

            index(i,j,3,1)=l+1
            index(i,j,3,2)=k

            index(i,j,4,1)=l+1
            index(i,j,4,2)=k+1

            swiss=0.

            do m=1,4
!
!              Spheric trigonometry
!
               arg=0.
               arg=(((sin(fr*dreal(lat(index(i,j,m,1),index(i,j,m,2)))))*(sin(fr*xlat)))+&
                    ((cos(fr*dreal(lat(index(i,j,m,1),index(i,j,m,2)))))*(cos(fr*xlat))*&
                     (cos(fr*((dreal(lon(index(i,j,m,1),index(i,j,m,2))))-xlon)))))
!               print*, lat(index(i,j,m,1),index(i,j,m,2)), xlat,&
!                       lon(index(i,j,m,1),index(i,j,m,2)), xlon, arg
!
!!               Plane
!
!               if (arg.ge.1) then
!!                  arg=pi*rearth*sqrt((lat(index(i,j,m,1),index(i,j,m,2))-xlat)**2+ &
!!                      (lon(index(i,j,m,1),index(i,j,m,2))-xlon)**2*  &
!!                      (cos((lat(index(i,j,m,1),index(i,j,m,2))+xlat)*pi/360.))**2)/180.
!!                  if ((arg).eq.0) arg=0.000001
!               else
                  arg=rearth*(acos(arg))
!               endif
               wiss(i,j,m)=sngl(1./arg)
               swiss=swiss+wiss(i,j,m)
            enddo

            do m=1,4
               wiss(i,j,m)=wiss(i,j,m)/swiss
!               print*, 'A', xlat, xlon, wiss(i,j,m)
            enddo

!            stop 121
  enddo
  enddo

  else
  do i=1,end_index_n(1)
    xlon=dreal(stlon+(i-1)*dlon)
  do j=1,end_index_n(2)
    xlat=dreal(stlat+(j-1)*dlat)

            call latlon_to_ij(proj, sngl(xlat), sngl(xlon), jll, ill)
            l=aint(ill)
            k=aint(jll)
!            print*, 'Titi_2', ill, jll
!            stop 122
!            print*, xlat, xlon, l,'<',end_index(1), k, '<', end_index(2)
            if (l.lt.1) then
               write(*,'(/,8x, "Error: start_lat must be greater",/)')
               stop
            endif
            if (l.ge.end_index(1)) then
                write(*,'(/,8x, "Error: end_lat must be lesser",/)')
                stop
            endif
            if (k.lt.1) then
                write(*,'(/,8x, "Error: start_lon must be greater",/)')
                stop
            endif
            if (k.ge.end_index(2)) then
                write(*,'(/,8x, "Error: end_lon must be lesser",/)')
                stop
            endif

            index(i,j,1,1)=l
            index(i,j,1,2)=k

            index(i,j,2,1)=l
            index(i,j,2,2)=k+1

            index(i,j,3,1)=l+1
            index(i,j,3,2)=k

            index(i,j,4,1)=l+1
            index(i,j,4,2)=k+1

            call latlon_to_ij(proj, 41.38, 2.12, jll, ill)
!            print*, 'Titi_1', ill, jll
!            print*, 'Titi_1', aint(ill), aint(jll)
!            stop 122

            swiss=0.

            do m=1,4
                arg=(((sin(fr*dreal(lat(index(i,j,m,1),index(i,j,m,2)))))*(sin(fr*xlat)))+&
                     ((cos(fr*dreal(lat(index(i,j,m,1),index(i,j,m,2)))))*(cos(fr*xlat))*&
                      (cos(fr*((dreal(lon(index(i,j,m,1),index(i,j,m,2))))-xlon)))))
!               print*, lat(index(i,j,m,1),index(i,j,m,2)), xlat,&
!                       lon(index(i,j,m,1),index(i,j,m,2)), xlon, arg
!               if (arg.ge.1) then
!!                  arg=pi*rearth*sqrt((lat(index(i,j,m,1),index(i,j,m,2))-xlat)**2+ &
!!                      (lon(index(i,j,m,1),index(i,j,m,2))-xlon)**2*  &
!!                      (cos((lat(index(i,j,m,1),index(i,j,m,2))+xlat)*pi/360.))**2)/180.
!!                  if ((arg).eq.0) arg=0.000001
!               else
                  arg=rearth*(acos(arg))
!               endif
                wiss(i,j,m)=sngl(1./arg)
                swiss=swiss+wiss(i,j,m)
! It computes the angles between the grids for the 4 points
!               call latlon_to_ij(proj, lat(index(i,j,m,1),index(i,j,m,2)),&
!                                 lon(index(i,j,m,1),index(i,j,m,2)), iall, jall)
!               call latlon_to_ij(proj, ((lat(index(i,j,m,1),index(i,j,m,2)))+0.01),&
!                                 lon(index(i,j,m,1),index(i,j,m,2)), ibll, jbll)
!               wagdot(i,j,m)=sngl(rf*datan(dreal((ibll-iall)/(jbll-jall))))

            enddo

            do m=1,4
               wiss(i,j,m)=wiss(i,j,m)/swiss
!               print*, 'A', xlat, xlon, wiss(i,j,m)
!               print*, 'Arga', wagdot(i,j,m)
            enddo

! It computes the angles between the grids for the point of lon/lat
            call latlon_to_ij(proj,sngl(xlat),sngl(xlon),jall,iall)
            call latlon_to_ij(proj,sngl(xlat+0.01),sngl(xlon),jbll,ibll)
            wagdot(i,j,1)=sngl(rf*datan(dreal((jbll-jall)/(ibll-iall))))
!            print*, 'Arga', i, j, jbll,'-',jall,ibll,'-',iall,wagdot(i,j,1)

! It computes the angles between the grids with spheric triangles

!            lata=(dreal(lat(index(i,j,1,1),index(i,j,1,2))))-&
!                ((dreal(lat(index(i,j,3,1),index(i,j,3,2))))-(xlat))
!            latb=xlat

!             arga=(((sin(fr*lata))*(sin(fr*latb)))+((cos(fr*lata))*(cos(fr*latb))*&
!                    (cos(fr*((dreal(lon(index(i,j,1,1),index(i,j,1,2))))-&
!                     dreal(lon(index(i,j,3,1),index(i,j,3,2))))))))
!             arga=rf*asin(((cos(fr*latb))*&
!                           (sin(fr*((dreal(lon(index(i,j,1,1),index(i,j,1,2))))-&
!                            dreal(lon(index(i,j,3,1),index(i,j,3,2)))))))/&
!                           (sin(acos(arga))))


!            wagdot=sngl(arga)

!            call latlon_to_ij(proj, sngl(xlat+0.001), sngl(xlon), all, bll)
!            argab=rf*datan(dreal((all-ill)/(bll-jll)))


!            print*, 'Arga', arga, wagdot(i,j)
!            print*, 'Arga', arga, argab

  enddo
  enddo

  endif

!    stop


  case (2)

  if (lcross) then
  do i=1,end_index_n(1)
    xlon=dreal(stlon+(i-1)*dlon)
  do j=1,end_index_n(2)
    xlat=dreal(stlat+(j-1)*dlat)

            call latlon_to_ij(proj, sngl(xlat), sngl(xlon), jll, ill)
!            print*, 'Titi_0', ill, jll

!            do m=1,4
!            call latlon_to_ij(proj, lat(index(i,j,m,1),index(i,j,m,2)),&
!                 lon(index(i,j,m,1),index(i,j,m,2)), jll, ill)
!            print*, 'Titi_1', ill, jll
!            xlat=41.38
!            xlon=2.12
!            call latlon_to_ij(proj, sngl(xlat), sngl(xlon), jll, ill)
!            print*, 'Titi_1', ill, jll
!            print*, 'Titi_1', aint(ill), aint(jll)
!            enddo
!            stop 121


            l=aint(ill)
            k=aint(jll)

            if (l.lt.1) then
               write(*,'(/,8x, "Error: start_lat must be greater",/)')
               stop
            endif
            if (l.ge.end_index(1)) then
               write(*,'(/,8x, "Error: end_lat must be lesser",/)')
               stop
            endif
            if (k.lt.1) then
               write(*,'(/,8x, "Error: start_lon must be greater",/)')
               stop
            endif
            if (k.ge.end_index(2)) then
               write(*,'(/,8x, "Error: end_lon must be lesser",/)')
               stop
            endif

            index(i,j,1,1)=l
            index(i,j,1,2)=k

            index(i,j,2,1)=l
            index(i,j,2,2)=k+1


            index(i,j,3,1)=l+1
            index(i,j,3,2)=k

            index(i,j,4,1)=l+1
            index(i,j,4,2)=k+1

            swiss=0.

            do m=1,4
!
!              Spheric trigonometry
!
               arg=0.
               arg=(((sin(fr*dreal(lat(index(i,j,m,1),index(i,j,m,2)))))*(sin(fr*xlat)))+&
                    ((cos(fr*dreal(lat(index(i,j,m,1),index(i,j,m,2)))))*(cos(fr*xlat))*&
                     (cos(fr*((dreal(lon(index(i,j,m,1),index(i,j,m,2))))-xlon)))))
!               print*, lat(index(i,j,m,1),index(i,j,m,2)), xlat,&
!                       lon(index(i,j,m,1),index(i,j,m,2)), xlon, arg
!
!!               Plane
!
!               if (arg.ge.1) then
!!                  arg=pi*rearth*sqrt((lat(index(i,j,m,1),index(i,j,m,2))-xlat)**2+ &
!!                      (lon(index(i,j,m,1),index(i,j,m,2))-xlon)**2*  &
!!                      (cos((lat(index(i,j,m,1),index(i,j,m,2))+xlat)*pi/360.))**2)/180.
!!                  if ((arg).eq.0) arg=0.000001
!               else
                  arg=rearth*(acos(arg))
!               endif
               wiss(i,j,m)=sngl(1./(arg*arg))
               swiss=swiss+wiss(i,j,m)
            enddo

            do m=1,4
               wiss(i,j,m)=wiss(i,j,m)/swiss
!               print*, 'A', xlat, xlon, wiss(i,j,m)
            enddo

!            stop 121
  enddo
  enddo

  else
  do i=1,end_index_n(1)
    xlon=dreal(stlon+(i-1)*dlon)
  do j=1,end_index_n(2)

    xlat=dreal(stlat+(j-1)*dlat)

            call latlon_to_ij(proj, sngl(xlat), sngl(xlon), jll, ill)
            l=aint(ill)
            k=aint(jll)
!            print*, 'Titi_2', ill, jll
!            stop 122
!            print*, xlat, xlon, l,'<',end_index(1), k, '<', end_index(2)
            if (l.lt.1) then
               write(*,'(/,8x, "Error: start_lat must be greater",/)')
               stop
            endif
            if (l.ge.end_index(1)) then
                write(*,'(/,8x, "Error: end_lat must be lesser",/)')
                stop
            endif
            if (k.lt.1) then
                write(*,'(/,8x, "Error: start_lon must be greater",/)')
                stop
            endif
            if (k.ge.end_index(2)) then
                write(*,'(/,8x, "Error: end_lon must be lesser",/)')
                stop
            endif

            index(i,j,1,1)=l
            index(i,j,1,2)=k

            index(i,j,2,1)=l
            index(i,j,2,2)=k+1

            index(i,j,3,1)=l+1
            index(i,j,3,2)=k

            index(i,j,4,1)=l+1
            index(i,j,4,2)=k+1

            call latlon_to_ij(proj, 41.38, 2.12, jll, ill)
!            print*, 'Titi_1', ill, jll
!            print*, 'Titi_1', aint(ill), aint(jll)
!            stop 122

            swiss=0.

            do m=1,4
                arg=(((sin(fr*dreal(lat(index(i,j,m,1),index(i,j,m,2)))))*(sin(fr*xlat)))+&
                     ((cos(fr*dreal(lat(index(i,j,m,1),index(i,j,m,2)))))*(cos(fr*xlat))*&
                     (cos(fr*((dreal(lon(index(i,j,m,1),index(i,j,m,2))))-xlon)))))
!               print*, lat(index(i,j,m,1),index(i,j,m,2)), xlat,&
!                       lon(index(i,j,m,1),index(i,j,m,2)), xlon, arg
!               if (arg.ge.1) then
!!                  arg=pi*rearth*sqrt((lat(index(i,j,m,1),index(i,j,m,2))-xlat)**2+ &
!!                      (lon(index(i,j,m,1),index(i,j,m,2))-xlon)**2*  &
!!                      (cos((lat(index(i,j,m,1),index(i,j,m,2))+xlat)*pi/360.))**2)/180.
!!                  if ((arg).eq.0) arg=0.000001
!               else
                  arg=rearth*(acos(arg))
!               endif
                wiss(i,j,m)=sngl(1./(arg*arg))
                swiss=swiss+wiss(i,j,m)
! It computes the angles between the grids for the 4 points
!               call latlon_to_ij(proj, lat(index(i,j,m,1),index(i,j,m,2)),&
!                                 lon(index(i,j,m,1),index(i,j,m,2)), iall, jall)
!               call latlon_to_ij(proj, ((lat(index(i,j,m,1),index(i,j,m,2)))+0.01),&
!                                 lon(index(i,j,m,1),index(i,j,m,2)), ibll, jbll)
!               wagdot(i,j,m)=sngl(rf*datan(dreal((ibll-iall)/(jbll-jall))))

            enddo

            do m=1,4
               wiss(i,j,m)=wiss(i,j,m)/swiss
!               print*, 'A', xlat, xlon, wiss(i,j,m)
!               print*, 'Arga', wagdot(i,j,m)
            enddo

! It computes the angles between the grids for the point of lon/lat
            call latlon_to_ij(proj,sngl(xlat),sngl(xlon),jall,iall)
            call latlon_to_ij(proj,sngl(xlat+0.01),sngl(xlon),jbll,ibll)
            wagdot(i,j,1)=sngl(rf*datan(dreal((jbll-jall)/(ibll-iall))))
!            print*, 'Arga', i, j, jbll,'-',jall,ibll,'-',iall,wagdot(i,j,1)

! It computes the angles between the grids with spheric triangles

!            lata=(dreal(lat(index(i,j,1,1),index(i,j,1,2))))-&
!                ((dreal(lat(index(i,j,3,1),index(i,j,3,2))))-(xlat))
!            latb=xlat

!             arga=(((sin(fr*lata))*(sin(fr*latb)))+((cos(fr*lata))*(cos(fr*latb))*&
!                    (cos(fr*((dreal(lon(index(i,j,1,1),index(i,j,1,2))))-&
!                     dreal(lon(index(i,j,3,1),index(i,j,3,2))))))))
!             arga=rf*asin(((cos(fr*latb))*&
!                           (sin(fr*((dreal(lon(index(i,j,1,1),index(i,j,1,2))))-&
!                            dreal(lon(index(i,j,3,1),index(i,j,3,2)))))))/&
!                           (sin(acos(arga))))


!            wagdot=sngl(arga)

!            call latlon_to_ij(proj, sngl(xlat+0.001), sngl(xlon), all, bll)
!            argab=rf*datan(dreal((all-ill)/(bll-jll)))


!            print*, 'Arga', arga, wagdot(i,j)
!            print*, 'Arga', arga, argab

  enddo
  enddo

  endif

!    stop

  case (1)

  if (lcross) then
  do i=1,end_index_n(1)
    xlon=dreal(stlon+(i-1)*dlon)
  do j=1,end_index_n(2)
    xlat=dreal(stlat+(j-1)*dlat)

            call latlon_to_ij(proj, sngl(xlat), sngl(xlon), jll, ill)
!            print*, 'Titi_0', ill, jll

            l=aint(ill)
            k=aint(jll)

            if (l.lt.1) then
               write(*,'(/,8x, "Error: start_lat must be greater",/)')
               stop
            endif
            if (l.ge.end_index(1)) then
               write(*,'(/,8x, "Error: end_lat must be lesser",/)')
               stop
            endif
            if (k.lt.1) then
               write(*,'(/,8x, "Error: start_lon must be greater",/)')
               stop
            endif
            if (k.ge.end_index(2)) then
               write(*,'(/,8x, "Error: end_lon must be lesser",/)')
               stop
            endif

!
!                    (l+1,k)     (l+1,k+1)
!               -        -------------
!      l+1-ill  |        |           |
!               |        |           |
!               -        |  *        |
!        ill-l  |        | (ill,jll) |
!               |        |           |
!               -        -------------
!                      (l,k)       (l,k+1)
!
!                        |--|--------|
!                      jll-k   k+1-jll 


            index(i,j,1,1)=l
            index(i,j,1,2)=k
            wiss(i,j,1)=(l+1-ill)*(k+1-jll)

            index(i,j,2,1)=l
            index(i,j,2,2)=k+1
            wiss(i,j,2)=(l+1-ill)*(jll-k)

            index(i,j,3,1)=l+1
            index(i,j,3,2)=k
            wiss(i,j,3)=(ill-l)*(k+1-jll)

            index(i,j,4,1)=l+1
            index(i,j,4,2)=k+1
            wiss(i,j,4)=(ill-l)*(jll-k)

!            print* , 'Hola', i, j, wiss(i,j,1), wiss(i,j,2), wiss(i,j,3), wiss(i,j,4)
!            print* , 'Hola2', wiss(i,j,1) + wiss(i,j,2) + wiss(i,j,3) + wiss(i,j,4)
!            stop 121
  enddo
  enddo
!  stop
  else
  do i=1,end_index_n(1)
    xlon=dreal(stlon+(i-1)*dlon)
  do j=1,end_index_n(2)
    xlat=dreal(stlat+(j-1)*dlat)

            call latlon_to_ij(proj, sngl(xlat), sngl(xlon), jll, ill)
            l=aint(ill)
            k=aint(jll)
!            print*, 'Titi_2', ill, jll
!            stop 122
!            print*, xlat, xlon, l,'<',end_index(1), k, '<', end_index(2)
            if (l.lt.1) then
               write(*,'(/,8x, "Error: start_lat must be greater",/)')
               stop
            endif
            if (l.ge.end_index(1)) then
                write(*,'(/,8x, "Error: end_lat must be lesser",/)')
                stop
            endif
            if (k.lt.1) then
                write(*,'(/,8x, "Error: start_lon must be greater",/)')
                stop
            endif
            if (k.ge.end_index(2)) then
                write(*,'(/,8x, "Error: end_lon must be lesser",/)')
                stop
            endif

            index(i,j,1,1)=l
            index(i,j,1,2)=k
            wiss(i,j,1)=(l+1-ill)*(k+1-jll)

            index(i,j,2,1)=l
            index(i,j,2,2)=k+1
            wiss(i,j,2)=(l+1-ill)*(jll-k)

            index(i,j,3,1)=l+1
            index(i,j,3,2)=k
            wiss(i,j,3)=(ill-l)*(k+1-jll)

            index(i,j,4,1)=l+1
            index(i,j,4,2)=k+1
            wiss(i,j,4)=(ill-l)*(jll-k)

! It computes the angles between the grids for the point of lon/lat
            call latlon_to_ij(proj,sngl(xlat),sngl(xlon),jall,iall)
            call latlon_to_ij(proj,sngl(xlat+0.01),sngl(xlon),jbll,ibll)
            wagdot(i,j,1)=sngl(rf*datan(dreal((jbll-jall)/(ibll-iall))))
!            print*, 'Arga', i, j, jbll,'-',jall,ibll,'-',iall,wagdot(i,j,1)

  enddo
  enddo
  endif

!    stop

  case default
    write(*,'(/,8x, "Subroutine weightsa: incorrect interpolation method",/)')
    stop
  end select

end subroutine weightsa

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE latlon_to_ij(proj, lat, lon, i, j)
! Converts input lat/lon values to the cartesian (i,j) value
! for the given projection.
! Source: /LITTLE_R/src/map_utils_module.F

IMPLICIT NONE
! TYPE(proj_info), INTENT(IN)          :: proj
REAL, INTENT(IN)                     :: lat
REAL, INTENT(IN)                     :: lon
REAL, INTENT(OUT)                    :: i
REAL, INTENT(OUT)                    :: j

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
    LOGICAL          :: init     ! Flag to indicate if this struct is
                                 ! ready for use
END TYPE proj_info

TYPE(proj_info), INTENT(IN)          :: proj

IF (.NOT.proj%init) THEN
   PRINT '(A)', 'You have not called map_set for this projection!'
   STOP 'LATLON_TO_IJ'
ENDIF

SELECT CASE(proj%code)

   CASE(5)
       CALL llij_ps(lat,lon,proj,i,j)

   CASE(3)
       CALL llij_lc(lat,lon,proj,i,j)

   CASE(1)
       CALL llij_merc(lat, lon, proj, i, j)

   CASE DEFAULT
       PRINT '(A,I2)', 'Unrecognized map projection code: ', proj%code
       STOP 'LATLON_TO_IJ'

 END SELECT
 RETURN
 END SUBROUTINE latlon_to_ij

SUBROUTINE llij_ps(lat,lon,proj,i,j)
! 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 i/j indices which
! if within the domain range from 1->nx and 1->ny, respectively.

 IMPLICIT NONE

 REAL, PARAMETER    :: pi = 3.14159265358
 REAL, PARAMETER    :: rad_per_deg = pi / 180.


! Delcare input arguments
 REAL, INTENT(IN)               :: lat
 REAL, INTENT(IN)               :: lon

 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
    LOGICAL          :: init     ! Flag to indicate if this struct is
                                 ! ready for use
 END TYPE proj_info

 TYPE(proj_info),INTENT(IN)     :: proj

! Declare output arguments
 REAL, INTENT(OUT)              :: i !(x-index)
 REAL, INTENT(OUT)              :: j !(y-index)

! Declare local variables
    
 REAL                           :: reflon
 REAL                           :: scale_top
 REAL                           :: ala
 REAL                           :: alo
 REAL                           :: rm

! BEGIN CODE

  
 reflon = proj%stdlon + 90.
   
! Compute numerator term of map scale factor

 scale_top = 1. + 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. + proj%hemi *SIN(ala))
 alo = (lon - reflon) * rad_per_deg
 i = proj%polei + rm * COS(alo)
 j = proj%polej + proj%hemi * rm * SIN(alo)

 RETURN

END SUBROUTINE llij_ps
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

SUBROUTINE llij_lc( lat, lon, proj, i, j)

! Subroutine to compute the geographical latitude and longitude values
! to the cartesian x/y on a Lambert Conformal projection.

  IMPLICIT NONE

  REAL, PARAMETER    :: pi = 3.14159265358
  REAL, PARAMETER    :: rad_per_deg = pi / 180.

! Input Args
  REAL, INTENT(IN)              :: lat      ! Latitude (-90->90 deg N)
  REAL, INTENT(IN)              :: lon      ! Longitude (-180->180 E)

  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
    LOGICAL          :: init     ! Flag to indicate if this struct is
                                 ! ready for use
  END TYPE proj_info

  TYPE(proj_info),INTENT(IN)      :: proj     ! Projection info structure

! Output Args
  REAL, INTENT(OUT)             :: i        ! Cartesian X coordinate
  REAL, INTENT(OUT)             :: j        ! Cartesian Y coordinate

! Locals
  REAL                          :: arg
  REAL                          :: deltalon
  REAL                          :: tl1r
  REAL                          :: rm
  REAL                          :: ctl1r


! BEGIN CODE

! Compute deltalon between known longitude and standard lon and ensure
! it is not in the cut zone
 deltalon = lon - proj%stdlon
 IF (deltalon .GT. +180.) deltalon = deltalon - 360.
 IF (deltalon .LT. -180.) deltalon = deltalon + 360.

! 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.*proj%hemi-lat)*rad_per_deg/2.) / &
 TAN((90.*proj%hemi-proj%truelat1)*rad_per_deg/2.))**proj%cone

 arg = proj%cone*(deltalon*rad_per_deg)
 i = proj%polei + proj%hemi * rm * SIN(arg)
 j = 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 .EQ. -1.) THEN
   i = 2. - i
   j = 2. - j
 ENDIF
 RETURN
END SUBROUTINE llij_lc

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
SUBROUTINE llij_merc(lat, lon, proj, i, j)

! Compute i/j coordinate from lat lon for mercator projection
 
  IMPLICIT NONE

  REAL, PARAMETER    :: pi = 3.14159265358
  REAL, PARAMETER    :: rad_per_deg = pi / 180.
  REAL, PARAMETER    :: deg_per_rad = 180./pi

  REAL, INTENT(IN)              :: lat
  REAL, INTENT(IN)              :: lon

  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
    LOGICAL          :: init     ! Flag to indicate if this struct is
                                 ! ready for use
  END TYPE proj_info

  TYPE(proj_info),INTENT(IN)    :: proj

  REAL,INTENT(OUT)              :: i
  REAL,INTENT(OUT)              :: j
  REAL                          :: deltalon

  deltalon = lon - proj%lon1
  IF (deltalon .LT. -180.) deltalon = deltalon + 360.
  IF (deltalon .GT. 180.) deltalon = deltalon - 360.
  i = 1. + (deltalon/(proj%dlon*deg_per_rad))
  j = 1. + (ALOG(TAN(0.5*((lat + 90.) * rad_per_deg)))) / &
         proj%dlon - proj%rsw

  RETURN
END SUBROUTINE llij_merc
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
