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

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 2330

Last change on this file since 2330 was 2287, checked in by smasson, 14 years ago
update licence of all NEMO files...
Property svn:keywords set to `Id`
File size: 23.4 KB

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

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

Download in other formats: