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_oce_interp.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_interp.F90 @ 11428

Last change on this file since 11428 was 11428, checked in by acc, 4 years ago
Branch 2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps. Prettification of NST routines to reduce long lines and incease code readability. No functional changes. SETTE tested
Property svn:keywords set to `Id`
File size: 65.1 KB

Line
1	MODULE agrif_oce_interp
2	!!======================================================================
3	!! * MODULE agrif_oce_interp *
4	!! AGRIF: interpolation package for the ocean dynamics (OPA)
5	!!======================================================================
6	!! History : 2.0 ! 2002-06 (L. Debreu) Original cade
7	!! 3.2 ! 2009-04 (R. Benshila)
8	!! 3.6 ! 2014-09 (R. Benshila)
9	!!----------------------------------------------------------------------
10	#if defined key_agrif
11	!!----------------------------------------------------------------------
12	!! 'key_agrif' AGRIF zoom
13	!!----------------------------------------------------------------------
14	!! Agrif_tra :
15	!! Agrif_dyn :
16	!! Agrif_ssh :
17	!! Agrif_dyn_ts :
18	!! Agrif_dta_ts :
19	!! Agrif_ssh_ts :
20	!! Agrif_avm :
21	!! interpu :
22	!! interpv :
23	!!----------------------------------------------------------------------
24	USE par_oce
25	USE oce
26	USE dom_oce
27	USE zdf_oce
28	USE agrif_oce
29	USE phycst
30	USE dynspg_ts, ONLY: un_adv, vn_adv
31	!
32	USE in_out_manager
33	USE agrif_oce_sponge
34	USE lib_mpp
35
36	IMPLICIT NONE
37	PRIVATE
38
39	PUBLIC Agrif_dyn, Agrif_ssh, Agrif_dyn_ts, Agrif_ssh_ts, Agrif_dta_ts
40	PUBLIC Agrif_tra, Agrif_avm
41	PUBLIC interpun , interpvn
42	PUBLIC interptsn, interpsshn, interpavm
43	PUBLIC interpunb, interpvnb , interpub2b, interpvb2b
44	PUBLIC interpe3t, interpumsk, interpvmsk
45
46	INTEGER :: bdy_tinterp = 0
47
48	# include "vectopt_loop_substitute.h90"
49	!!----------------------------------------------------------------------
50	!! NEMO/NST 4.0 , NEMO Consortium (2018)
51	!! $Id$
52	!! Software governed by the CeCILL license (see ./LICENSE)
53	!!----------------------------------------------------------------------
54	CONTAINS
55
56	SUBROUTINE Agrif_tra
57	!!----------------------------------------------------------------------
58	!! * ROUTINE Agrif_tra *
59	!!----------------------------------------------------------------------
60	!
61	IF( Agrif_Root() ) RETURN
62	!
63	Agrif_SpecialValue = 0._wp
64	Agrif_UseSpecialValue = .TRUE.
65	!
66	CALL Agrif_Bc_variable( tsn_id, procname=interptsn )
67	!
68	Agrif_UseSpecialValue = .FALSE.
69	!
70	END SUBROUTINE Agrif_tra
71
72
73	SUBROUTINE Agrif_dyn( kt )
74	!!----------------------------------------------------------------------
75	!! * ROUTINE Agrif_DYN *
76	!!----------------------------------------------------------------------
77	INTEGER, INTENT(in) :: kt
78	!
79	INTEGER :: ji, jj, jk ! dummy loop indices
80	INTEGER :: j1, j2, i1, i2
81	INTEGER :: ibdy1, jbdy1, ibdy2, jbdy2
82	REAL(wp), DIMENSION(jpi,jpj) :: zub, zvb
83	!!----------------------------------------------------------------------
84	!
85	IF( Agrif_Root() ) RETURN
86	!
87	Agrif_SpecialValue = 0._wp
88	Agrif_UseSpecialValue = ln_spc_dyn
89	!
90	CALL Agrif_Bc_variable( un_interp_id, procname=interpun )
91	CALL Agrif_Bc_variable( vn_interp_id, procname=interpvn )
92	!
93	Agrif_UseSpecialValue = .FALSE.
94	!
95	! prevent smoothing in ghost cells
96	i1 = 1 ; i2 = nlci
97	j1 = 1 ; j2 = nlcj
98	IF( nbondj == -1 .OR. nbondj == 2 ) j1 = 2 + nbghostcells
99	IF( nbondj == +1 .OR. nbondj == 2 ) j2 = nlcj - nbghostcells - 1
100	IF( nbondi == -1 .OR. nbondi == 2 ) i1 = 2 + nbghostcells
101	IF( nbondi == +1 .OR. nbondi == 2 ) i2 = nlci - nbghostcells - 1
102
103	! --- West --- !
104	IF( nbondi == -1 .OR. nbondi == 2 ) THEN
105	ibdy1 = 2
106	ibdy2 = 1+nbghostcells
107	!
108	IF( .NOT.ln_dynspg_ts ) THEN ! Store transport
109	uu_b(ibdy1:ibdy2,:,Krhs_a) = 0._wp
110	DO jk = 1, jpkm1
111	DO jj = 1, jpj
112	uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj, Krhs_a) &
113	& + e3u(ibdy1:ibdy2,jj,jk,Krhs_a) &
114	& * uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk)
115	END DO
116	END DO
117	DO jj = 1, jpj
118	uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) * r1_hu_a(ibdy1:ibdy2,jj)
119	END DO
120	ENDIF
121	!
122	IF( .NOT.lk_agrif_clp ) THEN
123	DO jk=1,jpkm1 ! Smooth
124	DO jj=j1,j2
125	uu(ibdy2,jj,jk,Krhs_a) = 0.25_wp( uu(ibdy2-1,jj,jk,Krhs_a)+2._wpuu(ibdy2,jj,jk,Krhs_a) &
126	& + uu(ibdy2+1,jj,jk,Krhs_a) )
127	END DO
128	END DO
129	ENDIF
130	!
131	zub(ibdy1:ibdy2,:) = 0._wp ! Correct transport
132	DO jk = 1, jpkm1
133	DO jj = 1, jpj
134	zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) + e3u(ibdy1:ibdy2,jj,jk,Krhs_a) &
135	& * uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk)
136	END DO
137	END DO
138	DO jj=1,jpj
139	zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)
140	END DO
141
142	DO jk = 1, jpkm1
143	DO jj = 1, jpj
144	uu(ibdy1:ibdy2,jj,jk,Krhs_a) = ( uu(ibdy1:ibdy2,jj,jk,Krhs_a) &
145	& + uu_b(ibdy1:ibdy2,jj, Krhs_a) &
146	& - zub(ibdy1:ibdy2,jj) ) &
147	& * umask(ibdy1:ibdy2,jj,jk)
148	END DO
149	END DO
150
151	IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate
152	zvb(ibdy1:ibdy2,:) = 0._wp
153	DO jk = 1, jpkm1
154	DO jj = 1, jpj
155	zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) + e3v(ibdy1:ibdy2,jj,jk,Krhs_a) &
156	& * vv(ibdy1:ibdy2,jj,jk,Krhs_a) * vmask(ibdy1:ibdy2,jj,jk)
157	END DO
158	END DO
159	DO jj = 1, jpj
160	zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv_a(ibdy1:ibdy2,jj)
161	END DO
162	DO jk = 1, jpkm1
163	DO jj = 1, jpj
164	vv(ibdy1:ibdy2,jj,jk,Krhs_a) = ( vv(ibdy1:ibdy2,jj,jk,Krhs_a) &
165	& + vv_b(ibdy1:ibdy2,jj, Krhs_a) &
166	& - zvb(ibdy1:ibdy2,jj) ) &
167	& * vmask(ibdy1:ibdy2,jj,jk)
168	END DO
169	END DO
170	ENDIF
171	!
172	DO jk = 1, jpkm1 ! Mask domain edges
173	DO jj = 1, jpj
174	uu(1,jj,jk,Krhs_a) = 0._wp
175	vv(1,jj,jk,Krhs_a) = 0._wp
176	END DO
177	END DO
178	ENDIF
179
180	! --- East --- !
181	IF( nbondi == 1 .OR. nbondi == 2 ) THEN
182	ibdy1 = nlci-1-nbghostcells
183	ibdy2 = nlci-2
184	!
185	IF( .NOT.ln_dynspg_ts ) THEN ! Store transport
186	uu_b(ibdy1:ibdy2,:,Krhs_a) = 0._wp
187	DO jk = 1, jpkm1
188	DO jj = 1, jpj
189	uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj, Krhs_a) &
190	& + e3u(ibdy1:ibdy2,jj,jk,Krhs_a) &
191	& * uu(ibdy1:ibdy2,jj,jk,Krhs_a) &
192	& * umask(ibdy1:ibdy2,jj,jk)
193	END DO
194	END DO
195	DO jj = 1, jpj
196	uu_b(ibdy1:ibdy2,jj,Krhs_a) = uu_b(ibdy1:ibdy2,jj,Krhs_a) * r1_hu_a(ibdy1:ibdy2,jj)
197	END DO
198	ENDIF
199	!
200	IF( .NOT.lk_agrif_clp ) THEN
201	DO jk=1,jpkm1 ! Smooth
202	DO jj=j1,j2
203	uu(ibdy1,jj,jk,Krhs_a) = 0.25_wp*( uu(ibdy1-1,jj,jk,Krhs_a) &
204	& + 2._wp*uu(ibdy1 ,jj,jk,Krhs_a) &
205	& + uu(ibdy1+1,jj,jk,Krhs_a) )
206	END DO
207	END DO
208	ENDIF
209	!
210	zub(ibdy1:ibdy2,:) = 0._wp ! Correct transport
211	DO jk = 1, jpkm1
212	DO jj = 1, jpj
213	zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) &
214	& + e3u(ibdy1:ibdy2,jj,jk,Krhs_a) &
215	& * uu(ibdy1:ibdy2,jj,jk,Krhs_a) * umask(ibdy1:ibdy2,jj,jk)
216	END DO
217	END DO
218	DO jj=1,jpj
219	zub(ibdy1:ibdy2,jj) = zub(ibdy1:ibdy2,jj) * r1_hu_a(ibdy1:ibdy2,jj)
220	END DO
221
222	DO jk = 1, jpkm1
223	DO jj = 1, jpj
224	uu(ibdy1:ibdy2,jj,jk,Krhs_a) = ( uu(ibdy1:ibdy2,jj,jk,Krhs_a) &
225	& + uu_b(ibdy1:ibdy2,jj, Krhs_a) &
226	& - zub(ibdy1:ibdy2,jj) &
227	& ) * umask(ibdy1:ibdy2,jj,jk)
228	END DO
229	END DO
230
231	IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate
232	ibdy1 = ibdy1 + 1
233	ibdy2 = ibdy2 + 1
234	zvb(ibdy1:ibdy2,:) = 0._wp
235	DO jk = 1, jpkm1
236	DO jj = 1, jpj
237	zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) &
238	& + e3v(ibdy1:ibdy2,jj,jk,Krhs_a) &
239	& * vv(ibdy1:ibdy2,jj,jk,Krhs_a) &
240	& * vmask(ibdy1:ibdy2,jj,jk)
241	END DO
242	END DO
243	DO jj = 1, jpj
244	zvb(ibdy1:ibdy2,jj) = zvb(ibdy1:ibdy2,jj) * r1_hv_a(ibdy1:ibdy2,jj)
245	END DO
246	DO jk = 1, jpkm1
247	DO jj = 1, jpj
248	vv(ibdy1:ibdy2,jj,jk,Krhs_a) = ( vv(ibdy1:ibdy2,jj,jk,Krhs_a) &
249	& + vv_b(ibdy1:ibdy2,jj, Krhs_a) &
250	& - zvb(ibdy1:ibdy2,jj) &
251	& ) * vmask(ibdy1:ibdy2,jj,jk)
252	END DO
253	END DO
254	ENDIF
255	!
256	DO jk = 1, jpkm1 ! Mask domain edges
257	DO jj = 1, jpj
258	uu(nlci-1,jj,jk,Krhs_a) = 0._wp
259	vv(nlci ,jj,jk,Krhs_a) = 0._wp
260	END DO
261	END DO
262	ENDIF
263
264	! --- South --- !
265	IF( nbondj == -1 .OR. nbondj == 2 ) THEN
266	jbdy1 = 2
267	jbdy2 = 1+nbghostcells
268	!
269	IF( .NOT.ln_dynspg_ts ) THEN ! Store transport
270	vv_b(:,jbdy1:jbdy2,Krhs_a) = 0._wp
271	DO jk = 1, jpkm1
272	DO ji = 1, jpi
273	vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2, Krhs_a) &
274	& + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) &
275	& * vv(ji,jbdy1:jbdy2,jk,Krhs_a) &
276	& * vmask(ji,jbdy1:jbdy2,jk)
277	END DO
278	END DO
279	DO ji=1,jpi
280	vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) * r1_hv_a(ji,jbdy1:jbdy2)
281	END DO
282	ENDIF
283	!
284	IF ( .NOT.lk_agrif_clp ) THEN
285	DO jk = 1, jpkm1 ! Smooth
286	DO ji = i1, i2
287	vv(ji,jbdy2,jk,Krhs_a) = 0.25_wp*( vv(ji,jbdy2-1,jk,Krhs_a) &
288	& + 2._wp*vv(ji,jbdy2 ,jk,Krhs_a) &
289	& + vv(ji,jbdy2+1,jk,Krhs_a) )
290	END DO
291	END DO
292	ENDIF
293	!
294	zvb(:,jbdy1:jbdy2) = 0._wp ! Correct transport
295	DO jk=1,jpkm1
296	DO ji=1,jpi
297	zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) &
298	& + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) &
299	& * vv(ji,jbdy1:jbdy2,jk,Krhs_a) &
300	& * vmask(ji,jbdy1:jbdy2,jk)
301	END DO
302	END DO
303	DO ji = 1, jpi
304	zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)
305	END DO
306
307	DO jk = 1, jpkm1
308	DO ji = 1, jpi
309	vv(ji,jbdy1:jbdy2,jk,Krhs_a) = ( vv(ji,jbdy1:jbdy2,jk,Krhs_a) &
310	& + vv_b(ji,jbdy1:jbdy2, Krhs_a) &
311	& - zvb(ji,jbdy1:jbdy2) &
312	& ) * vmask(ji,jbdy1:jbdy2,jk)
313	END DO
314	END DO
315
316	IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate
317	zub(:,jbdy1:jbdy2) = 0._wp
318	DO jk = 1, jpkm1
319	DO ji = 1, jpi
320	zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) &
321	& + e3u(ji,jbdy1:jbdy2,jk,Krhs_a) &
322	& * uu(ji,jbdy1:jbdy2,jk,Krhs_a) &
323	& * umask(ji,jbdy1:jbdy2,jk)
324	END DO
325	END DO
326	DO ji = 1, jpi
327	zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu_a(ji,jbdy1:jbdy2)
328	END DO
329
330	DO jk = 1, jpkm1
331	DO ji = 1, jpi
332	uu(ji,jbdy1:jbdy2,jk,Krhs_a) = ( uu(ji,jbdy1:jbdy2,jk,Krhs_a) &
333	& + uu_b(ji,jbdy1:jbdy2, Krhs_a) &
334	& - zub(ji,jbdy1:jbdy2) &
335	& ) * umask(ji,jbdy1:jbdy2,jk)
336	END DO
337	END DO
338	ENDIF
339	!
340	DO jk = 1, jpkm1 ! Mask domain edges
341	DO ji = 1, jpi
342	uu(ji,1,jk,Krhs_a) = 0._wp
343	vv(ji,1,jk,Krhs_a) = 0._wp
344	END DO
345	END DO
346	ENDIF
347
348	! --- North --- !
349	IF( nbondj == 1 .OR. nbondj == 2 ) THEN
350	jbdy1 = nlcj-1-nbghostcells
351	jbdy2 = nlcj-2
352	!
353	IF( .NOT.ln_dynspg_ts ) THEN ! Store transport
354	vv_b(:,jbdy1:jbdy2,Krhs_a) = 0._wp
355	DO jk = 1, jpkm1
356	DO ji = 1, jpi
357	vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2, Krhs_a) &
358	& + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) &
359	& * vv(ji,jbdy1:jbdy2,jk,Krhs_a) &
360	& * vmask(ji,jbdy1:jbdy2,jk)
361	END DO
362	END DO
363	DO ji=1,jpi
364	vv_b(ji,jbdy1:jbdy2,Krhs_a) = vv_b(ji,jbdy1:jbdy2,Krhs_a) * r1_hv_a(ji,jbdy1:jbdy2)
365	END DO
366	ENDIF
367	!
368	IF ( .NOT.lk_agrif_clp ) THEN
369	DO jk = 1, jpkm1 ! Smooth
370	DO ji = i1, i2
371	vv(ji,jbdy1,jk,Krhs_a) = 0.25_wp*( vv(ji,jbdy1-1,jk,Krhs_a) &
372	& + 2._wp*vv(ji,jbdy1 ,jk,Krhs_a) &
373	& + vv(ji,jbdy1+1,jk,Krhs_a) )
374	END DO
375	END DO
376	ENDIF
377	!
378	zvb(:,jbdy1:jbdy2) = 0._wp ! Correct transport
379	DO jk=1,jpkm1
380	DO ji=1,jpi
381	zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) &
382	& + e3v(ji,jbdy1:jbdy2,jk,Krhs_a) &
383	& * vv(ji,jbdy1:jbdy2,jk,Krhs_a) &
384	& * vmask(ji,jbdy1:jbdy2,jk)
385	END DO
386	END DO
387	DO ji = 1, jpi
388	zvb(ji,jbdy1:jbdy2) = zvb(ji,jbdy1:jbdy2) * r1_hv_a(ji,jbdy1:jbdy2)
389	END DO
390
391	DO jk = 1, jpkm1
392	DO ji = 1, jpi
393	vv(ji,jbdy1:jbdy2,jk,Krhs_a) = ( vv(ji,jbdy1:jbdy2,jk,Krhs_a) &
394	& + vv_b(ji,jbdy1:jbdy2, Krhs_a) &
395	& - zvb(ji,jbdy1:jbdy2) &
396	& ) * vmask(ji,jbdy1:jbdy2,jk)
397	END DO
398	END DO
399
400	IF( ln_dynspg_ts ) THEN ! Set tangential velocities to time splitting estimate
401	jbdy1 = jbdy1 + 1
402	jbdy2 = jbdy2 + 1
403	zub(:,jbdy1:jbdy2) = 0._wp
404	DO jk = 1, jpkm1
405	DO ji = 1, jpi
406	zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) &
407	& + e3u(ji,jbdy1:jbdy2,jk,Krhs_a) &
408	& * uu(ji,jbdy1:jbdy2,jk,Krhs_a) &
409	& * umask(ji,jbdy1:jbdy2,jk)
410	END DO
411	END DO
412	DO ji = 1, jpi
413	zub(ji,jbdy1:jbdy2) = zub(ji,jbdy1:jbdy2) * r1_hu_a(ji,jbdy1:jbdy2)
414	END DO
415
416	DO jk = 1, jpkm1
417	DO ji = 1, jpi
418	uu(ji,jbdy1:jbdy2,jk,Krhs_a) = ( uu(ji,jbdy1:jbdy2,jk,Krhs_a) &
419	& + uu_b(ji,jbdy1:jbdy2, Krhs_a) &
420	& - zub(ji,jbdy1:jbdy2) &
421	& ) * umask(ji,jbdy1:jbdy2,jk)
422	END DO
423	END DO
424	ENDIF
425	!
426	DO jk = 1, jpkm1 ! Mask domain edges
427	DO ji = 1, jpi
428	uu(ji,nlcj ,jk,Krhs_a) = 0._wp
429	vv(ji,nlcj-1,jk,Krhs_a) = 0._wp
430	END DO
431	END DO
432	ENDIF
433	!
434	END SUBROUTINE Agrif_dyn
435
436
437	SUBROUTINE Agrif_dyn_ts( jn )
438	!!----------------------------------------------------------------------
439	!! * ROUTINE Agrif_dyn_ts *
440	!!----------------------------------------------------------------------
441	INTEGER, INTENT(in) :: jn
442	!!
443	INTEGER :: ji, jj
444	!!----------------------------------------------------------------------
445	!
446	IF( Agrif_Root() ) RETURN
447	!
448	IF((nbondi == -1).OR.(nbondi == 2)) THEN
449	DO jj=1,jpj
450	va_e(2:nbghostcells+1,jj) = vbdy_w(1:nbghostcells,jj) * hvr_e(2:nbghostcells+1,jj)
451	! Specified fluxes:
452	ua_e(2:nbghostcells+1,jj) = ubdy_w(1:nbghostcells,jj) * hur_e(2:nbghostcells+1,jj)
453	! Characteristics method (only if ghostcells=1):
454	!alt ua_e(2,jj) = 0.5_wp * ( ubdy_w(jj) * hur_e(2,jj) + ua_e(3,jj) &
455	!alt & - sqrt(grav * hur_e(2,jj)) * (sshn_e(3,jj) - hbdy_w(jj)) )
456	END DO
457	ENDIF
458	!
459	IF((nbondi == 1).OR.(nbondi == 2)) THEN
460	DO jj=1,jpj
461	va_e(nlci-nbghostcells:nlci-1,jj) = vbdy_e(1:nbghostcells,jj) * hvr_e(nlci-nbghostcells:nlci-1,jj)
462	! Specified fluxes:
463	ua_e(nlci-nbghostcells-1:nlci-2,jj) = ubdy_e(1:nbghostcells,jj) * hur_e(nlci-nbghostcells-1:nlci-2,jj)
464	! Characteristics method (only if ghostcells=1):
465	!alt ua_e(nlci-2,jj) = 0.5_wp * ( ubdy_e(jj) * hur_e(nlci-2,jj) + ua_e(nlci-3,jj) &
466	!alt & + sqrt(grav * hur_e(nlci-2,jj)) * (sshn_e(nlci-2,jj) - hbdy_e(jj)) )
467	END DO
468	ENDIF
469	!
470	IF((nbondj == -1).OR.(nbondj == 2)) THEN
471	DO ji=1,jpi
472	ua_e(ji,2:nbghostcells+1) = ubdy_s(ji,1:nbghostcells) * hur_e(ji,2:nbghostcells+1)
473	! Specified fluxes:
474	va_e(ji,2:nbghostcells+1) = vbdy_s(ji,1:nbghostcells) * hvr_e(ji,2:nbghostcells+1)
475	! Characteristics method (only if ghostcells=1):
476	!alt va_e(ji,2) = 0.5_wp * ( vbdy_s(ji) * hvr_e(ji,2) + va_e(ji,3) &
477	!alt & - sqrt(grav * hvr_e(ji,2)) * (sshn_e(ji,3) - hbdy_s(ji)) )
478	END DO
479	ENDIF
480	!
481	IF((nbondj == 1).OR.(nbondj == 2)) THEN
482	DO ji=1,jpi
483	ua_e(ji,nlcj-nbghostcells:nlcj-1) = ubdy_n(ji,1:nbghostcells) * hur_e(ji,nlcj-nbghostcells:nlcj-1)
484	! Specified fluxes:
485	va_e(ji,nlcj-nbghostcells-1:nlcj-2) = vbdy_n(ji,1:nbghostcells) * hvr_e(ji,nlcj-nbghostcells-1:nlcj-2)
486	! Characteristics method (only if ghostcells=1):
487	!alt va_e(ji,nlcj-2) = 0.5_wp * ( vbdy_n(ji) * hvr_e(ji,nlcj-2) + va_e(ji,nlcj-3) &
488	!alt & + sqrt(grav * hvr_e(ji,nlcj-2)) * (sshn_e(ji,nlcj-2) - hbdy_n(ji)) )
489	END DO
490	ENDIF
491	!
492	END SUBROUTINE Agrif_dyn_ts
493
494
495	SUBROUTINE Agrif_dta_ts( kt )
496	!!----------------------------------------------------------------------
497	!! * ROUTINE Agrif_dta_ts *
498	!!----------------------------------------------------------------------
499	INTEGER, INTENT(in) :: kt
500	!!
501	INTEGER :: ji, jj
502	LOGICAL :: ll_int_cons
503	!!----------------------------------------------------------------------
504	!
505	IF( Agrif_Root() ) RETURN
506	!
507	ll_int_cons = ln_bt_fw ! Assume conservative temporal integration in the forward case only
508	!
509	! Enforce volume conservation if no time refinement:
510	IF ( Agrif_rhot()==1 ) ll_int_cons=.TRUE.
511	!
512	! Interpolate barotropic fluxes
513	Agrif_SpecialValue=0._wp
514	Agrif_UseSpecialValue = ln_spc_dyn
515	!
516	IF( ll_int_cons ) THEN ! Conservative interpolation
517	! order matters here !!!!!!
518	CALL Agrif_Bc_variable( ub2b_interp_id, calledweight=1._wp, procname=interpub2b ) ! Time integrated
519	CALL Agrif_Bc_variable( vb2b_interp_id, calledweight=1._wp, procname=interpvb2b )
520	bdy_tinterp = 1
521	CALL Agrif_Bc_variable( unb_id , calledweight=1._wp, procname=interpunb ) ! After
522	CALL Agrif_Bc_variable( vnb_id , calledweight=1._wp, procname=interpvnb )
523	bdy_tinterp = 2
524	CALL Agrif_Bc_variable( unb_id , calledweight=0._wp, procname=interpunb ) ! Before
525	CALL Agrif_Bc_variable( vnb_id , calledweight=0._wp, procname=interpvnb )
526	ELSE ! Linear interpolation
527	bdy_tinterp = 0
528	ubdy_w(:,:) = 0._wp ; vbdy_w(:,:) = 0._wp
529	ubdy_e(:,:) = 0._wp ; vbdy_e(:,:) = 0._wp
530	ubdy_n(:,:) = 0._wp ; vbdy_n(:,:) = 0._wp
531	ubdy_s(:,:) = 0._wp ; vbdy_s(:,:) = 0._wp
532	CALL Agrif_Bc_variable( unb_id, procname=interpunb )
533	CALL Agrif_Bc_variable( vnb_id, procname=interpvnb )
534	ENDIF
535	Agrif_UseSpecialValue = .FALSE.
536	!
537	END SUBROUTINE Agrif_dta_ts
538
539
540	SUBROUTINE Agrif_ssh( kt )
541	!!----------------------------------------------------------------------
542	!! * ROUTINE Agrif_ssh *
543	!!----------------------------------------------------------------------
544	INTEGER, INTENT(in) :: kt
545	!
546	INTEGER :: ji, jj, indx, indy
547	!!----------------------------------------------------------------------
548	!
549	IF( Agrif_Root() ) RETURN
550	!
551	! Linear time interpolation of sea level
552	!
553	Agrif_SpecialValue = 0._wp
554	Agrif_UseSpecialValue = .TRUE.
555	CALL Agrif_Bc_variable(sshn_id, procname=interpsshn )
556	Agrif_UseSpecialValue = .FALSE.
557	!
558	! --- West --- !
559	IF((nbondi == -1).OR.(nbondi == 2)) THEN
560	indx = 1+nbghostcells
561	DO jj = 1, jpj
562	DO ji = 2, indx
563	ssh(ji,jj,Krhs_a) = hbdy_w(ji-1,jj)
564	ENDDO
565	ENDDO
566	ENDIF
567	!
568	! --- East --- !
569	IF((nbondi == 1).OR.(nbondi == 2)) THEN
570	indx = nlci-nbghostcells
571	DO jj = 1, jpj
572	DO ji = indx, nlci-1
573	ssh(ji,jj,Krhs_a) = hbdy_e(ji-indx+1,jj)
574	ENDDO
575	ENDDO
576	ENDIF
577	!
578	! --- South --- !
579	IF((nbondj == -1).OR.(nbondj == 2)) THEN
580	indy = 1+nbghostcells
581	DO jj = 2, indy
582	DO ji = 1, jpi
583	ssh(ji,jj,Krhs_a) = hbdy_s(ji,jj-1)
584	ENDDO
585	ENDDO
586	ENDIF
587	!
588	! --- North --- !
589	IF((nbondj == 1).OR.(nbondj == 2)) THEN
590	indy = nlcj-nbghostcells
591	DO jj = indy, nlcj-1
592	DO ji = 1, jpi
593	ssh(ji,jj,Krhs_a) = hbdy_n(ji,jj-indy+1)
594	ENDDO
595	ENDDO
596	ENDIF
597	!
598	END SUBROUTINE Agrif_ssh
599
600
601	SUBROUTINE Agrif_ssh_ts( jn )
602	!!----------------------------------------------------------------------
603	!! * ROUTINE Agrif_ssh_ts *
604	!!----------------------------------------------------------------------
605	INTEGER, INTENT(in) :: jn
606	!!
607	INTEGER :: ji, jj, indx, indy
608	!!----------------------------------------------------------------------
609	!! clem ghost (starting at i,j=1 is important I think otherwise you introduce a grad(ssh)/=0 at point 2)
610	!
611	IF( Agrif_Root() ) RETURN
612	!
613	! --- West --- !
614	IF((nbondi == -1).OR.(nbondi == 2)) THEN
615	indx = 1+nbghostcells
616	DO jj = 1, jpj
617	DO ji = 2, indx
618	ssha_e(ji,jj) = hbdy_w(ji-1,jj)
619	ENDDO
620	ENDDO
621	ENDIF
622	!
623	! --- East --- !
624	IF((nbondi == 1).OR.(nbondi == 2)) THEN
625	indx = nlci-nbghostcells
626	DO jj = 1, jpj
627	DO ji = indx, nlci-1
628	ssha_e(ji,jj) = hbdy_e(ji-indx+1,jj)
629	ENDDO
630	ENDDO
631	ENDIF
632	!
633	! --- South --- !
634	IF((nbondj == -1).OR.(nbondj == 2)) THEN
635	indy = 1+nbghostcells
636	DO jj = 2, indy
637	DO ji = 1, jpi
638	ssha_e(ji,jj) = hbdy_s(ji,jj-1)
639	ENDDO
640	ENDDO
641	ENDIF
642	!
643	! --- North --- !
644	IF((nbondj == 1).OR.(nbondj == 2)) THEN
645	indy = nlcj-nbghostcells
646	DO jj = indy, nlcj-1
647	DO ji = 1, jpi
648	ssha_e(ji,jj) = hbdy_n(ji,jj-indy+1)
649	ENDDO
650	ENDDO
651	ENDIF
652	!
653	END SUBROUTINE Agrif_ssh_ts
654
655	SUBROUTINE Agrif_avm
656	!!----------------------------------------------------------------------
657	!! * ROUTINE Agrif_avm *
658	!!----------------------------------------------------------------------
659	REAL(wp) :: zalpha
660	!!----------------------------------------------------------------------
661	!
662	IF( Agrif_Root() ) RETURN
663	!
664	zalpha = 1._wp ! JC: proper time interpolation impossible
665	! => use last available value from parent
666	!
667	Agrif_SpecialValue = 0.e0
668	Agrif_UseSpecialValue = .TRUE.
669	!
670	CALL Agrif_Bc_variable( avm_id, calledweight=zalpha, procname=interpavm )
671	!
672	Agrif_UseSpecialValue = .FALSE.
673	!
674	END SUBROUTINE Agrif_avm
675
676
677	SUBROUTINE interptsn( ptab, i1, i2, j1, j2, k1, k2, n1, n2, before, nb, ndir )
678	!!----------------------------------------------------------------------
679	!! * ROUTINE interptsn *
680	!!----------------------------------------------------------------------
681	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: ptab
682	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
683	LOGICAL , INTENT(in ) :: before
684	INTEGER , INTENT(in ) :: nb , ndir
685	!
686	INTEGER :: ji, jj, jk, jn, iref, jref, ibdy, jbdy ! dummy loop indices
687	INTEGER :: imin, imax, jmin, jmax, N_in, N_out
688	REAL(wp) :: zrho, z1, z2, z3, z4, z5, z6, z7
689	LOGICAL :: western_side, eastern_side,northern_side,southern_side
690	! vertical interpolation:
691	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk,n1:n2) :: ptab_child
692	REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
693	REAL(wp), DIMENSION(k1:k2) :: h_in
694	REAL(wp), DIMENSION(1:jpk) :: h_out
695	REAL(wp) :: h_diff
696
697	IF( before ) THEN
698	DO jn = 1,jpts
699	DO jk=k1,k2
700	DO jj=j1,j2
701	DO ji=i1,i2
702	ptab(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kmm_a)
703	END DO
704	END DO
705	END DO
706	END DO
707
708	# if defined key_vertical
709	DO jk=k1,k2
710	DO jj=j1,j2
711	DO ji=i1,i2
712	ptab(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm_a)
713	END DO
714	END DO
715	END DO
716	# endif
717	ELSE
718
719	western_side = (nb == 1).AND.(ndir == 1) ; eastern_side = (nb == 1).AND.(ndir == 2)
720	southern_side = (nb == 2).AND.(ndir == 1) ; northern_side = (nb == 2).AND.(ndir == 2)
721
722	# if defined key_vertical
723	DO jj=j1,j2
724	DO ji=i1,i2
725	iref = ji
726	jref = jj
727	if(western_side) iref=MAX(2 ,ji)
728	if(eastern_side) iref=MIN(nlci-1,ji)
729	if(southern_side) jref=MAX(2 ,jj)
730	if(northern_side) jref=MIN(nlcj-1,jj)
731	N_in = 0
732	DO jk=k1,k2 !k2 = jpk of parent grid
733	IF (ptab(ji,jj,jk,n2) == 0) EXIT
734	N_in = N_in + 1
735	tabin(jk,:) = ptab(ji,jj,jk,n1:n2-1)
736	h_in(N_in) = ptab(ji,jj,jk,n2)
737	END DO
738	N_out = 0
739	DO jk=1,jpk ! jpk of child grid
740	IF (tmask(iref,jref,jk) == 0) EXIT
741	N_out = N_out + 1
742	h_out(jk) = e3t(iref,jref,jk,Kmm_a)
743	ENDDO
744	IF (N_in > 0) THEN
745	DO jn=1,jpts
746	call reconstructandremap(tabin(1:N_in,jn),h_in,ptab_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
747	ENDDO
748	ENDIF
749	ENDDO
750	ENDDO
751	# else
752	ptab_child(i1:i2,j1:j2,1:jpk,1:jpts) = ptab(i1:i2,j1:j2,1:jpk,1:jpts)
753	# endif
754	!
755	DO jn=1, jpts
756	ts(i1:i2,j1:j2,1:jpk,jn,Krhs_a)=ptab_child(i1:i2,j1:j2,1:jpk,jn)*tmask(i1:i2,j1:j2,1:jpk)
757	END DO
758
759	IF ( .NOT.lk_agrif_clp ) THEN
760	!
761	imin = i1 ; imax = i2
762	jmin = j1 ; jmax = j2
763	!
764	! Remove CORNERS
765	IF((nbondj == -1).OR.(nbondj == 2)) jmin = 2 + nbghostcells
766	IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj - nbghostcells - 1
767	IF((nbondi == -1).OR.(nbondi == 2)) imin = 2 + nbghostcells
768	IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci - nbghostcells - 1
769	!
770	IF( eastern_side ) THEN
771	zrho = Agrif_Rhox()
772	z1 = ( zrho - 1._wp ) * 0.5_wp
773	z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )
774	z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
775	z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp )
776	z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
777	!
778	ibdy = nlci-nbghostcells
779	DO jn = 1, jpts
780	ts(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) &
781	& + z2 * ptab_child(ibdy ,jmin:jmax,1:jpkm1,jn)
782	DO jk = 1, jpkm1
783	DO jj = jmin,jmax
784	IF( umask(ibdy-1,jj,jk) == 0._wp ) THEN
785	ts(ibdy,jj,jk,jn,Krhs_a) = ts(ibdy+1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk)
786	ELSE
787	ts(ibdy,jj,jk,jn,Krhs_a) = ( z4 * ts(ibdy+1,jj,jk,jn,Krhs_a) &
788	& + z3 * ts(ibdy-1,jj,jk,jn,Krhs_a) &
789	& ) * tmask(ibdy ,jj,jk)
790	IF( uu(ibdy-1,jj,jk,Kmm_a) > 0._wp ) THEN
791	ts(ibdy,jj,jk,jn,Krhs_a) = ( z6 * ts(ibdy-1,jj,jk,jn,Krhs_a) &
792	& + z5 * ts(ibdy+1,jj,jk,jn,Krhs_a) &
793	& + z7 * ts(ibdy-2,jj,jk,jn,Krhs_a) &
794	& ) * tmask(ibdy ,jj,jk)
795	ENDIF
796	ENDIF
797	END DO
798	END DO
799	! Restore ghost points:
800	ts(ibdy+1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy+1,jmin:jmax,1:jpkm1,jn) &
801	& * tmask(ibdy+1,jmin:jmax,1:jpkm1)
802	END DO
803	ENDIF
804	!
805	IF( northern_side ) THEN
806	zrho = Agrif_Rhoy()
807	z1 = ( zrho - 1._wp ) * 0.5_wp
808	z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )
809	z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
810	z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp )
811	z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
812	!
813	jbdy = nlcj-nbghostcells
814	DO jn = 1, jpts
815	ts(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) &
816	& + z2 * ptab_child(imin:imax,jbdy ,1:jpkm1,jn)
817	DO jk = 1, jpkm1
818	DO ji = imin,imax
819	IF( vmask(ji,jbdy-1,jk) == 0._wp ) THEN
820	ts(ji,jbdy,jk,jn,Krhs_a) = ts(ji,jbdy+1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)
821	ELSE
822	ts(ji,jbdy,jk,jn,Krhs_a)=( z4 * ts(ji,jbdy+1,jk,jn,Krhs_a) &
823	& + z3 * ts(ji,jbdy-1,jk,jn,Krhs_a) &
824	& ) * tmask(ji,jbdy ,jk)
825	IF (vv(ji,jbdy-1,jk,Kmm_a) > 0._wp ) THEN
826	ts(ji,jbdy,jk,jn,Krhs_a)=( z6 * ts(ji,jbdy-1,jk,jn,Krhs_a) &
827	& + z5 * ts(ji,jbdy+1,jk,jn,Krhs_a) &
828	& + z7 * ts(ji,jbdy-2,jk,jn,Krhs_a) &
829	& ) * tmask(ji,jbdy ,jk)
830	ENDIF
831	ENDIF
832	END DO
833	END DO
834	! Restore ghost points:
835	ts(imin:imax,jbdy+1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy+1,1:jpkm1,jn) &
836	& * tmask(imin:imax,jbdy+1,1:jpkm1)
837	END DO
838	ENDIF
839	!
840	IF( western_side ) THEN
841	zrho = Agrif_Rhox()
842	z1 = ( zrho - 1._wp ) * 0.5_wp
843	z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )
844	z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
845	z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp )
846	z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
847	!
848	ibdy = 1+nbghostcells
849	DO jn = 1, jpts
850	ts(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) &
851	& + z2 * ptab_child(ibdy ,jmin:jmax,1:jpkm1,jn)
852	DO jk = 1, jpkm1
853	DO jj = jmin,jmax
854	IF( umask(ibdy,jj,jk) == 0._wp ) THEN
855	ts(ibdy,jj,jk,jn,Krhs_a) = ts(ibdy-1,jj,jk,jn,Krhs_a) * tmask(ibdy,jj,jk)
856	ELSE
857	ts(ibdy,jj,jk,jn,Krhs_a) = ( z4 * ts(ibdy-1,jj,jk,jn,Krhs_a) &
858	& + z3 * ts(ibdy+1,jj,jk,jn,Krhs_a) &
859	& ) * tmask(ibdy ,jj,jk)
860	IF( uu(ibdy,jj,jk,Kmm_a) < 0._wp ) THEN
861	ts(ibdy,jj,jk,jn,Krhs_a) = ( z6 * ts(ibdy+1,jj,jk,jn,Krhs_a) &
862	& + z5 * ts(ibdy-1,jj,jk,jn,Krhs_a) &
863	& + z7 * ts(ibdy+2,jj,jk,jn,Krhs_a) &
864	& ) * tmask(ibdy,jj,jk)
865	ENDIF
866	ENDIF
867	END DO
868	END DO
869	! Restore ghost points:
870	ts(ibdy-1,jmin:jmax,1:jpkm1,jn,Krhs_a) = ptab_child(ibdy-1,jmin:jmax,1:jpkm1,jn) &
871	& * tmask(ibdy-1,jmin:jmax,1:jpkm1)
872	END DO
873	ENDIF
874	!
875	IF( southern_side ) THEN
876	zrho = Agrif_Rhoy()
877	z1 = ( zrho - 1._wp ) * 0.5_wp
878	z3 = ( zrho - 1._wp ) / ( zrho + 1._wp )
879	z6 = 2._wp * ( zrho - 1._wp ) / ( zrho + 1._wp )
880	z7 = - ( zrho - 1._wp ) / ( zrho + 3._wp )
881	z2 = 1._wp - z1 ; z4 = 1._wp - z3 ; z5 = 1._wp - z6 - z7
882	!
883	jbdy=1+nbghostcells
884	DO jn = 1, jpts
885	ts(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = z1 * ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) &
886	& + z2 * ptab_child(imin:imax,jbdy ,1:jpkm1,jn)
887	DO jk = 1, jpkm1
888	DO ji = imin,imax
889	IF( vmask(ji,jbdy,jk) == 0._wp ) THEN
890	ts(ji,jbdy,jk,jn,Krhs_a) = ts(ji,jbdy-1,jk,jn,Krhs_a) * tmask(ji,jbdy,jk)
891	ELSE
892	ts(ji,jbdy,jk,jn,Krhs_a) = ( z4 * ts(ji,jbdy-1,jk,jn,Krhs_a) &
893	& + z3 * ts(ji,jbdy+1,jk,jn,Krhs_a) &
894	& ) * tmask(ji,jbdy ,jk)
895	IF( vv(ji,jbdy,jk,Kmm_a) < 0._wp ) THEN
896	ts(ji,jbdy,jk,jn,Krhs_a) = ( z6 * ts(ji,jbdy+1,jk,jn,Krhs_a) &
897	& + z5 * ts(ji,jbdy-1,jk,jn,Krhs_a) &
898	& + z7 * ts(ji,jbdy+2,jk,jn,Krhs_a) &
899	& ) * tmask(ji,jbdy ,jk)
900	ENDIF
901	ENDIF
902	END DO
903	END DO
904	! Restore ghost points:
905	ts(imin:imax,jbdy-1,1:jpkm1,jn,Krhs_a) = ptab_child(imin:imax,jbdy-1,1:jpkm1,jn) &
906	& * tmask(imin:imax,jbdy-1,1:jpkm1)
907	END DO
908	ENDIF
909	!
910	ENDIF
911	ENDIF
912	!
913	END SUBROUTINE interptsn
914
915	SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir )
916	!!----------------------------------------------------------------------
917	!! * ROUTINE interpsshn *
918	!!----------------------------------------------------------------------
919	INTEGER , INTENT(in ) :: i1, i2, j1, j2
920	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
921	LOGICAL , INTENT(in ) :: before
922	INTEGER , INTENT(in ) :: nb , ndir
923	!
924	LOGICAL :: western_side, eastern_side,northern_side,southern_side
925	!!----------------------------------------------------------------------
926	!
927	IF( before) THEN
928	ptab(i1:i2,j1:j2) = ssh(i1:i2,j1:j2,Kmm_a)
929	ELSE
930	western_side = (nb == 1).AND.(ndir == 1)
931	eastern_side = (nb == 1).AND.(ndir == 2)
932	southern_side = (nb == 2).AND.(ndir == 1)
933	northern_side = (nb == 2).AND.(ndir == 2)
934	!! clem ghost
935	IF(western_side) hbdy_w(1:nbghostcells,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
936	IF(eastern_side) hbdy_e(1:nbghostcells,j1:j2) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
937	IF(southern_side) hbdy_s(i1:i2,1:nbghostcells) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
938	IF(northern_side) hbdy_n(i1:i2,1:nbghostcells) = ptab(i1:i2,j1:j2) * tmask(i1:i2,j1:j2,1)
939	ENDIF
940	!
941	END SUBROUTINE interpsshn
942
943	SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
944	!!----------------------------------------------------------------------
945	!! * ROUTINE interpun *
946	!!---------------------------------------------
947	!!
948	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
949	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
950	LOGICAL, INTENT(in) :: before
951	INTEGER, INTENT(in) :: nb , ndir
952	!!
953	INTEGER :: ji,jj,jk
954	REAL(wp) :: zrhoy
955	! vertical interpolation:
956	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
957	REAL(wp), DIMENSION(1:jpk) :: h_out
958	INTEGER :: N_in, N_out, iref
959	REAL(wp) :: h_diff
960	LOGICAL :: western_side, eastern_side
961	!!---------------------------------------------
962	!
963	IF (before) THEN
964	DO jk=1,jpk
965	DO jj=j1,j2
966	DO ji=i1,i2
967	ptab(ji,jj,jk,1) = ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) &
968	& * uu(ji,jj,jk,Kmm_a) * umask(ji,jj,jk) )
969	# if defined key_vertical
970	ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a))
971	# endif
972	END DO
973	END DO
974	END DO
975	ELSE
976	zrhoy = Agrif_rhoy()
977	# if defined key_vertical
978	! VERTICAL REFINEMENT BEGIN
979	western_side = (nb == 1).AND.(ndir == 1)
980	eastern_side = (nb == 1).AND.(ndir == 2)
981
982	DO ji=i1,i2
983	iref = ji
984	IF (western_side) iref = MAX(2,ji)
985	IF (eastern_side) iref = MIN(nlci-2,ji)
986	DO jj=j1,j2
987	N_in = 0
988	DO jk=k1,k2
989	IF (ptab(ji,jj,jk,2) == 0) EXIT
990	N_in = N_in + 1
991	tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
992	h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy)
993	ENDDO
994
995	IF (N_in == 0) THEN
996	uu(ji,jj,:,Krhs_a) = 0._wp
997	CYCLE
998	ENDIF
999
1000	N_out = 0
1001	DO jk=1,jpk
1002	if (umask(iref,jj,jk) == 0) EXIT
1003	N_out = N_out + 1
1004	h_out(N_out) = e3u(iref,jj,jk,Krhs_a)
1005	ENDDO
1006
1007	IF (N_out == 0) THEN
1008	uu(ji,jj,:,Krhs_a) = 0._wp
1009	CYCLE
1010	ENDIF
1011
1012	IF (N_in * N_out > 0) THEN
1013	h_diff = SUM( h_out(1:N_out) ) - SUM( h_in(1:N_in) )
1014	! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
1015	if (h_diff < -1.e4) then
1016	print *,'CHECK YOUR BATHY ...', h_diff, SUM( h_out(1:N_out) ), SUM( h_in(1:N_in) )
1017	! stop
1018	endif
1019	ENDIF
1020	call reconstructandremap( tabin(1:N_in) , h_in(1:N_in), uu(ji,jj,1:N_out,Krhs_a), &
1021	& h_out(1:N_out), N_in , N_out )
1022	ENDDO
1023	ENDDO
1024
1025	# else
1026	DO jk = 1, jpkm1
1027	DO jj=j1,j2
1028	uu(i1:i2,jj,jk,Krhs_a) = ptab(i1:i2,jj,jk,1) / ( zrhoy * e2u(i1:i2,jj) * e3u(i1:i2,jj,jk,Krhs_a) )
1029	END DO
1030	END DO
1031	# endif
1032
1033	ENDIF
1034	!
1035	END SUBROUTINE interpun
1036
1037	SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
1038	!!----------------------------------------------------------------------
1039	!! * ROUTINE interpvn *
1040	!!----------------------------------------------------------------------
1041	!
1042	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
1043	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
1044	LOGICAL, INTENT(in) :: before
1045	INTEGER, INTENT(in) :: nb , ndir
1046	!
1047	INTEGER :: ji,jj,jk
1048	REAL(wp) :: zrhox
1049	! vertical interpolation:
1050	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
1051	REAL(wp), DIMENSION(1:jpk) :: h_out
1052	INTEGER :: N_in, N_out, jref
1053	REAL(wp) :: h_diff
1054	LOGICAL :: northern_side,southern_side
1055	!!---------------------------------------------
1056	!
1057	IF (before) THEN
1058	DO jk=k1,k2
1059	DO jj=j1,j2
1060	DO ji=i1,i2
1061	ptab(ji,jj,jk,1) = ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) &
1062	& * vv(ji,jj,jk,Kmm_a) * vmask(ji,jj,jk) )
1063	# if defined key_vertical
1064	ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
1065	# endif
1066	END DO
1067	END DO
1068	END DO
1069	ELSE
1070	zrhox = Agrif_rhox()
1071	# if defined key_vertical
1072
1073	southern_side = (nb == 2).AND.(ndir == 1)
1074	northern_side = (nb == 2).AND.(ndir == 2)
1075
1076	DO jj=j1,j2
1077	jref = jj
1078	IF (southern_side) jref = MAX(2,jj)
1079	IF (northern_side) jref = MIN(nlcj-2,jj)
1080	DO ji=i1,i2
1081	N_in = 0
1082	DO jk=k1,k2
1083	if (ptab(ji,jj,jk,2) == 0) EXIT
1084	N_in = N_in + 1
1085	tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
1086	h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)
1087	END DO
1088	IF (N_in == 0) THEN
1089	vv(ji,jj,:,Krhs_a) = 0._wp
1090	CYCLE
1091	ENDIF
1092
1093	N_out = 0
1094	DO jk=1,jpk
1095	if (vmask(ji,jref,jk) == 0) EXIT
1096	N_out = N_out + 1
1097	h_out(N_out) = e3v(ji,jref,jk,Krhs_a)
1098	END DO
1099	IF (N_out == 0) THEN
1100	vv(ji,jj,:,Krhs_a) = 0._wp
1101	CYCLE
1102	ENDIF
1103	call reconstructandremap( tabin(1:N_in) , h_in(1:N_in), vv(ji,jj,1:N_out,Krhs_a), &
1104	& h_out(1:N_out), N_in , N_out )
1105	END DO
1106	END DO
1107	# else
1108	DO jk = 1, jpkm1
1109	vv(i1:i2,j1:j2,jk,Krhs_a) = ptab(i1:i2,j1:j2,jk,1) / ( zrhox * e1v(i1:i2,j1:j2) * e3v(i1:i2,j1:j2,jk,Krhs_a) )
1110	END DO
1111	# endif
1112	ENDIF
1113	!
1114	END SUBROUTINE interpvn
1115
1116	SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir )
1117	!!----------------------------------------------------------------------
1118	!! * ROUTINE interpunb *
1119	!!----------------------------------------------------------------------
1120	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1121	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1122	LOGICAL , INTENT(in ) :: before
1123	INTEGER , INTENT(in ) :: nb , ndir
1124	!
1125	INTEGER :: ji, jj
1126	REAL(wp) :: zrhoy, zrhot, zt0, zt1, ztcoeff
1127	LOGICAL :: western_side, eastern_side,northern_side,southern_side
1128	!!----------------------------------------------------------------------
1129	!
1130	IF( before ) THEN
1131	ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * uu_b(i1:i2,j1:j2,Kmm_a)
1132	ELSE
1133	western_side = (nb == 1).AND.(ndir == 1)
1134	eastern_side = (nb == 1).AND.(ndir == 2)
1135	southern_side = (nb == 2).AND.(ndir == 1)
1136	northern_side = (nb == 2).AND.(ndir == 2)
1137	zrhoy = Agrif_Rhoy()
1138	zrhot = Agrif_rhot()
1139	! Time indexes bounds for integration
1140	zt0 = REAL(Agrif_NbStepint() , wp) / zrhot
1141	zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1142	! Polynomial interpolation coefficients:
1143	IF( bdy_tinterp == 1 ) THEN
1144	ztcoeff = zrhot * ( zt1*2._wp ( zt1 - 1._wp) &
1145	& - zt0*2._wp ( zt0 - 1._wp) )
1146	ELSEIF( bdy_tinterp == 2 ) THEN
1147	ztcoeff = zrhot * ( zt1 * ( zt1 - 1._wp)**2._wp &
1148	& - zt0 * ( zt0 - 1._wp)**2._wp )
1149	ELSE
1150	ztcoeff = 1
1151	ENDIF
1152	!
1153	IF(western_side) ubdy_w(1:nbghostcells,j1:j2) = ubdy_w(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2)
1154	IF(eastern_side) ubdy_e(1:nbghostcells,j1:j2) = ubdy_e(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2)
1155	IF(southern_side) ubdy_s(i1:i2,1:nbghostcells) = ubdy_s(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2)
1156	IF(northern_side) ubdy_n(i1:i2,1:nbghostcells) = ubdy_n(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2)
1157	!
1158	IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1159	IF(western_side) ubdy_w(1:nbghostcells,j1:j2) = ubdy_w(1:nbghostcells,j1:j2) / (zrhoye2u(i1:i2,j1:j2)) umask(i1:i2,j1:j2,1)
1160	IF(eastern_side) ubdy_e(1:nbghostcells,j1:j2) = ubdy_e(1:nbghostcells,j1:j2) / (zrhoye2u(i1:i2,j1:j2)) umask(i1:i2,j1:j2,1)
1161	IF(southern_side) ubdy_s(i1:i2,1:nbghostcells) = ubdy_s(i1:i2,1:nbghostcells) / (zrhoye2u(i1:i2,j1:j2)) umask(i1:i2,j1:j2,1)
1162	IF(northern_side) ubdy_n(i1:i2,1:nbghostcells) = ubdy_n(i1:i2,1:nbghostcells) / (zrhoye2u(i1:i2,j1:j2)) umask(i1:i2,j1:j2,1)
1163	ENDIF
1164	ENDIF
1165	!
1166	END SUBROUTINE interpunb
1167
1168
1169	SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before, nb, ndir )
1170	!!----------------------------------------------------------------------
1171	!! * ROUTINE interpvnb *
1172	!!----------------------------------------------------------------------
1173	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1174	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1175	LOGICAL , INTENT(in ) :: before
1176	INTEGER , INTENT(in ) :: nb , ndir
1177	!
1178	INTEGER :: ji,jj
1179	REAL(wp) :: zrhox, zrhot, zt0, zt1, ztcoeff
1180	LOGICAL :: western_side, eastern_side,northern_side,southern_side
1181	!!----------------------------------------------------------------------
1182	!
1183	IF( before ) THEN
1184	ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vv_b(i1:i2,j1:j2,Kmm_a)
1185	ELSE
1186	western_side = (nb == 1).AND.(ndir == 1)
1187	eastern_side = (nb == 1).AND.(ndir == 2)
1188	southern_side = (nb == 2).AND.(ndir == 1)
1189	northern_side = (nb == 2).AND.(ndir == 2)
1190	zrhox = Agrif_Rhox()
1191	zrhot = Agrif_rhot()
1192	! Time indexes bounds for integration
1193	zt0 = REAL(Agrif_NbStepint() , wp) / zrhot
1194	zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1195	IF( bdy_tinterp == 1 ) THEN
1196	ztcoeff = zrhot * ( zt1*2._wp ( zt1 - 1._wp) &
1197	& - zt0*2._wp ( zt0 - 1._wp) )
1198	ELSEIF( bdy_tinterp == 2 ) THEN
1199	ztcoeff = zrhot * ( zt1 * ( zt1 - 1._wp)**2._wp &
1200	& - zt0 * ( zt0 - 1._wp)**2._wp )
1201	ELSE
1202	ztcoeff = 1
1203	ENDIF
1204	!! clem ghost
1205	IF(western_side) vbdy_w(1:nbghostcells,j1:j2) = vbdy_w(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2)
1206	IF(eastern_side) vbdy_e(1:nbghostcells,j1:j2) = vbdy_e(1:nbghostcells,j1:j2) + ztcoeff * ptab(i1:i2,j1:j2)
1207	IF(southern_side) vbdy_s(i1:i2,1:nbghostcells) = vbdy_s(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2)
1208	IF(northern_side) vbdy_n(i1:i2,1:nbghostcells) = vbdy_n(i1:i2,1:nbghostcells) + ztcoeff * ptab(i1:i2,j1:j2)
1209	!
1210	IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1211	IF(western_side) vbdy_w(1:nbghostcells,j1:j2) = vbdy_w(1:nbghostcells,j1:j2) / (zrhoxe1v(i1:i2,j1:j2)) vmask(i1:i2,j1:j2,1)
1212	IF(eastern_side) vbdy_e(1:nbghostcells,j1:j2) = vbdy_e(1:nbghostcells,j1:j2) / (zrhoxe1v(i1:i2,j1:j2)) vmask(i1:i2,j1:j2,1)
1213	IF(southern_side) vbdy_s(i1:i2,1:nbghostcells) = vbdy_s(i1:i2,1:nbghostcells) / (zrhoxe1v(i1:i2,j1:j2)) vmask(i1:i2,j1:j2,1)
1214	IF(northern_side) vbdy_n(i1:i2,1:nbghostcells) = vbdy_n(i1:i2,1:nbghostcells) / (zrhoxe1v(i1:i2,j1:j2)) vmask(i1:i2,j1:j2,1)
1215	ENDIF
1216	ENDIF
1217	!
1218	END SUBROUTINE interpvnb
1219
1220
1221	SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before, nb, ndir )
1222	!!----------------------------------------------------------------------
1223	!! * ROUTINE interpub2b *
1224	!!----------------------------------------------------------------------
1225	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1226	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1227	LOGICAL , INTENT(in ) :: before
1228	INTEGER , INTENT(in ) :: nb , ndir
1229	!
1230	INTEGER :: ji,jj
1231	REAL(wp) :: zrhot, zt0, zt1,zat
1232	LOGICAL :: western_side, eastern_side,northern_side,southern_side
1233	!!----------------------------------------------------------------------
1234	IF( before ) THEN
1235	IF ( ln_bt_fw ) THEN
1236	ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1237	ELSE
1238	ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * un_adv(i1:i2,j1:j2)
1239	ENDIF
1240	ELSE
1241	western_side = (nb == 1).AND.(ndir == 1)
1242	eastern_side = (nb == 1).AND.(ndir == 2)
1243	southern_side = (nb == 2).AND.(ndir == 1)
1244	northern_side = (nb == 2).AND.(ndir == 2)
1245	zrhot = Agrif_rhot()
1246	! Time indexes bounds for integration
1247	zt0 = REAL(Agrif_NbStepint() , wp) / zrhot
1248	zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1249	! Polynomial interpolation coefficients:
1250	zat = zrhot * ( zt1*2._wp (-2._wp*zt1 + 3._wp) &
1251	& - zt0*2._wp (-2._wp*zt0 + 3._wp) )
1252	!! clem ghost
1253	IF(western_side ) ubdy_w(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2)
1254	IF(eastern_side ) ubdy_e(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2)
1255	IF(southern_side) ubdy_s(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2)
1256	IF(northern_side) ubdy_n(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2)
1257	ENDIF
1258	!
1259	END SUBROUTINE interpub2b
1260
1261
1262	SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before, nb, ndir )
1263	!!----------------------------------------------------------------------
1264	!! * ROUTINE interpvb2b *
1265	!!----------------------------------------------------------------------
1266	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1267	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1268	LOGICAL , INTENT(in ) :: before
1269	INTEGER , INTENT(in ) :: nb , ndir
1270	!
1271	INTEGER :: ji,jj
1272	REAL(wp) :: zrhot, zt0, zt1,zat
1273	LOGICAL :: western_side, eastern_side,northern_side,southern_side
1274	!!----------------------------------------------------------------------
1275	!
1276	IF( before ) THEN
1277	IF ( ln_bt_fw ) THEN
1278	ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1279	ELSE
1280	ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vn_adv(i1:i2,j1:j2)
1281	ENDIF
1282	ELSE
1283	western_side = (nb == 1).AND.(ndir == 1)
1284	eastern_side = (nb == 1).AND.(ndir == 2)
1285	southern_side = (nb == 2).AND.(ndir == 1)
1286	northern_side = (nb == 2).AND.(ndir == 2)
1287	zrhot = Agrif_rhot()
1288	! Time indexes bounds for integration
1289	zt0 = REAL(Agrif_NbStepint() , wp) / zrhot
1290	zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1291	! Polynomial interpolation coefficients:
1292	zat = zrhot * ( zt1*2._wp (-2._wp*zt1 + 3._wp) &
1293	& - zt0*2._wp (-2._wp*zt0 + 3._wp) )
1294	!
1295	IF(western_side ) vbdy_w(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2)
1296	IF(eastern_side ) vbdy_e(1:nbghostcells,j1:j2) = zat * ptab(i1:i2,j1:j2)
1297	IF(southern_side) vbdy_s(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2)
1298	IF(northern_side) vbdy_n(i1:i2,1:nbghostcells) = zat * ptab(i1:i2,j1:j2)
1299	ENDIF
1300	!
1301	END SUBROUTINE interpvb2b
1302
1303
1304	SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1305	!!----------------------------------------------------------------------
1306	!! * ROUTINE interpe3t *
1307	!!----------------------------------------------------------------------
1308	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1309	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1310	LOGICAL , INTENT(in ) :: before
1311	INTEGER , INTENT(in ) :: nb , ndir
1312	!
1313	INTEGER :: ji, jj, jk
1314	LOGICAL :: western_side, eastern_side, northern_side, southern_side
1315	!!----------------------------------------------------------------------
1316	!
1317	IF( before ) THEN
1318	ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
1319	ELSE
1320	western_side = (nb == 1).AND.(ndir == 1)
1321	eastern_side = (nb == 1).AND.(ndir == 2)
1322	southern_side = (nb == 2).AND.(ndir == 1)
1323	northern_side = (nb == 2).AND.(ndir == 2)
1324	!
1325	DO jk = k1, k2
1326	DO jj = j1, j2
1327	DO ji = i1, i2
1328	!
1329	IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) ) > 1.D-2) THEN
1330	IF (western_side.AND.(ptab(i1+nbghostcells-1,jj,jk)>0._wp)) THEN
1331	WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1332	WRITE(numout,*) ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1333	kindic_agr = kindic_agr + 1
1334	ELSEIF (eastern_side.AND.(ptab(i2-nbghostcells+1,jj,jk)>0._wp)) THEN
1335	WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1336	WRITE(numout,*) ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1337	kindic_agr = kindic_agr + 1
1338	ELSEIF (southern_side.AND.(ptab(ji,j1+nbghostcells-1,jk)>0._wp)) THEN
1339	WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1340	WRITE(numout,*) ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1341	kindic_agr = kindic_agr + 1
1342	ELSEIF (northern_side.AND.(ptab(ji,j2-nbghostcells+1,jk)>0._wp)) THEN
1343	WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1344	WRITE(numout,*) ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1345	kindic_agr = kindic_agr + 1
1346	ENDIF
1347	ENDIF
1348	END DO
1349	END DO
1350	END DO
1351	!
1352	ENDIF
1353	!
1354	END SUBROUTINE interpe3t
1355
1356
1357	SUBROUTINE interpumsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1358	!!----------------------------------------------------------------------
1359	!! * ROUTINE interpumsk *
1360	!!----------------------------------------------------------------------
1361	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1362	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1363	LOGICAL , INTENT(in ) :: before
1364	INTEGER , INTENT(in ) :: nb , ndir
1365	!
1366	INTEGER :: ji, jj, jk
1367	LOGICAL :: western_side, eastern_side
1368	!!----------------------------------------------------------------------
1369	!
1370	IF( before ) THEN
1371	ptab(i1:i2,j1:j2,k1:k2) = umask(i1:i2,j1:j2,k1:k2)
1372	ELSE
1373	western_side = (nb == 1).AND.(ndir == 1)
1374	eastern_side = (nb == 1).AND.(ndir == 2)
1375	DO jk = k1, k2
1376	DO jj = j1, j2
1377	DO ji = i1, i2
1378	! Velocity mask at boundary edge points:
1379	IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN
1380	IF (western_side) THEN
1381	WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1382	WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1383	kindic_agr = kindic_agr + 1
1384	ELSEIF (eastern_side) THEN
1385	WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1386	WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1387	kindic_agr = kindic_agr + 1
1388	ENDIF
1389	ENDIF
1390	END DO
1391	END DO
1392	END DO
1393	!
1394	ENDIF
1395	!
1396	END SUBROUTINE interpumsk
1397
1398
1399	SUBROUTINE interpvmsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1400	!!----------------------------------------------------------------------
1401	!! * ROUTINE interpvmsk *
1402	!!----------------------------------------------------------------------
1403	INTEGER , INTENT(in ) :: i1,i2,j1,j2,k1,k2
1404	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1405	LOGICAL , INTENT(in ) :: before
1406	INTEGER , INTENT(in ) :: nb , ndir
1407	!
1408	INTEGER :: ji, jj, jk
1409	LOGICAL :: northern_side, southern_side
1410	!!----------------------------------------------------------------------
1411	!
1412	IF( before ) THEN
1413	ptab(i1:i2,j1:j2,k1:k2) = vmask(i1:i2,j1:j2,k1:k2)
1414	ELSE
1415	southern_side = (nb == 2).AND.(ndir == 1)
1416	northern_side = (nb == 2).AND.(ndir == 2)
1417	DO jk = k1, k2
1418	DO jj = j1, j2
1419	DO ji = i1, i2
1420	! Velocity mask at boundary edge points:
1421	IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN
1422	IF (southern_side) THEN
1423	WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1424	WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1425	kindic_agr = kindic_agr + 1
1426	ELSEIF (northern_side) THEN
1427	WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1428	WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1429	kindic_agr = kindic_agr + 1
1430	ENDIF
1431	ENDIF
1432	END DO
1433	END DO
1434	END DO
1435	!
1436	ENDIF
1437	!
1438	END SUBROUTINE interpvmsk
1439
1440
1441	SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before )
1442	!!----------------------------------------------------------------------
1443	!! * ROUTINE interavm *
1444	!!----------------------------------------------------------------------
1445	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, m1, m2
1446	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
1447	LOGICAL , INTENT(in ) :: before
1448	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
1449	REAL(wp), DIMENSION(1:jpk) :: h_out
1450	INTEGER :: N_in, N_out, ji, jj, jk
1451	!!----------------------------------------------------------------------
1452	!
1453	IF (before) THEN
1454	DO jk=k1,k2
1455	DO jj=j1,j2
1456	DO ji=i1,i2
1457	ptab(ji,jj,jk,1) = avm_k(ji,jj,jk)
1458	END DO
1459	END DO
1460	END DO
1461	#ifdef key_vertical
1462	DO jk=k1,k2
1463	DO jj=j1,j2
1464	DO ji=i1,i2
1465	ptab(ji,jj,jk,2) = wmask(ji,jj,jk) * e3w(ji,jj,jk,Kmm_a)
1466	END DO
1467	END DO
1468	END DO
1469	#endif
1470	ELSE
1471	#ifdef key_vertical
1472	avm_k(i1:i2,j1:j2,1:jpk) = 0.
1473	DO jj=j1,j2
1474	DO ji=i1,i2
1475	N_in = 0
1476	DO jk=k1,k2 !k2 = jpk of parent grid
1477	IF (ptab(ji,jj,jk,2) == 0) EXIT
1478	N_in = N_in + 1
1479	tabin(jk) = ptab(ji,jj,jk,1)
1480	h_in(N_in) = ptab(ji,jj,jk,2)
1481	END DO
1482	N_out = 0
1483	DO jk=1,jpk ! jpk of child grid
1484	IF (wmask(ji,jj,jk) == 0) EXIT
1485	N_out = N_out + 1
1486	h_out(jk) = e3t(ji,jj,jk,Kmm_a)
1487	ENDDO
1488	IF (N_in > 0) THEN
1489	CALL reconstructandremap(tabin(1:N_in),h_in,avm_k(ji,jj,1:N_out),h_out,N_in,N_out)
1490	ENDIF
1491	ENDDO
1492	ENDDO
1493	#else
1494	avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2,1)
1495	#endif
1496	ENDIF
1497	!
1498	END SUBROUTINE interpavm
1499
1500	#else
1501	!!----------------------------------------------------------------------
1502	!! Empty module no AGRIF zoom
1503	!!----------------------------------------------------------------------
1504	CONTAINS
1505	SUBROUTINE Agrif_OCE_Interp_empty
1506	WRITE(,) 'agrif_oce_interp : You should not have seen this print! error?'
1507	END SUBROUTINE Agrif_OCE_Interp_empty
1508	#endif
1509
1510	!!======================================================================
1511	END MODULE agrif_oce_interp

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

Download in other formats: