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

Context Navigation

source: trunk/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 1839

Last change on this file since 1839 was 1559, checked in by ctlod, 15 years ago
only cosmetic changes
File size: 23.1 KB

Rev	Line
[1231]	1	MODULE traadv_qck
	2	!!==============================================================================
	3	!! * MODULE traadv_qck *
	4	!! Ocean active tracers: horizontal & vertical advective trend
	5	!!==============================================================================
[1559]	6	!! History : 3.0 ! 2008-07 (G. Reffray) Original code
[1231]	7	!!----------------------------------------------------------------------
	8
	9	!!----------------------------------------------------------------------
	10	!! tra_adv_qck : update the tracer trend with the horizontal advection
	11	!! trends using a 3rd order finite difference scheme
[1559]	12	!! tra_adv_qck_i :
	13	!! tra_adv_qck_j :
	14	!! tra_adv_cen2_k : 2nd centered scheme for the vertical advection
[1231]	15	!!----------------------------------------------------------------------
	16	USE oce ! ocean dynamics and active tracers
	17	USE dom_oce ! ocean space and time domain
	18	USE trdmod ! ocean active tracers trends
	19	USE trdmod_oce ! ocean variables trends
	20	USE trabbl ! advective term in the BBL
	21	USE lib_mpp ! distribued memory computing
	22	USE lbclnk ! ocean lateral boundary condition (or mpp link)
	23	USE dynspg_oce ! surface pressure gradient variables
	24	USE in_out_manager ! I/O manager
	25	USE diaptr ! poleward transport diagnostics
	26	USE prtctl ! Print control
	27
	28	IMPLICIT NONE
	29	PRIVATE
	30
[1559]	31	PUBLIC tra_adv_qck ! routine called by step.F90
[1231]	32
[1559]	33	REAL(wp), DIMENSION(jpi,jpj) :: btr2
	34	REAL(wp) :: r1_6
	35
[1231]	36	!! * Substitutions
	37	# include "domzgr_substitute.h90"
	38	# include "vectopt_loop_substitute.h90"
	39	!!----------------------------------------------------------------------
[1559]	40	!! NEMO/OPA 3.2 , LOCEAN-IPSL (2009)
[1231]	41	!! $Id$
	42	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
	43	!!----------------------------------------------------------------------
	44
	45	CONTAINS
	46
	47	SUBROUTINE tra_adv_qck( kt, pun, pvn, pwn )
	48	!!----------------------------------------------------------------------
	49	!! * ROUTINE tra_adv_qck *
	50	!!
	51	!! ** Purpose : Compute the now trend due to the advection of tracers
	52	!! and add it to the general trend of passive tracer equations.
	53	!!
	54	!! ** Method : The advection is evaluated by a third order scheme
[1559]	55	!! For a positive velocity u : u(i)>0
	56	!! \|--FU--\|--FC--\|--FD--\|------\|
	57	!! i-1 i i+1 i+2
[1231]	58	!!
[1559]	59	!! For a negative velocity u : u(i)<0
	60	!! \|------\|--FD--\|--FC--\|--FU--\|
	61	!! i-1 i i+1 i+2
	62	!! where FU is the second upwind point
	63	!! FD is the first douwning point
	64	!! FC is the central point (or the first upwind point)
[1231]	65	!!
[1559]	66	!! Flux(i) = u(i) * { 0.5(FC+FD) -0.5C(i)(FD-FC) -((1-C(i))/6)(FU+FD-2FC) }
	67	!! with C(i)=\|u(i)\|dx(i)/dt (=Courant number)
[1231]	68	!!
	69	!! dt = 2*rdtra and the scalar values are tb and sb
	70	!!
	71	!! On the vertical, the simple centered scheme used tn and sn
	72	!!
[1559]	73	!! The fluxes are bounded by the ULTIMATE limiter to
	74	!! guarantee the monotonicity of the solution and to
[1231]	75	!! prevent the appearance of spurious numerical oscillations
	76	!!
	77	!! ** Action : - update (ta,sa) with the now advective tracer trends
	78	!! - save the trends ('key_trdtra')
	79	!!
	80	!! ** Reference : Leonard (1979, 1991)
	81	!!----------------------------------------------------------------------
	82	USE oce, ONLY : ztrdt => ua ! use ua as workspace
	83	USE oce, ONLY : ztrds => va ! use va as workspace
[1559]	84	!!
[1231]	85	INTEGER , INTENT(in) :: kt ! ocean time-step index
	86	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pun ! effective ocean velocity, u_component
	87	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pvn ! effective ocean velocity, v_component
	88	REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) :: pwn ! effective ocean velocity, w_component
	89	!!
	90	INTEGER :: ji, jj, jk ! dummy loop indices
	91	REAL(wp) :: z_hdivn_x, z_hdivn_y, z_hdivn ! temporary scalars
	92	REAL(wp) :: zbtr, z2 ! " "
	93	!!----------------------------------------------------------------------
	94
	95	IF( kt == nit000 ) THEN
	96	IF(lwp) WRITE(numout,*)
	97	IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme'
	98	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
	99	IF(lwp) WRITE(numout,*)
	100	btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
[1559]	101	r1_6 = 1. / 6.
[1231]	102	ENDIF
	103
	104	IF( neuler == 0 .AND. kt == nit000 ) THEN ; z2 = 1.
	105	ELSE ; z2 = 2.
	106	ENDIF
	107
	108	! I. The horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
	109	!---------------------------------------------------------------------------
	110
[1559]	111	CALL tra_adv_qck_i( pun, tb, tn, ta, ztrdt, z2)
	112	CALL tra_adv_qck_i( pun, sb, sn, sa, ztrds, z2)
[1231]	113
	114	IF( l_trdtra ) CALL trd_mod(ztrdt, ztrds, jptra_trd_xad, 'TRA', kt)
	115
[1559]	116	CALL tra_adv_qck_j( kt, pvn, tb, tn, ta, ztrdt, pht_adv, z2)
	117	CALL tra_adv_qck_j( kt, pvn, sb, sn, sa, ztrds, pst_adv, z2)
[1231]	118
	119	IF( l_trdtra ) THEN
	120	CALL trd_mod(ztrdt, ztrds, jptra_trd_yad, 'TRA', kt)
	121	!
[1559]	122	ztrdt(:,:,:) = ta(:,:,:) ! Save the horizontal up-to-date ta/sa trends
[1231]	123	ztrds(:,:,:) = sa(:,:,:)
	124	END IF
	125
	126	IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' qck had - Ta: ', mask1=tmask, &
	127	& tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
	128
	129	! II. The vertical fluxes are computed with the 2nd order centered scheme
	130	!-------------------------------------------------------------------------
	131	!
[1559]	132	CALL tra_adv_cen2_k( pwn, tn, ta )
	133	CALL tra_adv_cen2_k( pwn, sn, sa )
[1231]	134	!
	135	!Save the vertical advective trends for diagnostic
	136	! ----------------------------------------------------
	137	IF( l_trdtra ) THEN
	138	! Recompute the vertical advection zta & zsa trends computed
	139	! at the step 2. above in making the difference between the new
	140	! trends and the previous one: ta()/sa - ztrdt()/ztrds() and substract
	141	! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
	142
	143	DO jk = 1, jpkm1
	144	DO jj = 2, jpjm1
	145	DO ji = fs_2, fs_jpim1 ! vector opt.
	146	#if defined key_zco
	147	zbtr = btr2(ji,jj)
	148	z_hdivn_x = e2u(ji,jj)pun(ji,jj,jk) - e2u(ji-1,jj)pun(ji-1,jj,jk)
	149	z_hdivn_y = e1v(ji,jj)pvn(ji,jj,jk) - e1v(ji,jj-1)pvn(ji,jj-1,jk)
	150	#else
	151	zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
	152	z_hdivn_x = e2u(ji,jj)fse3u(ji,jj,jk)pun(ji,jj,jk) - e2u(ji-1,jj)fse3u(ji-1,jj,jk)pun(ji-1,jj,jk)
	153	z_hdivn_y = e1v(ji,jj)fse3v(ji,jj,jk)pvn(ji,jj,jk) - e1v(ji,jj-1)fse3v(ji,jj-1,jk)pvn(ji,jj-1,jk)
	154	#endif
	155	z_hdivn = (z_hdivn_x + z_hdivn_y) * zbtr
	156	ztrdt(ji,jj,jk) = ta(ji,jj,jk) - ztrdt(ji,jj,jk) - tn(ji,jj,jk) * z_hdivn
	157	ztrds(ji,jj,jk) = sa(ji,jj,jk) - ztrds(ji,jj,jk) - sn(ji,jj,jk) * z_hdivn
	158	END DO
	159	END DO
	160	END DO
	161	CALL trd_mod(ztrdt, ztrds, jptra_trd_zad, 'TRA', kt)
	162	ENDIF
	163
	164	IF(ln_ctl) CALL prt_ctl( tab3d_1=ta, clinfo1=' qck zad - Ta: ', mask1=tmask, &
	165	& tab3d_2=sa, clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
	166	!
	167	END SUBROUTINE tra_adv_qck
	168
	169
[1559]	170	SUBROUTINE tra_adv_qck_i ( pun, tra, tran, traa, ztrdtra, z2 )
[1231]	171	!!----------------------------------------------------------------------
	172	!!
	173	!!----------------------------------------------------------------------
	174	REAL, INTENT(in) :: z2
	175	REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pun, tra, tran ! horizontal effective velocity
	176	REAL(wp), INTENT(out) , DIMENSION(jpi,jpj,jpk) :: ztrdtra
	177	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: traa
	178	!
	179	INTEGER :: ji, jj, jk
	180	REAL(wp) :: za, zbtr, dir, dx, dt ! temporary scalars
	181	REAL(wp) :: z_hdivn_x
	182	REAL(wp), DIMENSION(jpi,jpj) :: zmask, zupst, zdwst, zc_cfl
	183	REAL(wp), DIMENSION(jpi,jpj) :: zfu, zfc, zfd, zfho, zmskl, zsc_e
	184	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux
	185	!----------------------------------------------------------------------
	186
	187	zfu (:,jpj) = 0.e0 ; zfc (:,jpj) = 0.e0
	188	zfd (:,jpj) = 0.e0 ; zc_cfl(:,jpj) = 0.e0
	189	zsc_e (:,jpj) = 0.e0 ; zmskl (:,jpj) = 0.e0
	190	zfho (:,jpj) = 0.e0
	191	! ===============
	192	DO jk = 1, jpkm1 ! Horizontal slab
	193	! ! ===============
	194	!--- Computation of the ustream and downstream value of the tracer and the mask
	195	DO jj = 2, jpjm1
	196	DO ji = 2, fs_jpim1 ! vector opt.
	197	! Upstream in the x-direction for the tracer
	198	zupst(ji,jj)=tra(ji-1,jj,jk)
	199	! Downstream in the x-direction for the tracer
	200	zdwst(ji,jj)=tra(ji+1,jj,jk)
	201	! Mask at the T-points in the x-direction (mask=0 or mask=1)
	202	zmask(ji,jj)=tmask(ji-1,jj,jk)+tmask(ji,jj,jk)+tmask(ji+1,jj,jk)-2
[1559]	203	END DO
	204	END DO
[1231]	205	!
	206	!--- Lateral boundary conditions
	207	CALL lbc_lnk( zupst(:,:), 'T', 1. )
	208	CALL lbc_lnk( zdwst(:,:), 'T', 1. )
	209	CALL lbc_lnk( zmask(:,:), 'T', 1. )
	210	!
	211	! Horizontal advective fluxes
	212	! ---------------------------
	213	!
	214	dt = z2 * rdttra(jk)
	215	!--- tracer flux at u-points
	216	DO jj = 1, jpjm1
	217	DO ji = 1, jpi
	218	#if defined key_zco
	219	zsc_e(ji,jj) = e2u(ji,jj)
	220	#else
	221	zsc_e(ji,jj) = e2u(ji,jj) * fse3u(ji,jj,jk)
	222	#endif
	223	dir = 0.5 + sign(0.5,pun(ji,jj,jk)) ! if pun>0 : dir = 1 otherwise dir = 0
	224	dx = dir * e1t(ji,jj) + (1-dir)* e1t(ji+1,jj)
	225	zc_cfl (ji,jj) = ABS(pun(ji,jj,jk))*dt/dx ! (0<zc_cfl<1 : Courant number on x-direction)
	226
	227	zfu(ji,jj) = dirzupst(ji ,jj )+(1-dir)zdwst(ji+1,jj ) ! FU in the x-direction for T
	228	zfc(ji,jj) = dirtra (ji ,jj,jk)+(1-dir)tra (ji+1,jj,jk) ! FC in the x-direction for T
	229	zfd(ji,jj) = dirtra (ji+1,jj,jk)+(1-dir)tra (ji ,jj,jk) ! FD in the x-direction for T
	230	zmskl(ji,jj) = dirzmask(ji ,jj) +(1-dir)zmask(ji+1,jj)
[1559]	231	END DO
	232	END DO
[1231]	233	!
	234	!--- QUICKEST scheme
	235	! Tracer flux on the x-direction
	236	CALL quickest(zfu,zfd,zfc,zfho,zc_cfl)
	237	!--- If the second ustream point is a land point
	238	!--- the flux is computed by the 1st order UPWIND scheme
	239	zfho(:,:) = zmskl(:,:)zfho(:,:) + (1.-zmskl(:,:))zfc(:,:)
	240	!--- Computation of fluxes
	241	zflux(:,:,jk) = zsc_e(:,:)pun(:,:,jk)zfho(:,:)
	242	!
	243	!--- Tracer flux divergence at t-point added to the general trend
	244	DO jj = 2, jpjm1
	245	DO ji = fs_2, fs_jpim1 ! vector opt.
	246	!--- horizontal advective trends
	247	#if defined key_zco
	248	zbtr = btr2(ji,jj)
	249	#else
	250	zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
	251	#endif
	252	za = - zbtr * ( zflux(ji,jj,jk) - zflux(ji-1,jj,jk) )
	253	!--- add it to the general tracer trends
	254	traa(ji,jj,jk) = traa(ji,jj,jk) + za
	255	END DO
	256	END DO
	257	! ! ===============
	258	END DO ! End of slab
	259	! ! ===============
	260	!
	261	! Save the horizontal advective trends for diagnostic
	262	! -----------------------------------------------------
	263	IF( l_trdtra ) THEN
	264	! T/S ZONAL advection trends
	265	ztrdtra(:,:,:) = 0.e0
	266	!
	267	DO jk = 1, jpkm1
	268	DO jj = 2, jpjm1
	269	DO ji = fs_2, fs_jpim1 ! vector opt.
	270	!-- Compute zonal divergence by splitting hdivn (see divcur.F90)
	271	! N.B. This computation is not valid along OBCs (if any)
	272	#if defined key_zco
	273	zbtr = btr2(ji,jj)
	274	z_hdivn_x = ( e2u(ji ,jj) * pun(ji ,jj,jk) &
	275	& - e2u(ji-1,jj) * pun(ji-1,jj,jk) ) * zbtr
	276	#else
	277	zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
	278	z_hdivn_x = ( e2u(ji ,jj) * fse3u(ji ,jj,jk) * pun(ji ,jj,jk) &
	279	& - e2u(ji-1,jj) * fse3u(ji-1,jj,jk) * pun(ji-1,jj,jk) ) * zbtr
	280	#endif
	281	ztrdtra(ji,jj,jk) = - zbtr * ( zflux(ji,jj,jk) - zflux(ji-1,jj,jk) ) + tran(ji,jj,jk) * z_hdivn_x
	282	END DO
	283	END DO
	284	END DO
	285	END IF
	286
[1559]	287	END SUBROUTINE tra_adv_qck_i
[1231]	288
	289
[1559]	290	SUBROUTINE tra_adv_qck_j ( kt, pvn, tra, tran, traa, ztrdtra, trd_adv, z2 )
[1231]	291	!!----------------------------------------------------------------------
	292	!!
	293	!!----------------------------------------------------------------------
	294	INTEGER, INTENT(in) :: kt ! ocean time-step index
	295	REAL, INTENT(in) :: z2
	296	REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pvn, tra, tran ! horizontal effective velocity
	297	REAL(wp), INTENT(out) , DIMENSION(jpj) :: trd_adv
	298	REAL(wp), INTENT(out) , DIMENSION(jpi,jpj,jpk) :: ztrdtra
	299	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: traa
[1559]	300	!!
[1231]	301	INTEGER :: ji, jj, jk
	302	REAL(wp) :: za, zbtr, dir, dx, dt ! temporary scalars
	303	REAL(wp) :: z_hdivn_y
	304	REAL(wp), DIMENSION(jpi,jpj) :: zmask, zupst, zdwst, zc_cfl
	305	REAL(wp), DIMENSION(jpi,jpj) :: zfu, zfc, zfd, zfho, zmskl, zsc_e
	306	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux
	307	!----------------------------------------------------------------------
	308	! II. Part 2 : y-direction
	309	!----------------------------------------------------------------------
	310
	311	zfu (:,jpj) = 0.e0 ; zfc (:,jpj) = 0.e0
	312	zfd (:,jpj) = 0.e0 ; zc_cfl(:,jpj) = 0.e0
	313	zsc_e (:,jpj) = 0.e0 ; zmskl (:,jpj) = 0.e0
	314	zfho (:,jpj) = 0.e0
	315
	316	! ===============
	317	DO jk = 1, jpkm1 ! Horizontal slab
	318	! ! ===============
	319	!--- Computation of the ustream and downstream value of the tracer and the mask
	320	DO jj = 2, jpjm1
	321	DO ji = 2, fs_jpim1 ! vector opt.
	322	! Upstream in the x-direction for the tracer
	323	zupst(ji,jj)=tra(ji,jj-1,jk)
	324	! Downstream in the x-direction for the tracer
	325	zdwst(ji,jj)=tra(ji,jj+1,jk)
	326	! Mask at the T-points in the x-direction (mask=0 or mask=1)
	327	zmask(ji,jj)=tmask(ji,jj-1,jk)+tmask(ji,jj,jk)+tmask(ji,jj+1,jk)-2
[1559]	328	END DO
	329	END DO
[1231]	330	!
	331	!--- Lateral boundary conditions
	332	CALL lbc_lnk( zupst(:,:), 'T', 1. )
	333	CALL lbc_lnk( zdwst(:,:), 'T', 1. )
	334	CALL lbc_lnk( zmask(:,:), 'T', 1. )
	335	!
	336	! Horizontal advective fluxes
	337	! ---------------------------
	338	!
	339	dt = z2 * rdttra(jk)
	340	!--- tracer flux at v-points
	341	DO jj = 1, jpjm1
	342	DO ji = 1, jpi
	343	#if defined key_zco
	344	zsc_e(ji,jj) = e1v(ji,jj)
	345	#else
	346	zsc_e(ji,jj) = e1v(ji,jj) * fse3v(ji,jj,jk)
	347	#endif
	348	dir = 0.5 + sign(0.5,pvn(ji,jj,jk)) ! if pvn>0 : dir = 1 otherwise dir = 0
	349	dx = dir * e2t(ji,jj) + (1-dir)* e2t(ji,jj+1)
	350	zc_cfl(ji,jj) = ABS(pvn(ji,jj,jk))*dt/dx ! (0<zc_cfl<1 : Courant number on y-direction)
	351
	352	zfu(ji,jj) = dirzupst(ji,jj )+(1-dir)zdwst(ji,jj+1 ) ! FU in the y-direction for T
	353	zfc(ji,jj) = dirtra (ji,jj ,jk)+(1-dir)tra (ji,jj+1,jk) ! FC in the y-direction for T
	354	zfd(ji,jj) = dirtra (ji,jj+1,jk)+(1-dir)tra (ji,jj ,jk) ! FD in the y-direction for T
	355	zmskl(ji,jj) = dirzmask(ji,jj )+(1-dir)zmask(ji,jj+1)
[1559]	356	END DO
	357	END DO
[1231]	358	!
	359	!--- QUICKEST scheme
	360	! Tracer flux on the y-direction
	361	CALL quickest(zfu,zfd,zfc,zfho,zc_cfl)
	362	!--- If the second ustream point is a land point
	363	!--- the flux is computed by the 1st order UPWIND scheme
	364	zfho(:,:) = zmskl(:,:)zfho(:,:) + (1.-zmskl(:,:))zfc(:,:)
	365	!--- Computation of fluxes
	366	zflux(:,:,jk) = zsc_e(:,:)pvn(:,:,jk)zfho(:,:)
	367	!
	368	!--- Tracer flux divergence at t-point added to the general trend
	369	DO jj = 2, jpjm1
	370	DO ji = fs_2, fs_jpim1 ! vector opt.
	371	!--- horizontal advective trends
	372	#if defined key_zco
	373	zbtr = btr2(ji,jj)
	374	#else
	375	zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
	376	#endif
	377	za = - zbtr * ( zflux(ji,jj,jk) - zflux(ji,jj-1,jk) )
	378	!--- add it to the general tracer trends
	379	traa(ji,jj,jk) = traa(ji,jj,jk) + za
	380	END DO
	381	END DO
	382	! ! ===============
	383	END DO ! End of slab
	384	! ! ===============
	385	!
	386	! Save the horizontal advective trends for diagnostic
	387	! -----------------------------------------------------
	388	IF( l_trdtra ) THEN
	389	! T/S MERIDIONAL advection trends
	390	DO jk = 1, jpkm1
	391	DO jj = 2, jpjm1
	392	DO ji = fs_2, fs_jpim1 ! vector opt.
	393	!-- Compute merid. divergence by splitting hdivn (see divcur.F90)
	394	! N.B. This computation is not valid along OBCs (if any)
	395	#if defined key_zco
	396	zbtr = btr2(ji,jj)
	397	z_hdivn_y = ( e1v(ji,jj ) * pvn(ji,jj ,jk) &
	398	& - e1v(ji,jj-1) * pvn(ji,jj-1,jk) ) * zbtr
	399	#else
	400	zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
	401	z_hdivn_y = ( e1v(ji, jj) * fse3v(ji,jj ,jk) * pvn(ji,jj ,jk) &
	402	& - e1v(ji,jj-1) * fse3v(ji,jj-1,jk) * pvn(ji,jj-1,jk) ) * zbtr
	403	#endif
	404	ztrdtra(ji,jj,jk) = - zbtr * ( zflux(ji,jj,jk) - zflux(ji,jj-1,jk) ) + tran(ji,jj,jk) * z_hdivn_y
	405	END DO
	406	END DO
	407	END DO
	408	END IF
	409
	410	! "zonal" mean advective heat and salt transport
	411	! ----------------------------------------------
	412
	413	IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
	414	IF( lk_zco ) THEN
	415	DO jk = 1, jpkm1
	416	DO jj = 2, jpjm1
	417	DO ji = fs_2, fs_jpim1 ! vector opt.
	418	zflux(ji,jj,jk) = zflux(ji,jj,jk) * fse3v(ji,jj,jk)
	419	END DO
	420	END DO
	421	END DO
	422	ENDIF
	423	trd_adv(:) = ptr_vj( zflux(:,:,:) )
	424	ENDIF
	425
[1559]	426	END SUBROUTINE tra_adv_qck_j
[1231]	427
	428
[1559]	429	SUBROUTINE tra_adv_cen2_k ( pwn, ptn, pta )
[1231]	430	!!----------------------------------------------------------------------
	431	!!
	432	!!----------------------------------------------------------------------
[1559]	433	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwn ! vertical effective velocity
	434	REAL(wp), INTENT(in ), DIMENSION(jpi,jpj,jpk) :: ptn ! now tracer
	435	REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: pta ! tracer general trend
	436	!!
	437	INTEGER :: ji, jj, jk ! dummy loop indices
	438	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zflux ! 3D workspace
	439	!!----------------------------------------------------------------------
[1231]	440	!
[1559]	441	! !== Vertical advective fluxes ==!
	442	zflux(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero
[1231]	443	!
[1559]	444	! ! Surface value
	445	IF( lk_vvl ) THEN ; zflux(:,:, 1 ) = 0.e0 ! Variable volume : flux set to zero
	446	ELSE ; zflux(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1) ! Constant volume : advective flux through the surface
[1231]	447	ENDIF
	448	!
[1559]	449	DO jk = 2, jpkm1 ! Interior point: second order centered tracer flux at w-point
[1231]	450	DO jj = 2, jpjm1
	451	DO ji = fs_2, fs_jpim1 ! vector opt.
[1559]	452	zflux(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1) + ptn(ji,jj,jk) )
[1231]	453	END DO
	454	END DO
	455	END DO
	456	!
[1559]	457	DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==!
[1231]	458	DO jj = 2, jpjm1
	459	DO ji = fs_2, fs_jpim1 ! vector opt.
[1559]	460	pta(ji,jj,jk) = pta(ji,jj,jk) - ( zflux(ji,jj,jk) - zflux(ji,jj,jk+1) ) &
	461	& / fse3t(ji,jj,jk)
[1231]	462	END DO
	463	END DO
	464	END DO
	465	!
[1559]	466	END SUBROUTINE tra_adv_cen2_k
[1231]	467
	468
[1559]	469	SUBROUTINE quickest( fu, fd, fc, fho, fc_cfl )
[1231]	470	!!----------------------------------------------------------------------
	471	!!
	472	!!----------------------------------------------------------------------
	473	REAL(wp), INTENT(in) , DIMENSION(jpi,jpj) :: fu, fd, fc, fc_cfl
	474	REAL(wp), INTENT(out) , DIMENSION(jpi,jpj) :: fho
	475	REAL(wp) , DIMENSION(jpi,jpj) :: zcurv, zcoef1, zcoef2, zcoef3 ! temporary scalars
	476	!
	477	zcurv (:,:) = fd(:,:) + fu(:,:) - 2.*fc(:,:)
	478	zcoef1(:,:) = 0.5*( fc(:,:) + fd(:,:) )
	479	zcoef2(:,:) = 0.5fc_cfl(:,:)( fd(:,:) - fc(:,:) )
[1559]	480	zcoef3(:,:) = ( ( 1. - ( fc_cfl(:,:)fc_cfl(:,:) ) )r1_6 )*zcurv(:,:)
[1231]	481	fho (:,:) = zcoef1(:,:) - zcoef2(:,:) - zcoef3(:,:) ! phi_f QUICKEST
	482	!
	483	zcoef1(:,:) = fd(:,:) - fu(:,:) ! DEL
	484	zcoef2(:,:) = ABS( zcoef1(:,:) ) ! ABS(DEL)
	485	zcoef3(:,:) = ABS( zcurv(:,:) ) ! ABS(CURV)
	486	!
	487	WHERE ( zcoef3(:,:) >= zcoef2(:,:) )
	488	fho(:,:) = fc(:,:)
	489	ELSEWHERE
	490	zcoef3(:,:) = fu(:,:) + ( ( fc(:,:) - fu(:,:) )/MAX(fc_cfl(:,:),1.e-9) ) ! phi_REF
	491	WHERE ( zcoef1(:,:) >= 0.e0 )
	492	fho(:,:) = MAX(fc(:,:),fho(:,:))
	493	fho(:,:) = MIN(fho(:,:),MIN(zcoef3(:,:),fd(:,:)))
	494	ELSEWHERE
	495	fho(:,:) = MIN(fc(:,:),fho(:,:))
	496	fho(:,:) = MAX(fho(:,:),MAX(zcoef3(:,:),fd(:,:)))
	497	ENDWHERE
	498	ENDWHERE
	499	!
	500	END SUBROUTINE quickest
	501
	502	!!======================================================================
	503	END MODULE traadv_qck

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

Download in other formats: