C----------------------------------------------------------------------- SUBROUTINE GDSWZD04(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT,LMAP,XLON,XLAT,YLON,YLAT,AREA) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: GDSWZD04 GDS WIZARD FOR GAUSSIAN CYLINDRICAL C PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10 C C ABSTRACT: THIS SUBPROGRAM DECODES THE GRIB GRID DESCRIPTION SECTION C (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63) C AND RETURNS ONE OF THE FOLLOWING: C (IOPT=+1) EARTH COORDINATES OF SELECTED GRID COORDINATES C (IOPT=-1) GRID COORDINATES OF SELECTED EARTH COORDINATES C FOR GAUSSIAN CYLINDRICAL PROJECTIONS. C IF THE SELECTED COORDINATES ARE MORE THAN ONE GRIDPOINT C BEYOND THE THE EDGES OF THE GRID DOMAIN, THEN THE RELEVANT C OUTPUT ELEMENTS ARE SET TO FILL VALUES. C THE ACTUAL NUMBER OF VALID POINTS COMPUTED IS RETURNED TOO. C OPTIONALLY, THE VECTOR ROTATIONS AND THE MAP JACOBIANS C FOR THIS GRID MAY BE RETURNED AS WELL. C C PROGRAM HISTORY LOG: C 96-04-10 IREDELL C 97-10-20 IREDELL INCLUDE MAP OPTIONS C 1999-04-08 IREDELL USE SUBROUTINE SPLAT C 2001-06-18 IREDELL CORRECT AREA COMPUTATION C C USAGE: CALL GDSWZD04(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, C & LROT,CROT,SROT,LMAP,XLON,XLAT,YLON,YLAT,AREA) C C INPUT ARGUMENT LIST: C KGDS - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63 C IOPT - INTEGER OPTION FLAG C (+1 TO COMPUTE EARTH COORDS OF SELECTED GRID COORDS) C (-1 TO COMPUTE GRID COORDS OF SELECTED EARTH COORDS) C NPTS - INTEGER MAXIMUM NUMBER OF COORDINATES C FILL - REAL FILL VALUE TO SET INVALID OUTPUT DATA C (MUST BE IMPOSSIBLE VALUE; SUGGESTED VALUE: -9999.) C XPTS - REAL (NPTS) GRID X POINT COORDINATES IF IOPT>0 C YPTS - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT>0 C RLON - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT<0 C (ACCEPTABLE RANGE: -360. TO 360.) C RLAT - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT<0 C (ACCEPTABLE RANGE: -90. TO 90.) C LROT - INTEGER FLAG TO RETURN VECTOR ROTATIONS IF 1 C LMAP - INTEGER FLAG TO RETURN MAP JACOBIANS IF 1 C C OUTPUT ARGUMENT LIST: C XPTS - REAL (NPTS) GRID X POINT COORDINATES IF IOPT<0 C YPTS - REAL (NPTS) GRID Y POINT COORDINATES IF IOPT<0 C RLON - REAL (NPTS) EARTH LONGITUDES IN DEGREES E IF IOPT>0 C RLAT - REAL (NPTS) EARTH LATITUDES IN DEGREES N IF IOPT>0 C NRET - INTEGER NUMBER OF VALID POINTS COMPUTED C CROT - REAL (NPTS) CLOCKWISE VECTOR ROTATION COSINES IF LROT=1 C SROT - REAL (NPTS) CLOCKWISE VECTOR ROTATION SINES IF LROT=1 C (UGRID=CROT*UEARTH-SROT*VEARTH; C VGRID=SROT*UEARTH+CROT*VEARTH) C XLON - REAL (NPTS) DX/DLON IN 1/DEGREES IF LMAP=1 C XLAT - REAL (NPTS) DX/DLAT IN 1/DEGREES IF LMAP=1 C YLON - REAL (NPTS) DY/DLON IN 1/DEGREES IF LMAP=1 C YLAT - REAL (NPTS) DY/DLAT IN 1/DEGREES IF LMAP=1 C AREA - REAL (NPTS) AREA WEIGHTS IN M**2 IF LMAP=1 C C SUBPROGRAMS CALLED: C SPLAT COMPUTE LATITUDE FUNCTIONS C C ATTRIBUTES: C LANGUAGE: FORTRAN 77 C C$$$ INTEGER KGDS(200) REAL XPTS(NPTS),YPTS(NPTS),RLON(NPTS),RLAT(NPTS) REAL CROT(NPTS),SROT(NPTS) REAL XLON(NPTS),XLAT(NPTS),YLON(NPTS),YLAT(NPTS),AREA(NPTS) PARAMETER(RERTH=6.3712E6) PARAMETER(PI=3.14159265358979,DPR=180./PI) PARAMETER(JGMAX=2000) REAL ALAT(0:JGMAX+1),BLAT(0:JGMAX+1) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF(KGDS(1).EQ.004.AND.KGDS(10)*2.LE.JGMAX) THEN IM=KGDS(2) JM=KGDS(3) RLAT1=KGDS(4)*1.E-3 RLON1=KGDS(5)*1.E-3 RLAT2=KGDS(7)*1.E-3 RLON2=KGDS(8)*1.E-3 JG=KGDS(10)*2 ISCAN=MOD(KGDS(11)/128,2) JSCAN=MOD(KGDS(11)/64,2) NSCAN=MOD(KGDS(11)/32,2) HI=(-1.)**ISCAN JH=(-1)**JSCAN DLON=HI*(MOD(HI*(RLON2-RLON1)-1+3600,360.)+1)/(IM-1) CALL SPLAT(4,JG,ALAT(1),BLAT) DO JA=1,JG ALAT(JA)=DPR*ASIN(ALAT(JA)) ENDDO ALAT(0)=180.-ALAT(1) ALAT(JG+1)=-ALAT(0) BLAT(0)=-BLAT(1) BLAT(JG+1)=BLAT(0) J1=1 DOWHILE(J1.LT.JG.AND.RLAT1.LT.(ALAT(J1)+ALAT(J1+1))/2) J1=J1+1 ENDDO J2=J1+JH*(JM-1) XMIN=0 XMAX=IM+1 IF(IM.EQ.NINT(360/ABS(DLON))) XMAX=IM+2 YMIN=0.5 YMAX=JM+0.5 NRET=0 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C TRANSLATE GRID COORDINATES TO EARTH COORDINATES IF(IOPT.EQ.0.OR.IOPT.EQ.1) THEN DO N=1,NPTS IF(XPTS(N).GE.XMIN.AND.XPTS(N).LE.XMAX.AND. & YPTS(N).GE.YMIN.AND.YPTS(N).LE.YMAX) THEN RLON(N)=MOD(RLON1+DLON*(XPTS(N)-1)+3600,360.) J=YPTS(N) WB=YPTS(N)-J RLATA=ALAT(J1+JH*(J-1)) RLATB=ALAT(J1+JH*J) RLAT(N)=RLATA+WB*(RLATB-RLATA) NRET=NRET+1 IF(LROT.EQ.1) THEN CROT(N)=1 SROT(N)=0 ENDIF IF(LMAP.EQ.1) THEN XLON(N)=1/DLON XLAT(N)=0. YLON(N)=0. YLAT(N)=1/(RLATB-RLATA) WLATA=BLAT(J1+JH*(J-1)) WLATB=BLAT(J1+JH*J) WLAT=WLATA+WB*(WLATB-WLATA) AREA(N)=RERTH**2*WLAT*DLON/DPR ENDIF ELSE RLON(N)=FILL RLAT(N)=FILL ENDIF ENDDO C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C TRANSLATE EARTH COORDINATES TO GRID COORDINATES ELSEIF(IOPT.EQ.-1) THEN DO N=1,NPTS XPTS(N)=FILL YPTS(N)=FILL IF(ABS(RLON(N)).LE.360.AND.ABS(RLAT(N)).LE.90) THEN XPTS(N)=1+HI*MOD(HI*(RLON(N)-RLON1)+3600,360.)/DLON JA=MIN(INT((JG+1)/180.*(90-RLAT(N))),JG) IF(RLAT(N).GT.ALAT(JA)) JA=MAX(JA-2,0) IF(RLAT(N).LT.ALAT(JA+1)) JA=MIN(JA+2,JG) IF(RLAT(N).GT.ALAT(JA)) JA=JA-1 IF(RLAT(N).LT.ALAT(JA+1)) JA=JA+1 YPTSA=1+JH*(JA-J1) YPTSB=1+JH*(JA+1-J1) WB=(ALAT(JA)-RLAT(N))/(ALAT(JA)-ALAT(JA+1)) YPTS(N)=YPTSA+WB*(YPTSB-YPTSA) IF(XPTS(N).GE.XMIN.AND.XPTS(N).LE.XMAX.AND. & YPTS(N).GE.YMIN.AND.YPTS(N).LE.YMAX) THEN NRET=NRET+1 IF(LROT.EQ.1) THEN CROT(N)=1 SROT(N)=0 ENDIF IF(LMAP.EQ.1) THEN XLON(N)=1/DLON XLAT(N)=0. YLON(N)=0. YLAT(N)=JH/(ALAT(JA)-ALAT(JA+1)) WLATA=BLAT(JA) WLATB=BLAT(JA+1) WLAT=WLATA+WB*(WLATB-WLATA) AREA(N)=RERTH**2*WLAT*DLON/DPR ENDIF ELSE XPTS(N)=FILL YPTS(N)=FILL ENDIF ENDIF ENDDO ENDIF C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C PROJECTION UNRECOGNIZED ELSE IRET=-1 IF(IOPT.GE.0) THEN DO N=1,NPTS RLON(N)=FILL RLAT(N)=FILL ENDDO ENDIF IF(IOPT.LE.0) THEN DO N=1,NPTS XPTS(N)=FILL YPTS(N)=FILL ENDDO ENDIF ENDIF C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - END