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.

ldfslp.F90 in NEMO/trunk/src/OCE/LDF – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/LDF/ldfslp.F90 @ 13295

Last change on this file since 13295 was 13295, checked in by acc, 4 years ago
Replace do-loop macros in the trunk with alternative forms with greater flexibility for extra halo applications. This alters a lot of routines but does not change any behaviour or results. do_loop_substitute.h90 is greatly simplified by this change. SETTE results are identical to those with the previous revision
Property svn:keywords set to `Id`
File size: 47.2 KB

Rev	Line
[3]	1	MODULE ldfslp
	2	!!======================================================================
	3	!! * MODULE ldfslp *
	4	!! Ocean physics: slopes of neutral surfaces
	5	!!======================================================================
[1515]	6	!! History : OPA ! 1994-12 (G. Madec, M. Imbard) Original code
	7	!! 8.0 ! 1997-06 (G. Madec) optimization, lbc
	8	!! 8.1 ! 1999-10 (A. Jouzeau) NEW profile in the mixed layer
[2528]	9	!! NEMO 1.0 ! 2002-10 (G. Madec) Free form, F90
	10	!! - ! 2005-10 (A. Beckmann) correction for s-coordinates
	11	!! 3.3 ! 2010-10 (G. Nurser, C. Harris, G. Madec) add Griffies operator
	12	!! - ! 2010-11 (F. Dupond, G. Madec) bug correction in slopes just below the ML
[5836]	13	!! 3.7 ! 2013-12 (F. Lemarie, G. Madec) add limiter on triad slopes
[1515]	14	!!----------------------------------------------------------------------
[5836]	15
[3]	16	!!----------------------------------------------------------------------
[3625]	17	!! ldf_slp : calculates the slopes of neutral surface (Madec operator)
[5836]	18	!! ldf_slp_triad : calculates the triads of isoneutral slopes (Griffies operator)
[3625]	19	!! ldf_slp_mxl : calculates the slopes at the base of the mixed layer (Madec operator)
	20	!! ldf_slp_init : initialization of the slopes computation
[3]	21	!!----------------------------------------------------------------------
[3848]	22	USE oce ! ocean dynamics and tracers
[12377]	23	USE isf_oce ! ice shelf
[3848]	24	USE dom_oce ! ocean space and time domain
[9490]	25	! USE ldfdyn ! lateral diffusion: eddy viscosity coef.
[3848]	26	USE phycst ! physical constants
	27	USE zdfmxl ! mixed layer depth
	28	USE eosbn2 ! equation of states
[4990]	29	!
	30	USE in_out_manager ! I/O manager
[5836]	31	USE prtctl ! Print control
[3848]	32	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[5836]	33	USE lib_mpp ! distribued memory computing library
	34	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
[3848]	35	USE timing ! Timing
[3]	36
	37	IMPLICIT NONE
	38	PRIVATE
	39
[5836]	40	PUBLIC ldf_slp ! routine called by step.F90
	41	PUBLIC ldf_slp_triad ! routine called by step.F90
	42	PUBLIC ldf_slp_init ! routine called by nemogcm.F90
[3]	43
[5836]	44	LOGICAL , PUBLIC :: l_ldfslp = .FALSE. !: slopes flag
	45
[9490]	46	LOGICAL , PUBLIC :: ln_traldf_iso = .TRUE. !: iso-neutral direction (nam_traldf namelist)
	47	LOGICAL , PUBLIC :: ln_traldf_triad = .FALSE. !: griffies triad scheme (nam_traldf namelist)
	48	LOGICAL , PUBLIC :: ln_dynldf_iso !: iso-neutral direction (nam_dynldf namelist)
[5836]	49
[9490]	50	LOGICAL , PUBLIC :: ln_triad_iso = .FALSE. !: pure horizontal mixing in ML (nam_traldf namelist)
	51	LOGICAL , PUBLIC :: ln_botmix_triad = .FALSE. !: mixing on bottom (nam_traldf namelist)
	52	REAL(wp), PUBLIC :: rn_sw_triad = 1._wp !: =1 switching triads ; =0 all four triads used (nam_traldf namelist)
	53	REAL(wp), PUBLIC :: rn_slpmax = 0.01_wp !: slope limit (nam_traldf namelist)
[5836]	54
	55	LOGICAL , PUBLIC :: l_grad_zps = .FALSE. !: special treatment for Horz Tgradients w partial steps (triad operator)
	56
	57	! !! Classic operator (Madec)
[2715]	58	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: uslp, wslpi !: i_slope at U- and W-points
	59	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vslp, wslpj !: j-slope at V- and W-points
[5836]	60	! !! triad operator (Griffies)
[2715]	61	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wslp2 !: wslp**2 from Griffies quarter cells
[3294]	62	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi_g, triadj_g !: skew flux slopes relative to geopotentials
[2715]	63	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi , triadj !: isoneutral slopes relative to model-coordinate
[5836]	64	! !! both operators
	65	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ah_wslp2 !: ah * slope^2 at w-point
	66	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: akz !: stabilizing vertical diffusivity
	67
	68	! !! Madec operator
[2715]	69	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: omlmask ! mask of the surface mixed layer at T-pt
	70	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: uslpml, wslpiml ! i_slope at U- and W-points just below the mixed layer
	71	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: vslpml, wslpjml ! j_slope at V- and W-points just below the mixed layer
[2528]	72
	73	REAL(wp) :: repsln = 1.e-25_wp ! tiny value used as minium of di(rho), dj(rho) and dk(rho)
	74
[3]	75	!! * Substitutions
[12377]	76	# include "do_loop_substitute.h90"
[13237]	77	# include "domzgr_substitute.h90"
[3]	78	!!----------------------------------------------------------------------
[9598]	79	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[1156]	80	!! $Id$
[10068]	81	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	82	!!----------------------------------------------------------------------
	83	CONTAINS
	84
[12377]	85	SUBROUTINE ldf_slp( kt, prd, pn2, Kbb, Kmm )
[3]	86	!!----------------------------------------------------------------------
	87	!! * ROUTINE ldf_slp *
[3294]	88	!!
[1515]	89	!! ** Purpose : Compute the slopes of neutral surface (slope of isopycnal
[2528]	90	!! surfaces referenced locally) (ln_traldf_iso=T).
[3294]	91	!!
	92	!! ** Method : The slope in the i-direction is computed at U- and
	93	!! W-points (uslp, wslpi) and the slope in the j-direction is
[3]	94	!! computed at V- and W-points (vslp, wslpj).
	95	!! They are bounded by 1/100 over the whole ocean, and within the
	96	!! surface layer they are bounded by the distance to the surface
	97	!! ( slope<= depth/l where l is the length scale of horizontal
	98	!! diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
	99	!! of 10cm/s)
	100	!! A horizontal shapiro filter is applied to the slopes
[461]	101	!! ln_sco=T, s-coordinate, add to the previously computed slopes
[3]	102	!! the slope of the model level surface.
	103	!! macro-tasked on horizontal slab (jk-loop) (2, jpk-1)
	104	!! [slopes already set to zero at level 1, and to zero or the ocean
[461]	105	!! bottom slope (ln_sco=T) at level jpk in inildf]
[3]	106	!!
[3294]	107	!! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and j-slopes
[3]	108	!! of now neutral surfaces at u-, w- and v- w-points, resp.
[1515]	109	!!----------------------------------------------------------------------
[2715]	110	INTEGER , INTENT(in) :: kt ! ocean time-step index
[12377]	111	INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices
[2715]	112	REAL(wp), INTENT(in), DIMENSION(:,:,:) :: prd ! in situ density
	113	REAL(wp), INTENT(in), DIMENSION(:,:,:) :: pn2 ! Brunt-Vaisala frequency (locally ref.)
[1515]	114	!!
	115	INTEGER :: ji , jj , jk ! dummy loop indices
[6140]	116	INTEGER :: ii0, ii1 ! temporary integer
	117	INTEGER :: ij0, ij1 ! temporary integer
[5836]	118	REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16, zcofw, z1_slpmax ! local scalars
[2528]	119	REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - -
	120	REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - -
	121	REAL(wp) :: zck, zfk, zbw ! - -
[6140]	122	REAL(wp) :: zdepu, zdepv ! - -
[9019]	123	REAL(wp), DIMENSION(jpi,jpj) :: zslpml_hmlpu, zslpml_hmlpv
	124	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgru, zwz, zdzr
	125	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zgrv, zww
[3]	126	!!----------------------------------------------------------------------
[3294]	127	!
[9019]	128	IF( ln_timing ) CALL timing_start('ldf_slp')
[3294]	129	!
[5836]	130	zeps = 1.e-20_wp !== Local constant initialization ==!
	131	z1_16 = 1.0_wp / 16._wp
	132	zm1_g = -1.0_wp / grav
	133	zm1_2g = -0.5_wp / grav
	134	z1_slpmax = 1._wp / rn_slpmax
	135	!
[7753]	136	zww(:,:,:) = 0._wp
	137	zwz(:,:,:) = 0._wp
[5836]	138	!
[13295]	139	DO_3D( 1, 0, 1, 0, 1, jpk )
[12377]	140	zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) )
	141	zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) )
	142	END_3D
[5836]	143	IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level
[13295]	144	DO_2D( 1, 0, 1, 0 )
[12377]	145	zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj)
	146	zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj)
	147	END_2D
[5836]	148	ENDIF
[6140]	149	IF( ln_zps .AND. ln_isfcav ) THEN ! partial steps correction at the bottom ocean level
[13295]	150	DO_2D( 1, 0, 1, 0 )
[12377]	151	IF( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj)
	152	IF( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj)
	153	END_2D
[6140]	154	ENDIF
[5836]	155	!
[7753]	156	zdzr(:,:,1) = 0._wp !== Local vertical density gradient at T-point == ! (evaluated from N^2)
[5836]	157	DO jk = 2, jpkm1
	158	! ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
	159	! ! trick: tmask(ik ) = 0 => all pn2 = 0 => zdzr = 0
	160	! ! else tmask(ik+1) = 0 => pn2(ik+1) = 0 => zdzr divides by 1
	161	! ! umask(ik+1) /= 0 => all pn2 /= 0 => zdzr divides by 2
	162	! ! NB: 1/(tmask+1) = (1-.5tmask) substitute a / by a ==> faster
[7753]	163	zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp ) &
	164	& * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) )
[5836]	165	END DO
	166	!
	167	! !== Slopes just below the mixed layer ==!
[12377]	168	CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr, Kmm ) ! output: uslpml, vslpml, wslpiml, wslpjml
[2389]	169
[3294]	170
[5836]	171	! I. slopes at u and v point \| uslp = d/di( prd ) / d/dz( prd )
	172	! =========================== \| vslp = d/dj( prd ) / d/dz( prd )
	173	!
[6140]	174	IF ( ln_isfcav ) THEN
[13295]	175	DO_2D( 0, 0, 0, 0 )
[12377]	176	zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / ( MAX(hmlpt (ji,jj), hmlpt (ji+1,jj ), 5._wp) &
	177	& - MAX(risfdep(ji,jj), risfdep(ji+1,jj ) ) )
	178	zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / ( MAX(hmlpt (ji,jj), hmlpt (ji ,jj+1), 5._wp) &
	179	& - MAX(risfdep(ji,jj), risfdep(ji ,jj+1) ) )
	180	END_2D
[6140]	181	ELSE
[13295]	182	DO_2D( 0, 0, 0, 0 )
[12377]	183	zslpml_hmlpu(ji,jj) = uslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp)
	184	zslpml_hmlpv(ji,jj) = vslpml(ji,jj) / MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp)
	185	END_2D
[6140]	186	END IF
	187
[13295]	188	DO_3D( 0, 0, 0, 0, 2, jpkm1 )
[12377]	189	! ! horizontal and vertical density gradient at u- and v-points
	190	zau = zgru(ji,jj,jk) * r1_e1u(ji,jj)
	191	zav = zgrv(ji,jj,jk) * r1_e2v(ji,jj)
	192	zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) )
	193	zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) )
	194	! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
	195	! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
	196	zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u(ji,jj,jk,Kmm)* ABS( zau ) )
	197	zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v(ji,jj,jk,Kmm)* ABS( zav ) )
	198	! ! Fred Dupont: add a correction for bottom partial steps:
	199	! ! max slope = 1/2 * e3 / e1
	200	IF (ln_zps .AND. jk==mbku(ji,jj)) &
[13237]	201	zbu = MIN( zbu, - z1_slpmax * ABS( zau ) , &
	202	& - 2._wp * e1u(ji,jj) / e3u(ji,jj,jk,Kmm)* ABS( zau ) )
[12377]	203	IF (ln_zps .AND. jk==mbkv(ji,jj)) &
[13237]	204	zbv = MIN( zbv, - z1_slpmax * ABS( zav ) , &
	205	& - 2._wp * e2v(ji,jj) / e3v(ji,jj,jk,Kmm)* ABS( zav ) )
[12377]	206	! ! uslp and vslp output in zwz and zww, resp.
	207	zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj,jk) )
	208	zfj = MAX( omlmask(ji,jj,jk), omlmask(ji,jj+1,jk) )
	209	! thickness of water column between surface and level k at u/v point
	210	zdepu = 0.5_wp * ( ( gdept (ji,jj,jk,Kmm) + gdept (ji+1,jj,jk,Kmm) ) &
[13237]	211	& - 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj) ) &
	212	& - e3u(ji,jj,miku(ji,jj),Kmm) )
[12377]	213	zdepv = 0.5_wp * ( ( gdept (ji,jj,jk,Kmm) + gdept (ji,jj+1,jk,Kmm) ) &
[13237]	214	& - 2 * MAX( risfdep(ji,jj), risfdep(ji,jj+1) ) &
	215	& - e3v(ji,jj,mikv(ji,jj),Kmm) )
[12377]	216	!
	217	zwz(ji,jj,jk) = ( ( 1._wp - zfi) * zau / ( zbu - zeps ) &
	218	& + zfi * zdepu * zslpml_hmlpu(ji,jj) ) * umask(ji,jj,jk)
	219	zww(ji,jj,jk) = ( ( 1._wp - zfj) * zav / ( zbv - zeps ) &
	220	& + zfj * zdepv * zslpml_hmlpv(ji,jj) ) * vmask(ji,jj,jk)
[2528]	221	!!gm modif to suppress omlmask.... (as in Griffies case)
[5836]	222	! ! ! jk must be >= ML level for zf=1. otherwise zf=0.
	223	! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp )
	224	! zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp )
[12377]	225	! zci = 0.5 * ( gdept(ji+1,jj,jk,Kmm)+gdept(ji,jj,jk,Kmm) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )
	226	! zcj = 0.5 * ( gdept(ji,jj+1,jk,Kmm)+gdept(ji,jj,jk,Kmm) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )
[5836]	227	! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk)
	228	! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk)
[2528]	229	!!gm end modif
[12377]	230	END_3D
[13226]	231	CALL lbc_lnk_multi( 'ldfslp', zwz, 'U', -1.0_wp, zww, 'V', -1.0_wp ) ! lateral boundary conditions
[5836]	232	!
	233	! !* horizontal Shapiro filter
	234	DO jk = 2, jpkm1
[13295]	235	DO_2D( 0, 0, 0, 0 )
[12377]	236	uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) &
	237	& + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &
	238	& + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &
	239	& + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &
	240	& + 4.* zwz(ji ,jj ,jk) )
	241	vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) &
	242	& + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &
	243	& + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &
	244	& + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &
	245	& + 4.* zww(ji,jj ,jk) )
	246	END_2D
[5836]	247	DO jj = 3, jpj-2 ! other rows
[12377]	248	DO ji = 2, jpim1 ! vector opt.
[5836]	249	uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) &
	250	& + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &
	251	& + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &
	252	& + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &
	253	& + 4.* zwz(ji ,jj ,jk) )
	254	vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) &
	255	& + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &
	256	& + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &
	257	& + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &
	258	& + 4.* zww(ji,jj ,jk) )
[3]	259	END DO
[5836]	260	END DO
	261	! !* decrease along coastal boundaries
[13295]	262	DO_2D( 0, 0, 0, 0 )
[12377]	263	uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp &
	264	& * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp
	265	vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp &
	266	& * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp
	267	END_2D
[5836]	268	END DO
[3]	269
	270
[5836]	271	! II. slopes at w point \| wslpi = mij( d/di( prd ) / d/dz( prd )
	272	! =========================== \| wslpj = mij( d/dj( prd ) / d/dz( prd )
	273	!
[13295]	274	DO_3D( 0, 0, 0, 0, 2, jpkm1 )
[12377]	275	! !* Local vertical density gradient evaluated from N^2
	276	zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. )
	277	! !* Slopes at w point
	278	! ! i- & j-gradient of density at w-points
	279	zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) &
	280	& + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj)
	281	zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) &
	282	& + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj)
	283	zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) &
	284	& + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci * wmask (ji,jj,jk)
	285	zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) &
	286	& + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj * wmask (ji,jj,jk)
	287	! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
	288	! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
	289	zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/e3w(ji,jj,jk,Kmm)* ABS( zai ) )
	290	zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w(ji,jj,jk,Kmm)* ABS( zaj ) )
	291	! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.)
	292	zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0
	293	zck = ( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mikt(ji,jj),Kmm) ) / MAX( hmlp(ji,jj) - gdepw(ji,jj,mikt(ji,jj),Kmm), 10._wp )
	294	zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk)
	295	zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) + zck * wslpjml(ji,jj) * zfk ) * wmask(ji,jj,jk)
[2528]	296
	297	!!gm modif to suppress omlmask.... (as in Griffies operator)
[5836]	298	! ! ! jk must be >= ML level for zfk=1. otherwise zfk=0.
	299	! zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp )
[13237]	300	! zck = gdepw(ji,jj,jk,Kmm) / MAX( hmlp(ji,jj), 10. )
[5836]	301	! zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk)
	302	! zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk)
[2528]	303	!!gm end modif
[12377]	304	END_3D
[13226]	305	CALL lbc_lnk_multi( 'ldfslp', zwz, 'T', -1.0_wp, zww, 'T', -1.0_wp ) ! lateral boundary conditions
[5836]	306	!
	307	! !* horizontal Shapiro filter
	308	DO jk = 2, jpkm1
[13295]	309	DO_2D( 0, 0, 0, 0 )
[12377]	310	zcofw = wmask(ji,jj,jk) * z1_16
	311	wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) &
	312	& + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &
	313	& + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &
	314	& + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &
	315	& + 4.* zwz(ji ,jj ,jk) ) * zcofw
[3]	316
[12377]	317	wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) &
	318	& + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &
	319	& + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &
	320	& + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &
	321	& + 4.* zww(ji ,jj ,jk) ) * zcofw
	322	END_2D
[5836]	323	DO jj = 3, jpj-2 ! other rows
[12377]	324	DO ji = 2, jpim1 ! vector opt.
[5836]	325	zcofw = wmask(ji,jj,jk) * z1_16
	326	wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) &
	327	& + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &
	328	& + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &
	329	& + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &
	330	& + 4.* zwz(ji ,jj ,jk) ) * zcofw
[3]	331
[5836]	332	wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) &
	333	& + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &
	334	& + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &
	335	& + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &
	336	& + 4.* zww(ji ,jj ,jk) ) * zcofw
[3]	337	END DO
[5836]	338	END DO
	339	! !* decrease in vicinity of topography
[13295]	340	DO_2D( 0, 0, 0, 0 )
[12377]	341	zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) &
	342	& * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25
	343	wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck
	344	wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck
	345	END_2D
[5836]	346	END DO
[3294]	347
[5836]	348	! IV. Lateral boundary conditions
	349	! ===============================
[13226]	350	CALL lbc_lnk_multi( 'ldfslp', uslp , 'U', -1.0_wp , vslp , 'V', -1.0_wp , wslpi, 'W', -1.0_wp, wslpj, 'W', -1.0_wp )
[3]	351
[12377]	352	IF(sn_cfctl%l_prtctl) THEN
[5836]	353	CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk)
	354	CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
[49]	355	ENDIF
[5836]	356	!
[9019]	357	IF( ln_timing ) CALL timing_stop('ldf_slp')
[2715]	358	!
[3]	359	END SUBROUTINE ldf_slp
	360
[3294]	361
[12377]	362	SUBROUTINE ldf_slp_triad ( kt, Kbb, Kmm )
[2528]	363	!!----------------------------------------------------------------------
[5836]	364	!! * ROUTINE ldf_slp_triad *
[2528]	365	!!
	366	!! ** Purpose : Compute the squared slopes of neutral surfaces (slope
[5836]	367	!! of iso-pycnal surfaces referenced locally) (ln_traldf_triad=T)
[3294]	368	!! at W-points using the Griffies quarter-cells.
[2528]	369	!!
[3294]	370	!! ** Method : calculates alpha and beta at T-points
[2528]	371	!!
	372	!! ** Action : - triadi_g, triadj_g T-pts i- and j-slope triads relative to geopot. (used for eiv)
	373	!! - triadi , triadj T-pts i- and j-slope triads relative to model-coordinate
	374	!! - wslp2 squared slope of neutral surfaces at w-points.
	375	!!----------------------------------------------------------------------
[3294]	376	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[12377]	377	INTEGER , INTENT(in) :: Kbb, Kmm ! ocean time level indices
[3294]	378	!!
[2528]	379	INTEGER :: ji, jj, jk, jl, ip, jp, kp ! dummy loop indices
[3294]	380	INTEGER :: iku, ikv ! local integer
	381	REAL(wp) :: zfacti, zfactj ! local scalars
	382	REAL(wp) :: znot_thru_surface ! local scalars
[5836]	383	REAL(wp) :: zdit, zdis, zdkt, zbu, zbti, zisw
	384	REAL(wp) :: zdjt, zdjs, zdks, zbv, zbtj, zjsw
[3294]	385	REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_g_raw, zti_g_lim
	386	REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim
[2528]	387	REAL(wp) :: zdzrho_raw
[5836]	388	REAL(wp) :: zbeta0, ze3_e1, ze3_e2
[9019]	389	REAL(wp), DIMENSION(jpi,jpj) :: z1_mlbw
	390	REAL(wp), DIMENSION(jpi,jpj,jpk,0:1) :: zdxrho , zdyrho, zdzrho ! Horizontal and vertical density gradients
	391	REAL(wp), DIMENSION(jpi,jpj,0:1,0:1) :: zti_mlb, ztj_mlb ! for Griffies operator only
[2528]	392	!!----------------------------------------------------------------------
[3294]	393	!
[9019]	394	IF( ln_timing ) CALL timing_start('ldf_slp_triad')
[3294]	395	!
[2528]	396	!--------------------------------!
	397	! Some preliminary calculation !
	398	!--------------------------------!
	399	!
[3294]	400	DO jl = 0, 1 !== unmasked before density i- j-, k-gradients ==!
	401	!
	402	ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln)
[13295]	403	DO_3D( 1, 0, 1, 0, 1, jpkm1 )
[12377]	404	zdit = ( ts(ji+1,jj,jk,jp_tem,Kbb) - ts(ji,jj,jk,jp_tem,Kbb) ) ! i-gradient of T & S at u-point
	405	zdis = ( ts(ji+1,jj,jk,jp_sal,Kbb) - ts(ji,jj,jk,jp_sal,Kbb) )
	406	zdjt = ( ts(ji,jj+1,jk,jp_tem,Kbb) - ts(ji,jj,jk,jp_tem,Kbb) ) ! j-gradient of T & S at v-point
	407	zdjs = ( ts(ji,jj+1,jk,jp_sal,Kbb) - ts(ji,jj,jk,jp_sal,Kbb) )
	408	zdxrho_raw = ( - rab_b(ji+ip,jj ,jk,jp_tem) * zdit + rab_b(ji+ip,jj ,jk,jp_sal) * zdis ) * r1_e1u(ji,jj)
	409	zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) * r1_e2v(ji,jj)
	410	zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign
	411	zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw )
	412	END_3D
[3294]	413	!
[4990]	414	IF( ln_zps .AND. l_grad_zps ) THEN ! partial steps: correction of i- & j-grad on bottom
[13295]	415	DO_2D( 1, 0, 1, 0 )
[12377]	416	iku = mbku(ji,jj) ; ikv = mbkv(ji,jj) ! last ocean level (u- & v-points)
	417	zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature
	418	zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity
	419	zdxrho_raw = ( - rab_b(ji+ip,jj ,iku,jp_tem) * zdit + rab_b(ji+ip,jj ,iku,jp_sal) * zdis ) * r1_e1u(ji,jj)
	420	zdyrho_raw = ( - rab_b(ji ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji ,jj+jp,ikv,jp_sal) * zdjs ) * r1_e2v(ji,jj)
	421	zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign
	422	zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw )
	423	END_2D
[3294]	424	ENDIF
	425	!
[2528]	426	END DO
[3294]	427
	428	DO kp = 0, 1 !== unmasked before density i- j-, k-gradients ==!
[13295]	429	DO_3D( 1, 1, 1, 1, 1, jpkm1 )
[12377]	430	IF( jk+kp > 1 ) THEN ! k-gradient of T & S a jk+kp
	431	zdkt = ( ts(ji,jj,jk+kp-1,jp_tem,Kbb) - ts(ji,jj,jk+kp,jp_tem,Kbb) )
	432	zdks = ( ts(ji,jj,jk+kp-1,jp_sal,Kbb) - ts(ji,jj,jk+kp,jp_sal,Kbb) )
	433	ELSE
	434	zdkt = 0._wp ! 1st level gradient set to zero
	435	zdks = 0._wp
	436	ENDIF
	437	zdzrho_raw = ( - rab_b(ji,jj,jk ,jp_tem) * zdkt &
	438	& + rab_b(ji,jj,jk ,jp_sal) * zdks &
	439	& ) / e3w(ji,jj,jk+kp,Kmm)
	440	zdzrho(ji,jj,jk,kp) = - MIN( - repsln , zdzrho_raw ) ! force zdzrho >= repsln
	441	END_3D
[2528]	442	END DO
	443	!
[13295]	444	DO_2D( 1, 1, 1, 1 )
[12377]	445	jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth
	446	z1_mlbw(ji,jj) = 1._wp / gdepw(ji,jj,jk,Kmm)
	447	END_2D
[2528]	448	!
[3294]	449	! !== intialisations to zero ==!
[2528]	450	!
[3294]	451	wslp2 (:,:,:) = 0._wp ! wslp2 will be cumulated 3D field set to zero
	452	triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero
[2528]	453	triadj_g(:,:,1,:,:) = 0._wp ; triadj_g(:,:,jpk,:,:) = 0._wp
[3294]	454	!!gm _iso set to zero missing
	455	triadi (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero
	456	triadj (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp
	457
[2528]	458	!-------------------------------------!
	459	! Triads just below the Mixed Layer !
	460	!-------------------------------------!
	461	!
[3294]	462	DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base
	463	DO kp = 0, 1 ! with only the slope-max limit and MASKED
[13295]	464	DO_2D( 1, 0, 1, 0 )
[12377]	465	ip = jl ; jp = jl
	466	!
	467	jk = nmln(ji+ip,jj) + 1
	468	IF( jk > mbkt(ji+ip,jj) ) THEN ! ML reaches bottom
	469	zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp
	470	ELSE
	471	! Add s-coordinate slope at t-points (do this by subtracting gradient of depth)
	472	zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) &
	473	& - ( gdept(ji+1,jj,jk-kp,Kmm) - gdept(ji,jj,jk-kp,Kmm) ) * r1_e1u(ji,jj) ) * umask(ji,jj,jk)
	474	ze3_e1 = e3w(ji+ip,jj,jk-kp,Kmm) * r1_e1u(ji,jj)
	475	zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1 , ABS( zti_g_raw ) ), zti_g_raw )
	476	ENDIF
	477	!
	478	jk = nmln(ji,jj+jp) + 1
	479	IF( jk > mbkt(ji,jj+jp) ) THEN !ML reaches bottom
	480	ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp
	481	ELSE
	482	ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) &
	483	& - ( gdept(ji,jj+1,jk-kp,Kmm) - gdept(ji,jj,jk-kp,Kmm) ) / e2v(ji,jj) ) * vmask(ji,jj,jk)
	484	ze3_e2 = e3w(ji,jj+jp,jk-kp,Kmm) / e2v(ji,jj)
	485	ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2 , ABS( ztj_g_raw ) ), ztj_g_raw )
	486	ENDIF
	487	END_2D
[2528]	488	END DO
	489	END DO
	490
	491	!-------------------------------------!
	492	! Triads with surface limits !
	493	!-------------------------------------!
	494	!
[3294]	495	DO kp = 0, 1 ! k-index of triads
[2528]	496	DO jl = 0, 1
[3294]	497	ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes)
[2528]	498	DO jk = 1, jpkm1
[3294]	499	! Must mask contribution to slope from dz/dx at constant s for triads jk=1,kp=0 that poke up though ocean surface
	500	znot_thru_surface = REAL( 1-1/(jk+kp), wp ) !jk+kp=1,=0.; otherwise=1.0
[13295]	501	DO_2D( 1, 0, 1, 0 )
[12377]	502	!
	503	! Calculate slope relative to geopotentials used for GM skew fluxes
	504	! Add s-coordinate slope at t-points (do this by subtracting gradient of depth)
	505	! Limit by slope relative to geopotentials by rn_slpmax, and mask by psi-point
	506	! masked by umask taken at the level of dz(rho)
	507	!
	508	! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti)
	509	!
	510	zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked
	511	ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp)
	512	!
	513	! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface
	514	zti_coord = znot_thru_surface * ( gdept(ji+1,jj ,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e1u(ji,jj)
	515	ztj_coord = znot_thru_surface * ( gdept(ji ,jj+1,jk,Kmm) - gdept(ji,jj,jk,Kmm) ) * r1_e2v(ji,jj) ! unmasked
	516	zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces
	517	ztj_g_raw = ztj_raw - ztj_coord
	518	! additional limit required in bilaplacian case
	519	ze3_e1 = e3w(ji+ip,jj ,jk+kp,Kmm) * r1_e1u(ji,jj)
	520	ze3_e2 = e3w(ji ,jj+jp,jk+kp,Kmm) * r1_e2v(ji,jj)
	521	! NB: hard coded factor 5 (can be a namelist parameter...)
	522	zti_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e1, ABS( zti_g_raw ) ), zti_g_raw )
	523	ztj_g_lim = SIGN( MIN( rn_slpmax, 5.0_wp * ze3_e2, ABS( ztj_g_raw ) ), ztj_g_raw )
	524	!
	525	! Below ML use limited zti_g as is & mask
	526	! Inside ML replace by linearly reducing sx_mlb towards surface & mask
	527	!
	528	zfacti = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji+ip,jj)), wp ) ! k index of uppermost point(s) of triad is jk+kp-1
	529	zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp ) ! must be .ge. nmln(ji,jj) for zfact=1
	530	! ! otherwise zfact=0
	531	zti_g_lim = ( zfacti * zti_g_lim &
	532	& + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) &
	533	& * gdepw(ji+ip,jj,jk+kp,Kmm) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp)
	534	ztj_g_lim = ( zfactj * ztj_g_lim &
	535	& + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) &
	536	& * gdepw(ji,jj+jp,jk+kp,Kmm) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp)
	537	!
	538	triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim
	539	triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim
	540	!
	541	! Get coefficients of isoneutral diffusion tensor
	542	! 1. Utilise gradients relative to s-coordinate, so add t-point slopes (subtract depth gradients)
	543	! 2. We require that isoneutral diffusion gives no vertical buoyancy flux
	544	! i.e. 33 term = (real slope* 31, 13 terms)
	545	! To do this, retain limited sx**2 in vertical flux, but divide by real slope for 13/31 terms
	546	! Equivalent to tapering A_iso = sx_limited2/(real slope)2
	547	!
	548	zti_lim = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp) ! remove coordinate slope => relative to coordinate surfaces
	549	ztj_lim = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp)
	550	!
	551	IF( ln_triad_iso ) THEN
	552	zti_raw = zti_lim*zti_lim / zti_raw
	553	ztj_raw = ztj_lim*ztj_lim / ztj_raw
	554	zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw )
	555	ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw )
	556	zti_lim = zfacti * zti_lim + ( 1._wp - zfacti ) * zti_raw
	557	ztj_lim = zfactj * ztj_lim + ( 1._wp - zfactj ) * ztj_raw
	558	ENDIF
	559	! ! switching triad scheme
	560	zisw = (1._wp - rn_sw_triad ) + rn_sw_triad &
	561	& * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - ip ) * SIGN( 1._wp , zdxrho(ji+ip,jj,jk,1-ip) ) )
	562	zjsw = (1._wp - rn_sw_triad ) + rn_sw_triad &
	563	& * 2._wp * ABS( 0.5_wp - kp - ( 0.5_wp - jp ) * SIGN( 1._wp , zdyrho(ji,jj+jp,jk,1-jp) ) )
	564	!
	565	triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim * zisw
	566	triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim * zjsw
	567	!
	568	zbu = e1e2u(ji ,jj ) * e3u(ji ,jj ,jk ,Kmm)
	569	zbv = e1e2v(ji ,jj ) * e3v(ji ,jj ,jk ,Kmm)
	570	zbti = e1e2t(ji+ip,jj ) * e3w(ji+ip,jj ,jk+kp,Kmm)
	571	zbtj = e1e2t(ji ,jj+jp) * e3w(ji ,jj+jp,jk+kp,Kmm)
	572	!
	573	wslp2(ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim*zti_g_lim ! masked
	574	wslp2(ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim*ztj_g_lim
	575	END_2D
[2528]	576	END DO
	577	END DO
	578	END DO
	579	!
	580	wslp2(:,:,1) = 0._wp ! force the surface wslp to zero
[3294]	581
[13226]	582	CALL lbc_lnk( 'ldfslp', wslp2, 'W', 1.0_wp ) ! lateral boundary confition on wslp2 only ==>>> gm : necessary ? to be checked
[2528]	583	!
[9019]	584	IF( ln_timing ) CALL timing_stop('ldf_slp_triad')
[2715]	585	!
[5836]	586	END SUBROUTINE ldf_slp_triad
[2528]	587
	588
[12377]	589	SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr, Kmm )
[3]	590	!!----------------------------------------------------------------------
	591	!! * ROUTINE ldf_slp_mxl *
	592	!!
[3294]	593	!! ** Purpose : Compute the slopes of iso-neutral surface just below
[1515]	594	!! the mixed layer.
	595	!!
[2528]	596	!! ** Method : The slope in the i-direction is computed at u- & w-points
	597	!! (uslpml, wslpiml) and the slope in the j-direction is computed
	598	!! at v- and w-points (vslpml, wslpjml) with the same bounds as
	599	!! in ldf_slp.
[3]	600	!!
[2389]	601	!! ** Action : uslpml, wslpiml : i- & j-slopes of neutral surfaces
[3294]	602	!! vslpml, wslpjml just below the mixed layer
[2389]	603	!! omlmask : mixed layer mask
[1515]	604	!!----------------------------------------------------------------------
[2715]	605	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: prd ! in situ density
	606	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: pn2 ! Brunt-Vaisala frequency (locally ref.)
	607	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_gru, p_grv ! i- & j-gradient of density (u- & v-pts)
	608	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_dzr ! z-gradient of density (T-point)
[12377]	609	INTEGER , INTENT(in) :: Kmm ! ocean time level indices
[3]	610	!!
[3294]	611	INTEGER :: ji , jj , jk ! dummy loop indices
	612	INTEGER :: iku, ikv, ik, ikm1 ! local integers
[5836]	613	REAL(wp) :: zeps, zm1_g, zm1_2g, z1_slpmax ! local scalars
[2528]	614	REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - -
	615	REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - -
	616	REAL(wp) :: zck, zfk, zbw ! - -
[3]	617	!!----------------------------------------------------------------------
[3294]	618	!
[2528]	619	zeps = 1.e-20_wp !== Local constant initialization ==!
	620	zm1_g = -1.0_wp / grav
	621	zm1_2g = -0.5_wp / grav
[5836]	622	z1_slpmax = 1._wp / rn_slpmax
[1515]	623	!
[7753]	624	uslpml (1,:) = 0._wp ; uslpml (jpi,:) = 0._wp
	625	vslpml (1,:) = 0._wp ; vslpml (jpi,:) = 0._wp
	626	wslpiml(1,:) = 0._wp ; wslpiml(jpi,:) = 0._wp
	627	wslpjml(1,:) = 0._wp ; wslpjml(jpi,:) = 0._wp
[2528]	628	!
[3294]	629	! !== surface mixed layer mask !
[13295]	630	DO_3D( 1, 1, 1, 1, 1, jpk )
[12377]	631	ik = nmln(ji,jj) - 1
	632	IF( jk <= ik ) THEN ; omlmask(ji,jj,jk) = 1._wp
	633	ELSE ; omlmask(ji,jj,jk) = 0._wp
	634	ENDIF
	635	END_3D
[3]	636
	637
	638	! Slopes of isopycnal surfaces just before bottom of mixed layer
	639	! --------------------------------------------------------------
[1515]	640	! The slope are computed as in the 3D case.
	641	! A key point here is the definition of the mixed layer at u- and v-points.
	642	! It is assumed to be the maximum of the two neighbouring T-point mixed layer depth.
	643	! Otherwise, a n2 value inside the mixed layer can be involved in the computation
	644	! of the slope, resulting in a too steep diagnosed slope and thus a spurious eddy
	645	! induce velocity field near the base of the mixed layer.
[3]	646	!-----------------------------------------------------------------------
[1515]	647	!
[13295]	648	DO_2D( 0, 0, 0, 0 )
[12377]	649	! !== Slope at u- & v-points just below the Mixed Layer ==!
	650	!
	651	! !- vertical density gradient for u- and v-slopes (from dzr at T-point)
	652	iku = MIN( MAX( 1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1)
	653	ikv = MIN( MAX( 1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) !
	654	zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) )
	655	zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) )
	656	! !- horizontal density gradient at u- & v-points
	657	zau = p_gru(ji,jj,iku) * r1_e1u(ji,jj)
	658	zav = p_grv(ji,jj,ikv) * r1_e2v(ji,jj)
	659	! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0
	660	! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
	661	zbu = MIN( zbu , - z1_slpmax * ABS( zau ) , -7.e+3_wp/e3u(ji,jj,iku,Kmm)* ABS( zau ) )
	662	zbv = MIN( zbv , - z1_slpmax * ABS( zav ) , -7.e+3_wp/e3v(ji,jj,ikv,Kmm)* ABS( zav ) )
	663	! !- Slope at u- & v-points (uslpml, vslpml)
	664	uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku)
	665	vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv)
	666	!
	667	! !== i- & j-slopes at w-points just below the Mixed Layer ==!
	668	!
	669	ik = MIN( nmln(ji,jj) + 1, jpk )
	670	ikm1 = MAX( 1, ik-1 )
	671	! !- vertical density gradient for w-slope (from N^2)
	672	zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
	673	! !- horizontal density i- & j-gradient at w-points
	674	zci = MAX( umask(ji-1,jj,ik ) + umask(ji,jj,ik ) &
	675	& + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj)
	676	zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) &
	677	& + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj)
	678	zai = ( p_gru(ji-1,jj,ik ) + p_gru(ji,jj,ik) &
	679	& + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1 ) ) / zci * tmask(ji,jj,ik)
	680	zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) &
	681	& + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1) ) / zcj * tmask(ji,jj,ik)
	682	! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0.
	683	! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
	684	zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/e3w(ji,jj,ik,Kmm)* ABS( zai ) )
	685	zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/e3w(ji,jj,ik,Kmm)* ABS( zaj ) )
	686	! !- i- & j-slope at w-points (wslpiml, wslpjml)
	687	wslpiml(ji,jj) = zai / ( zbi - zeps ) * tmask (ji,jj,ik)
	688	wslpjml(ji,jj) = zaj / ( zbj - zeps ) * tmask (ji,jj,ik)
	689	END_2D
[3294]	690	!!gm this lbc_lnk should be useless....
[13226]	691	CALL lbc_lnk_multi( 'ldfslp', uslpml , 'U', -1.0_wp , vslpml , 'V', -1.0_wp , wslpiml, 'W', -1.0_wp , wslpjml, 'W', -1.0_wp )
[1515]	692	!
[3]	693	END SUBROUTINE ldf_slp_mxl
	694
	695
	696	SUBROUTINE ldf_slp_init
	697	!!----------------------------------------------------------------------
	698	!! * ROUTINE ldf_slp_init *
	699	!!
	700	!! ** Purpose : Initialization for the isopycnal slopes computation
	701	!!
[5836]	702	!! ** Method :
[3]	703	!!----------------------------------------------------------------------
	704	INTEGER :: ji, jj, jk ! dummy loop indices
[2528]	705	INTEGER :: ierr ! local integer
[3]	706	!!----------------------------------------------------------------------
[3294]	707	!
	708	IF(lwp) THEN
[3]	709	WRITE(numout,*)
[2528]	710	WRITE(numout,*) 'ldf_slp_init : direction of lateral mixing'
	711	WRITE(numout,*) '~~~~~~~~~~~~'
[3]	712	ENDIF
[5836]	713	!
	714	ALLOCATE( ah_wslp2(jpi,jpj,jpk) , akz(jpi,jpj,jpk) , STAT=ierr )
	715	IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate ah_slp2 or akz' )
	716	!
	717	IF( ln_traldf_triad ) THEN ! Griffies operator : triad of slopes
[9490]	718	IF(lwp) WRITE(numout,*) ' ==>>> triad) operator (Griffies)'
[5836]	719	ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , &
	720	& triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , &
	721	& wslp2 (jpi,jpj,jpk) , STAT=ierr )
	722	IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' )
[2528]	723	IF( ln_dynldf_iso ) CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' )
	724	!
	725	ELSE ! Madec operator : slopes at u-, v-, and w-points
[9490]	726	IF(lwp) WRITE(numout,*) ' ==>>> iso operator (Madec)'
[5836]	727	ALLOCATE( omlmask(jpi,jpj,jpk) , &
	728	& uslp(jpi,jpj,jpk) , uslpml(jpi,jpj) , wslpi(jpi,jpj,jpk) , wslpiml(jpi,jpj) , &
	729	& vslp(jpi,jpj,jpk) , vslpml(jpi,jpj) , wslpj(jpi,jpj,jpk) , wslpjml(jpi,jpj) , STAT=ierr )
[2715]	730	IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Madec operator slope ' )
[3]	731
[2528]	732	! Direction of lateral diffusion (tracers and/or momentum)
	733	! ------------------------------
[7753]	734	uslp (:,:,:) = 0._wp ; uslpml (:,:) = 0._wp ! set the slope to zero (even in s-coordinates)
	735	vslp (:,:,:) = 0._wp ; vslpml (:,:) = 0._wp
	736	wslpi(:,:,:) = 0._wp ; wslpiml(:,:) = 0._wp
	737	wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp
[3]	738
[5836]	739	!!gm I no longer understand this.....
[6140]	740	!!gm IF( (ln_traldf_hor .OR. ln_dynldf_hor) .AND. .NOT. (.NOT.ln_linssh .AND. ln_rstart) ) THEN
[5836]	741	! IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
	742	!
	743	! ! geopotential diffusion in s-coordinates on tracers and/or momentum
	744	! ! The slopes of s-surfaces are computed once (no call to ldfslp in step)
	745	! ! The slopes for momentum diffusion are i- or j- averaged of those on tracers
	746	!
	747	! ! set the slope of diffusion to the slope of s-surfaces
[6140]	748	! ! ( c a u t i o n : minus sign as dep has positive value )
[5836]	749	! DO jk = 1, jpk
	750	! DO jj = 2, jpjm1
[12377]	751	! DO ji = 2, jpim1 ! vector opt.
	752	! uslp (ji,jj,jk) = - ( gdept(ji+1,jj,jk,Kmm) - gdept(ji ,jj ,jk,Kmm) ) * r1_e1u(ji,jj) * umask(ji,jj,jk)
	753	! vslp (ji,jj,jk) = - ( gdept(ji,jj+1,jk,Kmm) - gdept(ji ,jj ,jk,Kmm) ) * r1_e2v(ji,jj) * vmask(ji,jj,jk)
	754	! wslpi(ji,jj,jk) = - ( gdepw(ji+1,jj,jk,Kmm) - gdepw(ji-1,jj,jk,Kmm) ) * r1_e1t(ji,jj) * wmask(ji,jj,jk) * 0.5
	755	! wslpj(ji,jj,jk) = - ( gdepw(ji,jj+1,jk,Kmm) - gdepw(ji,jj-1,jk,Kmm) ) * r1_e2t(ji,jj) * wmask(ji,jj,jk) * 0.5
[5836]	756	! END DO
	757	! END DO
	758	! END DO
[10425]	759	! CALL lbc_lnk_multi( 'ldfslp', uslp , 'U', -1. ; CALL lbc_lnk( 'ldfslp', vslp , 'V', -1., wslpi, 'W', -1., wslpj, 'W', -1. )
[5836]	760	!!gm ENDIF
[2715]	761	ENDIF
	762	!
[3]	763	END SUBROUTINE ldf_slp_init
	764
	765	!!======================================================================
	766	END MODULE ldfslp

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

Download in other formats: