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 @ 1156

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

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

Download in other formats: