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_bilapg.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_bilapg.F90 @ 4432

Last change on this file since 4432 was 4432, checked in by trackstand2, 10 years ago
mbkmax changes in dynldf_bilapg
Property svn:keywords set to `Id`
File size: 40.2 KB

Rev	Line
[3]	1	MODULE dynldf_bilapg
	2	!!======================================================================
	3	!! * MODULE dynldf_bilapg *
	4	!! Ocean dynamics: lateral viscosity trend
	5	!!======================================================================
[2715]	6	!! History : OPA ! 1997-07 (G. Madec) Original code
	7	!! NEMO 1.0 ! 2002-08 (G. Madec) F90: Free form and module
	8	!! 2.0 ! 2004-08 (C. Talandier) New trends organization
	9	!!----------------------------------------------------------------------
[3]	10	#if defined key_ldfslp \|\| defined key_esopa
	11	!!----------------------------------------------------------------------
	12	!! 'key_ldfslp' Rotation of mixing tensor
	13	!!----------------------------------------------------------------------
	14	!! dyn_ldf_bilapg : update the momentum trend with the horizontal part
	15	!! of the horizontal s-coord. bilaplacian diffusion
	16	!! ldfguv :
	17	!!----------------------------------------------------------------------
	18	USE oce ! ocean dynamics and tracers
	19	USE dom_oce ! ocean space and time domain
	20	USE ldfdyn_oce ! ocean dynamics lateral physics
	21	USE zdf_oce ! ocean vertical physics
[216]	22	USE trdmod ! ocean dynamics trends
	23	USE trdmod_oce ! ocean variables trends
[3]	24	USE ldfslp ! iso-neutral slopes available
[2715]	25	USE in_out_manager ! I/O manager
	26	USE lib_mpp ! MPP library
[3]	27	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
[258]	28	USE prtctl ! Print control
[3]	29
	30	IMPLICIT NONE
	31	PRIVATE
	32
[2715]	33	PUBLIC dyn_ldf_bilapg ! called by step.F90
[3]	34
[3432]	35	!FTRANS zfuw zfvw zdiu zdiv :I :z
[2715]	36	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zfvw , zdiu, zdiv ! 2D workspace (ldfguv)
[3432]	37	!FTRANS zdju zdj1u zdjv zdj1v :I :z
[2715]	38	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdju, zdj1u, zdjv, zdj1v ! 2D workspace (ldfguv)
	39
[3211]	40	!! * Control permutation of array indices
	41	# include "oce_ftrans.h90"
	42	# include "dom_oce_ftrans.h90"
	43	# include "ldfdyn_oce_ftrans.h90"
	44	# include "zdf_oce_ftrans.h90"
	45	# include "ldfslp_ftrans.h90"
	46
[3]	47	!! * Substitutions
	48	# include "domzgr_substitute.h90"
	49	# include "ldfdyn_substitute.h90"
	50	!!----------------------------------------------------------------------
[2528]	51	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]	52	!! $Id$
[2715]	53	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	54	!!----------------------------------------------------------------------
	55	CONTAINS
	56
[2715]	57	INTEGER FUNCTION dyn_ldf_bilapg_alloc()
	58	!!----------------------------------------------------------------------
	59	!! * ROUTINE dyn_ldf_bilapg_alloc *
	60	!!----------------------------------------------------------------------
	61	ALLOCATE( zfuw(jpi,jpk) , zfvw (jpi,jpk) , zdiu(jpi,jpk) , zdiv (jpi,jpk) , &
	62	& zdju(jpi,jpk) , zdj1u(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_bilapg_alloc )
	63	!
	64	IF( dyn_ldf_bilapg_alloc /= 0 ) CALL ctl_warn('dyn_ldf_bilapg_alloc: failed to allocate arrays')
	65	END FUNCTION dyn_ldf_bilapg_alloc
	66
	67
[3]	68	SUBROUTINE dyn_ldf_bilapg( kt )
	69	!!----------------------------------------------------------------------
	70	!! * ROUTINE dyn_ldf_bilapg *
	71	!!
	72	!! ** Purpose : Compute the before trend of the horizontal momentum
	73	!! diffusion and add it to the general trend of momentum equation.
	74	!!
	75	!! ** Method : The lateral momentum diffusive trends is provided by a
	76	!! a 4th order operator rotated along geopotential surfaces. It is
	77	!! computed using before fields (forward in time) and geopotential
	78	!! slopes computed in routine inildf.
	79	!! -1- compute the geopotential harmonic operator applied to
	80	!! (ub,vb) and multiply it by the eddy diffusivity coefficient
	81	!! (done by a call to ldfgpu and ldfgpv routines) The result is in
[216]	82	!! (zwk1,zwk2) arrays. Applied the domain lateral boundary conditions
[3]	83	!! by call to lbc_lnk.
[216]	84	!! -2- applied to (zwk1,zwk2) the geopotential harmonic operator
[3]	85	!! by a second call to ldfgpu and ldfgpv routines respectively. The
[216]	86	!! result is in (zwk3,zwk4) arrays.
[3]	87	!! -3- Add this trend to the general trend (ta,sa):
[216]	88	!! (ua,va) = (ua,va) + (zwk3,zwk4)
[3]	89	!! 'key_trddyn' defined: the trend is saved for diagnostics.
	90	!!
	91	!! ** Action : - Update (ua,va) arrays with the before geopotential
	92	!! biharmonic mixing trend.
[216]	93	!! - save the trend in (zwk3,zwk4) ('key_trddyn')
[3]	94	!!----------------------------------------------------------------------
[2715]	95	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	96	USE wrk_nemo, ONLY: zwk1 => wrk_3d_3 , zwk2 => wrk_3d_4 ! 3D workspace
	97	USE oce , ONLY: zwk3 => ta , zwk4 => sa ! ta, sa used as 3D workspace
[3211]	98	!! DCSE_NEMO: need additional directives for renamed module variables
	99	!FTRANS zwk1 :I :I :z
	100	!FTRANS zwk2 :I :I :z
	101	!FTRANS zwk3 :I :I :z
	102	!FTRANS zwk4 :I :I :z
[2715]	103	!
[3]	104	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[2715]	105	!
[3]	106	INTEGER :: ji, jj, jk ! dummy loop indices
	107	!!----------------------------------------------------------------------
	108
[2715]	109	IF( wrk_in_use(3, 3,4) ) THEN
	110	CALL ctl_stop('dyn_ldf_bilapg: requested workspace arrays unavailable') ; RETURN
	111	ENDIF
	112
[3]	113	IF( kt == nit000 ) THEN
	114	IF(lwp) WRITE(numout,*)
	115	IF(lwp) WRITE(numout,*) 'dyn_ldf_bilapg : horizontal biharmonic operator in s-coordinate'
	116	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
[3432]	117	zwk1(:,:,:) = 0.e0_wp ; zwk3(:,:,:) = 0.e0_wp
	118	zwk2(:,:,:) = 0.e0_wp ; zwk4(:,:,:) = 0.e0_wp
[2715]	119	! ! allocate dyn_ldf_bilapg arrays
	120	IF( dyn_ldf_bilapg_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays')
[3]	121	ENDIF
	122
	123	! Laplacian of (ub,vb) multiplied by ahm
	124	! --------------------------------------
[2715]	125	CALL ldfguv( ub, vb, zwk1, zwk2, 1 ) ! rotated harmonic operator applied to (ub,vb)
	126	! ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) )
	127	CALL lbc_lnk( zwk1, 'U', -1. ) ; CALL lbc_lnk( zwk2, 'V', -1. ) ! Lateral boundary conditions
[3]	128
	129	! Bilaplacian of (ub,vb)
	130	! ----------------------
[2715]	131	CALL ldfguv( zwk1, zwk2, zwk3, zwk4, 2 ) ! rotated harmonic operator applied to (zwk1,zwk2)
	132	! ! (output in (zwk3,zwk4) )
[3]	133
[2715]	134	! Update the momentum trends
[3]	135	! --------------------------
[3211]	136	#if defined key_z_first
[2715]	137	DO jj = 2, jpjm1 ! add the diffusive trend to the general momentum trends
[3211]	138	DO ji = 2, jpim1
[4432]	139	DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
[3211]	140	#else
	141	DO jj = 2, jpjm1 ! add the diffusive trend to the general momentum trends
[3]	142	DO jk = 1, jpkm1
	143	DO ji = 2, jpim1
[3211]	144	#endif
[216]	145	ua(ji,jj,jk) = ua(ji,jj,jk) + zwk3(ji,jj,jk)
	146	va(ji,jj,jk) = va(ji,jj,jk) + zwk4(ji,jj,jk)
[3]	147	END DO
	148	END DO
[2715]	149	END DO
	150	!
	151	IF( wrk_not_released(3, 3,4) ) CALL ctl_stop('dyn_ldf_bilapg: failed to release workspace arrays')
	152	!
[3]	153	END SUBROUTINE dyn_ldf_bilapg
	154
	155
	156	SUBROUTINE ldfguv( pu, pv, plu, plv, kahm )
	157	!!----------------------------------------------------------------------
	158	!! * ROUTINE ldfguv *
	159	!!
	160	!! ** Purpose : Apply a geopotential harmonic operator to (pu,pv)
	161	!! (defined at u- and v-points) and multiply it by the eddy
	162	!! viscosity coefficient (if kahm=1).
	163	!!
	164	!! ** Method : The harmonic operator rotated along geopotential
	165	!! surfaces is applied to (pu,pv) using the slopes of geopotential
	166	!! surfaces computed in inildf routine. The result is provided in
	167	!! (plu,plv) arrays. It is computed in 2 stepv:
	168	!!
	169	!! First step: horizontal part of the operator. It is computed on
	170	!! ========== pu as follows (idem on pv)
	171	!! horizontal fluxes :
	172	!! zftu = e2ue3u/e1u di[ pu ] - e2uuslp dk[ mi(mk(pu)) ]
	173	!! zftv = e1ve3v/e2v dj[ pu ] - e1vvslp dk[ mj(mk(pu)) ]
	174	!! take the horizontal divergence of the fluxes (no divided by
	175	!! the volume element :
	176	!! plu = di-1[ zftu ] + dj-1[ zftv ]
	177	!!
	178	!! Second step: vertical part of the operator. It is computed on
	179	!! =========== pu as follows (idem on pv)
	180	!! vertical fluxes :
	181	!! zftw = e1te2t/e3w (wslpi^2+wslpj^2) dk-1[ pu ]
	182	!! - e2t * wslpi di[ mi(mk(pu)) ]
	183	!! - e1t * wslpj dj[ mj(mk(pu)) ]
	184	!! take the vertical divergence of the fluxes add it to the hori-
	185	!! zontal component, divide the result by the volume element and
	186	!! if kahm=1, multiply by the eddy diffusivity coefficient:
	187	!! plu = aht / (e1te2te3t) { plu + dk[ zftw ] }
	188	!! else:
	189	!! plu = 1 / (e1te2te3t) { plu + dk[ zftw ] }
	190	!!
	191	!! ** Action :
	192	!! plu, plv : partial harmonic operator applied to
	193	!! pu and pv (all the components except
	194	!! second order vertical derivative term)
	195	!! 'key_trddyn' defined: the trend is saved for diagnostics.
[2715]	196	!!----------------------------------------------------------------------
	197	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
	198	USE wrk_nemo, ONLY: ziut => wrk_2d_1 , zjuf => wrk_2d_2 , zjvt => wrk_2d_3
[3432]	199	USE wrk_nemo, ONLY: zivf => wrk_2d_4
	200	#if ! defined key_z_first
	201	USE wrk_nemo, ONLY: zdku => wrk_2d_5 , zdk1u => wrk_2d_6
[2715]	202	USE wrk_nemo, ONLY: zdkv => wrk_2d_7 , zdk1v => wrk_2d_8
[3432]	203	#endif
	204	USE timing , ONLY: timing_start, timing_stop
[3]	205	!!
[3211]	206	!FTRANS pu :I :I :z
	207	!FTRANS pv :I :I :z
	208	!FTRANS plu :I :I :z
	209	!FTRANS plv :I :I :z
	210	!! DCSE_NEMO: work around deficiency in ftrans
	211	! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu , pv ! 1st call: before horizontal velocity
[2715]	212	! ! 2nd call: ahm x these fields
[4409]	213	REAL(wp), INTENT(in ) :: pu(jpi,jpj,jpkorig) , pv(jpi,jpj,jpkorig)
[3211]	214	! REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: plu, plv ! partial harmonic operator applied to
[2715]	215	! ! pu and pv (all the components except
	216	! ! second order vertical derivative term)
[4409]	217	REAL(wp), INTENT( out) :: plu(jpi,jpj,jpkorig), plv(jpi,jpj,jpkorig) ! partial harmonic operator applied to
[2715]	218	INTEGER , INTENT(in ) :: kahm ! =1 1st call ; =2 2nd call
	219	!
	220	INTEGER :: ji, jj, jk ! dummy loop indices
[3432]	221	INTEGER :: jif, jjf ! dummy loop indices over full domain
[2715]	222	REAL(wp) :: zabe1 , zabe2 , zcof1 , zcof2 ! local scalar
	223	REAL(wp) :: zcoef0, zcoef3, zcoef4 ! - -
	224	REAL(wp) :: zbur, zbvr, zmkt, zmkf, zuav, zvav ! - -
	225	REAL(wp) :: zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - -
[3432]	226	#if defined key_z_first
	227	! Can use scalars instead of work arrays when built with z-first
	228	REAL(wp) :: zdku, zdkv, zdk1u, zdk1v
	229	#endif
[3]	230	!!----------------------------------------------------------------------
	231
[3432]	232	CALL timing_start('ldfguv')
	233
[2715]	234	IF( wrk_in_use(2, 1,2,3,4,5,6,7,8) ) THEN
	235	CALL ctl_stop('dyn:ldfguv: requested workspace arrays unavailable') ; RETURN
	236	END IF
[3432]	237	CALL timing_start('ldfguv_1st')
	238
	239	#if defined key_z_first
[4432]	240	! ! ********** !
	241	! First step !
	242	! ! ********** !
[3432]	243	DO jj = 2, jpjm1
	244	DO ji = 2, jpim1
	245
	246	! Treat jk = 1 separately as is special case
	247	jk = 1
	248
	249	plu(ji,jj,jk) = &
	250	! ------------- ziut (ji+1, jj) -
	251	tmask(ji+1,jj,jk) * &
	252	( (e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) &
	253	! + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	254	+ (-e2t(ji+1,jj) / MAX( umask(ji,jj,jk )+umask(ji+1,jj,jk+1) &
	255	+ umask(ji,jj,jk+1)+umask(ji+1,jj,jk ), 1._wp ) &
	256	* 0.5 * ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) ) ) * &
	257	( ( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + &
	258	( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + &
	259	! +zdk1u(ji,jj) + zdku (ji-1,jj) ) )
	260	( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + &
	261	( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) &
	262	) ) - &
	263	! ------------- ziut (ji,jj ) +
	264	tmask(ji,jj,jk) * &
	265	( (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) &
	266	! + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	267	+ (-e2t(ji,jj) / MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
	268	+ umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1._wp ) &
	269	* 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * &
	270	( ( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
	271	( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + &
	272	! +zdk1u(ji,jj) + zdku (ji-1,jj) ) )
	273	( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
	274	( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) &
	275	) ) + &
	276	! ------------- zjuf (ji ,jj) -
	277	fmask(ji,jj,jk) * &
	278	( (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) &
	279	! + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &
	280	+ (-e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1) &
	281	+ umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. ) &
	282	* 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) ) * &
	283	( &
	284	( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) &
	285	+ ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj ,jk+1) &
	286	! +zdk1u(ji,jj+1) + zdku (ji,jj) ) )
	287	+ ( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) &
	288	+ ( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) &
	289	) ) - &
	290	! ------------- zjuf (ji,jj-1)
	291	fmask(ji,jj-1,jk) * &
	292	( (e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1) ) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) &
	293	! + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &
	294	+ (-e1f(ji,jj-1) /MAX(umask(ji,jj,jk)+umask(ji,jj-1,jk+1) &
	295	+ umask(ji,jj,jk+1)+umask(ji,jj-1,jk), 1. ) &
	296	* 0.5 * ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) ) ) * &
	297	( &
	298	( pu(ji,jj,jk) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) &
	299	+ ( pu(ji,jj-1 ,jk) - pu(ji,jj-1 ,jk+1) ) * umask(ji,jj-1 ,jk+1) &
	300	! +zdk1u(ji,jj+1) + zdku (ji,jj) ) )
	301	+ ( pu(ji,jj,jk) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) &
	302	+ ( pu(ji,jj-1 ,jk) - pu(ji,jj-1 ,jk+1) ) * umask(ji,jj-1,jk+1) &
	303	) )
	304
	305
	306	plv(ji,jj,jk) = &
	307	! ------------- zivf (ji,jj ) -
	308	fmask(ji,jj,jk) * &
	309	( (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) &
	310	! + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) &
	311	+ ((-e2f(ji,jj) / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
	312	+ vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. )) &
	313	* 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( &
	314	( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + &
	315	( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + &
	316	! +zdk1u(ji,jj) + zdku (ji+1,jj) ) )
	317	( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + &
	318	( pu(ji+1,jj,jk) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) &
	319	) ) - &
	320	! ------------- zivf (ji-1,jj) +
	321	fmask(ji-1,jj,jk) * &
	322	( (e2f(ji-1,jj) * fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) &
	323	! + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) &
	324	+ ((-e2f(ji-1,jj) / MAX( vmask(ji,jj,jk )+vmask(ji-1,jj,jk+1) &
	325	+ vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk ), 1. )) &
	326	* 0.5 * ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * ( &
	327	( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + &
	328	( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + &
	329	! +zdk1u(ji,jj) + zdku (ji+1,jj) ) )
	330	( pu(ji-1,jj,jk) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + &
	331	( pu(ji ,jj,jk) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) &
	332	) ) + &
	333	! ------------- zjvt (ji,jj+1) -
	334	tmask(ji,jj+1,jk) * &
	335	( (e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / e2t(ji,jj+1)) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) &
	336	! + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) &
	337	+ ((-e1t(ji,jj+1)/MAX( vmask(ji,jj,jk )+vmask(ji,jj+1,jk+1) &
	338	+ vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk ), 1. ) ) &
	339	* 0.5 * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * ( &
	340	( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj ,jk+1) + &
	341	( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) + &
	342	! +zdk1u(ji,jj-1) + zdku (ji,jj) ) )
	343	( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj ,jk+1) + &
	344	( pu(ji,jj+1,jk) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) &
	345	) ) - &
	346	! ------------- zjvt (ji,jj )
	347	tmask(ji,jj,jk) * &
	348	( (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) &
	349	! + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) &
	350	+ ((-e1t(ji,jj)/MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	351	+ vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) ) &
	352	* 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * ( &
	353	( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
	354	( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) + &
	355	! +zdk1u(ji,jj-1) + zdku (ji,jj) ) )
	356	( pu(ji,jj-1,jk) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
	357	( pu(ji,jj ,jk) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) &
	358	) )
	359
	360
[4432]	361	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
[3432]	362
	363	! plu(ji,jj,jk) = ziut (ji+1,jj) - &
	364	! ziut (ji,jj ) + &
	365	! zjuf (ji ,jj) - &
	366	! zjuf (ji,jj-1)
	367	plu(ji,jj,jk) = &
	368	! ------------- ziut (ji+1, jj ) -
	369	tmask(ji+1,jj,jk) * &
	370	( (e2t(ji+1,jj) * fse3t(ji+1,jj,jk) / e1t(ji+1,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) &
	371	! + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	372	+ (-e2t(ji+1,jj) / MAX(umask(ji,jj,jk)+umask(ji+1,jj,jk+1) &
	373	+ umask(ji,jj,jk+1)+umask(ji+1,jj,jk), 1._wp ) &
	374	* 0.5 * ( uslp(ji,jj,jk) + uslp(ji+1,jj,jk) ) ) * &
	375	( ( pu(ji+1,jj,jk-1) - pu(ji+1,jj,jk ) ) * umask(ji+1,jj,jk) + &
	376	( pu(ji ,jj,jk ) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + &
	377	! +zdk1u(ji,jj) + zdku (ji-1,jj) ) )
	378	( pu(ji+1,jj,jk ) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + &
	379	( pu(ji ,jj,jk-1) - pu(ji ,jj,jk ) ) * umask(ji ,jj,jk) &
	380	) ) - &
	381	! ------------- ziut (ji , jj ) +
	382	tmask(ji,jj,jk) * &
	383	( (e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) &
	384	! + zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	385	+ (-e2t(ji,jj) / MAX(umask(ji-1,jj,jk)+umask(ji,jj,jk+1) &
	386	+ umask(ji-1,jj,jk+1)+umask(ji,jj,jk), 1._wp ) &
	387	* 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) ) ) * &
	388	( ( pu(ji ,jj,jk-1) - pu(ji ,jj,jk ) ) * umask(ji,jj,jk) + &
	389	( pu(ji-1,jj,jk ) - pu(ji-1,jj,jk+1) ) * umask(ji-1,jj,jk+1) + &
	390	! +zdk1u(ji,jj) + zdku (ji-1,jj) ) )
	391	( pu(ji ,jj,jk ) - pu(ji ,jj,jk+1) ) * umask(ji,jj,jk+1) + &
	392	( pu(ji-1,jj,jk-1) - pu(ji-1,jj,jk ) ) * umask(ji-1,jj,jk) &
	393	) ) + &
	394	! ------------- zjuf (ji , jj ) -
	395	fmask(ji,jj,jk) * &
	396	( (e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj) ) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) &
	397	! + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &
	398	+ (-e1f(ji,jj) /MAX(umask(ji,jj+1,jk)+umask(ji,jj,jk+1) &
	399	+ umask(ji,jj+1,jk+1)+umask(ji,jj,jk), 1. ) &
	400	* 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) ) ) * &
	401	( &
	402	(pu(ji,jj+1,jk-1) - pu(ji,jj+1,jk ) ) * umask(ji,jj+1,jk) + &
	403	(pu(ji,jj ,jk ) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) + &
	404	! +zdk1u(ji,jj+1) + zdku (ji,jj) ) )
	405	(pu(ji,jj+1,jk ) - pu(ji,jj+1,jk+1) ) * umask(ji,jj+1,jk+1) + &
	406	(pu(ji,jj ,jk-1) - pu(ji,jj ,jk ) ) * umask(ji,jj,jk) &
	407	) ) - &
	408	! ------------- zjuf (ji , jj-1)
	409	fmask(ji,jj-1,jk) * &
	410	( (e1f(ji,jj-1) * fse3f(ji,jj-1,jk) / e2f(ji,jj-1) ) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) &
	411	! + zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &
	412	+ (-e1f(ji,jj-1) /MAX(umask(ji,jj,jk)+umask(ji,jj-1,jk+1) &
	413	+ umask(ji,jj,jk+1)+umask(ji,jj-1,jk), 1. ) &
	414	* 0.5 * ( vslp(ji+1,jj-1,jk) + vslp(ji,jj-1,jk) ) ) * &
	415	( &
	416	(pu(ji,jj,jk-1) - pu(ji,jj,jk ) ) * umask(ji,jj,jk) + &
	417	(pu(ji,jj-1,jk ) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
	418	! +zdk1u(ji,jj+1) + zdku (ji,jj) ) )
	419	(pu(ji,jj,jk ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) + &
	420	(pu(ji,jj-1,jk-1) - pu(ji,jj-1 ,jk ) ) * umask(ji,jj-1,jk) &
	421	) )
	422
	423
	424	! plv(ji,jj,jk) = zivf (ji,jj ) - &
	425	! zivf (ji-1,jj) + &
	426	! zjvt (ji,jj+1) - &
	427	! zjvt (ji,jj )
	428	plv(ji,jj,jk) = &
	429	! ------------- zivf (ji,jj ) -
	430	fmask(ji,jj,jk) * &
	431	( (e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)) * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) &
	432	! + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) &
	433	+ ((-e2f(ji,jj) / MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
	434	+ vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. )) &
	435	* 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )) * ( &
	436	( pu(ji ,jj,jk-1) - pu(ji ,jj,jk ) ) * umask(ji ,jj,jk ) + &
	437	( pu(ji+1,jj,jk ) - pu(ji+1,jj,jk+1) ) * umask(ji+1,jj,jk+1) + &
	438	! +zdk1u(ji,jj) + zdku (ji+1,jj) ) )
	439	( pu(ji ,jj,jk ) - pu(ji ,jj,jk+1) ) * umask(ji ,jj,jk+1) + &
	440	( pu(ji+1,jj,jk-1) - pu(ji+1,jj,jk ) ) * umask(ji+1,jj,jk ) &
	441	) ) - &
	442	! ------------- zivf (ji-1,jj) +
	443	fmask(ji-1,jj,jk) * &
	444	( (e2f(ji-1,jj) * fse3f(ji-1,jj,jk) / e1f(ji-1,jj)) * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) &
	445	! + zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) &
	446	+ ((-e2f(ji-1,jj) / MAX( vmask(ji,jj,jk )+vmask(ji-1,jj,jk+1) &
	447	+ vmask(ji,jj,jk+1)+vmask(ji-1,jj,jk ), 1. )) &
	448	* 0.5 * ( uslp(ji-1,jj+1,jk) + uslp(ji-1,jj,jk) )) * ( &
	449	( pu(ji-1 ,jj,jk-1) - pu(ji-1 ,jj,jk ) ) * umask(ji-1 ,jj,jk ) + &
	450	( pu(ji,jj,jk ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) + &
	451	! +zdk1u(ji,jj) + zdku (ji+1,jj) ) )
	452	( pu(ji-1 ,jj,jk ) - pu(ji-1 ,jj,jk+1) ) * umask(ji-1 ,jj,jk+1) + &
	453	( pu(ji,jj,jk-1) - pu(ji,jj,jk ) ) * umask(ji,jj,jk ) &
	454	) ) + &
	455	! ------------- zjvt (ji,jj+1) -
	456	tmask(ji,jj+1,jk) * &
	457	( (e1t(ji,jj+1) * fse3t(ji,jj+1,jk) / e2t(ji,jj+1)) * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) &
	458	! + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) &
	459	+ ((-e1t(ji,jj+1)/MAX( vmask(ji,jj,jk )+vmask(ji,jj+1,jk+1) &
	460	+ vmask(ji,jj,jk+1)+vmask(ji,jj+1,jk ), 1. ) ) &
	461	* 0.5 * ( vslp(ji,jj,jk) + vslp(ji,jj+1,jk) )) * &
	462	( &
	463	( pu(ji,jj,jk-1) - pu(ji,jj,jk ) ) * umask(ji,jj,jk) + &
	464	( pu(ji,jj+1 ,jk ) - pu(ji,jj+1 ,jk+1) ) * umask(ji,jj+1,jk+1) + &
	465	! +zdk1u(ji,jj-1) + zdku (ji,jj) ) )
	466	( pu(ji,jj,jk ) - pu(ji,jj,jk+1) ) * umask(ji,jj,jk+1) + &
	467	( pu(ji,jj+1 ,jk-1) - pu(ji,jj+1 ,jk ) ) * umask(ji,jj+1,jk) &
	468	) ) - &
	469	! ------------- zjvt (ji,jj )
	470	tmask(ji,jj,jk) * &
	471	( (e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)) * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) &
	472	! + zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) &
	473	+ ((-e1t(ji,jj)/MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	474	+ vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. ) ) &
	475	* 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )) * &
	476	( &
	477	( pu(ji,jj-1,jk-1) - pu(ji,jj-1,jk ) ) * umask(ji,jj-1,jk) + &
	478	( pu(ji,jj ,jk ) - pu(ji,jj ,jk+1) ) * umask(ji,jj,jk+1) + &
	479	! +zdk1u(ji,jj-1) + zdku (ji,jj) ) )
	480	( pu(ji,jj-1,jk ) - pu(ji,jj-1,jk+1) ) * umask(ji,jj-1,jk+1) + &
	481	( pu(ji,jj ,jk-1) - pu(ji,jj ,jk ) ) * umask(ji,jj,jk) ) )
	482
	483	END DO
	484	END DO
	485
[4432]	486	!
	487	END DO
	488
[3432]	489	#else
	490	! ! ********** ! ! ===============
[3]	491	DO jk = 1, jpkm1 ! First step ! ! Horizontal slab
	492	! ! ********** ! ! ===============
	493
	494	! I.1 Vertical gradient of pu and pv at level jk and jk+1
	495	! -------------------------------------------------------
	496	! surface boundary condition: zdku(jk=1)=zdku(jk=2)
	497	! zdkv(jk=1)=zdkv(jk=2)
	498
	499	zdk1u(:,:) = ( pu(:,:,jk) - pu(:,:,jk+1) ) * umask(:,:,jk+1)
	500	zdk1v(:,:) = ( pv(:,:,jk) - pv(:,:,jk+1) ) * vmask(:,:,jk+1)
	501
	502	IF( jk == 1 ) THEN
	503	zdku(:,:) = zdk1u(:,:)
	504	zdkv(:,:) = zdk1v(:,:)
	505	ELSE
	506	zdku(:,:) = ( pu(:,:,jk-1) - pu(:,:,jk) ) * umask(:,:,jk)
	507	zdkv(:,:) = ( pv(:,:,jk-1) - pv(:,:,jk) ) * vmask(:,:,jk)
	508	ENDIF
	509
	510	! -----f-----
	511	! I.2 Horizontal fluxes on U \|
	512	! ------------------------=== t u t
	513	! \|
	514	! i-flux at t-point -----f-----
	515	DO jj = 1, jpjm1
	516	DO ji = 2, jpi
	517	zabe1 = e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
	518
	519	zmkt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
	520	+ umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )
	521
	522	zcof1 = -e2t(ji,jj) * zmkt &
	523	* 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
	524
	525	ziut(ji,jj) = tmask(ji,jj,jk) * &
	526	( zabe1 * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) &
	527	+ zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	528	+zdk1u(ji,jj) + zdku (ji-1,jj) ) )
	529	END DO
	530	END DO
	531
	532	! j-flux at f-point
	533	DO jj = 1, jpjm1
	534	DO ji = 1, jpim1
	535	zabe2 = e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
	536
	537	zmkf = 1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) &
	538	+ umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. )
	539
	540	zcof2 = -e1f(ji,jj) * zmkf &
	541	* 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
	542
	543	zjuf(ji,jj) = fmask(ji,jj,jk) * &
	544	( zabe2 * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) &
	545	+ zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &
	546	+zdk1u(ji,jj+1) + zdku (ji,jj) ) )
	547	END DO
	548	END DO
	549
	550	! \| t \|
	551	! I.3 Horizontal fluxes on V \| \|
	552	! ------------------------=== f---v---f
	553	! \| \|
	554	! i-flux at f-point \| t \|
	555	DO jj = 1, jpjm1
	556	DO ji = 1, jpim1
	557	zabe1 = e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
	558
	559	zmkf = 1./MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
	560	+ vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. )
	561
	562	zcof1 = -e2f(ji,jj) * zmkf &
	563	* 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
	564
	565	zivf(ji,jj) = fmask(ji,jj,jk) * &
	566	( zabe1 * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) &
	567	+ zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) &
	568	+zdk1u(ji,jj) + zdku (ji+1,jj) ) )
	569	END DO
	570	END DO
	571
	572	! j-flux at t-point
	573	DO jj = 2, jpj
	574	DO ji = 1, jpim1
	575	zabe2 = e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
	576
	577	zmkt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	578	+ vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. )
	579
	580	zcof2 = -e1t(ji,jj) * zmkt &
	581	* 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
	582
	583	zjvt(ji,jj) = tmask(ji,jj,jk) * &
	584	( zabe2 * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) &
	585	+ zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) &
	586	+zdk1u(ji,jj-1) + zdku (ji,jj) ) )
	587	END DO
	588	END DO
	589
	590
	591	! I.4 Second derivative (divergence) (not divided by the volume)
	592	! ---------------------
	593
	594	DO jj = 2, jpjm1
	595	DO ji = 2, jpim1
	596	plu(ji,jj,jk) = ziut (ji+1,jj) - ziut (ji,jj ) &
	597	+ zjuf (ji ,jj) - zjuf (ji,jj-1)
	598	plv(ji,jj,jk) = zivf (ji,jj ) - zivf (ji-1,jj) &
	599	+ zjvt (ji,jj+1) - zjvt (ji,jj )
	600	END DO
	601	END DO
	602
	603	! ! ===============
	604	END DO ! End of slab
	605	! ! ===============
[3432]	606	#endif
	607	CALL timing_stop('ldfguv_1st','section')
[3]	608
	609	!,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
	610
[3432]	611	CALL timing_start('ldfguv_2nd')
[3]	612	! ! ************ ! ! ===============
	613	DO jj = 2, jpjm1 ! Second step ! ! Horizontal slab
	614	! ! ************ ! ! ===============
	615
	616	! II.1 horizontal (pu,pv) gradients
	617	! ---------------------------------
	618
[3432]	619	#if defined key_z_first
	620	DO ji = 2, jpi
[4432]	621	DO jk = 1, mbkmax(ji,jj) ! jpk
[3432]	622	! i-gradient of u at jj
	623	zdiu (ji,jk) = tmask(ji,jj ,jk) * ( pu(ji,jj ,jk) - pu(ji-1,jj ,jk) )
	624	! j-gradient of u and v at jj
	625	zdju (ji,jk) = fmask(ji,jj ,jk) * ( pu(ji,jj+1,jk) - pu(ji ,jj ,jk) )
	626	zdjv (ji,jk) = tmask(ji,jj ,jk) * ( pv(ji,jj ,jk) - pv(ji ,jj-1,jk) )
	627	! j-gradient of u and v at jj+1
	628	zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj ,jk) - pu(ji ,jj-1,jk) )
	629	zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji ,jj ,jk) )
	630	END DO
	631	END DO
	632	DO ji = 1, jpim1
[4432]	633	DO jk = 1, mbkmax(ji,jj) ! jpk
[3432]	634	! i-gradient of v at jj
	635	zdiv (ji,jk) = fmask(ji,jj ,jk) * ( pv(ji+1,jj,jk) - pv(ji ,jj ,jk) )
	636	END DO
	637	END DO
	638	#else
[3]	639	DO jk = 1, jpk
	640	DO ji = 2, jpi
	641	! i-gradient of u at jj
	642	zdiu (ji,jk) = tmask(ji,jj ,jk) * ( pu(ji,jj ,jk) - pu(ji-1,jj ,jk) )
	643	! j-gradient of u and v at jj
	644	zdju (ji,jk) = fmask(ji,jj ,jk) * ( pu(ji,jj+1,jk) - pu(ji ,jj ,jk) )
	645	zdjv (ji,jk) = tmask(ji,jj ,jk) * ( pv(ji,jj ,jk) - pv(ji ,jj-1,jk) )
	646	! j-gradient of u and v at jj+1
	647	zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj ,jk) - pu(ji ,jj-1,jk) )
	648	zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji ,jj ,jk) )
	649	END DO
	650	END DO
	651	DO jk = 1, jpk
	652	DO ji = 1, jpim1
	653	! i-gradient of v at jj
	654	zdiv (ji,jk) = fmask(ji,jj ,jk) * ( pv(ji+1,jj,jk) - pv(ji ,jj ,jk) )
	655	END DO
	656	END DO
[3432]	657	#endif
[3]	658
	659	! II.2 Vertical fluxes
	660	! --------------------
	661
	662	! Surface and bottom vertical fluxes set to zero
	663
[4432]	664	#if defined key_z_first
	665	DO ji=1, jpi
	666	zfuw(ji, 1 ) = 0.e0
	667	zfvw(ji, 1 ) = 0.e0
	668	zfuw(ji,mbkmax(ji,jj):jpk) = 0.e0
	669	zfvw(ji,mbkmax(ji,jj):jpk) = 0.e0
	670	END DO
	671	#else
[3]	672	zfuw(:, 1 ) = 0.e0
	673	zfvw(:, 1 ) = 0.e0
	674	zfuw(:,jpk) = 0.e0
	675	zfvw(:,jpk) = 0.e0
[4432]	676	#endif
[3]	677	! interior (2=<jk=<jpk-1) on pu field
	678
[3432]	679	#if defined key_z_first
	680	DO ji = 2, jpim1
[4432]	681	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
[3432]	682	#else
[3]	683	DO jk = 2, jpkm1
	684	DO ji = 2, jpim1
[3432]	685	#endif
[3]	686	! i- and j-slopes at uw-point
	687	zuwslpi = 0.5 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
	688	zuwslpj = 0.5 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
	689	! coef. for the vertical dirative
	690	zcoef0 = e1u(ji,jj) * e2u(ji,jj) / fse3u(ji,jj,jk) &
	691	* ( zuwslpi * zuwslpi + zuwslpj * zuwslpj )
	692	! weights for the i-k, j-k averaging at t- and f-points, resp.
	693	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) &
	694	+ tmask(ji,jj,jk )+tmask(ji+1,jj,jk ), 1. )
	695	zmkf = 1./MAX( fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1) &
	696	+ fmask(ji,jj-1,jk )+fmask(ji,jj,jk ), 1. )
	697	! coef. for the horitontal derivative
	698	zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
	699	zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
	700	! vertical flux on u field
	701	zfuw(ji,jk) = umask(ji,jj,jk) * &
	702	( zcoef0 * ( pu (ji,jj,jk-1) - pu (ji,jj,jk) ) &
	703	+ zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) &
	704	+zdiu (ji,jk ) + zdiu (ji+1,jk ) ) &
	705	+ zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) &
	706	+zdj1u(ji,jk ) + zdju (ji ,jk ) ) )
	707	END DO
	708	END DO
	709
	710	! interior (2=<jk=<jpk-1) on pv field
	711
[3432]	712	#if defined key_z_first
	713	DO ji = 2, jpim1
[4432]	714	DO jk = 2, mbkmax(ji,jj)-1 ! jpkm1
[3432]	715	#else
[3]	716	DO jk = 2, jpkm1
	717	DO ji = 2, jpim1
[3432]	718	#endif
[3]	719	! i- and j-slopes at vw-point
	720	zvwslpi = 0.5 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
	721	zvwslpj = 0.5 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
	722	! coef. for the vertical derivative
	723	zcoef0 = e1v(ji,jj) * e2v(ji,jj) / fse3v(ji,jj,jk) &
	724	* ( zvwslpi * zvwslpi + zvwslpj * zvwslpj )
	725	! weights for the i-k, j-k averaging at f- and t-points, resp.
	726	zmkf = 1./MAX( fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1) &
	727	+ fmask(ji-1,jj,jk )+fmask(ji,jj,jk ), 1. )
	728	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1) &
	729	+ tmask(ji,jj,jk )+tmask(ji,jj+1,jk ), 1. )
	730	! coef. for the horizontal derivatives
	731	zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
	732	zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
	733	! vertical flux on pv field
	734	zfvw(ji,jk) = vmask(ji,jj,jk) * &
	735	( zcoef0 * ( pv (ji,jj,jk-1) - pv (ji,jj,jk) ) &
	736	+ zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) &
	737	+zdiv (ji,jk ) + zdiv (ji-1,jk ) ) &
	738	+ zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) &
	739	+zdjv (ji,jk ) + zdj1v(ji ,jk ) ) )
	740	END DO
	741	END DO
	742
	743
	744	! II.3 Divergence of vertical fluxes added to the horizontal divergence
	745	! ---------------------------------------------------------------------
	746
	747	IF( kahm == 1 ) THEN
	748	! multiply the laplacian by the eddy viscosity coefficient
[3432]	749	#if defined key_z_first
	750	DO ji = 2, jpim1
[4432]	751	DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
[3432]	752	#else
[3]	753	DO jk = 1, jpkm1
	754	DO ji = 2, jpim1
[3432]	755	#endif
[3]	756	! eddy coef. divided by the volume element
	757	zbur = fsahmu(ji,jj,jk) / ( e1u(ji,jj)e2u(ji,jj)fse3u(ji,jj,jk) )
	758	zbvr = fsahmv(ji,jj,jk) / ( e1v(ji,jj)e2v(ji,jj)fse3v(ji,jj,jk) )
	759	! vertical divergence
	760	zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
	761	zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
	762	! harmonic operator applied to (pu,pv) and multiply by ahm
	763	plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
	764	plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
	765	END DO
	766	END DO
	767	ELSEIF( kahm == 2 ) THEN
	768	! second call, no multiplication
[3432]	769	#if defined key_z_first
	770	DO ji = 2, jpim1
[4432]	771	DO jk = 1, mbkmax(ji,jj)-1 ! jpkm1
[3432]	772	#else
[3]	773	DO jk = 1, jpkm1
	774	DO ji = 2, jpim1
[3432]	775	#endif
[3]	776	! inverse of the volume element
	777	zbur = 1. / ( e1u(ji,jj)e2u(ji,jj)fse3u(ji,jj,jk) )
	778	zbvr = 1. / ( e1v(ji,jj)e2v(ji,jj)fse3v(ji,jj,jk) )
	779	! vertical divergence
	780	zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
	781	zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
	782	! harmonic operator applied to (pu,pv)
	783	plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
	784	plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
	785	END DO
	786	END DO
	787	ELSE
	788	IF(lwp)WRITE(numout,*) ' ldfguv: kahm= 1 or 2, here =', kahm
	789	IF(lwp)WRITE(numout,*) ' We stop'
	790	STOP 'ldfguv'
	791	ENDIF
	792	! ! ===============
	793	END DO ! End of slab
	794	! ! ===============
[3432]	795	CALL timing_stop('ldfguv_2nd','section')
[2715]	796
	797	IF( wrk_not_released(2, 1,2,3,4,5,6,7,8) ) CALL ctl_stop('dyn:ldfguv: failed to release workspace arrays')
	798	!
[3432]	799	CALL timing_stop('ldfguv','section')
	800
[3]	801	END SUBROUTINE ldfguv
	802
	803	#else
	804	!!----------------------------------------------------------------------
	805	!! Dummy module : NO rotation of mixing tensor
	806	!!----------------------------------------------------------------------
	807	CONTAINS
	808	SUBROUTINE dyn_ldf_bilapg( kt ) ! Dummy routine
[2715]	809	INTEGER, INTENT(in) :: kt
[32]	810	WRITE(,) 'dyn_ldf_bilapg: You should not have seen this print! error?', kt
[3]	811	END SUBROUTINE dyn_ldf_bilapg
	812	#endif
	813
	814	!!======================================================================
	815	END MODULE dynldf_bilapg

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

Download in other formats: