forked from erget/wgrib2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathipsector.f
executable file
·130 lines (130 loc) · 4.59 KB
/
ipsector.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
C-----------------------------------------------------------------------
SUBROUTINE IPSECTOR(I1,I2,J1,J2,NF,M,KGDS,L,F,MS,
& KGDSS,LS,FS,IRET)
C$$$ SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM: IPSECTOR CREATE A SUBSECTOR OF A GRID
C PRGMMR: IREDELL ORG: W/NMC23 DATE: 1998-04-08
C
C ABSTRACT: THIS SUBPROGRAM CREATES A SUBSECTOR OF A GRID,
C COMPUTING A NEW GRID DEFINITION SECTION AND COPYING DATA
C FROM THE ORIGINAL GRID DATA.
C
C PROGRAM HISTORY LOG:
C 1999-04-08 IREDELL
C
C USAGE: CALL IPSECTOR(I1,I2,J1,J2,NF,M,KGDS,L,F,MS,
C & KGDSS,LS,FS,IRET)
C
C INPUT ARGUMENT LIST:
C I1 - INTEGER FIRST X POINT OF THE SECTOR
C OR IF 0<=I2<I1,
C THE TOTAL NUMBER OF NON-OVERLAPPING X SECTORS
C OR IF 0<=I2<-I1,
C THE NEGATIVE TOTAL NUMBER OF OVERLAPPING X SECTORS.
C I2 - INTEGER LAST X POINT OF THE SECTOR
C OR IF 0<=I2<ABS(I1), THE SECTOR NUMBER.
C J1 - INTEGER FIRST Y POINT OF THE SECTOR
C OR IF 0<=J2<J1,
C THE TOTAL NUMBER OF NON-OVERLAPPING Y SECTORS
C OR IF 0<=J2<-J1,
C THE NEGATIVE TOTAL NUMBER OF OVERLAPPING Y SECTORS.
C J2 - INTEGER LAST Y POINT OF THE SECTOR
C OR IF 0<=J2<ABS(J1), THE SECTOR NUMBER.
C NF - INTEGER NUMBER OF FIELDS TO CUT
C M - INTEGER FIRST DIMENSION OF INPUT FIELD ARRAYS
C KGDS - INTEGER (200) INPUT GDS PARAMETERS AS DECODED BY W3FI63
C L - LOGICAL(1) (M,NF) BITMAP FOR INPUT FIELD
C F - REAL (M,NF) DATA FOR INPUT FIELD
C MS - INTEGER FIRST DIMENSION OF OUTPUT FIELD ARRAYS
C
C OUTPUT ARGUMENT LIST:
C KGDSS - INTEGER (200) OUTPUT GDS PARAMETERS
C LS - LOGICAL(1) (MS,NF) BITMAP FOR OUTPUT FIELD
C FS - REAL (MS,NF) DATA FOR OUTPUT FIELD
C IRET - INTEGER RETURN CODE
C (0 IF SUCCESSFUL;
C 71 IF UNSUPPORTED PROJECTION;
C 72 IF INVALID SECTOR SPECIFICATION)
C
C ATTRIBUTES:
C LANGUAGE: FORTRAN 90
C
C$$$
INTEGER I1,I2,J1,J2,NF,M,MS
INTEGER KGDS(200)
LOGICAL(1) L(M,NF)
REAL F(M,NF)
INTEGER KGDSS(200)
LOGICAL(1) LS(MS,NF)
REAL FS(MS,NF)
INTEGER IRET
REAL XPTS(2),YPTS(2),RLON(2),RLAT(2)
INTEGER I1A,I2A,INA,J1A,J2A,JNA,K,KS,NS,NSCAN
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C COMPUTE ACTUAL SECTOR BOUNDARIES
IF((KGDS(1).NE.0.AND.KGDS(1).NE.1.AND.KGDS(1).NE.3.AND.
& KGDS(1).NE.4.AND.KGDS(1).NE.5).OR.KGDS(20).NE.255) THEN
IRET=71
RETURN
ENDIF
I1A=I1
I2A=I2
J1A=J1
J2A=J2
IF(I2.GE.0.AND.I1.GT.I2) THEN
I1A=MIN(I2*((KGDS(2)-1)/I1+1)+1,KGDS(2)+1)
I2A=MIN((I2+1)*((KGDS(2)-1)/I1+1),KGDS(2))
ELSEIF(I2.GE.0.AND.-I1.GT.I2) THEN
I1A=MIN(I2*((KGDS(2)-2)/(-I1)+1)+1,KGDS(2)+1)
I2A=MIN((I2+1)*((KGDS(2)-2)/(-I1)+1)+1,KGDS(2))
ENDIF
IF(J2.GE.0.AND.J1.GT.J2) THEN
J1A=MIN(J2*((KGDS(3)-1)/J1+1)+1,KGDS(3)+1)
J2A=MIN((J2+1)*((KGDS(3)-1)/J1+1),KGDS(3))
ELSEIF(J2.GE.0.AND.-J1.GT.J2) THEN
J1A=MIN(J2*((KGDS(3)-2)/(-J1)+1)+1,KGDS(3)+1)
J2A=MIN((J2+1)*((KGDS(3)-2)/(-J1)+1)+1,KGDS(3))
ENDIF
IF(I1A.LT.1.OR.I2A.GT.KGDS(2).OR.I1A.GT.I2A.OR.
& J1A.LT.1.OR.J2A.GT.KGDS(3).OR.J1A.GT.J2A) THEN
IRET=72
RETURN
ENDIF
INA=I2A-I1A+1
JNA=J2A-J1A+1
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C COMPUTE NEW GRID DEFINITION SECTION PARAMETERS
KGDSS=KGDS
XPTS(1)=I1A
YPTS(1)=J1A
XPTS(2)=I2A
YPTS(2)=J2A
CALL GDSWIZ(KGDS,1,2,-9999.,XPTS,YPTS,RLON,RLAT,NRET,0,CROT,SROT)
KGDSS(2)=INA
KGDSS(3)=JNA
KGDSS(4)=NINT(1.E3*RLAT(1))
KGDSS(5)=NINT(1.E3*RLON(1))
IF(KGDS(1).EQ.0.OR.KGDS(1).EQ.1.OR.KGDS(1).EQ.4) THEN
IF(MOD(KGDS(6),2).EQ.0) KGDSS(6)=KGDS(6)+1
KGDSS(7)=NINT(1.E3*RLAT(2))
KGDSS(8)=NINT(1.E3*RLON(2))
ENDIF
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
C COPY BITMAPS AND DATA
NS=INA*JNA
NSCAN=MOD(KGDS(11)/32,2)
DO N=1,NF
DO KS=1,NS
IF(NSCAN.EQ.0) THEN
K=((KS-1)/INA+J1A-1)*KGDS(2)+MOD(KS-1,INA)+I1A-1
ELSE
K=((KS-1)/JNA+I1A-1)*KGDS(3)+MOD(KS-1,JNA)+J1A-1
ENDIF
LS(KS,N)=L(K,N)
FS(KS,N)=F(K,N)
ENDDO
ENDDO
IRET=0
C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
END