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 trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 3389

Last change on this file since 3389 was 3301, checked in by cbricaud, 12 years ago
undo a bug correction, for ticket 924
Property svn:keywords set to `Id`
File size: 24.1 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
[3294]	28	USE wrk_nemo ! Memory Allocation
	29	USE timing ! Timing
[1231]	30
	31	IMPLICIT NONE
	32	PRIVATE
	33
[1559]	34	PUBLIC tra_adv_qck ! routine called by step.F90
[1231]	35
[2528]	36	LOGICAL :: l_trd ! flag to compute trends
	37	REAL(wp) :: r1_6 = 1./ 6. ! 1/6 ratio
[1559]	38
[1231]	39	!! * Substitutions
	40	# include "domzgr_substitute.h90"
	41	# include "vectopt_loop_substitute.h90"
	42	!!----------------------------------------------------------------------
[2528]	43	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1231]	44	!! $Id$
[2528]	45	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1231]	46	!!----------------------------------------------------------------------
	47	CONTAINS
	48
[3294]	49	SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn, &
[2528]	50	& ptb, ptn, pta, kjpt )
[1231]	51	!!----------------------------------------------------------------------
	52	!! * ROUTINE tra_adv_qck *
	53	!!
	54	!! ** Purpose : Compute the now trend due to the advection of tracers
	55	!! and add it to the general trend of passive tracer equations.
	56	!!
	57	!! ** Method : The advection is evaluated by a third order scheme
[1559]	58	!! For a positive velocity u : u(i)>0
	59	!! \|--FU--\|--FC--\|--FD--\|------\|
	60	!! i-1 i i+1 i+2
[1231]	61	!!
[1559]	62	!! For a negative velocity u : u(i)<0
	63	!! \|------\|--FD--\|--FC--\|--FU--\|
	64	!! i-1 i i+1 i+2
	65	!! where FU is the second upwind point
	66	!! FD is the first douwning point
	67	!! FC is the central point (or the first upwind point)
[1231]	68	!!
[1559]	69	!! Flux(i) = u(i) * { 0.5(FC+FD) -0.5C(i)(FD-FC) -((1-C(i))/6)(FU+FD-2FC) }
	70	!! with C(i)=\|u(i)\|dx(i)/dt (=Courant number)
[1231]	71	!!
	72	!! dt = 2*rdtra and the scalar values are tb and sb
	73	!!
[2528]	74	!! On the vertical, the simple centered scheme used ptn
[1231]	75	!!
[1559]	76	!! The fluxes are bounded by the ULTIMATE limiter to
	77	!! guarantee the monotonicity of the solution and to
[1231]	78	!! prevent the appearance of spurious numerical oscillations
	79	!!
[2528]	80	!! ** Action : - update (pta) with the now advective tracer trends
	81	!! - save the trends
[1231]	82	!!
	83	!! ** Reference : Leonard (1979, 1991)
	84	!!----------------------------------------------------------------------
[2528]	85	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[3294]	86	INTEGER , INTENT(in ) :: kit000 ! first time step index
[2528]	87	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	88	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	89	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
	90	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
	91	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
	92	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[1231]	93	!!----------------------------------------------------------------------
	94
[3294]	95	!
	96	IF( nn_timing == 1 ) CALL timing_start('tra_adv_qck')
	97	!
	98	IF( kt == kit000 ) THEN
[1231]	99	IF(lwp) WRITE(numout,*)
[2528]	100	IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
[1231]	101	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
	102	IF(lwp) WRITE(numout,*)
[2528]	103	!
	104	l_trd = .FALSE.
	105	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
[1231]	106	ENDIF
	107
	108	! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
[2528]	109	CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt )
	110	CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt )
[1231]	111
	112	! II. The vertical fluxes are computed with the 2nd order centered scheme
[2528]	113	CALL tra_adv_cen2_k( kt, cdtype, pwn, ptn, pta, kjpt )
[1231]	114	!
[3294]	115	IF( nn_timing == 1 ) CALL timing_stop('tra_adv_qck')
	116	!
[1231]	117	END SUBROUTINE tra_adv_qck
	118
	119
[2528]	120	SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun, &
	121	& ptb, ptn, pta, kjpt )
[1231]	122	!!----------------------------------------------------------------------
	123	!!
	124	!!----------------------------------------------------------------------
[2715]	125	USE oce , ONLY: zwx => ua ! ua used as workspace
	126	!
	127	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	128	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	129	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	130	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
	131	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun ! i-velocity components
	132	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
	133	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[2528]	134	!!
[2715]	135	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	136	REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk ! local scalars
[3294]	137	REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu, zfc, zfd
[1231]	138	!----------------------------------------------------------------------
[2715]	139	!
[3294]	140	CALL wrk_alloc( jpi, jpj, jpk, zfu, zfc, zfd )
[2528]	141	! ! ===========
	142	DO jn = 1, kjpt ! tracer loop
	143	! ! ===========
	144	zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0
	145	zfd(:,:,:) = 0.0 ; zwx(:,:,:) = 0.0
	146	!
	147	DO jk = 1, jpkm1
	148	!
	149	!--- Computation of the ustream and downstream value of the tracer and the mask
	150	DO jj = 2, jpjm1
	151	DO ji = fs_2, fs_jpim1 ! vector opt.
	152	! Upstream in the x-direction for the tracer
	153	zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)
	154	! Downstream in the x-direction for the tracer
	155	zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)
	156	END DO
[1559]	157	END DO
	158	END DO
[2528]	159	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) ! Lateral boundary conditions
	160
[1231]	161	!
	162	! Horizontal advective fluxes
	163	! ---------------------------
	164	!
[2528]	165	DO jk = 1, jpkm1
	166	DO jj = 2, jpjm1
	167	DO ji = fs_2, fs_jpim1 ! vector opt.
	168	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	169	zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T
	170	END DO
	171	END DO
[1559]	172	END DO
[1231]	173	!
[2528]	174	DO jk = 1, jpkm1
	175	zdt = p2dt(jk)
	176	DO jj = 2, jpjm1
	177	DO ji = fs_2, fs_jpim1 ! vector opt.
	178	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
[3301]	179	zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)
[2528]	180	zwx(ji,jj,jk) = ABS( pun(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)
	181	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
	182	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
	183	END DO
	184	END DO
	185	END DO
	186	!--- Lateral boundary conditions
	187	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
	188	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
	189
[1231]	190	!--- QUICKEST scheme
[2528]	191	CALL quickest( zfu, zfd, zfc, zwx )
[1231]	192	!
[2528]	193	! Mask at the T-points in the x-direction (mask=0 or mask=1)
	194	DO jk = 1, jpkm1
	195	DO jj = 2, jpjm1
	196	DO ji = fs_2, fs_jpim1 ! vector opt.
	197	zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
[2715]	198	END DO
[1231]	199	END DO
	200	END DO
[2528]	201	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ! Lateral boundary conditions
	202
[1231]	203	!
[2528]	204	! Tracer flux on the x-direction
	205	DO jk = 1, jpkm1
	206	!
[1231]	207	DO jj = 2, jpjm1
[2528]	208	DO ji = fs_2, fs_jpim1 ! vector opt.
	209	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	210	!--- If the second ustream point is a land point
	211	!--- the flux is computed by the 1st order UPWIND scheme
	212	zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
	213	zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
	214	zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
[1231]	215	END DO
	216	END DO
[3300]	217	END DO
	218	!
	219	CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
	220	!
	221	! Computation of the trend
	222	DO jk = 1, jpkm1
[2528]	223	DO jj = 2, jpjm1
	224	DO ji = fs_2, fs_jpim1 ! vector opt.
	225	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	226	! horizontal advective trends
	227	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
	228	!--- add it to the general tracer trends
	229	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
	230	END DO
	231	END DO
[1231]	232	END DO
[2528]	233	! ! trend diagnostics (contribution of upstream fluxes)
	234	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
	235	!
	236	END DO
	237	!
[3294]	238	CALL wrk_dealloc( jpi, jpj, jpk, zfu, zfc, zfd )
[2715]	239	!
[1559]	240	END SUBROUTINE tra_adv_qck_i
[1231]	241
	242
[2528]	243	SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn, &
	244	& ptb, ptn, pta, kjpt )
[1231]	245	!!----------------------------------------------------------------------
	246	!!
	247	!!----------------------------------------------------------------------
[2715]	248	USE oce , ONLY: zwy => ua ! ua used as workspace
	249	!
	250	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	251	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	252	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	253	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
	254	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pvn ! j-velocity components
	255	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
	256	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[1559]	257	!!
[2715]	258	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	259	REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk ! local scalars
[3294]	260	REAL(wp), POINTER, DIMENSION(:,:,:) :: zfu, zfc, zfd
[1231]	261	!----------------------------------------------------------------------
[2715]	262	!
[3294]	263	CALL wrk_alloc( jpi, jpj, jpk, zfu, zfc, zfd )
	264	!
[2528]	265	! ! ===========
	266	DO jn = 1, kjpt ! tracer loop
	267	! ! ===========
	268	zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0
	269	zfd(:,:,:) = 0.0 ; zwy(:,:,:) = 0.0
	270	!
	271	DO jk = 1, jpkm1
	272	!
	273	!--- Computation of the ustream and downstream value of the tracer and the mask
	274	DO jj = 2, jpjm1
	275	DO ji = fs_2, fs_jpim1 ! vector opt.
	276	! Upstream in the x-direction for the tracer
	277	zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
	278	! Downstream in the x-direction for the tracer
	279	zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
	280	END DO
[1559]	281	END DO
	282	END DO
[2528]	283	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) ! Lateral boundary conditions
	284
	285
[1231]	286	!
	287	! Horizontal advective fluxes
	288	! ---------------------------
	289	!
[2528]	290	DO jk = 1, jpkm1
	291	DO jj = 2, jpjm1
	292	DO ji = fs_2, fs_jpim1 ! vector opt.
	293	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	294	zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T
	295	END DO
[1559]	296	END DO
	297	END DO
[1231]	298	!
[2528]	299	DO jk = 1, jpkm1
	300	zdt = p2dt(jk)
	301	DO jj = 2, jpjm1
	302	DO ji = fs_2, fs_jpim1 ! vector opt.
	303	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
[3301]	304	zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
[2528]	305	zwy(ji,jj,jk) = ABS( pvn(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)
	306	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
	307	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
	308	END DO
	309	END DO
	310	END DO
	311
	312	!--- Lateral boundary conditions
	313	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
	314	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
	315
[1231]	316	!--- QUICKEST scheme
[2528]	317	CALL quickest( zfu, zfd, zfc, zwy )
[1231]	318	!
[2528]	319	! Mask at the T-points in the x-direction (mask=0 or mask=1)
	320	DO jk = 1, jpkm1
	321	DO jj = 2, jpjm1
	322	DO ji = fs_2, fs_jpim1 ! vector opt.
	323	zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
	324	END DO
[1231]	325	END DO
	326	END DO
[2528]	327	!--- Lateral boundary conditions
	328	CALL lbc_lnk( zfu(:,:,:), 'T', 1. )
	329	!
	330	! Tracer flux on the x-direction
	331	DO jk = 1, jpkm1
	332	!
[1231]	333	DO jj = 2, jpjm1
[2528]	334	DO ji = fs_2, fs_jpim1 ! vector opt.
	335	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	336	!--- If the second ustream point is a land point
	337	!--- the flux is computed by the 1st order UPWIND scheme
	338	zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
	339	zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
	340	zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
[1231]	341	END DO
	342	END DO
[3300]	343	END DO
	344	!
	345	CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
	346	!
	347	! Computation of the trend
	348	DO jk = 1, jpkm1
[2528]	349	DO jj = 2, jpjm1
	350	DO ji = fs_2, fs_jpim1 ! vector opt.
	351	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	352	! horizontal advective trends
	353	ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
	354	!--- add it to the general tracer trends
	355	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[1231]	356	END DO
	357	END DO
[2528]	358	END DO
	359	! ! trend diagnostics (contribution of upstream fluxes)
	360	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
	361	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	362	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
	363	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) )
	364	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) )
[1231]	365	ENDIF
[2528]	366	!
	367	END DO
	368	!
[3294]	369	CALL wrk_dealloc( jpi, jpj, jpk, zfu, zfc, zfd )
[2715]	370	!
[1559]	371	END SUBROUTINE tra_adv_qck_j
[1231]	372
	373
[2528]	374	SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn, &
	375	& ptn, pta, kjpt )
[1231]	376	!!----------------------------------------------------------------------
	377	!!
	378	!!----------------------------------------------------------------------
[2715]	379	USE oce, ONLY: zwz => ua ! ua used as workspace
	380	!
	381	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	382	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	383	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	384	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pwn ! vertical velocity
	385	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptn ! before and now tracer fields
	386	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
	387	!
[2528]	388	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[2715]	389	REAL(wp) :: zbtr , ztra ! local scalars
[1559]	390	!!----------------------------------------------------------------------
[2528]	391
	392	! ! ===========
	393	DO jn = 1, kjpt ! tracer loop
	394	! ! ===========
	395	! 1. Bottom value : flux set to zero
	396	zwz(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero
	397	!
	398	! ! Surface value
	399	IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! Variable volume : flux set to zero
	400	ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn) ! Constant volume : advective flux through the surface
	401	ENDIF
	402	!
	403	DO jk = 2, jpkm1 ! Interior point: second order centered tracer flux at w-point
	404	DO jj = 2, jpjm1
	405	DO ji = fs_2, fs_jpim1 ! vector opt.
	406	zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
	407	END DO
[1231]	408	END DO
	409	END DO
[2528]	410	!
	411	DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==!
	412	DO jj = 2, jpjm1
	413	DO ji = fs_2, fs_jpim1 ! vector opt.
	414	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	415	! k- vertical advective trends
	416	ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )
	417	! added to the general tracer trends
	418	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
	419	END DO
[1231]	420	END DO
	421	END DO
[2528]	422	! ! Save the vertical advective trends for diagnostic
	423	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) )
	424	!
[1231]	425	END DO
	426	!
[1559]	427	END SUBROUTINE tra_adv_cen2_k
[1231]	428
	429
[2528]	430	SUBROUTINE quickest( pfu, pfd, pfc, puc )
[1231]	431	!!----------------------------------------------------------------------
	432	!!
[2528]	433	!! ** Purpose : Computation of advective flux with Quickest scheme
	434	!!
	435	!! ** Method :
[1231]	436	!!----------------------------------------------------------------------
[2528]	437	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfu ! second upwind point
	438	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfd ! first douwning point
	439	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfc ! the central point (or the first upwind point)
	440	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: puc ! input as Courant number ; output as flux
	441	!!
	442	INTEGER :: ji, jj, jk ! dummy loop indices
	443	REAL(wp) :: zcoef1, zcoef2, zcoef3 ! local scalars
	444	REAL(wp) :: zc, zcurv, zfho ! - -
	445	!----------------------------------------------------------------------
[3294]	446	!
	447	IF( nn_timing == 1 ) CALL timing_start('quickest')
	448	!
[2528]	449	DO jk = 1, jpkm1
	450	DO jj = 1, jpj
	451	DO ji = 1, jpi
	452	zc = puc(ji,jj,jk) ! Courant number
	453	zcurv = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
	454	zcoef1 = 0.5 * ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
	455	zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
	456	zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
	457	zfho = zcoef1 - zcoef2 - zcoef3 ! phi_f QUICKEST
	458	!
	459	zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
	460	zcoef2 = ABS( zcoef1 )
	461	zcoef3 = ABS( zcurv )
	462	IF( zcoef3 >= zcoef2 ) THEN
	463	zfho = pfc(ji,jj,jk)
	464	ELSE
	465	zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) ) ! phi_REF
	466	IF( zcoef1 >= 0. ) THEN
	467	zfho = MAX( pfc(ji,jj,jk), zfho )
	468	zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) )
	469	ELSE
	470	zfho = MIN( pfc(ji,jj,jk), zfho )
	471	zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) )
	472	ENDIF
	473	ENDIF
	474	puc(ji,jj,jk) = zfho
	475	END DO
	476	END DO
	477	END DO
[1231]	478	!
[3294]	479	IF( nn_timing == 1 ) CALL timing_stop('quickest')
	480	!
[1231]	481	END SUBROUTINE quickest
	482
	483	!!======================================================================
	484	END MODULE traadv_qck

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

Download in other formats: