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

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

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

Download in other formats: