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

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

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

Download in other formats:

Original Format