Context Navigation

agrif_opa_sponge.F90 @ 8965

Last change on this file since 8965 was 8965, checked in by jchanut, 6 years ago
Add zoom capability in default Gyre, remove bcs velocity update (non reproducible), add option to set clamped obcs - #1965
Property svn:keywords set to `Id`
File size: 16.5 KB

Rev	Line
[3680]	1	#define SPONGE && define SPONGE_TOP
[390]	2
[5656]	3	MODULE agrif_opa_sponge
[6140]	4	!!======================================================================
	5	!! * MODULE agrif_opa_update *
	6	!! AGRIF :
	7	!!======================================================================
	8	!! History :
	9	!!----------------------------------------------------------------------
[7646]	10	#if defined key_agrif
[636]	11	USE par_oce
	12	USE oce
	13	USE dom_oce
	14	USE in_out_manager
[782]	15	USE agrif_oce
[3294]	16	USE wrk_nemo
[5656]	17	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[390]	18
[636]	19	IMPLICIT NONE
	20	PRIVATE
[390]	21
[5656]	22	PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn
	23	PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge
[390]	24
[1156]	25	!!----------------------------------------------------------------------
[6140]	26	!! NEMO/NST 3.7 , NEMO Consortium (2015)
[1156]	27	!! $Id$
[2528]	28	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1156]	29	!!----------------------------------------------------------------------
[5656]	30	CONTAINS
[636]	31
[782]	32	SUBROUTINE Agrif_Sponge_Tra
[636]	33	!!---------------------------------------------
	34	!! * ROUTINE Agrif_Sponge_Tra *
	35	!!---------------------------------------------
	36	REAL(wp) :: timecoeff
[6140]	37	!!---------------------------------------------
	38	!
[390]	39	#if defined SPONGE
[8762]	40	!! timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
	41	!! Assume persistence:
	42	timecoeff = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
[636]	43
[5656]	44	CALL Agrif_Sponge
[390]	45	Agrif_SpecialValue=0.
	46	Agrif_UseSpecialValue = .TRUE.
[5656]	47	tabspongedone_tsn = .FALSE.
[390]	48
[5656]	49	CALL Agrif_Bc_Variable(tsn_sponge_id,calledweight=timecoeff,procname=interptsn_sponge)
[390]	50
[5656]	51	Agrif_UseSpecialValue = .FALSE.
[390]	52	#endif
[6140]	53	!
[636]	54	END SUBROUTINE Agrif_Sponge_Tra
[390]	55
[6140]	56
[782]	57	SUBROUTINE Agrif_Sponge_dyn
[636]	58	!!---------------------------------------------
	59	!! * ROUTINE Agrif_Sponge_dyn *
	60	!!---------------------------------------------
	61	REAL(wp) :: timecoeff
[6140]	62	!!---------------------------------------------
[390]	63
	64	#if defined SPONGE
[8762]	65	!! timecoeff = REAL(Agrif_NbStepint(),wp)/Agrif_rhot()
	66	!! Assume persistence:
	67	timecoeff = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
[390]	68
	69	Agrif_SpecialValue=0.
[782]	70	Agrif_UseSpecialValue = ln_spc_dyn
[390]	71
[5656]	72	tabspongedone_u = .FALSE.
	73	tabspongedone_v = .FALSE.
	74	CALL Agrif_Bc_Variable(un_sponge_id,calledweight=timecoeff,procname=interpun_sponge)
[390]	75
[5656]	76	tabspongedone_u = .FALSE.
	77	tabspongedone_v = .FALSE.
	78	CALL Agrif_Bc_Variable(vn_sponge_id,calledweight=timecoeff,procname=interpvn_sponge)
	79
[390]	80	Agrif_UseSpecialValue = .FALSE.
	81	#endif
[6140]	82	!
[636]	83	END SUBROUTINE Agrif_Sponge_dyn
[390]	84
[6140]	85
[3680]	86	SUBROUTINE Agrif_Sponge
	87	!!---------------------------------------------
	88	!! * ROUTINE Agrif_Sponge *
	89	!!---------------------------------------------
	90	INTEGER :: ji,jj,jk
	91	INTEGER :: ispongearea, ilci, ilcj
[4153]	92	LOGICAL :: ll_spdone
	93	REAL(wp) :: z1spongearea, zramp
	94	REAL(wp), POINTER, DIMENSION(:,:) :: ztabramp
[3680]	95
	96	#if defined SPONGE \|\| defined SPONGE_TOP
[4153]	97	ll_spdone=.TRUE.
	98	IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
	99	! Define ramp from boundaries towards domain interior
	100	! at T-points
	101	! Store it in ztabramp
	102	ll_spdone=.FALSE.
[3680]	103
[4153]	104	CALL wrk_alloc( jpi, jpj, ztabramp )
[3680]	105
[5656]	106	ispongearea = 2 + nn_sponge_len * Agrif_irhox()
[4153]	107	ilci = nlci - ispongearea
	108	ilcj = nlcj - ispongearea
	109	z1spongearea = 1._wp / REAL( ispongearea - 2 )
[3680]	110
[5656]	111	ztabramp(:,:) = 0._wp
[4153]	112
	113	IF( (nbondi == -1) .OR. (nbondi == 2) ) THEN
	114	DO jj = 1, jpj
	115	IF ( umask(2,jj,1) == 1._wp ) THEN
	116	DO ji = 2, ispongearea
	117	ztabramp(ji,jj) = ( ispongearea-ji ) * z1spongearea
	118	END DO
	119	ENDIF
	120	ENDDO
	121	ENDIF
	122
	123	IF( (nbondi == 1) .OR. (nbondi == 2) ) THEN
	124	DO jj = 1, jpj
	125	IF ( umask(nlci-2,jj,1) == 1._wp ) THEN
	126	DO ji = ilci+1,nlci-1
	127	zramp = (ji - (ilci+1) ) * z1spongearea
	128	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
	129	ENDDO
	130	ENDIF
	131	ENDDO
	132	ENDIF
	133
	134	IF( (nbondj == -1) .OR. (nbondj == 2) ) THEN
	135	DO ji = 1, jpi
	136	IF ( vmask(ji,2,1) == 1._wp ) THEN
	137	DO jj = 2, ispongearea
	138	zramp = ( ispongearea-jj ) * z1spongearea
	139	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
	140	END DO
	141	ENDIF
	142	ENDDO
	143	ENDIF
	144
	145	IF( (nbondj == 1) .OR. (nbondj == 2) ) THEN
	146	DO ji = 1, jpi
	147	IF ( vmask(ji,nlcj-2,1) == 1._wp ) THEN
	148	DO jj = ilcj+1,nlcj-1
	149	zramp = (jj - (ilcj+1) ) * z1spongearea
	150	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), zramp )
	151	END DO
	152	ENDIF
	153	ENDDO
	154	ENDIF
	155
	156	ENDIF
	157
[3680]	158	! Tracers
	159	IF( .NOT. spongedoneT ) THEN
[5656]	160	fsaht_spu(:,:) = 0._wp
	161	fsaht_spv(:,:) = 0._wp
	162	DO jj = 2, jpjm1
	163	DO ji = 2, jpim1 ! vector opt.
	164	fsaht_spu(ji,jj) = 0.5_wp * visc_tra * (ztabramp(ji,jj) + ztabramp(ji+1,jj ))
	165	fsaht_spv(ji,jj) = 0.5_wp * visc_tra * (ztabramp(ji,jj) + ztabramp(ji ,jj+1))
	166	END DO
	167	END DO
[3680]	168
[5656]	169	CALL lbc_lnk( fsaht_spu, 'U', 1. ) ! Lateral boundary conditions
	170	CALL lbc_lnk( fsaht_spv, 'V', 1. )
[3680]	171	spongedoneT = .TRUE.
	172	ENDIF
	173
	174	! Dynamics
	175	IF( .NOT. spongedoneU ) THEN
[5656]	176	fsahm_spt(:,:) = 0._wp
	177	fsahm_spf(:,:) = 0._wp
	178	DO jj = 2, jpjm1
	179	DO ji = 2, jpim1 ! vector opt.
	180	fsahm_spt(ji,jj) = visc_dyn * ztabramp(ji,jj)
	181	fsahm_spf(ji,jj) = 0.25_wp * visc_dyn * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1) &
	182	& +ztabramp(ji,jj) + ztabramp(ji+1,jj ) )
	183	END DO
	184	END DO
[3680]	185
[5656]	186	CALL lbc_lnk( fsahm_spt, 'T', 1. ) ! Lateral boundary conditions
	187	CALL lbc_lnk( fsahm_spf, 'F', 1. )
[3680]	188	spongedoneU = .TRUE.
	189	ENDIF
	190	!
[4153]	191	IF (.NOT.ll_spdone) CALL wrk_dealloc( jpi, jpj, ztabramp )
[3680]	192	!
	193	#endif
[6140]	194	!
[3680]	195	END SUBROUTINE Agrif_Sponge
	196
[6140]	197
[5656]	198	SUBROUTINE interptsn_sponge(tabres,i1,i2,j1,j2,k1,k2,n1,n2,before)
[636]	199	!!---------------------------------------------
[5656]	200	!! * ROUTINE interptsn_sponge *
[636]	201	!!---------------------------------------------
[3294]	202	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,n1,n2
	203	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
[5656]	204	LOGICAL, INTENT(in) :: before
[6140]	205	!
[5656]	206	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	207	INTEGER :: iku, ikv
	208	REAL(wp) :: ztsa, zabe1, zabe2, zbtr
	209	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ztu, ztv
	210	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2) ::tsbdiff
	211	!
[6140]	212	IF( before ) THEN
[8762]	213	tabres(i1:i2,j1:j2,k1:k2,n1:n2) = tsb(i1:i2,j1:j2,k1:k2,n1:n2)
[5656]	214	ELSE
[6140]	215	!
[5656]	216	tsbdiff(:,:,:,:) = tsb(i1:i2,j1:j2,:,:) - tabres(:,:,:,:)
	217	DO jn = 1, jpts
	218	DO jk = 1, jpkm1
	219	DO jj = j1,j2-1
	220	DO ji = i1,i2-1
[6140]	221	zabe1 = fsaht_spu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u_n(ji,jj,jk)
	222	zabe2 = fsaht_spv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v_n(ji,jj,jk)
[5656]	223	ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) )
	224	ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
[6140]	225	END DO
	226	END DO
	227	!
[5656]	228	IF( ln_zps ) THEN ! set gradient at partial step level
	229	DO jj = j1,j2-1
	230	DO ji = i1,i2-1
	231	! last level
	232	iku = mbku(ji,jj)
	233	ikv = mbkv(ji,jj)
[6140]	234	IF( iku == jk ) ztu(ji,jj,jk) = 0._wp
	235	IF( ikv == jk ) ztv(ji,jj,jk) = 0._wp
[5656]	236	END DO
	237	END DO
	238	ENDIF
[6140]	239	END DO
	240	!
[5656]	241	DO jk = 1, jpkm1
	242	DO jj = j1+1,j2-1
	243	DO ji = i1+1,i2-1
	244	IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
[6140]	245	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
[5656]	246	! horizontal diffusive trends
	247	ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) )
	248	! add it to the general tracer trends
	249	tsa(ji,jj,jk,jn) = tsa(ji,jj,jk,jn) + ztsa
	250	ENDIF
[6140]	251	END DO
	252	END DO
	253	END DO
	254	!
	255	END DO
	256	!
[5656]	257	tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
[6140]	258	!
[5656]	259	ENDIF
[6140]	260	!
[5656]	261	END SUBROUTINE interptsn_sponge
	262
[6140]	263
[5656]	264	SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2, before)
[636]	265	!!---------------------------------------------
[5656]	266	!! * ROUTINE interpun_sponge *
	267	!!---------------------------------------------
[636]	268	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
	269	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
[5656]	270	LOGICAL, INTENT(in) :: before
[390]	271
[5656]	272	INTEGER :: ji,jj,jk
[390]	273
[5656]	274	! sponge parameters
	275	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
	276	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: ubdiff
	277	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff
	278	INTEGER :: jmax
[6140]	279	!!---------------------------------------------
[5656]	280	!
[6140]	281	IF( before ) THEN
[8965]	282	tabres = ub(i1:i2,j1:j2,:)
[5656]	283	ELSE
	284	ubdiff(i1:i2,j1:j2,:) = (ub(i1:i2,j1:j2,:) - tabres(:,:,:))*umask(i1:i2,j1:j2,:)
[6140]	285	!
[5656]	286	DO jk = 1, jpkm1 ! Horizontal slab
	287	! ! ===============
	288
	289	! ! --------
	290	! Horizontal divergence ! div
	291	! ! --------
	292	DO jj = j1,j2
	293	DO ji = i1+1,i2 ! vector opt.
[6140]	294	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
	295	hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)e3u_n(ji ,jj,jk) ubdiff(ji ,jj,jk) &
	296	& -e2u(ji-1,jj)e3u_n(ji-1,jj,jk) ubdiff(ji-1,jj,jk) ) * zbtr
[5656]	297	END DO
	298	END DO
	299
	300	DO jj = j1,j2-1
	301	DO ji = i1,i2 ! vector opt.
[6140]	302	zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
[5656]	303	rotdiff(ji,jj,jk) = (-e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &
	304	+e1u(ji,jj ) * ubdiff(ji,jj ,jk) &
	305	& ) * fmask(ji,jj,jk) * zbtr
	306	END DO
	307	END DO
[6140]	308	END DO
[5656]	309	!
	310	DO jj = j1+1, j2-1
	311	DO ji = i1+1, i2-1 ! vector opt.
	312
	313	IF (.NOT. tabspongedone_u(ji,jj)) THEN
	314	DO jk = 1, jpkm1 ! Horizontal slab
	315	ze2u = rotdiff (ji,jj,jk)
	316	ze1v = hdivdiff(ji,jj,jk)
	317	! horizontal diffusive trends
[6140]	318	zua = - ( ze2u - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) ) &
[5656]	319	+ ( hdivdiff(ji+1,jj,jk) - ze1v ) / e1u(ji,jj)
	320
	321	! add it to the general momentum trends
	322	ua(ji,jj,jk) = ua(ji,jj,jk) + zua
	323
	324	END DO
	325	ENDIF
	326
	327	END DO
	328	END DO
	329
	330	tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
	331
	332	jmax = j2-1
	333	IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-3)
	334
	335	DO jj = j1+1, jmax
	336	DO ji = i1+1, i2 ! vector opt.
	337
	338	IF (.NOT. tabspongedone_v(ji,jj)) THEN
	339	DO jk = 1, jpkm1 ! Horizontal slab
	340	ze2u = rotdiff (ji,jj,jk)
	341	ze1v = hdivdiff(ji,jj,jk)
	342
	343	! horizontal diffusive trends
[6140]	344	zva = + ( ze2u - rotdiff (ji-1,jj,jk)) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) ) &
[5656]	345	+ ( hdivdiff(ji,jj+1,jk) - ze1v ) / e2v(ji,jj)
	346
	347	! add it to the general momentum trends
	348	va(ji,jj,jk) = va(ji,jj,jk) + zva
	349	END DO
	350	ENDIF
[6140]	351	!
[5656]	352	END DO
	353	END DO
[6140]	354	!
[5656]	355	tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
[6140]	356	!
[5656]	357	ENDIF
[6140]	358	!
[5656]	359	END SUBROUTINE interpun_sponge
	360
	361
	362	SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2, before,nb,ndir)
[636]	363	!!---------------------------------------------
[5656]	364	!! * ROUTINE interpvn_sponge *
	365	!!---------------------------------------------
[636]	366	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2
	367	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2), INTENT(inout) :: tabres
[5656]	368	LOGICAL, INTENT(in) :: before
	369	INTEGER, INTENT(in) :: nb , ndir
[6140]	370	!
	371	INTEGER :: ji, jj, jk
	372	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr
[5656]	373	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: vbdiff
	374	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2) :: rotdiff, hdivdiff
	375	INTEGER :: imax
[6140]	376	!!---------------------------------------------
[5656]	377
[6140]	378	IF( before ) THEN
[8965]	379	tabres = vb(i1:i2,j1:j2,:)
[5656]	380	ELSE
[6140]	381	!
[5656]	382	vbdiff(i1:i2,j1:j2,:) = (vb(i1:i2,j1:j2,:) - tabres(:,:,:))*vmask(i1:i2,j1:j2,:)
[6140]	383	!
[5656]	384	DO jk = 1, jpkm1 ! Horizontal slab
	385	! ! ===============
	386
	387	! ! --------
	388	! Horizontal divergence ! div
	389	! ! --------
	390	DO jj = j1+1,j2
	391	DO ji = i1,i2 ! vector opt.
[6140]	392	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk) * fsahm_spt(ji,jj)
	393	hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v_n(ji,jj ,jk) * vbdiff(ji,jj ,jk) &
	394	& -e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) * vbdiff(ji,jj-1,jk) ) * zbtr
[5656]	395	END DO
	396	END DO
	397	DO jj = j1,j2
	398	DO ji = i1,i2-1 ! vector opt.
[6140]	399	zbtr = r1_e1e2f(ji,jj) * e3f_n(ji,jj,jk) * fsahm_spf(ji,jj)
[5656]	400	rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) &
[6140]	401	& -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr
[5656]	402	END DO
	403	END DO
[6140]	404	END DO
[5656]	405
	406	! ! ===============
	407	!
	408
	409	imax = i2-1
[6140]	410	IF ((nbondi == 1).OR.(nbondi == 2)) imax = MIN(imax,nlci-3)
[5656]	411
	412	DO jj = j1+1, j2
	413	DO ji = i1+1, imax ! vector opt.
[6140]	414	IF( .NOT. tabspongedone_u(ji,jj) ) THEN
	415	DO jk = 1, jpkm1
	416	ua(ji,jj,jk) = ua(ji,jj,jk) &
	417	& - ( rotdiff (ji ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u_n(ji,jj,jk) ) &
	418	& + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj ,jk)) * r1_e1u(ji,jj)
[5656]	419	END DO
	420	ENDIF
	421	END DO
	422	END DO
[6140]	423	!
[5656]	424	tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
[6140]	425	!
[5656]	426	DO jj = j1+1, j2-1
	427	DO ji = i1+1, i2-1 ! vector opt.
[6140]	428	IF( .NOT. tabspongedone_v(ji,jj) ) THEN
	429	DO jk = 1, jpkm1
	430	va(ji,jj,jk) = va(ji,jj,jk) &
	431	& + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v_n(ji,jj,jk) ) &
	432	& + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj)
[5656]	433	END DO
	434	ENDIF
	435	END DO
	436	END DO
	437	tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
	438	ENDIF
[6140]	439	!
[5656]	440	END SUBROUTINE interpvn_sponge
	441
[390]	442	#else
[636]	443	CONTAINS
	444	SUBROUTINE agrif_opa_sponge_empty
	445	!!---------------------------------------------
	446	!! * ROUTINE agrif_OPA_sponge_empty *
	447	!!---------------------------------------------
	448	WRITE(,) 'agrif_opa_sponge : You should not have seen this print! error?'
	449	END SUBROUTINE agrif_opa_sponge_empty
[390]	450	#endif
[636]	451
[6140]	452	!!======================================================================
[636]	453	END MODULE agrif_opa_sponge

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/2017/dev_r8624_AGRIF3_VVL/NEMOGCM/NEMO/NST_SRC/agrif_opa_sponge.F90 @ 8965

Download in other formats: