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

Context Navigation

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 6004

Last change on this file since 6004 was 6004, checked in by gm, 8 years ago
#1613: vvl by default, step III: Merge with the trunk (free surface simplification) (see wiki)
Property svn:keywords set to `Id`
File size: 24.0 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
[4990]	19	USE trc_oce ! share passive tracers/Ocean variables
	20	USE trd_oce ! trends: ocean variables
	21	USE trdtra ! trends manager: tracers
	22	USE diaptr ! poleward transport diagnostics
	23	!
[1231]	24	USE lib_mpp ! distribued memory computing
	25	USE lbclnk ! ocean lateral boundary condition (or mpp link)
	26	USE in_out_manager ! I/O manager
[3294]	27	USE wrk_nemo ! Memory Allocation
	28	USE timing ! Timing
[3625]	29	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[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 "vectopt_loop_substitute.h90"
	41	!!----------------------------------------------------------------------
[2528]	42	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1231]	43	!! $Id$
[2528]	44	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[1231]	45	!!----------------------------------------------------------------------
	46	CONTAINS
	47
[3294]	48	SUBROUTINE tra_adv_qck ( kt, kit000, cdtype, p2dt, pun, pvn, pwn, &
[2528]	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	!!
[2528]	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	!!
[5883]	79	!! ** Action : - update pta with the now advective tracer trends
	80	!! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
	81	!! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
[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	!!----------------------------------------------------------------------
[3294]	94	!
	95	IF( nn_timing == 1 ) CALL timing_start('tra_adv_qck')
	96	!
	97	IF( kt == kit000 ) THEN
[1231]	98	IF(lwp) WRITE(numout,*)
[2528]	99	IF(lwp) WRITE(numout,*) 'tra_adv_qck : 3rd order quickest advection scheme on ', cdtype
[1231]	100	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
	101	IF(lwp) WRITE(numout,*)
	102	ENDIF
[5836]	103	!
[4990]	104	l_trd = .FALSE.
	105	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
[4499]	106	!
[5883]	107	! ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
[2528]	108	CALL tra_adv_qck_i( kt, cdtype, p2dt, pun, ptb, ptn, pta, kjpt )
	109	CALL tra_adv_qck_j( kt, cdtype, p2dt, pvn, ptb, ptn, pta, kjpt )
[1231]	110
[5883]	111	! ! vertical fluxes are computed with the 2nd order centered scheme
[2528]	112	CALL tra_adv_cen2_k( kt, cdtype, pwn, ptn, pta, kjpt )
[1231]	113	!
[3294]	114	IF( nn_timing == 1 ) CALL timing_stop('tra_adv_qck')
	115	!
[1231]	116	END SUBROUTINE tra_adv_qck
	117
	118
[2528]	119	SUBROUTINE tra_adv_qck_i( kt, cdtype, p2dt, pun, &
	120	& ptb, ptn, pta, kjpt )
[1231]	121	!!----------------------------------------------------------------------
	122	!!
	123	!!----------------------------------------------------------------------
[2715]	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
	127	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
	128	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun ! i-velocity components
	129	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
	130	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[2528]	131	!!
[5836]	132	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	133	REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk ! local scalars
	134	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwx, zfu, zfc, zfd
[1231]	135	!----------------------------------------------------------------------
[2715]	136	!
[4990]	137	CALL wrk_alloc( jpi, jpj, jpk, zwx, zfu, zfc, zfd )
[2528]	138	! ! ===========
	139	DO jn = 1, kjpt ! tracer loop
	140	! ! ===========
[5836]	141	zfu(:,:,:) = 0._wp ; zfc(:,:,:) = 0._wp
	142	zfd(:,:,:) = 0._wp ; zwx(:,:,:) = 0._wp
	143	!
	144	!!gm why not using a SHIFT instruction...
	145	DO jk = 1, jpkm1 !--- Computation of the ustream and downstream value of the tracer and the mask
[2528]	146	DO jj = 2, jpjm1
	147	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	148	zfc(ji,jj,jk) = ptb(ji-1,jj,jk,jn) ! Upstream in the x-direction for the tracer
	149	zfd(ji,jj,jk) = ptb(ji+1,jj,jk,jn) ! Downstream in the x-direction for the tracer
[2528]	150	END DO
[1559]	151	END DO
	152	END DO
[2528]	153	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) ! Lateral boundary conditions
	154
[1231]	155	!
	156	! Horizontal advective fluxes
	157	! ---------------------------
[2528]	158	DO jk = 1, jpkm1
	159	DO jj = 2, jpjm1
	160	DO ji = fs_2, fs_jpim1 ! vector opt.
	161	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	162	zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T
	163	END DO
	164	END DO
[1559]	165	END DO
[1231]	166	!
[2528]	167	DO jk = 1, jpkm1
	168	zdt = p2dt(jk)
	169	DO jj = 2, jpjm1
	170	DO ji = fs_2, fs_jpim1 ! vector opt.
	171	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
[5845]	172	zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u_n(ji,jj,jk)
[2528]	173	zwx(ji,jj,jk) = ABS( pun(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)
	174	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
	175	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
	176	END DO
	177	END DO
	178	END DO
	179	!--- Lateral boundary conditions
	180	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
	181	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zwx(:,:,:), 'T', 1. )
	182
[1231]	183	!--- QUICKEST scheme
[2528]	184	CALL quickest( zfu, zfd, zfc, zwx )
[1231]	185	!
[2528]	186	! Mask at the T-points in the x-direction (mask=0 or mask=1)
	187	DO jk = 1, jpkm1
	188	DO jj = 2, jpjm1
	189	DO ji = fs_2, fs_jpim1 ! vector opt.
	190	zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
[2715]	191	END DO
[1231]	192	END DO
	193	END DO
[2528]	194	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ! Lateral boundary conditions
	195
[1231]	196	!
[2528]	197	! Tracer flux on the x-direction
	198	DO jk = 1, jpkm1
	199	!
[1231]	200	DO jj = 2, jpjm1
[2528]	201	DO ji = fs_2, fs_jpim1 ! vector opt.
	202	zdir = 0.5 + SIGN( 0.5, pun(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	203	!--- If the second ustream point is a land point
	204	!--- the flux is computed by the 1st order UPWIND scheme
	205	zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
	206	zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
	207	zwx(ji,jj,jk) = zwx(ji,jj,jk) * pun(ji,jj,jk)
[1231]	208	END DO
	209	END DO
[3300]	210	END DO
	211	!
	212	CALL lbc_lnk( zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
	213	!
	214	! Computation of the trend
	215	DO jk = 1, jpkm1
[2528]	216	DO jj = 2, jpjm1
	217	DO ji = fs_2, fs_jpim1 ! vector opt.
[6004]	218	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
[2528]	219	! horizontal advective trends
	220	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
	221	!--- add it to the general tracer trends
	222	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
	223	END DO
	224	END DO
[1231]	225	END DO
[5883]	226	! ! trend diagnostics
[4990]	227	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptn(:,:,:,jn) )
[2528]	228	!
	229	END DO
	230	!
[4990]	231	CALL wrk_dealloc( jpi, jpj, jpk, zwx, zfu, zfc, zfd )
[2715]	232	!
[1559]	233	END SUBROUTINE tra_adv_qck_i
[1231]	234
	235
[2528]	236	SUBROUTINE tra_adv_qck_j( kt, cdtype, p2dt, pvn, &
	237	& ptb, ptn, pta, kjpt )
[1231]	238	!!----------------------------------------------------------------------
	239	!!
	240	!!----------------------------------------------------------------------
[2715]	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
[1559]	248	!!
[2715]	249	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	250	REAL(wp) :: ztra, zbtr, zdir, zdx, zdt, zmsk ! local scalars
[4990]	251	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwy, zfu, zfc, zfd
[1231]	252	!----------------------------------------------------------------------
[2715]	253	!
[4990]	254	CALL wrk_alloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd )
[3294]	255	!
[2528]	256	! ! ===========
	257	DO jn = 1, kjpt ! tracer loop
	258	! ! ===========
	259	zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0
	260	zfd(:,:,:) = 0.0 ; zwy(:,:,:) = 0.0
	261	!
	262	DO jk = 1, jpkm1
	263	!
	264	!--- Computation of the ustream and downstream value of the tracer and the mask
	265	DO jj = 2, jpjm1
	266	DO ji = fs_2, fs_jpim1 ! vector opt.
	267	! Upstream in the x-direction for the tracer
	268	zfc(ji,jj,jk) = ptb(ji,jj-1,jk,jn)
	269	! Downstream in the x-direction for the tracer
	270	zfd(ji,jj,jk) = ptb(ji,jj+1,jk,jn)
	271	END DO
[1559]	272	END DO
	273	END DO
[2528]	274	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. ) ! Lateral boundary conditions
	275
	276
[1231]	277	!
	278	! Horizontal advective fluxes
	279	! ---------------------------
	280	!
[2528]	281	DO jk = 1, jpkm1
	282	DO jj = 2, jpjm1
	283	DO ji = fs_2, fs_jpim1 ! vector opt.
	284	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	285	zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T
	286	END DO
[1559]	287	END DO
	288	END DO
[1231]	289	!
[2528]	290	DO jk = 1, jpkm1
	291	zdt = p2dt(jk)
	292	DO jj = 2, jpjm1
	293	DO ji = fs_2, fs_jpim1 ! vector opt.
	294	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
[5845]	295	zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v_n(ji,jj,jk)
[2528]	296	zwy(ji,jj,jk) = ABS( pvn(ji,jj,jk) ) * zdt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)
	297	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
	298	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
	299	END DO
	300	END DO
	301	END DO
	302
	303	!--- Lateral boundary conditions
	304	CALL lbc_lnk( zfu(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zfd(:,:,:), 'T', 1. )
	305	CALL lbc_lnk( zfc(:,:,:), 'T', 1. ) ; CALL lbc_lnk( zwy(:,:,:), 'T', 1. )
	306
[1231]	307	!--- QUICKEST scheme
[2528]	308	CALL quickest( zfu, zfd, zfc, zwy )
[1231]	309	!
[2528]	310	! Mask at the T-points in the x-direction (mask=0 or mask=1)
	311	DO jk = 1, jpkm1
	312	DO jj = 2, jpjm1
	313	DO ji = fs_2, fs_jpim1 ! vector opt.
	314	zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
	315	END DO
[1231]	316	END DO
	317	END DO
[2528]	318	!--- Lateral boundary conditions
	319	CALL lbc_lnk( zfu(:,:,:), 'T', 1. )
	320	!
	321	! Tracer flux on the x-direction
	322	DO jk = 1, jpkm1
	323	!
[1231]	324	DO jj = 2, jpjm1
[2528]	325	DO ji = fs_2, fs_jpim1 ! vector opt.
	326	zdir = 0.5 + SIGN( 0.5, pvn(ji,jj,jk) ) ! if pun > 0 : zdir = 1 otherwise zdir = 0
	327	!--- If the second ustream point is a land point
	328	!--- the flux is computed by the 1st order UPWIND scheme
	329	zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
	330	zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
	331	zwy(ji,jj,jk) = zwy(ji,jj,jk) * pvn(ji,jj,jk)
[1231]	332	END DO
	333	END DO
[3300]	334	END DO
	335	!
	336	CALL lbc_lnk( zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
	337	!
	338	! Computation of the trend
	339	DO jk = 1, jpkm1
[2528]	340	DO jj = 2, jpjm1
	341	DO ji = fs_2, fs_jpim1 ! vector opt.
[6004]	342	zbtr = r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
[2528]	343	! horizontal advective trends
	344	ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
	345	!--- add it to the general tracer trends
	346	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[1231]	347	END DO
	348	END DO
[2528]	349	END DO
[5883]	350	! ! trend diagnostics
[4990]	351	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptn(:,:,:,jn) )
[2528]	352	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
[5147]	353	IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
	354	IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) )
	355	IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) )
[1231]	356	ENDIF
[2528]	357	!
	358	END DO
	359	!
[4990]	360	CALL wrk_dealloc( jpi, jpj, jpk, zwy, zfu, zfc, zfd )
[2715]	361	!
[1559]	362	END SUBROUTINE tra_adv_qck_j
[1231]	363
	364
[2528]	365	SUBROUTINE tra_adv_cen2_k( kt, cdtype, pwn, &
	366	& ptn, pta, kjpt )
[1231]	367	!!----------------------------------------------------------------------
	368	!!
	369	!!----------------------------------------------------------------------
[2715]	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), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pwn ! vertical velocity
	374	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptn ! before and now tracer fields
	375	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
	376	!
[2528]	377	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[4990]	378	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz
[1559]	379	!!----------------------------------------------------------------------
[4990]	380	!
[5836]	381	CALL wrk_alloc( jpi,jpj,jpk, zwz )
	382	!
[5883]	383	zwz(:,:, 1 ) = 0._wp ! surface & bottom values set to zero for all tracers
	384	zwz(:,:,jpk) = 0._wp
[5836]	385	!
[2528]	386	! ! ===========
	387	DO jn = 1, kjpt ! tracer loop
	388	! ! ===========
	389	!
[5836]	390	DO jk = 2, jpkm1 !* Interior point (w-masked 2nd order centered flux)
[2528]	391	DO jj = 2, jpjm1
	392	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	393	zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
[2528]	394	END DO
[1231]	395	END DO
	396	END DO
[5866]	397	IF( ln_linssh ) THEN !* top value (only in linear free surf. as zwz is multiplied by wmask)
[5836]	398	IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean)
	399	DO jj = 1, jpj
	400	DO ji = 1, jpi
	401	zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptn(ji,jj,mikt(ji,jj),jn) ! linear free surface
	402	END DO
	403	END DO
[5883]	404	ELSE ! no ocean cavities (only ocean surface)
[5836]	405	zwz(:,:,1) = pwn(:,:,1) * ptn(:,:,1,jn)
	406	ENDIF
	407	ENDIF
[2528]	408	!
	409	DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==!
	410	DO jj = 2, jpjm1
	411	DO ji = fs_2, fs_jpim1 ! vector opt.
[5836]	412	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) &
[6004]	413	& * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
[2528]	414	END DO
[1231]	415	END DO
	416	END DO
[5883]	417	! ! Send trends for diagnostic
[4990]	418	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, ptn(:,:,:,jn) )
[2528]	419	!
[1231]	420	END DO
	421	!
[5866]	422	CALL wrk_dealloc( jpi,jpj,jpk, zwz )
[4990]	423	!
[1559]	424	END SUBROUTINE tra_adv_cen2_k
[1231]	425
	426
[2528]	427	SUBROUTINE quickest( pfu, pfd, pfc, puc )
[1231]	428	!!----------------------------------------------------------------------
	429	!!
[2528]	430	!! ** Purpose : Computation of advective flux with Quickest scheme
	431	!!
	432	!! ** Method :
[1231]	433	!!----------------------------------------------------------------------
[2528]	434	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfu ! second upwind point
	435	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfd ! first douwning point
	436	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfc ! the central point (or the first upwind point)
	437	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: puc ! input as Courant number ; output as flux
	438	!!
	439	INTEGER :: ji, jj, jk ! dummy loop indices
	440	REAL(wp) :: zcoef1, zcoef2, zcoef3 ! local scalars
	441	REAL(wp) :: zc, zcurv, zfho ! - -
	442	!----------------------------------------------------------------------
[3294]	443	!
	444	IF( nn_timing == 1 ) CALL timing_start('quickest')
	445	!
[2528]	446	DO jk = 1, jpkm1
	447	DO jj = 1, jpj
	448	DO ji = 1, jpi
	449	zc = puc(ji,jj,jk) ! Courant number
	450	zcurv = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
	451	zcoef1 = 0.5 * ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
	452	zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
	453	zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
	454	zfho = zcoef1 - zcoef2 - zcoef3 ! phi_f QUICKEST
	455	!
	456	zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
	457	zcoef2 = ABS( zcoef1 )
	458	zcoef3 = ABS( zcurv )
	459	IF( zcoef3 >= zcoef2 ) THEN
	460	zfho = pfc(ji,jj,jk)
	461	ELSE
	462	zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) ) ! phi_REF
	463	IF( zcoef1 >= 0. ) THEN
	464	zfho = MAX( pfc(ji,jj,jk), zfho )
	465	zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) )
	466	ELSE
	467	zfho = MIN( pfc(ji,jj,jk), zfho )
	468	zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) )
	469	ENDIF
	470	ENDIF
	471	puc(ji,jj,jk) = zfho
	472	END DO
	473	END DO
	474	END DO
[1231]	475	!
[3294]	476	IF( nn_timing == 1 ) CALL timing_stop('quickest')
	477	!
[1231]	478	END SUBROUTINE quickest
	479
	480	!!======================================================================
	481	END MODULE traadv_qck

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

Download in other formats: