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.

trabbl_adv.h90 in tags/nemo_v1_13_dev5/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: tags/nemo_v1_13_dev5/NEMO/OPA_SRC/TRA/trabbl_adv.h90 @ 5302

Last change on this file since 5302 was 481, checked in by opalod, 18 years ago
nemo_v1_bugfix_047 : CT : correction in order to take into account the partial steps for advective bbl
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 20.4 KB

Line
1	!!----------------------------------------------------------------------
2	!! * trabbl_adv.h90 *
3	!!----------------------------------------------------------------------
4
5	!!----------------------------------------------------------------------
6	!! OPA 9.0 , LOCEAN-IPSL (2005)
7	!! $Header$
8	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
9	!!----------------------------------------------------------------------
10
11	SUBROUTINE tra_bbl_adv( kt )
12	!!----------------------------------------------------------------------
13	!! * ROUTINE tra_bbl_adv *
14	!!
15	!! ** Purpose : Compute the before tracer (t & s) trend associated
16	!! with the bottom boundary layer and add it to the general trend
17	!! of tracer equations. The bottom boundary layer is supposed to be
18	!! both an advective and diffusive bottom boundary layer.
19	!!
20	!! ** Method : Computes the bottom boundary horizontal and vertical
21	!! advection terms. Add it to the general trend : ta =ta + adv_bbl.
22	!! When the product grad( rho) * grad(h) < 0 (where grad is a
23	!! along bottom slope gradient) an additional lateral 2nd order
24	!! diffusion along the bottom slope is added to the general
25	!! tracer trend, otherwise the additional trend is set to 0.
26	!! Second order operator (laplacian type) with variable coefficient
27	!! computed as follow for temperature (idem on s):
28	!! difft = 1/(e1te2te3t) { di-1[ ahbt e2u*e3u/e1u di[ztb] ]
29	!! + dj-1[ ahbt e1v*e3v/e2v dj[ztb] ] }
30	!! where ztb is a 2D array: the bottom ocean te;perature and ahtb
31	!! is a time and space varying diffusive coefficient defined by:
32	!! ahbt = zahbp if grad(rho).grad(h) < 0
33	!! = 0. otherwise.
34	!! Note that grad(.) is the along bottom slope gradient. grad(rho)
35	!! is evaluated using the local density (i.e. referenced at the
36	!! local depth). Typical value of ahbt is 2000 m2/s (equivalent to
37	!! a downslope velocity of 20 cm/s if the condition for slope
38	!! convection is satified)
39	!! Add this before trend to the general trend (ta,sa) of the
40	!! botton ocean tracer point:
41	!! ta = ta + difft
42	!!
43	!! ** Action : - update (ta,sa) at the bottom level with the bottom
44	!! boundary layer trend
45	!! - save the lateral diffusion trends in tldfbbl/sldfbbl ('key_trdtra')
46	!! - save the horizontal advection trends in tladbbl/sladbbl ('key_trdtra')
47	!!
48	!! References :
49	!! Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
50	!!
51	!! History :
52	!! 8.5 ! 02-12 (A. de Miranda, G. Madec) Original Code
53	!! 9.0 ! 04-01 (A. de Miranda, G. Madec, J.M. Molines )
54	!! 9.0 ! 04-08 (C. Talandier) New trends organization
55	!!----------------------------------------------------------------------
56	!! * Modules used
57	USE eosbn2 ! equation of state
58	USE oce, ONLY : ztdta => ua, & ! use ua as 3D workspace
59	ztdsa => va ! use va as 3D workspace
60
61	!! * Arguments
62	INTEGER, INTENT( in ) :: kt ! ocean time-step
63
64	!! * Local declarations
65	INTEGER :: ji, jj, jk ! dummy loop indices
66	INTEGER :: ik, iku, ikv ! temporary integers
67
68	REAL(wp) :: &
69	zsign, zt, zs, zh, zalbet, & ! temporary scalars
70	zgdrho, zbtr, zta, zsa ! " "
71	REAL(wp), DIMENSION(jpi,jpj) :: &
72	ztnb, zsnb, zdep, ztbb, zsbb
73	REAL(wp), DIMENSION(jpi,jpj) :: & ! temporary workspace arrays
74	zalphax, zwu, zunb, & ! " "
75	zalphay, zwv, zvnb, & ! " "
76	zwx, zwy, zww, zwz ! " "
77	REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
78	zhdivn ! temporary workspace arrays
79	REAL(wp) :: &
80	zfui, zfvj, zbt, zsigna, & ! temporary scalars
81	iku1, iku2, ikv1, ikv2, & ! temporary scalars
82	ze3u,ze3v
83	REAL(wp) :: &
84	fsalbt, pft, pfs, pfh ! statement function
85	!!----------------------------------------------------------------------
86	! ratio alpha/beta
87	! ================
88	! fsalbt: ratio of thermal over saline expension coefficients
89	! pft : potential temperature in degrees celcius
90	! pfs : salinity anomaly (s-35) in psu
91	! pfh : depth in meters
92
93	fsalbt( pft, pfs, pfh ) = &
94	( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft &
95	- 0.203814e-03 ) * pft &
96	+ 0.170907e-01 ) * pft &
97	+ 0.665157e-01 &
98	+(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs &
99	+ ( ( - 0.302285e-13 * pfh &
100	- 0.251520e-11 * pfs &
101	+ 0.512857e-12 * pft * pft ) * pfh &
102	- 0.164759e-06 * pfs &
103	+( 0.791325e-08 * pft - 0.933746e-06 ) * pft &
104	+ 0.380374e-04 ) * pfh
105	!!----------------------------------------------------------------------
106
107
108	IF( kt == nit000 ) CALL tra_bbl_init ! initialization at first time-step
109
110	! Save ta and sa trends
111	IF( l_trdtra ) THEN
112	ztdta(:,:,:) = ta(:,:,:)
113	ztdsa(:,:,:) = sa(:,:,:)
114	ENDIF
115
116	! 1. 2D fields of bottom temperature and salinity, and bottom slope
117	! -----------------------------------------------------------------
118	! mbathy= number of w-level, minimum value=1 (cf dommsk.F)
119
120	#if defined key_vectopt_loop && ! defined key_mpp_omp
121	jj = 1
122	DO ji = 1, jpij ! vector opt. (forced unrolling)
123	#else
124	DO jj = 1, jpj
125	DO ji = 1, jpi
126	#endif
127	ik = mbkt(ji,jj) ! index of the bottom ocean T-level
128	ztnb(ji,jj) = tn(ji,jj,ik)
129	zsnb(ji,jj) = sn(ji,jj,ik)
130	ztbb(ji,jj) = tb(ji,jj,ik)
131	zsbb(ji,jj) = sb(ji,jj,ik)
132	zdep(ji,jj) = fsdept(ji,jj,ik) ! depth of the ocean bottom T-level
133
134	zunb(ji,jj) = un(ji,jj,mbku(ji,jj))
135	zvnb(ji,jj) = vn(ji,jj,mbkv(ji,jj))
136	#if ! defined key_vectopt_loop \|\| defined key_mpp_omp
137	END DO
138	#endif
139	END DO
140
141
142	! 2. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0
143	! --------------------------------------------
144	! Sign of the local density gradient along the i- and j-slopes
145	! multiplied by the slope of the ocean bottom
146
147	SELECT CASE ( neos )
148
149	CASE ( 0 ) ! Jackett and McDougall (1994) formulation
150
151	DO jj = 1, jpjm1
152	DO ji = 1, fs_jpim1 ! vector opt.
153	! ... temperature, salinity anomalie and depth
154	zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
155	zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
156	zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
157	! ... masked ratio alpha/beta
158	zalbet = fsalbt( zt, zs, zh ) * umask(ji,jj,1)
159	! ... local density gradient along i-bathymetric slope
160	zgdrho = zalbet * ( ztnb(ji+1,jj) - ztnb(ji,jj) ) &
161	& - ( zsnb(ji+1,jj) - zsnb(ji,jj) )
162	zgdrho = zgdrho * umask(ji,jj,1)
163	! ... sign of local i-gradient of density multiplied by the i-slope
164	zsign = SIGN( 0.5, -zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
165	zsigna= SIGN( 0.5, zunb(ji,jj) * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
166	zalphax(ji,jj)=( 0.5 + zsigna ) * ( 0.5 - zsign ) * umask(ji,jj,1)
167	END DO
168	END DO
169
170	DO jj = 1, jpjm1
171	DO ji = 1, fs_jpim1 ! vector opt.
172	! ... temperature, salinity anomalie and depth
173	zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
174	zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
175	zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
176	! ... masked ratio alpha/beta
177	zalbet = fsalbt( zt, zs, zh ) * vmask(ji,jj,1)
178	! ... local density gradient along j-bathymetric slope
179	zgdrho = zalbet * ( ztnb(ji,jj+1) - ztnb(ji,jj) ) &
180	& - ( zsnb(ji,jj+1) - zsnb(ji,jj) )
181	zgdrho = zgdrho*vmask(ji,jj,1)
182	! ... sign of local j-gradient of density multiplied by the j-slope
183	zsign = SIGN( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
184	zsigna= SIGN( 0.5, zvnb(ji,jj) * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
185	zalphay(ji,jj)=( 0.5 + zsigna ) * ( 0.5 - zsign ) * vmask(ji,jj,1)
186	END DO
187	END DO
188
189
190	CASE ( 1 ) ! Linear formulation function of temperature only
191	!
192	DO jj = 1, jpjm1
193	DO ji = 1, fs_jpim1 ! vector opt.
194	! local 'density/temperature' gradient along i-bathymetric slope
195	zgdrho = ( ztnb(ji+1,jj) - ztnb(ji,jj) )
196	! sign of local i-gradient of density multiplied by the i-slope
197	zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
198	zsigna= SIGN( 0.5, zunb(ji,jj) * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
199	zalphax(ji,jj)=( 0.5 + zsigna ) * ( 0.5 - zsign ) * umask(ji,jj,1)
200
201	! local density gradient along j-bathymetric slope
202	zgdrho = ( ztnb(ji,jj+1) - ztnb(ji,jj) )
203	! sign of local j-gradient of density multiplied by the j-slope
204	zsign = SIGN( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
205	zsigna= SIGN( 0.5, zvnb(ji,jj) * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
206	zalphay(ji,jj)=( 0.5 + zsigna ) * ( 0.5 - zsign ) * vmask(ji,jj,1)
207	END DO
208	END DO
209
210	CASE ( 2 ) ! Linear formulation function of temperature and salinity
211	DO jj = 1, jpjm1
212	DO ji = 1, fs_jpim1 ! vector opt.
213	! local density gradient along i-bathymetric slope
214	zgdrho = - ( rbeta*( zsnb(ji+1,jj) - zsnb(ji,jj) ) &
215	& - ralpha*( ztnb(ji+1,jj) - ztnb(ji,jj) ) )
216	! sign of local i-gradient of density multiplied by the i-slope
217	zsign = SIGN( 0.5, - zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
218	zsigna= SIGN( 0.5, zunb(ji,jj) * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
219	zalphax(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * umask(ji,jj,1)
220
221	! local density gradient along j-bathymetric slope
222	zgdrho = - ( rbeta*( zsnb(ji,jj+1) - zsnb(ji,jj) ) &
223	- ralpha*( ztnb(ji,jj+1) - ztnb(ji,jj) ) )
224	! sign of local j-gradient of density multiplied by the j-slope
225	zsign = SIGN( 0.5, - zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
226	zsigna= SIGN( 0.5, zvnb(ji,jj) * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
227	zalphay(ji,jj) = ( 0.5 + zsigna ) * ( 0.5 - zsign ) * vmask(ji,jj,1)
228	END DO
229	END DO
230
231	CASE DEFAULT
232
233	WRITE(ctmp1,*) ' bad flag value for neos = ', neos
234	CALL ctl_stop( ctmp1 )
235
236	END SELECT
237
238	! lateral boundary conditions on zalphax and zalphay (unchanged sign)
239	CALL lbc_lnk( zalphax, 'U', 1. ) ; CALL lbc_lnk( zalphay, 'V', 1. )
240
241
242	! 3. Velocities that are exchanged between ajacent bottom boxes.
243	!---------------------------------------------------------------
244
245	! ... is equal to zero but where bbl will work.
246	u_bbl(:,:,:) = 0.e0
247	v_bbl(:,:,:) = 0.e0
248
249	IF( ln_zps ) THEN ! partial steps correction
250
251	# if defined key_vectopt_loop && ! defined key_mpp_omp
252	jj = 1
253	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
254	# else
255	DO jj = 1, jpjm1
256	DO ji = 1, jpim1
257	# endif
258	iku = mbku(ji ,jj )
259	ikv = mbkv(ji ,jj )
260	iku1 = mbkt(ji+1,jj )
261	iku2 = mbkt(ji ,jj )
262	ikv1 = mbkt(ji ,jj+1)
263	ikv2 = mbkt(ji ,jj )
264	ze3u = MIN( fse3u(ji,jj,iku1), fse3u(ji,jj,iku2) )
265	ze3v = MIN( fse3v(ji,jj,ikv1), fse3v(ji,jj,ikv2) )
266
267	IF( MAX(iku,ikv) > 1 ) THEN
268	u_bbl(ji,jj,iku) = zalphax(ji,jj) * un(ji,jj,iku) * ze3u / fse3u(ji,jj,iku)
269	v_bbl(ji,jj,ikv) = zalphay(ji,jj) * vn(ji,jj,ikv) * ze3v / fse3v(ji,jj,ikv)
270	ENDIF
271	# if ! defined key_vectopt_loop \|\| defined key_mpp_omp
272	END DO
273	# endif
274	END DO
275
276	! lateral boundary conditions on u_bbl and v_bbl (changed sign)
277	CALL lbc_lnk( u_bbl, 'U', -1. ) ; CALL lbc_lnk( v_bbl, 'V', -1. )
278
279	ELSE ! if not partial step loop over the whole domain no lbc call
280
281	#if defined key_vectopt_loop && ! defined key_mpp_omp
282	jj = 1
283	DO ji = 1, jpij ! vector opt. (forced unrolling)
284	#else
285	DO jj = 1, jpj
286	DO ji = 1, jpi
287	#endif
288	iku = mbku(ji,jj)
289	ikv = mbkv(ji,jj)
290	IF( MAX(iku,ikv) > 1 ) THEN
291	u_bbl(ji,jj,iku) = zalphax(ji,jj) * un(ji,jj,iku)
292	v_bbl(ji,jj,ikv) = zalphay(ji,jj) * vn(ji,jj,ikv)
293	ENDIF
294	#if ! defined key_vectopt_loop \|\| defined key_mpp_omp
295	END DO
296	#endif
297	END DO
298
299	ENDIF
300
301
302	! 5. Along sigma advective trend
303	! -------------------------------
304	! ... Second order centered tracer flux at u and v-points
305
306	# if defined key_vectopt_loop && ! defined key_mpp_omp
307	jj = 1
308	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
309	# else
310	DO jj = 1, jpjm1
311	DO ji = 1, jpim1
312	# endif
313	iku = mbku(ji,jj)
314	ikv = mbkv(ji,jj)
315	zfui = e2u(ji,jj) * fse3u(ji,jj,iku) * u_bbl(ji,jj,iku)
316	zfvj = e1v(ji,jj) * fse3v(ji,jj,ikv) * v_bbl(ji,jj,ikv)
317	! centered scheme
318	! zwx(ji,jj) = 0.5* zfui * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
319	! zwy(ji,jj) = 0.5* zfvj * ( ztnb(ji,jj) + ztnb(ji,jj+1) )
320	! zww(ji,jj) = 0.5* zfui * ( zsnb(ji,jj) + zsnb(ji+1,jj) )
321	! zwz(ji,jj) = 0.5* zfvj * ( zsnb(ji,jj) + zsnb(ji,jj+1) )
322	! upstream scheme
323	zwx(ji,jj) = ( ( zfui + ABS( zfui ) ) * ztbb(ji ,jj ) &
324	& +( zfui - ABS( zfui ) ) * ztbb(ji+1,jj ) ) * 0.5
325	zwy(ji,jj) = ( ( zfvj + ABS( zfvj ) ) * ztbb(ji ,jj ) &
326	& +( zfvj - ABS( zfvj ) ) * ztbb(ji ,jj+1) ) * 0.5
327	zww(ji,jj) = ( ( zfui + ABS( zfui ) ) * zsbb(ji ,jj ) &
328	& +( zfui - ABS( zfui ) ) * zsbb(ji+1,jj ) ) * 0.5
329	zwz(ji,jj) = ( ( zfvj + ABS( zfvj ) ) * zsbb(ji ,jj ) &
330	& +( zfvj - ABS( zfvj ) ) * zsbb(ji ,jj+1) ) * 0.5
331	#if ! defined key_vectopt_loop \|\| defined key_mpp_omp
332	END DO
333	#endif
334	END DO
335	# if defined key_vectopt_loop && ! defined key_mpp_omp
336	jj = 1
337	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
338	# else
339	DO jj = 2, jpjm1
340	DO ji = 2, jpim1
341	# endif
342	ik = mbkt(ji,jj)
343	zbtr = 1. / ( e1t(ji,jj)e2t(ji,jj)fse3t(ji,jj,ik) )
344	! horizontal advective trends
345	zta = - zbtr * ( zwx(ji,jj) - zwx(ji-1,jj ) &
346	& + zwy(ji,jj) - zwy(ji ,jj-1) )
347	zsa = - zbtr * ( zww(ji,jj) - zww(ji-1,jj ) &
348	& + zwz(ji,jj) - zwz(ji ,jj-1) )
349
350	! add it to the general tracer trends
351	ta(ji,jj,ik) = ta(ji,jj,ik) + zta
352	sa(ji,jj,ik) = sa(ji,jj,ik) + zsa
353	#if ! defined key_vectopt_loop \|\| defined key_mpp_omp
354	END DO
355	#endif
356	END DO
357
358	! save the trends for diagnostic
359	! BBL lateral advection tracers trends
360	IF( l_trdtra ) THEN
361	# if defined key_vectopt_loop && ! defined key_mpp_omp
362	jj = 1
363	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
364	# else
365	DO jj = 2, jpjm1
366	DO ji = 2, jpim1
367	# endif
368	ik = mbkt(ji,jj)
369	tladbbl(ji,jj) = ta(ji,jj,ik) - ztdta(ji,jj,ik)
370	sladbbl(ji,jj) = sa(ji,jj,ik) - ztdsa(ji,jj,ik)
371	# if ! defined key_vectopt_loop \|\| defined key_mpp_omp
372	END DO
373	# endif
374	END DO
375
376	ENDIF
377
378	IF(ln_ctl) THEN
379	CALL prt_ctl(tab3d_1=ta, clinfo1=' bbl - Ta: ', mask1=tmask, &
380	& tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
381	ENDIF
382
383
384	! 6. Vertical advection velocities
385	! --------------------------------
386	! ... computes divergence perturbation (velocties to be removed from upper t boxes :
387	DO jk= 1, jpkm1
388	DO jj=1, jpjm1
389	DO ji = 1, fs_jpim1 ! vector opt.
390	zwu(ji,jj) = -e2u(ji,jj) * u_bbl(ji,jj,jk) * fse3u(ji,jj,jk)
391	zwv(ji,jj) = -e1v(ji,jj) * v_bbl(ji,jj,jk) * fse3v(ji,jj,jk)
392	END DO
393	END DO
394
395	! ... horizontal divergence
396	DO jj = 2, jpjm1
397	DO ji = fs_2, fs_jpim1 ! vector opt.
398	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk)
399	zhdivn(ji,jj,jk) = ( zwu(ji,jj) - zwu(ji-1,jj ) &
400	+ zwv(ji,jj) - zwv(ji ,jj-1) ) / zbt
401	END DO
402	END DO
403	END DO
404
405
406	! ... horizontal bottom divergence
407
408	IF( ln_zps ) THEN
409
410	# if defined key_vectopt_loop && ! defined key_mpp_omp
411	jj = 1
412	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
413	# else
414	DO jj = 1, jpjm1
415	DO ji = 1, jpim1
416	# endif
417	iku = mbku(ji ,jj )
418	ikv = mbkv(ji ,jj )
419	iku1 = mbkt(ji+1,jj )
420	iku2 = mbkt(ji ,jj )
421	ikv1 = mbkt(ji ,jj+1)
422	ikv2 = mbkt(ji ,jj )
423	ze3u = MIN( fse3u(ji,jj,iku1), fse3u(ji,jj,iku2) )
424	ze3v = MIN( fse3v(ji,jj,ikv1), fse3v(ji,jj,ikv2) )
425
426	zwu(ji,jj) = zalphax(ji,jj) * e2u(ji,jj) * ze3u
427	zwv(ji,jj) = zalphay(ji,jj) * e1v(ji,jj) * ze3v
428	#if ! defined key_vectopt_loop \|\| defined key_mpp_omp
429	END DO
430	#endif
431	END DO
432
433	ELSE
434
435	# if defined key_vectopt_loop && ! defined key_mpp_omp
436	jj = 1
437	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
438	# else
439	DO jj = 1, jpjm1
440	DO ji = 1, jpim1
441	# endif
442	iku = mbku(ji,jj)
443	ikv = mbkv(ji,jj)
444	zwu(ji,jj) = zalphax(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,iku)
445	zwv(ji,jj) = zalphay(ji,jj) * e1v(ji,jj) * fse3v(ji,jj,ikv)
446	#if ! defined key_vectopt_loop \|\| defined key_mpp_omp
447	END DO
448	#endif
449	END DO
450
451	ENDIF
452
453
454	# if defined key_vectopt_loop && ! defined key_mpp_omp
455	jj = 1
456	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
457	# else
458	DO jj = 2, jpjm1
459	DO ji = 2, jpim1
460	# endif
461	ik = mbkt(ji,jj)
462	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik)
463	zhdivn(ji,jj,ik) = &
464	& ( zwu(ji ,jj ) * ( zunb(ji ,jj ) - un(ji ,jj ,ik) ) &
465	& - zwu(ji-1,jj ) * ( zunb(ji-1,jj ) - un(ji-1,jj ,ik) ) &
466	& + zwv(ji ,jj ) * ( zvnb(ji ,jj ) - vn(ji ,jj ,ik) ) &
467	& - zwv(ji ,jj-1) * ( zvnb(ji ,jj-1) - vn(ji ,jj-1,ik) ) &
468	& ) / zbt
469
470	# if ! defined key_vectopt_loop \|\| defined key_mpp_omp
471	END DO
472	# endif
473	END DO
474
475	! 7. compute additional vertical velocity to be used in t boxes
476	! -------------------------------------------------------------
477
478	! ... Computation from the bottom
479	! Note that w_bbl(:,:,jpk) has been set to 0 in tra_bbl_init
480	DO jk = jpkm1, 1, -1
481	DO jj= 2, jpjm1
482	DO ji = fs_2, fs_jpim1 ! vector opt.
483	w_bbl(ji,jj,jk) = w_bbl(ji,jj,jk+1) - fse3t(ji,jj,jk)*zhdivn(ji,jj,jk)
484	END DO
485	END DO
486	END DO
487
488	! Boundary condition on w_bbl (unchanged sign)
489	CALL lbc_lnk( w_bbl, 'W', 1. )
490
491	END SUBROUTINE tra_bbl_adv

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

Download in other formats:

Original Format