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 90-04-12 R.E.JONES CONVERT TO CFT77 FORTRAN
C 94-04-28 R.E.JONES ADD SAVE STATEMENT
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: CRAY CFT77 FORTRAN
C MACHINE: CRAY C916-128, CRAY Y-MP8/864, CRAY Y-MP EL2/256
C
C$$$
C
SAVE
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
SUBROUTINE W3FB12(XI,XJ,ALAT1,ELON1,DX,ELONV,ALATAN,ALAT,ELON,
& IERR)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: W3FB12 LAMBERT(I,J) TO LAT/LON FOR GRIB
C PRGMMR: STACKPOLE ORG: NMC42 DATE:88-11-28
C
C ABSTRACT: CONVERTS THE COORDINATES OF A LOCATION ON EARTH GIVEN IN A
C GRID COORDINATE SYSTEM OVERLAID ON A LAMBERT CONFORMAL TANGENT
C CONE PROJECTION TRUE AT A GIVEN N OR S LATITUDE TO THE
C NATURAL COORDINATE SYSTEM OF LATITUDE/LONGITUDE
C W3FB12 IS THE REVERSE OF W3FB11.
C USES GRIB SPECIFICATION OF THE LOCATION OF THE GRID
C
C PROGRAM HISTORY LOG:
C 88-11-25 ORIGINAL AUTHOR: STACKPOLE, W/NMC42
C 90-04-12 R.E.JONES CONVERT TO CFT77 FORTRAN
C 94-04-28 R.E.JONES ADD SAVE STATEMENT
C
C USAGE: CALL W3FB12(XI,XJ,ALAT1,ELON1,DX,ELONV,ALATAN,ALAT,ELON,IERR,
C IERR)
C INPUT ARGUMENT LIST:
C XI - I COORDINATE OF THE POINT REAL*4
C XJ - J COORDINATE OF THE POINT REAL*4
C ALAT1 - LATITUDE OF LOWER LEFT POINT OF GRID (POINT 1,1)
C LATITUDE <0 FOR SOUTHERN HEMISPHERE; REAL*4
C ELON1 - LONGITUDE OF LOWER LEFT POINT OF GRID (POINT 1,1)
C EAST LONGITUDE USED THROUGHOUT; 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 THE GRID) ALONG WHICH LATITUDE INCREASES AS
C THE Y-COORDINATE INCREASES. REAL*4
C THIS IS ALSO THE MERIDIAN (ON THE OTHER 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 (TOUCHES OR OSCULATES) THE SPHERICAL EARTH.
C SET NEGATIVE TO INDICATE A
C SOUTHERN HEMISPHERE PROJECTION; REAL*4
C
C OUTPUT ARGUMENT LIST:
C ALAT - LATITUDE IN DEGREES (NEGATIVE IN SOUTHERN HEMI.)
C ELON - EAST LONGITUDE IN DEGREES, REAL*4
C IERR - .EQ. 0 IF NO PROBLEM
C .GE. 1 IF THE REQUESTED XI,XJ POINT IS IN THE
C FORBIDDEN ZONE, I.E. OFF THE LAMBERT MAP
C IN THE OPEN SPACE WHERE THE CONE IS CUT.
C IF IERR.GE.1 THEN ALAT=999. AND ELON=999.
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: CRAY CFT77 FORTRAN
C MACHINE: CRAY C916-128, CRAY Y-MP8/864, CRAY Y-MP EL2/256
C
C$$$
C
LOGICAL NEWMAP
C
SAVE
C
DATA RERTH /6.3712E+6/, PI/3.14159/, OLDRML/99999./
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
PIBY2 = PI / 2.0
RADPD = PI / 180.0
DEGPRD = 1.0 / RADPD
REBYDX = RERTH / DX
ALATN1 = ALATAN * RADPD
AN = H * SIN(ALATN1)
COSLTN = COS(ALATN1)
C
C MAKE SURE THAT INPUT LONGITUDE DOES 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
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 RMLL TO TEST IF MAP AND GRID UNCHANGED FROM PREVIOUS
C CALL TO THIS CODE. THUS AVOID UNNEEDED RECOMPUTATIONS.
C
IF (RMLL.EQ.OLDRML) THEN
NEWMAP = .FALSE.
ELSE
NEWMAP = .TRUE.
OLDRML = RMLL
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)
ENDIF
C
C RADIUS TO THE I,J POINT (IN GRID UNITS)
C YY REVERSED SO POSITIVE IS DOWN
C
XX = XI - POLEI
YY = POLEJ - XJ
R2 = XX**2 + YY**2
C
C CHECK THAT THE REQUESTED I,J IS NOT IN THE FORBIDDEN ZONE
C YY MUST BE POSITIVE UP FOR THIS TEST
C
THETA = PI*(1.-AN)
BETA = ABS(ATAN2(XX,-YY))
IERR = 0
IF (BETA.LE.THETA) THEN
IERR = 1
ALAT = 999.
ELON = 999.
IF (.NOT.NEWMAP) RETURN
ENDIF
C
C NOW THE MAGIC FORMULAE
C
IF (R2.EQ.0) THEN
ALAT = H * 90.0
ELON = ELONV
ELSE
C
C FIRST THE LONGITUDE
C
ELON = ELONV + DEGPRD * ATAN2(H*XX,YY)/AN
ELON = AMOD(ELON+360., 360.)
C
C NOW THE LATITUDE
C RECALCULATE THE THING ONLY IF MAP IS NEW SINCE LAST TIME
C
IF (NEWMAP) THEN
ANINV = 1./AN
ANINV2 = ANINV/2.
THING = ((AN/REBYDX) ** ANINV)/
& ((COSLTN**((1.-AN)*ANINV))*(1.+ AN))
ENDIF
ALAT = H*(PIBY2 - 2.*ATAN(THING*(R2**ANINV2)))*DEGPRD
ENDIF
C
C FOLLOWING TO ASSURE ERROR VALUES IF FIRST TIME THRU
C IS OFF THE MAP
C
IF (IERR.NE.0) THEN
ALAT = 999.
ELON = 999.
IERR = 2
ENDIF
RETURN
END
SUBROUTINE W3FB08(ALAT,ALON,ALAT1,ALON1,ALATIN,DX,XI,XJ)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: W3FB08 LAT/LON TO MERC (I,J) FOR GRIB
C PRGMMR: STACKPOLE ORG: NMC42 DATE:88-04-05
C
C ABSTRACT: CONVERTS A LOCATION ON EARTH GIVEN IN
C THE COORDINATE SYSTEM OF LATITUDE/LONGITUDE TO AN (I,J)
C COORDINATE SYSTEM OVERLAID ON A MERCATOR MAP PROJECTION
C W3FB08 IS THE REVERSE OF W3FB09
C USES GRIB SPECIFICATION OF THE LOCATION OF THE GRID
C
C PROGRAM HISTORY LOG:
C 88-03-01 ORIGINAL AUTHOR: STACKPOLE, W/NMC42
C 90-04-12 R.E.JONES CONVERT TO CRAY CFT77 FORTRAN
C
C USAGE: CALL W3FB08 (ALAT,ALON,ALAT1,ALON1,ALATIN,DX,XI,XJ)
C INPUT ARGUMENT LIST:
C ALAT - LATITUDE IN DEGREES (NEGATIVE IN SOUTHERN HEMIS)
C ALON - EAST LONGITUDE IN DEGREES, REAL*4
C ALAT1 - LATITUDE OF LOWER LEFT CORNER OF GRID (POINT (1,1))
C ALON1 - LONGITUDE OF LOWER LEFT CORNER OF GRID (POINT (1,1))
C ALL REAL*4
C ALATIN - THE LATITUDE AT WHICH THE MERCATOR CYLINDER
C INTERSECTS THE EARTH
C DX - MESH LENGTH OF GRID IN METERS AT ALATIN
C
C OUTPUT ARGUMENT LIST:
C XI - I COORDINATE OF THE POINT SPECIFIED BY ALAT, ALON
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: CRAY CFT77 FORTRAN
C MACHINE: CRAY Y-MP8/832
C
C$$$
C
DATA RERTH /6.3712E+6/, PI/3.1416/
C
C PRELIMINARY VARIABLES AND REDIFINITIONS
C
RADPD = PI / 180.0
DEGPR = 180.0 / PI
CLAIN = COS(RADPD*ALATIN)
DELLON = DX / (RERTH*CLAIN)
C
C GET DISTANCE FROM EQUATOR TO ORIGIN ALAT1
C
DJEO = 0.
IF (ALAT1.NE.0.)
& DJEO = (ALOG(TAN(0.5*((ALAT1+90.0)*RADPD))))/DELLON
C
C NOW THE I AND J COORDINATES
C
XI = 1. + ((ALON - ALON1)/(DELLON*DEGPR))
XJ = 1. + (ALOG(TAN(0.5*((ALAT + 90.) * RADPD))))/
& DELLON
& - DJEO
C
RETURN
END
SUBROUTINE W3FB06(ALAT,ALON,ALAT1,ALON1,DX,ALONV,XI,XJ)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: W3FB06 LAT/LON TO POLA (I,J) FOR GRIB
C PRGMMR: STACKPOLE ORG: NMC42 DATE:88-04-05
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 POLAR STEREOGRAPHIC MAP PRO-
C JECTION TRUE AT 60 DEGREES N OR S LATITUDE. W3FB06 IS THE REVERSE
C OF W3FB07. USES GRIB SPECIFICATION OF THE LOCATION OF THE GRID
C
C PROGRAM HISTORY LOG:
C 88-01-01 ORIGINAL AUTHOR: STACKPOLE, W/NMC42
C 90-04-12 R.E.JONES CONVERT TO CRAY CFT77 FORTRAN
C
C USAGE: CALL W3FB06 (ALAT,ALON,ALAT1,ALON1,DX,ALONV,XI,XJ)
C INPUT ARGUMENT LIST:
C ALAT - LATITUDE IN DEGREES (NEGATIVE IN SOUTHERN HEMIS)
C ALON - EAST LONGITUDE IN DEGREES, REAL*4
C ALAT1 - LATITUDE OF LOWER LEFT POINT OF GRID (POINT (1,1))
C ALON1 - LONGITUDE OF LOWER LEFT POINT OF GRID (POINT (1,1))
C ALL REAL*4
C DX - MESH LENGTH OF GRID IN METERS AT 60 DEG LAT
C MUST BE SET NEGATIVE IF USING
C SOUTHERN HEMISPHERE PROJECTION.
C 190500.0 LFM GRID,
C 381000.0 NH PE GRID, -381000.0 SH PE GRID, ETC.
C ALONV - 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 FOR EXAMPLE:
C 255.0 FOR LFM GRID,
C 280.0 NH PE GRID, 100.0 SH PE GRID, ETC.
C
C OUTPUT ARGUMENT LIST:
C XI - I COORDINATE OF THE POINT SPECIFIED BY ALAT, ALON
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: CRAY CFT77 FORTRAN
C MACHINE: CRAY Y-MP8/832
C
C$$$
C
DATA RERTH /6.3712E+6/, PI/3.1416/
DATA SS60 /1.86603/
C
C PRELIMINARY VARIABLES AND REDIFINITIONS
C
C H = 1 FOR NORTHERN HEMISPHERE; = -1 FOR SOUTHERN
C
C REFLON IS LONGITUDE UPON WHICH THE POSITIVE X-COORDINATE
C DRAWN THROUGH THE POLE AND TO THE RIGHT LIES
C ROTATED AROUND FROM ORIENTATION (Y-COORDINATE) LONGITUDE
C DIFFERENTLY IN EACH HEMISPHERE
C
IF (DX.LT.0) THEN
H = -1.0
DXL = -DX
REFLON = ALONV - 90.0
ELSE
H = 1.0
DXL = DX
REFLON = ALONV - 270.0
ENDIF
C
RADPD = PI / 180.0
REBYDX = RERTH/DXL
C
C RADIUS TO LOWER LEFT HAND (LL) CORNER
C
ALA1 = ALAT1 * RADPD
RMLL = REBYDX * COS(ALA1) * SS60/(1. + H * SIN(ALA1))
C
C USE LL POINT INFO TO LOCATE POLE POINT
C
ALO1 = (ALON1 - REFLON) * RADPD
POLEI = 1. - RMLL * COS(ALO1)
POLEJ = 1. - H * RMLL * SIN(ALO1)
C
C RADIUS TO DESIRED POINT AND THE I J TOO
C
ALA = ALAT * RADPD
RM = REBYDX * COS(ALA) * SS60/(1. + H * SIN(ALA))
C
ALO = (ALON - REFLON) * RADPD
XI = POLEI + RM * COS(ALO)
XJ = POLEJ + H * RM * SIN(ALO)
C
RETURN
END
subroutine baciolgwv()
print *, ' GWVX REPLACE WITH A REAL BACIOL ROUTINE'
stop 99
end
subroutine dgeff()
print *, ' GWVX REPLACE WITH A REAL DGEF ROUTINE'
stop 99
end
subroutine dgesf()
print *, ' GWVX REPLACE WITH A REAL DGES ROUTINE'
stop 99
end
! subroutine ludcmp()
! print *, ' GWVX REPLACE WITH A REAL LUDCMP ROUTINE'
! stop 99
! end
! subroutine lubksb()
! print *, ' GWVX REPLACE WITH A REAL LUbksb ROUTINE'
! stop 99
! end
subroutine dgef()
print *,' GWVX REPLACE WITH A REAL DGEF'
end
subroutine dges()
print *,' GWVX REPLACE WITH A REAL DGEs'
end