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

Go to the source code of this file.

Functions/Subroutines

subroutine bal5 (NITER)
 

Function/Subroutine Documentation

subroutine bal5 (   NITER)

Definition at line 3 of file bal5.f.

Referenced by wocss().

3 C**********************************************************************
4 C
5 C THIS IS A MODIFIED VERSION OF ROY ENDLICH'S OCTOBER 1984 CODE. THE
6 C OPTIONS FOR VORTICITY CONSERVATION & NONZERO DIVERGENCE HAVE BEEN
7 C REMOVED. THE CODE WAS ALSO MODIFIED TO USE MORE
8 C IF-THEN-ELSE STRUCTURE BY LUDWIG DECEMBER 1987.
9 C
10 C THIS ROUTINE BALANCES DIVERGENCE TOWARD ZERO. DIV IS SCALED TO UNITS
11 C OF 10**-6/SECOND. THE METHOD USES DIRECT VECTOR ALTERATIONS.
12 C THIS FORM IS FOR A SQUARE GRID AND OMITS TRIGONOMETRIC FUNCTIONS.
13 C THE FLUX FORMULATION IS USED SO THAT WIND COMPONENTS ARE WEIGHTED BY
14 C THE THICKNESS OF THE LAYER WHEN THE FINITE DIFFERENCE SCHEME IS USED
15 C TO REDUCE THE DIVERGENCE. INDICES IN ARRAYS (I,J,K) ARE I=COLUMN,
16 C J=ROW, K=LEVEL; PT (1,1,1) IS AT SW CORNER AT GROUND. FOR
17 C COMPUTATION BOXES, INDICES REFER TO SW CORNER OF BOX.
18 C______________________________________________________________________
19 C
20 C MODIFIED MARCH 1996 SO THAT BELOW GROUND AND NEAR-STATION POINTS ARE
21 C TREATED SEPARATELY. ADJUSTMENTS AT GRID POINTS NEAR OBSERVATION
22 C SITES ARE ADJUSTED BY LESS (THE FACTOR ADJMAX) THAN OTHER POINTS.
23 C BELOW GROUND ARE LEFT UNADJUSTED AS BEFORE.
24 C
25 C F. L. LUDWIG 3/96
26 C______________________________________________________________________
27 C
28 C MODIFIFED SO THE ADJUSTMENTS OVER OVERRIDDEN AT THE END OF EACH
29 C ITERATION ALL POINTS ARE ADJUSTED, THEN THOSE NEAR OBSERVATIONS OR
30 C SURROUNDED BY 3 OR 4 SUBSURFACE POINT ARE SET BACK TO ORIGINAL, OR
31 C NEAR-ORIGINAL VALUES BEFORE STARTING NEXT ITERATION. ADDED SECTION
32 C TO GET VERTICAL MOTION ESTIMATES FROM SLOPE OF FLOW SURFACE AND
33 C HORIZONTAL COMPONENTS (FUNCTION DUBYOU)
34 C
35 C FLUDWIG 3/2000
36 C_____________________________________________________________________
37 C
38  include 'NGRIDS.PAR'
39 C
40  include 'FLOWER.CMM'
41  include 'LIMITS.CMM'
42  include 'STALOC.CMM'
43 C
44  LOGICAL doadj,ifxpt(nxgrd,nygrd),isdop
45  REAL di(nxgrd,nygrd),u1(nxgrd,nygrd),v1(nxgrd,nygrd)
46  REAL un(nxgrd,nygrd),vn(nxgrd,nygrd),thk(nxgrd,nygrd)
47  REAL ustart(nxgrd,nygrd),vstart(nxgrd,nygrd)
48  DATA doadj, enfrac /.false.,0.5/
49 C
50  gs=ds*1.0e-05
51 C
52 C USE GRID SPACING IN 100'S OF KM. DS IS IN M. FOR PROPER SCALING.
53 C
54  gsi=10.0/gs
55 C
56  DO 800 l=2,nlvl
57 C
58 C GET COMPONENTS FOR THIS LEVEL
59 C
60  DO 35 j=1,nrow
61  DO 33 i=1,ncol
62  un(i,j)=u(i,j,l)
63  vn(i,j)=v(i,j,l)
64 33 CONTINUE
65 35 CONTINUE
66 C
67 C IDENTIFY PTS NOT TO BE CHANGED
68 C
69  CALL fixwnd (ifxpt,l)
70 C
71  CALL setmat(0.0,di,ncol,nrow)
72 C
73 C COMPUTE LAYER THICKNESS AND MULTIPLY WIND COMPONENTS
74 C
75  DO 40 j=1,nrow
76  DO 40 i=1,ncol
77  la=l +1
78  IF (la.GT.nlvl) la=nlvl
79  hta=rhs(i,j,la)
80  lb=l-1
81  IF (lb.LT.1) lb=1
82  htb=rhs(i,j,lb)
83  IF (htb.LT.-1.0) htb=-1.0
84  thk(i,j)=0.5*(hta-htb)*0.01
85 C
86 C FOR NEG (OR VERY SMALL) RHS.
87 C
88  IF (thk(i,j).LE.0.01) thk(i,j)=0.01
89 C IF (THK(I,J).LE.0.1) THK(I,J)=0.1
90 C
91 C UNITS OF THICKNESS ARE HUNDREDS OF M FOR CONVENIENCE.
92 C SET INITIAL WINDS BEFORE ALTERATIONS.
93 C
94  u1(i,j)=un(i,j)
95  v1(i,j)=vn(i,j)
96 C
97 C WEIGHT WINDS WITH THICKNESS OF LAYER.
98 C
99  un(i,j)=u1(i,j)*thk(i,j)
100  vn(i,j)=v1(i,j)*thk(i,j)
101  40 CONTINUE
102 C
103 C COMPUTE DIVERGENCE (DI=DUE+DVN) FROM FLUX DIFFERENCES BETWEEN
104 C V COMPONENTS AT NORTH (VNO) & SOUTH (VSO) SIDES OF BOX & U
105 C COMPONENTS AT EAST (UE) & WEST (UW) SIDES OF BOX. DDIJ IS ZERO
106 C TO PRODUCE NONDIVERGENT WINDS; OTHER VALUES COULD BE USED TO
107 C ACCOMODATE AN AREA-WIDE GENERAL DIVERGENCE. RA=RELAXATION FACTOR.
108 C
109  ddij=0.0
110  ra=0.7
111  DO 240 lg=1,niter
112 C
113 C SAVE VALUES AT START OF ITERATION SO THE VALUES AT
114 C NEAR-OSERVATION AND OTHER FIXED POINTS CAN BE CHANGED BACK AT
115 C END OF ECH ITERATIVE STEP -- ADDED 3/2000 BY FLUDWIG
116 C
117  DO 70 j=1,nrow
118  DO 60 i=1,ncol
119  ustart(i,j)=un(i,j)
120  vstart(i,j)=vn(i,j)
121 60 CONTINUE
122 70 CONTINUE
123 C
124  DO 150 j=1,nrowm1
125  DO 140 i=1,ncolm1
126  ue=0.5*(un(i+1,j)+un(i+1,j+1))
127  uw=0.5*(un(i,j) +un(i,j+1))
128  vso=0.5*(vn(i+1,j)+vn(i,j))
129  vno=0.5*(vn(i,j+1)+vn(i+1,j+1))
130  due=gsi*(ue-uw)
131  dvn=gsi*(vno-vso)
132  di(i,j)=due+dvn
133  cuij=0.05*gs*(ddij-di(i,j))*ra
134  cvij=0.05*gs*(ddij-di(i,j))*ra
135 C
136 C LIMIT CHANGES TO LESS THAN 1 FOR NUMERICAL STABILITY.
137 C
138  IF (cuij .LT.-1.0) cuij=-1.0
139  IF (cuij .GT. 1.0) cuij=1.0
140  IF (cvij .LT.-1.0) cvij=-1.0
141  IF (cvij .GT. 1.0) cvij=1.0
142 C
143  un(i+1,j)=un(i+1,j)+cuij
144  un(i+1,j+1)=un(i+1,j+1) +cuij
145  un(i,j)=un(i,j) -cuij
146  un(i,j+1)=un(i,j+1) -cuij
147  vn(i+1,j)=vn(i+1,j)-cvij
148  vn(i,j)=vn(i,j)-cvij
149  vn(i,j+1)=vn(i,j+1)+cvij
150  vn(i+1,j+1)=vn(i+1,j+1)+cvij
151 140 CONTINUE
152 C
153 150 CONTINUE
154 C
155 C GO BACK AND SUBSTITUTE ORIGINAL VALUES, CORRECTED FOR MAXIMUM
156 C SPECIFIED ADJUSTMENT AT FIXED POINTS. INSERT ZEROS FOR SUBSURFACE
157 C POINTS. INCREASE ADJUSTMENTS FOR HIGHER LEVELS.
158 C
159 C ADDED BY FLUDWIG 3/2000
160 C
161  DO 170 j=1,nrow
162  DO 160 i=1,ncol
163 C
164  IF (ifxpt(i,j)) THEN
165 C
166 C FOR UPPER WIND SITE USE MINIMUM ADJUSTMENT AT ALL LEVELS
167 C
168  isdop=.false.
169  DO 153 jd=1,numdop
170  IF (i.EQ. nint(xdop(jd)) .AND.
171  $ j.EQ. nint(ydop(jd))) isdop=.true.
172 153 CONTINUE
173 C
174  IF (l.EQ.levbot(i,j) .OR. isdop) THEN
175  un(i,j)=ustart(i,j)+
176  $ adjmax*(un(i,j)-ustart(i,j))
177  vn(i,j)=vstart(i,j)+
178  $ adjmax*(vn(i,j)-vstart(i,j))
179 C
180  ELSE IF (l.GT.levbot(i,j)) THEN
181 C
182 C INCREASE ADJUSTMENT ALLOWED WITH HEIGHT AT OTHER FIXED POINTS
183 C
184  adjless=adjmax+(1.0-adjmax)*
185  $ (rhs(i,j,l)-rhs(i,j,levbot(i,j)))/
186  $ (rhs(i,j,nlvl)-rhs(i,j,levbot(i,j)))
187 C
188  un(i,j)=ustart(i,j)+
189  $ adjless*(un(i,j)-ustart(i,j))
190  vn(i,j)=vstart(i,j)+
191  $ adjless*(vn(i,j)-vstart(i,j))
192  END IF
193  END IF
194 160 CONTINUE
195 170 CONTINUE
196 C
197 C END OF ITERATION LOOP
198 C
199 240 CONTINUE
200 C
201  sum1=0.0
202  sum2=0.0
203  q1=0.0
204  DO 350 j=1,nrow
205  DO 340 i=1,ncol
206 C
207 C INCLUDE ONLY POINTS THAT ARE ABOVE THE SURFACE AND
208 C ARE NOT NEAREST TO AN OBSERVATION SITE.
209 C
210  IF (rhs(i,j,l).GT.0.0 ) THEN
211  un(i,j)=un(i,j)/thk(i,j)
212  vn(i,j)=vn(i,j)/thk(i,j)
213  u1(i,j)=u1(i,j)-un(i,j)
214  v1(i,j)=v1(i,j)-vn(i,j)
215  q1=q1+1.0
216  sum1=sum1+u1(i,j)
217  sum2=sum2+v1(i,j)
218  END IF
219 340 CONTINUE
220 350 CONTINUE
221 C
222  sum1=sum1/q1
223  sum2=sum2/q1
224 C
225 C NORMALIZE ORIGINAL AVERAGE VALUES FOR ABOVE GROUND POINTS
226 C THAT ARE NOT NEAREST OBSERVATION SITES.
227 C
228  DO 450 j=1,nrow
229  DO 445 i=1,ncol
230  IF (rhs(i,j,l) .GT.0.0) THEN
231  un(i,j)=un(i,j)+sum1
232  vn(i,j)=vn(i,j)+sum2
233  ELSE
234  un(i,j)=0.0
235  vn(i,j)=0.0
236  END IF
237 445 CONTINUE
238 450 CONTINUE
239 C
240 C CHANGE BACK TO 3D ARRAYS
241 C
242  DO 590 j=1,nrow
243  DO 580 i=1,ncol
244  u(i,j,l)=un(i,j)
245  v(i,j,l)=vn(i,j)
246  w(i,j,l)=dubyou(un(i,j),vn(i,j),i,j,l)
247 580 CONTINUE
248 590 CONTINUE
249 C
250 C END FLOW LEVEL LOOP
251 C
252 800 CONTINUE
253 C
254  RETURN
255 C
setmat
subroutine setmat(VALUE, ARRAY, NUM1, NUM2)
Definition: utils.f:258
fixwnd
subroutine fixwnd(IFXPT, LVL)
Definition: windest.f:3
dubyou
real function dubyou(UT, VT, I, J, K)
Definition: vertvel.f:3

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