Context Navigation

agrif_opa_sponge.F90 @ 8902

Last change on this file since 8902 was 8902, checked in by timgraham, 6 years ago
Bug fixes for non key_vertical case, reverted to interpolating tracer rather than tracer content
Property svn:keywords set to `Id`
File size: 22.5 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) = tsn(ji,jj,jk,jn)
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)
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
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	tsbdiff(ji,jj,jk,n1:n2-1) = tsb(ji,jj,jk,n1:n2) - tmask(ji,jj,jk)*tabres(ji,jj,jk,n1:n2-1)
267	enddo
268	enddo
269	enddo
270	#endif
271	DO jn = 1, jpts
272	DO jk = 1, jpkm1
273	DO jj = j1,j2-1
274	DO ji = i1,i2-1
275	zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
276	zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
277	ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) )
278	ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
279	END DO
280	END DO
281	!
282	IF( ln_zps ) THEN ! set gradient at partial step level
283	DO jj = j1,j2-1
284	DO ji = i1,i2-1
285	! last level
286	iku = mbku(ji,jj)
287	ikv = mbkv(ji,jj)
288	IF( iku == jk ) ztu(ji,jj,jk) = 0._wp
289	IF( ikv == jk ) ztv(ji,jj,jk) = 0._wp
290	END DO
291	END DO
292	ENDIF
293	END DO
294	!
295	DO jk = 1, jpkm1
296	DO jj = j1+1,j2-1
297	DO ji = i1+1,i2-1
298	IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
299	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
300	! horizontal diffusive trends
301	ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) )
302	! add it to the general tracer trends
303	tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa
304	ENDIF
305	END DO
306	END DO
307	END DO
308	!
309	END DO
310	!
311	tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
312	!
313	ENDIF
314	!
315	END SUBROUTINE interptsn_sponge
316
317
318	SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
319	!!---------------------------------------------
320	!! * ROUTINE interpun_sponge *
321	!!---------------------------------------------
322	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
323	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
324	LOGICAL, INTENT(in) :: before
325
326	INTEGER :: ji,jj,jk,N_in,N_out
327
328	! sponge parameters
329	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff, zrhoy
330	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
331	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
332	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
333	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
334	REAL(wp), DIMENSION(1:jpk) :: h_out
335	INTEGER :: jmax
336	!!---------------------------------------------
337	!
338	IF( before ) THEN
339	#ifdef key_vertical
340	DO jk=1,jpk
341	DO jj=j1,j2
342	DO ji=i1,i2
343	tabres(ji,jj,jk,m1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk))
344	tabres(ji,jj,jk,m2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk))
345	END DO
346	END DO
347	END DO
348
349	#else
350	tabres(i1:i2,j1:j2,:,m1) = un(i1:i2,j1:j2,:)
351	#endif
352	ELSE
353	#ifdef key_vertical
354	tabres_child(:,:,:) = 0._wp
355	zrhoy = Agrif_rhoy()
356	DO jj=j1,j2
357	DO ji=i1,i2
358	N_in = 0
359	DO jk=k1,k2
360	IF (tabres(ji,jj,jk,m2) == 0) EXIT
361	N_in = N_in + 1
362	tabin(jk) = tabres(ji,jj,jk,m1)/tabres(ji,jj,jk,m2)
363	h_in(N_in) = tabres(ji,jj,jk,m2)/(e2u(ji,jj)*zrhoy)
364	ENDDO
365	!
366	IF (N_in == 0) THEN
367	tabres_child(ji,jj,:) = 0.
368	CYCLE
369	ENDIF
370
371	N_out = 0
372	DO jk=1,jpk
373	if (umask(ji,jj,jk) == 0) EXIT
374	N_out = N_out + 1
375	h_out(N_out) = e3u_n(ji,jj,jk)
376	ENDDO
377
378	IF (N_out == 0) THEN
379	tabres_child(ji,jj,:) = 0.
380	CYCLE
381	ENDIF
382
383	IF (N_in * N_out > 0) THEN
384	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
385	if (h_diff < -1.e4) then
386	print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
387	endif
388	ENDIF
389	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)
390
391	ENDDO
392	ENDDO
393	DO jk = 1, jpkm1
394	DO jj=j1,j2
395	ubdiff(i1:i2,jj,jk) = (ub(i1:i2,jj,jk) - tabres_child(i1:i2,jj,jk))*umask(i1:i2,jj,jk)
396	END DO
397	END DO
398	#else
399	ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:,1))*umask(i1:i2,j1:j2,:)
400	#endif
401	!
402	DO jk = 1, jpkm1 ! Horizontal slab
403	! ! ===============
404
405	! ! --------
406	! Horizontal divergence ! div
407	! ! --------
408	DO jj = j1,j2
409	DO ji = i1+1,i2 ! vector opt.
410	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
411	hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)e3u_n(ji ,jj,jk) ubdiff(ji ,jj,jk) &
412	& -e2u(ji-1,jj)e3u_n(ji-1,jj,jk) ubdiff(ji-1,jj,jk) ) * zbtr
413	END DO
414	END DO
415
416	DO jj = j1,j2-1
417	DO ji = i1,i2 ! vector opt.
418	zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
419	rotdiff(ji,jj,jk) = (-e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &
420	+e1u(ji,jj ) * ubdiff(ji,jj ,jk) &
421	& ) * fmask(ji,jj,jk) * zbtr
422	END DO
423	END DO
424	END DO
425	!
426	DO jj = j1+1, j2-1
427	DO ji = i1+1, i2-1 ! vector opt.
428
429	IF (.NOT. tabspongedone_u(ji,jj)) THEN
430	DO jk = 1, jpkm1 ! Horizontal slab
431	ze2u = rotdiff (ji,jj,jk)
432	ze1v = hdivdiff(ji,jj,jk)
433	! horizontal diffusive trends
434	zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) ) &
435	+ ( hdivdiff(ji+1,jj,jk) - ze1v ) / e1u(ji,jj)
436
437	! add it to the general momentum trends
438	ua(ji,jj,jk) = ua(ji,jj,jk) + zua
439
440	END DO
441	ENDIF
442
443	END DO
444	END DO
445
446	tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
447
448	jmax = j2-1
449	IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-3)
450
451	DO jj = j1+1, jmax
452	DO ji = i1+1, i2 ! vector opt.
453
454	IF (.NOT. tabspongedone_v(ji,jj)) THEN
455	DO jk = 1, jpkm1 ! Horizontal slab
456	ze2u = rotdiff (ji,jj,jk)
457	ze1v = hdivdiff(ji,jj,jk)
458
459	! horizontal diffusive trends
460	zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) ) &
461	+ ( hdivdiff(ji,jj+1,jk) - ze1v ) / e2v(ji,jj)
462
463	! add it to the general momentum trends
464	va(ji,jj,jk) = va(ji,jj,jk) + zva
465	END DO
466	ENDIF
467	!
468	END DO
469	END DO
470	!
471	tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
472	!
473	ENDIF
474	!
475	END SUBROUTINE interpun_sponge
476
477
478	SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir)
479	!!---------------------------------------------
480	!! * ROUTINE interpvn_sponge *
481	!!---------------------------------------------
482	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
483	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
484	LOGICAL, INTENT(in) :: before
485	INTEGER, INTENT(in) :: nb , ndir
486	!
487	INTEGER :: ji, jj, jk, N_in, N_out
488	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff, zrhox
489	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
490	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
491	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
492	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
493	REAL(wp), DIMENSION(1:jpk) :: h_out
494	INTEGER :: imax
495	!!---------------------------------------------
496
497	IF( before ) THEN
498	#ifdef key_vertical
499	DO jk=1,jpk
500	DO jj=j1,j2
501	DO ji=i1,i2
502	tabres(ji,jj,jk,m1) = (e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)*vmask(ji,jj,jk))
503	tabres(ji,jj,jk,m2) = (vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk))
504	END DO
505	END DO
506	END DO
507	#else
508	tabres(i1:i2,j1:j2,:,m1) = vn(i1:i2,j1:j2,:)
509	#endif
510	ELSE
511	#ifdef key_vertical
512	zrhox = Agrif_rhox()
513	tabres_child(:,:,:) = 0._wp
514	DO jj=j1,j2
515	DO ji=i1,i2
516	N_in = 0
517	DO jk=k1,k2
518	IF (tabres(ji,jj,jk,m2) == 0) EXIT
519	N_in = N_in + 1
520	tabin(jk) = tabres(ji,jj,jk,m1)/tabres(ji,jj,jk,m2)
521	h_in(N_in) = tabres(ji,jj,jk,m2)/(e1v(ji,jj)*zrhox)
522	ENDDO
523
524	IF (N_in == 0) THEN
525	tabres_child(ji,jj,:) = 0.
526	CYCLE
527	ENDIF
528
529	N_out = 0
530	DO jk=1,jpk
531	if (umask(ji,jj,jk) == 0) EXIT
532	N_out = N_out + 1
533	h_out(N_out) = e3v_n(ji,jj,jk)
534	ENDDO
535
536	IF (N_in * N_out > 0) THEN
537	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
538	if (h_diff < -1.e4) then
539	print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
540	endif
541	ENDIF
542	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)
543	ENDDO
544	ENDDO
545	DO jk = 1, jpkm1
546	DO jj=j1,j2
547	vbdiff(i1:i2,jj,jk) = (vb(i1:i2,jj,jk) - tabres_child(i1:i2,jj,jk))*vmask(i1:i2,jj,jk)
548	END DO
549	END DO
550	#else
551	!
552	vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(:,:,:))*vmask(i1:i2,j1:j2,:)
553	#endif
554	!
555	DO jk = 1, jpkm1 ! Horizontal slab
556	! ! ===============
557
558	! ! --------
559	! Horizontal divergence ! div
560	! ! --------
561	DO jj = j1+1,j2
562	DO ji = i1,i2 ! vector opt.
563	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
564	hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v_n(ji,jj ,jk) * vbdiff(ji,jj ,jk) &
565	& -e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk) ) * zbtr
566	END DO
567	END DO
568	DO jj = j1,j2
569	DO ji = i1,i2-1 ! vector opt.
570	zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
571	rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) &
572	& -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr
573	END DO
574	END DO
575	END DO
576
577	! ! ===============
578	!
579
580	imax = i2-1
581	IF ((nbondi == 1).OR.(nbondi == 2)) imax = MIN(imax,nlci-3)
582
583	DO jj = j1+1, j2
584	DO ji = i1+1, imax ! vector opt.
585	IF( .NOT. tabspongedone_u(ji,jj) ) THEN
586	DO jk = 1, jpkm1
587	ua(ji,jj,jk) = ua(ji,jj,jk) &
588	& - ( rotdiff (ji ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) ) &
589	& + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj ,jk)) * r1_e1u(ji,jj)
590	END DO
591	ENDIF
592	END DO
593	END DO
594	!
595	tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
596	!
597	DO jj = j1+1, j2-1
598	DO ji = i1+1, i2-1 ! vector opt.
599	IF( .NOT. tabspongedone_v(ji,jj) ) THEN
600	DO jk = 1, jpkm1
601	va(ji,jj,jk) = va(ji,jj,jk) &
602	& + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) ) &
603	& + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj)
604	END DO
605	ENDIF
606	END DO
607	END DO
608	tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
609	ENDIF
610	!
611	END SUBROUTINE interpvn_sponge
612
613	#else
614	CONTAINS
615	SUBROUTINE agrif_opa_sponge_empty
616	!!---------------------------------------------
617	!! * ROUTINE agrif_OPA_sponge_empty *
618	!!---------------------------------------------
619	WRITE(,) 'agrif_opa_sponge : You should not have seen this print! error?'
620	END SUBROUTINE agrif_opa_sponge_empty
621	#endif
622
623	!!======================================================================
624	END MODULE agrif_opa_sponge

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

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

Context Navigation

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

Download in other formats: