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/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/NST – NEMO

Context Navigation

source: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/NST/agrif_oce_sponge.F90 @ 10419

Last change on this file since 10419 was 10170, checked in by smasson, 6 years ago
dev_r10164_HPC09_ESIWACE_PREP_MERGE: action 2a: add report calls of lbc_lnk, see #2133
Property svn:keywords set to `Id`
File size: 23.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)
[390]	24
[636]	25	IMPLICIT NONE
	26	PRIVATE
[390]	27
[5656]	28	PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn
	29	PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge
[390]	30
[1156]	31	!!----------------------------------------------------------------------
[9598]	32	!! NEMO/NST 4.0 , NEMO Consortium (2018)
[1156]	33	!! $Id$
[10068]	34	!! Software governed by the CeCILL license (see ./LICENSE)
[1156]	35	!!----------------------------------------------------------------------
[5656]	36	CONTAINS
[636]	37
[782]	38	SUBROUTINE Agrif_Sponge_Tra
[9019]	39	!!----------------------------------------------------------------------
	40	!! * ROUTINE Agrif_Sponge_Tra *
	41	!!----------------------------------------------------------------------
	42	REAL(wp) :: zcoef ! local scalar
[9056]	43
[9019]	44	!!----------------------------------------------------------------------
[6140]	45	!
[390]	46	#if defined SPONGE
[9031]	47	!! Assume persistence:
[9056]	48	zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
[9031]	49
[5656]	50	CALL Agrif_Sponge
[9019]	51	Agrif_SpecialValue = 0._wp
[390]	52	Agrif_UseSpecialValue = .TRUE.
[9019]	53	tabspongedone_tsn = .FALSE.
	54	!
	55	CALL Agrif_Bc_Variable( tsn_sponge_id, calledweight=zcoef, procname=interptsn_sponge )
	56	!
[5656]	57	Agrif_UseSpecialValue = .FALSE.
[390]	58	#endif
[6140]	59	!
[636]	60	END SUBROUTINE Agrif_Sponge_Tra
[390]	61
[6140]	62
[782]	63	SUBROUTINE Agrif_Sponge_dyn
[9019]	64	!!----------------------------------------------------------------------
	65	!! * ROUTINE Agrif_Sponge_dyn *
	66	!!----------------------------------------------------------------------
	67	REAL(wp) :: zcoef ! local scalar
	68	!!----------------------------------------------------------------------
	69	!
[390]	70	#if defined SPONGE
[9056]	71	zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
[9031]	72
	73	Agrif_SpecialValue=0.
[782]	74	Agrif_UseSpecialValue = ln_spc_dyn
[9019]	75	!
[5656]	76	tabspongedone_u = .FALSE.
	77	tabspongedone_v = .FALSE.
[9019]	78	CALL Agrif_Bc_Variable( un_sponge_id, calledweight=zcoef, procname=interpun_sponge )
	79	!
[5656]	80	tabspongedone_u = .FALSE.
	81	tabspongedone_v = .FALSE.
[9019]	82	CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge )
	83	!
[390]	84	Agrif_UseSpecialValue = .FALSE.
	85	#endif
[6140]	86	!
[636]	87	END SUBROUTINE Agrif_Sponge_dyn
[390]	88
[6140]	89
[3680]	90	SUBROUTINE Agrif_Sponge
[9019]	91	!!----------------------------------------------------------------------
	92	!! * ROUTINE Agrif_Sponge *
	93	!!----------------------------------------------------------------------
	94	INTEGER :: ji, jj, ind1, ind2
	95	INTEGER :: ispongearea
	96	REAL(wp) :: z1_spongearea
	97	REAL(wp), DIMENSION(jpi,jpj) :: ztabramp
	98	!!----------------------------------------------------------------------
	99	!
[3680]	100	#if defined SPONGE \|\| defined SPONGE_TOP
[4153]	101	IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
[9019]	102	! Define ramp from boundaries towards domain interior at T-points
[4153]	103	! Store it in ztabramp
[3680]	104
[9806]	105	ispongearea = 1 + nn_sponge_len * Agrif_irhox()
	106	z1_spongearea = 1._wp / REAL( ispongearea )
[9019]	107
[5656]	108	ztabramp(:,:) = 0._wp
[4153]	109
[9019]	110	! --- West --- !
[4153]	111	IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
[9019]	112	ind1 = 1+nbghostcells
[9806]	113	ind2 = 1+nbghostcells + ispongearea
[4153]	114	DO jj = 1, jpj
[9806]	115	DO ji = ind1, ind2
[9019]	116	ztabramp(ji,jj) = REAL( ind2 - ji ) * z1_spongearea * umask(ind1,jj,1)
	117	END DO
[4153]	118	ENDDO
	119	ENDIF
	120
[9019]	121	! --- East --- !
[4153]	122	IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
[9806]	123	ind1 = nlci - nbghostcells - ispongearea
	124	ind2 = nlci - nbghostcells
[4153]	125	DO jj = 1, jpj
[9019]	126	DO ji = ind1, ind2
[9786]	127	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ji - ind1 ) * z1_spongearea * umask(ind2-1,jj,1) )
[9019]	128	ENDDO
[4153]	129	ENDDO
	130	ENDIF
	131
[9019]	132	! --- South --- !
[4153]	133	IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
[9019]	134	ind1 = 1+nbghostcells
[9806]	135	ind2 = 1+nbghostcells + ispongearea
	136	DO jj = ind1, ind2
[9019]	137	DO ji = 1, jpi
	138	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - jj ) * z1_spongearea * vmask(ji,ind1,1) )
	139	END DO
[4153]	140	ENDDO
	141	ENDIF
	142
[9019]	143	! --- North --- !
[4153]	144	IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
[9806]	145	ind1 = nlcj - nbghostcells - ispongearea
	146	ind2 = nlcj - nbghostcells
[9019]	147	DO jj = ind1, ind2
	148	DO ji = 1, jpi
[9786]	149	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( jj - ind1 ) * z1_spongearea * vmask(ji,ind2-1,1) )
[9019]	150	END DO
[4153]	151	ENDDO
	152	ENDIF
	153
	154	ENDIF
	155
[3680]	156	! Tracers
	157	IF( .NOT. spongedoneT ) THEN
[5656]	158	fsaht_spu(:,:) = 0._wp
	159	fsaht_spv(:,:) = 0._wp
	160	DO jj = 2, jpjm1
	161	DO ji = 2, jpim1 ! vector opt.
[9019]	162	fsaht_spu(ji,jj) = 0.5_wp * visc_tra * ( ztabramp(ji,jj) + ztabramp(ji+1,jj ) )
	163	fsaht_spv(ji,jj) = 0.5_wp * visc_tra * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1) )
[5656]	164	END DO
	165	END DO
[10170]	166	CALL lbc_lnk( 'agrif_oce_sponge', fsaht_spu, 'U', 1. ) ! Lateral boundary conditions
	167	CALL lbc_lnk( 'agrif_oce_sponge', fsaht_spv, 'V', 1. )
[9019]	168
[3680]	169	spongedoneT = .TRUE.
	170	ENDIF
	171
	172	! Dynamics
	173	IF( .NOT. spongedoneU ) THEN
[5656]	174	fsahm_spt(:,:) = 0._wp
	175	fsahm_spf(:,:) = 0._wp
	176	DO jj = 2, jpjm1
	177	DO ji = 2, jpim1 ! vector opt.
	178	fsahm_spt(ji,jj) = visc_dyn * ztabramp(ji,jj)
[9806]	179	fsahm_spf(ji,jj) = 0.25_wp * visc_dyn * ( ztabramp(ji ,jj ) + ztabramp(ji ,jj+1) &
	180	& +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj ) )
[5656]	181	END DO
	182	END DO
[10170]	183	CALL lbc_lnk( 'agrif_oce_sponge', fsahm_spt, 'T', 1. ) ! Lateral boundary conditions
	184	CALL lbc_lnk( 'agrif_oce_sponge', fsahm_spf, 'F', 1. )
[9019]	185
[3680]	186	spongedoneU = .TRUE.
	187	ENDIF
	188	!
	189	#endif
[6140]	190	!
[3680]	191	END SUBROUTINE Agrif_Sponge
	192
[9019]	193	SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
	194	!!----------------------------------------------------------------------
	195	!! * ROUTINE interptsn_sponge *
	196	!!----------------------------------------------------------------------
	197	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
	198	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
	199	LOGICAL , INTENT(in ) :: before
[6140]	200	!
[5656]	201	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	202	INTEGER :: iku, ikv
	203	REAL(wp) :: ztsa, zabe1, zabe2, zbtr
[9031]	204	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv
	205	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
	206	! vertical interpolation:
	207	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
	208	REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
	209	REAL(wp), DIMENSION(k1:k2) :: h_in
	210	REAL(wp), DIMENSION(1:jpk) :: h_out
	211	INTEGER :: N_in, N_out
	212	REAL(wp) :: h_diff
[9019]	213	!!----------------------------------------------------------------------
[5656]	214	!
[6140]	215	IF( before ) THEN
[9031]	216	DO jn = 1, jpts
	217	DO jk=k1,k2
	218	DO jj=j1,j2
	219	DO ji=i1,i2
	220	tabres(ji,jj,jk,jn) = tsb(ji,jj,jk,jn)
	221	END DO
	222	END DO
	223	END DO
	224	END DO
	225
	226	# if defined key_vertical
	227	DO jk=k1,k2
	228	DO jj=j1,j2
	229	DO ji=i1,i2
	230	tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)
	231	END DO
	232	END DO
	233	END DO
	234	# endif
	235
[5656]	236	ELSE
[6140]	237	!
[9031]	238	# if defined key_vertical
	239	tabres_child(:,:,:,:) = 0.
	240	DO jj=j1,j2
	241	DO ji=i1,i2
	242	N_in = 0
	243	DO jk=k1,k2 !k2 = jpk of parent grid
	244	IF (tabres(ji,jj,jk,n2) == 0) EXIT
	245	N_in = N_in + 1
	246	tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)
	247	h_in(N_in) = tabres(ji,jj,jk,n2)
	248	END DO
	249	N_out = 0
	250	DO jk=1,jpk ! jpk of child grid
	251	IF (tmask(ji,jj,jk) == 0) EXIT
	252	N_out = N_out + 1
	253	h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above
	254	ENDDO
	255	IF (N_in > 0) THEN
	256	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
	257	tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for?????
	258	DO jn=1,jpts
	259	call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
	260	ENDDO
	261	ENDIF
	262	ENDDO
	263	ENDDO
	264	# endif
	265
	266	DO jj=j1,j2
	267	DO ji=i1,i2
	268	DO jk=1,jpkm1
	269	# if defined key_vertical
	270	tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)
	271	# else
	272	tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts)
	273	# endif
	274	ENDDO
	275	ENDDO
	276	ENDDO
	277
[5656]	278	DO jn = 1, jpts
	279	DO jk = 1, jpkm1
[9806]	280	ztu(i1:i2,j1:j2,jk) = 0._wp
	281	DO jj = j1,j2
[5656]	282	DO ji = i1,i2-1
[6140]	283	zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
[9806]	284	ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) )
	285	END DO
	286	END DO
	287	ztv(i1:i2,j1:j2,jk) = 0._wp
	288	DO ji = i1,i2
	289	DO jj = j1,j2-1
[6140]	290	zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
[5656]	291	ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
[6140]	292	END DO
	293	END DO
	294	!
[5656]	295	IF( ln_zps ) THEN ! set gradient at partial step level
[9806]	296	DO jj = j1,j2
	297	DO ji = i1,i2
[5656]	298	! last level
	299	iku = mbku(ji,jj)
	300	ikv = mbkv(ji,jj)
[6140]	301	IF( iku == jk ) ztu(ji,jj,jk) = 0._wp
	302	IF( ikv == jk ) ztv(ji,jj,jk) = 0._wp
[5656]	303	END DO
	304	END DO
	305	ENDIF
[6140]	306	END DO
	307	!
[5656]	308	DO jk = 1, jpkm1
	309	DO jj = j1+1,j2-1
	310	DO ji = i1+1,i2-1
	311	IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
[6140]	312	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
[5656]	313	! horizontal diffusive trends
[9019]	314	ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )
[5656]	315	! add it to the general tracer trends
	316	tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa
	317	ENDIF
[6140]	318	END DO
	319	END DO
	320	END DO
	321	!
	322	END DO
	323	!
[5656]	324	tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
[6140]	325	!
[5656]	326	ENDIF
[6140]	327	!
[5656]	328	END SUBROUTINE interptsn_sponge
	329
[9031]	330	SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
	331	!!---------------------------------------------
	332	!! * ROUTINE interpun_sponge *
	333	!!---------------------------------------------
	334	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
	335	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
	336	LOGICAL, INTENT(in) :: before
[6140]	337
[9031]	338	INTEGER :: ji,jj,jk,jmax
	339
[5656]	340	! sponge parameters
[9031]	341	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff
	342	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
	343	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
	344	! vertical interpolation:
	345	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
	346	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
	347	REAL(wp), DIMENSION(1:jpk) :: h_out
	348	INTEGER ::N_in,N_out
	349	!!---------------------------------------------
[5656]	350	!
[6140]	351	IF( before ) THEN
[9031]	352	DO jk=1,jpkm1
	353	DO jj=j1,j2
	354	DO ji=i1,i2
	355	tabres(ji,jj,jk,m1) = ub(ji,jj,jk)
	356	# if defined key_vertical
	357	tabres(ji,jj,jk,m2) = e3u_n(ji,jj,jk)*umask(ji,jj,jk)
	358	# endif
	359	END DO
	360	END DO
	361	END DO
	362
[5656]	363	ELSE
[9031]	364
	365	# if defined key_vertical
	366	tabres_child(:,:,:) = 0._wp
	367	DO jj=j1,j2
	368	DO ji=i1,i2
	369	N_in = 0
	370	DO jk=k1,k2
	371	IF (tabres(ji,jj,jk,m2) == 0) EXIT
	372	N_in = N_in + 1
	373	tabin(jk) = tabres(ji,jj,jk,m1)
	374	h_in(N_in) = tabres(ji,jj,jk,m2)
	375	ENDDO
	376	!
	377	IF (N_in == 0) THEN
	378	tabres_child(ji,jj,:) = 0.
	379	CYCLE
	380	ENDIF
	381
	382	N_out = 0
	383	DO jk=1,jpk
	384	if (umask(ji,jj,jk) == 0) EXIT
	385	N_out = N_out + 1
	386	h_out(N_out) = e3u_n(ji,jj,jk)
	387	ENDDO
	388
	389	IF (N_out == 0) THEN
	390	tabres_child(ji,jj,:) = 0.
	391	CYCLE
	392	ENDIF
	393
	394	IF (N_in * N_out > 0) THEN
	395	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
	396	if (h_diff < -1.e4) then
	397	print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
	398	endif
	399	ENDIF
	400	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)
	401
	402	ENDDO
	403	ENDDO
	404
	405	ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
	406	#else
	407	ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
	408	#endif
[6140]	409	!
[5656]	410	DO jk = 1, jpkm1 ! Horizontal slab
	411	! ! ===============
	412
	413	! ! --------
	414	! Horizontal divergence ! div
	415	! ! --------
	416	DO jj = j1,j2
	417	DO ji = i1+1,i2 ! vector opt.
[6140]	418	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
	419	hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)e3u_n(ji ,jj,jk) ubdiff(ji ,jj,jk) &
	420	& -e2u(ji-1,jj)e3u_n(ji-1,jj,jk) ubdiff(ji-1,jj,jk) ) * zbtr
[5656]	421	END DO
	422	END DO
	423
	424	DO jj = j1,j2-1
	425	DO ji = i1,i2 ! vector opt.
[6140]	426	zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
[9019]	427	rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &
	428	& +e1u(ji,jj ) * ubdiff(ji,jj ,jk) ) * fmask(ji,jj,jk) * zbtr
[5656]	429	END DO
	430	END DO
[6140]	431	END DO
[5656]	432	!
	433	DO jj = j1+1, j2-1
	434	DO ji = i1+1, i2-1 ! vector opt.
	435
	436	IF (.NOT. tabspongedone_u(ji,jj)) THEN
	437	DO jk = 1, jpkm1 ! Horizontal slab
	438	ze2u = rotdiff (ji,jj,jk)
	439	ze1v = hdivdiff(ji,jj,jk)
	440	! horizontal diffusive trends
[9019]	441	zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) ) &
	442	+ ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj)
[5656]	443
	444	! add it to the general momentum trends
	445	ua(ji,jj,jk) = ua(ji,jj,jk) + zua
	446
	447	END DO
	448	ENDIF
	449
	450	END DO
	451	END DO
	452
	453	tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
	454
	455	jmax = j2-1
[9019]	456	IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-nbghostcells-2) ! North
[5656]	457
	458	DO jj = j1+1, jmax
	459	DO ji = i1+1, i2 ! vector opt.
	460
	461	IF (.NOT. tabspongedone_v(ji,jj)) THEN
	462	DO jk = 1, jpkm1 ! Horizontal slab
	463	ze2u = rotdiff (ji,jj,jk)
	464	ze1v = hdivdiff(ji,jj,jk)
	465
	466	! horizontal diffusive trends
[9019]	467	zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) ) &
	468	+ ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
[5656]	469
	470	! add it to the general momentum trends
	471	va(ji,jj,jk) = va(ji,jj,jk) + zva
	472	END DO
	473	ENDIF
[6140]	474	!
[5656]	475	END DO
	476	END DO
[6140]	477	!
[5656]	478	tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
[6140]	479	!
[5656]	480	ENDIF
[6140]	481	!
[5656]	482	END SUBROUTINE interpun_sponge
	483
[9031]	484	SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir)
	485	!!---------------------------------------------
	486	!! * ROUTINE interpvn_sponge *
	487	!!---------------------------------------------
	488	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
	489	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
	490	LOGICAL, INTENT(in) :: before
	491	INTEGER, INTENT(in) :: nb , ndir
[6140]	492	!
[9031]	493	INTEGER :: ji, jj, jk, imax
	494	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff
	495	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
	496	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
	497	! vertical interpolation:
	498	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
	499	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
	500	REAL(wp), DIMENSION(1:jpk) :: h_out
	501	INTEGER :: N_in, N_out
	502	!!---------------------------------------------
	503
[6140]	504	IF( before ) THEN
[9031]	505	DO jk=1,jpkm1
	506	DO jj=j1,j2
	507	DO ji=i1,i2
	508	tabres(ji,jj,jk,m1) = vb(ji,jj,jk)
	509	# if defined key_vertical
	510	tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v_n(ji,jj,jk)
	511	# endif
	512	END DO
	513	END DO
	514	END DO
[5656]	515	ELSE
[9031]	516
	517	# if defined key_vertical
	518	tabres_child(:,:,:) = 0._wp
	519	DO jj=j1,j2
	520	DO ji=i1,i2
	521	N_in = 0
	522	DO jk=k1,k2
	523	IF (tabres(ji,jj,jk,m2) == 0) EXIT
	524	N_in = N_in + 1
	525	tabin(jk) = tabres(ji,jj,jk,m1)
	526	h_in(N_in) = tabres(ji,jj,jk,m2)
	527	ENDDO
	528
	529	IF (N_in == 0) THEN
	530	tabres_child(ji,jj,:) = 0.
	531	CYCLE
	532	ENDIF
	533
	534	N_out = 0
	535	DO jk=1,jpk
	536	if (vmask(ji,jj,jk) == 0) EXIT
	537	N_out = N_out + 1
	538	h_out(N_out) = e3v_n(ji,jj,jk)
	539	ENDDO
	540
	541	IF (N_in * N_out > 0) THEN
	542	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
	543	if (h_diff < -1.e4) then
	544	print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
	545	endif
	546	ENDIF
	547	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)
	548	ENDDO
	549	ENDDO
	550
	551	vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)
	552	# else
	553	vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
	554	# endif
[6140]	555	!
[5656]	556	DO jk = 1, jpkm1 ! Horizontal slab
	557	! ! ===============
	558
	559	! ! --------
	560	! Horizontal divergence ! div
	561	! ! --------
	562	DO jj = j1+1,j2
	563	DO ji = i1,i2 ! vector opt.
[6140]	564	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
	565	hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v_n(ji,jj ,jk) * vbdiff(ji,jj ,jk) &
	566	& -e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk) ) * zbtr
[5656]	567	END DO
	568	END DO
	569	DO jj = j1,j2
	570	DO ji = i1,i2-1 ! vector opt.
[6140]	571	zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
[5656]	572	rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) &
[6140]	573	& -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr
[5656]	574	END DO
	575	END DO
[6140]	576	END DO
[5656]	577
	578	! ! ===============
	579	!
	580
[9019]	581	imax = i2 - 1
	582	IF ((nbondi == 1).OR.(nbondi == 2)) imax = MIN(imax,nlci-nbghostcells-2) ! East
[5656]	583
	584	DO jj = j1+1, j2
	585	DO ji = i1+1, imax ! vector opt.
[6140]	586	IF( .NOT. tabspongedone_u(ji,jj) ) THEN
	587	DO jk = 1, jpkm1
	588	ua(ji,jj,jk) = ua(ji,jj,jk) &
	589	& - ( rotdiff (ji ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) ) &
	590	& + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj ,jk)) * r1_e1u(ji,jj)
[5656]	591	END DO
	592	ENDIF
	593	END DO
	594	END DO
[6140]	595	!
[5656]	596	tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
[6140]	597	!
[5656]	598	DO jj = j1+1, j2-1
	599	DO ji = i1+1, i2-1 ! vector opt.
[6140]	600	IF( .NOT. tabspongedone_v(ji,jj) ) THEN
	601	DO jk = 1, jpkm1
	602	va(ji,jj,jk) = va(ji,jj,jk) &
	603	& + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) ) &
	604	& + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj)
[5656]	605	END DO
	606	ENDIF
	607	END DO
	608	END DO
	609	tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
	610	ENDIF
[6140]	611	!
[5656]	612	END SUBROUTINE interpvn_sponge
	613
[390]	614	#else
[9019]	615	!!----------------------------------------------------------------------
	616	!! Empty module no AGRIF zoom
	617	!!----------------------------------------------------------------------
[636]	618	CONTAINS
[9570]	619	SUBROUTINE agrif_oce_sponge_empty
	620	WRITE(,) 'agrif_oce_sponge : You should not have seen this print! error?'
	621	END SUBROUTINE agrif_oce_sponge_empty
[390]	622	#endif
[636]	623
[6140]	624	!!======================================================================
[9570]	625	END MODULE agrif_oce_sponge

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

Download in other formats: