NCEPLIBS-w3emc  2.11.0
w3ft05v.f
Go to the documentation of this file.
1 C> @file
2 C> @brief Convert (145,37) grid to (65,65) n. hemi. grid
3 C> @author Ralph Jones @date 1985-04-10
4 
5 C> Convert a northern hemisphere 2.5 degree lat.,lon. 145 by
6 C> 37 grid to a polar stereographic 65 by 65 grid. The polar
7 C> stereographic map projection is true at 60 deg. n. , The mesh
8 C> length is 381 km. and the oriention is 80 deg. w.
9 C>
10 C> ### Program History Log:
11 C> Date | Programmer | Comment
12 C> -----|------------|--------
13 C> 1985-04-10 | Ralph Jones | Vectorized version of w3ft05().
14 C> 1989-10-21 | Ralph Jones | Changes to increase speed.
15 C> 1991-07-25 | Ralph Jones | Change to cray cft77 fortran.
16 C>
17 C> @param[in] ALOLA 145*37 gid 2.5 lat,lon grid n. hemisphere
18 C> 5365 point grid is o.n. 84 type 29 or 1d hex
19 C> interp - 1 linear interpolation , ne.1 biquadratic
20 C> @param[out] APOLA 65*65 grid of northern hemisphere
21 C> 4225 point grid is o.n.84 type 27 or 1b hex.
22 C> @param INTERP
23 C> @remark
24 C> - 1. W1 and w2 are used to store sets of constants which are
25 C> reusable for repeated calls to the subroutine.
26 C> - 2. Wind components are not rotated to the 65*65 grid orientation
27 C> after interpolation. you may use w3fc08 to do this.
28 C> - 3. The grid points values on the equator have been extrapolated
29 C> outward to all the grid points outside the equator on the 65*65
30 C> grid (about 1100 points).
31 C>
32 C> @author Ralph Jones @date 1985-04-10
33  SUBROUTINE w3ft05v(ALOLA,APOLA,INTERP)
34 C
35  REAL R2(4225), WLON(4225)
36  REAL XLAT(4225), XI(65,65), XJ(65,65)
37  REAL XII(4225), XJJ(4225), ANGLE(4225)
38  REAL ALOLA(145,37), APOLA(4225), ERAS(4225,4)
39  REAL W1(4225), W2(4225)
40  REAL XDELI(4225), XDELJ(4225)
41  REAL XI2TM(4225), XJ2TM(4225)
42 C
43  INTEGER IV(4225), JV(4225), JY(4225,4)
44  INTEGER IM1(4225), IP1(4225), IP2(4225)
45 C
46  LOGICAL LIN
47 C
48  SAVE
49 C
50  equivalence(xi(1,1),xii(1)),(xj(1,1),xjj(1))
51 C
52  DATA degprd/57.2957795/
53  DATA earthr/6371.2/
54  DATA intrpo/99/
55  DATA iswt /0/
56 C
57  lin = .false.
58  IF (interp.EQ.1) lin = .true.
59 C
60  IF (iswt.EQ.1) GO TO 900
61 C
62  orient = 80.0
63  deg = 2.5
64  xmesh = 381.0
65  gi2 = (1.86603 * earthr) / xmesh
66  gi2 = gi2 * gi2
67 C
68 C NEXT 32 LINES OF CODE PUTS SUBROUTINE W3FB01 IN LINE
69 C
70  DO 100 j = 1,65
71  xj1 = j - 33
72  DO 100 i = 1,65
73  xi(i,j) = i - 33
74  xj(i,j) = xj1
75  100 CONTINUE
76 C
77  DO 200 kk = 1,4225
78  r2(kk) = xjj(kk) * xjj(kk) + xii(kk) * xii(kk)
79  xlat(kk) = degprd *
80  & asin((gi2 - r2(kk)) / (gi2 + r2(kk)))
81  200 CONTINUE
82 C
83  xii(2113) = 1.0
84  DO 300 kk = 1,4225
85  angle(kk) = degprd * atan2(xjj(kk),xii(kk))
86  300 CONTINUE
87 C
88  DO 400 kk = 1,4225
89  IF (angle(kk).LT.0.0) angle(kk) = angle(kk) + 360.0
90  400 CONTINUE
91 C
92  DO 500 kk = 1,4225
93  wlon(kk) = 270.0 + orient - angle(kk)
94  500 CONTINUE
95 C
96  DO 600 kk = 1,4225
97  IF (wlon(kk).LT.0.0) wlon(kk) = wlon(kk) + 360.0
98  600 CONTINUE
99 C
100  DO 700 kk = 1,4225
101  IF (wlon(kk).GE.360.0) wlon(kk) = wlon(kk) - 360.0
102  700 CONTINUE
103 C
104  xlat(2113) = 90.0
105  wlon(2113) = 0.0
106 C
107  DO 800 kk = 1,4225
108  w1(kk) = (360.0 - wlon(kk)) / deg + 1.0
109  w2(kk) = xlat(kk) / deg + 1.0
110  800 CONTINUE
111 C
112  iswt = 1
113  intrpo = interp
114  GO TO 1000
115 C
116 C AFTER THE 1ST CALL TO W3FT05V TEST INTERP, IF IT HAS
117 C CHANGED RECOMPUTE SOME CONSTANTS
118 C
119  900 CONTINUE
120  IF (interp.EQ.intrpo) GO TO 2100
121  intrpo = interp
122 C
123  1000 CONTINUE
124  DO 1100 k = 1,4225
125  iv(k) = w1(k)
126  jv(k) = w2(k)
127  xdeli(k) = w1(k) - iv(k)
128  xdelj(k) = w2(k) - jv(k)
129  ip1(k) = iv(k) + 1
130  jy(k,3) = jv(k) + 1
131  jy(k,2) = jv(k)
132  1100 CONTINUE
133 C
134  IF (lin) GO TO 1400
135 C
136  DO 1200 k = 1,4225
137  ip2(k) = iv(k) + 2
138  im1(k) = iv(k) - 1
139  jy(k,1) = jv(k) - 1
140  jy(k,4) = jv(k) + 2
141  xi2tm(k) = xdeli(k) * (xdeli(k) - 1.0) * .25
142  xj2tm(k) = xdelj(k) * (xdelj(k) - 1.0) * .25
143  1200 CONTINUE
144 C
145  DO 1300 kk = 1,4225
146  IF (iv(kk).EQ.1) THEN
147  ip2(kk) = 3
148  im1(kk) = 144
149  ELSE IF (iv(kk).EQ.144) THEN
150  ip2(kk) = 2
151  im1(kk) = 143
152  ENDIF
153  1300 CONTINUE
154 C
155  1400 CONTINUE
156 C
157  IF (lin) GO TO 1700
158 C
159  DO 1500 kk = 1,4225
160  IF (jv(kk).LT.2.OR.jv(kk).GT.35) xj2tm(kk) = 0.0
161  1500 CONTINUE
162 C
163  DO 1600 kk = 1,4225
164  IF (ip2(kk).LT.1) ip2(kk) = 1
165  IF (im1(kk).LT.1) im1(kk) = 1
166  IF (ip2(kk).GT.145) ip2(kk) = 145
167  IF (im1(kk).GT.145) im1(kk) = 145
168  1600 CONTINUE
169 C
170  1700 CONTINUE
171  DO 1800 kk = 1,4225
172  IF (iv(kk).LT.1) iv(kk) = 1
173  IF (ip1(kk).LT.1) ip1(kk) = 1
174  IF (iv(kk).GT.145) iv(kk) = 145
175  IF (ip1(kk).GT.145) ip1(kk) = 145
176  1800 CONTINUE
177 C
178 C LINEAR INTERPOLATION
179 C
180  DO 1900 kk = 1,4225
181  IF (jy(kk,2).LT.1) jy(kk,2) = 1
182  IF (jy(kk,2).GT.37) jy(kk,2) = 37
183  IF (jy(kk,3).LT.1) jy(kk,3) = 1
184  IF (jy(kk,3).GT.37) jy(kk,3) = 37
185  1900 CONTINUE
186 C
187  IF (.NOT.lin) THEN
188  DO 2000 kk = 1,4225
189  IF (jy(kk,1).LT.1) jy(kk,1) = 1
190  IF (jy(kk,1).GT.37) jy(kk,1) = 37
191  IF (jy(kk,4).LT.1) jy(kk,4) = 1
192  IF (jy(kk,4).GT.37) jy(kk,4) = 37
193  2000 CONTINUE
194  ENDIF
195 C
196  2100 CONTINUE
197  IF (lin) THEN
198 C
199 C LINEAR INTERPOLATION
200 C
201  DO 2200 kk = 1,4225
202  eras(kk,2) = (alola(ip1(kk),jy(kk,2))-alola(iv(kk),jy(kk,2)))
203  & * xdeli(kk) + alola(iv(kk),jy(kk,2))
204  eras(kk,3) = (alola(ip1(kk),jy(kk,3))-alola(iv(kk),jy(kk,3)))
205  & * xdeli(kk) + alola(iv(kk),jy(kk,3))
206  2200 CONTINUE
207 C
208  DO 2300 kk = 1,4225
209  apola(kk) = eras(kk,2) + (eras(kk,3) - eras(kk,2))
210  & * xdelj(kk)
211  2300 CONTINUE
212 C
213  ELSE
214 C
215 C QUADRATIC INTERPOLATION
216 C
217  DO 2400 kk = 1,4225
218  eras(kk,1)=(alola(ip1(kk),jy(kk,1))-alola(iv(kk),jy(kk,1)))
219  & * xdeli(kk) + alola(iv(kk),jy(kk,1)) +
220  & ( alola(im1(kk),jy(kk,1)) - alola(iv(kk),jy(kk,1))
221  & - alola(ip1(kk),jy(kk,1))+alola(ip2(kk),jy(kk,1)))
222  & * xi2tm(kk)
223  eras(kk,2)=(alola(ip1(kk),jy(kk,2))-alola(iv(kk),jy(kk,2)))
224  & * xdeli(kk) + alola(iv(kk),jy(kk,2)) +
225  & ( alola(im1(kk),jy(kk,2)) - alola(iv(kk),jy(kk,2))
226  & - alola(ip1(kk),jy(kk,2))+alola(ip2(kk),jy(kk,2)))
227  & * xi2tm(kk)
228  eras(kk,3)=(alola(ip1(kk),jy(kk,3))-alola(iv(kk),jy(kk,3)))
229  & * xdeli(kk) + alola(iv(kk),jy(kk,3)) +
230  & ( alola(im1(kk),jy(kk,3)) - alola(iv(kk),jy(kk,3))
231  & - alola(ip1(kk),jy(kk,3))+alola(ip2(kk),jy(kk,3)))
232  & * xi2tm(kk)
233  eras(kk,4)=(alola(ip1(kk),jy(kk,4))-alola(iv(kk),jy(kk,4)))
234  & * xdeli(kk) + alola(iv(kk),jy(kk,4)) +
235  & ( alola(im1(kk),jy(kk,4)) - alola(iv(kk),jy(kk,4))
236  & - alola(ip1(kk),jy(kk,4))+alola(ip2(kk),jy(kk,4)))
237  & * xi2tm(kk)
238  2400 CONTINUE
239 C
240  DO 2500 kk = 1,4225
241  apola(kk) = eras(kk,2) + (eras(kk,3) - eras(kk,2))
242  & * xdelj(kk) + (eras(kk,1) - eras(kk,2)
243  & - eras(kk,3) + eras(kk,4)) * xj2tm(kk)
244  2500 CONTINUE
245 C
246  ENDIF
247 C
248 C SET POLE POINT , WMO STANDARD FOR U OR V
249 C
250  apola(2113) = alola(73,37)
251 C
252  RETURN
253  END
w3ft05v
subroutine w3ft05v(ALOLA, APOLA, INTERP)
Convert a northern hemisphere 2.5 degree lat.,lon.
Definition: w3ft05v.f:34