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_muscl2.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90 @ 2776

Last change on this file since 2776 was 2715, checked in by rblod, 13 years ago
First attempt to put dynamic allocation on the trunk
Property svn:keywords set to `Id`
File size: 15.0 KB

Rev	Line
[3]	1	MODULE traadv_muscl2
[2715]	2	!!======================================================================
[3]	3	!! * MODULE traadv_muscl2 *
[2528]	4	!! Ocean tracers: horizontal & vertical advective trend
[2715]	5	!!======================================================================
[2528]	6	!! History : 1.0 ! 2002-06 (G. Madec) from traadv_muscl
	7	!! 3.2 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
[503]	8	!!----------------------------------------------------------------------
[3]	9
	10	!!----------------------------------------------------------------------
	11	!! tra_adv_muscl2 : update the tracer trend with the horizontal
	12	!! and vertical advection trends using MUSCL2 scheme
	13	!!----------------------------------------------------------------------
	14	USE oce ! ocean dynamics and active tracers
	15	USE dom_oce ! ocean space and time domain
[2528]	16	USE trdmod_oce ! tracers trends
	17	USE trdtra ! tracers trends
[3]	18	USE in_out_manager ! I/O manager
[367]	19	USE dynspg_oce ! choice/control of key cpp for surface pressure gradient
[216]	20	USE trabbl ! tracers: bottom boundary layer
[2528]	21	USE lib_mpp ! distribued memory computing
[67]	22	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[132]	23	USE diaptr ! poleward transport diagnostics
[2528]	24	USE trc_oce ! share passive tracers/Ocean variables
[3]	25
[2528]	26
[3]	27	IMPLICIT NONE
	28	PRIVATE
	29
[2528]	30	PUBLIC tra_adv_muscl2 ! routine called by step.F90
[3]	31
[2528]	32	LOGICAL :: l_trd ! flag to compute trends
	33
[3]	34	!! * Substitutions
	35	# include "domzgr_substitute.h90"
	36	# include "vectopt_loop_substitute.h90"
	37	!!----------------------------------------------------------------------
[2528]	38	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]	39	!! $Id$
[2528]	40	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	41	!!----------------------------------------------------------------------
	42	CONTAINS
	43
[2528]	44	SUBROUTINE tra_adv_muscl2( kt, cdtype, p2dt, pun, pvn, pwn, &
	45	& ptb, ptn, pta, kjpt )
[3]	46	!!----------------------------------------------------------------------
	47	!! * ROUTINE tra_adv_muscl2 *
	48	!!
[216]	49	!! ** Purpose : Compute the now trend due to total advection of T and
	50	!! S using a MUSCL scheme (Monotone Upstream-centered Scheme for
	51	!! Conservation Laws) and add it to the general tracer trend.
[3]	52	!!
[216]	53	!! ** Method : MUSCL scheme plus centered scheme at ocean boundaries
[3]	54	!!
[2528]	55	!! ** Action : - update (pta) with the now advective tracer trends
	56	!! - save trends
[3]	57	!!
[503]	58	!! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
	59	!! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
	60	!!----------------------------------------------------------------------
[2715]	61	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	62	USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as 3D workspace
	63	USE wrk_nemo, ONLY: zslpx => wrk_3d_1 , zslpy => wrk_3d_2 ! 3D workspace
[3]	64	!!
[2528]	65	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	66	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	67	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	68	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
	69	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
	70	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before & now tracer fields
	71	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[503]	72	!!
[2528]	73	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[2715]	74	REAL(wp) :: zu, z0u, zzwx, zw ! local scalars
	75	REAL(wp) :: zv, z0v, zzwy, z0w ! - -
	76	REAL(wp) :: ztra, zbtr, zdt, zalpha ! - -
[3]	77	!!----------------------------------------------------------------------
	78
[2715]	79	IF( wrk_in_use(3, 1,2) ) THEN
	80	CALL ctl_stop('tra_adv_muscl2: requested workspace arrays are unavailable') ; RETURN
	81	ENDIF
	82
[2528]	83	IF( kt == nit000 ) THEN
	84	IF(lwp) WRITE(numout,*)
	85	IF(lwp) WRITE(numout,*) 'tra_adv_muscl2 : MUSCL2 advection scheme on ', cdtype
	86	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
	87	!
	88	l_trd = .FALSE.
[2715]	89	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
[3]	90	ENDIF
	91
[2528]	92	! ! ===========
	93	DO jn = 1, kjpt ! tracer loop
	94	! ! ===========
	95	! I. Horizontal advective fluxes
	96	! ------------------------------
	97	! first guess of the slopes
	98	zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0 ! bottom values
	99	! interior values
	100	DO jk = 1, jpkm1
	101	DO jj = 1, jpjm1
	102	DO ji = 1, fs_jpim1 ! vector opt.
	103	zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk,jn) - ptb(ji,jj,jk,jn) )
	104	zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
	105	END DO
	106	END DO
[3]	107	END DO
[2528]	108	!
	109	CALL lbc_lnk( zwx, 'U', -1. ) ! lateral boundary conditions on zwx, zwy (changed sign)
	110	CALL lbc_lnk( zwy, 'V', -1. )
	111	! !-- Slopes of tracer
	112	zslpx(:,:,jpk) = 0.e0 ; zslpy(:,:,jpk) = 0.e0 ! bottom values
	113	DO jk = 1, jpkm1 ! interior values
	114	DO jj = 2, jpj
	115	DO ji = fs_2, jpi ! vector opt.
	116	zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji-1,jj ,jk) ) &
	117	& * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj ,jk) ) )
	118	zslpy(ji,jj,jk) = ( zwy(ji,jj,jk) + zwy(ji ,jj-1,jk) ) &
	119	& * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji ,jj-1,jk) ) )
	120	END DO
[3]	121	END DO
	122	END DO
[2528]	123	!
	124	DO jk = 1, jpkm1 ! Slopes limitation
	125	DO jj = 2, jpj
	126	DO ji = fs_2, jpi ! vector opt.
	127	zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji ,jj,jk) ), &
	128	& 2.*ABS( zwx (ji-1,jj,jk) ), &
	129	& 2.*ABS( zwx (ji ,jj,jk) ) )
	130	zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN( ABS( zslpy(ji,jj ,jk) ), &
	131	& 2.*ABS( zwy (ji,jj-1,jk) ), &
	132	& 2.*ABS( zwy (ji,jj ,jk) ) )
	133	END DO
	134	END DO
	135	END DO ! interior values
[3]	136
[2528]	137	! !-- MUSCL horizontal advective fluxes
	138	DO jk = 1, jpkm1 ! interior values
	139	zdt = p2dt(jk)
	140	DO jj = 2, jpjm1
	141	DO ji = fs_2, fs_jpim1 ! vector opt.
	142	! MUSCL fluxes
	143	z0u = SIGN( 0.5, pun(ji,jj,jk) )
	144	zalpha = 0.5 - z0u
	145	zu = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
	146	zzwx = ptb(ji+1,jj,jk,jn) + zu * zslpx(ji+1,jj,jk)
	147	zzwy = ptb(ji ,jj,jk,jn) + zu * zslpx(ji ,jj,jk)
	148	zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
	149	!
	150	z0v = SIGN( 0.5, pvn(ji,jj,jk) )
	151	zalpha = 0.5 - z0v
	152	zv = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
	153	zzwx = ptb(ji,jj+1,jk,jn) + zv * zslpy(ji,jj+1,jk)
	154	zzwy = ptb(ji,jj ,jk,jn) + zv * zslpy(ji,jj ,jk)
	155	zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
	156	END DO
[3]	157	END DO
	158	END DO
	159
[2528]	160	!! centered scheme at lateral b.C. if off-shore velocity
[503]	161	DO jk = 1, jpkm1
	162	DO jj = 2, jpjm1
	163	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	164	IF( umask(ji,jj,jk) == 0. ) THEN
	165	IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
	166	zwx(ji+1,jj,jk) = 0.5 * pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk,jn) + ptn(ji+2,jj,jk,jn) )
	167	ENDIF
	168	IF( pun(ji-1,jj,jk) < 0. ) THEN
	169	zwx(ji-1,jj,jk) = 0.5 * pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk,jn) + ptn(ji,jj,jk,jn) )
	170	ENDIF
	171	ENDIF
	172	IF( vmask(ji,jj,jk) == 0. ) THEN
	173	IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
	174	zwy(ji,jj+1,jk) = 0.5 * pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk,jn) + ptn(ji,jj+2,jk,jn) )
	175	ENDIF
	176	IF( pvn(ji,jj-1,jk) < 0. ) THEN
	177	zwy(ji,jj-1,jk) = 0.5 * pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk,jn) + ptn(ji,jj,jk,jn) )
	178	ENDIF
	179	ENDIF
[503]	180	END DO
	181	END DO
	182	END DO
[2528]	183	CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! lateral boundary condition (changed sign)
[216]	184
[2528]	185	! Tracer flux divergence at t-point added to the general trend
[503]	186	DO jk = 1, jpkm1
	187	DO jj = 2, jpjm1
	188	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	189	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	190	! horizontal advective trends
	191	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	192	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) )
	193	! added to the general tracer trends
	194	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[503]	195	END DO
[2528]	196	END DO
[503]	197	END DO
[2528]	198	! ! trend diagnostics (contribution of upstream fluxes)
	199	IF( l_trd ) THEN
	200	CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptb(:,:,:,jn) )
	201	CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptb(:,:,:,jn) )
	202	END IF
[216]	203
[2528]	204	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	205	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
	206	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) )
	207	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) )
	208	ENDIF
[216]	209
[2528]	210	! II. Vertical advective fluxes
	211	! -----------------------------
	212	! !-- first guess of the slopes
	213	zwx (:,:, 1 ) = 0.e0 ; zwx (:,:,jpk) = 0.e0 ! surface & bottom boundary conditions
	214	DO jk = 2, jpkm1 ! interior values
	215	zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) )
	216	END DO
[97]	217
[2528]	218	! !-- Slopes of tracer
	219	zslpx(:,:,1) = 0.e0 ! surface values
	220	DO jk = 2, jpkm1 ! interior value
	221	DO jj = 1, jpj
	222	DO ji = 1, jpi
	223	zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) ) &
	224	& * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )
[3]	225	END DO
	226	END DO
	227	END DO
[2528]	228	! !-- Slopes limitation
	229	DO jk = 2, jpkm1 ! interior values
	230	DO jj = 1, jpj
	231	DO ji = 1, jpi
	232	zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji,jj,jk ) ), &
	233	& 2.*ABS( zwx (ji,jj,jk+1) ), &
	234	& 2.*ABS( zwx (ji,jj,jk ) ) )
	235	END DO
[3]	236	END DO
	237	END DO
[2528]	238	! !-- vertical advective flux
	239	! ! surface values (bottom already set to zero)
	240	IF( lk_vvl ) THEN ; zwx(:,:, 1 ) = 0.e0 ! variable volume
	241	ELSE ; zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface
	242	ENDIF
	243	!
	244	DO jk = 1, jpkm1 ! interior values
	245	zdt = p2dt(jk)
	246	DO jj = 2, jpjm1
	247	DO ji = fs_2, fs_jpim1 ! vector opt.
	248	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3w(ji,jj,jk+1) )
	249	z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
	250	zalpha = 0.5 + z0w
	251	zw = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr
	252	zzwx = ptb(ji,jj,jk+1,jn) + zw * zslpx(ji,jj,jk+1)
	253	zzwy = ptb(ji,jj,jk ,jn) + zw * zslpx(ji,jj,jk )
	254	zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
	255	END DO
[3]	256	END DO
	257	END DO
[2528]	258	!
	259	DO jk = 2, jpkm1 ! centered near the bottom
	260	DO jj = 2, jpjm1
	261	DO ji = fs_2, fs_jpim1 ! vector opt.
	262	IF( tmask(ji,jj,jk+1) == 0. ) THEN
	263	IF( pwn(ji,jj,jk) > 0. ) THEN
	264	zwx(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
	265	ENDIF
[3]	266	ENDIF
[2528]	267	END DO
[3]	268	END DO
	269	END DO
[2528]	270	!
	271	DO jk = 1, jpkm1 ! Compute & add the vertical advective trend
	272	DO jj = 2, jpjm1
[503]	273	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	274	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	275	! vertical advective trends
	276	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
	277	! added to the general tracer trends
	278	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[503]	279	END DO
	280	END DO
	281	END DO
[2528]	282	! ! trend diagnostics (contribution of upstream fluxes)
	283	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwx, pwn, ptb(:,:,:,jn) )
[503]	284	!
[2528]	285	END DO
[503]	286	!
[2715]	287	IF( wrk_not_released(3, 1,2) ) CALL ctl_stop('tra_adv_muscl2: failed to release workspace arrays')
	288	!
[3]	289	END SUBROUTINE tra_adv_muscl2
	290
	291	!!======================================================================
	292	END MODULE traadv_muscl2

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

Download in other formats: