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

Context Navigation

source: trunk/NEMO/OPA_SRC/TRA/trabbl_adv.h90 @ 3

Last change on this file since 3 was 3, checked in by opalod, 20 years ago
Initial revision
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 21.2 KB

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

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

Download in other formats:

Original Format