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/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps/src/OCE/TRA/traadv_ubs.F90 @ 10806

Last change on this file since 10806 was 10806, checked in by davestorkey, 5 years ago
2019/dev_r10721_KERNEL-02_Storkey_Coward_IMMERSE_first_steps branch: Latest updates. Make sure all time-dependent 3D variables are converted in previously modified modules.
Property svn:keywords set to `Id`
File size: 20.7 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
	12	!! advection trends using a third order biaised scheme
	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
[5836]	18	USE traadv_fct ! acces to routine interp_4th_cpt
[4990]	19	USE trdtra ! trends manager: tracers
	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)
	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
	40	# include "vectopt_loop_substitute.h90"
	41	!!----------------------------------------------------------------------
[9598]	42	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1152]	43	!! $Id$
[10068]	44	!! Software governed by the CeCILL license (see ./LICENSE)
[503]	45	!!----------------------------------------------------------------------
	46	CONTAINS
	47
[10806]	48	SUBROUTINE tra_adv_ubs( kt, kit000, ktlev, cdtype, p2dt, pu, pv, pw, &
[10802]	49	& pt_lev1, pt_lev2, pt_rhs, kjpt, kn_ubs_v )
[503]	50	!!----------------------------------------------------------------------
	51	!! * ROUTINE tra_adv_ubs *
	52	!!
	53	!! ** Purpose : Compute the now trend due to the advection of tracers
	54	!! and add it to the general trend of passive tracer equations.
	55	!!
[5836]	56	!! ** Method : The 3rd order Upstream Biased Scheme (UBS) is based on an
[3787]	57	!! upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005)
[519]	58	!! It is only used in the horizontal direction.
	59	!! For example the i-component of the advective fluxes are given by :
[10802]	60	!! ! e2u e3u uu ( mi(Tn) - zltu(i ) ,ktlev) if uu(i,ktlev) >= 0
[4990]	61	!! ztu = ! or
[10802]	62	!! ! e2u e3u uu ( mi(Tn) - zltu(i+1) ,ktlev) if uu(i,ktlev) < 0
[519]	63	!! where zltu is the second derivative of the before temperature field:
	64	!! zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
[5836]	65	!! This results in a dissipatively dominant (i.e. hyper-diffusive)
[519]	66	!! truncation error. The overall performance of the advection scheme
	67	!! is similar to that reported in (Farrow and Stevens, 1995).
[5836]	68	!! For stability reasons, the first term of the fluxes which corresponds
[519]	69	!! to a second order centered scheme is evaluated using the now velocity
	70	!! (centered in time) while the second term which is the diffusive part
	71	!! of the scheme, is evaluated using the before velocity (forward in time).
	72	!! Note that UBS is not positive. Do not use it on passive tracers.
[5836]	73	!! On the vertical, the advection is evaluated using a FCT scheme,
	74	!! as the UBS have been found to be too diffusive.
[6140]	75	!! kn_ubs_v argument controles whether the FCT is based on
	76	!! a 2nd order centrered scheme (kn_ubs_v=2) or on a 4th order compact
	77	!! scheme (kn_ubs_v=4).
[503]	78	!!
[10802]	79	!! ** Action : - update pt_rhs with the now advective tracer trends
[6140]	80	!! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
	81	!! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
[519]	82	!!
	83	!! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.
	84	!! Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741.
[503]	85	!!----------------------------------------------------------------------
[2528]	86	INTEGER , INTENT(in ) :: kt ! ocean time-step index
[3294]	87	INTEGER , INTENT(in ) :: kit000 ! first time step index
[10802]	88	INTEGER , INTENT(in ) :: ktlev ! time level index for source terms
[2528]	89	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	90	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[5836]	91	INTEGER , INTENT(in ) :: kn_ubs_v ! number of tracers
[6140]	92	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
[10806]	93	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pu, pv, pw ! 3 ocean transport components
[10802]	94	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: pt_lev1, pt_lev2 ! before and now tracer fields
	95	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pt_rhs ! tracer trend
[2715]	96	!
	97	INTEGER :: ji, jj, jk, jn ! dummy loop indices
[6140]	98	REAL(wp) :: ztra, zbtr, zcoef ! local scalars
[2715]	99	REAL(wp) :: zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk ! - -
	100	REAL(wp) :: zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn ! - -
[9019]	101	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztu, ztv, zltu, zltv, zti, ztw ! 3D workspace
[503]	102	!!----------------------------------------------------------------------
[3294]	103	!
	104	IF( kt == kit000 ) THEN
[503]	105	IF(lwp) WRITE(numout,*)
[2528]	106	IF(lwp) WRITE(numout,*) 'tra_adv_ubs : horizontal UBS advection scheme on ', cdtype
[503]	107	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
	108	ENDIF
[2528]	109	!
[4499]	110	l_trd = .FALSE.
[7646]	111	l_hst = .FALSE.
	112	l_ptr = .FALSE.
	113	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
	114	IF( cdtype == 'TRA' .AND. ln_diaptr ) l_ptr = .TRUE.
	115	IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
	116	& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
[4499]	117	!
[6140]	118	ztw (:,:, 1 ) = 0._wp ! surface & bottom value : set to zero for all tracers
	119	zltu(:,:,jpk) = 0._wp ; zltv(:,:,jpk) = 0._wp
[5836]	120	ztw (:,:,jpk) = 0._wp ; zti (:,:,jpk) = 0._wp
	121	!
[2528]	122	! ! ===========
	123	DO jn = 1, kjpt ! tracer loop
	124	! ! ===========
	125	!
[5836]	126	DO jk = 1, jpkm1 !== horizontal laplacian of before tracer ==!
	127	DO jj = 1, jpjm1 ! First derivative (masked gradient)
[2528]	128	DO ji = 1, fs_jpim1 ! vector opt.
[10802]	129	zeeu = e2_e1u(ji,jj) * e3u(ji,jj,jk,ktlev) * umask(ji,jj,jk)
	130	zeev = e1_e2v(ji,jj) * e3v(ji,jj,jk,ktlev) * vmask(ji,jj,jk)
	131	ztu(ji,jj,jk) = zeeu * ( pt_lev1(ji+1,jj ,jk,jn) - pt_lev1(ji,jj,jk,jn) )
	132	ztv(ji,jj,jk) = zeev * ( pt_lev1(ji ,jj+1,jk,jn) - pt_lev1(ji,jj,jk,jn) )
[2528]	133	END DO
[503]	134	END DO
[5836]	135	DO jj = 2, jpjm1 ! Second derivative (divergence)
[2528]	136	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	137	zcoef = 1._wp / ( 6._wp * e3t(ji,jj,jk,ktlev) )
[2528]	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 DO
[503]	141	END DO
[2528]	142	!
[5836]	143	END DO
[10425]	144	CALL lbc_lnk( 'traadv_ubs', zltu, 'T', 1. ) ; CALL lbc_lnk( 'traadv_ubs', zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn)
[2528]	145	!
[5836]	146	DO jk = 1, jpkm1 !== Horizontal advective fluxes ==! (UBS)
[2528]	147	DO jj = 1, jpjm1
	148	DO ji = 1, fs_jpim1 ! vector opt.
[10802]	149	zfp_ui = pu(ji,jj,jk) + ABS( pu(ji,jj,jk) ) ! upstream transport (x2)
	150	zfm_ui = pu(ji,jj,jk) - ABS( pu(ji,jj,jk) )
	151	zfp_vj = pv(ji,jj,jk) + ABS( pv(ji,jj,jk) )
	152	zfm_vj = pv(ji,jj,jk) - ABS( pv(ji,jj,jk) )
[5836]	153	! ! 2nd order centered advective fluxes (x2)
[10802]	154	zcenut = pu(ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji+1,jj ,jk,jn) )
	155	zcenvt = pv(ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji ,jj+1,jk,jn) )
[5836]	156	! ! UBS advective fluxes
[4990]	157	ztu(ji,jj,jk) = 0.5 * ( zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk) )
	158	ztv(ji,jj,jk) = 0.5 * ( zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk) )
[2528]	159	END DO
[503]	160	END DO
[5836]	161	END DO
	162	!
[10802]	163	zltu(:,:,:) = pt_rhs(:,:,:,jn) ! store the initial trends before its update
[5836]	164	!
	165	DO jk = 1, jpkm1 !== add the horizontal advective trend ==!
[503]	166	DO jj = 2, jpjm1
	167	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	168	pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) &
[4990]	169	& - ( ztu(ji,jj,jk) - ztu(ji-1,jj ,jk) &
[10802]	170	& + ztv(ji,jj,jk) - ztv(ji ,jj-1,jk) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
[503]	171	END DO
	172	END DO
[2528]	173	!
[5836]	174	END DO
	175	!
[10802]	176	zltu(:,:,:) = pt_rhs(:,:,:,jn) - zltu(:,:,:) ! Horizontal advective trend used in vertical 2nd order FCT case
[5836]	177	! ! and/or in trend diagnostic (l_trd=T)
[4990]	178	!
	179	IF( l_trd ) THEN ! trend diagnostics
[10802]	180	CALL trd_tra( kt, cdtype, jn, jptra_xad, ztu, pu, pt_lev2(:,:,:,jn) )
	181	CALL trd_tra( kt, cdtype, jn, jptra_yad, ztv, pv, pt_lev2(:,:,:,jn) )
[2528]	182	END IF
[7646]	183	!
	184	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	185	IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', ztv(:,:,:) )
	186	! ! heati/salt transport
	187	IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', ztu(:,:,:), ztv(:,:,:) )
[5836]	188	!
[7646]	189	!
[5836]	190	! !== vertical advective trend ==!
	191	!
	192	SELECT CASE( kn_ubs_v ) ! select the vertical advection scheme
	193	!
	194	CASE( 2 ) ! 2nd order FCT
	195	!
[10802]	196	IF( l_trd ) zltv(:,:,:) = pt_rhs(:,:,:,jn) ! store pt_rhs if trend diag.
[5836]	197	!
	198	! !* upstream advection with initial mass fluxes & intermediate update ==!
	199	DO jk = 2, jpkm1 ! Interior value (w-masked)
	200	DO jj = 1, jpj
	201	DO ji = 1, jpi
[10806]	202	zfp_wk = pw(ji,jj,jk) + ABS( pw(ji,jj,jk) )
	203	zfm_wk = pw(ji,jj,jk) - ABS( pw(ji,jj,jk) )
[10802]	204	ztw(ji,jj,jk) = 0.5_wp * ( zfp_wk * pt_lev1(ji,jj,jk,jn) + zfm_wk * pt_lev1(ji,jj,jk-1,jn) ) * wmask(ji,jj,jk)
[5836]	205	END DO
[2528]	206	END DO
[5836]	207	END DO
[6140]	208	IF( ln_linssh ) THEN ! top ocean value (only in linear free surface as ztw has been w-masked)
[5836]	209	IF( ln_isfcav ) THEN ! top of the ice-shelf cavities and at the ocean surface
	210	DO jj = 1, jpj
	211	DO ji = 1, jpi
[10806]	212	ztw(ji,jj, mikt(ji,jj) ) = pw(ji,jj,mikt(ji,jj)) * pt_lev1(ji,jj,mikt(ji,jj),jn) ! linear free surface
[5836]	213	END DO
	214	END DO
	215	ELSE ! no cavities: only at the ocean surface
[10806]	216	ztw(:,:,1) = pw(:,:,1) * pt_lev1(:,:,1,jn)
[5836]	217	ENDIF
	218	ENDIF
	219	!
	220	DO jk = 1, jpkm1 !* trend and after field with monotonic scheme
	221	DO jj = 2, jpjm1
	222	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	223	ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
	224	pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) + ztak
	225	zti(ji,jj,jk) = ( pt_lev1(ji,jj,jk,jn) + p2dt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
[5836]	226	END DO
	227	END DO
[2528]	228	END DO
[10425]	229	CALL lbc_lnk( 'traadv_ubs', zti, 'T', 1. ) ! Lateral boundary conditions on zti, zsi (unchanged sign)
[5836]	230	!
	231	! !* anti-diffusive flux : high order minus low order
	232	DO jk = 2, jpkm1 ! Interior value (w-masked)
	233	DO jj = 1, jpj
	234	DO ji = 1, jpi
[10806]	235	ztw(ji,jj,jk) = ( 0.5_wp * pw(ji,jj,jk) * ( pt_lev2(ji,jj,jk,jn) + pt_lev2(ji,jj,jk-1,jn) ) &
[5836]	236	& - ztw(ji,jj,jk) ) * wmask(ji,jj,jk)
	237	END DO
[503]	238	END DO
	239	END DO
[5836]	240	! ! top ocean value: high order == upstream ==>> zwz=0
[6140]	241	IF( ln_linssh ) ztw(:,:, 1 ) = 0._wp ! only ocean surface as interior zwz values have been w-masked
[5836]	242	!
[10802]	243	CALL nonosc_z( pt_lev1(:,:,:,jn), ztw, zti, p2dt, e3t(:,:,:,ktlev) ) ! monotonicity algorithm
[5836]	244	!
	245	CASE( 4 ) ! 4th order COMPACT
[10802]	246	CALL interp_4th_cpt( pt_lev2(:,:,:,jn) , ztw ) ! 4th order compact interpolation of T at w-point
[5836]	247	DO jk = 2, jpkm1
	248	DO jj = 2, jpjm1
	249	DO ji = fs_2, fs_jpim1
[10806]	250	ztw(ji,jj,jk) = pw(ji,jj,jk) * ztw(ji,jj,jk) * wmask(ji,jj,jk)
[5836]	251	END DO
[2528]	252	END DO
[503]	253	END DO
[10806]	254	IF( ln_linssh ) ztw(:,:, 1 ) = pw(:,:,1) * pt_lev2(:,:,1,jn) !!gm ISF & 4th COMPACT doesn't work
[5836]	255	!
	256	END SELECT
[2528]	257	!
[5836]	258	DO jk = 1, jpkm1 ! final trend with corrected fluxes
[2528]	259	DO jj = 2, jpjm1
	260	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	261	pt_rhs(ji,jj,jk,jn) = pt_rhs(ji,jj,jk,jn) - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
[2528]	262	END DO
[503]	263	END DO
	264	END DO
[5836]	265	!
	266	IF( l_trd ) THEN ! vertical advective trend diagnostics
[10802]	267	DO jk = 1, jpkm1 ! (compute -w.dk[ptn]= -dk[w.ptn] + pt_lev2.dk[w])
[2528]	268	DO jj = 2, jpjm1
	269	DO ji = fs_2, fs_jpim1 ! vector opt.
[10802]	270	zltv(ji,jj,jk) = pt_rhs(ji,jj,jk,jn) - zltv(ji,jj,jk) &
[10806]	271	& + pt_lev2(ji,jj,jk,jn) * ( pw(ji,jj,jk) - pw(ji,jj,jk+1) ) &
[10802]	272	& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,ktlev)
[2528]	273	END DO
	274	END DO
[503]	275	END DO
[4990]	276	CALL trd_tra( kt, cdtype, jn, jptra_zad, zltv )
[2528]	277	ENDIF
	278	!
[4990]	279	END DO
[503]	280	!
[2528]	281	END SUBROUTINE tra_adv_ubs
[503]	282
	283
[10802]	284	SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt, pe3t )
[503]	285	!!---------------------------------------------------------------------
	286	!! * ROUTINE nonosc_z *
	287	!!
	288	!! ** Purpose : compute monotonic tracer fluxes from the upstream
	289	!! scheme and the before field by a nonoscillatory algorithm
	290	!!
	291	!! ** Method : ... ???
	292	!! warning : pbef and paft must be masked, but the boundaries
	293	!! conditions on the fluxes are not necessary zalezak (1979)
	294	!! drange (1995) multi-dimensional forward-in-time and upstream-
	295	!! in-space based differencing for fluid
	296	!!----------------------------------------------------------------------
[6140]	297	REAL(wp), INTENT(in ) :: p2dt ! tracer time-step
[2528]	298	REAL(wp), DIMENSION (jpi,jpj,jpk) :: pbef ! before field
[10802]	299	REAL(wp), INTENT(in ), DIMENSION (jpi,jpj,jpk) :: pe3t ! now cell thickness field
[503]	300	REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: paft ! after field
	301	REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: pcc ! monotonic flux in the k direction
[2715]	302	!
	303	INTEGER :: ji, jj, jk ! dummy loop indices
	304	INTEGER :: ikm1 ! local integer
[6140]	305	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn ! local scalars
[9019]	306	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zbetup, zbetdo ! 3D workspace
[503]	307	!!----------------------------------------------------------------------
[3294]	308	!
[2715]	309	zbig = 1.e+40_wp
	310	zrtrn = 1.e-15_wp
	311	zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp
[5836]	312	!
[503]	313	! Search local extrema
	314	! --------------------
[5836]	315	! ! large negative value (-zbig) inside land
[503]	316	pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
	317	paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
[5836]	318	!
	319	DO jk = 1, jpkm1 ! search maximum in neighbourhood
[503]	320	ikm1 = MAX(jk-1,1)
	321	DO jj = 2, jpjm1
	322	DO ji = fs_2, fs_jpim1 ! vector opt.
	323	zbetup(ji,jj,jk) = MAX( pbef(ji ,jj ,jk ), paft(ji ,jj ,jk ), &
	324	& pbef(ji ,jj ,ikm1), pbef(ji ,jj ,jk+1), &
	325	& paft(ji ,jj ,ikm1), paft(ji ,jj ,jk+1) )
	326	END DO
	327	END DO
	328	END DO
[5836]	329	! ! large positive value (+zbig) inside land
[503]	330	pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
	331	paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
[5836]	332	!
	333	DO jk = 1, jpkm1 ! search minimum in neighbourhood
[503]	334	ikm1 = MAX(jk-1,1)
	335	DO jj = 2, jpjm1
	336	DO ji = fs_2, fs_jpim1 ! vector opt.
	337	zbetdo(ji,jj,jk) = MIN( pbef(ji ,jj ,jk ), paft(ji ,jj ,jk ), &
	338	& pbef(ji ,jj ,ikm1), pbef(ji ,jj ,jk+1), &
	339	& paft(ji ,jj ,ikm1), paft(ji ,jj ,jk+1) )
	340	END DO
	341	END DO
	342	END DO
[5836]	343	! ! restore masked values to zero
[503]	344	pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
	345	paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
[5836]	346	!
	347	! Positive and negative part of fluxes and beta terms
	348	! ---------------------------------------------------
[503]	349	DO jk = 1, jpkm1
	350	DO jj = 2, jpjm1
	351	DO ji = fs_2, fs_jpim1 ! vector opt.
	352	! positive & negative part of the flux
	353	zpos = MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
	354	zneg = MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
	355	! up & down beta terms
[10802]	356	zbt = e1e2t(ji,jj) * pe3t(ji,jj,jk) / p2dt
[503]	357	zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
	358	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
	359	END DO
	360	END DO
	361	END DO
[5836]	362	!
[503]	363	! monotonic flux in the k direction, i.e. pcc
	364	! -------------------------------------------
	365	DO jk = 2, jpkm1
	366	DO jj = 2, jpjm1
	367	DO ji = fs_2, fs_jpim1 ! vector opt.
	368	za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
	369	zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
	370	zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
	371	pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
	372	END DO
	373	END DO
	374	END DO
	375	!
	376	END SUBROUTINE nonosc_z
	377
	378	!!======================================================================
	379	END MODULE traadv_ubs

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

Download in other formats: