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 branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_tvd.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 20.8 KB

Rev	Line
[3]	1	MODULE traadv_tvd
	2	!!==============================================================================
	3	!! * MODULE traadv_tvd *
[2528]	4	!! Ocean tracers: horizontal & vertical advective trend
[3]	5	!!==============================================================================
[2528]	6	!! History : OPA ! 1995-12 (L. Mortier) Original code
	7	!! ! 2000-01 (H. Loukos) adapted to ORCA
	8	!! ! 2000-10 (MA Foujols E.Kestenare) include file not routine
	9	!! ! 2000-12 (E. Kestenare M. Levy) fix bug in trtrd indexes
	10	!! ! 2001-07 (E. Durand G. Madec) adaptation to ORCA config
	11	!! 8.5 ! 2002-06 (G. Madec) F90: Free form and module
	12	!! NEMO 1.0 ! 2004-01 (A. de Miranda, G. Madec, J.M. Molines ): advective bbl
	13	!! 2.0 ! 2008-04 (S. Cravatte) add the i-, j- & k- trends computation
	14	!! - ! 2009-11 (V. Garnier) Surface pressure gradient organization
	15	!! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
[503]	16	!!----------------------------------------------------------------------
[3]	17
	18	!!----------------------------------------------------------------------
	19	!! tra_adv_tvd : update the tracer trend with the horizontal
	20	!! and vertical advection trends using a TVD scheme
	21	!! nonosc : compute monotonic tracer fluxes by a nonoscillatory
	22	!! algorithm
	23	!!----------------------------------------------------------------------
	24	USE oce ! ocean dynamics and active tracers
	25	USE dom_oce ! ocean space and time domain
[2528]	26	USE trdmod_oce ! tracers trends
[2715]	27	USE trdtra ! tracers trends
[3]	28	USE in_out_manager ! I/O manager
[367]	29	USE dynspg_oce ! choice/control of key cpp for surface pressure gradient
[2715]	30	USE lib_mpp ! MPP library
[74]	31	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[132]	32	USE diaptr ! poleward transport diagnostics
[2528]	33	USE trc_oce ! share passive tracers/Ocean variables
[3]	34
[74]	35
[3]	36	IMPLICIT NONE
	37	PRIVATE
	38
[503]	39	PUBLIC tra_adv_tvd ! routine called by step.F90
[3]	40
[2715]	41	LOGICAL :: l_trd ! flag to compute trends
[2528]	42
[3211]	43	!! * Control permutation of array indices
	44	# include "oce_ftrans.h90"
	45	# include "dom_oce_ftrans.h90"
	46	# include "trc_oce_ftrans.h90"
	47
[3]	48	!! * Substitutions
	49	# include "domzgr_substitute.h90"
	50	# include "vectopt_loop_substitute.h90"
	51	!!----------------------------------------------------------------------
[2528]	52	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]	53	!! $Id$
[2528]	54	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	55	!!----------------------------------------------------------------------
	56	CONTAINS
	57
[2528]	58	SUBROUTINE tra_adv_tvd ( kt, cdtype, p2dt, pun, pvn, pwn, &
	59	& ptb, ptn, pta, kjpt )
[3]	60	!!----------------------------------------------------------------------
	61	!! * ROUTINE tra_adv_tvd *
	62	!!
	63	!! ** Purpose : Compute the now trend due to total advection of
	64	!! tracers and add it to the general trend of tracer equations
	65	!!
	66	!! ** Method : TVD scheme, i.e. 2nd order centered scheme with
	67	!! corrected flux (monotonic correction)
	68	!! note: - this advection scheme needs a leap-frog time scheme
	69	!!
[2528]	70	!! ** Action : - update (pta) with the now advective tracer trends
	71	!! - save the trends
[503]	72	!!----------------------------------------------------------------------
[2715]	73	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	74	USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace
	75	USE wrk_nemo, ONLY: zwi => wrk_3d_12 , zwz => wrk_3d_13 ! 3D workspace
[3211]	76
	77	!! DCSE_NEMO: need additional directives for renamed module variables
	78	!FTRANS zwx zwy zwi zwz :I :I :z
	79
[2715]	80	!
[2528]	81	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	82	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	83	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	84	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
[3211]	85
	86	!! DCSE_NEMO: This style defeats ftrans
	87	! REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
	88	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
	89	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
	90
	91	!FTRANS pun pvn pwn :I :I :z
	92	!FTRANS ptb ptn pta :I :I :z :
	93	REAL(wp), INTENT(in ) :: pun(jpi,jpj,jpk) ! ocean velocity component (u)
	94	REAL(wp), INTENT(in ) :: pvn(jpi,jpj,jpk) ! ocean velocity component (v)
	95	REAL(wp), INTENT(in ) :: pwn(jpi,jpj,jpk) ! ocean velocity component (w)
	96	REAL(wp), INTENT(in ) :: ptb(jpi,jpj,jpk,kjpt) ! tracer fields (before)
	97	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpk,kjpt) ! tracer fields (now)
	98	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpk,kjpt) ! tracer trend
	99
[2715]	100	!
[2528]	101	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	102	REAL(wp) :: z2dtt, zbtr, ztra ! local scalar
	103	REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - -
	104	REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - -
	105	REAL(wp), DIMENSION (:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz
[3211]	106	!FTRANS ztrdx ztrdy ztrdz :I :I :z
	107
[3]	108	!!----------------------------------------------------------------------
	109
[2715]	110	IF( wrk_in_use(3, 12,13) ) THEN
	111	CALL ctl_stop('tra_adv_tvd: requested workspace arrays unavailable') ; RETURN
	112	ENDIF
	113
[2528]	114	IF( kt == nit000 ) THEN
	115	IF(lwp) WRITE(numout,*)
	116	IF(lwp) WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype
	117	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	118	!
	119	l_trd = .FALSE.
	120	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
[3]	121	ENDIF
[2528]	122	!
	123	IF( l_trd ) THEN
	124	ALLOCATE( ztrdx(jpi,jpj,jpk) ) ; ztrdx(:,:,:) = 0.e0
	125	ALLOCATE( ztrdy(jpi,jpj,jpk) ) ; ztrdy(:,:,:) = 0.e0
	126	ALLOCATE( ztrdz(jpi,jpj,jpk) ) ; ztrdz(:,:,:) = 0.e0
	127	END IF
	128	!
	129	zwi(:,:,:) = 0.e0
	130	!
	131	! ! ===========
	132	DO jn = 1, kjpt ! tracer loop
	133	! ! ===========
	134	! 1. Bottom value : flux set to zero
	135	! ----------------------------------
	136	zwx(:,:,jpk) = 0.e0 ; zwz(:,:,jpk) = 0.e0
	137	zwy(:,:,jpk) = 0.e0 ; zwi(:,:,jpk) = 0.e0
[3]	138
[2528]	139	! 2. upstream advection with initial mass fluxes & intermediate update
	140	! --------------------------------------------------------------------
	141	! upstream tracer flux in the i and j direction
[3211]	142	#if defined key_z_first
	143	DO jj = 1, jpjm1
	144	DO ji = 1, jpim1
	145	DO jk = 1, jpkm1
	146	#else
[2528]	147	DO jk = 1, jpkm1
	148	DO jj = 1, jpjm1
	149	DO ji = 1, fs_jpim1 ! vector opt.
[3211]	150	#endif
[2528]	151	! upstream scheme
	152	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
	153	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
	154	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
	155	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
	156	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) )
	157	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) )
	158	END DO
[3]	159	END DO
	160	END DO
	161
[2528]	162	! upstream tracer flux in the k direction
	163	! Surface value
	164	IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! volume variable
	165	ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface
	166	ENDIF
	167	! Interior value
[3211]	168	#if defined key_z_first
	169	DO jj = 1, jpj
	170	DO ji = 1, jpi
	171	DO jk = 2, jpkm1
	172	#else
[2528]	173	DO jk = 2, jpkm1
	174	DO jj = 1, jpj
	175	DO ji = 1, jpi
[3211]	176	#endif
[2528]	177	zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
	178	zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
	179	zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )
	180	END DO
[3]	181	END DO
	182	END DO
	183
[2528]	184	! total advective trend
[3211]	185	#if defined key_z_first
	186	DO jj = 2, jpjm1
	187	DO ji = 2, jpim1
	188	DO jk = 1, jpkm1
	189	z2dtt = p2dt(jk)
	190	#else
[216]	191	DO jk = 1, jpkm1
[2528]	192	z2dtt = p2dt(jk)
[216]	193	DO jj = 2, jpjm1
	194	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	195	#endif
[503]	196	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
[2528]	197	! total intermediate advective trends
	198	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	199	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
	200	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
	201	! update and guess with monotonic sheme
	202	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
	203	zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)
[216]	204	END DO
	205	END DO
	206	END DO
[2528]	207	! ! Lateral boundary conditions on zwi (unchanged sign)
	208	CALL lbc_lnk( zwi, 'T', 1. )
	209
	210	! ! trend diagnostics (contribution of upstream fluxes)
	211	IF( l_trd ) THEN
	212	! store intermediate advective trends
	213	ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)
	214	END IF
	215	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	216	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
	217	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) )
	218	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) )
	219	ENDIF
	220
	221	! 3. antidiffusive flux : high order minus low order
	222	! --------------------------------------------------
	223	! antidiffusive flux on i and j
[3211]	224	#if defined key_z_first
	225	DO jj = 1, jpjm1
	226	DO ji = 1, jpim1
	227	DO jk = 1, jpkm1
	228	#else
[503]	229	DO jk = 1, jpkm1
[2528]	230	DO jj = 1, jpjm1
	231	DO ji = 1, fs_jpim1 ! vector opt.
[3211]	232	#endif
[2528]	233	zwx(ji,jj,jk) = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk)
	234	zwy(ji,jj,jk) = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk)
[503]	235	END DO
	236	END DO
	237	END DO
[2528]	238
	239	! antidiffusive flux on k
[3211]	240	#if defined key_z_first
	241	DO jj = 1, jpj
	242	DO ji = 1, jpi
	243	zwz(ji,jj,1) = 0.e0 ! Surface value
	244	DO jk = 2, jpkm1
	245	#else
	246	zwz(:,:,1) = 0.e0 ! Surface value
[503]	247	!
[3211]	248	DO jk = 2, jpkm1 ! Interior value
[2528]	249	DO jj = 1, jpj
	250	DO ji = 1, jpi
[3211]	251	#endif
[2528]	252	zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk)
[503]	253	END DO
	254	END DO
	255	END DO
[2528]	256	CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions
	257	CALL lbc_lnk( zwz, 'W', 1. )
[216]	258
[2528]	259	! 4. monotonicity algorithm
	260	! -------------------------
	261	CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
[3]	262
	263
[2528]	264	! 5. final trend with corrected fluxes
	265	! ------------------------------------
[3211]	266	#if defined key_z_first
	267	DO jj = 2, jpjm1
	268	DO ji = 2, jpim1
	269	DO jk = 1, jpkm1
	270	#else
[216]	271	DO jk = 1, jpkm1
	272	DO jj = 2, jpjm1
[2528]	273	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	274	#endif
[503]	275	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
[2528]	276	! total advective trends
	277	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	278	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
	279	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
	280	! add them to the general tracer trends
	281	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[216]	282	END DO
	283	END DO
	284	END DO
[2528]	285
	286	! ! trend diagnostics (contribution of upstream fluxes)
	287	IF( l_trd ) THEN
	288	ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed
	289	ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed
	290	ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed
	291
	292	CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, ztrdx, pun, ptn(:,:,:,jn) )
	293	CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, ztrdy, pvn, ptn(:,:,:,jn) )
	294	CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, ztrdz, pwn, ptn(:,:,:,jn) )
	295	END IF
	296	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	297	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
	298	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) ) + htr_adv(:)
	299	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) ) + str_adv(:)
	300	ENDIF
[503]	301	!
[2715]	302	END DO
[503]	303	!
[2528]	304	IF( l_trd ) THEN
	305	DEALLOCATE( ztrdx ) ; DEALLOCATE( ztrdy ) ; DEALLOCATE( ztrdz )
	306	END IF
	307	!
[2715]	308	IF( wrk_not_released(3, 12,13) ) CALL ctl_stop('tra_adv_tvd: failed to release workspace arrays')
	309	!
[3211]	310
	311	!! * Reset control of array index permutation
	312	!FTRANS CLEAR
	313	# include "oce_ftrans.h90"
	314	# include "dom_oce_ftrans.h90"
	315	# include "trc_oce_ftrans.h90"
	316
[3]	317	END SUBROUTINE tra_adv_tvd
	318
	319
[2528]	320	SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
[3]	321	!!---------------------------------------------------------------------
	322	!! * ROUTINE nonosc *
	323	!!
	324	!! ** Purpose : compute monotonic tracer fluxes from the upstream
	325	!! scheme and the before field by a nonoscillatory algorithm
	326	!!
	327	!! ** Method : ... ???
	328	!! warning : pbef and paft must be masked, but the boundaries
	329	!! conditions on the fluxes are not necessary zalezak (1979)
	330	!! drange (1995) multi-dimensional forward-in-time and upstream-
	331	!! in-space based differencing for fluid
	332	!!----------------------------------------------------------------------
[2715]	333	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	334	USE wrk_nemo, ONLY: zbetup => wrk_3d_8 , zbetdo => wrk_3d_9 ! 3D workspace
	335	USE wrk_nemo, ONLY: zbup => wrk_3d_10 , zbdo => wrk_3d_11 ! - -
[3211]	336
	337	!! DCSE_NEMO: need additional directives for renamed module variables
	338	!FTRANS zbetup zbetdo zbup zbdo :I :I :z
	339
[2715]	340	!
[2528]	341	REAL(wp), DIMENSION(jpk) , INTENT(in ) :: p2dt ! vertical profile of tracer time-step
[3211]	342
	343	!! DCSE_NEMO: This style defeats ftrans
	344	! REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field
	345	! REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions
	346
	347	!FTRANS pbef paft :I :I :z
	348	!FTRANS paa pbb pcc :I :I :z
	349	REAL(wp), INTENT(in ) :: pbef(jpi,jpj,jpk), paft(jpi,jpj,jpk) ! before & after field
	350	REAL(wp), INTENT(inout) :: paa(jpi,jpj,jpk) ! monotonic fluxes in the 1st direction
	351	REAL(wp), INTENT(inout) :: pbb(jpi,jpj,jpk) ! monotonic fluxes in the 2nd direction
	352	REAL(wp), INTENT(inout) :: pcc(jpi,jpj,jpk) ! monotonic fluxes in the 3rd direction
[2715]	353	!
	354	INTEGER :: ji, jj, jk ! dummy loop indices
	355	INTEGER :: ikm1 ! local integer
	356	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt ! local scalars
	357	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -
[3]	358	!!----------------------------------------------------------------------
	359
[2715]	360	IF( wrk_in_use(3, 8,9,10,11) ) THEN
	361	CALL ctl_stop('nonosc: requested workspace array unavailable') ; RETURN
	362	ENDIF
[3]	363
[2715]	364	zbig = 1.e+40_wp
	365	zrtrn = 1.e-15_wp
	366	zbetup(:,:,jpk) = 0._wp ; zbetdo(:,:,jpk) = 0._wp
[785]	367
[2715]	368
[3]	369	! Search local extrema
	370	! --------------------
[785]	371	! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
	372	zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ), &
	373	& paft * tmask - zbig * ( 1.e0 - tmask ) )
	374	zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ), &
	375	& paft * tmask + zbig * ( 1.e0 - tmask ) )
	376
[3211]	377	#if defined key_z_first
	378	DO jj = 2, jpjm1
	379	DO ji = 2, jpim1
	380	DO jk = 1, jpkm1
	381	ikm1 = MAX(jk-1,1)
	382	z2dtt = p2dt(jk)
	383	#else
[3]	384	DO jk = 1, jpkm1
	385	ikm1 = MAX(jk-1,1)
[2528]	386	z2dtt = p2dt(jk)
[3]	387	DO jj = 2, jpjm1
	388	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	389	#endif
[3]	390
[785]	391	! search maximum in neighbourhood
	392	zup = MAX( zbup(ji ,jj ,jk ), &
	393	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
	394	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), &
	395	& zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) )
[3]	396
[785]	397	! search minimum in neighbourhood
	398	zdo = MIN( zbdo(ji ,jj ,jk ), &
	399	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
	400	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), &
	401	& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
[3]	402
[785]	403	! positive part of the flux
[3]	404	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
	405	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
	406	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
[785]	407
	408	! negative part of the flux
[3]	409	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
	410	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
	411	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
[785]	412
[3]	413	! up & down beta terms
	414	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
[785]	415	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
	416	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
[3]	417	END DO
	418	END DO
	419	END DO
[2528]	420	CALL lbc_lnk( zbetup, 'T', 1. ) ; CALL lbc_lnk( zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
[3]	421
	422
	423
[237]	424	! 3. monotonic flux in the i & j direction (paa & pbb)
	425	! ----------------------------------------
[3211]	426	#if defined key_z_first
	427	DO jj = 2, jpjm1
	428	DO ji = 2, jpim1
	429	DO jk = 1, jpkm1
	430	#else
[3]	431	DO jk = 1, jpkm1
	432	DO jj = 2, jpjm1
	433	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	434	#endif
[785]	435	zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
	436	zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
	437	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
	438	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
[3]	439
[785]	440	zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
	441	zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
	442	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
	443	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv )
[3]	444
	445	! monotonic flux in the k direction, i.e. pcc
	446	! -------------------------------------------
[785]	447	za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
	448	zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
	449	zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
	450	pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1.e0 - zc) * zb )
[3]	451	END DO
	452	END DO
	453	END DO
[2528]	454	CALL lbc_lnk( paa, 'U', -1. ) ; CALL lbc_lnk( pbb, 'V', -1. ) ! lateral boundary condition (changed sign)
[503]	455	!
[2715]	456	IF( wrk_not_released(3, 8,9,10,11) ) CALL ctl_stop('nonosc: failed to release workspace arrays')
	457	!
[3]	458	END SUBROUTINE nonosc
	459
	460	!!======================================================================
	461	END MODULE traadv_tvd

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

Download in other formats: