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.

dynspg_ts.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: trunk/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 2724

Last change on this file since 2724 was 2724, checked in by rblod, 13 years ago
Correct dynspg_ts, see ticket #812
Property svn:keywords set to `Id`
File size: 36.7 KB

Line
1	MODULE dynspg_ts
2	!!======================================================================
3	!! History : 1.0 ! 2004-12 (L. Bessieres, G. Madec) Original code
4	!! - ! 2005-11 (V. Garnier, G. Madec) optimization
5	!! - ! 2006-08 (S. Masson) distributed restart using iom
6	!! 2.0 ! 2007-07 (D. Storkey) calls to BDY routines
7	!! - ! 2008-01 (R. Benshila) change averaging method
8	!! 3.2 ! 2009-07 (R. Benshila, G. Madec) Complete revisit associated to vvl reactivation
9	!! 3.3 ! 2010-09 (D. Storkey, E. O'Dea) update for BDY for Shelf configurations
10	!! 3.3 ! 2011-03 (R. Benshila, R. Hordoir, P. Oddo) update calculation of ub_b
11	!!---------------------------------------------------------------------
12	#if defined key_dynspg_ts \|\| defined key_esopa
13	!!----------------------------------------------------------------------
14	!! 'key_dynspg_ts' free surface cst volume with time splitting
15	!!----------------------------------------------------------------------
16	!! dyn_spg_ts : compute surface pressure gradient trend using a time-
17	!! splitting scheme and add to the general trend
18	!! ts_rst : read/write the time-splitting restart fields in the ocean restart file
19	!!----------------------------------------------------------------------
20	USE oce ! ocean dynamics and tracers
21	USE dom_oce ! ocean space and time domain
22	USE sbc_oce ! surface boundary condition: ocean
23	USE dynspg_oce ! surface pressure gradient variables
24	USE phycst ! physical constants
25	USE domvvl ! variable volume
26	USE zdfbfr ! bottom friction
27	USE obcdta ! open boundary condition data
28	USE obcfla ! Flather open boundary condition
29	USE dynvor ! vorticity term
30	USE obc_oce ! Lateral open boundary condition
31	USE obc_par ! open boundary condition parameters
32	USE bdy_oce ! unstructured open boundaries
33	USE bdy_par ! unstructured open boundaries
34	USE bdydta ! unstructured open boundaries
35	USE bdydyn ! unstructured open boundaries
36	USE bdytides ! tidal forcing at unstructured open boundaries.
37	USE lib_mpp ! distributed memory computing library
38	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
39	USE prtctl ! Print control
40	USE in_out_manager ! I/O manager
41	USE iom ! IOM library
42	USE restart ! only for lrst_oce
43	USE zdf_oce
44
45	IMPLICIT NONE
46	PRIVATE
47
48	PUBLIC dyn_spg_ts ! routine called by step.F90
49	PUBLIC ts_rst ! routine called by istate.F90
50	PUBLIC dyn_spg_ts_alloc ! routine called by dynspg.F90
51
52
53	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
54	REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
55
56	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_b, vn_b ! now averaged velocity
57	REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ub_b, vb_b ! before averaged velocity
58
59	!! * Substitutions
60	# include "domzgr_substitute.h90"
61	# include "vectopt_loop_substitute.h90"
62	!!----------------------------------------------------------------------
63	!! NEMO/OPA 4.0 , NEMO Consortium (2011)
64	!! $Id$
65	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
66	!!----------------------------------------------------------------------
67	CONTAINS
68
69	INTEGER FUNCTION dyn_spg_ts_alloc()
70	!!----------------------------------------------------------------------
71	!! * routine dyn_spg_ts_alloc *
72	!!----------------------------------------------------------------------
73	ALLOCATE( ftnw (jpi,jpj) , ftne(jpi,jpj) , un_b(jpi,jpj) , vn_b(jpi,jpj) , &
74	& ftsw (jpi,jpj) , ftse(jpi,jpj) , ub_b(jpi,jpj) , vb_b(jpi,jpj) , STAT= dyn_spg_ts_alloc )
75	!
76	IF( lk_mpp ) CALL mpp_sum( dyn_spg_ts_alloc )
77	IF( dyn_spg_ts_alloc /= 0 ) CALL ctl_warn('dynspg_oce_alloc: failed to allocate arrays')
78	!
79	END FUNCTION dyn_spg_ts_alloc
80
81
82	SUBROUTINE dyn_spg_ts( kt )
83	!!----------------------------------------------------------------------
84	!! * routine dyn_spg_ts *
85	!!
86	!! ** Purpose : Compute the now trend due to the surface pressure
87	!! gradient in case of free surface formulation with time-splitting.
88	!! Add it to the general trend of momentum equation.
89	!!
90	!! ** Method : Free surface formulation with time-splitting
91	!! -1- Save the vertically integrated trend. This general trend is
92	!! held constant over the barotropic integration.
93	!! The Coriolis force is removed from the general trend as the
94	!! surface gradient and the Coriolis force are updated within
95	!! the barotropic integration.
96	!! -2- Barotropic loop : updates of sea surface height (ssha_e) and
97	!! barotropic velocity (ua_e and va_e) through barotropic
98	!! momentum and continuity integration. Barotropic former
99	!! variables are time averaging over the full barotropic cycle
100	!! (= 2 * baroclinic time step) and saved in uX_b
101	!! and vX_b (X specifying after, now or before).
102	!! -3- The new general trend becomes :
103	!! ua = ua - sum_k(ua)/H + ( un_b - ub_b )
104	!!
105	!! ** Action : - Update (ua,va) with the surf. pressure gradient trend
106	!!
107	!! References : Griffies et al., (2003): A technical guide to MOM4. NOAA/GFDL
108	!!---------------------------------------------------------------------
109	USE wrk_nemo, ONLY: wrk_in_use, wrk_not_released
110	USE wrk_nemo, ONLY: zsshun_e => wrk_2d_1 , zsshb_e => wrk_2d_2 , zhdiv => wrk_2d_3
111	USE wrk_nemo, ONLY: zsshvn_e => wrk_2d_4 , zssh_sum => wrk_2d_5
112	USE wrk_nemo, ONLY: zcu => wrk_2d_6 , zwx => wrk_2d_7 , zua => wrk_2d_8 , zbfru => wrk_2d_9
113	USE wrk_nemo, ONLY: zcv => wrk_2d_10 , zwy => wrk_2d_11 , zva => wrk_2d_12 , zbfrv => wrk_2d_13
114	USE wrk_nemo, ONLY: zun => wrk_2d_14 , zun_e => wrk_2d_15 , zub_e => wrk_2d_16 , zu_sum => wrk_2d_17
115	USE wrk_nemo, ONLY: zvn => wrk_2d_18 , zvn_e => wrk_2d_19 , zvb_e => wrk_2d_20 , zv_sum => wrk_2d_21
116	!
117	INTEGER, INTENT(in) :: kt ! ocean time-step index
118	!
119	INTEGER :: ji, jj, jk, jn ! dummy loop indices
120	INTEGER :: icycle ! local scalar
121	REAL(wp) :: zraur, zcoef, z2dt_e, z2dt_b ! local scalars
122	REAL(wp) :: z1_8, zx1, zy1 ! - -
123	REAL(wp) :: z1_4, zx2, zy2 ! - -
124	REAL(wp) :: zu_spg, zu_cor, zu_sld, zu_asp ! - -
125	REAL(wp) :: zv_spg, zv_cor, zv_sld, zv_asp ! - -
126	!!----------------------------------------------------------------------
127
128	IF( wrk_in_use(2, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, &
129	& 11,12,13,14,15,16,17,18,19,20,21 ) ) THEN
130	CALL ctl_stop( 'dyn_spg_ts: requested workspace arrays unavailable' ) ; RETURN
131	ENDIF
132
133	IF( kt == nit000 ) THEN !* initialisation
134	!
135	IF(lwp) WRITE(numout,*)
136	IF(lwp) WRITE(numout,*) 'dyn_spg_ts : surface pressure gradient trend'
137	IF(lwp) WRITE(numout,*) '~~~~~~~~~~ free surface with time splitting'
138	IF(lwp) WRITE(numout,) ' Number of sub cycle in 1 time-step (2 rdt) : icycle = ', 2nn_baro
139	!
140	CALL ts_rst( nit000, 'READ' ) ! read or initialize the following fields: un_b, vn_b
141	!
142	ua_e (:,:) = un_b (:,:)
143	va_e (:,:) = vn_b (:,:)
144	hu_e (:,:) = hu (:,:)
145	hv_e (:,:) = hv (:,:)
146	hur_e (:,:) = hur (:,:)
147	hvr_e (:,:) = hvr (:,:)
148	IF( ln_dynvor_een ) THEN
149	ftne(1,:) = 0.e0 ; ftnw(1,:) = 0.e0 ; ftse(1,:) = 0.e0 ; ftsw(1,:) = 0.e0
150	DO jj = 2, jpj
151	DO ji = fs_2, jpi ! vector opt.
152	ftne(ji,jj) = ( ff(ji-1,jj ) + ff(ji ,jj ) + ff(ji ,jj-1) ) / 3.
153	ftnw(ji,jj) = ( ff(ji-1,jj-1) + ff(ji-1,jj ) + ff(ji ,jj ) ) / 3.
154	ftse(ji,jj) = ( ff(ji ,jj ) + ff(ji ,jj-1) + ff(ji-1,jj-1) ) / 3.
155	ftsw(ji,jj) = ( ff(ji ,jj-1) + ff(ji-1,jj-1) + ff(ji-1,jj ) ) / 3.
156	END DO
157	END DO
158	ENDIF
159	!
160	ENDIF
161
162	! !* Local constant initialization
163	z2dt_b = 2.0 * rdt ! baroclinic time step
164	z1_8 = 0.5 * 0.25 ! coefficient for vorticity estimates
165	z1_4 = 0.5 * 0.5
166	zraur = 1. / rau0 ! 1 / volumic mass
167	!
168	zhdiv(:,:) = 0.e0 ! barotropic divergence
169	zu_sld = 0.e0 ; zu_asp = 0.e0 ! tides trends (lk_tide=F)
170	zv_sld = 0.e0 ; zv_asp = 0.e0
171
172	! -----------------------------------------------------------------------------
173	! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
174	! -----------------------------------------------------------------------------
175	!
176	! !* e3d/dt(Ua), e3Ub, e3*Vn (Vertically integrated)
177	! ! --------------------------
178	zua(:,:) = 0.e0 ; zun(:,:) = 0.e0 ; ub_b(:,:) = 0.e0
179	zva(:,:) = 0.e0 ; zvn(:,:) = 0.e0 ; vb_b(:,:) = 0.e0
180	!
181	DO jk = 1, jpkm1
182	#if defined key_vectopt_loop
183	DO jj = 1, 1 !Vector opt. => forced unrolling
184	DO ji = 1, jpij
185	#else
186	DO jj = 1, jpj
187	DO ji = 1, jpi
188	#endif
189	! ! now trend
190	zua(ji,jj) = zua(ji,jj) + fse3u (ji,jj,jk) * ua(ji,jj,jk) * umask(ji,jj,jk)
191	zva(ji,jj) = zva(ji,jj) + fse3v (ji,jj,jk) * va(ji,jj,jk) * vmask(ji,jj,jk)
192	! ! now velocity
193	zun(ji,jj) = zun(ji,jj) + fse3u (ji,jj,jk) * un(ji,jj,jk)
194	zvn(ji,jj) = zvn(ji,jj) + fse3v (ji,jj,jk) * vn(ji,jj,jk)
195	!
196	#if defined key_vvl
197	ub_b(ji,jj) = ub_b(ji,jj) + (fse3u_0(ji,jj,jk)(1.+sshu_b(ji,jj)muu(ji,jj,jk)))* ub(ji,jj,jk)
198	vb_b(ji,jj) = vb_b(ji,jj) + (fse3v_0(ji,jj,jk)(1.+sshv_b(ji,jj)muv(ji,jj,jk)))* vb(ji,jj,jk)
199	#else
200	ub_b(ji,jj) = ub_b(ji,jj) + fse3u_0(ji,jj,jk) * ub(ji,jj,jk) * umask(ji,jj,jk)
201	vb_b(ji,jj) = vb_b(ji,jj) + fse3v_0(ji,jj,jk) * vb(ji,jj,jk) * vmask(ji,jj,jk)
202	#endif
203	END DO
204	END DO
205	END DO
206
207	! !* baroclinic momentum trend (remove the vertical mean trend)
208	DO jk = 1, jpkm1 ! --------------------------
209	DO jj = 2, jpjm1
210	DO ji = fs_2, fs_jpim1 ! vector opt.
211	ua(ji,jj,jk) = ua(ji,jj,jk) - zua(ji,jj) * hur(ji,jj)
212	va(ji,jj,jk) = va(ji,jj,jk) - zva(ji,jj) * hvr(ji,jj)
213	END DO
214	END DO
215	END DO
216
217	! !* barotropic Coriolis trends * H (vorticity scheme dependent)
218	! ! ---------------------------====
219	zwx(:,:) = zun(:,:) * e2u(:,:) ! now transport
220	zwy(:,:) = zvn(:,:) * e1v(:,:)
221	!
222	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN ! energy conserving or mixed scheme
223	DO jj = 2, jpjm1
224	DO ji = fs_2, fs_jpim1 ! vector opt.
225	zy1 = ( zwy(ji,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
226	zy2 = ( zwy(ji,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
227	zx1 = ( zwx(ji-1,jj) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
228	zx2 = ( zwx(ji ,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
229	! energy conserving formulation for planetary vorticity term
230	zcu(ji,jj) = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 )
231	zcv(ji,jj) =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 )
232	END DO
233	END DO
234	!
235	ELSEIF ( ln_dynvor_ens ) THEN ! enstrophy conserving scheme
236	DO jj = 2, jpjm1
237	DO ji = fs_2, fs_jpim1 ! vector opt.
238	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj)
239	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
240	zcu(ji,jj) = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) )
241	zcv(ji,jj) = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) )
242	END DO
243	END DO
244	!
245	ELSEIF ( ln_dynvor_een ) THEN ! enstrophy and energy conserving scheme
246	DO jj = 2, jpjm1
247	DO ji = fs_2, fs_jpim1 ! vector opt.
248	zcu(ji,jj) = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
249	& + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) )
250	zcv(ji,jj) = - z1_4 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) &
251	& + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) )
252	END DO
253	END DO
254	!
255	ENDIF
256
257	! !* Right-Hand-Side of the barotropic momentum equation
258	! ! ----------------------------------------------------
259	IF( lk_vvl ) THEN ! Variable volume : remove both Coriolis and Surface pressure gradient
260	DO jj = 2, jpjm1
261	DO ji = fs_2, fs_jpim1 ! vector opt.
262	zcu(ji,jj) = zcu(ji,jj) - grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn(ji+1,jj ) &
263	& - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hu(ji,jj) / e1u(ji,jj)
264	zcv(ji,jj) = zcv(ji,jj) - grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn(ji ,jj+1) &
265	& - ( rhd(ji ,jj ,1) + 1 ) * sshn(ji ,jj ) ) * hv(ji,jj) / e2v(ji,jj)
266	END DO
267	END DO
268	ENDIF
269
270	DO jj = 2, jpjm1 ! Remove coriolis term (and possibly spg) from barotropic trend
271	DO ji = fs_2, fs_jpim1
272	zua(ji,jj) = zua(ji,jj) - zcu(ji,jj)
273	zva(ji,jj) = zva(ji,jj) - zcv(ji,jj)
274	END DO
275	END DO
276
277
278	! ! Remove barotropic contribution of bottom friction
279	! ! from the barotropic transport trend
280	zcoef = -1. / z2dt_b
281	# if defined key_vectopt_loop
282	DO jj = 1, 1
283	DO ji = 1, jpij-jpi ! vector opt. (forced unrolling)
284	# else
285	DO jj = 2, jpjm1
286	DO ji = 2, jpim1
287	# endif
288	! Apply stability criteria for bottom friction
289	!RBbug for vvl and external mode we may need to use varying fse3
290	!!gm Rq: the bottom e3 present the smallest variation, the use of e3u_0 is not a big approx.
291	zbfru(ji,jj) = MAX( bfrua(ji,jj) , fse3u(ji,jj,mbku(ji,jj)) * zcoef )
292	zbfrv(ji,jj) = MAX( bfrva(ji,jj) , fse3v(ji,jj,mbkv(ji,jj)) * zcoef )
293	END DO
294	END DO
295
296	IF( lk_vvl ) THEN
297	DO jj = 2, jpjm1
298	DO ji = fs_2, fs_jpim1 ! vector opt.
299	zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) &
300	& / ( hu_0(ji,jj) + sshu_b(ji,jj) + 1.e0 - umask(ji,jj,1) )
301	zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) &
302	& / ( hv_0(ji,jj) + sshv_b(ji,jj) + 1.e0 - vmask(ji,jj,1) )
303	END DO
304	END DO
305	ELSE
306	DO jj = 2, jpjm1
307	DO ji = fs_2, fs_jpim1 ! vector opt.
308	zua(ji,jj) = zua(ji,jj) - zbfru(ji,jj) * ub_b(ji,jj) * hur(ji,jj)
309	zva(ji,jj) = zva(ji,jj) - zbfrv(ji,jj) * vb_b(ji,jj) * hvr(ji,jj)
310	END DO
311	END DO
312	ENDIF
313
314	! !* d/dt(Ua), Ub, Vn (Vertical mean velocity)
315	! ! --------------------------
316	zua(:,:) = zua(:,:) * hur(:,:)
317	zva(:,:) = zva(:,:) * hvr(:,:)
318	!
319	IF( lk_vvl ) THEN
320	ub_b(:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
321	vb_b(:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
322	ELSE
323	ub_b(:,:) = ub_b(:,:) * hur(:,:)
324	vb_b(:,:) = vb_b(:,:) * hvr(:,:)
325	ENDIF
326
327	! -----------------------------------------------------------------------
328	! Phase 2 : Integration of the barotropic equations with time splitting
329	! -----------------------------------------------------------------------
330	!
331	! ! ==================== !
332	! ! Initialisations !
333	! ! ==================== !
334	icycle = 2 * nn_baro ! Number of barotropic sub time-step
335
336	! ! Start from NOW field
337	hu_e (:,:) = hu (:,:) ! ocean depth at u- and v-points
338	hv_e (:,:) = hv (:,:)
339	hur_e (:,:) = hur (:,:) ! ocean depth inverted at u- and v-points
340	hvr_e (:,:) = hvr (:,:)
341	!RBbug zsshb_e(:,:) = sshn (:,:)
342	zsshb_e(:,:) = sshn_b(:,:) ! sea surface height (before and now)
343	sshn_e (:,:) = sshn (:,:)
344
345	zun_e (:,:) = un_b (:,:) ! barotropic velocity (external)
346	zvn_e (:,:) = vn_b (:,:)
347	zub_e (:,:) = un_b (:,:)
348	zvb_e (:,:) = vn_b (:,:)
349
350	zu_sum (:,:) = un_b (:,:) ! summation
351	zv_sum (:,:) = vn_b (:,:)
352	zssh_sum(:,:) = sshn (:,:)
353
354	#if defined key_obc
355	! set ssh corrections to 0
356	! ssh corrections are applied to normal velocities (Flather's algorithm) and averaged over the barotropic loop
357	IF( lp_obc_east ) sshfoe_b(:,:) = 0.e0
358	IF( lp_obc_west ) sshfow_b(:,:) = 0.e0
359	IF( lp_obc_south ) sshfos_b(:,:) = 0.e0
360	IF( lp_obc_north ) sshfon_b(:,:) = 0.e0
361	#endif
362
363	! ! ==================== !
364	DO jn = 1, icycle ! sub-time-step loop ! (from NOW to AFTER+1)
365	! ! ==================== !
366	z2dt_e = 2. * ( rdt / nn_baro )
367	IF( jn == 1 ) z2dt_e = rdt / nn_baro
368
369	! !* Update the forcing (OBC, BDY and tides)
370	! ! ------------------
371	IF( lk_obc ) CALL obc_dta_bt ( kt, jn )
372	IF( lk_bdy ) CALL bdy_dta_fla( kt, jn+1, icycle )
373
374	! !* after ssh_e
375	! ! -----------
376	DO jj = 2, jpjm1 ! Horizontal divergence of barotropic transports
377	DO ji = fs_2, fs_jpim1 ! vector opt.
378	zhdiv(ji,jj) = ( e2u(ji ,jj) * zun_e(ji ,jj) * hu_e(ji ,jj) &
379	& - e2u(ji-1,jj) * zun_e(ji-1,jj) * hu_e(ji-1,jj) &
380	& + e1v(ji,jj ) * zvn_e(ji,jj ) * hv_e(ji,jj ) &
381	& - e1v(ji,jj-1) * zvn_e(ji,jj-1) * hv_e(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) )
382	END DO
383	END DO
384	!
385	#if defined key_obc
386	! ! OBC : zhdiv must be zero behind the open boundary
387	!! mpp remark: The zeroing of hdiv can probably be extended to 1->jpi/jpj for the correct row/column
388	IF( lp_obc_east ) zhdiv(nie0p1:nie1p1,nje0 :nje1 ) = 0.e0 ! east
389	IF( lp_obc_west ) zhdiv(niw0 :niw1 ,njw0 :njw1 ) = 0.e0 ! west
390	IF( lp_obc_north ) zhdiv(nin0 :nin1 ,njn0p1:njn1p1) = 0.e0 ! north
391	IF( lp_obc_south ) zhdiv(nis0 :nis1 ,njs0 :njs1 ) = 0.e0 ! south
392	#endif
393	#if defined key_bdy
394	zhdiv(:,:) = zhdiv(:,:) * bdytmask(:,:) ! BDY mask
395	#endif
396	!
397	DO jj = 2, jpjm1 ! leap-frog on ssh_e
398	DO ji = fs_2, fs_jpim1 ! vector opt.
399	ssha_e(ji,jj) = ( zsshb_e(ji,jj) - z2dt_e * ( zraur * ( emp(ji,jj)-rnf(ji,jj) ) + zhdiv(ji,jj) ) ) * tmask(ji,jj,1)
400	END DO
401	END DO
402
403	! !* after barotropic velocities (vorticity scheme dependent)
404	! ! ---------------------------
405	zwx(:,:) = e2u(:,:) * zun_e(:,:) * hu_e(:,:) ! now_e transport
406	zwy(:,:) = e1v(:,:) * zvn_e(:,:) * hv_e(:,:)
407	!
408	IF( ln_dynvor_ene .OR. ln_dynvor_mix ) THEN !== energy conserving or mixed scheme ==!
409	DO jj = 2, jpjm1
410	DO ji = fs_2, fs_jpim1 ! vector opt.
411	! surface pressure gradient
412	IF( lk_vvl) THEN
413	zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) &
414	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj)
415	zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) &
416	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj)
417	ELSE
418	zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
419	zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
420	ENDIF
421	! energy conserving formulation for planetary vorticity term
422	zy1 = ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) ) / e1u(ji,jj)
423	zy2 = ( zwy(ji ,jj ) + zwy(ji+1,jj ) ) / e1u(ji,jj)
424	zx1 = ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) ) / e2v(ji,jj)
425	zx2 = ( zwx(ji ,jj ) + zwx(ji ,jj+1) ) / e2v(ji,jj)
426	zu_cor = z1_4 * ( ff(ji ,jj-1) * zy1 + ff(ji,jj) * zy2 ) * hur_e(ji,jj)
427	zv_cor =-z1_4 * ( ff(ji-1,jj ) * zx1 + ff(ji,jj) * zx2 ) * hvr_e(ji,jj)
428	! after velocities with implicit bottom friction
429	ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) &
430	& / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
431	va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) &
432	& / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
433	END DO
434	END DO
435	!
436	ELSEIF ( ln_dynvor_ens ) THEN !== enstrophy conserving scheme ==!
437	DO jj = 2, jpjm1
438	DO ji = fs_2, fs_jpim1 ! vector opt.
439	! surface pressure gradient
440	IF( lk_vvl) THEN
441	zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) &
442	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj)
443	zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) &
444	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj)
445	ELSE
446	zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
447	zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
448	ENDIF
449	! enstrophy conserving formulation for planetary vorticity term
450	zy1 = z1_8 * ( zwy(ji ,jj-1) + zwy(ji+1,jj-1) + zwy(ji,jj) + zwy(ji+1,jj ) ) / e1u(ji,jj)
451	zx1 = - z1_8 * ( zwx(ji-1,jj ) + zwx(ji-1,jj+1) + zwx(ji,jj) + zwx(ji ,jj+1) ) / e2v(ji,jj)
452	zu_cor = zy1 * ( ff(ji ,jj-1) + ff(ji,jj) ) * hur_e(ji,jj)
453	zv_cor = zx1 * ( ff(ji-1,jj ) + ff(ji,jj) ) * hvr_e(ji,jj)
454	! after velocities with implicit bottom friction
455	ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) &
456	& / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
457	va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) &
458	& / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
459	END DO
460	END DO
461	!
462	ELSEIF ( ln_dynvor_een ) THEN !== energy and enstrophy conserving scheme ==!
463	DO jj = 2, jpjm1
464	DO ji = fs_2, fs_jpim1 ! vector opt.
465	! surface pressure gradient
466	IF( lk_vvl) THEN
467	zu_spg = -grav * ( ( rhd(ji+1,jj ,1) + 1 ) * sshn_e(ji+1,jj ) &
468	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e1u(ji,jj)
469	zv_spg = -grav * ( ( rhd(ji ,jj+1,1) + 1 ) * sshn_e(ji ,jj+1) &
470	& - ( rhd(ji ,jj ,1) + 1 ) * sshn_e(ji ,jj ) ) / e2v(ji,jj)
471	ELSE
472	zu_spg = -grav * ( sshn_e(ji+1,jj) - sshn_e(ji,jj) ) / e1u(ji,jj)
473	zv_spg = -grav * ( sshn_e(ji,jj+1) - sshn_e(ji,jj) ) / e2v(ji,jj)
474	ENDIF
475	! energy/enstrophy conserving formulation for planetary vorticity term
476	zu_cor = + z1_4 / e1u(ji,jj) * ( ftne(ji,jj ) * zwy(ji ,jj ) + ftnw(ji+1,jj) * zwy(ji+1,jj ) &
477	& + ftse(ji,jj ) * zwy(ji ,jj-1) + ftsw(ji+1,jj) * zwy(ji+1,jj-1) ) * hur_e(ji,jj)
478	zv_cor = - z1_4 / e2v(ji,jj) * ( ftsw(ji,jj+1) * zwx(ji-1,jj+1) + ftse(ji,jj+1) * zwx(ji ,jj+1) &
479	& + ftnw(ji,jj ) * zwx(ji-1,jj ) + ftne(ji,jj ) * zwx(ji ,jj ) ) * hvr_e(ji,jj)
480	! after velocities with implicit bottom friction
481	ua_e(ji,jj) = ( zub_e(ji,jj) + z2dt_e * ( zu_cor + zu_spg + zu_sld + zu_asp + zua(ji,jj) ) ) * umask(ji,jj,1) &
482	& / ( 1.e0 - z2dt_e * bfrua(ji,jj) * hur_e(ji,jj) )
483	va_e(ji,jj) = ( zvb_e(ji,jj) + z2dt_e * ( zv_cor + zv_spg + zv_sld + zv_asp + zva(ji,jj) ) ) * vmask(ji,jj,1) &
484	& / ( 1.e0 - z2dt_e * bfrva(ji,jj) * hvr_e(ji,jj) )
485	END DO
486	END DO
487	!
488	ENDIF
489	! !* domain lateral boundary
490	! ! -----------------------
491	! ! Flather's boundary condition for the barotropic loop :
492	! ! - Update sea surface height on each open boundary
493	! ! - Correct the velocity
494
495	IF( lk_obc ) CALL obc_fla_ts ( ua_e, va_e, sshn_e, ssha_e )
496	IF( lk_bdy .OR. ln_tides ) CALL bdy_dyn_fla( sshn_e )
497	!
498	CALL lbc_lnk( ua_e , 'U', -1. ) ! local domain boundaries
499	CALL lbc_lnk( va_e , 'V', -1. )
500	CALL lbc_lnk( ssha_e, 'T', 1. )
501
502	zu_sum (:,:) = zu_sum (:,:) + ua_e (:,:) ! Sum over sub-time-steps
503	zv_sum (:,:) = zv_sum (:,:) + va_e (:,:)
504	zssh_sum(:,:) = zssh_sum(:,:) + ssha_e(:,:)
505
506	! !* Time filter and swap
507	! ! --------------------
508	IF( jn == 1 ) THEN ! Swap only (1st Euler time step)
509	zsshb_e(:,:) = sshn_e(:,:)
510	zub_e (:,:) = zun_e (:,:)
511	zvb_e (:,:) = zvn_e (:,:)
512	sshn_e (:,:) = ssha_e(:,:)
513	zun_e (:,:) = ua_e (:,:)
514	zvn_e (:,:) = va_e (:,:)
515	ELSE ! Swap + Filter
516	zsshb_e(:,:) = atfp * ( zsshb_e(:,:) + ssha_e(:,:) ) + atfp1 * sshn_e(:,:)
517	zub_e (:,:) = atfp * ( zub_e (:,:) + ua_e (:,:) ) + atfp1 * zun_e (:,:)
518	zvb_e (:,:) = atfp * ( zvb_e (:,:) + va_e (:,:) ) + atfp1 * zvn_e (:,:)
519	sshn_e (:,:) = ssha_e(:,:)
520	zun_e (:,:) = ua_e (:,:)
521	zvn_e (:,:) = va_e (:,:)
522	ENDIF
523
524	IF( lk_vvl ) THEN !* Update ocean depth (variable volume case only)
525	! ! ------------------
526	DO jj = 1, jpjm1 ! Sea Surface Height at u- & v-points
527	DO ji = 1, fs_jpim1 ! Vector opt.
528	zsshun_e(ji,jj) = 0.5 * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) ) &
529	& * ( e1t(ji ,jj) * e2t(ji ,jj) * sshn_e(ji ,jj) &
530	& + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn_e(ji+1,jj) )
531	zsshvn_e(ji,jj) = 0.5 * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) ) &
532	& * ( e1t(ji,jj ) * e2t(ji,jj ) * sshn_e(ji,jj ) &
533	& + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn_e(ji,jj+1) )
534	END DO
535	END DO
536	CALL lbc_lnk( zsshun_e, 'U', 1. ) ! lateral boundaries conditions
537	CALL lbc_lnk( zsshvn_e, 'V', 1. )
538	!
539	hu_e (:,:) = hu_0(:,:) + zsshun_e(:,:) ! Ocean depth at U- and V-points
540	hv_e (:,:) = hv_0(:,:) + zsshvn_e(:,:)
541	hur_e(:,:) = umask(:,:,1) / ( hu_e(:,:) + 1.e0 - umask(:,:,1) )
542	hvr_e(:,:) = vmask(:,:,1) / ( hv_e(:,:) + 1.e0 - vmask(:,:,1) )
543	!
544	ENDIF
545	! ! ==================== !
546	END DO ! end loop !
547	! ! ==================== !
548
549	#if defined key_obc
550	IF( lp_obc_east ) sshfoe_b(:,:) = zcoef * sshfoe_b(:,:) !!gm totally useless ?????
551	IF( lp_obc_west ) sshfow_b(:,:) = zcoef * sshfow_b(:,:)
552	IF( lp_obc_north ) sshfon_b(:,:) = zcoef * sshfon_b(:,:)
553	IF( lp_obc_south ) sshfos_b(:,:) = zcoef * sshfos_b(:,:)
554	#endif
555
556	! -----------------------------------------------------------------------------
557	! Phase 3. update the general trend with the barotropic trend
558	! -----------------------------------------------------------------------------
559	!
560	! !* Time average ==> after barotropic u, v, ssh
561	zcoef = 1.e0 / ( 2 * nn_baro + 1 )
562	zu_sum(:,:) = zcoef * zu_sum (:,:)
563	zv_sum(:,:) = zcoef * zv_sum (:,:)
564	!
565	! !* update the general momentum trend
566	DO jk=1,jpkm1
567	ua(:,:,jk) = ua(:,:,jk) + ( zu_sum(:,:) - ub_b(:,:) ) / z2dt_b
568	va(:,:,jk) = va(:,:,jk) + ( zv_sum(:,:) - vb_b(:,:) ) / z2dt_b
569	END DO
570	un_b (:,:) = zu_sum(:,:)
571	vn_b (:,:) = zv_sum(:,:)
572	sshn_b(:,:) = zcoef * zssh_sum(:,:)
573	!
574	! !* write time-spliting arrays in the restart
575	IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' )
576	!
577	!
578	IF( wrk_not_released(2, 1, 2, 3, 4, 5, 6, 7, 8, 9,10, &
579	& 11,12,13,14,15,16,17,18,19,20,21) ) &
580	CALL ctl_stop('dyn_spg_ts: failed to release workspace arrays')
581	!
582	END SUBROUTINE dyn_spg_ts
583
584
585	SUBROUTINE ts_rst( kt, cdrw )
586	!!---------------------------------------------------------------------
587	!! * ROUTINE ts_rst *
588	!!
589	!! ** Purpose : Read or write time-splitting arrays in restart file
590	!!----------------------------------------------------------------------
591	INTEGER , INTENT(in) :: kt ! ocean time-step
592	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
593	!
594	INTEGER :: ji, jk ! dummy loop indices
595	!!----------------------------------------------------------------------
596	!
597	IF( TRIM(cdrw) == 'READ' ) THEN
598	IF( iom_varid( numror, 'un_b', ldstop = .FALSE. ) > 0 ) THEN
599	CALL iom_get( numror, jpdom_autoglo, 'un_b' , un_b (:,:) ) ! external velocity issued
600	CALL iom_get( numror, jpdom_autoglo, 'vn_b' , vn_b (:,:) ) ! from barotropic loop
601	ELSE
602	un_b (:,:) = 0.e0
603	vn_b (:,:) = 0.e0
604	! vertical sum
605	IF( lk_vopt_loop ) THEN ! vector opt., forced unroll
606	DO jk = 1, jpkm1
607	DO ji = 1, jpij
608	un_b(ji,1) = un_b(ji,1) + fse3u(ji,1,jk) * un(ji,1,jk)
609	vn_b(ji,1) = vn_b(ji,1) + fse3v(ji,1,jk) * vn(ji,1,jk)
610	END DO
611	END DO
612	ELSE ! No vector opt.
613	DO jk = 1, jpkm1
614	un_b(:,:) = un_b(:,:) + fse3u(:,:,jk) * un(:,:,jk)
615	vn_b(:,:) = vn_b(:,:) + fse3v(:,:,jk) * vn(:,:,jk)
616	END DO
617	ENDIF
618	un_b (:,:) = un_b(:,:) * hur(:,:)
619	vn_b (:,:) = vn_b(:,:) * hvr(:,:)
620	ENDIF
621
622	! Vertically integrated velocity (before)
623	IF (neuler/=0) THEN
624	ub_b (:,:) = 0.e0
625	vb_b (:,:) = 0.e0
626
627	! vertical sum
628	IF( lk_vopt_loop ) THEN ! vector opt., forced unroll
629	DO jk = 1, jpkm1
630	DO ji = 1, jpij
631	ub_b(ji,1) = ub_b(ji,1) + fse3u_b(ji,1,jk) * ub(ji,1,jk)
632	vb_b(ji,1) = vb_b(ji,1) + fse3v_b(ji,1,jk) * vb(ji,1,jk)
633	END DO
634	END DO
635	ELSE ! No vector opt.
636	DO jk = 1, jpkm1
637	ub_b(:,:) = ub_b(:,:) + fse3u_b(:,:,jk) * ub(:,:,jk)
638	vb_b(:,:) = vb_b(:,:) + fse3v_b(:,:,jk) * vb(:,:,jk)
639	END DO
640	ENDIF
641
642	IF( lk_vvl ) THEN
643	ub_b (:,:) = ub_b(:,:) * umask(:,:,1) / ( hu_0(:,:) + sshu_b(:,:) + 1.e0 - umask(:,:,1) )
644	vb_b (:,:) = vb_b(:,:) * vmask(:,:,1) / ( hv_0(:,:) + sshv_b(:,:) + 1.e0 - vmask(:,:,1) )
645	ELSE
646	ub_b(:,:) = ub_b(:,:) * hur(:,:)
647	vb_b(:,:) = vb_b(:,:) * hvr(:,:)
648	ENDIF
649	ELSE ! neuler==0
650	ub_b (:,:) = un_b (:,:)
651	vb_b (:,:) = vn_b (:,:)
652	ENDIF
653
654	IF( iom_varid( numror, 'sshn_b', ldstop = .FALSE. ) > 0 ) THEN
655	CALL iom_get( numror, jpdom_autoglo, 'sshn_b' , sshn_b (:,:) ) ! filtered ssh
656	ELSE
657	sshn_b(:,:) = sshb(:,:) ! if not in restart set previous time mean to current baroclinic before value
658	ENDIF
659	ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN
660	CALL iom_rstput( kt, nitrst, numrow, 'un_b' , un_b (:,:) ) ! external velocity and ssh
661	CALL iom_rstput( kt, nitrst, numrow, 'vn_b' , vn_b (:,:) ) ! issued from barotropic loop
662	CALL iom_rstput( kt, nitrst, numrow, 'sshn_b' , sshn_b(:,:) ) !
663	ENDIF
664	!
665	END SUBROUTINE ts_rst
666
667	#else
668	!!----------------------------------------------------------------------
669	!! Default case : Empty module No standart free surface cst volume
670	!!----------------------------------------------------------------------
671	CONTAINS
672	INTEGER FUNCTION dyn_spg_ts_alloc() ! Dummy function
673	dyn_spg_ts_alloc = 0
674	END FUNCTION dyn_spg_ts_alloc
675	SUBROUTINE dyn_spg_ts( kt ) ! Empty routine
676	INTEGER, INTENT(in) :: kt
677	WRITE(,) 'dyn_spg_ts: You should not have seen this print! error?', kt
678	END SUBROUTINE dyn_spg_ts
679	SUBROUTINE ts_rst( kt, cdrw ) ! Empty routine
680	INTEGER , INTENT(in) :: kt ! ocean time-step
681	CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag
682	WRITE(,) 'ts_rst : You should not have seen this print! error?', kt, cdrw
683	END SUBROUTINE ts_rst
684	#endif
685
686	!!======================================================================
687	END MODULE dynspg_ts

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

Download in other formats: