Context Navigation

ldfslp.F90 @ 13095

Last change on this file since 13095 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

Property svn:keywords set to Id

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

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

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

Download in other formats: