Context Navigation

Back to Ticket #805

Ticket #805: u2dbc_im.F

File u2dbc_im.F, 51.8 KB (added by m.hadfield, 5 years ago)

Line
1	#include "cppdefs.h"
2	MODULE u2dbc_mod
3	!
4	!svn $Id$
5	!=======================================================================
6	! Copyright (c) 2002-2019 The ROMS/TOMS Group !
7	! Licensed under a MIT/X style license !
8	! See License_ROMS.txt Hernan G. Arango !
9	!========================================== Alexander F. Shchepetkin ===
10	! !
11	! This subroutine sets lateral boundary conditions for vertically !
12	! integrated U-velocity. !
13	! !
14	!=======================================================================
15	!
16	implicit none
17	!
18	PRIVATE
19	PUBLIC :: u2dbc, u2dbc_tile
20	!
21	CONTAINS
22	!
23	!***********************************************************************
24	SUBROUTINE u2dbc (ng, tile, kout)
25	!***********************************************************************
26	!
27	USE mod_param
28	USE mod_ocean
29	USE mod_stepping
30	!
31	! Imported variable declarations.
32	!
33	integer, intent(in) :: ng, tile, kout
34	!
35	! Local variable declarations.
36	!
37	#include "tile.h"
38	!
39	CALL u2dbc_tile (ng, tile, &
40	& LBi, UBi, LBj, UBj, &
41	& IminS, ImaxS, JminS, JmaxS, &
42	& krhs(ng), kstp(ng), kout, &
43	& OCEAN(ng) % ubar, &
44	& OCEAN(ng) % vbar, &
45	& OCEAN(ng) % zeta)
46
47	RETURN
48	END SUBROUTINE u2dbc
49	!
50	!***********************************************************************
51	SUBROUTINE u2dbc_tile (ng, tile, &
52	& LBi, UBi, LBj, UBj, &
53	& IminS, ImaxS, JminS, JmaxS, &
54	& krhs, kstp, kout, &
55	& ubar, vbar, zeta)
56	!***********************************************************************
57	!
58	USE mod_param
59	USE mod_boundary
60	USE mod_clima
61	USE mod_forces
62	USE mod_grid
63	USE mod_ncparam
64	#ifdef NESTING
65	USE mod_nesting
66	#endif
67	USE mod_scalars
68	!
69	! Imported variable declarations.
70	!
71	integer, intent(in) :: ng, tile
72	integer, intent(in) :: LBi, UBi, LBj, UBj
73	integer, intent(in) :: IminS, ImaxS, JminS, JmaxS
74	integer, intent(in) :: krhs, kstp, kout
75	!
76	#ifdef ASSUMED_SHAPE
77	real(r8), intent(in) :: vbar(LBi:,LBj:,:)
78	real(r8), intent(in) :: zeta(LBi:,LBj:,:)
79
80	real(r8), intent(inout) :: ubar(LBi:,LBj:,:)
81	#else
82	real(r8), intent(in) :: vbar(LBi:UBi,LBj:UBj,3)
83	real(r8), intent(in) :: zeta(LBi:UBi,LBj:UBj,3)
84
85	real(r8), intent(inout) :: ubar(LBi:UBi,LBj:UBj,3)
86	#endif
87	!
88	! Local variable declarations.
89	!
90	integer :: Imin, Imax
91	integer :: i, j, know
92	#ifdef NESTING
93	integer :: Idg, Jdg, cr, dg, m, rg, tnew, told
94	#endif
95
96	real(r8), parameter :: eps = 1.0E-20_r8
97
98	real(r8) :: Ce, Cx, Zx
99	real(r8) :: bry_pgr, bry_cor, bry_str, bry_val
100	real(r8) :: cff, cff1, cff2, cff3, dt2d, dUde, dUdt, dUdx
101	real(r8) :: obc_in, obc_out, phi, tau
102
103	# if defined ATM_PRESS && defined PRESS_COMPENSATE
104	real(r8) :: OneAtm, fac
105	# endif
106
107	real(r8), dimension(IminS:ImaxS,JminS:JmaxS) :: grad
108
109	#include "set_bounds.h"
110	!
111	!-----------------------------------------------------------------------
112	! Set time-indices
113	!-----------------------------------------------------------------------
114	!
115	IF (FIRST_2D_STEP) THEN
116	know=krhs
117	dt2d=dtfast(ng)
118	ELSE IF (PREDICTOR_2D_STEP(ng)) THEN
119	know=krhs
120	dt2d=2.0_r8*dtfast(ng)
121	ELSE
122	know=kstp
123	dt2d=dtfast(ng)
124	END IF
125	# if defined ATM_PRESS && defined PRESS_COMPENSATE
126	OneAtm=1013.25_r8 ! 1 atm = 1013.25 mb
127	fac=100.0_r8/(g*rho0)
128	# endif
129
130	!
131	!-----------------------------------------------------------------------
132	! Lateral boundary conditions at the western edge.
133	!-----------------------------------------------------------------------
134	!
135	IF (DOMAIN(ng)%Western_Edge(tile)) THEN
136	!
137	! Western edge, implicit upstream radiation condition.
138	!
139	IF (LBC(iwest,isUbar,ng)%radiation) THEN
140	DO j=Jstr,Jend+1
141	grad(Istr ,j)=ubar(Istr ,j ,know)- &
142	& ubar(Istr ,j-1,know)
143	grad(Istr+1,j)=ubar(Istr+1,j ,know)- &
144	& ubar(Istr+1,j-1,know)
145	END DO
146	DO j=Jstr,Jend
147	IF (LBC_apply(ng)%west(j)) THEN
148	dUdt=ubar(Istr+1,j,know)-ubar(Istr+1,j,kout)
149	dUdx=ubar(Istr+1,j,kout)-ubar(Istr+2,j,kout)
150
151	IF (LBC(iwest,isUbar,ng)%nudging) THEN
152	IF (LnudgeM2CLM(ng)) THEN
153	obc_out=0.5_r8* &
154	& (CLIMA(ng)%M2nudgcof(Istr-1,j)+ &
155	& CLIMA(ng)%M2nudgcof(Istr ,j))
156	obc_in =obcfac(ng)*obc_out
157	ELSE
158	obc_out=M2obc_out(ng,iwest)
159	obc_in =M2obc_in (ng,iwest)
160	END IF
161	IF ((dUdt*dUdx).lt.0.0_r8) THEN
162	tau=obc_in
163	ELSE
164	tau=obc_out
165	END IF
166	#ifdef IMPLICIT_NUDGING
167	IF (tau.gt.0.0_r8) tau=1.0_r8/tau
168	#else
169	tau=tau*dt2d
170	#endif
171	END IF
172
173	IF ((dUdt*dUdx).lt.0.0_r8) dUdt=0.0_r8
174	IF ((dUdt*(grad(Istr+1,j )+ &
175	& grad(Istr+1,j+1))).gt.0.0_r8) THEN
176	dUde=grad(Istr+1,j )
177	ELSE
178	dUde=grad(Istr+1,j+1)
179	END IF
180	cff=MAX(dUdxdUdx+dUdedUde,eps)
181	Cx=dUdt*dUdx
182	#ifdef RADIATION_2D
183	Ce=MIN(cff,MAX(dUdt*dUde,-cff))
184	#else
185	Ce=0.0_r8
186	#endif
187	#if defined CELERITY_WRITE && defined FORWARD_WRITE
188	BOUNDARY(ng)%ubar_west_Cx(j)=Cx
189	BOUNDARY(ng)%ubar_west_Ce(j)=Ce
190	BOUNDARY(ng)%ubar_west_C2(j)=cff
191	#endif
192	ubar(Istr,j,kout)=(cff*ubar(Istr ,j,know)+ &
193	& Cx *ubar(Istr+1,j,kout)- &
194	& MAX(Ce,0.0_r8)*grad(Istr,j )- &
195	& MIN(Ce,0.0_r8)*grad(Istr,j+1))/ &
196	& (cff+Cx)
197
198	IF (LBC(iwest,isUbar,ng)%nudging) THEN
199	#ifdef IMPLICIT_NUDGING
200	phi=dt(ng)/(tau+dt(ng))
201	ubar(Istr,j,kout)=(1.0_r8-phi)*ubar(Istr,j,kout)+ &
202	& phi*BOUNDARY(ng)%ubar_west(j)
203	#else
204	ubar(Istr,j,kout)=ubar(Istr,j,kout)+ &
205	& tau*(BOUNDARY(ng)%ubar_west(j)- &
206	& ubar(Istr,j,know))
207	#endif
208	END IF
209	#ifdef MASKING
210	ubar(Istr,j,kout)=ubar(Istr,j,kout)* &
211	& GRID(ng)%umask(Istr,j)
212	#endif
213	END IF
214	END DO
215	!
216	! Western edge, Flather boundary condition.
217	!
218	ELSE IF (LBC(iwest,isUbar,ng)%Flather) THEN
219	DO j=Jstr,Jend
220	IF (LBC_apply(ng)%west(j)) THEN
221	#if defined SSH_TIDES && !defined UV_TIDES
222	IF (LBC(iwest,isFsur,ng)%acquire) THEN
223	bry_pgr=-g*(zeta(Istr,j,know)- &
224	& BOUNDARY(ng)%zeta_west(j))* &
225	& 0.5_r8*GRID(ng)%pm(Istr,j)
226	ELSE
227	bry_pgr=-g*(zeta(Istr ,j,know)- &
228	& zeta(Istr-1,j,know))* &
229	& 0.5_r8*(GRID(ng)%pm(Istr-1,j)+ &
230	& GRID(ng)%pm(Istr ,j))
231	END IF
232	# ifdef UV_COR
233	bry_cor=0.125_r8*(vbar(Istr-1,j ,know)+ &
234	& vbar(Istr-1,j+1,know)+ &
235	& vbar(Istr ,j ,know)+ &
236	& vbar(Istr ,j+1,know))* &
237	& (GRID(ng)%f(Istr-1,j)+ &
238	& GRID(ng)%f(Istr ,j))
239	# else
240	bry_cor=0.0_r8
241	# endif
242	cff1=1.0_r8/(0.5_r8*(GRID(ng)%h(Istr-1,j)+ &
243	& zeta(Istr-1,j,know)+ &
244	& GRID(ng)%h(Istr ,j)+ &
245	& zeta(Istr ,j,know)))
246	bry_str=cff1*(FORCES(ng)%sustr(Istr,j)- &
247	& FORCES(ng)%bustr(Istr,j))
248	Cx=1.0_r8/SQRT(g0.5_r8(GRID(ng)%h(Istr-1,j)+ &
249	& zeta(Istr-1,j,know)+ &
250	& GRID(ng)%h(Istr ,j)+ &
251	& zeta(Istr ,j,know)))
252	cff2=GRID(ng)%om_u(Istr,j)*Cx
253	!! cff2=dt2d
254	bry_val=ubar(Istr+1,j,know)+ &
255	& cff2*(bry_pgr+ &
256	& bry_cor+ &
257	& bry_str)
258	#else
259	bry_val=BOUNDARY(ng)%ubar_west(j)
260	#endif
261	cff=1.0_r8/(0.5_r8*(GRID(ng)%h(Istr-1,j)+ &
262	& zeta(Istr-1,j,know)+ &
263	& GRID(ng)%h(Istr ,j)+ &
264	& zeta(Istr ,j,know)))
265	Cx=SQRT(g*cff)
266	# if defined ATM_PRESS && defined PRESS_COMPENSATE
267	ubar(Istr,j,kout)= &
268	& bry_val-Cx(0.5_r8(zeta(Istr-1,j,know)+ &
269	& zeta(Istr ,j,know)+ &
270	& fac*(FORCES(ng)%Pair(Istr-1,j)+ &
271	& FORCES(ng)%Pair(Istr ,j)- &
272	& 2.0_r8*OneAtm))- &
273	& BOUNDARY(ng)%zeta_west(j))
274	# else
275	ubar(Istr,j,kout)=bry_val- &
276	& Cx(0.5_r8(zeta(Istr-1,j,know)+ &
277	& zeta(Istr ,j,know))- &
278	& BOUNDARY(ng)%zeta_west(j))
279	# endif
280	#ifdef MASKING
281	ubar(Istr,j,kout)=ubar(Istr,j,kout)* &
282	& GRID(ng)%umask(Istr,j)
283	#endif
284	END IF
285	END DO
286	!
287	! Western edge, Shchepetkin boundary condition (Maison et al., 2010).
288	!
289	ELSE IF (LBC(iwest,isUbar,ng)%Shchepetkin) THEN
290	DO j=Jstr,Jend
291	IF (LBC_apply(ng)%west(j)) THEN
292	#if defined SSH_TIDES && !defined UV_TIDES
293	IF (LBC(iwest,isFsur,ng)%acquire) THEN
294	bry_pgr=-g*(zeta(Istr,j,know)- &
295	& BOUNDARY(ng)%zeta_west(j))* &
296	& 0.5_r8*GRID(ng)%pm(Istr,j)
297	ELSE
298	bry_pgr=-g*(zeta(Istr ,j,know)- &
299	& zeta(Istr-1,j,know))* &
300	& 0.5_r8*(GRID(ng)%pm(Istr-1,j)+ &
301	& GRID(ng)%pm(Istr ,j))
302	END IF
303	# ifdef UV_COR
304	bry_cor=0.125_r8*(vbar(Istr-1,j ,know)+ &
305	& vbar(Istr-1,j+1,know)+ &
306	& vbar(Istr ,j ,know)+ &
307	& vbar(Istr ,j+1,know))* &
308	& (GRID(ng)%f(Istr-1,j)+ &
309	& GRID(ng)%f(Istr ,j))
310	# else
311	bry_cor=0.0_r8
312	# endif
313	cff1=1.0_r8/(0.5_r8*(GRID(ng)%h(Istr-1,j)+ &
314	& zeta(Istr-1,j,know)+ &
315	& GRID(ng)%h(Istr ,j)+ &
316	& zeta(Istr ,j,know)))
317	bry_str=cff1*(FORCES(ng)%sustr(Istr,j)- &
318	& FORCES(ng)%bustr(Istr,j))
319	Cx=1.0_r8/SQRT(g0.5_r8(GRID(ng)%h(Istr-1,j)+ &
320	& zeta(Istr-1,j,know)+ &
321	& GRID(ng)%h(Istr ,j)+ &
322	& zeta(Istr ,j,know)))
323	cff2=GRID(ng)%om_u(Istr,j)*Cx
324	!! cff2=dt2d
325	bry_val=ubar(Istr+1,j,know)+ &
326	& cff2*(bry_pgr+ &
327	& bry_cor+ &
328	& bry_str)
329	#else
330	bry_val=BOUNDARY(ng)%ubar_west(j)
331	#endif
332	cff=0.5_r8*(GRID(ng)%h(Istr-1,j)+ &
333	& GRID(ng)%h(Istr ,j))
334	cff1=SQRT(g/cff)
335	Cx=dt2dcff1cff0.5_r8(GRID(ng)%pm(Istr-1,j)+ &
336	& GRID(ng)%pm(Istr ,j))
337	Zx=(0.5_r8+Cx)*zeta(Istr ,j,know)+ &
338	& (0.5_r8-Cx)*zeta(Istr-1,j,know)
339	IF (Cx.gt.Co) THEN
340	cff2=(1.0_r8-Co/Cx)**2
341	cff3=zeta(Istr,j,kout)+ &
342	& Cx*zeta(Istr-1,j,know)- &
343	& (1.0_r8+Cx)*zeta(Istr,j,know)
344	Zx=Zx+cff2*cff3
345	END IF
346	ubar(Istr,j,kout)=0.5_r8* &
347	& ((1.0_r8-Cx)*ubar(Istr,j,know)+ &
348	& Cx*ubar(Istr+1,j,know)+ &
349	& bry_val- &
350	& cff1*(Zx-BOUNDARY(ng)%zeta_west(j)))
351	#ifdef MASKING
352	ubar(Istr,j,kout)=ubar(Istr,j,kout)* &
353	& GRID(ng)%umask(Istr,j)
354	#endif
355	END IF
356	END DO
357	!
358	! Western edge, clamped boundary condition.
359	!
360	ELSE IF (LBC(iwest,isUbar,ng)%clamped) THEN
361	DO j=Jstr,Jend
362	IF (LBC_apply(ng)%west(j)) THEN
363	ubar(Istr,j,kout)=BOUNDARY(ng)%ubar_west(j)
364	#ifdef MASKING
365	ubar(Istr,j,kout)=ubar(Istr,j,kout)* &
366	& GRID(ng)%umask(Istr,j)
367	#endif
368	END IF
369	END DO
370	!
371	! Western edge, gradient boundary condition.
372	!
373	ELSE IF (LBC(iwest,isUbar,ng)%gradient) THEN
374	DO j=Jstr,Jend
375	IF (LBC_apply(ng)%west(j)) THEN
376	ubar(Istr,j,kout)=ubar(Istr+1,j,kout)
377	#ifdef MASKING
378	ubar(Istr,j,kout)=ubar(Istr,j,kout)* &
379	& GRID(ng)%umask(Istr,j)
380	#endif
381	END IF
382	END DO
383	!
384	! Western edge, reduced-physics boundary condition.
385	!
386	ELSE IF (LBC(iwest,isUbar,ng)%reduced) THEN
387	DO j=Jstr,Jend
388	IF (LBC_apply(ng)%west(j)) THEN
389	IF (LBC(iwest,isFsur,ng)%acquire) THEN
390	bry_pgr=-g*(zeta(Istr,j,know)- &
391	& BOUNDARY(ng)%zeta_west(j))* &
392	& 0.5_r8*GRID(ng)%pm(Istr,j)
393	ELSE
394	bry_pgr=-g*(zeta(Istr ,j,know)- &
395	& zeta(Istr-1,j,know))* &
396	& 0.5_r8*(GRID(ng)%pm(Istr-1,j)+ &
397	& GRID(ng)%pm(Istr ,j))
398	END IF
399	#ifdef UV_COR
400	bry_cor=0.125_r8*(vbar(Istr-1,j ,know)+ &
401	& vbar(Istr-1,j+1,know)+ &
402	& vbar(Istr ,j ,know)+ &
403	& vbar(Istr ,j+1,know))* &
404	& (GRID(ng)%f(Istr-1,j)+ &
405	& GRID(ng)%f(Istr ,j))
406	#else
407	bry_cor=0.0_r8
408	#endif
409	cff=1.0_r8/(0.5_r8*(GRID(ng)%h(Istr-1,j)+ &
410	& zeta(Istr-1,j,know)+ &
411	& GRID(ng)%h(Istr ,j)+ &
412	& zeta(Istr ,j,know)))
413	bry_str=cff*(FORCES(ng)%sustr(Istr,j)- &
414	& FORCES(ng)%bustr(Istr,j))
415	ubar(Istr,j,kout)=ubar(Istr,j,know)+ &
416	& dt2d*(bry_pgr+ &
417	& bry_cor+ &
418	& bry_str)
419	#ifdef MASKING
420	ubar(Istr,j,kout)=ubar(Istr,j,kout)* &
421	& GRID(ng)%umask(Istr,j)
422	#endif
423	END IF
424	END DO
425	!
426	! Western edge, closed boundary condition.
427	!
428	ELSE IF (LBC(iwest,isUbar,ng)%closed) THEN
429	DO j=Jstr,Jend
430	IF (LBC_apply(ng)%west(j)) THEN
431	ubar(Istr,j,kout)=0.0_r8
432	END IF
433	END DO
434
435	#ifdef NESTING
436	!
437	! If refinement grid and western edge, impose mass flux from donor
438	! coaser grid for volume and mass conservation.
439	!
440	ELSE IF (LBC(iwest,isUbar,ng)%nested) THEN
441	DO cr=1,Ncontact
442	dg=Ucontact(cr)%donor_grid
443	rg=Ucontact(cr)%receiver_grid
444	IF (RefinedGrid(ng).and. &
445	& (rg.eq.ng).and.(DXmax(dg).gt.DXmax(rg))) THEN
446	told=3-RollingIndex(cr)
447	tnew=RollingIndex(cr)
448	DO j=Jstr,Jend
449	m=BRY_CONTACT(iwest,cr)%C2Bindex(j)
450	Idg=Ucontact(cr)%Idg(m) ! for debugging
451	Jdg=Ucontact(cr)%Jdg(m) ! purposes
452	cff=0.5_r8GRID(ng)%on_u(Istr,j) &
453	& (GRID(ng)%h(Istr-1,j)+zeta(Istr-1,j,kout)+ &
454	& GRID(ng)%h(Istr ,j)+zeta(Istr ,j,kout))
455	cff1=GRID(ng)%on_u(Istr,j)/REFINED(cr)%on_u(m)
456	bry_val=cff1*REFINED(cr)%DU_avg2(1,m,tnew)/cff
457	# ifdef MASKING
458	bry_val=bry_val*GRID(ng)%umask(Istr,j)
459	# endif
460	# ifdef NESTING_DEBUG
461	BRY_CONTACT(iwest,cr)%Mflux(j)=cff*bry_val
462	# endif
463	ubar(Istr,j,kout)=bry_val
464	END DO
465	END IF
466	END DO
467	#endif
468	END IF
469	END IF
470	!
471	!-----------------------------------------------------------------------
472	! Lateral boundary conditions at the eastern edge.
473	!-----------------------------------------------------------------------
474	!
475	IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
476	!
477	! Eastern edge, implicit upstream radiation condition.
478	!
479	IF (LBC(ieast,isUbar,ng)%radiation) THEN
480	DO j=Jstr,Jend+1
481	grad(Iend ,j)=ubar(Iend ,j ,know)- &
482	& ubar(Iend ,j-1,know)
483	grad(Iend+1,j)=ubar(Iend+1,j ,know)- &
484	& ubar(Iend+1,j-1,know)
485	END DO
486	DO j=Jstr,Jend
487	IF (LBC_apply(ng)%east(j)) THEN
488	dUdt=ubar(Iend,j,know)-ubar(Iend ,j,kout)
489	dUdx=ubar(Iend,j,kout)-ubar(Iend-1,j,kout)
490
491	IF (LBC(ieast,isUbar,ng)%nudging) THEN
492	IF (LnudgeM2CLM(ng)) THEN
493	obc_out=0.5_r8* &
494	& (CLIMA(ng)%M2nudgcof(Iend ,j)+ &
495	& CLIMA(ng)%M2nudgcof(Iend+1,j))
496	obc_in =obcfac(ng)*obc_out
497	ELSE
498	obc_out=M2obc_out(ng,ieast)
499	obc_in =M2obc_in (ng,ieast)
500	END IF
501	IF ((dUdt*dUdx).lt.0.0_r8) THEN
502	tau=obc_in
503	ELSE
504	tau=obc_out
505	END IF
506	#ifdef IMPLICIT_NUDGING
507	IF (tau.gt.0.0_r8) tau=1.0_r8/tau
508	#else
509	tau=tau*dt2d
510	#endif
511	END IF
512
513	IF ((dUdt*dUdx).lt.0.0_r8) dUdt=0.0_r8
514	IF ((dUdt*(grad(Iend,j )+ &
515	& grad(Iend,j+1))).gt.0.0_r8) THEN
516	dUde=grad(Iend,j)
517	ELSE
518	dUde=grad(Iend,j+1)
519	END IF
520	cff=MAX(dUdxdUdx+dUdedUde,eps)
521	Cx=dUdt*dUdx
522	#ifdef RADIATION_2D
523	Ce=MIN(cff,MAX(dUdt*dUde,-cff))
524	#else
525	Ce=0.0_r8
526	#endif
527	#if defined CELERITY_WRITE && defined FORWARD_WRITE
528	BOUNDARY(ng)%ubar_east_Cx(j)=Cx
529	BOUNDARY(ng)%ubar_east_Ce(j)=Ce
530	BOUNDARY(ng)%ubar_east_C2(j)=cff
531	#endif
532	ubar(Iend+1,j,kout)=(cff*ubar(Iend+1,j,know)+ &
533	& Cx *ubar(Iend ,j,kout)- &
534	& MAX(Ce,0.0_r8)*grad(Iend+1,j )- &
535	& MIN(Ce,0.0_r8)*grad(Iend+1,j+1))/ &
536	& (cff+Cx)
537
538	IF (LBC(ieast,isUbar,ng)%nudging) THEN
539	#ifdef IMPLICIT_NUDGING
540	phi=dt(ng)/(tau+dt(ng))
541	ubar(Iend+1,j,kout)=(1.0_r8-phi)*ubar(Iend+1,j,kout)+ &
542	& phi*BOUNDARY(ng)%ubar_east(j)
543	#else
544	ubar(Iend+1,j,kout)=ubar(Iend+1,j,kout)+ &
545	& tau*(BOUNDARY(ng)%ubar_east(j)- &
546	& ubar(Iend+1,j,know))
547	#endif
548	END IF
549	#ifdef MASKING
550	ubar(Iend+1,j,kout)=ubar(Iend+1,j,kout)* &
551	& GRID(ng)%umask(Iend+1,j)
552	#endif
553	END IF
554	END DO
555	!
556	! Eastern edge, Flather boundary condition.
557	!
558	ELSE IF (LBC(ieast,isUbar,ng)%Flather) THEN
559	DO j=Jstr,Jend
560	IF (LBC_apply(ng)%east(j)) THEN
561	#if defined SSH_TIDES && !defined UV_TIDES
562	IF (LBC(ieast,isFsur,ng)%acquire) THEN
563	bry_pgr=-g*(BOUNDARY(ng)%zeta_east(j)- &
564	& zeta(Iend,j,know))* &
565	& 0.5_r8*GRID(ng)%pm(Iend,j)
566	ELSE
567	bry_pgr=-g*(zeta(Iend+1,j,know)- &
568	& zeta(Iend ,j,know))* &
569	& 0.5_r8*(GRID(ng)%pm(Iend ,j)+ &
570	& GRID(ng)%pm(Iend+1,j))
571	END IF
572	# ifdef UV_COR
573	bry_cor=0.125_r8*(vbar(Iend ,j ,know)+ &
574	& vbar(Iend ,j+1,know)+ &
575	& vbar(Iend+1,j ,know)+ &
576	& vbar(Iend+1,j+1,know))* &
577	& (GRID(ng)%f(Iend ,j)+ &
578	& GRID(ng)%f(Iend+1,j))
579	# else
580	bry_cor=0.0_r8
581	# endif
582	cff1=1.0_r8/(0.5_r8*(GRID(ng)%h(Iend ,j)+ &
583	& zeta(Iend ,j,know)+ &
584	& GRID(ng)%h(Iend+1,j)+ &
585	& zeta(Iend+1,j,know)))
586	bry_str=cff1*(FORCES(ng)%sustr(Iend+1,j)- &
587	& FORCES(ng)%bustr(Iend+1,j))
588	Cx=1.0_r8/SQRT(g0.5_r8(GRID(ng)%h(Iend+1,j)+ &
589	& zeta(Iend+1,j,know)+ &
590	& GRID(ng)%h(Iend ,j)+ &
591	& zeta(Iend ,j,know)))
592	cff2=GRID(ng)%om_u(Iend+1,j)*Cx
593	!! cff2=dt2d
594	bry_val=ubar(Iend,j,know)+ &
595	& cff2*(bry_pgr+ &
596	& bry_cor+ &
597	& bry_str)
598	#else
599	bry_val=BOUNDARY(ng)%ubar_east(j)
600	#endif
601	cff=1.0_r8/(0.5_r8*(GRID(ng)%h(Iend ,j)+ &
602	& zeta(Iend ,j,know)+ &
603	& GRID(ng)%h(Iend+1,j)+ &
604	& zeta(Iend+1,j,know)))
605	Cx=SQRT(g*cff)
606	# if defined ATM_PRESS && defined PRESS_COMPENSATE
607	ubar(Iend+1,j,kout)= &
608	& bry_val+Cx(0.5_r8(zeta(Iend ,j,know)+ &
609	& zeta(Iend+1,j,know)+ &
610	& fac*(FORCES(ng)%Pair(Iend ,j)+ &
611	& FORCES(ng)%Pair(Iend+1,j)- &
612	& 2.0_r8*OneAtm))- &
613	& BOUNDARY(ng)%zeta_east(j))
614	# else
615	ubar(Iend+1,j,kout)=bry_val+ &
616	& Cx(0.5_r8(zeta(Iend ,j,know)+ &
617	& zeta(Iend+1,j,know))- &
618	& BOUNDARY(ng)%zeta_east(j))
619	# endif
620	#ifdef MASKING
621	ubar(Iend+1,j,kout)=ubar(Iend+1,j,kout)* &
622	& GRID(ng)%umask(Iend+1,j)
623	#endif
624	END IF
625	END DO
626	!
627	! Eastern edge, Shchepetkin boundary condition (Maison et al., 2010).
628	!
629	ELSE IF (LBC(ieast,isUbar,ng)%Shchepetkin) THEN
630	DO j=Jstr,Jend
631	IF (LBC_apply(ng)%east(j)) THEN
632	#if defined SSH_TIDES && !defined UV_TIDES
633	IF (LBC(ieast,isFsur,ng)%acquire) THEN
634	bry_pgr=-g*(BOUNDARY(ng)%zeta_east(j)- &
635	& zeta(Iend,j,know))* &
636	& 0.5_r8*GRID(ng)%pm(Iend,j)
637	ELSE
638	bry_pgr=-g*(zeta(Iend+1,j,know)- &
639	& zeta(Iend ,j,know))* &
640	& 0.5_r8*(GRID(ng)%pm(Iend ,j)+ &
641	& GRID(ng)%pm(Iend+1,j))
642	END IF
643	# ifdef UV_COR
644	bry_cor=0.125_r8*(vbar(Iend ,j ,know)+ &
645	& vbar(Iend ,j+1,know)+ &
646	& vbar(Iend+1,j ,know)+ &
647	& vbar(Iend+1,j+1,know))* &
648	& (GRID(ng)%f(Iend ,j)+ &
649	& GRID(ng)%f(Iend+1,j))
650	# else
651	bry_cor=0.0_r8
652	# endif
653	cff1=1.0_r8/(0.5_r8*(GRID(ng)%h(Iend ,j)+ &
654	& zeta(Iend ,j,know)+ &
655	& GRID(ng)%h(Iend+1,j)+ &
656	& zeta(Iend+1,j,know)))
657	bry_str=cff1*(FORCES(ng)%sustr(Iend+1,j)- &
658	& FORCES(ng)%bustr(Iend+1,j))
659	Cx=1.0_r8/SQRT(g0.5_r8(GRID(ng)%h(Iend+1,j)+ &
660	& zeta(Iend+1,j,know)+ &
661	& GRID(ng)%h(Iend ,j)+ &
662	& zeta(Iend ,j,know)))
663	cff2=GRID(ng)%om_u(Iend+1,j)*Cx
664	!! cff2=dt2d
665	bry_val=ubar(Iend,j,know)+ &
666	& cff2*(bry_pgr+ &
667	& bry_cor+ &
668	& bry_str)
669	#else
670	bry_val=BOUNDARY(ng)%ubar_east(j)
671	#endif
672	cff=0.5_r8*(GRID(ng)%h(Iend ,j)+ &
673	& GRID(ng)%h(Iend+1,j))
674	cff1=SQRT(g/cff)
675	Cx=dt2dcff1cff0.5_r8(GRID(ng)%pm(Iend ,j)+ &
676	& GRID(ng)%pm(Iend+1,j))
677	Zx=(0.5_r8+Cx)*zeta(Iend ,j,know)+ &
678	& (0.5_r8-Cx)*zeta(Iend+1,j,know)
679	IF (Cx.gt.Co) THEN
680	cff2=(1.0_r8-Co/Cx)**2
681	cff3=zeta(Iend,j,kout)+ &
682	& Cx*zeta(Iend+1,j,know)- &
683	& (1.0_r8+Cx)*zeta(Iend,j,know)
684	Zx=Zx+cff2*cff3
685	END IF
686	ubar(Iend+1,j,kout)=0.5_r8* &
687	& ((1.0_r8-Cx)*ubar(Iend+1,j,know)+ &
688	& Cx*ubar(Iend,j,know)+ &
689	& bry_val+ &
690	& cff1*(Zx-BOUNDARY(ng)%zeta_east(j)))
691	#ifdef MASKING
692	ubar(Iend+1,j,kout)=ubar(Iend+1,j,kout)* &
693	& GRID(ng)%umask(Iend+1,j)
694	#endif
695	END IF
696	END DO
697	!
698	! Eastern edge, clamped boundary condition.
699	!
700	ELSE IF (LBC(ieast,isUbar,ng)%clamped) THEN
701	DO j=Jstr,Jend
702	IF (LBC_apply(ng)%east(j)) THEN
703	ubar(Iend+1,j,kout)=BOUNDARY(ng)%ubar_east(j)
704	#ifdef MASKING
705	ubar(Iend+1,j,kout)=ubar(Iend+1,j,kout)* &
706	& GRID(ng)%umask(Iend+1,j)
707	#endif
708	END IF
709	END DO
710	!
711	! Eastern edge, gradient boundary condition.
712	!
713	ELSE IF (LBC(ieast,isUbar,ng)%gradient) THEN
714	DO j=Jstr,Jend
715	IF (LBC_apply(ng)%east(j)) THEN
716	ubar(Iend+1,j,kout)=ubar(Iend,j,kout)
717	#ifdef MASKING
718	ubar(Iend+1,j,kout)=ubar(Iend+1,j,kout)* &
719	& GRID(ng)%umask(Iend+1,j)
720	#endif
721	END IF
722	END DO
723	!
724	! Eastern edge, reduced-physics boundary condition.
725	!
726	ELSE IF (LBC(ieast,isUbar,ng)%reduced) THEN
727	DO j=Jstr,Jend
728	IF (LBC_apply(ng)%east(j)) THEN
729	IF (LBC(ieast,isFsur,ng)%acquire) THEN
730	bry_pgr=-g*(BOUNDARY(ng)%zeta_east(j)- &
731	& zeta(Iend,j,know))* &
732	& 0.5_r8*GRID(ng)%pm(Iend,j)
733	ELSE
734	bry_pgr=-g*(zeta(Iend+1,j,know)- &
735	& zeta(Iend ,j,know))* &
736	& 0.5_r8*(GRID(ng)%pm(Iend ,j)+ &
737	& GRID(ng)%pm(Iend+1,j))
738	END IF
739	#ifdef UV_COR
740	bry_cor=0.125_r8*(vbar(Iend ,j ,know)+ &
741	& vbar(Iend ,j+1,know)+ &
742	& vbar(Iend+1,j ,know)+ &
743	& vbar(Iend+1,j+1,know))* &
744	& (GRID(ng)%f(Iend ,j)+ &
745	& GRID(ng)%f(Iend+1,j))
746	#else
747	bry_cor=0.0_r8
748	#endif
749	cff=1.0_r8/(0.5_r8*(GRID(ng)%h(Iend ,j)+ &
750	& zeta(Iend ,j,know)+ &
751	& GRID(ng)%h(Iend+1,j)+ &
752	& zeta(Iend+1,j,know)))
753	bry_str=cff*(FORCES(ng)%sustr(Iend+1,j)- &
754	& FORCES(ng)%bustr(Iend+1,j))
755	ubar(Iend+1,j,kout)=ubar(Iend+1,j,know)+ &
756	& dt2d*(bry_pgr+ &
757	& bry_cor+ &
758	& bry_str)
759	#ifdef MASKING
760	ubar(Iend+1,j,kout)=ubar(Iend+1,j,kout)* &
761	& GRID(ng)%umask(Iend+1,j)
762	#endif
763	END IF
764	END DO
765	!
766	! Eastern edge, closed boundary condition.
767	!
768	ELSE IF (LBC(ieast,isUbar,ng)%closed) THEN
769	DO j=Jstr,Jend
770	IF (LBC_apply(ng)%east(j)) THEN
771	ubar(Iend+1,j,kout)=0.0_r8
772	END IF
773	END DO
774
775	#ifdef NESTING
776	!
777	! If refinement grid and eastern edge, impose mass flux from donor
778	! coaser grid for volume and mass conservation.
779	!
780	ELSE IF (LBC(ieast,isUbar,ng)%nested) THEN
781	DO cr=1,Ncontact
782	dg=Ucontact(cr)%donor_grid
783	rg=Ucontact(cr)%receiver_grid
784	IF (RefinedGrid(ng).and. &
785	& (rg.eq.ng).and.(DXmax(dg).gt.DXmax(rg))) THEN
786	told=3-RollingIndex(cr)
787	tnew=RollingIndex(cr)
788	DO j=Jstr,Jend
789	m=BRY_CONTACT(ieast,cr)%C2Bindex(j)
790	Idg=Ucontact(cr)%Idg(m) ! for debugging
791	Jdg=Ucontact(cr)%Jdg(m) ! purposes
792	cff=0.5_r8GRID(ng)%on_u(Iend+1,j) &
793	& (GRID(ng)%h(Iend+1,j)+zeta(Iend+1,j,kout)+ &
794	& GRID(ng)%h(Iend ,j)+zeta(Iend ,j,kout))
795	cff1=GRID(ng)%on_u(Iend+1,j)/REFINED(cr)%on_u(m)
796	bry_val=cff1*REFINED(cr)%DU_avg2(1,m,tnew)/cff
797	# ifdef MASKING
798	bry_val=bry_val*GRID(ng)%umask(Iend+1,j)
799	# endif
800	# ifdef NESTING_DEBUG
801	BRY_CONTACT(ieast,cr)%Mflux(j)=cff*bry_val
802	# endif
803	ubar(Iend+1,j,kout)=bry_val
804	END DO
805	END IF
806	END DO
807	#endif
808	END IF
809	END IF
810	!
811	!-----------------------------------------------------------------------
812	! Lateral boundary conditions at the southern edge.
813	!-----------------------------------------------------------------------
814	!
815	IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
816	!
817	! Southern edge, implicit upstream radiation condition.
818	!
819	IF (LBC(isouth,isUbar,ng)%radiation) THEN
820	DO i=IstrU-1,Iend
821	grad(i,Jstr-1)=ubar(i+1,Jstr-1,know)- &
822	& ubar(i ,Jstr-1,know)
823	grad(i,Jstr )=ubar(i+1,Jstr ,know)- &
824	& ubar(i ,Jstr ,know)
825	END DO
826	DO i=IstrU,Iend
827	IF (LBC_apply(ng)%south(i)) THEN
828	dUdt=ubar(i,Jstr,know)-ubar(i,Jstr ,kout)
829	dUde=ubar(i,Jstr,kout)-ubar(i,Jstr+1,kout)
830
831	IF (LBC(isouth,isUbar,ng)%nudging) THEN
832	IF (LnudgeM2CLM(ng)) THEN
833	obc_out=0.5_r8* &
834	& (CLIMA(ng)%M2nudgcof(i-1,Jstr-1)+ &
835	& CLIMA(ng)%M2nudgcof(i ,Jstr-1))
836	obc_in =obcfac(ng)*obc_out
837	ELSE
838	obc_out=M2obc_out(ng,isouth)
839	obc_in =M2obc_in (ng,isouth)
840	END IF
841	IF ((dUdt*dUde).lt.0.0_r8) THEN
842	tau=obc_in
843	ELSE
844	tau=obc_out
845	END IF
846	#ifdef IMPLICIT_NUDGING
847	IF (tau.gt.0.0_r8) tau=1.0_r8/tau
848	#else
849	tau=tau*dt2d
850	#endif
851	END IF
852
853	IF ((dUdt*dUde).lt.0.0_r8) dUdt=0.0_r8
854	IF ((dUdt*(grad(i-1,Jstr)+ &
855	& grad(i ,Jstr))).gt.0.0_r8) THEN
856	dUdx=grad(i-1,Jstr)
857	ELSE
858	dUdx=grad(i ,Jstr)
859	END IF
860	cff=MAX(dUdxdUdx+dUdedUde,eps)
861	#ifdef RADIATION_2D
862	Cx=MIN(cff,MAX(dUdt*dUdx,-cff))
863	#else
864	Cx=0.0_r8
865	#endif
866	Ce=dUdt*dUde
867	#if defined CELERITY_WRITE && defined FORWARD_WRITE
868	BOUNDARY(ng)%ubar_south_Cx(i)=Cx
869	BOUNDARY(ng)%ubar_south_Ce(i)=Ce
870	BOUNDARY(ng)%ubar_south_C2(i)=cff
871	#endif
872	ubar(i,Jstr-1,kout)=(cff*ubar(i,Jstr-1,know)+ &
873	& Ce *ubar(i,Jstr ,kout)- &
874	& MAX(Cx,0.0_r8)*grad(i-1,Jstr-1)- &
875	& MIN(Cx,0.0_r8)*grad(i ,Jstr-1))/ &
876	& (cff+Ce)
877
878	IF (LBC(isouth,isUbar,ng)%nudging) THEN
879	#ifdef IMPLICIT_NUDGING
880	phi=dt(ng)/(tau+dt(ng))
881	ubar(i,Jstr-1,kout)=(1.0_r8-phi)*ubar(i,Jstr-1,kout)+ &
882	& phi*BOUNDARY(ng)%ubar_south(i)
883	#else
884	ubar(i,Jstr-1,kout)=ubar(i,Jstr-1,kout)+ &
885	& tau*(BOUNDARY(ng)%ubar_south(i)- &
886	& ubar(i,Jstr-1,know))
887	#endif
888	END IF
889	#ifdef MASKING
890	ubar(i,Jstr-1,kout)=ubar(i,Jstr-1,kout)* &
891	& GRID(ng)%umask(i,Jstr-1)
892	#endif
893	END IF
894	END DO
895	!
896	! Southern edge, Chapman boundary condition.
897	!
898	ELSE IF (LBC(isouth,isUbar,ng)%Flather.or. &
899	& LBC(isouth,isUbar,ng)%reduced.or. &
900	& LBC(isouth,isUbar,ng)%Shchepetkin) THEN
901	DO i=IstrU,Iend
902	IF (LBC_apply(ng)%south(i)) THEN
903	cff=dt2d0.5_r8(GRID(ng)%pn(i-1,Jstr)+ &
904	& GRID(ng)%pn(i ,Jstr))
905	cff1=SQRT(g0.5_r8(GRID(ng)%h(i-1,Jstr)+ &
906	& zeta(i-1,Jstr,know)+ &
907	& GRID(ng)%h(i ,Jstr)+ &
908	& zeta(i ,Jstr,know)))
909	Ce=cff*cff1
910	cff2=1.0_r8/(1.0_r8+Ce)
911	ubar(i,Jstr-1,kout)=cff2*(ubar(i,Jstr-1,know)+ &
912	& Ce*ubar(i,Jstr,kout))
913	#ifdef MASKING
914	ubar(i,Jstr-1,kout)=ubar(i,Jstr-1,kout)* &
915	& GRID(ng)%umask(i,Jstr-1)
916	#endif
917	END IF
918	END DO
919	!
920	! Southern edge, clamped boundary condition.
921	!
922	ELSE IF (LBC(isouth,isUbar,ng)%clamped) THEN
923	DO i=IstrU,Iend
924	IF (LBC_apply(ng)%south(i)) THEN
925	ubar(i,Jstr-1,kout)=BOUNDARY(ng)%ubar_south(i)
926	#ifdef MASKING
927	ubar(i,Jstr-1,kout)=ubar(i,Jstr-1,kout)* &
928	& GRID(ng)%umask(i,Jstr-1)
929	#endif
930	END IF
931	END DO
932	!
933	! Southern edge, gradient boundary condition.
934	!
935	ELSE IF (LBC(isouth,isUbar,ng)%gradient) THEN
936	DO i=IstrU,Iend
937	IF (LBC_apply(ng)%south(i)) THEN
938	ubar(i,Jstr-1,kout)=ubar(i,Jstr,kout)
939	#ifdef MASKING
940	ubar(i,Jstr-1,kout)=ubar(i,Jstr-1,kout)* &
941	& GRID(ng)%umask(i,Jstr-1)
942	#endif
943	END IF
944	END DO
945	!
946	! Southern edge, closed boundary condition: free slip (gamma2=1) or
947	! no slip (gamma2=-1).
948	!
949	ELSE IF (LBC(isouth,isUbar,ng)%closed) THEN
950	IF (EWperiodic(ng)) THEN
951	Imin=IstrU
952	Imax=Iend
953	ELSE
954	Imin=Istr
955	Imax=IendR
956	END IF
957	DO i=Imin,Imax
958	IF (LBC_apply(ng)%south(i)) THEN
959	ubar(i,Jstr-1,kout)=gamma2(ng)*ubar(i,Jstr,kout)
960	#ifdef MASKING
961	ubar(i,Jstr-1,kout)=ubar(i,Jstr-1,kout)* &
962	& GRID(ng)%umask(i,Jstr-1)
963	#endif
964	END IF
965	END DO
966	END IF
967	END IF
968	!
969	!-----------------------------------------------------------------------
970	! Lateral boundary conditions at the northern edge.
971	!-----------------------------------------------------------------------
972	!
973	IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
974	!
975	! Northern edge, implicit upstream radiation condition.
976	!
977	IF (LBC(inorth,isUbar,ng)%radiation) THEN
978	DO i=IstrU-1,Iend
979	grad(i,Jend )=ubar(i+1,Jend ,know)- &
980	& ubar(i ,Jend ,know)
981	grad(i,Jend+1)=ubar(i+1,Jend+1,know)- &
982	& ubar(i ,Jend+1,know)
983	END DO
984	DO i=IstrU,Iend
985	IF (LBC_apply(ng)%north(i)) THEN
986	dUdt=ubar(i,Jend,know)-ubar(i,Jend ,kout)
987	dUde=ubar(i,Jend,kout)-ubar(i,Jend-1,kout)
988
989	IF (LBC(inorth,isUbar,ng)%nudging) THEN
990	IF (LnudgeM2CLM(ng)) THEN
991	obc_out=0.5_r8* &
992	& (CLIMA(ng)%M2nudgcof(i-1,Jend+1)+ &
993	& CLIMA(ng)%M2nudgcof(i ,Jend+1))
994	obc_in =obcfac(ng)*obc_out
995	ELSE
996	obc_out=M2obc_out(ng,inorth)
997	obc_in =M2obc_in (ng,inorth)
998	END IF
999	IF ((dUdt*dUde).lt.0.0_r8) THEN
1000	tau=obc_in
1001	ELSE
1002	tau=obc_out
1003	END IF
1004	#ifdef IMPLICIT_NUDGING
1005	IF (tau.gt.0.0_r8) tau=1.0_r8/tau
1006	#else
1007	tau=tau*dt2d
1008	#endif
1009	END IF
1010
1011	IF ((dUdt*dUde).lt.0.0_r8) dUdt=0.0_r8
1012	IF ((dUdt*(grad(i-1,Jend)+ &
1013	& grad(i ,Jend))).gt.0.0_r8) THEN
1014	dUdx=grad(i-1,Jend)
1015	ELSE
1016	dUdx=grad(i ,Jend)
1017	END IF
1018	cff=MAX(dUdxdUdx+dUdedUde,eps)
1019	#ifdef RADIATION_2D
1020	Cx=MIN(cff,MAX(dUdt*dUdx,-cff))
1021	#else
1022	Cx=0.0_r8
1023	#endif
1024	Ce=dUdt*dUde
1025	#if defined CELERITY_WRITE && defined FORWARD_WRITE
1026	BOUNDARY(ng)%ubar_north_Cx(i)=Cx
1027	BOUNDARY(ng)%ubar_north_Ce(i)=Ce
1028	BOUNDARY(ng)%ubar_north_C2(i)=cff
1029	#endif
1030	ubar(i,Jend+1,kout)=(cff*ubar(i,Jend+1,know)+ &
1031	& Ce *ubar(i,Jend ,kout)- &
1032	& MAX(Cx,0.0_r8)*grad(i-1,Jend+1)- &
1033	& MIN(Cx,0.0_r8)*grad(i ,Jend+1))/ &
1034	& (cff+Ce)
1035
1036	IF (LBC(inorth,isUbar,ng)%nudging) THEN
1037	#ifdef IMPLICIT_NUDGING
1038	phi=dt(ng)/(tau+dt(ng))
1039	ubar(i,Jend+1,kout)=(1.0_r8-phi)*ubar(i,Jend+1,kout)+ &
1040	& phi*BOUNDARY(ng)%ubar_north(i)
1041	#else
1042	ubar(i,Jend+1,kout)=ubar(i,Jend+1,kout)+ &
1043	& tau*(BOUNDARY(ng)%ubar_north(i)- &
1044	& ubar(i,Jend+1,know))
1045	#endif
1046	END IF
1047	#ifdef MASKING
1048	ubar(i,Jend+1,kout)=ubar(i,Jend+1,kout)* &
1049	& GRID(ng)%umask(i,Jend+1)
1050	#endif
1051	END IF
1052	END DO
1053	!
1054	! Northern edge, Chapman boundary condition.
1055	!
1056	ELSE IF (LBC(inorth,isUbar,ng)%Flather.or. &
1057	& LBC(inorth,isUbar,ng)%reduced.or. &
1058	& LBC(inorth,isUbar,ng)%Shchepetkin) THEN
1059	DO i=IstrU,Iend
1060	IF (LBC_apply(ng)%north(i)) THEN
1061	cff=dt2d0.5_r8(GRID(ng)%pn(i-1,Jend)+ &
1062	& GRID(ng)%pn(i ,Jend))
1063	cff1=SQRT(g0.5_r8(GRID(ng)%h(i-1,Jend)+ &
1064	& zeta(i-1,Jend,know)+ &
1065	& GRID(ng)%h(i ,Jend)+ &
1066	& zeta(i ,Jend,know)))
1067	Ce=cff*cff1
1068	cff2=1.0_r8/(1.0_r8+Ce)
1069	ubar(i,Jend+1,kout)=cff2*(ubar(i,Jend+1,know)+ &
1070	& Ce*ubar(i,Jend,kout))
1071	#ifdef MASKING
1072	ubar(i,Jend+1,kout)=ubar(i,Jend+1,kout)* &
1073	& GRID(ng)%umask(i,Jend+1)
1074	#endif
1075	END IF
1076	END DO
1077	!
1078	! Northern edge, clamped boundary condition.
1079	!
1080	ELSE IF (LBC(inorth,isUbar,ng)%clamped) THEN
1081	DO i=IstrU,Iend
1082	IF (LBC_apply(ng)%north(i)) THEN
1083	ubar(i,Jend+1,kout)=BOUNDARY(ng)%ubar_north(i)
1084	#ifdef MASKING
1085	ubar(i,Jend+1,kout)=ubar(i,Jend+1,kout)* &
1086	& GRID(ng)%umask(i,Jend+1)
1087	#endif
1088	END IF
1089	END DO
1090	!
1091	! Northern edge, gradient boundary condition.
1092	!
1093	ELSE IF (LBC(inorth,isUbar,ng)%gradient) THEN
1094	DO i=IstrU,Iend
1095	IF (LBC_apply(ng)%north(i)) THEN
1096	ubar(i,Jend+1,kout)=ubar(i,Jend,kout)
1097	#ifdef MASKING
1098	ubar(i,Jend+1,kout)=ubar(i,Jend+1,kout)* &
1099	& GRID(ng)%umask(i,Jend+1)
1100	#endif
1101	END IF
1102	END DO
1103	!
1104	! Northern edge, closed boundary condition: free slip (gamma2=1) or
1105	! no slip (gamma2=-1).
1106	!
1107	ELSE IF (LBC(inorth,isUbar,ng)%closed) THEN
1108	IF (EWperiodic(ng)) THEN
1109	Imin=IstrU
1110	Imax=Iend
1111	ELSE
1112	Imin=Istr
1113	Imax=IendR
1114	END IF
1115	DO i=Imin,Imax
1116	IF (LBC_apply(ng)%north(i)) THEN
1117	ubar(i,Jend+1,kout)=gamma2(ng)*ubar(i,Jend,kout)
1118	#ifdef MASKING
1119	ubar(i,Jend+1,kout)=ubar(i,Jend+1,kout)* &
1120	& GRID(ng)%umask(i,Jend+1)
1121	#endif
1122	END IF
1123	END DO
1124	END IF
1125	END IF
1126	!
1127	!-----------------------------------------------------------------------
1128	! Boundary corners.
1129	!-----------------------------------------------------------------------
1130	!
1131	IF (.not.(EWperiodic(ng).or.NSperiodic(ng))) THEN
1132	IF (DOMAIN(ng)%SouthWest_Corner(tile)) THEN
1133	IF (LBC_apply(ng)%south(Istr ).and. &
1134	& LBC_apply(ng)%west (Jstr-1)) THEN
1135	ubar(Istr,Jstr-1,kout)=0.5_r8*(ubar(Istr+1,Jstr-1,kout)+ &
1136	& ubar(Istr ,Jstr ,kout))
1137	END IF
1138	END IF
1139	IF (DOMAIN(ng)%SouthEast_Corner(tile)) THEN
1140	IF (LBC_apply(ng)%south(Iend+1).and. &
1141	& LBC_apply(ng)%east (Jstr-1)) THEN
1142	ubar(Iend+1,Jstr-1,kout)=0.5_r8*(ubar(Iend ,Jstr-1,kout)+ &
1143	& ubar(Iend+1,Jstr ,kout))
1144	END IF
1145	END IF
1146	IF (DOMAIN(ng)%NorthWest_Corner(tile)) THEN
1147	IF (LBC_apply(ng)%north(Istr ).and. &
1148	& LBC_apply(ng)%west (Jend+1)) THEN
1149	ubar(Istr,Jend+1,kout)=0.5_r8*(ubar(Istr ,Jend ,kout)+ &
1150	& ubar(Istr+1,Jend+1,kout))
1151	END IF
1152	END IF
1153	IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN
1154	IF (LBC_apply(ng)%north(Iend+1).and. &
1155	& LBC_apply(ng)%east (Jend+1)) THEN
1156	ubar(Iend+1,Jend+1,kout)=0.5_r8*(ubar(Iend+1,Jend ,kout)+ &
1157	& ubar(Iend ,Jend+1,kout))
1158	END IF
1159	END IF
1160	END IF
1161
1162	#if defined WET_DRY
1163	!
1164	!-----------------------------------------------------------------------
1165	! Impose wetting and drying conditions.
1166	!-----------------------------------------------------------------------
1167	!
1168	IF (.not.EWperiodic(ng)) THEN
1169	IF (DOMAIN(ng)%Western_Edge(tile)) THEN
1170	DO j=Jstr,Jend
1171	IF (LBC_apply(ng)%west(j)) THEN
1172	cff1=ABS(ABS(GRID(ng)%umask_wet(Istr,j))-1.0_r8)
1173	cff2=0.5_r8+DSIGN(0.5_r8,ubar(Istr,j,kout))* &
1174	& GRID(ng)%umask_wet(Istr,j)
1175	cff=0.5_r8GRID(ng)%umask_wet(Istr,j)cff1+ &
1176	& cff2*(1.0_r8-cff1)
1177	ubar(Istr,j,kout)=ubar(Istr,j,kout)*cff
1178	END IF
1179	END DO
1180	END IF
1181	IF (DOMAIN(ng)%Eastern_Edge(tile)) THEN
1182	DO j=Jstr,Jend
1183	IF (LBC_apply(ng)%east(j)) THEN
1184	cff1=ABS(ABS(GRID(ng)%umask_wet(Iend+1,j))-1.0_r8)
1185	cff2=0.5_r8+DSIGN(0.5_r8,ubar(Iend+1,j,kout))* &
1186	& GRID(ng)%umask_wet(Iend+1,j)
1187	cff=0.5_r8GRID(ng)%umask_wet(Iend+1,j)cff1+ &
1188	& cff2*(1.0_r8-cff1)
1189	ubar(Iend+1,j,kout)=ubar(Iend+1,j,kout)*cff
1190	END IF
1191	END DO
1192	END IF
1193	END IF
1194
1195	IF (.not.NSperiodic(ng)) THEN
1196	IF (DOMAIN(ng)%Southern_Edge(tile)) THEN
1197	DO i=IstrU,Iend
1198	IF (LBC_apply(ng)%south(i)) THEN
1199	cff1=ABS(ABS(GRID(ng)%umask_wet(i,Jstr-1))-1.0_r8)
1200	cff2=0.5_r8+DSIGN(0.5_r8,ubar(i,Jstr-1,kout))* &
1201	& GRID(ng)%umask_wet(i,Jstr-1)
1202	cff=0.5_r8GRID(ng)%umask_wet(i,Jstr-1)cff1+ &
1203	& cff2*(1.0_r8-cff1)
1204	ubar(i,Jstr-1,kout)=ubar(i,Jstr-1,kout)*cff
1205	END IF
1206	END DO
1207	END IF
1208	IF (DOMAIN(ng)%Northern_Edge(tile)) THEN
1209	DO i=Istr,Iend
1210	IF (LBC_apply(ng)%north(i)) THEN
1211	cff1=ABS(ABS(GRID(ng)%umask_wet(i,Jend+1))-1.0_r8)
1212	cff2=0.5_r8+DSIGN(0.5_r8,ubar(i,Jend+1,kout))* &
1213	& GRID(ng)%umask_wet(i,Jend+1)
1214	cff=0.5_r8GRID(ng)%umask_wet(i,Jend+1)cff1+ &
1215	& cff2*(1.0_r8-cff1)
1216	ubar(i,Jend+1,kout)=ubar(i,Jend+1,kout)*cff
1217	END IF
1218	END DO
1219	END IF
1220	END IF
1221
1222	IF (.not.(EWperiodic(ng).or.NSperiodic(ng))) THEN
1223	IF (DOMAIN(ng)%SouthWest_Corner(tile)) THEN
1224	IF (LBC_apply(ng)%south(Istr ).and. &
1225	& LBC_apply(ng)%west (Jstr-1)) THEN
1226	cff1=ABS(ABS(GRID(ng)%umask_wet(Istr,Jstr-1))-1.0_r8)
1227	cff2=0.5_r8+DSIGN(0.5_r8,ubar(Istr,Jstr-1,kout))* &
1228	& GRID(ng)%umask_wet(Istr,Jstr-1)
1229	cff=0.5_r8GRID(ng)%umask_wet(Istr,Jstr-1)cff1+ &
1230	& cff2*(1.0_r8-cff1)
1231	ubar(Istr,Jstr-1,kout)=ubar(Istr,Jstr-1,kout)*cff
1232	END IF
1233	END IF
1234	IF (DOMAIN(ng)%SouthEast_Corner(tile)) THEN
1235	IF (LBC_apply(ng)%south(Iend+1).and. &
1236	& LBC_apply(ng)%east (Jstr-1)) THEN
1237	cff1=ABS(ABS(GRID(ng)%umask_wet(Iend+1,Jstr-1))-1.0_r8)
1238	cff2=0.5_r8+DSIGN(0.5_r8,ubar(Iend+1,Jstr-1,kout))* &
1239	& GRID(ng)%umask_wet(Iend+1,Jstr-1)
1240	cff=0.5_r8GRID(ng)%umask_wet(Iend+1,Jstr-1)cff1+ &
1241	& cff2*(1.0_r8-cff1)
1242	ubar(Iend+1,Jstr-1,kout)=ubar(Iend+1,Jstr-1,kout)*cff
1243	END IF
1244	END IF
1245	IF (DOMAIN(ng)%NorthWest_Corner(tile)) THEN
1246	IF (LBC_apply(ng)%north(Istr ).and. &
1247	& LBC_apply(ng)%west (Jend+1)) THEN
1248	cff1=ABS(ABS(GRID(ng)%umask_wet(Istr,Jend+1))-1.0_r8)
1249	cff2=0.5_r8+DSIGN(0.5_r8,ubar(Istr,Jend+1,kout))* &
1250	& GRID(ng)%umask_wet(Istr,Jend+1)
1251	cff=0.5_r8GRID(ng)%umask_wet(Istr,Jend+1)cff1+ &
1252	& cff2*(1.0_r8-cff1)
1253	ubar(Istr,Jend+1,kout)=ubar(Istr,Jend+1,kout)*cff
1254	END IF
1255	END IF
1256	IF (DOMAIN(ng)%NorthEast_Corner(tile)) THEN
1257	IF (LBC_apply(ng)%north(Iend+1).and. &
1258	& LBC_apply(ng)%east (Jend+1)) THEN
1259	cff1=ABS(ABS(GRID(ng)%umask_wet(Iend+1,Jend+1))-1.0_r8)
1260	cff2=0.5_r8+DSIGN(0.5_r8,ubar(Iend+1,Jend+1,kout))* &
1261	& GRID(ng)%umask_wet(Iend+1,Jend+1)
1262	cff=0.5_r8*GRID(ng)%umask_wet(Iend+1,Jend+1)+cff1+ &
1263	& cff2*(1.0_r8-cff1)
1264	ubar(Iend+1,Jend+1,kout)=ubar(Iend+1,Jend+1,kout)*cff
1265	END IF
1266	END IF
1267	END IF
1268	#endif
1269
1270	RETURN
1271	END SUBROUTINE u2dbc_tile
1272	END MODULE u2dbc_mod

Download in other formats:

Original Format