awips2/nativeLib/rary.ohd.pproc.gribit/TEXT/w3fb11.f
root 9bb8decbcf Initial revision of AWIPS2 11.9.0-7p5
Former-commit-id: 133dc97f67 [formerly a02aeb236c] [formerly 9f19e3f712] [formerly 06a8b51d6d [formerly 9f19e3f712 [formerly 64fa9254b946eae7e61bbc3f513b7c3696c4f54f]]]
Former-commit-id: 06a8b51d6d
Former-commit-id: 377dcd10b9 [formerly 3360eb6c5f]
Former-commit-id: 8e80217e59
2012-01-06 08:55:05 -06:00

132 lines
4.6 KiB
Fortran

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
C ================================= RCS keyword statements ==========
CHARACTER*68 RCSKW1,RCSKW2
DATA RCSKW1,RCSKW2 / '
.$Source: /fs/hseb/ob72/wfo_rfc/precip_proc/source/gribit/src/RCS/w3fb11.f,v $
. $', '
.$Id: w3fb11.f,v 1.1 2006/05/03 13:44:00 gsood Exp $
. $' /
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