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 branches/2011/dev_UKM0_2011/NEMOGCM/NEMO/OPA_SRC/DYN – NEMO

Context Navigation

source: branches/2011/dev_UKM0_2011/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg_ts.F90 @ 3062

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

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

Download in other formats: