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/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2014/dev_r4650_UKMO14.12_STAND_ALONE_OBSOPER/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl2.F90 @ 5600

Last change on this file since 5600 was 5600, checked in by andrewryan, 9 years ago
merged in latest version of trunk alongside changes to SAO_SRC to be compatible with latest OBS
Property svn:keywords set to `Id`
File size: 15.1 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
[5034]	15	USE trc_oce ! share passive tracers/Ocean variables
[3]	16	USE dom_oce ! ocean space and time domain
[5034]	17	USE trd_oce ! trends: ocean variables
	18	USE trdtra ! trends manager: tracers
[3]	19	USE in_out_manager ! I/O manager
[367]	20	USE dynspg_oce ! choice/control of key cpp for surface pressure gradient
[5034]	21	USE diaptr ! poleward transport diagnostics
	22	!
[2528]	23	USE lib_mpp ! distribued memory computing
[67]	24	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[3294]	25	USE wrk_nemo ! Memory Allocation
	26	USE timing ! Timing
[3625]	27	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[3]	28
	29	IMPLICIT NONE
	30	PRIVATE
	31
[2528]	32	PUBLIC tra_adv_muscl2 ! routine called by step.F90
[3]	33
	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
[3294]	44	SUBROUTINE tra_adv_muscl2( kt, kit000, cdtype, p2dt, pun, pvn, pwn, &
[2528]	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	!!----------------------------------------------------------------------
[2528]	61	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[3294]	62	INTEGER , INTENT(in ) :: kit000 ! first time step index
[2528]	63	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	64	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	65	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
	66	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
	67	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before & now tracer fields
	68	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[503]	69	!!
[2528]	70	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[2715]	71	REAL(wp) :: zu, z0u, zzwx, zw ! local scalars
	72	REAL(wp) :: zv, z0v, zzwy, z0w ! - -
	73	REAL(wp) :: ztra, zbtr, zdt, zalpha ! - -
[5034]	74	REAL(wp), POINTER, DIMENSION(:,:,:) :: zslpx, zslpy , zwx, zwy
[3]	75	!!----------------------------------------------------------------------
[3294]	76	!
	77	IF( nn_timing == 1 ) CALL timing_start('tra_adv_muscl2')
	78	!
[5034]	79	CALL wrk_alloc( jpi, jpj, jpk, zslpx, zslpy, zwx, zwy )
[3294]	80	!
[3]	81
[3294]	82	IF( kt == kit000 ) THEN
[2528]	83	IF(lwp) WRITE(numout,*)
	84	IF(lwp) WRITE(numout,*) 'tra_adv_muscl2 : MUSCL2 advection scheme on ', cdtype
	85	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
[3]	86	ENDIF
[4499]	87	!
[2528]	88	! ! ===========
	89	DO jn = 1, kjpt ! tracer loop
	90	! ! ===========
	91	! I. Horizontal advective fluxes
	92	! ------------------------------
	93	! first guess of the slopes
	94	zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0 ! bottom values
	95	! interior values
	96	DO jk = 1, jpkm1
	97	DO jj = 1, jpjm1
	98	DO ji = 1, fs_jpim1 ! vector opt.
	99	zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk,jn) - ptb(ji,jj,jk,jn) )
	100	zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
	101	END DO
	102	END DO
[3]	103	END DO
[2528]	104	!
	105	CALL lbc_lnk( zwx, 'U', -1. ) ! lateral boundary conditions on zwx, zwy (changed sign)
	106	CALL lbc_lnk( zwy, 'V', -1. )
	107	! !-- Slopes of tracer
	108	zslpx(:,:,jpk) = 0.e0 ; zslpy(:,:,jpk) = 0.e0 ! bottom values
	109	DO jk = 1, jpkm1 ! interior values
	110	DO jj = 2, jpj
	111	DO ji = fs_2, jpi ! vector opt.
	112	zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji-1,jj ,jk) ) &
	113	& * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj ,jk) ) )
	114	zslpy(ji,jj,jk) = ( zwy(ji,jj,jk) + zwy(ji ,jj-1,jk) ) &
	115	& * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji ,jj-1,jk) ) )
	116	END DO
[3]	117	END DO
	118	END DO
[2528]	119	!
	120	DO jk = 1, jpkm1 ! Slopes limitation
	121	DO jj = 2, jpj
	122	DO ji = fs_2, jpi ! vector opt.
	123	zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji ,jj,jk) ), &
	124	& 2.*ABS( zwx (ji-1,jj,jk) ), &
	125	& 2.*ABS( zwx (ji ,jj,jk) ) )
	126	zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN( ABS( zslpy(ji,jj ,jk) ), &
	127	& 2.*ABS( zwy (ji,jj-1,jk) ), &
	128	& 2.*ABS( zwy (ji,jj ,jk) ) )
	129	END DO
	130	END DO
	131	END DO ! interior values
[3]	132
[2528]	133	! !-- MUSCL horizontal advective fluxes
	134	DO jk = 1, jpkm1 ! interior values
	135	zdt = p2dt(jk)
	136	DO jj = 2, jpjm1
	137	DO ji = fs_2, fs_jpim1 ! vector opt.
	138	! MUSCL fluxes
	139	z0u = SIGN( 0.5, pun(ji,jj,jk) )
	140	zalpha = 0.5 - z0u
	141	zu = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
	142	zzwx = ptb(ji+1,jj,jk,jn) + zu * zslpx(ji+1,jj,jk)
	143	zzwy = ptb(ji ,jj,jk,jn) + zu * zslpx(ji ,jj,jk)
	144	zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
	145	!
	146	z0v = SIGN( 0.5, pvn(ji,jj,jk) )
	147	zalpha = 0.5 - z0v
	148	zv = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
	149	zzwx = ptb(ji,jj+1,jk,jn) + zv * zslpy(ji,jj+1,jk)
	150	zzwy = ptb(ji,jj ,jk,jn) + zv * zslpy(ji,jj ,jk)
	151	zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
	152	END DO
[3]	153	END DO
	154	END DO
	155
[2528]	156	!! centered scheme at lateral b.C. if off-shore velocity
[503]	157	DO jk = 1, jpkm1
	158	DO jj = 2, jpjm1
	159	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	160	IF( umask(ji,jj,jk) == 0. ) THEN
	161	IF( pun(ji+1,jj,jk) > 0. .AND. ji /= jpi ) THEN
	162	zwx(ji+1,jj,jk) = 0.5 * pun(ji+1,jj,jk) * ( ptn(ji+1,jj,jk,jn) + ptn(ji+2,jj,jk,jn) )
	163	ENDIF
	164	IF( pun(ji-1,jj,jk) < 0. ) THEN
	165	zwx(ji-1,jj,jk) = 0.5 * pun(ji-1,jj,jk) * ( ptn(ji-1,jj,jk,jn) + ptn(ji,jj,jk,jn) )
	166	ENDIF
	167	ENDIF
	168	IF( vmask(ji,jj,jk) == 0. ) THEN
	169	IF( pvn(ji,jj+1,jk) > 0. .AND. jj /= jpj ) THEN
	170	zwy(ji,jj+1,jk) = 0.5 * pvn(ji,jj+1,jk) * ( ptn(ji,jj+1,jk,jn) + ptn(ji,jj+2,jk,jn) )
	171	ENDIF
	172	IF( pvn(ji,jj-1,jk) < 0. ) THEN
	173	zwy(ji,jj-1,jk) = 0.5 * pvn(ji,jj-1,jk) * ( ptn(ji,jj-1,jk,jn) + ptn(ji,jj,jk,jn) )
	174	ENDIF
	175	ENDIF
[503]	176	END DO
	177	END DO
	178	END DO
[2528]	179	CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! lateral boundary condition (changed sign)
[216]	180
[2528]	181	! Tracer flux divergence at t-point added to the general trend
[503]	182	DO jk = 1, jpkm1
	183	DO jj = 2, jpjm1
	184	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	185	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	186	! horizontal advective trends
	187	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	188	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) )
	189	! added to the general tracer trends
	190	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[503]	191	END DO
[2528]	192	END DO
[503]	193	END DO
[2528]	194	! ! trend diagnostics (contribution of upstream fluxes)
[5034]	195	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. &
	196	&( cdtype == 'TRC' .AND. l_trdtrc ) ) THEN
	197	CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) )
	198	CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) )
[2528]	199	END IF
[216]	200
[2528]	201	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
[5600]	202	IF( cdtype == 'TRA' .AND. ln_diaptr ) THEN
	203	IF( jn == jp_tem ) htr_adv(:) = ptr_sj( zwy(:,:,:) )
	204	IF( jn == jp_sal ) str_adv(:) = ptr_sj( zwy(:,:,:) )
[2528]	205	ENDIF
[216]	206
[2528]	207	! II. Vertical advective fluxes
	208	! -----------------------------
	209	! !-- first guess of the slopes
	210	zwx (:,:, 1 ) = 0.e0 ; zwx (:,:,jpk) = 0.e0 ! surface & bottom boundary conditions
	211	DO jk = 2, jpkm1 ! interior values
	212	zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) )
	213	END DO
[97]	214
[2528]	215	! !-- Slopes of tracer
	216	zslpx(:,:,1) = 0.e0 ! surface values
	217	DO jk = 2, jpkm1 ! interior value
	218	DO jj = 1, jpj
	219	DO ji = 1, jpi
	220	zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) ) &
	221	& * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )
[3]	222	END DO
	223	END DO
	224	END DO
[2528]	225	! !-- Slopes limitation
	226	DO jk = 2, jpkm1 ! interior values
	227	DO jj = 1, jpj
	228	DO ji = 1, jpi
	229	zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji,jj,jk ) ), &
	230	& 2.*ABS( zwx (ji,jj,jk+1) ), &
	231	& 2.*ABS( zwx (ji,jj,jk ) ) )
	232	END DO
[3]	233	END DO
	234	END DO
[2528]	235	! !-- vertical advective flux
	236	! ! surface values (bottom already set to zero)
	237	IF( lk_vvl ) THEN ; zwx(:,:, 1 ) = 0.e0 ! variable volume
	238	ELSE ; zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface
	239	ENDIF
	240	!
	241	DO jk = 1, jpkm1 ! interior values
	242	zdt = p2dt(jk)
	243	DO jj = 2, jpjm1
	244	DO ji = fs_2, fs_jpim1 ! vector opt.
	245	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3w(ji,jj,jk+1) )
	246	z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
	247	zalpha = 0.5 + z0w
	248	zw = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr
	249	zzwx = ptb(ji,jj,jk+1,jn) + zw * zslpx(ji,jj,jk+1)
	250	zzwy = ptb(ji,jj,jk ,jn) + zw * zslpx(ji,jj,jk )
	251	zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
	252	END DO
[3]	253	END DO
	254	END DO
[2528]	255	!
	256	DO jk = 2, jpkm1 ! centered near the bottom
	257	DO jj = 2, jpjm1
	258	DO ji = fs_2, fs_jpim1 ! vector opt.
	259	IF( tmask(ji,jj,jk+1) == 0. ) THEN
	260	IF( pwn(ji,jj,jk) > 0. ) THEN
	261	zwx(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk-1,jn) + ptn(ji,jj,jk,jn) )
	262	ENDIF
[3]	263	ENDIF
[2528]	264	END DO
[3]	265	END DO
	266	END DO
[2528]	267	!
	268	DO jk = 1, jpkm1 ! Compute & add the vertical advective trend
	269	DO jj = 2, jpjm1
[503]	270	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	271	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	272	! vertical advective trends
	273	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
	274	! added to the general tracer trends
	275	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[503]	276	END DO
	277	END DO
	278	END DO
[2528]	279	! ! trend diagnostics (contribution of upstream fluxes)
[5034]	280	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. &
	281	&( cdtype == 'TRC' .AND. l_trdtrc ) ) &
	282	CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) )
[503]	283	!
[2528]	284	END DO
[503]	285	!
[5034]	286	CALL wrk_dealloc( jpi, jpj, jpk, zslpx, zslpy, zwx, zwy )
[2715]	287	!
[3294]	288	IF( nn_timing == 1 ) CALL timing_stop('tra_adv_muscl2')
	289	!
[3]	290	END SUBROUTINE tra_adv_muscl2
	291
	292	!!======================================================================
	293	END MODULE traadv_muscl2

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

Download in other formats: