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/2019/dev_r12072_MERGE_OPTION2_2019/src/NST – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/NST/agrif_oce_sponge.F90 @ 12256

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

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

Download in other formats: