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

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 9019

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