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

Last change on this file since 8861 was 8861, checked in by timgraham, 6 years ago
Added vertical interpolation of avm
Property svn:keywords set to `Id`
File size: 55.8 KB

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

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

Download in other formats: