Context Navigation

agrif_oce_sponge.F90 @ 12377

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Property svn:keywords set to Id

File size: 32.4 KB

Line
1	#define SPONGE && define SPONGE_TOP
2
3	MODULE agrif_oce_sponge
4	!!======================================================================
5	!! * MODULE agrif_oce_interp *
6	!! AGRIF: sponge package for the ocean dynamics (OPA)
7	!!======================================================================
8	!! History : 2.0 ! 2002-06 (XXX) Original cade
9	!! - ! 2005-11 (XXX)
10	!! 3.2 ! 2009-04 (R. Benshila)
11	!! 3.6 ! 2014-09 (R. Benshila)
12	!!----------------------------------------------------------------------
13	#if defined key_agrif
14	!!----------------------------------------------------------------------
15	!! 'key_agrif' AGRIF zoom
16	!!----------------------------------------------------------------------
17	USE par_oce
18	USE oce
19	USE dom_oce
20	!
21	USE in_out_manager
22	USE agrif_oce
23	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
24	USE iom
25	USE vremap
26
27	IMPLICIT NONE
28	PRIVATE
29
30	PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn
31	PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge
32
33	!! * Substitutions
34	# include "do_loop_substitute.h90"
35	!!----------------------------------------------------------------------
36	!! NEMO/NST 4.0 , NEMO Consortium (2018)
37	!! $Id$
38	!! Software governed by the CeCILL license (see ./LICENSE)
39	!!----------------------------------------------------------------------
40	CONTAINS
41
42	SUBROUTINE Agrif_Sponge_Tra
43	!!----------------------------------------------------------------------
44	!! * ROUTINE Agrif_Sponge_Tra *
45	!!----------------------------------------------------------------------
46	REAL(wp) :: zcoef ! local scalar
47
48	!!----------------------------------------------------------------------
49	!
50	#if defined SPONGE
51	!! Assume persistence:
52	zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
53
54	CALL Agrif_Sponge
55	Agrif_SpecialValue = 0._wp
56	Agrif_UseSpecialValue = .TRUE.
57	tabspongedone_tsn = .FALSE.
58	!
59	CALL Agrif_Bc_Variable( tsn_sponge_id, calledweight=zcoef, procname=interptsn_sponge )
60	!
61	Agrif_UseSpecialValue = .FALSE.
62	#endif
63	!
64	CALL iom_put( 'agrif_spu', fspu(:,:))
65	CALL iom_put( 'agrif_spv', fspv(:,:))
66	!
67	END SUBROUTINE Agrif_Sponge_Tra
68
69
70	SUBROUTINE Agrif_Sponge_dyn
71	!!----------------------------------------------------------------------
72	!! * ROUTINE Agrif_Sponge_dyn *
73	!!----------------------------------------------------------------------
74	REAL(wp) :: zcoef ! local scalar
75	!!----------------------------------------------------------------------
76	!
77	#if defined SPONGE
78	zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
79
80	Agrif_SpecialValue=0.
81	Agrif_UseSpecialValue = ln_spc_dyn
82	!
83	tabspongedone_u = .FALSE.
84	tabspongedone_v = .FALSE.
85	CALL Agrif_Bc_Variable( un_sponge_id, calledweight=zcoef, procname=interpun_sponge )
86	!
87	tabspongedone_u = .FALSE.
88	tabspongedone_v = .FALSE.
89	CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge )
90	!
91	Agrif_UseSpecialValue = .FALSE.
92	#endif
93	!
94	CALL iom_put( 'agrif_spt', fspt(:,:))
95	CALL iom_put( 'agrif_spf', fspf(:,:))
96	!
97	END SUBROUTINE Agrif_Sponge_dyn
98
99
100	SUBROUTINE Agrif_Sponge
101	!!----------------------------------------------------------------------
102	!! * ROUTINE Agrif_Sponge *
103	!!----------------------------------------------------------------------
104	INTEGER :: ji, jj, ind1, ind2
105	INTEGER :: ispongearea, jspongearea
106	REAL(wp) :: z1_ispongearea, z1_jspongearea
107	REAL(wp), DIMENSION(jpi,jpj) :: ztabramp
108	REAL(wp), DIMENSION(jpjmax) :: zmskwest, zmskeast
109	REAL(wp), DIMENSION(jpimax) :: zmsknorth, zmsksouth
110	!!----------------------------------------------------------------------
111	!
112	! Sponge 1d example with:
113	! iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2
114	!
115	!coarse : U T U T U T U
116	!\| \| \| \| \|
117	!fine : t u t u t u t u t u t u t u t u t u t u t
118	!sponge val:0 0 0 1 5/6 4/6 3/6 2/6 1/6 0 0
119	! \| ghost \| <-- sponge area -- > \|
120	! \| points \| \|
121	! \|--> dynamical interface
122
123	#if defined SPONGE \|\| defined SPONGE_TOP
124	IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
125	!
126	! Retrieve masks at open boundaries:
127
128	! --- West --- !
129	ztabramp(:,:) = 0._wp
130	ind1 = 1+nbghostcells
131	DO ji = mi0(ind1), mi1(ind1)
132	ztabramp(ji,:) = ssumask(ji,:)
133	END DO
134	!
135	zmskwest(:) = 0._wp
136	zmskwest(1:jpj) = MAXVAL(ztabramp(:,:), dim=1)
137
138	! --- East --- !
139	ztabramp(:,:) = 0._wp
140	ind1 = jpiglo - nbghostcells - 1
141	DO ji = mi0(ind1), mi1(ind1)
142	ztabramp(ji,:) = ssumask(ji,:)
143	END DO
144	!
145	zmskeast(:) = 0._wp
146	zmskeast(1:jpj) = MAXVAL(ztabramp(:,:), dim=1)
147
148	! --- South --- !
149	ztabramp(:,:) = 0._wp
150	ind1 = 1+nbghostcells
151	DO jj = mj0(ind1), mj1(ind1)
152	ztabramp(:,jj) = ssvmask(:,jj)
153	END DO
154	!
155	zmsksouth(:) = 0._wp
156	zmsksouth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2)
157
158	! --- North --- !
159	ztabramp(:,:) = 0._wp
160	ind1 = jpjglo - nbghostcells - 1
161	DO jj = mj0(ind1), mj1(ind1)
162	ztabramp(:,jj) = ssvmask(:,jj)
163	END DO
164	!
165	zmsknorth(:) = 0._wp
166	zmsknorth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2)
167	! JC: SPONGE MASKING TO BE SORTED OUT:
168	zmskwest(:) = 1._wp
169	zmskeast(:) = 1._wp
170	zmsknorth(:) = 1._wp
171	zmsksouth(:) = 1._wp
172	#if defined key_mpp_mpi
173	! CALL mpp_max( 'AGRIF_sponge', zmskwest(:) , jpjmax )
174	! CALL mpp_max( 'AGRIF_Sponge', zmskeast(:) , jpjmax )
175	! CALL mpp_max( 'AGRIF_Sponge', zmsksouth(:), jpimax )
176	! CALL mpp_max( 'AGRIF_Sponge', zmsknorth(:), jpimax )
177	#endif
178
179	! Define ramp from boundaries towards domain interior at T-points
180	! Store it in ztabramp
181
182	ispongearea = nn_sponge_len * Agrif_irhox()
183	z1_ispongearea = 1._wp / REAL( ispongearea )
184	jspongearea = nn_sponge_len * Agrif_irhoy()
185	z1_jspongearea = 1._wp / REAL( jspongearea )
186
187	ztabramp(:,:) = 0._wp
188
189	! Trick to remove sponge in 2DV domains:
190	IF ( nbcellsx <= 3 ) ispongearea = -1
191	IF ( nbcellsy <= 3 ) jspongearea = -1
192
193	! --- West --- !
194	ind1 = 1+nbghostcells
195	ind2 = 1+nbghostcells + ispongearea
196	DO ji = mi0(ind1), mi1(ind2)
197	DO jj = 1, jpj
198	ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj)
199	END DO
200	END DO
201
202	! ghost cells:
203	ind1 = 1
204	ind2 = nbghostcells + 1
205	DO ji = mi0(ind1), mi1(ind2)
206	DO jj = 1, jpj
207	ztabramp(ji,jj) = zmskwest(jj)
208	END DO
209	END DO
210
211	! --- East --- !
212	ind1 = jpiglo - nbghostcells - ispongearea
213	ind2 = jpiglo - nbghostcells
214	DO ji = mi0(ind1), mi1(ind2)
215	DO jj = 1, jpj
216	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj)
217	ENDDO
218	END DO
219
220	! ghost cells:
221	ind1 = jpiglo - nbghostcells
222	ind2 = jpiglo
223	DO ji = mi0(ind1), mi1(ind2)
224	DO jj = 1, jpj
225	ztabramp(ji,jj) = zmskeast(jj)
226	ENDDO
227	END DO
228
229	! --- South --- !
230	ind1 = 1+nbghostcells
231	ind2 = 1+nbghostcells + jspongearea
232	DO jj = mj0(ind1), mj1(ind2)
233	DO ji = 1, jpi
234	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea) * zmsksouth(ji)
235	END DO
236	END DO
237
238	! ghost cells:
239	ind1 = 1
240	ind2 = nbghostcells + 1
241	DO jj = mj0(ind1), mj1(ind2)
242	DO ji = 1, jpi
243	ztabramp(ji,jj) = zmsksouth(ji)
244	END DO
245	END DO
246
247	! --- North --- !
248	ind1 = jpjglo - nbghostcells - jspongearea
249	ind2 = jpjglo - nbghostcells
250	DO jj = mj0(ind1), mj1(ind2)
251	DO ji = 1, jpi
252	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mjg(jj) - ind1 ) * z1_jspongearea) * zmsknorth(ji)
253	END DO
254	END DO
255
256	! ghost cells:
257	ind1 = jpjglo - nbghostcells
258	ind2 = jpjglo
259	DO jj = mj0(ind1), mj1(ind2)
260	DO ji = 1, jpi
261	ztabramp(ji,jj) = zmsknorth(ji)
262	END DO
263	END DO
264
265	ENDIF
266
267	! Tracers
268	IF( .NOT. spongedoneT ) THEN
269	fspu(:,:) = 0._wp
270	fspv(:,:) = 0._wp
271	DO_2D_00_00
272	fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji+1,jj ) ) * ssumask(ji,jj)
273	fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1) ) * ssvmask(ji,jj)
274	END_2D
275	CALL lbc_lnk( 'agrif_Sponge', fspu, 'U', 1. ) ! Lateral boundary conditions
276	CALL lbc_lnk( 'agrif_Sponge', fspv, 'V', 1. )
277
278	spongedoneT = .TRUE.
279	ENDIF
280
281	! Dynamics
282	IF( .NOT. spongedoneU ) THEN
283	fspt(:,:) = 0._wp
284	fspf(:,:) = 0._wp
285	DO_2D_00_00
286	fspt(ji,jj) = ztabramp(ji,jj) * ssmask(ji,jj)
287	fspf(ji,jj) = 0.25_wp * ( ztabramp(ji ,jj ) + ztabramp(ji ,jj+1) &
288	& +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj ) ) &
289	& * ssvmask(ji,jj) * ssvmask(ji,jj+1)
290	END_2D
291	CALL lbc_lnk( 'agrif_Sponge', fspt, 'T', 1. ) ! Lateral boundary conditions
292	CALL lbc_lnk( 'agrif_Sponge', fspf, 'F', 1. )
293
294	spongedoneU = .TRUE.
295	ENDIF
296
297	#if defined key_vertical
298	! Remove vertical interpolation where not needed:
299	DO_2D_00_00
300	IF ((fspu(ji-1,jj)==0._wp).AND.(fspu(ji,jj)==0._wp).AND. &
301	& (fspv(ji,jj-1)==0._wp).AND.(fspv(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0
302	!
303	IF ((fspt(ji+1,jj)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. &
304	& (fspf(ji,jj-1)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbku_parent(ji,jj) = 0
305	!
306	IF ((fspt(ji,jj+1)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. &
307	& (fspf(ji-1,jj)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbkv_parent(ji,jj) = 0
308	!
309	IF ( ssmask(ji,jj) == 0._wp) mbkt_parent(ji,jj) = 0
310	IF (ssumask(ji,jj) == 0._wp) mbku_parent(ji,jj) = 0
311	IF (ssvmask(ji,jj) == 0._wp) mbkv_parent(ji,jj) = 0
312	END_2D
313	!
314	ztabramp(:,:) = REAL( mbkt_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'T', 1. )
315	mbkt_parent(:,:) = NINT( ztabramp(:,:) )
316	ztabramp(:,:) = REAL( mbku_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'U', 1. )
317	mbku_parent(:,:) = NINT( ztabramp(:,:) )
318	ztabramp(:,:) = REAL( mbkv_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'V', 1. )
319	mbkv_parent(:,:) = NINT( ztabramp(:,:) )
320	#endif
321	!
322	#endif
323	!
324	END SUBROUTINE Agrif_Sponge
325
326	SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
327	!!----------------------------------------------------------------------
328	!! * ROUTINE interptsn_sponge *
329	!!----------------------------------------------------------------------
330	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
331	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
332	LOGICAL , INTENT(in ) :: before
333	!
334	INTEGER :: ji, jj, jk, jn ! dummy loop indices
335	INTEGER :: iku, ikv
336	REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot, ztrelax
337	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv
338	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
339	! vertical interpolation:
340	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
341	REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
342	REAL(wp), DIMENSION(k1:k2) :: h_in
343	REAL(wp), DIMENSION(1:jpk) :: h_out
344	INTEGER :: N_in, N_out
345	!!----------------------------------------------------------------------
346	!
347	IF( before ) THEN
348	DO jn = 1, jpts
349	DO jk=k1,k2
350	DO jj=j1,j2
351	DO ji=i1,i2
352	tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a)
353	END DO
354	END DO
355	END DO
356	END DO
357
358	# if defined key_vertical
359	! Interpolate thicknesses
360	! Warning: these are masked, hence extrapolated prior interpolation.
361	DO jk=k1,k2
362	DO jj=j1,j2
363	DO ji=i1,i2
364	tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a)
365	END DO
366	END DO
367	END DO
368
369	! Extrapolate thicknesses in partial bottom cells:
370	! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
371	IF (ln_zps) THEN
372	DO jj=j1,j2
373	DO ji=i1,i2
374	jk = mbkt(ji,jj)
375	tabres(ji,jj,jk,jpts+1) = 0._wp
376	END DO
377	END DO
378	END IF
379
380	! Save ssh at last level:
381	IF (.NOT.ln_linssh) THEN
382	tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1)
383	ELSE
384	tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp
385	END IF
386	# endif
387
388	ELSE
389	!
390	# if defined key_vertical
391
392	IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp
393
394	DO jj=j1,j2
395	DO ji=i1,i2
396	tabres_child(ji,jj,:,:) = 0._wp
397	N_in = mbkt_parent(ji,jj)
398	zhtot = 0._wp
399	DO jk=1,N_in !k2 = jpk of parent grid
400	IF (jk==N_in) THEN
401	h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot
402	ELSE
403	h_in(jk) = tabres(ji,jj,jk,n2)
404	ENDIF
405	zhtot = zhtot + h_in(jk)
406	tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)
407	END DO
408	N_out = 0
409	DO jk=1,jpk ! jpk of child grid
410	IF (tmask(ji,jj,jk) == 0) EXIT
411	N_out = N_out + 1
412	h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above
413	ENDDO
414
415	! Account for small differences in free-surface
416	IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
417	h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
418	ELSE
419	h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
420	ENDIF
421	IF (N_in*N_out > 0) THEN
422	CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
423	ENDIF
424	ENDDO
425	ENDDO
426	# endif
427
428	DO jj=j1,j2
429	DO ji=i1,i2
430	DO jk=1,jpkm1
431	# if defined key_vertical
432	tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk)
433	# else
434	tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk)
435	# endif
436	ENDDO
437	ENDDO
438	ENDDO
439
440	!* set relaxation time scale
441	IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_tra / ( rdt )
442	ELSE ; ztrelax = rn_trelax_tra / (2._wp * rdt )
443	ENDIF
444
445	DO jn = 1, jpts
446	DO jk = 1, jpkm1
447	ztu(i1:i2,j1:j2,jk) = 0._wp
448	DO jj = j1,j2
449	DO ji = i1,i2-1
450	zabe1 = rn_sponge_tra * fspu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
451	ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) )
452	END DO
453	END DO
454	ztv(i1:i2,j1:j2,jk) = 0._wp
455	DO ji = i1,i2
456	DO jj = j1,j2-1
457	zabe2 = rn_sponge_tra * fspv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
458	ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
459	END DO
460	END DO
461	!
462	IF( ln_zps ) THEN ! set gradient at partial step level
463	DO jj = j1,j2
464	DO ji = i1,i2
465	! last level
466	iku = mbku(ji,jj)
467	ikv = mbkv(ji,jj)
468	IF( iku == jk ) ztu(ji,jj,jk) = 0._wp
469	IF( ikv == jk ) ztv(ji,jj,jk) = 0._wp
470	END DO
471	END DO
472	ENDIF
473	END DO
474	!
475	DO jk = 1, jpkm1
476	DO jj = j1+1,j2-1
477	DO ji = i1+1,i2-1
478	IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
479	zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
480	! horizontal diffusive trends
481	ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) &
482	& - ztrelax * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn)
483	! add it to the general tracer trends
484	ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa
485	ENDIF
486	END DO
487	END DO
488	END DO
489	!
490	END DO
491	!
492	tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
493	!
494	ENDIF
495	!
496	END SUBROUTINE interptsn_sponge
497
498	SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
499	!!---------------------------------------------
500	!! * ROUTINE interpun_sponge *
501	!!---------------------------------------------
502	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
503	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
504	LOGICAL, INTENT(in) :: before
505
506	INTEGER :: ji,jj,jk,jmax
507
508	! sponge parameters
509	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax
510	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
511	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
512	! vertical interpolation:
513	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
514	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
515	REAL(wp), DIMENSION(1:jpk) :: h_out
516	INTEGER ::N_in, N_out
517	!!---------------------------------------------
518	!
519	IF( before ) THEN
520	DO jk=k1,k2
521	DO jj=j1,j2
522	DO ji=i1,i2
523	tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a)
524	# if defined key_vertical
525	tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk)
526	# endif
527	END DO
528	END DO
529	END DO
530
531	# if defined key_vertical
532	! Extrapolate thicknesses in partial bottom cells:
533	! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
534	IF (ln_zps) THEN
535	DO jj=j1,j2
536	DO ji=i1,i2
537	jk = mbku(ji,jj)
538	tabres(ji,jj,jk,m2) = 0._wp
539	END DO
540	END DO
541	END IF
542	! Save ssh at last level:
543	tabres(i1:i2,j1:j2,k2,m2) = 0._wp
544	IF (.NOT.ln_linssh) THEN
545	! This vertical sum below should be replaced by the sea-level at U-points (optimization):
546	DO jk=1,jpk
547	tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u(i1:i2,j1:j2,jk,Kbb_a) * umask(i1:i2,j1:j2,jk)
548	END DO
549	tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2)
550	END IF
551	# endif
552
553	ELSE
554
555	# if defined key_vertical
556	IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
557
558	DO jj=j1,j2
559	DO ji=i1,i2
560	tabres_child(ji,jj,:) = 0._wp
561	N_in = mbku_parent(ji,jj)
562	zhtot = 0._wp
563	DO jk=1,N_in
564	IF (jk==N_in) THEN
565	h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
566	ELSE
567	h_in(jk) = tabres(ji,jj,jk,m2)
568	ENDIF
569	zhtot = zhtot + h_in(jk)
570	tabin(jk) = tabres(ji,jj,jk,m1)
571	ENDDO
572	!
573	N_out = 0
574	DO jk=1,jpk
575	IF (umask(ji,jj,jk) == 0) EXIT
576	N_out = N_out + 1
577	h_out(N_out) = e3u(ji,jj,jk,Kbb_a)
578	ENDDO
579
580	! Account for small differences in free-surface
581	IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
582	h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
583	ELSE
584	h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
585	ENDIF
586
587	IF (N_in * N_out > 0) THEN
588	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,1)
589	ENDIF
590	ENDDO
591	ENDDO
592
593	ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
594	#else
595	ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
596	#endif
597	!* set relaxation time scale
598	IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_dyn / ( rdt )
599	ELSE ; ztrelax = rn_trelax_dyn / (2._wp * rdt )
600	ENDIF
601	!
602	DO jk = 1, jpkm1 ! Horizontal slab
603	! ! ===============
604
605	! ! --------
606	! Horizontal divergence ! div
607	! ! --------
608	DO jj = j1,j2
609	DO ji = i1+1,i2 ! vector opt.
610	zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj)
611	hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)e3u(ji ,jj,jk,Kbb_a) ubdiff(ji ,jj,jk) &
612	& -e2u(ji-1,jj)e3u(ji-1,jj,jk,Kbb_a) ubdiff(ji-1,jj,jk) ) * zbtr
613	END DO
614	END DO
615
616	DO jj = j1,j2-1
617	DO ji = i1,i2 ! vector opt.
618	zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj)
619	rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &
620	& +e1u(ji,jj ) * ubdiff(ji,jj ,jk) ) * fmask(ji,jj,jk) * zbtr
621	END DO
622	END DO
623	END DO
624	!
625	DO jj = j1+1, j2-1
626	DO ji = i1+1, i2-1 ! vector opt.
627
628	IF (.NOT. tabspongedone_u(ji,jj)) THEN
629	DO jk = 1, jpkm1 ! Horizontal slab
630	ze2u = rotdiff (ji,jj,jk)
631	ze1v = hdivdiff(ji,jj,jk)
632	! horizontal diffusive trends
633	zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) &
634	& + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) &
635	& - ztrelax * fspu(ji,jj) * ubdiff(ji,jj,jk)
636
637	! add it to the general momentum trends
638	uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua
639	END DO
640	ENDIF
641
642	END DO
643	END DO
644
645	tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
646
647	jmax = j2-1
648	IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-nbghostcells-2) ! North
649
650	DO jj = j1+1, jmax
651	DO ji = i1+1, i2 ! vector opt.
652
653	IF (.NOT. tabspongedone_v(ji,jj)) THEN
654	DO jk = 1, jpkm1 ! Horizontal slab
655	ze2u = rotdiff (ji,jj,jk)
656	ze1v = hdivdiff(ji,jj,jk)
657
658	! horizontal diffusive trends
659	zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) ) &
660	+ ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
661
662	! add it to the general momentum trends
663	vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) + zva
664	END DO
665	ENDIF
666	!
667	END DO
668	END DO
669	!
670	tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
671	!
672	ENDIF
673	!
674	END SUBROUTINE interpun_sponge
675
676	SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir)
677	!!---------------------------------------------
678	!! * ROUTINE interpvn_sponge *
679	!!---------------------------------------------
680	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
681	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
682	LOGICAL, INTENT(in) :: before
683	INTEGER, INTENT(in) :: nb , ndir
684	!
685	INTEGER :: ji, jj, jk, imax
686	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax
687	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
688	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
689	! vertical interpolation:
690	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
691	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
692	REAL(wp), DIMENSION(1:jpk) :: h_out
693	INTEGER :: N_in, N_out
694	!!---------------------------------------------
695
696	IF( before ) THEN
697	DO jk=k1,k2
698	DO jj=j1,j2
699	DO ji=i1,i2
700	tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a)
701	# if defined key_vertical
702	tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_a)
703	# endif
704	END DO
705	END DO
706	END DO
707
708	# if defined key_vertical
709	! Extrapolate thicknesses in partial bottom cells:
710	! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
711	IF (ln_zps) THEN
712	DO jj=j1,j2
713	DO ji=i1,i2
714	jk = mbkv(ji,jj)
715	tabres(ji,jj,jk,m2) = 0._wp
716	END DO
717	END DO
718	END IF
719	! Save ssh at last level:
720	tabres(i1:i2,j1:j2,k2,m2) = 0._wp
721	IF (.NOT.ln_linssh) THEN
722	! This vertical sum below should be replaced by the sea-level at V-points (optimization):
723	DO jk=1,jpk
724	tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v(i1:i2,j1:j2,jk,Kbb_a) * vmask(i1:i2,j1:j2,jk)
725	END DO
726	tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2)
727	END IF
728	# endif
729
730	ELSE
731
732	# if defined key_vertical
733	IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
734	DO jj=j1,j2
735	DO ji=i1,i2
736	tabres_child(ji,jj,:) = 0._wp
737	N_in = mbkv_parent(ji,jj)
738	zhtot = 0._wp
739	DO jk=1,N_in
740	IF (jk==N_in) THEN
741	h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
742	ELSE
743	h_in(jk) = tabres(ji,jj,jk,m2)
744	ENDIF
745	zhtot = zhtot + h_in(jk)
746	tabin(jk) = tabres(ji,jj,jk,m1)
747	ENDDO
748	!
749	N_out = 0
750	DO jk=1,jpk
751	IF (vmask(ji,jj,jk) == 0) EXIT
752	N_out = N_out + 1
753	h_out(N_out) = e3v(ji,jj,jk,Kbb_a)
754	ENDDO
755
756	! Account for small differences in free-surface
757	IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
758	h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
759	ELSE
760	h_in(1) = h_in(1) - ( sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
761	ENDIF
762
763	IF (N_in * N_out > 0) THEN
764	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,1)
765	ENDIF
766	ENDDO
767	ENDDO
768
769	vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)
770	# else
771	vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
772	# endif
773	!* set relaxation time scale
774	IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_dyn / ( rdt )
775	ELSE ; ztrelax = rn_trelax_dyn / (2._wp * rdt )
776	ENDIF
777	!
778	DO jk = 1, jpkm1 ! Horizontal slab
779	! ! ===============
780
781	! ! --------
782	! Horizontal divergence ! div
783	! ! --------
784	DO jj = j1+1,j2
785	DO ji = i1,i2 ! vector opt.
786	zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj)
787	hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v(ji,jj ,jk,Kbb_a) * vbdiff(ji,jj ,jk) &
788	& -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kbb_a) * vbdiff(ji,jj-1,jk) ) * zbtr
789	END DO
790	END DO
791	DO jj = j1,j2
792	DO ji = i1,i2-1 ! vector opt.
793	zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj)
794	rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) &
795	& -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr
796	END DO
797	END DO
798	END DO
799
800	! ! ===============
801	!
802
803	imax = i2 - 1
804	IF ((nbondi == 1).OR.(nbondi == 2)) imax = MIN(imax,nlci-nbghostcells-2) ! East
805
806	DO jj = j1+1, j2
807	DO ji = i1+1, imax ! vector opt.
808	IF( .NOT. tabspongedone_u(ji,jj) ) THEN
809	DO jk = 1, jpkm1
810	uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) &
811	& - ( rotdiff (ji ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) &
812	& + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj ,jk)) * r1_e1u(ji,jj)
813	END DO
814	ENDIF
815	END DO
816	END DO
817	!
818	tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
819	!
820	DO jj = j1+1, j2-1
821	DO ji = i1+1, i2-1 ! vector opt.
822	IF( .NOT. tabspongedone_v(ji,jj) ) THEN
823	DO jk = 1, jpkm1
824	vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) &
825	& + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) ) &
826	& + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj) &
827	& - ztrelax * fspv(ji,jj) * vbdiff(ji,jj,jk)
828	END DO
829	ENDIF
830	END DO
831	END DO
832	tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
833	ENDIF
834	!
835	END SUBROUTINE interpvn_sponge
836
837	#else
838	!!----------------------------------------------------------------------
839	!! Empty module no AGRIF zoom
840	!!----------------------------------------------------------------------
841	CONTAINS
842	SUBROUTINE agrif_oce_sponge_empty
843	WRITE(,) 'agrif_oce_sponge : You should not have seen this print! error?'
844	END SUBROUTINE agrif_oce_sponge_empty
845	#endif
846
847	!!======================================================================
848	END MODULE agrif_oce_sponge

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

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

Context Navigation

source: NEMO/trunk/src/NST/agrif_oce_sponge.F90 @ 12377

Download in other formats: