Context Navigation

ldfslp.F90

Last change on this file was 11738, checked in by marc, 5 years ago
The Dr Hook changes from my perl code.
File size: 56.4 KB

Line
1	MODULE ldfslp
2	!!======================================================================
3	!! * MODULE ldfslp *
4	!! Ocean physics: slopes of neutral surfaces
5	!!======================================================================
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
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
13	!!----------------------------------------------------------------------
14	#if defined key_ldfslp \|\| defined key_esopa
15	!!----------------------------------------------------------------------
16	!! 'key_ldfslp' Rotation of lateral mixing tensor
17	!!----------------------------------------------------------------------
18	!! ldf_slp_grif : calculates the triads of isoneutral slopes (Griffies operator)
19	!! ldf_slp : calculates the slopes of neutral surface (Madec operator)
20	!! ldf_slp_mxl : calculates the slopes at the base of the mixed layer (Madec operator)
21	!! ldf_slp_init : initialization of the slopes computation
22	!!----------------------------------------------------------------------
23	USE oce ! ocean dynamics and tracers
24	USE dom_oce ! ocean space and time domain
25	USE ldftra_oce ! lateral diffusion: traceur
26	USE ldfdyn_oce ! lateral diffusion: dynamics
27	USE phycst ! physical constants
28	USE zdfmxl ! mixed layer depth
29	USE eosbn2 ! equation of states
30	!
31	USE in_out_manager ! I/O manager
32	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
33	USE prtctl ! Print control
34	USE wrk_nemo ! work arrays
35	USE timing ! Timing
36	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
37
38	USE yomhook, ONLY: lhook, dr_hook
39	USE parkind1, ONLY: jprb, jpim
40
41	IMPLICIT NONE
42	PRIVATE
43
44	PUBLIC ldf_slp ! routine called by step.F90
45	PUBLIC ldf_slp_grif ! routine called by step.F90
46	PUBLIC ldf_slp_init ! routine called by opa.F90
47
48	LOGICAL , PUBLIC, PARAMETER :: lk_ldfslp = .TRUE. !: slopes flag
49	! !! Madec operator
50	! Arrays allocated in ldf_slp_init() routine once we know whether we're using the Griffies or Madec operator
51	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: uslp, wslpi !: i_slope at U- and W-points
52	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: vslp, wslpj !: j-slope at V- and W-points
53	! !! Griffies operator
54	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: wslp2 !: wslp**2 from Griffies quarter cells
55	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi_g, triadj_g !: skew flux slopes relative to geopotentials
56	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:,:) :: triadi , triadj !: isoneutral slopes relative to model-coordinate
57
58	! !! Madec operator
59	! Arrays allocated in ldf_slp_init() routine once we know whether we're using the Griffies or Madec operator
60	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: omlmask ! mask of the surface mixed layer at T-pt
61	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: uslpml, wslpiml ! i_slope at U- and W-points just below the mixed layer
62	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: vslpml, wslpjml ! j_slope at V- and W-points just below the mixed layer
63
64	REAL(wp) :: repsln = 1.e-25_wp ! tiny value used as minium of di(rho), dj(rho) and dk(rho)
65
66	!! * Substitutions
67	# include "domzgr_substitute.h90"
68	# include "ldftra_substitute.h90"
69	# include "ldfeiv_substitute.h90"
70	# include "vectopt_loop_substitute.h90"
71	!!----------------------------------------------------------------------
72	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
73	!! $Id$
74	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
75	!!----------------------------------------------------------------------
76	CONTAINS
77
78	SUBROUTINE ldf_slp( kt, prd, pn2 )
79	!!----------------------------------------------------------------------
80	!! * ROUTINE ldf_slp *
81	!!
82	!! ** Purpose : Compute the slopes of neutral surface (slope of isopycnal
83	!! surfaces referenced locally) (ln_traldf_iso=T).
84	!!
85	!! ** Method : The slope in the i-direction is computed at U- and
86	!! W-points (uslp, wslpi) and the slope in the j-direction is
87	!! computed at V- and W-points (vslp, wslpj).
88	!! They are bounded by 1/100 over the whole ocean, and within the
89	!! surface layer they are bounded by the distance to the surface
90	!! ( slope<= depth/l where l is the length scale of horizontal
91	!! diffusion (here, aht=2000m2/s ==> l=20km with a typical velocity
92	!! of 10cm/s)
93	!! A horizontal shapiro filter is applied to the slopes
94	!! ln_sco=T, s-coordinate, add to the previously computed slopes
95	!! the slope of the model level surface.
96	!! macro-tasked on horizontal slab (jk-loop) (2, jpk-1)
97	!! [slopes already set to zero at level 1, and to zero or the ocean
98	!! bottom slope (ln_sco=T) at level jpk in inildf]
99	!!
100	!! ** Action : - uslp, wslpi, and vslp, wslpj, the i- and j-slopes
101	!! of now neutral surfaces at u-, w- and v- w-points, resp.
102	!!----------------------------------------------------------------------
103	INTEGER , INTENT(in) :: kt ! ocean time-step index
104	REAL(wp), INTENT(in), DIMENSION(:,:,:) :: prd ! in situ density
105	REAL(wp), INTENT(in), DIMENSION(:,:,:) :: pn2 ! Brunt-Vaisala frequency (locally ref.)
106	!!
107	INTEGER :: ji , jj , jk ! dummy loop indices
108	INTEGER :: ii0, ii1, iku ! temporary integer
109	INTEGER :: ij0, ij1, ikv ! temporary integer
110	REAL(wp) :: zeps, zm1_g, zm1_2g, z1_16, zcofw ! local scalars
111	REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - -
112	REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - -
113	REAL(wp) :: zck, zfk, zbw ! - -
114	REAL(wp) :: zdepv, zdepu ! - -
115	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwz, zww
116	REAL(wp), POINTER, DIMENSION(:,:,:) :: zdzr
117	REAL(wp), POINTER, DIMENSION(:,:,:) :: zgru, zgrv
118	REAL(wp), POINTER, DIMENSION(:,: ) :: zhmlpu, zhmlpv
119	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
120	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
121	REAL(KIND=jprb) :: zhook_handle
122
123	CHARACTER(LEN=*), PARAMETER :: RoutineName='LDF_SLP'
124
125	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
126
127	!!----------------------------------------------------------------------
128	!
129	IF( nn_timing == 1 ) CALL timing_start('ldf_slp')
130	!
131	CALL wrk_alloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv )
132	CALL wrk_alloc( jpi,jpj, zhmlpu, zhmlpv )
133
134	IF ( ln_traldf_iso .OR. ln_dynldf_iso ) THEN
135
136	zeps = 1.e-20_wp !== Local constant initialization ==!
137	z1_16 = 1.0_wp / 16._wp
138	zm1_g = -1.0_wp / grav
139	zm1_2g = -0.5_wp / grav
140	!
141	zww(:,:,:) = 0._wp
142	zwz(:,:,:) = 0._wp
143	!
144	DO jk = 1, jpk !== i- & j-gradient of density ==!
145	DO jj = 1, jpjm1
146	DO ji = 1, fs_jpim1 ! vector opt.
147	zgru(ji,jj,jk) = umask(ji,jj,jk) * ( prd(ji+1,jj ,jk) - prd(ji,jj,jk) )
148	zgrv(ji,jj,jk) = vmask(ji,jj,jk) * ( prd(ji ,jj+1,jk) - prd(ji,jj,jk) )
149	END DO
150	END DO
151	END DO
152	IF( ln_zps ) THEN ! partial steps correction at the bottom ocean level
153	DO jj = 1, jpjm1
154	DO ji = 1, jpim1
155	zgru(ji,jj,mbku(ji,jj)) = gru(ji,jj)
156	zgrv(ji,jj,mbkv(ji,jj)) = grv(ji,jj)
157	END DO
158	END DO
159	ENDIF
160	IF( ln_zps .AND. ln_isfcav ) THEN ! partial steps correction at the bottom ocean level
161	DO jj = 1, jpjm1
162	DO ji = 1, jpim1
163	IF ( miku(ji,jj) > 1 ) zgru(ji,jj,miku(ji,jj)) = grui(ji,jj)
164	IF ( mikv(ji,jj) > 1 ) zgrv(ji,jj,mikv(ji,jj)) = grvi(ji,jj)
165	END DO
166	END DO
167	ENDIF
168	!
169	!== Local vertical density gradient at T-point == ! (evaluated from N^2)
170	! interior value
171	DO jk = 2, jpkm1
172	! ! zdzr = d/dz(prd)= - ( prd ) / grav * mk(pn2) -- at t point
173	! ! trick: tmask(ik ) = 0 => all pn2 = 0 => zdzr = 0
174	! ! else tmask(ik+1) = 0 => pn2(ik+1) = 0 => zdzr divides by 1
175	! ! umask(ik+1) /= 0 => all pn2 /= 0 => zdzr divides by 2
176	! ! NB: 1/(tmask+1) = (1-.5tmask) substitute a / by a ==> faster
177	zdzr(:,:,jk) = zm1_g * ( prd(:,:,jk) + 1._wp ) &
178	& * ( pn2(:,:,jk) + pn2(:,:,jk+1) ) * ( 1._wp - 0.5_wp * tmask(:,:,jk+1) )
179	END DO
180	! surface initialisation
181	zdzr(:,:,1) = 0._wp
182	IF ( ln_isfcav ) THEN
183	! if isf need to overwrite the interior value at at the first ocean point
184	DO jj = 1, jpjm1
185	DO ji = 1, jpim1
186	zdzr(ji,jj,1:mikt(ji,jj)) = 0._wp
187	END DO
188	END DO
189	END IF
190	!
191	! !== Slopes just below the mixed layer ==!
192	CALL ldf_slp_mxl( prd, pn2, zgru, zgrv, zdzr ) ! output: uslpml, vslpml, wslpiml, wslpjml
193
194
195	! I. slopes at u and v point \| uslp = d/di( prd ) / d/dz( prd )
196	! =========================== \| vslp = d/dj( prd ) / d/dz( prd )
197	!
198	IF ( ln_isfcav ) THEN
199	DO jj = 2, jpjm1
200	DO ji = fs_2, fs_jpim1 ! vector opt.
201	zhmlpu(ji,jj) = ( MAX(hmlpt(ji,jj) , hmlpt (ji+1,jj ), 5._wp) &
202	& - MAX(risfdep(ji,jj), risfdep(ji+1,jj ) ) )
203	zhmlpv(ji,jj) = ( MAX(hmlpt (ji,jj), hmlpt (ji ,jj+1), 5._wp) &
204	& - MAX(risfdep(ji,jj), risfdep(ji ,jj+1) ) )
205	ENDDO
206	ENDDO
207	ELSE
208	DO jj = 2, jpjm1
209	DO ji = fs_2, fs_jpim1 ! vector opt.
210	zhmlpu(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji+1,jj ), 5._wp)
211	zhmlpv(ji,jj) = MAX(hmlpt(ji,jj), hmlpt(ji ,jj+1), 5._wp)
212	ENDDO
213	ENDDO
214	END IF
215	DO jk = 2, jpkm1 !* Slopes at u and v points
216	DO jj = 2, jpjm1
217	DO ji = fs_2, fs_jpim1 ! vector opt.
218	! ! horizontal and vertical density gradient at u- and v-points
219	zau = zgru(ji,jj,jk) / e1u(ji,jj)
220	zav = zgrv(ji,jj,jk) / e2v(ji,jj)
221	zbu = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji+1,jj ,jk) )
222	zbv = 0.5_wp * ( zdzr(ji,jj,jk) + zdzr(ji ,jj+1,jk) )
223	! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0
224	! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
225	zbu = MIN( zbu, -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,jk)* ABS( zau ) )
226	zbv = MIN( zbv, -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,jk)* ABS( zav ) )
227	! ! uslp and vslp output in zwz and zww, resp.
228	zfi = MAX( omlmask(ji,jj,jk), omlmask(ji+1,jj ,jk) )
229	zfj = MAX( omlmask(ji,jj,jk), omlmask(ji ,jj+1,jk) )
230	! thickness of water column between surface and level k at u/v point
231	zdepu = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji+1,jj ,jk) ) &
232	- 2 * MAX( risfdep(ji,jj), risfdep(ji+1,jj ) ) - fse3u(ji,jj,miku(ji,jj)) )
233	zdepv = 0.5_wp * ( ( fsdept(ji,jj,jk) + fsdept(ji ,jj+1,jk) ) &
234	- 2 * MAX( risfdep(ji,jj), risfdep(ji ,jj+1) ) - fse3v(ji,jj,mikv(ji,jj)) )
235	!
236	zwz(ji,jj,jk) = ( 1. - zfi) * zau / ( zbu - zeps ) &
237	& + zfi * uslpml(ji,jj) * zdepu / zhmlpu(ji,jj)
238	zwz(ji,jj,jk) = zwz(ji,jj,jk) * wumask(ji,jj,jk)
239	zww(ji,jj,jk) = ( 1. - zfj) * zav / ( zbv - zeps ) &
240	& + zfj * vslpml(ji,jj) * zdepv / zhmlpv(ji,jj)
241	zww(ji,jj,jk) = zww(ji,jj,jk) * wvmask(ji,jj,jk)
242
243
244	!!gm modif to suppress omlmask.... (as in Griffies case)
245	! ! ! jk must be >= ML level for zf=1. otherwise zf=0.
246	! zfi = REAL( 1 - 1/(1 + jk / MAX( nmln(ji+1,jj), nmln(ji,jj) ) ), wp )
247	! zfj = REAL( 1 - 1/(1 + jk / MAX( nmln(ji,jj+1), nmln(ji,jj) ) ), wp )
248	! zci = 0.5 * ( fsdept(ji+1,jj,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji+1,jj), 10. ) )
249	! zcj = 0.5 * ( fsdept(ji,jj+1,jk)+fsdept(ji,jj,jk) ) / MAX( hmlpt(ji,jj), hmlpt(ji,jj+1), 10. ) )
250	! zwz(ji,jj,jk) = ( zfi * zai / ( zbi - zeps ) + ( 1._wp - zfi ) * wslpiml(ji,jj) * zci ) * tmask(ji,jj,jk)
251	! zww(ji,jj,jk) = ( zfj * zaj / ( zbj - zeps ) + ( 1._wp - zfj ) * wslpjml(ji,jj) * zcj ) * tmask(ji,jj,jk)
252	!!gm end modif
253	END DO
254	END DO
255	END DO
256	CALL lbc_lnk( zwz, 'U', -1. ) ; CALL lbc_lnk( zww, 'V', -1. ) ! lateral boundary conditions
257	!
258	! !* horizontal Shapiro filter
259	DO jk = 2, jpkm1
260	DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only
261	DO ji = 2, jpim1
262	uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) &
263	& + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &
264	& + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &
265	& + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &
266	& + 4.* zwz(ji ,jj ,jk) )
267	vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) &
268	& + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &
269	& + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &
270	& + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &
271	& + 4.* zww(ji,jj ,jk) )
272	END DO
273	END DO
274	DO jj = 3, jpj-2 ! other rows
275	DO ji = fs_2, fs_jpim1 ! vector opt.
276	uslp(ji,jj,jk) = z1_16 * ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) &
277	& + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &
278	& + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &
279	& + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &
280	& + 4.* zwz(ji ,jj ,jk) )
281	vslp(ji,jj,jk) = z1_16 * ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) &
282	& + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &
283	& + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &
284	& + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &
285	& + 4.* zww(ji,jj ,jk) )
286	END DO
287	END DO
288	! !* decrease along coastal boundaries
289	DO jj = 2, jpjm1
290	DO ji = fs_2, fs_jpim1 ! vector opt.
291	uslp(ji,jj,jk) = uslp(ji,jj,jk) * ( umask(ji,jj+1,jk) + umask(ji,jj-1,jk ) ) * 0.5_wp &
292	& * ( umask(ji,jj ,jk) + umask(ji,jj ,jk+1) ) * 0.5_wp &
293	& * umask(ji,jj,jk-1)
294	vslp(ji,jj,jk) = vslp(ji,jj,jk) * ( vmask(ji+1,jj,jk) + vmask(ji-1,jj,jk ) ) * 0.5_wp &
295	& * ( vmask(ji ,jj,jk) + vmask(ji ,jj,jk+1) ) * 0.5_wp &
296	& * vmask(ji,jj,jk-1)
297	END DO
298	END DO
299	END DO
300
301
302	! II. slopes at w point \| wslpi = mij( d/di( prd ) / d/dz( prd )
303	! =========================== \| wslpj = mij( d/dj( prd ) / d/dz( prd )
304	!
305	DO jk = 2, jpkm1
306	DO jj = 2, jpjm1
307	DO ji = fs_2, fs_jpim1 ! vector opt.
308	! !* Local vertical density gradient evaluated from N^2
309	zbw = zm1_2g * pn2 (ji,jj,jk) * ( prd (ji,jj,jk) + prd (ji,jj,jk-1) + 2. ) * wmask(ji,jj,jk)
310	! !* Slopes at w point
311	! ! i- & j-gradient of density at w-points
312	zci = MAX( umask(ji-1,jj,jk ) + umask(ji,jj,jk ) &
313	& + umask(ji-1,jj,jk-1) + umask(ji,jj,jk-1) , zeps ) * e1t(ji,jj)
314	zcj = MAX( vmask(ji,jj-1,jk ) + vmask(ji,jj,jk-1) &
315	& + vmask(ji,jj-1,jk-1) + vmask(ji,jj,jk ) , zeps ) * e2t(ji,jj)
316	zai = ( zgru (ji-1,jj,jk ) + zgru (ji,jj,jk-1) &
317	& + zgru (ji-1,jj,jk-1) + zgru (ji,jj,jk ) ) / zci
318	zaj = ( zgrv (ji,jj-1,jk ) + zgrv (ji,jj,jk-1) &
319	& + zgrv (ji,jj-1,jk-1) + zgrv (ji,jj,jk ) ) / zcj
320	! ! bound the slopes: abs(zw.)<= 1/100 and zb..<0.
321	! ! + kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
322	zbi = MIN( zbw ,- 100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zai ) )
323	zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,jk)* ABS( zaj ) )
324	! ! wslpi and wslpj with ML flattening (output in zwz and zww, resp.)
325	zfk = MAX( omlmask(ji,jj,jk), omlmask(ji,jj,jk-1) ) ! zfk=1 in the ML otherwise zfk=0
326	zck = ( fsdepw(ji,jj,jk) - fsdepw(ji,jj,mikt(ji,jj) ) ) / MAX( hmlp(ji,jj), 10._wp )
327	zwz(ji,jj,jk) = ( zai / ( zbi - zeps ) * ( 1._wp - zfk ) &
328	& + zck * wslpiml(ji,jj) * zfk ) * wmask(ji,jj,jk)
329	zww(ji,jj,jk) = ( zaj / ( zbj - zeps ) * ( 1._wp - zfk ) &
330	& + zck * wslpjml(ji,jj) * zfk ) * wmask(ji,jj,jk)
331
332	!!gm modif to suppress omlmask.... (as in Griffies operator)
333	! ! ! jk must be >= ML level for zfk=1. otherwise zfk=0.
334	! zfk = REAL( 1 - 1/(1 + jk / nmln(ji+1,jj)), wp )
335	! zck = fsdepw(ji,jj,jk) / MAX( hmlp(ji,jj), 10. )
336	! zwz(ji,jj,jk) = ( zfk * zai / ( zbi - zeps ) + ( 1._wp - zfk ) * wslpiml(ji,jj) * zck ) * tmask(ji,jj,jk)
337	! zww(ji,jj,jk) = ( zfk * zaj / ( zbj - zeps ) + ( 1._wp - zfk ) * wslpjml(ji,jj) * zck ) * tmask(ji,jj,jk)
338	!!gm end modif
339	END DO
340	END DO
341	END DO
342	CALL lbc_lnk( zwz, 'T', -1. ) ; CALL lbc_lnk( zww, 'T', -1. ) ! lateral boundary conditions
343	!
344	! !* horizontal Shapiro filter
345	DO jk = 2, jpkm1
346	DO jj = 2, jpjm1, MAX(1, jpj-3) ! rows jj=2 and =jpjm1 only
347	DO ji = 2, jpim1
348	zcofw = tmask(ji,jj,jk) * z1_16
349	wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) &
350	& + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &
351	& + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &
352	& + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &
353	& + 4.* zwz(ji ,jj ,jk) ) * zcofw
354
355	wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) &
356	& + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &
357	& + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &
358	& + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &
359	& + 4.* zww(ji ,jj ,jk) ) * zcofw
360	END DO
361	END DO
362	DO jj = 3, jpj-2 ! other rows
363	DO ji = fs_2, fs_jpim1 ! vector opt.
364	zcofw = tmask(ji,jj,jk) * z1_16
365	wslpi(ji,jj,jk) = ( zwz(ji-1,jj-1,jk) + zwz(ji+1,jj-1,jk) &
366	& + zwz(ji-1,jj+1,jk) + zwz(ji+1,jj+1,jk) &
367	& + 2.*( zwz(ji ,jj-1,jk) + zwz(ji-1,jj ,jk) &
368	& + zwz(ji+1,jj ,jk) + zwz(ji ,jj+1,jk) ) &
369	& + 4.* zwz(ji ,jj ,jk) ) * zcofw
370
371	wslpj(ji,jj,jk) = ( zww(ji-1,jj-1,jk) + zww(ji+1,jj-1,jk) &
372	& + zww(ji-1,jj+1,jk) + zww(ji+1,jj+1,jk) &
373	& + 2.*( zww(ji ,jj-1,jk) + zww(ji-1,jj ,jk) &
374	& + zww(ji+1,jj ,jk) + zww(ji ,jj+1,jk) ) &
375	& + 4.* zww(ji ,jj ,jk) ) * zcofw
376	END DO
377	END DO
378	! !* decrease along coastal boundaries
379	DO jj = 2, jpjm1
380	DO ji = fs_2, fs_jpim1 ! vector opt.
381	zck = ( umask(ji,jj,jk) + umask(ji-1,jj,jk) ) &
382	& * ( vmask(ji,jj,jk) + vmask(ji,jj-1,jk) ) * 0.25
383	wslpi(ji,jj,jk) = wslpi(ji,jj,jk) * zck * wmask(ji,jj,jk)
384	wslpj(ji,jj,jk) = wslpj(ji,jj,jk) * zck * wmask(ji,jj,jk)
385	END DO
386	END DO
387	END DO
388
389	! III. Specific grid points
390	! ===========================
391	!
392	IF( cp_cfg == "orca" .AND. jp_cfg == 4 ) THEN ! ORCA_R4 configuration: horizontal diffusion in specific area
393	! ! Gibraltar Strait
394	ij0 = 50 ; ij1 = 53
395	ii0 = 69 ; ii1 = 71 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
396	ij0 = 51 ; ij1 = 53
397	ii0 = 68 ; ii1 = 71 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
398	ii0 = 69 ; ii1 = 71 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
399	ii0 = 69 ; ii1 = 71 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
400	!
401	! ! Mediterrannean Sea
402	ij0 = 49 ; ij1 = 56
403	ii0 = 71 ; ii1 = 90 ; uslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
404	ij0 = 50 ; ij1 = 56
405	ii0 = 70 ; ii1 = 90 ; vslp ( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
406	ii0 = 71 ; ii1 = 90 ; wslpi( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
407	ii0 = 71 ; ii1 = 90 ; wslpj( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) , : ) = 0._wp
408	ENDIF
409
410
411	! IV. Lateral boundary conditions
412	! ===============================
413	CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. )
414	CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. )
415
416
417	IF(ln_ctl) THEN
418	CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk)
419	CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
420	ENDIF
421	!
422
423	ELSEIF ( lk_vvl ) THEN
424
425	IF(lwp) THEN
426	WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
427	ENDIF
428
429	! geopotential diffusion in s-coordinates on tracers and/or momentum
430	! The slopes of s-surfaces are computed at each time step due to vvl
431	! The slopes for momentum diffusion are i- or j- averaged of those on tracers
432
433	! set the slope of diffusion to the slope of s-surfaces
434	! ( c a u t i o n : minus sign as fsdep has positive value )
435	DO jj = 2, jpjm1
436	DO ji = fs_2, fs_jpim1 ! vector opt.
437	uslp(ji,jj,1) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,1) - fsdept_b(ji ,jj ,1) ) * umask(ji,jj,1)
438	vslp(ji,jj,1) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,1) - fsdept_b(ji ,jj ,1) ) * vmask(ji,jj,1)
439	wslpi(ji,jj,1) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,1) - fsdepw_b(ji-1,jj,1) ) * tmask(ji,jj,1) * 0.5
440	wslpj(ji,jj,1) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,1) - fsdepw_b(ji,jj-1,1) ) * tmask(ji,jj,1) * 0.5
441	END DO
442	END DO
443
444	DO jk = 2, jpk
445	DO jj = 2, jpjm1
446	DO ji = fs_2, fs_jpim1 ! vector opt.
447	uslp(ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)
448	vslp(ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)
449	wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) &
450	& * wmask(ji,jj,jk) * 0.5
451	wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) &
452	& * wmask(ji,jj,jk) * 0.5
453	END DO
454	END DO
455	END DO
456
457	! Lateral boundary conditions on the slopes
458	CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. )
459	CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. )
460
461	if( kt == nit000 ) then
462	IF(lwp) WRITE(numout,) ' max slop: u',SQRT( MAXVAL(uslpuslp)), ' v ', SQRT(MAXVAL(vslp)), &
463	& ' wi', sqrt(MAXVAL(wslpi)), ' wj', sqrt(MAXVAL(wslpj))
464	endif
465
466	IF(ln_ctl) THEN
467	CALL prt_ctl(tab3d_1=uslp , clinfo1=' slp - u : ', tab3d_2=vslp, clinfo2=' v : ', kdim=jpk)
468	CALL prt_ctl(tab3d_1=wslpi, clinfo1=' slp - wi: ', tab3d_2=wslpj, clinfo2=' wj: ', kdim=jpk)
469	ENDIF
470
471	ENDIF
472
473	CALL wrk_dealloc( jpi,jpj,jpk, zwz, zww, zdzr, zgru, zgrv )
474	CALL wrk_dealloc( jpi,jpj, zhmlpu, zhmlpv)
475	!
476	IF( nn_timing == 1 ) CALL timing_stop('ldf_slp')
477	!
478	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
479	END SUBROUTINE ldf_slp
480
481
482	SUBROUTINE ldf_slp_grif ( kt )
483	!!----------------------------------------------------------------------
484	!! * ROUTINE ldf_slp_grif *
485	!!
486	!! ** Purpose : Compute the squared slopes of neutral surfaces (slope
487	!! of iso-pycnal surfaces referenced locally) (ln_traldf_grif=T)
488	!! at W-points using the Griffies quarter-cells.
489	!!
490	!! ** Method : calculates alpha and beta at T-points
491	!!
492	!! ** Action : - triadi_g, triadj_g T-pts i- and j-slope triads relative to geopot. (used for eiv)
493	!! - triadi , triadj T-pts i- and j-slope triads relative to model-coordinate
494	!! - wslp2 squared slope of neutral surfaces at w-points.
495	!!----------------------------------------------------------------------
496	INTEGER, INTENT( in ) :: kt ! ocean time-step index
497	!!
498	INTEGER :: ji, jj, jk, jl, ip, jp, kp ! dummy loop indices
499	INTEGER :: iku, ikv ! local integer
500	REAL(wp) :: zfacti, zfactj ! local scalars
501	REAL(wp) :: znot_thru_surface ! local scalars
502	REAL(wp) :: zdit, zdis, zdjt, zdjs, zdkt, zdks, zbu, zbv, zbti, zbtj
503	REAL(wp) :: zdxrho_raw, zti_coord, zti_raw, zti_lim, zti_g_raw, zti_g_lim
504	REAL(wp) :: zdyrho_raw, ztj_coord, ztj_raw, ztj_lim, ztj_g_raw, ztj_g_lim
505	REAL(wp) :: zdzrho_raw
506	REAL(wp), POINTER, DIMENSION(:,:) :: z1_mlbw
507	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zdxrho , zdyrho, zdzrho ! Horizontal and vertical density gradients
508	REAL(wp), POINTER, DIMENSION(:,:,:,:) :: zti_mlb, ztj_mlb ! for Griffies operator only
509	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
510	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
511	REAL(KIND=jprb) :: zhook_handle
512
513	CHARACTER(LEN=*), PARAMETER :: RoutineName='LDF_SLP_GRIF'
514
515	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
516
517	!!----------------------------------------------------------------------
518	!
519	IF( nn_timing == 1 ) CALL timing_start('ldf_slp_grif')
520	!
521	CALL wrk_alloc( jpi,jpj, z1_mlbw )
522	CALL wrk_alloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho, klstart = 0 )
523	CALL wrk_alloc( jpi,jpj, 2,2, zti_mlb, ztj_mlb, kkstart = 0, klstart = 0 )
524	!
525	!--------------------------------!
526	! Some preliminary calculation !
527	!--------------------------------!
528	!
529	DO jl = 0, 1 !== unmasked before density i- j-, k-gradients ==!
530	!
531	ip = jl ; jp = jl ! guaranteed nonzero gradients ( absolute value larger than repsln)
532	DO jk = 1, jpkm1 ! done each pair of triad
533	DO jj = 1, jpjm1 ! NB: not masked ==> a minimum value is set
534	DO ji = 1, fs_jpim1 ! vector opt.
535	zdit = ( tsb(ji+1,jj,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! i-gradient of T & S at u-point
536	zdis = ( tsb(ji+1,jj,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
537	zdjt = ( tsb(ji,jj+1,jk,jp_tem) - tsb(ji,jj,jk,jp_tem) ) ! j-gradient of T & S at v-point
538	zdjs = ( tsb(ji,jj+1,jk,jp_sal) - tsb(ji,jj,jk,jp_sal) )
539	zdxrho_raw = ( - rab_b(ji+ip,jj ,jk,jp_tem) * zdit + rab_b(ji+ip,jj ,jk,jp_sal) * zdis ) / e1u(ji,jj)
540	zdyrho_raw = ( - rab_b(ji ,jj+jp,jk,jp_tem) * zdjt + rab_b(ji ,jj+jp,jk,jp_sal) * zdjs ) / e2v(ji,jj)
541	zdxrho(ji+ip,jj ,jk,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign
542	zdyrho(ji ,jj+jp,jk,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw )
543	END DO
544	END DO
545	END DO
546	!
547	IF( ln_zps .AND. l_grad_zps ) THEN ! partial steps: correction of i- & j-grad on bottom
548	DO jj = 1, jpjm1
549	DO ji = 1, jpim1
550	iku = mbku(ji,jj) ; ikv = mbkv(ji,jj) ! last ocean level (u- & v-points)
551	zdit = gtsu(ji,jj,jp_tem) ; zdjt = gtsv(ji,jj,jp_tem) ! i- & j-gradient of Temperature
552	zdis = gtsu(ji,jj,jp_sal) ; zdjs = gtsv(ji,jj,jp_sal) ! i- & j-gradient of Salinity
553	zdxrho_raw = ( - rab_b(ji+ip,jj ,iku,jp_tem) * zdit + rab_b(ji+ip,jj ,iku,jp_sal) * zdis ) / e1u(ji,jj)
554	zdyrho_raw = ( - rab_b(ji ,jj+jp,ikv,jp_tem) * zdjt + rab_b(ji ,jj+jp,ikv,jp_sal) * zdjs ) / e2v(ji,jj)
555	zdxrho(ji+ip,jj ,iku,1-ip) = SIGN( MAX( repsln, ABS( zdxrho_raw ) ), zdxrho_raw ) ! keep the sign
556	zdyrho(ji ,jj+jp,ikv,1-jp) = SIGN( MAX( repsln, ABS( zdyrho_raw ) ), zdyrho_raw )
557	END DO
558	END DO
559	ENDIF
560	!
561	END DO
562
563	DO kp = 0, 1 !== unmasked before density i- j-, k-gradients ==!
564	DO jk = 1, jpkm1 ! done each pair of triad
565	DO jj = 1, jpj ! NB: not masked ==> a minimum value is set
566	DO ji = 1, jpi ! vector opt.
567	IF( jk+kp > 1 ) THEN ! k-gradient of T & S a jk+kp
568	zdkt = ( tsb(ji,jj,jk+kp-1,jp_tem) - tsb(ji,jj,jk+kp,jp_tem) )
569	zdks = ( tsb(ji,jj,jk+kp-1,jp_sal) - tsb(ji,jj,jk+kp,jp_sal) )
570	ELSE
571	zdkt = 0._wp ! 1st level gradient set to zero
572	zdks = 0._wp
573	ENDIF
574	zdzrho_raw = ( - rab_b(ji,jj,jk,jp_tem) * zdkt + rab_b(ji,jj,jk,jp_sal) * zdks ) / fse3w(ji,jj,jk+kp)
575	zdzrho(ji,jj,jk,kp) = - MIN( - repsln, zdzrho_raw ) ! force zdzrho >= repsln
576	END DO
577	END DO
578	END DO
579	END DO
580	!
581	DO jj = 1, jpj !== Reciprocal depth of the w-point below ML base ==!
582	DO ji = 1, jpi
583	jk = MIN( nmln(ji,jj), mbkt(ji,jj) ) + 1 ! MIN in case ML depth is the ocean depth
584	z1_mlbw(ji,jj) = 1._wp / fsdepw(ji,jj,jk)
585	END DO
586	END DO
587	!
588	! !== intialisations to zero ==!
589	!
590	wslp2 (:,:,:) = 0._wp ! wslp2 will be cumulated 3D field set to zero
591	triadi_g(:,:,1,:,:) = 0._wp ; triadi_g(:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero
592	triadj_g(:,:,1,:,:) = 0._wp ; triadj_g(:,:,jpk,:,:) = 0._wp
593	!!gm _iso set to zero missing
594	triadi (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp ! set surface and bottom slope to zero
595	triadj (:,:,1,:,:) = 0._wp ; triadj (:,:,jpk,:,:) = 0._wp
596
597	!-------------------------------------!
598	! Triads just below the Mixed Layer !
599	!-------------------------------------!
600	!
601	DO jl = 0, 1 ! calculate slope of the 4 triads immediately ONE level below mixed-layer base
602	DO kp = 0, 1 ! with only the slope-max limit and MASKED
603	DO jj = 1, jpjm1
604	DO ji = 1, fs_jpim1
605	ip = jl ; jp = jl
606	!
607	jk = nmln(ji+ip,jj) + 1
608	IF( jk .GT. mbkt(ji+ip,jj) ) THEN !ML reaches bottom
609	zti_mlb(ji+ip,jj ,1-ip,kp) = 0.0_wp
610	ELSE
611	! Add s-coordinate slope at t-points (do this by subtracting gradient of depth)
612	zti_g_raw = ( zdxrho(ji+ip,jj,jk-kp,1-ip) / zdzrho(ji+ip,jj,jk-kp,kp) &
613	& - ( fsdept(ji+1,jj,jk-kp) - fsdept(ji,jj,jk-kp) ) / e1u(ji,jj) ) * umask(ji,jj,jk)
614	zti_mlb(ji+ip,jj ,1-ip,kp) = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw )
615	ENDIF
616	!
617	jk = nmln(ji,jj+jp) + 1
618	IF( jk .GT. mbkt(ji,jj+jp) ) THEN !ML reaches bottom
619	ztj_mlb(ji ,jj+jp,1-jp,kp) = 0.0_wp
620	ELSE
621	ztj_g_raw = ( zdyrho(ji,jj+jp,jk-kp,1-jp) / zdzrho(ji,jj+jp,jk-kp,kp) &
622	& - ( fsdept(ji,jj+1,jk-kp) - fsdept(ji,jj,jk-kp) ) / e2v(ji,jj) ) * vmask(ji,jj,jk)
623	ztj_mlb(ji ,jj+jp,1-jp,kp) = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw )
624	ENDIF
625	END DO
626	END DO
627	END DO
628	END DO
629
630	!-------------------------------------!
631	! Triads with surface limits !
632	!-------------------------------------!
633	!
634	DO kp = 0, 1 ! k-index of triads
635	DO jl = 0, 1
636	ip = jl ; jp = jl ! i- and j-indices of triads (i-k and j-k planes)
637	DO jk = 1, jpkm1
638	! Must mask contribution to slope from dz/dx at constant s for triads jk=1,kp=0 that poke up though ocean surface
639	znot_thru_surface = REAL( 1-1/(jk+kp), wp ) !jk+kp=1,=0.; otherwise=1.0
640	DO jj = 1, jpjm1
641	DO ji = 1, fs_jpim1 ! vector opt.
642	!
643	! Calculate slope relative to geopotentials used for GM skew fluxes
644	! Add s-coordinate slope at t-points (do this by subtracting gradient of depth)
645	! Limit by slope relative to geopotentials by rn_slpmax, and mask by psi-point
646	! masked by umask taken at the level of dz(rho)
647	!
648	! raw slopes: unmasked unbounded slopes (relative to geopotential (zti_g) and model surface (zti)
649	!
650	zti_raw = zdxrho(ji+ip,jj ,jk,1-ip) / zdzrho(ji+ip,jj ,jk,kp) ! unmasked
651	ztj_raw = zdyrho(ji ,jj+jp,jk,1-jp) / zdzrho(ji ,jj+jp,jk,kp)
652
653	! Must mask contribution to slope for triad jk=1,kp=0 that poke up though ocean surface
654	zti_coord = znot_thru_surface * ( fsdept(ji+1,jj ,jk) - fsdept(ji,jj,jk) ) / e1u(ji,jj)
655	ztj_coord = znot_thru_surface * ( fsdept(ji ,jj+1,jk) - fsdept(ji,jj,jk) ) / e2v(ji,jj) ! unmasked
656	zti_g_raw = zti_raw - zti_coord ! ref to geopot surfaces
657	ztj_g_raw = ztj_raw - ztj_coord
658	zti_g_lim = SIGN( MIN( rn_slpmax, ABS( zti_g_raw ) ), zti_g_raw )
659	ztj_g_lim = SIGN( MIN( rn_slpmax, ABS( ztj_g_raw ) ), ztj_g_raw )
660	!
661	! Below ML use limited zti_g as is & mask
662	! Inside ML replace by linearly reducing sx_mlb towards surface & mask
663	!
664	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
665	zfactj = REAL( 1 - 1/(1 + (jk+kp-1)/nmln(ji,jj+jp)), wp ) ! must be .ge. nmln(ji,jj) for zfact=1
666	! ! otherwise zfact=0
667	zti_g_lim = ( zfacti * zti_g_lim &
668	& + ( 1._wp - zfacti ) * zti_mlb(ji+ip,jj,1-ip,kp) &
669	& * fsdepw(ji+ip,jj,jk+kp) * z1_mlbw(ji+ip,jj) ) * umask(ji,jj,jk+kp)
670	ztj_g_lim = ( zfactj * ztj_g_lim &
671	& + ( 1._wp - zfactj ) * ztj_mlb(ji,jj+jp,1-jp,kp) &
672	& * fsdepw(ji,jj+jp,jk+kp) * z1_mlbw(ji,jj+jp) ) * vmask(ji,jj,jk+kp)
673	!
674	triadi_g(ji+ip,jj ,jk,1-ip,kp) = zti_g_lim
675	triadj_g(ji ,jj+jp,jk,1-jp,kp) = ztj_g_lim
676	!
677	! Get coefficients of isoneutral diffusion tensor
678	! 1. Utilise gradients relative to s-coordinate, so add t-point slopes (subtract depth gradients)
679	! 2. We require that isoneutral diffusion gives no vertical buoyancy flux
680	! i.e. 33 term = (real slope* 31, 13 terms)
681	! To do this, retain limited sx**2 in vertical flux, but divide by real slope for 13/31 terms
682	! Equivalent to tapering A_iso = sx_limited2/(real slope)2
683	!
684	zti_lim = ( zti_g_lim + zti_coord ) * umask(ji,jj,jk+kp) ! remove coordinate slope => relative to coordinate surfaces
685	ztj_lim = ( ztj_g_lim + ztj_coord ) * vmask(ji,jj,jk+kp)
686	!
687	IF( ln_triad_iso ) THEN
688	zti_raw = zti_lim**2 / zti_raw
689	ztj_raw = ztj_lim**2 / ztj_raw
690	zti_raw = SIGN( MIN( ABS(zti_lim), ABS( zti_raw ) ), zti_raw )
691	ztj_raw = SIGN( MIN( ABS(ztj_lim), ABS( ztj_raw ) ), ztj_raw )
692	zti_lim = zfacti * zti_lim &
693	& + ( 1._wp - zfacti ) * zti_raw
694	ztj_lim = zfactj * ztj_lim &
695	& + ( 1._wp - zfactj ) * ztj_raw
696	ENDIF
697	triadi(ji+ip,jj ,jk,1-ip,kp) = zti_lim
698	triadj(ji ,jj+jp,jk,1-jp,kp) = ztj_lim
699	!
700	zbu = e1u(ji ,jj) * e2u(ji ,jj) * fse3u(ji ,jj,jk )
701	zbv = e1v(ji ,jj) * e2v(ji ,jj) * fse3v(ji ,jj,jk )
702	zbti = e1t(ji+ip,jj) * e2t(ji+ip,jj) * fse3w(ji+ip,jj,jk+kp)
703	zbtj = e1t(ji,jj+jp) * e2t(ji,jj+jp) * fse3w(ji,jj+jp,jk+kp)
704	!
705	!!gm this may inhibit vectorization on Vect Computers, and even on scalar computers.... ==> to be checked
706	wslp2 (ji+ip,jj,jk+kp) = wslp2(ji+ip,jj,jk+kp) + 0.25_wp * zbu / zbti * zti_g_lim**2 ! masked
707	wslp2 (ji,jj+jp,jk+kp) = wslp2(ji,jj+jp,jk+kp) + 0.25_wp * zbv / zbtj * ztj_g_lim**2
708	END DO
709	END DO
710	END DO
711	END DO
712	END DO
713	!
714	wslp2(:,:,1) = 0._wp ! force the surface wslp to zero
715
716	CALL lbc_lnk( wslp2, 'W', 1. ) ! lateral boundary confition on wslp2 only ==>>> gm : necessary ? to be checked
717	!
718	CALL wrk_dealloc( jpi,jpj, z1_mlbw )
719	CALL wrk_dealloc( jpi,jpj,jpk,2, zdxrho , zdyrho, zdzrho, klstart = 0 )
720	CALL wrk_dealloc( jpi,jpj, 2,2, zti_mlb, ztj_mlb, kkstart = 0, klstart = 0 )
721	!
722	IF( nn_timing == 1 ) CALL timing_stop('ldf_slp_grif')
723	!
724	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
725	END SUBROUTINE ldf_slp_grif
726
727
728	SUBROUTINE ldf_slp_mxl( prd, pn2, p_gru, p_grv, p_dzr )
729	!!----------------------------------------------------------------------
730	!! * ROUTINE ldf_slp_mxl *
731	!!
732	!! ** Purpose : Compute the slopes of iso-neutral surface just below
733	!! the mixed layer.
734	!!
735	!! ** Method : The slope in the i-direction is computed at u- & w-points
736	!! (uslpml, wslpiml) and the slope in the j-direction is computed
737	!! at v- and w-points (vslpml, wslpjml) with the same bounds as
738	!! in ldf_slp.
739	!!
740	!! ** Action : uslpml, wslpiml : i- & j-slopes of neutral surfaces
741	!! vslpml, wslpjml just below the mixed layer
742	!! omlmask : mixed layer mask
743	!!----------------------------------------------------------------------
744	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: prd ! in situ density
745	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: pn2 ! Brunt-Vaisala frequency (locally ref.)
746	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_gru, p_grv ! i- & j-gradient of density (u- & v-pts)
747	REAL(wp), DIMENSION(:,:,:), INTENT(in) :: p_dzr ! z-gradient of density (T-point)
748	!!
749	INTEGER :: ji , jj , jk ! dummy loop indices
750	INTEGER :: iku, ikv, ik, ikm1 ! local integers
751	REAL(wp) :: zeps, zm1_g, zm1_2g ! local scalars
752	REAL(wp) :: zci, zfi, zau, zbu, zai, zbi ! - -
753	REAL(wp) :: zcj, zfj, zav, zbv, zaj, zbj ! - -
754	REAL(wp) :: zck, zfk, zbw ! - -
755	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
756	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
757	REAL(KIND=jprb) :: zhook_handle
758
759	CHARACTER(LEN=*), PARAMETER :: RoutineName='LDF_SLP_MXL'
760
761	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
762
763	!!----------------------------------------------------------------------
764	!
765	IF( nn_timing == 1 ) CALL timing_start('ldf_slp_mxl')
766	!
767	zeps = 1.e-20_wp !== Local constant initialization ==!
768	zm1_g = -1.0_wp / grav
769	zm1_2g = -0.5_wp / grav
770	!
771	uslpml (1,:) = 0._wp ; uslpml (jpi,:) = 0._wp
772	vslpml (1,:) = 0._wp ; vslpml (jpi,:) = 0._wp
773	wslpiml(1,:) = 0._wp ; wslpiml(jpi,:) = 0._wp
774	wslpjml(1,:) = 0._wp ; wslpjml(jpi,:) = 0._wp
775	!
776	! !== surface mixed layer mask !
777	DO jk = 1, jpk ! =1 inside the mixed layer, =0 otherwise
778	DO jj = 1, jpj
779	DO ji = 1, jpi
780	ik = nmln(ji,jj) - 1
781	IF( jk <= ik .AND. jk >= mikt(ji,jj) ) THEN
782	omlmask(ji,jj,jk) = 1._wp
783	ELSE
784	omlmask(ji,jj,jk) = 0._wp
785	ENDIF
786	END DO
787	END DO
788	END DO
789
790
791	! Slopes of isopycnal surfaces just before bottom of mixed layer
792	! --------------------------------------------------------------
793	! The slope are computed as in the 3D case.
794	! A key point here is the definition of the mixed layer at u- and v-points.
795	! It is assumed to be the maximum of the two neighbouring T-point mixed layer depth.
796	! Otherwise, a n2 value inside the mixed layer can be involved in the computation
797	! of the slope, resulting in a too steep diagnosed slope and thus a spurious eddy
798	! induce velocity field near the base of the mixed layer.
799	!-----------------------------------------------------------------------
800	!
801	DO jj = 2, jpjm1
802	DO ji = 2, jpim1
803	! !== Slope at u- & v-points just below the Mixed Layer ==!
804	!
805	! !- vertical density gradient for u- and v-slopes (from dzr at T-point)
806	iku = MIN( MAX( miku(ji,jj)+1, nmln(ji,jj) , nmln(ji+1,jj) ) , jpkm1 ) ! ML (MAX of T-pts, bound by jpkm1)
807	ikv = MIN( MAX( mikv(ji,jj)+1, nmln(ji,jj) , nmln(ji,jj+1) ) , jpkm1 ) !
808	zbu = 0.5_wp * ( p_dzr(ji,jj,iku) + p_dzr(ji+1,jj ,iku) )
809	zbv = 0.5_wp * ( p_dzr(ji,jj,ikv) + p_dzr(ji ,jj+1,ikv) )
810	! !- horizontal density gradient at u- & v-points
811	zau = p_gru(ji,jj,iku) / e1u(ji,jj)
812	zav = p_grv(ji,jj,ikv) / e2v(ji,jj)
813	! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0
814	! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
815	zbu = MIN( zbu , -100._wp* ABS( zau ) , -7.e+3_wp/fse3u(ji,jj,iku)* ABS( zau ) )
816	zbv = MIN( zbv , -100._wp* ABS( zav ) , -7.e+3_wp/fse3v(ji,jj,ikv)* ABS( zav ) )
817	! !- Slope at u- & v-points (uslpml, vslpml)
818	uslpml(ji,jj) = zau / ( zbu - zeps ) * umask(ji,jj,iku)
819	vslpml(ji,jj) = zav / ( zbv - zeps ) * vmask(ji,jj,ikv)
820	!
821	! !== i- & j-slopes at w-points just below the Mixed Layer ==!
822	!
823	ik = MIN( nmln(ji,jj) + 1, jpk )
824	ikm1 = MAX( 1, ik-1 )
825	! !- vertical density gradient for w-slope (from N^2)
826	zbw = zm1_2g * pn2 (ji,jj,ik) * ( prd (ji,jj,ik) + prd (ji,jj,ikm1) + 2. )
827	! !- horizontal density i- & j-gradient at w-points
828	zci = MAX( umask(ji-1,jj,ik ) + umask(ji,jj,ik ) &
829	& + umask(ji-1,jj,ikm1) + umask(ji,jj,ikm1) , zeps ) * e1t(ji,jj)
830	zcj = MAX( vmask(ji,jj-1,ik ) + vmask(ji,jj,ik ) &
831	& + vmask(ji,jj-1,ikm1) + vmask(ji,jj,ikm1) , zeps ) * e2t(ji,jj)
832	zai = ( p_gru(ji-1,jj,ik ) + p_gru(ji,jj,ik) &
833	& + p_gru(ji-1,jj,ikm1) + p_gru(ji,jj,ikm1 ) ) / zci * tmask(ji,jj,ik)
834	zaj = ( p_grv(ji,jj-1,ik ) + p_grv(ji,jj,ik ) &
835	& + p_grv(ji,jj-1,ikm1) + p_grv(ji,jj,ikm1) ) / zcj * tmask(ji,jj,ik)
836	! !- bound the slopes: abs(zw.)<= 1/100 and zb..<0.
837	! kxz max= ah slope max =< e1 e3 /(pi**2 2 dt)
838	zbi = MIN( zbw , -100._wp* ABS( zai ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zai ) )
839	zbj = MIN( zbw , -100._wp* ABS( zaj ) , -7.e+3_wp/fse3w(ji,jj,ik)* ABS( zaj ) )
840	! !- i- & j-slope at w-points (wslpiml, wslpjml)
841	wslpiml(ji,jj) = zai / ( zbi - zeps ) * wmask (ji,jj,ik)
842	wslpjml(ji,jj) = zaj / ( zbj - zeps ) * wmask (ji,jj,ik)
843	END DO
844	END DO
845	!!gm this lbc_lnk should be useless....
846	CALL lbc_lnk( uslpml , 'U', -1. ) ; CALL lbc_lnk( vslpml , 'V', -1. ) ! lateral boundary cond. (sign change)
847	CALL lbc_lnk( wslpiml, 'W', -1. ) ; CALL lbc_lnk( wslpjml, 'W', -1. ) ! lateral boundary conditions
848	!
849	IF( nn_timing == 1 ) CALL timing_stop('ldf_slp_mxl')
850	!
851	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
852	END SUBROUTINE ldf_slp_mxl
853
854
855	SUBROUTINE ldf_slp_init
856	!!----------------------------------------------------------------------
857	!! * ROUTINE ldf_slp_init *
858	!!
859	!! ** Purpose : Initialization for the isopycnal slopes computation
860	!!
861	!! ** Method : read the nammbf namelist and check the parameter
862	!! values called by tra_dmp at the first timestep (nit000)
863	!!----------------------------------------------------------------------
864	INTEGER :: ji, jj, jk ! dummy loop indices
865	INTEGER :: ierr ! local integer
866	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
867	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
868	REAL(KIND=jprb) :: zhook_handle
869
870	CHARACTER(LEN=*), PARAMETER :: RoutineName='LDF_SLP_INIT'
871
872	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
873
874	!!----------------------------------------------------------------------
875	!
876	IF( nn_timing == 1 ) CALL timing_start('ldf_slp_init')
877	!
878	IF(lwp) THEN
879	WRITE(numout,*)
880	WRITE(numout,*) 'ldf_slp_init : direction of lateral mixing'
881	WRITE(numout,*) '~~~~~~~~~~~~'
882	ENDIF
883
884	IF( ln_traldf_grif ) THEN ! Griffies operator : triad of slopes
885	ALLOCATE( triadi_g(jpi,jpj,jpk,0:1,0:1) , triadj_g(jpi,jpj,jpk,0:1,0:1) , wslp2(jpi,jpj,jpk) , STAT=ierr )
886	ALLOCATE( triadi (jpi,jpj,jpk,0:1,0:1) , triadj (jpi,jpj,jpk,0:1,0:1) , STAT=ierr )
887	IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Griffies operator slope' )
888	!
889	IF( ln_dynldf_iso ) CALL ctl_stop( 'ldf_slp_init: Griffies operator on momentum not supported' )
890	!
891	ELSE ! Madec operator : slopes at u-, v-, and w-points
892	ALLOCATE( uslp(jpi,jpj,jpk) , vslp(jpi,jpj,jpk) , wslpi(jpi,jpj,jpk) , wslpj(jpi,jpj,jpk) , &
893	& omlmask(jpi,jpj,jpk) , uslpml(jpi,jpj) , vslpml(jpi,jpj) , wslpiml(jpi,jpj) , wslpjml(jpi,jpj) , STAT=ierr )
894	IF( ierr > 0 ) CALL ctl_stop( 'STOP', 'ldf_slp_init : unable to allocate Madec operator slope ' )
895
896	! Direction of lateral diffusion (tracers and/or momentum)
897	! ------------------------------
898	uslp (:,:,:) = 0._wp ; uslpml (:,:) = 0._wp ! set the slope to zero (even in s-coordinates)
899	vslp (:,:,:) = 0._wp ; vslpml (:,:) = 0._wp
900	wslpi(:,:,:) = 0._wp ; wslpiml(:,:) = 0._wp
901	wslpj(:,:,:) = 0._wp ; wslpjml(:,:) = 0._wp
902
903	IF(ln_sco .AND. (ln_traldf_hor .OR. ln_dynldf_hor )) THEN
904	IF(lwp) WRITE(numout,*) ' Horizontal mixing in s-coordinate: slope = slope of s-surfaces'
905
906	! geopotential diffusion in s-coordinates on tracers and/or momentum
907	! The slopes of s-surfaces are computed once (no call to ldfslp in step)
908	! The slopes for momentum diffusion are i- or j- averaged of those on tracers
909
910	! set the slope of diffusion to the slope of s-surfaces
911	! ( c a u t i o n : minus sign as fsdep has positive value )
912	DO jk = 1, jpk
913	DO jj = 2, jpjm1
914	DO ji = fs_2, fs_jpim1 ! vector opt.
915	uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept_b(ji+1,jj,jk) - fsdept_b(ji ,jj ,jk) ) * umask(ji,jj,jk)
916	vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept_b(ji,jj+1,jk) - fsdept_b(ji ,jj ,jk) ) * vmask(ji,jj,jk)
917	wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw_b(ji+1,jj,jk) - fsdepw_b(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
918	wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw_b(ji,jj+1,jk) - fsdepw_b(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
919	END DO
920	END DO
921	END DO
922	CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. ) ! Lateral boundary conditions
923	CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. )
924	ENDIF
925	ENDIF
926	!
927	IF( nn_timing == 1 ) CALL timing_stop('ldf_slp_init')
928	!
929	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
930	END SUBROUTINE ldf_slp_init
931
932	#else
933	!!------------------------------------------------------------------------
934	!! Dummy module : NO Rotation of lateral mixing tensor
935	!!------------------------------------------------------------------------
936	LOGICAL, PUBLIC, PARAMETER :: lk_ldfslp = .FALSE. !: slopes flag
937	CONTAINS
938	SUBROUTINE ldf_slp( kt, prd, pn2 ) ! Dummy routine
939	INTEGER, INTENT(in) :: kt
940	REAL, DIMENSION(:,:,:), INTENT(in) :: prd, pn2
941	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
942	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
943	REAL(KIND=jprb) :: zhook_handle
944
945	CHARACTER(LEN=*), PARAMETER :: RoutineName='LDF_SLP'
946
947	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
948
949	WRITE(,) 'ldf_slp: You should not have seen this print! error?', kt, prd(1,1,1), pn2(1,1,1)
950	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
951	END SUBROUTINE ldf_slp
952	SUBROUTINE ldf_slp_grif( kt ) ! Dummy routine
953	INTEGER, INTENT(in) :: kt
954	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
955	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
956	REAL(KIND=jprb) :: zhook_handle
957
958	CHARACTER(LEN=*), PARAMETER :: RoutineName='LDF_SLP_GRIF'
959
960	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
961
962	WRITE(,) 'ldf_slp_grif: You should not have seen this print! error?', kt
963	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
964	END SUBROUTINE ldf_slp_grif
965	SUBROUTINE ldf_slp_init ! Dummy routine
966	INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
967	INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
968	REAL(KIND=jprb) :: zhook_handle
969
970	CHARACTER(LEN=*), PARAMETER :: RoutineName='LDF_SLP_INIT'
971
972	IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
973
974	IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
975	END SUBROUTINE ldf_slp_init
976	#endif
977
978	!!======================================================================
979	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: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/LDF/ldfslp.F90

Download in other formats: