NCEPLIBS-w3emc  2.11.0
w3ft16.f
Go to the documentation of this file.
1 C> @file
2 C> @brief Convert (95,91) grid to (3447) grid
3 C> @author Ralph Jones @date 1994-05-03
4 
5 C> Convert a northern hemisphere 1.0 degree lat.,lon. 95 by
6 C> 91 grid to a wafs 1.25 degree thinned 3447 point grid.
7 C>
8 C> ### Program History Log:
9 C> Date | Programmer | Comment
10 C> -----|------------|--------
11 C> 1994-05-03 | Ralph Jones | Initial.
12 C>
13 C> @param[in] ALOLA 95 * 91 grid 1.0 deg. lat,lon grid northern hemisphere 8645 point grid.
14 C> @param[in] INTERP 1 linear interpolation , ne.1 biquadratic
15 C> @param[out] BTHIN 3447 point thinned grid of n. hemispere 3447 grid is for grib grids 37-40.
16 C>
17 C> @note
18 C> - W1 and w2 are used to store sets of constants which are
19 C> reusable for repeated calls to the subroutine. 10 other arrays
20 C> are saved and reused on the next call.
21 C>
22 C> @author Ralph Jones @date 1994-05-03
23  SUBROUTINE w3ft16(ALOLA,BTHIN,INTERP)
24 C
25  parameter(npts=3447)
26 C
27  REAL SEP(73)
28  REAL ALOLA(95,91), BTHIN(NPTS), ERAS(NPTS,4)
29  REAL W1(NPTS), W2(NPTS)
30  REAL XDELI(NPTS), XDELJ(NPTS)
31  REAL XI2TM(NPTS), XJ2TM(NPTS)
32 C
33  INTEGER NPT(73)
34  INTEGER IV(NPTS), JV(NPTS), JY(NPTS,4)
35  INTEGER IM1(NPTS), IP1(NPTS), IP2(NPTS)
36 C
37  LOGICAL LIN
38 C
39  SAVE
40 C
41  DATA intrpo/99/
42  DATA iswt /0/
43 C
44 C GRID POINT SEPARATION
45 C
46  DATA sep /1.250, 1.250, 1.250, 1.250, 1.250, 1.250,
47  & 1.250, 1.250, 1.268, 1.268, 1.268, 1.286,
48  & 1.286, 1.286, 1.304, 1.304, 1.324, 1.324,
49  & 1.343, 1.364, 1.364, 1.385, 1.406, 1.406,
50  & 1.429, 1.452, 1.475, 1.500, 1.525, 1.525,
51  & 1.552, 1.579, 1.607, 1.636, 1.667, 1.698,
52  & 1.765, 1.800, 1.837, 1.875, 1.915, 1.957,
53  & 2.045, 2.093, 2.143, 2.195, 2.308, 2.368,
54  & 2.432, 2.571, 2.647, 2.813, 2.903, 3.103,
55  & 3.214, 3.333, 3.600, 3.750, 4.091, 4.286,
56  & 4.737, 5.000, 5.625, 6.000, 6.923, 8.182,
57  & 9.000,11.250,12.857,18.000,22.500,45.000,
58  & 90.000/
59 C
60 C NUMBER OF POINTS ALONG LAT CIRCLE FOR ONE OCTANT
61 C
62  DATA npt / 73, 73, 73, 73, 73, 73,
63  & 73, 73, 72, 72, 72, 71,
64  & 71, 71, 70, 70, 69, 69,
65  & 68, 67, 67, 66, 65, 65,
66  & 64, 63, 62, 61, 60, 60,
67  & 59, 58, 57, 56, 55, 54,
68  & 52, 51, 50, 49, 48, 47,
69  & 45, 44, 43, 42, 40, 39,
70  & 38, 36, 35, 33, 32, 30,
71  & 29, 28, 26, 25, 23, 22,
72  & 20, 19, 17, 16, 14, 12,
73  & 11, 9, 8, 6, 5, 3,
74  & 2/
75 C
76  lin = .false.
77  IF (interp.EQ.1) lin = .true.
78 C
79  IF (iswt.EQ.1) GO TO 900
80 C
81  ijout = 0
82  DO 200 j = 1,73
83  xjou = (j-1) * 1.25 + 1.0
84  ii = npt(j)
85  rdglat = sep(j)
86  DO 100 i = 1,ii
87  ijout = ijout + 1
88  w1(ijout) = (i-1) * rdglat + 3.0
89  w2(ijout) = xjou
90  100 CONTINUE
91  200 CONTINUE
92 C
93  iswt = 1
94  intrpo = interp
95  GO TO 1000
96 C
97 C AFTER THE 1ST CALL TO W3FT16 TEST INTERP, IF IT HAS
98 C CHANGED RECOMPUTE SOME CONSTANTS
99 C
100  900 CONTINUE
101  IF (interp.EQ.intrpo) GO TO 2100
102  intrpo = interp
103 C
104  1000 CONTINUE
105  DO 1100 k = 1,npts
106  iv(k) = w1(k)
107  jv(k) = w2(k)
108  xdeli(k) = w1(k) - iv(k)
109  xdelj(k) = w2(k) - jv(k)
110  ip1(k) = iv(k) + 1
111  jy(k,3) = jv(k) + 1
112  jy(k,2) = jv(k)
113  1100 CONTINUE
114 C
115  IF (lin) GO TO 1400
116 C
117  DO 1200 k = 1,npts
118  ip2(k) = iv(k) + 2
119  im1(k) = iv(k) - 1
120  jy(k,1) = jv(k) - 1
121  jy(k,4) = jv(k) + 2
122  xi2tm(k) = xdeli(k) * (xdeli(k) - 1.0) * .25
123  xj2tm(k) = xdelj(k) * (xdelj(k) - 1.0) * .25
124  1200 CONTINUE
125 C
126  1400 CONTINUE
127 C
128  IF (lin) GO TO 1700
129 C
130  DO 1500 kk = 1,npts
131  IF (jv(kk).LT.2.OR.jv(kk).GE.90) xj2tm(kk) = 0.0
132  1500 CONTINUE
133 C
134 C LINEAR INTERPOLATION
135 C
136  1700 CONTINUE
137  DO 1900 kk = 1,npts
138  IF (jy(kk,3).GT.91) jy(kk,3) = 91
139  1900 CONTINUE
140 C
141  IF (.NOT.lin) THEN
142  DO 2000 kk = 1,npts
143  IF (jy(kk,1).LT.1) jy(kk,1) = 1
144  IF (jy(kk,4).GT.91) jy(kk,4) = 91
145  2000 CONTINUE
146  ENDIF
147 C
148  2100 CONTINUE
149  IF (lin) THEN
150 C
151 C LINEAR INTERPOLATION
152 C
153  DO 2200 kk = 1,npts
154  eras(kk,2) = (alola(ip1(kk),jy(kk,2))-alola(iv(kk),jy(kk,2)))
155  & * xdeli(kk) + alola(iv(kk),jy(kk,2))
156  eras(kk,3) = (alola(ip1(kk),jy(kk,3))-alola(iv(kk),jy(kk,3)))
157  & * xdeli(kk) + alola(iv(kk),jy(kk,3))
158  2200 CONTINUE
159 C
160  DO 2300 kk = 1,npts
161  bthin(kk) = eras(kk,2) + (eras(kk,3) - eras(kk,2))
162  & * xdelj(kk)
163  2300 CONTINUE
164 C
165  ELSE
166 C
167 C QUADRATIC INTERPOLATION
168 C
169  DO 2400 kk = 1,npts
170  eras(kk,1)=(alola(ip1(kk),jy(kk,1))-alola(iv(kk),jy(kk,1)))
171  & * xdeli(kk) + alola(iv(kk),jy(kk,1)) +
172  & ( alola(im1(kk),jy(kk,1)) - alola(iv(kk),jy(kk,1))
173  & - alola(ip1(kk),jy(kk,1))+alola(ip2(kk),jy(kk,1)))
174  & * xi2tm(kk)
175  eras(kk,2)=(alola(ip1(kk),jy(kk,2))-alola(iv(kk),jy(kk,2)))
176  & * xdeli(kk) + alola(iv(kk),jy(kk,2)) +
177  & ( alola(im1(kk),jy(kk,2)) - alola(iv(kk),jy(kk,2))
178  & - alola(ip1(kk),jy(kk,2))+alola(ip2(kk),jy(kk,2)))
179  & * xi2tm(kk)
180  eras(kk,3)=(alola(ip1(kk),jy(kk,3))-alola(iv(kk),jy(kk,3)))
181  & * xdeli(kk) + alola(iv(kk),jy(kk,3)) +
182  & ( alola(im1(kk),jy(kk,3)) - alola(iv(kk),jy(kk,3))
183  & - alola(ip1(kk),jy(kk,3))+alola(ip2(kk),jy(kk,3)))
184  & * xi2tm(kk)
185  eras(kk,4)=(alola(ip1(kk),jy(kk,4))-alola(iv(kk),jy(kk,4)))
186  & * xdeli(kk) + alola(iv(kk),jy(kk,4)) +
187  & ( alola(im1(kk),jy(kk,4)) - alola(iv(kk),jy(kk,4))
188  & - alola(ip1(kk),jy(kk,4))+alola(ip2(kk),jy(kk,4)))
189  & * xi2tm(kk)
190  2400 CONTINUE
191 C
192  DO 2500 kk = 1,npts
193  bthin(kk) = eras(kk,2) + (eras(kk,3) - eras(kk,2))
194  & * xdelj(kk) + (eras(kk,1) - eras(kk,2)
195  & - eras(kk,3) + eras(kk,4)) * xj2tm(kk)
196  2500 CONTINUE
197 C
198  ENDIF
199 C
200  RETURN
201  END
w3ft16
subroutine w3ft16(ALOLA, BTHIN, INTERP)
Convert a northern hemisphere 1.0 degree lat.,lon.
Definition: w3ft16.f:24