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.

traadv_tvd.F90 in trunk/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: trunk/NEMO/OPA_SRC/TRA/traadv_tvd.F90 @ 247

Last change on this file since 247 was 247, checked in by opalod, 20 years ago
CL : Add CVS Header and CeCILL licence information
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 22.1 KB

Rev	Line
[3]	1	MODULE traadv_tvd
	2	!!==============================================================================
	3	!! * MODULE traadv_tvd *
	4	!! Ocean active tracers: horizontal & vertical advective trend
	5	!!==============================================================================
	6
	7	!!----------------------------------------------------------------------
	8	!! tra_adv_tvd : update the tracer trend with the horizontal
	9	!! and vertical advection trends using a TVD scheme
	10	!! nonosc : compute monotonic tracer fluxes by a nonoscillatory
	11	!! algorithm
	12	!!----------------------------------------------------------------------
	13	!! * Modules used
	14	USE oce ! ocean dynamics and active tracers
	15	USE dom_oce ! ocean space and time domain
[216]	16	USE trdmod ! ocean active tracers trends
	17	USE trdmod_oce ! ocean variables trends
[3]	18	USE in_out_manager ! I/O manager
	19	USE dynspg_fsc ! surface pressure gradient
	20	USE dynspg_fsc_atsk ! autotasked surface pressure gradient
	21	USE trabbl ! Advective term of BBL
	22	USE lib_mpp
[74]	23	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[132]	24	USE diaptr ! poleward transport diagnostics
[3]	25
[74]	26
[3]	27	IMPLICIT NONE
	28	PRIVATE
	29
	30	!! * Accessibility
	31	PUBLIC tra_adv_tvd ! routine called by step.F90
	32
	33	!! * Substitutions
	34	# include "domzgr_substitute.h90"
	35	# include "vectopt_loop_substitute.h90"
	36	!!----------------------------------------------------------------------
[247]	37	!! OPA 9.0 , LOCEAN-IPSL (2005)
	38	!! $Header$
	39	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
[3]	40	!!----------------------------------------------------------------------
	41
	42	CONTAINS
	43
	44	SUBROUTINE tra_adv_tvd( kt )
	45	!!----------------------------------------------------------------------
	46	!! * ROUTINE tra_adv_tvd *
	47	!!
	48	!! ** Purpose : Compute the now trend due to total advection of
	49	!! tracers and add it to the general trend of tracer equations
	50	!!
	51	!! ** Method : TVD scheme, i.e. 2nd order centered scheme with
	52	!! corrected flux (monotonic correction)
	53	!! note: - this advection scheme needs a leap-frog time scheme
	54	!!
	55	!! ** Action : - update (ta,sa) with the now advective tracer trends
	56	!! - save the trends in (ttrdh,strdh) ('key_trdtra')
	57	!!
	58	!! History :
	59	!! ! 95-12 (L. Mortier) Original code
	60	!! ! 00-01 (H. Loukos) adapted to ORCA
	61	!! ! 00-10 (MA Foujols E.Kestenare) include file not routine
	62	!! ! 00-12 (E. Kestenare M. Levy) fix bug in trtrd indexes
	63	!! ! 01-07 (E. Durand G. Madec) adaptation to ORCA config
	64	!! 8.5 ! 02-06 (G. Madec) F90: Free form and module
	65	!! 9.0 ! 04-01 (A. de Miranda, G. Madec, J.M. Molines ): advective bbl
[216]	66	!! 9.0 ! 08-04 (S. Cravatte) add the i-, j- & k- trends computation
[3]	67	!!----------------------------------------------------------------------
	68	!! * Modules used
	69	#if defined key_trabbl_adv
	70	USE oce , zun => ua, & ! use ua as workspace
	71	& zvn => va ! use va as workspace
	72	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwn
	73	#else
	74	USE oce , zun => un, & ! When no bbl, zun == un
	75	zvn => vn, & ! zvn == vn
	76	zwn => wn ! zwn == wn
	77	#endif
[216]	78	USE trdmod_oce , ztay => tladj, & ! use tladj latter
	79	& zsay => sladj, & ! use sladj latter
	80	& ztaz => tladi, & ! use ua as workspace
	81	& zsaz => sladi ! use ua as workspace
	82
[3]	83	!! * Arguments
	84	INTEGER, INTENT( in ) :: kt ! ocean time-step
	85
	86	!! * Local declarations
	87	INTEGER :: ji, jj, jk ! dummy loop indices
[216]	88	REAL(wp) :: zta, zsa, & ! temporary scalar
	89	ztai, ztaj, ztak, & ! " "
	90	zsai, zsaj, zsak ! " "
[3]	91	REAL(wp), DIMENSION (jpi,jpj,jpk) :: &
	92	zti, ztu, ztv, ztw, & ! temporary workspace
[216]	93	zsi, zsu, zsv, zsw, & ! " "
	94	ztdta, ztdsa ! " "
[3]	95	REAL(wp) :: &
	96	z2dtt, zbtr, zeu, zev, zew, z2, & ! temporary scalar
	97	zfp_ui, zfp_vj, zfp_wk, & ! " "
	98	zfm_ui, zfm_vj, zfm_wk ! " "
	99	!!----------------------------------------------------------------------
	100
	101	IF( kt == nit000 .AND. lwp ) THEN
	102	WRITE(numout,*)
	103	WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme'
	104	WRITE(numout,*) '~~~~~~~~~~~'
	105	ENDIF
	106
	107	IF( neuler == 0 .AND. kt == nit000 ) THEN
	108	z2=1.
	109	ELSE
	110	z2=2.
	111	ENDIF
	112
[216]	113	! Save ta and sa trends
	114	IF( l_trdtra ) THEN
	115	ztdta(:,:,:) = ta(:,:,:)
	116	ztdsa(:,:,:) = sa(:,:,:)
	117	l_adv = 'tvd'
	118	ENDIF
	119
[3]	120	#if defined key_trabbl_adv
	121	! Advective Bottom boundary layer: add the velocity
	122	! -------------------------------------------------
	123	zun(:,:,:) = un (:,:,:) - u_bbl(:,:,:)
	124	zvn(:,:,:) = vn (:,:,:) - v_bbl(:,:,:)
	125	zwn(:,:,:) = wn (:,:,:) + w_bbl(:,:,:)
	126	#endif
	127
	128	! 1. Bottom value : flux set to zero
	129	! ---------------
	130	ztu(:,:,jpk) = 0.e0 ; zsu(:,:,jpk) = 0.e0
	131	ztv(:,:,jpk) = 0.e0 ; zsv(:,:,jpk) = 0.e0
	132	ztw(:,:,jpk) = 0.e0 ; zsw(:,:,jpk) = 0.e0
	133	zti(:,:,jpk) = 0.e0 ; zsi(:,:,jpk) = 0.e0
	134
	135
	136	! 2. upstream advection with initial mass fluxes & intermediate update
	137	! --------------------------------------------------------------------
	138	! upstream tracer flux in the i and j direction
	139	DO jk = 1, jpkm1
	140	DO jj = 1, jpjm1
	141	DO ji = 1, fs_jpim1 ! vector opt.
	142	zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
	143	zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
	144	! upstream scheme
	145	zfp_ui = zeu + ABS( zeu )
	146	zfm_ui = zeu - ABS( zeu )
	147	zfp_vj = zev + ABS( zev )
	148	zfm_vj = zev - ABS( zev )
	149	ztu(ji,jj,jk) = zfp_ui * tb(ji,jj,jk) + zfm_ui * tb(ji+1,jj ,jk)
	150	ztv(ji,jj,jk) = zfp_vj * tb(ji,jj,jk) + zfm_vj * tb(ji ,jj+1,jk)
	151	zsu(ji,jj,jk) = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj ,jk)
	152	zsv(ji,jj,jk) = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji ,jj+1,jk)
	153	END DO
	154	END DO
	155	END DO
	156
	157	! upstream tracer flux in the k direction
	158	! Surface value
	159	IF( lk_dynspg_fsc .OR. lk_dynspg_fsc_tsk ) THEN ! free surface-constant volume
	160	DO jj = 1, jpj
	161	DO ji = 1, jpi
	162	zew = e1t(ji,jj) * e2t(ji,jj) * zwn(ji,jj,1)
	163	ztw(ji,jj,1) = zew * tb(ji,jj,1)
	164	zsw(ji,jj,1) = zew * sb(ji,jj,1)
	165	END DO
	166	END DO
[216]	167	ELSE ! rigid lid : flux set to zero
[3]	168	ztw(:,:,1) = 0.e0
	169	zsw(:,:,1) = 0.e0
	170	ENDIF
	171
	172	! Interior value
	173	DO jk = 2, jpkm1
	174	DO jj = 1, jpj
	175	DO ji = 1, jpi
	176	zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * zwn(ji,jj,jk)
	177	zfp_wk = zew + ABS( zew )
	178	zfm_wk = zew - ABS( zew )
	179	ztw(ji,jj,jk) = zfp_wk * tb(ji,jj,jk) + zfm_wk * tb(ji,jj,jk-1)
	180	zsw(ji,jj,jk) = zfp_wk * sb(ji,jj,jk) + zfm_wk * sb(ji,jj,jk-1)
	181	END DO
	182	END DO
	183	END DO
	184
	185	! total advective trend
	186	DO jk = 1, jpkm1
	187	DO jj = 2, jpjm1
	188	DO ji = fs_2, fs_jpim1 ! vector opt.
	189	zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
[216]	190
	191	! i- j- horizontal & k- vertical advective trends
	192	ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk ) ) * zbtr
	193	ztaj = - ( ztv(ji,jj,jk) - ztv(ji ,jj-1,jk ) ) * zbtr
	194	ztak = - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) * zbtr
	195	zsai = - ( zsu(ji,jj,jk) - zsu(ji-1,jj ,jk ) ) * zbtr
	196	zsaj = - ( zsv(ji,jj,jk) - zsv(ji ,jj-1,jk ) ) * zbtr
	197	zsak = - ( zsw(ji,jj,jk) - zsw(ji ,jj ,jk+1) ) * zbtr
	198
	199	! total intermediate advective trends
	200	zti(ji,jj,jk) = ztai + ztaj + ztak
	201	zsi(ji,jj,jk) = zsai + zsaj + zsak
[3]	202	END DO
	203	END DO
	204	END DO
	205
[216]	206	! Save the intermediate vertical & j- horizontal advection trends
	207	IF( l_trdtra ) THEN
	208	DO jk = 1, jpkm1
	209	DO jj = 2, jpjm1
	210	DO ji = fs_2, fs_jpim1 ! vector opt.
	211	zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	212	ztay(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji ,jj-1,jk ) ) * zbtr
	213	zsay(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji ,jj-1,jk ) ) * zbtr
	214	ztaz(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) * zbtr
	215	zsaz(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji ,jj ,jk+1) ) * zbtr
	216	END DO
	217	END DO
	218	END DO
	219	ENDIF
	220
[3]	221	! update and guess with monotonic sheme
	222	DO jk = 1, jpkm1
	223	z2dtt = z2 * rdttra(jk)
	224	DO jj = 2, jpjm1
	225	DO ji = fs_2, fs_jpim1 ! vector opt.
	226	ta(ji,jj,jk) = ta(ji,jj,jk) + zti(ji,jj,jk)
	227	sa(ji,jj,jk) = sa(ji,jj,jk) + zsi(ji,jj,jk)
	228	zti (ji,jj,jk) = ( tb(ji,jj,jk) + z2dtt * zti(ji,jj,jk) ) * tmask(ji,jj,jk)
	229	zsi (ji,jj,jk) = ( sb(ji,jj,jk) + z2dtt * zsi(ji,jj,jk) ) * tmask(ji,jj,jk)
	230	END DO
	231	END DO
	232	END DO
	233
	234	! Lateral boundary conditions on zti, zsi (unchanged sign)
	235	CALL lbc_lnk( zti, 'T', 1. )
	236	CALL lbc_lnk( zsi, 'T', 1. )
	237
	238
	239	! 3. antidiffusive flux : high order minus low order
	240	! --------------------------------------------------
	241	! antidiffusive flux on i and j
	242	DO jk = 1, jpkm1
	243	DO jj = 1, jpjm1
	244	DO ji = 1, fs_jpim1 ! vector opt.
	245	zeu = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
	246	zev = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
	247	ztu(ji,jj,jk) = zeu * ( tn(ji,jj,jk) + tn(ji+1,jj,jk) ) - ztu(ji,jj,jk)
	248	zsu(ji,jj,jk) = zeu * ( sn(ji,jj,jk) + sn(ji+1,jj,jk) ) - zsu(ji,jj,jk)
	249	ztv(ji,jj,jk) = zev * ( tn(ji,jj,jk) + tn(ji,jj+1,jk) ) - ztv(ji,jj,jk)
	250	zsv(ji,jj,jk) = zev * ( sn(ji,jj,jk) + sn(ji,jj+1,jk) ) - zsv(ji,jj,jk)
	251	END DO
	252	END DO
	253	END DO
	254
	255	! antidiffusive flux on k
	256	! Surface value
	257	ztw(:,:,1) = 0.
	258	zsw(:,:,1) = 0.
	259
	260	! Interior value
	261	DO jk = 2, jpkm1
	262	DO jj = 1, jpj
	263	DO ji = 1, jpi
	264	zew = 0.5 * e1t(ji,jj) * e2t(ji,jj) * zwn(ji,jj,jk)
	265	ztw(ji,jj,jk) = zew * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) ) - ztw(ji,jj,jk)
	266	zsw(ji,jj,jk) = zew * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) ) - zsw(ji,jj,jk)
	267	END DO
	268	END DO
	269	END DO
	270
	271	! Lateral bondary conditions
	272	CALL lbc_lnk( ztu, 'U', -1. ) ; CALL lbc_lnk( zsu, 'U', -1. )
	273	CALL lbc_lnk( ztv, 'V', -1. ) ; CALL lbc_lnk( zsv, 'V', -1. )
	274	CALL lbc_lnk( ztw, 'W', 1. ) ; CALL lbc_lnk( zsw, 'W', 1. )
	275
	276	! 4. monotonicity algorithm
	277	! -------------------------
	278	CALL nonosc( tb, ztu, ztv, ztw, zti, z2 )
	279	CALL nonosc( sb, zsu, zsv, zsw, zsi, z2 )
	280
	281
	282	! 5. final trend with corrected fluxes
	283	! ------------------------------------
	284	DO jk = 1, jpkm1
	285	DO jj = 2, jpjm1
[216]	286	DO ji = fs_2, fs_jpim1 ! vector opt.
[3]	287	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
[216]	288	! i- j- horizontal & k- vertical advective trends
	289	ztai = - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk )) * zbtr
	290	ztaj = - ( ztv(ji,jj,jk) - ztv(ji ,jj-1,jk )) * zbtr
	291	ztak = - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1)) * zbtr
	292	zsai = - ( zsu(ji,jj,jk) - zsu(ji-1,jj ,jk )) * zbtr
	293	zsaj = - ( zsv(ji,jj,jk) - zsv(ji ,jj-1,jk )) * zbtr
	294	zsak = - ( zsw(ji,jj,jk) - zsw(ji ,jj ,jk+1)) * zbtr
	295
	296	! add them to the general tracer trends
	297	ta(ji,jj,jk) = ta(ji,jj,jk) + ztai + ztaj + ztak
	298	sa(ji,jj,jk) = sa(ji,jj,jk) + zsai + zsaj + zsak
[3]	299	END DO
	300	END DO
	301	END DO
	302
[216]	303	! save the advective trends for diagnostic
	304	! tracers trends
	305	IF( l_trdtra ) THEN
	306	! Compute the final vertical & j- horizontal advection trends
	307	DO jk = 1, jpkm1
	308	DO jj = 2, jpjm1
	309	DO ji = fs_2, fs_jpim1 ! vector opt.
	310	zbtr = 1./ ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	311	ztay(ji,jj,jk) = - ( ztv(ji,jj,jk) - ztv(ji ,jj-1,jk ) ) * zbtr &
	312	& + ztay(ji,jj,jk)
	313	zsay(ji,jj,jk) = - ( zsv(ji,jj,jk) - zsv(ji ,jj-1,jk ) ) * zbtr &
	314	& + zsay(ji,jj,jk)
	315	ztaz(ji,jj,jk) = - ( ztw(ji,jj,jk) - ztw(ji ,jj ,jk+1) ) * zbtr &
	316	& + ztaz(ji,jj,jk)
	317	zsaz(ji,jj,jk) = - ( zsw(ji,jj,jk) - zsw(ji ,jj ,jk+1) ) * zbtr &
	318	& + zsaz(ji,jj,jk)
	319	END DO
	320	END DO
	321	END DO
	322
	323	! horizontal advection:
	324	! make the difference between the new trends ta()/sa() and the
	325	! previous one ztdta()/ztdsa() to have the total advection trends
	326	! to which we substract the vertical trends ztaz()/zsaz()
	327	ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) - ztaz(:,:,:)
	328	ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) - zsaz(:,:,:)
	329
	330	! Add the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends
	331	ztdta(:,:,:) = ztdta(:,:,:) + tn(:,:,:) * hdivn(:,:,:)
	332	ztdsa(:,:,:) = ztdsa(:,:,:) + sn(:,:,:) * hdivn(:,:,:)
	333
	334	CALL trd_mod(ztdta, ztdsa, jpttdlad, 'TRA', kt)
	335
	336	! vertical advection:
	337	! Substract the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
	338	ztaz(:,:,:) = ztaz(:,:,:) - tn(:,:,:) * hdivn(:,:,:)
	339	zsaz(:,:,:) = zsaz(:,:,:) - sn(:,:,:) * hdivn(:,:,:)
	340
	341	CALL trd_mod(ztaz, zsaz, jpttdzad, 'TRA', kt)
	342
	343	ENDIF
	344
[106]	345	IF(l_ctl) THEN
	346	zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
	347	zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
[3]	348	WRITE(numout,*) ' zad - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl, ' tvd'
	349	t_ctl = zta ; s_ctl = zsa
	350	ENDIF
	351
[132]	352	! "zonal" mean advective heat and salt transport
	353	IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
	354	pht_adv(:) = ptr_vj( ztv(:,:,:) )
	355	pst_adv(:) = ptr_vj( zsv(:,:,:) )
[3]	356	ENDIF
	357
	358	END SUBROUTINE tra_adv_tvd
	359
	360
	361	SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, prdt )
	362	!!---------------------------------------------------------------------
	363	!! * ROUTINE nonosc *
	364	!!
	365	!! ** Purpose : compute monotonic tracer fluxes from the upstream
	366	!! scheme and the before field by a nonoscillatory algorithm
	367	!!
	368	!! ** Method : ... ???
	369	!! warning : pbef and paft must be masked, but the boundaries
	370	!! conditions on the fluxes are not necessary zalezak (1979)
	371	!! drange (1995) multi-dimensional forward-in-time and upstream-
	372	!! in-space based differencing for fluid
	373	!!
	374	!! History :
	375	!! ! 97-04 (L. Mortier) Original code
	376	!! ! 00-02 (H. Loukos) rewritting for opa8
	377	!! ! 00-10 (M.A Foujols, E. Kestenare) lateral b.c.
	378	!! ! 01-03 (E. Kestenare) add key_passivetrc
	379	!! ! 01-07 (E. Durand G. Madec) adapted for T & S
	380	!! 8.5 ! 02-06 (G. Madec) F90: Free form and module
	381	!!----------------------------------------------------------------------
	382	!! * Arguments
	383	REAL(wp), INTENT( in ) :: &
	384	prdt ! ???
	385	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT( inout ) :: &
	386	pbef, & ! before field
	387	paft, & ! after field
	388	paa, & ! monotonic flux in the i direction
	389	pbb, & ! monotonic flux in the j direction
	390	pcc ! monotonic flux in the k direction
	391
	392	!! * Local declarations
	393	INTEGER :: ji, jj, jk ! dummy loop indices
	394	INTEGER :: ikm1
	395	REAL(wp), DIMENSION (jpi,jpj,jpk) :: zbetup, zbetdo
	396	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt
	397	!!----------------------------------------------------------------------
	398
	399	zbig = 1.e+40
	400	zrtrn = 1.e-15
	401
	402	! Search local extrema
	403	! --------------------
	404	! large negative value (-zbig) inside land
[237]	405	pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
	406	paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
[3]	407	! search maximum in neighbourhood
	408	DO jk = 1, jpkm1
	409	ikm1 = MAX(jk-1,1)
	410	DO jj = 2, jpjm1
	411	DO ji = fs_2, fs_jpim1 ! vector opt.
	412	zbetup(ji,jj,jk) = MAX( pbef(ji ,jj ,jk ), paft(ji ,jj ,jk ), &
	413	& pbef(ji-1,jj ,jk ), pbef(ji+1,jj ,jk ), &
	414	& paft(ji-1,jj ,jk ), paft(ji+1,jj ,jk ), &
	415	& pbef(ji ,jj-1,jk ), pbef(ji ,jj+1,jk ), &
	416	& paft(ji ,jj-1,jk ), paft(ji ,jj+1,jk ), &
	417	& pbef(ji ,jj ,ikm1), pbef(ji ,jj ,jk+1), &
	418	& paft(ji ,jj ,ikm1), paft(ji ,jj ,jk+1) )
	419	END DO
	420	END DO
	421	END DO
	422	! large positive value (+zbig) inside land
[237]	423	pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
	424	paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
[3]	425	! search minimum in neighbourhood
	426	DO jk = 1, jpkm1
	427	ikm1 = MAX(jk-1,1)
	428	DO jj = 2, jpjm1
	429	DO ji = fs_2, fs_jpim1 ! vector opt.
	430	zbetdo(ji,jj,jk) = MIN( pbef(ji ,jj ,jk ), paft(ji ,jj ,jk ), &
	431	& pbef(ji-1,jj ,jk ), pbef(ji+1,jj ,jk ), &
	432	& paft(ji-1,jj ,jk ), paft(ji+1,jj ,jk ), &
	433	& pbef(ji ,jj-1,jk ), pbef(ji ,jj+1,jk ), &
	434	& paft(ji ,jj-1,jk ), paft(ji ,jj+1,jk ), &
	435	& pbef(ji ,jj ,ikm1), pbef(ji ,jj ,jk+1), &
	436	& paft(ji ,jj ,ikm1), paft(ji ,jj ,jk+1) )
	437	END DO
	438	END DO
	439	END DO
	440
	441	! restore masked values to zero
	442	pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
	443	paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
	444
	445
	446	! 2. Positive and negative part of fluxes and beta terms
	447	! ------------------------------------------------------
	448
	449	DO jk = 1, jpkm1
	450	z2dtt = prdt * rdttra(jk)
	451	DO jj = 2, jpjm1
	452	DO ji = fs_2, fs_jpim1 ! vector opt.
	453	! positive & negative part of the flux
	454	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
	455	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
	456	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
	457	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
	458	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
	459	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
	460	! up & down beta terms
	461	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
	462	zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
	463	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
	464	END DO
	465	END DO
	466	END DO
	467
	468	! lateral boundary condition on zbetup & zbetdo (unchanged sign)
	469	CALL lbc_lnk( zbetup, 'T', 1. )
	470	CALL lbc_lnk( zbetdo, 'T', 1. )
	471
	472
[237]	473	! 3. monotonic flux in the i & j direction (paa & pbb)
	474	! ----------------------------------------
[3]	475	DO jk = 1, jpkm1
	476	DO jj = 2, jpjm1
	477	DO ji = fs_2, fs_jpim1 ! vector opt.
[237]	478	za = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
	479	zb = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
	480	zc = 0.5 * ( 1.e0 + SIGN( 1.e0, paa(ji,jj,jk) ) )
	481	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
[3]	482
[237]	483	za = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
	484	zb = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
	485	zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pbb(ji,jj,jk) ) )
	486	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
[3]	487	END DO
	488	END DO
	489	END DO
	490
	491
	492	! monotonic flux in the k direction, i.e. pcc
	493	! -------------------------------------------
	494	DO jk = 2, jpkm1
	495	DO jj = 2, jpjm1
	496	DO ji = fs_2, fs_jpim1 ! vector opt.
[237]	497
	498	za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
	499	zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
	500	zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
	501	pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
[3]	502	END DO
	503	END DO
	504	END DO
	505
[237]	506	! lateral boundary condition on paa, pbb, pcc
	507	CALL lbc_lnk( paa, 'U', -1. ) ! changed sign
	508	CALL lbc_lnk( pbb, 'V', -1. ) ! changed sign
	509	CALL lbc_lnk( pcc, 'W', 1. ) ! NO changed sign
[3]	510
	511	END SUBROUTINE nonosc
	512
	513	!!======================================================================
	514	END MODULE traadv_tvd

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

Download in other formats: