GNU WOCSS (gwocss)  2.2.4-pre
GNU version of Winds On Critical Streamline Surfaces (WOCSS)
Functions/Subroutines
levwnd.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine levwnd
 

Function/Subroutine Documentation

subroutine levwnd ( )

Definition at line 3 of file levwnd.f.

Referenced by tstwnd().

3 C*********************************************************************
4 C
5 C INTERPOLATES MASS ADJUSTED FIELD TO ANEMOMETER HEIGHT
6 C (Z10, FOR Z INDEX=1)
7 C & TO NHORIZ-1 FLAT PLANES (FOR Z INDICES 2-NHORIZ) SET AT THE
8 C HEIGHTS ZCHOOZ ABOVE THE LOWEST TERRRAIN GRID POINT IN THE COARSE
9 C GRID, AS DEFINEDTHE VARIABLE ARRAY SIGMA. THIS SUBROUTINE ALSO
10 C CONVERTS INTERPOLATED WINDS BACK TO METEOROLOGICAL SPEEDS
11 C AND ANGLES THAT CAN BE PLOTTED IF DESIRED. LUDWIG--JANUARY 1988
12 C
13 C MODIFIED APRIL 2002 TO INTERPOLATE OTHER VARIABLES (RICHARDSON NO.,
14 C POT. TEMP. LAPSE RATE, PRESSURE, POT. TEMP, BRUNT VAISALA PERIOD)
15 C TO FLAT SURFACES AS WELL. ALSO SIMPLIFIED LOGIC.
16 C
17 C FLUDWIG 4/02
18 C
19  parameter(rad2d=180./3.14159)
20 C
21  include 'NGRIDS.PAR'
22 C
23  include 'ANCHOR.CMM'
24  include 'FLOWER.CMM'
25  include 'LIMITS.CMM'
26  include 'TSONDS.CMM'
27 C
28  LOGICAL ifxpt(nxgrd,nygrd),dotrpo,didfnd
29 C
30 C MAKE SURE EVERYTHING IS STILL HERE ON SUBSEQUENT CALLS
31 C
32  SAVE
33 C
34  k1=1
35  z0=zzero
36  z10=10.0
37  nbada=0
38 C
39  CALL fixwnd (ifxpt,1)
40 C
41  DO 400 ix=1,ncol
42  DO 390 iy=1,nrow
43 C
44 C 1ST LEVEL REPRESENTS A TERRAIN FOLLOWING SURFACE, Z10 M ABOVE SFC
45 C EVERYWHERE. IN ORDER TO REFLECT BOTH THE OBSERVATIONS USED TO
46 C GET FIRST ESTIMATES OF ANEMOMETER LEVEL WINDS AND THE ADJUSTED
47 C WINDS ON THE FIRST ABOVE GROUND FLOW SURFACE, WE USE FOR THE LOWER
48 C END OF THE INTERPOLATION, HALF THE INITIAL LEVEL 1 WIND ESTIMATE AT
49 C THE HT WHERE A LOG PROFILE WOULD HAVE THE WINDS BE 50% OF THEIR
50 C MAGNITUDE AT Z=Z10. VALUES AT THE 1ST LEVEL ABOVE GROUND ARE USED
51 C FOR THE UPPER END OF THE INTERPOLATION. FOR PTS NEAR OBS, THE
52 C WEIGHTING OF ANEMOMETER HT VALUES IS INCREASED BY USING Z FOR
53 C 90% OF ANEMOMETER HT VLAUE. LAPSE RATE UNRELIABLE AT THIS LEVEL, SO
54 C BULK RICHARDSON NUMBER AND BRUNT-VAISALA PERIOD NOT CALCULATED
55 C
56  zagl=z10
57  IF (ifxpt(ix,iy)) THEN
58  bottom=z0*((z10/z0)**0.9)
59  u0=0.9*u(ix,iy,1)
60  v0=0.9*v(ix,iy,1)
61  w0=0.9*w(ix,iy,1)
62  ELSE
63  bottom=z0*((z10/z0)**0.5)
64  u0=0.5*u(ix,iy,1)
65  v0=0.5*v(ix,iy,1)
66  w0=0.5*w(ix,iy,1)
67  END IF
68 C
69 C NOW GET VALUES TO USE FOR TOP END OF INTERPOLATION.
70 C
71  top=rhs(ix,iy,levbot(ix,iy))
72  u1=u(ix,iy,levbot(ix,iy))
73  v1=v(ix,iy,levbot(ix,iy))
74  w1=w(ix,iy,levbot(ix,iy))
75 C
76 C WE NOW HAVE VALUES ABOVE AND BELOW THE ANEMOMETER HEIGHT --
77 C INTERPOLATE
78 C
79  CALL lgntrp(uuuf,bottom,top,zagl,u0,u1)
80  u(ix,iy,1)=uuuf
81  CALL lgntrp(vvvf,bottom,top,zagl,v0,v1)
82  v(ix,iy,1)=vvvf
83  CALL lgntrp(wwwf,bottom,top,zagl,w0,w1)
84  w(ix,iy,1)=wwwf
85  iugraf(ix,iy,1)=nint(100.0*u(ix,iy,1))
86  ivgraf(ix,iy,1)=nint(100.0*v(ix,iy,1))
87  iwgraf(ix,iy,1)=nint(100.0*w(ix,iy,1))
88  iptgrf(ix,iy,1)=nint(10.0*tmpkel(ix,iy,1))
89 C IPRGRF(IX,IY,1)=NINT(100.0*SETALT(IX,IY,1))
90  patsfc=cvt2p(setalt(ix,iy,1),sfcht(ix,iy)+z10)
91  iprgrf(ix,iy,1)=nint(10.0*patsfc)
92 C
93 C DONE WITH ANEMOMETER HEIGHT VALUES, NOW INTERPOLATE TO HORIZONTAL
94 C PLANES (IZ>.
95 C
96  DO 300 iz=2,nflat
97 C
98 C GET HEIGHT OF HORIZONTAL PLANE ABOVE LOCAL GROUND LEVEL (ZAGL), THEN
99 C FIND 1ST FLOW SURFACE ABOVE THIS FLAT SURFACE, MEANWHILE FILLING
100 C BELOW SFC LEVELS WITH MISSING DATA VALUES. IF ALL SURFACES ARE
101 C LOWER THAN TERRAIN AT THIS GRID POINT SET ALL LEVELS TO MISSING
102 C VALUES AND PROCEED TO NEXT GRID PT.
103 C
104  IF (doflat) THEN
105  zagl=zchooz(iz)-sfcht(ix,iy)
106  zmsl=zchooz(iz)
107  ELSE
108  zagl=zchooz(iz)
109  zmsl=zchooz(iz)+sfcht(ix,iy)
110  END IF
111 C
112  IF (zagl.LE.z0) THEN
113  dotrpo=.false.
114 C
115 C THIS SURFACE BELOW GROUND; SET VALUES TO ZERO OR MISSING
116 C
117  iugraf(ix,iy,iz)=0
118  ivgraf(ix,iy,iz)=0
119  iwgraf(ix,iy,iz)=0
120  iptgrf(ix,iy,iz)=0
121  iprgrf(ix,iy,iz)=0
122 C
123  ELSE IF (zagl .GE. rhs(ix,iy,nlvl)) THEN
124  dotrpo=.false.
125 C
126 C EXRAPOLATION SURFACE ABOVE HIGHEST FLOW SURFACE. USE VALUES
127 C BASED ON THOSE ON HIGHEST SURFACE. ASSUME A LAPSE RATE OF
128 C
129  iugraf(ix,iy,iz)=nint(100.0*u(ix,iy,nlvl))
130  ivgraf(ix,iy,iz)=nint(100.0*v(ix,iy,nlvl))
131  iwgraf(ix,iy,iz)=nint(100.0*w(ix,iy,nlvl))
132  IF (dothet .OR. dopres .OR. dobri .OR. dobvpd) THEN
133  ttt=tmpkel(ix,iy,nlvl)+dthdzl(ix,iy,nlvl)*
134  $ (zagl-rhs(ix,iy,nlvl))
135  iptgrf(ix,iy,iz)=nint(10.0*ttt)
136  atop=setalt(ix,iy,nlvl)
137 C IPRGRF(IX,IY,IZ)=NINT(100.0*ATOP)
138  iprgrf(ix,iy,iz)=nint(10.0*cvt2p(atop,zmsl))
139  END IF
140 C
141  ELSE IF (zagl.LE. z10) THEN
142  dotrpo=.true.
143 C
144 C INTERPOLATION SURFACE IS BELOW ANEMOMETER HEIGHT. USE LOG PROFILE
145 C WIND VALUES FROM ANEMOMETER LEVEL. USE ANEMOMETER LEVEL VALUES FOR
146 C OTHER PARAMETERS.
147 C
148  ibr=4
149  bottom=z0
150  u0=0.0
151  v0=0.0
152  w0=0.0
153  IF (dothet .OR. dopres .OR. dobri .OR. dobvpd) THEN
154  t0=tmpkel(ix,iy,1)
155  t1=t0
156 C
157 C ASSUME AT LOWEST SURFACE APPLIES HERE.
158 C
159  alt0=setalt(ix,iy,levbot(ix,iy))
160 C IPRGRF(IX,IY,IZ)=NINT(100.0*ALT0)
161  iprgrf(ix,iy,iz)=nint(10.0*cvt2p(alt0,zmsl))
162  END IF
163 C
164 C NOW GET VALUES TO USE FOR TOP END OF INTERPOLATION. ASSUME NO
165 C CHANGE AND USE SAME VALUES FOR THE VARIABLES OTHER THAN WIND.
166 C
167  top=z10
168  u1=u(ix,iy,1)
169  v1=v(ix,iy,1)
170  w1=w(ix,iy,1)
171 C
172  ELSE IF (zagl .LT. rhs(ix,iy,levbot(ix,iy))
173  $ .AND. zagl.GT.z10) THEN
174  dotrpo=.true.
175 C
176 C THIS INTERPOLATION SFC MUST BE BETWEEN ANEMOMETER AND THE
177 C LOWEST ABOVE-GROUND FLOW SURFACE.
178 C
179  ibr=2
180 C
181 C GET VALUES TO USE FOR BOTTOM END OF INTERPOLATION.
182 C
183  bottom=z10
184  u0=u(ix,iy,1)
185  v0=v(ix,iy,1)
186  w0=w(ix,iy,1)
187  IF (dothet .OR. dopres .OR. dobri .OR. dobvpd) THEN
188  t0=tmpkel(ix,iy,1)
189  t1=tmpkel(ix,iy,levbot(ix,iy))
190  alt0=setalt(ix,iy,levbot(ix,iy))
191 C IPRGRF(IX,IY,IZ)=-NINT(100.0*ALT0)
192  iprgrf(ix,iy,iz)=nint(10.0*cvt2p(alt0,zmsl))
193  END IF
194 C
195 C NOW GET VALUES TO USE FOR TOP END OF INTERPOLATION.
196 C
197  top=rhs(ix,iy,levbot(ix,iy))
198  u1=u(ix,iy,levbot(ix,iy))
199  v1=v(ix,iy,levbot(ix,iy))
200  w1=w(ix,iy,levbot(ix,iy))
201 C
202  ELSE IF (zagl.GE. rhs(ix,iy,levbot(ix,iy))
203  $ .AND. zagl .LT. rhs(ix,iy,nlvl)) THEN
204  dotrpo=.true.
205 C
206 C GETTING VALUES WHEN INTERPOLATION SFC IS BETWEEN THE LOWEST AND
207 C THE HIGHEST FLOW SFCS
208 C
209  DO 130 ll=levbot(ix,iy),nlvl-1
210  didfnd=.false.
211  IF (zagl.GE.rhs(ix,iy,ll) .AND.
212  $ zagl.LT.rhs(ix,iy,ll+1)) THEN
213 C
214 C FOR INTERPOLATION SFC BETWEEN 2 FLOW SFCS
215 C
216  ibr=3
217  bottom=rhs(ix,iy,ll)
218  u0=u(ix,iy,ll)
219  v0=v(ix,iy,ll)
220  w0=w(ix,iy,ll)
221 C
222 C WIND VALUES TO USE FOR TOP END OF INTERPOLATION.
223 C
224  top=rhs(ix,iy,ll+1)
225  u1=u(ix,iy,ll+1)
226  v1=v(ix,iy,ll+1)
227  w1=w(ix,iy,ll+1)
228 C
229 C USE WEIGHTED ALTIMETER SETTING FROM FLOW SFCS.
230 C
231  IF (zagl.EQ.rhs(ix,iy,ll)) THEN
232  alt0=setalt(ix,iy,ll)
233  ELSE IF (rhs(ix,iy,ll+1) .EQ. zagl) THEN
234  alt0=setalt(ix,iy,ll+1)
235  ELSE
236  wt1=1.0/abs(zagl-rhs(ix,iy,ll))
237  wt2=1.0/abs(zagl-rhs(ix,iy,ll+1))
238  alt0=(wt1*setalt(ix,iy,ll)+
239  $ wt2*setalt(ix,iy,ll+1))/(wt1+wt2)
240  END IF
241 C IPRGRF(IX,IY,IZ)=NINT(100.0*ALT0)
242  iprgrf(ix,iy,iz)=nint(10.0*cvt2p(alt0,zmsl))
243  IF (dothet .OR. dopres .OR. dobri
244  $ .OR. dobvpd) THEN
245  t0=tmpkel(ix,iy,ll)
246  t1=tmpkel(ix,iy,ll+1)
247  END IF
248 C
249 C WE FOUND WHERE IT IS, WE CAN GO ON TO THE NEXT STEP
250 C
251  didfnd=.true.
252  END IF
253  IF (didfnd) GO TO 135
254 130 CONTINUE
255 135 CONTINUE
256 C
257  END IF
258 C
259 C IF WE HAVE NOT ALREADY DEFINED THE VALUES, THEN WE NEED
260 C TO INTERPLATE, USING VALUES DEFINED ABOVE.
261 C
262  IF (dotrpo) THEN
263 C
264 C IF ONE BVP VALUE, GOOD USE IT
265 C
266  CALL lgntrp(uuuf,bottom,top,zagl,u0,u1)
267  CALL lgntrp(vvvf,bottom,top,zagl,v0,v1)
268  CALL lgntrp(wwwf,bottom,top,zagl,w0,w1)
269  iugraf(ix,iy,iz)=nint(100.0*uuuf)
270  ivgraf(ix,iy,iz)=nint(100.0*vvvf)
271  iwgraf(ix,iy,iz)=nint(100.0*wwwf)
272 C
273 C VARIABLES ARE SCALED TO INTEGERS FOR WRITING FILES. UNITS ARE
274 C CM/SEC FOR WIND COMPONENTS, (DEG K)/10 FOR POTENTIAL TEMP., MB (HP)
275 C FOR PRESSURE, SECONDS FOR BRUNT VAISALA PERIOD, AND 100* BULK
276 C RICHARDSON NO.
277 C
278  aaset=xlintr(alt0,alt1,zagl,bottom,top)
279  IF (dothet .OR. dopres .OR. dobri.OR. dobvpd)
280  $ iptgrf(ix,iy,iz)=nint(10.0*
281  $ xlintr(t0,t1,zagl,bottom,top))
282 
283  END IF
284 C
285 300 CONTINUE
286 C
287 390 CONTINUE
288 400 CONTINUE
289 C
290 C CONVERT COMPONENTS IN HORIZONTAL OR TERRAIN FOLLOWING
291 C PLANES TO METEOROLOGICAL WINDS SPEED & DIRECTION
292 C
293  DO 500 lev=1,nflat
294  DO 490 ix=1,ncol
295  DO 480 iy=1,nrow
296 C
297 C CONVERT INTEGER CM/S BACK TO REAL M/S
298 C
299  ug=0.01*float(iugraf(ix,iy,lev))
300  vg=0.01*float(ivgraf(ix,iy,lev))
301 C
302  IF (ug.EQ.zero.AND.vg.EQ.zero) THEN
303  spdmet(ix,iy,lev)=zero
304  dirmet(ix,iy,lev)=zero
305  ELSE
306  spdmet(ix,iy,lev)=sp(ug,vg)
307  dirmet(ix,iy,lev)=dd(ug,vg)
308  END IF
309 480 CONTINUE
310 490 CONTINUE
311 500 CONTINUE
312 C
313  RETURN
314 C
lgntrp
subroutine lgntrp(Y, X0, X1, X, Y0, Y1)
Definition: utils.f:198
cvt2p
real function cvt2p(ALTIM, ELEV)
Definition: utils.f:129
fixwnd
subroutine fixwnd(IFXPT, LVL)
Definition: windest.f:3
xlintr
real function xlintr(X0, X1, Z, ZZ0, Z1)
Definition: utils.f:47
sp
real function sp(UU, VV)
Definition: utils.f:19
dd
real function dd(UU, VV)
Definition: utils.f:30

Here is the call graph for this function:

Here is the caller graph for this function:

Generated on Tue Sep 22 2015 23:54:04 for GNU WOCSS (gwocss) by   doxygen 1.8.10