NCEPLIBS-w3emc  2.11.0
w3ft207.f
Go to the documentation of this file.
1 C> @file
2 C> @brief Convert (361,91) grid to (49,35) n. hemi. grid
3 C> @author Ralph Jones @date 1993-10-19
4 
5 C> Convert a northern hemisphere 1.0 degree lat.,lon. 361 by
6 C> 91 grid to a polar stereographic 49 by 35 grid. The polar
7 C> stereographic map projection is true at 60 deg. n. , The mesh
8 C> length is 95.25 km. and the oriention is 150 deg. w.
9 C> awips grid 207 regional - alaska.
10 C>
11 C> ### Program History Log:
12 C> Date | Programmer | Comment
13 C> -----|------------|--------
14 C> 1993-10-19 | Ralph Jones | Initial.
15 C>
16 C> @param[in] ALOLA 361*91 grid 1.0 deg. lat,lon grid n. hemi.
17 C> 32851 point grid. 360 * 181 one degree grib grid 3 was flipped, greenwish added
18 C> to right side and cut to 361 * 91.
19 C> @param[in] INTERP 1 linear interpolation , ne.1 biquadratic
20 C> @param[out] APOLA 49*35 grid of northern hemisphere. 1715 point grid is
21 C> awips grid type 207
22 C>
23 C> @note
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 49*35 grid orientation
27 C> after interpolation. You may use w3fc08() to do this.
28 C>
29 C> @author Ralph Jones @date 1993-10-19
30  SUBROUTINE w3ft207(ALOLA,APOLA,INTERP)
31 C
32  parameter(npts=1715,ii=49,jj=35)
33  parameter(orient=150.0,ipole=25,jpole=51)
34  parameter(xmesh=95.250)
35 C
36  REAL R2(NPTS), WLON(NPTS)
37  REAL XLAT(NPTS), XI(II,JJ), XJ(II,JJ)
38  REAL XII(NPTS), XJJ(NPTS), ANGLE(NPTS)
39  REAL ALOLA(361,91), APOLA(NPTS), ERAS(NPTS,4)
40  REAL W1(NPTS), W2(NPTS)
41  REAL XDELI(NPTS), XDELJ(NPTS)
42  REAL XI2TM(NPTS), XJ2TM(NPTS)
43 C
44  INTEGER IV(NPTS), JV(NPTS), JY(NPTS,4)
45  INTEGER IM1(NPTS), IP1(NPTS), IP2(NPTS)
46 C
47  LOGICAL LIN
48 C
49  SAVE
50 C
51  equivalence(xi(1,1),xii(1)),(xj(1,1),xjj(1))
52 C
53  DATA degprd/57.2957795/
54  DATA earthr/6371.2/
55  DATA intrpo/99/
56  DATA iswt /0/
57 C
58  lin = .false.
59  IF (interp.EQ.1) lin = .true.
60 C
61  IF (iswt.EQ.1) GO TO 900
62 C
63  deg = 1.0
64  gi2 = (1.86603 * earthr) / xmesh
65  gi2 = gi2 * gi2
66 C
67 C NEXT 32 LINES OF CODE PUTS SUBROUTINE W3FB05 IN LINE
68 C
69  DO 100 j = 1,jj
70  xj1 = j - jpole
71  DO 100 i = 1,ii
72  xi(i,j) = i - ipole
73  xj(i,j) = xj1
74  100 CONTINUE
75 C
76  DO 200 kk = 1,npts
77  r2(kk) = xjj(kk) * xjj(kk) + xii(kk) * xii(kk)
78  xlat(kk) = degprd *
79  & asin((gi2 - r2(kk)) / (gi2 + r2(kk)))
80  200 CONTINUE
81 C
82  DO 300 kk = 1,npts
83  angle(kk) = degprd * atan2(xjj(kk),xii(kk))
84  300 CONTINUE
85 C
86  DO 400 kk = 1,npts
87  IF (angle(kk).LT.0.0) angle(kk) = angle(kk) + 360.0
88  400 CONTINUE
89 C
90  DO 500 kk = 1,npts
91  wlon(kk) = 270.0 + orient - angle(kk)
92  500 CONTINUE
93 C
94  DO 600 kk = 1,npts
95  IF (wlon(kk).LT.0.0) wlon(kk) = wlon(kk) + 360.0
96  600 CONTINUE
97 C
98  DO 700 kk = 1,npts
99  IF (wlon(kk).GE.360.0) wlon(kk) = wlon(kk) - 360.0
100  700 CONTINUE
101 C
102  DO 800 kk = 1,npts
103  w1(kk) = (360.0 - wlon(kk)) / deg + 1.0
104  w2(kk) = xlat(kk) / deg + 1.0
105  800 CONTINUE
106 C
107  iswt = 1
108  intrpo = interp
109  GO TO 1000
110 C
111 C AFTER THE 1ST CALL TO W3FT207 TEST INTERP, IF IT HAS
112 C CHANGED RECOMPUTE SOME CONSTANTS
113 C
114  900 CONTINUE
115  IF (interp.EQ.intrpo) GO TO 2100
116  intrpo = interp
117 C
118  1000 CONTINUE
119  DO 1100 k = 1,npts
120  iv(k) = w1(k)
121  jv(k) = w2(k)
122  xdeli(k) = w1(k) - iv(k)
123  xdelj(k) = w2(k) - jv(k)
124  ip1(k) = iv(k) + 1
125  jy(k,3) = jv(k) + 1
126  jy(k,2) = jv(k)
127  1100 CONTINUE
128 C
129  IF (lin) GO TO 1400
130 C
131  DO 1200 k = 1,npts
132  ip2(k) = iv(k) + 2
133  im1(k) = iv(k) - 1
134  jy(k,1) = jv(k) - 1
135  jy(k,4) = jv(k) + 2
136  xi2tm(k) = xdeli(k) * (xdeli(k) - 1.0) * .25
137  xj2tm(k) = xdelj(k) * (xdelj(k) - 1.0) * .25
138  1200 CONTINUE
139 C
140  DO 1300 kk = 1,npts
141  IF (iv(kk).EQ.1) THEN
142  ip2(kk) = 3
143  im1(kk) = 360
144  ELSE IF (iv(kk).EQ.360) THEN
145  ip2(kk) = 2
146  im1(kk) = 359
147  ENDIF
148  1300 CONTINUE
149 C
150  1400 CONTINUE
151 C
152  IF (lin) GO TO 1700
153 C
154  DO 1500 kk = 1,npts
155  IF (jv(kk).LT.2.OR.jv(kk).GT.89) xj2tm(kk) = 0.0
156  1500 CONTINUE
157 C
158  DO 1600 kk = 1,npts
159  IF (ip2(kk).LT.1) ip2(kk) = 1
160  IF (im1(kk).LT.1) im1(kk) = 1
161  IF (ip2(kk).GT.361) ip2(kk) = 361
162  IF (im1(kk).GT.361) im1(kk) = 361
163  1600 CONTINUE
164 C
165  1700 CONTINUE
166  DO 1800 kk = 1,npts
167  IF (iv(kk).LT.1) iv(kk) = 1
168  IF (ip1(kk).LT.1) ip1(kk) = 1
169  IF (iv(kk).GT.361) iv(kk) = 361
170  IF (ip1(kk).GT.361) ip1(kk) = 361
171  1800 CONTINUE
172 C
173 C LINEAR INTERPOLATION
174 C
175  DO 1900 kk = 1,npts
176  IF (jy(kk,2).LT.1) jy(kk,2) = 1
177  IF (jy(kk,2).GT.91) jy(kk,2) = 91
178  IF (jy(kk,3).LT.1) jy(kk,3) = 1
179  IF (jy(kk,3).GT.91) jy(kk,3) = 91
180  1900 CONTINUE
181 C
182  IF (.NOT.lin) THEN
183  DO 2000 kk = 1,npts
184  IF (jy(kk,1).LT.1) jy(kk,1) = 1
185  IF (jy(kk,1).GT.91) jy(kk,1) = 91
186  IF (jy(kk,4).LT.1) jy(kk,4) = 1
187  IF (jy(kk,4).GT.91) jy(kk,4) = 91
188  2000 CONTINUE
189  ENDIF
190 C
191  2100 CONTINUE
192  IF (lin) THEN
193 C
194 C LINEAR INTERPOLATION
195 C
196  DO 2200 kk = 1,npts
197  eras(kk,2) = (alola(ip1(kk),jy(kk,2))-alola(iv(kk),jy(kk,2)))
198  & * xdeli(kk) + alola(iv(kk),jy(kk,2))
199  eras(kk,3) = (alola(ip1(kk),jy(kk,3))-alola(iv(kk),jy(kk,3)))
200  & * xdeli(kk) + alola(iv(kk),jy(kk,3))
201  2200 CONTINUE
202 C
203  DO 2300 kk = 1,npts
204  apola(kk) = eras(kk,2) + (eras(kk,3) - eras(kk,2))
205  & * xdelj(kk)
206  2300 CONTINUE
207 C
208  ELSE
209 C
210 C QUADRATIC INTERPOLATION
211 C
212  DO 2400 kk = 1,npts
213  eras(kk,1)=(alola(ip1(kk),jy(kk,1))-alola(iv(kk),jy(kk,1)))
214  & * xdeli(kk) + alola(iv(kk),jy(kk,1)) +
215  & ( alola(im1(kk),jy(kk,1)) - alola(iv(kk),jy(kk,1))
216  & - alola(ip1(kk),jy(kk,1))+alola(ip2(kk),jy(kk,1)))
217  & * xi2tm(kk)
218  eras(kk,2)=(alola(ip1(kk),jy(kk,2))-alola(iv(kk),jy(kk,2)))
219  & * xdeli(kk) + alola(iv(kk),jy(kk,2)) +
220  & ( alola(im1(kk),jy(kk,2)) - alola(iv(kk),jy(kk,2))
221  & - alola(ip1(kk),jy(kk,2))+alola(ip2(kk),jy(kk,2)))
222  & * xi2tm(kk)
223  eras(kk,3)=(alola(ip1(kk),jy(kk,3))-alola(iv(kk),jy(kk,3)))
224  & * xdeli(kk) + alola(iv(kk),jy(kk,3)) +
225  & ( alola(im1(kk),jy(kk,3)) - alola(iv(kk),jy(kk,3))
226  & - alola(ip1(kk),jy(kk,3))+alola(ip2(kk),jy(kk,3)))
227  & * xi2tm(kk)
228  eras(kk,4)=(alola(ip1(kk),jy(kk,4))-alola(iv(kk),jy(kk,4)))
229  & * xdeli(kk) + alola(iv(kk),jy(kk,4)) +
230  & ( alola(im1(kk),jy(kk,4)) - alola(iv(kk),jy(kk,4))
231  & - alola(ip1(kk),jy(kk,4))+alola(ip2(kk),jy(kk,4)))
232  & * xi2tm(kk)
233  2400 CONTINUE
234 C
235  DO 2500 kk = 1,npts
236  apola(kk) = eras(kk,2) + (eras(kk,3) - eras(kk,2))
237  & * xdelj(kk) + (eras(kk,1) - eras(kk,2)
238  & - eras(kk,3) + eras(kk,4)) * xj2tm(kk)
239  2500 CONTINUE
240 C
241  ENDIF
242 C
243  RETURN
244  END
w3ft207
subroutine w3ft207(ALOLA, APOLA, INTERP)
Convert a northern hemisphere 1.0 degree lat.,lon.
Definition: w3ft207.f:31