Context Navigation

agrif_oce_sponge.F90 @ 15437

Last change on this file since 15437 was 15437, checked in by jchanut, 3 years ago
agrif fixes: i) Add masking prior interpolation on "before" tracers: these seem to be strangly unmasked which induces extrapolation issues. ii) Add flux limiter on tracer sponge (recall sponge is applied on anomalies) to prevent from generating spurious extrema.
Property svn:keywords set to `Id`
File size: 45.6 KB

Rev	Line
[3680]	1	#define SPONGE && define SPONGE_TOP
[390]	2
[9570]	3	MODULE agrif_oce_sponge
[6140]	4	!!======================================================================
[9570]	5	!! * MODULE agrif_oce_interp *
[14227]	6	!! AGRIF: sponge package for the ocean dynamics (OCE)
[6140]	7	!!======================================================================
[9019]	8	!! History : 2.0 ! 2002-06 (XXX) Original cade
	9	!! - ! 2005-11 (XXX)
	10	!! 3.2 ! 2009-04 (R. Benshila)
	11	!! 3.6 ! 2014-09 (R. Benshila)
[6140]	12	!!----------------------------------------------------------------------
[7646]	13	#if defined key_agrif
[9019]	14	!!----------------------------------------------------------------------
	15	!! 'key_agrif' AGRIF zoom
	16	!!----------------------------------------------------------------------
[636]	17	USE par_oce
	18	USE oce
	19	USE dom_oce
[9019]	20	!
[636]	21	USE in_out_manager
[782]	22	USE agrif_oce
[5656]	23	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[12377]	24	USE iom
	25	USE vremap
[390]	26
[636]	27	IMPLICIT NONE
	28	PRIVATE
[390]	29
[14086]	30	PUBLIC Agrif_Sponge, Agrif_Sponge_2d, Agrif_Sponge_Tra, Agrif_Sponge_Dyn
[5656]	31	PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge
[14086]	32	PUBLIC interpunb_sponge, interpvnb_sponge
[390]	33
[12377]	34	!! * Substitutions
[14053]	35	# include "domzgr_substitute.h90"
[12377]	36	# include "do_loop_substitute.h90"
[1156]	37	!!----------------------------------------------------------------------
[9598]	38	!! NEMO/NST 4.0 , NEMO Consortium (2018)
[1156]	39	!! $Id$
[10068]	40	!! Software governed by the CeCILL license (see ./LICENSE)
[1156]	41	!!----------------------------------------------------------------------
[5656]	42	CONTAINS
[636]	43
[782]	44	SUBROUTINE Agrif_Sponge_Tra
[9019]	45	!!----------------------------------------------------------------------
	46	!! * ROUTINE Agrif_Sponge_Tra *
	47	!!----------------------------------------------------------------------
	48	REAL(wp) :: zcoef ! local scalar
[14086]	49	INTEGER :: istart, iend, jstart, jend
[9019]	50	!!----------------------------------------------------------------------
[6140]	51	!
[390]	52	#if defined SPONGE
[9031]	53	!! Assume persistence:
[9056]	54	zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
[9031]	55
[9019]	56	Agrif_SpecialValue = 0._wp
[390]	57	Agrif_UseSpecialValue = .TRUE.
[14086]	58	l_vremap = ln_vert_remap
[9019]	59	tabspongedone_tsn = .FALSE.
	60	!
[14086]	61	CALL Agrif_Bc_Variable( ts_sponge_id, calledweight=zcoef, procname=interptsn_sponge )
[9019]	62	!
[5656]	63	Agrif_UseSpecialValue = .FALSE.
[14086]	64	l_vremap = .FALSE.
[390]	65	#endif
[6140]	66	!
[12377]	67	CALL iom_put( 'agrif_spu', fspu(:,:))
	68	CALL iom_put( 'agrif_spv', fspv(:,:))
	69	!
[636]	70	END SUBROUTINE Agrif_Sponge_Tra
[390]	71
[6140]	72
[782]	73	SUBROUTINE Agrif_Sponge_dyn
[9019]	74	!!----------------------------------------------------------------------
	75	!! * ROUTINE Agrif_Sponge_dyn *
	76	!!----------------------------------------------------------------------
	77	REAL(wp) :: zcoef ! local scalar
	78	!!----------------------------------------------------------------------
	79	!
[390]	80	#if defined SPONGE
[9056]	81	zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
[9031]	82
[13286]	83	Agrif_SpecialValue = 0._wp
[782]	84	Agrif_UseSpecialValue = ln_spc_dyn
[14086]	85	l_vremap = ln_vert_remap
[13286]	86	use_sign_north = .TRUE.
	87	sign_north = -1._wp
[9019]	88	!
[5656]	89	tabspongedone_u = .FALSE.
	90	tabspongedone_v = .FALSE.
[9019]	91	CALL Agrif_Bc_Variable( un_sponge_id, calledweight=zcoef, procname=interpun_sponge )
	92	!
[5656]	93	tabspongedone_u = .FALSE.
	94	tabspongedone_v = .FALSE.
[9019]	95	CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge )
[14086]	96
	97	IF ( nn_shift_bar>0 ) THEN ! then split sponge between 2d and 3d
	98	zcoef = REAL(Agrif_NbStepint(),wp)/REAL(Agrif_rhot()) ! forward tsplit
	99	tabspongedone_u = .FALSE.
	100	tabspongedone_v = .FALSE.
	101	CALL Agrif_Bc_Variable( unb_sponge_id, calledweight=zcoef, procname=interpunb_sponge )
	102	!
	103	tabspongedone_u = .FALSE.
	104	tabspongedone_v = .FALSE.
	105	CALL Agrif_Bc_Variable( vnb_sponge_id, calledweight=zcoef, procname=interpvnb_sponge )
	106	ENDIF
[9019]	107	!
[390]	108	Agrif_UseSpecialValue = .FALSE.
[13286]	109	use_sign_north = .FALSE.
[14086]	110	l_vremap = .FALSE.
	111	!
[390]	112	#endif
[6140]	113	!
[12377]	114	CALL iom_put( 'agrif_spt', fspt(:,:))
	115	CALL iom_put( 'agrif_spf', fspf(:,:))
	116	!
[636]	117	END SUBROUTINE Agrif_Sponge_dyn
[390]	118
[6140]	119
[3680]	120	SUBROUTINE Agrif_Sponge
[9019]	121	!!----------------------------------------------------------------------
	122	!! * ROUTINE Agrif_Sponge *
	123	!!----------------------------------------------------------------------
	124	INTEGER :: ji, jj, ind1, ind2
[12377]	125	INTEGER :: ispongearea, jspongearea
	126	REAL(wp) :: z1_ispongearea, z1_jspongearea
[9019]	127	REAL(wp), DIMENSION(jpi,jpj) :: ztabramp
	128	!!----------------------------------------------------------------------
	129	!
[12377]	130	! Sponge 1d example with:
	131	! iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2
	132	!
	133	!coarse : U T U T U T U
	134	!\| \| \| \| \|
	135	!fine : t u t u t u t u t u t u t u t u t u t u t
[14086]	136	!sponge val:0 1 1 1 1 5/6 4/6 3/6 2/6 1/6 0
[12377]	137	! \| ghost \| <-- sponge area -- > \|
	138	! \| points \| \|
	139	! \|--> dynamical interface
	140
[3680]	141	#if defined SPONGE \|\| defined SPONGE_TOP
[14086]	142	! Define ramp from boundaries towards domain interior at F-points
	143	! Store it in ztabramp
[12377]	144
[14086]	145	ispongearea = nn_sponge_len * Agrif_irhox()
	146	z1_ispongearea = 1._wp / REAL( MAX(ispongearea,1), wp )
	147	jspongearea = nn_sponge_len * Agrif_irhoy()
	148	z1_jspongearea = 1._wp / REAL( MAX(jspongearea,1), wp )
	149
	150	ztabramp(:,:) = 0._wp
	151
[14976]	152	IF( lk_west ) THEN ! --- West --- !
	153	ind1 = nn_hls + nbghostcells ! halo + nbghostcells
	154	ind2 = nn_hls + nbghostcells + ispongearea
[14086]	155	DO ji = mi0(ind1), mi1(ind2)
	156	DO jj = 1, jpj
	157	ztabramp(ji,jj) = REAL(ind2 - mig(ji), wp) * z1_ispongearea
[13216]	158	END DO
[14086]	159	END DO
	160	! ghost cells:
	161	ind1 = 1
[14976]	162	ind2 = nn_hls + nbghostcells ! halo + nbghostcells
[14086]	163	DO ji = mi0(ind1), mi1(ind2)
	164	DO jj = 1, jpj
	165	ztabramp(ji,jj) = 1._wp
[13216]	166	END DO
[14086]	167	END DO
	168	ENDIF
	169	IF( lk_east ) THEN ! --- East --- !
[14976]	170	ind1 = jpiglo - ( nn_hls + nbghostcells -1 ) - ispongearea - 1
	171	ind2 = jpiglo - ( nn_hls + nbghostcells -1 ) - 1 ! halo + land + nbghostcells - 1
[14086]	172	DO ji = mi0(ind1), mi1(ind2)
	173	DO jj = 1, jpj
	174	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea )
[13216]	175	END DO
[14086]	176	END DO
	177	! ghost cells:
[14976]	178	ind1 = jpiglo - ( nn_hls + nbghostcells -1 ) - 1 ! halo + land + nbghostcells - 1
[14086]	179	ind2 = jpiglo - 1
	180	DO ji = mi0(ind1), mi1(ind2)
	181	DO jj = 1, jpj
	182	ztabramp(ji,jj) = 1._wp
[13216]	183	END DO
[14086]	184	END DO
	185	ENDIF
	186	IF( lk_south ) THEN ! --- South --- !
[14976]	187	ind1 = nn_hls + nbghostcells ! halo + nbghostcells
	188	ind2 = nn_hls + nbghostcells + jspongearea
[14086]	189	DO jj = mj0(ind1), mj1(ind2)
	190	DO ji = 1, jpi
	191	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea )
	192	END DO
	193	END DO
	194	! ghost cells:
	195	ind1 = 1
[14976]	196	ind2 = nn_hls + nbghostcells ! halo + nbghostcells
[14086]	197	DO jj = mj0(ind1), mj1(ind2)
	198	DO ji = 1, jpi
	199	ztabramp(ji,jj) = 1._wp
	200	END DO
	201	END DO
	202	ENDIF
	203	IF( lk_north ) THEN ! --- North --- !
[14976]	204	ind1 = jpjglo - ( nn_hls + nbghostcells -1 ) - jspongearea - 1
	205	ind2 = jpjglo - ( nn_hls + nbghostcells -1 ) - 1 ! halo + nbghostcells - 1
[14086]	206	DO jj = mj0(ind1), mj1(ind2)
	207	DO ji = 1, jpi
	208	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea )
	209	END DO
	210	END DO
	211	! ghost cells:
[14976]	212	ind1 = jpjglo - ( nn_hls + nbghostcells -1 ) ! halo + land + nbghostcells - 1
[14086]	213	ind2 = jpjglo
	214	DO jj = mj0(ind1), mj1(ind2)
	215	DO ji = 1, jpi
	216	ztabramp(ji,jj) = 1._wp
	217	END DO
	218	END DO
	219	ENDIF
	220	!
	221	! Tracers
	222	fspu(:,:) = 0._wp
	223	fspv(:,:) = 0._wp
	224	DO_2D( 0, 0, 0, 0 )
	225	fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji,jj-1) ) * ssumask(ji,jj)
	226	fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji-1,jj) ) * ssvmask(ji,jj)
	227	END_2D
[13216]	228
[14086]	229	! Dynamics
	230	fspt(:,:) = 0._wp
	231	fspf(:,:) = 0._wp
	232	DO_2D( 0, 0, 0, 0 )
	233	fspt(ji,jj) = 0.25_wp * ( ztabramp(ji ,jj ) + ztabramp(ji-1,jj ) &
	234	& +ztabramp(ji ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj)
	235	fspf(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1)
	236	END_2D
	237
[14433]	238	CALL lbc_lnk( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp, fspt, 'T', 1._wp, fspf, 'F', 1._wp )
[14086]	239	!
	240	! Remove vertical interpolation where not needed:
	241	! (A null value in mbkx arrays does the job)
	242	WHERE (fspu(:,:) == 0._wp) mbku_parent(:,:) = 0
	243	WHERE (fspv(:,:) == 0._wp) mbkv_parent(:,:) = 0
	244	WHERE (fspt(:,:) == 0._wp) mbkt_parent(:,:) = 0
	245	!
[12377]	246	#endif
[14086]	247	!
	248	END SUBROUTINE Agrif_Sponge
[12377]	249
[3680]	250
[14086]	251	SUBROUTINE Agrif_Sponge_2d
	252	!!----------------------------------------------------------------------
	253	!! * ROUTINE Agrif_Sponge_2d *
	254	!!----------------------------------------------------------------------
	255	INTEGER :: ji, jj, ind1, ind2, ishift, jshift
	256	INTEGER :: ispongearea, jspongearea
	257	REAL(wp) :: z1_ispongearea, z1_jspongearea
	258	REAL(wp), DIMENSION(jpi,jpj) :: ztabramp
	259	!!----------------------------------------------------------------------
	260	!
	261	! Sponge 1d example with:
	262	! iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2
	263	!
	264	!coarse : U T U T U T U
	265	!\| \| \| \| \|
	266	!fine : t u t u t u t u t u t u t u t u t u t u t
	267	!sponge val:0 1 1 1 1 5/6 4/6 3/6 2/6 1/6 0
	268	! \| ghost \| <-- sponge area -- > \|
	269	! \| points \| \|
	270	! \|--> dynamical interface
[4153]	271
[14086]	272	#if defined SPONGE \|\| defined SPONGE_TOP
	273	! Define ramp from boundaries towards domain interior at F-points
	274	! Store it in ztabramp
[12377]	275
[14086]	276	ispongearea = nn_sponge_len * Agrif_irhox()
	277	z1_ispongearea = 1._wp / REAL( MAX(ispongearea,1), wp )
	278	jspongearea = nn_sponge_len * Agrif_irhoy()
	279	z1_jspongearea = 1._wp / REAL( MAX(jspongearea,1), wp )
	280	ishift = nn_shift_bar * Agrif_irhox()
	281	jshift = nn_shift_bar * Agrif_irhoy()
	282
	283	ztabramp(:,:) = 0._wp
	284
	285	IF( lk_west ) THEN ! --- West --- !
[14976]	286	ind1 = nn_hls + nbghostcells + ishift
	287	ind2 = nn_hls + nbghostcells + ishift + ispongearea
[14086]	288	DO ji = mi0(ind1), mi1(ind2)
	289	DO jj = 1, jpj
	290	ztabramp(ji,jj) = REAL(ind2 - mig(ji), wp) * z1_ispongearea
[13286]	291	END DO
[14086]	292	END DO
	293	! ghost cells:
	294	ind1 = 1
[14976]	295	ind2 = nn_hls + nbghostcells + ishift ! halo + nbghostcells
[14086]	296	DO ji = mi0(ind1), mi1(ind2)
	297	DO jj = 1, jpj
	298	ztabramp(ji,jj) = 1._wp
[12377]	299	END DO
[14086]	300	END DO
	301	ENDIF
	302	IF( lk_east ) THEN ! --- East --- !
[14976]	303	ind1 = jpiglo - ( nn_hls + nbghostcells -1 + ishift) - ispongearea - 1
	304	ind2 = jpiglo - ( nn_hls + nbghostcells -1 + ishift) - 1 ! halo + nbghostcells - 1
[14086]	305	DO ji = mi0(ind1), mi1(ind2)
	306	DO jj = 1, jpj
	307	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea )
[13216]	308	END DO
[14086]	309	END DO
	310	! ghost cells:
[14976]	311	ind1 = jpiglo - ( nn_hls + nbghostcells -1 + ishift) - 1 ! halo + nbghostcells - 1
[14086]	312	ind2 = jpiglo - 1
	313	DO ji = mi0(ind1), mi1(ind2)
	314	DO jj = 1, jpj
	315	ztabramp(ji,jj) = 1._wp
[13216]	316	END DO
[14086]	317	END DO
	318	ENDIF
	319	IF( lk_south ) THEN ! --- South --- !
[14976]	320	ind1 = nn_hls + nbghostcells + jshift ! halo + nbghostcells
	321	ind2 = nn_hls + nbghostcells + jshift + jspongearea
[14086]	322	DO jj = mj0(ind1), mj1(ind2)
	323	DO ji = 1, jpi
	324	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea )
[12377]	325	END DO
[14086]	326	END DO
	327	! ghost cells:
	328	ind1 = 1
[14976]	329	ind2 = nn_hls + nbghostcells + jshift ! halo + land + nbghostcells
[14086]	330	DO jj = mj0(ind1), mj1(ind2)
	331	DO ji = 1, jpi
	332	ztabramp(ji,jj) = 1._wp
[12377]	333	END DO
[14086]	334	END DO
	335	ENDIF
	336	IF( lk_north ) THEN ! --- North --- !
[14976]	337	ind1 = jpjglo - ( nn_hls + nbghostcells -1 + jshift) - jspongearea - 1
	338	ind2 = jpjglo - ( nn_hls + nbghostcells -1 + jshift) - 1 ! halo + land + nbghostcells - 1
[14086]	339	DO jj = mj0(ind1), mj1(ind2)
	340	DO ji = 1, jpi
	341	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea )
[12377]	342	END DO
[14086]	343	END DO
	344	! ghost cells:
[14976]	345	ind1 = jpjglo - ( nn_hls + nbghostcells -1 + jshift) ! halo + land + nbghostcells - 1
[14086]	346	ind2 = jpjglo
	347	DO jj = mj0(ind1), mj1(ind2)
	348	DO ji = 1, jpi
	349	ztabramp(ji,jj) = 1._wp
[12377]	350	END DO
[14086]	351	END DO
[4153]	352	ENDIF
[14086]	353	!
[3680]	354	! Tracers
[14086]	355	fspu_2d(:,:) = 0._wp
	356	fspv_2d(:,:) = 0._wp
	357	DO_2D( 0, 0, 0, 0 )
	358	fspu_2d(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji,jj-1) ) * ssumask(ji,jj)
	359	fspv_2d(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji-1,jj) ) * ssvmask(ji,jj)
	360	END_2D
[3680]	361
	362	! Dynamics
[14086]	363	fspt_2d(:,:) = 0._wp
	364	fspf_2d(:,:) = 0._wp
	365	DO_2D( 0, 0, 0, 0 )
	366	fspt_2d(ji,jj) = 0.25_wp * ( ztabramp(ji ,jj ) + ztabramp(ji-1,jj ) &
	367	& +ztabramp(ji ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj)
	368	fspf_2d(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1)
[12377]	369	END_2D
[14433]	370	CALL lbc_lnk( 'agrif_Sponge_2d', fspu_2d, 'U', 1._wp, fspv_2d, 'V', 1._wp, fspt_2d, 'T', 1._wp, fspf_2d, 'F', 1._wp )
[3680]	371	!
	372	#endif
[6140]	373	!
[14086]	374	END SUBROUTINE Agrif_Sponge_2d
[3680]	375
[14086]	376
	377	SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before)
[9019]	378	!!----------------------------------------------------------------------
	379	!! * ROUTINE interptsn_sponge *
	380	!!----------------------------------------------------------------------
	381	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
	382	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
	383	LOGICAL , INTENT(in ) :: before
[6140]	384	!
[5656]	385	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	386	INTEGER :: iku, ikv
[13216]	387	REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot
[15437]	388	REAl(wp) :: zflag, zdmod, zdtot
[14086]	389	REAL(wp), DIMENSION(i1-1:i2,j1-1:j2,jpk) :: ztu, ztv
[9031]	390	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
	391	! vertical interpolation:
	392	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
[14086]	393	REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin, tabin_i
	394	REAL(wp), DIMENSION(k1:k2) :: z_in, z_in_i, h_in_i
	395	REAL(wp), DIMENSION(1:jpk) :: h_out, z_out
[9031]	396	INTEGER :: N_in, N_out
[9019]	397	!!----------------------------------------------------------------------
[5656]	398	!
[6140]	399	IF( before ) THEN
[9031]	400	DO jn = 1, jpts
	401	DO jk=k1,k2
	402	DO jj=j1,j2
	403	DO ji=i1,i2
[15437]	404	! JC: masking is mandatory here: before tracer field seems
	405	! to hold non zero values where tmask=0
	406	tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a) * tmask(ji,jj,jk)
[9031]	407	END DO
	408	END DO
	409	END DO
	410	END DO
	411
[14086]	412	IF ( l_vremap.OR.ln_zps ) THEN
[12377]	413
[14086]	414	! Fill cell depths (i.e. gdept) to be interpolated
	415	! Warning: these are masked, hence extrapolated prior interpolation.
	416	DO jj=j1,j2
	417	DO ji=i1,i2
	418	tabres(ji,jj,k1,jpts+1) = 0.5_wp * tmask(ji,jj,k1) * e3t(ji,jj,k1,Kbb_a)
	419	DO jk=k1+1,k2
	420	tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * &
	421	& ( tabres(ji,jj,jk-1,jpts+1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a)+e3t(ji,jj,jk,Kbb_a)) )
	422	END DO
	423	END DO
	424	END DO
[9031]	425
[14086]	426	! Save ssh at last level:
	427	IF ( .NOT.ln_linssh ) THEN
	428	tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1)
	429	END IF
	430
	431	END IF
	432
	433	ELSE
[6140]	434	!
[14086]	435	IF ( l_vremap ) THEN
[12377]	436
[14086]	437	IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp
[12377]	438
[14086]	439	DO jj=j1,j2
	440	DO ji=i1,i2
	441
	442	tabres_child(ji,jj,:,:) = 0._wp
	443	! Build vertical grids:
	444	N_in = mbkt_parent(ji,jj)
	445	! Input grid (account for partial cells if any):
[14170]	446	IF ( N_in > 0 ) THEN
	447	DO jk=1,N_in
	448	z_in(jk) = tabres(ji,jj,jk,n2) - tabres(ji,jj,k2,n2)
	449	tabin(jk,1:jpts) = tabres(ji,jj,jk,1:jpts)
	450	END DO
[14086]	451
[14170]	452	! Intermediate grid:
	453	DO jk = 1, N_in
	454	h_in_i(jk) = e3t0_parent(ji,jj,jk) * &
	455	& (1._wp + tabres(ji,jj,k2,n2)/(ht0_parent(ji,jj)*ssmask(ji,jj) + 1._wp - ssmask(ji,jj)))
	456	END DO
	457	z_in_i(1) = 0.5_wp * h_in_i(1)
	458	DO jk=2,N_in
	459	z_in_i(jk) = z_in_i(jk-1) + 0.5_wp * ( h_in_i(jk) + h_in_i(jk-1) )
	460	END DO
	461	z_in_i(1:N_in) = z_in_i(1:N_in) - tabres(ji,jj,k2,n2)
	462	END IF
[14086]	463	! Output (Child) grid:
	464	N_out = mbkt(ji,jj)
	465	DO jk=1,N_out
	466	h_out(jk) = e3t(ji,jj,jk,Kbb_a)
	467	END DO
	468	z_out(1) = 0.5_wp * h_out(1)
	469	DO jk=2,N_out
	470	z_out(jk) = z_out(jk-1) + 0.5_wp * ( h_out(jk)+h_out(jk-1) )
	471	END DO
	472	IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Kbb_a)
	473
	474	! Account for small differences in the free-surface
	475	IF ( sum(h_out(1:N_out)) > sum(h_in_i(1:N_in) )) THEN
	476	h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in_i(1:N_in)) )
[12377]	477	ELSE
[14086]	478	h_in_i(1)= h_in_i(1) - ( sum(h_in_i(1:N_in))-sum(h_out(1:N_out)) )
	479	END IF
	480	IF (N_in*N_out > 0) THEN
	481	CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),tabin_i(1:N_in,1:jpts),z_in_i(1:N_in),N_in,N_in,jpts)
	482	CALL reconstructandremap(tabin_i(1:N_in,1:jpts),h_in_i(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
	483	! CALL remap_linear(tabin(1:N_in,1:jpts),z_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),z_out(1:N_in),N_in,N_out,jpts)
[12377]	484	ENDIF
[9031]	485	END DO
[14086]	486	END DO
	487
	488	DO jj=j1,j2
	489	DO ji=i1,i2
	490	DO jk=1,jpkm1
	491	tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk)
	492	END DO
[13286]	493	END DO
	494	END DO
[9031]	495
[14086]	496	ELSE
	497
	498	IF ( Agrif_Parent(ln_zps) ) THEN ! Account for partial cells
	499
	500	DO jj=j1,j2
	501	DO ji=i1,i2
	502	!
	503	N_in = mbkt(ji,jj)
	504	N_out = mbkt(ji,jj)
	505	z_in(1) = tabres(ji,jj,1,n2)
	506	tabin(1,1:jpts) = tabres(ji,jj,1,1:jpts)
	507	DO jk=2, N_in
	508	z_in(jk) = tabres(ji,jj,jk,n2)
	509	tabin(jk,1:jpts) = tabres(ji,jj,jk,1:jpts)
	510	END DO
	511	IF (.NOT.ln_linssh) z_in(1:N_in) = z_in(1:N_in) - tabres(ji,jj,k2,n2)
	512
	513	z_out(1) = 0.5_wp * e3t(ji,jj,1,Kbb_a)
	514	DO jk=2, N_out
	515	z_out(jk) = z_out(jk-1) + 0.5_wp * (e3t(ji,jj,jk-1,Kbb_a) + e3t(ji,jj,jk,Kbb_a))
	516	END DO
	517	IF (.NOT.ln_linssh) z_out(1:N_out) = z_out(1:N_out) - ssh(ji,jj,Kbb_a)
	518
	519	CALL remap_linear(tabin(1:N_in,1:jpts), z_in(1:N_in), tabres(ji,jj,1:N_out,1:jpts), &
	520	& z_out(1:N_out), N_in, N_out, jpts)
	521	END DO
[13286]	522	END DO
[14086]	523	ENDIF
	524
	525	DO jj=j1,j2
	526	DO ji=i1,i2
	527	DO jk=1,jpkm1
	528	tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk)
	529	END DO
	530	END DO
[13286]	531	END DO
[9031]	532
[14086]	533	END IF
	534
[5656]	535	DO jn = 1, jpts
	536	DO jk = 1, jpkm1
[14086]	537	ztu(i1-1:i2,j1-1:j2,jk) = 0._wp
[9806]	538	DO jj = j1,j2
[5656]	539	DO ji = i1,i2-1
[14086]	540	zabe1 = rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
[15437]	541	zdtot = tsbdiff(ji+1,jj,jk,jn) - tsbdiff(ji,jj,jk,jn)
	542	zdmod = ts(ji+1,jj,jk,jn,Kbb_a) - ts(ji,jj,jk,jn,Kbb_a)
	543	zflag = 0.5_wp + SIGN(0.5_wp, zdtot*zdmod)
	544	ztu(ji,jj,jk) = zabe1 * fspu(ji,jj) * ( zflag * zdtot + (1._wp - zflag) * zdmod )
[9806]	545	END DO
	546	END DO
[14086]	547	ztv(i1-1:i2,j1-1:j2,jk) = 0._wp
[9806]	548	DO ji = i1,i2
	549	DO jj = j1,j2-1
[14086]	550	zabe2 = rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
	551	ztv(ji,jj,jk) = zabe2 * fspv(ji,jj) * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
[15437]	552	zdtot = tsbdiff(ji,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn)
	553	zdmod = ts(ji,jj+1,jk,jn,Kbb_a) - ts(ji,jj,jk,jn,Kbb_a)
	554	zflag = 0.5_wp + SIGN(0.5_wp, zdtot*zdmod)
	555	ztv(ji,jj,jk) = zabe2 * fspv(ji,jj) * ( zflag * zdtot + (1._wp - zflag) * zdmod )
[6140]	556	END DO
	557	END DO
	558	!
[5656]	559	IF( ln_zps ) THEN ! set gradient at partial step level
[9806]	560	DO jj = j1,j2
	561	DO ji = i1,i2
[5656]	562	! last level
	563	iku = mbku(ji,jj)
	564	ikv = mbkv(ji,jj)
[6140]	565	IF( iku == jk ) ztu(ji,jj,jk) = 0._wp
	566	IF( ikv == jk ) ztv(ji,jj,jk) = 0._wp
[5656]	567	END DO
	568	END DO
	569	ENDIF
[6140]	570	END DO
	571	!
[14086]	572	! JC: there is something wrong with the Laplacian in corners
[5656]	573	DO jk = 1, jpkm1
[14086]	574	DO jj = j1,j2
	575	DO ji = i1,i2
[5656]	576	IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
[12377]	577	zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
[5656]	578	! horizontal diffusive trends
[14086]	579	ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) &
	580	& + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) &
	581	& - rn_trelax_tra * r1_Dt * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn)
[5656]	582	! add it to the general tracer trends
[12377]	583	ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa
[5656]	584	ENDIF
[6140]	585	END DO
	586	END DO
[14086]	587
[6140]	588	END DO
	589	!
	590	END DO
	591	!
[14086]	592	tabspongedone_tsn(i1:i2,j1:j2) = .TRUE.
[6140]	593	!
[5656]	594	ENDIF
[6140]	595	!
[5656]	596	END SUBROUTINE interptsn_sponge
	597
[13286]	598
[9031]	599	SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
	600	!!---------------------------------------------
	601	!! * ROUTINE interpun_sponge *
	602	!!---------------------------------------------
	603	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
	604	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
	605	LOGICAL, INTENT(in) :: before
[6140]	606
[13286]	607	INTEGER :: ji,jj,jk,jmax
	608	INTEGER :: ind1
[5656]	609	! sponge parameters
[13216]	610	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot
[9031]	611	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
	612	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
	613	! vertical interpolation:
	614	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
	615	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
	616	REAL(wp), DIMENSION(1:jpk) :: h_out
[12377]	617	INTEGER ::N_in, N_out
[9031]	618	!!---------------------------------------------
[5656]	619	!
[6140]	620	IF( before ) THEN
[12377]	621	DO jk=k1,k2
[9031]	622	DO jj=j1,j2
	623	DO ji=i1,i2
[15437]	624	tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a) * umask(ji,jj,jk)
[9031]	625	END DO
	626	END DO
	627	END DO
	628
[14086]	629	IF ( l_vremap ) THEN
	630
	631	DO jk=k1,k2
	632	DO jj=j1,j2
	633	DO ji=i1,i2
	634	tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk)
	635	END DO
[12377]	636	END DO
[14086]	637	END DO
	638
	639	! Extrapolate thicknesses in partial bottom cells:
	640	! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
	641	IF (ln_zps) THEN
	642	DO jj=j1,j2
	643	DO ji=i1,i2
	644	jk = mbku(ji,jj)
	645	tabres(ji,jj,jk,m2) = 0._wp
	646	END DO
	647	END DO
	648	END IF
	649	! Save ssh at last level:
	650	tabres(i1:i2,j1:j2,k2,m2) = 0._wp
	651	IF (.NOT.ln_linssh) THEN
	652	! This vertical sum below should be replaced by the sea-level at U-points (optimization):
	653	DO jk=1,jpk
	654	tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u(i1:i2,j1:j2,jk,Kbb_a) * umask(i1:i2,j1:j2,jk)
	655	END DO
	656	tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2)
	657	END IF
[12377]	658	END IF
	659
[5656]	660	ELSE
[9031]	661
[14086]	662	IF ( l_vremap ) THEN
[12377]	663
[14086]	664	IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
	665
	666	DO jj=j1,j2
	667	DO ji=i1,i2
	668	tabres_child(ji,jj,:) = 0._wp
	669	N_in = mbku_parent(ji,jj)
[14218]	670	N_out = mbku(ji,jj)
	671	IF (N_in * N_out > 0) THEN
	672	zhtot = 0._wp
	673	DO jk=1,N_in
	674	!IF (jk==N_in) THEN
	675	! h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
	676	!ELSE
	677	! h_in(jk) = tabres(ji,jj,jk,m2)
	678	!ENDIF
	679	h_in(jk) = e3u0_parent(ji,jj,jk)
	680	zhtot = zhtot + h_in(jk)
	681	tabin(jk) = tabres(ji,jj,jk,m1)
	682	END DO
	683	!
	684	DO jk=1,N_out
	685	h_out(jk) = e3u(ji,jj,jk,Kbb_a)
	686	END DO
[14086]	687
[14218]	688	! Account for small differences in free-surface
	689	IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
	690	h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
	691	ELSE
	692	h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
	693	ENDIF
[14086]	694
	695	CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
	696	ENDIF
[13286]	697	END DO
	698	END DO
[9031]	699
[15317]	700	ubdiff(i1:i2,j1:j2,1:jpk) = (uu(i1:i2,j1:j2,1:jpk,Kbb_a) - tabres_child(i1:i2,j1:j2,1:jpk))*umask(i1:i2,j1:j2,1:jpk)
[14086]	701	ELSE
	702
[15317]	703	ubdiff(i1:i2,j1:j2,1:jpk) = (uu(i1:i2,j1:j2,1:jpk,Kbb_a) - tabres(i1:i2,j1:j2,1:jpk,1))*umask(i1:i2,j1:j2,1:jpk)
[14086]	704
	705	ENDIF
[6140]	706	!
[5656]	707	DO jk = 1, jpkm1 ! Horizontal slab
	708	! ! ===============
	709
	710	! ! --------
	711	! Horizontal divergence ! div
	712	! ! --------
	713	DO jj = j1,j2
	714	DO ji = i1+1,i2 ! vector opt.
[13216]	715	zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a)
[12377]	716	hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)e3u(ji ,jj,jk,Kbb_a) ubdiff(ji ,jj,jk) &
	717	& -e2u(ji-1,jj)e3u(ji-1,jj,jk,Kbb_a) ubdiff(ji-1,jj,jk) ) * zbtr
[5656]	718	END DO
	719	END DO
	720
	721	DO jj = j1,j2-1
	722	DO ji = i1,i2 ! vector opt.
[13216]	723	zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk)
[9019]	724	rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &
	725	& +e1u(ji,jj ) * ubdiff(ji,jj ,jk) ) * fmask(ji,jj,jk) * zbtr
[5656]	726	END DO
	727	END DO
[6140]	728	END DO
[5656]	729	!
	730	DO jj = j1+1, j2-1
	731	DO ji = i1+1, i2-1 ! vector opt.
	732
	733	IF (.NOT. tabspongedone_u(ji,jj)) THEN
	734	DO jk = 1, jpkm1 ! Horizontal slab
	735	ze2u = rotdiff (ji,jj,jk)
	736	ze1v = hdivdiff(ji,jj,jk)
	737	! horizontal diffusive trends
[12377]	738	zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) &
	739	& + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) &
[13216]	740	& - rn_trelax_dyn * r1_Dt * fspu(ji,jj) * ubdiff(ji,jj,jk)
[5656]	741
	742	! add it to the general momentum trends
[12377]	743	uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua
[5656]	744	END DO
	745	ENDIF
	746
	747	END DO
	748	END DO
	749
	750	tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
	751
	752	jmax = j2-1
[14976]	753	ind1 = jpjglo - ( nn_hls + nbghostcells + 1 ) ! North
[13286]	754	DO jj = mj0(ind1), mj1(ind1)
	755	jmax = MIN(jmax,jj)
	756	END DO
[5656]	757
	758	DO jj = j1+1, jmax
	759	DO ji = i1+1, i2 ! vector opt.
	760
	761	IF (.NOT. tabspongedone_v(ji,jj)) THEN
	762	DO jk = 1, jpkm1 ! Horizontal slab
	763	ze2u = rotdiff (ji,jj,jk)
	764	ze1v = hdivdiff(ji,jj,jk)
	765
	766	! horizontal diffusive trends
[12377]	767	zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) ) &
[9019]	768	+ ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
[5656]	769
	770	! add it to the general momentum trends
[12377]	771	vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) + zva
[5656]	772	END DO
	773	ENDIF
[6140]	774	!
[5656]	775	END DO
	776	END DO
[6140]	777	!
[5656]	778	tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
[6140]	779	!
[5656]	780	ENDIF
[6140]	781	!
[5656]	782	END SUBROUTINE interpun_sponge
	783
[13286]	784
	785	SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
[9031]	786	!!---------------------------------------------
	787	!! * ROUTINE interpvn_sponge *
	788	!!---------------------------------------------
	789	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
	790	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
	791	LOGICAL, INTENT(in) :: before
[6140]	792	!
[9031]	793	INTEGER :: ji, jj, jk, imax
[13286]	794	INTEGER :: ind1
[13216]	795	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot
[9031]	796	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
	797	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
	798	! vertical interpolation:
	799	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
	800	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
	801	REAL(wp), DIMENSION(1:jpk) :: h_out
	802	INTEGER :: N_in, N_out
	803	!!---------------------------------------------
	804
[6140]	805	IF( before ) THEN
[12377]	806	DO jk=k1,k2
[9031]	807	DO jj=j1,j2
	808	DO ji=i1,i2
[15437]	809	tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a) * vmask(ji,jj,jk)
[9031]	810	END DO
	811	END DO
	812	END DO
[12377]	813
[14086]	814	IF ( l_vremap ) THEN
	815
	816	DO jk=k1,k2
	817	DO jj=j1,j2
	818	DO ji=i1,i2
	819	tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_a)
	820	END DO
[12377]	821	END DO
[14086]	822	END DO
	823	! Extrapolate thicknesses in partial bottom cells:
	824	! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
	825	IF (ln_zps) THEN
	826	DO jj=j1,j2
	827	DO ji=i1,i2
	828	jk = mbkv(ji,jj)
	829	tabres(ji,jj,jk,m2) = 0._wp
	830	END DO
	831	END DO
	832	END IF
	833	! Save ssh at last level:
	834	tabres(i1:i2,j1:j2,k2,m2) = 0._wp
	835	IF (.NOT.ln_linssh) THEN
	836	! This vertical sum below should be replaced by the sea-level at V-points (optimization):
	837	DO jk=1,jpk
	838	tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v(i1:i2,j1:j2,jk,Kbb_a) * vmask(i1:i2,j1:j2,jk)
	839	END DO
	840	tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2)
	841	END IF
[12377]	842
[14086]	843	END IF
	844
[5656]	845	ELSE
[9031]	846
[14086]	847	IF ( l_vremap ) THEN
	848	IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
	849	DO jj=j1,j2
	850	DO ji=i1,i2
	851	tabres_child(ji,jj,:) = 0._wp
	852	N_in = mbkv_parent(ji,jj)
[14218]	853	N_out = mbkv(ji,jj)
	854	IF (N_in * N_out > 0) THEN
	855	zhtot = 0._wp
	856	DO jk=1,N_in
	857	!IF (jk==N_in) THEN
	858	! h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
	859	!ELSE
	860	! h_in(jk) = tabres(ji,jj,jk,m2)
	861	!ENDIF
	862	h_in(jk) = e3v0_parent(ji,jj,jk)
	863	zhtot = zhtot + h_in(jk)
	864	tabin(jk) = tabres(ji,jj,jk,m1)
	865	END DO
	866	!
	867	DO jk=1,N_out
	868	h_out(jk) = e3v(ji,jj,jk,Kbb_a)
	869	END DO
[14086]	870
[14218]	871	! Account for small differences in free-surface
	872	IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
	873	h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
	874	ELSE
	875	h_in(1) = h_in(1) - ( sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
	876	ENDIF
[14086]	877
	878	CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
[14218]	879
[14086]	880	ENDIF
[13286]	881	END DO
	882	END DO
[9031]	883
[15317]	884	vbdiff(i1:i2,j1:j2,1:jpk) = (vv(i1:i2,j1:j2,1:jpk,Kbb_a) - tabres_child(i1:i2,j1:j2,1:jpk))*vmask(i1:i2,j1:j2,1:jpk)
[14086]	885	ELSE
	886
[15317]	887	vbdiff(i1:i2,j1:j2,1:jpk) = (vv(i1:i2,j1:j2,1:jpk,Kbb_a) - tabres(i1:i2,j1:j2,1:jpk,1))*vmask(i1:i2,j1:j2,1:jpk)
[14086]	888
	889	ENDIF
[6140]	890	!
[5656]	891	DO jk = 1, jpkm1 ! Horizontal slab
	892	! ! ===============
	893
	894	! ! --------
	895	! Horizontal divergence ! div
	896	! ! --------
	897	DO jj = j1+1,j2
	898	DO ji = i1,i2 ! vector opt.
[13216]	899	zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a)
[12377]	900	hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v(ji,jj ,jk,Kbb_a) * vbdiff(ji,jj ,jk) &
	901	& -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kbb_a) * vbdiff(ji,jj-1,jk) ) * zbtr
[5656]	902	END DO
	903	END DO
	904	DO jj = j1,j2
	905	DO ji = i1,i2-1 ! vector opt.
[13216]	906	zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk)
[5656]	907	rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) &
[6140]	908	& -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr
[5656]	909	END DO
	910	END DO
[6140]	911	END DO
[5656]	912
	913	! ! ===============
	914	!
	915
[9019]	916	imax = i2 - 1
[14976]	917	ind1 = jpiglo - ( nn_hls + nbghostcells + 1 ) ! East
[13286]	918	DO ji = mi0(ind1), mi1(ind1)
	919	imax = MIN(imax,ji)
	920	END DO
	921
[5656]	922	DO jj = j1+1, j2
	923	DO ji = i1+1, imax ! vector opt.
[6140]	924	IF( .NOT. tabspongedone_u(ji,jj) ) THEN
	925	DO jk = 1, jpkm1
[13216]	926	uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) &
[12377]	927	& - ( rotdiff (ji ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) &
[6140]	928	& + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj ,jk)) * r1_e1u(ji,jj)
[5656]	929	END DO
	930	ENDIF
	931	END DO
	932	END DO
[6140]	933	!
[5656]	934	tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
[6140]	935	!
[5656]	936	DO jj = j1+1, j2-1
	937	DO ji = i1+1, i2-1 ! vector opt.
[6140]	938	IF( .NOT. tabspongedone_v(ji,jj) ) THEN
	939	DO jk = 1, jpkm1
[13216]	940	vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) &
[12377]	941	& + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) ) &
[13216]	942	& + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj) &
	943	& - rn_trelax_dyn * r1_Dt * fspv(ji,jj) * vbdiff(ji,jj,jk)
[5656]	944	END DO
	945	ENDIF
	946	END DO
	947	END DO
	948	tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
	949	ENDIF
[6140]	950	!
[5656]	951	END SUBROUTINE interpvn_sponge
	952
[14086]	953	SUBROUTINE interpunb_sponge(tabres,i1,i2,j1,j2, before)
	954	!!---------------------------------------------
	955	!! * ROUTINE interpunb_sponge *
	956	!!---------------------------------------------
	957	INTEGER, INTENT(in) :: i1,i2,j1,j2
	958	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
	959	LOGICAL, INTENT(in) :: before
	960
	961	INTEGER :: ji, jj, ind1, jmax
	962	! sponge parameters
	963	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
	964	REAL(wp), DIMENSION(i1:i2,j1:j2) :: ubdiff
	965	REAL(wp), DIMENSION(i1:i2,j1:j2) :: rotdiff, hdivdiff
	966	!!---------------------------------------------
	967	!
	968	IF( before ) THEN
	969	DO jj=j1,j2
	970	DO ji=i1,i2
	971	tabres(ji,jj) = uu_b(ji,jj,Kmm_a)
	972	END DO
	973	END DO
	974
	975	ELSE
	976
	977	ubdiff(i1:i2,j1:j2) = (uu_b(i1:i2,j1:j2,Kmm_a) - tabres(i1:i2,j1:j2))*umask(i1:i2,j1:j2,1)
	978	!
	979	! ! --------
	980	! Horizontal divergence ! div
	981	! ! --------
	982	DO jj = j1,j2
	983	DO ji = i1+1,i2 ! vector opt.
	984	zbtr = rn_sponge_dyn * r1_Dt * fspt_2d(ji,jj) * r1_ht_0(ji,jj)
	985	hdivdiff(ji,jj) = ( e2u(ji ,jj)hu(ji ,jj,Kbb_a) ubdiff(ji ,jj) &
	986	&-e2u(ji-1,jj)hu(ji-1,jj,Kbb_a) ubdiff(ji-1,jj) ) * zbtr
	987	END DO
	988	END DO
	989
	990	DO jj = j1,j2-1
	991	DO ji = i1,i2 ! vector opt.
	992	zbtr = rn_sponge_dyn * r1_Dt * fspf_2d(ji,jj) * hf_0(ji,jj)
	993	rotdiff(ji,jj) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1) &
	994	& +e1u(ji,jj ) * ubdiff(ji,jj ) ) * fmask(ji,jj,1) * zbtr
	995	END DO
	996	END DO
	997	!
	998	DO jj = j1+1, j2-1
	999	DO ji = i1+1, i2-1 ! vector opt.
	1000	IF (.NOT. tabspongedone_u(ji,jj)) THEN
	1001	ze2u = rotdiff (ji,jj)
	1002	ze1v = hdivdiff(ji,jj)
	1003	! horizontal diffusive trends
	1004	zua = - ( ze2u - rotdiff (ji,jj-1) ) * r1_e2u(ji,jj) * r1_hu(ji,jj,Kmm_a) &
	1005	& + ( hdivdiff(ji+1,jj) - ze1v ) * r1_e1u(ji,jj) &
	1006	& - rn_trelax_dyn * r1_Dt * fspu_2d(ji,jj) * ubdiff(ji,jj)
	1007
	1008	! add it to the general momentum trends
	1009	uu(ji,jj,:,Krhs_a) = uu(ji,jj,:,Krhs_a) + zua
	1010	ENDIF
	1011	END DO
	1012	END DO
	1013
	1014	tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
	1015
	1016	jmax = j2-1
[14976]	1017	ind1 = jpjglo - ( nn_hls + nbghostcells + 1 ) ! North
[14086]	1018	DO jj = mj0(ind1), mj1(ind1)
	1019	jmax = MIN(jmax,jj)
	1020	END DO
	1021
	1022	DO jj = j1+1, jmax
	1023	DO ji = i1+1, i2 ! vector opt.
	1024	IF (.NOT. tabspongedone_v(ji,jj)) THEN
	1025	ze2u = rotdiff (ji,jj)
	1026	ze1v = hdivdiff(ji,jj)
	1027	zva = + ( ze2u - rotdiff (ji-1,jj) ) * r1_e1v(ji,jj) * r1_hv(ji,jj,Kmm_a) &
	1028	+ ( hdivdiff(ji,jj+1) - ze1v ) * r1_e2v(ji,jj)
	1029	vv(ji,jj,:,Krhs_a) = vv(ji,jj,:,Krhs_a) + zva
	1030	ENDIF
	1031	END DO
	1032	END DO
	1033	!
	1034	tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
	1035	!
	1036	ENDIF
	1037	!
	1038	END SUBROUTINE interpunb_sponge
	1039
	1040
	1041	SUBROUTINE interpvnb_sponge(tabres,i1,i2,j1,j2, before)
	1042	!!---------------------------------------------
	1043	!! * ROUTINE interpvnb_sponge *
	1044	!!---------------------------------------------
	1045	INTEGER, INTENT(in) :: i1,i2,j1,j2
	1046	REAL(wp), DIMENSION(i1:i2,j1:j2), INTENT(inout) :: tabres
	1047	LOGICAL, INTENT(in) :: before
	1048	!
	1049	INTEGER :: ji, jj, ind1, imax
	1050	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
	1051	REAL(wp), DIMENSION(i1:i2,j1:j2) :: vbdiff
	1052	REAL(wp), DIMENSION(i1:i2,j1:j2) :: rotdiff, hdivdiff
	1053	!!---------------------------------------------
	1054
	1055	IF( before ) THEN
	1056	DO jj=j1,j2
	1057	DO ji=i1,i2
	1058	tabres(ji,jj) = vv_b(ji,jj,Kmm_a)
	1059	END DO
	1060	END DO
	1061	ELSE
	1062	vbdiff(i1:i2,j1:j2) = (vv_b(i1:i2,j1:j2,Kmm_a) - tabres(i1:i2,j1:j2))*vmask(i1:i2,j1:j2,1)
	1063	! ! --------
	1064	! Horizontal divergence ! div
	1065	! ! --------
	1066	DO jj = j1+1,j2
	1067	DO ji = i1,i2 ! vector opt.
	1068	zbtr = rn_sponge_dyn * r1_Dt * fspt_2d(ji,jj) * r1_ht_0(ji,jj)
	1069	hdivdiff(ji,jj) = ( e1v(ji,jj ) * hv(ji,jj ,Kbb_a) * vbdiff(ji,jj ) &
	1070	& -e1v(ji,jj-1) * hv(ji,jj-1,Kbb_a) * vbdiff(ji,jj-1) ) * zbtr
	1071	END DO
	1072	END DO
	1073	DO jj = j1,j2
	1074	DO ji = i1,i2-1 ! vector opt.
	1075	zbtr = rn_sponge_dyn * r1_Dt * fspf_2d(ji,jj) * hf_0(ji,jj)
	1076	rotdiff(ji,jj) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj) &
	1077	& -e2v(ji ,jj) * vbdiff(ji ,jj) ) * fmask(ji,jj,1) * zbtr
	1078	END DO
	1079	END DO
	1080	! ! ===============
	1081	!
	1082
	1083	imax = i2 - 1
[14976]	1084	ind1 = jpiglo - ( nn_hls + nbghostcells + 1 ) ! East
[14086]	1085	DO ji = mi0(ind1), mi1(ind1)
	1086	imax = MIN(imax,ji)
	1087	END DO
	1088
	1089	DO jj = j1+1, j2
	1090	DO ji = i1+1, imax ! vector opt.
	1091	IF( .NOT. tabspongedone_u(ji,jj) ) THEN
	1092	zua = - ( rotdiff (ji ,jj) - rotdiff (ji,jj-1)) * r1_e2u(ji,jj) * r1_hu(ji,jj,Kmm_a) &
	1093	& + ( hdivdiff(ji+1,jj) - hdivdiff(ji,jj )) * r1_e1u(ji,jj)
	1094	uu(ji,jj,:,Krhs_a) = uu(ji,jj,:,Krhs_a) + zua
	1095	ENDIF
	1096	END DO
	1097	END DO
	1098	!
	1099	tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
	1100	!
	1101	DO jj = j1+1, j2-1
	1102	DO ji = i1+1, i2-1 ! vector opt.
	1103	IF( .NOT. tabspongedone_v(ji,jj) ) THEN
	1104	zva = ( rotdiff (ji,jj ) - rotdiff (ji-1,jj) ) * r1_e1v(ji,jj) *r1_hv(ji,jj,Kmm_a) &
	1105	& + ( hdivdiff(ji,jj+1) - hdivdiff(ji ,jj) ) * r1_e2v(ji,jj) &
	1106	& - rn_trelax_dyn * r1_Dt * fspv_2d(ji,jj) * vbdiff(ji,jj)
	1107	vv(ji,jj,:,Krhs_a) = vv(ji,jj,:,Krhs_a) + zva
	1108	ENDIF
	1109	END DO
	1110	END DO
	1111	tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
	1112	ENDIF
	1113	!
	1114	END SUBROUTINE interpvnb_sponge
	1115
	1116
[390]	1117	#else
[9019]	1118	!!----------------------------------------------------------------------
	1119	!! Empty module no AGRIF zoom
	1120	!!----------------------------------------------------------------------
[636]	1121	CONTAINS
[9570]	1122	SUBROUTINE agrif_oce_sponge_empty
	1123	WRITE(,) 'agrif_oce_sponge : You should not have seen this print! error?'
	1124	END SUBROUTINE agrif_oce_sponge_empty
[390]	1125	#endif
[636]	1126
[6140]	1127	!!======================================================================
[9570]	1128	END MODULE agrif_oce_sponge

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

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

Context Navigation

source: NEMO/trunk/src/NST/agrif_oce_sponge.F90 @ 15437

Download in other formats: