program main c c calculate pdef bilin file for ncep grid 218, a lambert c grid. the xdef/ydef that goes with this is: c xdef 700 linear 220 0.11 c ydef 500 linear 12 0.11 c real xi(700,500),xj(700,500),wrot(700,500) c do 10 j=1,500 rlat = 12.0+float(j-1)*0.11 do 10 i=1,700 rlon = 220.0+float(i-1)*0.11 call w3fb11(rlat,rlon,12.190,226.541,12191.,265.0,25.0, & xi(i,j),xj(i,j)) if (xi(i,j).lt.1.0 .or. xi(i,j).gt.614.0 .or. & xj(i,j).lt.1.0 .or. xj(i,j).gt.428.0) then xi(i,j) = -999.0 xj(i,j) = -999.0 wrot(i,j) = -999.0 else call w3fc03 (rlon,265.0,25.0,wrot(i,j)) endif 10 continue c open (8,file="grib218.pdef",form="unformatted") c c important to do seperate writes, so that grads c can figure it out c write(8) xi write(8) xj write(8) wrot stop end 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 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 c SUBROUTINE W3FC03(SLON,CENLON,XLAT1,angle) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: W3FC03 GRID U,V WIND COMPS. TO DIR. AND SPEED C AUTHOR: ROGERS/BRILL ORG: WD22 DATE: 00-03-27 C C ABSTRACT: GIVEN THE GRID-ORIENTED WIND COMPONENTS ON A C LAMBERT CONFORMAL GRID POINT, COMPUTE THE DIRECTION C AND SPEED OF THE WIND AT THAT POINT. INPUT WINDS AT THE NORTH C POLE POINT ARE ASSUMED TO HAVE THEIR COMPONENTS FOLLOW THE WMO C STANDARDS FOR REPORTING WINDS AT THE NORTH POLE. C (SEE OFFICE NOTE 241 FOR WMO DEFINITION). OUTPUT DIRECTION C WILL FOLLOW WMO CONVENTION. C C PROGRAM HISTORY LOG: C 00-03-27 BRILL/ROGERS c c 00-10-30 doty: modified to just return rotation angle C C USAGE: CALL W3FC03 (SLON,FGU,FGV,CENLON,XLAT1,DIR,SPD) C C INPUT VARIABLES: C NAMES INTERFACE DESCRIPTION OF VARIABLES AND TYPES C ------ --------- ----------------------------------------------- C SLON ARG LIST REAL*4 STATION LONGITUDE (-DEG W) C CENLON ARG LIST REAL*4 CENTRAL LONGITUDE C XLAT1 ARG LIST REAL*4 TRUE LATITUDE #1 C C OUTPUT VARIABLES: C NAMES INTERFACE DESCRIPTION OF VARIABLES AND TYPES C ------ --------- ----------------------------------------------- c angle arg list rotation angle C C SUBPROGRAMS CALLED: C NAMES LIBRARY C ------------------------------------------------------- -------- C ABS ACOS ATAN2 SQRT SYSTEM C C WARNING: THIS JOB WILL NOT VECTORIZE ON A CRAY C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C MACHINE: IBM SP C C$$$ PARAMETER (DTR=3.1415926/180.0) C C COMPUTE CONSTANT OF CONE C COCON = SIN(XLAT1*DTR) ANGLE = COCON * (cenlon-slon) * DTR c ANGLE = COCON * (SLON-CENLON) * DTR RETURN END