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