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 @ 258

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

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

Download in other formats: