CCC CONVERT MAPS LAMBERT CONFORMAL PROJECTION i,j TO LONG,LAT ccc using ncep routine w3fb12 real DX_P, LAT_LL_P, LON_LL_P, LON_XX_P, LAT_TAN_P CCC FOLLOWING FROM mapscon.f PARAMETER ( DX_P = 20317.63 ) PARAMETER ( LAT_LL_P = 16.2810 ) PARAMETER ( LON_LL_P = -126.1378 ) PARAMETER ( LON_XX_P = -95.0 ) C** LAT_TAN_P R LATITUDE AT LAMBERT CONFORMAL PROJECTION IS TRUE (DEG) PARAMETER ( LAT_TAN_P = 25.0 ) cunusedC** LAT_TRUE_P R LATITUDE AT WHICH X-Y SCALE IS TRUE (DEG) cunused real PARAMETER ( LAT_TRUE_P = 35.0 ) *** IF LCONVERT=1, CONVERT FROM IMAGE i,j CCC WITH SW,NE CORNERS AT (iimagemin,jimagemax),(imagemax,jimagemin) CCC CORRESPONDING TO MAPS CORNERS (imapsmin,jmapsmin),(mapsmax,jmapsmax) c--- data LCONVERT / 1 / data LCONVERT / 0 / ********** NB *** SHOULD CONVERT TO GRID 3 !!!! **************** ccc for ca/nv grid 3 c--- data imapsmin,jmapsmin, imapsmax,jmapsmax / 38,78, 98,136 / ccc for ca/nv grid 2b (expanded_to_CO) *** TEMPORARY data imapsmin,jmapsmin, imapsmax,jmapsmax / 37,86, 107,136 / ccc for ca/nv grid 1 c--- data imapsmin,jmapsmin, imapsmax,jmapsmax / 33,93, 73,143 / c-for_ghostview_image_coords (approx 2b) data iimagemin,jimagemin, iimagemax,jimagemax / 68,255, 538,588 / c-for_reversed_image_coords data iimagemin,jimagemax, iimagemax,jimagemin / 102,255, 538,588 / CCC FOR CALC OF KPT - MAPS 20km GRID data maxx / 301 / data deg2rad / 0.017453292 / character qlongdegminit*8, qlatdegminit*6 print *,'FOR 20km MAPS LAMBERT CONFORMAL GRID - (1,1)=SW_corner' if ( LCONVERT .eq. 1 ) then print *,'WARNING - LCONVERT=1 => images corners used =', & imapsmin,jmapsmin, imapsmax,jmapsmax endif maxiter=100 do i=1,maxiter if ( LCONVERT .ne. 1 ) then print *,'INPUT maps i,j (fp)(0,0=end)' read(5,*) ai,aj if( ai.eq.0.0 .and. aj.eq.0.0 ) stop else print *,'INPUT *image* i,j (fp)(0,0=end)' read(5,*) aiimage,ajimage if( aiimage.eq.0.0 .and. ajimage.eq.0.0 ) stop ai = imapsmin + float(imapsmax-imapsmin)* & (aiimage-iimagemin)/float(iimagemax-iimagemin) aj = jmapsmin + float(jmapsmax-jmapsmin)* & (ajimage-jimagemin)/float(jimagemax-jimagemin) print *,' =>ai,aj=',ai,aj endif ccc get flatfile pt number akpt = ai +(aj-1)*maxx call W3FB12( aI,aJ, LAT_LL_P,LON_LL_P,DX_P,LON_XX_P,LAT_TAN_P, & ALAT,aLON, IERR ) if( ierr.ne.0 ) then print *,'ERROR: ierr=', ierr cycle endif ccc convert to W longitude alon = alon - 360. ccc get long,lat in degrees and minutes ilongdeg = int( abs(alon) ) ilongminit = nint( ( abs(alon) - ilongdeg )*60.) write(qlongdegminit,'(i4,"+",i2,"m")') ilongdeg,ilongminit ilatdeg = int( alat ) ilatminit = nint( ( abs(alat) - ilatdeg )*60.) write(qlatdegminit,'(i2,"+",i2,"m")') ilatdeg,ilatminit CCC DETERMINE ANGLE OF TRUE NORTH RELATIVE TO Y AXIS (ccw=+) ccc (taken from blip.pl) ROTCON_P = 0.422618 angle2north = ROTCON_P * ( alon-LON_XX_P ) * deg2rad angle2northdeg = angle2north / deg2rad write(6,6100) akpt,alat,alon, qlongdegminit,qlatdegminit, & angle2northdeg 6100 format('KTP LAT,LONG=',f10.2,2x,2f12.4,3x,a8,1x,a6,f10.1 /) C-------------------------------------------------------------------- CCC USED THIS FOR LOOKING INTO ROTATING BLIPMAP TO INPUT INTO SEEYOU CCC DETERMINE COORDS IN ROTATED FRAME (cw=+) CCC airot,aijrot = rotated coords with airot=ajrot=0 at aicent,ajcent CCC rotate by rotdeg (deg) about aicent,ajcent = center of rotation (normal coords) CCC BELOW IS FOR CA-NV GRID CENTER c rotdeg = -10.0 c aicent = 60.0 c ajcent = 107.5 c rotrad = rotdeg * deg2rad c airot = (ai-aicent)*cos(rotrad) +(aj-ajcent)*sin(rotrad) c ajrot = -(ai-aicent)*sin(rotrad) +(aj-ajcent)*cos(rotrad) c write(6,6100) akpt,alon,alat, qlongdegminit,qlatdegminit, c & angle2northdeg, c & airot,ajrot c 6100 format('KTP LONG,LAT=',f10.2,2x,2f10.3,3x,a8,1x,a6,f10.1, c & 2x,2f7.1 /) C-------------------------------------------------------------------- enddo stop 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 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: IBM VS FORTRAN C MACHINE: NAS C C$$$ C LOGICAL NEWMAP DATA RERTH /6.3712E+6/, PI/3.1415926/, OLDRML/99999./ C C PRELIMINARY VARIABLES AND REDIFINITIONS C C H = 1 FOR NORTHERN HEMISPHERE; = -1 FOR SOUTHERN C SAVE BETA = 1. IERR = 0 IF(ALATAN.GT.0) THEN H = 1. ELSE H = -1. ENDIF C PIBY2 = PI/2. RADPD = PI/180.0 DEGPRD = 1./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. 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