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_r13383_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/traadv_qck.F90 @ 13516

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

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

Download in other formats: