source: trunk/00FlowSolve_PL/PROJECTS/NEW_SRC/apply_BCs.f90

Last change on this file was 11, checked in by xlvlod, 18 years ago

Import initial
File size: 12.0 KB
Line 
1subroutine apply_BCs
2#define DEBUG_LEVEL 1
3
4! apply BCs to intermediate velocities and pd
5 use grid_info
6 use dependent_variables
7 use intermediate_variables
8 use mpi_parameters
9 use boundary_information
10 use counters_flags_etc
11 use user_parameters   ! for surface_V_of_x
12 use dimensional_scales
13 implicit none
14 include '../input/problem_size.h'
15 include 'mpif.h'
16
17 real, dimension(:,:), allocatable            :: f0,f1,f2,f3,f4
18 real                                         :: h, xval,F_of_t,thresh
19
20#if DEBUG_LEVEL >= 1
21 if(myid==0) write(0,*) 'hello world from apply_BCs'
22#endif
23
24! FACE 1: xi=0
25#if DEBUG_LEVEL >= 2
26 write(0,*) 'apply_BCs, FACE 1',myid
27#endif
28
29 h=Grid(0)%dxi
30 allocate (f0(ny,locnz),f1(ny,locnz),f2(ny,locnz),f3(ny,locnz),f4(ny,locnz))
31 i=1
32 if( BC_vel1(1,1) == 'specified_value' ) then
33  U(i,1:ny,1:locnz)=BC_values1(1:ny,1:locnz,1)
34  v(i,1:ny,1:locnz)=BC_values1(1:ny,1:locnz,2)
35  W(i,1:ny,1:locnz)=BC_values1(1:ny,1:locnz,3)
36 elseif( BC_vel1(1,1) == 'specified_stress' ) then
37  if( boundary_type1(1,1) == 'solid' ) then
38   U(i,1:ny,1:locnz)=BC_values1(1:ny,1:locnz,1)
39  elseif( boundary_type1(1,1) == 'permeable' ) then
40   f1=U(2,:,:)
41   f2=U(3,:,:)
42   f3=U(4,:,:)
43   f4=U(5,:,:)
44   f0=(-24.*h*BC_values1(:,:,1) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
45   U(i,1:ny,1:locnz)=f0
46  endif
47
48  f1=v(2,:,:)
49  f2=v(3,:,:)
50  f3=v(4,:,:)
51  f4=v(5,:,:)
52  f0=(-24.*h*BC_values1(:,:,2) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
53  v(i,1:ny,1:locnz)=f0
54  f1=w(2,:,:)
55  f2=w(3,:,:)
56  f3=w(4,:,:)
57  f4=w(5,:,:)
58  f0=(-24.*h*BC_values1(:,:,3) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
59  w(i,1:ny,1:locnz)=f0
60 endif
61
62! bcs for scalar 1 at face 1
63 if( BC_scalar1(1,1,1) == 'specified_value' ) then
64  s1(i,1:ny,1:locnz)=BC_values1(1:ny,1:locnz,4)
65 elseif( BC_scalar1(1,1,1) == 'specified_gradient' ) then
66  f1=s1(2,:,:)
67  f2=s1(3,:,:)
68  f3=s1(4,:,:)
69  f4=s1(5,:,:)
70  f0=(-24.*h*BC_values1(:,:,4) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
71  s1(i,1:ny,1:locnz)=f0
72 endif
73
74! bcs for scalar 2 at face 1
75if(num_scalars>1) then
76 if( BC_scalar1(1,1,2) == 'specified_value' ) then
77  s2(i,1:ny,1:locnz)=BC_values1(1:ny,1:locnz,5)
78 elseif( BC_scalar1(1,1,2) == 'specified_gradient' ) then
79  f1=s2(2,:,:)
80  f2=s2(3,:,:)
81  f3=s2(4,:,:)
82  f4=s2(5,:,:)
83  f0=(-24.*h*BC_values1(:,:,5) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
84  s2(i,1:ny,1:locnz)=f0
85 endif
86endif
87
88
89! FACE 2: xi=1
90#if DEBUG_LEVEL >= 2
91 write(0,*) 'apply_BCs, FACE 2',myid
92#endif
93
94 i=nx
95 if( BC_vel2(1,1) == 'specified_value' ) then
96  U(i,1:ny,1:locnz)=BC_values2(1:ny,1:locnz,1)
97  v(i,1:ny,1:locnz)=BC_values2(1:ny,1:locnz,2)
98  W(i,1:ny,1:locnz)=BC_values2(1:ny,1:locnz,3)
99 elseif( BC_vel2(1,1) == 'specified_stress' ) then
100  if( boundary_type2(1,1) == 'solid' ) then
101   U(i,1:ny,1:locnz)=BC_values2(1:ny,1:locnz,1)
102  elseif( boundary_type2(1,1) == 'permeable' ) then
103   f0=U(nx-4,:,:)
104   f1=U(nx-3,:,:)
105   f2=U(nx-2,:,:)
106   f3=U(nx-1,:,:)
107   f4=(24.*h*BC_values2(:,:,1) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
108   U(i,1:ny,1:locnz)=f0
109  endif
110
111   f0=v(nx-4,:,:)
112   f1=v(nx-3,:,:)
113   f2=v(nx-2,:,:)
114   f3=v(nx-1,:,:)
115   f4=(24.*h*BC_values2(:,:,2) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
116   v(i,1:ny,1:locnz)=f4
117   f0=w(nx-4,:,:)
118   f1=w(nx-3,:,:)
119   f2=w(nx-2,:,:)
120   f3=w(nx-1,:,:)
121   f4=(24.*h*BC_values2(:,:,3) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
122   w(i,1:ny,1:locnz)=f4
123 endif
124
125 if( BC_scalar2(1,1,1) == 'specified_value' ) then
126  s1(i,1:ny,1:locnz)=BC_values2(1:ny,1:locnz,4)
127 elseif( BC_scalar2(1,1,1) == 'specified_gradient' ) then
128  f0=s1(nx-4,:,:)
129  f1=s1(nx-3,:,:)
130  f2=s1(nx-2,:,:)
131  f3=s1(nx-1,:,:)
132  f4=(24.*h*BC_values2(:,:,4) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
133  s1(i,1:ny,1:locnz)=f4
134 endif
135
136if(num_scalars>1) then
137 if( BC_scalar2(1,1,2) == 'specified_value' ) then
138  s2(i,1:ny,1:locnz)=BC_values2(1:ny,1:locnz,5)
139 elseif( BC_scalar2(1,1,2) == 'specified_gradient' ) then
140  f0=s2(nx-4,:,:)
141  f1=s2(nx-3,:,:)
142  f2=s2(nx-2,:,:)
143  f3=s2(nx-1,:,:)
144  f4=(24.*h*BC_values2(:,:,5) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
145  s2(i,1:ny,1:locnz)=f4
146 endif
147endif
148
149 deallocate (f0,f1,f2,f3,f4)
150
151
152! FACE 3: eta=0
153#if DEBUG_LEVEL >= 2
154 write(0,*) 'apply_BCs, FACE 3',myid
155#endif
156
157 h=Grid(0)%deta
158 allocate (f0(nx,locnz),f1(nx,locnz),f2(nx,locnz),f3(nx,locnz),f4(nx,locnz))
159 j=1
160 if( BC_vel3(1,1) == 'specified_value' ) then
161  U(1:nx,j,1:locnz)=BC_values3(1:nx,1:locnz,1)
162  v(1:nx,j,1:locnz)=BC_values3(1:nx,1:locnz,2)
163  W(1:nx,j,1:locnz)=BC_values3(1:nx,1:locnz,3)
164 elseif( BC_vel3(1,1) == 'specified_stress' ) then
165  if( boundary_type3(1,1) == 'solid' ) then
166   v(1:nx,j,1:locnz)=BC_values3(1:nx,1:locnz,2)
167  elseif( boundary_type3(1,1) == 'permeable' ) then
168   f1=v(1:nx,2,1:locnz)
169   f2=v(1:nx,3,1:locnz)
170   f3=v(1:nx,4,1:locnz)
171   f4=v(1:nx,5,1:locnz)
172   f0=(-24.*h*BC_values3(:,:,2) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
173   v(1:nx,j,1:locnz)=f0
174  endif
175
176   f1=U(1:nx,2,1:locnz)
177   f2=U(1:nx,3,1:locnz)
178   f3=U(1:nx,4,1:locnz)
179   f4=U(1:nx,5,1:locnz)
180   f0=(-24.*h*BC_values3(:,:,1) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
181   U(1:nx,j,1:locnz)=f0
182   f1=W(1:nx,2,1:locnz)
183   f2=W(1:nx,3,1:locnz)
184   f3=W(1:nx,4,1:locnz)
185   f4=W(1:nx,5,1:locnz)
186   f0=(-24.*h*BC_values3(:,:,3) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
187   W(1:nx,j,1:locnz)=f0
188 endif
189
190 if( BC_scalar3(1,1,1) == 'specified_value' ) then
191  s1(1:nx,j,1:locnz)=BC_values3(1:nx,1:locnz,4)
192 elseif( BC_scalar3(1,1,1) == 'specified_gradient' ) then
193  f1=s1(1:nx,2,1:locnz)
194  f2=s1(1:nx,3,1:locnz)
195  f3=s1(1:nx,4,1:locnz)
196  f4=s1(1:nx,5,1:locnz)
197  f0=(-24.*h*BC_values3(:,:,4) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
198  s1(1:nx,j,1:locnz)=f0
199 endif
200
201if(num_scalars>1) then
202 if( BC_scalar3(1,1,2) == 'specified_value' ) then
203  s2(1:nx,j,1:locnz)=BC_values3(1:nx,1:locnz,5)
204 elseif( BC_scalar3(1,1,2) == 'specified_gradient' ) then
205  f1=s2(1:nx,2,1:locnz)
206  f2=s2(1:nx,3,1:locnz)
207  f3=s2(1:nx,4,1:locnz)
208  f4=s2(1:nx,5,1:locnz)
209  f0=(-24.*h*BC_values3(:,:,5) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
210  s2(1:nx,j,1:locnz)=f0
211 endif
212endif
213
214! FACE 4: eta=1
215#if DEBUG_LEVEL >= 2
216 write(0,*) 'apply_BCs, FACE 4',myid
217#endif
218
219 j=ny
220 if( BC_vel4(1,1) == 'specified_value' ) then
221  U(1:nx,j,1:locnz)=BC_values4(1:nx,1:locnz,1)
222  v(1:nx,j,1:locnz)=BC_values4(1:nx,1:locnz,2)
223  W(1:nx,j,1:locnz)=BC_values4(1:nx,1:locnz,3)
224 elseif( BC_vel4(1,1) == 'specified_stress' ) then
225  if( boundary_type4(1,1) == 'solid' ) then
226   v(1:nx,j,1:locnz)=BC_values4(1:nx,1:locnz,2)
227  elseif( boundary_type4(1,1) == 'permeable' ) then
228   f0=v(1:nx,ny-4,1:locnz)
229   f1=v(1:nx,ny-3,1:locnz)
230   f2=v(1:nx,ny-2,1:locnz)
231   f3=v(1:nx,ny-1,1:locnz)
232   f4=(24.*h*BC_values4(:,:,2) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
233   v(1:nx,j,1:locnz)=f4
234  endif
235
236   f0=U(1:nx,ny-4,1:locnz)
237   f1=U(1:nx,ny-3,1:locnz)
238   f2=U(1:nx,ny-2,1:locnz)
239   f3=U(1:nx,ny-1,1:locnz)
240   f4=(24.*h*BC_values4(:,:,1) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
241   U(1:nx,j,1:locnz)=f4
242   f0=W(1:nx,ny-4,1:locnz)
243   f1=W(1:nx,ny-3,1:locnz)
244   f2=W(1:nx,ny-2,1:locnz)
245   f3=W(1:nx,ny-1,1:locnz)
246   f4=(24.*h*BC_values4(:,:,3) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
247   W(1:nx,j,1:locnz)=f4
248 endif
249
250 if( BC_scalar4(1,1,1) == 'specified_value' ) then
251  s1(1:nx,j,1:locnz)=BC_values4(1:nx,1:locnz,4)
252 elseif( BC_scalar4(1,1,1) == 'specified_gradient' ) then
253  f0=s1(1:nx,ny-4,1:locnz)
254  f1=s1(1:nx,ny-3,1:locnz)
255  f2=s1(1:nx,ny-2,1:locnz)
256  f3=s1(1:nx,ny-1,1:locnz)
257  f4=(24.*h*BC_values4(:,:,4) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
258  s1(1:nx,j,1:locnz)=f4
259 endif
260
261if(num_scalars>1) then
262 if( BC_scalar4(1,1,2) == 'specified_value' ) then
263  s2(1:nx,j,1:locnz)=BC_values4(1:nx,1:locnz,5)
264 elseif( BC_scalar4(1,1,2) == 'specified_gradient' ) then
265  f0=s2(1:nx,ny-4,1:locnz)
266  f1=s2(1:nx,ny-3,1:locnz)
267  f2=s2(1:nx,ny-2,1:locnz)
268  f3=s2(1:nx,ny-1,1:locnz)
269  f4=(24.*h*BC_values4(:,:,5) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
270  s2(1:nx,j,1:locnz)=f4
271 endif
272endif
273
274 deallocate (f0,f1,f2,f3,f4)
275
276
277! FACE 5: zeta=0
278#if DEBUG_LEVEL >= 2
279 write(0,*) 'apply_BCs, FACE 5',myid
280#endif
281
282 if( myid == 0 ) then
283  h=Grid(0)%dzeta
284  allocate (f0(nx,ny),f1(nx,ny),f2(nx,ny),f3(nx,ny),f4(nx,ny))
285  k=1
286  if( BC_vel5(1,1) == 'specified_value' ) then
287   U(1:nx,1:ny,k)=BC_values5(1:nx,1:ny,1)
288   v(1:nx,1:ny,k)=BC_values5(1:nx,1:ny,2)
289   W(1:nx,1:ny,k)=BC_values5(1:nx,1:ny,3)
290  elseif( BC_vel5(1,1) == 'specified_stress' ) then
291   if( boundary_type5(1,1) == 'solid' ) then
292    W(1:nx,1:ny,k)=BC_values5(1:nx,1:ny,3)
293   elseif( boundary_type5(1,1) == 'permeable' ) then
294!!!    f1=W(:,:,2)
295!!!    f2=W(:,:,3)
296!!!    f3=W(:,:,4)
297!!!    f4=W(:,:,5)
298!!!    f0=(-24.*h*BC_values5(:,:,3) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
299!!!    W(:,:,k)=f0
300      W(:,:,k)=-9999.
301   endif
302
303    f1=U(:,:,2)
304    f2=U(:,:,3)
305    f0 = (-2.*h*BC_values5(:,:,1) + 4.*f1 - f2) / 3.
306!!!    f3=U(:,:,4)
307!!!    f4=U(:,:,5)
308!!!    f0=(-24.*h*BC_values5(:,:,1) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
309    U(:,:,k)=f0
310
311    f1=v(:,:,2)
312    f2=v(:,:,3)
313    f0 = (-2.*h*BC_values5(:,:,2) + 4.*f1 - f2) / 3.
314!!!    f3=v(:,:,4)
315!!!    f4=v(:,:,5)
316!!!    f0=(-24.*h*BC_values5(:,:,2) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
317    v(:,:,k)=f0
318  endif
319
320  if( BC_scalar5(1,1,1) == 'specified_value' ) then
321   s1(1:nx,1:ny,k)=BC_values5(1:nx,1:ny,4)
322  elseif( BC_scalar5(1,1,1) == 'specified_gradient' ) then
323   f1=s1(:,:,2)
324   f2=s1(:,:,3)
325   f0 = (-2.*h*BC_values5(:,:,4) + 4.*f1 - f2) / 3.
326!!!   f3=s1(:,:,4)
327!!!   f4=s1(:,:,5)
328!!!   f0=(-24.*h*BC_values5(:,:,4) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
329   s1(:,:,k)=f0
330  endif
331
332 if(num_scalars>1) then
333  if( BC_scalar5(1,1,2) == 'specified_value' ) then
334   s2(1:nx,1:ny,k)=BC_values5(1:nx,1:ny,5)
335  elseif( BC_scalar5(1,1,2) == 'specified_gradient' ) then
336   f1=s2(:,:,2)
337   f2=s2(:,:,3)
338   f0 = (-2.*h*BC_values5(:,:,5) + 4.*f1 - f2) / 3.
339!!!   f3=s2(:,:,4)
340!!!   f4=s2(:,:,5)
341!!!   f0=(-24.*h*BC_values5(:,:,5) + 96.*f1 -72.*f2 +32.*f3 -6.*f4)/50.
342   s2(:,:,k)=f0
343  endif
344 endif
345
346  deallocate (f0,f1,f2,f3,f4)
347 endif ! myid==0 block
348
349! FACE 6: zeta=1
350
351#if DEBUG_LEVEL >= 2
352 write(0,*) 'apply_BCs, FACE 6',myid
353#endif
354
355 if( myid == numprocs-1 ) then
356
357  h=Grid(0)%dzeta
358  allocate (f0(nx,ny),f1(nx,ny),f2(nx,ny),f3(nx,ny),f4(nx,ny))
359  k=locnz
360
361!**** ENFORCE STRESS FREE CONDITION OVER ENTIRE UPPER SURFACE ****************
362  if( BC_vel6(1,1) == 'specified_value' ) then
363   U(1:nx,1:ny,k)=BC_values6(1:nx,1:ny,1)
364   v(1:nx,1:ny,k)=BC_values6(1:nx,1:ny,2)
365   W(1:nx,1:ny,k)=BC_values6(1:nx,1:ny,3)
366  elseif( BC_vel6(1,1) == 'specified_stress' ) then
367   if( boundary_type6(1,1) == 'solid' ) then
368    W(1:nx,1:ny,k)=BC_values6(1:nx,1:ny,3)
369   elseif( boundary_type6(1,1) == 'permeable' ) then
370!!!    f0=W(:,:,k-4)
371!!!    f1=W(:,:,k-3)
372!!!    f2=W(:,:,k-2)
373!!!    f3=W(:,:,k-1)
374!!!    f4=(24.*h*BC_values6(:,:,3) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
375!!!    W(:,:,k)=f4
376    W(:,:,k) = -9999.
377   endif
378
379!!!    f0=U(:,:,k-4)
380!!!    f1=U(:,:,k-3)
381!!!    f2=U(:,:,k-2)
382!!!    f3=U(:,:,k-1)
383!!!    f4=(24.*h*BC_values6(:,:,1) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
384!!!    U(:,:,k)=f4
385    f0=U(:,:,k-2)
386    f1=U(:,:,k-1)
387    f2 = (2.*h*BC_values6(:,:,1) + 4.*f1 - f0) / 3.
388    U(:,:,k)=f2
389
390!!!    f0=v(:,:,k-4)
391!!!    f1=v(:,:,k-3)
392!!!    f2=v(:,:,k-2)
393!!!    f3=v(:,:,k-1)
394!!!    f4=(24.*h*BC_values6(:,:,2) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
395!!!    v(:,:,k)=f4      ! set dV/dz everywhere on top surface
396    f0=v(:,:,k-2)
397    f1=v(:,:,k-1)
398    f2 = (2.*h*BC_values6(:,:,2) + 4.*f1 - f0) / 3.
399    v(:,:,k)=f2
400   endif
401!**********************************************************************
402
403  if( BC_scalar6(1,1,1) == 'specified_value' ) then
404   s1(1:nx,1:ny,k)=BC_values6(1:nx,1:ny,4)
405  elseif( BC_scalar6(1,1,1) == 'specified_gradient' ) then
406!!!   f0=s1(:,:,k-4)
407!!!   f1=s1(:,:,k-3)
408!!!   f2=s1(:,:,k-2)
409!!!   f3=s1(:,:,k-1)
410!!!   f4=(24.*h*BC_values6(:,:,4) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
411!!!   s1(:,:,k)=f4
412   f0=s1(:,:,k-2)
413   f1=s1(:,:,k-1)
414   f2 = (2.*h*BC_values6(:,:,4) + 4.*f1 - f0) / 3.
415   s1(:,:,k)=f2
416  endif
417
418 if(num_scalars>1) then
419  if( BC_scalar6(1,1,2) == 'specified_value' ) then
420   s2(1:nx,1:ny,k)=BC_values6(1:nx,1:ny,5)
421  elseif( BC_scalar6(1,1,2) == 'specified_gradient' ) then
422!!!   f0=s2(:,:,k-4)
423!!!   f1=s2(:,:,k-3)
424!!!   f2=s2(:,:,k-2)
425!!!   f3=s2(:,:,k-1)
426!!!   f4=(24.*h*BC_values6(:,:,5) -6.*f0 +32.*f1 -72.*f2 +96.*f3)/50.
427!!!   s2(:,:,k)=f4
428   f0=s2(:,:,k-2)
429   f1=s2(:,:,k-1)
430   f2 = (2.*h*BC_values6(:,:,5) + 4.*f1 - f0) / 3.
431   s2(:,:,k)=f2
432  endif
433 endif
434
435  deallocate (f0,f1,f2,f3,f4)
436 endif
437
438end subroutine apply_BCs
Note: See TracBrowser for help on using the repository browser.