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_ubs.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_ubs.F90 @ 2326

Last change on this file since 2326 was 2287, checked in by smasson, 14 years ago
update licence of all NEMO files...
Property svn:keywords set to `Id`
File size: 18.3 KB

Rev	Line
[503]	1	MODULE traadv_ubs
	2	!!==============================================================================
	3	!! * MODULE traadv_ubs *
	4	!! Ocean active tracers: horizontal & vertical advective trend
	5	!!==============================================================================
[2024]	6	!! History : 1.0 ! 2006-08 (L. Debreu, R. Benshila) Original code
[2104]	7	!! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
[503]	8	!!----------------------------------------------------------------------
	9
	10	!!----------------------------------------------------------------------
	11	!! tra_adv_ubs : update the tracer trend with the horizontal
	12	!! advection trends using a third order biaised scheme
	13	!!----------------------------------------------------------------------
	14	USE oce ! ocean dynamics and active tracers
	15	USE dom_oce ! ocean space and time domain
[2024]	16	USE trdmod_oce ! ocean space and time domain
	17	USE trdtra
[503]	18	USE lib_mpp
	19	USE lbclnk ! ocean lateral boundary condition (or mpp link)
	20	USE in_out_manager ! I/O manager
	21	USE diaptr ! poleward transport diagnostics
	22	USE dynspg_oce ! choice/control of key cpp for surface pressure gradient
[2082]	23	USE trc_oce ! share passive tracers/Ocean variables
[503]	24
	25	IMPLICIT NONE
	26	PRIVATE
	27
	28	PUBLIC tra_adv_ubs ! routine called by traadv module
	29
[2024]	30	LOGICAL :: l_trd ! flag to compute trends or not
[503]	31
	32	!! * Substitutions
	33	# include "domzgr_substitute.h90"
	34	# include "vectopt_loop_substitute.h90"
	35	!!----------------------------------------------------------------------
[2287]	36	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]	37	!! $Id$
[2287]	38	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[503]	39	!!----------------------------------------------------------------------
	40
	41	CONTAINS
	42
[2104]	43	SUBROUTINE tra_adv_ubs ( kt, cdtype, p2dt, pun, pvn, pwn, &
[2082]	44	& ptb, ptn, pta, kjpt )
[503]	45	!!----------------------------------------------------------------------
	46	!! * ROUTINE tra_adv_ubs *
	47	!!
	48	!! ** Purpose : Compute the now trend due to the advection of tracers
	49	!! and add it to the general trend of passive tracer equations.
	50	!!
[519]	51	!! ** Method : The upstream biased third (UBS) is order scheme based
	52	!! on an upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005)
	53	!! It is only used in the horizontal direction.
	54	!! For example the i-component of the advective fluxes are given by :
	55	!! ! e1u e3u un ( mi(Tn) - zltu(i ) ) if un(i) >= 0
	56	!! zwx = ! or
	57	!! ! e1u e3u un ( mi(Tn) - zltu(i+1) ) if un(i) < 0
	58	!! where zltu is the second derivative of the before temperature field:
	59	!! zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
	60	!! This results in a dissipatively dominant (i.e. hyper-diffusive)
	61	!! truncation error. The overall performance of the advection scheme
	62	!! is similar to that reported in (Farrow and Stevens, 1995).
	63	!! For stability reasons, the first term of the fluxes which corresponds
	64	!! to a second order centered scheme is evaluated using the now velocity
	65	!! (centered in time) while the second term which is the diffusive part
	66	!! of the scheme, is evaluated using the before velocity (forward in time).
	67	!! Note that UBS is not positive. Do not use it on passive tracers.
	68	!! On the vertical, the advection is evaluated using a TVD scheme, as
	69	!! the UBS have been found to be too diffusive.
[503]	70	!!
[2034]	71	!! ** Action : - update (pta) with the now advective tracer trends
[519]	72	!!
	73	!! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.
	74	!! Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741.
[503]	75	!!----------------------------------------------------------------------
[2024]	76	USE oce , zwx => ua ! use ua as workspace
	77	USE oce , zwy => va ! use va as workspace
[2104]	78	!!
	79	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	80	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	81	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	82	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
	83	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
	84	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
	85	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
	86	!!
[2024]	87	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[2104]	88	REAL(wp) :: ztra, zbtr, zcoef ! local scalars
	89	REAL(wp) :: zfp_ui, zfm_ui, zcenut ! - -
	90	REAL(wp) :: zfp_vj, zfm_vj, zcenvt ! - -
	91	REAL(wp) :: z2dtt ! - -
	92	REAL(wp) :: ztak, zfp_wk, zfm_wk ! - -
	93	REAL(wp) :: zeeu, zeev, z_hdivn ! - -
	94	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztu, ztv, zltu , zltv ! 3D workspace
	95	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zti, ztw ! - -
[503]	96	!!----------------------------------------------------------------------
	97
[2104]	98	IF( kt == nit000 ) THEN
[503]	99	IF(lwp) WRITE(numout,*)
[2082]	100	IF(lwp) WRITE(numout,*) 'tra_adv_ubs : horizontal UBS advection scheme on ', cdtype
[503]	101	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
	102	!
[2024]	103	l_trd = .FALSE.
	104	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
[503]	105	ENDIF
[2024]	106	!
	107	! ! ===========
	108	DO jn = 1, kjpt ! tracer loop
	109	! ! ===========
	110	! 1. Bottom value : flux set to zero
	111	! ----------------------------------
	112	zltu(:,:,jpk) = 0.e0 ; zltv(:,:,jpk) = 0.e0
[2104]	113	!
[2024]	114	DO jk = 1, jpkm1 ! Horizontal slab
[2104]	115	!
[2024]	116	! Laplacian
[2104]	117	DO jj = 1, jpjm1 ! First derivative (gradient)
[2024]	118	DO ji = 1, fs_jpim1 ! vector opt.
	119	zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
	120	zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
[2034]	121	ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) )
	122	ztv(ji,jj,jk) = zeev * ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
[2024]	123	END DO
[503]	124	END DO
[2104]	125	DO jj = 2, jpjm1 ! Second derivative (divergence)
[2024]	126	DO ji = fs_2, fs_jpim1 ! vector opt.
	127	zcoef = 1. / ( 6. * fse3t(ji,jj,jk) )
	128	zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef
	129	zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef
	130	END DO
[503]	131	END DO
[2104]	132	!
	133	END DO ! End of slab
	134	CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn)
[503]	135
[2024]	136	!
	137	! Horizontal advective fluxes
[2104]	138	DO jk = 1, jpkm1 ! Horizontal slab
[2024]	139	DO jj = 1, jpjm1
	140	DO ji = 1, fs_jpim1 ! vector opt.
	141	! upstream transport
	142	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
	143	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
	144	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
	145	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
	146	! centered scheme
[2034]	147	zcenut = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) )
	148	zcenvt = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) )
[2024]	149	! UBS scheme
	150	zwx(ji,jj,jk) = zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk)
	151	zwy(ji,jj,jk) = zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk)
	152	END DO
[503]	153	END DO
[2104]	154	END DO ! End of slab
[503]	155
[2034]	156	zltu(:,:,:) = pta(:,:,:,jn) ! store pta trends
[503]	157
[2024]	158	! Horizontal advective trends
[503]	159	DO jk = 1, jpkm1
[2024]	160	! Tracer flux divergence at t-point added to the general trend
[503]	161	DO jj = 2, jpjm1
	162	DO ji = fs_2, fs_jpim1 ! vector opt.
[2024]	163	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	164	! horizontal advective
	165	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) &
	166	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) )
	167	! add it to the general tracer trends
[2034]	168	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[503]	169	END DO
	170	END DO
[2104]	171	!
[2024]	172	END DO ! End of slab
[503]	173
[2024]	174	! Horizontal trend used in tra_adv_ztvd subroutine
[2034]	175	zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)
[503]	176
[2024]	177	! 3. Save the horizontal advective trends for diagnostic
	178	! ------------------------------------------------------
	179	! ! trend diagnostics (contribution of upstream fluxes)
	180	IF( l_trd ) THEN
[2083]	181	CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
	182	CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
[2024]	183	END IF
	184	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	185	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
	186	IF( jn == jp_tem ) pht_adv(:) = ptr_vj( zwy(:,:,:) )
	187	IF( jn == jp_sal ) pst_adv(:) = ptr_vj( zwy(:,:,:) )
[503]	188	ENDIF
[2024]	189
	190	! TVD scheme for the vertical direction
	191	! ----------------------
[2034]	192	IF( l_trd ) zltv(:,:,:) = pta(:,:,:,jn) ! store pta if trend diag.
[503]	193
[2024]	194	! Bottom value : flux set to zero
	195	ztw(:,:,jpk) = 0.e0 ; zti(:,:,jpk) = 0.e0
[503]	196
[2024]	197	! Surface value
	198	IF( lk_vvl ) THEN ; ztw(:,:,1) = 0.e0 ! variable volume : flux set to zero
[2034]	199	ELSE ; ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! free constant surface
[2024]	200	ENDIF
	201	! upstream advection with initial mass fluxes & intermediate update
	202	! -------------------------------------------------------------------
	203	! Interior value
	204	DO jk = 2, jpkm1
	205	DO jj = 1, jpj
	206	DO ji = 1, jpi
	207	zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
	208	zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
[2034]	209	ztw(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )
[2024]	210	END DO
	211	END DO
	212	END DO
	213	! update and guess with monotonic sheme
[503]	214	DO jk = 1, jpkm1
[2082]	215	z2dtt = p2dt(jk)
[503]	216	DO jj = 2, jpjm1
	217	DO ji = fs_2, fs_jpim1 ! vector opt.
[2024]	218	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	219	ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
[2034]	220	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztak
	221	zti(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
[503]	222	END DO
	223	END DO
	224	END DO
	225	!
[2024]	226	CALL lbc_lnk( zti, 'T', 1. ) ! Lateral boundary conditions on zti, zsi (unchanged sign)
[503]	227
[2024]	228	! antidiffusive flux : high order minus low order
	229	ztw(:,:,1) = 0.e0 ! Surface value
	230	DO jk = 2, jpkm1 ! Interior value
	231	DO jj = 1, jpj
	232	DO ji = 1, jpi
[2034]	233	ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - ztw(ji,jj,jk)
[2024]	234	END DO
[503]	235	END DO
	236	END DO
[2024]	237	!
[2082]	238	CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt ) ! monotonicity algorithm
[503]	239
[2024]	240	! final trend with corrected fluxes
	241	DO jk = 1, jpkm1
	242	DO jj = 2, jpjm1
	243	DO ji = fs_2, fs_jpim1 ! vector opt.
	244	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	245	! k- vertical advective trends
	246	ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )
	247	! added to the general tracer trends
[2034]	248	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[2024]	249	END DO
[503]	250	END DO
	251	END DO
	252
[2024]	253	! Save the final vertical advective trends
	254	IF( l_trd ) THEN ! vertical advective trend diagnostics
	255	DO jk = 1, jpkm1 ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
	256	DO jj = 2, jpjm1
	257	DO ji = fs_2, fs_jpim1 ! vector opt.
	258	zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	259	z_hdivn = ( pwn(ji,jj,jk) - pwn(ji,jj,jk+1) ) * zbtr
[2034]	260	zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn
[2024]	261	END DO
	262	END DO
[503]	263	END DO
[2083]	264	CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zltv )
[2024]	265	ENDIF
	266	!
	267	ENDDO
[503]	268	!
[2024]	269	END SUBROUTINE tra_adv_ubs
[503]	270
[2104]	271
[2082]	272	SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
[503]	273	!!---------------------------------------------------------------------
	274	!! * ROUTINE nonosc_z *
	275	!!
	276	!! ** Purpose : compute monotonic tracer fluxes from the upstream
	277	!! scheme and the before field by a nonoscillatory algorithm
	278	!!
	279	!! ** Method : ... ???
	280	!! warning : pbef and paft must be masked, but the boundaries
	281	!! conditions on the fluxes are not necessary zalezak (1979)
	282	!! drange (1995) multi-dimensional forward-in-time and upstream-
	283	!! in-space based differencing for fluid
	284	!!----------------------------------------------------------------------
[2104]	285	REAL(wp), INTENT(in ), DIMENSION(jpk) :: p2dt ! vertical profile of tracer time-step
[2024]	286	REAL(wp), DIMENSION (jpi,jpj,jpk) :: pbef ! before field
[503]	287	REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: paft ! after field
	288	REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: pcc ! monotonic flux in the k direction
	289	!!
	290	INTEGER :: ji, jj, jk ! dummy loop indices
	291	INTEGER :: ikm1
	292	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt
	293	REAL(wp), DIMENSION (jpi,jpj,jpk) :: zbetup, zbetdo
	294	!!----------------------------------------------------------------------
	295
	296	zbig = 1.e+40
	297	zrtrn = 1.e-15
	298	zbetup(:,:,:) = 0.e0 ; zbetdo(:,:,:) = 0.e0
	299
	300	! Search local extrema
	301	! --------------------
	302	! large negative value (-zbig) inside land
	303	pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
	304	paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
	305	! search maximum in neighbourhood
	306	DO jk = 1, jpkm1
	307	ikm1 = MAX(jk-1,1)
	308	DO jj = 2, jpjm1
	309	DO ji = fs_2, fs_jpim1 ! vector opt.
	310	zbetup(ji,jj,jk) = MAX( pbef(ji ,jj ,jk ), paft(ji ,jj ,jk ), &
	311	& pbef(ji ,jj ,ikm1), pbef(ji ,jj ,jk+1), &
	312	& paft(ji ,jj ,ikm1), paft(ji ,jj ,jk+1) )
	313	END DO
	314	END DO
	315	END DO
	316	! large positive value (+zbig) inside land
	317	pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
	318	paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
	319	! search minimum in neighbourhood
	320	DO jk = 1, jpkm1
	321	ikm1 = MAX(jk-1,1)
	322	DO jj = 2, jpjm1
	323	DO ji = fs_2, fs_jpim1 ! vector opt.
	324	zbetdo(ji,jj,jk) = MIN( pbef(ji ,jj ,jk ), paft(ji ,jj ,jk ), &
	325	& pbef(ji ,jj ,ikm1), pbef(ji ,jj ,jk+1), &
	326	& paft(ji ,jj ,ikm1), paft(ji ,jj ,jk+1) )
	327	END DO
	328	END DO
	329	END DO
	330
	331	! restore masked values to zero
	332	pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
	333	paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
	334
[2024]	335
[503]	336	! 2. Positive and negative part of fluxes and beta terms
	337	! ------------------------------------------------------
	338
	339	DO jk = 1, jpkm1
[2082]	340	z2dtt = p2dt(jk)
[503]	341	DO jj = 2, jpjm1
	342	DO ji = fs_2, fs_jpim1 ! vector opt.
	343	! positive & negative part of the flux
	344	zpos = MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
	345	zneg = MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
	346	! up & down beta terms
	347	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
	348	zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
	349	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
	350	END DO
	351	END DO
	352	END DO
	353	! monotonic flux in the k direction, i.e. pcc
	354	! -------------------------------------------
	355	DO jk = 2, jpkm1
	356	DO jj = 2, jpjm1
	357	DO ji = fs_2, fs_jpim1 ! vector opt.
	358	za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
	359	zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
	360	zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
	361	pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
	362	END DO
	363	END DO
	364	END DO
	365	!
	366	END SUBROUTINE nonosc_z
	367
	368	!!======================================================================
	369	END MODULE traadv_ubs

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

Download in other formats: