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 NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/TRA/traadv_ubs.F90 @ 14644

Last change on this file since 14644 was 14644, checked in by sparonuz, 3 years ago
Merge trunk -r14642:HEAD
Property svn:keywords set to `Id`
File size: 19.6 KB

Rev	Line
[503]	1	MODULE traadv_ubs
	2	!!==============================================================================
	3	!! * MODULE traadv_ubs *
	4	!! Ocean active tracers: horizontal & vertical advective trend
	5	!!==============================================================================
[2528]	6	!! History : 1.0 ! 2006-08 (L. Debreu, R. Benshila) Original code
	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
[14072]	12	!! advection trends using a third order biaised scheme
[503]	13	!!----------------------------------------------------------------------
[3625]	14	USE oce ! ocean dynamics and active tracers
	15	USE dom_oce ! ocean space and time domain
[4990]	16	USE trc_oce ! share passive tracers/Ocean variables
	17	USE trd_oce ! trends: ocean variables
[14072]	18	USE traadv_fct ! acces to routine interp_4th_cpt
	19	USE trdtra ! trends manager: tracers
[4990]	20	USE diaptr ! poleward transport diagnostics
[7646]	21	USE diaar5 ! AR5 diagnostics
[4990]	22	!
[9019]	23	USE iom ! I/O library
[9124]	24	USE in_out_manager ! I/O manager
[9019]	25	USE lib_mpp ! massively parallel library
[3625]	26	USE lbclnk ! ocean lateral boundary condition (or mpp link)
[14072]	27	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[503]	28
	29	IMPLICIT NONE
	30	PRIVATE
	31
	32	PUBLIC tra_adv_ubs ! routine called by traadv module
	33
[7646]	34	LOGICAL :: l_trd ! flag to compute trends
	35	LOGICAL :: l_ptr ! flag to compute poleward transport
	36	LOGICAL :: l_hst ! flag to compute heat transport
[503]	37
[7646]	38
[503]	39	!! * Substitutions
[12377]	40	# include "do_loop_substitute.h90"
[13237]	41	# include "domzgr_substitute.h90"
[14219]	42	# include "single_precision_substitute.h90"
[503]	43	!!----------------------------------------------------------------------
[9598]	44	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]	45	!! $Id$
[10068]	46	!! Software governed by the CeCILL license (see ./LICENSE)
[503]	47	!!----------------------------------------------------------------------
	48	CONTAINS
	49
[12377]	50	SUBROUTINE tra_adv_ubs( kt, kit000, cdtype, p2dt, pU, pV, pW, &
	51	& Kbb, Kmm, pt, kjpt, Krhs, kn_ubs_v )
[503]	52	!!----------------------------------------------------------------------
	53	!! * ROUTINE tra_adv_ubs *
[14072]	54	!!
[503]	55	!! ** Purpose : Compute the now trend due to the advection of tracers
	56	!! and add it to the general trend of passive tracer equations.
	57	!!
[5836]	58	!! ** Method : The 3rd order Upstream Biased Scheme (UBS) is based on an
[3787]	59	!! upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005)
[519]	60	!! It is only used in the horizontal direction.
	61	!! For example the i-component of the advective fluxes are given by :
[3787]	62	!! ! e2u e3u un ( mi(Tn) - zltu(i ) ) if un(i) >= 0
[14219]	63	!! ztu = ! or
[3787]	64	!! ! e2u e3u un ( mi(Tn) - zltu(i+1) ) if un(i) < 0
[519]	65	!! where zltu is the second derivative of the before temperature field:
	66	!! zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
[14072]	67	!! This results in a dissipatively dominant (i.e. hyper-diffusive)
	68	!! truncation error. The overall performance of the advection scheme
	69	!! is similar to that reported in (Farrow and Stevens, 1995).
[5836]	70	!! For stability reasons, the first term of the fluxes which corresponds
[14072]	71	!! to a second order centered scheme is evaluated using the now velocity
	72	!! (centered in time) while the second term which is the diffusive part
	73	!! of the scheme, is evaluated using the before velocity (forward in time).
[519]	74	!! Note that UBS is not positive. Do not use it on passive tracers.
[5836]	75	!! On the vertical, the advection is evaluated using a FCT scheme,
[14072]	76	!! as the UBS have been found to be too diffusive.
	77	!! kn_ubs_v argument controles whether the FCT is based on
	78	!! a 2nd order centrered scheme (kn_ubs_v=2) or on a 4th order compact
[6140]	79	!! scheme (kn_ubs_v=4).
[503]	80	!!
[12377]	81	!! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends
[6140]	82	!! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
[12377]	83	!! - poleward advective heat and salt transport (ln_diaptr=T)
[519]	84	!!
[14072]	85	!! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.
[13237]	86	!! Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731�1741.
[503]	87	!!----------------------------------------------------------------------
[12377]	88	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	89	INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices
	90	INTEGER , INTENT(in ) :: kit000 ! first time step index
	91	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	92	INTEGER , INTENT(in ) :: kjpt ! number of tracers
	93	INTEGER , INTENT(in ) :: kn_ubs_v ! number of tracers
	94	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
[13982]	95	! TEMP: [tiling] This can be A2D(nn_hls) if using XIOS (subdomain support)
[12377]	96	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume transport components
[14219]	97	REAL(dp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
[2715]	98	!
	99	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[6140]	100	REAL(wp) :: ztra, zbtr, zcoef ! local scalars
[2715]	101	REAL(wp) :: zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk ! - -
	102	REAL(wp) :: zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn ! - -
[13982]	103	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: ztu, ztv, zltu, zltv, zti, ztw ! 3D workspace
[503]	104	!!----------------------------------------------------------------------
[3294]	105	!
[13982]	106	IF( ntile == 0 .OR. ntile == 1 ) THEN ! Do only on the first tile
	107	IF( kt == kit000 ) THEN
	108	IF(lwp) WRITE(numout,*)
	109	IF(lwp) WRITE(numout,*) 'tra_adv_ubs : horizontal UBS advection scheme on ', cdtype
	110	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
	111	ENDIF
	112	!
	113	l_trd = .FALSE.
	114	l_hst = .FALSE.
	115	l_ptr = .FALSE.
	116	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
	117	IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
	118	IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
	119	& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
[503]	120	ENDIF
[2528]	121	!
[6140]	122	ztw (:,:, 1 ) = 0._wp ! surface & bottom value : set to zero for all tracers
	123	zltu(:,:,jpk) = 0._wp ; zltv(:,:,jpk) = 0._wp
[5836]	124	ztw (:,:,jpk) = 0._wp ; zti (:,:,jpk) = 0._wp
[2528]	125	! ! ===========
	126	DO jn = 1, kjpt ! tracer loop
	127	! ! ===========
[14072]	128	!
[13497]	129	DO jk = 1, jpkm1 !== horizontal laplacian of before tracer ==!
[13982]	130	DO_2D( nn_hls, nn_hls-1, nn_hls, nn_hls-1 ) ! First derivative (masked gradient)
[12377]	131	zeeu = e2_e1u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
	132	zeev = e1_e2v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
	133	ztu(ji,jj,jk) = zeeu * ( pt(ji+1,jj ,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
	134	ztv(ji,jj,jk) = zeev * ( pt(ji ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) )
	135	END_2D
[13982]	136	DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! Second derivative (divergence)
[12377]	137	zcoef = 1._wp / ( 6._wp * e3t(ji,jj,jk,Kmm) )
	138	zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef
	139	zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef
	140	END_2D
[14072]	141	!
	142	END DO
[14644]	143	IF (nn_hls.EQ.1) CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1.0_wp, zltv, 'T', 1.0_wp ) ! Lateral boundary cond. (unchanged sgn)
[14072]	144	!
[13497]	145	DO_3D( 1, 0, 1, 0, 1, jpkm1 ) !== Horizontal advective fluxes ==! (UBS)
	146	zfp_ui = pU(ji,jj,jk) + ABS( pU(ji,jj,jk) ) ! upstream transport (x2)
[12377]	147	zfm_ui = pU(ji,jj,jk) - ABS( pU(ji,jj,jk) )
	148	zfp_vj = pV(ji,jj,jk) + ABS( pV(ji,jj,jk) )
	149	zfm_vj = pV(ji,jj,jk) - ABS( pV(ji,jj,jk) )
	150	! ! 2nd order centered advective fluxes (x2)
	151	zcenut = pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) )
	152	zcenvt = pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) )
	153	! ! UBS advective fluxes
	154	ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) )
	155	ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) )
	156	END_3D
[5836]	157	!
[13982]	158	DO_3D( 1, 1, 1, 1, 1, jpk )
	159	zltu(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) ! store the initial trends before its update
	160	END_3D
[5836]	161	!
	162	DO jk = 1, jpkm1 !== add the horizontal advective trend ==!
[13295]	163	DO_2D( 0, 0, 0, 0 )
[12377]	164	pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) &
	165	& - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk) &
[13237]	166	& + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) ) &
	167	& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
[12377]	168	END_2D
[14072]	169	!
[5836]	170	END DO
	171	!
[13982]	172	DO_3D( 1, 1, 1, 1, 1, jpk )
	173	zltu(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltu(ji,jj,jk) ! Horizontal advective trend used in vertical 2nd order FCT case
	174	END_3D ! and/or in trend diagnostic (l_trd=T)
	175	!
[4990]	176	IF( l_trd ) THEN ! trend diagnostics
[14219]	177	CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, ztu, pU, CASTWP(pt(:,:,:,jn,Kmm)) )
	178	CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, ztv, pV, CASTWP(pt(:,:,:,jn,Kmm)) )
[2528]	179	END IF
[14072]	180	!
[7646]	181	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	182	IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) )
	183	! ! heati/salt transport
	184	IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) )
[5836]	185	!
[7646]	186	!
[5836]	187	! !== vertical advective trend ==!
	188	!
	189	SELECT CASE( kn_ubs_v ) ! select the vertical advection scheme
	190	!
[14072]	191	CASE( 2 ) ! 2nd order FCT
	192	!
[13982]	193	IF( l_trd ) THEN
	194	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	195	zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) ! store pt(:,:,:,:,Krhs) if trend diag.
	196	END_3D
	197	ENDIF
[5836]	198	!
[13497]	199	! !* upstream advection with initial mass fluxes & intermediate update ==!
[13295]	200	DO_3D( 1, 1, 1, 1, 2, jpkm1 )
[12377]	201	zfp_wk = pW(ji,jj,jk) + ABS( pW(ji,jj,jk) )
	202	zfm_wk = pW(ji,jj,jk) - ABS( pW(ji,jj,jk) )
	203	ztw(ji,jj,jk) = 0.5_wp * ( zfp_wk * pt(ji,jj,jk,jn,Kbb) + zfm_wk * pt(ji,jj,jk-1,jn,Kbb) ) * wmask(ji,jj,jk)
	204	END_3D
[13497]	205	IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as ztw has been w-masked)
	206	IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface
[13295]	207	DO_2D( 1, 1, 1, 1 )
[14072]	208	ztw(ji,jj, mikt(ji,jj) ) = pW(ji,jj,mikt(ji,jj)) * pt(ji,jj,mikt(ji,jj),jn,Kbb) ! linear free surface
[12377]	209	END_2D
[13497]	210	ELSE ! no cavities: only at the ocean surface
[13982]	211	DO_2D( 1, 1, 1, 1 )
	212	ztw(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb)
	213	END_2D
[5836]	214	ENDIF
	215	ENDIF
	216	!
[13497]	217	DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !* trend and after field with monotonic scheme
[13237]	218	ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) &
	219	& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
[14072]	220	pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + ztak
[12377]	221	zti(ji,jj,jk) = ( pt(ji,jj,jk,jn,Kbb) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
	222	END_3D
[5836]	223	!
	224	! !* anti-diffusive flux : high order minus low order
[13295]	225	DO_3D( 1, 1, 1, 1, 2, jpkm1 )
[12377]	226	ztw(ji,jj,jk) = ( 0.5_wp * pW(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji,jj,jk-1,jn,Kmm) ) &
	227	& - ztw(ji,jj,jk) ) * wmask(ji,jj,jk)
	228	END_3D
[5836]	229	! ! top ocean value: high order == upstream ==>> zwz=0
[6140]	230	IF( ln_linssh ) ztw(:,:, 1 ) = 0._wp ! only ocean surface as interior zwz values have been w-masked
[5836]	231	!
[12377]	232	CALL nonosc_z( Kmm, pt(:,:,:,jn,Kbb), ztw, zti, p2dt ) ! monotonicity algorithm
[5836]	233	!
	234	CASE( 4 ) ! 4th order COMPACT
[14219]	235	CALL interp_4th_cpt( CASTWP(pt(:,:,:,jn,Kmm)) , ztw )
[13295]	236	DO_3D( 0, 0, 0, 0, 2, jpkm1 )
[12377]	237	ztw(ji,jj,jk) = pW(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
	238	END_3D
[13982]	239	IF( ln_linssh ) THEN
	240	DO_2D( 1, 1, 1, 1 )
	241	ztw(ji,jj,1) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm) !!gm ISF & 4th COMPACT doesn't work
	242	END_2D
	243	ENDIF
[5836]	244	!
	245	END SELECT
[2528]	246	!
[13497]	247	DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! final trend with corrected fluxes
[13237]	248	pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) &
	249	& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
[12377]	250	END_3D
[5836]	251	!
[13497]	252	IF( l_trd ) THEN ! vertical advective trend diagnostics
	253	DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
[12377]	254	zltv(ji,jj,jk) = pt(ji,jj,jk,jn,Krhs) - zltv(ji,jj,jk) &
	255	& + pt(ji,jj,jk,jn,Kmm) * ( pW(ji,jj,jk) - pW(ji,jj,jk+1) ) &
	256	& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
	257	END_3D
	258	CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zltv )
[2528]	259	ENDIF
	260	!
[4990]	261	END DO
[503]	262	!
[2528]	263	END SUBROUTINE tra_adv_ubs
[503]	264
	265
[12377]	266	SUBROUTINE nonosc_z( Kmm, pbef, pcc, paft, p2dt )
[503]	267	!!---------------------------------------------------------------------
	268	!! * ROUTINE nonosc_z *
	269	!!
[14072]	270	!! ** Purpose : compute monotonic tracer fluxes from the upstream
	271	!! scheme and the before field by a nonoscillatory algorithm
	272	!!
[503]	273	!! ** Method : ... ???
	274	!! warning : pbef and paft must be masked, but the boundaries
	275	!! conditions on the fluxes are not necessary zalezak (1979)
	276	!! drange (1995) multi-dimensional forward-in-time and upstream-
	277	!! in-space based differencing for fluid
	278	!!----------------------------------------------------------------------
[13982]	279	INTEGER , INTENT(in ) :: Kmm ! time level index
	280	REAL(wp), INTENT(in ) :: p2dt ! tracer time-step
[14219]	281	REAL(dp), DIMENSION(jpi,jpj,jpk) :: pbef ! before field
[13982]	282	REAL(wp), INTENT(inout), DIMENSION(A2D(nn_hls) ,jpk) :: paft ! after field
	283	REAL(wp), INTENT(inout), DIMENSION(A2D(nn_hls) ,jpk) :: pcc ! monotonic flux in the k direction
[2715]	284	!
	285	INTEGER :: ji, jj, jk ! dummy loop indices
	286	INTEGER :: ikm1 ! local integer
[6140]	287	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars
[13982]	288	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zbetup, zbetdo ! 3D workspace
[503]	289	!!----------------------------------------------------------------------
[3294]	290	!
[13226]	291	zbig = 1.e+38_wp
[2715]	292	zrtrn = 1.e-15_wp
	293	zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp
[5836]	294	!
[503]	295	! Search local extrema
	296	! --------------------
[5836]	297	! ! large negative value (-zbig) inside land
[13982]	298	DO_3D( 0, 0, 0, 0, 1, jpk )
	299	pbef(ji,jj,jk) = pbef(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1.e0 - tmask(ji,jj,jk) )
	300	paft(ji,jj,jk) = paft(ji,jj,jk) * tmask(ji,jj,jk) - zbig * ( 1.e0 - tmask(ji,jj,jk) )
	301	END_3D
[5836]	302	!
	303	DO jk = 1, jpkm1 ! search maximum in neighbourhood
[503]	304	ikm1 = MAX(jk-1,1)
[13295]	305	DO_2D( 0, 0, 0, 0 )
[12377]	306	zbetup(ji,jj,jk) = MAX( pbef(ji ,jj ,jk ), paft(ji ,jj ,jk ), &
	307	& pbef(ji ,jj ,ikm1), pbef(ji ,jj ,jk+1), &
	308	& paft(ji ,jj ,ikm1), paft(ji ,jj ,jk+1) )
	309	END_2D
[503]	310	END DO
[5836]	311	! ! large positive value (+zbig) inside land
[13982]	312	DO_3D( 0, 0, 0, 0, 1, jpk )
	313	pbef(ji,jj,jk) = pbef(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1.e0 - tmask(ji,jj,jk) )
	314	paft(ji,jj,jk) = paft(ji,jj,jk) * tmask(ji,jj,jk) + zbig * ( 1.e0 - tmask(ji,jj,jk) )
	315	END_3D
[5836]	316	!
	317	DO jk = 1, jpkm1 ! search minimum in neighbourhood
[503]	318	ikm1 = MAX(jk-1,1)
[13295]	319	DO_2D( 0, 0, 0, 0 )
[12377]	320	zbetdo(ji,jj,jk) = MIN( pbef(ji ,jj ,jk ), paft(ji ,jj ,jk ), &
	321	& pbef(ji ,jj ,ikm1), pbef(ji ,jj ,jk+1), &
	322	& paft(ji ,jj ,ikm1), paft(ji ,jj ,jk+1) )
	323	END_2D
[503]	324	END DO
[5836]	325	! ! restore masked values to zero
[13982]	326	DO_3D( 0, 0, 0, 0, 1, jpk )
	327	pbef(ji,jj,jk) = pbef(ji,jj,jk) * tmask(ji,jj,jk)
	328	paft(ji,jj,jk) = paft(ji,jj,jk) * tmask(ji,jj,jk)
	329	END_3D
[5836]	330	!
	331	! Positive and negative part of fluxes and beta terms
	332	! ---------------------------------------------------
[13295]	333	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	334	! positive & negative part of the flux
	335	zpos = MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
	336	zneg = MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
	337	! up & down beta terms
	338	zbt = e1e2t(ji,jj) * e3t(ji,jj,jk,Kmm) / p2dt
	339	zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
	340	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
	341	END_3D
[5836]	342	!
[503]	343	! monotonic flux in the k direction, i.e. pcc
	344	! -------------------------------------------
[13295]	345	DO_3D( 0, 0, 0, 0, 2, jpkm1 )
[12377]	346	za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
	347	zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
[13226]	348	zc = 0.5 * ( 1.e0 + SIGN( 1.0_wp, pcc(ji,jj,jk) ) )
[12377]	349	pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
	350	END_3D
[503]	351	!
	352	END SUBROUTINE nonosc_z
	353
	354	!!======================================================================
	355	END MODULE traadv_ubs

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

Download in other formats: