Context Navigation

traadv_cen2.F90 @ 2696

Last change on this file since 2696 was 2690, checked in by gm, 13 years ago
dynamic mem: #785 ; homogeneization of the coding style associated with dyn allocation
Property svn:keywords set to `Id`
File size: 18.3 KB

Line
1	MODULE traadv_cen2
2	!!======================================================================
3	!! * MODULE traadv_cen2 *
4	!! Ocean tracers: horizontal & vertical advective trend
5	!!======================================================================
6	!! History : 8.2 ! 2001-08 (G. Madec, E. Durand) trahad+trazad=traadv
7	!! 1.0 ! 2002-06 (G. Madec) F90: Free form and module
8	!! 9.0 ! 2004-08 (C. Talandier) New trends organization
9	!! - ! 2005-11 (V. Garnier) Surface pressure gradient organization
10	!! 2.0 ! 2006-04 (R. Benshila, G. Madec) Step reorganization
11	!! - ! 2006-07 (G. madec) add ups_orca_set routine
12	!! 3.2 ! 2009-07 (G. Madec) add avmb, avtb in restart for cen2 advection
13	!! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
14	!!----------------------------------------------------------------------
15
16	!!----------------------------------------------------------------------
17	!! tra_adv_cen2 : update the tracer trend with the advection trends using a 2nd order centered scheme
18	!! ups_orca_set : allow mixed upstream/centered scheme in specific area (set for orca 2 and 4 only)
19	!!----------------------------------------------------------------------
20	USE oce, ONLY: tsn ! now ocean temperature and salinity
21	USE dom_oce ! ocean space and time domain
22	USE eosbn2 ! equation of state
23	USE trdmod_oce ! tracers trends
24	USE trdtra ! tracers trends
25	USE closea ! closed sea
26	USE sbcrnf ! river runoffs
27	USE in_out_manager ! I/O manager
28	USE iom ! IOM library
29	USE diaptr ! poleward transport diagnostics
30	USE zdf_oce ! ocean vertical physics
31	USE restart ! ocean restart
32	USE trc_oce ! share passive tracers/Ocean variables
33	USE lib_mpp ! MPP library
34
35	IMPLICIT NONE
36	PRIVATE
37
38	PUBLIC tra_adv_cen2 ! routine called by step.F90
39	PUBLIC ups_orca_set ! routine used by traadv_cen2_jki.F90
40
41	LOGICAL :: l_trd ! flag to compute trends
42
43	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits
44	! ! and in closed seas (orca 2 and 4 configurations)
45	!! * Substitutions
46	# include "domzgr_substitute.h90"
47	# include "vectopt_loop_substitute.h90"
48	!!----------------------------------------------------------------------
49	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
50	!! $Id$
51	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
52	!!----------------------------------------------------------------------
53	CONTAINS
54
55	SUBROUTINE tra_adv_cen2( kt, cdtype, pun, pvn, pwn, &
56	& ptb, ptn, pta, kjpt )
57	!!----------------------------------------------------------------------
58	!! * ROUTINE tra_adv_cen2 *
59	!!
60	!! ** Purpose : Compute the now trend due to the advection of tracers
61	!! and add it to the general trend of passive tracer equations.
62	!!
63	!! ** Method : The advection is evaluated by a second order centered
64	!! scheme using now fields (leap-frog scheme). In specific areas
65	!! (vicinity of major river mouths, some straits, or where tn is
66	!! approaching the freezing point) it is mixed with an upstream
67	!! scheme for stability reasons.
68	!! Part 0 : compute the upstream / centered flag
69	!! (3D array, zind, defined at T-point (0<zind<1))
70	!! Part I : horizontal advection
71	!! * centered flux:
72	!! zcenu = e2u*e3u un mi(ptn)
73	!! zcenv = e1v*e3v vn mj(ptn)
74	!! * upstream flux:
75	!! zupsu = e2u*e3u un (ptb(i) or ptb(i-1) ) [un>0 or <0]
76	!! zupsv = e1v*e3v vn (ptb(j) or ptb(j-1) ) [vn>0 or <0]
77	!! * mixed upstream / centered horizontal advection scheme
78	!! zcofi = max(zind(i+1), zind(i))
79	!! zcofj = max(zind(j+1), zind(j))
80	!! zwx = zcofi * zupsu + (1-zcofi) * zcenu
81	!! zwy = zcofj * zupsv + (1-zcofj) * zcenv
82	!! * horizontal advective trend (divergence of the fluxes)
83	!! ztra = 1/(e1te2te3t) { di-1[zwx] + dj-1[zwy] }
84	!! * Add this trend now to the general trend of tracer (ta,sa):
85	!! pta = pta + ztra
86	!! * trend diagnostic ('key_trdtra' defined): the trend is
87	!! saved for diagnostics. The trends saved is expressed as
88	!! Uh.gradh(T), i.e.
89	!! save trend = ztra + ptn divn
90	!!
91	!! Part II : vertical advection
92	!! For temperature (idem for salinity) the advective trend is com-
93	!! puted as follows :
94	!! ztra = 1/e3t dk+1[ zwz ]
95	!! where the vertical advective flux, zwz, is given by :
96	!! zwz = zcofk * zupst + (1-zcofk) * zcent
97	!! with
98	!! zupsv = upstream flux = wn * (ptb(k) or ptb(k-1) ) [wn>0 or <0]
99	!! zcenu = centered flux = wn * mk(tn)
100	!! The surface boundary condition is :
101	!! variable volume (lk_vvl = T) : zero advective flux
102	!! lin. free-surf (lk_vvl = F) : wn(:,:,1) * ptn(:,:,1)
103	!! Add this trend now to the general trend of tracer (ta,sa):
104	!! pta = pta + ztra
105	!! Trend diagnostic ('key_trdtra' defined): the trend is
106	!! saved for diagnostics. The trends saved is expressed as :
107	!! save trend = w.gradz(T) = ztra - ptn divn.
108	!!
109	!! ** Action : - update pta with the now advective tracer trends
110	!! - save trends if needed
111	!!----------------------------------------------------------------------
112	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
113	USE oce , ONLY: zwx => ua , zwy => va ! (ua,va) used as 3D workspace
114	USE wrk_nemo, ONLY: zwz => wrk_3d_1 , zind => wrk_3d_2 ! 3D workspace
115	USE wrk_nemo, ONLY: ztfreez => wrk_2d_1 ! 2D -
116	!
117	INTEGER , INTENT(in ) :: kt ! ocean time-step index
118	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
119	INTEGER , INTENT(in ) :: kjpt ! number of tracers
120	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
121	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb, ptn ! before and now tracer fields
122	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
123	!
124	INTEGER :: ji, jj, jk, jn ! dummy loop indices
125	INTEGER :: ierr ! local integer
126	REAL(wp) :: zbtr, ztra ! local scalars
127	REAL(wp) :: zfp_ui, zfp_vj, zfp_w, zcofi ! - -
128	REAL(wp) :: zfm_ui, zfm_vj, zfm_w, zcofj, zcofk ! - -
129	REAL(wp) :: zupsut, zcenut, zupst ! - -
130	REAL(wp) :: zupsvt, zcenvt, zcent, zice ! - -
131	!!----------------------------------------------------------------------
132
133	IF( wrk_in_use(2, 1) .OR. wrk_in_use(3, 1,2) ) THEN
134	CALL ctl_stop('tra_adv_cen2: requested workspace arrays unavailable') ; RETURN
135	ENDIF
136
137	IF( kt == nit000 ) THEN
138	IF(lwp) WRITE(numout,*)
139	IF(lwp) WRITE(numout,*) 'tra_adv_cen2 : 2nd order centered advection scheme on ', cdtype
140	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
141	IF(lwp) WRITE(numout,*)
142	!
143	ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
144	IF( ierr /= 0 ) CALL ctl_stop('STOP', 'tra_adv_cen2: unable to allocate array')
145	!
146	upsmsk(:,:) = 0._wp ! not upstream by default
147	!
148	IF( cp_cfg == "orca" ) CALL ups_orca_set ! set mixed Upstream/centered scheme near some straits
149	! ! and in closed seas (orca2 and orca4 only)
150	IF( jp_cfg == 2 .AND. .NOT. ln_rstart ) THEN ! Increase the background in the surface layers
151	avmb(1) = 10. * avmb(1) ; avtb(1) = 10. * avtb(1)
152	avmb(2) = 10. * avmb(2) ; avtb(2) = 10. * avtb(2)
153	avmb(3) = 5. * avmb(3) ; avtb(3) = 5. * avtb(3)
154	avmb(4) = 2.5 * avmb(4) ; avtb(4) = 2.5 * avtb(4)
155	ENDIF
156	!
157	l_trd = .FALSE.
158	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
159	ENDIF
160	!
161	! Upstream / centered scheme indicator
162	! ------------------------------------
163	!!gm not strickly exact : the freezing point should be computed at each ocean levels...
164	!!gm not a big deal since cen2 is no more used in global ice-ocean simulations
165	ztfreez(:,:) = tfreez( tsn(:,:,1,jp_sal) )
166	DO jk = 1, jpk
167	DO jj = 1, jpj
168	DO ji = 1, jpi
169	! ! below ice covered area (if tn < "freezing"+0.1 )
170	IF( tsn(ji,jj,jk,jp_tem) <= ztfreez(ji,jj) + 0.1 ) THEN ; zice = 1.e0
171	ELSE ; zice = 0.e0
172	ENDIF
173	zind(ji,jj,jk) = MAX ( &
174	rnfmsk(ji,jj) * rnfmsk_z(jk), & ! near runoff mouths (& closed sea outflows)
175	upsmsk(ji,jj) , & ! some of some straits
176	zice & ! below ice covered area (if tn < "freezing"+0.1 )
177	& ) * tmask(ji,jj,jk)
178	END DO
179	END DO
180	END DO
181
182	DO jn = 1, kjpt
183	!
184	! I. Horizontal advection
185	! ====================
186	!
187	DO jk = 1, jpkm1
188	! ! Second order centered tracer flux at u- and v-points
189	DO jj = 1, jpjm1
190	!
191	DO ji = 1, fs_jpim1 ! vector opt.
192	! upstream indicator
193	zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
194	zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
195	!
196	! upstream scheme
197	zfp_ui = pun(ji,jj,jk) + ABS( pun(ji,jj,jk) )
198	zfm_ui = pun(ji,jj,jk) - ABS( pun(ji,jj,jk) )
199	zfp_vj = pvn(ji,jj,jk) + ABS( pvn(ji,jj,jk) )
200	zfm_vj = pvn(ji,jj,jk) - ABS( pvn(ji,jj,jk) )
201	zupsut = zfp_ui * ptb(ji,jj,jk,jn) + zfm_ui * ptb(ji+1,jj ,jk,jn)
202	zupsvt = zfp_vj * ptb(ji,jj,jk,jn) + zfm_vj * ptb(ji ,jj+1,jk,jn)
203	! centered scheme
204	zcenut = pun(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji+1,jj ,jk,jn) )
205	zcenvt = pvn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji ,jj+1,jk,jn) )
206	! mixed centered / upstream scheme
207	zwx(ji,jj,jk) = 0.5 * ( zcofi * zupsut + (1.-zcofi) * zcenut )
208	zwy(ji,jj,jk) = 0.5 * ( zcofj * zupsvt + (1.-zcofj) * zcenvt )
209	END DO
210	END DO
211	END DO
212
213	! II. Vertical advection
214	! ==================
215	!
216	! ! Vertical advective fluxes
217	zwz(:,:,jpk) = 0.e0 ! Bottom value : flux set to zero
218	! ! Surface value :
219	IF( lk_vvl ) THEN ; zwz(:,:, 1 ) = 0.e0 ! volume variable
220	ELSE ; zwz(:,:, 1 ) = pwn(:,:,1) * ptn(:,:,1,jn) ! linear free surface
221	ENDIF
222	!
223	DO jk = 2, jpk ! Second order centered tracer flux at w-point
224	DO jj = 2, jpjm1
225	DO ji = fs_2, fs_jpim1 ! vector opt.
226	! upstream indicator
227	zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
228	! mixed centered / upstream scheme
229	zfp_w = pwn(ji,jj,jk) + ABS( pwn(ji,jj,jk) )
230	zfm_w = pwn(ji,jj,jk) - ABS( pwn(ji,jj,jk) )
231	zupst = zfp_w * ptb(ji,jj,jk,jn) + zfm_w * ptb(ji,jj,jk-1,jn)
232	! centered scheme
233	zcent = pwn(ji,jj,jk) * ( ptn(ji,jj,jk,jn) + ptn(ji,jj,jk-1,jn) )
234	! mixed centered / upstream scheme
235	zwz(ji,jj,jk) = 0.5 * ( zcofk * zupst + (1.-zcofk) * zcent )
236	END DO
237	END DO
238	END DO
239
240	! II. Divergence of advective fluxes
241	! ----------------------------------
242	DO jk = 1, jpkm1
243	DO jj = 2, jpjm1
244	DO ji = fs_2, fs_jpim1 ! vector opt.
245	zbtr = 1. / ( e1t(ji,jj) * e2t(ji,jj) * fse3t(ji,jj,jk) )
246	! advective trends
247	ztra = - zbtr * ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
248	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
249	& + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) )
250	! advective trends added to the general tracer trends
251	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) + ztra
252	END DO
253	END DO
254	END DO
255
256	! ! trend diagnostics (contribution of upstream fluxes)
257	IF( l_trd ) THEN
258	CALL trd_tra( kt, cdtype, jn, jptra_trd_xad, zwx, pun, ptn(:,:,:,jn) )
259	CALL trd_tra( kt, cdtype, jn, jptra_trd_yad, zwy, pvn, ptn(:,:,:,jn) )
260	CALL trd_tra( kt, cdtype, jn, jptra_trd_zad, zwz, pwn, ptn(:,:,:,jn) )
261	END IF
262	! ! "Poleward" heat and salt transports (contribution of upstream fluxes)
263	IF( cdtype == 'TRA' .AND. ln_diaptr .AND. ( MOD( kt, nn_fptr ) == 0 ) ) THEN
264	IF( jn == jp_tem ) htr_adv(:) = ptr_vj( zwy(:,:,:) )
265	IF( jn == jp_sal ) str_adv(:) = ptr_vj( zwy(:,:,:) )
266	ENDIF
267	!
268	ENDDO
269
270	! --------------------------- required in restart file to ensure restartability)
271	! avmb, avtb will be read in zdfini in restart case as they are used in zdftke, kpp etc...
272	IF( lrst_oce .AND. cdtype == 'TRA' ) THEN
273	CALL iom_rstput( kt, nitrst, numrow, 'avmb', avmb )
274	CALL iom_rstput( kt, nitrst, numrow, 'avtb', avtb )
275	ENDIF
276	!
277	IF( wrk_not_released(2, 1) .OR. &
278	wrk_not_released(3, 1,2) ) CALL ctl_stop('tra_adv_cen2: failed to release workspace arrays')
279	!
280	END SUBROUTINE tra_adv_cen2
281
282
283	SUBROUTINE ups_orca_set
284	!!----------------------------------------------------------------------
285	!! * ROUTINE ups_orca_set *
286	!!
287	!! ** Purpose : add a portion of upstream scheme in area where the
288	!! centered scheme generates too strong overshoot
289	!!
290	!! ** Method : orca (R4 and R2) confiiguration setting. Set upsmsk
291	!! array to nozero value in some straith.
292	!!
293	!! ** Action : - upsmsk set to 1 at some strait, 0 elsewhere for orca
294	!!----------------------------------------------------------------------
295	INTEGER :: ii0, ii1, ij0, ij1 ! temporary integers
296	!!----------------------------------------------------------------------
297
298	! mixed upstream/centered scheme near river mouths
299	! ------------------------------------------------
300	SELECT CASE ( jp_cfg )
301	! ! =======================
302	CASE ( 4 ) ! ORCA_R4 configuration
303	! ! =======================
304	! ! Gibraltar Strait
305	ii0 = 70 ; ii1 = 71
306	ij0 = 52 ; ij1 = 53 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
307	!
308	! ! =======================
309	CASE ( 2 ) ! ORCA_R2 configuration
310	! ! =======================
311	! ! Gibraltar Strait
312	ij0 = 102 ; ij1 = 102
313	ii0 = 138 ; ii1 = 138 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20
314	ii0 = 139 ; ii1 = 139 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
315	ii0 = 140 ; ii1 = 140 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
316	ij0 = 101 ; ij1 = 102
317	ii0 = 141 ; ii1 = 141 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
318	! ! Bab el Mandeb Strait
319	ij0 = 87 ; ij1 = 88
320	ii0 = 164 ; ii1 = 164 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10
321	ij0 = 88 ; ij1 = 88
322	ii0 = 163 ; ii1 = 163 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
323	ii0 = 162 ; ii1 = 162 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
324	ii0 = 160 ; ii1 = 161 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
325	ij0 = 89 ; ij1 = 89
326	ii0 = 158 ; ii1 = 160 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
327	ij0 = 90 ; ij1 = 90
328	ii0 = 160 ; ii1 = 160 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
329	! ! Sound Strait
330	ij0 = 116 ; ij1 = 116
331	ii0 = 144 ; ii1 = 144 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
332	ii0 = 145 ; ii1 = 147 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
333	ii0 = 148 ; ii1 = 148 ; upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
334	!
335	END SELECT
336
337	! mixed upstream/centered scheme over closed seas
338	! -----------------------------------------------
339	CALL clo_ups( upsmsk(:,:) )
340	!
341	END SUBROUTINE ups_orca_set
342
343	!!======================================================================
344	END MODULE traadv_cen2

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

New URL for NEMO forge! http://forge.nemo-ocean.eu

Context Navigation

source: branches/dev_r2586_dynamic_mem/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_cen2.F90 @ 2696

Download in other formats: