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_r12558_HPC-08_epico_Extra_Halo/src/NST – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12558_HPC-08_epico_Extra_Halo/src/NST/agrif_oce_sponge.F90 @ 13229

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

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

Download in other formats: