New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

agrif_oce_interp.F90 in NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap/src/NST/agrif_oce_interp.F90 @ 11463

Last change on this file since 11463 was 11463, checked in by acc, 5 years ago
Branch: dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps_rewrite_time_filterswap. Minor bugfix in step.F90 to enable AGRIF SETTE tests to run. Also merged prettification changes to NST routines from the dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps branch. Neither AGRIF_DEMO nor VORTEX restart perfectly (drifting after 8 and 121 timesteps, respectively).
Property svn:keywords set to `Id`
File size: 65.2 KB

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

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

Download in other formats: