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

Last change on this file since 4409 was 4409, checked in by trackstand2, 10 years ago
Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.
Property svn:keywords set to `Id`
File size: 29.6 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
[4409]	91	REAL(wp), DIMENSION( jpkorig ), 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 :
[4409]	101	REAL(wp), INTENT(in ) :: pun(jpi,jpj,jpkorig) ! ocean velocity component (u)
	102	REAL(wp), INTENT(in ) :: pvn(jpi,jpj,jpkorig) ! ocean velocity component (v)
	103	REAL(wp), INTENT(in ) :: pwn(jpi,jpj,jpkorig) ! ocean velocity component (w)
	104	REAL(wp), INTENT(in ) :: ptb(jpi,jpj,jpkorig,kjpt) ! tracer fields (before)
	105	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpkorig,kjpt) ! tracer fields (now)
	106	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpkorig,kjpt) ! tracer trend
[3211]	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
[4409]	153	REAL(wp), DIMENSION( jpkorig ), 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 :
[4409]	162	REAL(wp), INTENT(in ) :: pun(jpi,jpj,jpkorig) ! i-velocity component
	163	REAL(wp), INTENT(in ) :: ptb(jpi,jpj,jpkorig,kjpt) ! tracer field (before)
	164	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpkorig,kjpt) ! tracer field (now)
	165	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpkorig,kjpt) ! tracer trend
[3211]	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
[4409]	336	REAL(wp), DIMENSION( jpkorig ), 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 :
[4409]	345	REAL(wp), INTENT(in ) :: pvn(jpi,jpj,jpkorig) ! j-velocity component
	346	REAL(wp), INTENT(in ) :: ptb(jpi,jpj,jpkorig,kjpt) ! tracer field (before)
	347	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpkorig,kjpt) ! tracer field (now)
	348	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpkorig,kjpt) ! tracer trend
[3211]	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 :
[4409]	534	REAL(wp), INTENT(in ) :: pwn(jpi,jpj,jpkorig) ! vertical velocity
	535	REAL(wp), INTENT(in ) :: ptn(jpi,jpj,jpkorig,kjpt) ! tracer fields (now)
	536	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpkorig,kjpt) ! tracer trend
[3211]	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
[4409]	616	REAL(wp), INTENT(in ) :: pfu(jpi,jpj,jpkorig) ! second upwind point
	617	REAL(wp), INTENT(in ) :: pfd(jpi,jpj,jpkorig) ! first douwning point
	618	REAL(wp), INTENT(in ) :: pfc(jpi,jpj,jpkorig) ! the central point (or the first upwind point)
	619	REAL(wp), INTENT(inout) :: puc(jpi,jpj,jpkorig) ! input as Courant number ; output as flux
[3211]	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: