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_opa_interp.F90 in branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC – NEMO

Context Navigation

source: branches/2017/dev_METO_MERCATOR_2017_agrif/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 8999

Last change on this file since 8999 was 8999, checked in by timgraham, 6 years ago
Tidying of interp and added interpavm for vertical refinement case
Property svn:keywords set to `Id`
File size: 55.5 KB

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

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

Download in other formats: