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

Context Navigation

source: branches/UKMO/dev_r8864_restart_date/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_qck.F90 @ 9235

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

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

Download in other formats: