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

Context Navigation

source: tags/nemo_dev_x5/NEMO/OPA_SRC/TRA/trabbl_adv.h90 @ 1792

Last change on this file since 1792 was 106, checked in by opalod, 20 years ago
CT : UPDATE067 : Add control indices nictl, njctl used in SUM function output to compare mono versus multi procs runs
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 20.8 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, zwu, zunb, & ! " "
74	zalphay, zwv, zvnb, & ! " "
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	REAL(wp) :: &
81	fsalbt, pft, pfs, pfh ! statement function
82	!!----------------------------------------------------------------------
83	! ratio alpha/beta
84	! ================
85	! fsalbt: ratio of thermal over saline expension coefficients
86	! pft : potential temperature in degrees celcius
87	! pfs : salinity anomaly (s-35) in psu
88	! pfh : depth in meters
89
90	fsalbt( pft, pfs, pfh ) = &
91	( ( ( -0.255019e-07 * pft + 0.298357e-05 ) * pft &
92	- 0.203814e-03 ) * pft &
93	+ 0.170907e-01 ) * pft &
94	+ 0.665157e-01 &
95	+(-0.678662e-05 * pfs - 0.846960e-04 * pft + 0.378110e-02 ) * pfs &
96	+ ( ( - 0.302285e-13 * pfh &
97	- 0.251520e-11 * pfs &
98	+ 0.512857e-12 * pft * pft ) * pfh &
99	- 0.164759e-06 * pfs &
100	+( 0.791325e-08 * pft - 0.933746e-06 ) * pft &
101	+ 0.380374e-04 ) * pfh
102	!!----------------------------------------------------------------------
103
104
105	IF( kt == nit000 ) CALL tra_bbl_init ! initialization at first time-step
106
107
108	! 1. 2D fields of bottom temperature and salinity, and bottom slope
109	! -----------------------------------------------------------------
110	! mbathy= number of w-level, minimum value=1 (cf dommsk.F)
111
112	#if defined key_vectopt_loop && ! defined key_autotasking
113	jj = 1
114	DO ji = 1, jpij ! vector opt. (forced unrolling)
115	#else
116	DO jj = 1, jpj
117	DO ji = 1, jpi
118	#endif
119	ik = mbkt(ji,jj) ! index of the bottom ocean T-level
120	ztnb(ji,jj) = tn(ji,jj,ik) * tmask(ji,jj,1) ! masked now T at the ocean bottom
121	zsnb(ji,jj) = sn(ji,jj,ik) * tmask(ji,jj,1) ! masked now S at the ocean bottom
122	ztbb(ji,jj) = tb(ji,jj,ik) * tmask(ji,jj,1) ! masked before T at the ocean bottom
123	zsbb(ji,jj) = sb(ji,jj,ik) * tmask(ji,jj,1) ! masked before S at the ocean bottom
124	zdep(ji,jj) = fsdept(ji,jj,ik) ! depth of the ocean bottom T-level
125	#if ! defined key_vectopt_loop \|\| defined key_autotasking
126	END DO
127	#endif
128	END DO
129	#if defined key_vectopt_loop && ! defined key_autotasking
130	jj = 1
131	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
132	zunb(ji,jj) = un(ji,jj,mbku(ji,jj)) * umask(ji,jj,1)
133	zvnb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) * vmask(ji,jj,1) ! retirer le mask en u, v et t !
134	END DO
135	#else
136	DO jj = 1, jpjm1
137	DO ji = 1, jpim1
138	zunb(ji,jj) = un(ji,jj,mbku(ji,jj)) * umask(ji,jj,1)
139	zvnb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) * vmask(ji,jj,1)
140	END DO
141	END DO
142	#endif
143
144	! boundary conditions on zunb and zvnb (changed sign)
145	CALL lbc_lnk( zunb, 'U', -1. ) ; CALL lbc_lnk( zvnb, 'V', -1. )
146	! boundary condition on ztnb and znbb
147	CALL lbc_lnk( ztnb, 'T', 1. ) ; CALL lbc_lnk( ztbb, 'T', 1. )
148	! boundary condition on zsnb and zsbb
149	CALL lbc_lnk( zsnb, 'T', 1. ) ; CALL lbc_lnk( zsbb, 'T', 1. )
150
151	! Conditional diffusion along the slope in the bottom boundary layer
152
153	#if defined key_trabbl_dif
154	# if defined key_vectopt_loop && ! defined key_autotasking
155	jj = 1
156	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
157	# else
158	DO jj = 1, jpjm1
159	DO ji = 1, jpim1
160	# endif
161	iku = mbku(ji,jj)
162	ikv = mbkv(ji,jj)
163	zahu(ji,jj) = atrbble2u(ji,jj)fse3u(ji,jj,iku)/e1u(ji,jj) * umask(ji,jj,1)
164	zahv(ji,jj) = atrbble1v(ji,jj)fse3v(ji,jj,ikv)/e2v(ji,jj) * vmask(ji,jj,1)
165	# if ! defined key_vectopt_loop \|\| defined key_autotasking
166	END DO
167	# endif
168	END DO
169	#endif
170
171
172	! 2. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0
173	! --------------------------------------------
174	! Sign of the local density gradient along the i- and j-slopes
175	! multiplied by the slope of the ocean bottom
176
177	SELECT CASE ( neos )
178
179	CASE ( 0 ) ! Jackett and McDougall (1994) formulation
180
181	DO jj = 1, jpjm1
182	DO ji = 1, fs_jpim1 ! vector opt.
183	! ... temperature, salinity anomalie and depth
184	zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
185	zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
186	zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
187	! ... masked ratio alpha/beta
188	zalbet = fsalbt( zt, zs, zh )*umask(ji,jj,1)
189	! ... local density gradient along i-bathymetric slope
190	zgdrho = zalbet*( ztnb(ji+1,jj) - ztnb(ji,jj) ) &
191	- ( zsnb(ji+1,jj) - zsnb(ji,jj) )
192	zgdrho = zgdrho * umask(ji,jj,1)
193	! ... sign of local i-gradient of density multiplied by the i-slope
194	zsign = sign( 0.5, -zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
195	zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj)
196
197	zsigna= sign(0.5, zunb(ji,jj)*( zdep(ji+1,jj) - zdep(ji,jj) ))
198	zalphax(ji,jj)=(0.5+zsigna)(0.5-zsign)umask(ji,jj,1)
199	END DO
200	END DO
201
202	DO jj = 1, jpjm1
203	DO ji = 1, fs_jpim1 ! vector opt.
204	! ... temperature, salinity anomalie and depth
205	zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
206	zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
207	zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
208	! ... masked ratio alpha/beta
209	zalbet = fsalbt( zt, zs, zh )*vmask(ji,jj,1)
210	! ... local density gradient along j-bathymetric slope
211	zgdrho = zalbet*( ztnb(ji,jj+1) - ztnb(ji,jj) ) &
212	- ( zsnb(ji,jj+1) - zsnb(ji,jj) )
213	zgdrho = zgdrho*vmask(ji,jj,1)
214	! ... sign of local j-gradient of density multiplied by the j-slope
215	zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
216	zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj)
217
218	zsigna= sign(0.5, zvnb(ji,jj)*(zdep(ji,jj+1) - zdep(ji,jj) ) )
219	zalphay(ji,jj)=(0.5+zsigna)(0.5-zsign)vmask(ji,jj,1)
220	END DO
221	END DO
222
223
224	CASE ( 1 ) ! Linear formulation function of temperature only
225
226	IF(lwp) WRITE(numout,cform_err)
227	IF(lwp) WRITE(numout,) ' use of linear eos rho(T,S) = rau0 ( rbeta * S - ralpha * T )'
228	IF(lwp) WRITE(numout,*) ' bbl not implented: easy to do it '
229	nstop = nstop + 1
230
231	CASE ( 2 ) ! Linear formulation function of temperature and salinity
232
233	IF(lwp) WRITE(numout,cform_err)
234	IF(lwp) WRITE(numout,) ' use of linear eos rho(T,S) = rau0 ( rbeta * S - ralpha * T )'
235	IF(lwp) WRITE(numout,*) ' bbl not implented: easy to do it '
236	nstop = nstop + 1
237
238	CASE DEFAULT
239
240	IF(lwp) WRITE(numout,cform_err)
241	IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
242	nstop = nstop + 1
243
244	END SELECT
245
246	! lateral boundary conditions on zalphax and zalphay (unchanged sign)
247	CALL lbc_lnk( zalphax, 'U', 1. ) ; CALL lbc_lnk( zalphay, 'V', 1. )
248
249
250	! 3. Velocities that are exchanged between ajacent bottom boxes.
251	!---------------------------------------------------------------
252
253	! ... is equal to zero but where bbl will work.
254	u_bbl(:,:,:) = 0.e0
255	v_bbl(:,:,:) = 0.e0
256	# if defined key_vectopt_loop && ! defined key_autotasking
257	jj = 1
258	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
259	# else
260	DO jj = 1, jpjm1
261	DO ji = 1, jpim1
262	# endif
263	iku = mbku(ji,jj)
264	ikv = mbkv(ji,jj)
265	IF( MAX(iku,ikv) > 1 ) THEN
266	u_bbl(ji,jj,iku) = zalphax(ji,jj) * un(ji,jj,iku) * umask(ji,jj,1)
267	v_bbl(ji,jj,ikv) = zalphay(ji,jj) * vn(ji,jj,ikv) * vmask(ji,jj,1)
268	ENDIF
269	# if ! defined key_vectopt_loop \|\| defined key_autotasking
270	END DO
271	# endif
272	END DO
273
274	! lateral boundary conditions on u_bbl and v_bbl (changed sign)
275	CALL lbc_lnk( u_bbl, 'U', -1. ) ; CALL lbc_lnk( v_bbl, 'V', -1. )
276
277
278
279	#if defined key_trabbl_dif
280	! 4. Additional second order diffusive trends
281	! -------------------------------------------
282
283	! ... first derivative (gradient)
284	DO jj = 1, jpjm1
285	DO ji = 1, fs_jpim1 ! vertor opt.
286
287	zkx(ji,jj) = zki(ji,jj)*( ztbb(ji+1,jj) - ztbb(ji,jj) )
288	zkz(ji,jj) = zki(ji,jj)*( zsbb(ji+1,jj) - zsbb(ji,jj) )
289
290	zky(ji,jj) = zkj(ji,jj)*( ztbb(ji,jj+1) - ztbb(ji,jj) )
291	zkw(ji,jj) = zkj(ji,jj)*( zsbb(ji,jj+1) - zsbb(ji,jj) )
292	END DO
293	END DO
294
295	IF( cp_cfg == "orca" ) THEN
296	SELECT CASE ( jp_cfg )
297	! ! =======================
298	CASE ( 2 ) ! ORCA_R2 configuration
299	! ! =======================
300	! Gibraltar enhancement of BBL
301	zkx( mi0(139):mi1(140) , mj0(102):mj1(102) ) = 4.e0 * zkx( mi0(139):mi1(140) , mj0(102):mj1(102) )
302	zky( mi0(139):mi1(140) , mj0(102):mj1(102) ) = 4.e0 * zky( mi0(139):mi1(140) , mj0(102):mj1(102) )
303
304	! Red Sea enhancement of BBL
305	zkx( mi0(161):mi1(162) , mj0(88):mj1(88) ) = 10.e0 * zkx( mi0(161):mi1(162) , mj0(88):mj1(88) )
306	zky( mi0(161):mi1(162) , mj0(88):mj1(88) ) = 10.e0 * zky( mi0(161):mi1(162) , mj0(88):mj1(88) )
307
308	! ! =======================
309	CASE ( 4 ) ! ORCA_R4 configuration
310	! ! =======================
311	! Gibraltar enhancement of BBL
312	zkx( mi0(70):mi1(71) , mj0(52):mj1(52) ) = 4.e0 * zkx( mi0(70):mi1(71) , mj0(52):mj1(52) )
313	zky( mi0(70):mi1(71) , mj0(52):mj1(52) ) = 4.e0 * zky( mi0(70):mi1(71) , mj0(52):mj1(52) )
314
315	END SELECT
316
317	ENDIF
318
319	! ... second derivative (divergence) and add to the general tracer trend
320
321	# if defined key_vectopt_loop && ! defined key_autotasking
322	jj = 1
323	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
324	# else
325	DO jj = 2, jpjm1
326	DO ji = 2, jpim1
327	# endif
328	ik = mbkt(ji,jj)
329	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik) )
330	zta = ( zkx(ji,jj) - zkx(ji-1,jj ) &
331	& + zky(ji,jj) - zky(ji ,jj-1) ) * zbtr
332	zsa = ( zkz(ji,jj) - zkz(ji-1,jj ) &
333	& + zkw(ji,jj) - zkw(ji ,jj-1) ) * zbtr
334	ta(ji,jj,ik) = ta(ji,jj,ik) + zta
335	sa(ji,jj,ik) = sa(ji,jj,ik) + zsa
336	#if defined key_trdtra \|\| defined key_trdmld
337	bbltrd(ji,jj,1) = + zta
338	bbltrd(ji,jj,2) = + zsa
339	ttrd(ji,jj,ik,3) = ttrd(ji,jj,ik,3) + zta
340	strd(ji,jj,ik,3) = strd(ji,jj,ik,3) + zsa
341	#endif
342	#if ! defined key_vectopt_loop \|\| defined key_autotasking
343	END DO
344	#endif
345	END DO
346
347	#endif
348
349	! 5. Along sigma advective trend
350	! -------------------------------
351	! ... Second order centered tracer flux at u and v-points
352
353	# if defined key_vectopt_loop && ! defined key_autotasking
354	jj = 1
355	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
356	# else
357	DO jj = 1, jpjm1
358	DO ji = 1, jpim1
359	# endif
360	iku = mbku(ji,jj)
361	ikv = mbkv(ji,jj)
362	zfui = zalphax(ji,jj) e2u(ji,jj) fse3u(ji,jj,iku) * zunb(ji,jj)
363	zfvj = zalphay(ji,jj) e1v(ji,jj) fse3v(ji,jj,ikv) * zvnb(ji,jj)
364	! centered scheme
365	! zwx(ji,jj) = 0.5* zfui * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
366	! zwy(ji,jj) = 0.5* zfvj * ( ztnb(ji,jj) + ztnb(ji,jj+1) )
367	! zww(ji,jj) = 0.5* zfui * ( zsnb(ji,jj) + zsnb(ji+1,jj) )
368	! zwz(ji,jj) = 0.5* zfvj * ( zsnb(ji,jj) + zsnb(ji,jj+1) )
369	! upstream scheme
370	zwx(ji,jj) = ( ( zfui + ABS( zfui ) ) * ztbb(ji ,jj ) &
371	& +( zfui - ABS( zfui ) ) * ztbb(ji+1,jj ) ) * 0.5
372	zwy(ji,jj) = ( ( zfui + ABS( zfvj ) ) * ztbb(ji ,jj ) &
373	& +( zfui - ABS( zfvj ) ) * ztbb(ji ,jj+1) ) * 0.5
374	zww(ji,jj) = ( ( zfui + ABS( zfui ) ) * zsbb(ji ,jj ) &
375	& +( zfui - ABS( zfui ) ) * zsbb(ji+1,jj ) ) * 0.5
376	zwz(ji,jj) = ( ( zfui + ABS( zfvj ) ) * zsbb(ji ,jj ) &
377	& +( zfui - ABS( zfvj ) ) * zsbb(ji ,jj+1) ) * 0.5
378	#if ! defined key_vectopt_loop \|\| defined key_autotasking
379	END DO
380	#endif
381	END DO
382	# if defined key_vectopt_loop && ! defined key_autotasking
383	jj = 1
384	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
385	# else
386	DO jj = 2, jpjm1
387	DO ji = 2, jpim1
388	# endif
389	ik = mbkt(ji,jj)
390	zbtr = 1. / ( e1t(ji,jj)e2t(ji,jj)fse3t(ji,jj,ik) )
391	! horizontal advective trends
392	zta = - zbtr * ( zwx(ji,jj) - zwx(ji-1,jj ) &
393	& + zwy(ji,jj) - zwy(ji ,jj-1) )
394	zsa = - zbtr * ( zww(ji,jj) - zww(ji-1,jj ) &
395	& + zwz(ji,jj) - zwz(ji ,jj-1) )
396
397	! add it to the general tracer trends
398	ta(ji,jj,ik) = ta(ji,jj,ik) + zta
399	sa(ji,jj,ik) = sa(ji,jj,ik) + zsa
400	#if defined key_trdtra
401	ttrd(ji,jj,ik,1) = ttrd(ji,jj,ik,1) + zta
402	strd(ji,jj,ik,1) = strd(ji,jj,ik,1) + zsa
403	#endif
404	#if ! defined key_vectopt_loop \|\| defined key_autotasking
405	END DO
406	#endif
407	END DO
408
409	IF(l_ctl) THEN
410	zta = SUM( ta(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
411	zsa = SUM( sa(2:nictl,2:njctl,1:jpkm1) * tmask(2:nictl,2:njctl,1:jpkm1) )
412	WRITE(numout,*) ' bbl - Ta: ', zta-t_ctl, ' Sa: ', zsa-s_ctl
413	t_ctl = zta ; s_ctl = zsa
414	ENDIF
415
416
417	! 6. Vertical advection velocities
418	! --------------------------------
419	! ... computes divergence perturbation (velocties to be removed from upper t boxes :
420	DO jk= 1, jpkm1
421	DO jj=1, jpjm1
422	DO ji = 1, fs_jpim1 ! vertor opt.
423	zwu(ji,jj) = -e2u(ji,jj) * u_bbl(ji,jj,jk)
424	zwv(ji,jj) = -e1v(ji,jj) * v_bbl(ji,jj,jk)
425	END DO
426	END DO
427
428	! ... horizontal divergence
429	DO jj = 2, jpjm1
430	DO ji = fs_2, fs_jpim1 ! vector opt.
431	zbt = e1t(ji,jj) * e2t(ji,jj)
432	zhdivn(ji,jj,jk) = ( zwu(ji,jj) - zwu(ji-1,jj ) &
433	+ zwv(ji,jj) - zwv(ji ,jj-1) ) / zbt
434	END DO
435	END DO
436	END DO
437
438
439	! ... horizontal bottom divergence
440	# if defined key_vectopt_loop && ! defined key_autotasking
441	jj = 1
442	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
443	# else
444	DO jj = 1, jpjm1
445	DO ji = 1, jpim1
446	# endif
447	iku = mbku(ji,jj)
448	ikv = mbkv(ji,jj)
449	zwu(ji,jj) = zalphax(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,iku)
450	zwv(ji,jj) = zalphay(ji,jj) * e1v(ji,jj) * fse3v(ji,jj,ikv)
451	#if ! defined key_vectopt_loop \|\| defined key_autotasking
452	END DO
453	#endif
454	END DO
455
456	# if defined key_vectopt_loop && ! defined key_autotasking
457	jj = 1
458	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
459	# else
460	DO jj = 2, jpjm1
461	DO ji = 2, jpim1
462	# endif
463	ik = mbkt(ji,jj)
464	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik)
465	zhdivn(ji,jj,ik) = &
466	& ( zwu(ji ,jj ) * ( zunb(ji ,jj ) - un(ji ,jj ,ik) *umask(ji ,jj ,1) ) &
467	& - zwu(ji-1,jj ) * ( zunb(ji-1,jj ) - un(ji-1,jj ,ik) *umask(ji-1,jj ,1) ) &
468	& + zwv(ji ,jj ) * ( zvnb(ji ,jj ) - vn(ji ,jj ,ik) *vmask(ji ,jj ,1) ) &
469	& - zwv(ji ,jj-1) * ( zvnb(ji ,jj-1) - vn(ji ,jj-1,ik) *vmask(ji ,jj-1,1) ) &
470	& ) / zbt
471
472	# if ! defined key_vectopt_loop \|\| defined key_autotasking
473	END DO
474	# endif
475	END DO
476
477	! 7. compute additional vertical velocity to be used in t boxes
478	! -------------------------------------------------------------
479
480	! ... Computation from the bottom
481	! Note that w_bbl(:,:,jpk) has been set to 0 in tra_bbl_init
482	DO jk = jpkm1, 1, -1
483	DO jj= 2, jpjm1
484	DO ji = fs_2, fs_jpim1 ! vector opt.
485	w_bbl(ji,jj,jk) = w_bbl(ji,jj,jk+1) - fse3t(ji,jj,jk)*zhdivn(ji,jj,jk)
486	END DO
487	END DO
488	END DO
489
490	! Boundary condition on w_bbl (unchanged sign)
491	CALL lbc_lnk( w_bbl, 'W', 1. )
492
493	END SUBROUTINE tra_bbl_adv

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

Download in other formats:

Original Format