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 @ 10989

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

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

Download in other formats: