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_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11233_AGRIF-05_jchanut_vert_coord_interp/src/NST/agrif_oce_sponge.F90 @ 11802

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

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

Download in other formats: