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 branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90 @ 2333

Last change on this file since 2333 was 2333, checked in by rblod, 14 years ago
Add missing lwp in some advection routines
Property svn:keywords set to `Id`
File size: 14.7 KB

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

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

Download in other formats: