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 trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynldf_bilapg.F90 @ 3434

Last change on this file since 3434 was 3294, checked in by rblod, 12 years ago
Merge of 3.4beta into the trunk
Property svn:keywords set to `Id`
File size: 22.3 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
[3294]	29	USE wrk_nemo ! Memory Allocation
	30	USE timing ! Timing
[3]	31
[3294]	32
[3]	33	IMPLICIT NONE
	34	PRIVATE
	35
[2715]	36	PUBLIC dyn_ldf_bilapg ! called by step.F90
[3]	37
[2715]	38	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zfuw, zfvw , zdiu, zdiv ! 2D workspace (ldfguv)
	39	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: zdju, zdj1u, zdjv, zdj1v ! 2D workspace (ldfguv)
	40
[3]	41	!! * Substitutions
	42	# include "domzgr_substitute.h90"
	43	# include "ldfdyn_substitute.h90"
	44	!!----------------------------------------------------------------------
[2528]	45	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[1152]	46	!! $Id$
[2715]	47	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	48	!!----------------------------------------------------------------------
	49	CONTAINS
	50
[2715]	51	INTEGER FUNCTION dyn_ldf_bilapg_alloc()
	52	!!----------------------------------------------------------------------
	53	!! * ROUTINE dyn_ldf_bilapg_alloc *
	54	!!----------------------------------------------------------------------
	55	ALLOCATE( zfuw(jpi,jpk) , zfvw (jpi,jpk) , zdiu(jpi,jpk) , zdiv (jpi,jpk) , &
	56	& zdju(jpi,jpk) , zdj1u(jpi,jpk) , zdjv(jpi,jpk) , zdj1v(jpi,jpk) , STAT=dyn_ldf_bilapg_alloc )
	57	!
	58	IF( dyn_ldf_bilapg_alloc /= 0 ) CALL ctl_warn('dyn_ldf_bilapg_alloc: failed to allocate arrays')
	59	END FUNCTION dyn_ldf_bilapg_alloc
	60
	61
[3]	62	SUBROUTINE dyn_ldf_bilapg( kt )
	63	!!----------------------------------------------------------------------
	64	!! * ROUTINE dyn_ldf_bilapg *
	65	!!
	66	!! ** Purpose : Compute the before trend of the horizontal momentum
	67	!! diffusion and add it to the general trend of momentum equation.
	68	!!
	69	!! ** Method : The lateral momentum diffusive trends is provided by a
	70	!! a 4th order operator rotated along geopotential surfaces. It is
	71	!! computed using before fields (forward in time) and geopotential
	72	!! slopes computed in routine inildf.
	73	!! -1- compute the geopotential harmonic operator applied to
	74	!! (ub,vb) and multiply it by the eddy diffusivity coefficient
	75	!! (done by a call to ldfgpu and ldfgpv routines) The result is in
[216]	76	!! (zwk1,zwk2) arrays. Applied the domain lateral boundary conditions
[3]	77	!! by call to lbc_lnk.
[216]	78	!! -2- applied to (zwk1,zwk2) the geopotential harmonic operator
[3]	79	!! by a second call to ldfgpu and ldfgpv routines respectively. The
[216]	80	!! result is in (zwk3,zwk4) arrays.
[3]	81	!! -3- Add this trend to the general trend (ta,sa):
[216]	82	!! (ua,va) = (ua,va) + (zwk3,zwk4)
[3]	83	!! 'key_trddyn' defined: the trend is saved for diagnostics.
	84	!!
	85	!! ** Action : - Update (ua,va) arrays with the before geopotential
	86	!! biharmonic mixing trend.
[216]	87	!! - save the trend in (zwk3,zwk4) ('key_trddyn')
[3]	88	!!----------------------------------------------------------------------
	89	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[2715]	90	!
[3]	91	INTEGER :: ji, jj, jk ! dummy loop indices
[3294]	92	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwk1, zwk2, zwk3, zwk4
[3]	93	!!----------------------------------------------------------------------
[3294]	94	!
	95	IF( nn_timing == 1 ) CALL timing_start('dyn_ldf_bilapg')
	96	!
	97	CALL wrk_alloc( jpi, jpj, jpk, zwk1, zwk2, zwk3, zwk4 )
	98	!
[3]	99	IF( kt == nit000 ) THEN
	100	IF(lwp) WRITE(numout,*)
	101	IF(lwp) WRITE(numout,*) 'dyn_ldf_bilapg : horizontal biharmonic operator in s-coordinate'
	102	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~'
[2715]	103	! ! allocate dyn_ldf_bilapg arrays
	104	IF( dyn_ldf_bilapg_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_ldf_bilapg: failed to allocate arrays')
[3]	105	ENDIF
[3294]	106	!
	107	zwk1(:,:,:) = 0.e0 ; zwk3(:,:,:) = 0.e0
	108	zwk2(:,:,:) = 0.e0 ; zwk4(:,:,:) = 0.e0
[3]	109
	110	! Laplacian of (ub,vb) multiplied by ahm
	111	! --------------------------------------
[2715]	112	CALL ldfguv( ub, vb, zwk1, zwk2, 1 ) ! rotated harmonic operator applied to (ub,vb)
	113	! ! and multiply by ahmu, ahmv (output in (zwk1,zwk2) )
	114	CALL lbc_lnk( zwk1, 'U', -1. ) ; CALL lbc_lnk( zwk2, 'V', -1. ) ! Lateral boundary conditions
[3]	115
	116	! Bilaplacian of (ub,vb)
	117	! ----------------------
[2715]	118	CALL ldfguv( zwk1, zwk2, zwk3, zwk4, 2 ) ! rotated harmonic operator applied to (zwk1,zwk2)
	119	! ! (output in (zwk3,zwk4) )
[3]	120
[2715]	121	! Update the momentum trends
[3]	122	! --------------------------
[2715]	123	DO jj = 2, jpjm1 ! add the diffusive trend to the general momentum trends
[3]	124	DO jk = 1, jpkm1
	125	DO ji = 2, jpim1
[216]	126	ua(ji,jj,jk) = ua(ji,jj,jk) + zwk3(ji,jj,jk)
	127	va(ji,jj,jk) = va(ji,jj,jk) + zwk4(ji,jj,jk)
[3]	128	END DO
	129	END DO
[2715]	130	END DO
	131	!
[3294]	132	CALL wrk_dealloc( jpi, jpj, jpk, zwk1, zwk2, zwk3, zwk4 )
[2715]	133	!
[3294]	134	IF( nn_timing == 1 ) CALL timing_stop('dyn_ldf_bilapg')
	135	!
[3]	136	END SUBROUTINE dyn_ldf_bilapg
	137
	138
	139	SUBROUTINE ldfguv( pu, pv, plu, plv, kahm )
	140	!!----------------------------------------------------------------------
	141	!! * ROUTINE ldfguv *
	142	!!
	143	!! ** Purpose : Apply a geopotential harmonic operator to (pu,pv)
	144	!! (defined at u- and v-points) and multiply it by the eddy
	145	!! viscosity coefficient (if kahm=1).
	146	!!
	147	!! ** Method : The harmonic operator rotated along geopotential
	148	!! surfaces is applied to (pu,pv) using the slopes of geopotential
	149	!! surfaces computed in inildf routine. The result is provided in
	150	!! (plu,plv) arrays. It is computed in 2 stepv:
	151	!!
	152	!! First step: horizontal part of the operator. It is computed on
	153	!! ========== pu as follows (idem on pv)
	154	!! horizontal fluxes :
	155	!! zftu = e2ue3u/e1u di[ pu ] - e2uuslp dk[ mi(mk(pu)) ]
	156	!! zftv = e1ve3v/e2v dj[ pu ] - e1vvslp dk[ mj(mk(pu)) ]
	157	!! take the horizontal divergence of the fluxes (no divided by
	158	!! the volume element :
	159	!! plu = di-1[ zftu ] + dj-1[ zftv ]
	160	!!
	161	!! Second step: vertical part of the operator. It is computed on
	162	!! =========== pu as follows (idem on pv)
	163	!! vertical fluxes :
	164	!! zftw = e1te2t/e3w (wslpi^2+wslpj^2) dk-1[ pu ]
	165	!! - e2t * wslpi di[ mi(mk(pu)) ]
	166	!! - e1t * wslpj dj[ mj(mk(pu)) ]
	167	!! take the vertical divergence of the fluxes add it to the hori-
	168	!! zontal component, divide the result by the volume element and
	169	!! if kahm=1, multiply by the eddy diffusivity coefficient:
	170	!! plu = aht / (e1te2te3t) { plu + dk[ zftw ] }
	171	!! else:
	172	!! plu = 1 / (e1te2te3t) { plu + dk[ zftw ] }
	173	!!
	174	!! ** Action :
	175	!! plu, plv : partial harmonic operator applied to
	176	!! pu and pv (all the components except
	177	!! second order vertical derivative term)
	178	!! 'key_trddyn' defined: the trend is saved for diagnostics.
[2715]	179	!!----------------------------------------------------------------------
[3]	180	!!
[2715]	181	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: pu , pv ! 1st call: before horizontal velocity
	182	! ! 2nd call: ahm x these fields
	183	REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: plu, plv ! partial harmonic operator applied to
	184	! ! pu and pv (all the components except
	185	! ! second order vertical derivative term)
	186	INTEGER , INTENT(in ) :: kahm ! =1 1st call ; =2 2nd call
	187	!
	188	INTEGER :: ji, jj, jk ! dummy loop indices
	189	REAL(wp) :: zabe1 , zabe2 , zcof1 , zcof2 ! local scalar
	190	REAL(wp) :: zcoef0, zcoef3, zcoef4 ! - -
	191	REAL(wp) :: zbur, zbvr, zmkt, zmkf, zuav, zvav ! - -
	192	REAL(wp) :: zuwslpi, zuwslpj, zvwslpi, zvwslpj ! - -
[3294]	193	!
	194	REAL(wp), POINTER, DIMENSION(:,:) :: ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v
[3]	195	!!----------------------------------------------------------------------
[3294]	196	!
	197	IF( nn_timing == 1 ) CALL timing_start('ldfguv')
	198	!
	199	CALL wrk_alloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v )
	200	!
[3]	201	! ! ********** ! ! ===============
	202	DO jk = 1, jpkm1 ! First step ! ! Horizontal slab
	203	! ! ********** ! ! ===============
	204
	205	! I.1 Vertical gradient of pu and pv at level jk and jk+1
	206	! -------------------------------------------------------
	207	! surface boundary condition: zdku(jk=1)=zdku(jk=2)
	208	! zdkv(jk=1)=zdkv(jk=2)
	209
	210	zdk1u(:,:) = ( pu(:,:,jk) - pu(:,:,jk+1) ) * umask(:,:,jk+1)
	211	zdk1v(:,:) = ( pv(:,:,jk) - pv(:,:,jk+1) ) * vmask(:,:,jk+1)
	212
	213	IF( jk == 1 ) THEN
	214	zdku(:,:) = zdk1u(:,:)
	215	zdkv(:,:) = zdk1v(:,:)
	216	ELSE
	217	zdku(:,:) = ( pu(:,:,jk-1) - pu(:,:,jk) ) * umask(:,:,jk)
	218	zdkv(:,:) = ( pv(:,:,jk-1) - pv(:,:,jk) ) * vmask(:,:,jk)
	219	ENDIF
	220
	221	! -----f-----
	222	! I.2 Horizontal fluxes on U \|
	223	! ------------------------=== t u t
	224	! \|
	225	! i-flux at t-point -----f-----
	226	DO jj = 1, jpjm1
	227	DO ji = 2, jpi
	228	zabe1 = e2t(ji,jj) * fse3t(ji,jj,jk) / e1t(ji,jj)
	229
	230	zmkt = 1./MAX( umask(ji-1,jj,jk )+umask(ji,jj,jk+1) &
	231	+ umask(ji-1,jj,jk+1)+umask(ji,jj,jk ), 1. )
	232
	233	zcof1 = -e2t(ji,jj) * zmkt &
	234	* 0.5 * ( uslp(ji-1,jj,jk) + uslp(ji,jj,jk) )
	235
	236	ziut(ji,jj) = tmask(ji,jj,jk) * &
	237	( zabe1 * ( pu(ji,jj,jk) - pu(ji-1,jj,jk) ) &
	238	+ zcof1 * ( zdku (ji,jj) + zdk1u(ji-1,jj) &
	239	+zdk1u(ji,jj) + zdku (ji-1,jj) ) )
	240	END DO
	241	END DO
	242
	243	! j-flux at f-point
	244	DO jj = 1, jpjm1
	245	DO ji = 1, jpim1
	246	zabe2 = e1f(ji,jj) * fse3f(ji,jj,jk) / e2f(ji,jj)
	247
	248	zmkf = 1./MAX( umask(ji,jj+1,jk )+umask(ji,jj,jk+1) &
	249	+ umask(ji,jj+1,jk+1)+umask(ji,jj,jk ), 1. )
	250
	251	zcof2 = -e1f(ji,jj) * zmkf &
	252	* 0.5 * ( vslp(ji+1,jj,jk) + vslp(ji,jj,jk) )
	253
	254	zjuf(ji,jj) = fmask(ji,jj,jk) * &
	255	( zabe2 * ( pu(ji,jj+1,jk) - pu(ji,jj,jk) ) &
	256	+ zcof2 * ( zdku (ji,jj+1) + zdk1u(ji,jj) &
	257	+zdk1u(ji,jj+1) + zdku (ji,jj) ) )
	258	END DO
	259	END DO
	260
	261	! \| t \|
	262	! I.3 Horizontal fluxes on V \| \|
	263	! ------------------------=== f---v---f
	264	! \| \|
	265	! i-flux at f-point \| t \|
	266	DO jj = 1, jpjm1
	267	DO ji = 1, jpim1
	268	zabe1 = e2f(ji,jj) * fse3f(ji,jj,jk) / e1f(ji,jj)
	269
	270	zmkf = 1./MAX( vmask(ji+1,jj,jk )+vmask(ji,jj,jk+1) &
	271	+ vmask(ji+1,jj,jk+1)+vmask(ji,jj,jk ), 1. )
	272
	273	zcof1 = -e2f(ji,jj) * zmkf &
	274	* 0.5 * ( uslp(ji,jj+1,jk) + uslp(ji,jj,jk) )
	275
	276	zivf(ji,jj) = fmask(ji,jj,jk) * &
	277	( zabe1 * ( pu(ji+1,jj,jk) - pu(ji,jj,jk) ) &
	278	+ zcof1 * ( zdku (ji,jj) + zdk1u(ji+1,jj) &
	279	+zdk1u(ji,jj) + zdku (ji+1,jj) ) )
	280	END DO
	281	END DO
	282
	283	! j-flux at t-point
	284	DO jj = 2, jpj
	285	DO ji = 1, jpim1
	286	zabe2 = e1t(ji,jj) * fse3t(ji,jj,jk) / e2t(ji,jj)
	287
	288	zmkt = 1./MAX( vmask(ji,jj-1,jk )+vmask(ji,jj,jk+1) &
	289	+ vmask(ji,jj-1,jk+1)+vmask(ji,jj,jk ), 1. )
	290
	291	zcof2 = -e1t(ji,jj) * zmkt &
	292	* 0.5 * ( vslp(ji,jj-1,jk) + vslp(ji,jj,jk) )
	293
	294	zjvt(ji,jj) = tmask(ji,jj,jk) * &
	295	( zabe2 * ( pu(ji,jj,jk) - pu(ji,jj-1,jk) ) &
	296	+ zcof2 * ( zdku (ji,jj-1) + zdk1u(ji,jj) &
	297	+zdk1u(ji,jj-1) + zdku (ji,jj) ) )
	298	END DO
	299	END DO
	300
	301
	302	! I.4 Second derivative (divergence) (not divided by the volume)
	303	! ---------------------
	304
	305	DO jj = 2, jpjm1
	306	DO ji = 2, jpim1
	307	plu(ji,jj,jk) = ziut (ji+1,jj) - ziut (ji,jj ) &
	308	+ zjuf (ji ,jj) - zjuf (ji,jj-1)
	309	plv(ji,jj,jk) = zivf (ji,jj ) - zivf (ji-1,jj) &
	310	+ zjvt (ji,jj+1) - zjvt (ji,jj )
	311	END DO
	312	END DO
	313
	314	! ! ===============
	315	END DO ! End of slab
	316	! ! ===============
	317
	318	!,,,,,,,,,,,,,,,,,,,,,,,,,,,,,synchro,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
	319
	320	! ! ************ ! ! ===============
	321	DO jj = 2, jpjm1 ! Second step ! ! Horizontal slab
	322	! ! ************ ! ! ===============
	323
	324	! II.1 horizontal (pu,pv) gradients
	325	! ---------------------------------
	326
	327	DO jk = 1, jpk
	328	DO ji = 2, jpi
	329	! i-gradient of u at jj
	330	zdiu (ji,jk) = tmask(ji,jj ,jk) * ( pu(ji,jj ,jk) - pu(ji-1,jj ,jk) )
	331	! j-gradient of u and v at jj
	332	zdju (ji,jk) = fmask(ji,jj ,jk) * ( pu(ji,jj+1,jk) - pu(ji ,jj ,jk) )
	333	zdjv (ji,jk) = tmask(ji,jj ,jk) * ( pv(ji,jj ,jk) - pv(ji ,jj-1,jk) )
	334	! j-gradient of u and v at jj+1
	335	zdj1u(ji,jk) = fmask(ji,jj-1,jk) * ( pu(ji,jj ,jk) - pu(ji ,jj-1,jk) )
	336	zdj1v(ji,jk) = tmask(ji,jj+1,jk) * ( pv(ji,jj+1,jk) - pv(ji ,jj ,jk) )
	337	END DO
	338	END DO
	339	DO jk = 1, jpk
	340	DO ji = 1, jpim1
	341	! i-gradient of v at jj
	342	zdiv (ji,jk) = fmask(ji,jj ,jk) * ( pv(ji+1,jj,jk) - pv(ji ,jj ,jk) )
	343	END DO
	344	END DO
	345
	346
	347	! II.2 Vertical fluxes
	348	! --------------------
	349
	350	! Surface and bottom vertical fluxes set to zero
	351
	352	zfuw(:, 1 ) = 0.e0
	353	zfvw(:, 1 ) = 0.e0
	354	zfuw(:,jpk) = 0.e0
	355	zfvw(:,jpk) = 0.e0
	356
	357	! interior (2=<jk=<jpk-1) on pu field
	358
	359	DO jk = 2, jpkm1
	360	DO ji = 2, jpim1
	361	! i- and j-slopes at uw-point
	362	zuwslpi = 0.5 * ( wslpi(ji+1,jj,jk) + wslpi(ji,jj,jk) )
	363	zuwslpj = 0.5 * ( wslpj(ji+1,jj,jk) + wslpj(ji,jj,jk) )
	364	! coef. for the vertical dirative
	365	zcoef0 = e1u(ji,jj) * e2u(ji,jj) / fse3u(ji,jj,jk) &
	366	* ( zuwslpi * zuwslpi + zuwslpj * zuwslpj )
	367	! weights for the i-k, j-k averaging at t- and f-points, resp.
	368	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji+1,jj,jk-1) &
	369	+ tmask(ji,jj,jk )+tmask(ji+1,jj,jk ), 1. )
	370	zmkf = 1./MAX( fmask(ji,jj-1,jk-1)+fmask(ji,jj,jk-1) &
	371	+ fmask(ji,jj-1,jk )+fmask(ji,jj,jk ), 1. )
	372	! coef. for the horitontal derivative
	373	zcoef3 = - e2u(ji,jj) * zmkt * zuwslpi
	374	zcoef4 = - e1u(ji,jj) * zmkf * zuwslpj
	375	! vertical flux on u field
	376	zfuw(ji,jk) = umask(ji,jj,jk) * &
	377	( zcoef0 * ( pu (ji,jj,jk-1) - pu (ji,jj,jk) ) &
	378	+ zcoef3 * ( zdiu (ji,jk-1) + zdiu (ji+1,jk-1) &
	379	+zdiu (ji,jk ) + zdiu (ji+1,jk ) ) &
	380	+ zcoef4 * ( zdj1u(ji,jk-1) + zdju (ji ,jk-1) &
	381	+zdj1u(ji,jk ) + zdju (ji ,jk ) ) )
	382	END DO
	383	END DO
	384
	385	! interior (2=<jk=<jpk-1) on pv field
	386
	387	DO jk = 2, jpkm1
	388	DO ji = 2, jpim1
	389	! i- and j-slopes at vw-point
	390	zvwslpi = 0.5 * ( wslpi(ji,jj+1,jk) + wslpi(ji,jj,jk) )
	391	zvwslpj = 0.5 * ( wslpj(ji,jj+1,jk) + wslpj(ji,jj,jk) )
	392	! coef. for the vertical derivative
	393	zcoef0 = e1v(ji,jj) * e2v(ji,jj) / fse3v(ji,jj,jk) &
	394	* ( zvwslpi * zvwslpi + zvwslpj * zvwslpj )
	395	! weights for the i-k, j-k averaging at f- and t-points, resp.
	396	zmkf = 1./MAX( fmask(ji-1,jj,jk-1)+fmask(ji,jj,jk-1) &
	397	+ fmask(ji-1,jj,jk )+fmask(ji,jj,jk ), 1. )
	398	zmkt = 1./MAX( tmask(ji,jj,jk-1)+tmask(ji,jj+1,jk-1) &
	399	+ tmask(ji,jj,jk )+tmask(ji,jj+1,jk ), 1. )
	400	! coef. for the horizontal derivatives
	401	zcoef3 = - e2v(ji,jj) * zmkf * zvwslpi
	402	zcoef4 = - e1v(ji,jj) * zmkt * zvwslpj
	403	! vertical flux on pv field
	404	zfvw(ji,jk) = vmask(ji,jj,jk) * &
	405	( zcoef0 * ( pv (ji,jj,jk-1) - pv (ji,jj,jk) ) &
	406	+ zcoef3 * ( zdiv (ji,jk-1) + zdiv (ji-1,jk-1) &
	407	+zdiv (ji,jk ) + zdiv (ji-1,jk ) ) &
	408	+ zcoef4 * ( zdjv (ji,jk-1) + zdj1v(ji ,jk-1) &
	409	+zdjv (ji,jk ) + zdj1v(ji ,jk ) ) )
	410	END DO
	411	END DO
	412
	413
	414	! II.3 Divergence of vertical fluxes added to the horizontal divergence
	415	! ---------------------------------------------------------------------
	416
	417	IF( kahm == 1 ) THEN
	418	! multiply the laplacian by the eddy viscosity coefficient
	419	DO jk = 1, jpkm1
	420	DO ji = 2, jpim1
	421	! eddy coef. divided by the volume element
	422	zbur = fsahmu(ji,jj,jk) / ( e1u(ji,jj)e2u(ji,jj)fse3u(ji,jj,jk) )
	423	zbvr = fsahmv(ji,jj,jk) / ( e1v(ji,jj)e2v(ji,jj)fse3v(ji,jj,jk) )
	424	! vertical divergence
	425	zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
	426	zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
	427	! harmonic operator applied to (pu,pv) and multiply by ahm
	428	plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
	429	plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
	430	END DO
	431	END DO
	432	ELSEIF( kahm == 2 ) THEN
	433	! second call, no multiplication
	434	DO jk = 1, jpkm1
	435	DO ji = 2, jpim1
	436	! inverse of the volume element
	437	zbur = 1. / ( e1u(ji,jj)e2u(ji,jj)fse3u(ji,jj,jk) )
	438	zbvr = 1. / ( e1v(ji,jj)e2v(ji,jj)fse3v(ji,jj,jk) )
	439	! vertical divergence
	440	zuav = zfuw(ji,jk) - zfuw(ji,jk+1)
	441	zvav = zfvw(ji,jk) - zfvw(ji,jk+1)
	442	! harmonic operator applied to (pu,pv)
	443	plu(ji,jj,jk) = ( plu(ji,jj,jk) + zuav ) * zbur
	444	plv(ji,jj,jk) = ( plv(ji,jj,jk) + zvav ) * zbvr
	445	END DO
	446	END DO
	447	ELSE
	448	IF(lwp)WRITE(numout,*) ' ldfguv: kahm= 1 or 2, here =', kahm
	449	IF(lwp)WRITE(numout,*) ' We stop'
	450	STOP 'ldfguv'
	451	ENDIF
	452	! ! ===============
	453	END DO ! End of slab
	454	! ! ===============
[2715]	455
[3294]	456	CALL wrk_dealloc( jpi, jpj, ziut, zjuf, zjvt, zivf, zdku, zdk1u, zdkv, zdk1v )
[2715]	457	!
[3294]	458	IF( nn_timing == 1 ) CALL timing_stop('ldfguv')
	459	!
[3]	460	END SUBROUTINE ldfguv
	461
	462	#else
	463	!!----------------------------------------------------------------------
	464	!! Dummy module : NO rotation of mixing tensor
	465	!!----------------------------------------------------------------------
	466	CONTAINS
	467	SUBROUTINE dyn_ldf_bilapg( kt ) ! Dummy routine
[2715]	468	INTEGER, INTENT(in) :: kt
[32]	469	WRITE(,) 'dyn_ldf_bilapg: You should not have seen this print! error?', kt
[3]	470	END SUBROUTINE dyn_ldf_bilapg
	471	#endif
	472
	473	!!======================================================================
	474	END MODULE dynldf_bilapg

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

Download in other formats: