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 trunk/NEMO/NST_SRC – NEMO

Context Navigation

source: trunk/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 1146

Last change on this file since 1146 was 1146, checked in by rblod, 16 years ago
Add svn Id (first try), see ticket #210
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 15.1 KB

Rev	Line
[390]	1	#define SPONGE
	2
[636]	3	Module agrif_opa_sponge
[393]	4	#if defined key_agrif
[636]	5	USE par_oce
	6	USE oce
	7	USE dom_oce
	8	USE in_out_manager
[782]	9	USE agrif_oce
[390]	10
[636]	11	IMPLICIT NONE
	12	PRIVATE
[390]	13
[636]	14	PUBLIC Agrif_Sponge_Tra, Agrif_Sponge_Dyn, interptn, interpsn, interpun, interpvn
[390]	15
	16
[636]	17	CONTAINS
	18
[782]	19	SUBROUTINE Agrif_Sponge_Tra
[636]	20	!!---------------------------------------------
	21	!! * ROUTINE Agrif_Sponge_Tra *
	22	!!---------------------------------------------
	23	#include "domzgr_substitute.h90"
	24
[390]	25	INTEGER :: ji,jj,jk
	26	INTEGER :: spongearea
[636]	27	REAL(wp) :: timecoeff
	28	REAL(wp) :: zta, zsa, zabe1, zabe2, zbtr
[390]	29	REAL(wp), DIMENSION(jpi,jpj) :: localviscsponge
[782]	30	REAL(wp), DIMENSION(jpi,jpj,jpk) :: tbdiff, sbdiff
[390]	31	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztu ,ztv, zsu ,zsv
[782]	32	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztab
[390]	33
	34	#if defined SPONGE
	35
[636]	36	timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
	37
[390]	38	Agrif_SpecialValue=0.
	39	Agrif_UseSpecialValue = .TRUE.
[636]	40	ztab = 0.e0
	41	CALL Agrif_Bc_Variable(ztab, ta,calledweight=timecoeff,procname=interptn)
[390]	42	Agrif_UseSpecialValue = .FALSE.
	43
[636]	44	tbdiff(:,:,:) = tb(:,:,:) - ztab(:,:,:)
[390]	45
[636]	46	ztab = 0.e0
[390]	47	Agrif_SpecialValue=0.
	48	Agrif_UseSpecialValue = .TRUE.
[636]	49	CALL Agrif_Bc_Variable(ztab, sa,calledweight=timecoeff,procname=interpsn)
[390]	50	Agrif_UseSpecialValue = .FALSE.
	51
[636]	52	sbdiff(:,:,:) = sb(:,:,:) - ztab(:,:,:)
[390]	53
[782]	54
[390]	55	spongearea = 2 + 2 * Agrif_irhox()
	56
	57	localviscsponge = 0.
[782]	58
	59	IF (.NOT. spongedoneT) THEN
	60	spe1ur(:,:) = 0.
	61	spe2vr(:,:) = 0.
[390]	62
	63	IF ((nbondi == -1).OR.(nbondi == 2)) THEN
[636]	64	DO ji = 2, spongearea
	65	localviscsponge(ji,:) = visc_tra * (spongearea-ji)/real(spongearea-2)
	66	ENDDO
[782]	67
	68	spe1ur(2:spongearea-1,:)=0.5 * (localviscsponge(2:spongearea-1,:) + localviscsponge(3:spongearea,:)) &
	69	* e2u(2:spongearea-1,:) / e1u(2:spongearea-1,:)
	70
	71	spe2vr(2:spongearea,1:jpjm1) = 0.5 * (localviscsponge(2:spongearea,1:jpjm1) + &
	72	localviscsponge(2:spongearea,2:jpj)) &
	73	* e1v(2:spongearea,1:jpjm1) / e2v(2:spongearea,1:jpjm1)
[390]	74	ENDIF
	75
	76	IF ((nbondi == 1).OR.(nbondi == 2)) THEN
[636]	77	DO ji = nlci-spongearea + 1,nlci-1
	78	localviscsponge(ji,:) = visc_tra * (ji - (nlci-spongearea+1))/real(spongearea-2)
	79	ENDDO
[782]	80
	81	spe1ur(nlci-spongearea + 1:nlci-2,:)=0.5 * (localviscsponge(nlci-spongearea + 1:nlci-2,:) + &
	82	localviscsponge(nlci-spongearea + 2:nlci-1,:)) &
	83	* e2u(nlci-spongearea + 1:nlci-2,:) / e1u(nlci-spongearea + 1:nlci-2,:)
	84
	85	spe2vr(nlci-spongearea + 1:nlci-1,1:jpjm1) = 0.5 * (localviscsponge(nlci-spongearea + 1:nlci-1,1:jpjm1) &
	86	+ localviscsponge(nlci-spongearea + 1:nlci-1,2:jpj)) &
	87	* e1v(nlci-spongearea + 1:nlci-1,1:jpjm1) / e2v(nlci-spongearea + 1:nlci-1,1:jpjm1)
[390]	88	ENDIF
	89
	90
	91	IF ((nbondj == -1).OR.(nbondj == 2)) THEN
[636]	92	DO jj = 2, spongearea
	93	localviscsponge(:,jj) = visc_tra * (spongearea-jj)/real(spongearea-2)
	94	ENDDO
[782]	95
	96	spe1ur(1:jpim1,2:spongearea)=0.5 * (localviscsponge(1:jpim1,2:spongearea) + &
	97	localviscsponge(2:jpi,2:spongearea)) &
	98	* e2u(1:jpim1,2:spongearea) / e1u(1:jpim1,2:spongearea)
	99
	100	spe2vr(:,2:spongearea-1) = 0.5 * (localviscsponge(:,2:spongearea-1) + &
	101	localviscsponge(:,3:spongearea)) &
	102	* e1v(:,2:spongearea-1) / e2v(:,2:spongearea-1)
[390]	103	ENDIF
	104
	105	IF ((nbondj == 1).OR.(nbondj == 2)) THEN
[636]	106	DO jj = nlcj-spongearea + 1,nlcj-1
	107	localviscsponge(:,jj) = visc_tra * (jj - (nlcj-spongearea+1))/real(spongearea-2)
	108	ENDDO
[782]	109
	110	spe1ur(1:jpim1,nlcj-spongearea + 1:nlcj-1)=0.5 * (localviscsponge(1:jpim1,nlcj-spongearea + 1:nlcj-1) + &
	111	localviscsponge(2:jpi,nlcj-spongearea + 1:nlcj-1)) &
	112	* e2u(1:jpim1,nlcj-spongearea + 1:nlcj-1) / e1u(1:jpim1,nlcj-spongearea + 1:nlcj-1)
	113
	114	spe2vr(:,nlcj-spongearea + 1:nlcj-2) = 0.5 * (localviscsponge(:,nlcj-spongearea + 1:nlcj-2) + &
	115	localviscsponge(:,nlcj-spongearea + 2:nlcj-1)) &
	116	* e1v(:,nlcj-spongearea + 1:nlcj-2) / e2v(:,nlcj-spongearea + 1:nlcj-2)
[390]	117	ENDIF
[782]	118
	119	spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:))
[390]	120
	121	spongedoneT = .TRUE.
	122	ENDIF
	123
[636]	124	DO jk = 1, jpkm1
	125	DO jj = 1, jpjm1
[390]	126	DO ji = 1, jpim1
[469]	127	#if defined key_zco
[782]	128	zabe1 = umask(ji,jj,jk) * spe1ur(ji,jj)
	129	zabe2 = vmask(ji,jj,jk) * spe2vr(ji,jj)
[469]	130	#else
[782]	131	zabe1 = umask(ji,jj,jk) * spe1ur(ji,jj) * fse3u(ji,jj,jk)
	132	zabe2 = vmask(ji,jj,jk) * spe2vr(ji,jj) * fse3v(ji,jj,jk)
[390]	133	#endif
	134	ztu(ji,jj,jk) = zabe1 * ( tbdiff(ji+1,jj ,jk) - tbdiff(ji,jj,jk) )
	135	zsu(ji,jj,jk) = zabe1 * ( sbdiff(ji+1,jj ,jk) - sbdiff(ji,jj,jk) )
	136	ztv(ji,jj,jk) = zabe2 * ( tbdiff(ji ,jj+1,jk) - tbdiff(ji,jj,jk) )
	137	zsv(ji,jj,jk) = zabe2 * ( sbdiff(ji ,jj+1,jk) - sbdiff(ji,jj,jk) )
	138	ENDDO
[636]	139	ENDDO
[390]	140
	141	DO jj = 2,jpjm1
	142	DO ji = 2,jpim1
[469]	143	#if defined key_zco
[782]	144	zbtr = spbtr2(ji,jj)
[469]	145	#else
[782]	146	zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk)
[390]	147	#endif
	148	! horizontal diffusive trends
	149	zta = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) &
	150	& + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) )
	151	zsa = zbtr * ( zsu(ji,jj,jk) - zsu(ji-1,jj,jk) &
	152	& + zsv(ji,jj,jk) - zsv(ji,jj-1,jk) )
	153	! add it to the general tracer trends
	154	ta(ji,jj,jk) = (ta(ji,jj,jk) + zta)
	155	sa(ji,jj,jk) = (sa(ji,jj,jk) + zsa)
	156	END DO
	157	END DO
	158
[636]	159	ENDDO
[390]	160
	161	#endif
	162
[636]	163	END SUBROUTINE Agrif_Sponge_Tra
[390]	164
[782]	165	SUBROUTINE Agrif_Sponge_dyn
[636]	166	!!---------------------------------------------
	167	!! * ROUTINE Agrif_Sponge_dyn *
	168	!!---------------------------------------------
	169	#include "domzgr_substitute.h90"
[390]	170
	171	INTEGER :: ji,jj,jk
	172	INTEGER :: spongearea
[636]	173	REAL(wp) :: timecoeff
[782]	174	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
[390]	175	REAL(wp), DIMENSION(jpi,jpj) :: localviscsponge
[636]	176	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztab, ubdiff, vbdiff,rotdiff,hdivdiff
[390]	177
	178	#if defined SPONGE
	179
[636]	180	timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
[390]	181
	182	Agrif_SpecialValue=0.
[782]	183	Agrif_UseSpecialValue = ln_spc_dyn
[636]	184	ztab = 0.e0
	185	CALL Agrif_Bc_Variable(ztab, ua,calledweight=timecoeff,procname=interpun)
[390]	186	Agrif_UseSpecialValue = .FALSE.
	187
[782]	188	ubdiff(:,:,:) = (ub(:,:,:) - ztab(:,:,:))*umask(:,:,:)
[390]	189
[636]	190	ztab = 0.e0
[390]	191	Agrif_SpecialValue=0.
[782]	192	Agrif_UseSpecialValue = ln_spc_dyn
[636]	193	CALL Agrif_Bc_Variable(ztab, va,calledweight=timecoeff,procname=interpvn)
[390]	194	Agrif_UseSpecialValue = .FALSE.
	195
[782]	196	vbdiff(:,:,:) = (vb(:,:,:) - ztab(:,:,:))*vmask(:,:,:)
[390]	197
	198	spongearea = 2 + 2 * Agrif_irhox()
	199
	200	localviscsponge = 0.
	201
[782]	202	IF (.NOT. spongedoneU) THEN
	203	spe1ur2(:,:) = 0.
	204	spe2vr2(:,:) = 0.
	205
[390]	206	IF ((nbondi == -1).OR.(nbondi == 2)) THEN
[636]	207	DO ji = 2, spongearea
	208	localviscsponge(ji,:) = visc_dyn * (spongearea-ji)/real(spongearea-2)
	209	ENDDO
[782]	210
	211	spe1ur2(2:spongearea-1,:)=0.5 * (localviscsponge(2:spongearea-1,:) + localviscsponge(3:spongearea,:))
	212
	213	spe2vr2(2:spongearea,1:jpjm1) = 0.5 * (localviscsponge(2:spongearea,1:jpjm1) + &
	214	localviscsponge(2:spongearea,2:jpj))
[390]	215	ENDIF
	216
	217	IF ((nbondi == 1).OR.(nbondi == 2)) THEN
[636]	218	DO ji = nlci-spongearea + 1,nlci-1
	219	localviscsponge(ji,:) = visc_dyn * (ji - (nlci-spongearea+1))/real(spongearea-2)
	220	ENDDO
[782]	221
	222	spe1ur2(nlci-spongearea + 1:nlci-2,:)=0.5 * (localviscsponge(nlci-spongearea + 1:nlci-2,:) + &
	223	localviscsponge(nlci-spongearea + 2:nlci-1,:))
	224
	225	spe2vr2(nlci-spongearea + 1:nlci-1,1:jpjm1) = 0.5 * (localviscsponge(nlci-spongearea + 1:nlci-1,1:jpjm1) &
	226	+ localviscsponge(nlci-spongearea + 1:nlci-1,2:jpj))
[390]	227	ENDIF
	228
[782]	229
[390]	230	IF ((nbondj == -1).OR.(nbondj == 2)) THEN
[636]	231	DO jj = 2, spongearea
	232	localviscsponge(:,jj) = visc_dyn * (spongearea-jj)/real(spongearea-2)
	233	ENDDO
[782]	234
	235	spe1ur2(1:jpim1,2:spongearea)=0.5 * (localviscsponge(1:jpim1,2:spongearea) + &
	236	localviscsponge(2:jpi,2:spongearea))
	237
	238	spe2vr2(:,2:spongearea-1) = 0.5 * (localviscsponge(:,2:spongearea-1) + &
	239	localviscsponge(:,3:spongearea))
[390]	240	ENDIF
	241
	242	IF ((nbondj == 1).OR.(nbondj == 2)) THEN
[636]	243	DO jj = nlcj-spongearea + 1,nlcj-1
	244	localviscsponge(:,jj) = visc_dyn * (jj - (nlcj-spongearea+1))/real(spongearea-2)
	245	ENDDO
[782]	246
	247	spe1ur2(1:jpim1,nlcj-spongearea + 1:nlcj-1)=0.5 * (localviscsponge(1:jpim1,nlcj-spongearea + 1:nlcj-1) + &
	248	localviscsponge(2:jpi,nlcj-spongearea + 1:nlcj-1))
	249
	250	spe2vr2(:,nlcj-spongearea + 1:nlcj-2) = 0.5 * (localviscsponge(:,nlcj-spongearea + 1:nlcj-2) + &
	251	localviscsponge(:,nlcj-spongearea + 2:nlcj-1))
[636]	252	ENDIF
[390]	253
[782]	254	spongedoneU = .TRUE.
	255
	256	spbtr3(:,:) = 1./( e1f(:,:) * e2f(:,:))
	257	ENDIF
	258
	259	IF (.NOT. spongedoneT) THEN
	260	spbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:))
	261	ENDIF
	262
	263	DO jk=1,jpkm1
	264	ubdiff(:,:,jk) = ubdiff(:,:,jk) * spe1ur2(:,:)
	265	vbdiff(:,:,jk) = vbdiff(:,:,jk) * spe2vr2(:,:)
	266	ENDDO
	267
[390]	268	hdivdiff = 0.
	269	rotdiff = 0.
	270
	271	DO jk = 1, jpkm1 ! Horizontal slab
	272	! ! ===============
	273
	274	! ! --------
	275	! Horizontal divergence ! div
	276	! ! --------
	277	DO jj = 2, jpjm1
	278	DO ji = 2, jpim1 ! vector opt.
[469]	279	#if defined key_zco
[782]	280	zbtr = spbtr2(ji,jj)
[469]	281	hdivdiff(ji,jj,jk) = ( e2u(ji,jj) * ubdiff(ji,jj,jk) &
[636]	282	- e2u(ji-1,jj ) * ubdiff(ji-1,jj ,jk) &
	283	& + e1v(ji,jj) * vbdiff(ji,jj,jk) - &
[782]	284	& e1v(ji ,jj-1) * vbdiff(ji ,jj-1,jk) ) * zbtr
[469]	285	#else
[782]	286	zbtr = spbtr2(ji,jj) / fse3t(ji,jj,jk)
[390]	287	hdivdiff(ji,jj,jk) = &
	288	( e2u(ji,jj)fse3u(ji,jj,jk) &
	289	ubdiff(ji,jj,jk) - e2u(ji-1,jj )* &
	290	fse3u(ji-1,jj ,jk) * ubdiff(ji-1,jj ,jk) &
[636]	291	+ e1v(ji,jj)fse3v(ji,jj,jk) &
[390]	292	vbdiff(ji,jj,jk) - e1v(ji ,jj-1)* &
[782]	293	fse3v(ji ,jj-1,jk) * vbdiff(ji ,jj-1,jk) ) * zbtr
[390]	294	#endif
	295	END DO
	296	END DO
[636]	297
[390]	298	DO jj = 1, jpjm1
	299	DO ji = 1, jpim1 ! vector opt.
[782]	300	#if defined key_zco
	301	zbtr = spbtr3(ji,jj)
[390]	302	rotdiff(ji,jj,jk) = ( e2v(ji+1,jj ) * vbdiff(ji+1,jj ,jk) - e2v(ji,jj) * vbdiff(ji,jj,jk) &
	303	& - e1u(ji ,jj+1) * ubdiff(ji ,jj+1,jk) + e1u(ji,jj) * ubdiff(ji,jj,jk) ) &
[782]	304	& * fmask(ji,jj,jk) * zbtr
	305	#else
	306	zbtr = spbtr3(ji,jj) * fse3f(ji,jj,jk)
	307	rotdiff(ji,jj,jk) = ( e2v(ji+1,jj ) * vbdiff(ji+1,jj ,jk) - e2v(ji,jj) * vbdiff(ji,jj,jk) &
	308	& - e1u(ji ,jj+1) * ubdiff(ji ,jj+1,jk) + e1u(ji,jj) * ubdiff(ji,jj,jk) ) &
	309	& * fmask(ji,jj,jk) * zbtr
	310	#endif
[390]	311	END DO
	312	END DO
[636]	313
	314	ENDDO
	315
[390]	316	! ! ===============
	317	DO jk = 1, jpkm1 ! Horizontal slab
	318	! ! ===============
	319	DO jj = 2, jpjm1
	320	DO ji = 2, jpim1 ! vector opt.
[469]	321	#if defined key_zco
	322	! horizontal diffusive trends
	323	ze2u = rotdiff (ji,jj,jk)
	324	ze1v = hdivdiff(ji,jj,jk)
	325	zua = - ( ze2u - &
[636]	326	rotdiff (ji,jj-1,jk) ) / e2u(ji,jj) &
	327	+ ( hdivdiff(ji+1,jj,jk) - &
	328	ze1v ) / e1u(ji,jj)
[469]	329
	330	zva = + ( ze2u - &
[636]	331	rotdiff (ji-1,jj,jk) ) / e1v(ji,jj) &
	332	+ ( hdivdiff(ji,jj+1,jk) - &
	333	ze1v ) / e2v(ji,jj)
[469]	334	#else
[782]	335	ze2u = rotdiff (ji,jj,jk)
[390]	336	ze1v = hdivdiff(ji,jj,jk)
	337	! horizontal diffusive trends
[782]	338	zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * fse3u(ji,jj,jk) ) &
[636]	339	+ ( hdivdiff(ji+1,jj,jk) - ze1v &
	340	) / e1u(ji,jj)
[390]	341
[782]	342	zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * fse3v(ji,jj,jk) ) &
[636]	343	+ ( hdivdiff(ji,jj+1,jk) - ze1v &
	344	) / e2v(ji,jj)
[390]	345	#endif
	346
	347	! add it to the general momentum trends
	348	ua(ji,jj,jk) = ua(ji,jj,jk) + zua
	349	va(ji,jj,jk) = va(ji,jj,jk) + zva
	350	END DO
	351	END DO
	352	! ! ===============
	353	END DO ! End of slab
	354	! ! ===============
	355
	356	#endif
	357
[636]	358	END SUBROUTINE Agrif_Sponge_dyn
[390]	359
[636]	360	SUBROUTINE interptn(tabres,i1,i2,j1,j2,k1,k2)
	361	!!---------------------------------------------
	362	!! * ROUTINE interptn *
	363	!!---------------------------------------------
[390]	364	# include "domzgr_substitute.h90"
[636]	365
	366	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
	367	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
[390]	368
[636]	369	tabres(i1:i2,j1:j2,k1:k2) = tn(i1:i2,j1:j2,k1:k2)
[390]	370
[636]	371	END SUBROUTINE interptn
	372
	373	SUBROUTINE interpsn(tabres,i1,i2,j1,j2,k1,k2)
	374	!!---------------------------------------------
	375	!! * ROUTINE interpsn *
	376	!!---------------------------------------------
[390]	377	# include "domzgr_substitute.h90"
[636]	378
	379	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
	380	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
[390]	381
[636]	382	tabres(i1:i2,j1:j2,k1:k2) = sn(i1:i2,j1:j2,k1:k2)
[390]	383
[636]	384	END SUBROUTINE interpsn
	385
	386	SUBROUTINE interpun(tabres,i1,i2,j1,j2,k1,k2)
	387	!!---------------------------------------------
	388	!! * ROUTINE interpun *
	389	!!---------------------------------------------
[390]	390	# include "domzgr_substitute.h90"
[636]	391
	392	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
	393	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
[390]	394
[636]	395	tabres(i1:i2,j1:j2,k1:k2) = un(i1:i2,j1:j2,k1:k2)
[390]	396
[636]	397	END SUBROUTINE interpun
	398
	399	SUBROUTINE interpvn(tabres,i1,i2,j1,j2,k1,k2)
	400	!!---------------------------------------------
	401	!! * ROUTINE interpvn *
	402	!!---------------------------------------------
[390]	403	# include "domzgr_substitute.h90"
[636]	404
	405	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
	406	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
[390]	407
[636]	408	tabres(i1:i2,j1:j2,k1:k2) = vn(i1:i2,j1:j2,k1:k2)
[390]	409
[636]	410	END SUBROUTINE interpvn
[390]	411
	412	#else
[636]	413	CONTAINS
	414
	415	SUBROUTINE agrif_opa_sponge_empty
	416	!!---------------------------------------------
	417	!! * ROUTINE agrif_OPA_sponge_empty *
	418	!!---------------------------------------------
	419	WRITE(,) 'agrif_opa_sponge : You should not have seen this print! error?'
	420	END SUBROUTINE agrif_opa_sponge_empty
[390]	421	#endif
[636]	422
	423	END MODULE agrif_opa_sponge

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

Download in other formats: