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_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC – NEMO

Context Navigation

source: branches/2017/dev_r8126_UKMO_AGRIF_vert_interp/NEMOGCM/NEMO/NST_SRC/agrif_opa_interp.F90 @ 8690

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

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

Download in other formats: