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.

dynldf_iso.F90 in branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_iso.F90 @ 4460

Last change on this file since 4460 was 4460, checked in by trackstand2, 10 years ago
Addition of dbg SUM outputs and alteration of loop limits in dynldfiso and dynldf_bilapg
Property svn:keywords set to `Id`
File size: 56.5 KB

Rev	Line
[3]	1	MODULE dynldf_iso
	2	!!======================================================================
	3	!! * MODULE dynldf_iso *
	4	!! Ocean dynamics: lateral viscosity trend
	5	!!======================================================================
[2715]	6	!! History : OPA ! 97-07 (G. Madec) Original code
	7	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
	8	!! - ! 2004-08 (C. Talandier) New trends organization
	9	!! 2.0 ! 2005-11 (G. Madec) s-coordinate: horizontal diffusion
	10	!!----------------------------------------------------------------------
[3]	11	#if defined key_ldfslp \|\| defined key_esopa
	12	!!----------------------------------------------------------------------
	13	!! 'key_ldfslp' slopes of the direction of mixing
	14	!!----------------------------------------------------------------------
	15	!! dyn_ldf_iso : update the momentum trend with the horizontal part
	16	!! of the lateral diffusion using isopycnal or horizon-
	17	!! tal s-coordinate laplacian operator.
	18	!!----------------------------------------------------------------------
	19	USE oce ! ocean dynamics and tracers
	20	USE dom_oce ! ocean space and time domain
	21	USE ldfdyn_oce ! ocean dynamics lateral physics
	22	USE ldftra_oce ! ocean tracer lateral physics
	23	USE zdf_oce ! ocean vertical physics
[216]	24	USE trdmod ! ocean dynamics trends
	25	USE trdmod_oce ! ocean variables trends
[3]	26	USE ldfslp ! iso-neutral slopes
[455]	27	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[3]	28	USE in_out_manager ! I/O manager
[2715]	29	USE lib_mpp ! MPP library
[258]	30	USE prtctl ! Print control
[3]	31
	32	IMPLICIT NONE
	33	PRIVATE
	34
[2715]	35	PUBLIC dyn_ldf_iso ! called by step.F90
	36	PUBLIC dyn_ldf_iso_alloc ! called by nemogcm.F90
[3]	37
[3432]	38	!FTRANS zdiu zdju zdiv zdjv zdj1u zdj1v zfuw zfvw :I :z
[2715]	39	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zdiu, zdju, zdj1u ! 2D workspace (dyn_ldf_iso)
	40	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfvw, zdiv, zdjv, zdj1v ! - -
	41
[3211]	42	!! * Control permutation of array indices
	43	# include "oce_ftrans.h90"
	44	# include "dom_oce_ftrans.h90"
	45	# include "ldfdyn_oce_ftrans.h90"
	46	# include "ldftra_oce_ftrans.h90"
	47	# include "ldfslp_ftrans.h90"
	48	# include "zdf_oce_ftrans.h90"
	49
[3]	50	!! * Substitutions
	51	# include "domzgr_substitute.h90"
	52	# include "ldfdyn_substitute.h90"
	53	# include "vectopt_loop_substitute.h90"
	54	!!----------------------------------------------------------------------
[2715]	55	!! NEMO/OPA 3.3 , NEMO Consortium (2011)
[1156]	56	!! $Id$
[2715]	57	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	58	!!----------------------------------------------------------------------
	59	CONTAINS
	60
[2715]	61	INTEGER FUNCTION dyn_ldf_iso_alloc()
	62	!!----------------------------------------------------------------------
	63	!! * ROUTINE dyn_ldf_iso_alloc *
	64	!!----------------------------------------------------------------------
	65	ALLOCATE( zfuw(jpi,jpk) , zdiu(jpi,jpk) , zdju(jpi,jpk) , zdj1u(jpi,jpk) , &
	66	& zfvw(jpi,jpk) , zdiv(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_iso_alloc )
	67	!
	68	IF( dyn_ldf_iso_alloc /= 0 ) CALL ctl_warn('dyn_ldf_iso_alloc: array allocate failed.')
	69	END FUNCTION dyn_ldf_iso_alloc
	70
	71
[3]	72	SUBROUTINE dyn_ldf_iso( kt )
	73	!!----------------------------------------------------------------------
	74	!! * ROUTINE dyn_ldf_iso *
	75	!!
[455]	76	!! ** Purpose : Compute the before trend of the rotated laplacian
	77	!! operator of lateral momentum diffusion except the diagonal
	78	!! vertical term that will be computed in dynzdf module. Add it
	79	!! to the general trend of momentum equation.
[3]	80	!!
	81	!! ** Method :
[455]	82	!! The momentum lateral diffusive trend is provided by a 2nd
	83	!! order operator rotated along neutral or geopotential surfaces
	84	!! (in s-coordinates).
[3]	85	!! It is computed using before fields (forward in time) and isopyc-
[455]	86	!! nal or geopotential slopes computed in routine ldfslp.
[3]	87	!! Here, u and v components are considered as 2 independent scalar
	88	!! fields. Therefore, the property of splitting divergent and rota-
	89	!! tional part of the flow of the standard, z-coordinate laplacian
	90	!! momentum diffusion is lost.
	91	!! horizontal fluxes associated with the rotated lateral mixing:
	92	!! u-component:
	93	!! ziut = ( ahmt + ahmb0 ) e2t * e3t / e1t di[ ub ]
	94	!! - ahmt e2t * mi-1(uslp) dk[ mi(mk(ub)) ]
	95	!! zjuf = ( ahmf + ahmb0 ) e1f * e3f / e2f dj[ ub ]
	96	!! - ahmf e1f * mi(vslp) dk[ mj(mk(ub)) ]
	97	!! v-component:
	98	!! zivf = ( ahmf + ahmb0 ) e2t * e3t / e1t di[ vb ]
	99	!! - ahmf e2t * mj(uslp) dk[ mi(mk(vb)) ]
	100	!! zjvt = ( ahmt + ahmb0 ) e1f * e3f / e2f dj[ ub ]
	101	!! - ahmt e1f * mj-1(vslp) dk[ mj(mk(vb)) ]
	102	!! take the horizontal divergence of the fluxes:
	103	!! diffu = 1/(e1ue2ue3u) { di [ ziut ] + dj-1[ zjuf ] }
	104	!! diffv = 1/(e1ve2ve3v) { di-1[ zivf ] + dj [ zjvt ] }
	105	!! Add this trend to the general trend (ua,va):
	106	!! ua = ua + diffu
[455]	107	!! CAUTION: here the isopycnal part is with a coeff. of aht. This
	108	!! should be modified for applications others than orca_r2 (!!bug)
[3]	109	!!
	110	!! ** Action :
	111	!! Update (ua,va) arrays with the before geopotential biharmonic
	112	!! mixing trend.
[455]	113	!! Update (avmu,avmv) to accompt for the diagonal vertical component
	114	!! of the rotated operator in dynzdf module
[3]	115	!!----------------------------------------------------------------------
[2715]	116	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
[3432]	117	#if ! defined key_z_first
	118	! Then these workspace arrays can be left as 2D
	119	USE wrk_nemo, ONLY: zjvt => wrk_2d_3 ! 2D workspace
	120	USE wrk_nemo, ONLY: zivf => wrk_2d_4 ! 2D workspace
	121	USE wrk_nemo, ONLY: ziut => wrk_2d_1 , zjuf => wrk_2d_2
	122	USE wrk_nemo, ONLY: zdku => wrk_2d_5 , zdkv => wrk_2d_6
[2715]	123	USE wrk_nemo, ONLY: zdk1u => wrk_2d_7 , zdk1v => wrk_2d_8
[3432]	124	#endif
	125	USE timing, ONLY: timing_start, timing_stop
[2715]	126	!
	127	INTEGER, INTENT( in ) :: kt ! ocean time-step index
	128	!
	129	INTEGER :: ji, jj, jk ! dummy loop indices
	130	REAL(wp) :: zabe1, zabe2, zcof1, zcof2 ! local scalars
	131	REAL(wp) :: zmskt, zmskf, zbu, zbv, zuah, zvah ! - -
	132	REAL(wp) :: zcoef0, zcoef3, zcoef4, zmkt, zmkf ! - -
	133	REAL(wp) :: zuav, zvav, zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - -
[3432]	134	REAL(wp) :: zcof, zrecip
	135	#if defined key_z_first
	136	REAL(wp) :: zdku, zdk1u, zdki1u, zdk1i1u, zdkim1u, zdk1im1u
	137	REAL(wp) :: zdkj1u, zdk1j1u, zdkjm1u, zdk1jm1u
	138	REAL(wp) :: zdkv, zdk1v, zdki1v, zdk1i1v, zdkim1v, zdk1im1v
	139	REAL(wp) :: zdkj1v, zdk1j1v, zdkjm1v, zdk1jm1v
	140	REAL(wp) :: aziut, aziuti1, azjuf, azjufjm1
	141	REAL(wp) :: azivf, azivfim1, azjvtj1, azjvt
	142	#endif
[2715]	143	!!----------------------------------------------------------------------
[3]	144
[3432]	145	! ziut zjvt zjuf zivf :I :I :z
	146	!!$#if defined key_z_first
	147	!!$ REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ziut, zjuf, zivf, zjvt
	148	!!$ ALLOCATE( ziut(jpi,jpj,jpk), zjuf(jpi,jpj,jpk), &
	149	!!$ zivf(jpi,jpj,jpk), zjvt(jpi,jpj,jpk) )
	150	!!$#endif
	151	CALL timing_start('dyn_ldf_iso')
	152
[2715]	153	IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN
	154	CALL ctl_stop('dyn_ldf_iso: requested workspace arrays unavailable') ; RETURN
	155	END IF
[455]	156
[3]	157	IF( kt == nit000 ) THEN
	158	IF(lwp) WRITE(numout,*)
	159	IF(lwp) WRITE(numout,*) 'dyn_ldf_iso : iso-neutral laplacian diffusive operator or '
	160	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ s-coordinate horizontal diffusive operator'
[2715]	161	! ! allocate dyn_ldf_bilap arrays
	162	IF( dyn_ldf_iso_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_iso: failed to allocate arrays')
[3]	163	ENDIF
[216]	164
[2715]	165	! s-coordinate: Iso-level diffusion on momentum but not on tracer
[455]	166	IF( ln_dynldf_hor .AND. ln_traldf_iso ) THEN
[2715]	167	!
[3211]	168	#if defined key_z_first
	169	DO jj = 2, jpjm1 ! set the slopes of iso-level
	170	DO ji = fs_2, fs_jpim1
[4447]	171	DO jk = 1, mbkmax(ji,jj) ! jpk
[3211]	172	#else
[2715]	173	DO jk = 1, jpk ! set the slopes of iso-level
[455]	174	DO jj = 2, jpjm1
	175	DO ji = fs_2, fs_jpim1 ! vector opt.
[3211]	176	#endif
[455]	177	uslp (ji,jj,jk) = -1./e1u(ji,jj) * ( fsdept(ji+1,jj,jk) - fsdept(ji ,jj ,jk) ) * umask(ji,jj,jk)
	178	vslp (ji,jj,jk) = -1./e2v(ji,jj) * ( fsdept(ji,jj+1,jk) - fsdept(ji ,jj ,jk) ) * vmask(ji,jj,jk)
	179	wslpi(ji,jj,jk) = -1./e1t(ji,jj) * ( fsdepw(ji+1,jj,jk) - fsdepw(ji-1,jj,jk) ) * tmask(ji,jj,jk) * 0.5
	180	wslpj(ji,jj,jk) = -1./e2t(ji,jj) * ( fsdepw(ji,jj+1,jk) - fsdepw(ji,jj-1,jk) ) * tmask(ji,jj,jk) * 0.5
	181	END DO
	182	END DO
	183	END DO
	184	! Lateral boundary conditions on the slopes
	185	CALL lbc_lnk( uslp , 'U', -1. ) ; CALL lbc_lnk( vslp , 'V', -1. )
	186	CALL lbc_lnk( wslpi, 'W', -1. ) ; CALL lbc_lnk( wslpj, 'W', -1. )
	187
	188	!!bug
[2715]	189	IF( kt == nit000 ) then
	190	IF(lwp) WRITE(numout,) ' max slop: u', SQRT( MAXVAL(uslpuslp)), ' v ', SQRT(MAXVAL(vslp)), &
	191	& ' wi', sqrt(MAXVAL(wslpi)) , ' wj', sqrt(MAXVAL(wslpj))
[455]	192	endif
	193	!!end
[216]	194	ENDIF
[455]	195
[4453]	196	!CALL timing_start('dyn_ldf_iso_hslab')
[3432]	197
	198	#if defined key_z_first
	199
	200	! Vertical u- and v-shears at level jk and jk+1
	201	! ---------------------------------------------
	202	! surface boundary condition: zdku(jk=1)=zdku(jk=2)
	203	! zdkv(jk=1)=zdkv(jk=2)
	204	!!$ DO jj = 1, jpj, 1
	205	!!$ DO ji = 1, jpi, 1
	206	!!$
	207	!!$ ! jk=1 special case
	208	!!$ !zdk1u(ji,jj,1) = ( ub(ji,jj,1) -ub(ji,jj,2) ) * umask(ji,jj,2)
	209	!!$ zdk1v(ji,jj,1) = ( vb(ji,jj,1) -vb(ji,jj,2) ) * vmask(ji,jj,2)
	210	!!$ !zdku(ji,jj,1) = zdk1u(ji,jj,1)
	211	!!$ zdkv(ji,jj,1) = zdk1v(ji,jj,1)
	212	!!$
	213	!!$ DO jk = 2, jpkm1, 1
	214	!!$ !zdk1u(ji,jj,jk) = ( ub(ji,jj,jk) -ub(ji,jj,jk+1) ) * umask(ji,jj,jk+1)
	215	!!$ zdk1v(ji,jj,jk) = ( vb(ji,jj,jk) -vb(ji,jj,jk+1) ) * vmask(ji,jj,jk+1)
	216	!!$
	217	!!$ !zdku(ji,jj,jk) = ( ub(ji,jj,jk-1) - ub(ji,jj,jk) ) * umask(ji,jj,jk)
	218	!!$ zdkv(ji,jj,jk) = ( vb(ji,jj,jk-1) - vb(ji,jj,jk) ) * vmask(ji,jj,jk)
	219	!!$ END DO !jk
	220	!!$ END DO
	221	!!$ END DO
	222
	223	!!$ ! -----f-----
	224	!!$ ! Horizontal fluxes on U \|
	225	!!$ ! --------------------=== t u t
	226	!!$ ! \|
	227	!!$ ! i-flux at t-point -----f-----
	228	!!$
	229	!!$ IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
	230	!!$ DO jj = 2, jpjm1
	231	!!$ DO ji = 2, jpi
	232	!!$ DO jk = 1, jpkm1, 1
	233	!!$ zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj)
	234	!!$
	235	!!$ zmskt = 1._wp/MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
	236	!!$ & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1._wp )
	237	!!$
	238	!!$ zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5_wp * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
	239	!!$
	240	!!$ ziut(ji,jj,jk) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) &
	241	!!$ & + zcof1 * ( zdku + zdk1im1u &
	242	!!$ & +zdk1u + zdkim1u ) ) * tmask(ji,jj,jk)
	243	!!$ END DO
	244	!!$ END DO
	245	!!$ END DO
	246	!!$ ELSE ! other coordinate system (zco or sco) : e3t
	247	!!$ DO jj = 2, jpjm1
	248	!!$ DO ji = 2, jpi
	249	!!$ DO jk = 1, jpkm1, 1
	250	!!$ zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
	251	!!$
	252	!!$ zmskt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
	253	!!$ & + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )
	254	!!$
	255	!!$ zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
	256	!!$
	257	!!$ ziut(ji,jj,jk) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) &
	258	!!$ & + zcof1 * ( zdku(ji,jj,jk) + zdk1u(ji-1,jj,jk) &
	259	!!$ & +zdk1u(ji,jj,jk) + zdku (ji-1,jj,jk) ) ) * tmask(ji,jj,jk)
	260	!!$
	261	!!$ END DO
	262	!!$ END DO
	263	!!$ END DO
	264	!!$ ENDIF
	265
	266	! j-flux at f-point
	267	! BLOCKABLE(ji,jj,jk)
	268	! BLOCKING SIZE (0)
	269	!!$ DO jj = 1, jpjm1, 1
	270	!!$! BLOCKING SIZE (0)
	271	!!$ DO ji = 1, jpim1, 1
	272	!!$! BLOCKING SIZE (4)
	273	!!$ DO jk = 1, jpkm1, 1
	274	!!$ zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
	275	!!$
	276	!!$ zmskf = 1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) &
	277	!!$ & + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. )
	278	!!$
	279	!!$ zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
	280	!!$
	281	!!$ zjuf(ji,jj,jk) = ( zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) &
	282	!!$ & + zcof2 * ( zdku (ji,jj+1,jk) + zdk1u(ji,jj,jk) &
	283	!!$ & +zdk1u(ji,jj+1,jk) + zdku (ji,jj,jk) ) ) * fmask(ji,jj,jk)
	284	!!$ END DO
	285	!!$ END DO
	286	!!$ END DO
	287
	288	! \| t \|
	289	! Horizontal fluxes on V \| \|
	290	! --------------------=== f---v---f
	291	! \| \|
	292	! \| t \|
	293
	294	!!$ IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
	295	!!$ DO jj = 2, jpj
	296	!!$ DO ji = 1, jpim1
	297	!!$ DO jk = 1, jpkm1, 1
	298	!!$
	299	!!$ ! i-flux at f-point \| t \|
	300	!!$ zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
	301	!!$
	302	!!$ zmskf = 1._wp/MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
	303	!!$ + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1._wp )
	304	!!$
	305	!!$ zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
	306	!!$
	307	!!$ zivf(ji,jj,jk) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) &
	308	!!$ + zcof1 * ( zdkv(ji,jj,jk) + zdk1v(ji+1,jj,jk) &
	309	!!$ +zdk1v(ji,jj,jk) + zdkv(ji+1,jj,jk) ) ) * fmask(ji,jj,jk)
	310	!!$
	311	!!$
	312	!!$ ! j-flux at t-point
	313	!!$ zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj)
	314	!!$
	315	!!$ zmskt = 1._wp/MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	316	!!$ + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1._wp )
	317	!!$
	318	!!$ zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5_wp * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
	319	!!$
	320	!!$ zjvt(ji,jj,jk) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) &
	321	!!$ + zcof2 * ( zdkv(ji,jj-1,jk) + zdk1v(ji,jj,jk) &
	322	!!$ +zdk1v(ji,jj-1,jk) + zdkv (ji,jj,jk) ) ) * tmask(ji,jj,jk)
	323	!!$ END DO
	324	!!$ END DO
	325	!!$ END DO
	326	!!$ ELSE ! other coordinate system (zco or sco) : e3t
	327	!!$ DO jj = 2, jpj
	328	!!$ DO ji = 1, jpim1
	329	!!$ DO jk = 1, jpkm1
	330	!!$
	331	!!$ ! i-flux at f-point
	332	!!$ zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
	333	!!$
	334	!!$ zmskf = 1._wp/MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
	335	!!$ + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1._wp )
	336	!!$
	337	!!$ zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
	338	!!$
	339	!!$ zivf(ji,jj,jk) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) &
	340	!!$ + zcof1 * ( zdkv(ji,jj,jk) + zdk1v(ji+1,jj,jk) &
	341	!!$ +zdk1v(ji,jj,jk) + zdkv(ji+1,jj,jk) ) ) * fmask(ji,jj,jk)
	342	!!$ ! j-flux at t-point
	343	!!$ zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
	344	!!$
	345	!!$ zmskt = 1._wp/MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	346	!!$ + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1._wp )
	347	!!$
	348	!!$ zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5_wp * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
	349	!!$
	350	!!$ zjvt(ji,jj,jk) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) &
	351	!!$ + zcof2 * ( zdkv (ji,jj-1,jk) + zdk1v(ji,jj,jk) &
	352	!!$ +zdk1v(ji,jj-1,jk) + zdkv (ji,jj,jk) ) ) * tmask(ji,jj,jk)
	353	!!$ END DO
	354	!!$ END DO
	355	!!$ END DO
	356	!!$
	357	!!$ ENDIF
	358
	359
	360
	361	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
	362
	363	DO jj = 2, jpjm1
	364	DO ji = 2, jpim1
	365	!DIR$ SHORTLOOP
	366	! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the
	367	! If condition on jk>1. This does mean that there will be out-of-bounds reads
	368	! when jk is 1 but that doesn't matter.
	369	!DIR$ SAFE_ADDRESS
[4447]	370	DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1, 1
[3432]	371
	372	! Vertical u- and v-shears at level jk and jk+1
	373	! ---------------------------------------------
	374	! surface boundary condition: zdku(jk=1)=zdku(jk=2)
	375	! zdkv(jk=1)=zdkv(jk=2)
	376	zdk1u = ( ub(ji ,jj ,jk) -ub(ji ,jj ,jk+1) ) * umask(ji ,jj,jk+1)
	377	zdk1i1u = ( ub(ji+1,jj ,jk) -ub(ji+1,jj ,jk+1) ) * umask(ji+1,jj,jk+1)
	378	zdk1j1u = ( ub(ji ,jj+1,jk) -ub(ji ,jj+1,jk+1) ) * umask(ji ,jj+1,jk+1)
	379	zdk1im1u = ( ub(ji-1,jj ,jk) -ub(ji-1,jj ,jk+1) ) * umask(ji-1,jj,jk+1)
	380	zdk1jm1u = ( ub(ji ,jj-1,jk) -ub(ji ,jj-1,jk+1) ) * umask(ji ,jj-1,jk+1)
	381	IF(jk > 1)THEN
	382	zdku = ( ub(ji ,jj ,jk-1) - ub(ji ,jj ,jk) ) * umask(ji,jj,jk)
	383	zdki1u = ( ub(ji+1,jj ,jk-1) - ub(ji+1,jj ,jk) ) * umask(ji+1,jj,jk)
	384	zdkim1u= ( ub(ji-1,jj ,jk-1) - ub(ji-1,jj ,jk) ) * umask(ji-1,jj,jk)
	385	zdkj1u = ( ub(ji ,jj+1,jk-1) - ub(ji ,jj+1,jk) ) * umask(ji,jj+1,jk)
	386	zdkjm1u= ( ub(ji ,jj-1,jk-1) - ub(ji ,jj-1,jk) ) * umask(ji,jj-1,jk)
	387	ELSE
	388	zdku = zdk1u
	389	zdki1u = zdk1i1u
	390	zdkim1u= zdk1im1u
	391	zdkj1u = zdk1j1u
	392	zdkjm1u= zdk1jm1u
	393	END IF
	394
	395	zdku = zdku + zdk1u
	396	zdki1u = zdki1u + zdk1i1u
	397	zdkim1u = zdkim1u + zdk1im1u
	398	zdkj1u = zdkj1u +zdk1j1u
	399	zdkjm1u = zdk1jm1u + zdkjm1u
	400
	401	! volume elements
	402	zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
	403
	404	! horizontal component of isopycnal momentum diffusive trends
	405
	406	! -----f-----
	407	! Horizontal fluxes on U \|
	408	! --------------------=== t u t
	409	! \|
	410	! i-flux at t-point -----f-----
	411
	412	! z-coordinate - partial steps : min(e3u)
	413	aziuti1 = ( ((fsahmt(ji+1,jj,jk)+ahmb0) * e2t(ji+1,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji+1,jj)) * &
	414	( ub(ji+1,jj,jk) - ub(ji,jj,jk) ) + (- aht0 * e2t(ji+1,jj) * &
	415	(1._wp/MAX( umask(ji,jj,jk )+umask(ji+1,jj,jk+1) &
	416	+ umask(ji,jj,jk+1)+umask(ji+1,jj,jk ), 1._wp )) * 0.5 * &
	417	( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) )) * &
	418	( zdku + zdki1u ) ) * tmask(ji+1,jj,jk)
	419
	420	aziut = ( ((fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * &
	421	( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) + (- aht0 * e2t(ji,jj) * &
	422	(1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
	423	+ umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )) * 0.5 * &
	424	( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )) * &
	425	( zdku + zdkim1u ) ) * tmask(ji,jj,jk)
	426
	427	! j-flux at f-point - identical to non-ln_zps
	428
	429	azjuf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)) * &
	430	( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) + (- aht0 * e1f(ji,jj) * &
	431	(1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) &
	432	+ umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. )) * 0.5 * &
	433	( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )) * &
	434	( zdku + zdkj1u ) ) * fmask(ji,jj,jk)
	435
	436	azjufjm1 = ( (( fsahmf(ji,jj-1,jk) + ahmb0 ) * e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1)) * &
	437	( ub(ji,jj,jk) - ub(ji,jj-1,jk) ) + (- aht0 * e1f(ji,jj-1) * &
	438	(1./MAX( umask(ji,jj,jk )+umask(ji,jj-1,jk+1) &
	439	+ umask(ji,jj,jk+1)+umask(ji,jj-1,jk ), 1. )) * 0.5 * &
	440	( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) )) * &
	441	( zdku + zdkjm1u ) ) * fmask(ji,jj-1,jk)
	442
	443	! Second derivative (divergence) and add to the general trend
	444	! -----------------------------------------------------------
	445
	446	zuah =( &
	447	! ziut (ji+1,jj,jk) - &
	448	aziuti1 - &
	449	! ziut (ji ,jj,jk) + &
	450	aziut + &
	451	! zjuf (ji ,jj,jk) - &
	452	azjuf - &
	453	! zjuf (ji,jj-1,jk) &
	454	azjufjm1 &
	455	) / zbu
	456
	457	! add the trends to the general trends
	458	ua (ji,jj,jk) = ua (ji,jj,jk) + zuah
	459
	460	END DO
	461
	462	!DIR$ SHORTLOOP
	463	! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the
	464	! If condition on jk>1. This does mean that there will be out-of-bounds reads
	465	! when jk is 1 but that doesn't matter.
	466	!DIR$ SAFE_ADDRESS
[4447]	467	DO jk = 1, mbkmax(ji,jj)-1, 1 ! jpkm1, 1
[3432]	468
	469	zdk1v = ( vb(ji ,jj ,jk) -vb(ji ,jj ,jk+1) ) * vmask(ji ,jj,jk+1)
	470	zdk1i1v = ( vb(ji+1,jj ,jk) -vb(ji+1,jj ,jk+1) ) * vmask(ji+1,jj,jk+1)
	471	zdk1im1v = ( vb(ji-1,jj ,jk) -vb(ji-1,jj ,jk+1) ) * vmask(ji-1,jj,jk+1)
	472	zdk1j1v = ( vb(ji ,jj+1,jk) -vb(ji ,jj+1,jk+1) ) * vmask(ji ,jj+1,jk+1)
	473	zdk1jm1v = ( vb(ji ,jj-1,jk) -vb(ji ,jj-1,jk+1) ) * vmask(ji ,jj-1,jk+1)
	474	IF(jk > 1)THEN
	475	zdkv = ( vb(ji ,jj ,jk-1) - vb(ji ,jj ,jk) ) * vmask(ji,jj,jk)
	476	zdki1v = ( vb(ji+1,jj ,jk-1) - vb(ji+1,jj ,jk) ) * vmask(ji+1,jj,jk)
	477	zdkim1v= ( vb(ji-1,jj ,jk-1) - vb(ji-1,jj ,jk) ) * vmask(ji-1,jj,jk)
	478	zdkj1v = ( vb(ji ,jj+1,jk-1) - vb(ji ,jj+1,jk) ) * vmask(ji,jj+1,jk)
	479	zdkjm1v= ( vb(ji ,jj-1,jk-1) - vb(ji ,jj-1,jk) ) * vmask(ji,jj-1,jk)
	480	ELSE
	481	zdkv = zdk1v
	482	zdki1v = zdk1i1v
	483	zdkim1v= zdk1im1v
	484	zdkj1v = zdk1j1v
	485	zdkjm1v= zdk1jm1v
	486	END IF
	487
	488	zdkv = zdkv + zdk1v
	489	zdki1v = zdk1i1v + zdki1v
	490	zdkim1v = zdkim1v +zdk1im1v
	491	zdkj1v = zdk1j1v + zdkj1v
	492	zdkjm1v = zdkjm1v +zdk1jm1v
	493
	494	! \| t \|
	495	! Horizontal fluxes on V \| \|
	496	! --------------------=== f---v---f
	497	! \| \|
	498	! \| t \|
	499
	500	! i-flux at f-point - identical to non-ln_zps case
	501	azivf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / &
	502	e1f(ji,jj)) * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) &
	503	+ (- aht0 * e2f(ji,jj) / &
	504	MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
	505	+ vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1._wp ) * &
	506	0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * &
	507	( zdkv + zdki1v ) ) * fmask(ji,jj,jk)
	508
	509	azivfim1 = ( (( fsahmf(ji-1,jj,jk) + ahmb0 ) * e2f(ji-1,jj) * &
	510	fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( vb(ji,jj,jk) - vb(ji-1,jj,jk) ) &
	511	+ (- aht0 * e2f(ji-1,jj) / &
	512	MAX( vmask(ji,jj,jk )+vmask(ji-1,jj,jk+1) &
	513	+ vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk ), 1._wp ) * 0.5_wp * &
	514	( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * &
	515	( zdkim1v + zdkv ) ) * fmask(ji-1,jj,jk)
	516
	517	! j-flux at t-point - min(e3u) instead of e3t
	518	azjvtj1 = ( &
	519	((fsahmt(ji,jj+1,jk)+ahmb0) * e1t(ji,jj+1) * MIN( fse3v(ji,jj+1,jk), fse3v(ji,jj,jk) ) / &
	520	e2t(ji,jj+1)) * ( vb(ji,jj+1,jk) - vb(ji,jj,jk) ) &
	521	+ (- aht0 * e1t(ji,jj+1) / &
	522	MAX( vmask(ji,jj,jk )+vmask(ji,jj+1,jk+1) &
	523	+ vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk ), 1._wp ) * &
	524	0.5_wp * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * &
	525	( zdkv + zdkj1v ) &
	526	) * tmask(ji,jj+1,jk)
	527
	528	azjvt = ( &
	529	((fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / &
	530	e2t(ji,jj)) * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) &
	531	+ (- aht0 * e1t(ji,jj) /MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	532	+ vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1._wp ) * 0.5_wp * &
	533	( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * &
	534	( zdkjm1v + zdkv ) ) * tmask(ji,jj,jk)
	535
	536	! volume elements
	537	zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
	538
	539
	540	! Second derivative (divergence) and add to the general trend
	541	! -----------------------------------------------------------
	542	zvah =( &
	543	! zivf (ji ,jj ,jk) - &
	544	azivf - &
	545	! zivf (ji-1,jj ,jk) + &
	546	azivfim1 + &
	547	! zjvt (ji ,jj+1,jk) - &
	548	azjvtj1 - &
	549	! zjvt (ji ,jj ,jk) &
	550	azjvt &
	551	) / zbv
	552
	553	! add the trend to the general trend
	554	va (ji,jj,jk) = va (ji,jj,jk) + zvah
	555	END DO
	556	END DO
	557	END DO
	558
	559	ELSE ! other coordinate system (zco or sco) : e3t
	560
	561	DO jj = 2, jpjm1
	562	DO ji = 2, jpim1
	563	!DIR$ SHORTLOOP
	564	! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the
	565	! If condition on jk>1. This does mean that there will be out-of-bounds reads
	566	! when jk is 1 but that doesn't matter.
	567	!DIR$ SAFE_ADDRESS
[4447]	568	DO jk = 1, mbkmax(ji,jj)-1, 1 ! jpkm1, 1
[3432]	569
	570	! Vertical u- and v-shears at level jk and jk+1
	571	! ---------------------------------------------
	572	! surface boundary condition: zdku(jk=1)=zdku(jk=2)
	573	! zdkv(jk=1)=zdkv(jk=2)
	574	zdk1u = ( ub(ji ,jj ,jk) -ub(ji ,jj ,jk+1) ) * umask(ji ,jj,jk+1)
	575	zdk1i1u = ( ub(ji+1,jj ,jk) -ub(ji+1,jj ,jk+1) ) * umask(ji+1,jj,jk+1)
	576	zdk1j1u = ( ub(ji ,jj+1,jk) -ub(ji ,jj+1,jk+1) ) * umask(ji ,jj+1,jk+1)
	577	zdk1im1u = ( ub(ji-1,jj ,jk) -ub(ji-1,jj ,jk+1) ) * umask(ji-1,jj,jk+1)
	578	zdk1jm1u = ( ub(ji ,jj-1,jk) -ub(ji ,jj-1,jk+1) ) * umask(ji ,jj-1,jk+1)
	579	IF(jk > 1)THEN
	580	zdku = ( ub(ji ,jj ,jk-1) - ub(ji ,jj ,jk) ) * umask(ji,jj,jk)
	581	zdki1u = ( ub(ji+1,jj ,jk-1) - ub(ji+1,jj ,jk) ) * umask(ji+1,jj,jk)
	582	zdkim1u= ( ub(ji-1,jj ,jk-1) - ub(ji-1,jj ,jk) ) * umask(ji-1,jj,jk)
	583	zdkj1u = ( ub(ji ,jj+1,jk-1) - ub(ji ,jj+1,jk) ) * umask(ji,jj+1,jk)
	584	zdkjm1u= ( ub(ji ,jj-1,jk-1) - ub(ji ,jj-1,jk) ) * umask(ji,jj-1,jk)
	585	ELSE
	586	zdku = zdk1u
	587	zdki1u = zdk1i1u
	588	zdkim1u= zdk1im1u
	589	zdkj1u = zdk1j1u
	590	zdkjm1u= zdk1jm1u
	591	END IF
	592
	593	zdku = zdku + zdk1u
	594	zdki1u = zdki1u + zdk1i1u
	595	zdkim1u = zdkim1u + zdk1im1u
	596	zdkj1u = zdkj1u +zdk1j1u
	597	zdkjm1u = zdk1jm1u + zdkjm1u
	598
	599	! volume element
	600	zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
	601
	602	! horizontal component of isopycnal momentum diffusive trends
	603
	604	! -----f-----
	605	! Horizontal fluxes on U \|
	606	! --------------------=== t u t
	607	! \|
	608	! i-flux at t-point -----f-----
	609
	610	! other coordinate system (zco or sco) : e3t
	611	aziuti1 = ( ((fsahmt(ji+1,jj,jk)+ahmb0) * e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * &
	612	( ub(ji+1,jj,jk) - ub(ji,jj,jk) ) + (- aht0 * e2t(ji+1,jj) * &
	613	(1./MAX( umask(ji,jj,jk )+umask(ji+1,jj,jk+1) &
	614	+ umask(ji,jj,jk+1)+umask(ji+1,jj,jk ), 1. )) * 0.5 * &
	615	( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) )) * &
	616	( zdku + zdki1u ) ) * tmask(ji+1,jj,jk)
	617
	618	aziut = ( ((fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * &
	619	( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) + (- aht0 * e2t(ji,jj) * &
	620	(1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
	621	+ umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )) * 0.5 * &
	622	( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )) * &
	623	( zdku + zdkim1u ) ) * tmask(ji,jj,jk)
	624
	625	azjuf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)) * &
	626	( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) + (- aht0 * e1f(ji,jj) * &
	627	(1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) &
	628	+ umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. )) * 0.5 * &
	629	( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )) * &
	630	( zdku + zdkj1u ) ) * fmask(ji,jj,jk)
	631
	632	azjufjm1 = ( (( fsahmf(ji,jj-1,jk) + ahmb0 ) * e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1)) * &
	633	( ub(ji,jj,jk) - ub(ji,jj-1,jk) ) + (- aht0 * e1f(ji,jj-1) * &
	634	(1./MAX( umask(ji,jj,jk )+umask(ji,jj-1,jk+1) &
	635	+ umask(ji,jj,jk+1)+umask(ji,jj-1,jk ), 1. )) * 0.5 * &
	636	( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) )) * &
	637	( zdku + zdkjm1u ) ) * fmask(ji,jj-1,jk)
	638
	639
	640	! Second derivative (divergence) and add to the general trend
	641	! -----------------------------------------------------------
	642	zuah =( &
	643	! ziut (ji+1,jj,jk) - &
	644	aziuti1 - &
	645	! ziut (ji ,jj,jk) + &
	646	aziut + &
	647	! zjuf (ji ,jj,jk) - &
	648	azjuf - &
	649	! zjuf (ji,jj-1,jk) &
	650	azjufjm1 &
	651	) / zbu
	652
	653	! add the trend to the general trend
	654	ua (ji,jj,jk) = ua (ji,jj,jk) + zuah
	655
	656	END DO
	657
	658	!DIR$ SHORTLOOP
	659	! SAFE_ADDRESS allows the compiler to generate code that will vectorise despite the
	660	! If condition on jk>1. This does mean that there will be out-of-bounds reads
	661	! when jk is 1 but that doesn't matter.
	662	!DIR$ SAFE_ADDRESS
[4447]	663	DO jk = 1, mbkmax(ji,jj)-1, 1 ! jpkm1, 1
[3432]	664
	665	zdk1v = ( vb(ji ,jj ,jk) -vb(ji ,jj ,jk+1) ) * vmask(ji ,jj,jk+1)
	666	zdk1i1v = ( vb(ji+1,jj ,jk) -vb(ji+1,jj ,jk+1) ) * vmask(ji+1,jj,jk+1)
	667	zdk1im1v = ( vb(ji-1,jj ,jk) -vb(ji-1,jj ,jk+1) ) * vmask(ji-1,jj,jk+1)
	668	zdk1j1v = ( vb(ji ,jj+1,jk) -vb(ji ,jj+1,jk+1) ) * vmask(ji ,jj+1,jk+1)
	669	zdk1jm1v = ( vb(ji ,jj-1,jk) -vb(ji ,jj-1,jk+1) ) * vmask(ji ,jj-1,jk+1)
	670	IF(jk > 1)THEN
	671	zdkv = ( vb(ji ,jj ,jk-1) - vb(ji ,jj ,jk) ) * vmask(ji,jj,jk)
	672	zdki1v = ( vb(ji+1,jj ,jk-1) - vb(ji+1,jj ,jk) ) * vmask(ji+1,jj,jk)
	673	zdkim1v= ( vb(ji-1,jj ,jk-1) - vb(ji-1,jj ,jk) ) * vmask(ji-1,jj,jk)
	674	zdkj1v = ( vb(ji ,jj+1,jk-1) - vb(ji ,jj+1,jk) ) * vmask(ji,jj+1,jk)
	675	zdkjm1v= ( vb(ji ,jj-1,jk-1) - vb(ji ,jj-1,jk) ) * vmask(ji,jj-1,jk)
	676	ELSE
	677	zdkv = zdk1v
	678	zdki1v = zdk1i1v
	679	zdkim1v= zdk1im1v
	680	zdkj1v = zdk1j1v
	681	zdkjm1v= zdk1jm1v
	682	END IF
	683
	684	zdkv = zdkv + zdk1v
	685	zdki1v = zdk1i1v + zdki1v
	686	zdkim1v = zdkim1v +zdk1im1v
	687	zdkj1v = zdk1j1v + zdkj1v
	688	zdkjm1v = zdkjm1v +zdk1jm1v
	689
	690	azivf = ( (( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / &
	691	e1f(ji,jj)) * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) &
	692	+ (- aht0 * e2f(ji,jj) / &
	693	MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
	694	+ vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1._wp ) * &
	695	0.5_wp * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * &
	696	( zdkv + zdki1v ) ) * fmask(ji,jj,jk)
	697
	698	azivfim1 = ( (( fsahmf(ji-1,jj,jk) + ahmb0 ) * e2f(ji-1,jj) * &
	699	fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( vb(ji,jj,jk) - vb(ji-1,jj,jk) ) &
	700	+ (- aht0 * e2f(ji-1,jj) / &
	701	MAX( vmask(ji,jj,jk )+vmask(ji-1,jj,jk+1) &
	702	+ vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk ), 1._wp ) * 0.5_wp * &
	703	( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * &
	704	( zdkim1v + zdkv ) ) * fmask(ji-1,jj,jk)
	705
	706	azjvtj1 = ( &
	707	((fsahmt(ji,jj+1,jk)+ahmb0) * e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / &
	708	e2t(ji,jj+1)) * ( vb(ji,jj+1,jk) - vb(ji,jj,jk) ) &
	709	+ (- aht0 * e1t(ji,jj+1) / &
	710	MAX( vmask(ji,jj,jk )+vmask(ji,jj+1,jk+1) &
	711	+ vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk ), 1._wp ) * &
	712	0.5_wp * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * &
	713	( zdkv + zdkj1v ) &
	714	) * tmask(ji,jj+1,jk)
	715
	716	azjvt = ( ((fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * &
	717	( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) &
	718	+ (- aht0 * e1t(ji,jj) /MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	719	+ vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1._wp ) * 0.5_wp * &
	720	( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * &
	721	( zdkjm1v + zdkv ) ) * tmask(ji,jj,jk)
	722
	723	! volume elements
	724	zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
	725
	726	! Second derivative (divergence) and add to the general trend
	727	! -----------------------------------------------------------
	728	zvah =( &
	729	! zivf (ji ,jj ,jk) - &
	730	azivf - &
	731	! zivf (ji-1,jj ,jk) + &
	732	azivfim1 + &
	733	! zjvt (ji ,jj+1,jk) - &
	734	azjvtj1 - &
	735	! zjvt (ji ,jj ,jk) &
	736	azjvt &
	737	) / zbv
	738
	739	! add the trend to the general trend
	740	va (ji,jj,jk) = va (ji,jj,jk) + zvah
	741	END DO
	742	END DO
	743	END DO
	744
	745	END IF ! ln_zps
	746	#else
[3]	747	! ! ===============
	748	DO jk = 1, jpkm1 ! Horizontal slab
	749	! ! ===============
	750
	751	! Vertical u- and v-shears at level jk and jk+1
	752	! ---------------------------------------------
	753	! surface boundary condition: zdku(jk=1)=zdku(jk=2)
	754	! zdkv(jk=1)=zdkv(jk=2)
	755
	756	zdk1u(:,:) = ( ub(:,:,jk) -ub(:,:,jk+1) ) * umask(:,:,jk+1)
	757	zdk1v(:,:) = ( vb(:,:,jk) -vb(:,:,jk+1) ) * vmask(:,:,jk+1)
	758
	759	IF( jk == 1 ) THEN
	760	zdku(:,:) = zdk1u(:,:)
	761	zdkv(:,:) = zdk1v(:,:)
	762	ELSE
	763	zdku(:,:) = ( ub(:,:,jk-1) - ub(:,:,jk) ) * umask(:,:,jk)
	764	zdkv(:,:) = ( vb(:,:,jk-1) - vb(:,:,jk) ) * vmask(:,:,jk)
	765	ENDIF
	766
	767	! -----f-----
	768	! Horizontal fluxes on U \|
	769	! --------------------=== t u t
	770	! \|
	771	! i-flux at t-point -----f-----
	772
[455]	773	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
	774	DO jj = 2, jpjm1
	775	DO ji = fs_2, jpi ! vector opt.
	776	zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * MIN( fse3u(ji,jj,jk), fse3u(ji-1,jj,jk) ) / e1t(ji,jj)
[3]	777
[455]	778	zmskt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
	779	& + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )
[3]	780
[455]	781	zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
	782
	783	ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) &
	784	& + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	785	& +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
	786	END DO
	787	END DO
	788	ELSE ! other coordinate system (zco or sco) : e3t
	789	DO jj = 2, jpjm1
	790	DO ji = fs_2, jpi ! vector opt.
	791	zabe1 = (fsahmt(ji,jj,jk)+ahmb0) * e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
[3]	792
[455]	793	zmskt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
	794	& + umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )
	795
	796	zcof1 = - aht0 * e2t(ji,jj) * zmskt * 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
	797
	798	ziut(ji,jj) = ( zabe1 * ( ub(ji,jj,jk) - ub(ji-1,jj,jk) ) &
	799	& + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	800	& +zdk1u(ji,jj) + zdku (ji-1,jj) ) ) * tmask(ji,jj,jk)
	801	END DO
[3]	802	END DO
[455]	803	ENDIF
[3]	804
	805	! j-flux at f-point
	806	DO jj = 1, jpjm1
	807	DO ji = 1, fs_jpim1 ! vector opt.
[455]	808	zabe2 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
[3]	809
	810	zmskf = 1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) &
[455]	811	& + umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. )
[3]	812
[455]	813	zcof2 = - aht0 * e1f(ji,jj) * zmskf * 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
[3]	814
[455]	815	zjuf(ji,jj) = ( zabe2 * ( ub(ji,jj+1,jk) - ub(ji,jj,jk) ) &
	816	& + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &
	817	& +zdk1u(ji,jj+1) + zdku (ji,jj) ) ) * fmask(ji,jj,jk)
[3]	818	END DO
	819	END DO
	820
	821	! \| t \|
	822	! Horizontal fluxes on V \| \|
	823	! --------------------=== f---v---f
	824	! \| \|
	825	! i-flux at f-point \| t \|
	826
	827	DO jj = 2, jpjm1
	828	DO ji = 1, fs_jpim1 ! vector opt.
[455]	829	zabe1 = ( fsahmf(ji,jj,jk) + ahmb0 ) * e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
[3]	830
	831	zmskf = 1./MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
[455]	832	& + vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. )
[3]	833
[455]	834	zcof1 = - aht0 * e2f(ji,jj) * zmskf * 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
[3]	835
[455]	836	zivf(ji,jj) = ( zabe1 * ( vb(ji+1,jj,jk) - vb(ji,jj,jk) ) &
	837	& + zcof1 * ( zdkv (ji,jj) + zdk1v(ji+1,jj) &
	838	& +zdk1v(ji,jj) + zdkv (ji+1,jj) ) ) * fmask(ji,jj,jk)
[3]	839	END DO
	840	END DO
	841
	842	! j-flux at t-point
[455]	843	IF( ln_zps ) THEN ! z-coordinate - partial steps : min(e3u)
	844	DO jj = 2, jpj
	845	DO ji = 1, fs_jpim1 ! vector opt.
	846	zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * MIN( fse3v(ji,jj,jk), fse3v(ji,jj-1,jk) ) / e2t(ji,jj)
[3]	847
[455]	848	zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	849	& + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. )
[3]	850
[455]	851	zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
[3]	852
[455]	853	zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) &
	854	& + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
	855	& +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
	856	END DO
[3]	857	END DO
[455]	858	ELSE ! other coordinate system (zco or sco) : e3t
	859	DO jj = 2, jpj
	860	DO ji = 1, fs_jpim1 ! vector opt.
	861	zabe2 = (fsahmt(ji,jj,jk)+ahmb0) * e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
[3]	862
[455]	863	zmskt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	864	& + vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. )
[3]	865
[455]	866	zcof2 = - aht0 * e1t(ji,jj) * zmskt * 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
	867
	868	zjvt(ji,jj) = ( zabe2 * ( vb(ji,jj,jk) - vb(ji,jj-1,jk) ) &
	869	& + zcof2 * ( zdkv (ji,jj-1) + zdk1v(ji,jj) &
	870	& +zdk1v(ji,jj-1) + zdkv (ji,jj) ) ) * tmask(ji,jj,jk)
	871	END DO
	872	END DO
	873	ENDIF
	874
	875
[3]	876	! Second derivative (divergence) and add to the general trend
	877	! -----------------------------------------------------------
	878
	879	DO jj = 2, jpjm1
	880	DO ji = 2, jpim1 !! Question vectop possible??? !!bug
	881	! volume elements
	882	zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
	883	zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
	884	! horizontal component of isopycnal momentum diffusive trends
	885	zuah =( ziut (ji+1,jj) - ziut (ji,jj ) + &
[455]	886	& zjuf (ji ,jj) - zjuf (ji,jj-1) ) / zbu
[3]	887	zvah =( zivf (ji,jj ) - zivf (ji-1,jj) + &
[455]	888	& zjvt (ji,jj+1) - zjvt (ji,jj ) ) / zbv
[3]	889	! add the trends to the general trends
	890	ua (ji,jj,jk) = ua (ji,jj,jk) + zuah
	891	va (ji,jj,jk) = va (ji,jj,jk) + zvah
	892	END DO
	893	END DO
	894	! ! ===============
	895	END DO ! End of slab
	896	! ! ===============
[3432]	897	#endif
[4453]	898	!CALL timing_stop('dyn_ldf_iso_hslab','section')
[216]	899
[455]	900	! print sum trends (used for debugging)
	901	IF(ln_ctl) CALL prt_ctl( tab3d_1=ua, clinfo1=' ldfh - Ua: ', mask1=umask, &
	902	& tab3d_2=va, clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
[216]	903
[4453]	904	!CALL timing_start('dyn_ldf_iso_vslab')
[216]	905
[3432]	906	#if defined key_z_first
[455]	907	! ! ===============
	908	DO jj = 2, jpjm1 ! Vertical slab
	909	! ! ===============
	910
[3432]	911	! I. vertical trends associated with the lateral mixing
	912	! =====================================================
	913	! (excluding the vertical flux proportional to dk[t]
	914
	915
	916	! I.1 horizontal momentum gradient
	917	! --------------------------------
	918
	919	DO ji = 2, jpi
[4447]	920	DO jk = 1, mbkmax(ji,jj) ! jpk
[3432]	921	! i-gradient of u at jj
	922	zdiu (ji,jk) = tmask(ji,jj ,jk) * ( ub(ji,jj ,jk) - ub(ji-1,jj ,jk) )
	923	! j-gradient of u and v at jj
	924	zdju (ji,jk) = fmask(ji,jj ,jk) * ( ub(ji,jj+1,jk) - ub(ji ,jj ,jk) )
	925	zdjv (ji,jk) = tmask(ji,jj ,jk) * ( vb(ji,jj ,jk) - vb(ji ,jj-1,jk) )
	926	! j-gradient of u and v at jj+1
	927	zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj ,jk) - ub(ji ,jj-1,jk) )
	928	zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji ,jj ,jk) )
	929	END DO
	930	END DO
	931
	932	DO ji = 1, jpim1
[4447]	933	DO jk = 1, mbkmax(ji,jj) ! jpk
[3432]	934	! i-gradient of v at jj
	935	zdiv (ji,jk) = fmask(ji,jj ,jk) * ( vb(ji+1,jj,jk) - vb(ji ,jj ,jk) )
	936	END DO
	937	END DO
	938
	939
	940	! I.2 Vertical fluxes
	941	! -------------------
	942
	943	! Surface and bottom vertical fluxes set to zero
	944	!!$ DO ji = 1, jpi
	945	!!$ zfuw(ji, 1 ) = 0.e0
	946	!!$ zfvw(ji, 1 ) = 0.e0
	947	!!$ zfuw(ji,jpk) = 0.e0
	948	!!$ zfvw(ji,jpk) = 0.e0
	949	!!$ END DO
	950
	951	! interior (2=<jk=<jpk-1) on U field
	952	DO ji = 2, jpim1
[4447]	953	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
[3432]	954	zcoef0= 0.5 * aht0 * umask(ji,jj,jk)
	955
	956	zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
	957	zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
	958
	959	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) &
	960	+ tmask(ji,jj,jk )+tmask(ji+1,jj,jk ), 1. )
	961	zmkf = 1./MAX( fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1) &
	962	+ fmask(ji,jj-1,jk )+fmask(ji,jj,jk ), 1. )
	963
	964	zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
	965	zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
	966	! vertical flux on u field
	967	zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) &
	968	+zdiu (ji,jk ) + zdiu (ji+1,jk ) ) &
	969	+ zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) &
	970	+zdj1u(ji,jk ) + zdju (ji ,jk ) )
	971	! update avmu (add isopycnal vertical coefficient to avmu)
	972	! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
	973	avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0
	974	END DO
	975	END DO
	976
	977	! interior (2=<jk=<jpk-1) on V field
	978	DO ji = 2, jpim1
[4447]	979	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
[3432]	980	zcoef0= 0.5 * aht0 * vmask(ji,jj,jk)
	981
	982	zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
	983	zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
	984
	985	zmkf = 1./MAX( fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1) &
	986	+ fmask(ji-1,jj,jk )+fmask(ji,jj,jk ), 1. )
	987	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1) &
	988	+ tmask(ji,jj,jk )+tmask(ji,jj+1,jk ), 1. )
	989
	990	zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
	991	zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
	992	! vertical flux on v field
	993	zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) &
	994	+zdiv (ji,jk ) + zdiv (ji-1,jk ) ) &
	995	+ zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) &
	996	+zdjv (ji,jk ) + zdj1v(ji ,jk ) )
	997	! update avmv (add isopycnal vertical coefficient to avmv)
	998	! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
	999	avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0
	1000	END DO
	1001	! END DO
	1002
	1003
	1004	! I.3 Divergence of vertical fluxes added to the general tracer trend
	1005	! -------------------------------------------------------------------
	1006	! DO ji = 2, jpim1
[4447]	1007	zfuw(ji,1) = 0.0_wp
	1008	zfuw(ji,mbkmax(ji,jj)) = 0.0_wp
	1009	zfvw(ji,1) = 0.0_wp
	1010	zfvw(ji,mbkmax(ji,jj)) = 0.0_wp
[3432]	1011
[4447]	1012	DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
[3432]	1013	! volume elements
	1014	zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
	1015	zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
	1016	! part of the k-component of isopycnal momentum diffusive trends
	1017	zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu
	1018	zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv
	1019	! add the trends to the general trends
	1020	ua(ji,jj,jk) = ua(ji,jj,jk) + zuav
	1021	va(ji,jj,jk) = va(ji,jj,jk) + zvav
	1022	END DO
	1023	END DO
	1024	! ! ===============
	1025	END DO ! End of vertical slab
	1026	! ! ===============
	1027	#else
	1028	! ! ===============
	1029	DO jj = 2, jpjm1 ! Vertical slab
	1030	! ! ===============
	1031
[455]	1032
	1033	! I. vertical trends associated with the lateral mixing
	1034	! =====================================================
	1035	! (excluding the vertical flux proportional to dk[t]
	1036
	1037
	1038	! I.1 horizontal momentum gradient
	1039	! --------------------------------
	1040
	1041	DO jk = 1, jpk
	1042	DO ji = 2, jpi
	1043	! i-gradient of u at jj
	1044	zdiu (ji,jk) = tmask(ji,jj ,jk) * ( ub(ji,jj ,jk) - ub(ji-1,jj ,jk) )
	1045	! j-gradient of u and v at jj
	1046	zdju (ji,jk) = fmask(ji,jj ,jk) * ( ub(ji,jj+1,jk) - ub(ji ,jj ,jk) )
	1047	zdjv (ji,jk) = tmask(ji,jj ,jk) * ( vb(ji,jj ,jk) - vb(ji ,jj-1,jk) )
	1048	! j-gradient of u and v at jj+1
	1049	zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( ub(ji,jj ,jk) - ub(ji ,jj-1,jk) )
	1050	zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( vb(ji,jj+1,jk) - vb(ji ,jj ,jk) )
	1051	END DO
	1052	END DO
	1053	DO jk = 1, jpk
	1054	DO ji = 1, jpim1
	1055	! i-gradient of v at jj
	1056	zdiv (ji,jk) = fmask(ji,jj ,jk) * ( vb(ji+1,jj,jk) - vb(ji ,jj ,jk) )
	1057	END DO
	1058	END DO
	1059
	1060
	1061	! I.2 Vertical fluxes
	1062	! -------------------
	1063
	1064	! Surface and bottom vertical fluxes set to zero
	1065	DO ji = 1, jpi
	1066	zfuw(ji, 1 ) = 0.e0
	1067	zfvw(ji, 1 ) = 0.e0
	1068	zfuw(ji,jpk) = 0.e0
	1069	zfvw(ji,jpk) = 0.e0
	1070	END DO
	1071
	1072	! interior (2=<jk=<jpk-1) on U field
	1073	DO jk = 2, jpkm1
	1074	DO ji = 2, jpim1
	1075	zcoef0= 0.5 * aht0 * umask(ji,jj,jk)
	1076
	1077	zuwslpi = zcoef0 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
	1078	zuwslpj = zcoef0 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
	1079
	1080	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) &
	1081	+ tmask(ji,jj,jk )+tmask(ji+1,jj,jk ), 1. )
	1082	zmkf = 1./MAX( fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1) &
	1083	+ fmask(ji,jj-1,jk )+fmask(ji,jj,jk ), 1. )
	1084
	1085	zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
	1086	zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
	1087	! vertical flux on u field
	1088	zfuw(ji,jk) = zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) &
	1089	+zdiu (ji,jk ) + zdiu (ji+1,jk ) ) &
	1090	+ zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) &
	1091	+zdj1u(ji,jk ) + zdju (ji ,jk ) )
	1092	! update avmu (add isopycnal vertical coefficient to avmu)
	1093	! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
	1094	avmu(ji,jj,jk) = avmu(ji,jj,jk) + ( zuwslpi * zuwslpi + zuwslpj * zuwslpj ) / aht0
	1095	END DO
	1096	END DO
	1097
	1098	! interior (2=<jk=<jpk-1) on V field
	1099	DO jk = 2, jpkm1
	1100	DO ji = 2, jpim1
	1101	zcoef0= 0.5 * aht0 * vmask(ji,jj,jk)
	1102
	1103	zvwslpi = zcoef0 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
	1104	zvwslpj = zcoef0 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
	1105
	1106	zmkf = 1./MAX( fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1) &
	1107	+ fmask(ji-1,jj,jk )+fmask(ji,jj,jk ), 1. )
	1108	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1) &
	1109	+ tmask(ji,jj,jk )+tmask(ji,jj+1,jk ), 1. )
	1110
	1111	zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
	1112	zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
	1113	! vertical flux on v field
	1114	zfvw(ji,jk) = zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) &
	1115	+zdiv (ji,jk ) + zdiv (ji-1,jk ) ) &
	1116	+ zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) &
	1117	+zdjv (ji,jk ) + zdj1v(ji ,jk ) )
	1118	! update avmv (add isopycnal vertical coefficient to avmv)
	1119	! Caution: zcoef0 include aht0, so divided by aht0 to obtain slp^2 * aht0
	1120	avmv(ji,jj,jk) = avmv(ji,jj,jk) + ( zvwslpi * zvwslpi + zvwslpj * zvwslpj ) / aht0
	1121	END DO
	1122	END DO
	1123
	1124
	1125	! I.3 Divergence of vertical fluxes added to the general tracer trend
	1126	! -------------------------------------------------------------------
	1127	DO jk = 1, jpkm1
	1128	DO ji = 2, jpim1
	1129	! volume elements
	1130	zbu = e1u(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,jk)
	1131	zbv = e1v(ji,jj) * e2v(ji,jj) * fse3v(ji,jj,jk)
	1132	! part of the k-component of isopycnal momentum diffusive trends
	1133	zuav = ( zfuw(ji,jk) - zfuw(ji,jk+1) ) / zbu
	1134	zvav = ( zfvw(ji,jk) - zfvw(ji,jk+1) ) / zbv
	1135	! add the trends to the general trends
	1136	ua(ji,jj,jk) = ua(ji,jj,jk) + zuav
	1137	va(ji,jj,jk) = va(ji,jj,jk) + zvav
	1138	END DO
	1139	END DO
	1140	! ! ===============
	1141	END DO ! End of slab
	1142	! ! ===============
[3432]	1143	#endif
[4453]	1144	!CALL timing_stop('dyn_ldf_iso_vslab','section')
[4460]	1145	#if defined ARPDBGSUM
	1146	WRITE (,) narea,': ARPDBG: end dyn_ldf_iso SUM(ua)= ',SUM(ua),' at step=',kt
	1147	#endif
[455]	1148
[2715]	1149	IF( wrk_not_released(2, 1,2,3,4,5,6,7,8) ) CALL ctl_stop('dyn_ldf_iso: failed to release workspace arrays')
	1150	!
[3432]	1151	CALL timing_stop('dyn_ldf_iso','section')
	1152	!
[3]	1153	END SUBROUTINE dyn_ldf_iso
	1154
	1155	# else
	1156	!!----------------------------------------------------------------------
	1157	!! Dummy module NO rotation of mixing tensor
	1158	!!----------------------------------------------------------------------
	1159	CONTAINS
	1160	SUBROUTINE dyn_ldf_iso( kt ) ! Empty routine
[2715]	1161	INTEGER, INTENT(in) :: kt
[32]	1162	WRITE(,) 'dyn_ldf_iso: You should not have seen this print! error?', kt
[3]	1163	END SUBROUTINE dyn_ldf_iso
	1164	#endif
	1165
	1166	!!======================================================================
	1167	END MODULE dynldf_iso

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

Download in other formats: