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

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

agrif_oce_sponge.F90 in NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13312_AGRIF-03-04_jchanut_vinterp_tstep/src/NST/agrif_oce_sponge.F90 @ 13521

Last change on this file since 13521 was 13498, checked in by jchanut, 4 years ago
#2222 1) Define agrif sponge at f-points 2) Correct bathymetry check with vertical refinement just over needed points
Property svn:keywords set to `Id`
File size: 31.3 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 *
[9019]	6	!! AGRIF: sponge package for the ocean dynamics (OPA)
[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
[5656]	30	PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn
	31	PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge
[390]	32
[12377]	33	!! * Substitutions
	34	# include "do_loop_substitute.h90"
[1156]	35	!!----------------------------------------------------------------------
[9598]	36	!! NEMO/NST 4.0 , NEMO Consortium (2018)
[1156]	37	!! $Id$
[10068]	38	!! Software governed by the CeCILL license (see ./LICENSE)
[1156]	39	!!----------------------------------------------------------------------
[5656]	40	CONTAINS
[636]	41
[782]	42	SUBROUTINE Agrif_Sponge_Tra
[9019]	43	!!----------------------------------------------------------------------
	44	!! * ROUTINE Agrif_Sponge_Tra *
	45	!!----------------------------------------------------------------------
	46	REAL(wp) :: zcoef ! local scalar
[9056]	47
[9019]	48	!!----------------------------------------------------------------------
[6140]	49	!
[390]	50	#if defined SPONGE
[9031]	51	!! Assume persistence:
[9056]	52	zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
[9031]	53
[5656]	54	CALL Agrif_Sponge
[9019]	55	Agrif_SpecialValue = 0._wp
[390]	56	Agrif_UseSpecialValue = .TRUE.
[13351]	57	l_vremap = ln_vert_remap
[9019]	58	tabspongedone_tsn = .FALSE.
	59	!
[13351]	60	CALL Agrif_Bc_Variable( ts_sponge_id, calledweight=zcoef, procname=interptsn_sponge )
[9019]	61	!
[5656]	62	Agrif_UseSpecialValue = .FALSE.
[13337]	63	l_vremap = .FALSE.
[390]	64	#endif
[6140]	65	!
[12377]	66	CALL iom_put( 'agrif_spu', fspu(:,:))
	67	CALL iom_put( 'agrif_spv', fspv(:,:))
	68	!
[636]	69	END SUBROUTINE Agrif_Sponge_Tra
[390]	70
[6140]	71
[782]	72	SUBROUTINE Agrif_Sponge_dyn
[9019]	73	!!----------------------------------------------------------------------
	74	!! * ROUTINE Agrif_Sponge_dyn *
	75	!!----------------------------------------------------------------------
	76	REAL(wp) :: zcoef ! local scalar
	77	!!----------------------------------------------------------------------
	78	!
[390]	79	#if defined SPONGE
[9056]	80	zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
[9031]	81
[13286]	82	Agrif_SpecialValue = 0._wp
[782]	83	Agrif_UseSpecialValue = ln_spc_dyn
[13351]	84	l_vremap = ln_vert_remap
[13286]	85	use_sign_north = .TRUE.
	86	sign_north = -1._wp
[9019]	87	!
[5656]	88	tabspongedone_u = .FALSE.
	89	tabspongedone_v = .FALSE.
[9019]	90	CALL Agrif_Bc_Variable( un_sponge_id, calledweight=zcoef, procname=interpun_sponge )
	91	!
[5656]	92	tabspongedone_u = .FALSE.
	93	tabspongedone_v = .FALSE.
[9019]	94	CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge )
	95	!
[390]	96	Agrif_UseSpecialValue = .FALSE.
[13286]	97	use_sign_north = .FALSE.
[13337]	98	l_vremap = .FALSE.
[390]	99	#endif
[6140]	100	!
[12377]	101	CALL iom_put( 'agrif_spt', fspt(:,:))
	102	CALL iom_put( 'agrif_spf', fspf(:,:))
	103	!
[636]	104	END SUBROUTINE Agrif_Sponge_dyn
[390]	105
[6140]	106
[3680]	107	SUBROUTINE Agrif_Sponge
[9019]	108	!!----------------------------------------------------------------------
	109	!! * ROUTINE Agrif_Sponge *
	110	!!----------------------------------------------------------------------
	111	INTEGER :: ji, jj, ind1, ind2
[12377]	112	INTEGER :: ispongearea, jspongearea
	113	REAL(wp) :: z1_ispongearea, z1_jspongearea
[9019]	114	REAL(wp), DIMENSION(jpi,jpj) :: ztabramp
	115	!!----------------------------------------------------------------------
	116	!
[12377]	117	! Sponge 1d example with:
	118	! iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2
	119	!
	120	!coarse : U T U T U T U
	121	!\| \| \| \| \|
	122	!fine : t u t u t u t u t u t u t u t u t u t u t
[13498]	123	!sponge val:0 1 1 1 1 5/6 4/6 3/6 2/6 1/6 0
[12377]	124	! \| ghost \| <-- sponge area -- > \|
	125	! \| points \| \|
	126	! \|--> dynamical interface
	127
[3680]	128	#if defined SPONGE \|\| defined SPONGE_TOP
[4153]	129	IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
[13498]	130	! Define ramp from boundaries towards domain interior at F-points
[4153]	131	! Store it in ztabramp
[3680]	132
[13286]	133	ispongearea = nn_sponge_len * Agrif_irhox()
[13498]	134	z1_ispongearea = 1._wp / REAL( MAX(ispongearea,1), wp )
[13286]	135	jspongearea = nn_sponge_len * Agrif_irhoy()
[13498]	136	z1_jspongearea = 1._wp / REAL( MAX(jspongearea,1), wp )
[9019]	137
[5656]	138	ztabramp(:,:) = 0._wp
[4153]	139
[12377]	140	! Trick to remove sponge in 2DV domains:
	141	IF ( nbcellsx <= 3 ) ispongearea = -1
	142	IF ( nbcellsy <= 3 ) jspongearea = -1
	143
[13286]	144	IF( lk_west ) THEN ! --- West --- !
	145	ind1 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells
	146	ind2 = nn_hls + 1 + nbghostcells + ispongearea
[13216]	147	DO ji = mi0(ind1), mi1(ind2)
	148	DO jj = 1, jpj
[13337]	149	ztabramp(ji,jj) = REAL(ind2 - mig(ji), wp) * z1_ispongearea
[13216]	150	END DO
[13286]	151	END DO
[13216]	152	! ghost cells:
	153	ind1 = 1
[13286]	154	ind2 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells
[13216]	155	DO ji = mi0(ind1), mi1(ind2)
	156	DO jj = 1, jpj
[13337]	157	ztabramp(ji,jj) = 1._wp
[13216]	158	END DO
[12377]	159	END DO
[13216]	160	ENDIF
[13286]	161	IF( lk_east ) THEN ! --- East --- !
[13498]	162	ind1 = jpiglo - ( nn_hls + nbghostcells ) - ispongearea - 1
	163	ind2 = jpiglo - ( nn_hls + nbghostcells ) - 1 ! halo + land + nbghostcells - 1
[13216]	164	DO ji = mi0(ind1), mi1(ind2)
	165	DO jj = 1, jpj
[13337]	166	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mig(ji) - ind1, wp) * z1_ispongearea )
[13286]	167	END DO
[13216]	168	END DO
	169	! ghost cells:
[13498]	170	ind1 = jpiglo - ( nn_hls + nbghostcells ) - 1 ! halo + land + nbghostcells - 1
	171	ind2 = jpiglo - 1
[13216]	172	DO ji = mi0(ind1), mi1(ind2)
	173	DO jj = 1, jpj
[13337]	174	ztabramp(ji,jj) = 1._wp
[13286]	175	END DO
[13216]	176	END DO
[13286]	177	ENDIF
	178	IF( lk_south ) THEN ! --- South --- !
[13498]	179	ind1 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells
[13286]	180	ind2 = nn_hls + 1 + nbghostcells + jspongearea
[13216]	181	DO jj = mj0(ind1), mj1(ind2)
	182	DO ji = 1, jpi
[13337]	183	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(ind2 - mjg(jj), wp) * z1_jspongearea )
[13216]	184	END DO
[12377]	185	END DO
[13216]	186	! ghost cells:
	187	ind1 = 1
[13498]	188	ind2 = nn_hls + 1 + nbghostcells ! halo + land + nbghostcells
[13216]	189	DO jj = mj0(ind1), mj1(ind2)
	190	DO ji = 1, jpi
[13337]	191	ztabramp(ji,jj) = 1._wp
[13216]	192	END DO
[12377]	193	END DO
[13216]	194	ENDIF
[13286]	195	IF( lk_north ) THEN ! --- North --- !
[13498]	196	ind1 = jpjglo - ( nn_hls + nbghostcells ) - jspongearea - 1
	197	ind2 = jpjglo - ( nn_hls + nbghostcells ) - 1 ! halo + land + nbghostcells - 1
[13216]	198	DO jj = mj0(ind1), mj1(ind2)
	199	DO ji = 1, jpi
[13337]	200	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL(mjg(jj) - ind1, wp) * z1_jspongearea )
[13216]	201	END DO
[12377]	202	END DO
[13216]	203	! ghost cells:
[13286]	204	ind1 = jpjglo - ( nn_hls + nbghostcells ) ! halo + land + nbghostcells - 1
[13216]	205	ind2 = jpjglo
	206	DO jj = mj0(ind1), mj1(ind2)
	207	DO ji = 1, jpi
[13337]	208	ztabramp(ji,jj) = 1._wp
[13216]	209	END DO
[12377]	210	END DO
[13216]	211	ENDIF
[13286]	212	!
[4153]	213	ENDIF
	214
[3680]	215	! Tracers
	216	IF( .NOT. spongedoneT ) THEN
[12377]	217	fspu(:,:) = 0._wp
	218	fspv(:,:) = 0._wp
[13295]	219	DO_2D( 0, 0, 0, 0 )
[13498]	220	fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji,jj-1) ) * ssumask(ji,jj)
	221	fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji-1,jj) ) * ssvmask(ji,jj)
[12377]	222	END_2D
[3680]	223	ENDIF
	224
	225	! Dynamics
	226	IF( .NOT. spongedoneU ) THEN
[12377]	227	fspt(:,:) = 0._wp
	228	fspf(:,:) = 0._wp
[13295]	229	DO_2D( 0, 0, 0, 0 )
[13498]	230	fspt(ji,jj) = 0.25_wp * ( ztabramp(ji ,jj ) + ztabramp(ji-1,jj ) &
	231	& +ztabramp(ji ,jj-1) + ztabramp(ji-1,jj-1) ) * ssmask(ji,jj)
	232	fspf(ji,jj) = ztabramp(ji,jj) * ssvmask(ji,jj) * ssvmask(ji,jj+1)
[12377]	233	END_2D
[13286]	234	ENDIF
	235
	236	IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN
	237	CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp, fspt, 'T', 1._wp, fspf, 'F', 1._wp )
	238	spongedoneT = .TRUE.
[3680]	239	spongedoneU = .TRUE.
	240	ENDIF
[13286]	241	IF( .NOT. spongedoneT ) THEN
	242	CALL lbc_lnk_multi( 'agrif_Sponge', fspu, 'U', 1._wp, fspv, 'V', 1._wp )
	243	spongedoneT = .TRUE.
	244	ENDIF
	245	IF( .NOT. spongedoneT .AND. .NOT. spongedoneU ) THEN
	246	CALL lbc_lnk_multi( 'agrif_Sponge', fspt, 'T', 1._wp, fspf, 'F', 1._wp )
	247	spongedoneU = .TRUE.
	248	ENDIF
[13498]	249	!
[12377]	250	! Remove vertical interpolation where not needed:
[13498]	251	! (A null value in mbkx arrays does the job)
	252	WHERE (fspu(:,:) == 0._wp) mbku_parent(:,:) = 0
	253	WHERE (fspv(:,:) == 0._wp) mbkv_parent(:,:) = 0
	254	WHERE (fspt(:,:) == 0._wp) mbkt_parent(:,:) = 0
[3680]	255	!
[12377]	256	#endif
	257	!
[3680]	258	END SUBROUTINE Agrif_Sponge
	259
[13286]	260
[9019]	261	SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
	262	!!----------------------------------------------------------------------
	263	!! * ROUTINE interptsn_sponge *
	264	!!----------------------------------------------------------------------
	265	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
	266	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
	267	LOGICAL , INTENT(in ) :: before
[6140]	268	!
[5656]	269	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	270	INTEGER :: iku, ikv
[13216]	271	REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot
[9031]	272	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv
	273	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
	274	! vertical interpolation:
	275	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
	276	REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
	277	REAL(wp), DIMENSION(k1:k2) :: h_in
	278	REAL(wp), DIMENSION(1:jpk) :: h_out
	279	INTEGER :: N_in, N_out
[9019]	280	!!----------------------------------------------------------------------
[5656]	281	!
[6140]	282	IF( before ) THEN
[9031]	283	DO jn = 1, jpts
	284	DO jk=k1,k2
	285	DO jj=j1,j2
	286	DO ji=i1,i2
[12377]	287	tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a)
[9031]	288	END DO
	289	END DO
	290	END DO
	291	END DO
	292
[13337]	293	IF ( l_vremap ) THEN
	294
	295	! Interpolate thicknesses
	296	! Warning: these are masked, hence extrapolated prior interpolation.
	297	DO jk=k1,k2
	298	DO jj=j1,j2
	299	DO ji=i1,i2
	300	tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a)
	301	END DO
[12377]	302	END DO
	303	END DO
	304
[13337]	305	! Extrapolate thicknesses in partial bottom cells:
	306	! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
	307	IF (ln_zps) THEN
	308	DO jj=j1,j2
	309	DO ji=i1,i2
	310	jk = mbkt(ji,jj)
	311	tabres(ji,jj,jk,jpts+1) = 0._wp
	312	END DO
	313	END DO
	314	END IF
	315
	316	! Save ssh at last level:
	317	IF (.NOT.ln_linssh) THEN
	318	tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1)
	319	ELSE
	320	tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp
	321	END IF
[12377]	322	END IF
[9031]	323
[5656]	324	ELSE
[6140]	325	!
[13337]	326	IF ( l_vremap ) THEN
[12377]	327
[13337]	328	IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp
[12377]	329
[13337]	330	DO jj=j1,j2
	331	DO ji=i1,i2
	332	tabres_child(ji,jj,:,:) = 0._wp
	333	N_in = mbkt_parent(ji,jj)
	334	zhtot = 0._wp
	335	DO jk=1,N_in !k2 = jpk of parent grid
	336	IF (jk==N_in) THEN
	337	h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot
	338	ELSE
	339	h_in(jk) = tabres(ji,jj,jk,n2)
	340	END IF
	341	zhtot = zhtot + h_in(jk)
	342	tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)
	343	END DO
	344	N_out = 0
	345	DO jk=1,jpk ! jpk of child grid
	346	IF (tmask(ji,jj,jk) == 0) EXIT
	347	N_out = N_out + 1
	348	h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above
	349	END DO
	350
	351	! Account for small differences in the free-surface
	352	IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
	353	h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
[12377]	354	ELSE
[13337]	355	h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
	356	END IF
	357	IF (N_in*N_out > 0) THEN
	358	CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
[12377]	359	ENDIF
[9031]	360	END DO
[13337]	361	END DO
	362
	363	DO jj=j1,j2
	364	DO ji=i1,i2
	365	DO jk=1,jpkm1
	366	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)
	367	END DO
[13286]	368	END DO
	369	END DO
[9031]	370
[13337]	371	ELSE
	372
	373	DO jj=j1,j2
	374	DO ji=i1,i2
	375	DO jk=1,jpkm1
	376	tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk)
	377	END DO
[13286]	378	END DO
	379	END DO
[13337]	380	END IF
[9031]	381
[5656]	382	DO jn = 1, jpts
	383	DO jk = 1, jpkm1
[9806]	384	ztu(i1:i2,j1:j2,jk) = 0._wp
	385	DO jj = j1,j2
[5656]	386	DO ji = i1,i2-1
[13498]	387	zabe1 = fspu(ji,jj) * rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
[9806]	388	ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) )
	389	END DO
	390	END DO
	391	ztv(i1:i2,j1:j2,jk) = 0._wp
	392	DO ji = i1,i2
	393	DO jj = j1,j2-1
[13498]	394	zabe2 = fspv(ji,jj) * rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
[5656]	395	ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
[6140]	396	END DO
	397	END DO
	398	!
[5656]	399	IF( ln_zps ) THEN ! set gradient at partial step level
[9806]	400	DO jj = j1,j2
	401	DO ji = i1,i2
[5656]	402	! last level
	403	iku = mbku(ji,jj)
	404	ikv = mbkv(ji,jj)
[6140]	405	IF( iku == jk ) ztu(ji,jj,jk) = 0._wp
	406	IF( ikv == jk ) ztv(ji,jj,jk) = 0._wp
[5656]	407	END DO
	408	END DO
	409	ENDIF
[6140]	410	END DO
	411	!
[5656]	412	DO jk = 1, jpkm1
	413	DO jj = j1+1,j2-1
	414	DO ji = i1+1,i2-1
	415	IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
[12377]	416	zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
[5656]	417	! horizontal diffusive trends
[12377]	418	ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) &
[13216]	419	& - rn_trelax_tra * r1_Dt * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn)
[5656]	420	! add it to the general tracer trends
[12377]	421	ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa
[5656]	422	ENDIF
[6140]	423	END DO
	424	END DO
	425	END DO
	426	!
	427	END DO
	428	!
[5656]	429	tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
[6140]	430	!
[5656]	431	ENDIF
[6140]	432	!
[5656]	433	END SUBROUTINE interptsn_sponge
	434
[13286]	435
[9031]	436	SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
	437	!!---------------------------------------------
	438	!! * ROUTINE interpun_sponge *
	439	!!---------------------------------------------
	440	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
	441	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
	442	LOGICAL, INTENT(in) :: before
[6140]	443
[13286]	444	INTEGER :: ji,jj,jk,jmax
	445	INTEGER :: ind1
[5656]	446	! sponge parameters
[13216]	447	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot
[9031]	448	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
	449	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
	450	! vertical interpolation:
	451	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
	452	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
	453	REAL(wp), DIMENSION(1:jpk) :: h_out
[12377]	454	INTEGER ::N_in, N_out
[9031]	455	!!---------------------------------------------
[5656]	456	!
[6140]	457	IF( before ) THEN
[12377]	458	DO jk=k1,k2
[9031]	459	DO jj=j1,j2
	460	DO ji=i1,i2
[12377]	461	tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a)
[9031]	462	END DO
	463	END DO
	464	END DO
	465
[13337]	466	IF ( l_vremap ) THEN
	467
	468	DO jk=k1,k2
	469	DO jj=j1,j2
	470	DO ji=i1,i2
	471	tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk)
	472	END DO
[12377]	473	END DO
[13337]	474	END DO
	475
	476	! Extrapolate thicknesses in partial bottom cells:
	477	! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
	478	IF (ln_zps) THEN
	479	DO jj=j1,j2
	480	DO ji=i1,i2
	481	jk = mbku(ji,jj)
	482	tabres(ji,jj,jk,m2) = 0._wp
	483	END DO
	484	END DO
	485	END IF
	486	! Save ssh at last level:
	487	tabres(i1:i2,j1:j2,k2,m2) = 0._wp
	488	IF (.NOT.ln_linssh) THEN
	489	! This vertical sum below should be replaced by the sea-level at U-points (optimization):
	490	DO jk=1,jpk
	491	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)
	492	END DO
	493	tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2)
	494	END IF
[12377]	495	END IF
	496
[5656]	497	ELSE
[9031]	498
[13337]	499	IF ( l_vremap ) THEN
[12377]	500
[13337]	501	IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
	502
	503	DO jj=j1,j2
	504	DO ji=i1,i2
	505	tabres_child(ji,jj,:) = 0._wp
	506	N_in = mbku_parent(ji,jj)
	507	zhtot = 0._wp
	508	DO jk=1,N_in
	509	IF (jk==N_in) THEN
	510	h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
	511	ELSE
	512	h_in(jk) = tabres(ji,jj,jk,m2)
	513	ENDIF
	514	zhtot = zhtot + h_in(jk)
	515	tabin(jk) = tabres(ji,jj,jk,m1)
	516	END DO
	517	!
	518	N_out = 0
	519	DO jk=1,jpk
	520	IF (umask(ji,jj,jk) == 0) EXIT
	521	N_out = N_out + 1
	522	h_out(N_out) = e3u(ji,jj,jk,Kbb_a)
	523	END DO
	524
	525	! Account for small differences in free-surface
	526	IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
	527	h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
[12377]	528	ELSE
[13337]	529	h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
[12377]	530	ENDIF
[13337]	531
	532	IF (N_in * N_out > 0) THEN
	533	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)
	534	ENDIF
[13286]	535	END DO
	536	END DO
[9031]	537
[13337]	538	ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
	539
	540	ELSE
	541
	542	ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
	543
	544	ENDIF
[6140]	545	!
[5656]	546	DO jk = 1, jpkm1 ! Horizontal slab
	547	! ! ===============
	548
	549	! ! --------
	550	! Horizontal divergence ! div
	551	! ! --------
	552	DO jj = j1,j2
	553	DO ji = i1+1,i2 ! vector opt.
[13216]	554	zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a)
[12377]	555	hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)e3u(ji ,jj,jk,Kbb_a) ubdiff(ji ,jj,jk) &
	556	& -e2u(ji-1,jj)e3u(ji-1,jj,jk,Kbb_a) ubdiff(ji-1,jj,jk) ) * zbtr
[5656]	557	END DO
	558	END DO
	559
	560	DO jj = j1,j2-1
	561	DO ji = i1,i2 ! vector opt.
[13216]	562	zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk)
[9019]	563	rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &
	564	& +e1u(ji,jj ) * ubdiff(ji,jj ,jk) ) * fmask(ji,jj,jk) * zbtr
[5656]	565	END DO
	566	END DO
[6140]	567	END DO
[5656]	568	!
	569	DO jj = j1+1, j2-1
	570	DO ji = i1+1, i2-1 ! vector opt.
	571
	572	IF (.NOT. tabspongedone_u(ji,jj)) THEN
	573	DO jk = 1, jpkm1 ! Horizontal slab
	574	ze2u = rotdiff (ji,jj,jk)
	575	ze1v = hdivdiff(ji,jj,jk)
	576	! horizontal diffusive trends
[12377]	577	zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) &
	578	& + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) &
[13216]	579	& - rn_trelax_dyn * r1_Dt * fspu(ji,jj) * ubdiff(ji,jj,jk)
[5656]	580
	581	! add it to the general momentum trends
[12377]	582	uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua
[5656]	583	END DO
	584	ENDIF
	585
	586	END DO
	587	END DO
	588
	589	tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
	590
	591	jmax = j2-1
[13286]	592	ind1 = jpjglo - ( nn_hls + nbghostcells + 2 ) ! North
	593	DO jj = mj0(ind1), mj1(ind1)
	594	jmax = MIN(jmax,jj)
	595	END DO
[5656]	596
	597	DO jj = j1+1, jmax
	598	DO ji = i1+1, i2 ! vector opt.
	599
	600	IF (.NOT. tabspongedone_v(ji,jj)) THEN
	601	DO jk = 1, jpkm1 ! Horizontal slab
	602	ze2u = rotdiff (ji,jj,jk)
	603	ze1v = hdivdiff(ji,jj,jk)
	604
	605	! horizontal diffusive trends
[12377]	606	zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) ) &
[9019]	607	+ ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
[5656]	608
	609	! add it to the general momentum trends
[12377]	610	vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) + zva
[5656]	611	END DO
	612	ENDIF
[6140]	613	!
[5656]	614	END DO
	615	END DO
[6140]	616	!
[5656]	617	tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
[6140]	618	!
[5656]	619	ENDIF
[6140]	620	!
[5656]	621	END SUBROUTINE interpun_sponge
	622
[13286]	623
	624	SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
[9031]	625	!!---------------------------------------------
	626	!! * ROUTINE interpvn_sponge *
	627	!!---------------------------------------------
	628	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
	629	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
	630	LOGICAL, INTENT(in) :: before
[6140]	631	!
[9031]	632	INTEGER :: ji, jj, jk, imax
[13286]	633	INTEGER :: ind1
[13216]	634	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot
[9031]	635	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
	636	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
	637	! vertical interpolation:
	638	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
	639	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
	640	REAL(wp), DIMENSION(1:jpk) :: h_out
	641	INTEGER :: N_in, N_out
	642	!!---------------------------------------------
	643
[6140]	644	IF( before ) THEN
[12377]	645	DO jk=k1,k2
[9031]	646	DO jj=j1,j2
	647	DO ji=i1,i2
[12377]	648	tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a)
[9031]	649	END DO
	650	END DO
	651	END DO
[12377]	652
[13337]	653	IF ( l_vremap ) THEN
	654
	655	DO jk=k1,k2
	656	DO jj=j1,j2
	657	DO ji=i1,i2
	658	tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_a)
	659	END DO
[12377]	660	END DO
[13337]	661	END DO
	662	! Extrapolate thicknesses in partial bottom cells:
	663	! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
	664	IF (ln_zps) THEN
	665	DO jj=j1,j2
	666	DO ji=i1,i2
	667	jk = mbkv(ji,jj)
	668	tabres(ji,jj,jk,m2) = 0._wp
	669	END DO
	670	END DO
	671	END IF
	672	! Save ssh at last level:
	673	tabres(i1:i2,j1:j2,k2,m2) = 0._wp
	674	IF (.NOT.ln_linssh) THEN
	675	! This vertical sum below should be replaced by the sea-level at V-points (optimization):
	676	DO jk=1,jpk
	677	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)
	678	END DO
	679	tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2)
	680	END IF
[12377]	681
[13337]	682	END IF
	683
[5656]	684	ELSE
[9031]	685
[13337]	686	IF ( l_vremap ) THEN
	687	IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
	688	DO jj=j1,j2
	689	DO ji=i1,i2
	690	tabres_child(ji,jj,:) = 0._wp
	691	N_in = mbkv_parent(ji,jj)
	692	zhtot = 0._wp
	693	DO jk=1,N_in
	694	IF (jk==N_in) THEN
	695	h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
	696	ELSE
	697	h_in(jk) = tabres(ji,jj,jk,m2)
	698	ENDIF
	699	zhtot = zhtot + h_in(jk)
	700	tabin(jk) = tabres(ji,jj,jk,m1)
	701	END DO
	702	!
	703	N_out = 0
	704	DO jk=1,jpk
	705	IF (vmask(ji,jj,jk) == 0) EXIT
	706	N_out = N_out + 1
	707	h_out(N_out) = e3v(ji,jj,jk,Kbb_a)
	708	END DO
	709
	710	! Account for small differences in free-surface
	711	IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
	712	h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
[12377]	713	ELSE
[13337]	714	h_in(1) = h_in(1) - ( sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
[12377]	715	ENDIF
[13337]	716
	717	IF (N_in * N_out > 0) THEN
	718	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)
	719	ENDIF
[13286]	720	END DO
	721	END DO
[9031]	722
[13337]	723	vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)
	724
	725	ELSE
	726
	727	vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
	728
	729	ENDIF
[6140]	730	!
[5656]	731	DO jk = 1, jpkm1 ! Horizontal slab
	732	! ! ===============
	733
	734	! ! --------
	735	! Horizontal divergence ! div
	736	! ! --------
	737	DO jj = j1+1,j2
	738	DO ji = i1,i2 ! vector opt.
[13216]	739	zbtr = rn_sponge_dyn * r1_Dt * fspt(ji,jj) / e3t(ji,jj,jk,Kbb_a)
[12377]	740	hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v(ji,jj ,jk,Kbb_a) * vbdiff(ji,jj ,jk) &
	741	& -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kbb_a) * vbdiff(ji,jj-1,jk) ) * zbtr
[5656]	742	END DO
	743	END DO
	744	DO jj = j1,j2
	745	DO ji = i1,i2-1 ! vector opt.
[13216]	746	zbtr = rn_sponge_dyn * r1_Dt * fspf(ji,jj) * e3f(ji,jj,jk)
[5656]	747	rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) &
[6140]	748	& -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr
[5656]	749	END DO
	750	END DO
[6140]	751	END DO
[5656]	752
	753	! ! ===============
	754	!
	755
[9019]	756	imax = i2 - 1
[13286]	757	ind1 = jpiglo - ( nn_hls + nbghostcells + 2 ) ! East
	758	DO ji = mi0(ind1), mi1(ind1)
	759	imax = MIN(imax,ji)
	760	END DO
	761
[5656]	762	DO jj = j1+1, j2
	763	DO ji = i1+1, imax ! vector opt.
[6140]	764	IF( .NOT. tabspongedone_u(ji,jj) ) THEN
	765	DO jk = 1, jpkm1
[13216]	766	uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) &
[12377]	767	& - ( rotdiff (ji ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) &
[6140]	768	& + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj ,jk)) * r1_e1u(ji,jj)
[5656]	769	END DO
	770	ENDIF
	771	END DO
	772	END DO
[6140]	773	!
[5656]	774	tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
[6140]	775	!
[5656]	776	DO jj = j1+1, j2-1
	777	DO ji = i1+1, i2-1 ! vector opt.
[6140]	778	IF( .NOT. tabspongedone_v(ji,jj) ) THEN
	779	DO jk = 1, jpkm1
[13216]	780	vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) &
[12377]	781	& + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) ) &
[13216]	782	& + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj) &
	783	& - rn_trelax_dyn * r1_Dt * fspv(ji,jj) * vbdiff(ji,jj,jk)
[5656]	784	END DO
	785	ENDIF
	786	END DO
	787	END DO
	788	tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
	789	ENDIF
[6140]	790	!
[5656]	791	END SUBROUTINE interpvn_sponge
	792
[390]	793	#else
[9019]	794	!!----------------------------------------------------------------------
	795	!! Empty module no AGRIF zoom
	796	!!----------------------------------------------------------------------
[636]	797	CONTAINS
[9570]	798	SUBROUTINE agrif_oce_sponge_empty
	799	WRITE(,) 'agrif_oce_sponge : You should not have seen this print! error?'
	800	END SUBROUTINE agrif_oce_sponge_empty
[390]	801	#endif
[636]	802
[6140]	803	!!======================================================================
[9570]	804	END MODULE agrif_oce_sponge

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

Download in other formats: