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_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_qck.F90 @ 10802

Last change on this file since 10802 was 10802, checked in by davestorkey, 5 years ago
2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps : introduce new T/S variables and convert tracer advection routines (including calls from TOP).
Property svn:keywords set to `Id`
File size: 23.9 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	!
[9124]	24	USE in_out_manager ! I/O manager
[1231]	25	USE lib_mpp ! distribued memory computing
	26	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[3625]	27	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[1231]	28
	29	IMPLICIT NONE
	30	PRIVATE
	31
[1559]	32	PUBLIC tra_adv_qck ! routine called by step.F90
[1231]	33
[2528]	34	REAL(wp) :: r1_6 = 1./ 6. ! 1/6 ratio
[1559]	35
[7646]	36	LOGICAL :: l_trd ! flag to compute trends
	37	LOGICAL :: l_ptr ! flag to compute poleward transport
	38
	39
[1231]	40	!! * Substitutions
	41	# include "vectopt_loop_substitute.h90"
	42	!!----------------------------------------------------------------------
[9598]	43	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1231]	44	!! $Id$
[10068]	45	!! Software governed by the CeCILL license (see ./LICENSE)
[1231]	46	!!----------------------------------------------------------------------
	47	CONTAINS
	48
[10802]	49	SUBROUTINE tra_adv_qck ( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pwn, &
	50	& pt_lev1, pt_lev2, pt_rhs, kjpt )
[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	!!
[10802]	74	!! On the vertical, the simple centered scheme used pt_lev2
[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	!!
[10802]	80	!! ** Action : - update pt_rhs with the now advective tracer trends
[6140]	81	!! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
	82	!! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
[1231]	83	!!
	84	!! ** Reference : Leonard (1979, 1991)
	85	!!----------------------------------------------------------------------
[2528]	86	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[3294]	87	INTEGER , INTENT(in ) :: kit000 ! first time step index
[10802]	88	INTEGER , INTENT(in ) :: ktlev ! time level index for source terms
[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
[10802]	92	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pu, pv, pwn ! 3 ocean velocity components
	93	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt_lev1, pt_lev2 ! before and now tracer fields
	94	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend
[1231]	95	!!----------------------------------------------------------------------
[3294]	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.
[7646]	105	l_ptr = .FALSE.
	106	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
	107	IF( cdtype == 'TRA' .AND. ln_diaptr ) l_ptr = .TRUE.
[4499]	108	!
[7646]	109	!
[6140]	110	! ! horizontal fluxes are computed with the QUICKEST + ULTIMATE scheme
[10802]	111	CALL tra_adv_qck_i( kt, ktlev, cdtype, p2dt, pu, pt_lev1, pt_lev2, pt_rhs, kjpt )
	112	CALL tra_adv_qck_j( kt, ktlev, cdtype, p2dt, pv, pt_lev1, pt_lev2, pt_rhs, kjpt )
[1231]	113
[6140]	114	! ! vertical fluxes are computed with the 2nd order centered scheme
[10802]	115	CALL tra_adv_cen2_k( kt, ktlev, cdtype, pwn, pt_lev2, pt_rhs, kjpt )
[1231]	116	!
	117	END SUBROUTINE tra_adv_qck
	118
	119
[10802]	120	SUBROUTINE tra_adv_qck_i( kt, ktlev, cdtype, p2dt, pu, &
	121	& pt_lev1, pt_lev2, pt_rhs, kjpt )
[1231]	122	!!----------------------------------------------------------------------
	123	!!
	124	!!----------------------------------------------------------------------
[2715]	125	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[10802]	126	INTEGER , INTENT(in ) :: ktlev ! time level index for source terms
[2715]	127	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	128	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[6140]	129	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
[10802]	130	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pu ! i-velocity components
	131	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt_lev1, pt_lev2 ! before and now tracer fields
	132	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend
[2528]	133	!!
[5836]	134	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[6140]	135	REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk ! local scalars
[9019]	136	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx, zfu, zfc, zfd
[1231]	137	!----------------------------------------------------------------------
[2715]	138	!
[2528]	139	! ! ===========
	140	DO jn = 1, kjpt ! tracer loop
	141	! ! ===========
[5836]	142	zfu(:,:,:) = 0._wp ; zfc(:,:,:) = 0._wp
	143	zfd(:,:,:) = 0._wp ; zwx(:,:,:) = 0._wp
	144	!
	145	!!gm why not using a SHIFT instruction...
	146	DO jk = 1, jpkm1 !--- Computation of the ustream and downstream value of the tracer and the mask
[2528]	147	DO jj = 2, jpjm1
	148	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	149	zfc(ji,jj,jk) = pt_lev1(ji-1,jj,jk,jn) ! Upstream in the x-direction for the tracer
	150	zfd(ji,jj,jk) = pt_lev1(ji+1,jj,jk,jn) ! Downstream in the x-direction for the tracer
[2528]	151	END DO
[1559]	152	END DO
	153	END DO
[10425]	154	CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. ) ! Lateral boundary conditions
[2528]	155
[1231]	156	!
	157	! Horizontal advective fluxes
	158	! ---------------------------
[2528]	159	DO jk = 1, jpkm1
	160	DO jj = 2, jpjm1
	161	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	162	zdir = 0.5 + SIGN( 0.5, pu(ji,jj,jk) ) ! if pu > 0 : zdir = 1 otherwise zdir = 0
[2528]	163	zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji+1,jj,jk) ! FU in the x-direction for T
	164	END DO
	165	END DO
[1559]	166	END DO
[1231]	167	!
[2528]	168	DO jk = 1, jpkm1
	169	DO jj = 2, jpjm1
	170	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	171	zdir = 0.5 + SIGN( 0.5, pu(ji,jj,jk) ) ! if pu > 0 : zdir = 1 otherwise zdir = 0
	172	zdx = ( zdir * e1t(ji,jj) + ( 1. - zdir ) * e1t(ji+1,jj) ) * e2u(ji,jj) * e3u(ji,jj,jk,ktlev)
	173	zwx(ji,jj,jk) = ABS( pu(ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)
	174	zfc(ji,jj,jk) = zdir * pt_lev1(ji ,jj,jk,jn) + ( 1. - zdir ) * pt_lev1(ji+1,jj,jk,jn) ! FC in the x-direction for T
	175	zfd(ji,jj,jk) = zdir * pt_lev1(ji+1,jj,jk,jn) + ( 1. - zdir ) * pt_lev1(ji ,jj,jk,jn) ! FD in the x-direction for T
[2528]	176	END DO
	177	END DO
	178	END DO
	179	!--- Lateral boundary conditions
[10425]	180	CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1., zwx(:,:,:), 'T', 1. )
[2528]	181
[1231]	182	!--- QUICKEST scheme
[2528]	183	CALL quickest( zfu, zfd, zfc, zwx )
[1231]	184	!
[2528]	185	! Mask at the T-points in the x-direction (mask=0 or mask=1)
	186	DO jk = 1, jpkm1
	187	DO jj = 2, jpjm1
	188	DO ji = fs_2, fs_jpim1 ! vector opt.
	189	zfu(ji,jj,jk) = tmask(ji-1,jj,jk) + tmask(ji,jj,jk) + tmask(ji+1,jj,jk) - 2.
[2715]	190	END DO
[1231]	191	END DO
	192	END DO
[10425]	193	CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. ) ! Lateral boundary conditions
[2528]	194
[1231]	195	!
[2528]	196	! Tracer flux on the x-direction
	197	DO jk = 1, jpkm1
	198	!
[1231]	199	DO jj = 2, jpjm1
[2528]	200	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	201	zdir = 0.5 + SIGN( 0.5, pu(ji,jj,jk) ) ! if pu > 0 : zdir = 1 otherwise zdir = 0
[2528]	202	!--- If the second ustream point is a land point
	203	!--- the flux is computed by the 1st order UPWIND scheme
	204	zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji+1,jj,jk)
	205	zwx(ji,jj,jk) = zmsk * zwx(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
[10802]	206	zwx(ji,jj,jk) = zwx(ji,jj,jk) * pu(ji,jj,jk)
[1231]	207	END DO
	208	END DO
[3300]	209	END DO
	210	!
[10425]	211	CALL lbc_lnk( 'traadv_qck', zwx(:,:,:), 'T', 1. ) ! Lateral boundary conditions
[3300]	212	!
	213	! Computation of the trend
	214	DO jk = 1, jpkm1
[2528]	215	DO jj = 2, jpjm1
	216	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	217	zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
[2528]	218	! horizontal advective trends
	219	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj,jk) )
	220	!--- add it to the general tracer trends
[10802]	221	pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + ztra
[2528]	222	END DO
	223	END DO
[1231]	224	END DO
[6140]	225	! ! trend diagnostics
[10802]	226	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pu, pt_lev2(:,:,:,jn) )
[2528]	227	!
	228	END DO
	229	!
[1559]	230	END SUBROUTINE tra_adv_qck_i
[1231]	231
	232
[10802]	233	SUBROUTINE tra_adv_qck_j( kt, ktlev, cdtype, p2dt, pv, &
	234	& pt_lev1, pt_lev2, pt_rhs, kjpt )
[1231]	235	!!----------------------------------------------------------------------
	236	!!
	237	!!----------------------------------------------------------------------
[2715]	238	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[10802]	239	INTEGER , INTENT(in ) :: ktlev ! time level index for source terms
[2715]	240	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	241	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[6140]	242	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
[10802]	243	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pv ! j-velocity components
	244	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt_lev1, pt_lev2 ! before and now tracer fields
	245	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend
[1559]	246	!!
[9019]	247	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[6140]	248	REAL(wp) :: ztra, zbtr, zdir, zdx, zmsk ! local scalars
[9019]	249	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwy, zfu, zfc, zfd ! 3D workspace
[1231]	250	!----------------------------------------------------------------------
[2715]	251	!
[2528]	252	! ! ===========
	253	DO jn = 1, kjpt ! tracer loop
	254	! ! ===========
	255	zfu(:,:,:) = 0.0 ; zfc(:,:,:) = 0.0
	256	zfd(:,:,:) = 0.0 ; zwy(:,:,:) = 0.0
	257	!
	258	DO jk = 1, jpkm1
	259	!
	260	!--- Computation of the ustream and downstream value of the tracer and the mask
	261	DO jj = 2, jpjm1
	262	DO ji = fs_2, fs_jpim1 ! vector opt.
	263	! Upstream in the x-direction for the tracer
[10802]	264	zfc(ji,jj,jk) = pt_lev1(ji,jj-1,jk,jn)
[2528]	265	! Downstream in the x-direction for the tracer
[10802]	266	zfd(ji,jj,jk) = pt_lev1(ji,jj+1,jk,jn)
[2528]	267	END DO
[1559]	268	END DO
	269	END DO
[10425]	270	CALL lbc_lnk_multi( 'traadv_qck', zfc(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1. ) ! Lateral boundary conditions
[2528]	271
	272
[1231]	273	!
	274	! Horizontal advective fluxes
	275	! ---------------------------
	276	!
[2528]	277	DO jk = 1, jpkm1
	278	DO jj = 2, jpjm1
	279	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	280	zdir = 0.5 + SIGN( 0.5, pv(ji,jj,jk) ) ! if pu > 0 : zdir = 1 otherwise zdir = 0
[2528]	281	zfu(ji,jj,jk) = zdir * zfc(ji,jj,jk ) + ( 1. - zdir ) * zfd(ji,jj+1,jk) ! FU in the x-direction for T
	282	END DO
[1559]	283	END DO
	284	END DO
[1231]	285	!
[2528]	286	DO jk = 1, jpkm1
	287	DO jj = 2, jpjm1
	288	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	289	zdir = 0.5 + SIGN( 0.5, pv(ji,jj,jk) ) ! if pu > 0 : zdir = 1 otherwise zdir = 0
	290	zdx = ( zdir * e2t(ji,jj) + ( 1. - zdir ) * e2t(ji,jj+1) ) * e1v(ji,jj) * e3v(ji,jj,jk,ktlev)
	291	zwy(ji,jj,jk) = ABS( pv(ji,jj,jk) ) * p2dt / zdx ! (0<zc_cfl<1 : Courant number on x-direction)
	292	zfc(ji,jj,jk) = zdir * pt_lev1(ji,jj ,jk,jn) + ( 1. - zdir ) * pt_lev1(ji,jj+1,jk,jn) ! FC in the x-direction for T
	293	zfd(ji,jj,jk) = zdir * pt_lev1(ji,jj+1,jk,jn) + ( 1. - zdir ) * pt_lev1(ji,jj ,jk,jn) ! FD in the x-direction for T
[2528]	294	END DO
	295	END DO
	296	END DO
	297
	298	!--- Lateral boundary conditions
[10425]	299	CALL lbc_lnk_multi( 'traadv_qck', zfu(:,:,:), 'T', 1. , zfd(:,:,:), 'T', 1., zfc(:,:,:), 'T', 1., zwy(:,:,:), 'T', 1. )
[2528]	300
[1231]	301	!--- QUICKEST scheme
[2528]	302	CALL quickest( zfu, zfd, zfc, zwy )
[1231]	303	!
[2528]	304	! Mask at the T-points in the x-direction (mask=0 or mask=1)
	305	DO jk = 1, jpkm1
	306	DO jj = 2, jpjm1
	307	DO ji = fs_2, fs_jpim1 ! vector opt.
	308	zfu(ji,jj,jk) = tmask(ji,jj-1,jk) + tmask(ji,jj,jk) + tmask(ji,jj+1,jk) - 2.
	309	END DO
[1231]	310	END DO
	311	END DO
[10425]	312	CALL lbc_lnk( 'traadv_qck', zfu(:,:,:), 'T', 1. ) !--- Lateral boundary conditions
[2528]	313	!
	314	! Tracer flux on the x-direction
	315	DO jk = 1, jpkm1
	316	!
[1231]	317	DO jj = 2, jpjm1
[2528]	318	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	319	zdir = 0.5 + SIGN( 0.5, pv(ji,jj,jk) ) ! if pu > 0 : zdir = 1 otherwise zdir = 0
[2528]	320	!--- If the second ustream point is a land point
	321	!--- the flux is computed by the 1st order UPWIND scheme
	322	zmsk = zdir * zfu(ji,jj,jk) + ( 1. - zdir ) * zfu(ji,jj+1,jk)
	323	zwy(ji,jj,jk) = zmsk * zwy(ji,jj,jk) + ( 1. - zmsk ) * zfc(ji,jj,jk)
[10802]	324	zwy(ji,jj,jk) = zwy(ji,jj,jk) * pv(ji,jj,jk)
[1231]	325	END DO
	326	END DO
[3300]	327	END DO
	328	!
[10425]	329	CALL lbc_lnk( 'traadv_qck', zwy(:,:,:), 'T', 1. ) ! Lateral boundary conditions
[3300]	330	!
	331	! Computation of the trend
	332	DO jk = 1, jpkm1
[2528]	333	DO jj = 2, jpjm1
	334	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	335	zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
[2528]	336	! horizontal advective trends
	337	ztra = - zbtr * ( zwy(ji,jj,jk) - zwy(ji,jj-1,jk) )
	338	!--- add it to the general tracer trends
[10802]	339	pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + ztra
[1231]	340	END DO
	341	END DO
[2528]	342	END DO
[6140]	343	! ! trend diagnostics
[10802]	344	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pv, pt_lev2(:,:,:,jn) )
[2528]	345	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
[9019]	346	IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
[2528]	347	!
	348	END DO
	349	!
[1559]	350	END SUBROUTINE tra_adv_qck_j
[1231]	351
	352
[10802]	353	SUBROUTINE tra_adv_cen2_k( kt, ktlev, cdtype, pwn, &
	354	& pt_lev2, pt_rhs, kjpt )
[1231]	355	!!----------------------------------------------------------------------
	356	!!
	357	!!----------------------------------------------------------------------
[2715]	358	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[10802]	359	INTEGER , INTENT(in ) :: ktlev ! time level index for source terms
[2715]	360	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	361	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	362	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pwn ! vertical velocity
[10802]	363	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt_lev2 ! before and now tracer fields
	364	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend
[2715]	365	!
[2528]	366	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[9019]	367	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwz ! 3D workspace
[1559]	368	!!----------------------------------------------------------------------
[4990]	369	!
[6140]	370	zwz(:,:, 1 ) = 0._wp ! surface & bottom values set to zero for all tracers
	371	zwz(:,:,jpk) = 0._wp
[5836]	372	!
[2528]	373	! ! ===========
	374	DO jn = 1, kjpt ! tracer loop
	375	! ! ===========
	376	!
[5836]	377	DO jk = 2, jpkm1 !* Interior point (w-masked 2nd order centered flux)
[2528]	378	DO jj = 2, jpjm1
	379	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	380	zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( pt_lev2(ji,jj,jk-1,jn) + pt_lev2(ji,jj,jk,jn) ) * wmask(ji,jj,jk)
[2528]	381	END DO
[1231]	382	END DO
	383	END DO
[6140]	384	IF( ln_linssh ) THEN !* top value (only in linear free surf. as zwz is multiplied by wmask)
[5836]	385	IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean)
	386	DO jj = 1, jpj
	387	DO ji = 1, jpi
[10802]	388	zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * pt_lev2(ji,jj,mikt(ji,jj),jn) ! linear free surface
[5836]	389	END DO
	390	END DO
[6140]	391	ELSE ! no ocean cavities (only ocean surface)
[10802]	392	zwz(:,:,1) = pwn(:,:,1) * pt_lev2(:,:,1,jn)
[5836]	393	ENDIF
	394	ENDIF
[2528]	395	!
	396	DO jk = 1, jpkm1 !== Tracer flux divergence added to the general trend ==!
	397	DO jj = 2, jpjm1
	398	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	399	pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) - ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) &
	400	& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
[2528]	401	END DO
[1231]	402	END DO
	403	END DO
[6140]	404	! ! Send trends for diagnostic
[10802]	405	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwz, pwn, pt_lev2(:,:,:,jn) )
[2528]	406	!
[1231]	407	END DO
	408	!
[1559]	409	END SUBROUTINE tra_adv_cen2_k
[1231]	410
	411
[2528]	412	SUBROUTINE quickest( pfu, pfd, pfc, puc )
[1231]	413	!!----------------------------------------------------------------------
	414	!!
[2528]	415	!! ** Purpose : Computation of advective flux with Quickest scheme
	416	!!
	417	!! ** Method :
[1231]	418	!!----------------------------------------------------------------------
[2528]	419	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfu ! second upwind point
	420	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfd ! first douwning point
	421	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pfc ! the central point (or the first upwind point)
	422	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: puc ! input as Courant number ; output as flux
	423	!!
	424	INTEGER :: ji, jj, jk ! dummy loop indices
	425	REAL(wp) :: zcoef1, zcoef2, zcoef3 ! local scalars
	426	REAL(wp) :: zc, zcurv, zfho ! - -
	427	!----------------------------------------------------------------------
[3294]	428	!
[2528]	429	DO jk = 1, jpkm1
	430	DO jj = 1, jpj
	431	DO ji = 1, jpi
	432	zc = puc(ji,jj,jk) ! Courant number
	433	zcurv = pfd(ji,jj,jk) + pfu(ji,jj,jk) - 2. * pfc(ji,jj,jk)
	434	zcoef1 = 0.5 * ( pfc(ji,jj,jk) + pfd(ji,jj,jk) )
	435	zcoef2 = 0.5 * zc * ( pfd(ji,jj,jk) - pfc(ji,jj,jk) )
	436	zcoef3 = ( 1. - ( zc * zc ) ) * r1_6 * zcurv
	437	zfho = zcoef1 - zcoef2 - zcoef3 ! phi_f QUICKEST
	438	!
	439	zcoef1 = pfd(ji,jj,jk) - pfu(ji,jj,jk)
	440	zcoef2 = ABS( zcoef1 )
	441	zcoef3 = ABS( zcurv )
	442	IF( zcoef3 >= zcoef2 ) THEN
	443	zfho = pfc(ji,jj,jk)
	444	ELSE
	445	zcoef3 = pfu(ji,jj,jk) + ( ( pfc(ji,jj,jk) - pfu(ji,jj,jk) ) / MAX( zc, 1.e-9 ) ) ! phi_REF
	446	IF( zcoef1 >= 0. ) THEN
	447	zfho = MAX( pfc(ji,jj,jk), zfho )
	448	zfho = MIN( zfho, MIN( zcoef3, pfd(ji,jj,jk) ) )
	449	ELSE
	450	zfho = MIN( pfc(ji,jj,jk), zfho )
	451	zfho = MAX( zfho, MAX( zcoef3, pfd(ji,jj,jk) ) )
	452	ENDIF
	453	ENDIF
	454	puc(ji,jj,jk) = zfho
	455	END DO
	456	END DO
	457	END DO
[1231]	458	!
	459	END SUBROUTINE quickest
	460
	461	!!======================================================================
	462	END MODULE traadv_qck

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

Download in other formats: