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/2020/dev_r13898_Tiling_Cleanup_MPI3/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13898_Tiling_Cleanup_MPI3/src/OCE/TRA/traadv_qck.F90 @ 13906

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

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

Download in other formats: