Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

agrif_opa_sponge.F90 @ 9006

Last change on this file since 9006 was 8998, checked in by timgraham, 6 years ago
First commit of Jerome's modified versions of agrif_opa routines
Property svn:keywords set to `Id`
File size: 22.2 KB

Line
1	#define SPONGE && define 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	!! Assume persistence:
42	timecoeff = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
43
44	CALL Agrif_Sponge
45	Agrif_SpecialValue=0.
46	Agrif_UseSpecialValue = .TRUE.
47	tabspongedone_tsn = .FALSE.
48
49	CALL Agrif_Bc_Variable(tsn_sponge_id,calledweight=timecoeff,procname=interptsn_sponge)
50
51	Agrif_UseSpecialValue = .FALSE.
52	#endif
53	!
54	END SUBROUTINE Agrif_Sponge_Tra
55
56
57	SUBROUTINE Agrif_Sponge_dyn
58	!!---------------------------------------------
59	!! * ROUTINE Agrif_Sponge_dyn *
60	!!---------------------------------------------
61	REAL(wp) :: timecoeff
62	!!---------------------------------------------
63
64	#if defined SPONGE
65	!! timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
66	!! Assume persistence:
67	timecoeff = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
68
69	Agrif_SpecialValue=0.
70	Agrif_UseSpecialValue = ln_spc_dyn
71
72	tabspongedone_u = .FALSE.
73	tabspongedone_v = .FALSE.
74	CALL Agrif_Bc_Variable(un_sponge_id,calledweight=timecoeff,procname=interpun_sponge)
75
76	tabspongedone_u = .FALSE.
77	tabspongedone_v = .FALSE.
78	CALL Agrif_Bc_Variable(vn_sponge_id,calledweight=timecoeff,procname=interpvn_sponge)
79
80	Agrif_UseSpecialValue = .FALSE.
81	#endif
82	!
83	END SUBROUTINE Agrif_Sponge_dyn
84
85
86	SUBROUTINE Agrif_Sponge
87	!!---------------------------------------------
88	!! * ROUTINE Agrif_Sponge *
89	!!---------------------------------------------
90	INTEGER :: ji,jj,jk
91	INTEGER :: ispongearea, ilci, ilcj
92	LOGICAL :: ll_spdone
93	REAL(wp) :: z1spongearea, zramp
94	REAL(wp), POINTER, DIMENSION(:,:) :: ztabramp
95
96	#if defined SPONGE \|\| defined SPONGE_TOP
97	ll_spdone=.TRUE.
98	IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
99	! Define ramp from boundaries towards domain interior
100	! at T-points
101	! Store it in ztabramp
102	ll_spdone=.FALSE.
103
104	CALL wrk_alloc( jpi, jpj, ztabramp )
105
106	ispongearea = 2 + nn_sponge_len * Agrif_irhox()
107	ilci = nlci - ispongearea
108	ilcj = nlcj - ispongearea
109	z1spongearea = 1._wp / REAL( ispongearea - 2 )
110
111	ztabramp(:,:) = 0._wp
112
113	IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
114	DO jj = 1, jpj
115	IF ( umask(2,jj,1) == 1._wp ) THEN
116	DO ji = 2, ispongearea
117	ztabramp(ji,jj) = ( ispongearea-ji ) * z1spongearea
118	END DO
119	ENDIF
120	ENDDO
121	ENDIF
122
123	IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
124	DO jj = 1, jpj
125	IF ( umask(nlci-2,jj,1) == 1._wp ) THEN
126	DO ji = ilci+1,nlci-1
127	zramp = (ji - (ilci+1) ) * z1spongearea
128	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
129	ENDDO
130	ENDIF
131	ENDDO
132	ENDIF
133
134	IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
135	DO ji = 1, jpi
136	IF ( vmask(ji,2,1) == 1._wp ) THEN
137	DO jj = 2, ispongearea
138	zramp = ( ispongearea-jj ) * z1spongearea
139	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
140	END DO
141	ENDIF
142	ENDDO
143	ENDIF
144
145	IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
146	DO ji = 1, jpi
147	IF ( vmask(ji,nlcj-2,1) == 1._wp ) THEN
148	DO jj = ilcj+1,nlcj-1
149	zramp = (jj - (ilcj+1) ) * z1spongearea
150	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
151	END DO
152	ENDIF
153	ENDDO
154	ENDIF
155
156	ENDIF
157
158	! Tracers
159	IF( .NOT. spongedoneT ) THEN
160	fsaht_spu(:,:) = 0._wp
161	fsaht_spv(:,:) = 0._wp
162	DO jj = 2, jpjm1
163	DO ji = 2, jpim1 ! vector opt.
164	fsaht_spu(ji,jj) = 0.5_wp * visc_tra * (ztabramp(ji,jj) + ztabramp(ji+1,jj ))
165	fsaht_spv(ji,jj) = 0.5_wp * visc_tra * (ztabramp(ji,jj) + ztabramp(ji ,jj+1))
166	END DO
167	END DO
168
169	CALL lbc_lnk( fsaht_spu, 'U', 1. ) ! Lateral boundary conditions
170	CALL lbc_lnk( fsaht_spv, 'V', 1. )
171	spongedoneT = .TRUE.
172	ENDIF
173
174	! Dynamics
175	IF( .NOT. spongedoneU ) THEN
176	fsahm_spt(:,:) = 0._wp
177	fsahm_spf(:,:) = 0._wp
178	DO jj = 2, jpjm1
179	DO ji = 2, jpim1 ! vector opt.
180	fsahm_spt(ji,jj) = visc_dyn * ztabramp(ji,jj)
181	fsahm_spf(ji,jj) = 0.25_wp * visc_dyn * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1) &
182	& +ztabramp(ji,jj) + ztabramp(ji+1,jj ) )
183	END DO
184	END DO
185
186	CALL lbc_lnk( fsahm_spt, 'T', 1. ) ! Lateral boundary conditions
187	CALL lbc_lnk( fsahm_spf, 'F', 1. )
188	spongedoneU = .TRUE.
189	ENDIF
190	!
191	IF (.NOT.ll_spdone) CALL wrk_dealloc( jpi, jpj, ztabramp )
192	!
193	#endif
194	!
195	END SUBROUTINE Agrif_Sponge
196
197	SUBROUTINE interptsn_sponge(tabres,i1,i2,j1,j2,k1,k2,n1,n2,before)
198	!!---------------------------------------------
199	!! * ROUTINE interptsn_sponge *
200	!!---------------------------------------------
201	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
202	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
203	LOGICAL, INTENT(in) :: before
204	!
205	INTEGER :: ji, jj, jk, jn ! dummy loop indices
206	INTEGER :: iku, ikv
207	REAL(wp) :: ztsa, zabe1, zabe2, zbtr
208	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv
209	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
210	! vertical interpolation:
211	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
212	REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
213	REAL(wp), DIMENSION(k1:k2) :: h_in
214	REAL(wp), DIMENSION(1:jpk) :: h_out
215	INTEGER :: N_in, N_out
216	REAL(wp) :: h_diff
217	!
218	IF( before ) THEN
219	DO jn = 1, jpts
220	DO jk=k1,k2
221	DO jj=j1,j2
222	DO ji=i1,i2
223	tabres(ji,jj,jk,jn) = tsb(ji,jj,jk,jn)
224	END DO
225	END DO
226	END DO
227	END DO
228
229	# if defined key_vertical
230	DO jk=k1,k2
231	DO jj=j1,j2
232	DO ji=i1,i2
233	tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)
234	END DO
235	END DO
236	END DO
237	# endif
238
239	ELSE
240	!
241	# if defined key_vertical
242	tabres_child(:,:,:,:) = 0.
243	DO jj=j1,j2
244	DO ji=i1,i2
245	N_in = 0
246	DO jk=k1,k2 !k2 = jpk of parent grid
247	IF (tabres(ji,jj,jk,n2) == 0) EXIT
248	N_in = N_in + 1
249	tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)
250	h_in(N_in) = tabres(ji,jj,jk,n2)
251	END DO
252	N_out = 0
253	DO jk=1,jpk ! jpk of child grid
254	IF (tmask(ji,jj,jk) == 0) EXIT
255	N_out = N_out + 1
256	h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above
257	ENDDO
258	IF (N_in > 0) THEN
259	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
260	tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for?????
261	DO jn=1,jpts
262	call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
263	ENDDO
264	ENDIF
265	ENDDO
266	ENDDO
267	# endif
268
269	DO jj=j1,j2
270	DO ji=i1,i2
271	DO jk=1,jpkm1
272	# if defined key_vertical
273	tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)
274	# else
275	tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts)
276	# endif
277	ENDDO
278	ENDDO
279	ENDDO
280
281	DO jn = 1, jpts
282	DO jk = 1, jpkm1
283	DO jj = j1,j2-1
284	DO ji = i1,i2-1
285	zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
286	zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
287	ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) )
288	ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
289	END DO
290	END DO
291	!
292	IF( ln_zps ) THEN ! set gradient at partial step level
293	DO jj = j1,j2-1
294	DO ji = i1,i2-1
295	! last level
296	iku = mbku(ji,jj)
297	ikv = mbkv(ji,jj)
298	IF( iku == jk ) ztu(ji,jj,jk) = 0._wp
299	IF( ikv == jk ) ztv(ji,jj,jk) = 0._wp
300	END DO
301	END DO
302	ENDIF
303	END DO
304	!
305	DO jk = 1, jpkm1
306	DO jj = j1+1,j2-1
307	DO ji = i1+1,i2-1
308	IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
309	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
310	! horizontal diffusive trends
311	ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) )
312	! add it to the general tracer trends
313	tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa
314	ENDIF
315	END DO
316	END DO
317	END DO
318	!
319	END DO
320	!
321	tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
322	!
323	ENDIF
324	!
325	END SUBROUTINE interptsn_sponge
326
327	SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
328	!!---------------------------------------------
329	!! * ROUTINE interpun_sponge *
330	!!---------------------------------------------
331	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
332	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
333	LOGICAL, INTENT(in) :: before
334
335	INTEGER :: ji,jj,jk,jmax
336
337	! sponge parameters
338	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff
339	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
340	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
341	! vertical interpolation:
342	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
343	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
344	REAL(wp), DIMENSION(1:jpk) :: h_out
345	INTEGER ::N_in,N_out
346	!!---------------------------------------------
347	!
348	IF( before ) THEN
349	DO jk=1,jpkm1
350	DO jj=j1,j2
351	DO ji=i1,i2
352	tabres(ji,jj,jk,m1) = ub(ji,jj,jk)
353	# if defined key_vertical
354	tabres(ji,jj,jk,m2) = e3u_n(ji,jj,jk)*umask(ji,jj,jk)
355	# endif
356	END DO
357	END DO
358	END DO
359
360	ELSE
361
362	# if defined key_vertical
363	tabres_child(:,:,:) = 0._wp
364	DO jj=j1,j2
365	DO ji=i1,i2
366	N_in = 0
367	DO jk=k1,k2
368	IF (tabres(ji,jj,jk,m2) == 0) EXIT
369	N_in = N_in + 1
370	tabin(jk) = tabres(ji,jj,jk,m1)
371	h_in(N_in) = tabres(ji,jj,jk,m2)
372	ENDDO
373	!
374	IF (N_in == 0) THEN
375	tabres_child(ji,jj,:) = 0.
376	CYCLE
377	ENDIF
378
379	N_out = 0
380	DO jk=1,jpk
381	if (umask(ji,jj,jk) == 0) EXIT
382	N_out = N_out + 1
383	h_out(N_out) = e3u_n(ji,jj,jk)
384	ENDDO
385
386	IF (N_out == 0) THEN
387	tabres_child(ji,jj,:) = 0.
388	CYCLE
389	ENDIF
390
391	IF (N_in * N_out > 0) THEN
392	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
393	if (h_diff < -1.e4) then
394	print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
395	endif
396	ENDIF
397	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)
398
399	ENDDO
400	ENDDO
401
402	ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
403	#else
404	ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*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	SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir)
483	!!---------------------------------------------
484	!! * ROUTINE interpvn_sponge *
485	!!---------------------------------------------
486	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
487	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
488	LOGICAL, INTENT(in) :: before
489	INTEGER, INTENT(in) :: nb , ndir
490	!
491	INTEGER :: ji, jj, jk, imax
492	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff
493	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
494	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
495	! vertical interpolation:
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 :: N_in, N_out
500	!!---------------------------------------------
501
502	IF( before ) THEN
503	DO jk=1,jpkm1
504	DO jj=j1,j2
505	DO ji=i1,i2
506	tabres(ji,jj,jk,m1) = vb(ji,jj,jk)
507	# if defined key_vertical
508	tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_n(ji,jj,jk)
509	# endif
510	END DO
511	END DO
512	END DO
513	ELSE
514
515	# if defined key_vertical
516	tabres_child(:,:,:) = 0._wp
517	DO jj=j1,j2
518	DO ji=i1,i2
519	N_in = 0
520	DO jk=k1,k2
521	IF (tabres(ji,jj,jk,m2) == 0) EXIT
522	N_in = N_in + 1
523	tabin(jk) = tabres(ji,jj,jk,m1)
524	h_in(N_in) = tabres(ji,jj,jk,m2)
525	ENDDO
526
527	IF (N_in == 0) THEN
528	tabres_child(ji,jj,:) = 0.
529	CYCLE
530	ENDIF
531
532	N_out = 0
533	DO jk=1,jpk
534	if (vmask(ji,jj,jk) == 0) EXIT
535	N_out = N_out + 1
536	h_out(N_out) = e3v_n(ji,jj,jk)
537	ENDDO
538
539	IF (N_in * N_out > 0) THEN
540	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
541	if (h_diff < -1.e4) then
542	print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
543	endif
544	ENDIF
545	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)
546	ENDDO
547	ENDDO
548
549	vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)
550	# else
551	vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
552	# endif
553	!
554	DO jk = 1, jpkm1 ! Horizontal slab
555	! ! ===============
556
557	! ! --------
558	! Horizontal divergence ! div
559	! ! --------
560	DO jj = j1+1,j2
561	DO ji = i1,i2 ! vector opt.
562	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
563	hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v_n(ji,jj ,jk) * vbdiff(ji,jj ,jk) &
564	& -e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk) ) * zbtr
565	END DO
566	END DO
567	DO jj = j1,j2
568	DO ji = i1,i2-1 ! vector opt.
569	zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
570	rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) &
571	& -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr
572	END DO
573	END DO
574	END DO
575
576	! ! ===============
577	!
578
579	imax = i2-1
580	IF ((nbondi == 1).OR.(nbondi == 2)) imax = MIN(imax,nlci-3)
581
582	DO jj = j1+1, j2
583	DO ji = i1+1, imax ! vector opt.
584	IF( .NOT. tabspongedone_u(ji,jj) ) THEN
585	DO jk = 1, jpkm1
586	ua(ji,jj,jk) = ua(ji,jj,jk) &
587	& - ( rotdiff (ji ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) ) &
588	& + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj ,jk)) * r1_e1u(ji,jj)
589	END DO
590	ENDIF
591	END DO
592	END DO
593	!
594	tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
595	!
596	DO jj = j1+1, j2-1
597	DO ji = i1+1, i2-1 ! vector opt.
598	IF( .NOT. tabspongedone_v(ji,jj) ) THEN
599	DO jk = 1, jpkm1
600	va(ji,jj,jk) = va(ji,jj,jk) &
601	& + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) ) &
602	& + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj)
603	END DO
604	ENDIF
605	END DO
606	END DO
607	tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
608	ENDIF
609	!
610	END SUBROUTINE interpvn_sponge
611
612	#else
613	CONTAINS
614	SUBROUTINE agrif_opa_sponge_empty
615	!!---------------------------------------------
616	!! * ROUTINE agrif_OPA_sponge_empty *
617	!!---------------------------------------------
618	WRITE(,) 'agrif_opa_sponge : You should not have seen this print! error?'
619	END SUBROUTINE agrif_opa_sponge_empty
620	#endif
621
622	!!======================================================================
623	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_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 9006

Download in other formats: