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_qck.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_qck.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: 29.5 KB

Rev	Line
[1231]	1	MODULE traadv_qck
	2	!!==============================================================================
	3	!! * MODULE traadv_qck *
[2528]	4	!! Ocean tracers: horizontal & vertical advective trend
[1231]	5	!!==============================================================================
[1559]	6	!! History : 3.0 ! 2008-07 (G. Reffray) Original code
[2528]	7	!! 3.3 ! 2010-05 (C.Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
[1231]	8	!!----------------------------------------------------------------------
	9
	10	!!----------------------------------------------------------------------
[2528]	11	!! tra_adv_qck : update the tracer trend with the horizontal advection
	12	!! trends using a 3rd order finite difference scheme
	13	!! tra_adv_qck_i : apply QUICK scheme in i-direction
	14	!! tra_adv_qck_j : apply QUICK scheme in j-direction
[1559]	15	!! tra_adv_cen2_k : 2nd centered scheme for the vertical advection
[1231]	16	!!----------------------------------------------------------------------
	17	USE oce ! ocean dynamics and active tracers
	18	USE dom_oce ! ocean space and time domain
[2528]	19	USE trdmod_oce ! ocean space and time domain
	20	USE trdtra ! ocean tracers trends
[1231]	21	USE trabbl ! advective term in the BBL
	22	USE lib_mpp ! distribued memory computing
	23	USE lbclnk ! ocean lateral boundary condition (or mpp link)
	24	USE dynspg_oce ! surface pressure gradient variables
	25	USE in_out_manager ! I/O manager
	26	USE diaptr ! poleward transport diagnostics
[2528]	27	USE trc_oce ! share passive tracers/Ocean variables
[1231]	28
	29	IMPLICIT NONE
	30	PRIVATE
	31
[1559]	32	PUBLIC tra_adv_qck ! routine called by step.F90
[1231]	33
[2528]	34	LOGICAL :: l_trd ! flag to compute trends
	35	REAL(wp) :: r1_6 = 1./ 6. ! 1/6 ratio
[1559]	36
[3211]	37	!! * Control permutation of array indices
	38	# include "oce_ftrans.h90"
	39	# include "dom_oce_ftrans.h90"
	40	# include "trc_oce_ftrans.h90"
	41
[1231]	42	!! * Substitutions
	43	# include "domzgr_substitute.h90"
	44	# include "vectopt_loop_substitute.h90"
	45	!!----------------------------------------------------------------------
[2528]	46	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1231]	47	!! $Id$
[2528]	48	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1231]	49	!!----------------------------------------------------------------------
	50	CONTAINS
	51
[2528]	52	SUBROUTINE tra_adv_qck ( kt, cdtype, p2dt, pun, pvn, pwn, &
	53	& ptb, ptn, pta, kjpt )
[1231]	54	!!----------------------------------------------------------------------
	55	!! * ROUTINE tra_adv_qck *
	56	!!
	57	!! ** Purpose : Compute the now trend due to the advection of tracers
	58	!! and add it to the general trend of passive tracer equations.
	59	!!
	60	!! ** Method : The advection is evaluated by a third order scheme
[1559]	61	!! For a positive velocity u : u(i)>0
	62	!! \|--FU--\|--FC--\|--FD--\|------\|
	63	!! i-1 i i+1 i+2
[1231]	64	!!
[1559]	65	!! For a negative velocity u : u(i)<0
	66	!! \|------\|--FD--\|--FC--\|--FU--\|
	67	!! i-1 i i+1 i+2
	68	!! where FU is the second upwind point
	69	!! FD is the first douwning point
	70	!! FC is the central point (or the first upwind point)
[1231]	71	!!
[1559]	72	!! Flux(i) = u(i) * { 0.5(FC+FD) -0.5C(i)(FD-FC) -((1-C(i))/6)(FU+FD-2FC) }
	73	!! with C(i)=\|u(i)\|dx(i)/dt (=Courant number)
[1231]	74	!!
	75	!! dt = 2*rdtra and the scalar values are tb and sb
	76	!!
[2528]	77	!! On the vertical, the simple centered scheme used ptn
[1231]	78	!!
[1559]	79	!! The fluxes are bounded by the ULTIMATE limiter to
	80	!! guarantee the monotonicity of the solution and to
[1231]	81	!! prevent the appearance of spurious numerical oscillations
	82	!!
[2528]	83	!! ** Action : - update (pta) with the now advective tracer trends
	84	!! - save the trends
[1231]	85	!!
	86	!! ** Reference : Leonard (1979, 1991)
	87	!!----------------------------------------------------------------------
[2528]	88	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	89	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	90	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	91	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
[3211]	92
	93	!! DCSE_NEMO: This style defeats ftrans
	94	! REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
	95	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
	96	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
	97
	98	!FTRANS pun pvn pwn :I :I :z
	99	!FTRANS ptb ptn :I :I :z :
	100	!FTRANS pta :I :I :z :
	101	REAL(wp), INTENT(in ) :: pun(jpi,jpj,jpk) ! ocean velocity component (u)
	102	REAL(wp), INTENT(in ) :: pvn(jpi,jpj,jpk) ! ocean velocity component (v)
	103	REAL(wp), INTENT(in ) :: pwn(jpi,jpj,jpk) ! ocean velocity component (w)
	104	REAL(wp), INTENT(in ) :: ptb(jpi,jpj,jpk,kjpt) ! tracer fields (before)
	105	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpk,kjpt) ! tracer fields (now)
	106	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpk,kjpt) ! tracer trend
	107
[1231]	108	!!----------------------------------------------------------------------
	109
[2528]	110	IF( kt == nit000 ) THEN
[1231]	111	IF(lwp) WRITE(numout,*)
[2528]	112	IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
[1231]	113	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
	114	IF(lwp) WRITE(numout,*)
[2528]	115	!
	116	l_trd = .FALSE.
	117	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
[1231]	118	ENDIF
	119
	120	! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
[2528]	121	CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt )
	122	CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt )
[1231]	123
	124	! II. The vertical fluxes are computed with the 2nd order centered scheme
[2528]	125	CALL tra_adv_cen2_k( kt, cdtype, pwn, ptn, pta, kjpt )
[1231]	126	!
[3211]	127
	128	!! * Reset control of array index permutation
	129	!FTRANS CLEAR
	130	# include "oce_ftrans.h90"
	131	# include "dom_oce_ftrans.h90"
	132	# include "trc_oce_ftrans.h90"
	133
[1231]	134	END SUBROUTINE tra_adv_qck
	135
	136
[2528]	137	SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun, &
	138	& ptb, ptn, pta, kjpt )
[1231]	139	!!----------------------------------------------------------------------
	140	!!
	141	!!----------------------------------------------------------------------
[2715]	142	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	143	USE oce , ONLY: zwx => ua ! ua used as workspace
	144	USE wrk_nemo, ONLY: zfu => wrk_3d_1 , zfc => wrk_3d_2, zfd => wrk_3d_3 ! 3D workspace
[3211]	145
	146	!! DCSE_NEMO: need additional directives for renamed module variables
	147	!FTRANS zwx zfu zfc zfd :I :I :z
	148
[2715]	149	!
	150	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	151	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	152	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	153	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
[3211]	154
	155	!! DCSE_NEMO: This style defeats ftrans
	156	! REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun ! i-velocity components
	157	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
	158	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
	159
	160	!FTRANS pun :I :I :z
	161	!FTRANS ptb ptn pta :I :I :z :
	162	REAL(wp), INTENT(in ) :: pun(jpi,jpj,jpk) ! i-velocity component
	163	REAL(wp), INTENT(in ) :: ptb(jpi,jpj,jpk,kjpt) ! tracer field (before)
	164	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpk,kjpt) ! tracer field (now)
	165	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpk,kjpt) ! tracer trend
	166
[2528]	167	!!
[2715]	168	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	169	REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk ! local scalars
[1231]	170	!----------------------------------------------------------------------
[2715]	171	!
	172	IF( wrk_in_use(3, 1,2,3) ) THEN
	173	CALL ctl_stop('tra_adv_qck_i: requested workspace arrays unavailable') ; RETURN
	174	ENDIF
[2528]	175	! ! ===========
	176	DO jn = 1, kjpt ! tracer loop
	177	! ! ===========
	178	zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0
	179	zfd(:,:,:) = 0.0 ; zwx(:,:,:) = 0.0
	180	!
[3211]	181	#if defined key_z_first
	182	!--- Computation of the upstream and downstream value of the tracer and the mask
	183	DO jj = 2, jpjm1
	184	DO ji = 2, jpim1
	185	DO jk = 1, jpkm1
	186	#else
[2528]	187	DO jk = 1, jpkm1
	188	!
[3211]	189	!--- Computation of the upstream and downstream value of the tracer and the mask
[2528]	190	DO jj = 2, jpjm1
	191	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	192	#endif
[2528]	193	! Upstream in the x-direction for the tracer
	194	zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)
	195	! Downstream in the x-direction for the tracer
	196	zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)
	197	END DO
[1559]	198	END DO
	199	END DO
[2528]	200	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) ! Lateral boundary conditions
	201
[1231]	202	!
	203	! Horizontal advective fluxes
	204	! ---------------------------
	205	!
[3211]	206	#if defined key_z_first
	207	DO jj = 2, jpjm1
	208	DO ji = 2, jpim1
	209	DO jk = 1, jpkm1
	210	#else
[2528]	211	DO jk = 1, jpkm1
	212	DO jj = 2, jpjm1
	213	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	214	#endif
[2528]	215	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	216	zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T
	217	END DO
	218	END DO
[1559]	219	END DO
[1231]	220	!
[3211]	221	#if defined key_z_first
	222	DO jj = 2, jpjm1
	223	DO ji = 2, jpim1
	224	DO jk = 1, jpkm1
	225	zdt = p2dt(jk)
	226	#else
[2528]	227	DO jk = 1, jpkm1
	228	zdt = p2dt(jk)
	229	DO jj = 2, jpjm1
	230	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	231	#endif
[2528]	232	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	233	zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)
	234	zwx(ji,jj,jk) = ABS( pun(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)
	235	zfc(ji,jj,jk) = zdir * ptb(ji ,jj,jk,jn) + ( 1. - zdir ) * ptb(ji+1,jj,jk,jn) ! FC in the x-direction for T
	236	zfd(ji,jj,jk) = zdir * ptb(ji+1,jj,jk,jn) + ( 1. - zdir ) * ptb(ji ,jj,jk,jn) ! FD in the x-direction for T
	237	END DO
	238	END DO
	239	END DO
	240	!--- Lateral boundary conditions
	241	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
	242	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
	243
[1231]	244	!--- QUICKEST scheme
[2528]	245	CALL quickest( zfu, zfd, zfc, zwx )
[1231]	246	!
[2528]	247	! Mask at the T-points in the x-direction (mask=0 or mask=1)
[3211]	248	#if defined key_z_first
	249	DO jj = 2, jpjm1
	250	DO ji = 2, jpim1
	251	DO jk = 1, jpkm1
	252	#else
[2528]	253	DO jk = 1, jpkm1
	254	DO jj = 2, jpjm1
	255	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	256	#endif
[2528]	257	zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
[2715]	258	END DO
[1231]	259	END DO
	260	END DO
[2528]	261	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ! Lateral boundary conditions
	262
[1231]	263	!
[2528]	264	! Tracer flux on the x-direction
[3211]	265	#if defined key_z_first
	266	DO jj = 2, jpjm1
	267	DO ji = 2, jpim1
	268	DO jk = 1, jpkm1
	269	#else
[2528]	270	DO jk = 1, jpkm1
[1231]	271	DO jj = 2, jpjm1
[2528]	272	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	273	#endif
[2528]	274	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	275	!--- If the second ustream point is a land point
	276	!--- the flux is computed by the 1st order UPWIND scheme
	277	zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
	278	zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
	279	zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
[1231]	280	END DO
	281	END DO
[3211]	282	#if defined key_z_first
	283	END DO
	284	! Computation of the trend
	285	DO jj = 2, jpjm1
	286	DO ji = 2, jpim1
	287	DO jk = 1, jpkm1
	288	#else
[2528]	289	!
	290	! Computation of the trend
	291	DO jj = 2, jpjm1
	292	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	293	#endif
[2528]	294	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	295	! horizontal advective trends
	296	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
	297	!--- add it to the general tracer trends
	298	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
	299	END DO
	300	END DO
	301	!
[1231]	302	END DO
[2528]	303	! ! trend diagnostics (contribution of upstream fluxes)
	304	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
	305	!
	306	END DO
	307	!
[2715]	308	IF( wrk_not_released(3, 1,2,3) ) CALL ctl_stop('tra_adv_qck_i: 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
[1559]	317	END SUBROUTINE tra_adv_qck_i
[1231]	318
	319
[2528]	320	SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn, &
	321	& ptb, ptn, pta, kjpt )
[1231]	322	!!----------------------------------------------------------------------
	323	!!
	324	!!----------------------------------------------------------------------
[2715]	325	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	326	USE oce , ONLY: zwy => ua ! ua used as workspace
	327	USE wrk_nemo, ONLY: zfu => wrk_3d_1 , zfc => wrk_3d_2, zfd => wrk_3d_3 ! 3D workspace
[3211]	328
	329	!! DCSE_NEMO: need additional directives for renamed module variables
	330	!FTRANS zwy zfu zfc zfd :I :I :z
	331
[2715]	332	!
	333	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	334	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	335	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	336	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
[3211]	337
	338	!! DCSE_NEMO: This style defeats ftrans
	339	! REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pvn ! j-velocity components
	340	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
	341	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
	342
	343	!FTRANS pvn :I :I :z
	344	!FTRANS ptb ptn pta :I :I :z :
	345	REAL(wp), INTENT(in ) :: pvn(jpi,jpj,jpk) ! j-velocity component
	346	REAL(wp), INTENT(in ) :: ptb(jpi,jpj,jpk,kjpt) ! tracer field (before)
	347	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpk,kjpt) ! tracer field (now)
	348	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpk,kjpt) ! tracer trend
	349
[1559]	350	!!
[2715]	351	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	352	REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk ! local scalars
[1231]	353	!----------------------------------------------------------------------
[2715]	354	!
	355	IF(wrk_in_use(3, 1,2,3))THEN
	356	CALL ctl_stop('tra_adv_qck_j: ERROR: requested workspace arrays unavailable')
	357	RETURN
	358	END IF
[2528]	359	! ! ===========
	360	DO jn = 1, kjpt ! tracer loop
	361	! ! ===========
	362	zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0
	363	zfd(:,:,:) = 0.0 ; zwy(:,:,:) = 0.0
	364	!
[3211]	365	#if defined key_z_first
	366	!--- Computation of the ustream and downstream value of the tracer and the mask
	367	DO jj = 2, jpjm1
	368	DO ji = 2, jpim1
	369	DO jk = 1, jpkm1
	370	#else
[2528]	371	DO jk = 1, jpkm1
	372	!
	373	!--- Computation of the ustream and downstream value of the tracer and the mask
	374	DO jj = 2, jpjm1
	375	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	376	#endif
[2528]	377	! Upstream in the x-direction for the tracer
	378	zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
	379	! Downstream in the x-direction for the tracer
	380	zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
	381	END DO
[1559]	382	END DO
	383	END DO
[2528]	384	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) ! Lateral boundary conditions
	385
	386
[1231]	387	!
	388	! Horizontal advective fluxes
	389	! ---------------------------
	390	!
[3211]	391	#if defined key_z_first
	392	DO jj = 2, jpjm1
	393	DO ji = 2, jpim1
	394	DO jk = 1, jpkm1
	395	#else
[2528]	396	DO jk = 1, jpkm1
	397	DO jj = 2, jpjm1
	398	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	399	#endif
[2528]	400	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	401	zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T
	402	END DO
[1559]	403	END DO
	404	END DO
[1231]	405	!
[3211]	406	#if defined key_z_first
	407	DO jj = 2, jpjm1
	408	DO ji = 2, jpim1
	409	DO jk = 1, jpkm1
	410	zdt = p2dt(jk)
	411	#else
[2528]	412	DO jk = 1, jpkm1
	413	zdt = p2dt(jk)
	414	DO jj = 2, jpjm1
	415	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	416	#endif
[2528]	417	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	418	zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
	419	zwy(ji,jj,jk) = ABS( pvn(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)
	420	zfc(ji,jj,jk) = zdir * ptb(ji,jj ,jk,jn) + ( 1. - zdir ) * ptb(ji,jj+1,jk,jn) ! FC in the x-direction for T
	421	zfd(ji,jj,jk) = zdir * ptb(ji,jj+1,jk,jn) + ( 1. - zdir ) * ptb(ji,jj ,jk,jn) ! FD in the x-direction for T
	422	END DO
	423	END DO
	424	END DO
	425
	426	!--- Lateral boundary conditions
	427	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
	428	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
	429
[1231]	430	!--- QUICKEST scheme
[2528]	431	CALL quickest( zfu, zfd, zfc, zwy )
[1231]	432	!
[2528]	433	! Mask at the T-points in the x-direction (mask=0 or mask=1)
[3211]	434	#if defined key_z_first
	435	DO jj = 2, jpjm1
	436	DO ji = 2, jpim1
	437	DO jk = 1, jpkm1
	438	#else
[2528]	439	DO jk = 1, jpkm1
	440	DO jj = 2, jpjm1
	441	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	442	#endif
[2528]	443	zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
	444	END DO
[1231]	445	END DO
	446	END DO
[2528]	447	!--- Lateral boundary conditions
	448	CALL lbc_lnk( zfu(:,:,:), 'T', 1. )
	449	!
	450	! Tracer flux on the x-direction
[3211]	451	#if defined key_z_first
	452	DO jj = 2, jpjm1
	453	DO ji = 2, jpim1
	454	DO jk = 1, jpkm1
	455	#else
[2528]	456	DO jk = 1, jpkm1
	457	!
[1231]	458	DO jj = 2, jpjm1
[2528]	459	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	460	#endif
[2528]	461	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	462	!--- If the second ustream point is a land point
	463	!--- the flux is computed by the 1st order UPWIND scheme
	464	zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
	465	zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
	466	zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
[1231]	467	END DO
	468	END DO
[3211]	469	#if defined key_z_first
	470	END DO
	471	! Computation of the trend
	472	DO jj = 2, jpjm1
	473	DO ji = 2, jpim1
	474	DO jk = 1, jpkm1
	475	#else
[2528]	476	!
	477	! Computation of the trend
	478	DO jj = 2, jpjm1
	479	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	480	#endif
[2528]	481	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	482	! horizontal advective trends
	483	ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
	484	!--- add it to the general tracer trends
	485	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[1231]	486	END DO
	487	END DO
[2528]	488	!
	489	END DO
	490	! ! trend diagnostics (contribution of upstream fluxes)
	491	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
	492	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	493	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
	494	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) )
	495	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) )
[1231]	496	ENDIF
[2528]	497	!
	498	END DO
	499	!
[2715]	500	IF( wrk_not_released(3, 1,2,3) ) CALL ctl_stop('tra_adv_qck_j: failed to release workspace arrays')
	501	!
[3211]	502
	503	!! * Reset control of array index permutation
	504	!FTRANS CLEAR
	505	# include "oce_ftrans.h90"
	506	# include "dom_oce_ftrans.h90"
	507	# include "trc_oce_ftrans.h90"
	508
[1559]	509	END SUBROUTINE tra_adv_qck_j
[1231]	510
	511
[2528]	512	SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn, &
	513	& ptn, pta, kjpt )
[1231]	514	!!----------------------------------------------------------------------
	515	!!
	516	!!----------------------------------------------------------------------
[2715]	517	USE oce, ONLY: zwz => ua ! ua used as workspace
[3211]	518
	519	!! DCSE_NEMO: need additional directives for renamed module variables
	520	!FTRANS zwz :I :I :z
	521
[2715]	522	!
	523	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	524	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	525	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[3211]	526
	527	!! DCSE_NEMO: This style defeats ftrans
	528	! REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pwn ! vertical velocity
	529	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptn ! tracer fields (now)
	530	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
	531
	532	!FTRANS pwn :I :I :z
	533	!FTRANS ptn pta :I :I :z :
	534	REAL(wp), INTENT(in ) :: pwn(jpi,jpj,jpk) ! vertical velocity
	535	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpk,kjpt) ! tracer fields (now)
	536	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpk,kjpt) ! tracer trend
	537
[2715]	538	!
[2528]	539	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[2715]	540	REAL(wp) :: zbtr , ztra ! local scalars
[1559]	541	!!----------------------------------------------------------------------
[2528]	542
	543	! ! ===========
	544	DO jn = 1, kjpt ! tracer loop
	545	! ! ===========
	546	! 1. Bottom value : flux set to zero
	547	zwz(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero
	548	!
	549	! ! Surface value
	550	IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! Variable volume : flux set to zero
	551	ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn) ! Constant volume : advective flux through the surface
	552	ENDIF
	553	!
[3211]	554	#if defined key_z_first
	555	DO jj = 2, jpjm1
	556	DO ji = 2, jpim1
	557	DO jk = 2, jpkm1 ! Interior point: second order centered tracer flux at w-point
	558	#else
[2528]	559	DO jk = 2, jpkm1 ! Interior point: second order centered tracer flux at w-point
	560	DO jj = 2, jpjm1
	561	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	562	#endif
[2528]	563	zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
	564	END DO
[1231]	565	END DO
	566	END DO
[2528]	567	!
[3211]	568	#if defined key_z_first
	569	DO jj = 2, jpjm1
	570	DO ji = 2, jpim1
	571	DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==!
	572	#else
[2528]	573	DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==!
	574	DO jj = 2, jpjm1
	575	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	576	#endif
[2528]	577	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	578	! k- vertical advective trends
	579	ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )
	580	! added to the general tracer trends
	581	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
	582	END DO
[1231]	583	END DO
	584	END DO
[2528]	585	! ! Save the vertical advective trends for diagnostic
	586	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) )
	587	!
[1231]	588	END DO
	589	!
[3211]	590
	591	!! * Reset control of array index permutation
	592	!FTRANS CLEAR
	593	# include "oce_ftrans.h90"
	594	# include "dom_oce_ftrans.h90"
	595	# include "trc_oce_ftrans.h90"
	596
[1559]	597	END SUBROUTINE tra_adv_cen2_k
[1231]	598
	599
[2528]	600	SUBROUTINE quickest( pfu, pfd, pfc, puc )
[1231]	601	!!----------------------------------------------------------------------
	602	!!
[2528]	603	!! ** Purpose : Computation of advective flux with Quickest scheme
	604	!!
	605	!! ** Method :
[1231]	606	!!----------------------------------------------------------------------
[3211]	607
	608	!! DCSE_NEMO: This style defeats ftrans
	609
	610	! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfu ! second upwind point
	611	! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfd ! first douwning point
	612	! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfc ! the central point (or the first upwind point)
	613	! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: puc ! input as Courant number ; output as flux
	614
	615	!FTRANS pfu pfd pfc puc :I :I :z
	616	REAL(wp), INTENT(in ) :: pfu(jpi,jpj,jpk) ! second upwind point
	617	REAL(wp), INTENT(in ) :: pfd(jpi,jpj,jpk) ! first douwning point
	618	REAL(wp), INTENT(in ) :: pfc(jpi,jpj,jpk) ! the central point (or the first upwind point)
	619	REAL(wp), INTENT(inout) :: puc(jpi,jpj,jpk) ! input as Courant number ; output as flux
	620
[2528]	621	!!
	622	INTEGER :: ji, jj, jk ! dummy loop indices
	623	REAL(wp) :: zcoef1, zcoef2, zcoef3 ! local scalars
	624	REAL(wp) :: zc, zcurv, zfho ! - -
	625	!----------------------------------------------------------------------
	626
[3211]	627	#if defined key_z_first
	628	DO jj = 1, jpj
	629	DO ji = 1, jpi
	630	DO jk = 1, jpkm1
	631	#else
[2528]	632	DO jk = 1, jpkm1
	633	DO jj = 1, jpj
	634	DO ji = 1, jpi
[3211]	635	#endif
[2528]	636	zc = puc(ji,jj,jk) ! Courant number
	637	zcurv = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
	638	zcoef1 = 0.5 * ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
	639	zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
	640	zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
	641	zfho = zcoef1 - zcoef2 - zcoef3 ! phi_f QUICKEST
	642	!
	643	zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
	644	zcoef2 = ABS( zcoef1 )
	645	zcoef3 = ABS( zcurv )
	646	IF( zcoef3 >= zcoef2 ) THEN
	647	zfho = pfc(ji,jj,jk)
	648	ELSE
	649	zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) ) ! phi_REF
	650	IF( zcoef1 >= 0. ) THEN
	651	zfho = MAX( pfc(ji,jj,jk), zfho )
	652	zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) )
	653	ELSE
	654	zfho = MIN( pfc(ji,jj,jk), zfho )
	655	zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) )
	656	ENDIF
	657	ENDIF
	658	puc(ji,jj,jk) = zfho
	659	END DO
	660	END DO
	661	END DO
[1231]	662	!
	663	END SUBROUTINE quickest
	664
	665	!!======================================================================
	666	END MODULE traadv_qck

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

Download in other formats: