New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

agrif_opa_sponge.F90 in branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC – NEMO

Context Navigation

source: branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 8835

Last change on this file since 8835 was 8835, checked in by timgraham, 6 years ago
Finished sponge layer for tsn, U and V
Property svn:keywords set to `Id`
File size: 22.9 KB

Line
1	#undef SPONGE && undef SPONGE_TOP
2
3	MODULE agrif_opa_sponge
4	!!======================================================================
5	!! * MODULE agrif_opa_update *
6	!! AGRIF :
7	!!======================================================================
8	!! History :
9	!!----------------------------------------------------------------------
10	#if defined key_agrif
11	USE par_oce
12	USE oce
13	USE dom_oce
14	USE in_out_manager
15	USE agrif_oce
16	USE wrk_nemo
17	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
18
19	IMPLICIT NONE
20	PRIVATE
21
22	PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn
23	PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge
24
25	!!----------------------------------------------------------------------
26	!! NEMO/NST 3.7 , NEMO Consortium (2015)
27	!! $Id$
28	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
29	!!----------------------------------------------------------------------
30	CONTAINS
31
32	SUBROUTINE Agrif_Sponge_Tra
33	!!---------------------------------------------
34	!! * ROUTINE Agrif_Sponge_Tra *
35	!!---------------------------------------------
36	REAL(wp) :: timecoeff
37	!!---------------------------------------------
38	!
39	#if defined SPONGE
40	timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
41
42	CALL Agrif_Sponge
43	Agrif_SpecialValue=0.
44	Agrif_UseSpecialValue = .TRUE.
45	tabspongedone_tsn = .FALSE.
46
47	CALL Agrif_Bc_Variable(tsn_sponge_id,calledweight=timecoeff,procname=interptsn_sponge)
48
49	Agrif_UseSpecialValue = .FALSE.
50	#endif
51	!
52	END SUBROUTINE Agrif_Sponge_Tra
53
54
55	SUBROUTINE Agrif_Sponge_dyn
56	!!---------------------------------------------
57	!! * ROUTINE Agrif_Sponge_dyn *
58	!!---------------------------------------------
59	REAL(wp) :: timecoeff
60	!!---------------------------------------------
61	#if defined SPONGE
62	timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
63
64	Agrif_SpecialValue=0.
65	Agrif_UseSpecialValue = ln_spc_dyn
66
67	tabspongedone_u = .FALSE.
68	tabspongedone_v = .FALSE.
69	CALL Agrif_Bc_Variable(un_sponge_id,calledweight=timecoeff,procname=interpun_sponge)
70
71	tabspongedone_u = .FALSE.
72	tabspongedone_v = .FALSE.
73	CALL Agrif_Bc_Variable(vn_sponge_id,calledweight=timecoeff,procname=interpvn_sponge)
74
75	Agrif_UseSpecialValue = .FALSE.
76	#endif
77	!
78	END SUBROUTINE Agrif_Sponge_dyn
79
80
81	SUBROUTINE Agrif_Sponge
82	!!---------------------------------------------
83	!! * ROUTINE Agrif_Sponge *
84	!!---------------------------------------------
85	INTEGER :: ji,jj,jk
86	INTEGER :: ispongearea, ilci, ilcj
87	LOGICAL :: ll_spdone
88	REAL(wp) :: z1spongearea, zramp
89	REAL(wp), POINTER, DIMENSION(:,:) :: ztabramp
90
91	#if defined SPONGE \|\| defined SPONGE_TOP
92	ll_spdone=.TRUE.
93	IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
94	! Define ramp from boundaries towards domain interior
95	! at T-points
96	! Store it in ztabramp
97	ll_spdone=.FALSE.
98
99	CALL wrk_alloc( jpi, jpj, ztabramp )
100
101	ispongearea = 2 + nn_sponge_len * Agrif_irhox()
102	ilci = nlci - ispongearea
103	ilcj = nlcj - ispongearea
104	z1spongearea = 1._wp / REAL( ispongearea - 2 )
105
106	ztabramp(:,:) = 0._wp
107
108	IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
109	DO jj = 1, jpj
110	IF ( umask(2,jj,1) == 1._wp ) THEN
111	DO ji = 2, ispongearea
112	ztabramp(ji,jj) = ( ispongearea-ji ) * z1spongearea
113	END DO
114	ENDIF
115	ENDDO
116	ENDIF
117
118	IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
119	DO jj = 1, jpj
120	IF ( umask(nlci-2,jj,1) == 1._wp ) THEN
121	DO ji = ilci+1,nlci-1
122	zramp = (ji - (ilci+1) ) * z1spongearea
123	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
124	ENDDO
125	ENDIF
126	ENDDO
127	ENDIF
128
129	IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
130	DO ji = 1, jpi
131	IF ( vmask(ji,2,1) == 1._wp ) THEN
132	DO jj = 2, ispongearea
133	zramp = ( ispongearea-jj ) * z1spongearea
134	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
135	END DO
136	ENDIF
137	ENDDO
138	ENDIF
139
140	IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
141	DO ji = 1, jpi
142	IF ( vmask(ji,nlcj-2,1) == 1._wp ) THEN
143	DO jj = ilcj+1,nlcj-1
144	zramp = (jj - (ilcj+1) ) * z1spongearea
145	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
146	END DO
147	ENDIF
148	ENDDO
149	ENDIF
150
151	ENDIF
152
153	! Tracers
154	IF( .NOT. spongedoneT ) THEN
155	fsaht_spu(:,:) = 0._wp
156	fsaht_spv(:,:) = 0._wp
157	DO jj = 2, jpjm1
158	DO ji = 2, jpim1 ! vector opt.
159	fsaht_spu(ji,jj) = 0.5_wp * visc_tra * (ztabramp(ji,jj) + ztabramp(ji+1,jj ))
160	fsaht_spv(ji,jj) = 0.5_wp * visc_tra * (ztabramp(ji,jj) + ztabramp(ji ,jj+1))
161	END DO
162	END DO
163
164	CALL lbc_lnk( fsaht_spu, 'U', 1. ) ! Lateral boundary conditions
165	CALL lbc_lnk( fsaht_spv, 'V', 1. )
166	spongedoneT = .TRUE.
167	ENDIF
168
169	! Dynamics
170	IF( .NOT. spongedoneU ) THEN
171	fsahm_spt(:,:) = 0._wp
172	fsahm_spf(:,:) = 0._wp
173	DO jj = 2, jpjm1
174	DO ji = 2, jpim1 ! vector opt.
175	fsahm_spt(ji,jj) = visc_dyn * ztabramp(ji,jj)
176	fsahm_spf(ji,jj) = 0.25_wp * visc_dyn * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1) &
177	& +ztabramp(ji,jj) + ztabramp(ji+1,jj ) )
178	END DO
179	END DO
180
181	CALL lbc_lnk( fsahm_spt, 'T', 1. ) ! Lateral boundary conditions
182	CALL lbc_lnk( fsahm_spf, 'F', 1. )
183	spongedoneU = .TRUE.
184	ENDIF
185	!
186	IF (.NOT.ll_spdone) CALL wrk_dealloc( jpi, jpj, ztabramp )
187	!
188	#endif
189	!
190	END SUBROUTINE Agrif_Sponge
191
192
193	SUBROUTINE interptsn_sponge(tabres,i1,i2,j1,j2,k1,k2,n1,n2,before)
194	!!---------------------------------------------
195	!! * ROUTINE interptsn_sponge *
196	!!---------------------------------------------
197	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
198	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
199	LOGICAL, INTENT(in) :: before
200	!
201	INTEGER :: ji, jj, jk, jn ! dummy loop indices
202	INTEGER :: iku, ikv
203	REAL(wp) :: ztsa, zabe1, zabe2, zbtr
204	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv
205	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
206	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
207	REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
208	REAL(wp) :: h_in(k1:k2)
209	REAL(wp) :: h_out(1:jpk)
210	INTEGER :: N_in, N_out
211	REAL(wp) :: h_diff, zrhoxy
212	!
213	IF( before ) THEN
214	DO jn = n1,n2-1
215	DO jk=k1,k2
216	DO jj=j1,j2
217	DO ji=i1,i2
218	tabres(ji,jj,jk,jn) = tsb(ji,jj,jk,jn) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)
219	END DO
220	END DO
221	END DO
222	END DO
223	DO jk=k1,k2
224	DO jj=j1,j2
225	DO ji=i1,i2
226	tabres(ji,jj,jk,n2) = tmask(ji,jj,jk) * e1e2t(ji,jj) * e3t_n(ji,jj,jk)
227	END DO
228	END DO
229	END DO
230	ELSE
231	!
232	#ifdef key_vertical
233	tabres_child(:,:,:,:) = 0.
234	zrhoxy = Agrif_rhox()*Agrif_rhoy()
235	DO jj=j1,j2
236	DO ji=i1,i2
237	N_in = 0
238	DO jk=k1,k2 !k2 = jpk of parent grid
239	IF (tabres(ji,jj,jk,n2) == 0) EXIT
240	N_in = N_in + 1
241	tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
242	h_in(N_in) = tabres(ji,jj,jk,n2)/(e1e2t(ji,jj)*zrhoxy)
243	END DO
244	N_out = 0
245	DO jk=1,jpk ! jpk of child grid
246	IF (tmask(ji,jj,jk) == 0) EXIT ! TODO: Will not work with ISF. !This doesn't seem to work at the moment in GYRE but is OK in overflow model
247	N_out = N_out + 1
248	h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above
249	ENDDO
250	IF (N_in > 0) THEN
251	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
252	tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for?????
253	DO jn=1,jpts
254	call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
255	ENDDO
256	ENDIF
257	DO jk=1,jpk
258	tsbdiff(ji,jj,jk,n1:n2-1) = tsb(ji, jj,jk,n1:n2-1) - tabres_child(ji,jj,jk,n1:n2-1)
259	ENDDO
260	ENDDO
261	ENDDO
262	#else
263	do jk=k1,k2
264	do jj=j1,j2
265	do ji=i1,i2
266	! This will be slow - Is there a way to speed it up and avoid divide by zero?
267	IF (tabres(ji,jj,jk,n2) .NE. 0) THEN
268	tsbdiff(ji,jj,jk,n1:n2-1) = tsb(ji,jj,jk,n1:n2) - tabres(ji,jj,jk,n1:n2-1)/tabres(ji,jj,jk,n2)
269	ELSE
270	tsbdiff(ji,jj,jk,n1:n2-1) = 0._wp
271	ENDIF
272	enddo
273	enddo
274	enddo
275	#endif
276	DO jn = 1, jpts
277	DO jk = 1, jpkm1
278	DO jj = j1,j2-1
279	DO ji = i1,i2-1
280	zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
281	zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
282	ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) )
283	ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
284	END DO
285	END DO
286	!
287	IF( ln_zps ) THEN ! set gradient at partial step level
288	DO jj = j1,j2-1
289	DO ji = i1,i2-1
290	! last level
291	iku = mbku(ji,jj)
292	ikv = mbkv(ji,jj)
293	IF( iku == jk ) ztu(ji,jj,jk) = 0._wp
294	IF( ikv == jk ) ztv(ji,jj,jk) = 0._wp
295	END DO
296	END DO
297	ENDIF
298	END DO
299	!
300	DO jk = 1, jpkm1
301	DO jj = j1+1,j2-1
302	DO ji = i1+1,i2-1
303	IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
304	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
305	! horizontal diffusive trends
306	ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) )
307	! add it to the general tracer trends
308	tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa
309	ENDIF
310	END DO
311	END DO
312	END DO
313	!
314	END DO
315	!
316	tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
317	!
318	ENDIF
319	!
320	END SUBROUTINE interptsn_sponge
321
322
323	SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
324	!!---------------------------------------------
325	!! * ROUTINE interpun_sponge *
326	!!---------------------------------------------
327	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
328	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
329	LOGICAL, INTENT(in) :: before
330
331	INTEGER :: ji,jj,jk,N_in,N_out
332
333	! sponge parameters
334	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff, zrhoy
335	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
336	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
337	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
338	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
339	REAL(wp), DIMENSION(1:jpk) :: h_out
340	INTEGER :: jmax
341	!!---------------------------------------------
342	!
343	IF( before ) THEN
344	#ifdef key_vertical
345	DO jk=1,jpk
346	DO jj=j1,j2
347	DO ji=i1,i2
348	tabres(ji,jj,jk,m1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk))
349	tabres(ji,jj,jk,m2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk))
350	END DO
351	END DO
352	END DO
353
354	#else
355	tabres(i1:i2,j1:j2,:,m1) = un(i1:i2,j1:j2,:)
356	#endif
357	ELSE
358	#ifdef key_vertical
359	tabres_child(:,:,:) = 0._wp
360	zrhoy = Agrif_rhoy()
361	DO jj=j1,j2
362	DO ji=i1,i2
363	N_in = 0
364	DO jk=k1,k2
365	IF (tabres(ji,jj,jk,m2) == 0) EXIT
366	N_in = N_in + 1
367	tabin(jk) = tabres(ji,jj,jk,m1)/tabres(ji,jj,jk,m2)
368	h_in(N_in) = tabres(ji,jj,jk,m2)/(e2u(ji,jj)*zrhoy)
369	ENDDO
370	!
371	IF (N_in == 0) THEN
372	tabres_child(ji,jj,:) = 0.
373	CYCLE
374	ENDIF
375
376	N_out = 0
377	DO jk=1,jpk
378	if (umask(ji,jj,jk) == 0) EXIT
379	N_out = N_out + 1
380	h_out(N_out) = e3u_n(ji,jj,jk)
381	ENDDO
382
383	IF (N_out == 0) THEN
384	tabres_child(ji,jj,:) = 0.
385	CYCLE
386	ENDIF
387
388	IF (N_in * N_out > 0) THEN
389	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
390	if (h_diff < -1.e4) then
391	print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
392	endif
393	ENDIF
394	call reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
395
396	ENDDO
397	ENDDO
398	DO jk = 1, jpkm1
399	DO jj=j1,j2
400	ubdiff(i1:i2,jj,jk) = (ub(i1:i2,jj,jk) - tabres_child(i1:i2,jj,jk))*umask(i1:i2,jj,jk)
401	END DO
402	END DO
403	#else
404	ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:))*umask(i1:i2,j1:j2,:)
405	#endif
406	!
407	DO jk = 1, jpkm1 ! Horizontal slab
408	! ! ===============
409
410	! ! --------
411	! Horizontal divergence ! div
412	! ! --------
413	DO jj = j1,j2
414	DO ji = i1+1,i2 ! vector opt.
415	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
416	hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)e3u_n(ji ,jj,jk) ubdiff(ji ,jj,jk) &
417	& -e2u(ji-1,jj)e3u_n(ji-1,jj,jk) ubdiff(ji-1,jj,jk) ) * zbtr
418	END DO
419	END DO
420
421	DO jj = j1,j2-1
422	DO ji = i1,i2 ! vector opt.
423	zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
424	rotdiff(ji,jj,jk) = (-e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &
425	+e1u(ji,jj ) * ubdiff(ji,jj ,jk) &
426	& ) * fmask(ji,jj,jk) * zbtr
427	END DO
428	END DO
429	END DO
430	!
431	DO jj = j1+1, j2-1
432	DO ji = i1+1, i2-1 ! vector opt.
433
434	IF (.NOT. tabspongedone_u(ji,jj)) THEN
435	DO jk = 1, jpkm1 ! Horizontal slab
436	ze2u = rotdiff (ji,jj,jk)
437	ze1v = hdivdiff(ji,jj,jk)
438	! horizontal diffusive trends
439	zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) ) &
440	+ ( hdivdiff(ji+1,jj,jk) - ze1v ) / e1u(ji,jj)
441
442	! add it to the general momentum trends
443	ua(ji,jj,jk) = ua(ji,jj,jk) + zua
444
445	END DO
446	ENDIF
447
448	END DO
449	END DO
450
451	tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
452
453	jmax = j2-1
454	IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-3)
455
456	DO jj = j1+1, jmax
457	DO ji = i1+1, i2 ! vector opt.
458
459	IF (.NOT. tabspongedone_v(ji,jj)) THEN
460	DO jk = 1, jpkm1 ! Horizontal slab
461	ze2u = rotdiff (ji,jj,jk)
462	ze1v = hdivdiff(ji,jj,jk)
463
464	! horizontal diffusive trends
465	zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) ) &
466	+ ( hdivdiff(ji,jj+1,jk) - ze1v ) / e2v(ji,jj)
467
468	! add it to the general momentum trends
469	va(ji,jj,jk) = va(ji,jj,jk) + zva
470	END DO
471	ENDIF
472	!
473	END DO
474	END DO
475	!
476	tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
477	!
478	ENDIF
479	!
480	END SUBROUTINE interpun_sponge
481
482
483	SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir)
484	!!---------------------------------------------
485	!! * ROUTINE interpvn_sponge *
486	!!---------------------------------------------
487	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
488	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
489	LOGICAL, INTENT(in) :: before
490	INTEGER, INTENT(in) :: nb , ndir
491	!
492	INTEGER :: ji, jj, jk, N_in, N_out
493	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff, zrhox
494	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
495	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
496	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
497	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
498	REAL(wp), DIMENSION(1:jpk) :: h_out
499	INTEGER :: imax
500	!!---------------------------------------------
501
502	IF( before ) THEN
503	#ifdef key_vertical
504	DO jk=1,jpk
505	DO jj=j1,j2
506	DO ji=i1,i2
507	tabres(ji,jj,jk,m1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk))
508	tabres(ji,jj,jk,m2) = (vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk))
509	END DO
510	END DO
511	END DO
512	#else
513	tabres = vn(i1:i2,j1:j2,:)
514	#endif
515	ELSE
516	#ifdef key_vertical
517	zrhox = Agrif_rhox()
518	tabres_child(:,:,:) = 0._wp
519	DO jj=j1,j2
520	DO ji=i1,i2
521	N_in = 0
522	DO jk=k1,k2
523	IF (tabres(ji,jj,jk,m2) == 0) EXIT
524	N_in = N_in + 1
525	tabin(jk) = tabres(ji,jj,jk,m1)/tabres(ji,jj,jk,m2)
526	h_in(N_in) = tabres(ji,jj,jk,m2)/(e1v(ji,jj)*zrhox)
527	ENDDO
528
529	IF (N_in == 0) THEN
530	tabres_child(ji,jj,:) = 0.
531	CYCLE
532	ENDIF
533
534	N_out = 0
535	DO jk=1,jpk
536	if (umask(ji,jj,jk) == 0) EXIT
537	N_out = N_out + 1
538	h_out(N_out) = e3v_n(ji,jj,jk)
539	ENDDO
540
541	IF (N_in * N_out > 0) THEN
542	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
543	if (h_diff < -1.e4) then
544	print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
545	endif
546	ENDIF
547	call reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out)
548	ENDDO
549	ENDDO
550	DO jk = 1, jpkm1
551	DO jj=j1,j2
552	vbdiff(i1:i2,jj,jk) = (vb(i1:i2,jj,jk) - tabres_child(i1:i2,jj,jk))*vmask(i1:i2,jj,jk)
553	END DO
554	END DO
555	#else
556	!
557	vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(:,:,:))*vmask(i1:i2,j1:j2,:)
558	#endif
559	!
560	DO jk = 1, jpkm1 ! Horizontal slab
561	! ! ===============
562
563	! ! --------
564	! Horizontal divergence ! div
565	! ! --------
566	DO jj = j1+1,j2
567	DO ji = i1,i2 ! vector opt.
568	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
569	hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v_n(ji,jj ,jk) * vbdiff(ji,jj ,jk) &
570	& -e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk) ) * zbtr
571	END DO
572	END DO
573	DO jj = j1,j2
574	DO ji = i1,i2-1 ! vector opt.
575	zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
576	rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) &
577	& -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr
578	END DO
579	END DO
580	END DO
581
582	! ! ===============
583	!
584
585	imax = i2-1
586	IF ((nbondi == 1).OR.(nbondi == 2)) imax = MIN(imax,nlci-3)
587
588	DO jj = j1+1, j2
589	DO ji = i1+1, imax ! vector opt.
590	IF( .NOT. tabspongedone_u(ji,jj) ) THEN
591	DO jk = 1, jpkm1
592	ua(ji,jj,jk) = ua(ji,jj,jk) &
593	& - ( rotdiff (ji ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) ) &
594	& + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj ,jk)) * r1_e1u(ji,jj)
595	END DO
596	ENDIF
597	END DO
598	END DO
599	!
600	tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
601	!
602	DO jj = j1+1, j2-1
603	DO ji = i1+1, i2-1 ! vector opt.
604	IF( .NOT. tabspongedone_v(ji,jj) ) THEN
605	DO jk = 1, jpkm1
606	va(ji,jj,jk) = va(ji,jj,jk) &
607	& + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) ) &
608	& + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj)
609	END DO
610	ENDIF
611	END DO
612	END DO
613	tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
614	ENDIF
615	!
616	END SUBROUTINE interpvn_sponge
617
618	#else
619	CONTAINS
620	SUBROUTINE agrif_opa_sponge_empty
621	!!---------------------------------------------
622	!! * ROUTINE agrif_OPA_sponge_empty *
623	!!---------------------------------------------
624	WRITE(,) 'agrif_opa_sponge : You should not have seen this print! error?'
625	END SUBROUTINE agrif_opa_sponge_empty
626	#endif
627
628	!!======================================================================
629	END MODULE agrif_opa_sponge

Note: See TracBrowser for help on using the repository browser.

Download in other formats: