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

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

Download in other formats: