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 @ 258

Last change on this file since 258 was 258, checked in by opalod, 19 years ago
nemo_v1_update_004 : CT : Integration of the control print option for debugging work
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 22.2 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
58	USE flxrnf
59	USE ocfzpt
60	USE lbclnk
61	USE oce, ONLY : ztdta => ua, & ! use ua as 3D workspace
62	ztdsa => va ! use va as 3D workspace
63
64	!! * Arguments
65	INTEGER, INTENT( in ) :: kt ! ocean time-step
66
67	!! * Local declarations
68	INTEGER :: ji, jj, jk ! dummy loop indices
69	INTEGER :: ik, iku, ikv ! temporary integers
70
71	REAL(wp) :: &
72	zsign, zt, zs, zh, zalbet, & ! temporary scalars
73	zgdrho, zbtr, zta, zsa ! " "
74	REAL(wp), DIMENSION(jpi,jpj) :: &
75	zki, zkj, zkw, zkx, zky, zkz, & ! temporary workspace arrays
76	ztnb, zsnb, zdep, ztbb, zsbb, & ! " "
77	zahu, zahv ! " "
78	REAL(wp), DIMENSION(jpi,jpj) :: & ! temporary workspace arrays
79	zalphax, zwu, zunb, & ! " "
80	zalphay, zwv, zvnb, & ! " "
81	zwx, zwy, zww, zwz ! " "
82	REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
83	zhdivn ! temporary workspace arrays
84	REAL(wp) :: &
85	zfui, zfvj, zbt, zsigna ! temporary scalars
86	REAL(wp) :: &
87	fsalbt, pft, pfs, pfh ! statement function
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	!!----------------------------------------------------------------------
109
110
111	IF( kt == nit000 ) CALL tra_bbl_init ! initialization at first time-step
112
113	! Save ta and sa trends
114	IF( l_trdtra ) THEN
115	ztdta(:,:,:) = ta(:,:,:)
116	ztdsa(:,:,:) = sa(:,:,:)
117	ENDIF
118
119	! 1. 2D fields of bottom temperature and salinity, and bottom slope
120	! -----------------------------------------------------------------
121	! mbathy= number of w-level, minimum value=1 (cf dommsk.F)
122
123	#if defined key_vectopt_loop && ! defined key_autotasking
124	jj = 1
125	DO ji = 1, jpij ! vector opt. (forced unrolling)
126	#else
127	DO jj = 1, jpj
128	DO ji = 1, jpi
129	#endif
130	ik = mbkt(ji,jj) ! index of the bottom ocean T-level
131	ztnb(ji,jj) = tn(ji,jj,ik) * tmask(ji,jj,1) ! masked now T at the ocean bottom
132	zsnb(ji,jj) = sn(ji,jj,ik) * tmask(ji,jj,1) ! masked now S at the ocean bottom
133	ztbb(ji,jj) = tb(ji,jj,ik) * tmask(ji,jj,1) ! masked before T at the ocean bottom
134	zsbb(ji,jj) = sb(ji,jj,ik) * tmask(ji,jj,1) ! masked before S at the ocean bottom
135	zdep(ji,jj) = fsdept(ji,jj,ik) ! depth of the ocean bottom T-level
136	#if ! defined key_vectopt_loop \|\| defined key_autotasking
137	END DO
138	#endif
139	END DO
140	#if defined key_vectopt_loop && ! defined key_autotasking
141	jj = 1
142	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
143	zunb(ji,jj) = un(ji,jj,mbku(ji,jj)) * umask(ji,jj,1)
144	zvnb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) * vmask(ji,jj,1) ! retirer le mask en u, v et t !
145	END DO
146	#else
147	DO jj = 1, jpjm1
148	DO ji = 1, jpim1
149	zunb(ji,jj) = un(ji,jj,mbku(ji,jj)) * umask(ji,jj,1)
150	zvnb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) * vmask(ji,jj,1)
151	END DO
152	END DO
153	#endif
154
155	! boundary conditions on zunb and zvnb (changed sign)
156	CALL lbc_lnk( zunb, 'U', -1. ) ; CALL lbc_lnk( zvnb, 'V', -1. )
157	! boundary condition on ztnb and znbb
158	CALL lbc_lnk( ztnb, 'T', 1. ) ; CALL lbc_lnk( ztbb, 'T', 1. )
159	! boundary condition on zsnb and zsbb
160	CALL lbc_lnk( zsnb, 'T', 1. ) ; CALL lbc_lnk( zsbb, 'T', 1. )
161
162	! Conditional diffusion along the slope in the bottom boundary layer
163
164	#if defined key_trabbl_dif
165	# if defined key_vectopt_loop && ! defined key_autotasking
166	jj = 1
167	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
168	# else
169	DO jj = 1, jpjm1
170	DO ji = 1, jpim1
171	# endif
172	iku = mbku(ji,jj)
173	ikv = mbkv(ji,jj)
174	zahu(ji,jj) = atrbble2u(ji,jj)fse3u(ji,jj,iku)/e1u(ji,jj) * umask(ji,jj,1)
175	zahv(ji,jj) = atrbble1v(ji,jj)fse3v(ji,jj,ikv)/e2v(ji,jj) * vmask(ji,jj,1)
176	# if ! defined key_vectopt_loop \|\| defined key_autotasking
177	END DO
178	# endif
179	END DO
180	#endif
181
182
183	! 2. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0
184	! --------------------------------------------
185	! Sign of the local density gradient along the i- and j-slopes
186	! multiplied by the slope of the ocean bottom
187
188	SELECT CASE ( neos )
189
190	CASE ( 0 ) ! Jackett and McDougall (1994) formulation
191
192	DO jj = 1, jpjm1
193	DO ji = 1, fs_jpim1 ! vector opt.
194	! ... temperature, salinity anomalie and depth
195	zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
196	zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
197	zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
198	! ... masked ratio alpha/beta
199	zalbet = fsalbt( zt, zs, zh )*umask(ji,jj,1)
200	! ... local density gradient along i-bathymetric slope
201	zgdrho = zalbet*( ztnb(ji+1,jj) - ztnb(ji,jj) ) &
202	- ( zsnb(ji+1,jj) - zsnb(ji,jj) )
203	zgdrho = zgdrho * umask(ji,jj,1)
204	! ... sign of local i-gradient of density multiplied by the i-slope
205	zsign = sign( 0.5, -zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
206	zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj)
207
208	zsigna= sign(0.5, zunb(ji,jj)*( zdep(ji+1,jj) - zdep(ji,jj) ))
209	zalphax(ji,jj)=(0.5+zsigna)(0.5-zsign)umask(ji,jj,1)
210	END DO
211	END DO
212
213	DO jj = 1, jpjm1
214	DO ji = 1, fs_jpim1 ! vector opt.
215	! ... temperature, salinity anomalie and depth
216	zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
217	zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
218	zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
219	! ... masked ratio alpha/beta
220	zalbet = fsalbt( zt, zs, zh )*vmask(ji,jj,1)
221	! ... local density gradient along j-bathymetric slope
222	zgdrho = zalbet*( ztnb(ji,jj+1) - ztnb(ji,jj) ) &
223	- ( zsnb(ji,jj+1) - zsnb(ji,jj) )
224	zgdrho = zgdrho*vmask(ji,jj,1)
225	! ... sign of local j-gradient of density multiplied by the j-slope
226	zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
227	zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj)
228
229	zsigna= sign(0.5, zvnb(ji,jj)*(zdep(ji,jj+1) - zdep(ji,jj) ) )
230	zalphay(ji,jj)=(0.5+zsigna)(0.5-zsign)vmask(ji,jj,1)
231	END DO
232	END DO
233
234
235	CASE ( 1 ) ! Linear formulation function of temperature only
236
237	IF(lwp) WRITE(numout,cform_err)
238	IF(lwp) WRITE(numout,) ' use of linear eos rho(T,S) = rau0 ( rbeta * S - ralpha * T )'
239	IF(lwp) WRITE(numout,*) ' bbl not implented: easy to do it '
240	nstop = nstop + 1
241
242	CASE ( 2 ) ! Linear formulation function of temperature and salinity
243
244	IF(lwp) WRITE(numout,cform_err)
245	IF(lwp) WRITE(numout,) ' use of linear eos rho(T,S) = rau0 ( rbeta * S - ralpha * T )'
246	IF(lwp) WRITE(numout,*) ' bbl not implented: easy to do it '
247	nstop = nstop + 1
248
249	CASE DEFAULT
250
251	IF(lwp) WRITE(numout,cform_err)
252	IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
253	nstop = nstop + 1
254
255	END SELECT
256
257	! lateral boundary conditions on zalphax and zalphay (unchanged sign)
258	CALL lbc_lnk( zalphax, 'U', 1. ) ; CALL lbc_lnk( zalphay, 'V', 1. )
259
260
261	! 3. Velocities that are exchanged between ajacent bottom boxes.
262	!---------------------------------------------------------------
263
264	! ... is equal to zero but where bbl will work.
265	u_bbl(:,:,:) = 0.e0
266	v_bbl(:,:,:) = 0.e0
267	# if defined key_vectopt_loop && ! defined key_autotasking
268	jj = 1
269	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
270	# else
271	DO jj = 1, jpjm1
272	DO ji = 1, jpim1
273	# endif
274	iku = mbku(ji,jj)
275	ikv = mbkv(ji,jj)
276	IF( MAX(iku,ikv) > 1 ) THEN
277	u_bbl(ji,jj,iku) = zalphax(ji,jj) * un(ji,jj,iku) * umask(ji,jj,1)
278	v_bbl(ji,jj,ikv) = zalphay(ji,jj) * vn(ji,jj,ikv) * vmask(ji,jj,1)
279	ENDIF
280	# if ! defined key_vectopt_loop \|\| defined key_autotasking
281	END DO
282	# endif
283	END DO
284
285	! lateral boundary conditions on u_bbl and v_bbl (changed sign)
286	CALL lbc_lnk( u_bbl, 'U', -1. ) ; CALL lbc_lnk( v_bbl, 'V', -1. )
287
288
289
290	#if defined key_trabbl_dif
291	! 4. Additional second order diffusive trends
292	! -------------------------------------------
293
294	! ... first derivative (gradient)
295	DO jj = 1, jpjm1
296	DO ji = 1, fs_jpim1 ! vertor opt.
297	zkx(ji,jj) = zki(ji,jj)*( ztbb(ji+1,jj) - ztbb(ji,jj) )
298	zkz(ji,jj) = zki(ji,jj)*( zsbb(ji+1,jj) - zsbb(ji,jj) )
299
300	zky(ji,jj) = zkj(ji,jj)*( ztbb(ji,jj+1) - ztbb(ji,jj) )
301	zkw(ji,jj) = zkj(ji,jj)*( zsbb(ji,jj+1) - zsbb(ji,jj) )
302	END DO
303	END DO
304
305	IF( cp_cfg == "orca" ) THEN
306	SELECT CASE ( jp_cfg )
307	! ! =======================
308	CASE ( 2 ) ! ORCA_R2 configuration
309	! ! =======================
310	! Gibraltar enhancement of BBL
311	zkx( mi0(139):mi1(140) , mj0(102):mj1(102) ) = 4.e0 * zkx( mi0(139):mi1(140) , mj0(102):mj1(102) )
312	zky( mi0(139):mi1(140) , mj0(102):mj1(102) ) = 4.e0 * zky( mi0(139):mi1(140) , mj0(102):mj1(102) )
313
314	! Red Sea enhancement of BBL
315	zkx( mi0(161):mi1(162) , mj0(88):mj1(88) ) = 10.e0 * zkx( mi0(161):mi1(162) , mj0(88):mj1(88) )
316	zky( mi0(161):mi1(162) , mj0(88):mj1(88) ) = 10.e0 * zky( mi0(161):mi1(162) , mj0(88):mj1(88) )
317
318	! ! =======================
319	CASE ( 4 ) ! ORCA_R4 configuration
320	! ! =======================
321	! Gibraltar enhancement of BBL
322	zkx( mi0(70):mi1(71) , mj0(52):mj1(52) ) = 4.e0 * zkx( mi0(70):mi1(71) , mj0(52):mj1(52) )
323	zky( mi0(70):mi1(71) , mj0(52):mj1(52) ) = 4.e0 * zky( mi0(70):mi1(71) , mj0(52):mj1(52) )
324
325	END SELECT
326
327	ENDIF
328
329	! ... second derivative (divergence) and add to the general tracer trend
330
331	# if defined key_vectopt_loop && ! defined key_autotasking
332	jj = 1
333	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
334	# else
335	DO jj = 2, jpjm1
336	DO ji = 2, jpim1
337	# endif
338	ik = mbkt(ji,jj)
339	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik) )
340	zta = ( zkx(ji,jj) - zkx(ji-1,jj ) &
341	& + zky(ji,jj) - zky(ji ,jj-1) ) * zbtr
342	zsa = ( zkz(ji,jj) - zkz(ji-1,jj ) &
343	& + zkw(ji,jj) - zkw(ji ,jj-1) ) * zbtr
344	ta(ji,jj,ik) = ta(ji,jj,ik) + zta
345	sa(ji,jj,ik) = sa(ji,jj,ik) + zsa
346	#if ! defined key_vectopt_loop \|\| defined key_autotasking
347	END DO
348	#endif
349	END DO
350
351	! save the trends for diagnostic
352	! BBL lateral diffusion tracers trends
353	IF( l_trdtra ) THEN
354	# if defined key_vectopt_loop && ! defined key_autotasking
355	jj = 1
356	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
357	# else
358	DO jj = 2, jpjm1
359	DO ji = 2, jpim1
360	# endif
361	ik = mbkt(ji,jj)
362	tldfbbl(ji,jj) = ta(ji,jj,ik) - ztdta(ji,jj,ik)
363	sldfbbl(ji,jj) = sa(ji,jj,ik) - ztdsa(ji,jj,ik)
364	# if ! defined key_vectopt_loop \|\| defined key_autotasking
365	END DO
366	# endif
367	END DO
368
369	! save the new ta & sa trends
370	ztdta(:,:,:) = ta(:,:,:)
371	ztdsa(:,:,:) = sa(:,:,:)
372
373	ENDIF
374
375	#endif
376
377	! 5. Along sigma advective trend
378	! -------------------------------
379	! ... Second order centered tracer flux at u and v-points
380
381	# if defined key_vectopt_loop && ! defined key_autotasking
382	jj = 1
383	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
384	# else
385	DO jj = 1, jpjm1
386	DO ji = 1, jpim1
387	# endif
388	iku = mbku(ji,jj)
389	ikv = mbkv(ji,jj)
390	zfui = zalphax(ji,jj) e2u(ji,jj) fse3u(ji,jj,iku) * zunb(ji,jj)
391	zfvj = zalphay(ji,jj) e1v(ji,jj) fse3v(ji,jj,ikv) * zvnb(ji,jj)
392	! centered scheme
393	! zwx(ji,jj) = 0.5* zfui * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
394	! zwy(ji,jj) = 0.5* zfvj * ( ztnb(ji,jj) + ztnb(ji,jj+1) )
395	! zww(ji,jj) = 0.5* zfui * ( zsnb(ji,jj) + zsnb(ji+1,jj) )
396	! zwz(ji,jj) = 0.5* zfvj * ( zsnb(ji,jj) + zsnb(ji,jj+1) )
397	! upstream scheme
398	zwx(ji,jj) = ( ( zfui + ABS( zfui ) ) * ztbb(ji ,jj ) &
399	& +( zfui - ABS( zfui ) ) * ztbb(ji+1,jj ) ) * 0.5
400	zwy(ji,jj) = ( ( zfui + ABS( zfvj ) ) * ztbb(ji ,jj ) &
401	& +( zfui - ABS( zfvj ) ) * ztbb(ji ,jj+1) ) * 0.5
402	zww(ji,jj) = ( ( zfui + ABS( zfui ) ) * zsbb(ji ,jj ) &
403	& +( zfui - ABS( zfui ) ) * zsbb(ji+1,jj ) ) * 0.5
404	zwz(ji,jj) = ( ( zfui + ABS( zfvj ) ) * zsbb(ji ,jj ) &
405	& +( zfui - ABS( zfvj ) ) * zsbb(ji ,jj+1) ) * 0.5
406	#if ! defined key_vectopt_loop \|\| defined key_autotasking
407	END DO
408	#endif
409	END DO
410	# if defined key_vectopt_loop && ! defined key_autotasking
411	jj = 1
412	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
413	# else
414	DO jj = 2, jpjm1
415	DO ji = 2, jpim1
416	# endif
417	ik = mbkt(ji,jj)
418	zbtr = 1. / ( e1t(ji,jj)e2t(ji,jj)fse3t(ji,jj,ik) )
419	! horizontal advective trends
420	zta = - zbtr * ( zwx(ji,jj) - zwx(ji-1,jj ) &
421	& + zwy(ji,jj) - zwy(ji ,jj-1) )
422	zsa = - zbtr * ( zww(ji,jj) - zww(ji-1,jj ) &
423	& + zwz(ji,jj) - zwz(ji ,jj-1) )
424
425	! add it to the general tracer trends
426	ta(ji,jj,ik) = ta(ji,jj,ik) + zta
427	sa(ji,jj,ik) = sa(ji,jj,ik) + zsa
428	#if ! defined key_vectopt_loop \|\| defined key_autotasking
429	END DO
430	#endif
431	END DO
432
433	! save the trends for diagnostic
434	! BBL lateral advection tracers trends
435	IF( l_trdtra ) THEN
436	# if defined key_vectopt_loop && ! defined key_autotasking
437	jj = 1
438	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
439	# else
440	DO jj = 2, jpjm1
441	DO ji = 2, jpim1
442	# endif
443	ik = mbkt(ji,jj)
444	tladbbl(ji,jj) = ta(ji,jj,ik) - ztdta(ji,jj,ik)
445	sladbbl(ji,jj) = sa(ji,jj,ik) - ztdsa(ji,jj,ik)
446	# if ! defined key_vectopt_loop \|\| defined key_autotasking
447	END DO
448	# endif
449	END DO
450
451	ENDIF
452
453	IF(ln_ctl) THEN
454	CALL prt_ctl(tab3d_1=ta, clinfo1=' bbl - Ta: ', mask1=tmask, &
455	& tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
456	ENDIF
457
458
459	! 6. Vertical advection velocities
460	! --------------------------------
461	! ... computes divergence perturbation (velocties to be removed from upper t boxes :
462	DO jk= 1, jpkm1
463	DO jj=1, jpjm1
464	DO ji = 1, fs_jpim1 ! vertor opt.
465	zwu(ji,jj) = -e2u(ji,jj) * u_bbl(ji,jj,jk)
466	zwv(ji,jj) = -e1v(ji,jj) * v_bbl(ji,jj,jk)
467	END DO
468	END DO
469
470	! ... horizontal divergence
471	DO jj = 2, jpjm1
472	DO ji = fs_2, fs_jpim1 ! vector opt.
473	zbt = e1t(ji,jj) * e2t(ji,jj)
474	zhdivn(ji,jj,jk) = ( zwu(ji,jj) - zwu(ji-1,jj ) &
475	+ zwv(ji,jj) - zwv(ji ,jj-1) ) / zbt
476	END DO
477	END DO
478	END DO
479
480
481	! ... horizontal bottom divergence
482	# if defined key_vectopt_loop && ! defined key_autotasking
483	jj = 1
484	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
485	# else
486	DO jj = 1, jpjm1
487	DO ji = 1, jpim1
488	# endif
489	iku = mbku(ji,jj)
490	ikv = mbkv(ji,jj)
491	zwu(ji,jj) = zalphax(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,iku)
492	zwv(ji,jj) = zalphay(ji,jj) * e1v(ji,jj) * fse3v(ji,jj,ikv)
493	#if ! defined key_vectopt_loop \|\| defined key_autotasking
494	END DO
495	#endif
496	END DO
497
498	# if defined key_vectopt_loop && ! defined key_autotasking
499	jj = 1
500	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
501	# else
502	DO jj = 2, jpjm1
503	DO ji = 2, jpim1
504	# endif
505	ik = mbkt(ji,jj)
506	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik)
507	zhdivn(ji,jj,ik) = &
508	& ( zwu(ji ,jj ) * ( zunb(ji ,jj ) - un(ji ,jj ,ik) *umask(ji ,jj ,1) ) &
509	& - zwu(ji-1,jj ) * ( zunb(ji-1,jj ) - un(ji-1,jj ,ik) *umask(ji-1,jj ,1) ) &
510	& + zwv(ji ,jj ) * ( zvnb(ji ,jj ) - vn(ji ,jj ,ik) *vmask(ji ,jj ,1) ) &
511	& - zwv(ji ,jj-1) * ( zvnb(ji ,jj-1) - vn(ji ,jj-1,ik) *vmask(ji ,jj-1,1) ) &
512	& ) / zbt
513
514	# if ! defined key_vectopt_loop \|\| defined key_autotasking
515	END DO
516	# endif
517	END DO
518
519	! 7. compute additional vertical velocity to be used in t boxes
520	! -------------------------------------------------------------
521
522	! ... Computation from the bottom
523	! Note that w_bbl(:,:,jpk) has been set to 0 in tra_bbl_init
524	DO jk = jpkm1, 1, -1
525	DO jj= 2, jpjm1
526	DO ji = fs_2, fs_jpim1 ! vector opt.
527	w_bbl(ji,jj,jk) = w_bbl(ji,jj,jk+1) - fse3t(ji,jj,jk)*zhdivn(ji,jj,jk)
528	END DO
529	END DO
530	END DO
531
532	! Boundary condition on w_bbl (unchanged sign)
533	CALL lbc_lnk( w_bbl, 'W', 1. )
534
535	END SUBROUTINE tra_bbl_adv

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

Download in other formats:

Original Format