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_r13383_HPC-02_Daley_Tiling/src/OCE/TRA – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/OCE/TRA/traadv_ubs.F90 @ 13516

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

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

Download in other formats: