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

Last change on this file since 8822 was 8822, checked in by timgraham, 6 years ago
Bug fixes in interp and update. Modified update to deal with none matching depths at U/V points
Property svn:keywords set to `Id`
File size: 56.4 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)*zrhoxy)
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	zrhox = Agrif_Rhox()
719	!
720	zalpha1 = ( zrhox - 1. ) * 0.5
721	zalpha2 = 1. - zalpha1
722	!
723	zalpha3 = ( zrhox - 1. ) / ( zrhox + 1. )
724	zalpha4 = 1. - zalpha3
725	!
726	zalpha6 = 2. * ( zrhox - 1. ) / ( zrhox + 1. )
727	zalpha7 = - ( zrhox - 1. ) / ( zrhox + 3. )
728	zalpha5 = 1. - zalpha6 - zalpha7
729	!
730	imin = i1
731	imax = i2
732	jmin = j1
733	jmax = j2
734	!
735	! Remove CORNERS
736	IF((nbondj == -1).OR.(nbondj == 2)) jmin = 3
737	IF((nbondj == +1).OR.(nbondj == 2)) jmax = nlcj-2
738	IF((nbondi == -1).OR.(nbondi == 2)) imin = 3
739	IF((nbondi == +1).OR.(nbondi == 2)) imax = nlci-2
740	!
741	IF( eastern_side ) THEN
742	DO jn = 1, jpts
743	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)
744	DO jk = 1, jpkm1
745	DO jj = jmin,jmax
746	IF( umask(nlci-2,jj,jk) == 0._wp ) THEN
747	tsa(nlci-1,jj,jk,jn) = tsa(nlci,jj,jk,jn) * tmask(nlci-1,jj,jk)
748	ELSE
749	tsa(nlci-1,jj,jk,jn)=(zalpha4tsa(nlci,jj,jk,jn)+zalpha3tsa(nlci-2,jj,jk,jn))*tmask(nlci-1,jj,jk)
750	IF( un(nlci-2,jj,jk) > 0._wp ) THEN
751	tsa(nlci-1,jj,jk,jn)=( zalpha6tsa(nlci-2,jj,jk,jn)+zalpha5tsa(nlci,jj,jk,jn) &
752	+ zalpha7tsa(nlci-3,jj,jk,jn) ) tmask(nlci-1,jj,jk)
753	ENDIF
754	ENDIF
755	END DO
756	END DO
757	tsa(nlci,j1:j2,k1:k2,jn) = 0._wp
758	END DO
759	ENDIF
760	!
761	IF( northern_side ) THEN
762	DO jn = 1, jpts
763	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)
764	DO jk = 1, jpkm1
765	DO ji = imin,imax
766	IF( vmask(ji,nlcj-2,jk) == 0._wp ) THEN
767	tsa(ji,nlcj-1,jk,jn) = tsa(ji,nlcj,jk,jn) * tmask(ji,nlcj-1,jk)
768	ELSE
769	tsa(ji,nlcj-1,jk,jn)=(zalpha4tsa(ji,nlcj,jk,jn)+zalpha3tsa(ji,nlcj-2,jk,jn))*tmask(ji,nlcj-1,jk)
770	IF (vn(ji,nlcj-2,jk) > 0._wp ) THEN
771	tsa(ji,nlcj-1,jk,jn)=( zalpha6tsa(ji,nlcj-2,jk,jn)+zalpha5tsa(ji,nlcj,jk,jn) &
772	+ zalpha7tsa(ji,nlcj-3,jk,jn) ) tmask(ji,nlcj-1,jk)
773	ENDIF
774	ENDIF
775	END DO
776	END DO
777	tsa(i1:i2,nlcj,k1:k2,jn) = 0._wp
778	END DO
779	ENDIF
780	!
781	IF( western_side ) THEN
782	DO jn = 1, jpts
783	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)
784	DO jk = 1, jpkm1
785	DO jj = jmin,jmax
786	IF( umask(2,jj,jk) == 0._wp ) THEN
787	tsa(2,jj,jk,jn) = tsa(1,jj,jk,jn) * tmask(2,jj,jk)
788	ELSE
789	tsa(2,jj,jk,jn)=(zalpha4tsa(1,jj,jk,jn)+zalpha3tsa(3,jj,jk,jn))*tmask(2,jj,jk)
790	IF( un(2,jj,jk) < 0._wp ) THEN
791	tsa(2,jj,jk,jn)=(zalpha6tsa(3,jj,jk,jn)+zalpha5tsa(1,jj,jk,jn)+zalpha7tsa(4,jj,jk,jn))tmask(2,jj,jk)
792	ENDIF
793	ENDIF
794	END DO
795	END DO
796	tsa(1,j1:j2,k1:k2,jn) = 0._wp
797	END DO
798	ENDIF
799	!
800	IF( southern_side ) THEN
801	DO jn = 1, jpts
802	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)
803	DO jk = 1, jpk
804	DO ji=imin,imax
805	IF( vmask(ji,2,jk) == 0._wp ) THEN
806	tsa(ji,2,jk,jn)=tsa(ji,1,jk,jn) * tmask(ji,2,jk)
807	ELSE
808	tsa(ji,2,jk,jn)=(zalpha4tsa(ji,1,jk,jn)+zalpha3tsa(ji,3,jk,jn))*tmask(ji,2,jk)
809	IF( vn(ji,2,jk) < 0._wp ) THEN
810	tsa(ji,2,jk,jn)=(zalpha6tsa(ji,3,jk,jn)+zalpha5tsa(ji,1,jk,jn)+zalpha7tsa(ji,4,jk,jn))tmask(ji,2,jk)
811	ENDIF
812	ENDIF
813	END DO
814	END DO
815	tsa(i1:i2,1,k1:k2,jn) = 0._wp
816	END DO
817	ENDIF
818	!
819	! Treatment of corners
820	!
821	! East south
822	IF ((eastern_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
823	tsa(nlci-1,2,:,:) = ptab_child(nlci-1,2,:,1:jpts)
824	ENDIF
825	! East north
826	IF ((eastern_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
827	tsa(nlci-1,nlcj-1,:,:) = ptab_child(nlci-1,nlcj-1,:,1:jpts)
828	ENDIF
829	! West south
830	IF ((western_side).AND.((nbondj == -1).OR.(nbondj == 2))) THEN
831	tsa(2,2,:,:) = ptab_child(2,2,:,1:jpts)
832	ENDIF
833	! West north
834	IF ((western_side).AND.((nbondj == 1).OR.(nbondj == 2))) THEN
835	tsa(2,nlcj-1,:,:) = ptab_child(2,nlcj-1,:,1:jpts)
836	ENDIF
837	!
838	ENDIF
839	!
840	END SUBROUTINE interptsn
841
842
843	SUBROUTINE interpsshn( ptab, i1, i2, j1, j2, before, nb, ndir )
844	!!----------------------------------------------------------------------
845	!! * ROUTINE interpsshn *
846	!!----------------------------------------------------------------------
847	INTEGER , INTENT(in ) :: i1, i2, j1, j2
848	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
849	LOGICAL , INTENT(in ) :: before
850	INTEGER , INTENT(in ) :: nb , ndir
851	!
852	LOGICAL :: western_side, eastern_side,northern_side,southern_side
853	!!----------------------------------------------------------------------
854	!
855	IF( before) THEN
856	ptab(i1:i2,j1:j2) = sshn(i1:i2,j1:j2)
857	ELSE
858	western_side = (nb == 1).AND.(ndir == 1)
859	eastern_side = (nb == 1).AND.(ndir == 2)
860	southern_side = (nb == 2).AND.(ndir == 1)
861	northern_side = (nb == 2).AND.(ndir == 2)
862	IF(western_side) hbdy_w(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
863	IF(eastern_side) hbdy_e(j1:j2) = ptab(i1,j1:j2) * tmask(i1,j1:j2,1)
864	IF(southern_side) hbdy_s(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
865	IF(northern_side) hbdy_n(i1:i2) = ptab(i1:i2,j1) * tmask(i1:i2,j1,1)
866	ENDIF
867	!
868	END SUBROUTINE interpsshn
869
870
871	SUBROUTINE interpun( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
872	!!----------------------------------------------------------------------
873	!! * ROUTINE interpun *
874	!!---------------------------------------------
875	!!
876	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
877	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
878	LOGICAL, INTENT(in) :: before
879	INTEGER, INTENT(in) :: nb , ndir
880	!!
881	INTEGER :: ji,jj,jk
882	REAL(wp) :: zrhoy
883	! VERTICAL REFINEMENT BEGIN
884	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
885	REAL(wp), DIMENSION(k1:k2) :: tabin
886	REAL(wp) :: h_in(k1:k2)
887	REAL(wp) :: h_out(1:jpk)
888	INTEGER :: N_in, N_out
889	REAL(wp) :: h_diff
890	LOGICAL :: western_side, eastern_side
891	INTEGER :: iref
892
893	! VERTICAL REFINEMENT END
894	!!---------------------------------------------
895	!
896	zrhoy = Agrif_rhoy()
897	IF (before) THEN
898	!We can't use zero as the special value because we need to include zeros
899	!when interpolating the scale factors
900	IF(Agrif_UseSpecialValue) THEN
901	! Agrif_SpecialValue = -999._wp
902	Agrif_SpecialValue = 0._wp
903	ELSE
904	Agrif_SpecialValue = 0._wp
905	ENDIF
906	DO jk=1,jpk
907	DO jj=j1,j2
908	DO ji=i1,i2
909	ptab(ji,jj,jk,1) = (e2u(ji,jj) * e3u_n(ji,jj,jk) * un(ji,jj,jk)*umask(ji,jj,jk)) - &
910	& ((umask(ji,jj,jk)-1) * Agrif_SpecialValue)
911	ptab(ji,jj,jk,2) = (umask(ji,jj,jk) * e2u(ji,jj) * e3u_n(ji,jj,jk))
912	END DO
913	END DO
914	END DO
915	ELSE
916	! VERTICAL REFINEMENT BEGIN
917	western_side = (nb == 1).AND.(ndir == 1)
918	eastern_side = (nb == 1).AND.(ndir == 2)
919
920	Agrif_SpecialValue = 0._wp ! reset specialvalue to zero now interpolation completed
921
922	ptab_child(:,:,:) = 0.
923	#ifdef key_vertical
924	! VERTICAL REFINEMENT BEGIN
925	DO jj=j1,j2
926	DO ji=i1,i2
927	iref = ji
928	IF (western_side) iref = 2
929	IF (eastern_side) iref = nlci-2
930
931	N_in = 0
932	DO jk=k1,k2
933	IF (ptab(ji,jj,jk,2) == 0) EXIT
934	N_in = N_in + 1
935	tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
936	h_in(N_in) = ptab(ji,jj,jk,2)/(e2u(ji,jj)*zrhoy)
937	ENDDO
938
939	IF (N_in == 0) THEN
940	ptab_child(ji,jj,:) = 0.
941	CYCLE
942	ENDIF
943
944	N_out = 0
945	DO jk=1,jpk
946	if (umask(ji,jj,jk) == 0) EXIT
947	N_out = N_out + 1
948	h_out(N_out) = e3u_n(ji,jj,jk)
949	ENDDO
950
951	IF (N_out == 0) THEN
952	ptab_child(ji,jj,:) = 0.
953	CYCLE
954	ENDIF
955
956	IF (N_in * N_out > 0) THEN
957	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
958	! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
959	if (h_diff < -1.e4) then
960	print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
961	! stop
962	endif
963	ENDIF
964	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)
965	ENDDO
966	ENDDO
967
968	! VERTICAL REFINEMENT END
969	DO jk = 1, jpkm1
970	DO jj=j1,j2
971	ua(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk)
972	END DO
973	END DO
974	#else
975	DO jk = 1, jpkm1
976	DO jj=j1,j2
977	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) )
978	END DO
979	END DO
980	#endif
981	ENDIF
982	!
983	END SUBROUTINE interpun
984
985
986	SUBROUTINE interpvn( ptab, i1, i2, j1, j2, k1, k2, m1, m2, before, nb, ndir )
987	!!----------------------------------------------------------------------
988	!! * ROUTINE interpvn *
989	!!----------------------------------------------------------------------
990	!
991	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
992	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: ptab
993	LOGICAL, INTENT(in) :: before
994	INTEGER, INTENT(in) :: nb , ndir
995	!
996	INTEGER :: ji,jj,jk
997	REAL(wp) :: zrhox
998	! VERTICAL REFINEMENT BEGIN
999	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ptab_child
1000	REAL(wp), DIMENSION(k1:k2) :: tabin
1001	REAL(wp) :: h_in(k1:k2)
1002	REAL(wp) :: h_out(1:jpk)
1003	INTEGER :: N_in, N_out
1004	REAL(wp) :: h_diff
1005	LOGICAL :: northern_side,southern_side
1006	INTEGER :: jref
1007
1008	! VERTICAL REFINEMENT END
1009	!!---------------------------------------------
1010	!
1011	zrhox = Agrif_rhox()
1012	IF (before) THEN
1013	IF(Agrif_UseSpecialValue) THEN
1014	! Agrif_SpecialValue = -999._wp
1015	Agrif_SpecialValue = 0._wp
1016	ELSE
1017	Agrif_SpecialValue = 0._wp
1018	ENDIF
1019	DO jk=k1,k2
1020	DO jj=j1,j2
1021	DO ji=i1,i2
1022	ptab(ji,jj,jk,1) = e1v(ji,jj) * e3v_n(ji,jj,jk) * vn(ji,jj,jk)
1023	ptab(ji,jj,jk,2) = vmask(ji,jj,jk) * e1v(ji,jj) * e3v_n(ji,jj,jk)
1024	END DO
1025	END DO
1026	END DO
1027	ELSE
1028	Agrif_SpecialValue = 0._wp !Reset special value to zero now interpolation is done
1029	ptab_child(:,:,:) = 0.
1030	#ifdef key_vertical
1031	! VERTICAL REFINEMENT BEGIN
1032	southern_side = (nb == 2).AND.(ndir == 1)
1033	northern_side = (nb == 2).AND.(ndir == 2)
1034
1035	do jj=j1,j2
1036	jref = jj
1037	IF (southern_side) jref = 2
1038	IF (northern_side) jref = nlcj-2
1039	do ji=i1,i2
1040
1041	N_in = 0
1042	do jk=k1,k2
1043	if (ptab(ji,jj,jk,2) == 0) EXIT
1044	N_in = N_in + 1
1045	tabin(jk) = ptab(ji,jj,jk,1)/ptab(ji,jj,jk,2)
1046	h_in(N_in) = ptab(ji,jj,jk,2)/(e1v(ji,jj)*zrhox)
1047	enddo
1048	IF (N_in == 0) THEN
1049	ptab_child(ji,jj,:) = 0.
1050	CYCLE
1051	ENDIF
1052
1053	N_out = 0
1054	do jk=1,jpk
1055	if (vmask(ji,jref,jk) == 0) EXIT
1056	N_out = N_out + 1
1057	h_out(N_out) = e3v_n(ji,jj,jk)
1058	enddo
1059	IF (N_out == 0) THEN
1060	ptab_child(ji,jj,:) = 0.
1061	CYCLE
1062	ENDIF
1063
1064	! IF (N_in * N_out > 0) THEN
1065	! h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
1066	! Should be able to remove the next IF/ELSEIF statement once scale factors are dealt with properly
1067	! if (h_diff < 0.) then
1068	! print *,'CHECK YOUR BATHY interpvn...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
1069	! stop
1070	! endif
1071	! ENDIF
1072	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)
1073
1074	enddo
1075	enddo
1076	! in the following
1077	! remove division of va by fs e3v, zrhox and e1v (already done)
1078	DO jk=1,jpkm1
1079	DO jj=j1,j2
1080	va(i1:i2,jj,jk) = ptab_child(i1:i2,jj,jk)
1081	END DO
1082	END DO
1083	! VERTICAL REFINEMENT END
1084	#else
1085	DO jk=1,jpkm1
1086	va(i1:i2,j1:j2,jk) = vmask(i1:i2,j1:j2,jk) * ptab(i1:i2,j1:j2,jk,1) / &
1087	& ( zrhox * e1v(i1:i2,j1:j2) * e3v_n(i1:i2,j1:j2,jk) )
1088	END DO
1089	#endif
1090	ENDIF
1091	!
1092	END SUBROUTINE interpvn
1093
1094
1095	SUBROUTINE interpunb( ptab, i1, i2, j1, j2, before, nb, ndir )
1096	!!----------------------------------------------------------------------
1097	!! * ROUTINE interpunb *
1098	!!----------------------------------------------------------------------
1099	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1100	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1101	LOGICAL , INTENT(in ) :: before
1102	INTEGER , INTENT(in ) :: nb , ndir
1103	!
1104	INTEGER :: ji, jj
1105	REAL(wp) :: zrhoy, zrhot, zt0, zt1, ztcoeff
1106	LOGICAL :: western_side, eastern_side,northern_side,southern_side
1107	!!----------------------------------------------------------------------
1108	!
1109	IF( before ) THEN
1110	ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * hu_n(i1:i2,j1:j2) * un_b(i1:i2,j1:j2)
1111	ELSE
1112	western_side = (nb == 1).AND.(ndir == 1)
1113	eastern_side = (nb == 1).AND.(ndir == 2)
1114	southern_side = (nb == 2).AND.(ndir == 1)
1115	northern_side = (nb == 2).AND.(ndir == 2)
1116	zrhoy = Agrif_Rhoy()
1117	zrhot = Agrif_rhot()
1118	! Time indexes bounds for integration
1119	zt0 = REAL(Agrif_NbStepint() , wp) / zrhot
1120	zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1121	! Polynomial interpolation coefficients:
1122	IF( bdy_tinterp == 1 ) THEN
1123	ztcoeff = zrhot * ( zt1*2._wp ( zt1 - 1._wp) &
1124	& - zt0*2._wp ( zt0 - 1._wp) )
1125	ELSEIF( bdy_tinterp == 2 ) THEN
1126	ztcoeff = zrhot * ( zt1 * ( zt1 - 1._wp)**2._wp &
1127	& - zt0 * ( zt0 - 1._wp)**2._wp )
1128
1129	ELSE
1130	ztcoeff = 1
1131	ENDIF
1132	!
1133	IF(western_side) THEN
1134	ubdy_w(j1:j2) = ubdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2)
1135	ENDIF
1136	IF(eastern_side) THEN
1137	ubdy_e(j1:j2) = ubdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2)
1138	ENDIF
1139	IF(southern_side) THEN
1140	ubdy_s(i1:i2) = ubdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1)
1141	ENDIF
1142	IF(northern_side) THEN
1143	ubdy_n(i1:i2) = ubdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1)
1144	ENDIF
1145	!
1146	IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1147	IF(western_side) THEN
1148	ubdy_w(j1:j2) = ubdy_w(j1:j2) / (zrhoye2u(i1,j1:j2)) umask(i1,j1:j2,1)
1149	ENDIF
1150	IF(eastern_side) THEN
1151	ubdy_e(j1:j2) = ubdy_e(j1:j2) / (zrhoye2u(i1,j1:j2)) umask(i1,j1:j2,1)
1152	ENDIF
1153	IF(southern_side) THEN
1154	ubdy_s(i1:i2) = ubdy_s(i1:i2) / (zrhoye2u(i1:i2,j1)) umask(i1:i2,j1,1)
1155	ENDIF
1156	IF(northern_side) THEN
1157	ubdy_n(i1:i2) = ubdy_n(i1:i2) / (zrhoye2u(i1:i2,j1)) umask(i1:i2,j1,1)
1158	ENDIF
1159	ENDIF
1160	ENDIF
1161	!
1162	END SUBROUTINE interpunb
1163
1164
1165	SUBROUTINE interpvnb( ptab, i1, i2, j1, j2, before, nb, ndir )
1166	!!----------------------------------------------------------------------
1167	!! * ROUTINE interpvnb *
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) :: zrhox, zrhot, zt0, zt1, ztcoeff
1176	LOGICAL :: western_side, eastern_side,northern_side,southern_side
1177	!!----------------------------------------------------------------------
1178	!
1179	IF( before ) THEN
1180	ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * hv_n(i1:i2,j1:j2) * vn_b(i1:i2,j1:j2)
1181	ELSE
1182	western_side = (nb == 1).AND.(ndir == 1)
1183	eastern_side = (nb == 1).AND.(ndir == 2)
1184	southern_side = (nb == 2).AND.(ndir == 1)
1185	northern_side = (nb == 2).AND.(ndir == 2)
1186	zrhox = Agrif_Rhox()
1187	zrhot = Agrif_rhot()
1188	! Time indexes bounds for integration
1189	zt0 = REAL(Agrif_NbStepint() , wp) / zrhot
1190	zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1191	IF( bdy_tinterp == 1 ) THEN
1192	ztcoeff = zrhot * ( zt1*2._wp ( zt1 - 1._wp) &
1193	& - zt0*2._wp ( zt0 - 1._wp) )
1194	ELSEIF( bdy_tinterp == 2 ) THEN
1195	ztcoeff = zrhot * ( zt1 * ( zt1 - 1._wp)**2._wp &
1196	& - zt0 * ( zt0 - 1._wp)**2._wp )
1197	ELSE
1198	ztcoeff = 1
1199	ENDIF
1200	!
1201	IF(western_side) THEN
1202	vbdy_w(j1:j2) = vbdy_w(j1:j2) + ztcoeff * ptab(i1,j1:j2)
1203	ENDIF
1204	IF(eastern_side) THEN
1205	vbdy_e(j1:j2) = vbdy_e(j1:j2) + ztcoeff * ptab(i1,j1:j2)
1206	ENDIF
1207	IF(southern_side) THEN
1208	vbdy_s(i1:i2) = vbdy_s(i1:i2) + ztcoeff * ptab(i1:i2,j1)
1209	ENDIF
1210	IF(northern_side) THEN
1211	vbdy_n(i1:i2) = vbdy_n(i1:i2) + ztcoeff * ptab(i1:i2,j1)
1212	ENDIF
1213	!
1214	IF( bdy_tinterp == 0 .OR. bdy_tinterp == 2) THEN
1215	IF(western_side) THEN
1216	vbdy_w(j1:j2) = vbdy_w(j1:j2) / (zrhox*e1v(i1,j1:j2)) &
1217	& * vmask(i1,j1:j2,1)
1218	ENDIF
1219	IF(eastern_side) THEN
1220	vbdy_e(j1:j2) = vbdy_e(j1:j2) / (zrhox*e1v(i1,j1:j2)) &
1221	& * vmask(i1,j1:j2,1)
1222	ENDIF
1223	IF(southern_side) THEN
1224	vbdy_s(i1:i2) = vbdy_s(i1:i2) / (zrhox*e1v(i1:i2,j1)) &
1225	& * vmask(i1:i2,j1,1)
1226	ENDIF
1227	IF(northern_side) THEN
1228	vbdy_n(i1:i2) = vbdy_n(i1:i2) / (zrhox*e1v(i1:i2,j1)) &
1229	& * vmask(i1:i2,j1,1)
1230	ENDIF
1231	ENDIF
1232	ENDIF
1233	!
1234	END SUBROUTINE interpvnb
1235
1236
1237	SUBROUTINE interpub2b( ptab, i1, i2, j1, j2, before, nb, ndir )
1238	!!----------------------------------------------------------------------
1239	!! * ROUTINE interpub2b *
1240	!!----------------------------------------------------------------------
1241	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1242	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1243	LOGICAL , INTENT(in ) :: before
1244	INTEGER , INTENT(in ) :: nb , ndir
1245	!
1246	INTEGER :: ji,jj
1247	REAL(wp) :: zrhot, zt0, zt1,zat
1248	LOGICAL :: western_side, eastern_side,northern_side,southern_side
1249	!!----------------------------------------------------------------------
1250	IF( before ) THEN
1251	ptab(i1:i2,j1:j2) = e2u(i1:i2,j1:j2) * ub2_b(i1:i2,j1:j2)
1252	ELSE
1253	western_side = (nb == 1).AND.(ndir == 1)
1254	eastern_side = (nb == 1).AND.(ndir == 2)
1255	southern_side = (nb == 2).AND.(ndir == 1)
1256	northern_side = (nb == 2).AND.(ndir == 2)
1257	zrhot = Agrif_rhot()
1258	! Time indexes bounds for integration
1259	zt0 = REAL(Agrif_NbStepint() , wp) / zrhot
1260	zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1261	! Polynomial interpolation coefficients:
1262	zat = zrhot * ( zt1*2._wp (-2._wp*zt1 + 3._wp) &
1263	& - zt0*2._wp (-2._wp*zt0 + 3._wp) )
1264	!
1265	IF(western_side ) ubdy_w(j1:j2) = zat * ptab(i1,j1:j2)
1266	IF(eastern_side ) ubdy_e(j1:j2) = zat * ptab(i1,j1:j2)
1267	IF(southern_side) ubdy_s(i1:i2) = zat * ptab(i1:i2,j1)
1268	IF(northern_side) ubdy_n(i1:i2) = zat * ptab(i1:i2,j1)
1269	ENDIF
1270	!
1271	END SUBROUTINE interpub2b
1272
1273
1274	SUBROUTINE interpvb2b( ptab, i1, i2, j1, j2, before, nb, ndir )
1275	!!----------------------------------------------------------------------
1276	!! * ROUTINE interpvb2b *
1277	!!----------------------------------------------------------------------
1278	INTEGER , INTENT(in ) :: i1, i2, j1, j2
1279	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: ptab
1280	LOGICAL , INTENT(in ) :: before
1281	INTEGER , INTENT(in ) :: nb , ndir
1282	!
1283	INTEGER :: ji,jj
1284	REAL(wp) :: zrhot, zt0, zt1,zat
1285	LOGICAL :: western_side, eastern_side,northern_side,southern_side
1286	!!----------------------------------------------------------------------
1287	!
1288	IF( before ) THEN
1289	ptab(i1:i2,j1:j2) = e1v(i1:i2,j1:j2) * vb2_b(i1:i2,j1:j2)
1290	ELSE
1291	western_side = (nb == 1).AND.(ndir == 1)
1292	eastern_side = (nb == 1).AND.(ndir == 2)
1293	southern_side = (nb == 2).AND.(ndir == 1)
1294	northern_side = (nb == 2).AND.(ndir == 2)
1295	zrhot = Agrif_rhot()
1296	! Time indexes bounds for integration
1297	zt0 = REAL(Agrif_NbStepint() , wp) / zrhot
1298	zt1 = REAL(Agrif_NbStepint()+1, wp) / zrhot
1299	! Polynomial interpolation coefficients:
1300	zat = zrhot * ( zt1*2._wp (-2._wp*zt1 + 3._wp) &
1301	& - zt0*2._wp (-2._wp*zt0 + 3._wp) )
1302	!
1303	IF(western_side ) vbdy_w(j1:j2) = zat * ptab(i1,j1:j2)
1304	IF(eastern_side ) vbdy_e(j1:j2) = zat * ptab(i1,j1:j2)
1305	IF(southern_side) vbdy_s(i1:i2) = zat * ptab(i1:i2,j1)
1306	IF(northern_side) vbdy_n(i1:i2) = zat * ptab(i1:i2,j1)
1307	ENDIF
1308	!
1309	END SUBROUTINE interpvb2b
1310
1311
1312	SUBROUTINE interpe3t( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1313	!!----------------------------------------------------------------------
1314	!! * ROUTINE interpe3t *
1315	!!----------------------------------------------------------------------
1316	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1317	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1318	LOGICAL , INTENT(in ) :: before
1319	INTEGER , INTENT(in ) :: nb , ndir
1320	!
1321	INTEGER :: ji, jj, jk
1322	LOGICAL :: western_side, eastern_side, northern_side, southern_side
1323	REAL(wp) :: ztmpmsk
1324	!!----------------------------------------------------------------------
1325	!
1326	IF( before ) THEN
1327	ptab(i1:i2,j1:j2,k1:k2) = tmask(i1:i2,j1:j2,k1:k2) * e3t_0(i1:i2,j1:j2,k1:k2)
1328	ELSE
1329	western_side = (nb == 1).AND.(ndir == 1)
1330	eastern_side = (nb == 1).AND.(ndir == 2)
1331	southern_side = (nb == 2).AND.(ndir == 1)
1332	northern_side = (nb == 2).AND.(ndir == 2)
1333
1334	DO jk = k1, k2
1335	DO jj = j1, j2
1336	DO ji = i1, i2
1337	! Get velocity mask at boundary edge points:
1338	IF( western_side ) ztmpmsk = umask(ji ,jj ,1)
1339	IF( eastern_side ) ztmpmsk = umask(nlci-2,jj ,1)
1340	IF( northern_side) ztmpmsk = vmask(ji ,nlcj-2,1)
1341	IF( southern_side) ztmpmsk = vmask(ji ,2 ,1)
1342	!
1343	IF( ABS( ptab(ji,jj,jk) - tmask(ji,jj,jk) * e3t_0(ji,jj,jk) )*ztmpmsk > 1.D-2) THEN
1344	IF (western_side) THEN
1345	WRITE(numout,*) 'ERROR bathymetry merge at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1346	ELSEIF (eastern_side) THEN
1347	WRITE(numout,*) 'ERROR bathymetry merge at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1348	ELSEIF (southern_side) THEN
1349	WRITE(numout,*) 'ERROR bathymetry merge at the southern border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1350	ELSEIF (northern_side) THEN
1351	WRITE(numout,*) 'ERROR bathymetry merge at the northen border ji,jj,jk', ji+nimpp-1,jj+njmpp-1,jk
1352	ENDIF
1353	WRITE(numout,*) ' ptab(ji,jj,jk), e3t(ji,jj,jk) ', ptab(ji,jj,jk), e3t_0(ji,jj,jk)
1354	kindic_agr = kindic_agr + 1
1355	ENDIF
1356	END DO
1357	END DO
1358	END DO
1359	!
1360	ENDIF
1361	!
1362	END SUBROUTINE interpe3t
1363
1364
1365	SUBROUTINE interpumsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1366	!!----------------------------------------------------------------------
1367	!! * ROUTINE interpumsk *
1368	!!----------------------------------------------------------------------
1369	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1370	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1371	LOGICAL , INTENT(in ) :: before
1372	INTEGER , INTENT(in ) :: nb , ndir
1373	!
1374	INTEGER :: ji, jj, jk
1375	LOGICAL :: western_side, eastern_side
1376	!!----------------------------------------------------------------------
1377	!
1378	IF( before ) THEN
1379	ptab(i1:i2,j1:j2,k1:k2) = umask(i1:i2,j1:j2,k1:k2)
1380	ELSE
1381	western_side = (nb == 1).AND.(ndir == 1)
1382	eastern_side = (nb == 1).AND.(ndir == 2)
1383	DO jk = k1, k2
1384	DO jj = j1, j2
1385	DO ji = i1, i2
1386	! Velocity mask at boundary edge points:
1387	IF (ABS(ptab(ji,jj,jk) - umask(ji,jj,jk)) > 1.D-2) THEN
1388	IF (western_side) THEN
1389	WRITE(numout,*) 'ERROR with umask at the western border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1390	WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1391	kindic_agr = kindic_agr + 1
1392	ELSEIF (eastern_side) THEN
1393	WRITE(numout,*) 'ERROR with umask at the eastern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1394	WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), umask(ji,jj,jk)
1395	kindic_agr = kindic_agr + 1
1396	ENDIF
1397	ENDIF
1398	END DO
1399	END DO
1400	END DO
1401	!
1402	ENDIF
1403	!
1404	END SUBROUTINE interpumsk
1405
1406
1407	SUBROUTINE interpvmsk( ptab, i1, i2, j1, j2, k1, k2, before, nb, ndir )
1408	!!----------------------------------------------------------------------
1409	!! * ROUTINE interpvmsk *
1410	!!----------------------------------------------------------------------
1411	INTEGER , INTENT(in ) :: i1,i2,j1,j2,k1,k2
1412	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1413	LOGICAL , INTENT(in ) :: before
1414	INTEGER , INTENT(in ) :: nb , ndir
1415	!
1416	INTEGER :: ji, jj, jk
1417	LOGICAL :: northern_side, southern_side
1418	!!----------------------------------------------------------------------
1419	!
1420	IF( before ) THEN
1421	ptab(i1:i2,j1:j2,k1:k2) = vmask(i1:i2,j1:j2,k1:k2)
1422	ELSE
1423	southern_side = (nb == 2).AND.(ndir == 1)
1424	northern_side = (nb == 2).AND.(ndir == 2)
1425	DO jk = k1, k2
1426	DO jj = j1, j2
1427	DO ji = i1, i2
1428	! Velocity mask at boundary edge points:
1429	IF (ABS(ptab(ji,jj,jk) - vmask(ji,jj,jk)) > 1.D-2) THEN
1430	IF (southern_side) THEN
1431	WRITE(numout,*) 'ERROR with vmask at the southern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1432	WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1433	kindic_agr = kindic_agr + 1
1434	ELSEIF (northern_side) THEN
1435	WRITE(numout,*) 'ERROR with vmask at the northern border ji,jj,jk ', ji+nimpp-1,jj+njmpp-1,jk
1436	WRITE(numout,*) ' masks: parent, child ', ptab(ji,jj,jk), vmask(ji,jj,jk)
1437	kindic_agr = kindic_agr + 1
1438	ENDIF
1439	ENDIF
1440	END DO
1441	END DO
1442	END DO
1443	!
1444	ENDIF
1445	!
1446	END SUBROUTINE interpvmsk
1447
1448	# if defined key_zdftke
1449
1450	SUBROUTINE interpavm( ptab, i1, i2, j1, j2, k1, k2, before )
1451	!!----------------------------------------------------------------------
1452	!! * ROUTINE interavm *
1453	!!----------------------------------------------------------------------
1454	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2
1455	REAL(wp),DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: ptab
1456	LOGICAL , INTENT(in ) :: before
1457	!!----------------------------------------------------------------------
1458	!
1459	IF( before ) THEN
1460	ptab (i1:i2,j1:j2,k1:k2) = avm_k(i1:i2,j1:j2,k1:k2)
1461	ELSE
1462	avm_k(i1:i2,j1:j2,k1:k2) = ptab (i1:i2,j1:j2,k1:k2)
1463	ENDIF
1464	!
1465	END SUBROUTINE interpavm
1466
1467	# endif /* key_zdftke */
1468
1469	#else
1470	!!----------------------------------------------------------------------
1471	!! Empty module no AGRIF zoom
1472	!!----------------------------------------------------------------------
1473	CONTAINS
1474	SUBROUTINE Agrif_OPA_Interp_empty
1475	WRITE(,) 'agrif_opa_interp : You should not have seen this print! error?'
1476	END SUBROUTINE Agrif_OPA_Interp_empty
1477	#endif
1478
1479	!!======================================================================
1480	END MODULE agrif_opa_interp

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

Download in other formats: