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 NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/TRA/traadv_qck.F90 @ 12269

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

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

Download in other formats: