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

Context Navigation

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 2082

Last change on this file since 2082 was 2082, checked in by cetlod, 14 years ago
Improve the merge of TRA-TRC, see ticket #717
File size: 23.6 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	!!----------------------------------------------------------------------
	11	!! tra_adv_qck : update the tracer trend with the horizontal advection
	12	!! trends using a 3rd order finite difference scheme
[1559]	13	!! tra_adv_qck_i :
	14	!! tra_adv_qck_j :
	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
[2024]	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
[2034]	34	REAL(wp) :: r1_6 = 1./ 6.
[2024]	35	LOGICAL :: l_trd ! flag to compute trends
[1559]	36
[1231]	37	!! * Substitutions
	38	# include "domzgr_substitute.h90"
	39	# include "vectopt_loop_substitute.h90"
	40	!!----------------------------------------------------------------------
[2034]	41	!! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
	42	!! $Id: traadv_qck.F90 2024 2010-07-29 10:57:35Z cetlod $
[1231]	43	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
	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	!!
[2024]	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
[2082]	88	REAL(wp) , INTENT(in ), DIMENSION(jpk) :: p2dt ! vertical profile of tracer time-step
[2024]	89	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pun, pvn, pwn ! 3 ocean velocity components
[2034]	90	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb, ptn ! before and now tracer fields
	91	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend
[1231]	92	!!----------------------------------------------------------------------
	93
[2082]	94	IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) ) 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
	105	!---------------------------------------------------------------------------
	106
[2082]	107	CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt )
	108	CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt )
[1231]	109
	110	! II. The vertical fluxes are computed with the 2nd order centered scheme
	111	!-------------------------------------------------------------------------
	112	!
[2034]	113	CALL tra_adv_cen2_k( kt, cdtype, pwn, ptn, pta, kjpt )
[1231]	114	!
	115	END SUBROUTINE tra_adv_qck
	116
[2082]	117	SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun, &
	118	& ptb, ptn, pta, kjpt )
[1231]	119	!!----------------------------------------------------------------------
	120	!!
	121	!!----------------------------------------------------------------------
[2024]	122	USE oce , zwx => ua ! use ua as workspace
[2034]	123	!!
[2024]	124	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	125	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	126	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[2082]	127	REAL(wp) , INTENT(in ), DIMENSION(jpk) :: p2dt ! vertical profile of tracer time-step
[2024]	128	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pun ! zonal velocity component
[2034]	129	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb, ptn ! before tracer fields
	130	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend
	131	!!
[2024]	132	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	133	REAL(wp) :: ztra, zbtr ! temporary scalars
	134	REAL(wp) :: zdir, zdx, zdt, zmsk ! temporary scalars
	135	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfu, zfc, zfd
[1231]	136	!----------------------------------------------------------------------
	137
[2024]	138
	139	DO jn = 1, kjpt ! tracer loop
	140	! ! ===========
	141	zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0
	142	zfd(:,:,:) = 0.0 ; zwx(:,:,:) = 0.0
	143	!
	144	DO jk = 1, jpkm1
	145	!
	146	!--- Computation of the ustream and downstream value of the tracer and the mask
	147	DO jj = 2, jpjm1
	148	DO ji = fs_2, fs_jpim1 ! vector opt.
	149	! Upstream in the x-direction for the tracer
[2034]	150	zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn)
[2024]	151	! Downstream in the x-direction for the tracer
[2034]	152	zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn)
[2024]	153	END DO
[1559]	154	END DO
	155	END DO
[1231]	156	!
	157	!--- Lateral boundary conditions
[2024]	158	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
	159
[1231]	160	!
	161	! Horizontal advective fluxes
	162	! ---------------------------
	163	!
[2024]	164	DO jk = 1, jpkm1
	165	DO jj = 2, jpjm1
	166	DO ji = fs_2, fs_jpim1 ! vector opt.
	167	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	168	zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T
	169	END DO
	170	END DO
[1559]	171	END DO
[1231]	172	!
[2024]	173	DO jk = 1, jpkm1
[2082]	174	zdt = p2dt(jk)
[2024]	175	DO jj = 2, jpjm1
	176	DO ji = fs_2, fs_jpim1 ! vector opt.
	177	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	178	zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * fse3u(ji,jj,jk)
	179	zwx(ji,jj,jk) = ABS( pun(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)
[2034]	180	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
	181	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]	182	END DO
	183	END DO
	184	END DO !
	185
	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
[2024]	191	CALL quickest( zfu, zfd, zfc, zwx )
[1231]	192	!
[2024]	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.
	198	ENDDO
[1231]	199	END DO
	200	END DO
[2024]	201	!--- Lateral boundary conditions
	202	CALL lbc_lnk( zfu(:,:,:), 'T', 1. )
[1231]	203	!
[2024]	204	! Tracer flux on the x-direction
	205	DO jk = 1, jpkm1
	206	!
[1231]	207	DO jj = 2, jpjm1
[2024]	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
[2024]	217	!
	218	! Computation of the trend
	219	DO jj = 2, jpjm1
	220	DO ji = fs_2, fs_jpim1 ! vector opt.
	221	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	222	! horizontal advective trends
	223	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
	224	!--- add it to the general tracer trends
[2034]	225	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[2024]	226	END DO
	227	END DO
	228	!
[1231]	229	END DO
[2024]	230	! ! trend diagnostics (contribution of upstream fluxes)
[2034]	231	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jpt_trd_xad, zwx, pun, ptn(:,:,:,jn) )
[2024]	232	!
	233	END DO
	234	!
[1559]	235	END SUBROUTINE tra_adv_qck_i
[1231]	236
[2082]	237	SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn, &
	238	& ptb, ptn, pta, kjpt )
[1231]	239	!!----------------------------------------------------------------------
	240	!!
	241	!!----------------------------------------------------------------------
[2034]	242	!!
[2024]	243	USE oce , zwy => ua ! use ua as workspace
[2034]	244	!!
[2024]	245	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	246	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	247	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[2082]	248	REAL(wp) , INTENT(in ), DIMENSION(jpk) :: p2dt ! vertical profile of tracer time-step
[2024]	249	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pvn ! meridional velocity component
[2034]	250	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb, ptn ! before tracer fields
	251	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend
	252	!!
[2024]	253	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	254	REAL(wp) :: ztra, zbtr ! temporary scalars
	255	REAL(wp) :: zdir, zdx, zdt, zmsk ! temporary scalars
	256	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zfu, zfc, zfd
[1231]	257	!----------------------------------------------------------------------
	258
[2024]	259	DO jn = 1, kjpt ! tracer loop
	260	! ! ===========
	261	zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0
	262	zfd(:,:,:) = 0.0 ; zwy(:,:,:) = 0.0
	263	!
	264	DO jk = 1, jpkm1
	265	!
	266	!--- Computation of the ustream and downstream value of the tracer and the mask
	267	DO jj = 2, jpjm1
	268	DO ji = fs_2, fs_jpim1 ! vector opt.
	269	! Upstream in the x-direction for the tracer
[2034]	270	zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
[2024]	271	! Downstream in the x-direction for the tracer
[2034]	272	zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
[2024]	273	END DO
[1559]	274	END DO
	275	END DO
[1231]	276	!
[2024]	277	!--- Lateral boundary conditions
	278	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
	279
[1231]	280	!
	281	! Horizontal advective fluxes
	282	! ---------------------------
	283	!
[2024]	284	DO jk = 1, jpkm1
	285	DO jj = 2, jpjm1
	286	DO ji = fs_2, fs_jpim1 ! vector opt.
	287	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	288	zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T
	289	END DO
[1559]	290	END DO
	291	END DO
[1231]	292	!
[2024]	293	DO jk = 1, jpkm1
[2082]	294	zdt = p2dt(jk)
[2024]	295	DO jj = 2, jpjm1
	296	DO ji = fs_2, fs_jpim1 ! vector opt.
	297	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	298	zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * fse3v(ji,jj,jk)
	299	zwy(ji,jj,jk) = ABS( pvn(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)
[2034]	300	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
	301	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]	302	END DO
	303	END DO
	304	END DO !
	305
	306	!--- Lateral boundary conditions
	307	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
	308	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
	309
[1231]	310	!--- QUICKEST scheme
[2024]	311	CALL quickest( zfu, zfd, zfc, zwy )
[1231]	312	!
[2024]	313	! Mask at the T-points in the x-direction (mask=0 or mask=1)
	314	DO jk = 1, jpkm1
	315	DO jj = 2, jpjm1
	316	DO ji = fs_2, fs_jpim1 ! vector opt.
	317	zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
	318	ENDDO
[1231]	319	END DO
	320	END DO
[2024]	321	!--- Lateral boundary conditions
	322	CALL lbc_lnk( zfu(:,:,:), 'T', 1. )
	323	!
	324	! Tracer flux on the x-direction
	325	DO jk = 1, jpkm1
	326	!
[1231]	327	DO jj = 2, jpjm1
[2024]	328	DO ji = fs_2, fs_jpim1 ! vector opt.
	329	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	330	!--- If the second ustream point is a land point
	331	!--- the flux is computed by the 1st order UPWIND scheme
	332	zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
	333	zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
	334	zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
[1231]	335	END DO
	336	END DO
[2024]	337	!
	338	! Computation of the trend
	339	DO jj = 2, jpjm1
	340	DO ji = fs_2, fs_jpim1 ! vector opt.
	341	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	342	! horizontal advective trends
	343	ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
	344	!--- add it to the general tracer trends
[2034]	345	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[1231]	346	END DO
	347	END DO
[2024]	348	!
	349	END DO
	350	! ! trend diagnostics (contribution of upstream fluxes)
[2034]	351	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jpt_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
[2024]	352	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	353	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
	354	IF( jn == jp_tem ) pht_adv(:) = ptr_vj( zwy(:,:,:) )
	355	IF( jn == jp_sal ) pst_adv(:) = ptr_vj( zwy(:,:,:) )
[1231]	356	ENDIF
[2024]	357	!
	358	END DO
[1231]	359
[1559]	360	END SUBROUTINE tra_adv_qck_j
[1231]	361
[2034]	362	SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn, &
	363	& ptn, pta, kjpt )
[1231]	364	!!----------------------------------------------------------------------
	365	!!
	366	!!----------------------------------------------------------------------
[2034]	367	!!
[2024]	368	USE oce , zwz => ua ! use ua as workspace
[2034]	369	!!
[2024]	370	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	371	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	372	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	373	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pwn ! vertical velocity component
[2034]	374	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk,kjpt) :: ptn ! now tracer field
	375	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend
	376	!!
[2024]	377	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	378	REAL(wp) :: zbtr , ztra ! temporary scalars
[1559]	379	!!----------------------------------------------------------------------
[2024]	380
[1231]	381	!
[2024]	382	DO jn = 1, kjpt ! tracer loop
	383	! ! ===========
	384	! 1. Bottom value : flux set to zero
	385	zwz(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero
	386	!
	387	! ! Surface value
	388	IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! Variable volume : flux set to zero
[2034]	389	ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn) ! Constant volume : advective flux through the surface
[2024]	390	ENDIF
	391	!
	392	DO jk = 2, jpkm1 ! Interior point: second order centered tracer flux at w-point
	393	DO jj = 2, jpjm1
	394	DO ji = fs_2, fs_jpim1 ! vector opt.
[2034]	395	zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
[2024]	396	END DO
[1231]	397	END DO
	398	END DO
[2024]	399	!
	400	DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==!
	401	DO jj = 2, jpjm1
	402	DO ji = fs_2, fs_jpim1 ! vector opt.
	403	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	404	! k- vertical advective trends
	405	ztra = - zbtr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )
	406	! added to the general tracer trends
[2034]	407	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[2024]	408	END DO
[1231]	409	END DO
	410	END DO
[2024]	411	! ! Save the vertical advective trends for diagnostic
[2034]	412	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jpt_trd_zad, zwz, pwn, ptn(:,:,:,jn) )
[2024]	413	!
[1231]	414	END DO
	415	!
[1559]	416	END SUBROUTINE tra_adv_cen2_k
[1231]	417
	418
[2024]	419	SUBROUTINE quickest( pfu, pfd, pfc, puc )
[1231]	420	!!----------------------------------------------------------------------
	421	!!
[2024]	422	!! ** Purpose : Computation of advective flux with Quickest scheme
	423	!!
	424	!! ** Method :
[1231]	425	!!----------------------------------------------------------------------
[2024]	426	REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pfu ! second upwind point
	427	REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pfd ! first douwning point
	428	REAL(wp), INTENT(in) , DIMENSION(jpi,jpj,jpk) :: pfc ! the central point (or the first upwind point)
	429	REAL(wp), INTENT(inout) , DIMENSION(jpi,jpj,jpk) :: puc ! input as Courant number ; output as flux
	430	!!
	431	INTEGER :: ji, jj, jk ! dummy loop indices
	432	REAL(wp) :: zcoef1, zcoef2, zcoef3 ! temporary scalars
	433	REAL(wp) :: zc, zcurv, zfho !
	434	!----------------------------------------------------------------------
	435
	436	DO jk = 1, jpkm1
	437	DO jj = 1, jpj
	438	DO ji = 1, jpi
	439	zc = puc(ji,jj,jk) ! Courant number
	440	zcurv = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
	441	zcoef1 = 0.5 * ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
	442	zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
	443	zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
	444	zfho = zcoef1 - zcoef2 - zcoef3 ! phi_f QUICKEST
	445	!
	446	zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
	447	zcoef2 = ABS( zcoef1 )
	448	zcoef3 = ABS( zcurv )
	449	IF( zcoef3 >= zcoef2 ) THEN
	450	zfho = pfc(ji,jj,jk)
	451	ELSE
	452	zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) ) ! phi_REF
	453	IF( zcoef1 >= 0. ) THEN
	454	zfho = MAX( pfc(ji,jj,jk), zfho )
	455	zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) )
	456	ELSE
	457	zfho = MIN( pfc(ji,jj,jk), zfho )
	458	zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) )
	459	ENDIF
	460	ENDIF
	461	puc(ji,jj,jk) = zfho
	462	ENDDO
	463	ENDDO
	464	ENDDO
[1231]	465	!
	466	END SUBROUTINE quickest
	467
	468	!!======================================================================
	469	END MODULE traadv_qck

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

Download in other formats: