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_muscl.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90 @ 3211

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

Property svn:keywords set to Id

File size: 16.4 KB

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

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

Download in other formats: