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

Last change on this file since 202 was 202, checked in by opalod, 19 years ago
CT : UPDATE142 : Check the consistency between passive tracers transport modules (in TRP directory) and those used for the active tracers
Property svn:eol-style set to `native` Property svn:executable set to ``* Property svn:keywords set to `Author Date Id Revision`
File size: 20.0 KB

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

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

Download in other formats:

Original Format