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.

trcbbl_adv.h90 in trunk/NEMO/TOP_SRC/TRP – NEMO

Context Navigation

source: trunk/NEMO/TOP_SRC/TRP/trcbbl_adv.h90 @ 349

Last change on this file since 349 was 349, checked in by opalod, 18 years ago
nemo_v1_update_031 : CT : change header names
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 20.1 KB

Line
1	!!----------------------------------------------------------------------
2	!! * trcbbl_adv.h90 *
3	!!----------------------------------------------------------------------
4
5	!!----------------------------------------------------------------------
6	!! TOP 1.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 trc_bbl_adv( kt )
12	!!----------------------------------------------------------------------
13	!! * ROUTINE trc_bbl_adv *
14	!!
15	!! ** Purpose : Compute the before tracer 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 : tra =tra + 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 tra of the
40	!! botton ocean tracer point:
41	!! tra = tra + difft
42	!!
43	!! ** Action : - update tra at the bottom level with the bottom
44	!! boundary layer trend
45	!!
46	!! References :
47	!! Beckmann, A., and R. Doscher, 1997, J. Phys.Oceanogr., 581-591.
48	!!
49	!! History :
50	!! 8.5 ! 02-12 (A. de Miranda, G. Madec) Original Code
51	!! 9.0 ! 04-01 (A. de Miranda, G. Madec, J.M. Molines )
52	!! 9.0 ! 04-03 (C. Ethe) Adaptation for Passive tracers
53	!!----------------------------------------------------------------------
54	!! * Modules used
55	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
56
57	!! * Arguments
58	INTEGER, INTENT( in ) :: kt ! ocean time-step
59
60	!! * Local declarations
61	INTEGER :: ji, jj, jk, jn ! dummy loop indices
62	INTEGER :: ik, iku, ikv ! temporary integers
63
64	REAL(wp) :: &
65	zsign, zt, zs, zh, zalbet, & ! temporary scalars
66	zgdrho, zbtr, ztra ! " "
67	REAL(wp), DIMENSION(jpi,jpj) :: &
68	zki, zkj, zkw, zkx, zky, zkz, & ! temporary workspace arrays
69	ztnb, zsnb, zdep, ztrb, & ! " "
70	zahu, zahv ! " "
71	REAL(wp), DIMENSION(jpi,jpj) :: & ! temporary workspace arrays
72	zalphax, zwu, zunb, & ! " "
73	zalphay, zwv, zvnb, & ! " "
74	zwx, zwy ! " "
75	REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
76	zhdivn ! temporary workspace arrays
77	REAL(wp) :: &
78	zfui, zfvj, zbt, zsigna ! temporary scalars
79	REAL(wp) :: &
80	fsalbt, pft, pfs, pfh ! statement function
81	CHARACTER (len=22) :: charout
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	IF( kt == nittrc000 ) CALL trc_bbl_init ! initialization at first time-step
105
106	! 1. 2D fields of bottom temperature and salinity, and bottom slope
107	! -----------------------------------------------------------------
108	! mbathy= number of w-level, minimum value=1 (cf dommsk.F)
109
110	#if defined key_vectopt_loop && ! defined key_autotasking
111	jj = 1
112	DO ji = 1, jpij ! vector opt. (forced unrolling)
113	#else
114	DO jj = 1, jpj
115	DO ji = 1, jpi
116	#endif
117	ik = mbkt(ji,jj) ! index of the bottom ocean T-level
118	ztnb(ji,jj) = tn(ji,jj,ik) * tmask(ji,jj,1) ! masked now T at the ocean bottom
119	zsnb(ji,jj) = sn(ji,jj,ik) * tmask(ji,jj,1) ! masked now S at the ocean bottom
120	zdep(ji,jj) = fsdept(ji,jj,ik) ! depth of the ocean bottom T-level
121	#if ! defined key_vectopt_loop \|\| defined key_autotasking
122	END DO
123	#endif
124	END DO
125	#if defined key_vectopt_loop && ! defined key_autotasking
126	jj = 1
127	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
128	zunb(ji,jj) = un(ji,jj,mbku(ji,jj)) * umask(ji,jj,1)
129	zvnb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) * vmask(ji,jj,1) ! retirer le mask en u, v et t !
130	END DO
131	#else
132	DO jj = 1, jpjm1
133	DO ji = 1, jpim1
134	zunb(ji,jj) = un(ji,jj,mbku(ji,jj)) * umask(ji,jj,1)
135	zvnb(ji,jj) = vn(ji,jj,mbkv(ji,jj)) * vmask(ji,jj,1)
136	END DO
137	END DO
138	#endif
139
140	! boundary conditions on zunb and zvnb (changed sign)
141	CALL lbc_lnk( zunb, 'U', -1. ) ; CALL lbc_lnk( zvnb, 'V', -1. )
142
143	! Conditional diffusion along the slope in the bottom boundary layer
144
145	#if defined key_trcbbl_dif
146	# if defined key_vectopt_loop && ! defined key_autotasking
147	jj = 1
148	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
149	# else
150	DO jj = 1, jpjm1
151	DO ji = 1, jpim1
152	# endif
153	iku = mbku(ji,jj)
154	ikv = mbkv(ji,jj)
155	zahu(ji,jj) = atrbble2u(ji,jj)fse3u(ji,jj,iku)/e1u(ji,jj) * umask(ji,jj,1)
156	zahv(ji,jj) = atrbble1v(ji,jj)fse3v(ji,jj,ikv)/e2v(ji,jj) * vmask(ji,jj,1)
157	# if ! defined key_vectopt_loop \|\| defined key_autotasking
158	END DO
159	# endif
160	END DO
161	#endif
162
163
164	! 2. Criteria of additional bottom diffusivity: grad(rho).grad(h)<0
165	! --------------------------------------------
166	! Sign of the local density gradient along the i- and j-slopes
167	! multiplied by the slope of the ocean bottom
168
169	SELECT CASE ( neos )
170
171	CASE ( 0 ) ! Jackett and McDougall (1994) formulation
172
173	DO jj = 1, jpjm1
174	DO ji = 1, fs_jpim1 ! vector opt.
175	! ... temperature, salinity anomalie and depth
176	zt = 0.5 * ( ztnb(ji,jj) + ztnb(ji+1,jj) )
177	zs = 0.5 * ( zsnb(ji,jj) + zsnb(ji+1,jj) ) - 35.0
178	zh = 0.5 * ( zdep(ji,jj) + zdep(ji+1,jj) )
179	! ... masked ratio alpha/beta
180	zalbet = fsalbt( zt, zs, zh )*umask(ji,jj,1)
181	! ... local density gradient along i-bathymetric slope
182	zgdrho = zalbet*( ztnb(ji+1,jj) - ztnb(ji,jj) ) &
183	- ( zsnb(ji+1,jj) - zsnb(ji,jj) )
184	zgdrho = zgdrho * umask(ji,jj,1)
185	! ... sign of local i-gradient of density multiplied by the i-slope
186	zsign = sign( 0.5, -zgdrho * ( zdep(ji+1,jj) - zdep(ji,jj) ) )
187	zki(ji,jj) = ( 0.5 - zsign ) * zahu(ji,jj)
188
189	zsigna= sign(0.5, zunb(ji,jj)*( zdep(ji+1,jj) - zdep(ji,jj) ))
190	zalphax(ji,jj)=(0.5+zsigna)(0.5-zsign)umask(ji,jj,1)
191	END DO
192	END DO
193
194	DO jj = 1, jpjm1
195	DO ji = 1, fs_jpim1 ! vector opt.
196	! ... temperature, salinity anomalie and depth
197	zt = 0.5 * ( ztnb(ji,jj+1) + ztnb(ji,jj) )
198	zs = 0.5 * ( zsnb(ji,jj+1) + zsnb(ji,jj) ) - 35.0
199	zh = 0.5 * ( zdep(ji,jj+1) + zdep(ji,jj) )
200	! ... masked ratio alpha/beta
201	zalbet = fsalbt( zt, zs, zh )*vmask(ji,jj,1)
202	! ... local density gradient along j-bathymetric slope
203	zgdrho = zalbet*( ztnb(ji,jj+1) - ztnb(ji,jj) ) &
204	- ( zsnb(ji,jj+1) - zsnb(ji,jj) )
205	zgdrho = zgdrho*vmask(ji,jj,1)
206	! ... sign of local j-gradient of density multiplied by the j-slope
207	zsign = sign( 0.5, -zgdrho * ( zdep(ji,jj+1) - zdep(ji,jj) ) )
208	zkj(ji,jj) = ( 0.5 - zsign ) * zahv(ji,jj)
209
210	zsigna= sign(0.5, zvnb(ji,jj)*(zdep(ji,jj+1) - zdep(ji,jj) ) )
211	zalphay(ji,jj)=(0.5+zsigna)(0.5-zsign)vmask(ji,jj,1)
212	END DO
213	END DO
214
215
216	CASE ( 1 ) ! Linear formulation function of temperature only
217
218	IF(lwp) WRITE(numout,cform_err)
219	IF(lwp) WRITE(numout,) ' use of linear eos rho(T,S) = rau0 ( rbeta * S - ralpha * T )'
220	IF(lwp) WRITE(numout,*) ' bbl not implented: easy to do it '
221	nstop = nstop + 1
222
223	CASE ( 2 ) ! Linear formulation function of temperature and salinity
224
225	IF(lwp) WRITE(numout,cform_err)
226	IF(lwp) WRITE(numout,) ' use of linear eos rho(T,S) = rau0 ( rbeta * S - ralpha * T )'
227	IF(lwp) WRITE(numout,*) ' bbl not implented: easy to do it '
228	nstop = nstop + 1
229
230	CASE DEFAULT
231
232	IF(lwp) WRITE(numout,cform_err)
233	IF(lwp) WRITE(numout,*) ' bad flag value for neos = ', neos
234	nstop = nstop + 1
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	! ... is equal to zero but where bbl will work.
245	u_trc_bbl(:,:,:) = 0.e0
246	v_trc_bbl(:,:,:) = 0.e0
247	# if defined key_vectopt_loop && ! defined key_autotasking
248	jj = 1
249	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
250	# else
251	DO jj = 1, jpjm1
252	DO ji = 1, jpim1
253	# endif
254	iku = mbku(ji,jj)
255	ikv = mbkv(ji,jj)
256	IF( MAX(iku,ikv) > 1 ) THEN
257	u_trc_bbl(ji,jj,iku) = zalphax(ji,jj) * un(ji,jj,iku) * umask(ji,jj,1)
258	v_trc_bbl(ji,jj,ikv) = zalphay(ji,jj) * vn(ji,jj,ikv) * vmask(ji,jj,1)
259	ENDIF
260	# if ! defined key_vectopt_loop \|\| defined key_autotasking
261	END DO
262	# endif
263	END DO
264
265	! lateral boundary conditions on u_trc_bbl and v_trc_bbl (changed sign)
266	CALL lbc_lnk( u_trc_bbl, 'U', -1. ) ; CALL lbc_lnk( v_trc_bbl, 'V', -1. )
267
268
269
270	DO jn = 1, jptra
271
272	#if defined key_vectopt_loop && ! defined key_autotasking
273	jj = 1
274	DO ji = 1, jpij ! vector opt. (forced unrolling)
275	#else
276	DO jj = 1, jpj
277	DO ji = 1, jpi
278	#endif
279	ik = mbkt(ji,jj) ! index of the bottom ocean T-level
280	ztrb(ji,jj) = trb(ji,jj,ik,jn) * tmask(ji,jj,1) ! masked now T at the ocean bottom
281	#if ! defined key_vectopt_loop \|\| defined key_autotasking
282	END DO
283	#endif
284	END DO
285
286	#if defined key_trcbbl_dif
287	! 4. Additional second order diffusive trends
288	! -------------------------------------------
289
290	! ... first derivative (gradient)
291	DO jj = 1, jpjm1
292	DO ji = 1, fs_jpim1 ! vertor opt.
293	zkx(ji,jj) = zki(ji,jj)*( ztrb(ji+1,jj) - ztrb(ji,jj) )
294	zky(ji,jj) = zkj(ji,jj)*( ztrb(ji,jj+1) - ztrb(ji,jj) )
295	END DO
296	END DO
297
298	IF( cp_cfg == "orca" ) THEN
299	SELECT CASE ( jp_cfg )
300	! ! =======================
301	CASE ( 2 ) ! ORCA_R2 configuration
302	! ! =======================
303	! Gibraltar enhancement of BBL
304	zkx( mi0(139):mi1(140) , mj0(102):mj1(102) ) = 4.e0 * zkx( mi0(139):mi1(140) , mj0(102):mj1(102) )
305	zky( mi0(139):mi1(140) , mj0(102):mj1(102) ) = 4.e0 * zky( mi0(139):mi1(140) , mj0(102):mj1(102) )
306
307	! Red Sea enhancement of BBL
308	zkx( mi0(161):mi1(162) , mj0(88):mj1(88) ) = 10.e0 * zkx( mi0(161):mi1(162) , mj0(88):mj1(88) )
309	zky( mi0(161):mi1(162) , mj0(88):mj1(88) ) = 10.e0 * zky( mi0(161):mi1(162) , mj0(88):mj1(88) )
310
311	! ! =======================
312	CASE ( 4 ) ! ORCA_R4 configuration
313	! ! =======================
314	! Gibraltar enhancement of BBL
315	zkx( mi0(70):mi1(71) , mj0(52):mj1(52) ) = 4.e0 * zkx( mi0(70):mi1(71) , mj0(52):mj1(52) )
316	zky( mi0(70):mi1(71) , mj0(52):mj1(52) ) = 4.e0 * zky( mi0(70):mi1(71) , mj0(52):mj1(52) )
317
318	END SELECT
319
320	ENDIF
321
322	! ... second derivative (divergence) and add to the general tracer trend
323
324	# if defined key_vectopt_loop && ! defined key_autotasking
325	jj = 1
326	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
327	# else
328	DO jj = 2, jpjm1
329	DO ji = 2, jpim1
330	# endif
331	ik = mbkt(ji,jj)
332	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik) )
333	ztra = ( zkx(ji,jj) - zkx(ji-1,jj ) &
334	& + zky(ji,jj) - zky(ji ,jj-1) ) * zbtr
335	tra(ji,jj,ik,jn) = tra(ji,jj,ik,jn) + ztra
336	#if ! defined key_vectopt_loop \|\| defined key_autotasking
337	END DO
338	#endif
339	END DO
340
341	#endif
342
343
344	! 5. Along sigma advective trend
345	! -------------------------------
346	! ... Second order centered tracer flux at u and v-points
347
348	# if defined key_vectopt_loop && ! defined key_autotasking
349	jj = 1
350	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
351	# else
352	DO jj = 1, jpjm1
353	DO ji = 1, jpim1
354	# endif
355	iku = mbku(ji,jj)
356	ikv = mbkv(ji,jj)
357	zfui = zalphax(ji,jj) e2u(ji,jj) fse3u(ji,jj,iku) * zunb(ji,jj)
358	zfvj = zalphay(ji,jj) e1v(ji,jj) fse3v(ji,jj,ikv) * zvnb(ji,jj)
359	! upstream scheme
360	zwx(ji,jj) = ( ( zfui + ABS( zfui ) ) * ztrb(ji ,jj ) &
361	& +( zfui - ABS( zfui ) ) * ztrb(ji+1,jj ) ) * 0.5
362	zwy(ji,jj) = ( ( zfui + ABS( zfvj ) ) * ztrb(ji ,jj ) &
363	& +( zfui - ABS( zfvj ) ) * ztrb(ji ,jj+1) ) * 0.5
364	#if ! defined key_vectopt_loop \|\| defined key_autotasking
365	END DO
366	#endif
367	END DO
368	# if defined key_vectopt_loop && ! defined key_autotasking
369	jj = 1
370	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
371	# else
372	DO jj = 2, jpjm1
373	DO ji = 2, jpim1
374	# endif
375	ik = mbkt(ji,jj)
376	zbtr = 1. / ( e1t(ji,jj)e2t(ji,jj)fse3t(ji,jj,ik) )
377	! horizontal advective trends
378	ztra = - zbtr * ( zwx(ji,jj) - zwx(ji-1,jj ) &
379	& + zwy(ji,jj) - zwy(ji ,jj-1) )
380
381	! add it to the general tracer trends
382	tra(ji,jj,ik,jn) = tra(ji,jj,ik,jn) + ztra
383	#if ! defined key_vectopt_loop \|\| defined key_autotasking
384	END DO
385	#endif
386	END DO
387
388	END DO
389
390	IF(ln_ctl) THEN ! print mean trends (used for debugging)
391	WRITE(charout, FMT="('bbl - adv')")
392	CALL prt_ctl_trc_info(charout)
393	CALL prt_ctl_trc(tab4d=tra, mask=tmask, clinfo=ctrcnm,clinfo2='trd')
394	ENDIF
395	! 6. Vertical advection velocities
396	! --------------------------------
397	! ... computes divergence perturbation (velocties to be removed from upper t boxes :
398	DO jk= 1, jpkm1
399
400	DO jj=1, jpjm1
401	DO ji = 1, fs_jpim1 ! vertor opt.
402	zwu(ji,jj) = -e2u(ji,jj) * u_trc_bbl(ji,jj,jk)
403	zwv(ji,jj) = -e1v(ji,jj) * v_trc_bbl(ji,jj,jk)
404	END DO
405	END DO
406
407	! ... horizontal divergence
408	DO jj = 2, jpjm1
409	DO ji = fs_2, fs_jpim1 ! vector opt.
410	zbt = e1t(ji,jj) * e2t(ji,jj)
411	zhdivn(ji,jj,jk) = ( zwu(ji,jj) - zwu(ji-1,jj ) &
412	& + zwv(ji,jj) - zwv(ji ,jj-1) ) / zbt
413	END DO
414	END DO
415	END DO
416
417
418	! ... horizontal bottom divergence
419	# if defined key_vectopt_loop && ! defined key_autotasking
420	jj = 1
421	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
422	# else
423	DO jj = 1, jpjm1
424	DO ji = 1, jpim1
425	# endif
426	iku = mbku(ji,jj)
427	ikv = mbkv(ji,jj)
428	zwu(ji,jj) = zalphax(ji,jj) * e2u(ji,jj) * fse3u(ji,jj,iku)
429	zwv(ji,jj) = zalphay(ji,jj) * e1v(ji,jj) * fse3v(ji,jj,ikv)
430	#if ! defined key_vectopt_loop \|\| defined key_autotasking
431	END DO
432	#endif
433	END DO
434
435	# if defined key_vectopt_loop && ! defined key_autotasking
436	jj = 1
437	DO ji = jpi+2, jpij-jpi-1 ! vector opt. (forced unrolling)
438	# else
439	DO jj = 2, jpjm1
440	DO ji = 2, jpim1
441	# endif
442	ik = mbkt(ji,jj)
443	zbt = e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,ik)
444	zhdivn(ji,jj,ik) = &
445	& ( zwu(ji ,jj ) * ( zunb(ji ,jj ) - un(ji ,jj ,ik) *umask(ji ,jj ,1) ) &
446	& - zwu(ji-1,jj ) * ( zunb(ji-1,jj ) - un(ji-1,jj ,ik) *umask(ji-1,jj ,1) ) &
447	& + zwv(ji ,jj ) * ( zvnb(ji ,jj ) - vn(ji ,jj ,ik) *vmask(ji ,jj ,1) ) &
448	& - zwv(ji ,jj-1) * ( zvnb(ji ,jj-1) - vn(ji ,jj-1,ik) *vmask(ji ,jj-1,1) ) &
449	& ) / zbt
450
451	# if ! defined key_vectopt_loop \|\| defined key_autotasking
452	END DO
453	# endif
454	END DO
455
456	! 7. compute additional vertical velocity to be used in t boxes
457	! -------------------------------------------------------------
458
459	! ... Computation from the bottom
460	! Note that w_trc_bbl(:,:,jpk) has been set to 0 in trc_bbl_init
461	DO jk = jpkm1, 1, -1
462	DO jj= 2, jpjm1
463	DO ji = fs_2, fs_jpim1 ! vector opt.
464	w_trc_bbl(ji,jj,jk) = w_trc_bbl(ji,jj,jk+1) - fse3t(ji,jj,jk)*zhdivn(ji,jj,jk)
465	END DO
466	END DO
467	END DO
468
469	! Boundary condition on w_bbl (unchanged sign)
470	CALL lbc_lnk( w_trc_bbl, 'W', 1. )
471
472	END SUBROUTINE trc_bbl_adv

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

Download in other formats:

Original Format