C----------------------------------------------------------------------- C PROGRAM BILL_PREC c--------------------------------------------------------------------- c This program performs areal average of Stage IV 4km observed precip c data onto model's own grid. It reads in the degribbed data c having name degrib1d.dat that is created by bill_org.f. c--------------------------------------------------------------------- c SUBROUTINE IPOLATES(IP,IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI, c & NO,RLAT,RLON,IBO,LO,GO,IRET) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: IPOLATES IREDELL'S POLATE FOR SCALAR FIELDS C PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10 C C ABSTRACT: THIS SUBPROGRAM INTERPOLATES SCALAR FIELDS C FROM ANY GRID TO ANY GRID (JOE IRWIN'S DREAM). C ONLY HORIZONTAL INTERPOLATION IS PERFORMED. C THE FOLLOWING INTERPOLATION METHODS ARE POSSIBLE: C (IP=0) BILINEAR C (IP=1) BICUBIC C (IP=2) NEIGHBOR C (IP=3) BUDGET C (IP=4) SPECTRAL C (IP=6) NEIGHBOR-BUDGET C SOME OF THESE METHODS HAVE INTERPOLATION OPTIONS AND/OR C RESTRICTIONS ON THE INPUT OR OUTPUT GRIDS, BOTH OF WHICH C ARE DOCUMENTED MORE FULLY IN THEIR RESPECTIVE SUBPROGRAMS. C THE GRIDS ARE DEFINED BY THEIR GRID DESCRIPTION SECTIONS C (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63). C THE CURRENT CODE RECOGNIZES THE FOLLOWING PROJECTIONS: C (KGDS(1)=000) EQUIDISTANT CYLINDRICAL C (KGDS(1)=001) MERCATOR CYLINDRICAL C (KGDS(1)=003) LAMBERT CONFORMAL CONICAL C (KGDS(1)=004) GAUSSIAN CYLINDRICAL C (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL C (KGDS(1)=201) ROTATED EQUIDISTANT CYLINDRICAL C (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL C WHERE KGDS COULD BE EITHER INPUT KGDSI OR OUTPUT KGDSO. C AS AN ADDED BONUS THE NUMBER OF OUTPUT GRID POINTS C AND THEIR LATITUDES AND LONGITUDES ARE ALSO RETURNED. C ON THE OTHER HAND, THE OUTPUT CAN BE A SET OF STATION POINTS C IF KGDSO(1)<0, IN WHICH CASE THE NUMBER OF POINTS C AND THEIR LATITUDES AND LONGITUDES MUST BE INPUT. C INPUT BITMAPS WILL BE INTERPOLATED TO OUTPUT BITMAPS. C OUTPUT BITMAPS WILL ALSO BE CREATED WHEN THE OUTPUT GRID C EXTENDS OUTSIDE OF THE DOMAIN OF THE INPUT GRID. C THE OUTPUT FIELD IS SET TO 0 WHERE THE OUTPUT BITMAP IS OFF. C C PROGRAM HISTORY LOG: C 96-04-10 IREDELL C C USAGE: CALL IPOLATES(IP,IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI, C & NO,RLAT,RLON,IBO,LO,GO,IRET) C C INPUT ARGUMENT LIST: C IP - INTEGER INTERPOLATION METHOD C (IP=0 FOR BILINEAR; C IP=1 FOR BICUBIC; C IP=2 FOR NEIGHBOR; C IP=3 FOR BUDGET; C IP=4 FOR SPECTRAL; C IP=6 FOR NEIGHBOR-BUDGET) C IPOPT - INTEGER (20) INTERPOLATION OPTIONS C (IP=0: (NO OPTIONS) C IP=1: CONSTRAINT OPTION C IP=2: (NO OPTIONS) C IP=3: NUMBER IN RADIUS, RADIUS WEIGHTS ... C IP=4: SPECTRAL SHAPE, SPECTRAL TRUNCATION C IP=6: NUMBER IN RADIUS, RADIUS WEIGHTS ...) C KGDSI - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63 C KGDSO - INTEGER (200) OUTPUT GDS PARAMETERS C MI - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1 C OR DIMENSION OF INPUT GRID FIELDS IF KM=1 C MO - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1 C OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1 C KM - INTEGER NUMBER OF FIELDS TO INTERPOLATE C IBI - INTEGER (KM) INPUT BITMAP FLAGS C LI - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF RESPECTIVE IBI(K)=1) C GI - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE C NO - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)<0) C RLAT - REAL (NO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)<0) C RLON - REAL (NO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)<0) C C OUTPUT ARGUMENT LIST: C NO - INTEGER NUMBER OF OUTPUT POINTS (ONLY IF KGDSO(1)>=0) C RLAT - REAL (MO) OUTPUT LATITUDES IN DEGREES (IF KGDSO(1)>=0) C RLON - REAL (MO) OUTPUT LONGITUDES IN DEGREES (IF KGDSO(1)>=0) C IBO - INTEGER (KM) OUTPUT BITMAP FLAGS C LO - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT) C GO - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED C IRET - INTEGER RETURN CODE C 0 SUCCESSFUL INTERPOLATION C 1 UNRECOGNIZED INTERPOLATION METHOD C 2 UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP C 3 UNRECOGNIZED OUTPUT GRID C 1X INVALID BICUBIC METHOD PARAMETERS C 3X INVALID BUDGET METHOD PARAMETERS C 4X INVALID SPECTRAL METHOD PARAMETERS C C SUBPROGRAMS CALLED: C POLATES0 INTERPOLATE SCALAR FIELDS (BILINEAR) C POLATES1 INTERPOLATE SCALAR FIELDS (BICUBIC) C POLATES2 INTERPOLATE SCALAR FIELDS (NEIGHBOR) C POLATES3 INTERPOLATE SCALAR FIELDS (BUDGET) C POLATES4 INTERPOLATE SCALAR FIELDS (SPECTRAL) C C REMARKS: EXAMPLES DEMONSTRATING RELATIVE CPU COSTS. C THIS EXAMPLE IS INTERPOLATING 12 LEVELS OF TEMPERATURES C FROM THE 360 X 181 GLOBAL GRID (NCEP GRID 3) C TO THE 93 X 68 HAWAIIAN MERCATOR GRID (NCEP GRID 204). C THE EXAMPLE TIMES ARE FOR THE C90. AS A REFERENCE, THE CP TIME C FOR UNPACKING THE GLOBAL 12 TEMPERATURE FIELDS IS 0.04 SECONDS. C C METHOD IP IPOPT CP SECONDS C -------- -- ------------- ---------- C BILINEAR 0 0.03 C BICUBIC 1 0 0.07 C BICUBIC 1 1 0.07 C NEIGHBOR 2 0.01 C BUDGET 3 -1,-1 0.48 C SPECTRAL 4 0,40 0.22 C SPECTRAL 4 1,40 0.24 C SPECTRAL 4 0,-1 0.42 C N-BUDGET 6 -1,-1 0.15 C C THE SPECTRAL INTERPOLATION IS FAST FOR THE MERCATOR GRID. C HOWEVER, FOR SOME GRIDS THE SPECTRAL INTERPOLATION IS SLOW. C THE FOLLOWING EXAMPLE IS INTERPOLATING 12 LEVELS OF TEMPERATURES C FROM THE 360 X 181 GLOBAL GRID (NCEP GRID 3) C TO THE 93 X 65 CONUS LAMBERT CONFORMAL GRID (NCEP GRID 211). C C METHOD IP IPOPT CP SECONDS C -------- -- ------------- ---------- C BILINEAR 0 0.03 C BICUBIC 1 0 0.07 C BICUBIC 1 1 0.07 C NEIGHBOR 2 0.01 C BUDGET 3 -1,-1 0.51 C SPECTRAL 4 0,40 3.94 C SPECTRAL 4 1,40 5.02 C SPECTRAL 4 0,-1 11.36 C N-BUDGET 6 -1,-1 0.18 C C ATTRIBUTES: C LANGUAGE: FORTRAN 77 C C$$$ c ****mi=wrf grib specs, mo=eta grid specs****** PARAMETER(KM=1,MI=1121*881,MO=502*543) INTEGER IPOPT(20) INTEGER KGDSI(200),KGDSO(200) INTEGER IBI(KM),IBO(KM) LOGICAL*1 LI(MI,KM),LO(MO,KM) REAL GI(MI,KM),GO(MO,KM) REAL RLAT(MO),RLON(MO) cgal c these parameters are valid for 4 km stage 4 data c*********** find these for wrf grid ******************** IP=6 DO I=1,20 IPOPT(I)=-1 END DO KGDSI(1)=5 KGDSI(2)=1121 KGDSI(3)=881 KGDSI(4)=23117 KGDSI(5)=240977 KGDSI(6)=8 KGDSI(7)=255000 KGDSI(8)=4763 KGDSI(9)=4763 KGDSI(10)=0 KGDSI(11)=64 DO I=12,19 KGDSI(I)=0 END DO KGDSI(20)=255 DO I=21,200 KGDSI(I)=0 END DO c projection type (2D Eta) KGDSO(1)=003 c (x) i-pts KGDSO(2)=502 c (y) j-pts KGDSO(3)=543 c lat of 1st pt KGDSO(4)=28760 c lon of 1st pt KGDSO(5)=255010 c KGDSO(5)=253500 c flag KGDSO(6)=8 c center lon KGDSO(7)=263000 c center lat KGDSO(8)=4000 c d-lam*1000 KGDSO(9)=4000 c d-lon*1000 KGDSO(10)=0 c flag KGDSO(11)=64 KGDSO(12)=32000 KGDSO(13)=45000 DO I=14,19 KGDSO(I)=0 END DO KGDSO(20)=255 DO I=21,200 KGDSO(I)=0 END DO DO I=1,KM IBI(I)=1 END DO OPEN(9,file='degrib1d.dat',form='unformatted') READ(9) GI CLOSE(unit=9) pmax=0. do k=1,MI pmax=amax1(pmax,gi(k,1)) LI(k,1)=.true. if(GI(k,1).gt.1.0E10) LI(k,1)=.false. end do write(6,*) 'pmax is ',pmax C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C BILINEAR INTERPOLATION IF(IP.EQ.0) THEN cgal CALL POLATES0(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI, cgal & NO,RLAT,RLON,IBO,LO,GO,IRET) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C BICUBIC INTERPOLATION ELSEIF(IP.EQ.1) THEN cgal CALL POLATES1(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI, cgal & NO,RLAT,RLON,IBO,LO,GO,IRET) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C NEIGHBOR INTERPOLATION ELSEIF(IP.EQ.2) THEN cgal CALL POLATES2(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI, cgal & NO,RLAT,RLON,IBO,LO,GO,IRET) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C BUDGET INTERPOLATION ELSEIF(IP.EQ.3) THEN cgal CALL POLATES3(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI, cgal & NO,RLAT,RLON,IBO,LO,GO,IRET) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C SPECTRAL INTERPOLATION ELSEIF(IP.EQ.4) THEN cgal CALL POLATES4(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI, cgal & NO,RLAT,RLON,IBO,LO,GO,IRET) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C NEIGHBOR-BUDGET INTERPOLATION ELSEIF(IP.EQ.6) THEN CALL POLATES6(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI, & NO,RLAT,RLON,IBO,LO,GO,IRET) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C UNRECOGNIZED INTERPOLATION METHOD ELSE IF(KGDSO(1).GE.0) NO=0 CMIC$ DO ALL AUTOSCOPE DO K=1,KM IBO(K)=1 DO N=1,NO LO(N,K)=.FALSE. GO(N,K)=0. ENDDO ENDDO IRET=1 ENDIF c open(10,file='obsprec.dat',form='unformatted') write(10) go close(unit=10) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - END C----------------------------------------------------------------------- SUBROUTINE POLATES6(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI, & NO,RLAT,RLON,IBO,LO,GO,IRET) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: POLATES6 INTERPOLATE SCALAR FIELDS (BUDGET) C PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10 C C ABSTRACT: THIS SUBPROGRAM PERFORMS BUDGET INTERPOLATION C FROM ANY GRID TO ANY GRID FOR SCALAR FIELDS. C IT REQUIRES A GRID FOR THE OUTPUT FIELDS (KGDSO(1)>=0). C THE ALGORITHM SIMPLY COMPUTES (WEIGHTED) AVERAGES C OF NEIGHBOR POINTS ARRANGED IN A SQUARE BOX C CENTERED AROUND EACH OUTPUT GRID POINT AND STRETCHING C NEARLY HALFWAY TO EACH OF THE NEIGHBORING GRID POINTS. C OPTIONS ALLOW CHOICES OF NUMBER OF POINTS IN EACH RADIUS C FROM THE CENTER POINT (IPOPT(1)) WHICH DEFAULTS TO 2 C (IF IPOPT(1)=-1) MEANING THAT 25 POINTS WILL BE AVERAGED; C FURTHER OPTIONS ARE THE RESPECTIVE WEIGHTS FOR THE RADIUS C POINTS STARTING AT THE CENTER POINT (IPOPT(2:2+IPOPT(1)) C WHICH DEFAULTS TO ALL 1 (IF IPOPT(2)=-1.). C ONLY HORIZONTAL INTERPOLATION IS PERFORMED. C THE GRIDS ARE DEFINED BY THEIR GRID DESCRIPTION SECTIONS C (PASSED IN INTEGER FORM AS DECODED BY SUBPROGRAM W3FI63). C THE CURRENT CODE RECOGNIZES THE FOLLOWING PROJECTIONS: C (KGDS(1)=000) EQUIDISTANT CYLINDRICAL C (KGDS(1)=001) MERCATOR CYLINDRICAL C (KGDS(1)=003) LAMBERT CONFORMAL CONICAL C (KGDS(1)=004) GAUSSIAN CYLINDRICAL (SPECTRAL NATIVE) C (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL C (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL (ETA NATIVE) C WHERE KGDS COULD BE EITHER INPUT KGDSI OR OUTPUT KGDSO. C AS AN ADDED BONUS THE NUMBER OF OUTPUT GRID POINTS C AND THEIR LATITUDES AND LONGITUDES ARE ALSO RETURNED. C INPUT BITMAPS WILL BE INTERPOLATED TO OUTPUT BITMAPS. C OUTPUT BITMAPS WILL ALSO BE CREATED WHEN THE OUTPUT GRID C EXTENDS OUTSIDE OF THE DOMAIN OF THE INPUT GRID. C THE OUTPUT FIELD IS SET TO 0 WHERE THE OUTPUT BITMAP IS OFF. C C PROGRAM HISTORY LOG: C 96-04-10 IREDELL C 96-10-04 IREDELL NEIGHBOR POINTS NOT BILINEAR INTERPOLATION C 1999-04-08 IREDELL SPLIT IJKGDS INTO TWO PIECES C C USAGE: CALL POLATES6(IPOPT,KGDSI,KGDSO,MI,MO,KM,IBI,LI,GI, C & NO,RLAT,RLON,IBO,LO,GO,IRET) C C INPUT ARGUMENT LIST: C IPOPT - INTEGER (20) INTERPOLATION OPTIONS C IPOPT(1) IS NUMBER OF RADIUS POINTS C (DEFAULTS TO 2 IF IPOPT(1)=-1); C IPOPT(2:2+IPOPT(1)) ARE RESPECTIVE WEIGHTS C (DEFAULTS TO ALL 1 IF IPOPT(2)=-1). C KGDSI - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63 C KGDSO - INTEGER (200) OUTPUT GDS PARAMETERS C MI - INTEGER SKIP NUMBER BETWEEN INPUT GRID FIELDS IF KM>1 C OR DIMENSION OF INPUT GRID FIELDS IF KM=1 C MO - INTEGER SKIP NUMBER BETWEEN OUTPUT GRID FIELDS IF KM>1 C OR DIMENSION OF OUTPUT GRID FIELDS IF KM=1 C KM - INTEGER NUMBER OF FIELDS TO INTERPOLATE C IBI - INTEGER (KM) INPUT BITMAP FLAGS C LI - LOGICAL*1 (MI,KM) INPUT BITMAPS (IF SOME IBI(K)=1) C GI - REAL (MI,KM) INPUT FIELDS TO INTERPOLATE C C OUTPUT ARGUMENT LIST: C NO - INTEGER NUMBER OF OUTPUT POINTS C RLAT - REAL (MO) OUTPUT LATITUDES IN DEGREES C RLON - REAL (MO) OUTPUT LONGITUDES IN DEGREES C IBO - INTEGER (KM) OUTPUT BITMAP FLAGS C LO - LOGICAL*1 (MO,KM) OUTPUT BITMAPS (ALWAYS OUTPUT) C GO - REAL (MO,KM) OUTPUT FIELDS INTERPOLATED C IRET - INTEGER RETURN CODE C 0 SUCCESSFUL INTERPOLATION C 2 UNRECOGNIZED INPUT GRID OR NO GRID OVERLAP C 3 UNRECOGNIZED OUTPUT GRID C 31 INVALID UNDEFINED OUTPUT GRID C 32 INVALID BUDGET METHOD PARAMETERS C C SUBPROGRAMS CALLED: C GDSWIZ GRID DESCRIPTION SECTION WIZARD C IJKGDS0 SET UP PARAMETERS FOR IJKGDS1 C (IJKGDS1) RETURN FIELD POSITION FOR A GIVEN GRID POINT C POLFIXS MAKE MULTIPLE POLE SCALAR VALUES CONSISTENT C C ATTRIBUTES: C LANGUAGE: FORTRAN 77 C C$$$ CFPP$ EXPAND(IJKGDS1) INTEGER IPOPT(20) INTEGER KGDSI(200),KGDSO(200) INTEGER IBI(KM),IBO(KM) LOGICAL*1 LI(MI,KM),LO(MO,KM) REAL GI(MI,KM),GO(MO,KM) REAL RLAT(MO),RLON(MO) REAL XPTS(MO),YPTS(MO) REAL XPTB(MO),YPTB(MO),RLOB(MO),RLAB(MO) INTEGER N11(MO) REAL WO(MO,KM) INTEGER IJKGDSA(20) PARAMETER(FILL=-9999.) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C COMPUTE NUMBER OF OUTPUT POINTS AND THEIR LATITUDES AND LONGITUDES. IRET=0 IF(KGDSO(1).GE.0) THEN CALL GDSWIZ(KGDSO, 0,MO,FILL,XPTS,YPTS,RLON,RLAT,NO,0,DUM,DUM) IF(NO.EQ.0) IRET=3 ELSE IRET=31 ENDIF write(6,*) 'in polates lat and lon ',rlon(30),rlat(30),rlon(300) 2,rlat(300),rlon(3000),rlat(3000) pmax=0. do i=1,mi pmax=amax1(pmax,gi(i,1)) end do do i=1,mi if(li(i,1).eqv..true.) then write(6,*) 'true bitmap' go to 193 endif end do 193 continue write(6,*) 'in polates gi max ',pmax,no C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C SET PARAMETERS NB1=IPOPT(1) IF(NB1.EQ.-1) NB1=2 IF(IRET.EQ.0.AND.NB1.LT.0.) IRET=32 IF(IRET.EQ.0.AND.NB1.GE.20.AND.IPOPT(2).NE.-1) IRET=32 IF(IRET.EQ.0) THEN NB2=2*NB1+1 NB3=NB2*NB2 NB4=NB3 IF(IPOPT(2).NE.-1) THEN NB4=IPOPT(2) DO IB=1,NB1 NB4=NB4+8*IB*IPOPT(2+IB) ENDDO ENDIF ELSE NB2=0 NB3=0 NB4=0 ENDIF CMIC$ DO ALL AUTOSCOPE DO K=1,KM DO N=1,NO GO(N,K)=0. WO(N,K)=0. ENDDO ENDDO C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C LOOP OVER SAMPLE POINTS IN OUTPUT GRID BOX CALL IJKGDS0(KGDSI,IJKGDSA) write(6,*) 'nb3 is ',nb3,nb2,nb1 DO NB=1,NB3 C LOCATE INPUT POINTS AND COMPUTE THEIR WEIGHTS JB=(NB-1)/NB2-NB1 IB=NB-(JB+NB1)*NB2-NB1-1 LB=MAX(ABS(IB),ABS(JB)) WB=1 IF(IPOPT(2).NE.-1) WB=IPOPT(2+LB) IF(WB.NE.0) THEN DO N=1,NO XPTB(N)=XPTS(N)+IB/REAL(NB2) YPTB(N)=YPTS(N)+JB/REAL(NB2) ENDDO CALL GDSWIZ(KGDSO, 1,NO,FILL,XPTB,YPTB,RLOB,RLAB,NV,0,DUM,DUM) CALL GDSWIZ(KGDSI,-1,NO,FILL,XPTB,YPTB,RLOB,RLAB,NV,0,DUM,DUM) IF(IRET.EQ.0.AND.NV.EQ.0.AND.LB.EQ.0) IRET=2 DO N=1,NO XI=XPTB(N) YI=YPTB(N) IF(XI.NE.FILL.AND.YI.NE.FILL) THEN I1=NINT(XI) J1=NINT(YI) N11(N)=IJKGDS1(I1,J1,IJKGDSA) ELSE N11(N)=0 ENDIF ENDDO C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C INTERPOLATE WITH OR WITHOUT BITMAPS CMIC$ DO ALL AUTOSCOPE DO K=1,KM DO N=1,NO IF(N11(N).GT.0) THEN IF(IBI(K).EQ.0.OR.LI(N11(N),K)) THEN GO(N,K)=GO(N,K)+WB*GI(N11(N),K) WO(N,K)=WO(N,K)+WB ENDIF ENDIF ENDDO ENDDO ENDIF ENDDO C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C COMPUTE OUTPUT BITMAPS AND FIELDS CMIC$ DO ALL AUTOSCOPE DO K=1,KM IBO(K)=IBI(K) DO N=1,NO LO(N,K)=WO(N,K).GE.0.5*NB4 c if(lo(n,k)) then c write(6,*) 'true points at ',n c endif IF(LO(N,K)) THEN GO(N,K)=GO(N,K)/WO(N,K) ELSE IBO(K)=1 GO(N,K)=0. ENDIF ENDDO ENDDO ppmax=0. do j=1,mo ppmax=amax1(ppmax,go(j,1)) end do write(6,*) 'ppmax in ipolates is ',ppmax,iret IF(KGDSO(1).EQ.0) CALL POLFIXS(NO,MO,KM,RLAT,RLON,IBO,LO,GO) ppmax=0. do j=1,mo ppmax=amax1(ppmax,go(j,1)) end do write(6,*) 'ppmax in ipolates is ',ppmax,no C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - END C----------------------------------------------------------------------- SUBROUTINE GDSWIZ(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: GDSWIZ GRID DESCRIPTION SECTION WIZARD 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= 0) GRID AND EARTH COORDINATES OF ALL GRID POINTS C (IOPT=+1) EARTH COORDINATES OF SELECTED GRID COORDINATES C (IOPT=-1) GRID COORDINATES OF SELECTED EARTH COORDINATES C THE CURRENT CODE RECOGNIZES THE FOLLOWING PROJECTIONS: C (KGDS(1)=000) EQUIDISTANT CYLINDRICAL C (KGDS(1)=001) MERCATOR CYLINDRICAL C (KGDS(1)=003) LAMBERT CONFORMAL CONICAL C (KGDS(1)=004) GAUSSIAN CYLINDRICAL C (KGDS(1)=005) POLAR STEREOGRAPHIC AZIMUTHAL C (KGDS(1)=201) STAGGERED ROTATED EQUIDISTANT CYLINDRICAL C (KGDS(1)=202) ROTATED EQUIDISTANT CYLINDRICAL C (KGDS(1)=203) STAGGERED ROTATED EQUIDISTANT CYLINDRICAL 2-D 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. ALSO IF IOPT=0, C IF THE NUMBER OF GRID POINTS EXCEEDS THE NUMBER ALLOTTED, C THEN ALL THE OUTPUT ELEMENTS ARE SET TO FILL VALUES. C THE ACTUAL NUMBER OF VALID POINTS COMPUTED IS RETURNED TOO. C C PROGRAM HISTORY LOG: C 96-04-10 IREDELL C 98-08-20 BALDWIN ADD TYPE 203 STAGGERED 2-D ETA GRIDS C C USAGE: CALL GDSWIZ(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, C & LROT,CROT,SROT) C C INPUT ARGUMENT LIST: C KGDS - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63 C IOPT - INTEGER OPTION FLAG C ( 0 TO COMPUTE EARTH COORDS OF ALL THE GRID POINTS) 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 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 (-1 IF PROJECTION UNRECOGNIZED) 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 C SUBPROGRAMS CALLED: C GDSWIZ00 GDS WIZARD FOR EQUIDISTANT CYLINDRICAL C GDSWIZ01 GDS WIZARD FOR MERCATOR CYLINDRICAL C GDSWIZ03 GDS WIZARD FOR LAMBERT CONFORMAL CONICAL C GDSWIZ04 GDS WIZARD FOR GAUSSIAN CYLINDRICAL C GDSWIZ05 GDS WIZARD FOR POLAR STEREOGRAPHIC AZIMUTHAL C GDSWIZC9 GDS WIZARD FOR ROTATED EQUIDISTANT CYLINDRICAL C GDSWIZCA GDS WIZARD FOR ROTATED EQUIDISTANT CYLINDRICAL C GDSWIZCB GDS WIZARD FOR ROTATED EQUIDISTANT CYLINDRICAL 2-D 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) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C COMPUTE GRID COORDINATES FOR ALL GRID POINTS IF(IOPT.EQ.0) THEN IF(KGDS(1).EQ.201) THEN IM=KGDS(7)*2-1 JM=KGDS(8) KSCAN=MOD(KGDS(11)/256,2) IF(KSCAN.EQ.0) THEN IS1=(JM+1)/2 NM=(IM/2+1)*JM-JM/2 ELSE IS1=JM/2 NM=IM/2*JM+JM/2 ENDIF ELSEIF(KGDS(1).EQ.202) THEN IM=KGDS(7) JM=KGDS(8) NM=IM*JM ELSEIF(KGDS(1).EQ.203) THEN IM=KGDS(2) JM=KGDS(3) NM=IM*JM KSCAN=MOD(KGDS(11)/256,2) IF(KSCAN.EQ.0) THEN IS1=(JM+1)/2 ELSE IS1=JM/2 ENDIF ELSE IM=KGDS(2) JM=KGDS(3) NM=IM*JM ENDIF NSCAN=MOD(KGDS(11)/32,2) IF(NM.LE.NPTS) THEN IF(KGDS(1).EQ.201) THEN DO N=1,NM NN=2*N-1+KSCAN IF(NSCAN.EQ.0) THEN J=(NN-1)/IM+1 I=NN-IM*(J-1) ELSE I=(NN-1)/JM+1 J=NN-JM*(I-1) ENDIF XPTS(N)=IS1+(I-(J-KSCAN))/2 YPTS(N)=(I+(J-KSCAN))/2 ENDDO ELSEIF(KGDS(1).EQ.203) THEN DO N=1,NM IF(NSCAN.EQ.0) THEN J=(N-1)/IM+1 I=(N-IM*(J-1))*2-MOD(J+KSCAN,2) ELSE I=(N-1)/JM+1 J=(N-JM*(I-1))*2-MOD(I+KSCAN,2) ENDIF XPTS(N)=IS1+(I-(J-KSCAN))/2 YPTS(N)=(I+(J-KSCAN))/2 ENDDO ELSE DO N=1,NM IF(NSCAN.EQ.0) THEN J=(N-1)/IM+1 I=N-IM*(J-1) ELSE I=(N-1)/JM+1 J=N-JM*(I-1) ENDIF XPTS(N)=I YPTS(N)=J ENDDO ENDIF DO N=NM+1,NPTS XPTS(N)=FILL YPTS(N)=FILL ENDDO ELSE DO N=1,NPTS XPTS(N)=FILL YPTS(N)=FILL ENDDO ENDIF IOPF=1 ELSE IOPF=IOPT ENDIF C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C EQUIDISTANT CYLINDRICAL IF(KGDS(1).EQ.000) THEN CALL GDSWIZ00(KGDS,IOPF,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C MERCATOR CYLINDRICAL ELSEIF(KGDS(1).EQ.001) THEN CALL GDSWIZ01(KGDS,IOPF,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C LAMBERT CONFORMAL CONICAL ELSEIF(KGDS(1).EQ.003) THEN CALL GDSWIZ03(KGDS,IOPF,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C GAUSSIAN CYLINDRICAL c ELSEIF(KGDS(1).EQ.004) THEN c CALL GDSWIZ04(KGDS,IOPF,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, c & LROT,CROT,SROT) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C POLAR STEREOGRAPHIC AZIMUTHAL ELSEIF(KGDS(1).EQ.005) THEN CALL GDSWIZ05(KGDS,IOPF,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C STAGGERED ROTATED EQUIDISTANT CYLINDRICAL ELSEIF(KGDS(1).EQ.201) THEN CALL GDSWIZC9(KGDS,IOPF,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C ROTATED EQUIDISTANT CYLINDRICAL ELSEIF(KGDS(1).EQ.202) THEN CALL GDSWIZCA(KGDS,IOPF,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C STAGGERED ROTATED EQUIDISTANT CYLINDRICAL ELSEIF(KGDS(1).EQ.203) THEN CALL GDSWIZCB(KGDS,IOPF,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) 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 C----------------------------------------------------------------------- SUBROUTINE IJKGDS0(KGDS,IJKGDSA) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: IJKGDS0 SET UP PARAMETERS FOR FUNCTION IJKGDS1 C PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10 C C ABSTRACT: THIS SUBPROGRAM DECODES THE GRIB GRID DESCRIPTION SECTION C AND RETURNS A NAVIGATION PARAMETER ARRAY TO ALLOW FUNCTION C IJKGDS1 TO DECODE THE FIELD POSITION FOR A GIVEN GRID POINT. C C PROGRAM HISTORY LOG: C 96-04-10 IREDELL C 97-03-11 IREDELL ALLOWED HEMISPHERIC GRIDS TO WRAP OVER ONE POLE C 98-07-13 BALDWIN ADD 2D STAGGERED ETA GRID INDEXING (203) C 1999-04-08 IREDELL SPLIT IJKGDS INTO TWO C C USAGE: CALL IJKGDS0(KGDS,IJKGDSA) C C INPUT ARGUMENT LIST: C KGDS - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63 C C OUTPUT ARGUMENT LIST: C IJKGDSA - INTEGER (20) NAVIGATION PARAMETER ARRAY C IJKGDSA(1) IS NUMBER OF X POINTS C IJKGDSA(2) IS NUMBER OF Y POINTS C IJKGDSA(3) IS X WRAPAROUND INCREMENT C (0 IF NO WRAPAROUND) C IJKGDSA(4) IS Y WRAPAROUND LOWER PIVOT POINT C (0 IF NO WRAPAROUND) C IJKGDSA(5) IS Y WRAPAROUND UPPER PIVOT POINT C (0 IF NO WRAPAROUND) C IJKGDSA(6) IS SCANNING MODE C (0 IF X FIRST THEN Y; 1 IF Y FIRST THEN X; C 2 IF STAGGERED DIAGONAL LIKE PROJECTION 201; C 3 IF STAGGERED DIAGONAL LIKE PROJECTION 203) C IJKGDSA(7) IS MASS/WIND FLAG FOR STAGGERED DIAGONAL C (0 IF MASS; 1 IF WIND) C IJKGDSA(8:20) ARE UNUSED AT THE MOMENT C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C C$$$ INTEGER KGDS(200) INTEGER IJKGDSA(20) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C SET USUAL VALUES IM=KGDS(2) JM=KGDS(3) IWRAP=0 JWRAP1=0 JWRAP2=0 NSCAN=MOD(KGDS(11)/32,2) KSCAN=0 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C SET EXCEPTIONAL VALUES FOR PROJECTION 000 IF(KGDS(1).EQ.0) THEN RLON1=KGDS(5)*1.E-3 RLON2=KGDS(8)*1.E-3 ISCAN=MOD(KGDS(11)/128,2) IF(ISCAN.EQ.0) THEN DLON=(MOD(RLON2-RLON1-1+3600,360.)+1)/(IM-1) ELSE DLON=-(MOD(RLON1-RLON2-1+3600,360.)+1)/(IM-1) ENDIF IWRAP=NINT(360/ABS(DLON)) IF(IM.LT.IWRAP) IWRAP=0 IF(IWRAP.GT.0.AND.MOD(IWRAP,2).EQ.0) THEN RLAT1=KGDS(4)*1.E-3 RLAT2=KGDS(7)*1.E-3 DLAT=ABS(RLAT2-RLAT1)/(JM-1) IF(ABS(RLAT1).GT.90-0.25*DLAT) THEN JWRAP1=2 ELSEIF(ABS(RLAT1).GT.90-0.75*DLAT) THEN JWRAP1=1 ENDIF IF(ABS(RLAT2).GT.90-0.25*DLAT) THEN JWRAP2=2*JM ELSEIF(ABS(RLAT2).GT.90-0.75*DLAT) THEN JWRAP2=2*JM+1 ENDIF ENDIF C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C SET EXCEPTIONAL VALUES FOR PROJECTION 001 ELSEIF(KGDS(1).EQ.1) THEN RLON1=KGDS(5)*1.E-3 RLON2=KGDS(8)*1.E-3 ISCAN=MOD(KGDS(11)/128,2) IF(ISCAN.EQ.0) THEN DLON=(MOD(RLON2-RLON1-1+3600,360.)+1)/(IM-1) ELSE DLON=-(MOD(RLON1-RLON2-1+3600,360.)+1)/(IM-1) ENDIF IWRAP=NINT(360/ABS(DLON)) IF(IM.LT.IWRAP) IWRAP=0 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C SET EXCEPTIONAL VALUES FOR PROJECTION 004 ELSEIF(KGDS(1).EQ.4) THEN RLON1=KGDS(5)*1.E-3 RLON2=KGDS(8)*1.E-3 ISCAN=MOD(KGDS(11)/128,2) IF(ISCAN.EQ.0) THEN DLON=(MOD(RLON2-RLON1-1+3600,360.)+1)/(IM-1) ELSE DLON=-(MOD(RLON1-RLON2-1+3600,360.)+1)/(IM-1) ENDIF IWRAP=NINT(360/ABS(DLON)) IF(IM.LT.IWRAP) IWRAP=0 IF(IWRAP.GT.0.AND.MOD(IWRAP,2).EQ.0) THEN JG=KGDS(10)*2 IF(JM.EQ.JG) THEN JWRAP1=1 JWRAP2=2*JM+1 ENDIF ENDIF C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C SET EXCEPTIONAL VALUES FOR PROJECTION 201 ELSEIF(KGDS(1).EQ.201) THEN IM=KGDS(7) JM=KGDS(8) NSCAN=2 KSCAN=MOD(KGDS(11)/256,2) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C SET EXCEPTIONAL VALUES FOR PROJECTION 202 ELSEIF(KGDS(1).EQ.202) THEN IM=KGDS(7) JM=KGDS(8) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C SET EXCEPTIONAL VALUES FOR PROJECTION 203 ELSEIF(KGDS(1).EQ.203) THEN NSCAN=3 KSCAN=MOD(KGDS(11)/256,2) ENDIF C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C FILL NAVIGATION PARAMETER ARRAY IJKGDSA(1)=IM IJKGDSA(2)=JM IJKGDSA(3)=IWRAP IJKGDSA(4)=JWRAP1 IJKGDSA(5)=JWRAP2 IJKGDSA(6)=NSCAN IJKGDSA(7)=KSCAN C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - END C----------------------------------------------------------------------- SUBROUTINE POLFIXS(NM,NX,KM,RLAT,RLON,IB,LO,GO) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: POLFIXS MAKE MULTIPLE POLE SCALAR VALUES CONSISTENT C PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10 C C ABSTRACT: THIS SUBPROGRAM AVERAGES MULTIPLE POLE SCALAR VALUES C ON A LATITUDE/LONGITUDE GRID. BITMAPS MAY BE AVERAGED TOO. C C PROGRAM HISTORY LOG: C 96-04-10 IREDELL C C USAGE: CALL POLFIXS(NM,NX,KM,RLAT,RLON,IB,LO,GO) C C INPUT ARGUMENT LIST: C NO - INTEGER NUMBER OF GRID POINTS C NX - INTEGER LEADING DIMENSION OF FIELDS C KM - INTEGER NUMBER OF FIELDS C RLAT - REAL (NO) LATITUDES IN DEGREES C RLON - REAL (NO) LONGITUDES IN DEGREES C IB - INTEGER (KM) BITMAP FLAGS C LO - LOGICAL*1 (NX,KM) BITMAPS (IF SOME IB(K)=1) C GO - REAL (NX,KM) FIELDS C C OUTPUT ARGUMENT LIST: C LO - LOGICAL*1 (NX,KM) BITMAPS (IF SOME IB(K)=1) C GO - REAL (NX,KM) FIELDS C C ATTRIBUTES: C LANGUAGE: FORTRAN 77 C C$$$ PARAMETER(RLATNP=89.9995,RLATSP=-RLATNP) REAL RLAT(NM),RLON(NM) INTEGER IB(KM) REAL GO(NX,KM) LOGICAL*1 LO(NX,KM) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - CMIC$ DO ALL AUTOSCOPE DO K=1,KM WNP=0. GNP=0. TNP=0. WSP=0. GSP=0. TSP=0. C AVERAGE MULTIPLE POLE VALUES DO N=1,NM IF(RLAT(N).GE.RLATNP) THEN WNP=WNP+1 IF(IB(K).EQ.0.OR.LO(N,K)) THEN GNP=GNP+GO(N,K) TNP=TNP+1 ENDIF ELSEIF(RLAT(N).LE.RLATSP) THEN WSP=WSP+1 IF(IB(K).EQ.0.OR.LO(N,K)) THEN GSP=GSP+GO(N,K) TSP=TSP+1 ENDIF ENDIF ENDDO C DISTRIBUTE AVERAGE VALUES BACK TO MULTIPLE POLES IF(WNP.GT.1) THEN IF(TNP.GE.WNP/2) THEN GNP=GNP/TNP ELSE GNP=0. ENDIF DO N=1,NM IF(RLAT(N).GE.RLATNP) THEN IF(IB(K).NE.0) LO(N,K)=TNP.GE.WNP/2 GO(N,K)=GNP ENDIF ENDDO ENDIF IF(WSP.GT.1) THEN IF(TSP.GE.WSP/2) THEN GSP=GSP/TSP ELSE GSP=0. ENDIF DO N=1,NM IF(RLAT(N).LE.RLATSP) THEN IF(IB(K).NE.0) LO(N,K)=TSP.GE.WSP/2 GO(N,K)=GSP ENDIF ENDDO ENDIF ENDDO C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - END C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SUBROUTINE GDSWIZ00(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: GDSWIZ00 GDS WIZARD FOR EQUIDISTANT 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 EQUIDISTANT 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 C PROGRAM HISTORY LOG: C 96-04-10 IREDELL C C USAGE: CALL GDSWIZ00(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, C & LROT,CROT,SROT) 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 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 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) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF(KGDS(1).EQ.000) 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 ISCAN=MOD(KGDS(11)/128,2) JSCAN=MOD(KGDS(11)/64,2) NSCAN=MOD(KGDS(11)/32,2) HI=(-1.)**ISCAN HJ=(-1.)**(1-JSCAN) DLON=HI*(MOD(HI*(RLON2-RLON1)-1+3600,360.)+1)/(IM-1) DLAT=(RLAT2-RLAT1)/(JM-1) XMIN=0 XMAX=IM+1 IF(IM.EQ.NINT(360/ABS(DLON))) XMAX=IM+2 YMIN=0 YMAX=JM+1 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.) RLAT(N)=MIN(MAX(RLAT1+DLAT*(YPTS(N)-1),-90.),90.) NRET=NRET+1 IF(LROT.EQ.1) THEN CROT(N)=1 SROT(N)=0 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 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 YPTS(N)=1+(RLAT(N)-RLAT1)/DLAT 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 ELSE XPTS(N)=FILL YPTS(N)=FILL ENDIF ELSE XPTS(N)=FILL YPTS(N)=FILL 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 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SUBROUTINE GDSWIZ01(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: GDSWIZ01 GDS WIZARD FOR MERCATOR 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 MERCATOR 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 C PROGRAM HISTORY LOG: C 96-04-10 IREDELL C 96-10-01 IREDELL PROTECTED AGAINST UNRESOLVABLE POINTS C C USAGE: CALL GDSWIZ01(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, C & LROT,CROT,SROT) 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 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 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) PARAMETER(RERTH=6.3712E6) PARAMETER(PI=3.14159265358979,DPR=180./PI) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF(KGDS(1).EQ.001) 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 RLATI=KGDS(9)*1.E-3 ISCAN=MOD(KGDS(11)/128,2) JSCAN=MOD(KGDS(11)/64,2) NSCAN=MOD(KGDS(11)/32,2) DX=KGDS(12) DY=KGDS(13) HI=(-1.)**ISCAN HJ=(-1.)**(1-JSCAN) DLON=HI*(MOD(HI*(RLON2-RLON1)-1+3600,360.)+1)/(IM-1) DLAT=HJ*DY/(RERTH*COS(RLATI/DPR)) YE=1-LOG(TAN((RLAT1+90)/2/DPR))/DLAT XMIN=0 XMAX=IM+1 IF(IM.EQ.NINT(360/ABS(DLON))) XMAX=IM+2 YMIN=0 YMAX=JM+1 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.) RLAT(N)=2*ATAN(EXP(DLAT*(YPTS(N)-YE)))*DPR-90 NRET=NRET+1 IF(LROT.EQ.1) THEN CROT(N)=1 SROT(N)=0 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 IF(ABS(RLON(N)).LE.360.AND.ABS(RLAT(N)).LT.90) THEN XPTS(N)=1+HI*MOD(HI*(RLON(N)-RLON1)+3600,360.)/DLON YPTS(N)=YE+LOG(TAN((RLAT(N)+90)/2/DPR))/DLAT 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 ELSE XPTS(N)=FILL YPTS(N)=FILL ENDIF ELSE XPTS(N)=FILL YPTS(N)=FILL 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 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SUBROUTINE GDSWIZ03(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: GDSWIZ03 GDS WIZARD FOR LAMBERT CONFORMAL CONICAL 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 LAMBERT CONFORMAL CONICAL 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 C PROGRAM HISTORY LOG: C 96-04-10 IREDELL C 96-10-01 IREDELL PROTECTED AGAINST UNRESOLVABLE POINTS C 1999-04-27 GILBERT CORRECTED MINOR ERROR CALCULATING VARIABLE AN C FOR THE SECANT PROJECTION CASE (RLATI1.NE.RLATI2). C C USAGE: CALL GDSWIZ03(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, C & LROT,CROT,SROT) 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 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 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) PARAMETER(RERTH=6.3712E6) PARAMETER(PI=3.14159265358979,DPR=180./PI) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF(KGDS(1).EQ.003) THEN IM=KGDS(2) JM=KGDS(3) RLAT1=KGDS(4)*1.E-3 RLON1=KGDS(5)*1.E-3 IROT=MOD(KGDS(6)/8,2) ORIENT=KGDS(7)*1.E-3 DX=KGDS(8) DY=KGDS(9) IPROJ=MOD(KGDS(10)/128,2) ISCAN=MOD(KGDS(11)/128,2) JSCAN=MOD(KGDS(11)/64,2) NSCAN=MOD(KGDS(11)/32,2) RLATI1=KGDS(12)*1.E-3 RLATI2=KGDS(13)*1.E-3 H=(-1.)**IPROJ HI=(-1.)**ISCAN HJ=(-1.)**(1-JSCAN) DXS=DX*HI DYS=DY*HJ IF(RLATI1.EQ.RLATI2) THEN AN=SIN(H*RLATI1/DPR) ELSE AN=LOG(COS(RLATI1/DPR)/COS(RLATI2/DPR))/ & LOG(TAN((H*90-RLATI1)/2/DPR)/TAN((H*90-RLATI2)/2/DPR)) ENDIF DE=RERTH*COS(RLATI1/DPR)*TAN((H*RLATI1+90)/2/DPR)**AN/AN IF(H*RLAT1.EQ.90) THEN XP=1 YP=1 ELSE DR=DE/TAN((H*RLAT1+90)/2/DPR)**AN DLON1=MOD(RLON1-ORIENT+180+3600,360.)-180 XP=1-H*SIN(AN*DLON1/DPR)*DR/DXS YP=1+COS(AN*DLON1/DPR)*DR/DYS ENDIF ANTR=1/(2*AN) DE2=DE**2 XMIN=0 XMAX=IM+1 YMIN=0 YMAX=JM+1 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 DI=(XPTS(N)-XP)*DXS DJ=(YPTS(N)-YP)*DYS DR2=DI**2+DJ**2 IF(DR2.LT.DE2*1.E-6) THEN RLON(N)=0. RLAT(N)=H*90. ELSE RLON(N)=MOD(ORIENT+H/AN*DPR*ATAN2(DI,-DJ)+3600,360.) RLAT(N)=H*(2*DPR*ATAN((DE2/DR2)**ANTR)-90) ENDIF NRET=NRET+1 IF(LROT.EQ.1) THEN IF(IROT.EQ.1) THEN DLON=MOD(RLON(N)-ORIENT+180+3600,360.)-180 CROT(N)=H*COS(AN*DLON/DPR) SROT(N)=SIN(AN*DLON/DPR) ELSE CROT(N)=1 SROT(N)=0 ENDIF 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 IF(ABS(RLON(N)).LE.360.AND.ABS(RLAT(N)).LE.90.AND. & H*RLAT(N).NE.-90) THEN DR=DE*TAN((90-H*RLAT(N))/2/DPR)**AN DLON=MOD(RLON(N)-ORIENT+180+3600,360.)-180 XPTS(N)=XP+H*SIN(AN*DLON/DPR)*DR/DXS YPTS(N)=YP-COS(AN*DLON/DPR)*DR/DYS 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 IF(IROT.EQ.1) THEN CROT(N)=H*COS(AN*DLON/DPR) SROT(N)=SIN(AN*DLON/DPR) ELSE CROT(N)=1 SROT(N)=0 ENDIF ENDIF ELSE XPTS(N)=FILL YPTS(N)=FILL ENDIF ELSE XPTS(N)=FILL YPTS(N)=FILL 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 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SUBROUTINE GDSWIZ05(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: GDSWIZ05 GDS WIZARD FOR POLAR STEREOGRAPHIC AZIMUTHAL 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 POLAR STEREOGRAPHIC AZIMUTHAL 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 C PROGRAM HISTORY LOG: C 96-04-10 IREDELL C C USAGE: CALL GDSWIZ05(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, C & LROT,CROT,SROT) 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 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 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) PARAMETER(RERTH=6.3712E6) PARAMETER(PI=3.14159265358979,DPR=180./PI) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF(KGDS(1).EQ.005) THEN IM=KGDS(2) JM=KGDS(3) RLAT1=KGDS(4)*1.E-3 RLON1=KGDS(5)*1.E-3 IROT=MOD(KGDS(6)/8,2) ORIENT=KGDS(7)*1.E-3 DX=KGDS(8) DY=KGDS(9) IPROJ=MOD(KGDS(10)/128,2) ISCAN=MOD(KGDS(11)/128,2) JSCAN=MOD(KGDS(11)/64,2) NSCAN=MOD(KGDS(11)/32,2) H=(-1.)**IPROJ HI=(-1.)**ISCAN HJ=(-1.)**(1-JSCAN) DXS=DX*HI DYS=DY*HJ DE=(1.+SIN(60./DPR))*RERTH DR=DE*COS(RLAT1/DPR)/(1+H*SIN(RLAT1/DPR)) XP=1-H*SIN((RLON1-ORIENT)/DPR)*DR/DXS YP=1+COS((RLON1-ORIENT)/DPR)*DR/DYS DE2=DE**2 XMIN=0 XMAX=IM+1 YMIN=0 YMAX=JM+1 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 DI=(XPTS(N)-XP)*DXS DJ=(YPTS(N)-YP)*DYS DR2=DI**2+DJ**2 IF(DR2.LT.DE2*1.E-6) THEN RLON(N)=0. RLAT(N)=H*90. ELSE RLON(N)=MOD(ORIENT+H*DPR*ATAN2(DI,-DJ)+3600,360.) RLAT(N)=H*DPR*ASIN((DE2-DR2)/(DE2+DR2)) ENDIF NRET=NRET+1 IF(LROT.EQ.1) THEN IF(IROT.EQ.1) THEN CROT(N)=H*COS((RLON(N)-ORIENT)/DPR) SROT(N)=SIN((RLON(N)-ORIENT)/DPR) ELSE CROT(N)=1 SROT(N)=0 ENDIF 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 IF(ABS(RLON(N)).LE.360.AND.ABS(RLAT(N)).LE.90.AND. & H*RLAT(N).NE.-90) THEN DR=DE*TAN((90-H*RLAT(N))/2/DPR) XPTS(N)=XP+H*SIN((RLON(N)-ORIENT)/DPR)*DR/DXS YPTS(N)=YP-COS((RLON(N)-ORIENT)/DPR)*DR/DYS 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 IF(IROT.EQ.1) THEN CROT(N)=H*COS((RLON(N)-ORIENT)/DPR) SROT(N)=SIN((RLON(N)-ORIENT)/DPR) ELSE CROT(N)=1 SROT(N)=0 ENDIF ENDIF ELSE XPTS(N)=FILL YPTS(N)=FILL ENDIF ELSE XPTS(N)=FILL YPTS(N)=FILL 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 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SUBROUTINE GDSWIZC9(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: GDSWIZC9 GDS WIZARD FOR ROTATED EQUIDISTANT 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 STAGGERED ROTATED EQUIDISTANT CYLINDRICAL PROJECTIONS. C (SEE UNDER THE DESCRIPTION OF KGDS TO DETERMINE WHETHER C TO COMPUTE A STAGGERED WIND GRID OR A STAGGERED MASS GRID.) 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 C PROGRAM HISTORY LOG: C 96-04-10 IREDELL C C USAGE: CALL GDSWIZC9(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, C & LROT,CROT,SROT) C C INPUT ARGUMENT LIST: C KGDS - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63 C IMPORTANT NOTE: IF THE 9TH BIT (FROM RIGHT) OF KGDS(11) C (SCANNING MODE FLAG) IS 1, THEN THIS C THE GRID IS COMPUTED FOR A WIND FIELD; C OTHERWISE IT IS FOR A MASS FIELD. THUS C MOD(KGDS(11)/256,2)=0 FOR MASS GRID. 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 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 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) PARAMETER(RERTH=6.3712E6) PARAMETER(PI=3.14159265358979,DPR=180./PI) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF(KGDS(1).EQ.201) THEN RLAT1=KGDS(4)*1.E-3 RLON1=KGDS(5)*1.E-3 IROT=MOD(KGDS(6)/8,2) IM=KGDS(7)*2-1 JM=KGDS(8) DLON=KGDS(9)*1.E-3 DLAT=KGDS(10)*1.E-3 KSCAN=MOD(KGDS(11)/256,2) ISCAN=MOD(KGDS(11)/128,2) JSCAN=MOD(KGDS(11)/64,2) NSCAN=MOD(KGDS(11)/32,2) HI=(-1.)**ISCAN HJ=(-1.)**(1-JSCAN) DLONS=HI*DLON DLATS=HJ*DLAT RLONR=-(IM-1)/2*DLONS RLATR=-(JM-1)/2*DLATS SLAT1=SIN(RLAT1/DPR) CLAT1=COS(RLAT1/DPR) SLONR=SIN(RLONR/DPR) CLONR=COS(RLONR/DPR) SLATR=SIN(RLATR/DPR) CLATR=COS(RLATR/DPR) DENOM=1-(CLATR*SLONR)**2 SLAT0=(SLAT1*CLATR*CLONR-SLATR*SQRT(DENOM-SLAT1**2))/DENOM CLAT0=SQRT(1-SLAT0**2) RLON0=RLON1+HI*DPR*ACOS((CLAT0*CLATR*CLONR-SLAT0*SLATR)/CLAT1) C THE FOLLOWING INDENTED LINES ARE A TEMPORARY FIX OF IMPRECISE GRID. C CAUTION: CENTRAL LATITUDE AND LONGITUDE ARE ASSUMED TO BE INTEGERS. SLAT0=SIN(NINT(ASIN(SLAT0)*DPR)/DPR) CLAT0=SQRT(1-SLAT0**2) RLON0=NINT(RLON0) HS=SIGN(1.,MOD(RLON1-RLON0+180+3600,360.)-180) CLON1=COS((RLON1-RLON0)/DPR) SLATR=CLAT0*SLAT1-SLAT0*CLAT1*CLON1 CLATR=SQRT(1-SLATR**2) CLONR=(CLAT0*CLAT1*CLON1+SLAT0*SLAT1)/CLATR RLATR=DPR*ASIN(SLATR) RLONR=HS*DPR*ACOS(CLONR) DLATS=RLATR/(-(JM-1)/2) DLONS=RLONR/(-(IM-1)/2) IF(KSCAN.EQ.0) THEN IS1=(JM+1)/2 ELSE IS1=JM/2 ENDIF XMIN=0 XMAX=IM+1 IF(IM.EQ.NINT(360/ABS(DLONS))) XMAX=IM+2 YMIN=0 YMAX=JM+1 NRET=0 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C TRANSLATE GRID COORDINATES TO EARTH COORDINATES IF(IOPT.EQ.0.OR.IOPT.EQ.1) THEN DO N=1,NPTS XPTF=YPTS(N)+(XPTS(N)-IS1) YPTF=YPTS(N)-(XPTS(N)-IS1)+KSCAN IF(XPTF.GE.XMIN.AND.XPTF.LE.XMAX.AND. & YPTF.GE.YMIN.AND.YPTF.LE.YMAX) THEN HS=HI*SIGN(1.,XPTF-(IM+1)/2) RLONR=(XPTF-(IM+1)/2)*DLONS RLATR=(YPTF-(JM+1)/2)*DLATS CLONR=COS(RLONR/DPR) SLATR=SIN(RLATR/DPR) CLATR=COS(RLATR/DPR) SLAT=CLAT0*SLATR+SLAT0*CLATR*CLONR IF(SLAT.LE.-1) THEN CLAT=0. CLON=COS(RLON0/DPR) RLON(N)=0 RLAT(N)=-90 ELSEIF(SLAT.GE.1) THEN CLAT=0. CLON=COS(RLON0/DPR) RLON(N)=0 RLAT(N)=90 ELSE CLAT=SQRT(1-SLAT**2) CLON=(CLAT0*CLATR*CLONR-SLAT0*SLATR)/CLAT CLON=MIN(MAX(CLON,-1.),1.) RLON(N)=MOD(RLON0+HS*DPR*ACOS(CLON)+3600,360.) RLAT(N)=DPR*ASIN(SLAT) ENDIF NRET=NRET+1 IF(LROT.EQ.1) THEN IF(IROT.EQ.1) THEN IF(CLATR.LE.0) THEN CROT(N)=-SIGN(1.,SLATR*SLAT0) SROT(N)=0 ELSE SLON=SIN((RLON(N)-RLON0)/DPR) CROT(N)=(CLAT0*CLAT+SLAT0*SLAT*CLON)/CLATR SROT(N)=SLAT0*SLON/CLATR ENDIF ELSE CROT(N)=1 SROT(N)=0 ENDIF 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 IF(ABS(RLON(N)).LE.360.AND.ABS(RLAT(N)).LE.90) THEN HS=SIGN(1.,MOD(RLON(N)-RLON0+180+3600,360.)-180) CLON=COS((RLON(N)-RLON0)/DPR) SLAT=SIN(RLAT(N)/DPR) CLAT=COS(RLAT(N)/DPR) SLATR=CLAT0*SLAT-SLAT0*CLAT*CLON IF(SLATR.LE.-1) THEN CLATR=0. RLONR=0 RLATR=-90 ELSEIF(SLATR.GE.1) THEN CLATR=0. RLONR=0 RLATR=90 ELSE CLATR=SQRT(1-SLATR**2) CLONR=(CLAT0*CLAT*CLON+SLAT0*SLAT)/CLATR CLONR=MIN(MAX(CLONR,-1.),1.) RLONR=HS*DPR*ACOS(CLONR) RLATR=DPR*ASIN(SLATR) ENDIF XPTF=(IM+1)/2+RLONR/DLONS YPTF=(JM+1)/2+RLATR/DLATS IF(XPTF.GE.XMIN.AND.XPTF.LE.XMAX.AND. & YPTF.GE.YMIN.AND.YPTF.LE.YMAX) THEN XPTS(N)=IS1+(XPTF-(YPTF-KSCAN))/2 YPTS(N)=(XPTF+(YPTF-KSCAN))/2 NRET=NRET+1 IF(LROT.EQ.1) THEN IF(IROT.EQ.1) THEN IF(CLATR.LE.0) THEN CROT(N)=-SIGN(1.,SLATR*SLAT0) SROT(N)=0 ELSE SLON=SIN((RLON(N)-RLON0)/DPR) CROT(N)=(CLAT0*CLAT+SLAT0*SLAT*CLON)/CLATR SROT(N)=SLAT0*SLON/CLATR ENDIF ELSE CROT(N)=1 SROT(N)=0 ENDIF ENDIF ELSE XPTS(N)=FILL YPTS(N)=FILL ENDIF ELSE XPTS(N)=FILL YPTS(N)=FILL 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 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SUBROUTINE GDSWIZCA(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: GDSWIZCA GDS WIZARD FOR ROTATED EQUIDISTANT 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 ROTATED EQUIDISTANT 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 C PROGRAM HISTORY LOG: C 96-04-10 IREDELL C C USAGE: CALL GDSWIZCA(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, C & LROT,CROT,SROT) 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 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 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) PARAMETER(RERTH=6.3712E6) PARAMETER(PI=3.14159265358979,DPR=180./PI) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF(KGDS(1).EQ.202) THEN RLAT1=KGDS(4)*1.E-3 RLON1=KGDS(5)*1.E-3 IROT=MOD(KGDS(6)/8,2) IM=KGDS(7) JM=KGDS(8) DLON=KGDS(9)*1.E-3 DLAT=KGDS(10)*1.E-3 ISCAN=MOD(KGDS(11)/128,2) JSCAN=MOD(KGDS(11)/64,2) NSCAN=MOD(KGDS(11)/32,2) HI=(-1.)**ISCAN HJ=(-1.)**(1-JSCAN) DLONS=HI*DLON DLATS=HJ*DLAT RLONR=-(IM-1)/2*DLONS RLATR=-(JM-1)/2*DLATS SLAT1=SIN(RLAT1/DPR) CLAT1=COS(RLAT1/DPR) SLONR=SIN(RLONR/DPR) CLONR=COS(RLONR/DPR) SLATR=SIN(RLATR/DPR) CLATR=COS(RLATR/DPR) DENOM=1-(CLATR*SLONR)**2 SLAT0=(SLAT1*CLATR*CLONR-SLATR*SQRT(DENOM-SLAT1**2))/DENOM CLAT0=SQRT(1-SLAT0**2) RLON0=RLON1+HI*DPR*ACOS((CLAT0*CLATR*CLONR-SLAT0*SLATR)/CLAT1) C THE FOLLOWING INDENTED LINES ARE A TEMPORARY FIX OF IMPRECISE GRID. C CAUTION: CENTRAL LATITUDE AND LONGITUDE ARE ASSUMED TO BE INTEGERS. SLAT0=SIN(NINT(ASIN(SLAT0)*DPR)/DPR) CLAT0=SQRT(1-SLAT0**2) RLON0=NINT(RLON0) HS=SIGN(1.,MOD(RLON1-RLON0+180+3600,360.)-180) CLON1=COS((RLON1-RLON0)/DPR) SLATR=CLAT0*SLAT1-SLAT0*CLAT1*CLON1 CLATR=SQRT(1-SLATR**2) CLONR=(CLAT0*CLAT1*CLON1+SLAT0*SLAT1)/CLATR RLATR=DPR*ASIN(SLATR) RLONR=HS*DPR*ACOS(CLONR) DLATS=RLATR/(-(JM-1)/2) DLONS=RLONR/(-(IM-1)/2) XMIN=0 XMAX=IM+1 IF(IM.EQ.NINT(360/ABS(DLONS))) XMAX=IM+2 YMIN=0 YMAX=JM+1 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 HS=HI*SIGN(1.,XPTS(N)-(IM+1)/2) RLONR=(XPTS(N)-(IM+1)/2)*DLONS RLATR=(YPTS(N)-(JM+1)/2)*DLATS CLONR=COS(RLONR/DPR) SLATR=SIN(RLATR/DPR) CLATR=COS(RLATR/DPR) SLAT=CLAT0*SLATR+SLAT0*CLATR*CLONR IF(SLAT.LE.-1) THEN CLAT=0. CLON=COS(RLON0/DPR) RLON(N)=0 RLAT(N)=-90 ELSEIF(SLAT.GE.1) THEN CLAT=0. CLON=COS(RLON0/DPR) RLON(N)=0 RLAT(N)=90 ELSE CLAT=SQRT(1-SLAT**2) CLON=(CLAT0*CLATR*CLONR-SLAT0*SLATR)/CLAT CLON=MIN(MAX(CLON,-1.),1.) RLON(N)=MOD(RLON0+HS*DPR*ACOS(CLON)+3600,360.) RLAT(N)=DPR*ASIN(SLAT) ENDIF NRET=NRET+1 IF(LROT.EQ.1) THEN IF(IROT.EQ.1) THEN IF(CLATR.LE.0) THEN CROT(N)=-SIGN(1.,SLATR*SLAT0) SROT(N)=0 ELSE SLON=SIN((RLON(N)-RLON0)/DPR) CROT(N)=(CLAT0*CLAT+SLAT0*SLAT*CLON)/CLATR SROT(N)=SLAT0*SLON/CLATR ENDIF ELSE CROT(N)=1 SROT(N)=0 ENDIF 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 IF(ABS(RLON(N)).LE.360.AND.ABS(RLAT(N)).LE.90) THEN HS=SIGN(1.,MOD(RLON(N)-RLON0+180+3600,360.)-180) CLON=COS((RLON(N)-RLON0)/DPR) SLAT=SIN(RLAT(N)/DPR) CLAT=COS(RLAT(N)/DPR) SLATR=CLAT0*SLAT-SLAT0*CLAT*CLON IF(SLATR.LE.-1) THEN CLATR=0. RLONR=0 RLATR=-90 ELSEIF(SLATR.GE.1) THEN CLATR=0. RLONR=0 RLATR=90 ELSE CLATR=SQRT(1-SLATR**2) CLONR=(CLAT0*CLAT*CLON+SLAT0*SLAT)/CLATR CLONR=MIN(MAX(CLONR,-1.),1.) RLONR=HS*DPR*ACOS(CLONR) RLATR=DPR*ASIN(SLATR) ENDIF XPTS(N)=(IM+1)/2+RLONR/DLONS YPTS(N)=(JM+1)/2+RLATR/DLATS 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 IF(IROT.EQ.1) THEN IF(CLATR.LE.0) THEN CROT(N)=-SIGN(1.,SLATR*SLAT0) SROT(N)=0 ELSE SLON=SIN((RLON(N)-RLON0)/DPR) CROT(N)=(CLAT0*CLAT+SLAT0*SLAT*CLON)/CLATR SROT(N)=SLAT0*SLON/CLATR ENDIF ELSE CROT(N)=1 SROT(N)=0 ENDIF ENDIF ELSE XPTS(N)=FILL YPTS(N)=FILL ENDIF ELSE XPTS(N)=FILL YPTS(N)=FILL 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 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - SUBROUTINE GDSWIZCB(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, & LROT,CROT,SROT) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: GDSWIZCB GDS WIZARD FOR ROTATED EQUIDISTANT 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 STAGGERED ROTATED EQUIDISTANT CYLINDRICAL PROJECTIONS. C (SEE UNDER THE DESCRIPTION OF KGDS TO DETERMINE WHETHER C TO COMPUTE A STAGGERED WIND GRID OR A STAGGERED MASS GRID.) 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 C PROGRAM HISTORY LOG: C 96-04-10 IREDELL C 98-08-19 BALDWIN MODIFY GDSWIZC9 FOR TYPE 203 ETA GRIDS C C USAGE: CALL GDSWIZCB(KGDS,IOPT,NPTS,FILL,XPTS,YPTS,RLON,RLAT,NRET, C & LROT,CROT,SROT) C C INPUT ARGUMENT LIST: C KGDS - INTEGER (200) GDS PARAMETERS AS DECODED BY W3FI63 C IMPORTANT NOTE: IF THE 9TH BIT (FROM RIGHT) OF KGDS(11) C (SCANNING MODE FLAG) IS 1, THEN THIS C THE GRID IS COMPUTED FOR A WIND FIELD; C OTHERWISE IT IS FOR A MASS FIELD. THUS C MOD(KGDS(11)/256,2)=0 FOR MASS GRID. 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 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 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) PARAMETER(RERTH=6.3712E6) PARAMETER(PI=3.14159265358979,DPR=180./PI) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - IF(KGDS(1).EQ.203) THEN RLAT1=KGDS(4)*1.E-3 RLON1=KGDS(5)*1.E-3 RLAT0=KGDS(7)*1.E-3 RLON0=KGDS(8)*1.E-3 IROT=MOD(KGDS(6)/8,2) IM=KGDS(2)*2-1 JM=KGDS(3) KSCAN=MOD(KGDS(11)/256,2) ISCAN=MOD(KGDS(11)/128,2) JSCAN=MOD(KGDS(11)/64,2) NSCAN=MOD(KGDS(11)/32,2) HI=(-1.)**ISCAN HJ=(-1.)**(1-JSCAN) SLAT1=SIN(RLAT1/DPR) CLAT1=COS(RLAT1/DPR) SLAT0=SIN(RLAT0/DPR) CLAT0=COS(RLAT0/DPR) HS=SIGN(1.,MOD(RLON1-RLON0+180+3600,360.)-180) CLON1=COS((RLON1-RLON0)/DPR) SLATR=CLAT0*SLAT1-SLAT0*CLAT1*CLON1 CLATR=SQRT(1-SLATR**2) CLONR=(CLAT0*CLAT1*CLON1+SLAT0*SLAT1)/CLATR RLATR=DPR*ASIN(SLATR) RLONR=HS*DPR*ACOS(CLONR) DLATS=RLATR/(-(JM-1)/2) DLONS=RLONR/(-(IM-1)/2) IF(KSCAN.EQ.0) THEN IS1=(JM+1)/2 ELSE IS1=JM/2 ENDIF XMIN=0 XMAX=IM+1 IF(IM.EQ.NINT(360/ABS(DLONS))) XMAX=IM+2 YMIN=0 YMAX=JM+1 NRET=0 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C TRANSLATE GRID COORDINATES TO EARTH COORDINATES IF(IOPT.EQ.0.OR.IOPT.EQ.1) THEN DO N=1,NPTS XPTF=YPTS(N)+(XPTS(N)-IS1) YPTF=YPTS(N)-(XPTS(N)-IS1)+KSCAN IF(XPTF.GE.XMIN.AND.XPTF.LE.XMAX.AND. & YPTF.GE.YMIN.AND.YPTF.LE.YMAX) THEN HS=HI*SIGN(1.,XPTF-(IM+1)/2) RLONR=(XPTF-(IM+1)/2)*DLONS RLATR=(YPTF-(JM+1)/2)*DLATS CLONR=COS(RLONR/DPR) SLATR=SIN(RLATR/DPR) CLATR=COS(RLATR/DPR) SLAT=CLAT0*SLATR+SLAT0*CLATR*CLONR IF(SLAT.LE.-1) THEN CLAT=0. CLON=COS(RLON0/DPR) RLON(N)=0 RLAT(N)=-90 ELSEIF(SLAT.GE.1) THEN CLAT=0. CLON=COS(RLON0/DPR) RLON(N)=0 RLAT(N)=90 ELSE CLAT=SQRT(1-SLAT**2) CLON=(CLAT0*CLATR*CLONR-SLAT0*SLATR)/CLAT CLON=MIN(MAX(CLON,-1.),1.) RLON(N)=MOD(RLON0+HS*DPR*ACOS(CLON)+3600,360.) RLAT(N)=DPR*ASIN(SLAT) ENDIF NRET=NRET+1 IF(LROT.EQ.1) THEN IF(IROT.EQ.1) THEN IF(CLATR.LE.0) THEN CROT(N)=-SIGN(1.,SLATR*SLAT0) SROT(N)=0 ELSE SLON=SIN((RLON(N)-RLON0)/DPR) CROT(N)=(CLAT0*CLAT+SLAT0*SLAT*CLON)/CLATR SROT(N)=SLAT0*SLON/CLATR ENDIF ELSE CROT(N)=1 SROT(N)=0 ENDIF 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 IF(ABS(RLON(N)).LE.360.AND.ABS(RLAT(N)).LE.90) THEN HS=SIGN(1.,MOD(RLON(N)-RLON0+180+3600,360.)-180) CLON=COS((RLON(N)-RLON0)/DPR) SLAT=SIN(RLAT(N)/DPR) CLAT=COS(RLAT(N)/DPR) SLATR=CLAT0*SLAT-SLAT0*CLAT*CLON IF(SLATR.LE.-1) THEN CLATR=0. RLONR=0 RLATR=-90 ELSEIF(SLATR.GE.1) THEN CLATR=0. RLONR=0 RLATR=90 ELSE CLATR=SQRT(1-SLATR**2) CLONR=(CLAT0*CLAT*CLON+SLAT0*SLAT)/CLATR CLONR=MIN(MAX(CLONR,-1.),1.) RLONR=HS*DPR*ACOS(CLONR) RLATR=DPR*ASIN(SLATR) ENDIF XPTF=(IM+1)/2+RLONR/DLONS YPTF=(JM+1)/2+RLATR/DLATS IF(XPTF.GE.XMIN.AND.XPTF.LE.XMAX.AND. & YPTF.GE.YMIN.AND.YPTF.LE.YMAX) THEN XPTS(N)=IS1+(XPTF-(YPTF-KSCAN))/2 YPTS(N)=(XPTF+(YPTF-KSCAN))/2 NRET=NRET+1 IF(LROT.EQ.1) THEN IF(IROT.EQ.1) THEN IF(CLATR.LE.0) THEN CROT(N)=-SIGN(1.,SLATR*SLAT0) SROT(N)=0 ELSE SLON=SIN((RLON(N)-RLON0)/DPR) CROT(N)=(CLAT0*CLAT+SLAT0*SLAT*CLON)/CLATR SROT(N)=SLAT0*SLON/CLATR ENDIF ELSE CROT(N)=1 SROT(N)=0 ENDIF ENDIF ELSE XPTS(N)=FILL YPTS(N)=FILL ENDIF ELSE XPTS(N)=FILL YPTS(N)=FILL 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 C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C----------------------------------------------------------------------- FUNCTION IJKGDS1(I,J,IJKGDSA) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C C SUBPROGRAM: IJKGDS1 RETURN FIELD POSITION FOR A GIVEN GRID POINT C PRGMMR: IREDELL ORG: W/NMC23 DATE: 96-04-10 C C ABSTRACT: THIS SUBPROGRAM DECODES THE GRIB GRID DESCRIPTION SECTION C AND RETURNS THE FIELD POSITION FOR A GIVEN GRID POINT. C CALL IJKGDS0 TO SET UP THE NAVIGATION PARAMETER ARRAY. C C PROGRAM HISTORY LOG: C 96-04-10 IREDELL C 97-03-11 IREDELL ALLOWED HEMISPHERIC GRIDS TO WRAP OVER ONE POLE C 98-07-13 BALDWIN ADD 2D STAGGERED ETA GRID INDEXING (203) C 1999-04-08 IREDELL SPLIT IJKGDS INTO TWO C C USAGE: ...IJKGDS1(I,J,IJKGDSA) C C INPUT ARGUMENT LIST: C I - INTEGER X GRID POINT C J - INTEGER Y GRID POINT C IJKGDSA - INTEGER (20) NAVIGATION PARAMETER ARRAY C IJKGDSA(1) IS NUMBER OF X POINTS C IJKGDSA(2) IS NUMBER OF Y POINTS C IJKGDSA(3) IS X WRAPAROUND INCREMENT C (0 IF NO WRAPAROUND) C IJKGDSA(4) IS Y WRAPAROUND LOWER PIVOT POINT C (0 IF NO WRAPAROUND) C IJKGDSA(5) IS Y WRAPAROUND UPPER PIVOT POINT C (0 IF NO WRAPAROUND) C IJKGDSA(6) IS SCANNING MODE C (0 IF X FIRST THEN Y; 1 IF Y FIRST THEN X; C 2 IF STAGGERED DIAGONAL LIKE PROJECTION 201; C 3 IF STAGGERED DIAGONAL LIKE PROJECTION 203) C IJKGDSA(7) IS MASS/WIND FLAG FOR STAGGERED DIAGONAL C (0 IF MASS; 1 IF WIND) C IJKGDSA(8:20) ARE UNUSED AT THE MOMENT C C OUTPUT ARGUMENT LIST: C IJKGDS - INTEGER POSITION IN GRIB FIELD TO LOCATE GRID POINT C (0 IF OUT OF BOUNDS) C C ATTRIBUTES: C LANGUAGE: FORTRAN 90 C C$$$ INTEGER I,J INTEGER IJKGDSA(20) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C EXTRACT FROM NAVIGATION PARAMETER ARRAY IM=IJKGDSA(1) JM=IJKGDSA(2) IWRAP=IJKGDSA(3) JWRAP1=IJKGDSA(4) JWRAP2=IJKGDSA(5) NSCAN=IJKGDSA(6) KSCAN=IJKGDSA(7) C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C COMPUTE WRAPAROUNDS IN X AND Y IF NECESSARY AND POSSIBLE II=I JJ=J IF(IWRAP.GT.0) THEN II=MOD(I-1+IWRAP,IWRAP)+1 IF(J.LT.1.AND.JWRAP1.GT.0) THEN JJ=JWRAP1-J II=MOD(II-1+IWRAP/2,IWRAP)+1 ELSEIF(J.GT.JM.AND.JWRAP2.GT.0) THEN JJ=JWRAP2-J II=MOD(II-1+IWRAP/2,IWRAP)+1 ENDIF ENDIF C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C COMPUTE POSITION FOR THE APPROPRIATE SCANNING MODE IJKGDS1=0 IF(NSCAN.EQ.0) THEN IF(II.GE.1.AND.II.LE.IM.AND.JJ.GE.1.AND.JJ.LE.JM) & IJKGDS1=II+(JJ-1)*IM ELSEIF(NSCAN.EQ.1) THEN IF(II.GE.1.AND.II.LE.IM.AND.JJ.GE.1.AND.JJ.LE.JM) & IJKGDS1=JJ+(II-1)*JM ELSEIF(NSCAN.EQ.2) THEN IS1=(JM+1-KSCAN)/2 IIF=JJ+(II-IS1) JJF=JJ-(II-IS1)+KSCAN IF(IIF.GE.1.AND.IIF.LE.2*IM-1.AND.JJF.GE.1.AND.JJF.LE.JM) & IJKGDS1=(IIF+(JJF-1)*(2*IM-1)+1-KSCAN)/2 ELSEIF(NSCAN.EQ.3) THEN IS1=(JM+1-KSCAN)/2 IIF=JJ+(II-IS1) JJF=JJ-(II-IS1)+KSCAN IF(IIF.GE.1.AND.IIF.LE.2*IM-1.AND.JJF.GE.1.AND.JJF.LE.JM) & IJKGDS1=(IIF+1)/2+(JJF-1)*IM ENDIF C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - END C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -