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_tvd.F90 in branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/DEV_r2006_merge_TRA_TRC/NEMO/OPA_SRC/TRA/traadv_tvd.F90 @ 2083

Last change on this file since 2083 was 2083, checked in by cetlod, 14 years ago
correction of minor bugs
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 17.7 KB

Rev	Line
[3]	1	MODULE traadv_tvd
	2	!!==============================================================================
	3	!! * MODULE traadv_tvd *
[2024]	4	!! Ocean tracers: horizontal & vertical advective trend
[3]	5	!!==============================================================================
[503]	6	!! History : ! 95-12 (L. Mortier) Original code
	7	!! ! 00-01 (H. Loukos) adapted to ORCA
	8	!! ! 00-10 (MA Foujols E.Kestenare) include file not routine
	9	!! ! 00-12 (E. Kestenare M. Levy) fix bug in trtrd indexes
[2034]	10	!! ! 01-07 (E. Durand G. Madec) adaptation to ORCA config
[503]	11	!! 8.5 ! 02-06 (G. Madec) F90: Free form and module
	12	!! 9.0 ! 04-01 (A. de Miranda, G. Madec, J.M. Molines ): advective bbl
	13	!! 9.0 ! 08-04 (S. Cravatte) add the i-, j- & k- trends computation
[2024]	14	!! " " ! 09-11 (V. Garnier) Surface pressure gradient organization
[2034]	15	!! 3.3 ! 10-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
[503]	16	!!----------------------------------------------------------------------
[3]	17
[503]	18
[3]	19	!!----------------------------------------------------------------------
	20	!! tra_adv_tvd : update the tracer trend with the horizontal
	21	!! and vertical advection trends using a TVD scheme
	22	!! nonosc : compute monotonic tracer fluxes by a nonoscillatory
	23	!! algorithm
	24	!!----------------------------------------------------------------------
	25	USE oce ! ocean dynamics and active tracers
	26	USE dom_oce ! ocean space and time domain
[2024]	27	USE trdmod_oce ! tracers trends
	28	USE trdtra ! tracers trends
[3]	29	USE in_out_manager ! I/O manager
[367]	30	USE dynspg_oce ! choice/control of key cpp for surface pressure gradient
[3]	31	USE lib_mpp
[74]	32	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[132]	33	USE diaptr ! poleward transport diagnostics
[2082]	34	USE trc_oce ! share passive tracers/Ocean variables
[3]	35
[74]	36
[3]	37	IMPLICIT NONE
	38	PRIVATE
	39
[503]	40	PUBLIC tra_adv_tvd ! routine called by step.F90
[3]	41
[2024]	42	LOGICAL :: l_trd ! flag to compute trends
	43
[3]	44	!! * Substitutions
	45	# include "domzgr_substitute.h90"
	46	# include "vectopt_loop_substitute.h90"
	47	!!----------------------------------------------------------------------
[2034]	48	!! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
[1152]	49	!! $Id$
[503]	50	!! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
[3]	51	!!----------------------------------------------------------------------
	52
	53	CONTAINS
	54
[2082]	55	SUBROUTINE tra_adv_tvd ( kt, cdtype, p2dt, pun, pvn, pwn, &
	56	& ptb, ptn, pta, kjpt )
[3]	57	!!----------------------------------------------------------------------
	58	!! * ROUTINE tra_adv_tvd *
	59	!!
	60	!! ** Purpose : Compute the now trend due to total advection of
	61	!! tracers and add it to the general trend of tracer equations
	62	!!
	63	!! ** Method : TVD scheme, i.e. 2nd order centered scheme with
	64	!! corrected flux (monotonic correction)
	65	!! note: - this advection scheme needs a leap-frog time scheme
	66	!!
[2034]	67	!! ** Action : - update (pta) with the now advective tracer trends
[2024]	68	!! - save the trends
[503]	69	!!----------------------------------------------------------------------
[2024]	70	USE oce , zwx => ua ! use ua as workspace
	71	USE oce , zwy => va ! use va as workspace
[2034]	72	!!
[2024]	73	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	74	CHARACTER(len=3), INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	75	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[2082]	76	REAL(wp) , INTENT(in ), DIMENSION(jpk) :: p2dt ! vertical profile of tracer time-step
[2024]	77	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk) :: pun, pvn, pwn ! 3 ocean velocity components
[2034]	78	REAL(wp) , INTENT(in ), DIMENSION(jpi,jpj,jpk,kjpt) :: ptb, ptn ! before and now tracer fields
	79	REAL(wp) , INTENT(inout), DIMENSION(jpi,jpj,jpk,kjpt) :: pta ! tracer trend
	80	!!
[2024]	81	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[3]	82	REAL(wp) :: &
[2082]	83	z2dtt, zbtr, ztra, & ! temporary scalar
[2024]	84	zfp_ui, zfp_vj, zfp_wk, & ! " "
	85	zfm_ui, zfm_vj, zfm_wk ! " "
	86	REAL(wp), DIMENSION (jpi,jpj,jpk) :: zwi, zwz ! temporary workspace
	87	REAL(wp), DIMENSION (:,:,:), ALLOCATABLE :: ztrdx, ztrdy, ztrdz
[3]	88	!!----------------------------------------------------------------------
	89
[408]	90
[2082]	91	IF( ( cdtype == 'TRA' .AND. kt == nit000 ) .OR. ( cdtype == 'TRC' .AND. kt == nittrc000 ) ) THEN
[3]	92	WRITE(numout,*)
[2082]	93	WRITE(numout,*) 'tra_adv_tvd : TVD advection scheme on ', cdtype
[3]	94	WRITE(numout,*) '~~~~~~~~~~~'
[2024]	95	!
	96	l_trd = .FALSE.
	97	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
[3]	98	ENDIF
[2024]	99	!
	100	IF( l_trd ) THEN
	101	ALLOCATE( ztrdx(jpi,jpj,jpk) ) ; ztrdx(:,:,:) = 0.
	102	ALLOCATE( ztrdy(jpi,jpj,jpk) ) ; ztrdy(:,:,:) = 0.
	103	ALLOCATE( ztrdz(jpi,jpj,jpk) ) ; ztrdz(:,:,:) = 0.
	104	END IF
	105	!
[2082]	106	zwi(:,:,:) = 0.e0
[2024]	107	!
	108	! ! ===========
	109	DO jn = 1, kjpt ! tracer loop
	110	! ! ===========
	111	! 1. Bottom value : flux set to zero
	112	! ----------------------------------
	113	zwx(:,:,jpk) = 0.e0 ; zwz(:,:,jpk) = 0.e0
	114	zwy(:,:,jpk) = 0.e0 ; zwi(:,:,jpk) = 0.e0
[3]	115
[2024]	116	! 2. upstream advection with initial mass fluxes & intermediate update
	117	! --------------------------------------------------------------------
	118	! upstream tracer flux in the i and j direction
	119	DO jk = 1, jpkm1
	120	DO jj = 1, jpjm1
	121	DO ji = 1, fs_jpim1 ! vector opt.
	122	! upstream scheme
	123	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
	124	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
	125	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
	126	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
[2034]	127	zwx(ji,jj,jk) = 0.5 * ( zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn) )
	128	zwy(ji,jj,jk) = 0.5 * ( zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn) )
[2024]	129	END DO
[3]	130	END DO
	131	END DO
	132
[2024]	133	! upstream tracer flux in the k direction
	134	! Surface value
	135	IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! volume variable
[2034]	136	ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface
[2024]	137	ENDIF
	138	! Interior value
	139	DO jk = 2, jpkm1
	140	DO jj = 1, jpj
	141	DO ji = 1, jpi
	142	zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
	143	zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
[2034]	144	zwz(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )
[2024]	145	END DO
[3]	146	END DO
	147	END DO
	148
[2024]	149	! total advective trend
[216]	150	DO jk = 1, jpkm1
[2082]	151	z2dtt = p2dt(jk)
[216]	152	DO jj = 2, jpjm1
	153	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	154	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
[2024]	155	! total intermediate advective trends
	156	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	157	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
	158	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
	159	! update and guess with monotonic sheme
[2034]	160	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
	161	zwi(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ztra ) * tmask(ji,jj,jk)
[216]	162	END DO
	163	END DO
	164	END DO
[2024]	165	! ! Lateral boundary conditions on zwi (unchanged sign)
	166	CALL lbc_lnk( zwi, 'T', 1. )
	167
	168	! ! trend diagnostics (contribution of upstream fluxes)
	169	IF( l_trd ) THEN
	170	! store intermediate advective trends
	171	ztrdx(:,:,:) = zwx(:,:,:) ; ztrdy(:,:,:) = zwy(:,:,:) ; ztrdz(:,:,:) = zwz(:,:,:)
	172	END IF
	173	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	174	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
	175	IF( jn == jp_tem ) pht_adv(:) = ptr_vj( zwy(:,:,:) )
	176	IF( jn == jp_sal ) pst_adv(:) = ptr_vj( zwy(:,:,:) )
	177	ENDIF
	178
	179	! 3. antidiffusive flux : high order minus low order
	180	! --------------------------------------------------
	181	! antidiffusive flux on i and j
[503]	182	DO jk = 1, jpkm1
[2024]	183	DO jj = 1, jpjm1
	184	DO ji = 1, fs_jpim1 ! vector opt.
[2034]	185	zwx(ji,jj,jk) = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj,jk,jn) ) - zwx(ji,jj,jk)
	186	zwy(ji,jj,jk) = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj+1,jk,jn) ) - zwy(ji,jj,jk)
[503]	187	END DO
	188	END DO
	189	END DO
[2024]	190
	191	! antidiffusive flux on k
	192	! Surface value
	193	zwz(:,:,1) = 0.e0
	194	! Interior value
	195	DO jk = 2, jpkm1
	196	DO jj = 1, jpj
	197	DO ji = 1, jpi
[2034]	198	zwz(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - zwz(ji,jj,jk)
[503]	199	END DO
	200	END DO
	201	END DO
[216]	202
[2024]	203	! Lateral bondary conditions
	204	CALL lbc_lnk( zwx, 'U', -1. )
	205	CALL lbc_lnk( zwy, 'V', -1. )
	206	CALL lbc_lnk( zwz, 'W', 1. )
[3]	207
[2024]	208	! 4. monotonicity algorithm
	209	! -------------------------
[2082]	210	CALL nonosc( ptb(:,:,:,jn), zwx, zwy, zwz, zwi, p2dt )
[3]	211
	212
[2024]	213	! 5. final trend with corrected fluxes
	214	! ------------------------------------
[216]	215	DO jk = 1, jpkm1
	216	DO jj = 2, jpjm1
[2024]	217	DO ji = fs_2, fs_jpim1 ! vector opt.
[503]	218	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
[2024]	219	! total advective trends
	220	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	221	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
	222	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
	223	! add them to the general tracer trends
[2034]	224	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[216]	225	END DO
	226	END DO
	227	END DO
[2024]	228
	229	! ! trend diagnostics (contribution of upstream fluxes)
	230	IF( l_trd ) THEN
	231	ztrdx(:,:,:) = ztrdx(:,:,:) + zwx(:,:,:) ! <<< Add to previously computed
	232	ztrdy(:,:,:) = ztrdy(:,:,:) + zwy(:,:,:) ! <<< Add to previously computed
	233	ztrdz(:,:,:) = ztrdz(:,:,:) + zwz(:,:,:) ! <<< Add to previously computed
	234
[2083]	235	CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, ztrdx, pun, ptn(:,:,:,jn) )
	236	CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, ztrdy, pvn, ptn(:,:,:,jn) )
	237	CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, ztrdz, pwn, ptn(:,:,:,jn) )
[2024]	238	END IF
	239	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	240	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
	241	IF( jn == jp_tem ) pht_adv(:) = ptr_vj( zwy(:,:,:) ) + pht_adv(:)
	242	IF( jn == jp_sal ) pst_adv(:) = ptr_vj( zwy(:,:,:) ) + pst_adv(:)
	243	ENDIF
[503]	244	!
[2024]	245	ENDDO
[503]	246	!
[2024]	247	IF( l_trd ) THEN
	248	DEALLOCATE( ztrdx ) ; DEALLOCATE( ztrdy ) ; DEALLOCATE( ztrdz )
	249	END IF
	250	!
[3]	251	END SUBROUTINE tra_adv_tvd
	252
	253
[2082]	254	SUBROUTINE nonosc( pbef, paa, pbb, pcc, paft, p2dt )
[3]	255	!!---------------------------------------------------------------------
	256	!! * ROUTINE nonosc *
	257	!!
	258	!! ** Purpose : compute monotonic tracer fluxes from the upstream
	259	!! scheme and the before field by a nonoscillatory algorithm
	260	!!
	261	!! ** Method : ... ???
	262	!! warning : pbef and paft must be masked, but the boundaries
	263	!! conditions on the fluxes are not necessary zalezak (1979)
	264	!! drange (1995) multi-dimensional forward-in-time and upstream-
	265	!! in-space based differencing for fluid
	266	!!----------------------------------------------------------------------
[2082]	267	REAL(wp), DIMENSION(jpk) , INTENT( in ) :: &
	268	p2dt ! vertical profile of tracer time-step
[2024]	269	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT( in ) :: &
	270	pbef, & ! before field
	271	paft ! after field
[3]	272	REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT( inout ) :: &
	273	paa, & ! monotonic flux in the i direction
	274	pbb, & ! monotonic flux in the j direction
	275	pcc ! monotonic flux in the k direction
[503]	276	!!
[3]	277	INTEGER :: ji, jj, jk ! dummy loop indices
	278	INTEGER :: ikm1
	279	REAL(wp), DIMENSION (jpi,jpj,jpk) :: zbetup, zbetdo
[785]	280	REAL(wp), DIMENSION (jpi,jpj,jpk) :: zbup, zbdo
[3]	281	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt
[785]	282	REAL(wp) :: zau, zbu, zcu, zav, zbv, zcv
	283	REAL(wp) :: zup, zdo
[3]	284	!!----------------------------------------------------------------------
	285
	286	zbig = 1.e+40
	287	zrtrn = 1.e-15
[785]	288	zbetup(:,:,jpk) = 0.e0 ; zbetdo(:,:,jpk) = 0.e0
[3]	289
[785]	290
[3]	291	! Search local extrema
	292	! --------------------
[785]	293	! max/min of pbef & paft with large negative/positive value (-/+zbig) inside land
	294	zbup = MAX( pbef * tmask - zbig * ( 1.e0 - tmask ), &
	295	& paft * tmask - zbig * ( 1.e0 - tmask ) )
	296	zbdo = MIN( pbef * tmask + zbig * ( 1.e0 - tmask ), &
	297	& paft * tmask + zbig * ( 1.e0 - tmask ) )
	298
[3]	299	DO jk = 1, jpkm1
	300	ikm1 = MAX(jk-1,1)
[2082]	301	z2dtt = p2dt(jk)
[3]	302	DO jj = 2, jpjm1
	303	DO ji = fs_2, fs_jpim1 ! vector opt.
	304
[785]	305	! search maximum in neighbourhood
	306	zup = MAX( zbup(ji ,jj ,jk ), &
	307	& zbup(ji-1,jj ,jk ), zbup(ji+1,jj ,jk ), &
	308	& zbup(ji ,jj-1,jk ), zbup(ji ,jj+1,jk ), &
	309	& zbup(ji ,jj ,ikm1), zbup(ji ,jj ,jk+1) )
[3]	310
[785]	311	! search minimum in neighbourhood
	312	zdo = MIN( zbdo(ji ,jj ,jk ), &
	313	& zbdo(ji-1,jj ,jk ), zbdo(ji+1,jj ,jk ), &
	314	& zbdo(ji ,jj-1,jk ), zbdo(ji ,jj+1,jk ), &
	315	& zbdo(ji ,jj ,ikm1), zbdo(ji ,jj ,jk+1) )
[3]	316
[785]	317	! positive part of the flux
[3]	318	zpos = MAX( 0., paa(ji-1,jj ,jk ) ) - MIN( 0., paa(ji ,jj ,jk ) ) &
	319	& + MAX( 0., pbb(ji ,jj-1,jk ) ) - MIN( 0., pbb(ji ,jj ,jk ) ) &
	320	& + MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
[785]	321
	322	! negative part of the flux
[3]	323	zneg = MAX( 0., paa(ji ,jj ,jk ) ) - MIN( 0., paa(ji-1,jj ,jk ) ) &
	324	& + MAX( 0., pbb(ji ,jj ,jk ) ) - MIN( 0., pbb(ji ,jj-1,jk ) ) &
	325	& + MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
[785]	326
[3]	327	! up & down beta terms
	328	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
[785]	329	zbetup(ji,jj,jk) = ( zup - paft(ji,jj,jk) ) / ( zpos + zrtrn ) * zbt
	330	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zdo ) / ( zneg + zrtrn ) * zbt
[3]	331	END DO
	332	END DO
	333	END DO
	334
	335	! lateral boundary condition on zbetup & zbetdo (unchanged sign)
	336	CALL lbc_lnk( zbetup, 'T', 1. )
	337	CALL lbc_lnk( zbetdo, 'T', 1. )
	338
	339
[237]	340	! 3. monotonic flux in the i & j direction (paa & pbb)
	341	! ----------------------------------------
[3]	342	DO jk = 1, jpkm1
	343	DO jj = 2, jpjm1
	344	DO ji = fs_2, fs_jpim1 ! vector opt.
[785]	345	zau = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji+1,jj,jk) )
	346	zbu = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji+1,jj,jk) )
	347	zcu = ( 0.5 + SIGN( 0.5 , paa(ji,jj,jk) ) )
	348	paa(ji,jj,jk) = paa(ji,jj,jk) * ( zcu * zau + ( 1.e0 - zcu) * zbu )
[3]	349
[785]	350	zav = MIN( 1.e0, zbetdo(ji,jj,jk), zbetup(ji,jj+1,jk) )
	351	zbv = MIN( 1.e0, zbetup(ji,jj,jk), zbetdo(ji,jj+1,jk) )
	352	zcv = ( 0.5 + SIGN( 0.5 , pbb(ji,jj,jk) ) )
	353	pbb(ji,jj,jk) = pbb(ji,jj,jk) * ( zcv * zav + ( 1.e0 - zcv) * zbv )
[3]	354
	355	! monotonic flux in the k direction, i.e. pcc
	356	! -------------------------------------------
[785]	357	za = MIN( 1., zbetdo(ji,jj,jk+1), zbetup(ji,jj,jk) )
	358	zb = MIN( 1., zbetup(ji,jj,jk+1), zbetdo(ji,jj,jk) )
	359	zc = ( 0.5 + SIGN( 0.5 , pcc(ji,jj,jk+1) ) )
	360	pcc(ji,jj,jk+1) = pcc(ji,jj,jk+1) * ( zc * za + ( 1.e0 - zc) * zb )
[3]	361	END DO
	362	END DO
	363	END DO
	364
[237]	365	! lateral boundary condition on paa, pbb, pcc
	366	CALL lbc_lnk( paa, 'U', -1. ) ! changed sign
	367	CALL lbc_lnk( pbb, 'V', -1. ) ! changed sign
[503]	368	!
[3]	369	END SUBROUTINE nonosc
	370
	371	!!======================================================================
	372	END MODULE traadv_tvd

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

Download in other formats: