Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

agrif_opa_sponge.F90 @ 8306

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

Download in other formats: