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/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/NST_SRC – NEMO

Context Navigation

source: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 5006

Last change on this file since 5006 was 4153, checked in by cetlod, 10 years ago
dev_LOCEAN_2013: merge in trunk changes between r3940 and r4028, see ticket #1169
Property svn:keywords set to `Id`
File size: 16.7 KB

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

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

Download in other formats: