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_opa_sponge.F90 in branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 9031

Last change on this file since 9031 was 9031, checked in by timgraham, 6 years ago
Resolved AGRIF conflicts
Property svn:keywords set to `Id`
File size: 23.4 KB

Rev	Line
[3680]	1	#define SPONGE && define SPONGE_TOP
[390]	2
[5656]	3	MODULE agrif_opa_sponge
[6140]	4	!!======================================================================
[9019]	5	!! * MODULE agrif_opa_interp *
	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	!!----------------------------------------------------------------------
[9019]	32	!! NEMO/NST 4.0 , NEMO Consortium (2017)
[1156]	33	!! $Id$
[2528]	34	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[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
	43	!!----------------------------------------------------------------------
[6140]	44	!
[390]	45	#if defined SPONGE
[9031]	46	!! Assume persistence:
	47	timecoeff = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
	48
[5656]	49	CALL Agrif_Sponge
[9019]	50	Agrif_SpecialValue = 0._wp
[390]	51	Agrif_UseSpecialValue = .TRUE.
[9019]	52	tabspongedone_tsn = .FALSE.
	53	!
	54	CALL Agrif_Bc_Variable( tsn_sponge_id, calledweight=zcoef, procname=interptsn_sponge )
	55	!
[5656]	56	Agrif_UseSpecialValue = .FALSE.
[390]	57	#endif
[6140]	58	!
[636]	59	END SUBROUTINE Agrif_Sponge_Tra
[390]	60
[6140]	61
[782]	62	SUBROUTINE Agrif_Sponge_dyn
[9019]	63	!!----------------------------------------------------------------------
	64	!! * ROUTINE Agrif_Sponge_dyn *
	65	!!----------------------------------------------------------------------
	66	REAL(wp) :: zcoef ! local scalar
	67	!!----------------------------------------------------------------------
	68	!
[390]	69	#if defined SPONGE
[9031]	70	timecoeff = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
	71
	72	Agrif_SpecialValue=0.
[782]	73	Agrif_UseSpecialValue = ln_spc_dyn
[9019]	74	!
[5656]	75	tabspongedone_u = .FALSE.
	76	tabspongedone_v = .FALSE.
[9019]	77	CALL Agrif_Bc_Variable( un_sponge_id, calledweight=zcoef, procname=interpun_sponge )
	78	!
[5656]	79	tabspongedone_u = .FALSE.
	80	tabspongedone_v = .FALSE.
[9019]	81	CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge )
	82	!
[390]	83	Agrif_UseSpecialValue = .FALSE.
	84	#endif
[6140]	85	!
[636]	86	END SUBROUTINE Agrif_Sponge_dyn
[390]	87
[6140]	88
[3680]	89	SUBROUTINE Agrif_Sponge
[9019]	90	!!----------------------------------------------------------------------
	91	!! * ROUTINE Agrif_Sponge *
	92	!!----------------------------------------------------------------------
	93	INTEGER :: ji, jj, ind1, ind2
	94	INTEGER :: ispongearea
	95	REAL(wp) :: z1_spongearea
	96	REAL(wp), DIMENSION(jpi,jpj) :: ztabramp
	97	!!----------------------------------------------------------------------
	98	!
[3680]	99	#if defined SPONGE \|\| defined SPONGE_TOP
[4153]	100	IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
[9019]	101	! Define ramp from boundaries towards domain interior at T-points
[4153]	102	! Store it in ztabramp
[3680]	103
[5656]	104	ispongearea = 2 + nn_sponge_len * Agrif_irhox()
[9019]	105	z1_spongearea = 1._wp / REAL( ispongearea - 1 )
	106
[5656]	107	ztabramp(:,:) = 0._wp
[4153]	108
[9019]	109	! --- West --- !
[4153]	110	IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
[9019]	111	ind1 = 1+nbghostcells
	112	ind2 = 1+nbghostcells + (ispongearea-1)
[4153]	113	DO jj = 1, jpj
[9019]	114	DO ji = ind1, ind2
	115	ztabramp(ji,jj) = REAL( ind2 - ji ) * z1_spongearea * umask(ind1,jj,1)
	116	END DO
[4153]	117	ENDDO
	118	ENDIF
	119
[9019]	120	! --- East --- !
[4153]	121	IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
[9019]	122	ind1 = nlci - (1+nbghostcells) - (ispongearea-1)
	123	ind2 = nlci - (1+nbghostcells)
[4153]	124	DO jj = 1, jpj
[9019]	125	DO ji = ind1, ind2
	126	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ji - ind2 ) * z1_spongearea * umask(ind2-1,jj,1) )
	127	ENDDO
[4153]	128	ENDDO
	129	ENDIF
	130
[9019]	131	! --- South --- !
[4153]	132	IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
[9019]	133	ind1 = 1+nbghostcells
	134	ind2 = 1+nbghostcells + (ispongearea-1)
	135	DO jj = ind1, ind2
	136	DO ji = 1, jpi
	137	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - jj ) * z1_spongearea * vmask(ji,ind1,1) )
	138	END DO
[4153]	139	ENDDO
	140	ENDIF
	141
[9019]	142	! --- North --- !
[4153]	143	IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
[9019]	144	ind1 = nlcj - (1+nbghostcells) - (ispongearea-1)
	145	ind2 = nlcj - (1+nbghostcells)
	146	DO jj = ind1, ind2
	147	DO ji = 1, jpi
	148	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( jj - ind2 ) * z1_spongearea * vmask(ji,ind2-1,1) )
	149	END DO
[4153]	150	ENDDO
	151	ENDIF
	152
	153	ENDIF
	154
[3680]	155	! Tracers
	156	IF( .NOT. spongedoneT ) THEN
[5656]	157	fsaht_spu(:,:) = 0._wp
	158	fsaht_spv(:,:) = 0._wp
	159	DO jj = 2, jpjm1
	160	DO ji = 2, jpim1 ! vector opt.
[9019]	161	fsaht_spu(ji,jj) = 0.5_wp * visc_tra * ( ztabramp(ji,jj) + ztabramp(ji+1,jj ) )
	162	fsaht_spv(ji,jj) = 0.5_wp * visc_tra * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1) )
[5656]	163	END DO
	164	END DO
	165	CALL lbc_lnk( fsaht_spu, 'U', 1. ) ! Lateral boundary conditions
	166	CALL lbc_lnk( fsaht_spv, 'V', 1. )
[9019]	167
[3680]	168	spongedoneT = .TRUE.
	169	ENDIF
	170
	171	! Dynamics
	172	IF( .NOT. spongedoneU ) THEN
[5656]	173	fsahm_spt(:,:) = 0._wp
	174	fsahm_spf(:,:) = 0._wp
	175	DO jj = 2, jpjm1
	176	DO ji = 2, jpim1 ! vector opt.
	177	fsahm_spt(ji,jj) = visc_dyn * ztabramp(ji,jj)
	178	fsahm_spf(ji,jj) = 0.25_wp * visc_dyn * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1) &
	179	& +ztabramp(ji,jj) + ztabramp(ji+1,jj ) )
	180	END DO
	181	END DO
	182	CALL lbc_lnk( fsahm_spt, 'T', 1. ) ! Lateral boundary conditions
	183	CALL lbc_lnk( fsahm_spf, 'F', 1. )
[9019]	184
[3680]	185	spongedoneU = .TRUE.
	186	ENDIF
	187	!
	188	#endif
[6140]	189	!
[3680]	190	END SUBROUTINE Agrif_Sponge
	191
[9019]	192	SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
	193	!!----------------------------------------------------------------------
	194	!! * ROUTINE interptsn_sponge *
	195	!!----------------------------------------------------------------------
	196	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
	197	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
	198	LOGICAL , INTENT(in ) :: before
[6140]	199	!
[5656]	200	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	201	INTEGER :: iku, ikv
	202	REAL(wp) :: ztsa, zabe1, zabe2, zbtr
[9031]	203	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv
	204	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
	205	! vertical interpolation:
	206	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
	207	REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
	208	REAL(wp), DIMENSION(k1:k2) :: h_in
	209	REAL(wp), DIMENSION(1:jpk) :: h_out
	210	INTEGER :: N_in, N_out
	211	REAL(wp) :: h_diff
[9019]	212	!!----------------------------------------------------------------------
[5656]	213	!
[6140]	214	IF( before ) THEN
[9031]	215	DO jn = 1, jpts
	216	DO jk=k1,k2
	217	DO jj=j1,j2
	218	DO ji=i1,i2
	219	tabres(ji,jj,jk,jn) = tsb(ji,jj,jk,jn)
	220	END DO
	221	END DO
	222	END DO
	223	END DO
	224
	225	# if defined key_vertical
	226	DO jk=k1,k2
	227	DO jj=j1,j2
	228	DO ji=i1,i2
	229	tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t_n(ji,jj,jk)
	230	END DO
	231	END DO
	232	END DO
	233	# endif
	234
[5656]	235	ELSE
[6140]	236	!
[9031]	237	# if defined key_vertical
	238	tabres_child(:,:,:,:) = 0.
	239	DO jj=j1,j2
	240	DO ji=i1,i2
	241	N_in = 0
	242	DO jk=k1,k2 !k2 = jpk of parent grid
	243	IF (tabres(ji,jj,jk,n2) == 0) EXIT
	244	N_in = N_in + 1
	245	tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)
	246	h_in(N_in) = tabres(ji,jj,jk,n2)
	247	END DO
	248	N_out = 0
	249	DO jk=1,jpk ! jpk of child grid
	250	IF (tmask(ji,jj,jk) == 0) EXIT
	251	N_out = N_out + 1
	252	h_out(jk) = e3t_n(ji,jj,jk) !Child grid scale factors. Could multiply by e1e2t here instead of division above
	253	ENDDO
	254	IF (N_in > 0) THEN
	255	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
	256	tabres(ji,jj,k2,:) = tabres(ji,jj,k2-1,:) !what is this line for?????
	257	DO jn=1,jpts
	258	call reconstructandremap(tabin(1:N_in,jn),h_in,tabres_child(ji,jj,1:N_out,jn),h_out,N_in,N_out)
	259	ENDDO
	260	ENDIF
	261	ENDDO
	262	ENDDO
	263	# endif
	264
	265	DO jj=j1,j2
	266	DO ji=i1,i2
	267	DO jk=1,jpkm1
	268	# if defined key_vertical
	269	tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres_child(ji,jj,jk,1:jpts)
	270	# else
	271	tsbdiff(ji,jj,jk,1:jpts) = tsb(ji,jj,jk,1:jpts) - tabres(ji,jj,jk,1:jpts)
	272	# endif
	273	ENDDO
	274	ENDDO
	275	ENDDO
	276
[5656]	277	DO jn = 1, jpts
	278	DO jk = 1, jpkm1
	279	DO jj = j1,j2-1
	280	DO ji = i1,i2-1
[6140]	281	zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
	282	zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
[5656]	283	ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) )
	284	ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
[6140]	285	END DO
	286	END DO
	287	!
[5656]	288	IF( ln_zps ) THEN ! set gradient at partial step level
	289	DO jj = j1,j2-1
	290	DO ji = i1,i2-1
	291	! last level
	292	iku = mbku(ji,jj)
	293	ikv = mbkv(ji,jj)
[6140]	294	IF( iku == jk ) ztu(ji,jj,jk) = 0._wp
	295	IF( ikv == jk ) ztv(ji,jj,jk) = 0._wp
[5656]	296	END DO
	297	END DO
	298	ENDIF
[6140]	299	END DO
	300	!
[5656]	301	DO jk = 1, jpkm1
	302	DO jj = j1+1,j2-1
	303	DO ji = i1+1,i2-1
	304	IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
[6140]	305	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
[5656]	306	! horizontal diffusive trends
[9019]	307	ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )
[5656]	308	! add it to the general tracer trends
	309	tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa
	310	ENDIF
[6140]	311	END DO
	312	END DO
	313	END DO
	314	!
	315	END DO
	316	!
[5656]	317	tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
[6140]	318	!
[5656]	319	ENDIF
[6140]	320	!
[5656]	321	END SUBROUTINE interptsn_sponge
	322
[9031]	323	SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
	324	!!---------------------------------------------
	325	!! * ROUTINE interpun_sponge *
	326	!!---------------------------------------------
	327	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
	328	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
	329	LOGICAL, INTENT(in) :: before
[6140]	330
[9031]	331	INTEGER :: ji,jj,jk,jmax
	332
[5656]	333	! sponge parameters
[9031]	334	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff
	335	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
	336	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
	337	! vertical interpolation:
	338	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
	339	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
	340	REAL(wp), DIMENSION(1:jpk) :: h_out
	341	INTEGER ::N_in,N_out
	342	!!---------------------------------------------
[5656]	343	!
[6140]	344	IF( before ) THEN
[9031]	345	DO jk=1,jpkm1
	346	DO jj=j1,j2
	347	DO ji=i1,i2
	348	tabres(ji,jj,jk,m1) = ub(ji,jj,jk)
	349	# if defined key_vertical
	350	tabres(ji,jj,jk,m2) = e3u_n(ji,jj,jk)*umask(ji,jj,jk)
	351	# endif
	352	END DO
	353	END DO
	354	END DO
	355
[5656]	356	ELSE
[9031]	357
	358	# if defined key_vertical
	359	tabres_child(:,:,:) = 0._wp
	360	DO jj=j1,j2
	361	DO ji=i1,i2
	362	N_in = 0
	363	DO jk=k1,k2
	364	IF (tabres(ji,jj,jk,m2) == 0) EXIT
	365	N_in = N_in + 1
	366	tabin(jk) = tabres(ji,jj,jk,m1)
	367	h_in(N_in) = tabres(ji,jj,jk,m2)
	368	ENDDO
	369	!
	370	IF (N_in == 0) THEN
	371	tabres_child(ji,jj,:) = 0.
	372	CYCLE
	373	ENDIF
	374
	375	N_out = 0
	376	DO jk=1,jpk
	377	if (umask(ji,jj,jk) == 0) EXIT
	378	N_out = N_out + 1
	379	h_out(N_out) = e3u_n(ji,jj,jk)
	380	ENDDO
	381
	382	IF (N_out == 0) THEN
	383	tabres_child(ji,jj,:) = 0.
	384	CYCLE
	385	ENDIF
	386
	387	IF (N_in * N_out > 0) THEN
	388	h_diff = sum(h_out(1:N_out))-sum(h_in(1:N_in))
	389	if (h_diff < -1.e4) then
	390	print *,'CHECK YOUR BATHY ...', h_diff, sum(h_out(1:N_out)), sum(h_in(1:N_in))
	391	endif
	392	ENDIF
	393	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)
	394
	395	ENDDO
	396	ENDDO
	397
	398	ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
	399	#else
	400	ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
	401	#endif
	402	>>>>>>> .merge-right.r9019
[6140]	403	!
[5656]	404	DO jk = 1, jpkm1 ! Horizontal slab
	405	! ! ===============
	406
	407	! ! --------
	408	! Horizontal divergence ! div
	409	! ! --------
	410	DO jj = j1,j2
	411	DO ji = i1+1,i2 ! vector opt.
[6140]	412	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
	413	hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)e3u_n(ji ,jj,jk) ubdiff(ji ,jj,jk) &
	414	& -e2u(ji-1,jj)e3u_n(ji-1,jj,jk) ubdiff(ji-1,jj,jk) ) * zbtr
[5656]	415	END DO
	416	END DO
	417
	418	DO jj = j1,j2-1
	419	DO ji = i1,i2 ! vector opt.
[6140]	420	zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
[9019]	421	rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &
	422	& +e1u(ji,jj ) * ubdiff(ji,jj ,jk) ) * fmask(ji,jj,jk) * zbtr
[5656]	423	END DO
	424	END DO
[6140]	425	END DO
[5656]	426	!
	427	DO jj = j1+1, j2-1
	428	DO ji = i1+1, i2-1 ! vector opt.
	429
	430	IF (.NOT. tabspongedone_u(ji,jj)) THEN
	431	DO jk = 1, jpkm1 ! Horizontal slab
	432	ze2u = rotdiff (ji,jj,jk)
	433	ze1v = hdivdiff(ji,jj,jk)
	434	! horizontal diffusive trends
[9019]	435	zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) ) &
	436	+ ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj)
[5656]	437
	438	! add it to the general momentum trends
	439	ua(ji,jj,jk) = ua(ji,jj,jk) + zua
	440
	441	END DO
	442	ENDIF
	443
	444	END DO
	445	END DO
	446
	447	tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
	448
	449	jmax = j2-1
[9019]	450	IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-nbghostcells-2) ! North
[5656]	451
	452	DO jj = j1+1, jmax
	453	DO ji = i1+1, i2 ! vector opt.
	454
	455	IF (.NOT. tabspongedone_v(ji,jj)) THEN
	456	DO jk = 1, jpkm1 ! Horizontal slab
	457	ze2u = rotdiff (ji,jj,jk)
	458	ze1v = hdivdiff(ji,jj,jk)
	459
	460	! horizontal diffusive trends
[9019]	461	zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) ) &
	462	+ ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
[5656]	463
	464	! add it to the general momentum trends
	465	va(ji,jj,jk) = va(ji,jj,jk) + zva
	466	END DO
	467	ENDIF
[6140]	468	!
[5656]	469	END DO
	470	END DO
[6140]	471	!
[5656]	472	tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
[6140]	473	!
[5656]	474	ENDIF
[6140]	475	!
[5656]	476	END SUBROUTINE interpun_sponge
	477
[9031]	478	SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir)
	479	!!---------------------------------------------
	480	!! * ROUTINE interpvn_sponge *
	481	!!---------------------------------------------
	482	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
	483	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
	484	LOGICAL, INTENT(in) :: before
	485	INTEGER, INTENT(in) :: nb , ndir
[6140]	486	!
[9019]	487	INTEGER :: ji, jj, jk
	488	INTEGER :: imax
	489	REAL(wp):: ze2u, ze1v, zua, zva, zbtr
	490	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: vbdiff, rotdiff, hdivdiff
	491	!!----------------------------------------------------------------------
[9031]	492	INTEGER :: ji, jj, jk, imax
	493	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, h_diff
	494	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
	495	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
	496	! vertical interpolation:
	497	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
	498	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
	499	REAL(wp), DIMENSION(1:jpk) :: h_out
	500	INTEGER :: N_in, N_out
	501	!!---------------------------------------------
	502	>>>>>>> .merge-right.r9019
	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
	619	SUBROUTINE agrif_opa_sponge_empty
	620	WRITE(,) 'agrif_opa_sponge : You should not have seen this print! error?'
	621	END SUBROUTINE agrif_opa_sponge_empty
[390]	622	#endif
[636]	623
[6140]	624	!!======================================================================
[636]	625	END MODULE agrif_opa_sponge

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

Download in other formats: