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_r13327_KERNEL-06_2_techene_e3/src/NST – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13327_KERNEL-06_2_techene_e3/src/NST/agrif_oce_sponge.F90 @ 14010

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

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

Download in other formats: