Context Navigation

agrif_oce_sponge.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

Property svn:keywords set to Id

File size: 32.4 KB

Rev	Line
[3680]	1	#define SPONGE && define SPONGE_TOP
[390]	2
[9570]	3	MODULE agrif_oce_sponge
[6140]	4	!!======================================================================
[9570]	5	!! * MODULE agrif_oce_interp *
[9019]	6	!! AGRIF: sponge package for the ocean dynamics (OPA)
[6140]	7	!!======================================================================
[9019]	8	!! History : 2.0 ! 2002-06 (XXX) Original cade
	9	!! - ! 2005-11 (XXX)
	10	!! 3.2 ! 2009-04 (R. Benshila)
	11	!! 3.6 ! 2014-09 (R. Benshila)
[6140]	12	!!----------------------------------------------------------------------
[7646]	13	#if defined key_agrif
[9019]	14	!!----------------------------------------------------------------------
	15	!! 'key_agrif' AGRIF zoom
	16	!!----------------------------------------------------------------------
[636]	17	USE par_oce
	18	USE oce
	19	USE dom_oce
[9019]	20	!
[636]	21	USE in_out_manager
[782]	22	USE agrif_oce
[5656]	23	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[12229]	24	USE iom
	25	USE vremap
[390]	26
[636]	27	IMPLICIT NONE
	28	PRIVATE
[390]	29
[5656]	30	PUBLIC Agrif_Sponge, Agrif_Sponge_Tra, Agrif_Sponge_Dyn
	31	PUBLIC interptsn_sponge, interpun_sponge, interpvn_sponge
[390]	32
[12340]	33	!! * Substitutions
	34	# include "do_loop_substitute.h90"
[1156]	35	!!----------------------------------------------------------------------
[9598]	36	!! NEMO/NST 4.0 , NEMO Consortium (2018)
[1156]	37	!! $Id$
[10068]	38	!! Software governed by the CeCILL license (see ./LICENSE)
[1156]	39	!!----------------------------------------------------------------------
[5656]	40	CONTAINS
[636]	41
[782]	42	SUBROUTINE Agrif_Sponge_Tra
[9019]	43	!!----------------------------------------------------------------------
	44	!! * ROUTINE Agrif_Sponge_Tra *
	45	!!----------------------------------------------------------------------
	46	REAL(wp) :: zcoef ! local scalar
[9056]	47
[9019]	48	!!----------------------------------------------------------------------
[6140]	49	!
[390]	50	#if defined SPONGE
[9031]	51	!! Assume persistence:
[9056]	52	zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
[9031]	53
[5656]	54	CALL Agrif_Sponge
[9019]	55	Agrif_SpecialValue = 0._wp
[390]	56	Agrif_UseSpecialValue = .TRUE.
[9019]	57	tabspongedone_tsn = .FALSE.
	58	!
	59	CALL Agrif_Bc_Variable( tsn_sponge_id, calledweight=zcoef, procname=interptsn_sponge )
	60	!
[5656]	61	Agrif_UseSpecialValue = .FALSE.
[390]	62	#endif
[6140]	63	!
[12229]	64	CALL iom_put( 'agrif_spu', fspu(:,:))
	65	CALL iom_put( 'agrif_spv', fspv(:,:))
	66	!
[636]	67	END SUBROUTINE Agrif_Sponge_Tra
[390]	68
[6140]	69
[782]	70	SUBROUTINE Agrif_Sponge_dyn
[9019]	71	!!----------------------------------------------------------------------
	72	!! * ROUTINE Agrif_Sponge_dyn *
	73	!!----------------------------------------------------------------------
	74	REAL(wp) :: zcoef ! local scalar
	75	!!----------------------------------------------------------------------
	76	!
[390]	77	#if defined SPONGE
[9056]	78	zcoef = REAL(Agrif_rhot()-1,wp)/REAL(Agrif_rhot())
[9031]	79
	80	Agrif_SpecialValue=0.
[782]	81	Agrif_UseSpecialValue = ln_spc_dyn
[9019]	82	!
[5656]	83	tabspongedone_u = .FALSE.
	84	tabspongedone_v = .FALSE.
[9019]	85	CALL Agrif_Bc_Variable( un_sponge_id, calledweight=zcoef, procname=interpun_sponge )
	86	!
[5656]	87	tabspongedone_u = .FALSE.
	88	tabspongedone_v = .FALSE.
[9019]	89	CALL Agrif_Bc_Variable( vn_sponge_id, calledweight=zcoef, procname=interpvn_sponge )
	90	!
[390]	91	Agrif_UseSpecialValue = .FALSE.
	92	#endif
[6140]	93	!
[12229]	94	CALL iom_put( 'agrif_spt', fspt(:,:))
	95	CALL iom_put( 'agrif_spf', fspf(:,:))
	96	!
[636]	97	END SUBROUTINE Agrif_Sponge_dyn
[390]	98
[6140]	99
[3680]	100	SUBROUTINE Agrif_Sponge
[9019]	101	!!----------------------------------------------------------------------
	102	!! * ROUTINE Agrif_Sponge *
	103	!!----------------------------------------------------------------------
	104	INTEGER :: ji, jj, ind1, ind2
[12229]	105	INTEGER :: ispongearea, jspongearea
	106	REAL(wp) :: z1_ispongearea, z1_jspongearea
[9019]	107	REAL(wp), DIMENSION(jpi,jpj) :: ztabramp
[12229]	108	REAL(wp), DIMENSION(jpjmax) :: zmskwest, zmskeast
	109	REAL(wp), DIMENSION(jpimax) :: zmsknorth, zmsksouth
[9019]	110	!!----------------------------------------------------------------------
	111	!
[12229]	112	! Sponge 1d example with:
	113	! iraf = 3 ; nbghost = 3 ; nn_sponge_len = 2
	114	!
	115	!coarse : U T U T U T U
	116	!\| \| \| \| \|
	117	!fine : t u t u t u t u t u t u t u t u t u t u t
	118	!sponge val:0 0 0 1 5/6 4/6 3/6 2/6 1/6 0 0
	119	! \| ghost \| <-- sponge area -- > \|
	120	! \| points \| \|
	121	! \|--> dynamical interface
	122
[3680]	123	#if defined SPONGE \|\| defined SPONGE_TOP
[4153]	124	IF (( .NOT. spongedoneT ).OR.( .NOT. spongedoneU )) THEN
[12229]	125	!
	126	! Retrieve masks at open boundaries:
	127
	128	! --- West --- !
	129	ztabramp(:,:) = 0._wp
	130	ind1 = 1+nbghostcells
	131	DO ji = mi0(ind1), mi1(ind1)
	132	ztabramp(ji,:) = ssumask(ji,:)
	133	END DO
	134	!
	135	zmskwest(:) = 0._wp
	136	zmskwest(1:jpj) = MAXVAL(ztabramp(:,:), dim=1)
	137
	138	! --- East --- !
	139	ztabramp(:,:) = 0._wp
	140	ind1 = jpiglo - nbghostcells - 1
	141	DO ji = mi0(ind1), mi1(ind1)
	142	ztabramp(ji,:) = ssumask(ji,:)
	143	END DO
	144	!
	145	zmskeast(:) = 0._wp
	146	zmskeast(1:jpj) = MAXVAL(ztabramp(:,:), dim=1)
	147
	148	! --- South --- !
	149	ztabramp(:,:) = 0._wp
	150	ind1 = 1+nbghostcells
	151	DO jj = mj0(ind1), mj1(ind1)
	152	ztabramp(:,jj) = ssvmask(:,jj)
	153	END DO
	154	!
	155	zmsksouth(:) = 0._wp
	156	zmsksouth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2)
	157
	158	! --- North --- !
	159	ztabramp(:,:) = 0._wp
	160	ind1 = jpjglo - nbghostcells - 1
	161	DO jj = mj0(ind1), mj1(ind1)
	162	ztabramp(:,jj) = ssvmask(:,jj)
	163	END DO
	164	!
	165	zmsknorth(:) = 0._wp
	166	zmsknorth(1:jpi) = MAXVAL(ztabramp(:,:), dim=2)
	167	! JC: SPONGE MASKING TO BE SORTED OUT:
	168	zmskwest(:) = 1._wp
	169	zmskeast(:) = 1._wp
	170	zmsknorth(:) = 1._wp
	171	zmsksouth(:) = 1._wp
	172	#if defined key_mpp_mpi
	173	! CALL mpp_max( 'AGRIF_sponge', zmskwest(:) , jpjmax )
	174	! CALL mpp_max( 'AGRIF_Sponge', zmskeast(:) , jpjmax )
	175	! CALL mpp_max( 'AGRIF_Sponge', zmsksouth(:), jpimax )
	176	! CALL mpp_max( 'AGRIF_Sponge', zmsknorth(:), jpimax )
	177	#endif
	178
[9019]	179	! Define ramp from boundaries towards domain interior at T-points
[4153]	180	! Store it in ztabramp
[3680]	181
[12229]	182	ispongearea = nn_sponge_len * Agrif_irhox()
	183	z1_ispongearea = 1._wp / REAL( ispongearea )
	184	jspongearea = nn_sponge_len * Agrif_irhoy()
	185	z1_jspongearea = 1._wp / REAL( jspongearea )
[9019]	186
[5656]	187	ztabramp(:,:) = 0._wp
[4153]	188
[12229]	189	! Trick to remove sponge in 2DV domains:
	190	IF ( nbcellsx <= 3 ) ispongearea = -1
	191	IF ( nbcellsy <= 3 ) jspongearea = -1
	192
[9019]	193	! --- West --- !
[12229]	194	ind1 = 1+nbghostcells
	195	ind2 = 1+nbghostcells + ispongearea
	196	DO ji = mi0(ind1), mi1(ind2)
	197	DO jj = 1, jpj
	198	ztabramp(ji,jj) = REAL( ind2 - mig(ji) ) * z1_ispongearea * zmskwest(jj)
	199	END DO
	200	END DO
	201
	202	! ghost cells:
	203	ind1 = 1
	204	ind2 = nbghostcells + 1
	205	DO ji = mi0(ind1), mi1(ind2)
	206	DO jj = 1, jpj
	207	ztabramp(ji,jj) = zmskwest(jj)
	208	END DO
	209	END DO
	210
	211	! --- East --- !
	212	ind1 = jpiglo - nbghostcells - ispongearea
	213	ind2 = jpiglo - nbghostcells
	214	DO ji = mi0(ind1), mi1(ind2)
[4153]	215	DO jj = 1, jpj
[12229]	216	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mig(ji) - ind1 ) * z1_ispongearea) * zmskeast(jj)
[4153]	217	ENDDO
[12229]	218	END DO
[4153]	219
[12229]	220	! ghost cells:
	221	ind1 = jpiglo - nbghostcells
	222	ind2 = jpiglo
	223	DO ji = mi0(ind1), mi1(ind2)
[4153]	224	DO jj = 1, jpj
[12229]	225	ztabramp(ji,jj) = zmskeast(jj)
[4153]	226	ENDDO
[12229]	227	END DO
[4153]	228
[9019]	229	! --- South --- !
[12229]	230	ind1 = 1+nbghostcells
	231	ind2 = 1+nbghostcells + jspongearea
	232	DO jj = mj0(ind1), mj1(ind2)
	233	DO ji = 1, jpi
	234	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( ind2 - mjg(jj) ) * z1_jspongearea) * zmsksouth(ji)
	235	END DO
	236	END DO
[4153]	237
[12229]	238	! ghost cells:
	239	ind1 = 1
	240	ind2 = nbghostcells + 1
	241	DO jj = mj0(ind1), mj1(ind2)
	242	DO ji = 1, jpi
	243	ztabramp(ji,jj) = zmsksouth(ji)
	244	END DO
	245	END DO
	246
[9019]	247	! --- North --- !
[12229]	248	ind1 = jpjglo - nbghostcells - jspongearea
	249	ind2 = jpjglo - nbghostcells
	250	DO jj = mj0(ind1), mj1(ind2)
	251	DO ji = 1, jpi
	252	ztabramp(ji,jj) = MAX( ztabramp(ji,jj), REAL( mjg(jj) - ind1 ) * z1_jspongearea) * zmsknorth(ji)
	253	END DO
	254	END DO
[4153]	255
[12229]	256	! ghost cells:
	257	ind1 = jpjglo - nbghostcells
	258	ind2 = jpjglo
	259	DO jj = mj0(ind1), mj1(ind2)
	260	DO ji = 1, jpi
	261	ztabramp(ji,jj) = zmsknorth(ji)
	262	END DO
	263	END DO
	264
[4153]	265	ENDIF
	266
[3680]	267	! Tracers
	268	IF( .NOT. spongedoneT ) THEN
[12229]	269	fspu(:,:) = 0._wp
	270	fspv(:,:) = 0._wp
[12340]	271	DO_2D_00_00
	272	fspu(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji+1,jj ) ) * ssumask(ji,jj)
	273	fspv(ji,jj) = 0.5_wp * ( ztabramp(ji,jj) + ztabramp(ji ,jj+1) ) * ssvmask(ji,jj)
	274	END_2D
[12229]	275	CALL lbc_lnk( 'agrif_Sponge', fspu, 'U', 1. ) ! Lateral boundary conditions
	276	CALL lbc_lnk( 'agrif_Sponge', fspv, 'V', 1. )
	277
[3680]	278	spongedoneT = .TRUE.
	279	ENDIF
	280
	281	! Dynamics
	282	IF( .NOT. spongedoneU ) THEN
[12229]	283	fspt(:,:) = 0._wp
	284	fspf(:,:) = 0._wp
[12340]	285	DO_2D_00_00
	286	fspt(ji,jj) = ztabramp(ji,jj) * ssmask(ji,jj)
	287	fspf(ji,jj) = 0.25_wp * ( ztabramp(ji ,jj ) + ztabramp(ji ,jj+1) &
	288	& +ztabramp(ji+1,jj+1) + ztabramp(ji+1,jj ) ) &
	289	& * ssvmask(ji,jj) * ssvmask(ji,jj+1)
	290	END_2D
[12229]	291	CALL lbc_lnk( 'agrif_Sponge', fspt, 'T', 1. ) ! Lateral boundary conditions
	292	CALL lbc_lnk( 'agrif_Sponge', fspf, 'F', 1. )
[9019]	293
[3680]	294	spongedoneU = .TRUE.
	295	ENDIF
[12229]	296
	297	#if defined key_vertical
	298	! Remove vertical interpolation where not needed:
[12340]	299	DO_2D_00_00
	300	IF ((fspu(ji-1,jj)==0._wp).AND.(fspu(ji,jj)==0._wp).AND. &
	301	& (fspv(ji,jj-1)==0._wp).AND.(fspv(ji,jj)==0._wp)) mbkt_parent(ji,jj) = 0
[12229]	302	!
[12340]	303	IF ((fspt(ji+1,jj)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. &
	304	& (fspf(ji,jj-1)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbku_parent(ji,jj) = 0
[12229]	305	!
[12340]	306	IF ((fspt(ji,jj+1)==0._wp).AND.(fspt(ji,jj)==0._wp).AND. &
	307	& (fspf(ji-1,jj)==0._wp).AND.(fspf(ji,jj)==0._wp)) mbkv_parent(ji,jj) = 0
[12229]	308	!
[12340]	309	IF ( ssmask(ji,jj) == 0._wp) mbkt_parent(ji,jj) = 0
	310	IF (ssumask(ji,jj) == 0._wp) mbku_parent(ji,jj) = 0
	311	IF (ssvmask(ji,jj) == 0._wp) mbkv_parent(ji,jj) = 0
	312	END_2D
[3680]	313	!
[12229]	314	ztabramp(:,:) = REAL( mbkt_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'T', 1. )
	315	mbkt_parent(:,:) = NINT( ztabramp(:,:) )
	316	ztabramp(:,:) = REAL( mbku_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'U', 1. )
	317	mbku_parent(:,:) = NINT( ztabramp(:,:) )
	318	ztabramp(:,:) = REAL( mbkv_parent(:,:), wp ) ; CALL lbc_lnk( 'Agrif_Sponge', ztabramp, 'V', 1. )
	319	mbkv_parent(:,:) = NINT( ztabramp(:,:) )
[3680]	320	#endif
[6140]	321	!
[12229]	322	#endif
	323	!
[3680]	324	END SUBROUTINE Agrif_Sponge
	325
[12229]	326	SUBROUTINE interptsn_sponge( tabres, i1, i2, j1, j2, k1, k2, n1, n2, before )
[9019]	327	!!----------------------------------------------------------------------
	328	!! * ROUTINE interptsn_sponge *
	329	!!----------------------------------------------------------------------
	330	INTEGER , INTENT(in ) :: i1, i2, j1, j2, k1, k2, n1, n2
	331	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,n1:n2), INTENT(inout) :: tabres
	332	LOGICAL , INTENT(in ) :: before
[6140]	333	!
[5656]	334	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	335	INTEGER :: iku, ikv
[12229]	336	REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot, ztrelax
[9031]	337	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv
	338	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
	339	! vertical interpolation:
	340	REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
	341	REAL(wp), DIMENSION(k1:k2,n1:n2-1) :: tabin
	342	REAL(wp), DIMENSION(k1:k2) :: h_in
	343	REAL(wp), DIMENSION(1:jpk) :: h_out
	344	INTEGER :: N_in, N_out
[9019]	345	!!----------------------------------------------------------------------
[5656]	346	!
[6140]	347	IF( before ) THEN
[9031]	348	DO jn = 1, jpts
	349	DO jk=k1,k2
	350	DO jj=j1,j2
	351	DO ji=i1,i2
[11949]	352	tabres(ji,jj,jk,jn) = ts(ji,jj,jk,jn,Kbb_a)
[9031]	353	END DO
	354	END DO
	355	END DO
	356	END DO
	357
	358	# if defined key_vertical
[12229]	359	! Interpolate thicknesses
	360	! Warning: these are masked, hence extrapolated prior interpolation.
	361	DO jk=k1,k2
	362	DO jj=j1,j2
	363	DO ji=i1,i2
	364	tabres(ji,jj,jk,jpts+1) = tmask(ji,jj,jk) * e3t(ji,jj,jk,Kbb_a)
	365	END DO
	366	END DO
	367	END DO
	368
	369	! Extrapolate thicknesses in partial bottom cells:
	370	! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
	371	IF (ln_zps) THEN
	372	DO jj=j1,j2
	373	DO ji=i1,i2
	374	jk = mbkt(ji,jj)
	375	tabres(ji,jj,jk,jpts+1) = 0._wp
	376	END DO
	377	END DO
	378	END IF
	379
	380	! Save ssh at last level:
	381	IF (.NOT.ln_linssh) THEN
	382	tabres(i1:i2,j1:j2,k2,jpts+1) = ssh(i1:i2,j1:j2,Kbb_a)*tmask(i1:i2,j1:j2,1)
	383	ELSE
	384	tabres(i1:i2,j1:j2,k2,jpts+1) = 0._wp
	385	END IF
[9031]	386	# endif
	387
[5656]	388	ELSE
[6140]	389	!
[9031]	390	# if defined key_vertical
[12229]	391
	392	IF (ln_linssh) tabres(i1:i2,j1:j2,k2,n2) = 0._wp
	393
[9031]	394	DO jj=j1,j2
	395	DO ji=i1,i2
[12229]	396	tabres_child(ji,jj,:,:) = 0._wp
	397	N_in = mbkt_parent(ji,jj)
	398	zhtot = 0._wp
	399	DO jk=1,N_in !k2 = jpk of parent grid
	400	IF (jk==N_in) THEN
	401	h_in(jk) = ht0_parent(ji,jj) + tabres(ji,jj,k2,n2) - zhtot
	402	ELSE
	403	h_in(jk) = tabres(ji,jj,jk,n2)
	404	ENDIF
	405	zhtot = zhtot + h_in(jk)
[9031]	406	tabin(jk,:) = tabres(ji,jj,jk,n1:n2-1)
	407	END DO
	408	N_out = 0
	409	DO jk=1,jpk ! jpk of child grid
	410	IF (tmask(ji,jj,jk) == 0) EXIT
	411	N_out = N_out + 1
[12229]	412	h_out(jk) = e3t(ji,jj,jk,Kbb_a) !Child grid scale factors. Could multiply by e1e2t here instead of division above
[9031]	413	ENDDO
[12229]	414
	415	! Account for small differences in free-surface
	416	IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
	417	h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
	418	ELSE
	419	h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
[9031]	420	ENDIF
[12229]	421	IF (N_in*N_out > 0) THEN
	422	CALL reconstructandremap(tabin(1:N_in,1:jpts),h_in(1:N_in),tabres_child(ji,jj,1:N_out,1:jpts),h_out(1:N_out),N_in,N_out,jpts)
	423	ENDIF
[9031]	424	ENDDO
	425	ENDDO
	426	# endif
	427
	428	DO jj=j1,j2
	429	DO ji=i1,i2
	430	DO jk=1,jpkm1
	431	# if defined key_vertical
[12229]	432	tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres_child(ji,jj,jk,1:jpts)) * tmask(ji,jj,jk)
[9031]	433	# else
[12229]	434	tsbdiff(ji,jj,jk,1:jpts) = (ts(ji,jj,jk,1:jpts,Kbb_a) - tabres(ji,jj,jk,1:jpts))*tmask(ji,jj,jk)
[9031]	435	# endif
	436	ENDDO
	437	ENDDO
	438	ENDDO
	439
[12229]	440	!* set relaxation time scale
	441	IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_tra / ( rdt )
	442	ELSE ; ztrelax = rn_trelax_tra / (2._wp * rdt )
	443	ENDIF
	444
[5656]	445	DO jn = 1, jpts
	446	DO jk = 1, jpkm1
[9806]	447	ztu(i1:i2,j1:j2,jk) = 0._wp
	448	DO jj = j1,j2
[5656]	449	DO ji = i1,i2-1
[12229]	450	zabe1 = rn_sponge_tra * fspu(ji,jj) * umask(ji,jj,jk) * e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
[9806]	451	ztu(ji,jj,jk) = zabe1 * ( tsbdiff(ji+1,jj ,jk,jn) - tsbdiff(ji,jj,jk,jn) )
	452	END DO
	453	END DO
	454	ztv(i1:i2,j1:j2,jk) = 0._wp
	455	DO ji = i1,i2
	456	DO jj = j1,j2-1
[12229]	457	zabe2 = rn_sponge_tra * fspv(ji,jj) * vmask(ji,jj,jk) * e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
[5656]	458	ztv(ji,jj,jk) = zabe2 * ( tsbdiff(ji ,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn) )
[6140]	459	END DO
	460	END DO
	461	!
[5656]	462	IF( ln_zps ) THEN ! set gradient at partial step level
[9806]	463	DO jj = j1,j2
	464	DO ji = i1,i2
[5656]	465	! last level
	466	iku = mbku(ji,jj)
	467	ikv = mbkv(ji,jj)
[6140]	468	IF( iku == jk ) ztu(ji,jj,jk) = 0._wp
	469	IF( ikv == jk ) ztv(ji,jj,jk) = 0._wp
[5656]	470	END DO
	471	END DO
	472	ENDIF
[6140]	473	END DO
	474	!
[5656]	475	DO jk = 1, jpkm1
	476	DO jj = j1+1,j2-1
	477	DO ji = i1+1,i2-1
	478	IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
[11949]	479	zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
[5656]	480	! horizontal diffusive trends
[12229]	481	ztsa = zbtr * ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) + ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) &
	482	& - ztrelax * fspt(ji,jj) * tsbdiff(ji,jj,jk,jn)
[5656]	483	! add it to the general tracer trends
[11949]	484	ts(ji,jj,jk,jn,Krhs_a) = ts(ji,jj,jk,jn,Krhs_a) + ztsa
[5656]	485	ENDIF
[6140]	486	END DO
	487	END DO
	488	END DO
	489	!
	490	END DO
	491	!
[5656]	492	tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
[6140]	493	!
[5656]	494	ENDIF
[6140]	495	!
[5656]	496	END SUBROUTINE interptsn_sponge
	497
[9031]	498	SUBROUTINE interpun_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before)
	499	!!---------------------------------------------
	500	!! * ROUTINE interpun_sponge *
	501	!!---------------------------------------------
	502	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
	503	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
	504	LOGICAL, INTENT(in) :: before
[6140]	505
[9031]	506	INTEGER :: ji,jj,jk,jmax
	507
[5656]	508	! sponge parameters
[12229]	509	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax
[9031]	510	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: ubdiff
	511	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
	512	! vertical interpolation:
	513	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
	514	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
	515	REAL(wp), DIMENSION(1:jpk) :: h_out
[12229]	516	INTEGER ::N_in, N_out
[9031]	517	!!---------------------------------------------
[5656]	518	!
[6140]	519	IF( before ) THEN
[12229]	520	DO jk=k1,k2
[9031]	521	DO jj=j1,j2
	522	DO ji=i1,i2
[11949]	523	tabres(ji,jj,jk,m1) = uu(ji,jj,jk,Kbb_a)
[9031]	524	# if defined key_vertical
[12229]	525	tabres(ji,jj,jk,m2) = e3u(ji,jj,jk,Kbb_a)*umask(ji,jj,jk)
[9031]	526	# endif
	527	END DO
	528	END DO
	529	END DO
	530
[12229]	531	# if defined key_vertical
	532	! Extrapolate thicknesses in partial bottom cells:
	533	! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
	534	IF (ln_zps) THEN
	535	DO jj=j1,j2
	536	DO ji=i1,i2
	537	jk = mbku(ji,jj)
	538	tabres(ji,jj,jk,m2) = 0._wp
	539	END DO
	540	END DO
	541	END IF
	542	! Save ssh at last level:
	543	tabres(i1:i2,j1:j2,k2,m2) = 0._wp
	544	IF (.NOT.ln_linssh) THEN
	545	! This vertical sum below should be replaced by the sea-level at U-points (optimization):
	546	DO jk=1,jpk
	547	tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3u(i1:i2,j1:j2,jk,Kbb_a) * umask(i1:i2,j1:j2,jk)
	548	END DO
	549	tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hu_0(i1:i2,j1:j2)
	550	END IF
	551	# endif
	552
[5656]	553	ELSE
[9031]	554
	555	# if defined key_vertical
[12229]	556	IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
	557
[9031]	558	DO jj=j1,j2
	559	DO ji=i1,i2
[12229]	560	tabres_child(ji,jj,:) = 0._wp
	561	N_in = mbku_parent(ji,jj)
	562	zhtot = 0._wp
	563	DO jk=1,N_in
	564	IF (jk==N_in) THEN
	565	h_in(jk) = hu0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
	566	ELSE
	567	h_in(jk) = tabres(ji,jj,jk,m2)
	568	ENDIF
	569	zhtot = zhtot + h_in(jk)
[9031]	570	tabin(jk) = tabres(ji,jj,jk,m1)
[12229]	571	ENDDO
	572	!
	573	N_out = 0
	574	DO jk=1,jpk
	575	IF (umask(ji,jj,jk) == 0) EXIT
	576	N_out = N_out + 1
	577	h_out(N_out) = e3u(ji,jj,jk,Kbb_a)
	578	ENDDO
	579
	580	! Account for small differences in free-surface
	581	IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
	582	h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
	583	ELSE
	584	h_in(1) = h_in(1) - (sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
	585	ENDIF
	586
	587	IF (N_in * N_out > 0) THEN
	588	CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
	589	ENDIF
[9031]	590	ENDDO
	591	ENDDO
	592
[12229]	593	ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*umask(i1:i2,j1:j2,:)
[9031]	594	#else
[12229]	595	ubdiff(i1:i2,j1:j2,:) = (uu(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*umask(i1:i2,j1:j2,:)
[9031]	596	#endif
[12229]	597	!* set relaxation time scale
	598	IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_dyn / ( rdt )
	599	ELSE ; ztrelax = rn_trelax_dyn / (2._wp * rdt )
	600	ENDIF
[6140]	601	!
[5656]	602	DO jk = 1, jpkm1 ! Horizontal slab
	603	! ! ===============
	604
	605	! ! --------
	606	! Horizontal divergence ! div
	607	! ! --------
	608	DO jj = j1,j2
	609	DO ji = i1+1,i2 ! vector opt.
[12229]	610	zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj)
	611	hdivdiff(ji,jj,jk) = ( e2u(ji ,jj)e3u(ji ,jj,jk,Kbb_a) ubdiff(ji ,jj,jk) &
	612	& -e2u(ji-1,jj)e3u(ji-1,jj,jk,Kbb_a) ubdiff(ji-1,jj,jk) ) * zbtr
[5656]	613	END DO
	614	END DO
	615
	616	DO jj = j1,j2-1
	617	DO ji = i1,i2 ! vector opt.
[12229]	618	zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj)
[9019]	619	rotdiff(ji,jj,jk) = ( -e1u(ji,jj+1) * ubdiff(ji,jj+1,jk) &
	620	& +e1u(ji,jj ) * ubdiff(ji,jj ,jk) ) * fmask(ji,jj,jk) * zbtr
[5656]	621	END DO
	622	END DO
[6140]	623	END DO
[5656]	624	!
	625	DO jj = j1+1, j2-1
	626	DO ji = i1+1, i2-1 ! vector opt.
	627
	628	IF (.NOT. tabspongedone_u(ji,jj)) THEN
	629	DO jk = 1, jpkm1 ! Horizontal slab
	630	ze2u = rotdiff (ji,jj,jk)
	631	ze1v = hdivdiff(ji,jj,jk)
	632	! horizontal diffusive trends
[11949]	633	zua = - ( ze2u - rotdiff (ji,jj-1,jk) ) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) &
[12229]	634	& + ( hdivdiff(ji+1,jj,jk) - ze1v ) * r1_e1u(ji,jj) &
	635	& - ztrelax * fspu(ji,jj) * ubdiff(ji,jj,jk)
[5656]	636
	637	! add it to the general momentum trends
[12229]	638	uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) + zua
[5656]	639	END DO
	640	ENDIF
	641
	642	END DO
	643	END DO
	644
	645	tabspongedone_u(i1+1:i2-1,j1+1:j2-1) = .TRUE.
	646
	647	jmax = j2-1
[9019]	648	IF ((nbondj == 1).OR.(nbondj == 2)) jmax = MIN(jmax,nlcj-nbghostcells-2) ! North
[5656]	649
	650	DO jj = j1+1, jmax
	651	DO ji = i1+1, i2 ! vector opt.
	652
	653	IF (.NOT. tabspongedone_v(ji,jj)) THEN
	654	DO jk = 1, jpkm1 ! Horizontal slab
	655	ze2u = rotdiff (ji,jj,jk)
	656	ze1v = hdivdiff(ji,jj,jk)
	657
	658	! horizontal diffusive trends
[11949]	659	zva = + ( ze2u - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) ) &
[9019]	660	+ ( hdivdiff(ji,jj+1,jk) - ze1v ) * r1_e2v(ji,jj)
[5656]	661
	662	! add it to the general momentum trends
[11949]	663	vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) + zva
[5656]	664	END DO
	665	ENDIF
[6140]	666	!
[5656]	667	END DO
	668	END DO
[6140]	669	!
[5656]	670	tabspongedone_v(i1+1:i2,j1+1:jmax) = .TRUE.
[6140]	671	!
[5656]	672	ENDIF
[6140]	673	!
[5656]	674	END SUBROUTINE interpun_sponge
	675
[9031]	676	SUBROUTINE interpvn_sponge(tabres,i1,i2,j1,j2,k1,k2,m1,m2, before,nb,ndir)
	677	!!---------------------------------------------
	678	!! * ROUTINE interpvn_sponge *
	679	!!---------------------------------------------
	680	INTEGER, INTENT(in) :: i1,i2,j1,j2,k1,k2,m1,m2
	681	REAL(wp), DIMENSION(i1:i2,j1:j2,k1:k2,m1:m2), INTENT(inout) :: tabres
	682	LOGICAL, INTENT(in) :: before
	683	INTEGER, INTENT(in) :: nb , ndir
[6140]	684	!
[9031]	685	INTEGER :: ji, jj, jk, imax
[12229]	686	REAL(wp) :: ze2u, ze1v, zua, zva, zbtr, zhtot, ztrelax
[9031]	687	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: vbdiff
	688	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: rotdiff, hdivdiff
	689	! vertical interpolation:
	690	REAL(wp), DIMENSION(i1:i2,j1:j2,1:jpk) :: tabres_child
	691	REAL(wp), DIMENSION(k1:k2) :: tabin, h_in
	692	REAL(wp), DIMENSION(1:jpk) :: h_out
	693	INTEGER :: N_in, N_out
	694	!!---------------------------------------------
	695
[6140]	696	IF( before ) THEN
[12229]	697	DO jk=k1,k2
[9031]	698	DO jj=j1,j2
	699	DO ji=i1,i2
[11949]	700	tabres(ji,jj,jk,m1) = vv(ji,jj,jk,Kbb_a)
[9031]	701	# if defined key_vertical
[12229]	702	tabres(ji,jj,jk,m2) = vmask(ji,jj,jk) * e3v(ji,jj,jk,Kbb_a)
[9031]	703	# endif
	704	END DO
	705	END DO
	706	END DO
[12229]	707
	708	# if defined key_vertical
	709	! Extrapolate thicknesses in partial bottom cells:
	710	! Set them to Agrif_SpecialValue (0.). Correct bottom thicknesses are retrieved later on
	711	IF (ln_zps) THEN
	712	DO jj=j1,j2
	713	DO ji=i1,i2
	714	jk = mbkv(ji,jj)
	715	tabres(ji,jj,jk,m2) = 0._wp
	716	END DO
	717	END DO
	718	END IF
	719	! Save ssh at last level:
	720	tabres(i1:i2,j1:j2,k2,m2) = 0._wp
	721	IF (.NOT.ln_linssh) THEN
	722	! This vertical sum below should be replaced by the sea-level at V-points (optimization):
	723	DO jk=1,jpk
	724	tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) + e3v(i1:i2,j1:j2,jk,Kbb_a) * vmask(i1:i2,j1:j2,jk)
	725	END DO
	726	tabres(i1:i2,j1:j2,k2,m2) = tabres(i1:i2,j1:j2,k2,m2) - hv_0(i1:i2,j1:j2)
	727	END IF
	728	# endif
	729
[5656]	730	ELSE
[9031]	731
	732	# if defined key_vertical
[12229]	733	IF (ln_linssh) tabres(i1:i2,j1:j2,k2,m2) = 0._wp
[9031]	734	DO jj=j1,j2
	735	DO ji=i1,i2
[12229]	736	tabres_child(ji,jj,:) = 0._wp
	737	N_in = mbkv_parent(ji,jj)
	738	zhtot = 0._wp
	739	DO jk=1,N_in
	740	IF (jk==N_in) THEN
	741	h_in(jk) = hv0_parent(ji,jj) + tabres(ji,jj,k2,m2) - zhtot
	742	ELSE
	743	h_in(jk) = tabres(ji,jj,jk,m2)
	744	ENDIF
	745	zhtot = zhtot + h_in(jk)
[9031]	746	tabin(jk) = tabres(ji,jj,jk,m1)
[12229]	747	ENDDO
	748	!
	749	N_out = 0
	750	DO jk=1,jpk
	751	IF (vmask(ji,jj,jk) == 0) EXIT
	752	N_out = N_out + 1
	753	h_out(N_out) = e3v(ji,jj,jk,Kbb_a)
	754	ENDDO
	755
	756	! Account for small differences in free-surface
	757	IF ( sum(h_out(1:N_out)) > sum(h_in(1:N_in) )) THEN
	758	h_out(1) = h_out(1) - ( sum(h_out(1:N_out))-sum(h_in(1:N_in)) )
	759	ELSE
	760	h_in(1) = h_in(1) - ( sum(h_in(1:N_in))-sum(h_out(1:N_out)) )
	761	ENDIF
[9031]	762
[12229]	763	IF (N_in * N_out > 0) THEN
	764	CALL reconstructandremap(tabin(1:N_in),h_in(1:N_in),tabres_child(ji,jj,1:N_out),h_out(1:N_out),N_in,N_out,1)
	765	ENDIF
[9031]	766	ENDDO
	767	ENDDO
	768
[12229]	769	vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres_child(i1:i2,j1:j2,:))*vmask(i1:i2,j1:j2,:)
[9031]	770	# else
[12229]	771	vbdiff(i1:i2,j1:j2,:) = (vv(i1:i2,j1:j2,:,Kbb_a) - tabres(i1:i2,j1:j2,:,1))*vmask(i1:i2,j1:j2,:)
[9031]	772	# endif
[12229]	773	!* set relaxation time scale
	774	IF( neuler == 0 .AND. lk_agrif_fstep ) THEN ; ztrelax = rn_trelax_dyn / ( rdt )
	775	ELSE ; ztrelax = rn_trelax_dyn / (2._wp * rdt )
	776	ENDIF
[6140]	777	!
[5656]	778	DO jk = 1, jpkm1 ! Horizontal slab
	779	! ! ===============
	780
	781	! ! --------
	782	! Horizontal divergence ! div
	783	! ! --------
	784	DO jj = j1+1,j2
	785	DO ji = i1,i2 ! vector opt.
[12229]	786	zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kbb_a) * rn_sponge_dyn * fspt(ji,jj)
	787	hdivdiff(ji,jj,jk) = ( e1v(ji,jj ) * e3v(ji,jj ,jk,Kbb_a) * vbdiff(ji,jj ,jk) &
	788	& -e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kbb_a) * vbdiff(ji,jj-1,jk) ) * zbtr
[5656]	789	END DO
	790	END DO
	791	DO jj = j1,j2
	792	DO ji = i1,i2-1 ! vector opt.
[12229]	793	zbtr = r1_e1e2f(ji,jj) * e3f(ji,jj,jk) * rn_sponge_dyn * fspf(ji,jj)
	794	rotdiff(ji,jj,jk) = ( e2v(ji+1,jj) * vbdiff(ji+1,jj,jk) &
	795	& -e2v(ji ,jj) * vbdiff(ji ,jj,jk) ) * fmask(ji,jj,jk) * zbtr
[5656]	796	END DO
	797	END DO
[6140]	798	END DO
[5656]	799
	800	! ! ===============
	801	!
	802
[9019]	803	imax = i2 - 1
	804	IF ((nbondi == 1).OR.(nbondi == 2)) imax = MIN(imax,nlci-nbghostcells-2) ! East
[5656]	805
	806	DO jj = j1+1, j2
	807	DO ji = i1+1, imax ! vector opt.
[6140]	808	IF( .NOT. tabspongedone_u(ji,jj) ) THEN
	809	DO jk = 1, jpkm1
[12229]	810	uu(ji,jj,jk,Krhs_a) = uu(ji,jj,jk,Krhs_a) &
	811	& - ( rotdiff (ji ,jj,jk) - rotdiff (ji,jj-1,jk)) / ( e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a) ) &
	812	& + ( hdivdiff(ji+1,jj,jk) - hdivdiff(ji,jj ,jk)) * r1_e1u(ji,jj)
[5656]	813	END DO
	814	ENDIF
	815	END DO
	816	END DO
[6140]	817	!
[5656]	818	tabspongedone_u(i1+1:imax,j1+1:j2) = .TRUE.
[6140]	819	!
[5656]	820	DO jj = j1+1, j2-1
	821	DO ji = i1+1, i2-1 ! vector opt.
[6140]	822	IF( .NOT. tabspongedone_v(ji,jj) ) THEN
	823	DO jk = 1, jpkm1
[12229]	824	vv(ji,jj,jk,Krhs_a) = vv(ji,jj,jk,Krhs_a) &
	825	& + ( rotdiff (ji,jj ,jk) - rotdiff (ji-1,jj,jk) ) / ( e1v(ji,jj) * e3v(ji,jj,jk,Kmm_a) ) &
	826	& + ( hdivdiff(ji,jj+1,jk) - hdivdiff(ji ,jj,jk) ) * r1_e2v(ji,jj) &
	827	& - ztrelax * fspv(ji,jj) * vbdiff(ji,jj,jk)
[5656]	828	END DO
	829	ENDIF
	830	END DO
	831	END DO
	832	tabspongedone_v(i1+1:i2-1,j1+1:j2-1) = .TRUE.
	833	ENDIF
[6140]	834	!
[5656]	835	END SUBROUTINE interpvn_sponge
	836
[390]	837	#else
[9019]	838	!!----------------------------------------------------------------------
	839	!! Empty module no AGRIF zoom
	840	!!----------------------------------------------------------------------
[636]	841	CONTAINS
[9570]	842	SUBROUTINE agrif_oce_sponge_empty
	843	WRITE(,) 'agrif_oce_sponge : You should not have seen this print! error?'
	844	END SUBROUTINE agrif_oce_sponge_empty
[390]	845	#endif
[636]	846
[6140]	847	!!======================================================================
[9570]	848	END MODULE agrif_oce_sponge

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

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

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/NST/agrif_oce_sponge.F90 @ 12340

Download in other formats: