C
C     ***please read the REAME file i40_60.README before you try this***

C     This program serves as an example for converting the 40km RUC-2
C     output toi the 60-km RUC-1 grid. It essentially ONLY does the horizontal
C     remapping...that is, if you want to do a vertical interpolation for 3-D
C     grids, you'll need to do that on your own. If you want to do a coordinate
C     transformation from hybrid-B to pressure, see the info on the web site
C     for a routine called "hb2p."
C
C     Keep in mind that what we are doing is to take output on a Lambert 
C     Conformal grid(RUC-2) and translate it onto a polar sterographic
C     projection grid (RUC-1). The main routine to do this is an NMC routine
C     called W3FB11 which is included in this file. You will need the latitude
C     and longitudes of the 60km grid points which we have out on the web site.
C     Basically you can read all about the NMC routine; but given the correct
C     parameters about the 40-km RUIC-2 grid (see below), it will take
C     latitudes and longitudes of grid points and convert them to (i,j)'s
C     in a lambert conformal grid.
C
C     This program serves as an EXAMPLE of how you would do the RUC-2 to
C     RUC-1 remapping. The example used here is for precipitation data,
C     but you could easily substitute any 2-D field (or any 1 level of a 3-D
C     field instead.
C
C     ***please read the file i40_60.README before trying this routine***
C
C
C                                                          Barry Schwartz
C                                                          NOAA/FSL
C                                                          March 1998
C_____________________________________________________________________________
C
C     Program starts here.....
C
      parameter(nx=81,ny=62) !60 km grid dimensions
c 
c     parameters for W3fb11 (specifying the Conus 212R grid)
c                           (***longitudes are East not west!!***)
c
      parameter (alat1=16.281,elon1=233.862) !!lon = -126.1378
      parameter(dx=40635.25) !grid mesh
      parameter(elonv=265.00) !this is -95W + 360
      parameter(alatan=25.00) !true
c
      real lat(nx,ny),lon(nx,ny)
      real ix(nx,ny),jy(nx,ny)
      real maps_precip60(nx,ny) !precip interpolated to 60km maps grid
      real precip(151,113) !precip on the 40km (RUC-2) grid
c
c     start by reading in the lat/lons of the MAPS grid points
C
C     you would change the file name to your file (received from us)
C     of the file containing the lat/lons of the 60km grid points
C
      open(unit=10,file='latlon.mps',status='old',form='unformatted')
C
      read(10) lat,lon
      close(10)
c
c     now get each grid point's (i,j) for each 60-km (RUC-1) grid point
c     using NMC routine W3FB11
c
      do 100 j=1,ny
      do  75 i=1,nx
C
      lon(i,j)=lon(i,j) + 360. !this convert'swest longitude to east longitude
c
      call w3fb11(lat(i,j),lon(i,j),alat1,elon1,dx,elonv,alatan,
     *ix(i,j),jy(i,j))
c
c
C USAGE:  CALL W3FB11 (ALAT,ELON,ALAT1,ELON1,DX,ELONV,ALATAN,XI,XJ)
C  
C     see W3FB11 routine velow for all definitions
C
 75   continue
 100  continue
c
c     now bring in some data you want to remap. In this case we try precip
c     data...you need to change the input file to read your (151,113) RUC-2
c     2-D array!
c
      open(unit=25,'my_40km_data',status='old',form='unformatted')
c
c     read in the data any way you like; array must be (151,113)
c     ours is unformattted but it doesn't matter
c
      read(25) precip
      close(25)

c     now find the precip values for the 60-km (RUC-1) grid points
   
      do 300 i=1,nx
      do 250 j=1,ny
c
c
      ii=ix(i,j) + 0.5 !get the nearest grid point precip value
      jj=jy(i,j) + 0.5
c
      nvalue = (precip(ii,jj) + 0.005)*100
      maps_precip60(i,j) = nvalue/100. !!precip interpolated to the RUC-1 grid
c
c     whatever.....

 250  continue
 300  continue
c
c
c     etc etc....up to you what to do..with your data...output to plotting
c     routine..write a file...your deal not mine....
c
      stop
      end
c
c**********************************************************************
       SUBROUTINE W3FB11(ALAT,ELON,ALAT1,ELON1,DX,ELONV,ALATAN,XI,XJ)
C$$$   SUBPROGRAM  DOCUMENTATION  BLOCK
C
C SUBPROGRAM:  W3FB11        LAT/LON TO LAMBERT(I,J) FOR GRIB
C   PRGMMR: STACKPOLE        ORG: NMC42       DATE:88-11-28
C
C ABSTRACT: CONVERTS THE COORDINATES OF A LOCATION ON EARTH GIVEN IN
C   THE NATURAL COORDINATE SYSTEM OF LATITUDE/LONGITUDE TO A GRID
C   COORDINATE SYSTEM OVERLAID ON A LAMBERT CONFORMAL TANGENT CONE
C   PROJECTION TRUE AT A GIVEN N OR S LATITUDE. W3FB11 IS THE REVERSE
C   OF W3FB12. USES GRIB SPECIFICATION OF THE LOCATION OF THE GRID
C
C PROGRAM HISTORY LOG:
C   88-11-25  ORIGINAL AUTHOR:  STACKPOLE, W/NMC42
C
C USAGE:  CALL W3FB11 (ALAT,ELON,ALAT1,ELON1,DX,ELONV,ALATAN,XI,XJ)
C   INPUT ARGUMENT LIST:
C     ALAT     - LATITUDE IN DEGREES (NEGATIVE IN SOUTHERN HEMIS)
C     ELON     - EAST LONGITUDE IN DEGREES, REAL*4
C     ALAT1    - LATITUDE  OF LOWER LEFT POINT OF GRID (POINT (1,1))
C     ELON1    - LONGITUDE OF LOWER LEFT POINT OF GRID (POINT (1,1))
C                ALL REAL*4
C     DX       - MESH LENGTH OF GRID IN METERS AT TANGENT LATITUDE
C     ELONV    - THE ORIENTATION OF THE GRID.  I.E.,
C                THE EAST LONGITUDE VALUE OF THE VERTICAL MERIDIAN
C                WHICH IS PARALLEL TO THE Y-AXIS (OR COLUMNS OF
C                OF THE GRID) ALONG WHICH LATITUDE INCREASES AS
C                THE Y-COORDINATE INCREASES.  REAL*4
C                THIS IS ALSO THE MERIDIAN (ON THE BACK SIDE OF THE
C                TANGENT CONE) ALONG WHICH THE CUT IS MADE TO LAY
C                THE CONE FLAT.
C     ALATAN   - THE LATITUDE AT WHICH THE LAMBERT CONE IS TANGENT TO
C                (TOUCHING) THE SPHERICAL EARTH.
C                 SET NEGATIVE TO INDICATE A
C                 SOUTHERN HEMISPHERE PROJECTION.
C
C   OUTPUT ARGUMENT LIST:
C     XI       - I COORDINATE OF THE POINT SPECIFIED BY ALAT, ELON
C     XJ       - J COORDINATE OF THE POINT; BOTH REAL*4
C
C   REMARKS: FORMULAE AND NOTATION LOOSELY BASED ON HOKE, HAYES,
C     AND RENNINGER'S "MAP PROJECTIONS AND GRID SYSTEMS...", MARCH 1981
C     AFGWC/TN-79/003
C
C ATTRIBUTES:
C   LANGUAGE: IBM VS FORTRAN
C   MACHINE:  NAS
C
C$$$
C
         DATA  RERTH /6.3712E+6/, PI/3.14159/
C
C        PRELIMINARY VARIABLES AND REDIFINITIONS
C
C        H = 1 FOR NORTHERN HEMISPHERE; = -1 FOR SOUTHERN
C
         IF(ALATAN.GT.0) THEN
           H = 1.
         ELSE
           H = -1.
         ENDIF
C
         RADPD = PI/180.0
         REBYDX = RERTH/DX
         ALATN1 = ALATAN * RADPD
         AN = H * SIN(ALATN1)
         COSLTN = COS(ALATN1)
C
C        MAKE SURE THAT INPUT LONGITUDES DO NOT PASS THROUGH
C        THE CUT ZONE (FORBIDDEN TERRITORY) OF THE FLAT MAP
C        AS MEASURED FROM THE VERTICAL (REFERENCE) LONGITUDE.
C
         ELON1L = ELON1
         IF((ELON1 - ELONV).GT.180.)
     &     ELON1L = ELON1 - 360.
         IF((ELON1 - ELONV).LT.(-180.))
     &     ELON1L = ELON1 + 360.
C
         ELONL = ELON
         IF((ELON  - ELONV).GT.180.)
     &     ELONL  = ELON  - 360.
         IF((ELON - ELONV).LT.(-180.))
     &     ELONL = ELON + 360.
C
         ELONVR = ELONV *RADPD
C
C        RADIUS TO LOWER LEFT HAND (LL) CORNER
C
         ALA1 =  ALAT1 * RADPD
         RMLL = REBYDX * (((COSLTN)**(1.-AN))*(1.+AN)**AN) *
     &           (((COS(ALA1))/(1.+H*SIN(ALA1)))**AN)/AN
C
C        USE LL POINT INFO TO LOCATE POLE POINT
C
         ELO1 = ELON1L * RADPD
         ARG = AN * (ELO1-ELONVR)
         POLEI = 1. - H * RMLL * SIN(ARG)
         POLEJ = 1. + RMLL * COS(ARG)
C
C        RADIUS TO DESIRED POINT AND THE I J TOO
C
         ALA =  ALAT * RADPD
         RM = REBYDX * ((COSLTN**(1.-AN))*(1.+AN)**AN) *
     &           (((COS(ALA))/(1.+H*SIN(ALA)))**AN)/AN
C
         ELO = ELONL * RADPD
         ARG = AN*(ELO-ELONVR)
         XI = POLEI + H * RM * SIN(ARG)
         XJ = POLEJ - RM * COS(ARG)
C
C        IF COORDINATE LESS THAN 1
C        COMPENSATE FOR ORIGIN AT (1,1)
C
         IF(XI.LT.1.)  XI = XI - 1.
         IF(XJ.LT.1.)  XJ = XJ - 1.
C
      RETURN
      END
c************************************************************************