Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

traadv_mus.F90 @ 7910

Last change on this file since 7910 was 7910, checked in by timgraham, 7 years ago
All wrk_alloc removed
Property svn:keywords set to `Id`
File size: 15.2 KB

Line
1	MODULE traadv_mus
2	!!======================================================================
3	!! * MODULE traadv_mus *
4	!! Ocean tracers: horizontal & vertical advective trend
5	!!======================================================================
6	!! History : ! 2000-06 (A.Estublier) for passive tracers
7	!! ! 2001-08 (E.Durand, G.Madec) adapted for T & S
8	!! NEMO 1.0 ! 2002-06 (G. Madec) F90: Free form and module
9	!! 3.2 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
10	!! 3.4 ! 2012-06 (P. Oddo, M. Vichi) include the upstream where needed
11	!! 3.7 ! 2015-09 (G. Madec) add the ice-shelf cavities boundary condition
12	!!----------------------------------------------------------------------
13
14	!!----------------------------------------------------------------------
15	!! tra_adv_mus : update the tracer trend with the horizontal
16	!! and vertical advection trends using MUSCL scheme
17	!!----------------------------------------------------------------------
18	USE oce ! ocean dynamics and active tracers
19	USE trc_oce ! share passive tracers/Ocean variables
20	USE dom_oce ! ocean space and time domain
21	USE trd_oce ! trends: ocean variables
22	USE trdtra ! tracers trends manager
23	USE sbcrnf ! river runoffs
24	USE diaptr ! poleward transport diagnostics
25	USE diaar5 ! AR5 diagnostics
26
27	!
28	USE iom
29	USE timing ! Timing
30	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
31	USE in_out_manager ! I/O manager
32	USE lib_mpp ! distribued memory computing
33	USE lbclnk ! ocean lateral boundary condition (or mpp link)
34
35	IMPLICIT NONE
36	PRIVATE
37
38	PUBLIC tra_adv_mus ! routine called by traadv.F90
39
40	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits
41	! ! and in closed seas (orca 2 and 1 configurations)
42	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xind !: mixed upstream/centered index
43
44	LOGICAL :: l_trd ! flag to compute trends
45	LOGICAL :: l_ptr ! flag to compute poleward transport
46	LOGICAL :: l_hst ! flag to compute heat/salt transport
47
48	!! * Substitutions
49	# include "vectopt_loop_substitute.h90"
50	!!----------------------------------------------------------------------
51	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
52	!! $Id$
53	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
54	!!----------------------------------------------------------------------
55	CONTAINS
56
57	SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pun, pvn, pwn, &
58	& ptb, pta, kjpt, ld_msc_ups )
59	!!----------------------------------------------------------------------
60	!! * ROUTINE tra_adv_mus *
61	!!
62	!! ** Purpose : Compute the now trend due to total advection of tracers
63	!! using a MUSCL scheme (Monotone Upstream-centered Scheme for
64	!! Conservation Laws) and add it to the general tracer trend.
65	!!
66	!! ** Method : MUSCL scheme plus centered scheme at ocean boundaries
67	!! ld_msc_ups=T :
68	!!
69	!! ** Action : - update pta with the now advective tracer trends
70	!! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
71	!! - htr_adv, str_adv : poleward advective heat and salt transport (ln_diaptr=T)
72	!!
73	!! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
74	!! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
75	!!----------------------------------------------------------------------
76	INTEGER , INTENT(in ) :: kt ! ocean time-step index
77	INTEGER , INTENT(in ) :: kit000 ! first time step index
78	CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
79	INTEGER , INTENT(in ) :: kjpt ! number of tracers
80	LOGICAL , INTENT(in ) :: ld_msc_ups ! use upstream scheme within muscl
81	REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
82	REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pun, pvn, pwn ! 3 ocean velocity components
83	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(in ) :: ptb ! before tracer field
84	REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt), INTENT(inout) :: pta ! tracer trend
85	!
86	INTEGER :: ji, jj, jk, jn ! dummy loop indices
87	INTEGER :: ierr ! local integer
88	REAL(wp) :: zu, z0u, zzwx, zw ! local scalars
89	REAL(wp) :: zv, z0v, zzwy, z0w ! - -
90	REAL(wp) :: zalpha ! - -
91	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zslpx, zslpy ! 3D workspace
92	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwx , zwy ! - -
93	!!----------------------------------------------------------------------
94	!
95	IF( nn_timing == 1 ) CALL timing_start('tra_adv_mus')
96	!
97	!
98	IF( kt == kit000 ) THEN
99	IF(lwp) WRITE(numout,*)
100	IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
101	IF(lwp) WRITE(numout,*) ' : mixed up-stream ', ld_msc_ups
102	IF(lwp) WRITE(numout,*) '~~~~~~~'
103	IF(lwp) WRITE(numout,*)
104	!
105	! Upstream / MUSCL scheme indicator
106	!
107	ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
108	xind(:,:,:) = 1._wp ! set equal to 1 where up-stream is not needed
109	!
110	IF( ld_msc_ups ) THEN ! define the upstream indicator (if asked)
111	ALLOCATE( upsmsk(jpi,jpj), STAT=ierr )
112	upsmsk(:,:) = 0._wp ! not upstream by default
113	!
114	DO jk = 1, jpkm1
115	xind(:,:,jk) = 1._wp & ! =>1 where up-stream is not needed
116	& - MAX ( rnfmsk(:,:) * rnfmsk_z(jk), & ! =>0 near runoff mouths (& closed sea outflows)
117	& upsmsk(:,:) ) * tmask(:,:,jk) ! =>0 in some user defined area
118	END DO
119	ENDIF
120	!
121	ENDIF
122	!
123	l_trd = .FALSE.
124	l_hst = .FALSE.
125	l_ptr = .FALSE.
126	IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
127	IF( cdtype == 'TRA' .AND. ln_diaptr ) l_ptr = .TRUE.
128	IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
129	& iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
130	!
131	DO jn = 1, kjpt !== loop over the tracers ==!
132	!
133	! !* Horizontal advective fluxes
134	!
135	! !-- first guess of the slopes
136	zwx(:,:,jpk) = 0._wp ! bottom values
137	zwy(:,:,jpk) = 0._wp
138	DO jk = 1, jpkm1 ! interior values
139	DO jj = 1, jpjm1
140	DO ji = 1, fs_jpim1 ! vector opt.
141	zwx(ji,jj,jk) = umask(ji,jj,jk) * ( ptb(ji+1,jj,jk,jn) - ptb(ji,jj,jk,jn) )
142	zwy(ji,jj,jk) = vmask(ji,jj,jk) * ( ptb(ji,jj+1,jk,jn) - ptb(ji,jj,jk,jn) )
143	END DO
144	END DO
145	END DO
146	CALL lbc_lnk( zwx, 'U', -1. ) ! lateral boundary conditions (changed sign)
147	CALL lbc_lnk( zwy, 'V', -1. )
148	! !-- Slopes of tracer
149	zslpx(:,:,jpk) = 0._wp ! bottom values
150	zslpy(:,:,jpk) = 0._wp
151	DO jk = 1, jpkm1 ! interior values
152	DO jj = 2, jpj
153	DO ji = fs_2, jpi ! vector opt.
154	zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji-1,jj ,jk) ) &
155	& * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji-1,jj ,jk) ) )
156	zslpy(ji,jj,jk) = ( zwy(ji,jj,jk) + zwy(ji ,jj-1,jk) ) &
157	& * ( 0.25 + SIGN( 0.25, zwy(ji,jj,jk) * zwy(ji ,jj-1,jk) ) )
158	END DO
159	END DO
160	END DO
161	!
162	DO jk = 1, jpkm1 !-- Slopes limitation
163	DO jj = 2, jpj
164	DO ji = fs_2, jpi ! vector opt.
165	zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji ,jj,jk) ), &
166	& 2.*ABS( zwx (ji-1,jj,jk) ), &
167	& 2.*ABS( zwx (ji ,jj,jk) ) )
168	zslpy(ji,jj,jk) = SIGN( 1., zslpy(ji,jj,jk) ) * MIN( ABS( zslpy(ji,jj ,jk) ), &
169	& 2.*ABS( zwy (ji,jj-1,jk) ), &
170	& 2.*ABS( zwy (ji,jj ,jk) ) )
171	END DO
172	END DO
173	END DO
174	!
175	DO jk = 1, jpkm1 !-- MUSCL horizontal advective fluxes
176	DO jj = 2, jpjm1
177	DO ji = fs_2, fs_jpim1 ! vector opt.
178	! MUSCL fluxes
179	z0u = SIGN( 0.5, pun(ji,jj,jk) )
180	zalpha = 0.5 - z0u
181	zu = z0u - 0.5 * pun(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u_n(ji,jj,jk)
182	zzwx = ptb(ji+1,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj,jk)
183	zzwy = ptb(ji ,jj,jk,jn) + xind(ji,jj,jk) * zu * zslpx(ji ,jj,jk)
184	zwx(ji,jj,jk) = pun(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
185	!
186	z0v = SIGN( 0.5, pvn(ji,jj,jk) )
187	zalpha = 0.5 - z0v
188	zv = z0v - 0.5 * pvn(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v_n(ji,jj,jk)
189	zzwx = ptb(ji,jj+1,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1,jk)
190	zzwy = ptb(ji,jj ,jk,jn) + xind(ji,jj,jk) * zv * zslpy(ji,jj ,jk)
191	zwy(ji,jj,jk) = pvn(ji,jj,jk) * ( zalpha * zzwx + (1.-zalpha) * zzwy )
192	END DO
193	END DO
194	END DO
195	CALL lbc_lnk( zwx, 'U', -1. ) ; CALL lbc_lnk( zwy, 'V', -1. ) ! lateral boundary conditions (changed sign)
196	!
197	DO jk = 1, jpkm1 !-- Tracer advective trend
198	DO jj = 2, jpjm1
199	DO ji = fs_2, fs_jpim1 ! vector opt.
200	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
201	& + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
202	& * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
203	END DO
204	END DO
205	END DO
206	! ! trend diagnostics
207	IF( l_trd ) THEN
208	CALL trd_tra( kt, cdtype, jn, jptra_xad, zwx, pun, ptb(:,:,:,jn) )
209	CALL trd_tra( kt, cdtype, jn, jptra_yad, zwy, pvn, ptb(:,:,:,jn) )
210	END IF
211	! ! "Poleward" heat and salt transports
212	IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
213	! ! heat transport
214	IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
215	!
216	! !* Vertical advective fluxes
217	!
218	! !-- first guess of the slopes
219	zwx(:,:, 1 ) = 0._wp ! surface & bottom boundary conditions
220	zwx(:,:,jpk) = 0._wp
221	DO jk = 2, jpkm1 ! interior values
222	zwx(:,:,jk) = tmask(:,:,jk) * ( ptb(:,:,jk-1,jn) - ptb(:,:,jk,jn) )
223	END DO
224	! !-- Slopes of tracer
225	zslpx(:,:,1) = 0._wp ! surface values
226	DO jk = 2, jpkm1 ! interior value
227	DO jj = 1, jpj
228	DO ji = 1, jpi
229	zslpx(ji,jj,jk) = ( zwx(ji,jj,jk) + zwx(ji,jj,jk+1) ) &
230	& * ( 0.25 + SIGN( 0.25, zwx(ji,jj,jk) * zwx(ji,jj,jk+1) ) )
231	END DO
232	END DO
233	END DO
234	DO jk = 2, jpkm1 !-- Slopes limitation
235	DO jj = 1, jpj ! interior values
236	DO ji = 1, jpi
237	zslpx(ji,jj,jk) = SIGN( 1., zslpx(ji,jj,jk) ) * MIN( ABS( zslpx(ji,jj,jk ) ), &
238	& 2.*ABS( zwx (ji,jj,jk+1) ), &
239	& 2.*ABS( zwx (ji,jj,jk ) ) )
240	END DO
241	END DO
242	END DO
243	DO jk = 1, jpk-2 !-- vertical advective flux
244	DO jj = 2, jpjm1
245	DO ji = fs_2, fs_jpim1 ! vector opt.
246	z0w = SIGN( 0.5, pwn(ji,jj,jk+1) )
247	zalpha = 0.5 + z0w
248	zw = z0w - 0.5 * pwn(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w_n(ji,jj,jk+1)
249	zzwx = ptb(ji,jj,jk+1,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk+1)
250	zzwy = ptb(ji,jj,jk ,jn) + xind(ji,jj,jk) * zw * zslpx(ji,jj,jk )
251	zwx(ji,jj,jk+1) = pwn(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)
252	END DO
253	END DO
254	END DO
255	IF( ln_linssh ) THEN ! top values, linear free surface only
256	IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean)
257	DO jj = 1, jpj
258	DO ji = 1, jpi
259	zwx(ji,jj, mikt(ji,jj) ) = pwn(ji,jj,mikt(ji,jj)) * ptb(ji,jj,mikt(ji,jj),jn)
260	END DO
261	END DO
262	ELSE ! no cavities: only at the ocean surface
263	zwx(:,:,1) = pwn(:,:,1) * ptb(:,:,1,jn)
264	ENDIF
265	ENDIF
266	!
267	DO jk = 1, jpkm1 !-- vertical advective trend
268	DO jj = 2, jpjm1
269	DO ji = fs_2, fs_jpim1 ! vector opt.
270	pta(ji,jj,jk,jn) = pta(ji,jj,jk,jn) - ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) ) * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
271	END DO
272	END DO
273	END DO
274	! ! send trends for diagnostic
275	IF( l_trd ) CALL trd_tra( kt, cdtype, jn, jptra_zad, zwx, pwn, ptb(:,:,:,jn) )
276	!
277	END DO ! end of tracer loop
278	!
279	!
280	IF( nn_timing == 1 ) CALL timing_stop('tra_adv_mus')
281	!
282	END SUBROUTINE tra_adv_mus
283
284	!!======================================================================
285	END MODULE traadv_mus

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

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

Context Navigation

source: branches/2017/dev_r7881_no_wrk_alloc/NEMOGCM/NEMO/OPA_SRC/TRA/traadv_mus.F90 @ 7910

Download in other formats: