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

Go to the source code of this file.

Functions/Subroutines

subroutine strat
 

Function/Subroutine Documentation

subroutine strat ( )

Definition at line 3 of file strat.f.

Referenced by geosig().

3 C**********************************************************************
4 C
5 C READS TEMPERATURE SOUNDING INFORMATION,CALCULATES LAPSE RATES &
6 C OTHER PARAMETERS FOR VARIOUS FLOW LEVELS.
7 C MODIFIED BY LUDWIG, JAN 2002.
8 C
9  include 'NGRIDS.PAR'
10 C
11  include 'ANCHOR.CMM'
12  include 'FLOWER.CMM'
13  include 'LIMITS.CMM'
14  include 'STALOC.CMM'
15  include 'TSONDS.CMM'
16 C
17  REAL rhs1(ntsite,nzgrd),sythyt(ntsite),pmids(nsndht,ntsite)
18  INTEGER lowzee(ntsite,nzgrd),nxtzee(ntsite,nzgrd)
19  INTEGER levhi(ntsite),levtop,ixrass(2),iyrass(2)
20 C
21 C INTEGER COORDINATES OF RASS SITES -- NOT USED FOR UPPER PRESSURE
22 C
23  DATA ixrass,iyrass /423,432,4497,4499/
24 C
25  IF(numtmp.LE.0) THEN
26  WRITE (*,*) 'NO TEMP SONDE DATA -- CANNOT RUN'
27  stop
28  END IF
29 C
30  levtop=0
31 C
32  DO 100 jtsond=1,numtmp
33 C
34 C FILL ARRAYS WITH MSG DATA IDENTIFIER
35 C
36  DO 8 i=1,nsndht
37  zsnd(i,jtsond)=-9999.
38  tsnd(i,jtsond)=-9999.
39  psnd(i,jtsond)=-9999.
40  altmtr(i,jtsond)=-9999.
41  zmids(i,jtsond)=-9999.
42  pmids(i,jtsond)=-9999.
43  dptdzs(i,jtsond)=-9999.
44  potemp(i,jtsond)=-9999.
45 8 CONTINUE
46  ii=jtsond
47 C
48 C READ STATION NO., NO. OF HEIGHTS (NHITES) IN SOUNDING & DATA INPUT
49 C TYPE (ITYP) DATA ARE READ AS FOLLOWS FOR ITYP= :
50 C 1--HEIGHT(M--MSL),TEMPERATURE(C), POTENTIAL TEMPERATURE
51 C LAPSE RATE(K/M) -- THIS OPTION NOT AVAILABLE HERE
52 C 2--HEIGHT(M--MSL), TEMP(C) & PRESSURE(MB)
53 C 3--HEIGHT(IN FT. MSL), TEMP(C), PRESSURE (MB)
54 C
55 C SKIP BLANK LINE
56 C
57  READ(12,*)
58  READ(12,*) xs,ys,sythyt(jtsond)
59  IF ((nint(xs) .EQ. ixrass(1) .AND. nint(ys).EQ.iyrass(1))
60  $ .OR. (nint(xs).EQ.ixrass(2)
61  $ .AND. nint(ys).EQ.iyrass(2)) ) THEN
62  goodt(jtsond)=.false.
63  ELSE
64  goodt(jtsond)=.true.
65  END IF
66  xtmp(jtsond)=1.0+(xs-xorig)/dscrs
67  ytmp(jtsond)=1.0+(ys-yorig)/dscrs
68  READ(12,*) nthts(jtsond),ityp
69 C
70 C SKIP COLUMN HEADERS
71 C
72  READ(12,*)
73 C
74  nhites=nthts(jtsond)
75  IF (nhites .GT.0) THEN
76 C
77  i=0
78 C
79  DO 10 j=1,nhites
80  READ(12,*,end=186) z,t,p
81  IF(nint(z) .LE.-9999 .OR. nint(t).LE.-9999
82  $ .OR. nint(p) .LE.-9999 .OR.
83  $ z-sythyt(jtsond).LT. -1.0) GO TO 10
84  i=i+1
85  IF (ityp.GE.2) THEN
86 C
87 C CONVERT FEET TO METERS & T TO POT TMP, IF REQUIRED (ITYP=3)
88 C CONVERT TO HEIGHT ABOVE SFC--
89 C
90  IF (ityp .EQ. 3 ) z=0.3048*z
91  IF (i .LE. nsndht) THEN
92  zsnd(i,jtsond)=z-sythyt(jtsond)
93  potemp(i,jtsond)=tpot(p,t)
94  tsnd(i,jtsond)=t+273.13
95  psnd(i,jtsond)=p
96  altmtr(i,jtsond)=altset(z,p)
97  IF (altmtr(i,jtsond).LT. 28.) WRITE(*,*)
98  $ i,' CHECK T-SONDE ',jtsond,z,p,altset(z,p)
99  ELSE
100  zsnd(nsndht,jtsond)=z-sythyt(jtsond)
101  potemp(nsndht,jtsond)=tpot(p,t)
102  tsnd(nsndht,jtsond)=t+273.13
103  psnd(nsndht,jtsond)=p
104  altmtr(i,jtsond)=altset(z,p)
105  IF (altmtr(i,jtsond).LT. 28.) WRITE(*,*)
106  $ i,' CHECK T-SONDE ',jtsond,z,p,altset(z,p)
107  END IF
108  ELSE IF (ityp.EQ.1) THEN
109 C
110 C REMOVED LAPSE RATE INPUT OPTION FOR SLC VTMX APPLICATION
111 C FLUDWIG 2/02
112 C*********************************************************
113  WRITE (*,*) 'TEMP SOUNDINGS MUST PROVIDE HEIGHT ',
114  $ '(M, TYP2 OR FT, TYP3), TEMP (DEG C) & ',
115  $ 'PRESS (MB/HP)'
116  stop
117  END IF
118 C
119 10 CONTINUE
120 C
121 C EMPTY SOUNDING
122 C
123  ELSE
124  WRITE (*,*) 'TEMPERATURE SONDE ', jtsond, ' HAS NO OBS.'
125  stop
126  END IF
127 186 nhites=i
128  nthts(jtsond)=i
129 C
130 C SOUNDING READ; FIND HEIGHTS OF FLOW SURFACES
131 C FIRST, FIND NEAREST LOW REFERENCE POINT
132 C
133  idis2=(nint(xtmp(jtsond))-lowix(1))**2+
134  $ (nint(ytmp(jtsond))-lowiy(1))**2
135  jnear(jtsond)=1
136  DO 17 jd=2,5
137  id2here=(nint(xtmp(jtsond))-lowix(jd))**2+
138  $ (nint(ytmp(jtsond))-lowiy(jd))**2
139  IF (id2here .LT. idis2) THEN
140  idis2=id2here
141  jnear(jtsond)=jd
142  END IF
143 17 CONTINUE
144 C
145  lllx=lowix(jnear(jtsond))
146  llly=lowiy(jnear(jtsond))
147 C
148 C CHECK TO SEE IF TOP OBSERVATION IS ABOVE THIS SURFACE. RECORD
149 C HIGHEST LEVELS REACHED FOR INDIVIDUAL SONDES AND ALL SONDES.
150 C
151  stoplev=zsnd(nthts(jtsond),jtsond)
152  DO 25 kl=1,nlvl
153  zlevkl=rhs(lllx,llly,kl)
154  rhs1(jtsond,kl)=zlevkl
155 C
156 C FINDING HIGHEST FLOW SURFACE REACHED BY THIS SONDE
157 C
158  IF (stoplev .GE. zlevkl) THEN
159  levhi(jtsond)=kl
160  IF (kl .GT. levtop) levtop=kl
161  END IF
162 25 CONTINUE
163 C
164 C GETTING INDICES OF OBS THAT ARE NEAREST MIDPOINTS ABOVE & BELOW
165 C THE FLOW SURFACES.
166 C
167  lowzee(jtsond,1)=1
168  DO 45 kl=1,levhi(jtsond)
169  zlevkl=rhs1(jtsond,kl)
170  IF (kl.LT.nlvl) THEN
171  zbar=0.5*(zlevkl+rhs1(jtsond,kl+1))
172  ELSE IF (kl.EQ.nlvl) THEN
173  zbar=zlevkl+0.5*(zlevkl-rhs1(jtsond,kl-1))
174  END IF
175  close1=2000000.
176  DO 40 iz=1,nthts(jtsond)-1
177  IF (abs(zbar-zsnd(iz,jtsond)).LT.close1) THEN
178  close1=abs(zsnd(iz,jtsond)-zbar)
179  nxtzee(jtsond,kl)=iz
180  IF (nxtzee(jtsond,kl).EQ.
181  $ lowzee(jtsond,kl)) nxtzee(jtsond,kl)=iz+1
182 C
183 C LOWER POINT FOR NEXT SURFACE WILL BE SAME AS UPPER POINT FOR THIS
184 C ONE
185 C
186  IF (kl .LT. nlvl) THEN
187  lowzee(jtsond,kl+1)=iz
188  END IF
189  END IF
190 40 CONTINUE
191 45 CONTINUE
192 C
193 C CHECKING TO SEE THAT FLOW SFCS ARE BETWEEN BOUNDS WHEN POSSIBLE.
194 C
195  DO 70 kl=2,levhi(jtsond)
196  IF (zsnd(lowzee(jtsond,kl),jtsond) .GE.
197  $ rhs1(jtsond,kl) .AND. lowzee(jtsond,kl) .GT. 1) THEN
198 C
199 C PICK 1ST PLACE YOU COME TO BELOW THE PROPER MIDPOINT
200 C
201  DO 60 iz=lowzee(jtsond,kl),1,-1
202  zbar=0.5*(rhs1(jtsond,kl-1)+rhs1(jtsond,kl))
203  IF (zsnd(iz,jtsond) .LT. zbar) THEN
204  lowzee(jtsond,kl)=iz
205  GO TO 70
206  END IF
207 60 CONTINUE
208  END IF
209 70 CONTINUE
210 C
211 C FINALLY, CHECK TO SEE IF TOP PT IS HIGH ENOUGH (ABOVE MIDPT TO
212 C NEXT SURFACE) TO EXTRAPOLATE TO ONE MORE LEVEL
213 C
214  IF (levhi(jtsond).LT.nlvl) THEN
215  zbar=0.5*(rhs1(jtsond,levhi(jtsond))+
216  $ rhs1(jtsond,levhi(jtsond)+1))
217  IF (zsnd(nthts(jtsond),jtsond) .GE. zbar) THEN
218  levhi(jtsond)=1+levhi(jtsond)
219  IF (levhi(jtsond).GT.levtop) levtop=levhi(jtsond)
220  lowzee(jtsond,levhi(jtsond))=
221  $ (lowzee(jtsond,levhi(jtsond)-1) +
222  $ nxtzee(jtsond,levhi(jtsond)-1))/2
223  nxtzee(jtsond,levhi(jtsond))=nthts(jtsond)
224  END IF
225  END IF
226 100 CONTINUE
227 C
228 C ALL GOOD DATA READ; CONVERT WHERE NECESSARY & FIND HEIGHTS NEAREST
229 C THE POINTS TO USE FOR ESTIMATING GRADIENTSAT FLOW SURFACE HEIGHTS
230 C HAVE BEEN FOUND. NOW, GET POTENTIAL TEMP LAPSE RATES TO TOP OF
231 C EACH SOUNDING AND ESTIMATE LAPSE RATES ON THE FLOW SFCS.
232 C
233  DO 170 jtsond=1,numtmp
234  DO 160 kl=levhi(jtsond),1,-1
235 C
236 C USE OBS UP TO THE HIGHEST OB. INTERPOLATE STD ATMOS HTS
237 C CORRESPONDING TO OBSERVED PRESSURES -- THEN GET PRESSURE
238 C CORRSPONDING TO THE INTERPOLATED HT. OTHER VARIABLES ARE OBTAINED
239 C BY SIMPLE INTERP, OR FINITE DIFF ESIMATE.
240 C
241  IF (kl.GT.1) THEN
242  zz0=zsnd(lowzee(jtsond,kl),jtsond)+sythyt(jtsond)
243  z1=zsnd(nxtzee(jtsond,kl),jtsond)+sythyt(jtsond)
244  zlvl=rhs1(jtsond,kl)+sythyt(jtsond)
245  aaa1=altmtr(nxtzee(jtsond,kl),jtsond)
246  aaa0=altmtr(lowzee(jtsond,kl),jtsond)
247  tt0=tsnd(lowzee(jtsond,kl),jtsond)
248  tt1=tsnd(nxtzee(jtsond,kl),jtsond)
249  pt0=potemp(lowzee(jtsond,kl),jtsond)
250  pt1=potemp(nxtzee(jtsond,kl),jtsond)
251  ptlaps(jtsond,kl)=(pt1-pt0)/(z1-zz0)
252  altsig(jtsond,kl)=xlintr(aaa0,aaa1,zlvl,zz0,z1)
253  psigl(jtsond,kl)=cvt2p(altsig(jtsond,kl),zlvl)
254  ptsigl(jtsond,kl)=xlintr(pt0,pt1,zlvl,zz0,z1)
255  tsigl(jtsond,kl)=
256  $ pt2t(psigl(jtsond,kl),ptsigl(jtsond,kl))
257 C
258  ELSE IF (kl.EQ.1) THEN
259 C
260 C EXTRAPOLATE TO 10 M LEVEL
261 C
262  z1=zsnd(2,jtsond)+sythyt(jtsond)
263  zz0=zsnd(1,jtsond)+sythyt(jtsond)
264  zlvl=sythyt(jtsond)
265  ptlaps(jtsond,kl)=
266  $ (potemp(2,jtsond)-potemp(1,jtsond))/(z1-zz0)
267  IF (abs(zlvl-zz0).LT.0.1) THEN
268  altsig(jtsond,kl)=altmtr(1,jtsond)
269  ptsigl(jtsond,kl)=potemp(1,jtsond)
270  tsigl(jtsond,kl)=tsnd(1,jtsond)
271  psigl(jtsond,kl)=cvt2p(altsig(jtsond,kl),zlvl)
272  ELSE IF (abs(zlvl-z1) .LT. 0.1) THEN
273  altsig(jtsond,kl)=altmtr(2,jtsond)
274  ptsigl(jtsond,kl)=potemp(2,jtsond)
275  psigl(jtsond,kl)=cvt2p(altsig(jtsond,kl),zlvl)
276  tsigl(jtsond,kl)=
277  $ pt2t(psigl(jtsond,kl),ptsigl(jtsond,kl))
278  ELSE
279  pt1=potemp(2,jtsond)
280  pt0=potemp(1,jtsond)
281  tt1=tsnd(2,jtsond)
282  tt0=tsnd(1,jtsond)
283  aaa1=altmtr(2,jtsond)
284  aaa0=altmtr(1,jtsond)
285  altsig(jtsond,kl)=xlintr(aaa0,aaa1,zlvl,zz0,z1)
286  ptsigl(jtsond,kl)=xlintr(pt0,pt1,zlvl,zz0,z1)
287  psigl(jtsond,kl)=cvt2p(altsig(jtsond,kl),zlvl)
288  tsigl(jtsond,kl)=
289  $ pt2t(psigl(jtsond,kl),ptsigl(jtsond,kl))
290  END IF
291  END IF
292 C
293 160 CONTINUE
294 170 CONTINUE
295 C
296 C NOW GET THE VALUES ABOVE THOSE OBSERVED AT EACH SITE, BUT BELOW
297 C THE HIGHEST LEVEL OBSERVED ANYWHERE.
298 C
299  DO 250 jtsond=1,numtmp
300  IF (levhi(jtsond)+1.LE.levtop) THEN
301  DO 240 kl=levhi(jtsond)+1,levtop
302 C
303 C INVERSE DISTANCE HORIZ INTERP FROM THE HIGHEST OB HERE TO
304 C HIGHEST OB ANYWHERE
305 C
306  xx=xtmp(jtsond)
307  yy=ytmp(jtsond)
308  sumpt=0.0
309  sumdpt=0.0
310  sumwt=0.0
311  sumalt=0.0
312  DO 210 jd=1,numtmp
313 C
314 C SKIP WHEN LOOKING AT SAME SOUNDING BEING INTERPOLATED FOR.
315 C
316  IF (jd.NE.jtsond .AND. levhi(jd).GE.kl) THEN
317  weight=wndwt(xx,yy,xtmp(jd),ytmp(jd))
318  sumwt=sumwt+weight
319  sumdpt=sumdpt+weight*ptlaps(jd,kl)
320  sumpt=sumt+weight*ptsigl(jd,kl)
321 C
322 C USE INTERPOLATED ALTIMETER SETTING TO GET ESTIMATE TO BE USED
323 C FOR PRESSURE.
324 C
325  sumalt=sumalt+weight*altsig(jd,kl)
326  END IF
327 210 CONTINUE
328 C
329  IF (sumwt.GT.0.0) THEN
330  ptlaps(jtsond,kl)=sumdpt/sumwt
331  altsig(jtsond,kl)=sumalt/sumwt
332  ptsigl(jtsond,kl)=sumpt/sumwt
333  psigl(jtsond,kl)=cvt2p(altsig(jtsond,kl),
334  $ rhs1(jtsond,kl)+sythyt(jtsond))
335  tsigl(jtsond,kl)=pt2t(psigl(jtsond,kl),
336  $ ptsigl(jtsond,kl))
337 C
338 C CHECK FOR SUPERADIABATIC LAPSE RATE
339 C
340  IF (ptsigl(jtsond,kl).LT.ptsigl(jtsond,kl-1)) THEN
341  ptsigl(jtsond,kl)=ptsigl(jtsond,kl-1)
342  tsigl(jtsond,kl)=
343  $ pt2t(psigl(jtsond,kl),ptsigl(jtsond,kl))
344  END IF
345 C
346  ELSE
347  WRITE (*,*) 'CHECK STRAT NEAR 210'
348  stop
349  END IF
350 240 CONTINUE
351  END IF
352 C
353 250 CONTINUE
354 C
355 C FOR LEVELS ABOVE THE HIGHEST OB, USE LAST VALUE OF LAPSE
356 C BASED ON OBS AND EXTRPOLATE TEMPERATURE.
357 C
358  IF (levtop.LT.nlvl) THEN
359  DO 275 jtsond=1,numtmp
360 C
361 C IF NO SONDE REACHES LOWEST LEVEL, USE A MODERATELY STABLE LAPSE
362 C RATE OF 0.005 DEGK/M AND SFC ALTIMETER SETTING
363 C
364  IF (levtop.LT.2) THEN
365  ptlaps(jtsond,1)=0.005
366  levtop=1
367  ELSE
368  padjust=altsig(jtsond,levtop)
369  END IF
370 C
371  DO 265 kl=levtop+1,nlvl
372  altsig(jtsond,kl)=altsig(jtsond,levtop)
373  psigl(jtsond,kl)=cvt2p(altsig(jtsond,levtop),
374  $ rhs1(jtsond,kl)+sythyt(jtsond))
375  ptlaps(jtsond,kl)=ptlaps(jtsond,levtop)
376 C
377 C CHECK FOR SUPERADIABATIC LAPSE RATE
378 C
379  IF (ptlaps(jtsond,kl).LT.0.0) ptlaps(jtsond,kl)=0.0
380 C
381 C ESTIMATE PRESSURE AT LAST LEVEL FROM STAND ATMOS HEIGHT
382 C USE IT TO ESTIMATE PTENTIAL TEMP AT THIS LEVEL BY EXTRAPOL
383 C FROM PREVIOUS LEVEL.
384 C
385  ptsigl(jtsond,kl)=ptsigl(jtsond,levtop) +
386  $ ((rhs1(jtsond,kl)-rhs1(jtsond,levtop))*
387  $ ptlaps(jtsond,levtop))
388 C
389 C NOW CONVERT POTENT TEMP BACK TO TEMP
390 C
391  tsigl(jtsond,kl)=
392  $ pt2t(psigl(jtsond,kl),ptsigl(jtsond,kl))
393 265 CONTINUE
394 275 CONTINUE
395  END IF
396 C
397  RETURN
398 C
cvt2p
real function cvt2p(ALTIM, ELEV)
Definition: utils.f:129
wndwt
real function wndwt(X, Y, XOBS, YOBS)
Definition: utils.f:216
xlintr
real function xlintr(X0, X1, Z, ZZ0, Z1)
Definition: utils.f:47
pt2t
real function pt2t(PP, TT)
Definition: utils.f:80
altset
real function altset(Z, P)
Definition: utils.f:94
tpot
real function tpot(PP, TT)
Definition: utils.f:66

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