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.

traadv_cen2_atsk.h90 in tags/nemo_v1_03/NEMO/OPA_SRC/TRA – NEMO

Context Navigation

source: tags/nemo_v1_03/NEMO/OPA_SRC/TRA/traadv_cen2_atsk.h90 @ 6556

Last change on this file since 6556 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: 16.7 KB

Line
1	!!----------------------------------------------------------------------
2	!! * traadv_cen2_atsk.h90 *
3	!!----------------------------------------------------------------------
4	!! tra_adv_cen2 : update the tracer trend with the horizontal and
5	!! vertical advection trends using a seconder order
6	!! centered scheme. Auto-tasking case, k-slab for
7	!! hor. adv., j-slab for vert. adv.
8	!!----------------------------------------------------------------------
9	!! OPA 9.0 , LOCEAN-IPSL (2005)
10	!! $Header$
11	!! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
12	!!----------------------------------------------------------------------
13
14	SUBROUTINE tra_adv_cen2( kt )
15	!!----------------------------------------------------------------------
16	!! * ROUTINE tra_adv_cen2 *
17	!!
18	!! ** Purpose : Compute the now trend due to the advection of tracers
19	!! and add it to the general trend of passive tracer equations.
20	!!
21	!! ** Method : The advection is evaluated by a second order centered
22	!! scheme using now fields (leap-frog scheme). In specific areas
23	!! (vicinity of major river mouths, some straits, or where tn is
24	!! approaching the freezing point) it is mixed with an upstream
25	!! scheme for stability reasons.
26	!! Part 0 : compute the upstream / centered flag
27	!! (3D array, zind, defined at T-point (0<zind<1))
28	!! Part I : horizontal advection
29	!! * centered flux:
30	!! * s-coordinate (lk_sco=T) or
31	!! * z-coordinate with partial steps (lk_zps=T),
32	!! the vertical scale factors e3. are inside the derivatives:
33	!! zcenu = e2u*e3u un mi(tn)
34	!! zcenv = e1v*e3v vn mj(tn)
35	!! * z-coordinate (lk_zco=T), e3t=e3u=e3v:
36	!! zcenu = e2u un mi(tn)
37	!! zcenv = e1v vn mj(tn)
38	!! * upstream flux:
39	!! * s-coordinate (lk_sco=T) or
40	!! * z-coordinate with partial steps (lk_zps=T)
41	!! zupsu = e2u*e3u un (tb(i) or tb(i-1) ) [un>0 or <0]
42	!! zupsv = e1v*e3v vn (tb(j) or tb(j-1) ) [vn>0 or <0]
43	!! * z-coordinate (default key)
44	!! zupsu = e2u*e3u un (tb(i) or tb(i-1) ) [un>0 or <0]
45	!! zupsv = e1v*e3v vn (tb(j) or tb(j-1) ) [vn>0 or <0]
46	!! * mixed upstream / centered horizontal advection scheme
47	!! zcofi = max(zind(i+1), zind(i))
48	!! zcofj = max(zind(j+1), zind(j))
49	!! zwx = zcofi * zupsu + (1-zcofi) * zcenu
50	!! zwy = zcofj * zupsv + (1-zcofj) * zcenv
51	!! * horizontal advective trend (divergence of the fluxes)
52	!! * s-coordinate (lk_sco=T)
53	!! or z-coordinate with partial steps (lk_zps=T)
54	!! zta = 1/(e1te2te3t) { di-1[zwx] + dj-1[zwy] }
55	!! * z-coordinate (default key), e3t=e3u=e3v:
56	!! zta = 1/(e1t*e2t) { di-1[zwx] + dj-1[zwy] }
57	!! * Add this trend now to the general trend of tracer (ta,sa):
58	!! (ta,sa) = (ta,sa) + ( zta , zsa )
59	!! * trend diagnostic ('key_trdtra' defined): the trend is
60	!! saved for diagnostics. The trends saved is expressed as
61	!! Uh.gradh(T), i.e.
62	!! save trend = zta + tn divn
63	!! In addition, the advective trend in the two horizontal direc-
64	!! tion is also re-computed as Uh gradh(T). Indeed hadt+tn divn is
65	!! equal to (in s-coordinates, and similarly in z-coord.):
66	!! zta+tndivn=1/(e1te2te3t) { mi-1( e2ue3u un di[tn] )
67	!! +mj-1( e1v*e3v vn mj[tn] ) }
68	!! C A U T I O N : the trend saved is the centered trend only.
69	!! It doesn't take into account the upstream part of the scheme.
70	!!
71	!! Part II : vertical advection
72	!! For temperature (idem for salinity) the advective trend is com-
73	!! puted as follows :
74	!! zta = 1/e3t dk+1[ zwz ]
75	!! where the vertical advective flux, zwz, is given by :
76	!! zwz = zcofk * zupst + (1-zcofk) * zcent
77	!! with
78	!! zupsv = upstream flux = wn * (tb(k) or tb(k-1) ) [wn>0 or <0]
79	!! zcenu = centered flux = wn * mk(tn)
80	!! The surface boundary condition is :
81	!! rigid-lid (lk_dynspg_frd = T) : zero advective flux
82	!! free-surf (lk_dynspg_fsc = T) : wn(:,:,1) * tn(:,:,1)
83	!! Add this trend now to the general trend of tracer (ta,sa):
84	!! (ta,sa) = (ta,sa) + ( zta , zsa )
85	!! Trend diagnostic ('key_trdtra' defined): the trend is
86	!! saved for diagnostics. The trends saved is expressed as :
87	!! save trend = w.gradz(T) = zta - tn divn.
88	!!
89	!! ** Action : - update (ta,sa) with the now advective tracer trends
90	!! - save trends in (ttrdh,ttrd,strdhi,strd) ('key_trdtra')
91	!!
92	!! History :
93	!! 8.2 ! 01-08 (G. Madec, E. Durand) trahad+trazad = traadv
94	!! 8.5 ! 02-06 (G. Madec) F90: Free form and module
95	!! 9.0 ! 04-08 (C. Talandier) New trends organization
96	!!----------------------------------------------------------------------
97	!! * Modules used
98	USE oce , zwx => ua, & ! use ua as workspace
99	& zwy => va ! use va as workspace
100	#if defined key_trabbl_adv
101	REAL(wp), DIMENSION(jpi,jpj,jpk) :: & ! temporary arrays
102	& zun, zvn, zwn
103	#else
104	USE oce , zun => un, & ! When no bbl, zun == un
105	& zvn => vn, & ! When no bbl, zvn == vn
106	& zwn => wn ! When no bbl, zwn == wn
107	#endif
108	!! * Arguments
109	INTEGER, INTENT( in ) :: kt ! ocean time-step index
110
111	!! * Local save
112	REAL(wp), DIMENSION(jpi,jpj), SAVE :: &
113	zbtr2
114
115	!! * Local declarations
116	INTEGER :: ji, jj, jk ! dummy loop indices
117	REAL(wp) :: &
118	zbtr, zta, zsa, zfui, zfvj, & ! temporary scalars
119	zhw, ze3tr, zcofi, zcofj, & ! " "
120	zupsut, zupsvt, zupsus, zupsvs, & ! " "
121	zfp_ui, zfp_vj, zfm_ui, zfm_vj, & ! " "
122	zcofk, zupst, zupss, zcent, & ! " "
123	zcens, zfp_w, zfm_w, & ! " "
124	zcenut, zcenvt, zcenus, zcenvs, & ! " "
125	zfui1, zfvj1 ! " "
126	REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
127	zwz, zww, zind, & ! temporary workspace arrays
128	ztdta, ztdsa ! " "
129	!!----------------------------------------------------------------------
130
131	IF( kt == nit000 ) THEN
132	IF(lwp) WRITE(numout,*)
133	IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme'
134	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ Auto-tasking case'
135	IF(lwp) WRITE(numout,*)
136
137	zbtr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
138	ENDIF
139
140	! Save ta and sa trends
141	IF( l_trdtra ) THEN
142	ztdta(:,:,:) = ta(:,:,:)
143	ztdsa(:,:,:) = sa(:,:,:)
144	l_adv = 'ce2'
145	ENDIF
146
147	! ! ===============
148	DO jk = 1, jpkm1 ! Horizontal slab
149	! ! ===============
150	#if defined key_trabbl_adv
151	! Advective bottom boundary layer
152	! -------------------------------
153	zun(:,:,jk) = un (:,:,jk) - u_bbl(:,:,jk)
154	zvn(:,:,jk) = vn (:,:,jk) - v_bbl(:,:,jk)
155	zwn(:,:,jk) = wn (:,:,jk) + w_bbl(:,:,jk)
156	#endif
157
158	! 0. Upstream / centered scheme indicator
159	! ---------------------------------------
160	DO jj = 1, jpj
161	DO ji = 1, jpi
162	zind(ji,jj,jk) = MAX ( &
163	upsrnfh(ji,jj) * upsrnfz(jk), & ! changing advection scheme near runoff
164	upsadv(ji,jj) & ! in the vicinity of some straits
165	#if defined key_ice_lim
166	, tmask(ji,jj,jk) & ! half upstream tracer fluxes if tn < ("freezing"+0.1 )
167	* MAX( 0., SIGN( 1., fzptn(ji,jj)+0.1-tn(ji,jj,jk) ) ) &
168	#endif
169	)
170	END DO
171	END DO
172
173
174	! I. Horizontal advective fluxes
175	! ------------------------------
176	! Second order centered tracer flux at u and v-points
177	DO jj = 1, jpjm1
178	DO ji = 1, fs_jpim1 ! vector opt.
179	! upstream indicator
180	zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
181	zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
182	! volume fluxes * 1/2
183	#if defined key_s_coord \|\| defined key_partial_steps
184	zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * zun(ji,jj,jk)
185	zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * zvn(ji,jj,jk)
186	#else
187	zfui = 0.5 * e2u(ji,jj) * zun(ji,jj,jk)
188	zfvj = 0.5 * e1v(ji,jj) * zvn(ji,jj,jk)
189	#endif
190	! upstream scheme
191	zfp_ui = zfui + ABS( zfui )
192	zfp_vj = zfvj + ABS( zfvj )
193	zfm_ui = zfui - ABS( zfui )
194	zfm_vj = zfvj - ABS( zfvj )
195	zupsut = zfp_ui * tb(ji,jj,jk) + zfm_ui * tb(ji+1,jj ,jk)
196	zupsvt = zfp_vj * tb(ji,jj,jk) + zfm_vj * tb(ji ,jj+1,jk)
197	zupsus = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj ,jk)
198	zupsvs = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji ,jj+1,jk)
199	! centered scheme
200	zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj ,jk) )
201	zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji ,jj+1,jk) )
202	zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj ,jk) )
203	zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji ,jj+1,jk) )
204	! mixed centered / upstream scheme
205	zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut
206	zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt
207	zww(ji,jj,jk) = zcofi * zupsus + (1.-zcofi) * zcenus
208	zwz(ji,jj,jk) = zcofj * zupsvs + (1.-zcofj) * zcenvs
209	END DO
210	END DO
211
212	! 2. Tracer flux divergence at t-point added to the general trend
213	! ---------------------------------------------------------------
214
215	DO jj = 2, jpjm1
216	DO ji = fs_2, fs_jpim1 ! vector opt.
217	#if defined key_s_coord \|\| defined key_partial_steps
218	zbtr = zbtr2(ji,jj) / fse3t(ji,jj,jk)
219	#else
220	zbtr = zbtr2(ji,jj)
221	#endif
222	! horizontal advective trends
223	zta = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk) &
224	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk) )
225	zsa = - zbtr * ( zww(ji,jj,jk) - zww(ji-1,jj ,jk) &
226	& + zwz(ji,jj,jk) - zwz(ji ,jj-1,jk) )
227	! add it to the general tracer trends
228	ta(ji,jj,jk) = ta(ji,jj,jk) + zta
229	sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
230	END DO
231	END DO
232	! ! ===============
233	END DO ! End of slab
234	! ! ===============
235
236	! 3. Save the horizontal advective trends for diagnostic
237	! ------------------------------------------------------
238	IF( l_trdtra ) THEN
239	! Recompute the hoizontal advection zta & zsa trends computed
240	! at the step 2. above in making the difference between the new
241	! trends and the previous one ta()/sa - ztdta()/ztdsa() and add
242	! the term tn()/sn()*hdivn() to recover the Uh gradh(T/S) trends
243	ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) + tn(:,:,:) * hdivn(:,:,:)
244	ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) + sn(:,:,:) * hdivn(:,:,:)
245
246	CALL trd_mod(ztdta, ztdsa, jpttdlad, 'TRA', kt)
247
248	! Save the new ta and sa trends
249	ztdta(:,:,:) = ta(:,:,:)
250	ztdsa(:,:,:) = sa(:,:,:)
251
252	ENDIF
253
254	IF(ln_ctl) THEN
255	CALL prt_ctl(tab3d_1=ta, clinfo1=' centered2 had - Ta: ', mask1=tmask, &
256	& tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
257	ENDIF
258
259	! "zonal" mean advective heat and salt transport
260	IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
261	# if defined key_s_coord \|\| defined key_partial_steps
262	pht_adv(:) = ptr_vj( zwy(:,:,:) )
263	pst_adv(:) = ptr_vj( zwz(:,:,:) )
264	# else
265	DO jk = 1, jpkm1
266	DO jj = 2, jpjm1
267	DO ji = fs_2, fs_jpim1 ! vector opt.
268	zwy(ji,jj,jk) = zwy(ji,jj,jk) * fse3v(ji,jj,jk)
269	zwz(ji,jj,jk) = zwz(ji,jj,jk) * fse3v(ji,jj,jk)
270	END DO
271	END DO
272	END DO
273	pht_adv(:) = ptr_vj( zwy(:,:,:) )
274	pst_adv(:) = ptr_vj( zwz(:,:,:) )
275	# endif
276	ENDIF
277
278	! II. Vertical advection
279	! ----------------------
280	! ! ===============
281	DO jj = 2, jpjm1 ! Vertical slab
282	! ! ===============
283	! Bottom value : flux and indicator set to zero
284	zwz (:,jj,jpk) = 0.e0 ; zww(:,jj,jpk) = 0.e0
285	zind(:,jj,jpk) = 0.e0
286
287	! Surface value
288	IF( lk_dynspg_fsc_tsk ) THEN
289	! free surface-constant volume
290	zwz(:,jj, 1 ) = zwn(:,jj,1) * tn(:,jj,1)
291	zww(:,jj, 1 ) = zwn(:,jj,1) * sn(:,jj,1)
292	ELSE
293	! rigid lid : flux set to zero
294	zwz(:,jj, 1 ) = 0.e0 ; zww(:,jj, 1 ) = 0.e0
295	ENDIF
296
297	! 1. Vertical advective fluxes
298	! ----------------------------
299	! Second order centered tracer flux at w-point
300	DO jk = 2, jpk
301	DO ji = 2, jpim1
302	! upstream indicator
303	zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
304	! velocity * 1/2
305	zhw = 0.5 * zwn(ji,jj,jk)
306	! upstream scheme
307	zfp_w = zhw + ABS( zhw )
308	zfm_w = zhw - ABS( zhw )
309	zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1)
310	zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1)
311	! centered scheme
312	zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) )
313	zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) )
314	! mixed centered / upstream scheme
315	zwz(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent
316	zww(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens
317	END DO
318	END DO
319
320	! 2. Tracer flux divergence at t-point added to the general trend
321	! -------------------------
322	DO jk = 1, jpkm1
323	DO ji = 2, jpim1
324	ze3tr = 1. / fse3t(ji,jj,jk)
325	! vertical advective trends
326	zta = - ze3tr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )
327	zsa = - ze3tr * ( zww(ji,jj,jk) - zww(ji,jj,jk+1) )
328	! add it to the general tracer trends
329	ta(ji,jj,jk) = ta(ji,jj,jk) + zta
330	sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
331	END DO
332	END DO
333	! ! ===============
334	END DO ! End of slab
335	! ! ===============
336
337	! 3. Save the vertical advective trends for diagnostic
338	! ----------------------------------------------------
339	IF( l_trdtra ) THEN
340	! Recompute the vertical advection zta & zsa trends computed
341	! at the step 2. above in making the difference between the new
342	! trends and the previous one: ta()/sa - ztdta()/ztdsa() and substract
343	! the term tn()/sn()*hdivn() to recover the W gradz(T/S) trends
344	ztdta(:,:,:) = ta(:,:,:) - ztdta(:,:,:) - tn(:,:,:) * hdivn(:,:,:)
345	ztdsa(:,:,:) = sa(:,:,:) - ztdsa(:,:,:) - sn(:,:,:) * hdivn(:,:,:)
346
347	CALL trd_mod(ztdta, ztdsa, jpttdzad, 'TRA', kt)
348	ENDIF
349
350	IF(ln_ctl) THEN
351	CALL prt_ctl(tab3d_1=ta, clinfo1=' centered2 zad - Ta: ', mask1=tmask, &
352	& tab3d_2=sa, clinfo2=' Sa: ', mask2=tmask, clinfo3='tra')
353	ENDIF
354
355	END SUBROUTINE tra_adv_cen2

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

Download in other formats:

Original Format