NCEPLIBS-w3emc 2.12.0
Loading...
Searching...
No Matches
w3ft05.f
Go to the documentation of this file.
1C> @file
2C> @brief Convert (145,37) to (65,65) n. hemi. grid
3C> @author Ralph Jones @date 1985-04-08
4
5C> Convert a northern hemisphere 2.5 degree lat.,lon. 145 by
6C> 37 grid to a polar stereographic 65 by 65 grid. The polar
7C> stereographic map projection is true at 60 deg. n. , The mesh
8C> length is 381 km. and the oriention is 80 deg. w.
9C>
10C> ### Program History Log:
11C> Date | Programmer | Comment
12C> -----|------------|--------
13C> 1985-04-08 | Ralph Jones | Initial.
14C> 1991-07-30 | Ralph Jones | convert to cray cft77 fortran.
15C> 1992-05-02 | Ralph Jones | add save.
16C>
17C> @param[in] ALOLA 145*37 grid 2.5 lat,lon grid n. hemi.
18C> 5365 point grid is type 29 or 1d hex o.n. 84
19C> @param[in] LINEAR 1 linear interpolation , ne.1 biquadratic
20C> @param[out] APOLA 65*65 grid of northern hemi.
21C> 4225 point grid is type 27 or 1b hex o.n. 84
22C> @param[out] W1 65*65 scratch field
23C> @param[out] W2 65*65 scratch field
24C>
25C> @remark
26C> - 1. W1 and w2 are used to store sets of constants which are
27C> reusable for repeated calls to the subroutine. If they are
28C> over written by the user, a warning message will be printed
29C> and w1 and w2 will be recomputed.
30C> - 2. Wind components are not rotated to the 65*65 grid orientation
31C> after interpolation. You may use w3fc08() to do this.
32C> - 3. The grid points values on the equator have been extrapolated
33C> outward to all the grid points outside the equator on the 65*65
34C> grid (about 1100 points).
35C> - 4. You should use the cray vectorized version w3ft05v on the cray
36C> it has 3 parameters in the call, runs about 10 times faster. Uses
37C> more memory.
38C>
39C> @author Ralph Jones @date 1985-04-08
40 SUBROUTINE w3ft05(ALOLA,APOLA,W1,W2,LINEAR)
41C
42 REAL ALOLA(145,37)
43 REAL APOLA(4225)
44 REAL ERAS(4)
45 REAL SAVEW1(10)
46 REAL SAVEW2(10)
47 REAL W1(4225)
48 REAL W2(4225)
49C
50 INTEGER JY(4)
51 INTEGER OUT
52C
53 LOGICAL LIN
54C
55 SAVE
56C
57 DATA degprd/57.2957795/
58 DATA earthr/6371.2/
59 DATA iswt /0/
60 DATA out /6/
61C
62 4000 FORMAT ( 52h *** warning , w1 or w2 scratch files over written ,,
63 & 43h i will restore them , burning up cpu time,,
64 & 14h in w3ft05 ***)
65C
66 lin = .false.
67 IF (linear.EQ.1) lin = .true.
68C
69 IF (iswt.EQ.0) GO TO 300
70C
71C TEST W1 AND W2 TO SEE IF THEY WERE WRITTEN OVER
72C
73 DO 100 kk=1,10
74 IF (savew1(kk).NE.w1(kk)) GO TO 200
75 IF (savew2(kk).NE.w2(kk)) GO TO 200
76 100 CONTINUE
77 GOTO 1000
78C
79 200 CONTINUE
80 WRITE (out,4000)
81C
82 300 CONTINUE
83 deg = 2.5
84 nn = 0
85 xmesh = 381.0
86 gi2 = (1.86603*earthr) / xmesh
87 gi2 = gi2 * gi2
88C
89C DO LOOP 800 PUTS SUBROUTINE W3FB01 IN LINE
90C
91 DO 800 j = 1,65
92 xj = j - 33
93 xj2 = xj * xj
94 DO 800 i=1,65
95 xi = i - 33
96 r2 = xi*xi + xj2
97 IF (r2.NE.0.0) GO TO 400
98 wlon = 0.0
99 xlat = 90.0
100 GO TO 700
101 400 CONTINUE
102 xlong = degprd * atan2(xj,xi)
103 IF (xlong.GE.0.0) GO TO 500
104 wlon = -10.0 - xlong
105 IF (wlon.LT.0.0) wlon = wlon + 360.0
106 GO TO 600
107C
108 500 CONTINUE
109 wlon = 350.0 - xlong
110 600 CONTINUE
111 xlat = asin((gi2-r2)/(gi2+r2))*degprd
112 700 CONTINUE
113 IF (wlon.GT.360.0) wlon = wlon - 360.0
114 IF (wlon.LT.0.0) wlon = wlon + 360.0
115 nn = nn + 1
116 w1(nn) = ( 360.0 - wlon ) / deg + 1.0
117 w2(nn) = xlat / deg + 1.0
118 800 CONTINUE
119C
120 DO 900 kk = 1,10
121 savew1(kk) = w1(kk)
122 savew2(kk) = w2(kk)
123 900 CONTINUE
124C
125 iswt = 1
126C
127 1000 CONTINUE
128C
129 DO 2100 kk = 1,4225
130 i = w1(kk)
131 j = w2(kk)
132 fi = i
133 fj = j
134 xdeli = w1(kk) - fi
135 xdelj = w2(kk) - fj
136 ip1 = i + 1
137 jy(3) = j + 1
138 jy(2) = j
139 IF (lin) GO TO 1100
140 ip2 = i + 2
141 im1 = i - 1
142 jy(4) = j + 2
143 jy(1) = j - 1
144 xi2tm = xdeli * (xdeli-1.) * 0.25
145 xj2tm = xdelj * (xdelj-1.) * 0.25
146C
147 1100 CONTINUE
148 IF ((i.LT.2).OR.(j.LT.2)) GO TO 1200
149 IF ((i.GT.142).OR.(j.GT.34)) GO TO 1200
150C
151C QUADRATIC (LINEAR TOO) OK W/O FURTHER ADO SO GO TO 1700
152C
153 GO TO 1700
154C
155 1200 CONTINUE
156 IF (i.EQ.1) GO TO 1300
157 IF (i.EQ.144) GO TO 1400
158 ip2 = i + 2
159 im1 = i - 1
160 GO TO 1500
161C
162 1300 CONTINUE
163 ip2 = 3
164 im1 = 144
165 GO TO 1500
166C
167 1400 CONTINUE
168 ip2 = 2
169 im1 = 143
170C
171 1500 CONTINUE
172 ip1 = i + 1
173 IF (lin) GO TO 1600
174 IF ((j.LT.2).OR.(j.GE.36)) xj2tm=0.
175C.....DO NOT ALLOW POINT OFF GRID
176 IF (ip2.LT.1) ip2 = 1
177 IF (im1.LT.1) im1 = 1
178 IF (ip2.GT.145) ip2 = 145
179 IF (im1.GT.145) im1 = 145
180C
181 1600 CONTINUE
182C.....DO NOT ALLOW POINT OFF GRID
183 IF (i.LT.1) i = 1
184 IF (ip1.LT.1) ip1 = 1
185 IF (i.GT.145) i = 145
186 IF (ip1.GT.145) ip1 = 145
187C
188 1700 CONTINUE
189 IF (.NOT.lin) GO TO 1900
190C
191C LINEAR INTERPLOATION
192C
193 DO 1800 k = 2,3
194 j1 = jy(k)
195 IF (j1.LT.1) j1 = 1
196 IF (j1.GT.37) j1 = 37
197 eras(k) = (alola(ip1,j1) - alola(i,j1)) * xdeli + alola(i,j1)
198 1800 CONTINUE
199C
200 apola(kk) = eras(2) + (eras(3) - eras(2)) * xdelj
201 GO TO 2100
202C
203 1900 CONTINUE
204C
205C QUADRATIC INTERPOLATION
206C
207 DO 2000 k = 1,4
208 j1 = jy(k)
209C.....DO NOT ALLOW POINT OFF GRID
210 IF (j1.LT.1) j1 = 1
211 IF (j1.GT.37) j1 = 37
212 eras(k) = (alola(ip1,j1)-alola(i,j1))*xdeli+alola(i,j1)+
213 & (alola(im1,j1)-alola(i,j1)-alola(ip1,j1)+
214 & alola(ip2,j1))*xi2tm
215 2000 CONTINUE
216C
217 apola(kk) = eras(2)+(eras(3)-eras(2))*xdelj+(eras(1)-
218 & eras(2)-eras(3)+eras(4)) * xj2tm
219C
220 2100 CONTINUE
221C
222C SET POLE POINT , WMO STANDARD FOR U OR V
223C
224 apola(2113) = alola(73,37)
225C
226 RETURN
227 END
w3ft05
subroutine w3ft05(alola, apola, w1, w2, linear)
Convert a northern hemisphere 2.5 degree lat.,lon.
Definition w3ft05.f:41