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

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_muscl.F90 @ 2804

Last change on this file since 2804 was 2715, checked in by rblod, 13 years ago
First attempt to put dynamic allocation on the trunk
Property svn:keywords set to `Id`
File size: 13.6 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
[3]	36	!! * Substitutions
	37	# include "domzgr_substitute.h90"
	38	# include "vectopt_loop_substitute.h90"
	39	!!----------------------------------------------------------------------
[2528]	40	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]	41	!! $Id$
[2528]	42	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	43	!!----------------------------------------------------------------------
	44	CONTAINS
	45
[2528]	46	SUBROUTINE tra_adv_muscl( kt, cdtype, p2dt, pun, pvn, pwn, &
	47	& ptb, pta, kjpt )
[3]	48	!!----------------------------------------------------------------------
	49	!! * ROUTINE tra_adv_muscl *
[216]	50	!!
[3]	51	!! ** Purpose : Compute the now trend due to total advection of T and
	52	!! S using a MUSCL scheme (Monotone Upstream-centered Scheme for
	53	!! Conservation Laws) and add it to the general tracer trend.
	54	!!
[216]	55	!! ** Method : MUSCL scheme plus centered scheme at ocean boundaries
[3]	56	!!
	57	!! ** Action : - update (ta,sa) with the now advective tracer trends
[2528]	58	!! - save trends
[3]	59	!!
[503]	60	!! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
	61	!! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
	62	!!----------------------------------------------------------------------
[2715]	63	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	64	USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace
	65	USE wrk_nemo, ONLY: zslpx => wrk_3d_1 , zslpy => wrk_3d_2 ! 3D workspace
	66	!
[2528]	67	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	68	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	69	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	70	REAL(wp), DIMENSION( jpk ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
	71	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
	72	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before tracer field
	73	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
[2715]	74	!
[2528]	75	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[2715]	76	REAL(wp) :: zu, z0u, zzwx, zw ! local scalars
	77	REAL(wp) :: zv, z0v, zzwy, z0w ! - -
	78	REAL(wp) :: ztra, zbtr, zdt, zalpha ! - -
[3]	79	!!----------------------------------------------------------------------
	80
[2715]	81	IF( wrk_in_use(3, 1,2) ) THEN
	82	CALL ctl_stop('tra_adv_muscl: requested workspace arrays unavailable') ; RETURN
	83	ENDIF
	84
[2528]	85	IF( kt == nit000 ) THEN
	86	IF(lwp) WRITE(numout,*)
	87	IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
	88	IF(lwp) WRITE(numout,*) '~~~~~~~'
	89	!
	90	l_trd = .FALSE.
	91	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
[3]	92	ENDIF
	93
[2528]	94	! ! ===========
	95	DO jn = 1, kjpt ! tracer loop
	96	! ! ===========
	97	! I. Horizontal advective fluxes
	98	! ------------------------------
	99	! first guess of the slopes
	100	zwx(:,:,jpk) = 0.e0 ; zwy(:,:,jpk) = 0.e0 ! bottom values
	101	! interior values
	102	DO jk = 1, jpkm1
	103	DO jj = 1, jpjm1
	104	DO ji = 1, fs_jpim1 ! vector opt.
	105	zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk,jn) - ptb(ji,jj,jk,jn) )
	106	zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
	107	END DO
	108	END DO
[3]	109	END DO
[2528]	110	!
	111	CALL lbc_lnk( zwx, 'U', -1. ) ! lateral boundary conditions on zwx, zwy (changed sign)
	112	CALL lbc_lnk( zwy, 'V', -1. )
	113	! !-- Slopes of tracer
	114	zslpx(:,:,jpk) = 0.e0 ; zslpy(:,:,jpk) = 0.e0 ! bottom values
	115	DO jk = 1, jpkm1 ! interior values
	116	DO jj = 2, jpj
	117	DO ji = fs_2, jpi ! vector opt.
	118	zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji-1,jj ,jk) ) &
	119	& * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj ,jk) ) )
	120	zslpy(ji,jj,jk) = ( zwy(ji,jj,jk) + zwy(ji ,jj-1,jk) ) &
	121	& * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji ,jj-1,jk) ) )
	122	END DO
[3]	123	END DO
	124	END DO
[503]	125	!
[2528]	126	DO jk = 1, jpkm1 ! Slopes limitation
	127	DO jj = 2, jpj
	128	DO ji = fs_2, jpi ! vector opt.
	129	zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji ,jj,jk) ), &
	130	& 2.*ABS( zwx (ji-1,jj,jk) ), &
	131	& 2.*ABS( zwx (ji ,jj,jk) ) )
	132	zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN( ABS( zslpy(ji,jj ,jk) ), &
	133	& 2.*ABS( zwy (ji,jj-1,jk) ), &
	134	& 2.*ABS( zwy (ji,jj ,jk) ) )
[503]	135	END DO
[2528]	136	END DO
	137	END DO ! interior values
[216]	138
[2528]	139	! !-- MUSCL horizontal advective fluxes
	140	DO jk = 1, jpkm1 ! interior values
	141	zdt = p2dt(jk)
[503]	142	DO jj = 2, jpjm1
	143	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	144	! MUSCL fluxes
	145	z0u = SIGN( 0.5, pun(ji,jj,jk) )
	146	zalpha = 0.5 - z0u
	147	zu = z0u - 0.5 * pun(ji,jj,jk) * zdt / ( e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk) )
	148	zzwx = ptb(ji+1,jj,jk,jn) + zu * zslpx(ji+1,jj,jk)
	149	zzwy = ptb(ji ,jj,jk,jn) + zu * zslpx(ji ,jj,jk)
	150	zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
	151	!
	152	z0v = SIGN( 0.5, pvn(ji,jj,jk) )
	153	zalpha = 0.5 - z0v
	154	zv = z0v - 0.5 * pvn(ji,jj,jk) * zdt / ( e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk) )
	155	zzwx = ptb(ji,jj+1,jk,jn) + zv * zslpy(ji,jj+1,jk)
	156	zzwy = ptb(ji,jj ,jk,jn) + zv * zslpy(ji,jj ,jk)
	157	zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
[503]	158	END DO
	159	END DO
	160	END DO
[2528]	161	! ! lateral boundary conditions on zwx, zwy (changed sign)
	162	CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. )
[503]	163	!
[2528]	164	! Tracer flux divergence at t-point added to the general trend
	165	DO jk = 1, jpkm1
	166	DO jj = 2, jpjm1
	167	DO ji = fs_2, fs_jpim1 ! vector opt.
	168	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	169	! horizontal advective trends
	170	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
	171	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) )
	172	! add it to the general tracer trends
	173	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[3]	174	END DO
[2528]	175	END DO
	176	END DO
	177	! ! trend diagnostics (contribution of upstream fluxes)
	178	IF( l_trd ) THEN
	179	CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptb(:,:,:,jn) )
	180	CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptb(:,:,:,jn) )
	181	END IF
	182	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	183	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
	184	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) )
	185	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) )
[457]	186	ENDIF
[3]	187
[2528]	188	! II. Vertical advective fluxes
	189	! -----------------------------
	190	! !-- first guess of the slopes
	191	zwx (:,:, 1 ) = 0.e0 ; zwx (:,:,jpk) = 0.e0 ! surface & bottom boundary conditions
	192	DO jk = 2, jpkm1 ! interior values
	193	zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) )
[3]	194	END DO
	195
[2528]	196	! !-- Slopes of tracer
	197	zslpx(:,:,1) = 0.e0 ! surface values
	198	DO jk = 2, jpkm1 ! interior value
	199	DO jj = 1, jpj
	200	DO ji = 1, jpi
	201	zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) ) &
	202	& * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )
	203	END DO
[3]	204	END DO
	205	END DO
[2528]	206	! !-- Slopes limitation
	207	DO jk = 2, jpkm1 ! interior values
	208	DO jj = 1, jpj
	209	DO ji = 1, jpi
	210	zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji,jj,jk ) ), &
	211	& 2.*ABS( zwx (ji,jj,jk+1) ), &
	212	& 2.*ABS( zwx (ji,jj,jk ) ) )
	213	END DO
[3]	214	END DO
	215	END DO
[2528]	216	! !-- vertical advective flux
	217	! ! surface values (bottom already set to zero)
	218	IF( lk_vvl ) THEN ; zwx(:,:, 1 ) = 0.e0 ! variable volume
	219	ELSE ; zwx(:,:, 1 ) = pwn(:,:,1) * ptb(:,:,1,jn) ! linear free surface
	220	ENDIF
	221	!
	222	DO jk = 1, jpkm1 ! interior values
	223	zdt = p2dt(jk)
	224	DO jj = 2, jpjm1
	225	DO ji = fs_2, fs_jpim1 ! vector opt.
	226	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3w(ji,jj,jk+1) )
	227	z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
	228	zalpha = 0.5 + z0w
	229	zw = z0w - 0.5 * pwn(ji,jj,jk+1) * zdt * zbtr
	230	zzwx = ptb(ji,jj,jk+1,jn) + zw * zslpx(ji,jj,jk+1)
	231	zzwy = ptb(ji,jj,jk ,jn) + zw * zslpx(ji,jj,jk )
	232	zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
	233	END DO
[3]	234	END DO
	235	END DO
	236
[2528]	237	! Compute & add the vertical advective trend
[503]	238	DO jk = 1, jpkm1
[2528]	239	DO jj = 2, jpjm1
[503]	240	DO ji = fs_2, fs_jpim1 ! vector opt.
[2528]	241	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	242	! vertical advective trends
	243	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
	244	! add it to the general tracer trends
	245	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[503]	246	END DO
	247	END DO
	248	END DO
[2528]	249	! ! Save the vertical advective trends for diagnostic
	250	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwx, pwn, ptb(:,:,:,jn) )
[503]	251	!
[2528]	252	ENDDO
[503]	253	!
[2715]	254	IF( wrk_not_released(3, 1,2) ) CALL ctl_stop('tra_adv_muscl: requested workspace arrays unavailable')
	255	!
[3]	256	END SUBROUTINE tra_adv_muscl
	257
	258	!!======================================================================
	259	END MODULE traadv_muscl

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

Download in other formats: