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_fct.F90 in NEMO/branches/UKMO/r8395_coupling_sequence/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: NEMO/branches/UKMO/r8395_coupling_sequence/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_fct.F90 @ 10761

Last change on this file since 10761 was 10761, checked in by jcastill, 5 years ago
Remove svn keys
File size: 49.3 KB

Rev	Line
[5770]	1	MODULE traadv_fct
[3]	2	!!==============================================================================
[5770]	3	!! * MODULE traadv_fct *
	4	!! Ocean tracers: horizontal & vertical advective trend (2nd/4th order Flux Corrected Transport method)
[3]	5	!!==============================================================================
[5770]	6	!! History : 3.7 ! 2015-09 (L. Debreu, G. Madec) original code (inspired from traadv_tvd.F90)
[503]	7	!!----------------------------------------------------------------------
[3]	8
	9	!!----------------------------------------------------------------------
[5770]	10	!! tra_adv_fct : update the tracer trend with a 3D advective trends using a 2nd or 4th order FCT scheme
	11	!! tra_adv_fct_zts: update the tracer trend with a 3D advective trends using a 2nd order FCT scheme
	12	!! with sub-time-stepping in the vertical direction
	13	!! nonosc : compute monotonic tracer fluxes by a non-oscillatory algorithm
	14	!! interp_4th_cpt : 4th order compact scheme for the vertical component of the advection
[3]	15	!!----------------------------------------------------------------------
[3625]	16	USE oce ! ocean dynamics and active tracers
	17	USE dom_oce ! ocean space and time domain
[4990]	18	USE trc_oce ! share passive tracers/Ocean variables
	19	USE trd_oce ! trends: ocean variables
[3625]	20	USE trdtra ! tracers trends
[4990]	21	USE diaptr ! poleward transport diagnostics
[7646]	22	USE diaar5 ! AR5 diagnostics
	23	USE phycst, ONLY: rau0_rcp
[4990]	24	!
[5770]	25	USE in_out_manager ! I/O manager
[7646]	26	USE iom
[3625]	27	USE lib_mpp ! MPP library
	28	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[5770]	29	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[3625]	30	USE wrk_nemo ! Memory Allocation
	31	USE timing ! Timing
[3]	32
	33	IMPLICIT NONE
	34	PRIVATE
	35
[5770]	36	PUBLIC tra_adv_fct ! routine called by traadv.F90
	37	PUBLIC tra_adv_fct_zts ! routine called by traadv.F90
	38	PUBLIC interp_4th_cpt ! routine called by traadv_cen.F90
[3]	39
[5770]	40	LOGICAL :: l_trd ! flag to compute trends
[7646]	41	LOGICAL :: l_ptr ! flag to compute poleward transport
	42	LOGICAL :: l_hst ! flag to compute heat/salt transport
[5770]	43	REAL(wp) :: r1_6 = 1._wp / 6._wp ! =1/6
[2528]	44
[7646]	45	! ! tridiag solver associated indices:
	46	INTEGER, PARAMETER :: np_NH = 0 ! Neumann homogeneous boundary condition
	47	INTEGER, PARAMETER :: np_CEN2 = 1 ! 2nd order centered boundary condition
	48
[3]	49	!! * Substitutions
	50	# include "vectopt_loop_substitute.h90"
	51	!!----------------------------------------------------------------------
[5770]	52	!! NEMO/OPA 3.7 , NEMO Consortium (2014)
[1152]	53	!! $Id$
[2528]	54	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	55	!!----------------------------------------------------------------------
	56	CONTAINS
	57
[5770]	58	SUBROUTINE tra_adv_fct( kt, kit000, cdtype, p2dt, pun, pvn, pwn, &
	59	& ptb, ptn, pta, kjpt, kn_fct_h, kn_fct_v )
[3]	60	!!----------------------------------------------------------------------
[5770]	61	!! * ROUTINE tra_adv_fct *
[3]	62	!!
[6140]	63	!! ** Purpose : Compute the now trend due to total advection of tracers
	64	!! and add it to the general trend of tracer equations
[3]	65	!!
[5770]	66	!! ** Method : - 2nd or 4th FCT scheme on the horizontal direction
	67	!! (choice through the value of kn_fct)
[6140]	68	!! - on the vertical the 4th order is a compact scheme
[5770]	69	!! - corrected flux (monotonic correction)
[3]	70	!!
[6140]	71	!! ** Action : - update pta with the now advective tracer trends
	72	!! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
	73	!! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
[503]	74	!!----------------------------------------------------------------------
[2528]	75	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[3294]	76	INTEGER , INTENT(in ) :: kit000 ! first time step index
[2528]	77	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	78	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[5770]	79	INTEGER , INTENT(in ) :: kn_fct_h ! order of the FCT scheme (=2 or 4)
	80	INTEGER , INTENT(in ) :: kn_fct_v ! order of the FCT scheme (=2 or 4)
[6140]	81	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
[2528]	82	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
	83	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
	84	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[2715]	85	!
[5770]	86	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[6140]	87	REAL(wp) :: ztra ! local scalar
[5770]	88	REAL(wp) :: zfp_ui, zfp_vj, zfp_wk, zC2t_u, zC4t_u ! - -
	89	REAL(wp) :: zfm_ui, zfm_vj, zfm_wk, zC2t_v, zC4t_v ! - -
	90	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw
[7646]	91	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz, zptry
	92	REAL(wp), POINTER, DIMENSION(:,:) :: z2d
[3]	93	!!----------------------------------------------------------------------
[3294]	94	!
[5770]	95	IF( nn_timing == 1 ) CALL timing_start('tra_adv_fct')
[3294]	96	!
[5770]	97	CALL wrk_alloc( jpi,jpj,jpk, zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw )
[3294]	98	!
	99	IF( kt == kit000 ) THEN
[2528]	100	IF(lwp) WRITE(numout,*)
[5770]	101	IF(lwp) WRITE(numout,*) 'tra_adv_fct : FCT advection scheme on ', cdtype
[2528]	102	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
[3]	103	ENDIF
[2528]	104	!
[5770]	105	l_trd = .FALSE.
[7646]	106	l_hst = .FALSE.
	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.
	110	IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
	111	& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
[5770]	112	!
[7646]	113	IF( l_trd .OR. l_hst ) THEN
[3294]	114	CALL wrk_alloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
[7753]	115	ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp
[3294]	116	ENDIF
[2528]	117	!
[7646]	118	IF( l_ptr ) THEN
	119	CALL wrk_alloc( jpi, jpj, jpk, zptry )
[7753]	120	zptry(:,:,:) = 0._wp
[7646]	121	ENDIF
[6140]	122	! ! surface & bottom value : flux set to zero one for all
[7753]	123	zwz(:,:, 1 ) = 0._wp
	124	zwx(:,:,jpk) = 0._wp ; zwy(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp
[2528]	125	!
[7753]	126	zwi(:,:,:) = 0._wp
	127	!
[6140]	128	DO jn = 1, kjpt !== loop over the tracers ==!
[5770]	129	!
	130	! !== upstream advection with initial mass fluxes & intermediate update ==!
	131	! !* upstream tracer flux in the i and j direction
[2528]	132	DO jk = 1, jpkm1
	133	DO jj = 1, jpjm1
	134	DO ji = 1, fs_jpim1 ! vector opt.
	135	! upstream scheme
	136	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
	137	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
	138	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
	139	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
	140	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) )
	141	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) )
	142	END DO
[3]	143	END DO
	144	END DO
[5770]	145	! !* upstream tracer flux in the k direction *!
[6140]	146	DO jk = 2, jpkm1 ! Interior value ( multiplied by wmask)
[4990]	147	DO jj = 1, jpj
	148	DO ji = 1, jpi
[2528]	149	zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
	150	zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
[5120]	151	zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
[2528]	152	END DO
[3]	153	END DO
	154	END DO
[6140]	155	IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as zwz has been w-masked)
[5770]	156	IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface
[5120]	157	DO jj = 1, jpj
	158	DO ji = 1, jpi
	159	zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface
	160	END DO
	161	END DO
[5770]	162	ELSE ! no cavities: only at the ocean surface
[7753]	163	zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
[5770]	164	ENDIF
[5120]	165	ENDIF
[5770]	166	!
	167	DO jk = 1, jpkm1 !* trend and after field with monotonic scheme
[216]	168	DO jj = 2, jpjm1
	169	DO ji = fs_2, fs_jpim1 ! vector opt.
[6691]	170	! ! total intermediate advective trends
[5770]	171	ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	172	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
[6691]	173	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
	174	! ! update and guess with monotonic sheme
	175	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
	176	zwi(ji,jj,jk) = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
[216]	177	END DO
	178	END DO
	179	END DO
[5770]	180	CALL lbc_lnk( zwi, 'T', 1. ) ! Lateral boundary conditions on zwi (unchanged sign)
	181	!
[7646]	182	IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes)
[7753]	183	ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)
[2528]	184	END IF
[5770]	185	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
[7753]	186	IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:)
[5770]	187	!
	188	! !== anti-diffusive flux : high order minus low order ==!
	189	!
[6140]	190	SELECT CASE( kn_fct_h ) !* horizontal anti-diffusive fluxes
[5770]	191	!
[6140]	192	CASE( 2 ) !- 2nd order centered
[5770]	193	DO jk = 1, jpkm1
	194	DO jj = 1, jpjm1
	195	DO ji = 1, fs_jpim1 ! vector opt.
	196	zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk)
	197	zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk)
	198	END DO
[503]	199	END DO
	200	END DO
[5770]	201	!
[6140]	202	CASE( 4 ) !- 4th order centered
[7753]	203	zltu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero
	204	zltv(:,:,jpk) = 0._wp
[6140]	205	DO jk = 1, jpkm1 ! Laplacian
	206	DO jj = 1, jpjm1 ! 1st derivative (gradient)
[5770]	207	DO ji = 1, fs_jpim1 ! vector opt.
	208	ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
	209	ztv(ji,jj,jk) = ( ptn(ji ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
	210	END DO
[503]	211	END DO
[6140]	212	DO jj = 2, jpjm1 ! 2nd derivative * 1/ 6
[5770]	213	DO ji = fs_2, fs_jpim1 ! vector opt.
	214	zltu(ji,jj,jk) = ( ztu(ji,jj,jk) + ztu(ji-1,jj,jk) ) * r1_6
	215	zltv(ji,jj,jk) = ( ztv(ji,jj,jk) + ztv(ji,jj-1,jk) ) * r1_6
	216	END DO
	217	END DO
[503]	218	END DO
[5770]	219	CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn)
	220	!
[6140]	221	DO jk = 1, jpkm1 ! Horizontal advective fluxes
[5770]	222	DO jj = 1, jpjm1
	223	DO ji = 1, fs_jpim1 ! vector opt.
	224	zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ! 2 x C2 interpolation of T at u- & v-points
	225	zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn)
	226	! ! C4 minus upstream advective fluxes
	227	zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( zC2t_u + zltu(ji,jj,jk) - zltu(ji+1,jj,jk) ) - zwx(ji,jj,jk)
	228	zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( zC2t_v + zltv(ji,jj,jk) - zltv(ji,jj+1,jk) ) - zwy(ji,jj,jk)
	229	END DO
[5120]	230	END DO
[5770]	231	END DO
	232	!
[6140]	233	CASE( 41 ) !- 4th order centered ==>> !!gm coding attempt need to be tested
[7753]	234	ztu(:,:,jpk) = 0._wp ! Bottom value : flux set to zero
	235	ztv(:,:,jpk) = 0._wp
[6140]	236	DO jk = 1, jpkm1 ! 1st derivative (gradient)
	237	DO jj = 1, jpjm1
[5770]	238	DO ji = 1, fs_jpim1 ! vector opt.
	239	ztu(ji,jj,jk) = ( ptn(ji+1,jj ,jk,jn) - ptn(ji,jj,jk,jn) ) * umask(ji,jj,jk)
	240	ztv(ji,jj,jk) = ( ptn(ji ,jj+1,jk,jn) - ptn(ji,jj,jk,jn) ) * vmask(ji,jj,jk)
	241	END DO
	242	END DO
[5120]	243	END DO
[5770]	244	CALL lbc_lnk( ztu, 'U', -1. ) ; CALL lbc_lnk( ztv, 'V', -1. ) ! Lateral boundary cond. (unchanged sgn)
	245	!
[6140]	246	DO jk = 1, jpkm1 ! Horizontal advective fluxes
[5770]	247	DO jj = 2, jpjm1
	248	DO ji = 2, fs_jpim1 ! vector opt.
	249	zC2t_u = ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) ! 2 x C2 interpolation of T at u- & v-points (x2)
	250	zC2t_v = ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn)
	251	! ! C4 interpolation of T at u- & v-points (x2)
	252	zC4t_u = zC2t_u + r1_6 * ( ztu(ji-1,jj ,jk) - ztu(ji+1,jj ,jk) )
	253	zC4t_v = zC2t_v + r1_6 * ( ztv(ji ,jj-1,jk) - ztv(ji ,jj+1,jk) )
	254	! ! C4 minus upstream advective fluxes
	255	zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * zC4t_u - zwx(ji,jj,jk)
	256	zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * zC4t_v - zwy(ji,jj,jk)
	257	END DO
	258	END DO
	259	END DO
	260	!
	261	END SELECT
[6140]	262	!
	263	SELECT CASE( kn_fct_v ) !* vertical anti-diffusive fluxes (w-masked interior values)
[5770]	264	!
[6140]	265	CASE( 2 ) !- 2nd order centered
[5770]	266	DO jk = 2, jpkm1
	267	DO jj = 2, jpjm1
	268	DO ji = fs_2, fs_jpim1
[6140]	269	zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * 0.5_wp * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) &
	270	& - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
[5770]	271	END DO
	272	END DO
	273	END DO
	274	!
[6140]	275	CASE( 4 ) !- 4th order COMPACT
	276	CALL interp_4th_cpt( ptn(:,:,:,jn) , ztw ) ! zwt = COMPACT interpolation of T at w-point
[5770]	277	DO jk = 2, jpkm1
	278	DO jj = 2, jpjm1
	279	DO ji = fs_2, fs_jpim1
	280	zwz(ji,jj,jk) = ( pwn(ji,jj,jk) * ztw(ji,jj,jk) - zwz(ji,jj,jk) ) * wmask(ji,jj,jk)
	281	END DO
	282	END DO
	283	END DO
	284	!
	285	END SELECT
[6140]	286	IF( ln_linssh ) THEN ! top ocean value: high order = upstream ==>> zwz=0
[7753]	287	zwz(:,:,1) = 0._wp ! only ocean surface as interior zwz values have been w-masked
[6140]	288	ENDIF
[5770]	289	!
[2528]	290	CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions
	291	CALL lbc_lnk( zwz, 'W', 1. )
[6140]	292	!
[5770]	293	! !== monotonicity algorithm ==!
	294	!
[2528]	295	CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
[6140]	296	!
[5770]	297	! !== final trend with corrected fluxes ==!
	298	!
[216]	299	DO jk = 1, jpkm1
	300	DO jj = 2, jpjm1
[2528]	301	DO ji = fs_2, fs_jpim1 ! vector opt.
[5770]	302	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	303	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
	304	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) &
[6140]	305	& * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
[216]	306	END DO
	307	END DO
	308	END DO
[5770]	309	!
[7646]	310	IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes)
[7753]	311	ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed
	312	ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed
	313	ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed
[7646]	314	ENDIF
[5770]	315	!
[7646]	316	IF( l_trd ) THEN
[5770]	317	CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
	318	CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
	319	CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
	320	!
[2528]	321	END IF
[7646]	322	! ! heat/salt transport
	323	IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
	324
	325	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	326	IF( l_ptr ) THEN
[7753]	327	zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed
[7646]	328	CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
[2528]	329	ENDIF
[503]	330	!
[6140]	331	END DO ! end of tracer loop
[503]	332	!
[7646]	333	CALL wrk_dealloc( jpi,jpj,jpk, zwi, zwx, zwy, zwz, ztu, ztv, zltu, zltv, ztw )
	334	IF( l_trd .OR. l_hst ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
	335	IF( l_ptr ) CALL wrk_dealloc( jpi, jpj, jpk, zptry )
[2528]	336	!
[5770]	337	IF( nn_timing == 1 ) CALL timing_stop('tra_adv_fct')
[2715]	338	!
[5770]	339	END SUBROUTINE tra_adv_fct
[3]	340
[5737]	341
[5770]	342	SUBROUTINE tra_adv_fct_zts( kt, kit000, cdtype, p2dt, pun, pvn, pwn, &
	343	& ptb, ptn, pta, kjpt, kn_fct_zts )
[4990]	344	!!----------------------------------------------------------------------
[5770]	345	!! * ROUTINE tra_adv_fct_zts *
[4990]	346	!!
	347	!! ** Purpose : Compute the now trend due to total advection of
	348	!! tracers and add it to the general trend of tracer equations
	349	!!
	350	!! ** Method : TVD ZTS scheme, i.e. 2nd order centered scheme with
	351	!! corrected flux (monotonic correction). This version use sub-
	352	!! timestepping for the vertical advection which increases stability
	353	!! when vertical metrics are small.
	354	!! note: - this advection scheme needs a leap-frog time scheme
	355	!!
	356	!! ** Action : - update (pta) with the now advective tracer trends
	357	!! - save the trends
	358	!!----------------------------------------------------------------------
	359	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	360	INTEGER , INTENT(in ) :: kit000 ! first time step index
	361	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	362	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[5770]	363	INTEGER , INTENT(in ) :: kn_fct_zts ! number of number of vertical sub-timesteps
[6140]	364	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
[4990]	365	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
	366	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
	367	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
	368	!
	369	REAL(wp), DIMENSION( jpk ) :: zts ! length of sub-timestep for vertical advection
[6140]	370	REAL(wp) :: zr_p2dt ! reciprocal of tracer timestep
[4990]	371	INTEGER :: ji, jj, jk, jl, jn ! dummy loop indices
	372	INTEGER :: jtb, jtn, jta ! sub timestep pointers for leap-frog/euler forward steps
	373	INTEGER :: jtaken ! toggle for collecting appropriate fluxes from sub timesteps
	374	REAL(wp) :: z_rzts ! Fractional length of Euler forward sub-timestep for vertical advection
[6140]	375	REAL(wp) :: ztra ! local scalar
[4990]	376	REAL(wp) :: zfp_ui, zfp_vj, zfp_wk ! - -
	377	REAL(wp) :: zfm_ui, zfm_vj, zfm_wk ! - -
[5770]	378	REAL(wp), POINTER, DIMENSION(:,: ) :: zwx_sav , zwy_sav
	379	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwi, zwx, zwy, zwz, zhdiv, zwzts, zwz_sav
	380	REAL(wp), POINTER, DIMENSION(:,:,:) :: ztrdx, ztrdy, ztrdz
[7646]	381	REAL(wp), POINTER, DIMENSION(:,:,:) :: zptry
[5770]	382	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: ztrs
[4990]	383	!!----------------------------------------------------------------------
	384	!
[5770]	385	IF( nn_timing == 1 ) CALL timing_start('tra_adv_fct_zts')
[4990]	386	!
[5770]	387	CALL wrk_alloc( jpi,jpj, zwx_sav, zwy_sav )
[6691]	388	CALL wrk_alloc( jpi,jpj,jpk, zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav )
[5770]	389	CALL wrk_alloc( jpi,jpj,jpk,kjpt+1, ztrs )
[4990]	390	!
	391	IF( kt == kit000 ) THEN
	392	IF(lwp) WRITE(numout,*)
[5770]	393	IF(lwp) WRITE(numout,*) 'tra_adv_fct_zts : 2nd order FCT scheme with ', kn_fct_zts, ' vertical sub-timestep on ', cdtype
[4990]	394	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
	395	ENDIF
	396	!
	397	l_trd = .FALSE.
[7646]	398	l_hst = .FALSE.
	399	l_ptr = .FALSE.
	400	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
	401	IF( cdtype == 'TRA' .AND. ln_diaptr ) l_ptr = .TRUE.
	402	IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
	403	& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
[4990]	404	!
[7646]	405	IF( l_trd .OR. l_hst ) THEN
[5770]	406	CALL wrk_alloc( jpi,jpj,jpk, ztrdx, ztrdy, ztrdz )
[4990]	407	ztrdx(:,:,:) = 0._wp ; ztrdy(:,:,:) = 0._wp ; ztrdz(:,:,:) = 0._wp
	408	ENDIF
	409	!
[7646]	410	IF( l_ptr ) THEN
	411	CALL wrk_alloc( jpi, jpj,jpk, zptry )
	412	zptry(:,:,:) = 0._wp
	413	ENDIF
[4990]	414	zwi(:,:,:) = 0._wp
[5770]	415	z_rzts = 1._wp / REAL( kn_fct_zts, wp )
[6140]	416	zr_p2dt = 1._wp / p2dt
[4990]	417	!
[6140]	418	! surface & Bottom value : flux set to zero for all tracers
	419	zwz(:,:, 1 ) = 0._wp
	420	zwx(:,:,jpk) = 0._wp ; zwz(:,:,jpk) = 0._wp
	421	zwy(:,:,jpk) = 0._wp ; zwi(:,:,jpk) = 0._wp
	422	!
[4990]	423	! ! ===========
	424	DO jn = 1, kjpt ! tracer loop
	425	! ! ===========
[6140]	426	!
	427	! Upstream advection with initial mass fluxes & intermediate update
	428	DO jk = 1, jpkm1 ! upstream tracer flux in the i and j direction
[4990]	429	DO jj = 1, jpjm1
	430	DO ji = 1, fs_jpim1 ! vector opt.
	431	! upstream scheme
	432	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
	433	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
	434	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
	435	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
	436	zwx(ji,jj,jk) = 0.5_wp * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) )
	437	zwy(ji,jj,jk) = 0.5_wp * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) )
	438	END DO
	439	END DO
	440	END DO
[6140]	441	! ! upstream tracer flux in the k direction
	442	DO jk = 2, jpkm1 ! Interior value
[4990]	443	DO jj = 1, jpj
	444	DO ji = 1, jpi
	445	zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
	446	zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
[5770]	447	zwz(ji,jj,jk) = 0.5_wp * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
[4990]	448	END DO
	449	END DO
	450	END DO
[6140]	451	IF( ln_linssh ) THEN ! top value : linear free surface case only (as zwz is multiplied by wmask)
	452	IF( ln_isfcav ) THEN ! ice-shelf cavities: top value
[5120]	453	DO jj = 1, jpj
	454	DO ji = 1, jpi
[5770]	455	zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)
[5120]	456	END DO
[5770]	457	END DO
[6140]	458	ELSE ! no cavities, surface value
[5770]	459	zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
	460	ENDIF
[5120]	461	ENDIF
[5770]	462	!
	463	DO jk = 1, jpkm1 ! total advective trend
[4990]	464	DO jj = 2, jpjm1
	465	DO ji = fs_2, fs_jpim1 ! vector opt.
[6140]	466	! ! total intermediate advective trends
[5770]	467	ztra = - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	468	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
[6691]	469	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) * r1_e1e2t(ji,jj)
[6140]	470	! ! update and guess with monotonic sheme
[6691]	471	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra / e3t_n(ji,jj,jk) * tmask(ji,jj,jk)
	472	zwi(ji,jj,jk) = ( e3t_b(ji,jj,jk) * ptb(ji,jj,jk,jn) + p2dt * ztra ) / e3t_a(ji,jj,jk) * tmask(ji,jj,jk)
[4990]	473	END DO
	474	END DO
	475	END DO
[5770]	476	!
	477	CALL lbc_lnk( zwi, 'T', 1. ) ! Lateral boundary conditions on zwi (unchanged sign)
	478	!
[7646]	479	IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes)
[4990]	480	ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)
	481	END IF
[5770]	482	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
[7646]	483	IF( l_ptr ) zptry(:,:,:) = zwy(:,:,:)
[4990]	484
[5770]	485	! 3. anti-diffusive flux : high order minus low order
	486	! ---------------------------------------------------
[4990]	487
[5770]	488	DO jk = 1, jpkm1 !* horizontal anti-diffusive fluxes
	489	!
[4990]	490	DO jj = 1, jpjm1
	491	DO ji = 1, fs_jpim1 ! vector opt.
	492	zwx_sav(ji,jj) = zwx(ji,jj,jk)
	493	zwy_sav(ji,jj) = zwy(ji,jj,jk)
[5770]	494	!
[4990]	495	zwx(ji,jj,jk) = 0.5_wp * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) )
	496	zwy(ji,jj,jk) = 0.5_wp * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) )
	497	END DO
	498	END DO
[5770]	499	!
	500	DO jj = 2, jpjm1 ! partial horizontal divergence
[4990]	501	DO ji = fs_2, fs_jpim1
	502	zhdiv(ji,jj,jk) = ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) &
	503	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) )
	504	END DO
	505	END DO
[5770]	506	!
[4990]	507	DO jj = 1, jpjm1
	508	DO ji = 1, fs_jpim1 ! vector opt.
[5770]	509	zwx(ji,jj,jk) = zwx(ji,jj,jk) - zwx_sav(ji,jj)
	510	zwy(ji,jj,jk) = zwy(ji,jj,jk) - zwy_sav(ji,jj)
[4990]	511	END DO
	512	END DO
	513	END DO
[6140]	514	!
[5770]	515	! !* vertical anti-diffusive flux
	516	zwz_sav(:,:,:) = zwz(:,:,:)
	517	ztrs (:,:,:,1) = ptb(:,:,:,jn)
[6771]	518	ztrs (:,:,1,2) = ptb(:,:,1,jn)
	519	ztrs (:,:,1,3) = ptb(:,:,1,jn)
[5770]	520	zwzts (:,:,:) = 0._wp
[4990]	521	!
[5770]	522	DO jl = 1, kn_fct_zts ! Start of sub timestepping loop
	523	!
	524	IF( jl == 1 ) THEN ! Euler forward to kick things off
	525	jtb = 1 ; jtn = 1 ; jta = 2
[6140]	526	zts(:) = p2dt * z_rzts
[5770]	527	jtaken = MOD( kn_fct_zts + 1 , 2) ! Toggle to collect every second flux
	528	! ! starting at jl =1 if kn_fct_zts is odd;
	529	! ! starting at jl =2 otherwise
	530	ELSEIF( jl == 2 ) THEN ! First leapfrog step
	531	jtb = 1 ; jtn = 2 ; jta = 3
[6140]	532	zts(:) = 2._wp * p2dt * z_rzts
[5770]	533	ELSE ! Shuffle pointers for subsequent leapfrog steps
	534	jtb = MOD(jtb,3) + 1
	535	jtn = MOD(jtn,3) + 1
	536	jta = MOD(jta,3) + 1
[4990]	537	ENDIF
[5770]	538	DO jk = 2, jpkm1 ! interior value
[4990]	539	DO jj = 2, jpjm1
	540	DO ji = fs_2, fs_jpim1
[5770]	541	zwz(ji,jj,jk) = 0.5_wp * pwn(ji,jj,jk) * ( ztrs(ji,jj,jk,jtn) + ztrs(ji,jj,jk-1,jtn) ) * wmask(ji,jj,jk)
	542	IF( jtaken == 0 ) zwzts(ji,jj,jk) = zwzts(ji,jj,jk) + zwz(ji,jj,jk) * zts(jk) ! Accumulate time-weighted vertcal flux
[4990]	543	END DO
	544	END DO
	545	END DO
[6140]	546	IF( ln_linssh ) THEN ! top value (only in linear free surface case)
[5770]	547	IF( ln_isfcav ) THEN ! ice-shelf cavities
	548	DO jj = 1, jpj
	549	DO ji = 1, jpi
	550	zwz(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn) ! linear free surface
	551	END DO
	552	END DO
[6140]	553	ELSE ! no ocean cavities
[5770]	554	zwz(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
	555	ENDIF
	556	ENDIF
	557	!
[4990]	558	jtaken = MOD( jtaken + 1 , 2 )
[5770]	559	!
	560	DO jk = 2, jpkm1 ! total advective trends
[4990]	561	DO jj = 2, jpjm1
	562	DO ji = fs_2, fs_jpim1
[5770]	563	ztrs(ji,jj,jk,jta) = ztrs(ji,jj,jk,jtb) &
	564	& - zts(jk) * ( zhdiv(ji,jj,jk) + zwz(ji,jj,jk) - zwz(ji,jj,jk+1) ) &
[6140]	565	& * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
[4990]	566	END DO
	567	END DO
	568	END DO
[5770]	569	!
[4990]	570	END DO
	571
	572	DO jk = 2, jpkm1 ! Anti-diffusive vertical flux using average flux from the sub-timestepping
	573	DO jj = 2, jpjm1
	574	DO ji = fs_2, fs_jpim1
[6140]	575	zwz(ji,jj,jk) = ( zwzts(ji,jj,jk) * zr_p2dt - zwz_sav(ji,jj,jk) ) * wmask(ji,jj,jk)
[4990]	576	END DO
	577	END DO
	578	END DO
	579	CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! Lateral bondary conditions
	580	CALL lbc_lnk( zwz, 'W', 1. )
	581
	582	! 4. monotonicity algorithm
	583	! -------------------------
	584	CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
	585
	586
	587	! 5. final trend with corrected fluxes
	588	! ------------------------------------
	589	DO jk = 1, jpkm1
	590	DO jj = 2, jpjm1
	591	DO ji = fs_2, fs_jpim1 ! vector opt.
[5770]	592	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ( zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
	593	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) &
[6140]	594	& * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
[4990]	595	END DO
	596	END DO
	597	END DO
	598
[7646]	599	!
	600	IF( l_trd .OR. l_hst ) THEN ! trend diagnostics (contribution of upstream fluxes)
[4990]	601	ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed
	602	ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed
	603	ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed
[7646]	604	ENDIF
[5770]	605	!
[7646]	606	IF( l_trd ) THEN
	607	CALL trd_tra( kt, cdtype, jn, jptra_xad, ztrdx, pun, ptn(:,:,:,jn) )
	608	CALL trd_tra( kt, cdtype, jn, jptra_yad, ztrdy, pvn, ptn(:,:,:,jn) )
	609	CALL trd_tra( kt, cdtype, jn, jptra_zad, ztrdz, pwn, ptn(:,:,:,jn) )
[5770]	610	!
[4990]	611	END IF
[7646]	612	! ! heat/salt transport
	613	IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztrdx(:,:,:), ztrdy(:,:,:) )
	614
	615	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	616	IF( l_ptr ) THEN
	617	zptry(:,:,:) = zptry(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed
	618	CALL dia_ptr_hst( jn, 'adv', zptry(:,:,:) )
[4990]	619	ENDIF
	620	!
	621	END DO
	622	!
[7646]	623	CALL wrk_alloc( jpi,jpj, zwx_sav, zwy_sav )
	624	CALL wrk_alloc( jpi,jpj, jpk, zwx, zwy, zwz, zwi, zhdiv, zwzts, zwz_sav )
	625	CALL wrk_alloc( jpi,jpj,jpk,kjpt+1, ztrs )
	626	IF( l_trd .OR. l_hst ) CALL wrk_dealloc( jpi, jpj, jpk, ztrdx, ztrdy, ztrdz )
	627	IF( l_ptr ) CALL wrk_dealloc( jpi, jpj, jpk, zptry )
[4990]	628	!
[5770]	629	IF( nn_timing == 1 ) CALL timing_stop('tra_adv_fct_zts')
[4990]	630	!
[5770]	631	END SUBROUTINE tra_adv_fct_zts
[4990]	632
[5737]	633
[2528]	634	SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
[3]	635	!!---------------------------------------------------------------------
	636	!! * ROUTINE nonosc *
	637	!!
	638	!! ** Purpose : compute monotonic tracer fluxes from the upstream
	639	!! scheme and the before field by a nonoscillatory algorithm
	640	!!
	641	!! ** Method : ... ???
	642	!! warning : pbef and paft must be masked, but the boundaries
	643	!! conditions on the fluxes are not necessary zalezak (1979)
	644	!! drange (1995) multi-dimensional forward-in-time and upstream-
	645	!! in-space based differencing for fluid
	646	!!----------------------------------------------------------------------
[6140]	647	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
[2528]	648	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in ) :: pbef, paft ! before & after field
	649	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(inout) :: paa, pbb, pcc ! monotonic fluxes in the 3 directions
[2715]	650	!
[4990]	651	INTEGER :: ji, jj, jk ! dummy loop indices
	652	INTEGER :: ikm1 ! local integer
[6140]	653	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars
[2715]	654	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv, zup, zdo ! - -
[3294]	655	REAL(wp), POINTER, DIMENSION(:,:,:) :: zbetup, zbetdo, zbup, zbdo
[3]	656	!!----------------------------------------------------------------------
[3294]	657	!
	658	IF( nn_timing == 1 ) CALL timing_start('nonosc')
	659	!
	660	CALL wrk_alloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
	661	!
[2715]	662	zbig = 1.e+40_wp
	663	zrtrn = 1.e-15_wp
[7753]	664	zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp
[785]	665
[3]	666	! Search local extrema
	667	! --------------------
[785]	668	! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
[4990]	669	zbup = MAX( pbef * tmask - zbig * ( 1._wp - tmask ), &
	670	& paft * tmask - zbig * ( 1._wp - tmask ) )
	671	zbdo = MIN( pbef * tmask + zbig * ( 1._wp - tmask ), &
	672	& paft * tmask + zbig * ( 1._wp - tmask ) )
[785]	673
[5120]	674	DO jk = 1, jpkm1
	675	ikm1 = MAX(jk-1,1)
	676	DO jj = 2, jpjm1
	677	DO ji = fs_2, fs_jpim1 ! vector opt.
	678
[785]	679	! search maximum in neighbourhood
	680	zup = MAX( zbup(ji ,jj ,jk ), &
	681	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
	682	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), &
	683	& zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) )
[3]	684
[785]	685	! search minimum in neighbourhood
	686	zdo = MIN( zbdo(ji ,jj ,jk ), &
	687	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
	688	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), &
	689	& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
[3]	690
[785]	691	! positive part of the flux
[3]	692	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
	693	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
	694	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
[785]	695
	696	! negative part of the flux
[3]	697	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
	698	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
	699	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
[785]	700
[3]	701	! up & down beta terms
[6140]	702	zbt = e1e2t(ji,jj) * e3t_n(ji,jj,jk) / p2dt
[785]	703	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
	704	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
[3]	705	END DO
	706	END DO
	707	END DO
[2528]	708	CALL lbc_lnk( zbetup, 'T', 1. ) ; CALL lbc_lnk( zbetdo, 'T', 1. ) ! lateral boundary cond. (unchanged sign)
[3]	709
[237]	710	! 3. monotonic flux in the i & j direction (paa & pbb)
	711	! ----------------------------------------
[3]	712	DO jk = 1, jpkm1
	713	DO jj = 2, jpjm1
	714	DO ji = fs_2, fs_jpim1 ! vector opt.
[4990]	715	zau = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
	716	zbu = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
[785]	717	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
[4990]	718	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1._wp - zcu) * zbu )
[3]	719
[4990]	720	zav = MIN( 1._wp, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
	721	zbv = MIN( 1._wp, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
[785]	722	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
[4990]	723	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1._wp - zcv) * zbv )
[3]	724
	725	! monotonic flux in the k direction, i.e. pcc
	726	! -------------------------------------------
[785]	727	za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
	728	zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
	729	zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
[4990]	730	pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1._wp - zc) * zb )
[3]	731	END DO
	732	END DO
	733	END DO
[2528]	734	CALL lbc_lnk( paa, 'U', -1. ) ; CALL lbc_lnk( pbb, 'V', -1. ) ! lateral boundary condition (changed sign)
[503]	735	!
[3294]	736	CALL wrk_dealloc( jpi, jpj, jpk, zbetup, zbetdo, zbup, zbdo )
[2715]	737	!
[3294]	738	IF( nn_timing == 1 ) CALL timing_stop('nonosc')
	739	!
[3]	740	END SUBROUTINE nonosc
	741
[5770]	742
[7646]	743	SUBROUTINE interp_4th_cpt_org( pt_in, pt_out )
[5770]	744	!!----------------------------------------------------------------------
[7646]	745	!! * ROUTINE interp_4th_cpt_org *
[5770]	746	!!
	747	!! ** Purpose : Compute the interpolation of tracer at w-point
	748	!!
	749	!! ** Method : 4th order compact interpolation
	750	!!----------------------------------------------------------------------
	751	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! now tracer fields
	752	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! now tracer field interpolated at w-pts
	753	!
	754	INTEGER :: ji, jj, jk ! dummy loop integers
	755	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
	756	!!----------------------------------------------------------------------
	757
	758	DO jk = 3, jpkm1 !== build the three diagonal matrix ==!
	759	DO jj = 1, jpj
	760	DO ji = 1, jpi
	761	zwd (ji,jj,jk) = 4._wp
	762	zwi (ji,jj,jk) = 1._wp
	763	zws (ji,jj,jk) = 1._wp
	764	zwrm(ji,jj,jk) = 3._wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
	765	!
	766	IF( tmask(ji,jj,jk+1) == 0._wp) THEN ! Switch to second order centered at bottom
	767	zwd (ji,jj,jk) = 1._wp
	768	zwi (ji,jj,jk) = 0._wp
	769	zws (ji,jj,jk) = 0._wp
	770	zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
	771	ENDIF
	772	END DO
	773	END DO
	774	END DO
	775	!
[7646]	776	jk = 2 ! Switch to second order centered at top
	777	DO jj = 1, jpj
	778	DO ji = 1, jpi
[5770]	779	zwd (ji,jj,jk) = 1._wp
	780	zwi (ji,jj,jk) = 0._wp
	781	zws (ji,jj,jk) = 0._wp
	782	zwrm(ji,jj,jk) = 0.5 * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
	783	END DO
	784	END DO
	785	!
	786	! !== tridiagonal solve ==!
	787	DO jj = 1, jpj ! first recurrence
	788	DO ji = 1, jpi
	789	zwt(ji,jj,2) = zwd(ji,jj,2)
	790	END DO
	791	END DO
	792	DO jk = 3, jpkm1
	793	DO jj = 1, jpj
	794	DO ji = 1, jpi
	795	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
	796	END DO
	797	END DO
	798	END DO
	799	!
	800	DO jj = 1, jpj ! second recurrence: Zk = Yk - Ik / Tk-1 Zk-1
	801	DO ji = 1, jpi
	802	pt_out(ji,jj,2) = zwrm(ji,jj,2)
	803	END DO
	804	END DO
	805	DO jk = 3, jpkm1
	806	DO jj = 1, jpj
	807	DO ji = 1, jpi
	808	pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
	809	END DO
	810	END DO
	811	END DO
	812
	813	DO jj = 1, jpj ! third recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
	814	DO ji = 1, jpi
	815	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
	816	END DO
	817	END DO
	818	DO jk = jpk-2, 2, -1
	819	DO jj = 1, jpj
	820	DO ji = 1, jpi
	821	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
	822	END DO
	823	END DO
	824	END DO
	825	!
[7646]	826	END SUBROUTINE interp_4th_cpt_org
	827
	828
	829	SUBROUTINE interp_4th_cpt( pt_in, pt_out )
	830	!!----------------------------------------------------------------------
	831	!! * ROUTINE interp_4th_cpt *
	832	!!
	833	!! ** Purpose : Compute the interpolation of tracer at w-point
	834	!!
	835	!! ** Method : 4th order compact interpolation
	836	!!----------------------------------------------------------------------
	837	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pt_in ! field at t-point
	838	REAL(wp),DIMENSION(jpi,jpj,jpk), INTENT( out) :: pt_out ! field interpolated at w-point
	839	!
	840	INTEGER :: ji, jj, jk ! dummy loop integers
	841	INTEGER :: ikt, ikb ! local integers
	842	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwd, zwi, zws, zwrm, zwt
	843	!!----------------------------------------------------------------------
	844	!
	845	! !== build the three diagonal matrix & the RHS ==!
	846	!
	847	DO jk = 3, jpkm1 ! interior (from jk=3 to jpk-1)
	848	DO jj = 2, jpjm1
	849	DO ji = fs_2, fs_jpim1
	850	zwd (ji,jj,jk) = 3._wp * wmask(ji,jj,jk) + 1._wp ! diagonal
	851	zwi (ji,jj,jk) = wmask(ji,jj,jk) ! lower diagonal
	852	zws (ji,jj,jk) = wmask(ji,jj,jk) ! upper diagonal
	853	zwrm(ji,jj,jk) = 3._wp * wmask(ji,jj,jk) & ! RHS
	854	& * ( pt_in(ji,jj,jk) + pt_in(ji,jj,jk-1) )
	855	END DO
	856	END DO
	857	END DO
	858	!
	859	!!gm
	860	! SELECT CASE( kbc ) !* boundary condition
	861	! CASE( np_NH ) ! Neumann homogeneous at top & bottom
	862	! CASE( np_CEN2 ) ! 2nd order centered at top & bottom
	863	! END SELECT
	864	!!gm
	865	!
	866	DO jj = 2, jpjm1 ! 2nd order centered at top & bottom
	867	DO ji = fs_2, fs_jpim1
	868	ikt = mikt(ji,jj) + 1 ! w-point below the 1st wet point
	869	ikb = mbkt(ji,jj) ! - above the last wet point
	870	!
	871	zwd (ji,jj,ikt) = 1._wp ! top
	872	zwi (ji,jj,ikt) = 0._wp
	873	zws (ji,jj,ikt) = 0._wp
	874	zwrm(ji,jj,ikt) = 0.5_wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
	875	!
	876	zwd (ji,jj,ikb) = 1._wp ! bottom
	877	zwi (ji,jj,ikb) = 0._wp
	878	zws (ji,jj,ikb) = 0._wp
	879	zwrm(ji,jj,ikb) = 0.5_wp * ( pt_in(ji,jj,jk-1) + pt_in(ji,jj,jk) )
	880	END DO
	881	END DO
	882	!
	883	! !== tridiagonal solver ==!
	884	!
	885	DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1
	886	DO ji = fs_2, fs_jpim1
	887	zwt(ji,jj,2) = zwd(ji,jj,2)
	888	END DO
	889	END DO
	890	DO jk = 3, jpkm1
	891	DO jj = 2, jpjm1
	892	DO ji = fs_2, fs_jpim1
	893	zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) /zwt(ji,jj,jk-1)
	894	END DO
	895	END DO
	896	END DO
	897	!
	898	DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1
	899	DO ji = fs_2, fs_jpim1
	900	pt_out(ji,jj,2) = zwrm(ji,jj,2)
	901	END DO
	902	END DO
	903	DO jk = 3, jpkm1
	904	DO jj = 2, jpjm1
	905	DO ji = fs_2, fs_jpim1
	906	pt_out(ji,jj,jk) = zwrm(ji,jj,jk) - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
	907	END DO
	908	END DO
	909	END DO
	910
	911	DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
	912	DO ji = fs_2, fs_jpim1
	913	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
	914	END DO
	915	END DO
	916	DO jk = jpk-2, 2, -1
	917	DO jj = 2, jpjm1
	918	DO ji = fs_2, fs_jpim1
	919	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - zws(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
	920	END DO
	921	END DO
	922	END DO
	923	!
[5770]	924	END SUBROUTINE interp_4th_cpt
[7646]	925
	926
	927	SUBROUTINE tridia_solver( pD, pU, pL, pRHS, pt_out , klev )
	928	!!----------------------------------------------------------------------
	929	!! * ROUTINE tridia_solver *
	930	!!
	931	!! ** Purpose : solve a symmetric 3diagonal system
	932	!!
	933	!! ** Method : solve M.t_out = RHS(t) where M is a tri diagonal matrix ( jpk*jpk )
	934	!!
	935	!! ( D_1 U_1 0 0 0 )( t_1 ) ( RHS_1 )
	936	!! ( L_2 D_2 U_2 0 0 )( t_2 ) ( RHS_2 )
	937	!! ( 0 L_3 D_3 U_3 0 )( t_3 ) = ( RHS_3 )
	938	!! ( ... )( ... ) ( ... )
	939	!! ( 0 0 0 L_k D_k )( t_k ) ( RHS_k )
	940	!!
	941	!! M is decomposed in the product of an upper and lower triangular matrix.
	942	!! The tri-diagonals matrix is given as input 3D arrays: pD, pU, pL
	943	!! (i.e. the Diagonal, the Upper diagonal, and the Lower diagonal).
	944	!! The solution is pta.
	945	!! The 3d array zwt is used as a work space array.
	946	!!----------------------------------------------------------------------
	947	REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pD, pU, PL ! 3-diagonal matrix
	948	REAL(wp),DIMENSION(:,:,:), INTENT(in ) :: pRHS ! Right-Hand-Side
	949	REAL(wp),DIMENSION(:,:,:), INTENT( out) :: pt_out !!gm field at level=F(klev)
	950	INTEGER , INTENT(in ) :: klev ! =1 pt_out at w-level
	951	! ! =0 pt at t-level
	952	INTEGER :: ji, jj, jk ! dummy loop integers
	953	INTEGER :: kstart ! local indices
	954	REAL(wp),DIMENSION(jpi,jpj,jpk) :: zwt ! 3D work array
	955	!!----------------------------------------------------------------------
	956	!
	957	kstart = 1 + klev
	958	!
	959	DO jj = 2, jpjm1 !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1
	960	DO ji = fs_2, fs_jpim1
	961	zwt(ji,jj,kstart) = pD(ji,jj,kstart)
	962	END DO
	963	END DO
	964	DO jk = kstart+1, jpkm1
	965	DO jj = 2, jpjm1
	966	DO ji = fs_2, fs_jpim1
	967	zwt(ji,jj,jk) = pD(ji,jj,jk) - pL(ji,jj,jk) * pU(ji,jj,jk-1) /zwt(ji,jj,jk-1)
	968	END DO
	969	END DO
	970	END DO
	971	!
	972	DO jj = 2, jpjm1 !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1
	973	DO ji = fs_2, fs_jpim1
	974	pt_out(ji,jj,kstart) = pRHS(ji,jj,kstart)
	975	END DO
	976	END DO
	977	DO jk = kstart+1, jpkm1
	978	DO jj = 2, jpjm1
	979	DO ji = fs_2, fs_jpim1
	980	pt_out(ji,jj,jk) = pRHS(ji,jj,jk) - pL(ji,jj,jk) / zwt(ji,jj,jk-1) *pt_out(ji,jj,jk-1)
	981	END DO
	982	END DO
	983	END DO
	984
	985	DO jj = 2, jpjm1 !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk
	986	DO ji = fs_2, fs_jpim1
	987	pt_out(ji,jj,jpkm1) = pt_out(ji,jj,jpkm1) / zwt(ji,jj,jpkm1)
	988	END DO
	989	END DO
	990	DO jk = jpk-2, kstart, -1
	991	DO jj = 2, jpjm1
	992	DO ji = fs_2, fs_jpim1
	993	pt_out(ji,jj,jk) = ( pt_out(ji,jj,jk) - pU(ji,jj,jk) * pt_out(ji,jj,jk+1) ) / zwt(ji,jj,jk)
	994	END DO
	995	END DO
	996	END DO
	997	!
	998	END SUBROUTINE tridia_solver
	999
[3]	1000	!!======================================================================
[5770]	1001	END MODULE traadv_fct

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

Download in other formats: