subroutine grid_lonlat_input(MY,MX,start_lon,start_lat,interval_lon,interval_lat,inter,lmore)

 implicit none 

 character(len=20) :: gridfile='grid_lonlat.input'
 character(len=100) :: gridline
 character(len=20) :: parametre

 integer :: iunit=12
 integer :: ierr

 real :: start_lon,end_lon,interval_lon,start_lat,end_lat,interval_lat
 integer :: MX,MY
 integer :: inter

 logical :: lmore

 open(iunit, file=gridfile, err=299, form='formatted', status='old',&
      action='read')

       if (lmore) write(*,'("Subroutine grid_lonlat_input: Reading grid_lonlat.input file ...")')
       read(iunit, '(A100)', iostat=ierr, end=200) gridline
       if (lmore) write(*,'(A100)') gridline
       read(iunit, '(A100)', iostat=ierr, end=200) gridline
       if (lmore) write(*,'(A100)') gridline
       read(iunit, '(A100)', iostat=ierr, end=200) gridline
       if (gridline(1:8).ne.'&record1') then
         parametre='&record1'
         goto 297
       endif
       if (lmore) write(*,'(A100)') gridline
       read(iunit, 210, iostat=ierr, end=200) gridline, start_lon
       if (gridline(1:15).ne.' start_lon    =') then
         parametre='start_lon'
         goto 297
       endif
       if (lmore) write(*,210) gridline, start_lon
       read(iunit, 210, iostat=ierr, end=200) gridline, end_lon
       if (gridline(1:15).ne.' end_lon      =') then
         parametre='end_lon'
         goto 297
       endif
       if (lmore) write(*,210) gridline, end_lon
       read(iunit, 210, iostat=ierr, end=200) gridline, interval_lon
       if (gridline(1:15).ne.' interval_lon =') then
        parametre='interval_lon'
        goto 297
       endif
       if (lmore) write(*,210) gridline, interval_lon
       read(iunit, 210, iostat=ierr, end=200) gridline, start_lat
       if (gridline(1:15).ne.' start_lat    =') then
         parametre='start_lat'
         goto 297
       endif
       if (lmore) write(*,210) gridline, start_lat
       read(iunit, 210, iostat=ierr, end=200) gridline, end_lat
       if (gridline(1:15).ne.' end_lat      =') then
         parametre='end_lat'
         goto 297
       endif
       if (lmore) write(*,210) gridline, end_lat
       read(iunit, 210, iostat=ierr, end=200) gridline, interval_lat
       if (gridline(1:15).ne.' interval_lat =') then
         parametre='interval_lat'
         goto 297
       endif
       if (lmore) write(*,210) gridline, interval_lat

       read(iunit, '(A100)', iostat=ierr, end=200) gridline
       if (lmore) write(*,'(A100)') gridline
       read(iunit, '(A100)', iostat=ierr, end=200) gridline
       if (lmore) write(*,'(A100)') gridline
       read(iunit, '(A100)', iostat=ierr, end=200) gridline
       if (lmore) write(*,'(A100)') gridline
       read(iunit, '(A100)', iostat=ierr, end=200) gridline
       if (lmore) write(*,'(A100)') gridline
       read(iunit, '(A100)', iostat=ierr, end=200) gridline
       if (lmore) write(*,'(A100)') gridline
       read(iunit, '(A100)', iostat=ierr, end=200) gridline
       if (lmore) write(*,'(A100)') gridline
       read(iunit, '(A100)', iostat=ierr, end=200) gridline
       if (lmore) write(*,'(A100)') gridline
       read(iunit, '(A100)', iostat=ierr, end=200) gridline
       if (lmore) write(*,'(A100)') gridline
       read(iunit, 211, iostat=ierr, end=200) gridline, inter
       if (gridline(1:15).ne.' inter        =') then
         parametre='inter'
         goto 297
       endif
       if (lmore) write(*,211) gridline, inter
       read(iunit, '(A100)', iostat=ierr, end=200) gridline
       if (lmore) write(*,'(A100)') gridline
       read(iunit, '(A100)', iostat=ierr, end=200) gridline
       if (lmore) write(*,'(A100)') gridline

210 format(A15,F9.3)
211 format(A15,I7)
     goto 200
297 write(*,'(/,8x, "Subroutine grid_lonlat_input: Error reading in", A, 1x, "the parameter ", A)')&
          trim(gridfile), trim(parametre)
    stop 297
200 close(iunit, err=298)

!    stop
! First check of limits
!    write(*,'(F8.2)') amod(end_lon-start_lon,interval_lon)
!!    write(*,'(F8.2)') ((end_lon-start_lon)-int((end_lon-start_lon)/interval_lon)*interval_lon)
!    if ((amod(end_lon-start_lon,interval_lon)).ne.0) then
    if (abs(((end_lon-start_lon)*1000000)-int(((end_lon-start_lon)/interval_lon)*1000000)*interval_lon).ge.1) then
      write(*,'(/,8x, "Subroutine grid_lonlat_input: Error: end_lon - start_lon must be multiple of interval_lon",/)')
      stop
    else
       MY=nint((end_lon-start_lon)/interval_lon)
       MY=MY+1
       if (lmore) write(*,'("Subroutine grid_lonlat_input: Number of grid points MY in longitude:", I5)') MY
    endif
!    if ((amod(end_lat-start_lat,interval_lat)).ne.0) then
    if (abs(((end_lat-start_lat)*1000000)-int(((end_lat-start_lat)/interval_lat)*1000000)*interval_lat).ge.1) then
      write(*,'(/,8x, "Subroutine grid_lonlat_input: Error: end_lat - start_lat must be multiple by interval_lat",/)')
      stop
    else
      MX=nint((end_lat-start_lat)/interval_lat)
      MX=MX+1
      if (lmore) write(*,'("Subroutine grid_lonlat_input: Number of grid points MX in  latitude:", I5)') MX
    endif

!
! Check interpolation method
!
    if ((inter.ne.0).and.(inter.ne.1).and.(inter.ne.2).and.(inter.ne.3)) then
      write(*,'(/,8x, "Subroutine grid_lonlat_input: Error: inter must be 1, 2 or 3",/)')
      stop
    endif

    goto 300
298   write(*,'(/,8x, "Subroutine grid_lonlat_input: Error closing the ", A, 1x, "file",/)') trim(gridfile)
       stop 298
299   write(*,'(/,8x, "Subroutine grid_lonlat_input: Cannot open the ", A, 1x, "file",/)') trim(gridfile)
       stop 299
300   continue

! stop 0

return

end subroutine grid_lonlat_input

subroutine interpol_lonlat(bhi,bhr,loncrs,latcrs,londot,latdot,data,datan,datanUV,name,&
                           index,indexUV,wiscrs,wisdot,wagdot,lmore)

!
! This utility is experimental
! This utility interpolate (bilinial) 
! grids from MM5 (Lambert Conformal
! and Polar Stereographic) to regular
!  longitude and latitude grid
!
! Author: Alfons Callado Pallares
! E-mail: acallado@inm.es
!

implicit none

common /interpol/ lweight,lcross,lUV,end_index,end_index_n,&
                  inter,start_lon,start_lat,interval_lon,interval_lat

logical :: lweight
logical :: lcross
logical :: lmore
integer :: lUV
integer, dimension(4) :: end_index, end_index_n
real, dimension(end_index(1),end_index(2)) :: loncrs,latcrs
real, dimension(end_index(1),end_index(2)) :: londot,latdot
real, dimension(end_index_n(1),end_index_n(2),4) :: wiscrs
real, dimension(end_index_n(1),end_index_n(2),4) :: wisdot
real, dimension(end_index_n(1),end_index_n(2),1) :: wagdot
real, dimension(end_index(1),end_index(2),end_index(3),1) :: data
real, dimension(end_index_n(1), end_index_n(2), end_index_n(3), end_index_n(4)) :: datan
real, dimension(end_index_n(1), end_index_n(2), end_index_n(3), end_index_n(4), 2) :: datanUV
real, dimension(end_index(1),end_index(2)) :: databint
real :: bint
integer, dimension(end_index_n(1),end_index_n(2),4,2) :: index
integer, dimension(end_index_n(1),end_index_n(2),4,2) :: indexUV

integer :: inter
real :: start_lon,start_lat,interval_lon,interval_lat

character (len= 9) :: name

real :: pi, fr

integer :: bhi
real, dimension(20,20) :: bhr

integer :: i, j, k, l

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

TYPE(proj_info) :: proj,projUV

pi=3.14159265358
fr=pi/180.

! Read the grid_lonlat.input parameters

if (lweight) then
!    print*, bhi,bhr(3,1)
    call set_proj( bhi, bhr, projUV, londot(1,1), latdot(1,1))
!    print*,'Hola2', londot(1,1), latdot(1,1)
!    print*,'Hola2',projUV%code,projUV%lat1,projUV%lon1,projUV%polei,projUV%polej
    call set_proj( bhi, bhr, proj, loncrs(1,1), latcrs(1,1))

    if (lmore) write(*,'(/,8x, "Computing weights for cross points (interpolation step)",/)')
    call weightsa(proj,inter,start_lon,start_lat,interval_lon,interval_lat,end_index,end_index_n,&
                 latcrs,loncrs,wiscrs,index,.true.,wagdot,lmore)
    if (lmore) write(*,'(/,8x, "Computing weights for dot points (interpolation step)",/)')
    call weightsa(projUV,inter,start_lon,start_lat,interval_lon,interval_lat,end_index,end_index_n,&
                 latdot,londot,wisdot,indexUV,.false.,wagdot,lmore)

    lweight = .false.
!    stop

 endif
! stop

! Interpolate

if (lmore) write(*,'(8x, "Subroutine interpol_lonlat: Interpolating ... ",/)')
if (inter.eq.0) then
if (lcross) then
 if (lmore) print*, '       Subroutine interpol_lonlat: lUV(0 -> No UV; 1,2 -> UV) and lcross=', lUV, lcross
 datan=0.
! print*, 'Hello case (0)', end_index_n(3)
 do k=1,end_index_n(3)
 databint=0.
 databint=data(1:end_index(1),1:end_index(2),k,1)
 do i=1,end_index_n(1)
 do j=1,end_index_n(2)
    datan(i,j,k,1)=bint(wiscrs(i,j,1),wiscrs(i,j,2),databint,end_index(1),end_index(2),0)
!    print*, 'Hello0 case(0)', i,j,k,'1', datan(i,j,k,1), data(index(i,j,1,1),index(i,j,1,2),k,1),&
!                                         wiscrs(i,j,1), wiscrs(i,j,2)
! stop 0
 enddo
 enddo
 enddo
else
 if (lmore) print*, '       Subroutine interpol_lonlat: lUV(0 -> No UV; 1,2 -> UV)= and lcross', lUV, lcross
 if(lUV.eq.1) then
    datanUV= 0.
    datanUV= 0.
 endif
 if ((name.eq.'U        ').or.(name.eq.'U10      ')) then
 do k=1,end_index_n(3)
 databint=0.
 databint=data(1:end_index(1),1:end_index(2),k,1)
 do i=1,end_index_n(1)
 do j=1,end_index_n(2)
    datan(i,j,k,1)=bint(wiscrs(i,j,1),wiscrs(i,j,2),databint,end_index(1),end_index(2),0)
    datanUV(i,j,k,1,1) = datanUV(i,j,k,1,1) + (datan(i,j,k,1))*(cos(-1.*fr*(wagdot(i,j,1))))
    datanUV(i,j,k,1,2) = datanUV(i,j,k,1,2) + (datan(i,j,k,1))*(sin(-1.*fr*(wagdot(i,j,1))))
!    print*, 'Hello11 case (0)', i,j,k,'1', datan(i,j,k,1), index(i,j,1,1), index(i,j,1,2),&
!                                           data(index(i,j,1,1),index(i,j,1,2),k,1),&
!                                           wiscrs(i,j,1), wiscrs(i,j,2)
!    print*, 'Hello11 case (0)', i,j,k,'1', wagdot(i,j,1), datanUV(i,j,k,1,1), datanUV(i,j,k,1,2)
!    stop 0
 enddo
 enddo
 enddo
 elseif ((name.eq.'V        ').or.(name.eq.'V10      ')) then
 do k=1,end_index_n(3)
 databint=0.
 databint=data(1:end_index(1),1:end_index(2),k,1)
 do i=1,end_index_n(1)
 do j=1,end_index_n(2)
    datan(i,j,k,1)=bint(wiscrs(i,j,1),wiscrs(i,j,2),databint,end_index(1),end_index(2),0)
    datanUV(i,j,k,1,1) = datanUV(i,j,k,1,1) + (datan(i,j,k,1))*(sin(fr*(wagdot(i,j,1))))
    datanUV(i,j,k,1,2) = datanUV(i,j,k,1,2) + (datan(i,j,k,1))*(cos(fr*(wagdot(i,j,1))))
!    print*, 'Hello11 case (0)', i,j,k,'1', datan(i,j,k,1), index(i,j,1,1), index(i,j,1,2),&
!                                           data(index(i,j,1,1),index(i,j,1,2),k,1),&
!                                           wiscrs(i,j,1), wiscrs(i,j,2)
!    print*, 'Hello11 case (0)', i,j,k,'1', wagdot(i,j,1), datanUV(i,j,k,1,1), datanUV(i,j,k,1,2)
!    stop 0
 enddo
 enddo
 enddo
 endif
 endif
else
if (lcross) then
 if (lmore) print*, '       Subroutine interpol_lonlat: lUV(0 -> No UV; 1,2 -> UV) and lcross=', lUV, lcross
 datan=0.
! print*, 'Hello-0', end_index_n(3)
 do k=1,end_index_n(3)
 do i=1,end_index_n(1)
 do j=1,end_index_n(2)
    datan(i,j,k,1) = 0.
!    print*, 'Hello-1'
!    print*, 'Hello0', i,j,k,'1', datan(i,j,k,1), wiscrs(i,j,1), wiscrs(i,j,2), wiscrs(i,j,3),&
!                                 wiscrs(i,j,4)
!    print*, 'Hello0', i,j,k,'1', datan(i,j,k,1), wisdot(i,j,1), wisdot(i,j,2), wisdot(i,j,3),&
!                                 wisdot(i,j,4), wagdot(i,j,1)
!    print*, 'Hello0', i,j,k,'1', index(i,j,1,1), index(i,j,2,1), index(i,j,3,1), index(i,j,4,1)
!    print*, 'Hello0', i,j,k,'1', index(i,j,1,2), index(i,j,2,2), index(i,j,3,2), index(i,j,4,2)
!    stop
    do l=1,4
       datan(i,j,k,1) = datan(i,j,k,1) + wiscrs(i,j,l) * data(index(i,j,l,1),index(i,j,l,2),k,1)
    enddo
!    print*, 'Hello0', i,j,k,'1', datan(i,j,k,1), wiscrs(i,j,1), wiscrs(i,j,2), wiscrs(i,j,3),&
!                                 wiscrs(i,j,4)
 enddo
 enddo
 enddo
! print*, 'Hello01', i-1,j-1,k-1,'1', datan(i-1,j-1,k-1,1), wiscrs(i-1,j-1,1), wiscrs(i-1,j-1,2),&
!                                     wiscrs(i-1,j-1,3), wiscrs(i-1,j-1,4)
else
 if (lmore) print*, '       Subroutine interpol_lonlat: lUV(0 -> No UV; 1,2 -> UV)= and lcross', lUV, lcross
 if(lUV.eq.1) then
    datanUV= 0.
    datanUV= 0.
 endif
 if ((name.eq.'U        ').or.(name.eq.'U10      ')) then
 do k=1,end_index_n(3)
 do i=1,end_index_n(1)
 do j=1,end_index_n(2)
    datan(i,j,k,1) = 0.
    do l=1,4

       datan(i,j,k,1) = datan(i,j,k,1) + wisdot(i,j,l) * data(indexUV(i,j,l,1),indexUV(i,j,l,2),k,1)
!       datanUV(i,j,k,1,1)= datanUV(i,j,k,1,1) + (wisdot(i,j,l))*&
!                         (cos(-1.*fr*(wagdot(i,j,1))))*(data(indexUV(i,j,l,1),indexUV(i,j,l,2),k,1))
!       datanUV(i,j,k,1,2)= datanUV(i,j,k,1,2) + (wisdot(i,j,l))*&
!                         (sin(-1.*fr*(wagdot(i,j,1))))*(data(indexUV(i,j,l,1),indexUV(i,j,l,2),k,1))

    enddo
    datanUV(i,j,k,1,1) = datanUV(i,j,k,1,1) + (datan(i,j,k,1))*(cos(-1.*fr*(wagdot(i,j,1))))
    datanUV(i,j,k,1,2) = datanUV(i,j,k,1,2) + (datan(i,j,k,1))*(sin(-1.*fr*(wagdot(i,j,1)))) 

!    print*, 'Hello11', i,j,k,'1', datan(i,j,k,1), wagdot(i,j,1), datanUV(i,j,k,1,1), datanUV(i,j,k,1,2)
 enddo
 enddo
 enddo
! print*, 'Hello11', i-1,j-1,k-1,'1', datan(i-1,j-1,k-1,1), wagdot(i-1,j-1,1), datanUV(i-1,j-1,k-1,1,1),&
!                                    datanUV(i-1,j-1,k-1,1,2)
 elseif ((name.eq.'V        ').or.(name.eq.'V10      ')) then
 do k=1,end_index_n(3)
 do i=1,end_index_n(1)
 do j=1,end_index_n(2)
    datan(i,j,k,1) = 0.
    do l=1,4

       datan(i,j,k,1) = datan(i,j,k,1) + wisdot(i,j,l) * data(indexUV(i,j,l,1),indexUV(i,j,l,2),k,1)
!       datanUV(i,j,k,1,1)= datanUV(i,j,k,1,1) + (wisdot(i,j,l))*&
!                         (sin(fr*(wagdot(i,j,k))))*(data(indexUV(i,j,l,1),indexUV(i,j,l,2),k,1))
!       datanUV(i,j,k,1,2)= datanUV(i,j,k,1,2) + (wisdot(i,j,l))*&
!                         (cos(fr*(wagdot(i,j,k))))*(data(indexUV(i,j,l,1),indexUV(i,j,l,2),k,1))

    enddo
    datanUV(i,j,k,1,1) = datanUV(i,j,k,1,1) + (datan(i,j,k,1))*(sin(fr*(wagdot(i,j,1))))
    datanUV(i,j,k,1,2) = datanUV(i,j,k,1,2) + (datan(i,j,k,1))*(cos(fr*(wagdot(i,j,1))))
!    print*, 'Hello2', i,j,k,'1', datan(i,j,k,1), wagdot(i,j,1), datanUV(i,j,k,1,1), datanUV(i,j,k,1,2)
 enddo
 enddo
 enddo
! print*, 'Hello21', i-1,j-1,k-1,'1', datan(i-1,j-1,k-1,1), wagdot(i-1,j-1,1), datanUV(i-1,j-1,k-1,1,1),&
!                                     datanUV(i-1,j-1,k-1,1,2)

 endif
endif
endif

    return
!    stop
end subroutine interpol_lonlat  

SUBROUTINE set_proj( bhi, bhr, proj, LONGI,LATIT) 
! Initializes the map projection structure to missing values
! source: /LITTLE_R/src/map_utils_module.F 

   IMPLICIT NONE

   integer :: bhi 
   real, dimension(20,20) :: bhr

   real :: LONGI,LATIT

   REAL, PARAMETER :: earth_radius_m = 6370000.

   REAL :: ala1,alo1,reflon,scale_top

   REAL :: arg,deltalon1,tl1r,ctl1r

   REAL :: clain

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

   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

 !     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%rebydx   = -999.9
   proj%init     = .FALSE.

   proj%lat1=LATIT
   proj%lon1=LONGI

   proj%truelat1=BHR(5,1)
   proj%truelat2=BHR(6,1)

   proj%stdlon=BHR(3,1) 

   proj%dx=BHR(9,1)
   IF (proj%truelat1 .LT. 0.) THEN
     proj%hemi = -1.0
   ELSE
     proj%hemi = 1.0
   ENDIF
   proj%rebydx = earth_radius_m / proj%dx 

   proj%cone=BHR(4,1)

   select case (bhi) 
! Lambert Conformal MM5 projection
     case (1)
       proj%code=3

       deltalon1 = proj%lon1 - proj%stdlon
       IF (deltalon1 .GT. +180.) deltalon1 = deltalon1 - 360.
       IF (deltalon1 .LT. -180.) deltalon1 = deltalon1 + 360.

       tl1r = proj%truelat1 * rad_per_deg
       ctl1r = COS(tl1r)

       proj%rsw = proj%rebydx * ctl1r/proj%cone * &
              (TAN((90.*proj%hemi-proj%lat1)*rad_per_deg/2.) / &
               TAN((90.*proj%hemi-proj%truelat1)*rad_per_deg/2.))**proj%cone

       arg = proj%cone*(deltalon1*rad_per_deg)
       proj%polei = 1. - proj%hemi * proj%rsw * SIN(arg)
       proj%polej = 1. + proj%rsw * COS(arg)
!       PRINT '(A,2F10.3)', 'Computed pole (x,y) = ', proj%polei, proj%polej

! Stereographic polar MM5 projection
     case (2)
       proj%code=5

       reflon = proj%stdlon + 90.

       scale_top = 1. + proj%hemi * SIN(proj%truelat1 * rad_per_deg)

       ala1 = proj%lat1 * rad_per_deg
       proj%rsw = proj%rebydx*COS(ala1)*scale_top/(1.+proj%hemi*SIN(ala1))

       alo1 = (proj%lon1 - reflon) * rad_per_deg
       proj%polei = 1. - proj%rsw * COS(alo1)
       proj%polej = 1. - proj%hemi * proj%rsw * SIN(alo1)
!      PRINT '(A,2F10.1)', 'Computed (I,J) of pole point: ',proj%polei,proj%polej

! Mercator MM5 projection
     case (3)
       proj%code=1

!      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.
       IF (proj%lat1 .NE. 0.) THEN
         proj%rsw = (ALOG(TAN(0.5*((proj%lat1+90.)*rad_per_deg))))/proj%dlon
       ENDIF
   end select

   proj%init= .true.

!   print*, 'Set_proj', bhi,bhr(3,1)
!   print*, 'Set_proj', proj%code,proj%lat1,proj%lon1,proj%polei,proj%polej

END SUBROUTINE set_proj


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   FUNCTION bint(xx,yy,list,iii,jjj,ibint)
   
      !  Bi-linear interpolation among four grid values
   
      IMPLICIT NONE
   
      REAL :: xx , yy
      INTEGER :: ibint , iii, jjj
      REAL list(iii,jjj),stl(4,4)
   
      INTEGER :: ib , jb, n , i , j , k , kk , l , ll
      REAL :: bint , x , y , a , b , c , d , e , f , g , h

      real :: oned
   
      ib=iii-ibint
      jb=jjj-ibint
      bint = 0.0
      n = 0
      i = INT(xx+0.00001)
      j = INT(yy+0.00001)
      x = xx - i
      y = yy-j
   
      IF ( ( ABS(x).GT.0.0001 ) .OR. ( abs(y).gt.0.0001 ) ) THEN
         loop_1 : DO k = 1,4
            kk = i + k - 2
            IF ( kk .LT. 1) then 
               kk = 1
            END IF
            IF ( kk .GT. ib) THEN
               kk = ib
            END IF
            loop_2 : DO l = 1,4
               stl(k,l) = 0.
               ll = j + l - 2
               IF ( ll .LT. 1) THEN
                  ll = 1
               END IF
               IF ( ll .GT. jb) THEN
                  ll = jb
               END IF
               stl(k,l) = list(kk,ll)
               n = n + 1
               IF ( stl(k,l) .EQ. 0. ) THEN
                  stl(k,l) = 1.E-20
               END IF
            END DO loop_2
         END DO loop_1
   
         a = oned(x,stl(1,1),stl(2,1),stl(3,1),stl(4,1))
         b = oned(x,stl(1,2),stl(2,2),stl(3,2),stl(4,2))
         c = oned(x,stl(1,3),stl(2,3),stl(3,3),stl(4,3))
         d = oned(x,stl(1,4),stl(2,4),stl(3,4),stl(4,4))
         bint = oned(y,a,b,c,d)
   
         IF(n.NE.16) THEN
            e = oned(y,stl(1,1),stl(1,2),stl(1,3),stl(1,4))
            f = oned(y,stl(2,1),stl(2,2),stl(2,3),stl(2,4))
            g = oned(y,stl(3,1),stl(3,2),stl(3,3),stl(3,4))
            h = oned(y,stl(4,1),stl(4,2),stl(4,3),stl(4,4))
            bint = (bint+oned(x,e,f,g,h)) * 0.5
         END IF
   
      ELSE
         bint = list(i,j)
      END IF
   
   END FUNCTION bint

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

   real FUNCTION oned(x,a,b,c,d) 
   
      IMPLICIT NONE
   
      REAL :: x,a,b,c,d,oned
   
      oned = 0.                
   
      IF      ( x .EQ. 0. ) THEN
         oned = b      
      ELSE IF ( x .EQ. 1. ) THEN
         oned = c      
      END IF
   
      IF(b*c.NE.0.) THEN
         IF ( a*d .EQ. 0. ) THEN
            IF      ( ( a .EQ. 0 ) .AND. ( d .EQ. 0 ) ) THEN
               oned = b*(1.0-x)+c*x                                        
            ELSE IF ( a .NE. 0. ) THEN
               oned = b+x*(0.5*(c-a)+x*(0.5*(c+a)-b))            
            ELSE IF ( d .NE. 0. ) THEN
               oned = c+(1.0-x)*(0.5*(b-d)+(1.0-x)*(0.5*(b+d)-c)) 
            END IF
         ELSE
            oned = (1.0-x)*(b+x*(0.5*(c-a)+x*(0.5*(c+a)-b)))+x*(c+(1.0-x)*(0.5*(b-d)+(1.0-x)*(0.5*(b+d)-c)))                                   
         END IF
      END IF
   
   END FUNCTION oned                                                       

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
