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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/NST/agrif_oce_sponge.F90 @ 10989

Last change on this file since 10989 was 10989, checked in by acc, 5 years ago
2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : Convert NST routines in preparation for getting AGRIF back up and running. AGRIF conv stage now works but requires some renaming of recently changes DIU modules (included in this commit). AGRIF compile and link stage not yet working (agrif routines need to be passed the time-level indices) but non-AGRIF SETTE tests are all OK
Property svn:keywords set to `Id`
File size: 23.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)
[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
[10425]	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
[10425]	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
[10989]	220	tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb)
[9031]	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
[10989]	230	tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kmm)
[9031]	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
[10989]	253	h_out(jk) = e3t(ji,jj,jk,Kmm) !Child grid scale factors. Could multiply by e1e2t here instead of division above
[9031]	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
[10989]	270	tsbdiff(ji,jj,jk,1:jpts) = ts(ji,jj,jk,1:jpts,Kbb) - tabres_child(ji,jj,jk,1:jpts)
[9031]	271	# else
[10989]	272	tsbdiff(ji,jj,jk,1:jpts) = ts(ji,jj,jk,1:jpts,Kbb) - tabres(ji,jj,jk,1:jpts)
[9031]	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
[10989]	283	zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm)
[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
[10989]	290	zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm)
[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
[10989]	312	zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
[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
[10989]	316	ts(ji,jj,jk,jn,Krhs) = ts(ji,jj,jk,jn,Krhs) + ztsa
[5656]	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
[10989]	355	tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb)
[9031]	356	# if defined key_vertical
[10989]	357	tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kmm)*umask(ji,jj,jk)
[9031]	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
[10989]	386	h_out(N_out) = e3u(ji,jj,jk,Kmm)
[9031]	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
[10989]	405	ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
[9031]	406	#else
[10989]	407	ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
[9031]	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.
[10989]	418	zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * fsahm_spt(ji,jj)
	419	hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)e3u(ji ,jj,jk,Kmm) ubdiff(ji ,jj,jk) &
	420	& -e2u(ji-1,jj)e3u(ji-1,jj,jk,Kmm) 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.
[10989]	426	zbtr = r1_e1e2f(ji,jj) * e3f(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
[10989]	441	zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm) ) &
[9019]	442	+ ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj)
[5656]	443
	444	! add it to the general momentum trends
[10989]	445	uu(ji,jj,jk,Krhs) = uu(ji,jj,jk,Krhs) + zua
[5656]	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
[10989]	467	zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm) ) &
[9019]	468	+ ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
[5656]	469
	470	! add it to the general momentum trends
[10989]	471	vv(ji,jj,jk,Krhs) = vv(ji,jj,jk,Krhs) + zva
[5656]	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
[10989]	508	tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb)
[9031]	509	# if defined key_vertical
[10989]	510	tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kmm)
[9031]	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
[10989]	538	h_out(N_out) = e3v(ji,jj,jk,Kmm)
[9031]	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
[10989]	551	vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)
[9031]	552	# else
[10989]	553	vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
[9031]	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.
[10989]	564	zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm) * fsahm_spt(ji,jj)
	565	hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * vbdiff(ji,jj ,jk) &
	566	& -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * 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.
[10989]	571	zbtr = r1_e1e2f(ji,jj) * e3f(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
[10989]	588	uu(ji,jj,jk,Krhs) = uu(ji,jj,jk,Krhs) &
	589	& - ( rotdiff (ji ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm) ) &
[6140]	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
[10989]	602	vv(ji,jj,jk,Krhs) = vv(ji,jj,jk,Krhs) &
	603	& + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm) ) &
[6140]	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: