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 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_ubs.F90 @ 4409

Last change on this file since 4409 was 4409, checked in by trackstand2, 10 years ago
Changes to allow jpk to be modified to deepest level within a subdomain. jpkorig holds original value.
Property svn:keywords set to `Id`
File size: 22.3 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	!!----------------------------------------------------------------------
	14	USE oce ! ocean dynamics and active tracers
	15	USE dom_oce ! ocean space and time domain
[2528]	16	USE trdmod_oce ! ocean space and time domain
	17	USE trdtra
[503]	18	USE lib_mpp
	19	USE lbclnk ! ocean lateral boundary condition (or mpp link)
	20	USE in_out_manager ! I/O manager
	21	USE diaptr ! poleward transport diagnostics
	22	USE dynspg_oce ! choice/control of key cpp for surface pressure gradient
[2528]	23	USE trc_oce ! share passive tracers/Ocean variables
[503]	24
	25	IMPLICIT NONE
	26	PRIVATE
	27
	28	PUBLIC tra_adv_ubs ! routine called by traadv module
	29
[2528]	30	LOGICAL :: l_trd ! flag to compute trends or not
[503]	31
[3211]	32	!! * Control permutation of array indices
	33	# include "oce_ftrans.h90"
	34	# include "dom_oce_ftrans.h90"
	35	# include "trc_oce_ftrans.h90"
	36
[503]	37	!! * Substitutions
	38	# include "domzgr_substitute.h90"
	39	# include "vectopt_loop_substitute.h90"
	40	!!----------------------------------------------------------------------
[2528]	41	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]	42	!! $Id$
[2528]	43	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[503]	44	!!----------------------------------------------------------------------
	45	CONTAINS
	46
[2528]	47	SUBROUTINE tra_adv_ubs ( kt, cdtype, p2dt, pun, pvn, pwn, &
	48	& ptb, ptn, pta, kjpt )
[503]	49	!!----------------------------------------------------------------------
	50	!! * ROUTINE tra_adv_ubs *
	51	!!
	52	!! ** Purpose : Compute the now trend due to the advection of tracers
	53	!! and add it to the general trend of passive tracer equations.
	54	!!
[519]	55	!! ** Method : The upstream biased third (UBS) is order scheme based
	56	!! on an upstream-biased parabolic interpolation (Shchepetkin and McWilliams 2005)
	57	!! It is only used in the horizontal direction.
	58	!! For example the i-component of the advective fluxes are given by :
	59	!! ! e1u e3u un ( mi(Tn) - zltu(i ) ) if un(i) >= 0
	60	!! zwx = ! or
	61	!! ! e1u e3u un ( mi(Tn) - zltu(i+1) ) if un(i) < 0
	62	!! where zltu is the second derivative of the before temperature field:
	63	!! zltu = 1/e3t di[ e2u e3u / e1u di[Tb] ]
	64	!! This results in a dissipatively dominant (i.e. hyper-diffusive)
	65	!! truncation error. The overall performance of the advection scheme
	66	!! is similar to that reported in (Farrow and Stevens, 1995).
	67	!! For stability reasons, the first term of the fluxes which corresponds
	68	!! to a second order centered scheme is evaluated using the now velocity
	69	!! (centered in time) while the second term which is the diffusive part
	70	!! of the scheme, is evaluated using the before velocity (forward in time).
	71	!! Note that UBS is not positive. Do not use it on passive tracers.
	72	!! On the vertical, the advection is evaluated using a TVD scheme, as
	73	!! the UBS have been found to be too diffusive.
[503]	74	!!
[2528]	75	!! ** Action : - update (pta) with the now advective tracer trends
[519]	76	!!
	77	!! Reference : Shchepetkin, A. F., J. C. McWilliams, 2005, Ocean Modelling, 9, 347-404.
	78	!! Farrow, D.E., Stevens, D.P., 1995, J. Phys. Ocean. 25, 1731Ð1741.
[503]	79	!!----------------------------------------------------------------------
[2715]	80	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	81	USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as workspace
	82	USE wrk_nemo, ONLY: ztu => wrk_3d_1 , ztv => wrk_3d_2 ! 3D workspace
	83	USE wrk_nemo, ONLY: zltu => wrk_3d_3 , zltv => wrk_3d_4 ! - -
	84	USE wrk_nemo, ONLY: zti => wrk_3d_5 , ztw => wrk_3d_6 ! - -
[3211]	85
	86	!! DCSE_NEMO: need additional directives for renamed module variables
	87	!FTRANS zwx zwy ztu ztv zltu zltv zti ztw :I :I :z
	88
[2715]	89	!
[2528]	90	INTEGER , INTENT(in ) :: kt ! ocean time-step index
	91	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
	92	INTEGER , INTENT(in ) :: kjpt ! number of tracers
[4409]	93	REAL(wp), DIMENSION( jpkorig ), INTENT(in ) :: p2dt ! vertical profile of tracer time-step
[3211]	94
	95	! REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
	96	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
	97	! REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
	98
	99	!FTRANS pun pvn pwn :I :I :z
	100	!FTRANS ptb ptn pta :I :I :z :
[4409]	101	REAL(wp), INTENT(in ) :: pun(jpi,jpj,jpkorig) ! ocean velocity component (u)
	102	REAL(wp), INTENT(in ) :: pvn(jpi,jpj,jpkorig) ! ocean velocity component (v)
	103	REAL(wp), INTENT(in ) :: pwn(jpi,jpj,jpkorig) ! ocean velocity component (w)
[3211]	104	!! DCSE_NEMO: Next two arguments made inout to silence the cray compile,
	105	!! which rightly complains about the call to nonosc_v (which also has them
	106	!! as inout)
[4409]	107	REAL(wp), INTENT(inout) :: ptb(jpi,jpj,jpkorig,kjpt) ! tracer fields (before)
	108	REAL(wp), INTENT(inout) :: ptn(jpi,jpj,jpkorig,kjpt) ! tracer fields (now)
	109	REAL(wp), INTENT(inout) :: pta(jpi,jpj,jpkorig,kjpt) ! tracer trend
[3211]	110
[2715]	111	!
	112	INTEGER :: ji, jj, jk, jn ! dummy loop indices
	113	REAL(wp) :: ztra, zbtr, zcoef, z2dtt ! local scalars
	114	REAL(wp) :: zfp_ui, zfm_ui, zcenut, ztak, zfp_wk, zfm_wk ! - -
	115	REAL(wp) :: zfp_vj, zfm_vj, zcenvt, zeeu, zeev, z_hdivn ! - -
[503]	116	!!----------------------------------------------------------------------
	117
[2715]	118	IF( wrk_in_use(3, 1,2,3,4,5,6) )THEN
	119	CALL ctl_stop('tra_adv_ubs: requested workspace arrays unavailable') ; RETURN
	120	ENDIF
	121
[2528]	122	IF( kt == nit000 ) THEN
[503]	123	IF(lwp) WRITE(numout,*)
[2528]	124	IF(lwp) WRITE(numout,*) 'tra_adv_ubs : horizontal UBS advection scheme on ', cdtype
[503]	125	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
	126	!
[2528]	127	l_trd = .FALSE.
	128	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
[503]	129	ENDIF
[2528]	130	!
	131	! ! ===========
	132	DO jn = 1, kjpt ! tracer loop
	133	! ! ===========
	134	! 1. Bottom value : flux set to zero
	135	! ----------------------------------
	136	zltu(:,:,jpk) = 0.e0 ; zltv(:,:,jpk) = 0.e0
	137	!
[3211]	138	#if defined key_z_first
	139	DO jj = 1, jpjm1
	140	DO ji = 1, jpim1
	141	DO jk = 1, jpkm1
	142	#else
[2528]	143	DO jk = 1, jpkm1 ! Horizontal slab
	144	!
	145	! Laplacian
	146	DO jj = 1, jpjm1 ! First derivative (gradient)
	147	DO ji = 1, fs_jpim1 ! vector opt.
[3211]	148	#endif
[2528]	149	zeeu = e2u(ji,jj) * fse3u(ji,jj,jk) / e1u(ji,jj) * umask(ji,jj,jk)
	150	zeev = e1v(ji,jj) * fse3v(ji,jj,jk) / e2v(ji,jj) * vmask(ji,jj,jk)
	151	ztu(ji,jj,jk) = zeeu * ( ptb(ji+1,jj ,jk,jn) - ptb(ji,jj,jk,jn) )
	152	ztv(ji,jj,jk) = zeev * ( ptb(ji ,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
	153	END DO
[503]	154	END DO
[3211]	155	#if defined key_z_first
	156	END DO
	157	DO jj = 2, jpjm1 ! Second derivative (divergence)
	158	DO ji = 2, jpim1
	159	DO jk = 1, jpkm1
	160	#else
[2528]	161	DO jj = 2, jpjm1 ! Second derivative (divergence)
	162	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	163	#endif
[2528]	164	zcoef = 1. / ( 6. * fse3t(ji,jj,jk) )
	165	zltu(ji,jj,jk) = ( ztu(ji,jj,jk) - ztu(ji-1,jj,jk) ) * zcoef
	166	zltv(ji,jj,jk) = ( ztv(ji,jj,jk) - ztv(ji,jj-1,jk) ) * zcoef
	167	END DO
[503]	168	END DO
[2528]	169	!
	170	END DO ! End of slab
	171	CALL lbc_lnk( zltu, 'T', 1. ) ; CALL lbc_lnk( zltv, 'T', 1. ) ! Lateral boundary cond. (unchanged sgn)
[503]	172
[2528]	173	!
	174	! Horizontal advective fluxes
[3211]	175	#if defined key_z_first
	176	DO jj = 1, jpjm1
	177	DO ji = 1, jpim1
	178	DO jk = 1, jpkm1
	179	#else
[2528]	180	DO jk = 1, jpkm1 ! Horizontal slab
	181	DO jj = 1, jpjm1
	182	DO ji = 1, fs_jpim1 ! vector opt.
[3211]	183	#endif
[2528]	184	! upstream transport
	185	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
	186	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
	187	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
	188	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
	189	! centered scheme
	190	zcenut = 0.5 * pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) )
	191	zcenvt = 0.5 * pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) )
	192	! UBS scheme
	193	zwx(ji,jj,jk) = zcenut - zfp_ui * zltu(ji,jj,jk) - zfm_ui * zltu(ji+1,jj,jk)
	194	zwy(ji,jj,jk) = zcenvt - zfp_vj * zltv(ji,jj,jk) - zfm_vj * zltv(ji,jj+1,jk)
	195	END DO
[503]	196	END DO
[2528]	197	END DO ! End of slab
[503]	198
[2528]	199	zltu(:,:,:) = pta(:,:,:,jn) ! store pta trends
[503]	200
[2528]	201	! Horizontal advective trends
[3211]	202	#if defined key_z_first
	203	DO jj = 2, jpjm1
	204	DO ji = 2, jpim1
	205	DO jk = 1, jpkm1
	206	#else
[503]	207	DO jk = 1, jpkm1
[2528]	208	! Tracer flux divergence at t-point added to the general trend
[503]	209	DO jj = 2, jpjm1
	210	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	211	#endif
[2528]	212	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	213	! horizontal advective
	214	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) &
	215	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) )
	216	! add it to the general tracer trends
	217	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
[503]	218	END DO
	219	END DO
[2528]	220	!
	221	END DO ! End of slab
[503]	222
[2528]	223	! Horizontal trend used in tra_adv_ztvd subroutine
	224	zltu(:,:,:) = pta(:,:,:,jn) - zltu(:,:,:)
[503]	225
[2528]	226	! 3. Save the horizontal advective trends for diagnostic
	227	! ------------------------------------------------------
	228	! ! trend diagnostics (contribution of upstream fluxes)
	229	IF( l_trd ) THEN
	230	CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
	231	CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
	232	END IF
	233	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
	234	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
	235	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) )
	236	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) )
[503]	237	ENDIF
[2528]	238
	239	! TVD scheme for the vertical direction
	240	! ----------------------
	241	IF( l_trd ) zltv(:,:,:) = pta(:,:,:,jn) ! store pta if trend diag.
[503]	242
[2528]	243	! Bottom value : flux set to zero
	244	ztw(:,:,jpk) = 0.e0 ; zti(:,:,jpk) = 0.e0
[503]	245
[2528]	246	! Surface value
	247	IF( lk_vvl ) THEN ; ztw(:,:,1) = 0.e0 ! variable volume : flux set to zero
	248	ELSE ; ztw(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn) ! free constant surface
	249	ENDIF
	250	! upstream advection with initial mass fluxes & intermediate update
	251	! -------------------------------------------------------------------
	252	! Interior value
[3211]	253	#if defined key_z_first
	254	DO jj = 1, jpj
	255	DO ji = 1, jpi
	256	DO jk = 2, jpkm1
	257	#else
[2528]	258	DO jk = 2, jpkm1
	259	DO jj = 1, jpj
	260	DO ji = 1, jpi
[3211]	261	#endif
[2528]	262	zfp_wk = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
	263	zfm_wk = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
	264	ztw(ji,jj,jk) = 0.5 * ( zfp_wk * ptb(ji,jj,jk,jn) + zfm_wk * ptb(ji,jj,jk-1,jn) )
	265	END DO
	266	END DO
	267	END DO
	268	! update and guess with monotonic sheme
[3211]	269	#if defined key_z_first
	270	DO jj = 2, jpjm1
	271	DO ji = 2, jpim1
	272	DO jk = 1, jpkm1
	273	z2dtt = p2dt(jk)
	274	#else
[503]	275	DO jk = 1, jpkm1
[2528]	276	z2dtt = p2dt(jk)
[503]	277	DO jj = 2, jpjm1
	278	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	279	#endif
[2528]	280	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	281	ztak = - ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) ) * zbtr
	282	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztak
	283	zti(ji,jj,jk) = ( ptb(ji,jj,jk,jn) + z2dtt * ( ztak + zltu(ji,jj,jk) ) ) * tmask(ji,jj,jk)
[503]	284	END DO
	285	END DO
	286	END DO
	287	!
[2528]	288	CALL lbc_lnk( zti, 'T', 1. ) ! Lateral boundary conditions on zti, zsi (unchanged sign)
[503]	289
[2528]	290	! antidiffusive flux : high order minus low order
[3211]	291	#if defined key_z_first
	292	DO jj = 1, jpj
	293	DO ji = 1, jpi
	294	ztw(ji,jj,1) = 0.e0 ! Surface value
	295	DO jk = 2, jpkm1 ! Interior value
	296	#else
[2528]	297	ztw(:,:,1) = 0.e0 ! Surface value
	298	DO jk = 2, jpkm1 ! Interior value
	299	DO jj = 1, jpj
	300	DO ji = 1, jpi
[3211]	301	#endif
[2528]	302	ztw(ji,jj,jk) = 0.5 * pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) ) - ztw(ji,jj,jk)
	303	END DO
[503]	304	END DO
	305	END DO
[2528]	306	!
	307	CALL nonosc_z( ptb(:,:,:,jn), ztw, zti, p2dt ) ! monotonicity algorithm
[503]	308
[2528]	309	! final trend with corrected fluxes
[3211]	310	#if defined key_z_first
	311	DO jj = 2, jpjm1
	312	DO ji = 2, jpim1
	313	DO jk = 1, jpkm1
	314	#else
[2528]	315	DO jk = 1, jpkm1
	316	DO jj = 2, jpjm1
	317	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	318	#endif
[2528]	319	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	320	! k- vertical advective trends
	321	ztra = - zbtr * ( ztw(ji,jj,jk) - ztw(ji,jj,jk+1) )
	322	! added to the general tracer trends
	323	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
	324	END DO
[503]	325	END DO
	326	END DO
	327
[2528]	328	! Save the final vertical advective trends
	329	IF( l_trd ) THEN ! vertical advective trend diagnostics
[3211]	330	! (compute -w.dk[ptn]= -dk[w.ptn] + ptn.dk[w])
	331	#if defined key_z_first
	332	DO jj = 2, jpjm1
	333	DO ji = 2, jpim1
	334	DO jk = 1, jpkm1
	335	#else
	336	DO jk = 1, jpkm1
[2528]	337	DO jj = 2, jpjm1
	338	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	339	#endif
[2528]	340	zbtr = 1.e0 / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
	341	z_hdivn = ( pwn(ji,jj,jk) - pwn(ji,jj,jk+1) ) * zbtr
	342	zltv(ji,jj,jk) = pta(ji,jj,jk,jn) - zltv(ji,jj,jk) + ptn(ji,jj,jk,jn) * z_hdivn
	343	END DO
	344	END DO
[503]	345	END DO
[2528]	346	CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zltv )
	347	ENDIF
	348	!
	349	ENDDO
[503]	350	!
[2715]	351	IF( wrk_not_released(3, 1,2,3,4,5,6) ) CALL ctl_stop('tra_adv_ubs: failed to release workspace arrays')
	352	!
[3211]	353
	354	!! * Reset control of array index permutation
	355	!FTRANS CLEAR
	356	# include "oce_ftrans.h90"
	357	# include "dom_oce_ftrans.h90"
	358	# include "trc_oce_ftrans.h90"
	359
[2528]	360	END SUBROUTINE tra_adv_ubs
[503]	361
	362
[2528]	363	SUBROUTINE nonosc_z( pbef, pcc, paft, p2dt )
[503]	364	!!---------------------------------------------------------------------
	365	!! * ROUTINE nonosc_z *
	366	!!
	367	!! ** Purpose : compute monotonic tracer fluxes from the upstream
	368	!! scheme and the before field by a nonoscillatory algorithm
	369	!!
	370	!! ** Method : ... ???
	371	!! warning : pbef and paft must be masked, but the boundaries
	372	!! conditions on the fluxes are not necessary zalezak (1979)
	373	!! drange (1995) multi-dimensional forward-in-time and upstream-
	374	!! in-space based differencing for fluid
	375	!!----------------------------------------------------------------------
[2715]	376	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	377	USE wrk_nemo, ONLY: zbetup => wrk_3d_1, zbetdo => wrk_3d_2 ! 3D workspace
[3211]	378
	379	!! DCSE_NEMO: need additional directives for renamed module variables
	380	!FTRANS zbetup zbetdo :I :I :z
	381
[2715]	382	!
[4409]	383	REAL(wp), INTENT(in ), DIMENSION(jpkorig) :: p2dt ! vertical profile of tracer time-step
[3211]	384
	385	!! DCSE_NEMO: This style defeats ftrans
	386	! REAL(wp), DIMENSION (jpi,jpj,jpk) :: pbef ! before field
	387	! REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: paft ! after field
	388	! REAL(wp), INTENT(inout), DIMENSION (jpi,jpj,jpk) :: pcc ! monotonic flux in the k direction
	389
	390	!FTRANS pbef paft pcc :I :I :z
[4409]	391	REAL(wp), INTENT(inout) :: pbef(jpi,jpj,jpkorig) ! before field
	392	REAL(wp), INTENT(inout) :: paft(jpi,jpj,jpkorig) ! after field
	393	REAL(wp), INTENT(inout) :: pcc(jpi,jpj,jpkorig) ! monotonic flux in the k direction
[2715]	394	!
	395	INTEGER :: ji, jj, jk ! dummy loop indices
	396	INTEGER :: ikm1 ! local integer
	397	REAL(wp) :: zpos, zneg, zbt, za, zb, zc, zbig, zrtrn, z2dtt ! local scalars
[503]	398	!!----------------------------------------------------------------------
	399
[2715]	400	IF( wrk_in_use(3, 1,2) ) THEN
	401	CALL ctl_stop('nonosc_z: requested workspace arrays unavailable') ; RETURN
	402	ENDIF
[503]	403
[2715]	404	zbig = 1.e+40_wp
	405	zrtrn = 1.e-15_wp
	406	zbetup(:,:,:) = 0._wp ; zbetdo(:,:,:) = 0._wp
	407
[503]	408	! Search local extrema
	409	! --------------------
	410	! large negative value (-zbig) inside land
	411	pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
	412	paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) - zbig * ( 1.e0 - tmask(:,:,:) )
	413	! search maximum in neighbourhood
[3211]	414	#if defined key_z_first
	415	DO jj = 2, jpjm1
	416	DO ji = 2, jpim1
	417	DO jk = 1, jpkm1
	418	ikm1 = MAX(jk-1,1)
	419	#else
[503]	420	DO jk = 1, jpkm1
	421	ikm1 = MAX(jk-1,1)
	422	DO jj = 2, jpjm1
	423	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	424	#endif
[503]	425	zbetup(ji,jj,jk) = MAX( pbef(ji ,jj ,jk ), paft(ji ,jj ,jk ), &
	426	& pbef(ji ,jj ,ikm1), pbef(ji ,jj ,jk+1), &
	427	& paft(ji ,jj ,ikm1), paft(ji ,jj ,jk+1) )
	428	END DO
	429	END DO
	430	END DO
	431	! large positive value (+zbig) inside land
	432	pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
	433	paft(:,:,:) = paft(:,:,:) * tmask(:,:,:) + zbig * ( 1.e0 - tmask(:,:,:) )
	434	! search minimum in neighbourhood
[3211]	435	#if defined key_z_first
	436	DO jj = 2, jpjm1
	437	DO ji = 2, jpim1
	438	DO jk = 1, jpkm1
	439	ikm1 = MAX(jk-1,1)
	440	#else
[503]	441	DO jk = 1, jpkm1
	442	ikm1 = MAX(jk-1,1)
	443	DO jj = 2, jpjm1
	444	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	445	#endif
[503]	446	zbetdo(ji,jj,jk) = MIN( pbef(ji ,jj ,jk ), paft(ji ,jj ,jk ), &
	447	& pbef(ji ,jj ,ikm1), pbef(ji ,jj ,jk+1), &
	448	& paft(ji ,jj ,ikm1), paft(ji ,jj ,jk+1) )
	449	END DO
	450	END DO
	451	END DO
	452
	453	! restore masked values to zero
	454	pbef(:,:,:) = pbef(:,:,:) * tmask(:,:,:)
	455	paft(:,:,:) = paft(:,:,:) * tmask(:,:,:)
	456
[2528]	457
[503]	458	! 2. Positive and negative part of fluxes and beta terms
	459	! ------------------------------------------------------
	460
[3211]	461	#if defined key_z_first
	462	DO jj = 2, jpjm1
	463	DO ji = 2, jpim1
	464	DO jk = 1, jpkm1
	465	z2dtt = p2dt(jk)
	466	#else
[503]	467	DO jk = 1, jpkm1
[2528]	468	z2dtt = p2dt(jk)
[503]	469	DO jj = 2, jpjm1
	470	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	471	#endif
[503]	472	! positive & negative part of the flux
	473	zpos = MAX( 0., pcc(ji ,jj ,jk+1) ) - MIN( 0., pcc(ji ,jj ,jk ) )
	474	zneg = MAX( 0., pcc(ji ,jj ,jk ) ) - MIN( 0., pcc(ji ,jj ,jk+1) )
	475	! up & down beta terms
	476	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) / z2dtt
	477	zbetup(ji,jj,jk) = ( zbetup(ji,jj,jk) - paft(ji,jj,jk) ) / (zpos+zrtrn) * zbt
	478	zbetdo(ji,jj,jk) = ( paft(ji,jj,jk) - zbetdo(ji,jj,jk) ) / (zneg+zrtrn) * zbt
	479	END DO
	480	END DO
	481	END DO
	482	! monotonic flux in the k direction, i.e. pcc
	483	! -------------------------------------------
[3211]	484	#if defined key_z_first
	485	DO jj = 2, jpjm1
	486	DO ji = 2, jpim1
	487	DO jk = 2, jpkm1
	488	#else
[503]	489	DO jk = 2, jpkm1
	490	DO jj = 2, jpjm1
	491	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	492	#endif
[503]	493	za = MIN( 1., zbetdo(ji,jj,jk), zbetup(ji,jj,jk-1) )
	494	zb = MIN( 1., zbetup(ji,jj,jk), zbetdo(ji,jj,jk-1) )
	495	zc = 0.5 * ( 1.e0 + SIGN( 1.e0, pcc(ji,jj,jk) ) )
	496	pcc(ji,jj,jk) = pcc(ji,jj,jk) * ( zc * za + ( 1.e0 - zc) * zb )
	497	END DO
	498	END DO
	499	END DO
	500	!
[2715]	501	IF( wrk_not_released(3, 1,2) ) CALL ctl_stop('nonosc_z: failed to release workspace arrays')
	502	!
[503]	503	END SUBROUTINE nonosc_z
	504
	505	!!======================================================================
	506	END MODULE traadv_ubs

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

Download in other formats: